time_ave Program

Uses

  • program~~time_ave~~UsesGraph program~time_ave time_ave ESMF ESMF program~time_ave->ESMF gFTL_StringVector gFTL_StringVector program~time_ave->gFTL_StringVector ieee_arithmetic ieee_arithmetic program~time_ave->ieee_arithmetic iso_fortran_env iso_fortran_env program~time_ave->iso_fortran_env module~mapl MAPL program~time_ave->module~mapl module~mapl_filemetadatautilsmod MAPL_FileMetadataUtilsMod program~time_ave->module~mapl_filemetadatautilsmod mpi mpi program~time_ave->mpi module~esmf_cfiomod ESMF_CFIOMod module~mapl->module~esmf_cfiomod module~mapl_esmffieldbundleread MAPL_ESMFFieldBundleRead module~mapl->module~mapl_esmffieldbundleread module~mapl_esmffieldbundlewrite MAPL_ESMFFieldBundleWrite module~mapl->module~mapl_esmffieldbundlewrite module~mapl_fieldutils MAPL_FieldUtils module~mapl->module~mapl_fieldutils module~mapl_genericmod MAPL_GenericMod module~mapl->module~mapl_genericmod module~mapl_gridcompsmod MAPL_GridCompsMod module~mapl->module~mapl_gridcompsmod module~mapl_openmp_support MAPL_OpenMP_Support module~mapl->module~mapl_openmp_support module~mapl_profiler mapl_Profiler module~mapl->module~mapl_profiler module~mapl_stubcomponent MAPL_StubComponent module~mapl->module~mapl_stubcomponent module~mapl_varspecmiscmod MAPL_VarSpecMiscMod module~mapl->module~mapl_varspecmiscmod module~maplbase_mod MAPLBase_Mod module~mapl->module~maplbase_mod module~pfio pFIO module~mapl->module~pfio module~mapl_filemetadatautilsmod->ESMF module~mapl_filemetadatautilsmod->iso_fortran_env gFTL2_StringIntegerMap gFTL2_StringIntegerMap module~mapl_filemetadatautilsmod->gFTL2_StringIntegerMap module~mapl_abstractgridfactorymod MAPL_AbstractGridFactoryMod module~mapl_filemetadatautilsmod->module~mapl_abstractgridfactorymod module~mapl_exceptionhandling MAPL_ExceptionHandling module~mapl_filemetadatautilsmod->module~mapl_exceptionhandling module~mapl_gridmanagermod MAPL_GridManagerMod module~mapl_filemetadatautilsmod->module~mapl_gridmanagermod module~mapl_keywordenforcermod mapl_KeywordEnforcerMod module~mapl_filemetadatautilsmod->module~mapl_keywordenforcermod module~mapl_filemetadatautilsmod->module~pfio

Calls

time_avewadd_new_field_to_bundle
w
wcheck_quad
w
wcompute_nhmsf
w
wcompute_nsecf
w
wcopy_bundle_to_bundle
w
wcreate_output_grid
w
wdefined
w
wDistributedProfiler%reduce
w
wDistributedProfiler%start
w
wESMF_CalendarSetDefault
w
wESMF_ClockCreate
w
wESMF_ClockSet
w
wESMF_ConfigCreate
w
wESMF_ConfigFindLabel
w
wESMF_ConfigGetAttribute
w
wESMF_ConfigLoadFile
w
wESMF_ConfigNextLine
w
wesmf_fieldbundlecreate
w
wesmf_fieldbundleget
w
wesmf_fieldbundleset
w
wesmf_fieldget
w
wESMF_InfoGet
w
wESMF_InfoGetFromHost
w
wESMF_Initialize
w
wESMF_TimeIntervalSet
w
wFieldBundleWriter%create_from_bundle
w
wFieldBundleWriter%start_new_file
w
wFieldBundleWriter%write_to_file
w
wFileMetadataUtils%create
w
wFileMetadataUtils%get_coordinate_info
w
wFileMetadataUtils%get_level_name
w
wFileMetadataUtils%var_get_missing_value
w
wgenerate_report
w
wget_file_levels
w
wget_file_times
w
wget_level_info
w
wget_long_names
w
wget_units
w
wGridManager%make_factory
w
wGridManager%make_grid
w
whas_level
w
wlatlon_zstar
w
wlocal_esmf_timeset
w
wMAPL_abort
w
wMAPL_Assert
w
wMAPL_FieldBundleAdd
w
wMAPL_Finalize
w
wMAPL_GridGet
w
wMAPL_Initialize
w
wMAPL_read_bundle
w
wMAPL_Verify
w
wmpi_barrier
w
wmpi_comm_rank
w
wmpi_comm_size
w
wmpi_finalize
w
wmpi_init
w
wmpi_reduce
w
wncagt
w
wServerManager%finalize
w
wServerManager%initialize
w
wtick
w
wusage
w

Variables

Type Attributes Name Initial
character(len=1) :: char
character(len=2) :: time0
character(len=8) :: date0
character(len=:), allocatable :: lev_name
character(len=:), allocatable :: output_grid_name
character(len=ESMF_MAXPATHLEN) :: arg_str
character(len=ESMF_MAXSTR), allocatable :: aliases(:,:)
character(len=ESMF_MAXSTR), allocatable :: aliastmp(:,:)
character(len=ESMF_MAXSTR) :: doutput
character(len=ESMF_MAXSTR) :: dummy
character(len=ESMF_MAXSTR) :: ext
character(len=ESMF_MAXSTR), allocatable :: fname(:)
character(len=ESMF_MAXSTR) :: grid_type
character(len=ESMF_MAXSTR) :: hdfile
character(len=ESMF_MAXSTR) :: lev_units
character(len=ESMF_MAXSTR) :: name
character(len=ESMF_MAXSTR) :: name1
character(len=ESMF_MAXSTR) :: name2
character(len=ESMF_MAXSTR) :: name3
character(len=ESMF_MAXSTR) :: output
character(len=ESMF_MAXSTR), allocatable :: quadratics(:,:)
character(len=ESMF_MAXSTR), allocatable :: quadtmp(:,:)
character(len=ESMF_MAXSTR) :: rcfile
character(len=ESMF_MAXSTR) :: template
character(len=ESMF_MAXSTR), allocatable :: vname(:)
character(len=ESMF_MAXSTR), allocatable :: vname2(:)
character(len=ESMF_MAXSTR), allocatable :: vtitle(:)
character(len=ESMF_MAXSTR), allocatable :: vtitle2(:)
character(len=ESMF_MAXSTR), allocatable :: vunits(:)
character(len=ESMF_MAXSTR), allocatable :: vunits2(:)
integer :: L
integer :: begdate
integer :: begtime
integer :: comm
integer :: enddate
integer :: endtime
integer :: global_dims(3)
integer, allocatable :: hhmmss(:)
integer :: i
integer :: id
integer :: ierror
integer, allocatable :: iloc(:)
integer :: im
integer :: imglobal
integer :: j
integer :: jm
integer :: jmglobal
integer :: k
integer :: kbeg
integer :: kend
integer :: km
integer, allocatable :: kmvar(:)
integer, allocatable :: kmvar2(:)
integer :: lm
integer :: loc1
integer :: loc2
integer :: local_dims(3)
integer :: m
integer :: month
integer :: monthm
integer :: monthp
integer :: mv
integer :: mvars
integer :: myid
integer :: n
integer :: n_times
integer :: nalias
integer :: nargs
integer :: nc
integer :: ncvid
integer :: ndt
integer :: ndvars
integer :: nfiles
integer :: nhms
integer :: nhms0
integer :: nhmsm
integer :: nhmsp
integer, allocatable :: nloc(:)
integer :: nmax
integer :: npes
integer :: nquad
integer :: nstar
integer :: ntcrit
integer :: ntime
integer, allocatable :: ntimes(:,:,:,:)
integer :: ntmin
integer :: ntod
integer :: ntods
integer :: nv
integer :: nvars
integer :: nvars2
integer :: nymd
integer :: nymd0
integer :: nymdm
integer :: nymdp
integer, allocatable :: qloc(:,:)
integer :: rc
integer :: status
integer :: timeid
integer :: timeinc
integer :: timinc
integer :: year
integer :: yearm
integer :: yearp
integer, allocatable :: yymmdd(:)
logical :: allow_zonal_means
logical :: diurnal
logical :: file_has_lev
logical :: first
logical :: ignore_nan
logical :: ldquad
logical :: lquad
logical, allocatable :: lzstar(:)
logical :: mdiurnal
logical :: root
logical :: strict
logical :: tend
real(kind=ESMF_KIND_R4), pointer :: ptr2d(:,:)
real(kind=ESMF_KIND_R4), pointer :: ptr3d(:,:,:)
real(kind=REAL64), allocatable :: q(:,:,:,:)
real, allocatable :: dum(:,:,:)
real, allocatable :: dumz1(:,:)
real, allocatable :: dumz2(:,:)
real, allocatable :: lev(:)
real :: plev
real :: previous_undef
real, allocatable :: qmax(:)
real :: qmaxg
real, allocatable :: qmin(:)
real :: qming
real :: undef
type(ESMF_Clock) :: clock
type(ESMF_Config) :: config
type(ESMF_Field) :: field
type(ESMF_Grid) :: input_grid
type(ESMF_Grid) :: output_grid
type(ESMF_Info) :: infoh
type(ESMF_Time) :: etime
type(ESMF_Time), allocatable :: time_series(:)
type(ESMF_TimeInterval) :: time_interval
class(AbstractGridFactory), allocatable :: factory
type(DistributedProfiler), target :: t_prof
type(ESMF_FieldBundle) :: diurnal_bundle
type(ESMF_FieldBundle) :: final_bundle
type(ESMF_FieldBundle) :: primary_bundle
type(FieldBundleWriter) :: diurnal_writer
type(FieldBundleWriter) :: standard_writer
type(FileMetadata) :: basic_metadata
type(FileMetadataUtils) :: file_metadata
type(NetCDF4_FileFormatter) :: file_handle
type(ProfileReporter) :: reporter
type(ServerManager) :: io_server
type(verticalData) :: vertical_data

Functions

function get_long_names(bundle, rc) result(long_names)

Arguments

Type IntentOptional Attributes Name
type(ESMF_FieldBundle), intent(in) :: bundle
integer, intent(out), optional :: rc

Return Value character(len=ESMF_MAXSTR), allocatable, (:)

function get_units(bundle, rc) result(units)

Arguments

Type IntentOptional Attributes Name
type(ESMF_FieldBundle), intent(in) :: bundle
integer, intent(out), optional :: rc

Return Value character(len=ESMF_MAXSTR), allocatable, (:)

function compute_incymd(nymd, m) result(incymd)

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nymd
integer, intent(in) :: m

Return Value integer

function compute_nhmsf(nsec) result(nhmsf)

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nsec

Return Value integer

function compute_nsecf(nhms) result(seconds)

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nhms

Return Value integer

function get_level_info(bundle, rc) result(kmvar)

Arguments

Type IntentOptional Attributes Name
type(ESMF_FieldBundle), intent(in) :: bundle
integer, intent(out), optional :: rc

Return Value integer, allocatable, (:)

function defined(q, undef)

Arguments

Type IntentOptional Attributes Name
real :: q
real :: undef

Return Value logical

function has_level(grid, rc) result(grid_has_level)

Arguments

Type IntentOptional Attributes Name
type(ESMF_Grid), intent(in) :: grid
integer, intent(out), optional :: rc

Return Value logical

function is_leap_year(year)

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: year

Return Value logical

function create_output_grid(grid_name, lm, rc) result(new_grid)

Arguments

Type IntentOptional Attributes Name
character(len=*), intent(inout) :: grid_name
integer, intent(in) :: lm
integer, intent(out), optional :: rc

Return Value type(ESMF_Grid)

function local_esmf_timeset(yymmdd, hhmmss, rc) result(etime)

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: yymmdd
integer, intent(in) :: hhmmss
integer, intent(out), optional :: rc

Return Value type(ESMF_Time)


Subroutines

subroutine add_new_field_to_bundle(bundle, grid, lm, field_name, long_name, units, rc)

Arguments

Type IntentOptional Attributes Name
type(ESMF_FieldBundle), intent(inout) :: bundle
type(ESMF_Grid), intent(in) :: grid
integer, intent(in) :: lm
character(len=*), intent(in) :: field_name
character(len=*), intent(in) :: long_name
character(len=*), intent(in) :: units
integer, intent(out), optional :: rc

subroutine check_quad(quad, vname, nvars, aliases, nalias, qloc)

Arguments

Type IntentOptional Attributes Name
character(len=ESMF_MAXSTR) :: quad(2)
character(len=ESMF_MAXSTR) :: vname(nvars)
integer :: nvars
character(len=ESMF_MAXSTR) :: aliases(2,nalias)
integer :: nalias
integer :: qloc(2)

subroutine copy_bundle_to_bundle(input_bundle, output_bundle, rc)

Arguments

Type IntentOptional Attributes Name
type(ESMF_FieldBundle), intent(inout) :: input_bundle
type(ESMF_FieldBundle), intent(inout) :: output_bundle
integer, intent(out), optional :: rc

subroutine generate_report()

Arguments

None

subroutine get_esmf_grid_layout(grid, nx, ny, ims_out, jms_out, rc)

Arguments

Type IntentOptional Attributes Name
type(ESMF_Grid), intent(inout) :: grid
integer, intent(out) :: nx
integer, intent(out) :: ny
integer, intent(inout), allocatable :: ims_out(:)
integer, intent(inout), allocatable :: jms_out(:)
integer, intent(out), optional :: rc

subroutine get_file_levels(filename, vertical_data, rc)

Arguments

Type IntentOptional Attributes Name
character(len=*), intent(in) :: filename
type(verticalData), intent(inout) :: vertical_data
integer, intent(out), optional :: rc

subroutine get_file_times(file_metadata, num_times, time_series, time_interval, yymmdd, hhmmss, rc)

Arguments

Type IntentOptional Attributes Name
type(FileMetadataUtils), intent(inout) :: file_metadata
integer, intent(out) :: num_times
type(ESMF_Time), intent(inout), allocatable :: time_series(:)
integer, intent(out) :: time_interval
integer, intent(inout), allocatable :: yymmdd(:)
integer, intent(inout), allocatable :: hhmmss(:)
integer, intent(out), optional :: rc

subroutine latlon_zmean(q, qz, undef, grid, rc)

Arguments

Type IntentOptional Attributes Name
real, intent(inout) :: q(:,:)
real, intent(inout) :: qz(:)
real, intent(in) :: undef
type(ESMF_Grid), intent(inout) :: grid
integer, intent(out), optional :: rc

subroutine latlon_zstar(q, qp, undef, grid, rc)

Arguments

Type IntentOptional Attributes Name
real, intent(inout) :: q(:,:)
real, intent(out) :: qp(:,:)
real, intent(in) :: undef
type(ESMF_Grid), intent(inout) :: grid
integer, intent(out), optional :: rc

subroutine tick(nymd, nhms, ndt)

Arguments

Type IntentOptional Attributes Name
integer, intent(inout) :: nymd
integer, intent(inout) :: nhms
integer, intent(in) :: ndt

subroutine usage(root)

Arguments

Type IntentOptional Attributes Name
logical, intent(in) :: root

Source Code

program  time_ave

   use ESMF
   use MAPL
   use MAPL_FileMetadataUtilsMod
   use gFTL_StringVector
   use MPI
   use, intrinsic :: iso_fortran_env, only: int32, int64, int16, real32, real64
   use ieee_arithmetic, only: isnan => ieee_is_nan

   implicit none

   integer  comm,myid,npes,ierror
   integer  imglobal
   integer  jmglobal
   logical  root

! **********************************************************************
! **********************************************************************
! ****                                                              ****
! ****        Program to create time-averaged HDF files             ****
! ****                                                              ****
! **********************************************************************
! **********************************************************************

   integer  im,jm,lm

   integer  nymd, nhms
   integer  nymd0,nhms0
   integer  nymdp,nhmsp
   integer  nymdm,nhmsm
   integer  ntod, ndt, ntods
   integer  month,  year
   integer  monthp, yearp
   integer  monthm, yearm
   integer  begdate, begtime
   integer  enddate, endtime

   integer id,rc,timeinc,timeid
   integer ntime,nvars,ncvid,nvars2

   character(len=ESMF_MAXSTR), allocatable :: fname(:)
   character(len=ESMF_MAXSTR)  template
   character(len=ESMF_MAXSTR)  name
   character(len=ESMF_MAXSTR)  ext
   character(len=ESMF_MAXSTR)  output, doutput, hdfile, rcfile
   character(len=8)    date0
   character(len=2)    time0
   character(len=1)    char
   data output  /'monthly_ave'/
   data rcfile  /'NULL'/
   data doutput /'NULL'/
   data template/'NULL'/

   integer n,m,nargs,L,nfiles,nv,km,mvars,mv,ndvars

   real     plev,qming,qmaxg
   real     previous_undef,undef
   real,    allocatable ::    lev(:)
   integer, allocatable ::  kmvar(:)  ,  kmvar2(:)
   integer, allocatable :: yymmdd(:)
   integer, allocatable :: hhmmss(:)
   integer, allocatable ::   nloc(:)
   integer, allocatable ::   iloc(:)

   character(len=ESMF_MAXSTR), allocatable ::  vname(:),  vname2(:)
   character(len=ESMF_MAXSTR), allocatable :: vtitle(:), vtitle2(:)
   character(len=ESMF_MAXSTR), allocatable :: vunits(:), vunits2(:)

   real,    allocatable ::   qmin(:)
   real,    allocatable ::   qmax(:)
   real,    allocatable ::  dumz1(:,:)
   real,    allocatable ::  dumz2(:,:)
   real,    allocatable ::    dum(:,:,:)
   real(REAL64),  allocatable ::  q(:,:,:,:)
   integer, allocatable :: ntimes(:,:,:,:)

   integer timinc,i,j,k,nmax,kbeg,kend,loc1,loc2
   integer nstar
   logical tend, first, strict, diurnal, mdiurnal, lquad, ldquad
   logical ignore_nan
   data first  /.true./
   data strict /.true./

   type(ESMF_Config)   :: config

   integer,       allocatable ::       qloc(:,:)
   character(len=ESMF_MAXSTR), allocatable :: quadratics(:,:)
   character(len=ESMF_MAXSTR), allocatable ::    quadtmp(:,:)
   character(len=ESMF_MAXSTR), allocatable ::    aliases(:,:)
   character(len=ESMF_MAXSTR), allocatable ::   aliastmp(:,:)
   character(len=ESMF_MAXSTR)  name1, name2, name3, dummy
   integer        nquad
   integer        nalias
   logical,       allocatable :: lzstar(:)

   integer ntmin, ntcrit, nc

   type(FileMetadata) :: basic_metadata
   type(FileMetadataUtils) :: file_metadata
   type(NetCDF4_FileFormatter) :: file_handle
   integer :: status
   class(AbstractGridfactory), allocatable :: factory
   type(ESMF_Grid) :: output_grid,input_grid
   character(len=:), allocatable :: output_grid_name
   integer :: global_dims(3), local_dims(3)
   type(ESMF_Time), allocatable :: time_series(:)
   type(ESMF_TIme) :: etime
   type(ESMF_Clock) :: clock
   type(ESMF_TimeInterval) :: time_interval
   type(ESMF_FieldBundle) :: primary_bundle,final_bundle,diurnal_bundle
   type(ESMF_Field) :: field
   type(ServerManager) :: io_server
   type(FieldBundleWriter) :: standard_writer, diurnal_writer
   real(ESMF_KIND_R4), pointer :: ptr2d(:,:),ptr3d(:,:,:)
   character(len=ESMF_MAXSTR) :: grid_type
   logical :: allow_zonal_means
   character(len=ESMF_MAXPATHLEN) :: arg_str
   character(len=:), allocatable :: lev_name
   character(len=ESMF_MAXSTR) :: lev_units
   integer :: n_times
   type(verticalData) :: vertical_data
   logical :: file_has_lev
   type(DistributedProfiler), target :: t_prof
   type(ProfileReporter) :: reporter
   type(ESMF_Info) :: infoh

! **********************************************************************
! ****                       Initialization                         ****
! **********************************************************************

!call timebeg ('main')

   call mpi_init                ( ierror ) 
   _VERIFY(ierror)
   comm = mpi_comm_world
   call mpi_comm_rank ( comm,myid,ierror )
   _VERIFY(ierror)
   call mpi_comm_size ( comm,npes,ierror )
   _VERIFY(ierror)
   call ESMF_Initialize(logKindFlag=ESMF_LOGKIND_NONE,mpiCommunicator=MPI_COMM_WORLD, _RC)
   call MAPL_Initialize(_RC)
   t_prof = DistributedProfiler('time_ave_util',MpiTImerGauge(),MPI_COMM_WORLD)
   call t_prof%start(_RC)
   call io_server%initialize(MPI_COMM_WORLD,_RC)
   root = myid.eq.0
   call ESMF_CalendarSetDefault(ESMF_CALKIND_GREGORIAN,_RC)

! Read Command Line Arguments
! ---------------------------
   begdate = -999
   begtime = -999
   enddate = -999
   endtime = -999
   ndt = -999
   ntod = -999
   ntmin = -999
   nargs = command_argument_count()
   if( nargs.eq.0 ) then
      call usage(root)
   else
      lquad     = .TRUE.
      ldquad    = .FALSE.
      diurnal   = .FALSE.
      mdiurnal   = .FALSE.
      ignore_nan = .FALSE.
      do n=1,nargs
         call get_command_argument(n,arg_str)
         select case(trim(arg_str))
         case('-template')
            call get_command_argument(n+1,template)
         case('-tag')
            call get_command_argument(n+1,output)
         case('-rc')
            call get_command_argument(n+1,rcfile)
         case('-begdate')
            call get_command_argument(n+1,arg_str)
            read(arg_str,*)begdate
         case('-begtime')
            call get_command_argument(n+1,arg_str)
            read(arg_str,*)begtime
         case('-enddate')
            call get_command_argument(n+1,arg_str)
            read(arg_str,*)enddate
         case('-endtime')
            call get_command_argument(n+1,arg_str)
            read(arg_str,*)endtime
         case('-ntmin')
            call get_command_argument(n+1,arg_str)
            read(arg_str,*)ntmin
         case('-ntod')
            call get_command_argument(n+1,arg_str)
            read(arg_str,*)ntod
         case('-ndt')
            call get_command_argument(n+1,arg_str)
            read(arg_str,*)ndt
         case('-strict')
            call get_command_argument(n+1,arg_str)
            read(arg_str,*)strict
         case('-ogrid')
            call get_command_argument(n+1,arg_str)
            output_grid_name = trim(arg_str)
         case('-noquad')
            lquad = .FALSE.
         case('-ignore_nan')
            ignore_nan = .TRUE.
         case('-d')
            diurnal = .true.
            if (n+1 .le. nargs) then
               call get_command_argument(n+1,arg_str)
               read(arg_str,fmt='(a1)') char
               if (char.ne.'-') doutput=arg_str
            end if
         case('-md')
            mdiurnal = .true.
            if (n+1 .le. nargs) then
               call get_command_argument(n+1,arg_str)
               read(arg_str,fmt='(a1)') char
               if (char.ne.'-') doutput=arg_str
            end if
         case('-dv')
            ldquad = .true.
            diurnal = .true.
            if (n+1 .le. nargs) then
               call get_command_argument(n+1,arg_str)
               read(arg_str,fmt='(a1)') char
               if (char.ne.'-') doutput=arg_str
            end if
         case('-mdv')
            ldquad = .true.
            mdiurnal = .true.
            if (n+1 .le. nargs) then
               call get_command_argument(n+1,arg_str)
               read(arg_str,fmt='(a1)') char
               if (char.ne.'-') doutput=arg_str
            end if
         case('-eta')
            nfiles = 1
            call get_command_argument(n+nfiles,arg_str)
            read(arg_str,fmt='(a1)') char
            do while (char .ne. '-' .and. n+nfiles.ne.nargs)
               nfiles = nfiles + 1
               call get_command_argument(n+nfiles,arg_str)
               read(arg_str,fmt='(a1)') char
            enddo
            if (char.eq.'-') nfiles = nfiles-1
            allocate(fname(nfiles))
            do m=1,nfiles
               call get_command_argument(n+m,fname(m))
            enddo
         case('-hdf')
            nfiles = 1
            call get_command_argument(n+nfiles,arg_str)
            read(arg_str,fmt='(a1)') char
            do while (char .ne. '-' .and. n+nfiles.ne.nargs)
               nfiles = nfiles + 1
               call get_command_argument(n+nfiles,arg_str)
               read(arg_str,fmt='(a1)') char
            enddo
            if (char.eq.'-') nfiles = nfiles-1
            allocate(fname(nfiles))
            do m=1,nfiles
               call get_command_argument(n+m,fname(m))
            enddo
         end select
      enddo
   end if

   if( (diurnal.or.mdiurnal) .and. trim(doutput).eq.'NULL' ) then
      doutput = trim(output) // "_diurnal"
      if( mdiurnal ) diurnal = .FALSE.
   endif

   if (root .and. ignore_nan) print *,' ignore nan is true'


! Read RC Quadratics
! ------------------
   if( trim(rcfile).eq.'NULL' ) then
      nquad  = 0
      nalias = 0
   else
      config = ESMF_ConfigCreate    ( rc=rc )
      call   ESMF_ConfigLoadFile  ( config, trim(rcfile),  rc=rc )
      call   ESMF_ConfigFindLabel ( config, 'QUADRATICS:', rc=rc )
      tend = .false.
      m = 0
      do while (.not.tend)
         m = m+1
         allocate( quadtmp(3,m) )
         call ESMF_ConfigGetAttribute ( config,value=name1,   rc=rc )
         call ESMF_ConfigGetAttribute ( config,value=dummy,   rc=rc )
         call ESMF_ConfigGetAttribute ( config,value=name2,   rc=rc )
         call ESMF_ConfigGetAttribute ( config,value=dummy,   rc=rc )
         call ESMF_ConfigGetAttribute ( config,value=name3,default='XXX',rc=rc )
         call ESMF_ConfigNextLine     ( config,tableEnd=tend, rc=rc )
         if( m==1 ) then
            quadtmp(1,m)     = name1
            quadtmp(2,m)     = name2
            quadtmp(3,m)     = name3
            allocate( quadratics(3,m) )
            quadratics = quadtmp
         else
            quadtmp(1,1:m-1) = quadratics(1,:)
            quadtmp(2,1:m-1) = quadratics(2,:)
            quadtmp(3,1:m-1) = quadratics(3,:)
            quadtmp(1,m)     = name1
            quadtmp(2,m)     = name2
            quadtmp(3,m)     = name3
            deallocate( quadratics      )
            allocate( quadratics(3,m) )
            quadratics = quadtmp
         endif
         deallocate (quadtmp)
      enddo
      nquad = m

! Read RC Aliases
! ---------------
      call   ESMF_ConfigFindLabel ( config, 'ALIASES:', rc=rc )
      tend = .false.
      m = 0
      do while (.not.tend)
         m = m+1
         allocate( aliastmp(2,m) )
         call ESMF_ConfigGetAttribute ( config,value=name1,   rc=rc )
         call ESMF_ConfigGetAttribute ( config,value=dummy,   rc=rc )
         call ESMF_ConfigGetAttribute ( config,value=name2,   rc=rc )
         call ESMF_ConfigNextLine     ( config,tableEnd=tend  ,rc=rc )
         if( m==1 ) then
            aliastmp(1,m)     = name1
            aliastmp(2,m)     = name2
            allocate( aliases(2,m) )
            aliases = aliastmp
         else
            aliastmp(1,1:m-1) = aliases(1,:)
            aliastmp(2,1:m-1) = aliases(2,:)
            aliastmp(1,m)     = name1
            aliastmp(2,m)     = name2
            deallocate( aliases      )
            allocate( aliases(2,m) )
            aliases = aliastmp
         endif
         deallocate (aliastmp)
      enddo
      nalias = m
   endif
   if (.not. allocated(aliases)) allocate(aliases(0,0))

! **********************************************************************
! ****                        Read HDF File                         ****
! **********************************************************************

   call t_prof%start('initialize')

   if( trim(template).ne.'NULL' ) then
      name = template
   else
      name = fname(1)
   endif

   n = index(trim(name),'.',back=.true.)
   ext = trim(name(n+1:))

   call file_handle%open(trim(name),PFIO_READ,_RC)
   basic_metadata = file_handle%read(_RC)
   call file_handle%close(_RC)

   allocate(factory, source=grid_manager%make_factory(trim(name)))
   input_grid = grid_manager%make_grid(factory)
   file_has_lev = has_level(input_grid,_RC)
   call MAPL_GridGet(input_grid,globalCellCountPerDim=global_dims,_RC)
   lm = global_dims(3)

   if (file_has_lev) then
      call get_file_levels(trim(name),vertical_data,_RC)
   end if

   if (allocated(output_grid_name)) then
      output_grid = create_output_grid(output_grid_name,lm,_RC)
   else
      output_grid = input_grid
   end if
   call ESMF_InfoGetFromHost(output_grid,infoh,_RC)
   call ESMF_InfoGet(infoh,'GridType',grid_type,_RC)
   allow_zonal_means = trim(grid_type) == 'LatLon'
   if (trim(grid_type) == "Cubed-Sphere") then
      _ASSERT(mod(npes,6)==0,"If input files are Cubed-Sphere, must be run on multiple of 6 proccessors")
   end if
   call MAPL_GridGet(output_grid,localCellCountPerDim=local_dims,globalCellCountPerDim=global_dims,_RC)
   im = local_dims(1)
   jm = local_dims(2)
   lm = local_dims(3)
   imglobal = global_dims(1)
   jmglobal = global_dims(2)

   call file_metadata%create(basic_metadata,trim(name))
   call get_file_times(file_metadata,ntime,time_series,timinc,yymmdd,hhmmss,_RC)
   primary_bundle = ESMF_FieldBundleCreate(name="first_file",_RC)
   call ESMF_FieldBundleSet(primary_bundle,grid=output_grid,_RC)
   call MAPL_Read_Bundle(primary_bundle,trim(name),time=time_series(1),_RC)
   call ESMF_FieldBundleGet(primary_bundle,fieldCount=nvars,_RC)
   allocate(vname(nvars))
   call ESMF_FieldBundleGet(primary_bundle,fieldNameList=vname,_RC)
   kmvar = get_level_info(primary_bundle,_RC)
   vtitle = get_long_names(primary_bundle,_RC)
   vunits = get_units(primary_bundle,_RC)

   final_bundle = ESMF_FieldBundleCreate(name="first_file",_RC)
   call ESMF_FieldBundleSet(final_bundle,grid=output_grid,_RC)
   diurnal_bundle = ESMF_FieldBundleCreate(name="first_file",_RC)
   call ESMF_FieldBundleSet(diurnal_bundle,grid=output_grid,_RC)
   call copy_bundle_to_bundle(primary_bundle,final_bundle,_RC)

   if (size(time_series)>1) then
      time_interval = time_series(2) - time_series(1)
   else if (size(time_series)==1) then
      call ESMF_TimeIntervalSet(time_interval,h=6,_RC)
   end if
   clock = ESMF_ClockCreate(startTime=time_series(1),timeStep=time_interval,_RC)

   nvars2 = nvars

   if (file_has_lev) then
      lev_name = file_metadata%get_level_name(_RC)
      call file_metadata%get_coordinate_info(lev_name,coords=lev,coordUnits=lev_units,_RC)
   end if

   previous_undef = file_metadata%var_get_missing_value(trim(vname(1)),_RC)
   do i=2,size(vname)
      undef = file_metadata%var_get_missing_value(trim(vname(i)),_RC)
      _ASSERT(undef == previous_undef,"conflicting undefined values in your variables")
      previous_undef = undef
   enddo
   undef = previous_undef


! Set NDT for Strict Time Testing
! -------------------------------
   if( ntod.ne.-999 ) ndt = 86400
   if( ndt .eq.-999 ) ndt = compute_nsecf (timinc)
   if( timinc .eq. 0 ) then
      timeId = ncvid (id, 'time', rc)
      call ncagt     (id, timeId, 'time_increment', timinc, rc)
      if( timinc .eq. 0 ) then
         if( root ) then
            print *
            print *, 'Warning, GFIO Inquire states TIMINC = ',timinc
            print *, '         This will be reset to 060000 '
            print *, '         Use -ndt NNN (in seconds) to overide this'
         endif
         timinc = 060000
      endif
      ndt = compute_nsecf (timinc)
   endif

! Determine Number of Time Periods within 1-Day
! ---------------------------------------------
   ntods = 0
   if( diurnal .or. mdiurnal ) then
      if( ndt.lt.86400 ) ntods = 86400/ndt
   endif

! Set Minimum Required Times for Time Average (Default: 10 Days for Monthly Mean)
! -------------------------------------------------------------------------------
   if( ntmin.eq.-999 ) then
      if( ntod.eq.-999 ) then
         ntcrit = 10 * ( 86400.0/real(compute_nsecf(timinc)) )
      else
         ntcrit = 10
      endif
   else
      ntcrit = ntmin
   endif

! Determine Location Index for Each Variable in File
! --------------------------------------------------
   if( root ) print *
   allocate ( nloc(nvars) )
   nloc(1) = 1
   if( root ) write(6,7000) 1,trim(vname(1)),nloc(1),trim(vtitle(1)),max(1,kmvar(1))
   do n=2,nvars
      nloc(n) = nloc(n-1)+max(1,kmvar(n-1))
      if( root ) write(6,7000) n,trim(vname(n)),nloc(n),trim(vtitle(n)),max(1,kmvar(n))
7000  format(1x,'Primary Field:  ',i4,'  Name: ',a12,'  at location: ',i4,3x,a40,2x,i2,3x,i2,3x,i2)
   enddo

   nmax  = nloc(nvars)+max(1,kmvar(nvars))-1
   allocate( dum  (im,jm,nmax) )
   allocate( dumz1(im,jm) )
   allocate( dumz2(im,jm) )

! Append Default Quadratics to User-Supplied List
! -----------------------------------------------
   if( lquad ) then
      if( nquad.eq.0 ) then
         allocate( quadratics(3,nvars) )
         do n=1,nvars
            quadratics(1,n) = trim( vname(n) )
            quadratics(2,n) = trim( vname(n) )
            quadratics(3,n) = 'XXX'
         enddo
         nquad = nvars
      else
         allocate( quadtmp(3,nquad+nvars) )
         quadtmp(1,1:nquad) = quadratics(1,:)
         quadtmp(2,1:nquad) = quadratics(2,:)
         quadtmp(3,1:nquad) = quadratics(3,:)
         do n=1,nvars
            quadtmp(1,nquad+n) = trim( vname(n) )
            quadtmp(2,nquad+n) = trim( vname(n) )
            quadtmp(3,nquad+n) = 'XXX'
         enddo
         nquad = nquad + nvars
         deallocate( quadratics   )
         allocate( quadratics(3,nquad) )
         quadratics = quadtmp
         deallocate( quadtmp )
      endif
   endif

   allocate ( qloc(2,nquad) )
   allocate ( lzstar(nquad) ) ; lzstar = .FALSE.

! Determine Possible Quadratics
! -----------------------------
   km=kmvar(nvars)
   m=      nvars
   do n=1,nquad
      call check_quad ( quadratics(1,n),vname,nvars,aliases,nalias,qloc(1,n) )
      if( qloc(1,n)*qloc(2,n).ne.0 ) then
         m=m+1
         allocate ( iloc(m) )
         iloc(1:m-1) = nloc
         iloc(m) = iloc(m-1)+max(1,km)
         deallocate ( nloc )
         allocate ( nloc(m) )
         nloc = iloc
         deallocate ( iloc )
         km=kmvar( qloc(1,n) )
      endif
   enddo

   mvars = m
   nmax  = nloc(m)+max(1,km)-1

   allocate (  vname2(  mvars) )
   allocate ( vtitle2(  mvars) )
   allocate ( vunits2(  mvars) )
   allocate (  kmvar2(  mvars) )

   vname2(  1:nvars) = vname
   vtitle2(  1:nvars) = vtitle
   vunits2(  1:nvars) = vunits
   kmvar2(  1:nvars) = kmvar

   if( root .and. mvars.gt.nvars ) print *
   mv=  nvars
   do nv=1,nquad
      if( qloc(1,nv)*qloc(2,nv).ne.0 ) then
         mv  = mv+1

         if( trim(quadratics(1,nv)).eq.trim(quadratics(2,nv)) ) then
            vname2(mv) = "Var_"         //  trim(vname(qloc(1,nv)))
            vtitle2(mv) = "Variance_of_" //  trim(vname(qloc(1,nv)))
         else
            vname2(mv) = "Cov_"            // trim(vname(qloc(1,nv))) // "_"     // trim(vname(qloc(2,nv)))
            vtitle2(mv) = "Covariance_of_"  // trim(vname(qloc(1,nv))) // "_and_" // trim(vname(qloc(2,nv)))
         endif

         if( trim(quadratics(3,nv)).ne.'XXX' ) vname2(mv) = trim(quadratics(3,nv))

         nstar = index( trim(quadratics(1,nv)),'star',back=.true. )
         if( nstar.ne.0 ) then
            _ASSERT(allow_zonal_means,"grid is not lat-lon so cannot compute zonal means")
            lzstar(nv) = .TRUE.
            vtitle2(mv) = "Product_of_Zonal_Mean_Deviations_of_" // trim(vname(qloc(1,nv))) // "_and_" // trim(vname(qloc(2,nv)))
         endif

         vunits2(mv) = trim(vunits(qloc(1,nv))) // " " // trim(vunits(qloc(2,nv)))
         kmvar2(mv) =  kmvar(qloc(1,nv))

         call add_new_field_to_bundle(final_bundle,output_grid,kmvar(qloc(1,nv)),vname2(mv),vtitle2(mv),vunits2(mv),_RC)

         if( root ) write(6,7001) mv,trim(vname2(mv)),nloc(mv),trim(vtitle2(mv)),max(1,kmvar(qloc(1,nv))),qloc(1,nv),qloc(2,nv)
7001     format(1x,'   Quad Field:  ',i4,'  Name: ',a12,'  at location: ',i4,3x,a50,2x,i2,3x,i3,3x,i3)
      endif
   enddo

!deallocate ( lev )
   deallocate ( yymmdd )
   deallocate ( hhmmss )
   deallocate (  vname )
   deallocate ( vtitle )
   deallocate ( vunits )
   deallocate (  kmvar )

   allocate(   qmin(nmax) )
   allocate(   qmax(nmax) )
   allocate(      q(im,jm,nmax,0:ntods) )
   allocate( ntimes(im,jm,nmax,0:ntods) )
   ntimes = 0
   q = 0
   qmin =  abs(undef)
   qmax = -abs(undef)

   if( root ) then
      print *
      write(6,7002) mvars,nmax,im,jm,nmax,ntods
7002  format(1x,'Total Number of Variables: ',i3,/ &
            1x,'Total Size: ',i5,/ &
            1x,'Allocating q(',i4,',',i3,',',i5,',0:',i2.2,')')
      print *
      print *, 'Files: '
      do n=1,nfiles
         print *, n,trim(fname(n))
      enddo
      print *
      if( ntod.eq.-999 ) then
         print *, 'Averging Time-Period NHMS: ',ntod,' (ALL Possible Time Periods Used)'
      else
         print *, 'Averging Time-Period NHMS: ',ntod
      endif
      if( begdate.ne.-999 .or. begtime.ne.-999 ) print *, 'Beginning Date for Averaging: ',begdate,begtime
      if( enddate.ne.-999 .or. endtime.ne.-999 ) print *, '   Ending Date for Averaging: ',enddate,endtime
      if( strict ) then
         print *, 'Every Time Period Required for Averaging, STRICT = ',strict
      else
         print *, 'Only Averaging Time Periods Supplied, STRICT = ',strict
      endif
      write(6,7003) ntcrit
7003  format(1x,'Required Minimum Number of Defined Time Periods: ',i3,' (Otherwise, UNDEF)')
      print *
   endif

  call t_prof%stop('initialize')

! **********************************************************************
! ****                      Read HDF Files                          ****
! **********************************************************************

   k = 0

   do n=1,nfiles

      if (allocated(time_series)) deallocate(time_series)
      if (allocated(yymmdd)) deallocate(yymmdd)
      if (allocated(hhmmss)) deallocate(hhmmss)
      call file_handle%open(trim(fname(n)),PFIO_READ,_RC)
      basic_metadata = file_handle%read(_RC)
      call file_handle%close(_RC)
      call file_metadata%create(basic_metadata,trim(fname(n)))
      call get_file_times(file_metadata,ntime,time_series,timinc,yymmdd,hhmmss,_RC)


      do m=1,ntime
         nymd = yymmdd(m)
         nhms = hhmmss(m)
         if( nhms<0 ) then
            nhms = compute_nhmsf( compute_nsecf(nhms) + 86400 )
            call tick (nymd,nhms,-86400)
         endif

         if( ( begdate.ne.-999 .and. begtime.ne.-999 ) .and. &
               ( begdate.gt.nymd .or. &
               ( begdate.eq.nymd.and.begtime.gt.nhms ) ) ) cycle

         if( ( enddate.ne.-999 .and. endtime.ne.-999 ) .and. &
               ( enddate.lt.nymd .or. &
               ( enddate.eq.nymd.and.endtime.lt.nhms ) ) ) cycle

         k = k+1
         if( k.gt.ntods ) k = 1
         if( ntod.eq.-999 .or. ntod.eq.nhms ) then
            if( root ) write(6,3000) nymd,nhms,timinc,trim(fname(n)),k
3000        format(1x,'Reading nymd: ',i8.8,' nhms: ',i6.6,' TimInc: ',i6.6,'  from File: ',a,'  tod = ',i2)
            year  =     nymd/10000
            month = mod(nymd,10000)/100

! Check for Correct First Dataset
! -------------------------------
            if( strict .and. first ) then
               nymdm = nymd
               nhmsm = nhms
               call tick (nymdm,nhmsm,-ndt)
               yearm  =     nymdm/10000
               monthm = mod(nymdm,10000)/100
               if( year.eq.yearm .and. month.eq.monthm ) then
                  if( root ) print *, 'Date: ',nymd,' Time: ',nhms,' is NOT correct First Time Period!'
                  _FAIL("error processing dataset")
               endif
            endif

! Check Date and Time for STRICT Time Testing
! -------------------------------------------
            if( strict .and. .not.first ) then
               if( nymd.ne.nymdp .or. nhms.ne.nhmsp ) then
                  if( root ) print *, 'Date: ',nymdp,' Time: ',nhmsp,' not found!'
                  _FAIL("error processing dataset")
               endif
            endif
            nymdp = nymd
            nhmsp = nhms

! Primary Fields
! --------------

            etime = local_esmf_timeset(nymd,nhms,_RC)
            call MAPL_Read_Bundle(primary_bundle,trim(fname(1)),time=etime,file_override=trim(fname(n)),_RC)
            do nv=1,nvars2
               call ESMF_FieldBundleGet(primary_bundle,trim(vname2(nv)),field=field,_RC)
               call t_prof%start('PRIME')
               if( kmvar2(nv).eq.0 ) then
                  kbeg = 0
                  kend = 1
                  call ESMF_FieldGet(field,0,farrayPtr=ptr2d,_RC)
                  dum(:,:,nloc(nv))=ptr2d
               else
                  kbeg = 1
                  kend = kmvar2(nv)

                  call ESMF_FieldGet(field,0,farrayPtr=ptr3d,_RC)
                  dum(:,:,nloc(nv):nloc(nv)+kmvar2(nv)-1) = ptr3d
               endif

               rc = 0
               do L=1,max(1,kmvar2(nv))
                  do j=1,jm
                     do i=1,im
                        if( isnan( dum(i,j,nloc(nv)+L-1) ) .or. ( dum(i,j,nloc(nv)+L-1).gt.HUGE(dum(i,j,nloc(nv)+L-1)) ) ) then
!print *, 'Warning!  Nan or Infinity detected for ',trim(vname2(nv)),' at lat: ',lattice%jglobal(j),' lon: ',lattice%iglobal(i)
                           if( root .and. ignore_nan ) then
                              print *, 'Setting Nan or Infinity to UNDEF'
                              print *
                           else
                              rc = 1
                           endif
                           dum(i,j,nloc(nv)+L-1) = undef
                        endif
                        if( defined(dum(i,j,nloc(nv)+L-1),undef) ) then
                           q(i,j,nloc(nv)+L-1,0) =      q(i,j,nloc(nv)+L-1,0) + dum(i,j,nloc(nv)+L-1)
                           ntimes(i,j,nloc(nv)+L-1,0) = ntimes(i,j,nloc(nv)+L-1,0) + 1
                           if(   qmin(nloc(nv)+L-1).gt.dum(i,j,nloc(nv)+L-1) ) qmin(nloc(nv)+L-1) = dum(i,j,nloc(nv)+L-1)
                           if(   qmax(nloc(nv)+L-1).lt.dum(i,j,nloc(nv)+L-1) ) qmax(nloc(nv)+L-1) = dum(i,j,nloc(nv)+L-1)
                           if( ntods.ne.0 ) then
                              q(i,j,nloc(nv)+L-1,k) =      q(i,j,nloc(nv)+L-1,k) + dum(i,j,nloc(nv)+L-1)
                              ntimes(i,j,nloc(nv)+L-1,k) = ntimes(i,j,nloc(nv)+L-1,k) + 1
                           endif
                        endif
                     enddo
                  enddo
               enddo
               call t_prof%stop('PRIME')

            enddo

! Quadratics
! ----------
            call t_prof%start('QUAD')
            mv=  nvars2
            do nv=1,nquad
               if( qloc(1,nv)*qloc(2,nv).ne.0 ) then
                  mv=mv+1
                  do L=1,max(1,kmvar2(qloc(1,nv)))
                     if( lzstar(nv) ) then
                        call latlon_zstar (dum(:,:,nloc(qloc(1,nv))+L-1),dumz1,undef,output_grid,_RC)
                        call latlon_zstar (dum(:,:,nloc(qloc(2,nv))+L-1),dumz2,undef,output_grid,_RC)
                        do j=1,jm
                           do i=1,im
                              if( defined(dumz1(i,j),undef)  .and. &
                                    defined(dumz2(i,j),undef) ) then
                                 q(i,j,nloc(mv)+L-1,0) =      q(i,j,nloc(mv)+L-1,0) + dumz1(i,j)*dumz2(i,j)
                                 ntimes(i,j,nloc(mv)+L-1,0) = ntimes(i,j,nloc(mv)+L-1,0) + 1
                                 if( ntods.ne.0 ) then
                                    q(i,j,nloc(mv)+L-1,k) =      q(i,j,nloc(mv)+L-1,k) + dumz1(i,j)*dumz2(i,j)
                                    ntimes(i,j,nloc(mv)+L-1,k) = ntimes(i,j,nloc(mv)+L-1,k) + 1
                                 endif
                              endif
                           enddo
                        enddo
                     else
                        do j=1,jm
                           do i=1,im
                              if( defined(dum(i,j,nloc(qloc(1,nv))+L-1),undef)  .and. &
                                    defined(dum(i,j,nloc(qloc(2,nv))+L-1),undef) ) then
                                 q(i,j,nloc(mv)+L-1,0) = q(i,j,nloc(mv)+L-1,0) + dum(i,j,nloc(qloc(1,nv))+L-1) &
                                       * dum(i,j,nloc(qloc(2,nv))+L-1)
                                 ntimes(i,j,nloc(mv)+L-1,0) = ntimes(i,j,nloc(mv)+L-1,0) + 1
                                 if( ntods.ne.0 ) then
                                    q(i,j,nloc(mv)+L-1,k) = q(i,j,nloc(mv)+L-1,k) + dum(i,j,nloc(qloc(1,nv))+L-1) &
                                          * dum(i,j,nloc(qloc(2,nv))+L-1)
                                    ntimes(i,j,nloc(mv)+L-1,k) = ntimes(i,j,nloc(mv)+L-1,k) + 1
                                 endif
                              endif
                           enddo
                        enddo
                     endif
                  enddo
               endif
            enddo
            call t_prof%stop('QUAD')

            if( first ) then
               nymd0 = nymd
               nhms0 = nhms
               first = .false.
            endif

! Update Date and Time for Strict Test
! ------------------------------------
            call tick (nymdp,nhmsp,ndt)
            yearp  =     nymdp/10000
            monthp = mod(nymdp,10000)/100

         endif ! End ntod  Test
      enddo    ! End ntime Loop within file

      call MPI_BARRIER(comm,status)
      _VERIFY(status)
   enddo

   do k=0,ntods
      if( k.eq.0 ) then
         nc = ntcrit
      else
         nc = max( 1,ntcrit/ntods )
      endif
      do n=1,nmax
         do j=1,jm
            do i=1,im
               if( ntimes(i,j,n,k).lt.nc ) then
                  q(i,j,n,k) = undef
               else
                  q(i,j,n,k) = q(i,j,n,k)/ntimes(i,j,n,k)
               endif
            enddo
         enddo
      enddo
   enddo

! **********************************************************************
! ****               Write HDF Monthly Output File                  ****
! **********************************************************************

call t_prof%start('Write_AVE')

! Check for Correct Last Dataset
! ------------------------------
   if( strict .and. ( year.eq.yearp .and. month.eq.monthp ) ) then
      if( root ) print *, 'Date: ',nymd,' Time: ',nhms,' is NOT correct Last Time Period!'
      _FAIL("Error processing dataset")
   endif

   write(date0,4000) nymd0/100
   write(time0,2000) nhms0/10000

   hdfile = trim(output) // "." // trim(date0) // "." // trim(ext)

1000 format(i8.8)
2000 format(i2.2)
4000 format(i6.6)

   timeinc   = 060000

! Primary Fields
! --------------
   if( root ) print *
   do n=1,nvars2
      call ESMF_FieldBundleGet(final_bundle,trim(vname2(n)),field=field,_RC)
      if( kmvar2(n).eq.0 ) then
         kbeg = 0
         kend = 1
         call ESMF_FieldGet(field,0,farrayPtr=ptr2d,_RC)
         ptr2d = q(:,:,nloc(n),0)
      else
         kbeg = 1
         kend = kmvar2(n)
         call ESMF_FieldGet(field,0,farrayPtr=ptr3d,_RC)
         ptr3d = q(:,:,nloc(n):nloc(n)+kend-1,0)
      endif
      if( root ) write(6,3001) trim(vname2(n)),nloc(n),trim(hdfile)
3001  format(1x,'Writing ',a,' at location ',i6,' into File: ',a)
      dum(:,:,1:kend) = q(:,:,nloc(n):nloc(n)+kend-1,0)
   enddo

! Quadratics
! ----------
   mv=  nvars2
   do nv=1,nquad
      if( qloc(1,nv)*qloc(2,nv).ne.0 ) then
         mv=mv+1
         if( root ) write(6,3001) trim(vname2(mv)),nloc(mv),trim(hdfile)
         call ESMF_FieldBundleGet(final_bundle,trim(vname2(mv)),field=field,_RC)

         if( kmvar2(qloc(1,nv)).eq.0 ) then
            kbeg = 0
            kend = 1
         else
            kbeg = 1
            kend = kmvar2(qloc(1,nv))
         endif
         loc1 = nloc( qloc(1,nv) )
         loc2 = nloc( qloc(2,nv) )
         if( .not.lzstar(nv) ) then
            where( q(:,:,nloc(mv):nloc(mv)+kend-1,0).ne.undef )
               dum(:,:,1:kend) = q(:,:,nloc(mv):nloc(mv)+kend-1,0) - q(:,:,loc1:loc1+kend-1,0) &
                     * q(:,:,loc2:loc2+kend-1,0)
            elsewhere
               dum(:,:,1:kend) = undef
            endwhere
         else
            dum(:,:,1:kend) = q(:,:,nloc(mv):nloc(mv)+kend-1,0)
         endif
         if( kmvar2(qloc(1,nv)).eq.0 ) then
            call ESMF_FieldGet(field,0,farrayPtr=ptr2d,_RC)
            ptr2d = dum(:,:,1)
         else
            kend = kmvar2(qloc(1,nv))
            call ESMF_FieldGet(field,0,farrayPtr=ptr3d,_RC)
            ptr3d = dum(:,:,1:kend)
         endif
      endif
   enddo

   if( root ) then
      print *
      print *, 'Created: ',trim(hdfile)
      print *
   endif
   call t_prof%stop('Write_AVE')
   etime = local_esmf_timeset(nymd0,nhms0,_RC)
   call ESMF_ClockSet(clock,currTime=etime, _RC)
   call standard_writer%create_from_bundle(final_bundle,clock,n_steps=1,time_interval=timeinc,vertical_data=vertical_data,_RC)
   call standard_writer%start_new_file(trim(hdfile),_RC)
   call standard_writer%write_to_file(_RC)

! **********************************************************************
! ****               Write HDF Monthly Diurnal Output File          ****
! **********************************************************************

   if( ntods.ne.0 ) then
      call t_prof%start('Write_Diurnal')
      timeinc   = compute_nhmsf( 86400/ntods )

      do k=1,ntods

         if( k.eq.1 .or. mdiurnal ) then

            write(date0,4000) nymd0/100
            write(time0,2000) nhms0/10000

            if(  diurnal ) hdfile = trim(doutput) // "." // trim(date0) // "." // trim(ext)
            if( mdiurnal ) hdfile = trim(doutput) // "." // trim(date0) // "_" // trim(time0) // "z." // trim(ext)

            if( ldquad ) then
               ndvars = mvars  ! Include Quadratics in Diurnal Files
               if (k==1) then
                  call copy_bundle_to_bundle(final_bundle,diurnal_bundle,_RC)
               end if
            else
               ndvars = nvars2 ! Only Include Primary Fields in Diurnal Files (Default)
               if (k==1) then
                  do n=1,nvars
                     call ESMF_FieldBundleGet(final_bundle,trim(vname2(n)),field=field,_RC)
                     call MAPL_FieldBundleAdd(diurnal_bundle,field,_RC)
                  enddo
               endif
            endif
         endif

! Primary Fields
! --------------
         do n=1,nvars2
            call ESMF_FieldBundleGet(diurnal_bundle,trim(vname2(n)),field=field,_RC)
            if( kmvar2(n).eq.0 ) then
               kbeg = 0
               kend = 1
               call ESMF_FieldGet(field,0,farrayPtr=ptr2d,_RC)
               ptr2d = q(:,:,nloc(n),k)
            else
               kbeg = 1
               kend = kmvar2(n)
               call ESMF_FieldGet(field,0,farrayPtr=ptr3d,_RC)
               ptr3d = q(:,:,nloc(n):nloc(n)+kend-1,k)
            endif
            dum(:,:,1:kend) = q(:,:,nloc(n):nloc(n)+kend-1,k)
         enddo

! Quadratics
! ----------
         if( ndvars.eq.mvars ) then
            mv=  nvars2
            do nv=1,nquad
               if( qloc(1,nv)*qloc(2,nv).ne.0 ) then
                  mv=mv+1
                  call ESMF_FieldBundleGet(diurnal_bundle,trim(vname2(mv)),field=field,_RC)
                  if( kmvar2(qloc(1,nv)).eq.0 ) then
                     kbeg = 0
                     kend = 1
                  else
                     kbeg = 1
                     kend = kmvar2(qloc(1,nv))
                  endif
                  loc1 = nloc( qloc(1,nv) )
                  loc2 = nloc( qloc(2,nv) )
                  if( .not.lzstar(nv) ) then
                     where( q(:,:,nloc(mv):nloc(mv)+kend-1,0).ne.undef )
                        dum(:,:,1:kend) = q(:,:,nloc(mv):nloc(mv)+kend-1,k) - q(:,:,loc1:loc1+kend-1,k) &
                              * q(:,:,loc2:loc2+kend-1,k)
                     elsewhere
                        dum(:,:,1:kend) = undef
                     endwhere
                  else
                     dum(:,:,1:kend) = q(:,:,nloc(mv):nloc(mv)+kend-1,k)
                  endif
                  if( kmvar2(qloc(1,nv)).eq.0 ) then
                     call ESMF_FieldGet(field,0,farrayPtr=ptr2d,_RC)
                     ptr2d = dum(:,:,1)
                  else
                     kend = kmvar2(qloc(1,nv))
                     call ESMF_FieldGet(field,0,farrayPtr=ptr3d,_RC)
                     ptr3d = dum(:,:,1:kend)
                  endif
               endif
            enddo
         endif


         etime = local_esmf_timeset(nymd0,nhms0,_RC)
         call ESMF_ClockSet(clock,currTime=etime, _RC)
         if (k==1 .or. mdiurnal) then
            if (mdiurnal) then
               n_times = 1
            else
               n_times = ntods
            end if
            if (k==1) then
               call diurnal_writer%create_from_bundle(diurnal_bundle,clock,n_steps=n_times,time_interval=timeinc,vertical_data=vertical_data)
            end if
            call diurnal_writer%start_new_file(trim(hdfile),_RC)
         end if
         call diurnal_writer%write_to_file(_RC)
         if( root .and. mdiurnal ) then
            print *, 'Created: ',trim(hdfile)
         endif
         call tick (nymd0,nhms0,ndt)
      enddo

      if( root .and. diurnal ) then
         print *, 'Created: ',trim(hdfile)
      endif
      if( root ) print *

      call t_prof%stop('Write_Diurnal')
   endif

! **********************************************************************
! ****                 Write Min/Max Information                    ****
! **********************************************************************

   if( root ) print *
   do n=1,nvars2
      do L=1,max(1,kmvar2(n))
         if( kmvar2(n).eq.0 ) then
            plev = 0
         else
            plev = lev(L)
         endif

         call mpi_reduce( qmin(nloc(n)+L-1),qming,1,mpi_real,mpi_min,0,comm,ierror )
         _VERIFY(ierror)
         call mpi_reduce( qmax(nloc(n)+L-1),qmaxg,1,mpi_real,mpi_max,0,comm,ierror )
         _VERIFY(ierror)
         if( root ) then
            if(L.eq.1) then
               write(6,3101) trim(vname2(n)),plev,qming,qmaxg
            else
               write(6,3102) trim(vname2(n)),plev,qming,qmaxg
            endif
         endif
3101     format(1x,'Primary Field: ',a20,' Level: ',f9.3,'  Min: ',g15.8,'  Max: ',g15.8)
3102     format(1x,'               ',a20,' Level: ',f9.3,'  Min: ',g15.8,'  Max: ',g15.8)
      enddo
      call MPI_BARRIER(comm,status)
      _VERIFY(status)
      if( root ) print *
   enddo
   if( root ) print *

! **********************************************************************
! ****                     Timing Information                       ****
! **********************************************************************

   call io_server%finalize()
   call t_prof%stop()
   call t_prof%reduce()
   call t_prof%finalize()
   call generate_report()
   call MAPL_Finalize()
   call MPI_Finalize(status)
   stop

contains

   function create_output_grid(grid_name,lm,rc) result(new_grid)
      type(ESMF_Grid) :: new_grid
      character(len=*), intent(inout) :: grid_name
      integer, intent(in) :: lm
      integer, optional, intent(out) :: rc

      type(ESMF_Config) :: cf
      integer :: nn,im_world,jm_world,nx, ny
      character(len=5) :: imsz,jmsz
      character(len=2) :: pole,dateline

      nn   = len_trim(grid_name)
      imsz = grid_name(3:index(grid_name,'x')-1)
      jmsz = grid_name(index(grid_name,'x')+1:nn-3)
      pole = grid_name(1:2)
      dateline = grid_name(nn-1:nn)
      read(IMSZ,*) im_world
      read(JMSZ,*) jm_world

      cf = MAPL_ConfigCreate(_RC)
      call MAPL_ConfigSetAttribute(cf,value=lm, label=trim(grid_name)//".LM:",_RC)
      if (dateline=='CF') then
         call MAPL_MakeDecomposition(nx,ny,reduceFactor=6,_RC)
         call MAPL_ConfigSetAttribute(cf,value=NX, label=trim(grid_name)//".NX:",_RC)
         call MAPL_ConfigSetAttribute(cf,value="Cubed-Sphere", label=trim(grid_name)//".GRID_TYPE:",_RC)
         call MAPL_ConfigSetAttribute(cf,value=6, label=trim(grid_name)//".NF:",_RC)
         call MAPL_ConfigSetAttribute(cf,value=im_world,label=trim(grid_name)//".IM_WORLD:",_RC)
         call MAPL_ConfigSetAttribute(cf,value=ny, label=trim(grid_name)//".NY:",_RC)
      else if (dateline=='TM') then
         _FAIL("Tripolar not yet implemented for outpout")
      else
         call MAPL_MakeDecomposition(nx,ny,_RC)
         call MAPL_ConfigSetAttribute(cf,value=NX, label=trim(grid_name)//".NX:",_RC)
         call MAPL_ConfigSetAttribute(cf,value="LatLon", label=trim(grid_name)//".GRID_TYPE:",_RC)
         call MAPL_ConfigSetAttribute(cf,value=im_world,label=trim(grid_name)//".IM_WORLD:",_RC)
         call MAPL_ConfigSetAttribute(cf,value=jm_world,label=trim(grid_name)//".JM_WORLD:",_RC)
         call MAPL_ConfigSetAttribute(cf,value=ny, label=trim(grid_name)//".NY:",_RC)
         call MAPL_ConfigSetAttribute(cf,value=pole, label=trim(grid_name)//".POLE:",_RC)
         call MAPL_ConfigSetAttribute(cf,value=dateline, label=trim(grid_name)//".DATELINE:",_RC)
         if (pole=='XY' .and. dateline=='XY') then
            _FAIL("regional lat-lon output not supported")
         end if
      end if

      new_grid = grid_manager%make_grid(cf,prefix=trim(grid_name)//".",_RC)
      if (present(rc)) then
         rc=_SUCCESS
      end if
   end function create_output_grid

   subroutine get_file_levels(filename,vertical_data,rc)
      character(len=*), intent(in) :: filename
      type(VerticalData), intent(inout) :: vertical_data
      integer, intent(out), optional :: rc

      integer :: status
      type(NetCDF4_fileFormatter) :: formatter
      type(FileMetadata) :: basic_metadata
      type(FileMetadataUtils) :: metadata
      character(len=:), allocatable :: lev_name
      character(len=ESMF_MAXSTR) :: long_name
      character(len=ESMF_MAXSTR) :: standard_name
      character(len=ESMF_MAXSTR) :: vcoord
      character(len=ESMF_MAXSTR) :: lev_units
      real, allocatable, target :: levs(:)
      real, pointer :: plevs(:)

      call formatter%open(trim(filename),pFIO_Read,_RC)
      basic_metadata=formatter%read(_RC)
      call metadata%create(basic_metadata,trim(filename))
      lev_name = metadata%get_level_name(_RC)
      if (lev_name /= '') then
         call metadata%get_coordinate_info(lev_name,coords=levs,coordUnits=lev_units,long_name=long_name,&
               standard_name=standard_name,coordinate_attr=vcoord,_RC)
         plevs => levs
         vertical_data = VerticalData(levels=plevs,vunit=lev_units,vcoord=vcoord,standard_name=standard_name,long_name=long_name, &
               force_no_regrid=.true.,_RC)
         nullify(plevs)
      end if

      if (present(rc)) then
         rc=_SUCCESS
      end if

   end subroutine get_file_levels

   function has_level(grid,rc) result(grid_has_level)
      logical :: grid_has_level
      type(ESMF_Grid), intent(in) :: grid
      integer, intent(out), optional :: rc
      integer :: status, global_dims(3)
      call MAPL_GridGet(grid,globalCellCountPerDim=global_dims,_RC)
      grid_has_level = (global_dims(3)>0)
      if (present(rc)) then
         RC=_SUCCESS
      end if
   end function has_level

   subroutine copy_bundle_to_bundle(input_bundle,output_bundle,rc)
      type(ESMF_FieldBundle), intent(inout) :: input_bundle
      type(ESMF_FieldBundle), intent(inout) :: output_bundle
      integer, intent(out), optional :: rc
      integer :: status
      character(len=ESMF_MAXSTR), allocatable :: field_list(:)
      type(ESMF_Field) :: field
      integer :: i,num_fields
      call ESMF_FieldBundleGet(input_bundle,fieldCount=num_fields,_RC)
      allocate(field_list(num_fields))
      call ESMF_FieldBundleGet(input_bundle,fieldNameList=field_list,_RC)
      do i=1,num_fields
         call ESMF_FieldBundleGet(input_bundle,field_list(i),field=field,_RC)
         call MAPL_FieldBundleAdd(output_bundle,field,_RC)
      enddo
      if (present(rc)) then
         RC=_SUCCESS
      end if
   end subroutine copy_bundle_to_bundle

   subroutine add_new_field_to_bundle(bundle,grid,lm,field_name,long_name,units,rc)
      type(ESMF_FieldBundle), intent(inout) :: bundle
      type(ESMF_Grid), intent(in) :: grid
      integer, intent(in) :: lm
      character(len=*), intent(in) :: field_name
      character(len=*), intent(in) :: long_name
      character(len=*), intent(in) :: units
      integer, intent(out), optional :: rc

      integer :: status
      type(ESMF_Field) :: field
      type(ESMF_Info) :: infoh

      if (lm == 0) then
         field = ESMF_FieldCreate(grid,name=trim(field_name),typekind=ESMF_TYPEKIND_R4,_RC)
      else if (lm > 0) then
         field = ESMF_FieldCreate(grid,name=trim(field_name),typekind=ESMF_TYPEKIND_R4, &
               ungriddedLBound=[1],ungriddedUBound=[lm],_RC)
      end if
      call ESMF_InfoGetFromHost(field,infoh,_RC)
      call ESMF_InfoSet(infoh,key='LONG_NAME',value=trim(long_name),_RC)
      call ESMF_InfoSet(infoh,key='UNITS',value=trim(units),_RC)
      if (lm == 0) then
         call ESMF_InfoSet(infoh,key='DIMS',value=MAPL_DimsHorzOnly,_RC)
         call ESMF_InfoSet(infoh,key='VLOCATION',value=MAPL_VLocationNone,_RC)
      else if (lm > 0) then
         call ESMF_InfoSet(infoh,key='DIMS',value=MAPL_DimsHorzVert,_RC)
         call ESMF_InfoSet(infoh,key='VLOCATION',value=MAPL_VLocationCenter,_RC)
      end if
      call MAPL_FieldBundleAdd(bundle,field,_RC)
      if (present(rc)) then
         RC=_SUCCESS
      end if
   end subroutine add_new_field_to_bundle

   subroutine get_file_times(file_metadata,num_times,time_series,time_interval,yymmdd,hhmmss,rc)
      type(FileMetadataUtils), intent(inout) :: file_metadata
      integer, intent(out) :: num_times
      type(ESMF_Time), allocatable, intent(inout) :: time_series(:)
      integer, intent(inout), allocatable :: yymmdd(:)
      integer, intent(inout), allocatable :: hhmmss(:)
      integer, intent(out) :: time_interval
      integer, intent(out), optional :: rc

      integer :: status
      type(ESMF_TimeInterval) :: esmf_time_interval
      integer :: hour, minute, second, year, month, day, i

      num_times = file_metadata%get_dimension('time',_RC)
      call file_metadata%get_time_info(timeVector=time_series,_RC)
      if (num_times == 1) then
         time_interval = file_metadata%get_var_attr_int32('time','time_increment',_RC)
      else if (num_times > 1) then
         esmf_time_interval = time_series(2)-time_series(1)
         call ESMF_TimeIntervalGet(esmf_time_interval,h=hour,m=minute,s=second,_RC)
         time_interval = hour*10000+minute*100+second
      end if

      allocate(yymmdd(num_times),hhmmss(num_times))
      do i = 1,num_times
         call ESMF_TimeGet(time_series(i),yy=year,mm=month,dd=day,h=hour,m=minute,s=second,_RC)
         yymmdd(i)=year*10000+month*100+day
         hhmmss(i)=hour*10000+minute*100+second
      enddo
      if (present(rc)) then
         rc=_SUCCESS
      end if
   end subroutine get_file_times

   function get_level_info(bundle,rc) result(kmvar)
      integer, allocatable :: kmvar(:)
      type(ESMF_FieldBundle), intent(in) :: bundle
      integer, optional, intent(out) :: rc

      integer :: status
      character(len=ESMF_MAXSTR), allocatable :: field_list(:)
      type(ESMF_Field) :: field
      integer :: rank,i,num_fields,lb(1),ub(1)
      call ESMF_FieldBundleGet(bundle,fieldCount=num_fields,_RC)
      allocate(field_list(num_fields))
      allocate(kmvar(num_fields))
      call ESMF_FieldBundleGet(bundle,fieldNameList=field_list,_RC)
      do i=1,num_fields
         call ESMF_FieldBundleGet(bundle,field_list(i),field=field,_RC)
         call ESMF_FieldGet(field,rank=rank,_RC)
         if (rank==2) then
            kmvar(i)=0
         else if (rank==3) then
            call ESMF_FieldGet(field,ungriddedLBound=lb,ungriddedUBound=ub,_RC)
            kmvar(i)=ub(1)-lb(1)+1
         else
            _FAIL("Unsupported rank")
         end if
      end do
      if (present(rc)) then
         RC=_SUCCESS
      end if
   end function get_level_info

   function get_long_names(bundle,rc) result(long_names)
      character(len=ESMF_MAXSTR), allocatable :: long_names(:)
      type(ESMF_FieldBundle), intent(in) :: bundle
      integer, optional, intent(out) :: rc

      integer :: status
      character(len=ESMF_MAXSTR), allocatable :: field_list(:)
      type(ESMF_Field) :: field
      type(ESMF_Info) :: infoh
      integer :: i,num_fields

      call ESMF_FieldBundleGet(bundle,fieldCount=num_fields,_RC)
      allocate(field_list(num_fields))
      allocate(long_names(num_fields))
      call ESMF_FieldBundleGet(bundle,fieldNameList=field_list,_RC)
      do i=1,num_fields
         call ESMF_FieldBundleGet(bundle,field_list(i),field=field,_RC)
         call ESMF_InfoGetFromHost(field,infoh,_RC)
         call ESMF_InfoGet(infoh,key='LONG_NAME',value=long_names(i),_RC)
      enddo
      if (present(rc)) then
         RC=_SUCCESS
      end if
   end function get_long_names

   function get_units(bundle,rc) result(units)
      character(len=ESMF_MAXSTR), allocatable :: units(:)
      type(ESMF_FieldBundle), intent(in) :: bundle
      integer, optional, intent(out) :: rc

      integer :: status
      character(len=ESMF_MAXSTR), allocatable :: field_list(:)
      type(ESMF_Field) :: field
      type(ESMF_Info) :: infoh
      integer :: i,num_fields

      call ESMF_FieldBundleGet(bundle,fieldCount=num_fields,_RC)
      allocate(field_list(num_fields))
      allocate(units(num_fields))
      call ESMF_FieldBundleGet(bundle,fieldNameList=field_list,_RC)
      do i=1,num_fields
         call ESMF_FieldBundleGet(bundle,field_list(i),field=field,_RC)
         call ESMF_InfoGetFromHost(field,infoh,_RC)
         call ESMF_InfoGet(infoh,key='UNITS',value=units(i),_RC)
      enddo
      if (present(rc)) then
         RC=_SUCCESS
      end if
   end function get_units

   function local_esmf_timeset(yymmdd,hhmmss,rc) result(etime)
      type(ESMF_Time) :: etime
      integer, intent(in) :: yymmdd
      integer, intent(in) :: hhmmss
      integer, intent(out), optional :: rc

      integer :: year,month,day,hour,minute,second,status
      year = yymmdd/10000
      month = mod(yymmdd/100,100)
      day = mod(yymmdd,100)

      hour = hhmmss/10000
      minute = mod(hhmmss/100,100)
      second = mod(hhmmss,100)

      call ESMF_TimeSet(etime,yy=year,mm=month,dd=day,h=hour,m=minute,s=second,_RC)
      if (present(rc)) then
         rc=_SUCCESS
      endif
   end function local_esmf_timeset

   function defined ( q,undef )
      implicit none
      logical  defined
      real     q,undef
      defined = q /= undef
   end function defined

   subroutine latlon_zstar (q,qp,undef,grid,rc)
      real, intent(inout) :: q(:,:)
      real, intent(out) :: qp(:,:)
      real, intent(in) :: undef
      type (ESMF_Grid), intent(inout) :: grid
      integer, optional, intent(out) :: rc

      integer :: local_dims(3)
      integer im,jm,i,j,status
      real, allocatable :: qz(:)

      call MAPL_GridGet(grid,localCellCountPerDim=local_dims,_RC)
      im = local_dims(1)
      jm = local_dims(2)
      allocate(qz(jm))

      call latlon_zmean ( q,qz,undef,grid )
      do j=1,jm
         if( qz(j).eq. undef ) then
            qp(:,j) = undef
         else
            do i=1,im
               if( defined( q(i,j),undef) ) then
                  qp(i,j) = q(i,j) - qz(j)
               else
                  qp(i,j) = undef
               endif
            enddo
         endif
      enddo
      if (present(rc)) then
         rc=_SUCCESS
      endif
   end subroutine latlon_zstar

   subroutine latlon_zmean ( q,qz,undef,grid,rc)
      real, intent(inout) ::  q(:,:)
      real, intent(inout) ::  qz(:)
      real, intent(in) :: undef
      type(ESMF_Grid), intent(inout) :: grid
      integer, optional, intent(out) :: rc

      integer :: im,jm,im_global,jm_global,local_dims(3),global_dims(3),status,nx,ny
      real, allocatable ::   qg(:,:)
      real, allocatable :: buf(:,:)
      real :: qsum
      integer :: mpistatus(mpi_status_size)
      integer, allocatable :: ims(:),jms(:)
      integer j,n,peid,peid0,i1,j1,in,jn,mypet,i_start,i_end,isum
      type(ESMF_VM) :: vm

      call ESMF_VMGetCurrent(vm,_RC)
      call ESMF_VMGet(vm,localPet=mypet,_RC)
      call MAPL_GridGet(grid,localCellCountPerDim=local_dims,globalCellCountPerDim=global_dims,_RC)
      im = local_dims(1)
      jm = local_dims(2)
      im_global = global_dims(1)
      jm_global = global_dims(2)
      call get_esmf_grid_layout(grid,nx,ny,ims,jms,_RC)
      call mapl_grid_interior(grid,i1,in,j1,jn)

      qz = 0.0
      allocate( qg(im_global,jm) )
      peid0 = (mypet/nx)*ny
      if (i1==1) then
         i_start = 1
         i_end = ims(1)
         qg(i_start:i_end,:)=q
         do n=1,nx-1
            allocate(buf(ims(n+1),jm))
            peid = mypet + n
            call mpi_recv(buf,ims(n+1)*jm,MPI_FLOAT,peid,peid,MPI_COMM_WORLD,mpistatus,status)
            _VERIFY(status)
            i_start=i_end+1
            i_end = i_start+ims(n)-1
            qg(i_start:i_end,:)=buf
            deallocate(buf)
         enddo
      else
         call mpi_send(q,im*jm,MPI_FLOAT,peid0,mypet,MPI_COMM_WORLD,status)
         _VERIFY(status)
      end if

! compute zonal mean
      if (i1 == 1) then
         do j=1,jm
            isum = count(qg(:,j) /= undef)
            qsum = sum(qg(:,j),mask=qg(:,j)/=undef)
            if (isum == 0) then
               qz(j)=undef
            else
               qz(j)=qsum/real(isum)
            end if
         enddo

! send mean back to other ranks
         do n=1,nx-1
            peid = peid0+n
            call mpi_send(qz,jm,MPI_FLOAT,peid,peid0,MPI_COMM_WORLD,status)
            _VERIFY(status)
         enddo
      else
         call mpi_recv(qz,jm,MPI_FLOAT,peid0,peid0,MPI_COMM_WORLD,mpistatus,status)
         _VERIFY(status)
      end if

      if (present(rc)) then
         rc=_SUCCESS
      endif

   end subroutine latlon_zmean

   subroutine get_esmf_grid_layout(grid,nx,ny,ims_out,jms_out,rc)
      type(ESMF_Grid), intent(inout) :: grid
      integer, intent(out) :: nx
      integer, intent(out) :: ny
      integer, intent(inout), allocatable :: ims_out(:)
      integer, intent(inout), allocatable :: jms_out(:)
      integer, optional, intent(out) :: rc

      type(ESMF_VM) :: vm
      integer :: status
      type(ESMF_DistGrid) :: dist_grid
      integer, allocatable :: minindex(:,:),maxindex(:,:)
      integer :: dim_count, ndes
      integer, pointer :: ims(:),jms(:)

      call ESMF_VMGetCurrent(vm,_RC)
      call ESMF_VMGet(vm,petCount=ndes,_RC)
      call ESMF_GridGet(grid,distgrid=dist_grid,dimCOunt=dim_count,_RC)
      allocate(minindex(dim_count,ndes),maxindex(dim_count,ndes))
      call MAPL_DistGridGet(dist_grid,minIndex=minindex,maxIndex=maxindex,_RC)
      call MAPL_GetImsJms(minindex(1,:),maxindex(1,:),minindex(2,:),maxindex(2,:),ims,jms,_RC)
      nx = size(ims)
      ny = size(jms)
      allocate(ims_out(nx),jms_out(ny))
      ims_out = ims
      jms_out = jms

      if (present(rc)) then
         rc=_SUCCESS
      endif

   end subroutine get_esmf_grid_layout

   subroutine check_quad ( quad,vname,nvars,aliases,nalias,qloc )
      integer :: nvars, nalias
      character(len=ESMF_MAXSTR)  quad(2), aliases(2,nalias), vname(nvars)
      integer  qloc(2)
      integer  m,n

! Initialize Location of Quadratics
! ---------------------------------
      qloc = 0

! Check Quadratic Name against HDF Variable Names
! -----------------------------------------------
      do n=1,nvars
         if( trim(vname(n)).eq.trim(quad(1)) ) qloc(1) = n
         if( trim(vname(n)).eq.trim(quad(2)) ) qloc(2) = n
      enddo

! Check Quadratic Name against Aliases
! ------------------------------------
      do m=1,nalias
         if( trim(quad(1)).eq.trim(aliases(1,m)) ) then
            do n=1,nvars
               if( trim(vname(n)).eq.trim(quad(1)) .or. &
                     trim(vname(n)).eq.trim(aliases(2,m)) ) then
                  qloc(1) = n
                  exit
               endif
            enddo
         endif
         if( trim(quad(2)).eq.trim(aliases(1,m)) ) then
            do n=1,nvars
               if( trim(vname(n)).eq.trim(quad(2)) .or. &
                     trim(vname(n)).eq.trim(aliases(2,m)) ) then
                  qloc(2) = n
                  exit
               endif
            enddo
         endif
      enddo

   end subroutine check_quad

   function compute_nsecf (nhms) result(seconds)
      integer :: seconds
      integer, intent(in) :: nhms
      seconds =  nhms/10000*3600 + mod(nhms,10000)/100*60 + mod(nhms,100)
   end function compute_nsecf

   function compute_nhmsf (nsec) result(nhmsf)
      integer :: nhmsf
      integer, intent(in) :: nsec
      nhmsf =  nsec/3600*10000 + mod(nsec,3600)/60*100 + mod(nsec,60)
   end function compute_nhmsf

   subroutine tick (nymd,nhms,ndt)
      integer, intent(inout) :: nymd
      integer, intent(inout) :: nhms
      integer, intent(in) :: ndt

      integer :: nsec

      if(ndt.ne.0) then
         nsec = compute_nsecf(nhms) + ndt

         if (nsec.gt.86400)  then
            do while (nsec.gt.86400)
               nsec = nsec - 86400
               nymd = compute_incymd (nymd,1)
            enddo
         endif

         if (nsec.eq.86400)  then
            nsec = 0
            nymd = compute_incymd (nymd,1)
         endif

         if (nsec.lt.00000)  then
            do while (nsec.lt.0)
               nsec = 86400 + nsec
               nymd = compute_incymd (nymd,-1)
            enddo
         endif

         nhms = compute_nhmsf (nsec)
      endif

   end subroutine tick

   function compute_incymd (nymd,m) result(incymd)
      integer :: incymd
      integer, intent(in) :: nymd
      integer, intent(in) :: m
!***********************************************************************
!  purpose
!     incymd:  nymd changed by one day
!     modymd:  nymd converted to julian date
!  description of parameters
!     nymd     current date in yymmdd format
!     m        +/- 1 (day adjustment)
!
!***********************************************************************
!*                  goddard laboratory for atmospheres                 *
!***********************************************************************

      integer ndpm(12)
      data    ndpm /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/
      integer :: ny,nm,nd
!***********************************************************************
!
      ny = nymd / 10000
      nm = mod(nymd,10000) / 100
      nd = mod(nymd,100) + m

      if (nd.eq.0) then
         nm = nm - 1
         if (nm.eq.0) then
            nm = 12
            ny = ny - 1
         endif
         nd = ndpm(nm)
         if (nm.eq.2 .and. is_leap_year(ny))  nd = 29
      endif

      if (nd.eq.29 .and. nm.eq.2 .and. is_leap_year(ny))  go to 20

      if (nd.gt.ndpm(nm)) then
         nd = 1
         nm = nm + 1
         if (nm.gt.12) then
            nm = 1
            ny = ny + 1
         endif
      endif

20    continue
      incymd = ny*10000 + nm*100 + nd
      return

   end function compute_incymd

   logical function is_leap_year(year)
      integer, intent(in) :: year
      is_leap_year = (mod(year,4) == 0) .and. (mod(year,100) == 0 .or. mod(year,400) == 0)
   end function is_leap_year

   subroutine usage(root)
      logical, intent(in) :: root
      integer :: status,errorcode,rc
      if(root) then
         write(6,100)
100      format(  "usage:  ",/,/ &
               " time_ave.x -hdf      filenames (in hdf format)",/ &
               "           <-template template>"    ,/ &
               "           <-tag      tag>"    ,/ &
               "           <-rc       rcfile>" ,/ &
               "           <-ntod     ntod>"   ,/ &
               "           <-ntmin    ntmin>"  ,/ &
               "           <-strict   strict>" ,/ &
               "           <-d>"               ,/ &
               "           <-md>"              ,/,/ &
               "where:",/,/ &
               "  -hdf      filenames:  filenames (in hdf format) to average",/ &
               "  -template  template:  filename to use as template if hdf files differ (default: 1st filename)",/ &
               "  -begdate   yyyymmdd:  optional parameter for date to begin averaging",/ &
               "  -begtime     hhmmss:  optional parameter for time to begin averaging",/ &
               "  -enddate   yyyymmdd:  optional parameter for date to   end averaging",/ &
               "  -endtime     hhmmss:  optional parameter for time to   end averaging",/ &
               "  -tag            tag:  optional tag for output file (default: monthly_ave)",/ &
               "  -rc          rcfile:  optional resource filename for quadratics (default: no quadratics)",/ &
               "  -ntod          ntod:  optional time-of-day (hhmmss) to average (default: all time periods)",/ &
               "  -ntmin        ntmin:  optional parameter for required min. timeperiods (default: 10 days equiv)",/ &
               "  -strict      strict:  optional logical parameter for strict time testing (default: .true.)",/ &
               "  -d             dtag:  optional parameter to create & tag          monthly mean diurnal file  ", &
               "(all times included)",/ &
               "  -md            dtag:  optional parameter to create & tag multiple monthly mean diurnal files ", &
               "(one time  per file)",/ &
               "  -dv            dtag:  like -d  but includes diurnal variances",/ &
               "  -mdv           dtag:  like -md but includes diurnal variances",/ &
               )
      endif
      call MPI_Abort(MPI_COMM_WORLD,errorcode,status)
      _VERIFY(status)
   end subroutine usage

    subroutine generate_report()

         character(:), allocatable :: report_lines(:)
         integer :: i
         character(1) :: empty(0)

         reporter = ProfileReporter(empty)
         call reporter%add_column(NameColumn(20))
         call reporter%add_column(FormattedTextColumn('Inclusive','(f12.2)', 12, InclusiveColumn('MEAN')))
         call reporter%add_column(FormattedTextColumn('% Incl','(f6.2)', 6, PercentageColumn(InclusiveColumn('MEAN'),'MAX')))
         call reporter%add_column(FormattedTextColumn('Exclusive','(f12.2)', 12, ExclusiveColumn('MEAN')))
         call reporter%add_column(FormattedTextColumn('% Excl','(f6.2)', 6, PercentageColumn(ExclusiveColumn('MEAN'))))
         call reporter%add_column(FormattedTextColumn(' Max Excl)','(f12.2)', 12, ExclusiveColumn('MAX')))
         call reporter%add_column(FormattedTextColumn(' Min Excl)','(f12.2)', 12, ExclusiveColumn('MIN')))
         call reporter%add_column(FormattedTextColumn('Max PE)','(1x,i5.5,1x)', 7, ExclusiveColumn('MAX_PE')))
         call reporter%add_column(FormattedTextColumn('Min PE)','(1x,i5.5,1x)', 7, ExclusiveColumn('MIN_PE')))
        report_lines = reporter%generate_report(t_prof)
         if (mapl_am_I_root()) then
            write(*,'(a)')'Final profile'
            write(*,'(a)')'============='
            do i = 1, size(report_lines)
               write(*,'(a)') report_lines(i)
            end do
            write(*,'(a)') ''
         end if
    end subroutine generate_report


end program time_ave