#define I_AM_MAIN #include "MAPL_Generic.h" 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