File: C:\NOAA\NEMS_11731\src\chem\gocart\src\GMAO_Shared\GMAO_gfio\gfioutil.f

1     !-------------------------------------------------------------------------
2     !         NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
3     !-------------------------------------------------------------------------
4     !BOP
5     !
6     ! !ROUTINE:  IdentifyDim - Identify a cooridate variable
7     !
8     ! !DESCRIPTION: This function attempts to identify a coordiante variable
9     !               from the name or units of the variable.  It does so by 
10     !               attempting to match the units specified in the COARDS 
11     !               conventions or by checking the name against commonly used
12     !               names.
13     ! !INTERFACE:
14     !
15           integer function IdentifyDim (dimName, dimUnits)
16     !
17     ! !USES:
18     !
19           Implicit NONE
20     !
21     ! !INPUT PARAMETERS:
22     !
23           character*(*) dimName  ! Name of the coordinate variable
24           character*(*) dimUnits ! Units of the coordinate variable
25     !
26     ! !RETURN VALUES:
27     !
28     !     0 = X dimension (longitude)
29     !     1 = Y dimension (latitude)
30     !     2 = Z dimension (level)
31     !     3 = Time
32     !    -1 = Unable to determine dimension
33     !
34     ! !REVISION HISTORY:
35     !
36     !  1998.12.22  Lucchesi       Initial coding.
37     !  1999.11.02  da Silva       Made LATS4D compatible,
38     !  2001.01.02  da Silva       Cimcurventing PGI bugs.
39     !
40     !EOP
41     !-------------------------------------------------------------------------
42     
43             if (TRIM(dimUnits) .EQ. "degrees_north" ) then
44                 IdentifyDim = 1
45                 return
46              end if
47     
48             if (TRIM(dimUnits) .EQ. "hPa" ) then
49                 IdentifyDim = 2
50                 return
51             end if
52     
53             if ( trim(dimName) .eq. "time" ) then 
54                 IdentifyDim = 3
55                 return
56             end if
57     
58     
59             if (TRIM(dimUnits) .EQ. "degrees_east" .OR.
60          .       trim(dimName)  .eq. "longitude"    .OR.
61          .       trim(dimName)  .eq. "lon"  ) then
62                 IdentifyDim = 0
63             else if (TRIM(dimUnits) .EQ. "degrees_north" ) then
64                 IdentifyDim = 1
65             else if (  trim(dimName)  .eq. "latitude"    .OR.
66          .             trim(dimName)  .eq. "lat"  ) then
67                 IdentifyDim = 1
68             else if (INDEX(dimName,"lev") .NE. 0 .OR.
69          .           INDEX(dimName,"Height") .NE. 0) then
70               IdentifyDim = 2
71             else if (TRIM(dimUnits) .EQ. "mb" .OR.
72          .           TRIM(dimUnits) .EQ. "millibar" .OR.
73          .           TRIM(dimUnits) .EQ. "sigma_level" .OR.
74          .           TRIM(dimUnits) .EQ. "hPa") then
75               IdentifyDim = 2
76             else if (INDEX(dimName,"TIME") .EQ. 0 .OR.
77          .           INDEX(dimName,"TIME:EOSGRID") .EQ. 0 .OR.
78          .           INDEX(dimName,"time") .EQ. 0 .OR.
79          .           INDEX(dimName,"Time") .EQ. 0) then
80               IdentifyDim = 3
81             else
82               IdentifyDim = -1
83             endif
84     
85             return
86             end
87     
88     
89     !-------------------------------------------------------------------------
90     !         NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
91     !-------------------------------------------------------------------------
92     !BOP
93     !
94     ! !ROUTINE:  GetBegDateTime - Get begin date/time of file
95     !
96     ! !DESCRIPTION: This routine returns the begin date/begin time on file.
97     !               For a native GFIO file, it simply returns the value of
98     !  global attributes begin_date/begin_time. If these do no exist then
99     !  it attempts to parse the COARDS compliant time unit.
100     !
101     ! !INTERFACE:
102     !
103           subroutine GetBegDateTime ( fid, begDate, begTime, incSecs, rc )
104     !
105     ! !USES:
106     !
107           Implicit NONE
108           include "netcdf.inc"
109           include "gfio.h"
110     !
111     ! !INPUT PARAMETERS:
112     !
113           integer fid      ! file ID
114     !
115     ! !OUTPUT PARAMETERS:
116     !
117           integer begDate  ! beginning date
118           integer begTime  ! beginning time
119           integer incSecs  ! time increment in secs
120           integer rc       ! error return code
121     !
122     ! !REVISION HISTORY:
123     !
124     !  1999.11.01  da Silva  Initial code.
125     !  1999.11.08  da Silva  Generic time coordinate variable (no name assumed).
126     !  2000.10.18  Lucchesi  Added ParseTimeUnits subroutine to handle a wider
127     !                        variety of Units formats.
128     !  2003.11.02  da Silva  On some lats4d file with 1 time, it was returning
129     !                        incSecs = 0; now it returns incSecs = 1, so no
130     !                        division by zero occurs later on.
131     !
132     !EOP
133     !-------------------------------------------------------------------------
134     
135           integer i, timeId, hour, min, sec, corner(1), err, timInc
136           integer year, month, day
137           character(len=MAXCHR) timeUnits, strTmp, dimUnits
138     
139           character*(MAXCHR) varName, dimName
140           integer type, nvDims, vdims(MAXVDIMS), nvAtts, dimSize
141           integer nDims, nvars, ngatts, dimId, IdentifyDim
142     
143     !     Time conversion local variables
144           real*4    rtime
145           real*8    dtime
146           integer*2 itime
147           integer*4 ltime
148           integer   t1, t2
149     
150     
151     !     Start by determing the ID of the time coordinate variable
152     !     ---------------------------------------------------------
153           timeId = -1
154           call ncinq (fid, nDims, nvars, ngatts, dimId, rc)
155           if (err("GetBegDateTime: ncinq failed",rc,-48) .NE. 0)return
156           do i=1,nDims
157             call ncdinq (fid, i, dimName, dimSize, rc)  
158             if (err("GetBegDateTime: can't get dim info",rc,-41) .NE. 0) return
159             dimId = ncvid (fid, dimName, rc)
160             if (err("GetBegDateTime: ncvid failed",rc,-40) .NE. 0) return
161             call ncagtc (fid, dimId, 'units', dimUnits, MAXCHR, rc)
162             if (err("GetBegDateTime: could not get units for dimension",rc,-53)
163          .       .NE. 0) return
164             if ( IdentifyDim (dimName, dimUnits) .eq. 3 ) then
165                  timeId = dimId
166                  timeUnits = dimUnits
167                  exit
168             end if
169           end do
170     
171           if ( timeId .lt. 0 ) then
172              rc = -43
173              print *, "GetBegDateTime: could not find time coord"
174              return
175           end if
176     
177           call ncpopt(0)  ! tell netcdf to shut up
178     
179     !     Try assuming this file has been written with GFIO
180     !     -------------------------------------------------
181           call ncagt (fid, timeId, 'begin_date', begDate, rc)
182           if ( rc .eq. 0 ) then
183                call ncagt (fid, timeId, 'begin_time', begTime, rc)
184           end if
185     
186     !     Well, it must be a native GFIO file
187     !     -----------------------------------
188           if ( rc .eq. 0 ) then
189              call ncagt (fid, timeId, 'time_increment', timInc, rc)
190              if (err("GetBegDateTime: missing time increment",rc,-44) .NE. 0) return
191      !ams        write (strTmp,'(i6)') timinc
192      !ams        read (strTmp,'(3I2)') hour, min, sec
193     
194              call GFIO_parseIntTime ( timinc, hour, min, sec )
195     
196              incSecs = hour*3600 + min*60 + sec
197     !ams     print *, 'begdate, begtime, incsecs: ',begdate, begtime, incsecs
198              return                               ! all done.
199           end if
200     
201     !     If could not find begin_date/begin_time attributes
202     !     then this is not a native GFIO file. In this case
203     !     attempt to parse the COARDS compliant time units
204     !     --------------------------------------------------
205     !ams      call ncagtc (fid, timeId, 'units', timeUnits, MAXCHR, rc)
206     !ams      if (err("GetBegDateTime: missing time.units",rc,-44) .NE. 0) return
207           i = index(timeUnits,'since')
208           if ( i .le. 0 ) then
209               if (err("GetBegDateTime: invalid time units",1,-44) .NE. 0) return
210           endif
211     
212     !     Call to ParseTimeUnits replaces an internal read, that made assumptions
213     !     about the format of the Time Units string that were not always true.  
214     !     (RL: 10/2000)
215     
216           call ParseTimeUnits ( timeUnits, year, month, day, hour, min, sec, rc )
217           begDate = year*10000 + month*100 + day
218           begTime = hour*10000 + min*100   + sec
219     
220     !     Determine time increment.
221     !     -------------------------
222           call ncvinq (fid, timeID, varName, type, nvDims, vDims, 
223          &     nvAtts, rc)
224           if (err("GetBegDateTime: error in time variable inquire",
225          &    rc,-52) .NE. 0) return
226          
227           if ( type .eq. NCFLOAT )  then
228                corner(1) = 1
229                call ncvgt1(fid,timeID,corner,rtime,rc)
230                t1 = int(rtime) 
231                corner(1) = 2
232                call ncvgt1(fid,timeID,corner,rtime,rc)
233                t2 = int(rtime)
234           else if ( type .eq. NCDOUBLE ) then
235                corner(1) = 1
236                call ncvgt1(fid,timeID,corner,dtime,rc)
237                t1 = int(dtime)
238     !ams       print *, t1, dtime, rc
239                corner(1) = 2
240                call ncvgt1(fid,timeID,corner,dtime,rc)
241                t2 = int(dtime)
242     !ams       print *, t2, dtime, rc
243           else if ( type .eq. NCSHORT  ) then
244                corner(1) = 1
245                call ncvgt1(fid,timeID,corner,itime,rc)
246                t1 = itime
247                corner(1) = 2
248                call ncvgt1(fid,timeID,corner,itime,rc)
249                t2 = itime
250           else if ( type .eq. NCLONG   ) then
251                corner(1) = 1
252                call ncvgt1(fid,timeID,corner,ltime,rc)
253                t1 = ltime
254                corner(1) = 2
255                call ncvgt1(fid,timeID,corner,ltime,rc)
256                t2 = ltime
257           else
258                if (err("GetBegDateTime: invalid time data type",
259          &         1,-44) .NE. 0) return
260           endif
261     
262     
263     !     Convert time increment to seconds if necessary
264     !     ----------------------------------------------
265           incSecs = t2 - t1
266           if ( timeUnits(1:6) .eq.  'minute' ) then
267                incSecs = incSecs * 60 
268           else if ( timeUnits(1:4) .eq. 'hour'   ) then
269                incSecs = incSecs * 60 * 60 
270           else if ( timeUnits(1:3) .eq.  'day' ) then
271                incSecs = incSecs * 60 * 60 * 24
272           else
273                if (err("GetBegDateTime: invalid time unit name",
274          &         1,-44) .NE. 0) return
275           endif
276     
277           incSecs = max ( 1, incSecs )
278     
279     !ams      print *, 'begdate, begtime, incsecs: ',begdate, begtime, incsecs
280     
281           rc = 0 ! all done
282           call ncpopt(NCVERBOS)
283     
284           return
285           end
286     
287     
288     !-------------------------------------------------------------------------
289     !         NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
290     !-------------------------------------------------------------------------
291     !BOP
292     !
293     ! !ROUTINE:  GFIO_GetMissing --- Return missing value
294     !
295     ! !DESCRIPTION: Returns missing value on file
296     !
297     ! !INTERFACE:
298     
299           real function GFIO_GetMissing ( fid, rc )
300     
301     !
302     ! !USES:
303     !
304           Implicit NONE
305           include "netcdf.inc"
306           include "gfio.h"
307     !
308     ! !INPUT PARAMETERS:
309     
310           integer fid      ! file id
311           integer rc       ! Error code
312     
313     !
314     ! !REVISION HISTORY:
315     !
316     !  1999.11.01  da Silva  Initial code.
317     !
318     !EOP
319     !-------------------------------------------------------------------------
320     
321           integer nDims, recdim, ngatts, seconds
322           integer varType, nvDims, vDims(MAXVDIMS), nvAtts
323           character*8 strBuf
324           character*(MAXCHR) dimName
325           character*(MAXCHR) dimUnits
326           character*(MAXCHR) vnameTemp
327           integer dimSize
328           integer i
329           logical surfaceOnly
330           logical noTimeInfo
331           integer attType, attLen
332           integer allVars            ! all variables - includes dimension vars
333     
334           real*4 amiss_32
335           integer err
336     
337     ! Get basic information from the file
338     
339           call ncpopt(0)
340     
341           call ncinq (fid,nDims,allVars,ngatts,recdim,rc)
342           if (err("Inqure: ncinq failed",rc,-48) .NE. 0) return
343     
344           if (nDims .EQ. 3) then
345             surfaceOnly = .TRUE.
346           endif
347     
348           do i= 1, allVars
349             call ncvinq (fid,i,vnameTemp,varType,nvDims,vDims,nvAtts,rc)
350             if (err("GFIO_GetMissing: variable inquire error",rc,-52) .NE. 0) return
351             if (nvDims .EQ. 1) then   ! coord variable
352               cycle
353             else                      ! noon-coord variable
354                call ncagt (fid, i,'fmissing_value',amiss_32,rc)
355                if (rc .NE. 0) then
356                    call ncainq (fid, i, 'missing_value', attType, attLen, rc)
357                   if (rc.eq.0 .and. attType .EQ. NCFLOAT) then
358                      call ncagt (fid, allVars, 'missing_value', amiss_32, rc)
359                      if (err("GFIO_GetMissing: error getting missing value",rc,-53) 
360          .                .NE. 0) return
361                   else
362                         print *, 
363          .              'GFIO_GetMissing: Cannot find missing value, assuming 1E+15'
364                         amiss_32 = 1.0E+15
365                   end if
366                endif
367                exit    ! just check first non-ccordinate variable
368             endif
369           end do
370     
371           GFIO_GetMissing = amiss_32
372     
373           call ncpopt(NCVERBOS)
374     
375           rc = 0
376           end
377     
378     !-------------------------------------------------------------------------
379     !         NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
380     !-------------------------------------------------------------------------
381     !BOP
382     !
383     ! !ROUTINE:  ParseTimeUnits -- Parse the COARDS time units string 
384     !
385     ! !DESCRIPTION: This subroutine takes as input the COARDS metadata time 
386     !               units string and parses out the date included in the string.
387     !               This date typically represents the first time in a COARDS
388     !               HDF files.
389     ! !INTERFACE:
390     !
391           subroutine ParseTimeUnits ( TimeUnits, year, month, day, hour, min, sec, rc )
392     !
393     ! !USES:
394     !
395           implicit none
396           include "netcdf.inc"
397           include "gfio.h"
398     !
399     ! !INPUT PARAMETERS:
400     !
401           character*(MAXCHR) TimeUnits      ! Units metadata string from the Time coord var
402     !
403     ! !OUTPUT PARAMETERS:
404     !
405           integer        year               ! 4-digit year
406           integer        month              ! month
407           integer        day                ! day
408           integer        hour               ! hour
409           integer        min                ! minute
410           integer        sec                ! second
411           integer        rc                 ! return code
412                                             !  0 = no error
413                                             ! -1 = problem parsing string
414     
415     ! !REVISION HISTORY:
416     !
417     !  2000.10.18  Lucchesi       Initial coding.
418     !  2001.08.13  Lucchesi       Modified to better parse time:units string (needed for lats4d support)
419     !
420     !EOP
421     !-------------------------------------------------------------------------
422     
423     
424           ! Local variables
425     
426           character*(MAXCHR)  NewUnits
427           integer ypos(2), mpos(2), dpos(2), hpos(2), minpos(2), spos(2)
428           integer i, j, inew, strlen
429           integer firstdash, lastdash
430           integer firstcolon, lastcolon
431           integer lastspace
432     
433           strlen = LEN_TRIM (TimeUnits)
434           
435           firstdash = index(TimeUnits, '-')
436           lastdash  = index(TimeUnits, '-', BACK=.TRUE.)
437     
438           if (firstdash .LE. 0 .OR. lastdash .LE. 0) then
439             rc = -1
440             return
441           endif
442           
443           ypos(2) = firstdash - 1
444           mpos(1) = firstdash + 1
445           ypos(1) = ypos(2) - 4
446     
447           mpos(2) = lastdash - 1
448           dpos(1) = lastdash + 1
449           dpos(2) = dpos(1) + 2
450     
451           read ( TimeUnits(ypos(1):ypos(2)), * ) year
452           read ( TimeUnits(mpos(1):mpos(2)), * ) month
453           read ( TimeUnits(dpos(1):dpos(2)), * ) day
454     
455           firstcolon = index(TimeUnits, ':')
456     
457           if (firstcolon .LE. 0) then
458             
459             ! If no colons, check for hour.
460     
461             lastspace = index(TRIM(TimeUnits), ' ', BACK=.TRUE.)
462             if ((strlen-lastspace).eq.2 .or. (strlen-lastspace).eq.3) then
463               hpos(1) = lastspace+1
464               hpos(2) = strlen-1
465               read (TimeUnits(hpos(1):hpos(2)), * ) hour
466               min  = 0
467               sec  = 0
468             else
469               print *, 'ParseTimeUnits: Assuming a starting time of 00z'
470               hour = 0
471               min  = 0
472               sec  = 0
473             endif
474     
475           else
476             hpos(1) = firstcolon - 2
477             hpos(2) = firstcolon - 1
478             lastcolon =  index(TimeUnits, ':', BACK=.TRUE.)
479             if ( lastcolon .EQ. firstcolon ) then
480               mpos(1) = firstcolon + 1
481               mpos(2) = firstcolon + 2
482               read (TimeUnits(hpos(1):hpos(2)), * ) hour
483               read (TimeUnits(mpos(1):mpos(2)), * ) min
484               sec = 0
485             else 
486               mpos(1) = firstcolon + 1
487               mpos(2) = lastcolon - 1
488               spos(1) = lastcolon + 1
489               spos(2) = lastcolon + 2
490               read (TimeUnits(hpos(1):hpos(2)), * ) hour
491               read (TimeUnits(mpos(1):mpos(2)), * ) min
492               read (TimeUnits(spos(1):spos(2)), * ) sec
493             endif
494           endif
495            
496           rc = 0
497           return
498           end
499     
500           subroutine GFIO_parseIntTime ( hhmmss, hour, min, sec )      
501           integer, intent(in)  :: hhmmss
502           integer, intent(out) :: hour, min, sec 
503           hour = hhmmss / 10000
504           min  = mod(hhmmss,10000)/100
505           sec  = mod(hhmmss,100)
506           end
507