/[lmdze]/trunk/libf/IOIPSL/flincom.f90
ViewVC logotype

Contents of /trunk/libf/IOIPSL/flincom.f90

Parent Directory Parent Directory | Revision Log Revision Log


Revision 42 - (show annotations)
Thu Mar 24 11:52:41 2011 UTC (13 years, 2 months ago) by guez
File size: 14047 byte(s)
Removed programs "test_inter_barxy" and "test_disvert".

Added option "read" for "s_sampling" in "disvert".

1 MODULE flincom
2
3 ! From flincom.f90, version 2.2 2006/03/07 09:21:51
4
5 IMPLICIT NONE
6
7 PRIVATE
8 PUBLIC flinclo, flinopen_nozoom, flininfo, ncids
9
10 ! This is the data we keep on each file we open:
11 INTEGER, PARAMETER:: nbfile_max = 200
12 INTEGER, SAVE:: ncids(nbfile_max)
13 INTEGER, SAVE:: ncnbva(nbfile_max)
14 LOGICAL, SAVE:: ncfileopen(nbfile_max)=.FALSE.
15
16 CONTAINS
17
18 SUBROUTINE flinopen_nozoom(iim, jjm, llm, lon, lat, lev, ttm, itaus, date0, &
19 dt, fid_out)
20
21 ! This procedure opens an input file. There is no test of the
22 ! content of the file against the input from the model.
23
24 ! The user should check that the dimensions of lon lat and lev are
25 ! correct when passed to flinopen. This can be done after the call
26 ! when iim and jjm have been retrieved from the netCDF file. In
27 ! Fortran 90 this problem will be solved with an internal assign
28 ! IF iim, jjm, llm or ttm are parameters in the calling program it
29 ! will create a segmentation fault
30
31 USE calendar, ONLY: ymds2ju, ioconf_calendar
32 USE errioipsl, ONLY: histerr
33 USE netcdf, ONLY: nf90_get_att, nf90_get_var, nf90_global, &
34 nf90_inquire_variable
35
36 INTEGER, intent(in):: iim ! size in the x direction in the file (longitude)
37 INTEGER, intent(in):: jjm ! size in the y direction
38
39 INTEGER, intent(in):: llm ! number of levels
40 ! (llm = 0 means no axis to be expected)
41
42 INTEGER, intent(in):: ttm ! size of time axis
43 real, intent(out):: lon(iim, jjm) ! longitude
44 real, intent(out):: lat(iim, jjm) ! latitude
45 real, intent(out):: lev(llm)
46 INTEGER, intent(out):: itaus(ttm) ! time steps within this file
47 REAL, intent(out):: date0 ! Julian date at which itau = 0
48 REAL, intent(out):: dt ! length of the time steps of the data
49
50 INTEGER, intent(in):: fid_out
51 ! (file ID which is later used to read the data)
52
53 ! Variables local to the procedure:
54 INTEGER iret, vid, fid, nbdim, i
55 INTEGER gdtt_id, old_id, iv, gdtmaf_id
56 CHARACTER(LEN=250) name
57 CHARACTER(LEN=80) units, my_calendar
58 INTEGER year, month, day
59 REAL r_year, r_month, r_day
60 INTEGER year0, month0, day0, hours0, minutes0, seci
61 REAL sec, sec0
62 CHARACTER strc
63 REAL, DIMENSION(:), ALLOCATABLE:: vec_tmp
64
65 !---------------------------------------------------------------------
66
67 IF ((fid_out < 1) .OR. (fid_out > nbfile_max)) THEN
68 ! Either the fid_out has not been initialized (0 or very large)
69 ! then we have to open anyway. Else we only need to open the file
70 ! if it has not been opened before.
71 print *, "Call flinfo before flinopen"
72 stop 1
73 end IF
74 IF (.NOT. ncfileopen(fid_out)) THEN
75 print *, "Call flinfo before flinopen"
76 stop 1
77 end IF
78
79 ! The user has already opened the file
80 ! and we trust that he knows the dimensions
81
82 fid = ncids(fid_out)
83
84 ! 4.0 extracting the coordinates
85
86 CALL flinfindcood (fid_out, 'lon', vid, nbdim)
87 IF (nbdim == 2) THEN
88 iret = NF90_GET_VAR (fid, vid, lon, &
89 start=(/ 1, 1 /), count=(/ iim, jjm /))
90 ELSE
91 ALLOCATE(vec_tmp(iim))
92 iret = NF90_GET_VAR (fid, vid, vec_tmp, &
93 start=(/ 1 /), count=(/ iim /))
94 DO i=1, jjm
95 lon(:, i) = vec_tmp(:)
96 ENDDO
97 DEALLOCATE(vec_tmp)
98 ENDIF
99
100 CALL flinfindcood (fid_out, 'lat', vid, nbdim)
101 IF (nbdim == 2) THEN
102 iret = NF90_GET_VAR (fid, vid, lat, start=(/ 1, 1 /), &
103 count=(/ iim, jjm /))
104 ELSE
105 ALLOCATE(vec_tmp(jjm))
106 iret = NF90_GET_VAR (fid, vid, vec_tmp, start=(/ 1 /), count=(/ jjm /))
107 DO i=1, iim
108 lat(i, :) = vec_tmp(:)
109 ENDDO
110 DEALLOCATE(vec_tmp)
111 ENDIF
112
113 IF (llm > 0) THEN
114 CALL flinfindcood (fid_out, 'lev', vid, nbdim)
115 IF (nbdim == 1) THEN
116 iret = NF90_GET_VAR (fid, vid, lev, start=(/ 1 /), count=(/ llm /))
117 ELSE
118 CALL histerr (3, 'flinopen', &
119 'Can not handle vertical coordinates that have more', &
120 'than 1 dimension', ' ')
121 ENDIF
122 ENDIF
123
124 ! 5.0 Get all the details for the time if possible needed
125
126 IF (ttm > 0) THEN
127
128 ! 5.1 Find the time axis. Prefered method is the 'timestep since'
129
130 gdtmaf_id = -1
131 gdtt_id = -1
132 old_id = -1
133 DO iv=1, ncnbva(fid_out)
134 name=''
135 iret = NF90_INQUIRE_VARIABLE (fid, iv, name=name)
136 units=''
137 iret = NF90_GET_ATT (fid, iv, 'units', units)
138 IF (INDEX(units, 'seconds since') > 0) gdtmaf_id = iv
139 IF (INDEX(units, 'timesteps since') > 0) gdtt_id = iv
140 IF (INDEX(name, 'tstep') > 0) old_id = iv
141 ENDDO
142
143 IF (gdtt_id > 0) THEN
144 vid = gdtt_id
145 ELSE IF (gdtmaf_id > 0) THEN
146 vid = gdtmaf_id
147 ELSE IF (old_id > 0) THEN
148 vid = old_id
149 ELSE
150 CALL histerr (3, 'flinopen', 'No time axis found', ' ', ' ')
151 ENDIF
152
153 ALLOCATE(vec_tmp(ttm))
154 iret = NF90_GET_VAR (fid, vid, vec_tmp, &
155 start=(/ 1 /), count=(/ ttm /))
156 itaus(1:ttm) = NINT(vec_tmp(1:ttm))
157 DEALLOCATE(vec_tmp)
158
159 ! Getting all the details for the time axis
160
161 ! Find the calendar
162 my_calendar='XXXX'
163 iret = NF90_GET_ATT (fid, gdtmaf_id, 'calendar', my_calendar)
164 IF (INDEX(my_calendar, 'XXXX') < 1) THEN
165 CALL ioconf_calendar(my_calendar)
166 ENDIF
167
168 units = ''
169 iret = NF90_GET_ATT (fid, vid, 'units', units)
170 IF (gdtt_id > 0) THEN
171 units = units(INDEX(units, 'since')+6:LEN_TRIM(units))
172 READ (units, '(I4.4, 5(a, I2.2))') &
173 year0, strc, month0, strc, day0, &
174 strc, hours0, strc, minutes0, strc, seci
175 sec0 = hours0*3600. + minutes0*60. + seci
176 CALL ymds2ju (year0, month0, day0, sec0, date0)
177 iret = NF90_GET_ATT (fid, gdtt_id, 'tstep_sec', dt)
178 ELSE IF (gdtmaf_id > 0) THEN
179 units = units(INDEX(units, 'since')+6:LEN_TRIM(units))
180 READ (units, '(I4.4, 5(a, I2.2))') &
181 year0, strc, month0, strc, day0, &
182 strc, hours0, strc, minutes0, strc, seci
183 sec0 = hours0*3600. + minutes0*60. + seci
184 CALL ymds2ju (year0, month0, day0, sec0, date0)
185
186 ELSE IF (old_id > 0) THEN
187 iret = NF90_GET_ATT (fid, NF90_GLOBAL, 'delta_tstep_sec', dt)
188 iret = NF90_GET_ATT (fid, NF90_GLOBAL, 'day0', r_day)
189 iret = NF90_GET_ATT (fid, NF90_GLOBAL, 'sec0', sec)
190 iret = NF90_GET_ATT (fid, NF90_GLOBAL, 'year0', r_year)
191 iret = NF90_GET_ATT (fid, NF90_GLOBAL, 'month0', r_month)
192
193 day = INT(r_day)
194 month = INT(r_month)
195 year = INT(r_year)
196
197 CALL ymds2ju (year, month, day, sec, date0)
198 ENDIF
199 ENDIF
200
201 END SUBROUTINE flinopen_nozoom
202
203 !***************************************************************
204
205 SUBROUTINE flininfo(filename, iim, jjm, llm, ttm, fid_out)
206
207 ! This subroutine allows to get some information.
208 ! It is usually called by "flinopen" but the user may want to call
209 ! it before in order to allocate the space needed to extract the
210 ! data from the file.
211
212 USE strlowercase_m, ONLY: strlowercase
213 USE errioipsl, ONLY: histerr
214 USE netcdf, ONLY: nf90_inquire, nf90_inquire_dimension, nf90_noerr, &
215 nf90_nowrite
216 use netcdf95, only: nf95_open
217
218 CHARACTER(LEN=*), intent(in):: filename
219 INTEGER, intent(out):: iim, jjm, llm, ttm, fid_out
220
221 ! Variables local to the procedure:
222 INTEGER, SAVE:: nbfiles = 0
223 INTEGER, SAVE:: ncdims(nbfile_max, 4)
224 INTEGER iret, fid, ndims, nvars, nb_atts, id_unlim
225 INTEGER iv, lll
226 CHARACTER(LEN=80) name
227 CHARACTER(LEN=30) axname
228
229 !---------------------------------------------------------------------
230
231 lll = LEN_TRIM(filename)
232 IF (filename(lll-2:lll) /= '.nc') THEN
233 name = filename(1:lll)//'.nc'
234 ELSE
235 name = filename(1:lll)
236 ENDIF
237
238 call NF95_OPEN(name, NF90_NOWRITE, fid)
239 iret = NF90_INQUIRE(fid, nDimensions=ndims, nVariables=nvars, &
240 nAttributes=nb_atts, unlimitedDimId=id_unlim)
241
242 iim = 0
243 jjm = 0
244 llm = 0
245 ttm = 0
246
247 DO iv=1, ndims
248 iret = NF90_INQUIRE_DIMENSION(fid, iv, name=axname, len=lll)
249 CALL strlowercase(axname)
250 axname = ADJUSTL(axname)
251
252 IF ((INDEX(axname, 'x') == 1) .OR. (INDEX(axname, 'lon') == 1)) THEN
253 iim = lll
254 ELSE IF ((INDEX(axname, 'y') == 1) &
255 .OR. (INDEX(axname, 'lat') == 1)) THEN
256 jjm = lll
257 ELSE IF ((INDEX(axname, 'lev') == 1) .OR. (INDEX(axname, 'plev') == 1) &
258 .OR. (INDEX(axname, 'z') == 1) &
259 .OR. (INDEX(axname, 'depth') == 1)) THEN
260 llm = lll
261 ELSE IF ((INDEX(axname, 'tstep') == 1) &
262 .OR. (INDEX(axname, 'time_counter') == 1)) THEN
263 ! For the time we certainly need to allow for other names
264 ttm = lll
265 ELSE IF (ndims == 1) THEN
266 ! Nothing was found and ndims=1 then we have a vector of data
267 iim = lll
268 ENDIF
269 ENDDO
270
271 ! Keep all this information
272
273 nbfiles = nbfiles+1
274
275 IF (nbfiles > nbfile_max) THEN
276 CALL histerr(3, 'flininfo', &
277 'Too many files. Please increase nbfil_max', &
278 'in program flincom.F90.', ' ')
279 ENDIF
280
281 ncids(nbfiles) = fid
282 ncdims(nbfiles, :) = (/ iim, jjm, llm, ttm /)
283 ncnbva(nbfiles) = nvars
284 ncfileopen(nbfiles) = .TRUE.
285 fid_out = nbfiles
286
287 END SUBROUTINE flininfo
288
289 !***************************************************************
290
291 SUBROUTINE flinfindcood (fid_in, axtype, vid, ndim)
292
293 ! This subroutine explores the file in order to find
294 ! the coordinate according to a number of rules
295
296 USE strlowercase_m, ONLY: strlowercase
297 USE errioipsl, ONLY: histerr
298 USE netcdf, ONLY: nf90_get_att, nf90_inquire_dimension, &
299 nf90_inquire_variable, nf90_noerr
300
301 ! ARGUMENTS
302
303 INTEGER, intent(in):: fid_in
304 integer vid, ndim
305 CHARACTER(LEN=3):: axtype
306
307 ! LOCAL
308
309 INTEGER:: iv, iret, dimnb
310 CHARACTER(LEN=40):: dimname, dimuni1, dimuni2, dimuni3
311 CHARACTER(LEN=30):: str1
312 LOGICAL:: found_rule = .FALSE.
313 !---------------------------------------------------------------------
314 vid = -1
315
316 ! Make sure all strings are invalid
317
318 dimname = '?-?'
319 dimuni1 = '?-?'
320 dimuni2 = '?-?'
321 dimuni3 = '?-?'
322
323 ! First rule: we look for the correct units
324 ! lon: east
325 ! lat: north
326 ! We make an exact check as it would be too easy to mistake
327 ! some units by just comparing the substrings.
328
329 SELECTCASE(axtype)
330 CASE ('lon')
331 dimuni1 = 'degree_e'
332 dimuni2 = 'degrees_e'
333 found_rule = .TRUE.
334 CASE('lat')
335 dimuni1 = 'degree_n'
336 dimuni2 = 'degrees_n'
337 found_rule = .TRUE.
338 CASE('lev')
339 dimuni1 = 'm'
340 dimuni2 = 'km'
341 dimuni3 = 'hpa'
342 found_rule = .TRUE.
343 CASE DEFAULT
344 found_rule = .FALSE.
345 END SELECT
346
347 IF (found_rule) THEN
348 iv = 0
349 DO WHILE ((vid < 0).AND.(iv < ncnbva(fid_in)))
350 iv = iv+1
351 str1 = ''
352 iret = NF90_GET_ATT (ncids(fid_in), iv, 'units', str1)
353 IF (iret == NF90_NOERR) THEN
354 CALL strlowercase (str1)
355 IF ((INDEX(str1, TRIM(dimuni1)) == 1) &
356 .OR.(INDEX(str1, TRIM(dimuni2)) == 1) &
357 .OR.(INDEX(str1, TRIM(dimuni3)) == 1)) THEN
358 vid = iv
359 iret = NF90_INQUIRE_VARIABLE (ncids(fid_in), iv, ndims=ndim)
360 ENDIF
361 ENDIF
362 ENDDO
363 ENDIF
364
365 ! Second rule: we find specific names:
366 ! lon: nav_lon
367 ! lat: nav_lat
368 ! Here we can check if we find the substring as the
369 ! names are more specific.
370
371 SELECTCASE(axtype)
372 CASE ('lon')
373 dimname = 'nav_lon lon longitude'
374 found_rule = .TRUE.
375 CASE('lat')
376 dimname = 'nav_lat lat latitude'
377 found_rule = .TRUE.
378 CASE('lev')
379 dimname = 'plev level depth deptht'
380 found_rule = .TRUE.
381 CASE DEFAULT
382 found_rule = .FALSE.
383 END SELECT
384
385 IF (found_rule) THEN
386 iv = 0
387 DO WHILE ((vid < 0).AND.(iv < ncnbva(fid_in)))
388 iv = iv+1
389 str1=''
390 iret = NF90_INQUIRE_VARIABLE (ncids(fid_in), iv, &
391 name=str1, ndims=ndim)
392 IF (INDEX(dimname, TRIM(str1)) >= 1) THEN
393 vid = iv
394 ENDIF
395 ENDDO
396 ENDIF
397
398 ! Third rule: we find a variable with the same name as the dimension
399 ! lon = 1
400 ! lat = 2
401 ! lev = 3
402
403 IF (vid < 0) THEN
404 SELECTCASE(axtype)
405 CASE ('lon')
406 dimnb = 1
407 found_rule = .TRUE.
408 CASE('lat')
409 dimnb = 2
410 found_rule = .TRUE.
411 CASE('lev')
412 dimnb = 3
413 found_rule = .TRUE.
414 CASE DEFAULT
415 found_rule = .FALSE.
416 END SELECT
417
418 IF (found_rule) THEN
419 iret = NF90_INQUIRE_DIMENSION (ncids(fid_in), dimnb, name=dimname)
420 iv = 0
421 DO WHILE ((vid < 0).AND.(iv < ncnbva(fid_in)))
422 iv = iv+1
423 str1=''
424 iret = NF90_INQUIRE_VARIABLE (ncids(fid_in), iv, &
425 name=str1, ndims=ndim)
426 IF (INDEX(dimname, TRIM(str1)) == 1) THEN
427 vid = iv
428 ENDIF
429 ENDDO
430 ENDIF
431 ENDIF
432
433 ! Stop the program if no coordinate was found
434
435 IF (vid < 0) THEN
436 CALL histerr (3, 'flinfindcood', &
437 'No coordinate axis was found in the file', &
438 'The data in this file can not be used', axtype)
439 ENDIF
440
441 END SUBROUTINE flinfindcood
442
443 !***************************************************************
444
445 SUBROUTINE flinclo (fid_in)
446
447 USE netcdf, ONLY: nf90_close
448
449 INTEGER:: fid_in
450
451 INTEGER:: iret
452
453 !---------------------------------------------------------------------
454
455 iret = NF90_CLOSE (ncids(fid_in))
456 ncfileopen(fid_in) = .FALSE.
457
458 END SUBROUTINE flinclo
459
460 END MODULE flincom

  ViewVC Help
Powered by ViewVC 1.1.21