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
! **********************************************************************
! **** 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_AttributeGet(output_grid,'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
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_AttributeSet(field,name='LONG_NAME',value=trim(long_name),_RC)
call ESMF_AttributeSet(field,name='UNITS',value=trim(units),_RC)
if (lm == 0) then
call ESMF_AttributeSet(field,name='DIMS',value=MAPL_DimsHorzOnly,_RC)
call ESMF_AttributeSet(field,name='VLOCATION',value=MAPL_VLocationNone,_RC)
else if (lm > 0) then
call ESMF_AttributeSet(field,name='DIMS',value=MAPL_DimsHorzVert,_RC)
call ESMF_AttributeSet(field,name='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
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_AttributeGet(field,name='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
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_AttributeGet(field,name='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