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 |