File: C:\NOAA\NEMS_11731\src\chem\gocart\src\GMAO_Shared\Chem_Shared\Chem_UtilMod.F90

1     #include "MAPL_Generic.h"
2     #if 0
3     #define min(x,y) amin(real(x),real(y))
4     #define MIN(x,y) AMIN(real(x),real(y))
5     #define max(x,y) amax(real(x),real(y))
6     #define MAX(x,y) AMAX(real(x),real(y))
7     #endif
8     
9     !-------------------------------------------------------------------------
10     !         NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
11     !-------------------------------------------------------------------------
12     !BOP
13     !
14     ! !MODULE:  Chem_UtilMod --- Assorted Utilities for fvChem
15     !
16     ! !INTERFACE:
17     !
18     
19        module  Chem_UtilMod
20     
21     ! !USES:
22     
23        use ESMF_Mod
24        use MAPL_Mod
25     
26        use MAPL_CommsMod
27        use MAPL_CFIOMod
28     !  use HorzBinMod
29        use MAPL_HorzTransformMod
30     
31        use Chem_Mod                  ! Chemistry Base Class
32        use mod_diag                  ! fvGCM diagnostics
33        use m_die
34        use m_StrTemplate             ! string templates
35        use m_chars, only: uppercase
36     
37        implicit NONE
38     
39     !
40     ! !PUBLIIC MEMBER FUNCTIONS:
41     !
42     
43        PRIVATE
44        PUBLIC  Chem_UtilMPread          ! Array reader
45        PUBLIC  Chem_UtilNegFiller       ! Fills negative values in a column
46        PUBLIC  Chem_UtilTroppFixer      ! Repairs tropopause pressure bad values
47        PUBLIC  Chem_UtilGetTimeInfo     ! Time info on file
48        PUBLIC  Chem_UtilExtractIntegers ! Extract integers from a character-delimited string
49        PUBLIC  Chem_BiomassDiurnal      ! Biomass burning diurnal cycle
50     
51        PUBLIC tick      ! GEOS-4 stub
52        PUBLIC mcalday   ! GEOS-4 stub
53        PUBLIC pmaxmin   ! functional
54        PUBLIC zenith    ! GEOS-4 stub
55     
56     !
57     ! !DESCRIPTION:
58     !
59     !  This module implements assorted odds & ends for fvChem.
60     !
61     ! !REVISION HISTORY:
62     !
63     !  29oct2003  da Silva  First crack.
64     !  16aug2005  da Silva  Introduced scatter from MAPL_CommsMod.
65     !
66     !EOP
67     !-------------------------------------------------------------------------
68     
69        interface pmaxmin
70          module procedure pmaxmin2d
71          module procedure pmaxmin3d
72        end interface
73     
74     CONTAINS
75     
76     #ifdef USE_MAPL_MPREAD
77     
78     !-------------------------------------------------------------------------
79     !     NASA/GSFC, Global Modeling and Assimilation Office, Code 900.3     !
80     !-------------------------------------------------------------------------
81     !BOP
82     !
83     ! !IROUTINE:  Chem_UtilMPread --- Reads fields from file and distribute
84     !
85     ! !INTERFACE:
86     !
87        subroutine Chem_UtilMPread_g5 ( filen, varn, nymd, nhms, &
88                                     i1, i2, ig, im, j1, j2, jg, jm, km, &
89                                     var2d, var3d, cyclic, grid )
90     
91     ! !USES:
92     
93       implicit NONE
94     
95     ! !INPUT PARAMETERS:
96     
97       character(len=*), intent(in) :: filen   ! GFIO compatible file name
98       character(len=*), intent(in) :: varn    ! variable name
99       integer, intent(in)          :: nymd, nhms ! date/time
100     
101                                               ! Distributed grid info:
102       integer,      intent(in)     :: i1, i2  !   local  zonal indices
103       integer,      intent(in)     :: ig      !   zonal ghosting
104       integer,      intent(in)     :: im      !   global zonal dimension
105       integer,      intent(in)     :: j1, j2  !   local meridional indices
106       integer,      intent(in)     :: jg      !   meridional ghosting
107       integer,      intent(in)     :: jm      !   global zonal dimension
108       integer,      intent(in)     :: km      !   vertical dimension
109     
110       logical, OPTIONAL, intent(in) :: cyclic ! whether time dimension is periodic
111      
112                                               ! ESMF Grid; this is required
113                                               !  in GEOS-5 under ESMF
114       type(ESMF_Grid), OPTIONAL, intent(in) :: grid 
115     
116                                                       
117     
118     
119     ! !OUTPUT PARAMETERS:
120     
121       real, OPTIONAL, intent(out), target :: var2d(i1-ig:i2+ig,j1-jg:j2+jg)
122       real, OPTIONAL, intent(out), target :: var3d(i1-ig:i2+ig,j1-jg:j2+jg,km)
123     
124     ! !DESCRIPTION: 
125     !
126     ! !REVISION HISTORY:
127     !
128     !  15Aug2006  da Silva  It is now a simple wrap around CFIOReadArray().
129     !
130     !EOP
131     !-------------------------------------------------------------------------
132     
133         character(len=*), parameter ::  Iam = 'Chem_UtilMPread_g5'
134         type(ESMF_TIME)             :: time
135         integer                     :: yy, mm, dd, h, m, s, rc, STATUS
136     
137         real, pointer               :: ptr2(:,:), ptr3(:,:,:)
138     
139         
140     !   Convert time to ESMF format
141     !   ---------------------------
142         call parseIntTime_ ( nymd, yy, mm, dd )
143         call parseIntTime_ ( nhms,  h,  m,  s )
144         call ESMF_TimeSet(time, yy=yy, mm=mm, dd=dd,  h=h,  m=m, s=s, rc=status )
145         if ( status /= 0 ) call die(Iam,'failed to convert time')
146     
147     !   Read either 2D or 3D array
148     !   --------------------------
149         if ( .not. present(grid) ) &
150            call die ( Iam,'when running under the ESMF "grid" must be specified' )
151     
152         if ( present(var2d) ) then
153     
154              ptr2 => var2d
155              call MAPL_CFIORead ( varn, filen, time, grid, ptr2, rc=STATUS,  &
156                                   verbose = .true., time_is_cyclic=cyclic,   &
157                                   time_interp = .true. )
158     
159         else if ( present(var3d) ) then
160     
161              ptr3 => var3d
162              call MAPL_CFIORead ( varn, filen, time, grid, ptr3, rc=STATUS,  &
163                                   verbose = .true., time_is_cyclic=cyclic,   &
164                                   time_interp = .true. )
165     
166         else
167     
168            call die ( Iam,'either "var2d" or "var3d" must be specified' )
169     
170         end if
171     
172         if ( status /= 0 ) call die(Iam,'cannot read '//trim(varn))
173     
174     
175     CONTAINS
176         subroutine parseIntTime_ ( hhmmss, hour, min, sec )      
177           integer, intent(in)  :: hhmmss
178           integer, intent(out) :: hour, min, sec 
179           hour = hhmmss / 10000
180           min  = mod(hhmmss,10000)/100
181           sec  = mod(hhmmss,100)
182         end subroutine parseIntTime_
183     
184       end subroutine Chem_UtilMPread_g5
185     
186     #endif
187     
188     !-------------------------------------------------------------------------
189     !     NASA/GSFC, Global Modeling and Assimilation Office, Code 900.3     !
190     !-------------------------------------------------------------------------
191     !BOP
192     !
193     ! !IROUTINE:  Chem_UtilMPread_g4 --- Reads fields from file and distribute
194     !
195     ! !INTERFACE:
196     !
197        subroutine Chem_UtilMPread ( filen, varn, nymd, nhms,             &
198                                     i1, i2, ig, im, j1, j2, jg, jm, km,  &
199                                     var2d, var3d, cyclic, grid,          &
200                                     ForceBinning,                        &
201                                     maskString, gridMask, instanceNumber )
202     
203     ! !USES:
204     
205       implicit NONE
206     
207     ! !INPUT PARAMETERS:
208     
209       character(len=*), intent(in) :: filen   ! GFIO compatible file name
210       character(len=*), intent(in) :: varn    ! variable name
211       integer, intent(in)          :: nymd, nhms ! date/time
212     
213                                               ! Distributed grid info:
214       integer,      intent(in)     :: i1, i2  !   local  zonal indices
215       integer,      intent(in)     :: ig      !   zonal ghosting
216       integer,      intent(in)     :: im      !   global zonal dimension
217       integer,      intent(in)     :: j1, j2  !   local meridional indices
218       integer,      intent(in)     :: jg      !   meridional ghosting
219       integer,      intent(in)     :: jm      !   global zonal dimension
220       integer,      intent(in)     :: km      !   vertical dimension
221     
222       logical, OPTIONAL, intent(in) :: cyclic ! whether time dimension is periodic
223                                               ! ESMF Grid; this is required
224                                               !  in GEOS-5 under ESMF
225     
226                                            
227            
228       logical, OPTIONAL, intent(in) :: ForceBinning ! Whether remapping uses
229                                                     !  binning or interpolation
230                                                     ! Starting with ARCTAS default
231                                                     ! is .TRUE. so binning always
232                                                     ! takes place.
233     
234       type(ESMF_Grid), OPTIONAL, intent(in) :: grid
235     
236       character(len=*), OPTIONAL, intent(in) :: maskString !Delimited string of integers
237       real, OPTIONAL, intent(in) :: gridMask(i1:i2,j1:j2)  !Grid mask (NOTE: No ghosting) 
238       integer, OPTIONAL, intent(in) :: instanceNumber      !Instantiation, for debugging.
239     
240     ! !OUTPUT PARAMETERS:
241     
242       real, OPTIONAL, intent(out)   :: var2d(i1-ig:i2+ig,j1-jg:j2+jg)
243       real, OPTIONAL, intent(out)   :: var3d(i1-ig:i2+ig,j1-jg:j2+jg,km)
244     
245     ! !DESCRIPTION: 
246     !
247     ! !REVISION HISTORY:
248     !
249     !  28oct2003 da Silva  First crack.
250     !  03ayg2004 da Silva  Uses GetVarT for time interpolation
251     !  18nov2004 da Silva  Added cyclic option for climatological files.
252     !  30may2005 da Silva  Introduced template expansing in file names.
253     !  14aug2006 da Silva  Swap longitudes of input arrays (always), and
254     !                      made it read only one layer at a time. Notice that
255     !                      the lon-swap is dangerous and it is a temporary
256     !                      device until re retire this routine in favor of
257     !                      CFIOArrayRead().
258     !  15Aug2006 da Silva  Renamed with _g4 as this is now obsolete.
259     !  21Aug2006 da Silva  Back to GFIO compatible version; renamed CFIO one g5;
260     !                      now the GFIO one does horizontal interp/binning
261     !                      and swap only if necessary.
262     !  29Feb2008 Nielsen   Masking
263     !
264     !EOP
265     !-------------------------------------------------------------------------
266     
267         character(len=*), parameter ::  myname = 'Chem_UtilMPread_g4'
268         character(len=*), parameter ::  Iam=myname 
269         logical :: tcyclic, verb
270     
271         integer  :: READ_ONLY=1, nokm=0
272         integer :: fid, rc, ios, k
273         character(len=257) :: fname, vname
274         real :: const
275         integer :: imf, jmf, kmf                  ! dimensions on file
276         real,  pointer :: latf(:), lonf(:), levf(:)
277         real,  pointer :: lat(:),  lon(:)
278     
279         real, pointer :: ptr2(:,:), ptr2f(:,:)
280     !
281     !** added for the NEMS/GFS-GOCART                           ! Sarah Lu
282         real :: local(i1-ig:i2+ig,j1-jg:j2+jg)                  ! Sarah Lu
283     
284         logical :: doingMasking, ForceBinning_
285     
286         INTEGER, ALLOCATABLE :: regionNumbers(:),flag(:)
287         INTEGER, ALLOCATABLE :: mask(:,:)
288     
289     !   type(HorzBinTransform)     :: Trans
290         type(MAPL_HorzTransform)   :: Trans
291     
292         vname = trim(varn) ! make a copy
293     
294     !   Consistency check
295     !   -----------------
296         if ( .not. ( present(var2d) .or. present(var3d) ) ) then
297            call die ( myname, 'missing var2d or var3d' )
298         else if ( present(var2d) .and. present(var3d) ) then
299            call die ( myname, 'either var2d or var3d, but not both' )
300         end if
301     
302         if ( present(ForceBinning) ) then
303              ForceBinning_ = ForceBinning
304         else
305              ForceBinning_ = .true.
306         endif
307     
308         if ( present(cyclic) ) then
309             tcyclic = cyclic
310         else
311             tcyclic = .false. ! by default time dimension is not periodic
312         end if
313     
314         if ( .not. present(grid) ) then
315            call die ( myname, 'ESMF is a must when running under the ESMF' )
316         end if
317     
318     !   When masking, both the mask and the string
319     !   of integers (region numbers) must be present
320     !   --------------------------------------------
321         IF ( (PRESENT(maskString) .AND. .NOT. PRESENT(gridMask)) .OR.   &
322              (PRESENT(gridMask) .AND. .NOT. PRESENT(maskString))      ) &
323            CALL die ( myname, ": Both gridMask and maskString must be specified." )
324     
325         IF(PRESENT(gridMask)) THEN
326            IF(TRIM(maskString) == "-1") THEN
327             doingMasking = .FALSE.
328            ELSE
329             doingMasking = .TRUE.
330            END IF
331         ELSE
332          doingMasking = .FALSE.
333         END IF
334     
335     !   Masking initialization
336     !   ----------------------
337         SetMask: IF(doingMasking) THEN
338     
339          k = 32
340          ALLOCATE(regionNumbers(k),flag(k),mask(i1:i2,j1:j2),STAT=ios)
341          IF ( ios /= 0 ) CALL die ( myname, ": Cannot allocate for masking.")
342     
343     !   Obtain region numbers from delimited list of integers
344     !   -----------------------------------------------------
345          regionNumbers(:) = 0
346          CALL Chem_UtilExtractIntegers(maskString,k,regionNumbers,RC=ios)
347          IF ( ios /= 0 ) CALL die ( myname, ": Unable to extract integers for regionNumbers.")
348     
349     !   How many integers were found?
350     !   -----------------------------
351          flag(:) = 1
352          WHERE(regionNumbers(:) == 0) flag(:) = 0
353          k = SUM(flag)
354          DEALLOCATE(flag,STAT=ios)
355          IF ( ios /= 0 ) CALL die ( myname, ": Cannot dallocate flag.")
356     
357     !   Set local mask to 1 where gridMask matches each integer (within precision!) 
358     !   ---------------------------------------------------------------------------
359          mask(i1:i2,j1:j2) = 0
360          DO ios=1,k
361           WHERE(regionNumbers(ios)-0.01 <= gridMask .AND. &
362                 gridMask <= regionNumbers(ios)+0.01) mask = 1
363          END DO
364     
365         END IF SetMask
366     
367     !   Expand templates
368     !   ----------------
369         if ( index(filen,'%') .gt. 0 ) then
370              call StrTemplate ( fname, filen, xid='unknown', &
371                                 nymd=nymd, nhms=nhms )
372         else
373              fname = filen
374         end if
375     
376     !   Trick: when filename is "/dev/null" it is assumed the user wants to
377     !   set the veriable to a constant
378     !   -------------------------------------------------------------------
379         if ( fname(1:9) == '/dev/null' ) then    
380              ios = -1
381              k = index(fname,':')
382              if ( k > 9 ) then
383                   read(fname(k+1:),*,iostat=ios) const
384              end if
385              if ( ios /= 0 ) const = 0.0
386     
387              IF ( PRESENT(var2d) ) THEN
388               var2d = const
389               IF(doingMasking .AND. const /= 0.0) THEN
390                WHERE(mask == 0) var2d = 0.0
391               END IF
392              END IF
393     
394              IF ( PRESENT(var3d) ) THEN
395               var3d = const
396               IF(doingMasking .AND. const /= 0.0) THEN
397                DO k=1,km
398                 WHERE(mask == 0) var3d(:,:,k) = 0.0
399                END DO
400               END IF
401              END IF
402     
403              if ( MAPL_am_I_root() ) then
404                 print *, myname // ': input file is /dev/null'
405                 print *, myname // ':    setting variable ' // trim(vname) // &
406                                    ' to constant ', const
407              end if
408              return    ! exit here
409         end if
410     
411     !   Get the lat/lons from the input grid
412     !   ------------------------------------
413         call GridGetLatLons_ ( grid, lon, lat )
414     
415         local(:,:) = 0.0                                                    ! Sarah Lu
416     
417     !
418     !   For GFS, all PEs owns its own global array and do it's own scatter  ! Sarah Lu
419     
420     !   Read file
421     !   ---------
422     !*  if ( MAPL_am_I_root() ) then                                        ! Sarah Lu
423     
424     !      Allocate work space for scatter
425     !      -------------------------------
426            allocate(ptr2(im,jm),stat=ios)
427            if ( ios /= 0 ) call die ( myname, 'cannot allocate ptr2' )
428     
429     
430            print *, myname // ': Reading ' // trim(vname) // ' from GFIO file ' &
431                     // trim(fname) // ' at ', nymd, nhms
432     
433     !      Open file
434     !      ---------
435            call GFIO_Open ( fname, READ_ONLY, fid, rc )
436            if ( rc .ne. 0 ) then
437               call die(myname,'cannot open GFIO file '//trim(fname))
438            end if
439     
440     !      Query file metadata
441     !      -------------------
442            call queryFile_ ( fid, imf, jmf, kmf, lonf, latf, levf, vname )
443     
444     !      If binning, define transform
445     !      ----------------------------
446            if (im < imf .OR. jm < jmf .OR. ForceBinning_ ) then
447               call MAPL_HorzTransformCreate (Trans, imf, jmf, im, jm, rc=rc)
448               if ( rc /= 0 ) call die(myname,'cannot create transform',rc)
449            end if
450     
451     !      Allocate work space to reada data in
452     !      ------------------------------------
453            allocate(ptr2f(imf,jmf),stat=ios)
454            if ( ios /= 0 ) call die ( myname, 'cannot allocate ptr2f' )
455     
456     !*  end if ! masterproc                                                  ! Sarah Lu
457     
458     !   2D Variables
459     !   ------------
460         verb = .true.
461         if ( present(var2d) ) then
462     
463     !        Read global array and swap longitudes
464     !        -------------------------------------
465     !*       if ( MAPL_am_I_root() ) then                                    ! Sarah Lu
466     
467     !           Read the file
468     !           -------------
469                 call GFIO_GetVarT1 ( fid, trim(vname), nymd, nhms, imf, jmf, &
470                                      nokm, 1, ptr2f, rc, tcyclic, fid )
471                 if ( rc .ne. 0 ) call die(myname,'cannot read '//trim(vname) )
472     
473     !           Interpolate/bin if necessary
474     !           ----------------------------
475                 call Regrid_ ( ptr2f, imf, jmf, ptr2, im, jm, lonf, lon, verb )
476     
477     !*       end if ! masterproc                                            ! Sarah Lu
478     
479     !        Scatter the array
480     !        -----------------           
481     !*       call ArrayScatter ( var2d, ptr2, grid, rc=ios )                ! Sarah Lu
482     !*        if ( ios /= 0 ) call die ( myname, 'cannot scatter '//trim(vname) ) !Sarah Lu
483     
484     !* 
485     !* NOTE: emissions are S_to_N while GFS is N_to_S                        ! Sarah Lu
486              call GFS_Simple_Scatter ( ptr2f(:,jmf:1:-1), local )            ! Sarah lu
487              var2d(:,:) = local(:,:)                                         ! Sarah Lu
488     
489     !        Apply mask when present
490     !        -----------------------
491              IF(doingMasking) THEN
492               WHERE(mask(i1:i2,j1:j2) == 0) var2d(i1:i2,j1:j2) = 0.00
493              END IF
494     
495     !   3D Variables
496     !   ------------
497         else
498     
499           do k = 1, km
500     
501     !        Read 1 level of global array and swap longitudes
502     !        ------------------------------------------------
503     !*       if ( MAPL_am_I_root() ) then                                    ! Sarah Lu
504         
505     !            Read the file
506     !            -------------
507                  call GFIO_GetVarT1 ( fid, trim(vname), nymd, nhms, imf, jmf, k, &
508                                       1, ptr2f, rc, tcyclic, fid )
509                  if ( rc .ne. 0 ) call die(myname,'cannot read '//trim(vname) )
510     
511     !            Interpolate/bin if necessary
512     !            ----------------------------
513                  call Regrid_ ( ptr2f, imf, jmf, ptr2, im, jm, lonf, lon, verb ) 
514     
515     !*       end if ! masterproc                                             ! Sarah Lu
516     
517     !        Scatter the array with 1 level
518     !        ------------------------------           
519     !*       call ArrayScatter ( var3d(:,:,k), ptr2, grid, rc=ios )          ! Sarah Lu
520     !*       if ( ios /= 0 ) call die ( myname, 'cannot scatter v3d'//trim(vname)) !Sarah Lu
521     
522              call GFS_Simple_Scatter ( ptr2f(:,jmf:1:-1), local )            ! Sarah lu
523              var3d(:,:,k) = local(:,:)                                       ! Sarah Lu
524     
525              verb = .false. ! true only for k=1
526     
527     !        Apply mask when present
528     !        -----------------------
529              IF(doingMasking) THEN
530               WHERE(mask(i1:i2,j1:j2) == 0) var3d(i1:i2,j1:j2,k) = 0.00
531              END IF
532     
533             end do ! over k
534     
535           end if ! 2D or 3D
536     
537     !   Close file
538     !   ----------
539         if ( MAPL_am_I_root() ) then
540     
541     !      If binning, destroy transform
542     !      -----------------------------
543            if (im < imf .OR. jm < jmf .or. ForceBinning_ ) then
544               call MAPL_HorzTransformDestroy (Trans )
545            end if
546     
547            call GFIO_Close ( fid, rc )
548     !       print *, myname // ': Closing GFIO file ' // trim(fname)
549     
550            deallocate(ptr2, ptr2f, lonf, latf, levf, stat=ios)
551            if ( ios /= 0 ) call die ( myname, 'cannot deallocate ptr2/ptr2f' )
552     
553         end if ! masterproc
554     
555     !   All done
556     !   --------
557         deallocate(lat, lon, stat=ios)
558         if ( ios /= 0 ) call die ( myname, 'cannot deallocate lat/lon/lev')
559         IF(doingMasking) THEN
560          DEALLOCATE(regionNumbers, mask, STAT=ios)
561          IF ( ios /= 0 ) CALL die ( myname, ': Cannot deallocate masking tape.')
562         END IF
563         return
564     
565     CONTAINS
566     
567         subroutine Regrid_ (Gptr2file, im, jm, Gptr2bundle, im0, jm0, &
568                                 LONSfile, LONSbundle, verb ) 
569         integer, intent(in) :: im, jm, im0, jm0
570         real, pointer :: Gptr2bundle(:,:), Gptr2file(:,:)
571         real, pointer :: LONSfile(:), LONSbundle(:)
572         logical, intent(inout) :: verb
573     
574     !   Local buffers for cases compiled with -r8
575     !   -----------------------------------------
576         real(ESMF_KIND_R4) :: r4_Gptr2file(im,jm)
577         real(ESMF_KIND_R4) :: r4_Gptr2bundle(im0,jm0)
578     
579     !   Code borrowed from MAPL_CFIO - Note (im,jm) here differs from above
580     !   Notice that since ARCTAS unless ForceBinning has been specified as
581     !   .FALSE. *binning* is always the regridding  method.
582     !                             ---
583         integer :: STATUS
584     !ams    print *, 'Regrid_:  im, jm  (file)   = ', im, jm
585     !ams    print *, 'Regrid_: im0, jm0 (bundle) = ', im0, jm0
586         if ( im /= im0 .OR. jm /= JM0 ) then
587            r4_Gptr2file = Gptr2file
588            if (IM0 <  IM .or. JM0 < JM .or. ForceBinning_ ) then
589               if ( verb ) print *, myname // ': Binning... '
590               call MAPL_HorzTransformRun(Trans, r4_Gptr2file, r4_Gptr2bundle, MAPL_undef, rc=STATUS )
591               VERIFY_(STATUS)
592            else
593               if ( verb ) print *, myname // ': Interpolating... '
594               call hinterp ( r4_Gptr2file,im,jm, r4_Gptr2bundle,im0,jm0,1,MAPL_Undef )
595            end if ! coarsening
596            Gptr2bundle = r4_Gptr2bundle
597         else
598            Gptr2bundle = Gptr2file
599         end if ! change resolution
600         if ( abs(LONSbundle(1)-LONSfile(1)+180.) <= 10.*tiny(1.) ) then
601            if ( verb ) print *, myname // &
602                 ': shifting input longitudes by 180 degrees'
603            call shift180Lon2D_ ( Gptr2bundle, im0, jm0 )
604         end if
605         end subroutine Regrid_
606     
607     
608         subroutine shift180Lon2D_ ( c, im, jm )
609         integer, intent(in) :: im, jm
610         real, intent(inout) :: c(im,jm)
611         real :: cj(im)
612         integer :: m(4), n(4), imh, j
613         imh = nint(im/2.)
614         m = (/ 1,      imh, 1+imh,    im   /)
615         n = (/ 1,   im-imh, 1+im-imh, im   /)
616         do j = 1, jm
617            cj(n(1):n(2)) = c(m(3):m(4),j)
618            cj(n(3):n(4)) = c(m(1):m(2),j)
619            c(:,j) = cj
620         end do
621         return
622         end subroutine shift180Lon2D_
623     
624         subroutine queryFile_ ( fid, im, jm, km, lon, lat, lev, varn )
625           integer, intent(in)  :: fid
626           integer, intent(out) :: im, jm, km
627           character(len=255), intent(inout) :: varn
628           real, pointer :: lon(:), lat(:), lev(:)
629     !                    ----
630           character(len=255)              :: title, source, contact, levunits
631           character(len=255), allocatable :: vtitle(:), vunits(:)
632           character(len=257), allocatable :: vname(:)
633         
634           real,    allocatable :: valid_range(:,:), packing_range(:,:)
635           integer, allocatable :: kmvar(:), yyyymmdd(:), hhmmss(:)
636           integer :: lm, nvars, timinc, ngatts, n
637           real    :: amiss
638     
639           call GFIO_DimInquire ( fid, im, jm, km, lm, nvars, ngatts, rc )
640           if ( rc /= 0 ) call die(myname,'cannot get file dimensions')
641         
642           allocate ( lat(jm), lon(im), lev(km), yyyymmdd(lm), hhmmss(lm),    &
643                   vname(nvars), vunits(nvars), vtitle(nvars), kmvar(nvars), &
644                   valid_range(2,nvars), packing_range(2,nvars),             &
645                   stat=rc )
646           if ( rc /= 0 ) call die(myname,'cannot allocate file attribs' )
647     
648           call GFIO_Inquire ( fid, im, jm, km, lm, nvars,     &
649              title, source, contact, amiss,  &
650              lon, lat, lev, levunits,        &
651              yyyymmdd, hhmmss, timinc,       &
652              vname, vtitle, vunits, kmvar,   &
653              valid_range , packing_range, rc )
654           if ( rc /= 0 ) call die(myname,'cannot get file dimensions')
655         
656     !     Make variable names case insensitive
657     !     ------------------------------------
658           do n = 1, nvars
659              if ( uppercase(trim(varn)) .EQ. uppercase(trim(vname(n))) ) then
660                   varn = trim(vname(n))
661                   exit
662              end if
663           end do 
664     
665           deallocate ( yyyymmdd, hhmmss,                  &
666                       vname, vunits, vtitle, kmvar,      &
667                       valid_range, packing_range,        &
668                       stat=rc )
669           
670         end subroutine queryFile_
671     
672     end subroutine Chem_UtilMPread
673     
674     subroutine Chem_UtilGetTimeInfo ( fname, begDate, begTime, nTimes, incSecs )
675     
676       implicit NONE
677       character(len=*), intent(in) :: fname    ! GFIO/CFIO filename
678       integer, intent(out) :: begDate, begTime ! initial time/date on file
679                                                ! given as YYYYMMDD and HHMMSS
680       integer, intent(out) :: nTimes           ! number of time steps on file
681       integer, intent(out) :: incSecs          ! time steps in seconds
682      
683     
684       integer :: READ_ONLY=1
685       integer :: fid, rc, im,jm,km,nvars,ngatts
686       character(len=*), parameter ::  myname = 'Chem_UtilGetTimeInfo'
687       
688     ! Special case
689     ! -----------
690       if ( fname(1:9) == '/dev/null' ) then    
691          begDate=0; begTime=0;  nTimes=0; incSecs=0
692          return
693       end if
694     
695     ! Open file
696     ! ---------
697       call GFIO_Open ( fname, READ_ONLY, fid, rc )
698       if ( rc /= 0 ) call die(myname, 'Unable to open '// trim(fname))
699       if ( rc .ne. 0 ) then
700          begDate=-1; begTime=-1;  nTimes=-1; incSecs=-1
701          return
702       endif
703     
704     ! Get dimension sizes
705     ! -------------------
706       call GFIO_DimInquire (fid,im,jm,km,nTimes,nvars,ngatts,rc)
707       if ( rc /= 0 ) call die(myname, 'Unable to inquire about dimension sizes for file '// trim(fname))
708       if ( rc .ne. 0 ) then
709          begDate=-1; begTime=-1;  nTimes=-1; incSecs=-1
710          return
711       endif
712     
713     ! Get initial time/timestep
714     ! -------------------------
715       call GetBegDateTime ( fid, begDate, begTime, incSecs, rc )
716       if ( rc .ne. 0 ) then
717          begDate=-1; begTime=-1;  nTimes=-1; incSecs=-1
718          return
719       endif
720     
721       call GFIO_close(fid,rc)
722     
723     end subroutine Chem_UtilGetTimeInfo
724     
725     !............................... geos4 stubs  ..........................
726     
727     ! Parallelized utility routine for computing/printing
728     ! max/min of an input array
729     !
730           subroutine pmaxmin3d ( qname, a, pmin, pmax, im, jt, fac )
731           implicit none
732           character*(*)  qname
733           integer im, jt
734           real :: a(:,:,:)
735           real pmax, pmin
736           real fac                     ! multiplication factor
737           call pmaxmin2d ( qname, reshape(a,(/ im, jt /)), &
738                                  pmin, pmax, im, jt, fac )
739           
740           end subroutine pmaxmin3d
741     
742           subroutine pmaxmin2d ( qname, a, pmin, pmax, im, jt, fac )
743     
744           use parutilitiesmodule, only : commglobal, gid, maxop, parcollective
745     
746           implicit none
747     
748           character*(*)  qname
749           integer im, jt
750           real a(im,jt)
751           real pmax, pmin
752           real fac                     ! multiplication factor
753     
754           integer :: i, j, maxop_, two=2
755     
756           real qmin(jt), qmax(jt)
757           real pm1(2)
758     
759           character(len=16) :: name
760     
761     !$omp parallel do private(i, j, pmax, pmin)
762     
763           do j=1,jt
764              pmax = a(1,j)
765              pmin = a(1,j)
766              do i=2,im
767                 pmax = max(pmax, a(i,j))
768                 pmin = min(pmin, a(i,j))
769              enddo
770              qmax(j) = pmax
771              qmin(j) = pmin
772           enddo
773     !
774     ! Now find max/min of amax/amin
775     !
776           pmax = qmax(1)
777           pmin = qmin(1)
778           do j=2,jt
779              pmax = max(pmax, qmax(j))
780              pmin = min(pmin, qmin(j))
781           enddo
782     
783           pm1(1) = pmax
784           pm1(2) = -pmin
785           maxop_ = maxop
786           call parcollective(commglobal, maxop_, two, pm1  )
787           pmax=pm1(1)
788           pmin=-pm1(2)
789          
790           if ( fac /= 0.0 ) then  ! trick to prevent printing
791           if ( MAPL_am_I_root() ) then
792                name = '            '
793                name(1:len(qname)) = qname
794                write(*,*) name, ' max = ', pmax*fac, ' min = ', pmin*fac
795                return
796           end if
797           end if
798     
799         end subroutine pmaxmin2d
800     
801     
802           function leap_year(ny)
803     !
804     ! Determine if year ny is a leap year
805     !
806     ! Author: S.-J. Lin
807           implicit none
808           logical leap_year
809           integer ny
810           integer ny00
811     
812     !
813     ! No leap years prior to 0000
814     !
815           parameter ( ny00 = 0000 )   ! The threshold for starting leap-year 
816     
817           if( ny >= ny00 ) then
818              if( mod(ny,100) == 0. .and. mod(ny,400) == 0. ) then
819                  leap_year = .true.
820              elseif( mod(ny,4) == 0. .and. mod(ny,100) /= 0.  ) then
821                  leap_year = .true.
822              else
823                  leap_year = .false.
824              endif
825           else
826               leap_year = .false.
827           endif
828     
829           return 
830         end function leap_year
831     
832     
833           integer FUNCTION INCYMD (NYMD,M)
834     
835     !  PURPOSE
836     !     INCYMD:  NYMD CHANGED BY ONE DAY
837     !     MODYMD:  NYMD CONVERTED TO JULIAN DATE
838     !  DESCRIPTION OF PARAMETERS
839     !     NYMD     CURRENT DATE IN YYMMDD FORMAT
840     !     M        +/- 1 (DAY ADJUSTMENT)
841     
842           integer nymd, m, ny00, ny, nm, nd
843     
844           INTEGER NDPM(12)
845           DATA    NDPM /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/
846     !!!      logical leap_year
847           DATA    NY00     / 1900 /
848     
849           NY = NYMD / 10000
850           NM = MOD(NYMD,10000) / 100
851           ND = MOD(NYMD,100) + M
852     
853           IF (ND.EQ.0) THEN
854           NM = NM - 1
855           IF (NM.EQ.0) THEN
856               NM = 12
857               NY = NY - 1
858           ENDIF
859           ND = NDPM(NM)
860           IF (NM.EQ.2 .AND. leap_year(NY))  ND = 29
861           ENDIF
862     
863           IF (ND.EQ.29 .AND. NM.EQ.2 .AND. leap_year(ny))  GO TO 20
864     
865           IF (ND.GT.NDPM(NM)) THEN
866           ND = 1
867           NM = NM + 1
868           IF (NM.GT.12) THEN
869               NM = 1
870               NY = NY + 1
871           ENDIF
872           ENDIF
873     
874        20 CONTINUE
875           INCYMD = NY*10000 + NM*100 + ND
876           RETURN
877         END FUNCTION INCYMD
878     
879           subroutine tick (nymd, nhms, ndt)
880     
881     ! Input:
882           integer ndt                     ! TIME-STEP
883     ! Inpuit/Output:
884           integer nymd                    ! CURRENT YYYYMMDD
885           integer nhms                    ! CURRENT HHMMSS
886     !!!      integer incymd
887     
888     ! Revision:   S.-J. Lin Mar 2000
889            integer nsecf, nhmsf, n, nsec
890     
891            NSECF(N)   = N/10000*3600 + MOD(N,10000)/100* 60 + MOD(N,100)
892            NHMSF(N)   = N/3600*10000 + MOD(N,3600 )/ 60*100 + MOD(N, 60)
893     
894            NSEC = NSECF(NHMS) + ndt
895     
896            IF (NSEC.GT.86400)  THEN
897                DO WHILE (NSEC.GT.86400)
898                   NSEC = NSEC - 86400
899                   NYMD = INCYMD (NYMD,1)
900                ENDDO
901            ENDIF
902     
903            IF (NSEC.EQ.86400)  THEN
904                NSEC = 0
905                NYMD = INCYMD (NYMD,1)
906            ENDIF
907     
908            IF (NSEC .LT. 0)  THEN
909                DO WHILE (NSEC .LT. 0)
910                    NSEC = 86400 + NSEC
911                    NYMD = INCYMD (NYMD,-1)
912                ENDDO
913             ENDIF
914     
915               NHMS = NHMSF (NSEC)
916           return
917         end subroutine tick
918     
919           subroutine mcalday(nymd, nhms, calday)
920           implicit none
921     
922     ! input:
923           integer nymd
924           integer nhms
925     ! Output:
926           real calday                    ! Julian day (1 to 366 for non-leap year)
927                                          ! Julian day (-1 to -367 for   leap year)
928     ! Local:
929           logical leapyr                 ! leap_year?
930     !!!      logical leap_year
931           real tsec
932           integer n, nsecf, m, mm
933           integer dd, ds
934           integer days(12)
935           integer ny
936     
937           data days /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/
938           nsecf(n)  = n/10000*3600 + mod(n,10000)/100* 60 + mod(n,100)
939     
940           ny = nymd / 10000
941           mm = mod(nymd, 10000) / 100 
942           dd = mod(nymd,   100)
943     
944           ds = dd -1
945     
946           if( mm .ne. 1) then
947           do m=1, mm-1
948              if( m.eq.2  .and. leap_year(ny) ) then 
949                  ds = ds + 29
950              else
951                  ds = ds + days(m)
952              endif
953           enddo
954           endif
955     
956           tsec = ds * 86400 + nsecf(nhms)
957     
958           calday = tsec / 86400.  + 1.
959           if( leap_year(ny) ) calday = -calday
960     
961           return
962         end subroutine mcalday
963     
964           subroutine zenith(calday  ,dodiavg ,clat    ,coszrs  )
965     
966     !
967     ! Input arguments
968     !
969           real calday              ! Calendar day, including fraction
970           logical dodiavg          ! true => do diurnal averaging
971           real clat                ! Current latitude (radians)
972     !
973     ! Output arguments
974     !
975           real coszrs(*)       ! Cosine solar zenith angle
976     !
977     !---------------------------Local variables-----------------------------
978     !
979     
980           call die ('zenith','stub only, please do not call zenith_()' )
981     
982         end subroutine zenith
983     
984     !-------------------------------------------------------------------------
985     !     NASA/GSFC, Global Modeling and Assimilation Office, Code 900.3     !
986     !-------------------------------------------------------------------------
987     !BOP
988     !
989     ! !IROUTINE:  Chem_UtilNegFiller --- Negative Filler
990     !
991     ! !INTERFACE:
992     !
993        subroutine Chem_UtilNegFiller ( q, delp, in, jn, qmin )
994     
995     ! !USES:
996     
997       implicit NONE
998     
999     ! !INPUT/OUTPUT PARAMETERS:
1000     
1001       integer                    :: in, jn             ! number of local lon/lat 
1002       real, pointer              :: delp(:,:,:)
1003       real, OPTIONAL, intent(in) :: qmin
1004     
1005     ! !OUTPUT PARAMETERS:
1006     
1007       real, pointer :: q(:,:,:)     ! 3D tracer
1008     
1009     ! !DESCRIPTION: 
1010     !
1011     ! !REVISION HISTORY: Makes sure tracer has no negative values. This is
1012     !                    a "flat tax" algorithm: first negative values are
1013     !  replaced with tiny() or user specified value. Then profiles are rescaled 
1014     !  to preserve column mass, whenever possible. No mass conservation is
1015     !  imposed when the initial column mass is negative or zero.
1016     !
1017     !  18May2007  da Silva  First crack.
1018     !
1019     !EOP
1020     !-------------------------------------------------------------------------
1021     
1022       real :: mass_1(in,jn), mass_2(in,jn), ratio(in,jn)
1023       real :: qmin_, vmax, vmin
1024       integer :: k, k1, k2, km
1025     
1026       k1 = lbound(q,3)
1027       k2 = ubound(q,3)
1028       km = k2 - k1 + 1
1029     
1030     ! Unless specified, minimum is the smallest positive float, not zero
1031     ! ------------------------------------------------------------------
1032       if ( present(qmin) ) then
1033          qmin_ = qmin
1034       else  
1035          qmin_ = tiny(1.0)
1036       end if
1037     
1038     #ifdef DEBUG
1039       call pmaxmin ( 'NegFill:  q_beg', q, vmax, vmin, in*jn, km, 1. )
1040     #endif
1041     
1042     ! Column mass before fixer
1043     ! ------------------------
1044       mass_1 = sum ( delp * q, 3 )
1045     
1046     ! Cap q
1047     ! -----
1048       where ( q < qmin_ ) q = qmin_
1049     
1050     ! Enforce conservation of column mass
1051     ! -----------------------------------
1052       mass_2 = sum ( delp * q, 3 )
1053       where ( (mass_2 /= mass_1) .AND. (mass_1 > 0.0) )
1054               ratio = mass_1 / mass_2
1055       elsewhere
1056               ratio = 1.0
1057       end where
1058     
1059     ! Next correct q in each layer
1060     ! ----------------------------
1061       do k = k1, k2
1062          where ( ratio /= 1.0 )
1063                  q(:,:,k) = ratio * q(:,:,k)
1064          end where
1065       end do
1066     
1067     #ifdef DEBUG
1068       call pmaxmin ( 'NegFill: mass_1', mass_1, vmax, vmin, in*jn,1, 1. )
1069       call pmaxmin ( 'NegFill: mass_2', mass_2, vmax, vmin, in*jn,1, 1. )
1070       call pmaxmin ( 'NegFill:  ratio',  ratio, vmax, vmin, in*jn,1, 1. )
1071       call pmaxmin ( 'NegFill:  q_end', q, vmax, vmin, in*jn, km, 1. )
1072     #endif
1073     
1074     end subroutine Chem_UtilNegFiller
1075     
1076     !-------------------------------------------------------------------------
1077     !     NASA/GSFC, Global Modeling and Assimilation Office, Code 610.1     !
1078     !-------------------------------------------------------------------------
1079     !BOP
1080     !
1081     ! !IROUTINE:  Chem_UtilTroppFixer - Repair tropopause pressure bad values
1082     !
1083     ! !INTERFACE:
1084     !
1085       SUBROUTINE Chem_UtilTroppFixer(im, jm, tropp, threshold, verbose, newTropp, rc)
1086     
1087     ! !USES:
1088     
1089     ! !USES:
1090     
1091       IMPLICIT NONE
1092     
1093     ! !INPUT/OUTPUT PARAMETERS:
1094     
1095       INTEGER, INTENT(IN)		:: im          !Index range, longitude
1096       INTEGER, INTENT(IN)		:: jm          !Index range, latitude
1097       REAL, POINTER , INTENT(INOUT)	:: tropp(:,:)  !Tropopause pressure (Pa)
1098       REAL, OPTIONAL, INTENT(IN)    :: threshold   !User-supplied threshold (Pa).  If
1099                                                    ! not present, default is 110000 Pa.
1100       LOGICAL, OPTIONAL, INTENT(IN) :: verbose     !Write message when bad value is 
1101                                                    ! encountered. DEBUG directive turns
1102                                                    ! on the message even if verbose is not
1103                                                    ! present or if verbose = .FALSE.
1104       INTEGER, INTENT(OUT)		:: rc          !Return code: 0=OK, 1=Unable to repair
1105     
1106     
1107     ! !OUTPUT PARAMETERS:
1108     
1109       REAL, OPTIONAL, INTENT(OUT)   :: newTropp(1:im,1:jm)
1110                                                    !If present, fill with repaired
1111                                                    ! tropopause pressures (Pa).  Otherwise
1112                                                    ! tropp is overwritten.
1113     
1114     ! !DESCRIPTION: 
1115     !
1116     ! Replace bad values in tropopause pressure array with an average of nearby good
1117     ! values.  Working outward from the afflicted cell, nearby cells are scanned for 
1118     ! the presence of good values. In the first scan, all adjacent cells, which can
1119     ! number up to eight, are considered.  If at least one of the cells has a valid
1120     ! pressure, the scanning is terminated.  Otherwise, the "radius" of the scan is
1121     ! increased by one, and up to 24 cells are considered, and so on, until, at the
1122     ! extreme, all cells on the current processor fall under consideration.
1123     !
1124     ! After the scanning is done, the bad value is replaced with the average the 
1125     ! valid pressures found by the scan.  Thus, the accuracy of the replaced value
1126     ! drops rapidly as the scan expands outward.  At the same time, the only case
1127     ! in which a valid pressure will not be found to replace a bad one is if all 
1128     ! pressures on the current processor are invalid.  In this case the return 
1129     ! code is set to 1.
1130     !
1131     ! If newTropp is not present, then the input array, tropp, is overwritten.  If
1132     ! it is present, it is filled, even if there are no invalid pressures in tropp.
1133     !
1134     ! No scanning is done if all tropopause pressures are valid.
1135     !
1136     ! Assumptions/bugs:
1137     !
1138     ! Bad values are high, but not Indef, and the primary purpose of this routine
1139     ! is to repair the case where GEOS-5 fails to find the tropopause and assigns
1140     ! MAPL_UNDEF as the tropopause pressure.  We recommend using the blended 
1141     ! tropopause values for tropp, because the frequency of bad values is quite
1142     ! rare compared to the unblended case.
1143     !
1144     ! !REVISION HISTORY: 
1145     !
1146     !  25Jan2008  Nielsen  Initial code and testing.
1147     !
1148     !EOP
1149     !-------------------------------------------------------------------------
1150       CHARACTER(LEN=*), PARAMETER :: myName = 'Chem_UtilTroppFixer'
1151     
1152       LOGICAL :: tellMe
1153     
1154       INTEGER :: i,ier,j,l,m
1155       INTEGER :: ie,iw,jn,js
1156     
1157       INTEGER, ALLOCATABLE :: mask(:,:)
1158       INTEGER, ALLOCATABLE :: mx(:)
1159       REAL, ALLOCATABLE :: p(:,:)
1160       
1161       REAL :: badValue, r
1162     
1163       rc = 0
1164     
1165     ! Determine verbosity, letting the DEBUG 
1166     ! directive override local specification
1167     ! --------------------------------------
1168       tellMe = .FALSE.
1169       IF(PRESENT(verbose)) THEN
1170        IF(verbose) tellMe = .TRUE.
1171       END IF
1172     #ifdef DEBUG
1173       tellMe = .TRUE.
1174     #endif
1175     
1176     ! Set the bad value to 110000 Pa (1100 hPa)
1177     ! -----------------------------------------
1178       IF(PRESENT(threshold)) THEN
1179        badValue = threshold
1180       ELSE
1181        badValue = 1.10E+05 !Pa, 1100 hPa
1182       END IF
1183     
1184     ! There may be no bad values ...
1185     ! ------------------------------
1186       IF( ALL( tropp(1:im,1:jm) < badValue ) ) THEN
1187        IF(PRESENT(newTropp)) newTropp(1:im,1:jm) = tropp(1:im,1:jm)
1188        RETURN
1189       END IF
1190     
1191     ! ... or there is at least one bad value
1192     ! --------------------------------------
1193       ALLOCATE(mask(1:im,1:jm),STAT=ier)
1194       ALLOCATE(p(1:im,1:jm),STAT=ier)
1195       ALLOCATE(mx(4),STAT=ier)
1196     
1197     ! Loop over each cell
1198     ! -------------------
1199       DO j=1,jm
1200        DO i=1,im
1201     
1202     ! Invalid pressure found at cell(i,j)
1203     ! -----------------------------------
1204         IF(tropp(i,j) >= badValue) THEN
1205     
1206     ! Determine maximum "radius" of search
1207     ! ------------------------------------
1208          mx(1) = im-i
1209          mx(2) = i-1
1210          mx(3) = jm-j
1211          mx(4) = j-1
1212     
1213     ! Start search
1214     ! ------------
1215          DO m=1,MAXVAL(mx)
1216     
1217     ! Clear the mask
1218     ! --------------
1219           mask(1:im,1:jm) = 0
1220     
1221     ! Range of search
1222     ! ---------------
1223           iw = MAX( 1,i-m)
1224           ie = MIN(im,i+m)
1225           js = MAX( 1,j-m)
1226           jn = MIN(jm,j+m)
1227     
1228     ! Set mask to one for cells in range of search
1229     ! --------------------------------------------
1230            mask(iw:ie,js:jn) = 1
1231     
1232     ! Set mask back to zero for cells in range
1233     ! of search that have invalid pressures.
1234     ! ----------------------------------------
1235            WHERE(tropp(iw:ie,js:jn) >= badValue) mask(iw:ie,js:jn) = 0
1236     
1237     ! One valid pressure is enough ...
1238     ! --------------------------------
1239            IF(SUM(MASK) >= 1) EXIT
1240     
1241     ! ... or "radius" of search needs to be extended
1242     ! ----------------------------------------------
1243          END DO
1244     
1245     ! Repair bad value at cell(i,j) with average 
1246     ! of valid pressures found in range of search
1247     ! -------------------------------------------
1248          r = SUM(tropp,mask == 1)
1249          p(i,j) = r/(1.00*SUM(mask))
1250     
1251     ! For debugging
1252     ! -------------
1253          IF(tellMe) THEN
1254           WRITE(*,FMT="(A,': ',ES12.5,' becomes ',ES12.5,' Pa [',I4,2X,I4,']')") &
1255                 TRIM(myName),tropp(i,j),p(i,j),m,SUM(mask)
1256          END IF
1257     
1258         ELSE
1259     
1260     ! Input pressure at cell(i,j) was valid
1261     ! -------------------------------------
1262          p(i,j) = tropp(i,j)
1263     
1264         END IF
1265     
1266     ! Next cell
1267     ! ---------
1268        END DO
1269       END DO
1270     
1271     ! Clean up
1272     ! --------
1273       DEALLOCATE(mask,STAT=ier)
1274       DEALLOCATE(mx,STAT=ier)
1275     
1276     ! If all cells have bad values, then 
1277     ! register a failure, but continue.
1278     ! ----------------------------------
1279       IF( ANY( p(1:im,1:jm) >= badValue ) ) THEN
1280        PRINT *, myName,": WARNING Unable to fix bad tropopause pressure(s)"
1281        rc = 1
1282       END IF
1283       
1284     ! Overwrite input or fill output array
1285     ! ------------------------------------
1286       IF(PRESENT(newTropp)) THEN
1287        newTropp(1:im,1:jm) = p(1:im,1:jm)
1288       ELSE
1289        tropp(1:im,1:jm) = p(1:im,1:jm)
1290       END IF
1291     
1292     ! Clean up some more
1293     ! ------------------
1294       DEALLOCATE(p,STAT=ier)
1295     
1296       RETURN
1297       END SUBROUTINE Chem_UtilTroppFixer
1298     
1299     !-------------------------------------------------------------------------
1300     !     NASA/GSFC, Global Modeling and Assimilation Office, Code 610.1     !
1301     !-------------------------------------------------------------------------
1302     !BOP
1303     !
1304     ! !ROUTINE:  Chem_UtilExtractIntegers - Extract integers from a delimited string
1305     !
1306     ! !INTERFACE:
1307     !
1308       SUBROUTINE Chem_UtilExtractIntegers(string,iSize,iValues,delimiter,verbose,rc)
1309     
1310     ! !USES:
1311     
1312       IMPLICIT NONE
1313     
1314     ! !INPUT/OUTPUT PARAMETERS:
1315     
1316       CHARACTER(LEN=*), INTENT(IN)   :: string	  ! Character-delimited string of integers
1317       INTEGER, INTENT(IN)            :: iSize
1318       INTEGER, INTENT(INOUT)         :: iValues(iSize)! Space allocated for extracted integers
1319       CHARACTER(LEN=*), OPTIONAL     :: delimiter     ! 1-character delimiter
1320       LOGICAL, OPTIONAL, INTENT(IN)  :: verbose	  ! Let me know iValues as they are found. 
1321                                     		  ! DEBUG directive turns on the message even 
1322                                     		  ! if verbose is not present or if 
1323                                     		  ! verbose = .FALSE.
1324       INTEGER, OPTIONAL, INTENT(OUT) :: rc            ! Return code
1325     
1326     ! !DESCRIPTION: 
1327     !
1328     !  Extract integers from a character-delimited string, for example, "-1,45,256,7,10".  In the context
1329     !  of Chem_Util, this is provided for determining the numerically indexed regions over which an 
1330     !  emission might be applied.
1331     !
1332     !  In multiple passes, the string is parsed for the delimiter, and the characters up to, but not
1333     !  including the delimiter are taken as consecutive digits of an integer.  A negative sign ("-") is
1334     !  allowed.  After the first pass, each integer and its trailing delimiter are lopped of the head of
1335     !  the (local copy of the) string, and the process is started over.
1336     !
1337     !  The default delimiter is a comma (",").
1338     !
1339     !  "Unfilled" iValues are zero.
1340     !  
1341     !  Return codes:
1342     !  1 Zero-length string.
1343     !  2 iSize needs to be increased.
1344     !
1345     !  Assumptions/bugs:
1346     !
1347     !  A non-zero return code does not stop execution.
1348     !  Allowed numerals are: 0,1,2,3,4,5,6,7,8,9.
1349     !  A delimiter must be separated from another delimiter by at least one numeral.
1350     !  The delimiter cannot be a numeral or a negative sign.
1351     !  The character following a negative sign must be an allowed numeral.
1352     !  The first character must be an allowed numeral or a negative sign.
1353     !  The last character must be an allowed numeral.
1354     !  The blank character (" ") cannot serve as a delimiter.
1355     !
1356     !  Examples of strings that will work:
1357     !  "1"
1358     !  "-1"
1359     !  "-1,2004,-3"
1360     !  "1+-2+3"
1361     !  "-1A100A5"
1362     !
1363     !  Examples of strings that will not work:
1364     !  "1,--2,3"
1365     !  "1,,2,3"
1366     !  "1,A,3"
1367     !  "1,-,2"
1368     !  "1,2,3,4,"
1369     !  "+1"
1370     !  "1 3 6"
1371     !
1372     ! !REVISION HISTORY: 
1373     !
1374     !  29Feb2008  Nielsen  Initial code and testing.
1375     !
1376     !EOP
1377     !-------------------------------------------------------------------------
1378      CHARACTER(LEN=*), PARAMETER :: myName = 'Chem_UtilExtractIntegers'
1379     
1380      INTEGER :: base,count,i,iDash,last,lenStr
1381      INTEGER :: multiplier,pos,posDelim,sign
1382      CHARACTER(LEN=255) :: str
1383      CHARACTER(LEN=1) :: char,delimChar
1384      LOGICAL :: Done
1385      LOGICAL :: tellMe
1386     
1387     ! Initializations
1388     ! ---------------
1389      rc = 0
1390      count = 1
1391      Done = .FALSE.
1392      iValues(:) = 0
1393      base = ICHAR("0")
1394      iDash = ICHAR("-")
1395     
1396     ! Determine verbosity, letting the DEBUG 
1397     ! directive override local specification
1398     ! --------------------------------------
1399       tellMe = .FALSE.
1400       IF(PRESENT(verbose)) THEN
1401        IF(verbose) tellMe = .TRUE.
1402      END IF
1403     #ifdef DEBUG
1404       tellMe = .TRUE.
1405     #endif
1406     
1407     ! Check for zero-length string
1408     ! ----------------------------
1409      lenStr = LEN_TRIM(string)
1410      IF(lenStr == 0) THEN
1411       rc = 1
1412       PRINT *,myname,": ERROR - Found zero-length string."
1413       RETURN
1414      END IF
1415     
1416     ! Default delimiter is a comma
1417     ! ----------------------------
1418      delimChar = ","
1419      IF(PRESENT(delimiter)) delimChar(1:1) = delimiter(1:1)
1420     
1421     ! Work on a local copy
1422     ! --------------------
1423      str = TRIM(string)
1424     
1425     ! One pass for each delimited integer
1426     ! -----------------------------------
1427      Parse: DO
1428     
1429       lenStr = LEN_TRIM(str)
1430     
1431     ! Parse the string for the delimiter
1432     ! ----------------------------------
1433       posDelim = INDEX(TRIM(str),TRIM(delimChar))
1434       IF(tellMe) PRINT *,myname,": Input string is >",TRIM(string),"<"
1435     
1436     ! If the delimiter does not exist,
1437     ! one integer remains to be extracted.
1438     ! ------------------------------------
1439       IF(posDelim == 0) THEN
1440        Done = .TRUE.
1441        last = lenStr
1442       ELSE
1443        last = posDelim-1
1444       END IF
1445       multiplier = 10**last
1446     
1447     ! Examine the characters of this integer
1448     ! --------------------------------------
1449       Extract: DO pos=1,last
1450     
1451        char = str(pos:pos)
1452        i = ICHAR(char)
1453     
1454     ! Account for a leading "-"
1455     ! -------------------------
1456        IF(pos == 1) THEN 
1457         IF(i == iDash) THEN
1458          sign = -1
1459         ELSE
1460          sign = 1
1461         END IF
1462        END IF
1463     
1464     ! "Power" of 10 for this character
1465     ! --------------------------------
1466        multiplier = multiplier/10
1467     
1468        IF(pos == 1 .AND. sign == -1) CYCLE Extract
1469     
1470     ! Integer comes from remaining characters
1471     ! ---------------------------------------
1472        i = (i-base)*multiplier
1473        iValues(count) = iValues(count)+i
1474        IF(pos == last) THEN
1475         iValues(count) = iValues(count)*sign
1476         IF(tellMe) PRINT *,myname,":Integer number ",count," is ",iValues(count)
1477        END IF
1478       
1479       END DO Extract
1480     
1481       IF(Done) EXIT
1482     
1483     ! Lop off the leading integer and try again
1484     ! -----------------------------------------
1485       str(1:lenStr-posDelim) = str(posDelim+1:lenStr)
1486       str(lenStr-posDelim+1:255) = " "
1487       count = count+1
1488     
1489     ! Check size
1490     ! ----------
1491       IF(count > iSize) THEN
1492        rc = 2
1493        PRINT *,myname,": ERROR - iValues does not have enough elements."
1494       END IF
1495     
1496      END DO Parse
1497     
1498      RETURN
1499     END SUBROUTINE Chem_UtilExtractIntegers
1500     
1501     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1502     !
1503     ! This is a candidate for ESMFL, here for dependency reasons
1504     !  
1505     
1506       subroutine GridGetLatLons_ ( grid, lons, lats )
1507     
1508         implicit NONE
1509         type(ESMF_Grid) :: grid
1510     
1511         real, pointer   :: lons(:), lats(:)
1512     
1513     !                     ---
1514     
1515         character(len=*), parameter :: Iam = 'GridGetLatLons'
1516     
1517         real(KIND=8), pointer  :: R8D2(:,:)
1518         real, pointer          :: lons2d(:,:), lats2d(:,:)
1519         real, pointer          :: LONSLocal(:,:), LATSlocal(:,:)
1520         integer                :: IM_WORLD, JM_WORLD, dims(3), STATUS, RC
1521     
1522     !                          ----
1523     
1524     !      Get world dimensions
1525     !      --------------------
1526            call MAPL_GridGet(GRID, globalCellCountPerDim=DIMS, RC=STATUS)
1527     
1528            IM_WORLD = dims(1)
1529            JM_WORLD = dims(2)
1530     
1531     !      Allocate memory for output if necessary
1532     !      ---------------------------------------
1533            if ( .not. associated(lons) ) then
1534                 allocate(lons(IM_WORLD), stat=STATUS)
1535            else
1536                 if(size(LONS,1) /= IM_WORLD) STATUS = 1
1537            end if
1538            VERIFY_(status)
1539            if ( .not. associated(lats) ) then
1540                 allocate(lats(JM_WORLD), stat=STATUS)
1541            else
1542                 if(size(LATS,1) /= JM_WORLD) STATUS = 1
1543            end if
1544            VERIFY_(status)
1545     
1546     !      Local work space
1547     !      ----------------
1548            allocate(LONS2d(IM_WORLD,JM_WORLD), LATS2d(IM_WORLD,JM_WORLD), &
1549                     STAT=status)             
1550            VERIFY_(status)
1551            LONS2d=0
1552            LATS2d=0
1553     
1554     !      Get the local longitudes and gather them into a global array
1555     !      ------------------------------------------------------------
1556            call ESMF_GridGetCoord(grid, localDE=0, coordDim=1, &
1557                  staggerloc=ESMF_STAGGERLOC_CENTER, doCopy=ESMF_DATA_REF, &
1558                  fptr=R8D2, rc=status)
1559     
1560            allocate(LONSLOCAL(size(R8D2,1),size(R8D2,2)), STAT=status)             
1561            VERIFY_(status)
1562     
1563            LONSLOCAL = R8D2*(180/MAPL_PI)
1564     
1565            call ArrayGather(LONSLOCAL, LONS2D, GRID, RC=STATUS)
1566     
1567     !      Get the local longitudes and gather them into a global array
1568     !      ------------------------------------------------------------
1569            call ESMF_GridGetCoord(grid, localDE=0, coordDim=2, &
1570                  staggerloc=ESMF_STAGGERLOC_CENTER, doCopy=ESMF_DATA_REF, &
1571                  fptr=R8D2, rc=status)
1572     
1573            allocate(LATSLOCAL(size(R8D2,1),size(R8D2,2)), STAT=status)             
1574            VERIFY_(status)
1575     
1576            LATSlocal = R8D2*(180/MAPL_PI)
1577     
1578            call ArrayGather(LATSLOCAL, LATS2D, GRID, RC=STATUS)
1579            VERIFY_(STATUS)
1580     
1581     !      Return 1D arrays
1582     !      ----------------
1583            LONS = LONS2D(:,1)
1584            LATS = LATS2D(1,:)
1585     
1586            DEALLOCATE(LONSLOCAL, LATSLOCAL, LONS2d, LATS2d )
1587             
1588          end subroutine GridGetLatLons_
1589     
1590     !..........................................................................................
1591     
1592     
1593     !-------------------------------------------------------------------------
1594     !     NASA/GSFC, Global Modeling and Assimilation Office, Code 610.1     !
1595     !-------------------------------------------------------------------------
1596     !BOP
1597     !
1598     ! !ROUTINE:  Chem_UtilExtractIntegers - Extract integers from a delimited string
1599     !
1600     ! !INTERFACE:
1601     !
1602     
1603          subroutine Chem_BiomassDiurnal ( Eout, Ein, lons, lats, nhms, cdt)
1604     
1605     ! !USES:
1606     
1607       IMPLICIT NONE
1608     
1609     ! !ARGUMENTS:
1610     
1611            real, intent(out)   :: Eout(:,:) ! Emissions valid at NHMS
1612            real, intent(in)    :: Ein(:,:)  ! Daily-mean emissions
1613            real, intent(in)    :: lons(:)   ! Latitudes in degrees
1614            real, intent(in)    :: lats(:)   ! Latitudes in degrees
1615            integer, intent(in) :: nhms
1616            real, intent(in)    :: cdt       ! time step in seconds
1617     
1618     ! !DESCRIPTION: 
1619     !
1620     !      Applies diurnal cycle to biomass emissions.       
1621     !
1622     ! !DESCRIPTION:
1623     !
1624     !  This module implements assorted odds & ends for fvChem.
1625     !
1626     ! !REVISION HISTORY:
1627     !
1628     !  13nov2009  da Silva  First crack.
1629     !
1630     !EOP
1631     !-------------------------------------------------------------------------
1632     
1633     !      Hardwired diurnal cycle (multiplied by 100)
1634     !      These numbers were derived from GOES-12
1635     !      fire counts for 2003-2007.
1636     !      -------------------------------------------
1637            integer, parameter :: N = 240
1638            real,    parameter :: DT = 86400. / N
1639            real,    parameter :: Boreal(N) = &
1640            (/ 0.0277, 0.0292, 0.0306, 0.0318, 0.0327, 0.0335, &
1641               0.0340, 0.0342, 0.0341, 0.0338, 0.0333, 0.0326, &
1642               0.0316, 0.0305, 0.0292, 0.0278, 0.0263, 0.0248, &
1643               0.0233, 0.0217, 0.0202, 0.0187, 0.0172, 0.0158, &
1644               0.0145, 0.0133, 0.0121, 0.0110, 0.0100, 0.0091, &
1645               0.0083, 0.0075, 0.0068, 0.0062, 0.0056, 0.0051, &
1646               0.0046, 0.0042, 0.0038, 0.0035, 0.0032, 0.0030, &
1647               0.0028, 0.0026, 0.0025, 0.0024, 0.0024, 0.0024, &
1648               0.0024, 0.0026, 0.0027, 0.0030, 0.0033, 0.0036, &
1649               0.0041, 0.0046, 0.0052, 0.0060, 0.0069, 0.0079, &
1650               0.0090, 0.0104, 0.0119, 0.0137, 0.0157, 0.0180, &
1651               0.0205, 0.0235, 0.0268, 0.0305, 0.0346, 0.0393, &
1652               0.0444, 0.0502, 0.0565, 0.0634, 0.0711, 0.0794, &
1653               0.0884, 0.0982, 0.1087, 0.1201, 0.1323, 0.1453, &
1654               0.1593, 0.1742, 0.1900, 0.2069, 0.2249, 0.2439, &
1655               0.2642, 0.2858, 0.3086, 0.3329, 0.3587, 0.3860, &
1656               0.4149, 0.4455, 0.4776, 0.5115, 0.5470, 0.5840, &
1657               0.6227, 0.6628, 0.7043, 0.7470, 0.7908, 0.8355, &
1658               0.8810, 0.9271, 0.9735, 1.0200, 1.0665, 1.1126, &
1659               1.1580, 1.2026, 1.2460, 1.2880, 1.3282, 1.3664, &
1660               1.4023, 1.4356, 1.4660, 1.4933, 1.5174, 1.5379, &
1661               1.5548, 1.5679, 1.5772, 1.5826, 1.5841, 1.5818, &
1662               1.5758, 1.5661, 1.5529, 1.5365, 1.5169, 1.4944, &
1663               1.4693, 1.4417, 1.4119, 1.3801, 1.3467, 1.3117, &
1664               1.2755, 1.2383, 1.2003, 1.1616, 1.1225, 1.0832, &
1665               1.0437, 1.0044, 0.9653, 0.9265, 0.8882, 0.8504, &
1666               0.8134, 0.7771, 0.7416, 0.7070, 0.6734, 0.6407, &
1667               0.6092, 0.5787, 0.5493, 0.5210, 0.4939, 0.4680, &
1668               0.4433, 0.4197, 0.3974, 0.3763, 0.3565, 0.3380, &
1669               0.3209, 0.3051, 0.2907, 0.2777, 0.2662, 0.2561, &
1670               0.2476, 0.2407, 0.2352, 0.2313, 0.2289, 0.2279, &
1671               0.2283, 0.2300, 0.2329, 0.2369, 0.2417, 0.2474, &
1672               0.2536, 0.2602, 0.2670, 0.2738, 0.2805, 0.2869, &
1673               0.2927, 0.2979, 0.3024, 0.3059, 0.3085, 0.3101, &
1674               0.3107, 0.3102, 0.3087, 0.3061, 0.3026, 0.2983, &
1675               0.2931, 0.2871, 0.2806, 0.2735, 0.2659, 0.2579, &
1676               0.2497, 0.2412, 0.2326, 0.2240, 0.2153, 0.2066, &
1677               0.1979, 0.1894, 0.1809, 0.1726, 0.1643, 0.1562, &
1678               0.1482, 0.1404, 0.1326, 0.1250, 0.1175, 0.1101, &
1679               0.1028, 0.0956, 0.0886, 0.0818, 0.0751, 0.0687 /)       
1680            real,    parameter :: NonBoreal(N) = &
1681            (/ 0.0121, 0.0150, 0.0172, 0.0185, 0.0189, 0.0184, &
1682               0.0174, 0.0162, 0.0151, 0.0141, 0.0133, 0.0126, &
1683               0.0121, 0.0117, 0.0115, 0.0114, 0.0114, 0.0116, &
1684               0.0120, 0.0126, 0.0133, 0.0142, 0.0151, 0.0159, &
1685               0.0167, 0.0174, 0.0180, 0.0184, 0.0187, 0.0189, &
1686               0.0190, 0.0190, 0.0191, 0.0192, 0.0192, 0.0193, &
1687               0.0194, 0.0194, 0.0193, 0.0192, 0.0190, 0.0187, &
1688               0.0185, 0.0182, 0.0180, 0.0178, 0.0177, 0.0176, &
1689               0.0174, 0.0172, 0.0169, 0.0166, 0.0162, 0.0158, &
1690               0.0153, 0.0149, 0.0144, 0.0138, 0.0132, 0.0126, &
1691               0.0118, 0.0109, 0.0101, 0.0092, 0.0085, 0.0081, &
1692               0.0080, 0.0083, 0.0091, 0.0102, 0.0117, 0.0135, &
1693               0.0157, 0.0182, 0.0210, 0.0240, 0.0273, 0.0308, &
1694               0.0345, 0.0387, 0.0432, 0.0483, 0.0540, 0.0606, &
1695               0.0683, 0.0775, 0.0886, 0.1022, 0.1188, 0.1388, &
1696               0.1625, 0.1905, 0.2229, 0.2602, 0.3025, 0.3500, &
1697               0.4031, 0.4623, 0.5283, 0.6016, 0.6824, 0.7705, &
1698               0.8650, 0.9646, 1.0676, 1.1713, 1.2722, 1.3662, &
1699               1.4491, 1.5174, 1.5685, 1.6014, 1.6173, 1.6200, &
1700               1.6150, 1.6082, 1.6040, 1.6058, 1.6157, 1.6353, &
1701               1.6651, 1.7045, 1.7513, 1.8024, 1.8541, 1.9022, &
1702               1.9429, 1.9738, 1.9947, 2.0072, 2.0132, 2.0141, &
1703               2.0096, 1.9994, 1.9829, 1.9604, 1.9321, 1.8977, &
1704               1.8562, 1.8052, 1.7419, 1.6646, 1.5738, 1.4734, &
1705               1.3693, 1.2676, 1.1724, 1.0851, 1.0052, 0.9317, &
1706               0.8637, 0.8004, 0.7414, 0.6862, 0.6348, 0.5871, &
1707               0.5434, 0.5037, 0.4682, 0.4368, 0.4097, 0.3864, &
1708               0.3667, 0.3499, 0.3355, 0.3231, 0.3123, 0.3029, &
1709               0.2944, 0.2862, 0.2773, 0.2670, 0.2547, 0.2402, &
1710               0.2238, 0.2061, 0.1882, 0.1712, 0.1562, 0.1434, &
1711               0.1332, 0.1251, 0.1189, 0.1141, 0.1103, 0.1071, &
1712               0.1043, 0.1018, 0.0996, 0.0979, 0.0968, 0.0964, &
1713               0.0966, 0.0970, 0.0973, 0.0970, 0.0959, 0.0938, &
1714               0.0909, 0.0873, 0.0831, 0.0784, 0.0732, 0.0676, &
1715               0.0618, 0.0565, 0.0521, 0.0491, 0.0475, 0.0473, &
1716               0.0480, 0.0492, 0.0504, 0.0514, 0.0519, 0.0521, &
1717               0.0520, 0.0517, 0.0513, 0.0510, 0.0507, 0.0507, &
1718               0.0508, 0.0512, 0.0515, 0.0518, 0.0519, 0.0518, &
1719               0.0513, 0.0506, 0.0496, 0.0482, 0.0465, 0.0443, &
1720               0.0418, 0.0387, 0.0351, 0.0310, 0.0263, 0.0214 /)
1721     
1722     !      Fixed normalization factors; a more accurate normalization would take
1723     !      in consideration longitude and time step
1724     !      ---------------------------------------------------------------------
1725            real*8, save :: fBoreal = -1., fNonBoreal = -1
1726            real,   save :: fDT=-1
1727     
1728            integer :: hh, mm, ss, ndt, i, j, k
1729            integer :: i1, i2, j1, j2
1730            real :: secs, secs_local, aBoreal, aNonBoreal, alpha
1731     
1732     !                              -----
1733     
1734     !      Normalization factor depends on timestep
1735     !      ----------------------------------------
1736            if ( fDT /= cdt ) then
1737                 fBoreal = 0.0
1738                 fNonBoreal = 0.0
1739                 ndt = max(1,nint(cdt/DT))
1740                 do k = 1, N, ndt
1741                    fBoreal    = fboreal    + Boreal(k)
1742                    fNonBoreal = fNonBoreal + NonBoreal(k)
1743                 end do
1744                 fDT = cdt ! so it recalculates only if necessary
1745            end if
1746     
1747     
1748     !      Find number of secs since begining of the day (GMT)
1749     !      ---------------------------------------------------
1750            hh = nhms/10000
1751            mm = (nhms - 10000*hh) / 100
1752            ss = nhms - 10000*hh - 100*ss
1753            secs = 3600.*hh + 60.*mm + ss
1754     
1755     !      Apply factors depending on latitude
1756     !      -----------------------------------
1757            do i = lbound(Ein,1), ubound(Ein,1)
1758     
1759     !            Find corresponding index in hardwired diurnal cycle
1760     !            240 = 24 * 60 * 60 secs / 360 deg
1761     !            ---------------------------------------------------
1762                  secs_local = secs + 240. * lons(i)
1763                  k = 1 + mod(nint(secs_local/DT),N)
1764                  if ( k < 1 ) k = N + k
1765     
1766     !            Apply diurnal cycle
1767     !            -------------------
1768                  aBoreal = Boreal(k) / fBoreal
1769                  aNonBoreal = NonBoreal(k) / fNonBoreal
1770                  do j = lbound(Ein,2), ubound(Ein,2)
1771                     if ( lats(j) >= 50. ) then
1772                        Eout(i,j) = aBoreal    * Ein(i,j)
1773                     else if ( lats(j) >= 45. ) then
1774                        alpha = (lats(j) - 45. ) / 5.
1775                        Eout(i,j) = (1-alpha) * aNonBoreal * Ein(i,j) + &
1776                                       alpha  * aBoreal    * Ein(i,j)
1777                     else                  
1778                        Eout(i,j) = aNonBoreal * Ein(i,j)
1779                     end if
1780               end do
1781            end do
1782     
1783          end subroutine Chem_BiomassDiurnal
1784     
1785      end module Chem_UtilMod
1786     
1787