File: C:\NOAA\NEMS_11731\src\chem\gocart\src\GMAO_Shared\GMAO_gfio\GFIO_mean.f90

1        Program GFIO_mean
2     
3     !-------------------------------------------------------------------------
4     !         NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
5     !-------------------------------------------------------------------------
6     !BOP
7     !
8     ! !ROUTINE:  GFIO_mean --- Computing monthly means.
9     !
10     !
11     ! !USAGE: see the routine usage() below
12     !
13     ! !USES:
14     !
15     
16        Implicit NONE
17     
18     ! !DESCRIPTION: Computing fvdas monthly means.
19     !
20     ! !REVISION HISTORY:
21     !
22     !  ??Jan2005  Govindaraju  Initial design.
23     !  25Jan2006  da Silva     Removed inVars(:) as it was not used and did not
24     !                          compile on Absoft.
25     !  24Mar2006  da Silva     Implemented linear combination of files.
26     !  15Dec2009  Redder       Bug fix by changing the array, field, to a 
27     !                          scalar (and the name to field_val) which also
28     !                          removed unnecessary use of memory.(Bug affected
29     !                          only the program efficiency, not the results)
30     !-------------------------------------------------------------------------
31     !EOP
32     
33     
34        character(len=*), parameter :: myname = 'GFIO_mean'
35     
36       
37     !                              -----------------------
38     !                               Hardwired Parameters
39     !                              -----------------------
40     
41           integer, parameter :: mFiles = 256       ! Max. number of input files
42           integer, parameter :: mVars  = 256       ! Max. number of variables
43           integer, parameter :: mLevs  = 256       ! Max. number of levels    
44     
45     
46     !                              -----------------------
47     !                              User Defined Parameters
48     !                              -----------------------
49     
50     
51           integer            :: nFiles             ! Actual number of input files
52           character(len=255) :: inFiles(mFiles)    ! Input file names
53           character(len=255) :: out_total_file     ! Output running total file name
54           character(len=255) :: out_counter_file   ! Output running counter file name
55           character(len=255) :: outFile            ! Output file name
56           integer  :: inc_hhmmss                   ! increment hours specified from command line
57     
58           real              :: alpha(mFiles)   ! linear combination factors
59           logical           :: linear_comb     ! whether doing linear combination
60           logical           :: xswap           ! whether to swap lon-dimension: [0,360)->[-180,180)
61     
62     
63           integer  :: im                           ! zonal dimension
64           integer  :: jm                           ! meridional dimension
65           integer  :: km                           ! vertical dimension
66     
67           integer  :: iflag                        ! Initial flag
68           integer  :: irflag                       ! Monthly Mean flag
69     
70     
71           real, pointer     :: lon(:)              ! longitudes in deg (im)
72           real, pointer     :: lat(:)              ! latitudes in deg (jm)
73           real, pointer     :: lev(:)              ! levels in hPa (km)
74      
75           integer           :: nLevs = 0           ! total number of levels
76           real, pointer     :: Levs(:)             ! vertical levels
77           character(len=25) :: cLevs(mLevs)        ! Character reprsentation of levels
78     
79           integer           :: nVars               ! Actual number of variables
80           character(len=64) :: outVars(mVars)      ! output variable names (nVars)
81           character(len=64) :: outUnits(mVars)     ! Units of output variables (nVars)
82         
83           integer          :: outPrec              ! Output file precision:
84                                                    ! 0 = 32 bits,  1 = 64bits
85     
86     !                              -----------------------
87     !                                Variable Work Space
88     !                              -----------------------
89     
90           real, pointer ::  inField(:,:,:)         ! Input variable
91           real, pointer ::  write_out(:,:,:)       ! Input variable
92           real, pointer :: outField(:,:,:)         ! Ouput variable
93           real, pointer :: cntField(:,:,:)         ! Ouput variable
94     
95           real, pointer   :: total3d(:,:,:,:)          ! Monthly total for 3d
96           real, pointer   :: total1d(:,:,:)          ! Monthly total for 1d
97           real,pointer :: kount3d(:,:,:,:)        ! Kounter for 3d
98           real,pointer :: kount1d(:,:,:)          ! Kounter for 1d
99     
100     !                                  Local Work Space
101     !                              -----------------------
102     
103           integer count, iff, it, iv, lm, itest, ii, i, j, k
104           real dx,dx1, dx2,dy, dy1, dy2, field_val
105           double precision pi
106     
107     
108     !                              -----------------------
109     !                                  Output Meta Data
110     !                              -----------------------
111     
112           character(len=255) :: title              ! meta data title
113           character(len=255) :: source             ! data source
114           character(len=255) :: contact            ! contact org.   
115           character(len=255) :: levunits           ! Vertical levels
116           character(len=25)  :: append             ! im*jm
117           real               :: missing_val
118     
119           integer, pointer :: yyyymmdd(:)          ! Date
120           integer, pointer :: hhmmss(:)            !
121           integer          :: yyyymmdd_new          ! Date
122           integer          :: hhmmss_new            !
123           integer          :: ndate   ! Date
124           integer          :: yyyymmdd1,yyyymmdd2  ! Date
125           integer          :: yyyymmddp,hhmmssp    ! previous Date & time
126           integer          :: ntimep               ! counter for total number of times previously accumulated.
127           integer          :: yyyymmdd3,hhmmss3    ! previous Date & time
128           integer          :: yyyymmddp1           ! previous Date + 1
129           integer          :: hhmmss1              ! Time
130           integer          :: ntime 
131           integer          :: timinc,timinc_new    ! Time increment
132           integer          :: timinc_save          ! Time increment
133     
134           integer          :: in_fmode = 1         ! non-zero for READ-ONLY
135           integer          :: out_fmode = 0        ! 0 for READ-WRITE 
136           integer          :: fid                  ! input file ID
137           integer          :: out_fid              ! output file ID
138           integer          :: fidt                 ! output running total file ID
139           integer          :: fidc                 ! output running counter file ID
140           integer          :: nkount
141           integer          :: rc, rc1              ! return error code
142     
143           character(len=255) :: vtitle(mVars)      ! output title
144           character(len=255) :: vunits(mVars)      ! output title
145           character(len=257) :: vName(mVars)       ! output variable names (nVars)
146           integer            :: outKm(mVars)       ! number of levels for variables;
147           real              :: valid_range_prs(2, mVars)
148           real              :: packing_range_prs(2, mVars)
149     
150     
151     !                              -----------------------
152     !                                  eta information 
153     !                              -----------------------
154     
155           integer           :: im_e                ! input zonal dimension       
156           integer           :: jm_e                ! input meridional dimension       
157           integer           :: km_e                ! input vertical dimension    
158           integer           :: lm_e                ! input time dimension    
159           integer           :: nVars_e             ! input number of variables   
160           integer           :: buf(3)
161           real              :: undef               ! Missing value
162           real, pointer     :: lon_e(:)            ! longitudes in deg (im)
163           real, pointer     :: lat_e(:)            ! latitudes in deg (jm)
164           real, pointer     :: lev_e(:)            ! levels in hPa (km)
165           integer, pointer  :: kmVar_e(:)          ! Number of vertical levels for variables
166     
167           character(len=255) :: vtitle_in(mVars)   ! output title
168           real              :: valid_range(2, mVars)
169           real              :: packing_range(2, mVars)
170           integer           :: ngatts              ! Number of attributes for GFIO
171           integer           :: imin,jmin,xmin,imax,jmax,xmax
172           logical              initial,file_exist,rms
173     !.................................................................................
174     
175     
176        initial = .true.
177        nkount = 0
178     
179     !  Get user input
180     !  --------------
181        call  Init_ ( mFiles, nFiles, inFiles, outFile,                            &
182                      out_total_file,out_counter_file,iflag,irflag,                &
183                      im, jm, km, nLevs, Levs,                                     &
184                      mVars, nVars, outVars, outPrec, append,inc_hhmmss,           &
185                      alpha, rms, linear_comb, xswap ) 
186     
187     !  Need to check whether this will work if iflag /= 1; for now disable it
188     !  ----------------------------------------------------------------------
189        if ( linear_comb.and. iflag /= 1 ) &
190             call die ( myname, 'for now, linear combination mode only with iflag=1' )
191     
192     !  --------------------------------------------------------------------------------
193     !      Set the totals and the counter initialization flag.
194     !      iflag   0     Starting date of the month, initialize the totals and counters.
195     !              1     Compute monthly means.
196     !              2     donot open the monthly mean dataset.
197     !              3     Compute monthly means using the accumulated totals and counters.
198     
199     !      irflag  0     Continue saving the running totals.
200     !              1     Accumulate the totals using the previous accumulated totals
201     !                    and compute monthly means.
202     !              2     Donot accumulate the totals, compute the monthly means using
203     !                    the previously accumulated totals and coutners.
204     !              3     No more input data to be added.
205     !  --------------------------------------------------------------------------------
206     
207         if(iflag == 3 .and. irflag == 3) then
208     
209     !      -----------------------------------------
210     !       Compute monthly means from the existing
211     !       totals and the counters.
212     !      -----------------------------------------
213     
214            call tot2mean(out_total_file,out_counter_file,outFile,outPrec,mVars)
215            stop
216         endif
217     
218     !  -------------------------
219     !  Loop over input files ...
220     !  -------------------------
221        do iff = 1, nFiles
222     
223     !    Open GFIO file
224     !    --------------
225          call GFIO_Open ( inFiles(iff), in_fmode, fid, rc )
226          if ( rc /= 0 )  call die (myname, 'can not open input files')
227     
228     !    Determine on file
229     !    ------------------
230          call GFIO_DimInquire ( fid, im_e, jm_e, km_e, lm_e, nvars_e, ngatts, rc)
231          if ( rc /= 0 )  call die (myname, 'can not do GFIO_DimInquire')
232     
233     !     -------------------------------------------------------------------
234     !     Allocate and initialize the 3d total and counters to accumulate the
235     !     monthly means.
236     !     -------------------------------------------------------------------
237     
238          if(iff == 1) then
239           print *,' im_e,jm_e,km_e,nvars_e = ',im_e,jm_e,km_e,nvars_e
240           allocate ( total3d(im_e,jm_e,km_e,nvars_e), total1d(im_e,jm_e,nvars_e),stat=rc )
241           allocate ( kount3d(im_e,jm_e,km_e,nvars_e), kount1d(im_e,jm_e,nvars_e),stat=rc )
242           total3d = 0.0
243           total1d = 0.0
244           kount3d = 0.0
245           kount1d = 0.0
246           initial = .false.
247     !
248     !    ------------------------------------------------------------
249     !      If the irflag is 1 or 2, read the running totals and counters
250     !      from the existing running total and counter files.
251     !    ------------------------------------------------------------
252     
253           if(iflag /= 1 .and. (irflag == 1 .or. irflag == 2)) then
254            print *, ' Going in to get_rtotals *** '
255            call get_rtotals(out_total_file,out_counter_file,buf,mVars,  &
256                             yyyymmdd3,hhmmss3,total3d, kount3d,total1d, &
257                             kount1d,fidt,fidc)
258            print *, ' Out from  get_rtotals *** '
259            yyyymmddp = buf(1)
260            hhmmssp = buf(2)
261            nkount = buf(3)
262           endif
263          endif
264     
265     !    Allocate memory for meta data
266          allocate ( yyyymmdd(lm_e),hhmmss(lm_e),lon_e(im_e),lat_e(jm_e),lev_e(km_e), &
267                     kmVar_e(mVars), lev(km_e), stat = rc )
268          if ( rc /= 0 )  call die (myname, 'can not allocate yyyymmdd,hhmmss,lon_e,lat_e,lev')
269     
270     !    Get meta data
271          call GFIO_Inquire ( fid, im_e, jm_e, km_e, lm_e, nVars_e,  &
272                                    title, source, contact, undef,   &
273                                    lon_e, lat_e, lev_e, levunits,   &
274                                    yyyymmdd, hhmmss, timinc,        &
275                                    vname, vtitle, vunits, kmVar_e,  &
276                                    valid_range , packing_range, rc)
277     
278     !    print *,' timinc ',timinc
279          timinc_save = timinc
280     
281          if ( rc /= 0 )  call die (myname, 'can not do GFIO_Inquire')
282     
283     !    Find long names and number of levels for meta data
284     !    -------------------------------
285          if (nVars .le. 0) then
286             nVars = nVars_e
287          end if
288     
289          lev = lev_e
290          if (nLevs .le. 0) then 
291     !       print *, "setting level information."
292             nLevs = km_e
293             km = km_e
294             allocate(Levs(nLevs), stat=rc)
295             Levs = lev_e
296          end if
297             
298          do iv = 1, nVars 
299             if (nVars .eq. nVars_e) then
300                outVars(iv) = vname(iv)
301                outUnits(iv) = vunits(iv)
302                vtitle_in(iv) = vtitle(iv)
303                outKm(iv) = kmVar_e(iv)
304                vtitle_in(iv) = vtitle(iv)
305                valid_range_prs(1,iv) = valid_range(1,iv)
306                packing_range_prs(1,iv) = packing_range(1,iv)
307                valid_range_prs(2,iv) = valid_range(2,iv)
308                packing_range_prs(2,iv) = packing_range(2,iv)
309             else
310                do itest = 1, nVars_e
311                   if ( outVars(iv) .eq. vname(itest) ) then
312                      vtitle_in(iv) = vtitle(itest)
313                      outUnits(iv) = vunits(itest)
314                      vtitle_in(iv) = vtitle(itest)
315                      outKm(iv) = kmVar_e(itest)
316                      valid_range_prs(1,iv) = valid_range(1,itest)
317                      packing_range_prs(1,iv) = packing_range(1,itest)
318                      valid_range_prs(2,iv) = valid_range(2,itest)
319                      packing_range_prs(2,iv) = packing_range(2,itest)
320                   end if
321                end do
322             end if
323          end do
324     
325     
326     !     ------------------------------------------------------------------
327     !        Acquire the date and time only if the monthly mean flag is set.
328     !     ------------------------------------------------------------------
329     
330          if(initial) then
331           yyyymmdd1    = yyyymmdd(1)
332           hhmmss1    = 120000
333           timinc_new = 060000
334          endif
335     
336           if ( iff == 1 ) then
337             if((iflag /= 1 .and. iflag /= 3) .and. (irflag == 1 .or. irflag == 2)) then
338            
339              yyyymmddp1 = yyyymmddp + 1
340              if(yyyymmdd(1) /= yyyymmddp1 ) then
341                print *,' Previous totaled date ==> ',yyyymmddp
342                print *,' Current date ==> ',yyyymmdd(1)
343                print *,' Current time ==> ',hhmmss(1)
344                print *,' Previously stored time periods ==> ',nkount
345      
346                call die (myname, ' Check the continuation date and time ')
347               endif
348              endif
349     
350     !        Create output GFIO file
351     !        -----------------------
352              if ( nLevs .le. 0 )  then
353                 nLevs = km
354                 allocate(Levs(km), stat = rc)
355                 if ( rc /= 0 )  call die (myname, 'cannot allocate Levs')
356                 Levs = lev_e
357              end if
358     !       print *, "Levs=",Levs
359              
360     
361     
362              yyyymmdd1    = yyyymmdd(1)
363              timinc_new = 060000
364     
365           end if
366     !        yyyymmdd1    = yyyymmdd(1)
367     !        hhmmss1      = hhmmss(1)
368     
369           if(iff == nFiles) yyyymmdd2 = yyyymmdd(1)
370     
371     !     Loop over times on file
372     !     -----------------------
373     
374            print *, ' Reading ',trim( inFiles(iff)),' lm_e ',lm_e, ', alpha = ', alpha(iff), linear_comb
375     
376           do it = 1, lm_e
377     
378              nkount = nkount + 1
379     
380     !        Loop over variables
381     !        -------------------
382              do iv = 1, nVars 
383           
384     !           Read variable from GFIO file
385     !           ----------------------------
386                 if ( outKm(iv) > 0 ) then             ! 3D file
387     !           Allocated memory 
388                    allocate ( inField(im_e,jm_e,km_e),stat=rc )
389                    if ( rc /= 0 )  call die (myname, 'can not allocate inField  ')
390                    call GFIO_GetVar ( fid, outVars(iv), yyyymmdd(it), hhmmss(it), &
391                                       im_e, jm_e, 1, km_e, inField, rc )
392                    if ( rc /= 0 )  call die (myname, 'something wrong in GFIO_GetVarT for 3D file')
393                 else                                       ! 2D file
394                    allocate ( inField(im_e, jm_e, 1),stat = rc )
395                    if ( rc /= 0 )  call die (myname, 'can not allocate inField ')
396                    call GFIO_GetVar ( fid, outVars(iv), yyyymmdd(it), hhmmss(it), &
397                                       im_e, jm_e, 0, 1, inField, rc )
398                    if ( rc /= 0 )  call die (myname, 'something wrong in GFIO_GetVarT for 2D file')
399                 end if
400     
401     
402     !           -------------------------------------------------
403     !           Accumulate the totals and the counters to compute
404     !           monthly means.
405     !           -------------------------------------------------
406     !            print *,' nLevs ==> ',nLevs,outKm(iv),trim(outVars(iv))
407     
408                 if ( outKm(iv) > 0 ) then
409                    do k = 1, nLevs    
410                      do  j = 1,jm_e
411                       do  i = 1,im_e
412                        if(inField(i,j,k) < 1.e+10) then
413                         field_val = alpha(iff) * inField(i,j,k)
414                         if(rms) then
415                          field_val = field_val * field_val
416                         endif
417                         total3d(i,j,k,iv) = total3d(i,j,k,iv) + field_val
418                         kount3d(i,j,k,iv) = kount3d(i,j,k,iv) + 1.0
419                        endif
420                       end do
421                      end do
422                    end do
423                 else
424         
425     !               print *, ' Computing totals for ',trim(outVars(iv)),' in ',trim(inFiles(iff))
426                      do  j = 1,jm_e
427                       do  i = 1,im_e
428                        if(inField(i,j,1) < 1.e+10) then
429                         field_val = alpha(iff) * inField(i,j,1)
430                         if(rms) then
431                          field_val = field_val * field_val
432                         endif
433                         total1d(i,j,iv) = total1d(i,j,iv) + field_val
434                         kount1d(i,j,iv) = kount1d(i,j,iv) + 1.0
435                        endif
436                       end do
437                      end do
438                 endif
439                 deallocate( inField)
440              end do  ! variables
441           end do ! times
442     
443     !     Close input file
444     !     ----------------
445           call GFIO_Close ( fid, rc )
446     
447        end do ! input files
448     !
449          if ( inc_hhmmss .ne. 999999 ) then
450              timinc = inc_hhmmss
451          else
452              timinc_new = timinc_save
453          end if
454     
455         hhmmss1    = 120000
456         if(iflag == 1 .or. iflag == 3 .or. irflag == 2 ) then
457           if(yyyymmdd_new == 999999) then
458             yyyymmdd_new = (yyyymmdd1+yyyymmdd2)/2
459             hhmmss_new = 120000
460             print *,' yyyymmdd_new, hhmmss_new ,timinc_new ==> ',yyyymmdd_new, hhmmss_new,timinc_new
461           endif
462     
463         if((iflag == 3 .and. irflag == 2) ) then
464             yyyymmdd_new = (yyyymmdd3+yyyymmdd2)/2
465             hhmmss_new = 120000
466             print *,' yyyymmdd_new, hhmmss_new ,timinc_new ==> ',yyyymmdd_new, hhmmss_new,timinc_new
467           endif
468     
469            
470              inquire(file=trim(outFile),exist=file_exist)
471              if(file_exist) then
472                call GFIO_Open ( outFile, out_fmode, out_fid, rc )
473              else
474                if(rms) then
475                   title = 'RMS'
476                else
477                   title = 'MEAN'
478                endif
479     
480                call GFIO_Create ( outFile, title, source, contact, undef,         &
481                                   im_e, jm_e, km_e, lon_e, lat_e, Levs, levunits, &
482                                   yyyymmdd_new,hhmmss_new,timinc_new,             &
483                                   nVars, outVars, vtitle_in, outUnits, outKm,     &
484                                   valid_range_prs, packing_range_prs, outPrec,    &
485                                   out_fid, rc )
486              endif
487     
488     !        print *, ' rc ==> ',rc
489              if ( rc /= 0 )  call die (myname, 'wrong in GFIO_Create')
490     
491              print *,' Computing the monthly means '
492         elseif(iflag == 0) then
493              print *,' GFIO_Create ',trim(out_total_file)
494              if(ndate /= 999999) then
495               yyyymmdd3 = ndate
496              else
497               yyyymmdd3 = yyyymmdd1
498              endif
499              if(ntime /= 999999) then
500               hhmmss3   = ntime
501              else
502               hhmmss3   = hhmmss1
503              endif
504              print *,' yyyymmdd1,hhmmss1,timinc_new ',yyyymmdd3,hhmmss3,timinc_new
505              call GFIO_Create ( out_total_file, title, source, contact, undef,  &
506                                 im_e, jm_e, km_e, lon_e, lat_e, Levs, levunits, &
507                                 yyyymmdd3,hhmmss3,timinc_new,                   &
508                                 nVars, outVars, vtitle_in, outUnits, outKm,     &
509                                 valid_range_prs, packing_range_prs, outPrec,    &
510                                 fidt, rc )
511              if ( rc /= 0 )  call die (myname, 'wrong in GFIO_Create')
512     
513              print *,'  GFIO_Create ',trim(out_counter_file)
514              call GFIO_Create ( out_counter_file, title, source, contact, undef, &
515                                 im_e, jm_e, km_e, lon_e, lat_e, Levs, levunits,  &
516                                 yyyymmdd3,hhmmss3,timinc_new,                    &
517                                 nVars, outVars, vtitle_in, outUnits, outKm,      &
518                                 valid_range_prs, packing_range_prs, outPrec,     &
519                                 fidc, rc )
520              if ( rc /= 0 )  call die (myname, 'wrong in GFIO_Create')
521              print *,' Creation is completed '
522         endif
523     
524     
525     !   In linear combination mode, all counts reset to 1 as there is no 
526     !   normalization
527     !   -----------------------------------------------------------------
528         if ( linear_comb ) then
529              kount3d = 1.0
530              kount1d = 1.0 
531         endif
532     
533     !        Loop over variables
534     !        -------------------
535     
536             if(irflag <= 1) then
537               print *,' Adding the global attributes ', yyyymmdd(lm_e),hhmmss(lm_e),nkount
538              buf(1) = yyyymmdd(lm_e)
539              buf(2) = hhmmss(lm_e)
540              buf(3) = nkount
541              call GFIO_PutIntAtt ( fidt, 'yymmddp', 1, buf(1),outPrec, rc )
542              call GFIO_PutIntAtt ( fidt, 'hhmmssp', 1, buf(2),outPrec, rc )
543              call GFIO_PutIntAtt ( fidt, 'ntimep',  1, buf(3),outPrec, rc )
544              call GFIO_PutIntAtt ( fidc, 'yymmddp', 1, buf(1),outPrec, rc )
545              call GFIO_PutIntAtt ( fidc, 'hhmmssp', 1, buf(2),outPrec, rc )
546              call GFIO_PutIntAtt ( fidc, 'ntimep',  1, buf(3),outPrec, rc )
547             endif
548     
549             do iv = 1, nVars 
550              if(ndate /= 999999) yyyymmdd_new = ndate
551              if(ntime /= 999999) hhmmss_new = ntime
552              if(iflag == 1 .or. irflag == 2) then
553                 if ( outKm(iv) > 0 ) then
554                   allocate (outField(im_e,jm_e,km_e),stat=rc )
555                   do k = 1, km_e 
556                    do j = 1,jm_e
557                     do i = 1,im_e
558                      if(kount3d(i,j,k,iv) > 0) then
559                       outField(i,j,k) = total3d(i,j,k,iv)/kount3d(i,j,k,iv)
560                       if(rms) then
561                        outField(i,j,k) = sqrt(outField(i,j,k))
562                       endif
563                      else
564                       outField(i,j,k) = undef
565                      endif
566                     end do
567                    end do
568                   end do
569     
570     !           Write monthly mean variable to output file
571     !           -----------------------------------------           
572                    allocate ( write_out(im_e,jm_e,nLevs), stat = rc)
573                    if ( rc /= 0 )  call die (myname, 'cannot allocate write_out')
574     
575                    write_out = outField
576                    call GFIO_PutVar (out_fid,outVars(iv),yyyymmdd_new,hhmmss_new,  &
577                                      im_e, jm_e, 1, nLevs, write_out, rc )
578                    if ( rc /= 0 )  call die (myname, 'can not write')
579                    deallocate(outField,write_out)
580                 else
581                    allocate ( outField(im_e, jm_e, 1),  stat = rc )
582                    do j = 1,jm_e
583                     do i = 1,im_e
584                      if(kount1d(i,j,iv) > 0) then
585                       outField(i,j,1) = total1d(i,j,iv)/kount1d(i,j,iv)
586                       if(rms) then
587                        outField(i,j,1) = sqrt(outField(i,j,1))
588                       endif
589                      else
590                       outField(i,j,1) = undef
591                      endif
592                     end do
593                    end do
594     
595                    call GFIO_PutVar (out_fid,outVars(iv),yyyymmdd_new,hhmmss_new,  &
596                                      im_e, jm_e, 0, 1, outField, rc )
597                    if ( rc /= 0 )  call die(myname, 'can not write 2D file' )
598                    deallocate(outField)
599                 end if
600                else
601                  if(ndate /= 999999) yyyymmdd3 = ndate
602                  if(ntime /= 999999) hhmmss3   = ntime
603                  if(irflag <= 1) then
604     
605                   if ( outKm(iv) > 0 ) then
606                    allocate (outField(im_e,jm_e,km_e),stat=rc )
607                    allocate (cntField(im_e,jm_e,km_e),stat=rc )
608     
609                    do k = 1, km_e 
610                      do j = 1,jm_e
611                       do i = 1,im_e
612                          outField(i,j,k) = total3d(i,j,k,iv)
613                          cntField(i,j,k) = kount3d(i,j,k,iv)
614                       end do
615                      end do
616                    end do
617                    call GFIO_PutVar (fidt,outVars(iv),yyyymmdd3,hhmmss3,  &
618                                      im_e, jm_e, 1, nLevs, outField, rc )
619     
620                    call GFIO_PutVar (fidc,outVars(iv),yyyymmdd3,hhmmss3,  &
621                                      im_e, jm_e, 1, nLevs, cntField, rc )
622                   else
623                    allocate ( outField(im_e, jm_e, 1),  stat = rc )
624                    allocate ( cntField(im_e, jm_e, 1),  stat = rc )
625                    do j = 1,jm_e
626                     do i = 1,im_e
627                       outField(i,j,1) = total1d(i,j,iv)
628                       cntField(i,j,1) = kount1d(i,j,iv)
629                     end do
630                    end do
631      
632                    call GFIO_PutVar (fidt,outVars(iv),yyyymmdd3,hhmmss3,  &
633                                      im_e, jm_e, 0, 1, outField, rc )
634     
635                    call GFIO_PutVar (fidc,outVars(iv),yyyymmdd3,hhmmss3,  &
636                                      im_e, jm_e, 0, 1, cntField, rc )
637                    if ( rc /= 0 )  call die(myname, 'can not write 2D file' )
638                   endif
639                   deallocate(outField,cntField)
640     
641                 endif
642                endif
643              end do  ! variables
644     
645           deallocate( total3d, total1d,kount3d,kount1d)
646     !     if ( outKm(iv) .gt. 0 ) deallocate( write_out )
647     !     deallocate( outField )
648     
649     
650     !  Close output file
651     !  ----------------
652        if(iflag == 1 .or. irflag == 2) then
653         call GFIO_Close ( out_fid, rc )
654        endif
655     
656        if( iflag /= 1 .and. irflag <= 1) then
657         call GFIO_Close ( fidt, rc )
658         call GFIO_Close ( fidc, rc )
659        endif
660     
661     
662     !  All done
663     !  --------
664        call exit(0)
665     
666     CONTAINS
667     
668     !-------------------------------------------------------------------------
669     !         NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
670     !-------------------------------------------------------------------------
671     !BOP
672     !
673     ! !IROUTINE:  Init_ --- Parses command line 
674     ! 
675     ! !INTERFACE:
676     !
677     
678        subroutine Init_ ( mFiles, nFiles, inFiles, outFile,                          &
679                           out_total_file, out_counter_file,iflag,irflag,             &
680                           im, jm, km, nLevs, Levs,                                   &
681                           mVars, nVars, outVars, outPrec, append,inc_hhmmss, &
682                           alpha, rms, linear_comb, xswap ) 
683     
684     !
685     ! !USES:
686     !
687        Implicit NONE
688     
689     !
690     ! !INPUT PARAMETERS: 
691     !
692     
693           integer, intent(in)  :: mFiles           !  Max. number of input files as
694                                                    !   declared in calling program
695           integer, intent(in)  :: mVars            !  Max. number of variables
696     
697     
698     !
699     ! !OUTPUT PARAMETERS:
700     !
701     
702           integer, intent(out)          :: nFiles           !  Actual number of input files
703           character(len=*), intent(out) :: inFiles(:)       !  Input file names
704           character(len=*), intent(out) :: out_total_file   !  Input/output total files
705           character(len=*), intent(out) :: out_counter_file !  Input/output counter file
706           character(len=*), intent(out) :: outFile          !  Output file name 
707           character(len=*), intent(out) :: append       !  im*jm
708           integer  :: inc_hhmmss                        ! increment hours specified from command line
709     
710     
711           integer, intent(out)  :: im              !  zonal dimension
712           integer, intent(out)  :: jm              !  meridional dimension
713           integer, intent(out)  :: km              !  vertical dimension
714           integer, intent(out)  :: iflag           !  Initial flag
715           integer, intent(out)  :: irflag          !  Running total flag.
716     
717           real, pointer         :: lev(:)          ! levels in hPa (km)
718           real, pointer         :: Levs(:)         ! actual levels
719           integer, intent(out)  :: nLevs           ! actual number of levels
720           
721     
722           integer,          intent(out) :: nVars        ! Actual number of variables
723           character(len=*), intent(out) :: outVars(:)   ! output variable names (nVars)
724     
725           real,             intent(out) :: alpha(mFiles) ! coeficients for linear combination
726           logical,          intent(out) :: linear_comb  ! whether or not doing linear combination
727           logical,          intent(out) :: xswap  ! whether or not to swap lon-
728                                                   !  dimension; use this to change
729                                                   !  from [0,36) to [-180,180).
730     
731           logical,             intent(out) :: rms ! Flag for rms (true compute the rms
732                                                   ! false compute mean, default)
733     
734           integer, intent(out)          :: outPrec ! Output file precision:
735                                                    ! 0 = 32 bits,  1 = 6 4bits
736     
737     ! !DESCRIPTION: This routine initializes {\tt fvReduce}. It parses the command.
738     !
739     ! !REVISION HISTORY: 
740     !
741     ! Oct2002  Baoyu Yin Initial design and prologue.
742     !
743     !EOP
744     !-------------------------------------------------------------------------
745     
746        integer             iarg, argc
747        integer :: iargc
748        character(len=2048)  argv
749     
750        character(len=255)   rcfile, label, var, Vars(mVars)
751     
752        integer, parameter :: mKm = 256  ! max no of levels
753     
754        integer i, j, n, nVars0, rc, ios
755        real    xWest, p
756        logical :: debug = .false.
757        character(len=10) nLx, nLy
758     
759     print *
760     print *, "-------------------------------------------------------------------"
761     print *, "GFIO_mean - Compute monthly means.                                  "
762     print *, "-------------------------------------------------------------------"
763     print *
764     
765        argc = iargc()
766        if ( argc < 1 ) call usage_()
767     
768     !  Defaults
769     !  --------
770        nFiles = 0
771        nVars = 0
772        outFile = 'GFIO_mean.prs.nc4'
773        out_total_file = 'GFIO_total.prs.nc4'
774        out_counter_file = 'GFIO_counter.prs.nc4'
775     
776        outPrec = 0
777        km = -1
778        im =  288 
779        jm =  181 
780        nLx = '288'
781        nLy = '181'
782        iflag  = 1
783        irflag = 2
784        rms    = .false.
785        inc_hhmmss = 99999999
786        append = '288x181'
787        yyyymmdd_new = 999999
788        hhmmss_new = 999999
789        ndate      = 999999
790        ntime      = 999999
791     
792        iarg = 0
793        do i = 1, 32767
794           iarg = iarg + 1
795           if ( iarg .gt. argc ) then
796                exit
797           endif
798           call GetArg ( iArg, argv )
799           if(index(argv,'-o') .gt. 0 ) then
800              if ( iarg+1 .gt. argc ) call usage_()
801              iarg = iarg + 1
802              call GetArg ( iArg, outFile )
803           elseif(index(argv,'-tfile') .gt. 0 ) then
804              if ( iarg+1 .gt. argc ) call usage_()
805              iarg = iarg + 1
806              call GetArg ( iArg, out_total_file )
807           elseif(index(argv,'-cfile') .gt. 0 ) then
808              if ( iarg+1 .gt. argc ) call usage_()
809              iarg = iarg + 1
810              call GetArg ( iArg, out_counter_file )
811           else if(index(argv,'-inc') .gt. 0 ) then
812              if ( iarg+1 .gt. argc ) call usage_()
813              iarg = iarg + 1
814              call GetArg ( iArg, argv )
815              read(argv,*) inc_hhmmss
816           else if(index(argv,'-date') .gt. 0 ) then
817              if ( iarg+1 .gt. argc ) call usage_()
818              iarg = iarg + 1
819              call GetArg ( iArg, argv )
820              read(argv,*) ndate
821              yyyymmdd_new = ndate
822           else if(index(argv,'-time') .gt. 0 ) then
823              if ( iarg+1 .gt. argc ) call usage_()
824              iarg = iarg + 1
825              call GetArg ( iArg, argv )
826              read(argv,*) ntime
827              hhmmss_new = ntime
828           else if(index(argv,'-xswap') .gt. 0 ) then
829                xswap = .true.
830           else if(index(argv,'-rms') .gt. 0 ) then
831                rms = .true.
832           else if(index(argv,'-iflag') .gt. 0 ) then
833              if ( iarg+1 .gt. argc ) call usage_()
834              iarg = iarg + 1
835              call GetArg ( iArg, argv )
836              read(argv,*) iflag
837           else if(index(argv,'-irflag') .gt. 0 ) then
838              if ( iarg+1 .gt. argc ) call usage_()
839              iarg = iarg + 1
840              call GetArg ( iArg, argv )
841              read(argv,*) irflag
842           else if(index(argv,'-vars') .gt. 0 ) then
843              if ( iarg+1 .gt. argc ) call usage_()
844              iarg = iarg + 1
845              call GetArg ( iArg, argv )
846              call split_ ( ',', argv, mVars, Vars, nVars )
847           else if(index(argv,'-levels') .gt. 0 ) then
848              if ( iarg+1 .gt. argc ) call usage_()
849              iarg = iarg + 1
850              call GetArg ( iArg, argv )
851              call split_ ( ',', argv, mLevs, cLevs, nLevs )
852              allocate( Levs(nLevs), stat = rc)
853              if ( rc /= 0 )  call die (myname, 'wrong in allocating nLevs')
854     
855              km = nLevs
856              do k = 1, nLevs
857                 read(cLevs(k),*) Levs(k)
858              enddo
859           else if(index(argv,'-prec') .gt. 0 ) then
860              if ( iarg+1 .gt. argc ) call usage_()
861              iarg = iarg + 1
862              call GetArg ( iArg, argv )
863              read(argv,*) outPrec
864           else if(index(argv,'-d') .gt. 0 ) then
865              debug = .true.
866           else
867              nFiles = nFiles + 1
868              inFiles(nFiles) = argv
869           end if
870     
871        end do
872     
873        if ( outPrec .eq. 32 ) outPrec = 0   
874        if ( outPrec .eq. 64 ) outPrec = 1   
875     
876        if(iflag < 3 .and. irflag < 3) then
877         if ( nFiles .eq. 0 ) call die (myname, 'no files specified')
878        endif
879     
880     !  Resolve variable classes
881     !  ------------------------
882        do n = 1, nVars
883            outVars(n) = Vars(n)
884        end do 
885     
886     !  Handle linear combination specification. This is specificed as part 
887     !  of the file name: Example"
888     !                   GFIO_mean:  1:filename1.nc4 -1:filename2.nc4 will
889     !  will subtract filename2 from filename2. The format for each file
890     !  name is "alpha:filename". If alpha is omitted it is assumed 1.
891     !  This utility will work in linear combination mode if at least one 
892     !  alpha is specfied
893     !......................................................................
894        alpha = 1.0
895        linear_comb = .false.
896        xswap = .false.
897        do n = 1, nFiles
898           i = index(inFiles(n),',') 
899           if ( i > 1 ) then
900                read(inFiles(n)(1:i-1),*,iostat=ios) alpha(n)
901                if ( ios /= 0 ) then
902                     print *, 'GFIO_mean: invalid alpha in '//trim(inFiles(n))
903                     call exit(1)
904                else
905                    infiles(n) =  inFiles(n)(i+1:)
906                    linear_comb = .true.
907                end if
908           end if
909        end do
910     
911        print *, 'Input   Files: ', nFiles
912        if ( linear_comb ) then
913           do i = 1, nFiles
914              print *, "               ", trim(inFiles(i)), ', alpha = ', alpha(i)
915           end do
916        else
917           do i = 1, nFiles
918              print *, "               ", trim(inFiles(i)), alpha(i)
919           end do
920        end if
921     
922     !   print *
923     !   print *, 'Output   File: ', trim(outFile), ', prec = ', outPrec
924     !   print *
925     
926        if ( nLevs .gt. 0 ) then
927           write(*,'(a,i3,/(10x,6f10.2))') '        Levels: ', nLevs,Levs
928        end if 
929        print *
930        if ( nVars .gt. 0 ) then
931           write(*,'(a,i3,/(10x,6a10))') '     Variables: ', nVars, outVars(1:nVars)
932        end if 
933     
934        end subroutine Init_
935     
936     
937         subroutine Usage_()
938        
939     print *, "NAME"
940     print *, "   GFIO_mean  Computes mean or linear combination of files"
941     print *
942     print *, "SYNOPYSIS"
943     print *, "   GFIO_mean  [options] input_fname(s)"
944     print *
945     print *, "OPTIONS"
946     print *
947     print *, "  -o         ofname    output file name (default: GFIO_mean.prs.nc4)"
948     print *, "  -tfile     tfname    output Running total file (default: GFIO_total.prs.nc4)"
949     print *, "  -cfile     cfname    output file name (default: GFIO_counter.prs.nc4)" 
950     print *, "  -prec n     precision: "
951     print *, "                    n=0 for 32 bits (default) "
952     print *, "                    n=1 for 64 bits"
953     print *, "  -vars      varn     actual variable names, e.g., -vars hght,uwnd"
954     print *, "                      (default: all variables in the input file)"
955     print *, "  -date      ndate    Date to be intialized for the monthly mean (if not given will be computed)."
956     print *, "  -time      ntime    Time to be intialized for the monthly mean (if not given will be computed)."
957     print *, "  -iflag     flag     Initial flag. 0  Starting point for the running total"
958     print *, "                                       Running total and counter files will be created. "
959     print *, "                                    1  Compute the monthly means  (DEFAULT).            "
960     print *, "                                    2  Donot open the monthly mean dataset.             "
961     print *, "                                    3  Compute monthly means using the accumulated      "
962     print *, "                                       totals and counters.                             "
963     print *
964     print *, "  -irflag    flag                   0  Continue to accumulate the running totals and    "
965     print *, "                                       counters.                                        "
966     print *, "                                    1   Accumulate the totals using the previously      "
967     print *, "                                        accumulated totals                              "
968     print *, "                                    2  Save the running totals or compute               "
969     print *, "                                       monthly means (DEFAULT).                         "
970     print *, "  -rms       flag                   true, Computes rms                                   "
971     print *, "                                    false, Computes mean (default)                       "
972     print * 
973     print *, "  -levels vertical levels (default: all levels in the input file)"
974     print *, "  -xswap  swap longitudinal dimension; useful to change from [0,360) to"
975     print *, "          -[180,180) or vice versa"
976     print *
977     print *, "DESCRIPTION"
978     print *, "  Computes FVDAS monthly means or accumulated totals and counters."
979     print *
980     print *, "    Start accumulate the totals and the counters from the starting date and time.  "
981     print *
982     print *, "   ex:  GFIO_mean.x -irflag 0 -iflag 0 -tfile totals.nc4 -cfile counters.nc4 input_files"
983     print * 
984     print *, "   Continue to accumulate the totals and counters using next set of data. "
985     print * 
986     print *, "   ex:  GFIO_mean.x -iflag 2 -irflag 1 -tfile totals.nc4 -cfile counters.nc4 input_files"
987     print *
988     print *, "   Compute the monthly means already accumulated totals and counters without any input files."
989     print * 
990     print *, "   ex:  GFIO_mean.x -iflag 3 -irflag 3 -o out_file -tfile totals.nc4 -cfile counters.nc4 input_files"
991     print *
992     print *, "   Compute monthly means using the previously accumulated totals and counters. "
993     print * 
994     print *, "   ex:  GFIO_mean.x -iflag 3 -irflag 2  -o out_file -tfile totals.nc4 -cfile counters.nc4 input_files"
995     print *
996     print *, "   Compute monthly means. "
997     print * 
998     print *, "   ex:  GFIO_mean.x -iflag 1 -irflag 2 -o output_file  input_files"
999     print * 
1000     print *, "   Compute linear combination of files: result = 1.2*file1 - 3.2 * file2 + file3"
1001     print * 
1002     print *, "   ex:  GFIO_mean.x -o output_file  1.2,file1.nc4 -3.2,file2.nc4 file3.nc4"
1003     print * 
1004     print *, "   This feature requires iflag=1 for now. You need at least one occurence of"
1005     print *, "   ',' in the file name to trigger the linear combination mode."
1006     
1007     
1008         call die ( myname, 'exiting' )
1009     
1010         end subroutine Usage_
1011     
1012     !............................................................................
1013     
1014         subroutine split_ ( tok, str, mList, List, nList )
1015         implicit NONE
1016         character(len=1), intent(in)  :: tok  ! delimitter
1017         character(len=*), intent(in)  :: str  ! string
1018         integer,          intent(in)  :: mList
1019         character(len=*), intent(out) :: List(mList)
1020         integer,          intent(out) :: nList
1021     
1022         integer i, l, n, j
1023     
1024         i = 1
1025         l = len_trim(str)
1026         nList = 0
1027         do n = 1, mList
1028            j = index(trim(str(i:)),tok) - 1
1029            if ( j < 1 ) then
1030                 nList = nList + 1
1031                 List(nList) = str(i:)
1032                 exit
1033            end if
1034            nList = nList + 1
1035            List(nList) = str(i:i+j-1)
1036            i = i + j + 1
1037            if ( i > l ) then
1038                 exit
1039            endif
1040         end do
1041     
1042         end subroutine split_
1043     
1044     !
1045           subroutine die(myname,string)
1046           character(len=*) myname
1047           character(len=*) string
1048     !
1049           print *, ' --------------------------------'
1050           print *, '        ',myname
1051           print *, '   ',string
1052           print *, ' --------------------------------'
1053           call exit(1)
1054           return
1055           end subroutine die
1056     !
1057          subroutine get_rtotals(out_total_file,out_counter_file,buf,mVars,    &
1058                                 yyyymmdd3,hhmmss3,total3d, kount3d,total1d, &
1059                                 kount1d,fidt,fidc)
1060                                                                                                                               
1061           character(len=255) :: title              ! meta data title
1062           character(len=255) :: source             ! data source
1063           character(len=255) :: contact            ! contact org.
1064           character(len=255) :: levunits           ! Vertical levels
1065           character(len=25)  :: append             ! im*jm
1066           real               :: missing_val,undef
1067                                                                                                                                    
1068           character(len=*), intent(out) :: out_total_file   !  Input/output total files
1069           character(len=*), intent(out) :: out_counter_file !  Input/output counter file
1070           integer          :: fidt               ! Output running total file   ID
1071           integer          :: fidc               ! Output running counter file ID
1072           integer          :: in_fmode = 1       ! non-zero for READ-ONLY
1073           integer          :: out_fmode = 0      ! 0 for READ-WRITE
1074           integer          :: yyyymmdd3,hhmmss3  ! date & time
1075           integer          :: im_t,jm_t,km_t,lm_t,it
1076           integer          :: nvars_t,ngattst,rc,mVars
1077           integer          :: buf(3)
1078           integer, pointer :: yyyymmdd(:)          ! Date
1079           integer, pointer :: hhmmss(:)            ! time
1080           real, pointer     :: lon_t(:)            ! longitudes in deg (im)
1081           real, pointer     :: lat_t(:)            ! latitudes in deg (jm)
1082           real*8, pointer     :: lev_t(:)            ! levels in hPa (km)
1083           integer, pointer  :: kmVar_t(:)          ! Number of vertical levels for variables
1084     
1085           real, pointer   :: total3d(:,:,:,:)        ! Monthly total for 3d
1086           real, pointer   :: total1d(:,:,:)          ! Monthly total for 1d
1087           real, pointer   :: InField(:,:,:)          ! Monthly total for 1d
1088           real,pointer :: kount3d(:,:,:,:)        ! Kounter for 3d
1089           real,pointer :: kount1d(:,:,:)          ! Kounter for 1d
1090           character(len=255) :: vtitle(mVars)      ! output title
1091           character(len=255) :: vunits(mVars)      ! output title
1092           character(len=257) :: vName(mVars)       ! output variable names (nVars)
1093           integer            :: outKm(mVars)       ! number of levels for variables;
1094           integer            :: timinc
1095           real              :: valid_range_prs(2, mVars)
1096           real              :: packing_range_prs(2, mVars)
1097     
1098            
1099     !    Open GFIO running total and counter files
1100     !    --------------
1101          call GFIO_Open ( out_total_file, out_fmode, fidt, rc )
1102          if ( rc /= 0 )  call die (myname, 'can not open total file')
1103            
1104          call GFIO_Open ( out_counter_file, out_fmode, fidc, rc )
1105          if ( rc /= 0 )  call die (myname, 'can not open counter file')
1106     !    Determine on file
1107     !    ------------------
1108          call GFIO_DimInquire ( fidt, im_t, jm_t, km_t, lm_t, nvars_t, ngattst, rc)
1109          if ( rc /= 0 )  call die (myname, 'can not do GFIO_DimInquire')
1110     
1111            
1112     !    Allocate memory for meta data
1113     
1114          allocate ( yyyymmdd(lm_t),hhmmss(lm_t),lon_t(im_t),lat_t(jm_t),lev_t(km_t), &
1115                     kmVar_t(nVars_t), lev(km_t), stat = rc )
1116     !    Get meta data
1117            
1118          call GFIO_Inquire ( fidt, im_t, jm_t, km_t, lm_t, nVars_t,  &
1119                                    title, source, contact, undef,   &
1120                                    lon_t, lat_t, lev_t, levunits,   &
1121                                    yyyymmdd, hhmmss, timinc,        &
1122                                    vname, vtitle, vunits, kmVar_t,  &
1123                                    valid_range_prs , packing_range_prs, rc)
1124            
1125          print *,' yyyymmdd, hhmmss, timinc ',yyyymmdd, hhmmss, timinc
1126     
1127          if ( rc /= 0 )  call die (myname, 'can not allocate yyyymmdd,hhmmss,lon,lat,lev')
1128            
1129              it = 1
1130     
1131              call GFIO_GetIntAtt ( fidt, 'yymmddp', 1, buf(1), rc )
1132              call GFIO_GetIntAtt ( fidt, 'hhmmssp', 1, buf(2), rc )
1133              call GFIO_GetIntAtt ( fidt, 'ntimep',  1, buf(3), rc )
1134              call GFIO_GetIntAtt ( fidc, 'yymmddp', 1, buf(1), rc )
1135              call GFIO_GetIntAtt ( fidt, 'hhmmssp', 1, buf(2), rc )
1136              call GFIO_GetIntAtt ( fidt, 'ntimep',  1, buf(3), rc )
1137            
1138     !        Loop over variables
1139     !        -------------------
1140              do iv = 1, nVars_t
1141            
1142     !           Read variable from GFIO file
1143            
1144                 if ( kmVar_t(iv) > 0 ) then             ! 3D file
1145     !           Allocated memory
1146                    allocate ( inField(im_t,jm_t,km_t),stat=rc )
1147                    if ( rc /= 0 )  call die (myname, 'can not allocate  ')
1148            
1149                    call GFIO_GetVar ( fidt, vname(iv), yyyymmdd(it), hhmmss(it), &
1150                                       im_t, jm_t, 1, km_t, inField, rc )
1151                    if ( rc /= 0 )  call die (myname, 'something wrong in GFIO_GetVarT for 3D total file')
1152                    total3d(1:im_t,1:jm_t,1:km_t,iv) = inField(1:im_t,1:jm_t,1:km_t)
1153            
1154                    call GFIO_GetVar ( fidc, vname(iv), yyyymmdd(it), hhmmss(it), &
1155                                       im_t, jm_t, 1, km_t, inField, rc )
1156                    if ( rc /= 0 )  call die (myname, 'something wrong in GFIO_GetVarT for 3D counter file')
1157     
1158                    kount3d(1:im_t,1:jm_t,1:km_t,iv) = int(inField(1:im_t,1:jm_t,1:km_t))
1159                 else                                       ! 2D file
1160                    allocate ( inField(im_t, jm_t, 1),stat = rc )
1161                    if ( rc /= 0 )  call die (myname, 'can not allocate ')
1162                    call GFIO_GetVar ( fidt, vname(iv), yyyymmdd(it), hhmmss(it), &
1163                                       im_t, jm_t, 0, 1, inField, rc )
1164                    if ( rc /= 0 )  call die (myname, 'something wrong in GFIO_GetVarT for 2D file')
1165     
1166                    total1d(1:im_t,1:jm_t,iv) = inField(1:im_t,1:jm_t,1)
1167                    call GFIO_GetVar ( fidc, vname(iv), yyyymmdd(it), hhmmss(it), &
1168                                       im_t, jm_t, 0, 1, inField, rc )
1169                    kount1d(1:im_t,1:jm_t,iv) = int(inField(1:im_t,1:jm_t,1))
1170                 end if
1171               end do
1172               yyyymmdd3 = yyyymmdd(1)
1173               hhmmss3   = hhmmss(1)
1174            
1175               deallocate ( inField,stat=rc )
1176               deallocate ( yyyymmdd,hhmmss,lon_t,lat_t,lev_t, &
1177                            kmVar_t, lev)
1178     !         call GFIO_Close ( fidt, rc )
1179     !         call GFIO_Close ( fidc, rc )
1180             end subroutine  get_rtotals
1181     !
1182          subroutine tot2mean(out_total_file,out_counter_file,outFile,outPrec,mVars)
1183                                                                                                                               
1184           character(len=255) :: title              ! meta data title
1185           character(len=255) :: source             ! data source
1186           character(len=255) :: contact            ! contact org.
1187           character(len=255) :: levunits           ! Vertical levels
1188           character(len=25)  :: append             ! im*jm
1189           real               :: missing_val
1190                                                                                                                                    
1191           character(len=*), intent(inout) :: outFile          !  Input/output total files
1192           character(len=*), intent(inout) :: out_total_file   !  Input/output total files
1193           character(len=*), intent(inout) :: out_counter_file !  Input/output counter file
1194           integer          :: out_fid              ! monthly mean output file ID
1195           integer          :: fidt               ! Output running total file   ID
1196           integer          :: fidc               ! Output running counter file ID
1197           integer          :: in_fmode = 1       ! non-zero for READ-ONLY
1198           integer          :: out_fmode = 0      ! 0 for READ-WRITE
1199           integer          :: im_t,jm_t,km_t,lm_t,it
1200           integer          :: nvars_t,ngattst,rc,mVars
1201           integer          :: buf(3)
1202           integer, pointer :: yyyymmdd(:)          ! Date
1203           integer, pointer :: hhmmss(:)            ! time
1204           integer          :: yyyymmdd_new          ! Date
1205           integer          :: hhmmss_new            ! time
1206           real, pointer     :: lon_t(:)            ! longitudes in deg (im)
1207           real, pointer     :: lat_t(:)            ! latitudes in deg (jm)
1208           real*8, pointer     :: lev_t(:)            ! levels in hPa (km)
1209           integer, pointer  :: kmVar_t(:)          ! Number of vertical levels for variables
1210           integer           :: timinc_new,yyyymmddn
1211           real, pointer   :: totField(:,:,:)       !  totals
1212           real, pointer   :: kntField(:,:,:)       ! Counters    
1213           real, pointer   :: outField(:,:,:)       ! Monthly 
1214           character(len=255) :: vtitle(mVars)      ! output title
1215           character(len=255) :: vunits(mVars)      ! output title
1216           character(len=257) :: vName(mVars)       ! output variable names (nVars)
1217           integer            :: outKm(mVars),dd    ! number of levels for variables;
1218           integer            :: timinc
1219           real              :: valid_range_prs(2, mVars)
1220           real              :: packing_range_prs(2, mVars),undef
1221           integer          :: outPrec              ! Output file precision:
1222     
1223            
1224     !    Open GFIO running total and counter files
1225     !    --------------
1226          call GFIO_Open ( out_total_file, out_fmode, fidt, rc )
1227          if ( rc /= 0 )  call die (myname, 'can not open total file')
1228            
1229          call GFIO_Open ( out_counter_file, out_fmode, fidc, rc )
1230          if ( rc /= 0 )  call die (myname, 'can not open counter file')
1231     !    Determine on file
1232     !    ------------------
1233          call GFIO_DimInquire ( fidt, im_t, jm_t, km_t, lm_t, nvars_t, ngattst, rc)
1234          if ( rc /= 0 )  call die (myname, 'can not do GFIO_DimInquire')
1235     
1236     !    Allocate memory for meta data
1237     
1238          allocate ( yyyymmdd(lm_t),hhmmss(lm_t),lon_t(im_t),lat_t(jm_t),lev_t(km_t), &
1239                     kmVar_t(nVars_t), lev(km_t), stat = rc )
1240     !    Get meta data
1241            
1242          call GFIO_Inquire ( fidt, im_t, jm_t, km_t, lm_t, nVars_t,  &
1243                                    title, source, contact, undef,   &
1244                                    lon_t, lat_t, lev_t, levunits,   &
1245                                    yyyymmdd, hhmmss, timinc,        &
1246                                    vname, vtitle, vunits, kmVar_t,  &
1247                                    valid_range_prs , packing_range_prs, rc)
1248            
1249          if ( rc /= 0 )  call die (myname, 'can not allocate yyyymmdd,hhmmss,lon,lat,lev')
1250     
1251            
1252              it = 1
1253     
1254              call GFIO_GetIntAtt ( fidt, 'yymmddp', 1, buf(1), rc )
1255              call GFIO_GetIntAtt ( fidt, 'hhmmssp', 1, buf(2), rc )
1256              call GFIO_GetIntAtt ( fidt, 'ntimep',  1, buf(3), rc )
1257              call GFIO_GetIntAtt ( fidc, 'yymmddp', 1, buf(1), rc )
1258              call GFIO_GetIntAtt ( fidt, 'hhmmssp', 1, buf(2), rc )
1259              call GFIO_GetIntAtt ( fidt, 'ntimep',  1, buf(3), rc )
1260     
1261              yyyymmdd_new = (buf(1) + yyyymmdd(1))/2
1262              hhmmss_new = 120000
1263              timinc_new = timinc
1264              print *,' yyyymmdd,hhmmss ',yyyymmdd,hhmmss
1265              print *,'  yyyymmdd_new hhmmss_new timinc_new ', yyyymmdd_new,hhmmss_new,timinc_new
1266              
1267     ! 
1268     !            ----------------------------------------------
1269     !              Create the monthly mean output file.
1270     !            ----------------------------------------------
1271     
1272              call GFIO_Create ( outFile, title, source, contact, undef,           &
1273                                 im_t, jm_t, km_t, lon_t, lat_t, lev_t, levunits,  &
1274                                 yyyymmdd_new,hhmmss_new,timinc_new,               &
1275                                 nVars_t, vname, vtitle, vunits, kmVar_t,          &
1276                                 valid_range_prs, packing_range_prs, outPrec,      &
1277                                 out_fid, rc )
1278     !        print *, ' rc ==> ',rc
1279              if ( rc /= 0 )  call die (myname, 'wrong in GFIO_Create')
1280     
1281     !        Loop over variables
1282     !        -------------------
1283              do iv = 1, nVars_t
1284            
1285     !           Read variable from GFIO file
1286            
1287                 if ( kmVar_t(iv) > 0 ) then             ! 3D file
1288     !           Allocated memory
1289                    allocate ( totField(im_t,jm_t,km_t),stat=rc )
1290                    if ( rc /= 0 )  call die (myname, 'can not allocate  ')
1291                    allocate ( kntField(im_t,jm_t,km_t),stat=rc )
1292                    if ( rc /= 0 )  call die (myname, 'can not allocate  ')
1293                    allocate ( outField(im_t,jm_t,km_t),stat=rc )
1294                    if ( rc /= 0 )  call die (myname, 'can not allocate  ')
1295            
1296                    call GFIO_GetVar ( fidt, vname(iv), yyyymmdd(it), hhmmss(it), &
1297                                       im_t, jm_t, 1, km_t, totField, rc )
1298                    if ( rc /= 0 )  call die (myname, 'something wrong in GFIO_GetVarT for 3D total file')
1299            
1300                    call GFIO_GetVar ( fidc, vname(iv), yyyymmdd(it), hhmmss(it), &
1301                                       im_t, jm_t, 1, km_t, kntField, rc )
1302                    if ( rc /= 0 )  call die (myname, 'something wrong in GFIO_GetVarT for 3D counter file')
1303      
1304                    print *,' Field ',vname(iv)
1305                    do  k = 1,km_t
1306                     do j = 1,jm_t
1307                      do  i = 1,im_t
1308                        if(kntField(i,j,k) > 0) then
1309                          outField(i,j,k) = totField(i,j,k)/kntField(i,j,k)
1310                        else
1311                          outField(i,j,k) = undef
1312                        endif
1313                      end do
1314                     end do
1315                    end do
1316     
1317                    call GFIO_PutVar (out_fid,vname(iv),yyyymmdd_new,hhmmss_new,  &
1318                                      im_t, jm_t, 1, km_t, outField, rc )
1319     
1320                 else                                       ! 2D file
1321                    allocate ( totField(im_t,jm_t,1),stat=rc )
1322                    allocate ( kntField(im_t,jm_t,1),stat=rc )
1323                    if ( rc /= 0 )  call die (myname, 'can not allocate ')
1324                    call GFIO_GetVar ( fidt, vname(iv), yyyymmdd(it), hhmmss(it), &
1325                                       im_t, jm_t, 0, 1, totField, rc )
1326                    if ( rc /= 0 )  call die (myname, 'something wrong in GFIO_GetVarT for 2D file')
1327     
1328                    call GFIO_GetVar ( fidc, vname(iv), yyyymmdd(it), hhmmss(it), &
1329                                       im_t, jm_t, 0, 1, kntField, rc )
1330     
1331                     print *,' Field ',vname(iv)
1332                     do j = 1,jm_t
1333                      do  i = 1,im_t
1334                        if(kntField(i,j,k) > 0) then
1335                          outField(i,j,1) = totField(i,j,1)/kntField(i,j,1)
1336                        else
1337                          outField(i,j,1) = undef
1338                        endif
1339                      end do
1340                     end do
1341     
1342                    call GFIO_PutVar (out_fid,vname(iv),yyyymmdd_new,hhmmss_new,  &
1343                                      im_t, jm_t, 0, 1, outField, rc )
1344                 end if
1345                 deallocate(outField,totField,kntField)
1346               end do
1347            
1348               deallocate ( yyyymmdd,hhmmss,lon_t,lat_t,lev_t, &
1349                            kmVar_t, lev)
1350               call GFIO_Close ( out_fid, rc )
1351               call GFIO_Close ( fidt, rc )
1352               call GFIO_Close ( fidc, rc )
1353             end subroutine  tot2mean
1354     end Program GFIO_mean
1355