forked from metno/NORTRIP_multiroad
-
Notifications
You must be signed in to change notification settings - Fork 0
/
NORTRIP_multiroad_read_meteo_netcdf3.f90
142 lines (114 loc) · 6.8 KB
/
NORTRIP_multiroad_read_meteo_netcdf3.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
subroutine NORTRIP_read_meteo_netcdf3
use NORTRIP_multiroad_index_definitions
implicit none
include 'netcdf.inc'
!Local variables
integer status_nc !Error message
integer id_nc
integer dim_id_nc(num_dims_nc)
integer xtype_nc(num_var_nc)
integer natts_nc(num_var_nc)
integer var_id_nc(num_var_nc)
character(256) dimname_temp
integer i
integer i_grid_mid,j_grid_mid
real dlat_nc
integer exists
double precision, allocatable :: var2d_nc_dp(:,:)
double precision, allocatable :: var3d_nc_dp(:,:,:)
write(unit_logfile,'(A)') '================================================================'
write(unit_logfile,'(A)') 'Reading meteorological data (NORTRIP_read_meteo_netcdf)'
write(unit_logfile,'(A)') '================================================================'
!pathname_nc='C:\BEDRE BYLUFT\NORTRIP implementation\test\';
!filename_nc='AROME_1KM_OSLO_20141028_EPI.nc'
pathfilename_nc=trim(pathname_nc)//trim(filename_nc)
!Test existence of the filename. If does not exist then use default
inquire(file=trim(pathfilename_nc),exist=exists)
if (.not.exists) then
write(unit_logfile,'(A,A)') ' ERROR: Meteo netcdf file does not exist: ', trim(pathfilename_nc)
write(unit_logfile,'(A)') ' STOPPING'
!write(*,'(A,A)') ' ERROR: Meteo netcdf file does not exist. Stopping: ', trim(pathfilename_nc)
stop
endif
!Open the netcdf file for reading
write(unit_logfile,'(2A)') ' Opening netcdf meteo file: ',trim(pathfilename_nc)
status_nc = NF_OPEN (pathfilename_nc, NF_NOWRITE, id_nc)
if (status_nc .NE. NF_NOERR) write(unit_logfile,'(A,I)') 'ERROR opening netcdf file: ',status_nc
!Find out the x,y and time dimmensions of the file by looking at pressure variable
status_nc = NF_INQ_DIMID (id_nc,dim_name_nc(x_index),dim_id_nc(x_index))
status_nc = NF_INQ_DIM (id_nc,dim_id_nc(x_index),dimname_temp,dim_length_nc(x_index))
status_nc = NF_INQ_DIMID (id_nc,dim_name_nc(y_index),dim_id_nc(y_index))
status_nc = NF_INQ_DIM (id_nc,dim_id_nc(y_index),dimname_temp,dim_length_nc(y_index))
status_nc = NF_INQ_DIMID (id_nc,dim_name_nc(time_index),dim_id_nc(time_index))
status_nc = NF_INQ_DIM (id_nc,dim_id_nc(time_index),dimname_temp,dim_length_nc(time_index))
write(unit_logfile,'(A,3I)') ' Size of dimensions (x,y,t): ',dim_length_nc
!Allocate the nc arrays for reading
allocate (var1d_nc(num_dims_nc,maxval(dim_length_nc))) !x and y and time maximum dimmensions
allocate (var3d_nc(num_var_nc,dim_length_nc(x_index),dim_length_nc(y_index),dim_length_nc(time_index)))
allocate (var2d_nc(2,dim_length_nc(x_index),dim_length_nc(y_index))) !Lat and lon
allocate (var3d_nc_dp(dim_length_nc(x_index),dim_length_nc(y_index),dim_length_nc(time_index)))
allocate (var2d_nc_dp(dim_length_nc(x_index),dim_length_nc(y_index))) !Lat and lon
!Set the number of hours to be read
!Read the x, y and time values
do i=1,num_dims_nc
status_nc = NF_INQ_VARID (id_nc, trim(dim_name_nc(i)), var_id_nc(i))
status_nc = NF_GET_VARA_REAL (id_nc, var_id_nc(i), dim_start_nc(i), dim_length_nc(i), var1d_nc(i,:))
if (i.eq.time_index) then
write(unit_logfile,'(3A,2i12)') ' ',trim(dim_name_nc(i)),' (min, max in hours): ' &
,minval(int((var1d_nc(i,1:dim_length_nc(i))-var1d_nc(i,dim_start_nc(i)))/3600.+.5)+1) &
,maxval(int((var1d_nc(i,1:dim_length_nc(i))-var1d_nc(i,dim_start_nc(i)))/3600.+.5)+1)
else
write(unit_logfile,'(3A,2f12.2)') ' ',trim(dim_name_nc(i)),' (min, max in km): ' &
,minval(var1d_nc(i,1:dim_length_nc(i))/1000.),maxval(var1d_nc(i,1:dim_length_nc(i))/1000.)
endif
enddo
!Read through the variables in a loop
do i=1,num_var_nc
!write(*,*) i,trim(var_name_nc(i))
status_nc = NF_INQ_VARID (id_nc, trim(var_name_nc(i)), var_id_nc(i))
!write(*,*) 'Status1: ',status_nc,var_id_nc(i),trim(var_name_nc(i))
!write(*,*) 'Status1: ',dim_start_nc
!write(*,*) 'Status1: ',dim_length_nc
if (status_nc.ge.0) then
if (i.eq.lat_index.or.i.eq.lon_index) then
if (index(meteo_data_type,'arome').gt.0) then
status_nc = NF_GET_VARA_DOUBLE (id_nc, var_id_nc(i), dim_start_nc(1:2), dim_length_nc(1:2), var2d_nc_dp);var2d_nc(i,:,:)=real(var2d_nc_dp)
else
status_nc = NF_GET_VARA_REAL (id_nc, var_id_nc(i), dim_start_nc(1:2), dim_length_nc(1:2), var2d_nc(i,:,:))
endif
write(unit_logfile,'(3A,2f16.4)') ' ',trim(var_name_nc(i)),' (min, max): ' &
,minval(var2d_nc(i,:,:)),maxval(var2d_nc(i,:,:))
else
if (index(meteo_data_type,'arome').gt.0) then
status_nc = NF_GET_VARA_DOUBLE (id_nc, var_id_nc(i), dim_start_nc, dim_length_nc, var3d_nc_dp);var3d_nc(i,:,:,:)=real(var3d_nc_dp)
else
status_nc = NF_GET_VARA_REAL (id_nc, var_id_nc(i), dim_start_nc, dim_length_nc, var3d_nc(i,:,:,:))
endif
write(unit_logfile,'(3A,2f16.2)') ' ',trim(var_name_nc(i)),' (min, max): ' &
,minval(var3d_nc(i,:,:,:)),maxval(var3d_nc(i,:,:,:))
endif
var_available_nc(i)=.true.
else
write(unit_logfile,'(8A,8A)') ' Cannot read ',trim(var_name_nc(i))
var_available_nc(i)=.false.
endif
enddo
!NOTE: round off errors in precipitation. Need to include a 0 minimum.
status_nc = NF_CLOSE (id_nc)
!Calculate angle difference between North and the Model Y direction based on the middle grids
i_grid_mid=int(dim_length_nc(x_index)/2)
j_grid_mid=int(dim_length_nc(y_index)/2)
dgrid_nc(x_index)=var1d_nc(x_index,i_grid_mid)-var1d_nc(x_index,i_grid_mid-1)
dgrid_nc(y_index)=var1d_nc(y_index,j_grid_mid)-var1d_nc(y_index,j_grid_mid-1)
dlat_nc=var2d_nc(lat_index,i_grid_mid,j_grid_mid)-var2d_nc(lat_index,i_grid_mid,j_grid_mid-1)
!If the coordinates are in km instead of metres then change to metres (assuming the difference is not going to be > 100 km
if (dgrid_nc(x_index).lt.100) dgrid_nc=dgrid_nc*1000.
angle_nc=180./3.14159*acos(dlat_nc*3.14159/180.*6.37e6/dgrid_nc(x_index))
write(unit_logfile,'(A,2f12.1)') ' Grid spacing X and Y (m): ', dgrid_nc(x_index),dgrid_nc(y_index)
write(unit_logfile,'(A,2i,f12.4)') ' Angle difference between grid and geo North (i,j,deg): ', i_grid_mid,j_grid_mid,angle_nc
!Set the array dimensions to the available ones. Can be changed later based on input information, particularly for time
end_dim_nc=dim_length_nc
start_dim_nc=dim_start_nc
deallocate (var3d_nc_dp)
deallocate (var2d_nc_dp)
end subroutine NORTRIP_read_meteo_netcdf3