MI - Fimex
fortran_test.f90

Example on using the fimex interface

1 
3 PROGRAM fortran_test
4  USE fimex, ONLY : fimexio, set_filetype, axis_geox, axis_geoy, axis_lon, axis_lat,interpol_bilinear,&
5  filetype_rw
6  IMPLICIT NONE
7  TYPE(fimexio) :: fio, finter, finter2, frw
8  INTEGER :: ierr,i
9  CHARACTER(LEN=80) :: input_file
10  CHARACTER(LEN=80) :: config_file
11  CHARACTER(LEN=80) :: varname
12  REAL(KIND=8),DIMENSION(:),ALLOCATABLE,TARGET :: field
13  REAL(KIND=8),DIMENSION(:,:,:,:),POINTER :: field4d
14  REAL(KIND=8),DIMENSION(:,:,:),POINTER :: field3d
15  REAL(KIND=8), DIMENSION(:),TARGET,ALLOCATABLE :: lonvals
16  REAL(KIND=8), DIMENSION(:),TARGET,ALLOCATABLE :: latvals
17  INTEGER :: dataread
18  INTEGER :: nx,ny
19  INTEGER :: ndims
20  INTEGER(KIND=4), ALLOCATABLE, DIMENSION(:) :: start, vsize, atypes
21  INTEGER(KIND=8) :: vars
22  CHARACTER(LEN=10) :: cunit,cfiletype
23  CHARACTER(LEN=1024) :: dimname
24 
25 
26  IF (( command_argument_count() /= 4 ) .AND. ( command_argument_count() /= 5 )) THEN
27  WRITE(*,*) "Usage: ./fortran_test input-file file-type variable unit [config]"
28  WRITE(*,*) "Example ./fortran_test /opdata/arome_norway25/AROME_Norway25_00.nc netcdf surface_air_pressure hPa"
29  WRITE(*,*) " ./fortran_test /opdata/arome_norway25/preprod/AROME_Norway25_00.dat felt surface_air_pressure "//&
30  "hPa ~/metop/arome_norway25/etc/AromeFeltReaderConfig.xml"
31  ELSE
32  CALL get_command_argument(1,input_file)
33  CALL get_command_argument(2,cfiletype)
34  CALL get_command_argument(3,varname)
35  CALL get_command_argument(4,cunit)
36  IF ( command_argument_count() == 5 ) THEN
37  CALL get_command_argument(5,config_file)
38  ELSE
39  config_file = ""
40  ENDIF
41 
42  ! Open file
43  ierr=fio%open(input_file,config_file,set_filetype(cfiletype))
44  IF ( ierr /= 0 ) CALL error("Can't make io-object with file:"//trim(input_file)//" config: "//config_file)
45  WRITE(0,*) "open_file: success"
46 
47  write(*,*) "refTime = ", fio%get_refTime('days since 2013-01-01 00:00:00 +0000'), " days since 2013-01-01"
48 
49  write(*,*) "unlimited dimension = ", trim(fio%file_ulim_dimname())
50 
51  write(*,*) "dimensions"
52  DO i = 1, fio%dimensions_size()
53  write(*,*) i, " ", trim(fio%file_dimname(i)), fio%get_dimsize(fio%file_dimname(i))
54  END DO
55 
56  write(*,*) "variables"
57  DO i = 1, fio%variables_size()
58  write(*,*) i, " ", trim(fio%get_varname(i))," (", trim(fio%get_var_longitude(fio%get_varname(i))), ",",&
59  trim(fio%get_var_latitude(fio%get_varname(i))), ")", &
60  fio%get_vartype(fio%get_varname(i))
61  END DO
62 
63  ! Get dimensions
64  ndims=fio%get_dimensions(varname)
65  IF ( ndims <= 0 ) CALL error("Can't make slicebuilder")
66  WRITE(0,*) "get_dimensions: ", ndims
67 
68  ALLOCATE(start(ndims))
69  ALLOCATE(vsize(ndims))
70  ierr = fio%get_dimension_start_size(start, vsize)
71  DO i = 1, ndims
72  WRITE(*,*) "dimname ", i, ": ", trim(fio%get_dimname(i)), vsize(i)
73  END DO
74  ALLOCATE(atypes(ndims))
75  ierr = fio%get_axistypes(atypes)
76 
77  DO i = 1, ndims
78  !WRITE (*,*) i, " axistype: ", atypes(i)
79  !WRITE (*,*) AXIS_GeoX, AXIS_GeoY, AXIS_Lon, AXIS_Lat
80  SELECT CASE (atypes(i))
81  CASE(axis_geox, axis_lon)
82  nx = vsize(i)
83  CASE(axis_geoy, axis_lat)
84  ny = vsize(i)
85  CASE DEFAULT
86  WRITE(*,*) "reducind dimension ", i, " ",trim(fio%get_dimname(i))
87  ierr=fio%reduce_dimension(fio%get_dimname(i), 0, 1)
88  END SELECT
89  END DO
90  WRITE(*,*) "end reduce"
91  ALLOCATE(field(nx*ny))
92  ierr=fio%read(varname,field,cunit)
93  field4d(1:nx,1:ny,1:1,1:1) => field
94  IF ( ierr /= 0 ) THEN
95  CALL error("Can't read field")
96  ELSE
97  DO i = 1, 10
98  WRITE(*,*) field4d(i,1,1,1)
99  END DO
100  ENDIF
101  DEALLOCATE(field)
102 
103  ALLOCATE(field(vsize(1)*vsize(2)))
104  field3d(1:vsize(1),1:vsize(2),1:1) => field
105  ierr=fio%read(varname,field,cunit)
106  IF ( ierr /= 0 ) THEN
107  CALL error("Can't read field")
108  ELSE
109  DO i = 1, 10
110  WRITE(*,*) field3d(i,1,1)
111  END DO
112  ENDIF
113  DEALLOCATE(field)
114  DEALLOCATE(start)
115  DEALLOCATE(vsize)
116  DEALLOCATE(atypes)
117 
118  ! interpolate to 1x1 lat lon around oslo
119  WRITE(*,*) "method ", interpol_bilinear
120  ierr = finter%interpolate(fio, interpol_bilinear, "+proj=latlon +datum=WGS84", "8,9,...,12", "58,59,...,62", .true.)
121  IF ( ierr /= 0 ) THEN
122  CALL error("Can't interpolate file")
123  END IF
124  ! Get dimensions
125  ndims=finter%get_dimensions(varname)
126  IF ( ndims <= 0 ) CALL error("Can't make slicebuilder for interpol")
127  WRITE(0,*) "inter-get_dimensions: ", ndims
128 
129  ALLOCATE(start(ndims))
130  ALLOCATE(vsize(ndims))
131  ierr = finter%get_dimension_start_size(start, vsize)
132  ALLOCATE(atypes(ndims))
133  ierr = finter%get_axistypes(atypes)
134 
135  DO i = 1, ndims
136  !WRITE (*,*) i, " axistype: ", atypes(i)
137  !WRITE (*,*) AXIS_GeoX, AXIS_GeoY, AXIS_Lon, AXIS_Lat
138  SELECT CASE (atypes(i))
139  CASE(axis_geox, axis_lon)
140  nx = vsize(i)
141  CASE(axis_geoy, axis_lat)
142  ny = vsize(i)
143  CASE DEFAULT
144  WRITE(*,*) "reducind dimension ", i, " ",trim(finter%get_dimname(i))
145  ierr=finter%reduce_dimension(finter%get_dimname(i), 0, 1)
146  END SELECT
147  END DO
148  WRITE(*,*) "end reduce"
149  ALLOCATE(field(nx*ny))
150  ierr=finter%read(varname,field,cunit)
151  field4d(1:nx,1:ny,1:1,1:1) => field
152  IF ( ierr /= 0 ) THEN
153  CALL error("Can't read field")
154  ELSE
155  WRITE(*,*) field4d
156  ENDIF
157 
158  DEALLOCATE(field)
159  ndims=finter%get_dimensions('x')
160  ALLOCATE(field(5))
161  ierr=finter%read('x',field)
162  IF ( ierr /= 0 ) THEN
163  CALL error("Can't read field x/latitude")
164  ELSE
165  WRITE(*,*) field
166  ENDIF
167 
168  DEALLOCATE(field)
169  DEALLOCATE(start)
170  DEALLOCATE(vsize)
171  DEALLOCATE(atypes)
172  ierr=finter%close()
173 
174  ! interpolate to lat lon values
175  ALLOCATE(lonvals(2))
176  ALLOCATE(latvals(2))
177  lonvals(1) = 10
178  latvals(1) = 60
179  lonvals(2) = 10
180  latvals(2) = 55
181  WRITE(*,*) "method ", interpol_bilinear
182  ierr = finter2%interpolate_lonlat(fio, interpol_bilinear, lonvals, latvals)
183  IF ( ierr /= 0 ) THEN
184  CALL error("Can't interpolate file to lonlat")
185  END IF
186  ! Get dimensions
187  ndims=finter2%get_dimensions(varname)
188  IF ( ndims <= 0 ) CALL error("Can't make slicebuilder for interpol")
189  WRITE(0,*) "inter-get_dimensions: ", ndims
190 
191  ALLOCATE(start(ndims))
192  ALLOCATE(vsize(ndims))
193  ierr = finter2%get_dimension_start_size(start, vsize)
194  ALLOCATE(atypes(ndims))
195  ierr = finter2%get_axistypes(atypes)
196 
197  DO i = 1, ndims
198  !WRITE (*,*) i, " axistype: ", atypes(i)
199  !WRITE (*,*) AXIS_GeoX, AXIS_GeoY, AXIS_Lon, AXIS_Lat
200  SELECT CASE (atypes(i))
201  CASE(axis_geox, axis_lon)
202  nx = vsize(i)
203  CASE(axis_geoy, axis_lat)
204  ny = vsize(i)
205  CASE DEFAULT
206  WRITE(*,*) "reducind dimension ", i, " ",trim(finter2%get_dimname(i))
207  ierr=finter2%reduce_dimension(finter2%get_dimname(i), 0, 1)
208  END SELECT
209  END DO
210  WRITE(*,*) "end reduce"
211  ALLOCATE(field(nx*ny))
212  ierr=finter2%read(varname,field,cunit)
213  field4d(1:nx,1:ny,1:1,1:1) => field
214  IF ( ierr /= 0 ) THEN
215  CALL error("Can't read field")
216  ELSE
217  WRITE(*,*) field4d
218  ENDIF
219 
220  ndims=finter2%get_dimensions(finter2%get_var_longitude(varname))
221  IF ( ndims == 1 ) then
222  WRITE(*,*) trim(varname), " has longitude as dimension"
223  ELSE
224  WRITE(*,*) trim(varname), " has longitude as coordinates"
225  ierr=finter2%read(finter2%get_var_longitude(varname),field)
226  field4d(1:nx,1:ny,1:1,1:1) => field
227  IF ( ierr /= 0 ) THEN
228  CALL error("Can't read field longitude")
229  ELSE
230  WRITE(*,*) field4d
231  ENDIF
232  ENDIF
233  ! write the data to testOut.nc
234  ! Open file
235  ierr=frw%open("testOut.nc","",set_filetype("netcdf",filetype_rw))
236  IF ( ierr /= 0 ) CALL error("Can't make rw-object with file: testOut.nc")
237  ndims=frw%get_dimensions("pressure")
238  IF ( ndims <= 0 ) CALL error("Can't make slicebuilder for pressure in testOut.nc")
239  WRITE(0,*) "frw-get_dimensions: ", ndims
240  !
241  ! resize the slicebuilder as above, not really needed here since output-size known
242  dimname = "time"
243  ierr=frw%reduce_dimension(dimname, 0, 1)
244  ! write the 1-d field at time 0
245  write(*,*) "writing data in cunit to t=0: ", cunit
246  ierr=frw%write("pressure", field, cunit)
247  IF ( ierr /= 0 ) THEN
248  CALL error("Can't write field at t=0")
249  ENDIF
250 
251  ! write the 1-d field at time-position 1 with 10hPa more
252  ierr=frw%reduce_dimension(dimname, 1, 1)
253  field = field + 10
254  write(*,*) "writing data in cunit to t=1: ", cunit
255  ierr=frw%write("pressure", field, cunit)
256  IF ( ierr /= 0 ) THEN
257  CALL error("Can't write field at t=1")
258  ENDIF
259 
260  ! write lat/lon axes
261  ndims=frw%get_dimensions("lon")
262  ! both lonvals are 10, so I can make a 1x2 matrix, but need to reallocate
263  ! lonvals
264  DEALLOCATE(lonvals)
265  ALLOCATE(lonvals(1))
266  lonvals(1) = 10
267  ierr=frw%write("lon", lonvals, "degrees_east")
268  IF ( ierr /= 0 ) CALL error("Can't write lon")
269  ndims=frw%get_dimensions("lat")
270  ierr=frw%write("lat", latvals, "degrees_north")
271  IF ( ierr /= 0 ) CALL error("Can't write lat")
272 
273  ierr = frw%close()
274 
275 
276  DEALLOCATE(field)
277 
278  ierr=finter2%close()
279  ! Close file (free memory)
280  ierr=fio%close()
281  IF ( ierr /= 0 ) CALL error("Can't free memory")
282 
283  END IF
284 END PROGRAM fortran_test
285 
286 SUBROUTINE error(error_string)
287  IMPLICIT NONE
288  CHARACTER(*),INTENT(IN) :: error_string
289 
290  WRITE(*,*) '******** ERROR ***********'
291  WRITE(*,*) error_string
292  WRITE(*,*) '**************************'
293 
294  CALL abort()
295 END SUBROUTINE error
Generated on Mon Mar 20 2017 13:30:16 for MI - Fimex by   doxygen 1.8.11