!----------------------------------------------------------------------- ! module module_wrt_grid_comp ! !----------------------------------------------------------------------- !*** This module includes the functionality of write gridded component. !----------------------------------------------------------------------- !*** At initialization step, write grid is defined. The forecast field !*** bundle is mirrored and output field information inside the field !*** bundle is used to create ESMF field on the write grid and added in !*** the output field bundle on write grid component. Also the IO_BaseTime !*** is set to the initial clock time. !*** At the run step, output time is set from the write grid comp clock !*** the ESMF field bundles that contains the data on write grid are !*** written out through ESMF field bundle write to netcdf files. !*** The ESMF field bundle write uses parallel write, so if output grid !*** is cubed sphere grid, the six tiles file will be written out at !*** same time. !----------------------------------------------------------------------- !*** !*** Revision history !*** ! Jul 2017: J. Wang/G. Theurich - initial code for fv3 write grid component ! Mar 2018: S Moorthi - changing cfhour to accommodate up to 99999 hours ! Aug 2019: J. Wang - add inline post ! !--------------------------------------------------------------------------------- ! use mpi use esmf use write_internal_state use module_fv3_io_def, only : num_pes_fcst, & n_group, num_files, & filename_base, output_grid, output_file, & imo,jmo,ichunk2d,jchunk2d, & ichunk3d,jchunk3d,kchunk3d,nbits, & nsout => nsout_io, & cen_lon, cen_lat, & lon1, lat1, lon2, lat2, dlon, dlat, & stdlat1, stdlat2, dx, dy, iau_offset, & ideflate, lflname_fulltime use module_write_netcdf, only : write_netcdf use physcons, only : pi => con_pi #ifdef INLINE_POST use post_fv3, only : post_run_fv3 #endif ! !----------------------------------------------------------------------- ! implicit none ! !----------------------------------------------------------------------- private ! !----------------------------------------------------------------------- ! ! integer,save :: lead_write_task !<-- Rank of the first write task in the write group integer,save :: last_write_task !<-- Rank of the last write task in the write group integer,save :: ntasks !<-- # of write tasks in the current group integer,save :: itasks, jtasks !<-- # of write tasks in i/j direction in the current group integer,save :: ngrids integer,save :: wrt_mpi_comm !<-- the mpi communicator in the write comp integer,save :: idate(7) logical,save :: write_nsflip logical,save :: change_wrtidate=.false. ! !----------------------------------------------------------------------- ! type(ESMF_FieldBundle) :: gridFB integer :: FBcount character(len=esmf_maxstr),allocatable :: fcstItemNameList(:) ! !----------------------------------------------------------------------- REAL(KIND=8) :: btim,btim0 REAL(KIND=8),PUBLIC,SAVE :: write_init_tim, write_run_tim REAL(KIND=8), parameter :: radi=180.0d0/pi !----------------------------------------------------------------------- ! public SetServices ! interface splat module procedure splat4 module procedure splat8 end interface splat ! type optimizeT type(ESMF_State) :: state type(ESMF_GridComp), allocatable :: comps(:) end type contains ! !----------------------------------------------------------------------- !####################################################################### !----------------------------------------------------------------------- ! subroutine SetServices(wrt_comp, rc) type(ESMF_GridComp) :: wrt_comp integer, intent(out) :: rc rc = ESMF_SUCCESS call ESMF_GridCompSetEntryPoint(wrt_comp, ESMF_METHOD_INITIALIZE, phase=1, & userRoutine=wrt_initialize_p1, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_GridCompSetEntryPoint(wrt_comp, ESMF_METHOD_INITIALIZE, phase=2, & userRoutine=wrt_initialize_p2, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_GridCompSetEntryPoint(wrt_comp, ESMF_METHOD_INITIALIZE, phase=3, & userRoutine=wrt_initialize_p3, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_GridCompSetEntryPoint(wrt_comp, ESMF_METHOD_RUN, & userRoutine=wrt_run, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_GridCompSetEntryPoint(wrt_comp, ESMF_METHOD_FINALIZE, & userRoutine=wrt_finalize, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return end subroutine SetServices ! !----------------------------------------------------------------------- !####################################################################### !----------------------------------------------------------------------- ! subroutine wrt_initialize_p1(wrt_comp, imp_state_write, exp_state_write, clock, rc) ! !----------------------------------------------------------------------- !*** INITIALIZE THE WRITE GRIDDED COMPONENT. !----------------------------------------------------------------------- ! type(esmf_GridComp) :: wrt_comp type(ESMF_State) :: imp_state_write, exp_state_write type(esmf_Clock) :: clock integer,intent(out) :: rc ! !*** LOCAL VARIABLES ! TYPE(ESMF_VM) :: VM type(write_wrap) :: WRAP type(wrt_internal_state),pointer :: wrt_int_state integer :: tl, i, j, n, k integer,dimension(2,6) :: decomptile integer,dimension(2) :: regDecomp !define delayout for the nest grid integer :: fieldCount integer :: vm_mpi_comm character(40) :: fieldName type(ESMF_Config) :: cf, cf_output_grid type(ESMF_Info) :: info type(ESMF_DELayout) :: delayout type(ESMF_Grid) :: fcstGrid type(ESMF_Grid), allocatable :: wrtGrid(:) type(ESMF_Array) :: array type(ESMF_Field) :: field_work, field type(ESMF_Decomp_Flag) :: decompflagPTile(2,6) type(ESMF_StateItem_Flag), allocatable :: fcstItemTypeList(:) type(ESMF_FieldBundle) :: fcstFB, fieldbundle, mirrorFB type(ESMF_Field), allocatable :: fcstField(:) type(ESMF_TypeKind_Flag) :: typekind character(len=80), allocatable :: fieldnamelist(:) integer :: fieldDimCount, gridDimCount, tk integer, allocatable :: petMap(:) integer, allocatable :: gridToFieldMap(:) integer, allocatable :: ungriddedLBound(:) integer, allocatable :: ungriddedUBound(:) character(len=80) :: attName character(len=80), allocatable :: attNameList(:),attNameList2(:) type(ESMF_TypeKind_Flag), allocatable :: typekindList(:) character(len=80) :: valueS integer :: valueI4 real(ESMF_KIND_R4) :: valueR4 real(ESMF_KIND_R8) :: valueR8 logical, allocatable :: is_moving(:) integer :: attCount, jidx, idx, noutfile character(19) :: newdate character(128) :: FBlist_outfilename(100), outfile_name character(128),dimension(:,:), allocatable :: outfilename real(8), dimension(:), allocatable :: slat real(8), dimension(:), allocatable :: lat, lon real(ESMF_KIND_R8), dimension(:,:), pointer :: lonPtr, latPtr real(ESMF_KIND_R8) :: rot_lon, rot_lat real(ESMF_KIND_R8) :: geo_lon, geo_lat real(ESMF_KIND_R8) :: lon1_r8, lat1_r8 real(ESMF_KIND_R8) :: x1, y1, x, y, delat, delon type(ESMF_TimeInterval) :: IAU_offsetTI character(256) :: cf_open, cf_close character(256) :: gridfile integer :: num_output_file ! logical :: lprnt integer :: grid_id logical :: top_parent_is_global ! !----------------------------------------------------------------------- !*********************************************************************** !----------------------------------------------------------------------- ! rc = ESMF_SUCCESS ! !----------------------------------------------------------------------- !*** initialize the write component timers. !----------------------------------------------------------------------- ! write_init_tim = 0. write_run_tim = 0. btim0 = MPI_Wtime() ! !----------------------------------------------------------------------- !*** set the write component's internal state. !----------------------------------------------------------------------- ! allocate(wrt_int_state,stat=RC) wrap%write_int_state => wrt_int_state call ESMF_GridCompSetInternalState(wrt_comp, wrap, rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! call ESMF_VMGetCurrent(vm=VM,rc=RC) call ESMF_VMGet(vm=VM, localPet=wrt_int_state%mype, & petCount=wrt_int_state%petcount,mpiCommunicator=vm_mpi_comm,rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call mpi_comm_dup(vm_mpi_comm, wrt_mpi_comm, rc) ntasks = wrt_int_state%petcount jidx = wrt_int_state%petcount/6 lead_write_task = 0 last_write_task = ntasks -1 lprnt = lead_write_task == wrt_int_state%mype ! print *,'in wrt, lead_write_task=', & ! lead_write_task,'last_write_task=',last_write_task, & ! 'mype=',wrt_int_state%mype,'jidx=',jidx,' comm=',wrt_mpi_comm ! !----------------------------------------------------------------------- !*** get configuration variables !----------------------------------------------------------------------- ! call ESMF_GridCompGet(gridcomp=wrt_comp,config=CF,rc=RC) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__)) return ! bail out ! variables for post call ESMF_ConfigGetAttribute(config=CF,value=wrt_int_state%output_history,default=.true., & label='output_history:',rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__)) return call ESMF_ConfigGetAttribute(config=CF,value=wrt_int_state%write_dopost,default=.false., & label='write_dopost:',rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__)) return call ESMF_ConfigGetAttribute(config=CF,value=write_nsflip,default=.false., & label='write_nsflip:',rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__)) return if( wrt_int_state%write_dopost ) then #ifdef INLINE_POST call ESMF_ConfigGetAttribute(config=CF,value=wrt_int_state%post_namelist,default='itag', & label ='post_namelist:',rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__)) return #else rc = ESMF_RC_NOT_IMPL print *,'inline post not available on this machine' if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__)) return #endif endif allocate(output_file(num_files)) num_output_file = ESMF_ConfigGetLen(config=CF, label ='output_file:',rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if (num_files == num_output_file) then call ESMF_ConfigGetAttribute(CF,valueList=output_file,label='output_file:', & count=num_files, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return do i = 1, num_files if(output_file(i) /= "netcdf" .and. output_file(i) /= "netcdf_parallel") then write(0,*)"Only netcdf and netcdf_parallel are allowed for multiple values of output_file" call ESMF_Finalize(endflag=ESMF_END_ABORT) endif enddo else if ( num_output_file == 1) then call ESMF_ConfigGetAttribute(CF,valuelist=output_file,label='output_file:', count=1, rc=rc) output_file(1:num_files) = output_file(1) else output_file(1:num_files) = 'netcdf' endif if(lprnt) then print *,'num_files=',num_files do i=1,num_files print *,'num_file=',i,'filename_base= ',trim(filename_base(i)),' output_file= ',trim(output_file(i)) enddo endif call ESMF_InfoGetFromHost(imp_state_write, info=info, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_InfoGetAlloc(info, key="is_moving", values=is_moving, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeGet(imp_state_write, convention="NetCDF", purpose="FV3", & name="ngrids", value=ngrids, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeGet(imp_state_write, convention="NetCDF", purpose="FV3", & name="top_parent_is_global", value=top_parent_is_global, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return allocate(wrtGrid(ngrids)) allocate(output_grid(ngrids)) allocate(imo(ngrids)) allocate(jmo(ngrids)) allocate(cen_lon(ngrids)) allocate(cen_lat(ngrids)) allocate(lon1(ngrids)) allocate(lat1(ngrids)) allocate(lon2(ngrids)) allocate(lat2(ngrids)) allocate(dlon(ngrids)) allocate(dlat(ngrids)) allocate(stdlat1(ngrids)) allocate(stdlat2(ngrids)) allocate(dx(ngrids)) allocate(dy(ngrids)) allocate(ichunk2d(ngrids)) allocate(jchunk2d(ngrids)) allocate(ichunk3d(ngrids)) allocate(jchunk3d(ngrids)) allocate(kchunk3d(ngrids)) allocate(ideflate(ngrids)) allocate(nbits(ngrids)) allocate(wrt_int_state%out_grid_info(ngrids)) do n=1, ngrids if (n == 1) then ! for top level domain look directly in cf cf_output_grid = cf else ! for nest domains, look under specific section write(cf_open,'("")') n write(cf_close,'("")') n cf_output_grid = ESMF_ConfigCreate(cf, openLabel=trim(cf_open), closeLabel=trim(cf_close), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return end if call ESMF_ConfigGetAttribute(config=cf_output_grid, value=output_grid(n), label ='output_grid:',rc=rc) if (lprnt) then print *,'grid_id= ', n, ' output_grid= ', trim(output_grid(n)) end if call ESMF_ConfigGetAttribute(config=CF, value=itasks,default=1,label ='itasks:',rc=rc) jtasks = ntasks if(itasks > 0 ) jtasks = ntasks/itasks if( itasks*jtasks /= ntasks ) then itasks = 1 jtasks = ntasks endif if (trim(output_grid(n)) == 'gaussian_grid' .or. trim(output_grid(n)) == 'global_latlon') then call ESMF_ConfigGetAttribute(config=cf_output_grid, value=imo(n), label ='imo:',rc=rc) call ESMF_ConfigGetAttribute(config=cf_output_grid, value=jmo(n), label ='jmo:',rc=rc) if (lprnt) then print *,'imo=',imo(n),'jmo=',jmo(n) end if else if (trim(output_grid(n)) == 'regional_latlon') then call ESMF_ConfigGetAttribute(config=cf_output_grid, value=lon1(n), label ='lon1:',rc=rc) call ESMF_ConfigGetAttribute(config=cf_output_grid, value=lat1(n), label ='lat1:',rc=rc) call ESMF_ConfigGetAttribute(config=cf_output_grid, value=lon2(n), label ='lon2:',rc=rc) call ESMF_ConfigGetAttribute(config=cf_output_grid, value=lat2(n), label ='lat2:',rc=rc) call ESMF_ConfigGetAttribute(config=cf_output_grid, value=dlon(n), label ='dlon:',rc=rc) call ESMF_ConfigGetAttribute(config=cf_output_grid, value=dlat(n), label ='dlat:',rc=rc) imo(n) = (lon2(n)-lon1(n))/dlon(n) + 1 jmo(n) = (lat2(n)-lat1(n))/dlat(n) + 1 if (lprnt) then print *,'lon1=',lon1(n),' lat1=',lat1(n) print *,'lon2=',lon2(n),' lat2=',lat2(n) print *,'dlon=',dlon(n),' dlat=',dlat(n) print *,'imo =',imo(n), ' jmo =',jmo(n) end if else if (trim(output_grid(n)) == 'rotated_latlon') then call ESMF_ConfigGetAttribute(config=cf_output_grid, value=cen_lon(n), label ='cen_lon:',rc=rc) call ESMF_ConfigGetAttribute(config=cf_output_grid, value=cen_lat(n), label ='cen_lat:',rc=rc) call ESMF_ConfigGetAttribute(config=cf_output_grid, value=lon1(n), label ='lon1:', rc=rc) call ESMF_ConfigGetAttribute(config=cf_output_grid, value=lat1(n), label ='lat1:', rc=rc) call ESMF_ConfigGetAttribute(config=cf_output_grid, value=lon2(n), label ='lon2:', rc=rc) call ESMF_ConfigGetAttribute(config=cf_output_grid, value=lat2(n), label ='lat2:', rc=rc) call ESMF_ConfigGetAttribute(config=cf_output_grid, value=dlon(n), label ='dlon:', rc=rc) call ESMF_ConfigGetAttribute(config=cf_output_grid, value=dlat(n), label ='dlat:', rc=rc) imo(n) = (lon2(n)-lon1(n))/dlon(n) + 1 jmo(n) = (lat2(n)-lat1(n))/dlat(n) + 1 if (lprnt) then print *,'cen_lon=',cen_lon(n),' cen_lat=',cen_lat(n) print *,'lon1 =',lon1(n), ' lat1 =',lat1(n) print *,'lon2 =',lon2(n), ' lat2 =',lat2(n) print *,'dlon =',dlon(n), ' dlat =',dlat(n) print *,'imo =',imo(n), ' jmo =',jmo(n) end if else if (trim(output_grid(n)) == 'rotated_latlon_moving' .or. & trim(output_grid(n)) == 'regional_latlon_moving') then call ESMF_ConfigGetAttribute(config=cf_output_grid, value=imo(n), label ='imo:', rc=rc) call ESMF_ConfigGetAttribute(config=cf_output_grid, value=jmo(n), label ='jmo:', rc=rc) call ESMF_ConfigGetAttribute(config=cf_output_grid, value=dlon(n), label ='dlon:',rc=rc) call ESMF_ConfigGetAttribute(config=cf_output_grid, value=dlat(n), label ='dlat:',rc=rc) if (lprnt) then print *,'imo =',imo(n), ' jmo =',jmo(n) print *,'dlon=',dlon(n),' dlat=',dlat(n) end if else if (trim(output_grid(n)) == 'lambert_conformal') then call ESMF_ConfigGetAttribute(config=cf_output_grid, value=cen_lon(n), label ='cen_lon:',rc=rc) call ESMF_ConfigGetAttribute(config=cf_output_grid, value=cen_lat(n), label ='cen_lat:',rc=rc) call ESMF_ConfigGetAttribute(config=cf_output_grid, value=stdlat1(n), label ='stdlat1:',rc=rc) call ESMF_ConfigGetAttribute(config=cf_output_grid, value=stdlat2(n), label ='stdlat2:',rc=rc) call ESMF_ConfigGetAttribute(config=cf_output_grid, value=imo(n), label ='nx:', rc=rc) call ESMF_ConfigGetAttribute(config=cf_output_grid, value=jmo(n), label ='ny:', rc=rc) call ESMF_ConfigGetAttribute(config=cf_output_grid, value=lon1(n), label ='lon1:', rc=rc) call ESMF_ConfigGetAttribute(config=cf_output_grid, value=lat1(n), label ='lat1:', rc=rc) call ESMF_ConfigGetAttribute(config=cf_output_grid, value=dx(n), label ='dx:', rc=rc) call ESMF_ConfigGetAttribute(config=cf_output_grid, value=dy(n), label ='dy:', rc=rc) if (lprnt) then print *,'cen_lon=',cen_lon(n),' cen_lat=',cen_lat(n) print *,'stdlat1=',stdlat1(n),' stdlat2=',stdlat2(n) print *,'lon1=',lon1(n),' lat1=',lat1(n) print *,'nx=',imo(n), ' ny=',jmo(n) print *,'dx=',dx(n),' dy=',dy(n) endif endif ! output_grid ! chunksizes for netcdf_parallel call ESMF_ConfigGetAttribute(config=CF,value=ichunk2d(n),default=0,label ='ichunk2d:',rc=rc) call ESMF_ConfigGetAttribute(config=CF,value=jchunk2d(n),default=0,label ='jchunk2d:',rc=rc) call ESMF_ConfigGetAttribute(config=CF,value=ichunk3d(n),default=0,label ='ichunk3d:',rc=rc) call ESMF_ConfigGetAttribute(config=CF,value=jchunk3d(n),default=0,label ='jchunk3d:',rc=rc) call ESMF_ConfigGetAttribute(config=CF,value=kchunk3d(n),default=0,label ='kchunk3d:',rc=rc) ! zlib compression flag call ESMF_ConfigGetAttribute(config=CF,value=ideflate(n),default=0,label ='ideflate:',rc=rc) if (ideflate(n) < 0) ideflate(n)=0 call ESMF_ConfigGetAttribute(config=CF,value=nbits(n),default=0,label ='nbits:',rc=rc) if (lprnt) then print *,'ideflate=',ideflate(n),' nbits=',nbits(n) end if ! nbits quantization level for lossy compression (must be between 1 and 31) ! 1 is most compression, 31 is least. If outside this range, set to zero ! which means use lossless compression. if (nbits(n) < 1 .or. nbits(n) > 31) nbits(n)=0 ! lossless compression (no quantization) if (cf_output_grid /= cf) then ! destroy the temporary config object created for nest domains call ESMF_ConfigDestroy(config=cf_output_grid, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return endif if ( trim(output_grid(n)) == 'cubed_sphere_grid' ) then !*** Create cubed sphere grid from file if (top_parent_is_global .and. n==1) then gridfile = 'grid_spec.nc' ! global top-level parent do tl=1,6 decomptile(1,tl) = 1 decomptile(2,tl) = jidx decompflagPTile(:,tl) = (/ESMF_DECOMP_SYMMEDGEMAX,ESMF_DECOMP_SYMMEDGEMAX/) enddo call ESMF_AttributeGet(imp_state_write, convention="NetCDF", purpose="FV3", & name="gridfile", value=gridfile, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_LogWrite("wrtComp: gridfile:"//trim(gridfile),ESMF_LOGMSG_INFO,rc=rc) wrtGrid(n) = ESMF_GridCreateMosaic(filename="INPUT/"//trim(gridfile), & regDecompPTile=decomptile,tileFilePath="INPUT/", & decompflagPTile=decompflagPTile, & staggerlocList=(/ESMF_STAGGERLOC_CENTER, ESMF_STAGGERLOC_CORNER/), & name='wrt_grid', rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return else if (top_parent_is_global) then write(gridfile,'(A,I2.2,A,I1,A)') 'grid.nest', n, '.tile', n+5, '.nc' else if (n == 1) then gridfile='grid.tile7.halo0.nc' ! regional top-level parent else write(gridfile,'(A,I2.2,A,I1,A)') 'grid.nest', n, '.tile', n, '.nc' endif end if regDecomp(1) = 1 regDecomp(2) = ntasks allocate(petMap(ntasks)) do i=1, ntasks petMap(i) = i-1 enddo delayout = ESMF_DELayoutCreate(petMap=petMap, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! create the nest Grid by reading it from file but use DELayout call ESMF_LogWrite("wrtComp: gridfile:"//trim(gridfile),ESMF_LOGMSG_INFO,rc=rc) wrtGrid(n) = ESMF_GridCreate(filename="INPUT/"//trim(gridfile), & fileformat=ESMF_FILEFORMAT_GRIDSPEC, regDecomp=regDecomp, & decompflag=(/ESMF_DECOMP_SYMMEDGEMAX,ESMF_DECOMP_SYMMEDGEMAX/), & delayout=delayout, isSphere=.false., indexflag=ESMF_INDEX_DELOCAL, & rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if (lprnt) print *,'in nested/regional cubed_sphere grid, regDecomp=',regDecomp,' PetMap=',petMap(1),petMap(ntasks), & 'gridfile=',trim(gridfile) deallocate(petMap) endif else ! non 'cubed_sphere_grid' if ( trim(output_grid(n)) == 'gaussian_grid') then wrtGrid(n) = ESMF_GridCreate1PeriDim(minIndex=(/1,1/), & maxIndex=(/imo(n),jmo(n)/), regDecomp=(/itasks,jtasks/), & indexflag=ESMF_INDEX_GLOBAL, & name='wrt_grid',rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_GridAddCoord(wrtGrid(n), staggerLoc=ESMF_STAGGERLOC_CENTER, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_GridGetCoord(wrtGrid(n), coordDim=1, farrayPtr=lonPtr, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_GridGetCoord(wrtGrid(n), coordDim=2, farrayPtr=latPtr, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return allocate(slat(jmo(n)), lat(jmo(n)), lon(imo(n))) call splat(4, jmo(n), slat) if(write_nsflip) then do j=1,jmo(n) lat(j) = asin(slat(j)) * radi enddo else do j=1,jmo(n) lat(jmo(n)-j+1) = asin(slat(j)) * radi enddo endif do j=1,imo(n) lon(j) = 360.d0/real(imo(n),8) *real(j-1,8) enddo do j=lbound(latPtr,2),ubound(latPtr,2) do i=lbound(lonPtr,1),ubound(lonPtr,1) lonPtr(i,j) = 360.d0/real(imo(n),8) * real(i-1,8) latPtr(i,j) = lat(j) enddo enddo lon1(n) = lon(1) lon2(n) = lon(imo(n)) lat1(n) = lat(1) lat2(n) = lat(jmo(n)) dlon(n) = 360.d0/real(imo(n),8) dlat(n) = 180.d0/real(jmo(n),8) deallocate(slat, lat, lon) else if ( trim(output_grid(n)) == 'global_latlon') then wrtGrid(n) = ESMF_GridCreate1PeriDim(minIndex=(/1,1/), & maxIndex=(/imo(n),jmo(n)/), regDecomp=(/itasks,jtasks/), & indexflag=ESMF_INDEX_GLOBAL, & name='wrt_grid',rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_GridAddCoord(wrtGrid(n), staggerLoc=ESMF_STAGGERLOC_CENTER, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_GridGetCoord(wrtGrid(n), coordDim=1, farrayPtr=lonPtr, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_GridGetCoord(wrtGrid(n), coordDim=2, farrayPtr=latPtr, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return allocate(lat(jmo(n)), lon(imo(n))) if (mod(jmo(n),2) == 0) then ! if jmo even, lats do not include poles and equator delat = 180.d0/real(jmo(n),8) if(write_nsflip) then do j=1,jmo(n) lat(j) = 90.d0 - 0.5*delat - real(j-1,8)*delat enddo else do j=1,jmo(n) lat(j) = -90.d0 + 0.5*delat + real(j-1,8)*delat enddo endif else ! if jmo odd, lats include poles and equator delat = 180.d0/real(jmo(n)-1,8) if(write_nsflip) then do j=1,jmo(n) lat(j) = 90.d0 - real(j-1,8)*delat enddo else do j=1,jmo(n) lat(j) = -90.d0 + real(j-1,8)*delat enddo endif endif delon = 360.d0/real(imo(n),8) do i=1,imo(n) lon(i) = real(i-1,8)*delon enddo do j=lbound(latPtr,2),ubound(latPtr,2) do i=lbound(lonPtr,1),ubound(lonPtr,1) lonPtr(i,j) = lon(i) latPtr(i,j) = lat(j) enddo enddo lon1(n) = lon(1) lon2(n) = lon(imo(n)) lat1(n) = lat(1) lat2(n) = lat(jmo(n)) dlon(n) = delon dlat(n) = delat deallocate(lat, lon) else if ( trim(output_grid(n)) == 'regional_latlon' .or. & trim(output_grid(n)) == 'regional_latlon_moving' .or. & trim(output_grid(n)) == 'rotated_latlon' .or. & trim(output_grid(n)) == 'rotated_latlon_moving' .or. & trim(output_grid(n)) == 'lambert_conformal' ) then wrtGrid(n) = ESMF_GridCreateNoPeriDim(minIndex=(/1,1/), & maxIndex=(/imo(n),jmo(n)/), regDecomp=(/itasks,jtasks/), & indexflag=ESMF_INDEX_GLOBAL, & name='wrt_grid',rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_GridAddCoord(wrtGrid(n), staggerLoc=ESMF_STAGGERLOC_CENTER, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_GridGetCoord(wrtGrid(n), coordDim=1, farrayPtr=lonPtr, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_GridGetCoord(wrtGrid(n), coordDim=2, farrayPtr=latPtr, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if ( trim(output_grid(n)) == 'regional_latlon' ) then do j=lbound(lonPtr,2),ubound(lonPtr,2) do i=lbound(lonPtr,1),ubound(lonPtr,1) lonPtr(i,j) = lon1(n) + (lon2(n)-lon1(n))/(imo(n)-1) * (i-1) latPtr(i,j) = lat1(n) + (lat2(n)-lat1(n))/(jmo(n)-1) * (j-1) enddo enddo else if ( trim(output_grid(n)) == 'regional_latlon_moving' ) then ! Do not compute lonPtr, latPtr here. Will be done in the run phase else if ( trim(output_grid(n)) == 'rotated_latlon' ) then do j=lbound(lonPtr,2),ubound(lonPtr,2) do i=lbound(lonPtr,1),ubound(lonPtr,1) rot_lon = lon1(n) + (lon2(n)-lon1(n))/(imo(n)-1) * (i-1) rot_lat = lat1(n) + (lat2(n)-lat1(n))/(jmo(n)-1) * (j-1) call rtll(rot_lon, rot_lat, geo_lon, geo_lat, dble(cen_lon(n)), dble(cen_lat(n))) if (geo_lon < 0.0) geo_lon = geo_lon + 360.0 lonPtr(i,j) = geo_lon latPtr(i,j) = geo_lat enddo enddo rot_lon = lon1(n) rot_lat = lat1(n) call rtll(rot_lon, rot_lat, geo_lon, geo_lat, dble(cen_lon(n)), dble(cen_lat(n))) if (geo_lon < 0.0) geo_lon = geo_lon + 360.0 wrt_int_state%out_grid_info(n)%lonstart = geo_lon wrt_int_state%out_grid_info(n)%latstart = geo_lat rot_lon = lon2(n) rot_lat = lat1(n) call rtll(rot_lon, rot_lat, geo_lon, geo_lat, dble(cen_lon(n)), dble(cen_lat(n))) if (geo_lon < 0.0) geo_lon = geo_lon + 360.0 wrt_int_state%out_grid_info(n)%lonse = geo_lon wrt_int_state%out_grid_info(n)%latse = geo_lat rot_lon = lon1(n) rot_lat = lat2(n) call rtll(rot_lon, rot_lat, geo_lon, geo_lat, dble(cen_lon(n)), dble(cen_lat(n))) if (geo_lon < 0.0) geo_lon = geo_lon + 360.0 wrt_int_state%out_grid_info(n)%lonnw = geo_lon wrt_int_state%out_grid_info(n)%latnw = geo_lat rot_lon = lon2(n) rot_lat = lat2(n) call rtll(rot_lon, rot_lat, geo_lon, geo_lat, dble(cen_lon(n)), dble(cen_lat(n))) if (geo_lon < 0.0) geo_lon = geo_lon + 360.0 wrt_int_state%out_grid_info(n)%lonlast = geo_lon wrt_int_state%out_grid_info(n)%latlast = geo_lat else if ( trim(output_grid(n)) == 'rotated_latlon_moving' ) then ! Do not compute lonPtr, latPtr here. Will be done in the run phase else if ( trim(output_grid(n)) == 'lambert_conformal' ) then lon1_r8 = dble(lon1(n)) lat1_r8 = dble(lat1(n)) call lambert(dble(stdlat1(n)),dble(stdlat2(n)),dble(cen_lat(n)),dble(cen_lon(n)), & lon1_r8,lat1_r8,x1,y1, 1) do j=lbound(lonPtr,2),ubound(lonPtr,2) do i=lbound(lonPtr,1),ubound(lonPtr,1) x = x1 + dx(n) * (i-1) y = y1 + dy(n) * (j-1) call lambert(dble(stdlat1(n)),dble(stdlat2(n)),dble(cen_lat(n)),dble(cen_lon(n)), & geo_lon,geo_lat,x,y,-1) if (geo_lon <0.0) geo_lon = geo_lon + 360.0 lonPtr(i,j) = geo_lon latPtr(i,j) = geo_lat enddo enddo endif else write(0,*)"wrt_initialize_p1: Unknown output_grid ", trim(output_grid(n)) call ESMF_LogWrite("wrt_initialize_p1: Unknown output_grid "//trim(output_grid(n)),ESMF_LOGMSG_ERROR,rc=RC) call ESMF_Finalize(endflag=ESMF_END_ABORT) endif wrt_int_state%out_grid_info(n)%i_start = lbound(lonPtr,1) wrt_int_state%out_grid_info(n)%i_end = ubound(lonPtr,1) wrt_int_state%out_grid_info(n)%j_start = lbound(latPtr,2) wrt_int_state%out_grid_info(n)%j_end = ubound(latPtr,2) allocate( wrt_int_state%out_grid_info(n)%i_start_wrtgrp(wrt_int_state%petcount) ) allocate( wrt_int_state%out_grid_info(n)%i_end_wrtgrp (wrt_int_state%petcount) ) allocate( wrt_int_state%out_grid_info(n)%j_start_wrtgrp(wrt_int_state%petcount) ) allocate( wrt_int_state%out_grid_info(n)%j_end_wrtgrp (wrt_int_state%petcount) ) call mpi_allgather(wrt_int_state%out_grid_info(n)%i_start, 1, MPI_INTEGER, & wrt_int_state%out_grid_info(n)%i_start_wrtgrp, 1, MPI_INTEGER, wrt_mpi_comm, rc) call mpi_allgather(wrt_int_state%out_grid_info(n)%i_end, 1, MPI_INTEGER, & wrt_int_state%out_grid_info(n)%i_end_wrtgrp, 1, MPI_INTEGER, wrt_mpi_comm, rc) call mpi_allgather(wrt_int_state%out_grid_info(n)%j_start, 1, MPI_INTEGER, & wrt_int_state%out_grid_info(n)%j_start_wrtgrp, 1, MPI_INTEGER, wrt_mpi_comm, rc) call mpi_allgather(wrt_int_state%out_grid_info(n)%j_end, 1, MPI_INTEGER, & wrt_int_state%out_grid_info(n)%j_end_wrtgrp, 1, MPI_INTEGER, wrt_mpi_comm, rc) allocate( wrt_int_state%out_grid_info(n)%lonPtr(wrt_int_state%out_grid_info(n)%i_start:wrt_int_state%out_grid_info(n)%i_end, & wrt_int_state%out_grid_info(n)%j_start:wrt_int_state%out_grid_info(n)%j_end) ) allocate( wrt_int_state%out_grid_info(n)%latPtr(wrt_int_state%out_grid_info(n)%i_start:wrt_int_state%out_grid_info(n)%i_end, & wrt_int_state%out_grid_info(n)%j_start:wrt_int_state%out_grid_info(n)%j_end) ) do j=wrt_int_state%out_grid_info(n)%j_start, wrt_int_state%out_grid_info(n)%j_end do i=wrt_int_state%out_grid_info(n)%i_start, wrt_int_state%out_grid_info(n)%i_end wrt_int_state%out_grid_info(n)%latPtr(i,j) = latPtr(i,j) wrt_int_state%out_grid_info(n)%lonPtr(i,j) = lonPtr(i,j) enddo enddo wrt_int_state%out_grid_info(n)%im = imo(n) wrt_int_state%out_grid_info(n)%jm = jmo(n) end if ! non 'cubed_sphere_grid' end do ! n = 1, ngrids ! !----------------------------------------------------------------------- !*** get write grid component initial time from clock !----------------------------------------------------------------------- ! call ESMF_ClockGet(clock =CLOCK & !<-- The ESMF Clock ,startTime=wrt_int_state%IO_BASETIME & !<-- The Clock's starting time ,rc =RC) call ESMF_TimeGet(time=wrt_int_state%IO_BASETIME,yy=idate(1),mm=idate(2),dd=idate(3),h=idate(4), & m=idate(5),s=idate(6),rc=rc) ! if (lprnt) write(0,*) 'in wrt initial, io_baseline time=',idate,'rc=',rc idate(7) = 1 wrt_int_state%idate = idate wrt_int_state%fdate = idate ! update IO-BASETIME and idate on write grid comp when IAU is enabled if(iau_offset > 0 ) then call ESMF_TimeIntervalSet(IAU_offsetTI, h=iau_offset, rc=rc) wrt_int_state%IO_BASETIME = wrt_int_state%IO_BASETIME + IAU_offsetTI call ESMF_TimeGet(time=wrt_int_state%IO_BASETIME,yy=idate(1),mm=idate(2),dd=idate(3),h=idate(4), & m=idate(5),s=idate(6),rc=rc) wrt_int_state%idate = idate change_wrtidate = .true. if (lprnt) print *,'in wrt initial, with iau, io_baseline time=',idate,'rc=',rc endif ! !--- Look at the incoming FieldBundles in the imp_state_write, and mirror them as 'output_' bundles ! call ESMF_StateGet(imp_state_write, itemCount=FBCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return wrt_int_state%FBCount = FBCount ! if (lprnt) write(0,*) 'in wrt,fcst FBCount=',FBCount ! if (lprnt) write(0,*) 'in wrt,fcst wrt_int_state%FBCount=',wrt_int_state%FBCount allocate(fcstItemNameList(FBCount), fcstItemTypeList(FBCount)) allocate(wrt_int_state%wrtFB_names(FBCount)) allocate(wrt_int_state%ncount_fields(FBCount),wrt_int_state%ncount_attribs(FBCount)) allocate(wrt_int_state%field_names(2000,FBCount)) allocate(outfilename(2000,FBcount)) outfilename = '' call ESMF_StateGet(imp_state_write, itemNameList=fcstItemNameList, & itemTypeList=fcstItemTypeList, & !itemorderflag=ESMF_ITEMORDER_ADDORDER, & rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return !loop over all items in the imp_state_write and collect all FieldBundles do i=1, FBCount if (fcstItemTypeList(i) == ESMF_STATEITEM_FIELDBUNDLE) then call ESMF_StateGet(imp_state_write, itemName=fcstItemNameList(i), & fieldbundle=fcstFB, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeGet(fcstFB, convention="NetCDF", purpose="FV3", & name="grid_id", value=grid_id, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return !--- get grid dim count call ESMF_GridGet(wrtGrid(grid_id), dimCount=gridDimCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! create a mirrored 'output_' FieldBundle and add it to importState fieldbundle = ESMF_FieldBundleCreate(name="output_"//trim(fcstItemNameList(i)), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_StateAdd(imp_state_write, (/fieldbundle/), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! copy the fcstFB Attributes to the 'output_' FieldBundle call ESMF_AttributeCopy(fcstFB, fieldbundle, attcopy=ESMF_ATTCOPY_REFERENCE, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! grids in fcstFB for which 'is_moving' is .true. must provide a first level mirror for the Redist() target if (is_moving(grid_id)) then ! create a mirrored 'mirror_' FieldBundle and add it to importState mirrorFB = ESMF_FieldBundleCreate(name="mirror_"//trim(fcstItemNameList(i)), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_StateAdd(imp_state_write, (/mirrorFB/), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! copy the fcstFB Attributes to the 'mirror_' FieldBundle call ESMF_AttributeCopy(fcstFB, mirrorFB, attcopy=ESMF_ATTCOPY_REFERENCE, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return endif ! deal with all of the Fields inside this fcstFB call ESMF_FieldBundleGet(fcstFB, fieldCount=fieldCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if (fieldCount > 0) then wrt_int_state%ncount_fields(i) = fieldCount allocate(fcstField(fieldCount)) call ESMF_FieldBundleGet(fcstFB, fieldList=fcstField, & itemorderflag=ESMF_ITEMORDER_ADDORDER, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return do j=1, fieldCount call ESMF_FieldGet(fcstField(j), typekind=typekind, & dimCount=fieldDimCount, name=fieldName, grid=fcstGrid, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return allocate(gridToFieldMap(gridDimCount)) allocate(ungriddedLBound(fieldDimCount-gridDimCount)) allocate(ungriddedUBound(fieldDimCount-gridDimCount)) call ESMF_FieldGet(fcstField(j), gridToFieldMap=gridToFieldMap, & ungriddedLBound=ungriddedLBound, ungriddedUBound=ungriddedUBound, & rc=rc) call ESMF_LogWrite("after field create on wrt comp",ESMF_LOGMSG_INFO,rc=RC) ! if (lprnt) print *,'in wrt,fcstfld,fieldname=', & ! trim(fieldname),'fieldDimCount=',fieldDimCount,'gridDimCount=',gridDimCount, & ! 'gridToFieldMap=',gridToFieldMap,'ungriddedLBound=',ungriddedLBound, & ! 'ungriddedUBound=',ungriddedUBound,'rc=',rc if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! create the output field call ESMF_LogWrite("call field create on wrt comp",ESMF_LOGMSG_INFO,rc=RC) field_work = ESMF_FieldCreate(wrtGrid(grid_id), typekind, name=fieldName, & gridToFieldMap=gridToFieldMap, & ungriddedLBound=ungriddedLBound, & ungriddedUBound=ungriddedUBound, rc=rc) call ESMF_LogWrite("aft call field create on wrt comp",ESMF_LOGMSG_INFO,rc=RC) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return wrt_int_state%field_names(j,i) = trim(fieldName) call ESMF_AttributeCopy(fcstField(j), field_work, attcopy=ESMF_ATTCOPY_REFERENCE, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! ! get output file name call ESMF_AttributeGet(fcstField(j), convention="NetCDF", purpose="FV3", & name="output_file", value=outfile_name, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_LogWrite("bf fcstfield, get output_file "//trim(outfile_name)//" "//trim(fieldName),ESMF_LOGMSG_INFO,rc=RC) if (trim(outfile_name) /= '') then outfilename(j,i) = trim(outfile_name) endif call ESMF_LogWrite("af fcstfield, get output_file",ESMF_LOGMSG_INFO,rc=RC) ! if (lprnt) print *,' i=',i,' j=',j,' outfilename=',trim(outfilename(j,i)) ! add the output field to the 'output_' FieldBundle call ESMF_FieldBundleAdd(fieldbundle, (/field_work/), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! deal with grids for which 'is_moving' is .true. if (is_moving(grid_id)) then ! create an empty field that will serve as acceptor for GridTransfer of fcstGrid field_work = ESMF_FieldEmptyCreate(name=fieldName, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! use attributes to carry information for later FieldEmptyComplete() call ESMF_InfoGetFromHost(field_work, info=info, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return tk = typekind ! convert TypeKind_Flag to integer call ESMF_InfoSet(info, key="typekind", value=tk, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_InfoSet(info, key="gridToFieldMap", values=gridToFieldMap, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_InfoSet(info, key="ungriddedLBound", values=ungriddedLBound, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_InfoSet(info, key="ungriddedUBound", values=ungriddedUBound, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! add to 'mirror_' FieldBundle call ESMF_FieldBundleAdd(mirrorFB, (/field_work/), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return endif ! local garbage collection deallocate(gridToFieldMap, ungriddedLBound, ungriddedUBound) enddo ! call ESMF_AttributeCopy(fcstGrid, wrtGrid(grid_id), & attcopy=ESMF_ATTCOPY_REFERENCE, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return deallocate(fcstField) endif !if (fieldCount > 0) then else ! anything but a FieldBundle in the state is unexpected here call ESMF_LogSetError(ESMF_RC_ARG_BAD, & msg="Only FieldBundles supported in fcstState.", line=__LINE__, file=__FILE__) return endif !end FBCount enddo ! !loop over all items in the imp_state_write and count output FieldBundles call get_outfile(FBcount, outfilename,FBlist_outfilename,noutfile) wrt_int_state%FBCount = noutfile ! !create output field bundles allocate(wrt_int_state%wrtFB(wrt_int_state%FBcount)) do i=1, wrt_int_state%FBcount wrt_int_state%wrtFB_names(i) = trim(FBlist_outfilename(i)) wrt_int_state%wrtFB(i) = ESMF_FieldBundleCreate(name=trim(wrt_int_state%wrtFB_names(i)), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return do n=1, FBcount call ESMF_StateGet(imp_state_write, itemName="output_"//trim(fcstItemNameList(n)), & fieldbundle=fcstFB, rc=rc) if( index(trim(fcstItemNameList(n)),trim(FBlist_outfilename(i))) == 1 ) then ! ! copy the fcstfield bundle Attributes to the output field bundle call ESMF_AttributeCopy(fcstFB, wrt_int_state%wrtFB(i), & attcopy=ESMF_ATTCOPY_REFERENCE, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeGet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="grid_id", value=grid_id, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! call ESMF_FieldBundleGet(fcstFB, fieldCount=fieldCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return allocate(fcstField(fieldCount),fieldnamelist(fieldCount)) call ESMF_FieldBundleGet(fcstFB, fieldList=fcstField, fieldNameList=fieldnamelist, & itemorderflag=ESMF_ITEMORDER_ADDORDER, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return do j=1, fieldCount call ESMF_AttributeGet(fcstField(j),convention="NetCDF", purpose="FV3", & name='output_file',value=outfile_name, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! if (lprnt) print *,'in wrt,add field,i=',i,'n=',n,' j=',j, & ! 'fieldname=',trim(fieldnamelist(j)), ' outfile_name=',trim(outfile_name), & ! ' field bundle name, FBlist_outfilename(i)=',trim(FBlist_outfilename(i)) if( trim(outfile_name) == trim(FBlist_outfilename(i))) then call ESMF_FieldBundleAdd(wrt_int_state%wrtFB(i), (/fcstField(j)/), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return endif enddo deallocate(fcstField, fieldnamelist) endif ! add output grid related attributes call ESMF_AttributeAdd(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & attrList=(/"source","grid "/), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="source", value="FV3GFS", rc=rc) if (trim(output_grid(grid_id)) == 'cubed_sphere_grid') then call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="grid", value="cubed_sphere", rc=rc) else if (trim(output_grid(grid_id)) == 'gaussian_grid') then call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="grid", value="gaussian", rc=rc) call ESMF_AttributeAdd(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & attrList=(/"im","jm"/), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="im", value=imo(grid_id), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="jm", value=jmo(grid_id), rc=rc) else if (trim(output_grid(grid_id)) == 'regional_latlon' & .or. trim(output_grid(grid_id)) == 'regional_latlon_moving' & .or. trim(output_grid(grid_id)) == 'global_latlon') then ! for 'regional_latlon_moving' lon1/2 and lat1/2 will be overwritten in run phase call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="grid", value="latlon", rc=rc) call ESMF_AttributeAdd(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & attrList=(/"lon1","lat1","lon2","lat2","dlon","dlat"/), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="lon1", value=lon1(grid_id), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="lat1", value=lat1(grid_id), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="lon2", value=lon2(grid_id), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="lat2", value=lat2(grid_id), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="dlon", value=dlon(grid_id), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="dlat", value=dlat(grid_id), rc=rc) else if (trim(output_grid(grid_id)) == 'rotated_latlon' & .or. trim(output_grid(grid_id)) == 'rotated_latlon_moving') then call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="grid", value="rotated_latlon", rc=rc) call ESMF_AttributeAdd(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & attrList=(/"cen_lon",& "cen_lat",& "lon1 ",& "lat1 ",& "lon2 ",& "lat2 ",& "dlon ",& "dlat "/), rc=rc) ! for 'rotated_latlon_moving' cen_lon and cen_lat will be overwritten in run phase call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="cen_lon", value=cen_lon(grid_id), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="cen_lat", value=cen_lat(grid_id), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="lon1", value=lon1(grid_id), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="lat1", value=lat1(grid_id), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="lon2", value=lon2(grid_id), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="lat2", value=lat2(grid_id), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="dlon", value=dlon(grid_id), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="dlat", value=dlat(grid_id), rc=rc) else if (trim(output_grid(grid_id)) == 'lambert_conformal') then call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="grid", value="lambert_conformal", rc=rc) call ESMF_AttributeAdd(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & attrList=(/"cen_lon",& "cen_lat",& "stdlat1",& "stdlat2",& "nx ",& "ny ",& "lon1 ",& "lat1 ",& "dx ",& "dy "/), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="cen_lon", value=cen_lon(grid_id), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="cen_lat", value=cen_lat(grid_id), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="stdlat1", value=stdlat1(grid_id), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="stdlat2", value=stdlat2(grid_id), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="nx", value=imo(grid_id), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="ny", value=jmo(grid_id), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="lat1", value=lat1(grid_id), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="lon1", value=lon1(grid_id), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="dx", value=dx(grid_id), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="dy", value=dy(grid_id), rc=rc) end if enddo ! end FBcount enddo ! end wrt_int_state%FBcount ! ! add time Attribute ! look at the importState attributes and copy those starting with "time" call ESMF_AttributeGet(imp_state_write, convention="NetCDF", purpose="FV3", & attnestflag=ESMF_ATTNEST_OFF, count=attCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_LogWrite("Write component AttributeGet, attCount ", ESMF_LOGMSG_INFO, rc=rc) ! prepare the lists needed to transfer attributes allocate(attNameList(attCount), attNameList2(attCount)) allocate(typekindList(attCount)) ! loop over all the attributes on importState within AttPack j = 1 k = 1 do i=1, attCount call ESMF_AttributeGet(imp_state_write, convention="NetCDF", purpose="FV3", & attnestflag=ESMF_ATTNEST_OFF, attributeIndex=i, name=attName, & typekind=typekind, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! test for name starting with "time" if (index(trim(attName), "time") == 1) then ! add this attribute to the list of transfers attNameList(j) = attName typekindList(j) = typekind j = j + 1 if (index(trim(attName), "time:") == 1) then ! store names of attributes starting with "time:" for later use attNameList2(k) = attName k = k+1 endif endif enddo do n = 1, ngrids ! add the transfer attributes from importState to grid call ESMF_AttributeAdd(wrtGrid(n), convention="NetCDF", purpose="FV3", & attrList=attNameList(1:j-1), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! loop over the added attributes, access the value (only scalar allowed), ! and set them on the grid do i=1, j-1 if (typekindList(i) == ESMF_TYPEKIND_CHARACTER) then call ESMF_AttributeGet(imp_state_write, & convention="NetCDF", purpose="FV3", & name=trim(attNameList(i)), value=valueS, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! update the time:units when idate on write grid component is changed if ( index(trim(attNameList(i)),'time:units')>0) then if ( change_wrtidate ) then idx = index(trim(valueS),' since ') if(lprnt) print *,'in write grid comp, time:unit=',trim(valueS) write(newdate,'(I4.4,a,I2.2,a,I2.2,a,I2.2,a,I2.2,a,I2.2)') idate(1),'-', & idate(2),'-',idate(3),' ',idate(4),':',idate(5),':',idate(6) valueS = valueS(1:idx+6)//newdate if(lprnt) print *,'in write grid comp, new time:unit=',trim(valueS) endif endif call ESMF_AttributeSet(wrtGrid(n), convention="NetCDF", purpose="FV3", & name=trim(attNameList(i)), value=valueS, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return else if (typekindList(i) == ESMF_TYPEKIND_I4) then call ESMF_AttributeGet(imp_state_write, & convention="NetCDF", purpose="FV3", & name=trim(attNameList(i)), value=valueI4, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeSet(wrtGrid(n), convention="NetCDF", purpose="FV3", & name=trim(attNameList(i)), value=valueI4, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return else if (typekindList(i) == ESMF_TYPEKIND_R4) then call ESMF_AttributeGet(imp_state_write, & convention="NetCDF", purpose="FV3", & name=trim(attNameList(i)), value=valueR4, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeSet(wrtGrid(n), convention="NetCDF", purpose="FV3", & name=trim(attNameList(i)), value=valueR4, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return else if (typekindList(i) == ESMF_TYPEKIND_R8) then call ESMF_AttributeGet(imp_state_write, & convention="NetCDF", purpose="FV3", & name=trim(attNameList(i)), value=valueR8, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeSet(wrtGrid(n), convention="NetCDF", purpose="FV3", & name=trim(attNameList(i)), value=valueR8, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return endif enddo ! Add special attribute that holds names of "time" related attributes ! for faster access during Run(). call ESMF_AttributeAdd(wrtGrid(n), convention="NetCDF", purpose="FV3", & attrList=(/"TimeAttributes"/), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeSet(wrtGrid(n), convention="NetCDF", purpose="FV3", & name="TimeAttributes", valueList=attNameList2(1:k-1), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! !*** create temporary field bundle for axes information ! write the Grid coordinate arrays into the output files via temporary FB gridFB = ESMF_FieldBundleCreate(rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_GridGetCoord(wrtGrid(n), coordDim=1, & staggerloc=ESMF_STAGGERLOC_CENTER, array=array, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return field = ESMF_FieldCreate(wrtGrid(n), array, name="grid_xt", rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return !add attribute info ! long name call ESMF_AttributeAdd(field,convention="NetCDF",purpose="FV3", & attrList=(/'long_name'/), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeSet(field,convention="NetCDF",purpose="FV3",name='long_name', & value="T-cell longitude", rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! units call ESMF_AttributeAdd(field,convention="NetCDF",purpose="FV3", & attrList=(/'units'/), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeSet(field,convention="NetCDF",purpose="FV3",name='units', & value="degrees_E", rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! cartesian_axis call ESMF_AttributeAdd(field,convention="NetCDF",purpose="FV3", & attrList=(/'cartesian_axis'/), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeSet(field,convention="NetCDF",purpose="FV3",name='cartesian_axis', & value="X", rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! add field to bundle call ESMF_FieldBundleAdd(gridFB, (/field/), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! ! get 2nd dimension call ESMF_GridGetCoord(wrtGrid(n), coordDim=2, & staggerloc=ESMF_STAGGERLOC_CENTER, array=array, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return field = ESMF_FieldCreate(wrtGrid(n), array, name="grid_yt", rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return !add attribute info ! long name call ESMF_AttributeAdd(field,convention="NetCDF",purpose="FV3", & attrList=(/'long_name'/), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeSet(field,convention="NetCDF",purpose="FV3",name='long_name', & value="T-cell latitude", rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! units call ESMF_AttributeAdd(field,convention="NetCDF",purpose="FV3", & attrList=(/'units'/), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeSet(field,convention="NetCDF",purpose="FV3",name='units', & value="degrees_N", rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! cartesian_axis call ESMF_AttributeAdd(field,convention="NetCDF",purpose="FV3", & attrList=(/'cartesian_axis'/), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeSet(field,convention="NetCDF",purpose="FV3",name='cartesian_axis', & value="Y", rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_FieldBundleAdd(gridFB, (/field/), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return end do ! n=1, ngrids deallocate(attNameList, attNameList2, typekindList) ! ! write_init_tim = MPI_Wtime() - btim0 ! !----------------------------------------------------------------------- ! end subroutine wrt_initialize_p1 ! !----------------------------------------------------------------------- !####################################################################### !----------------------------------------------------------------------- ! subroutine wrt_initialize_p2(wrt_comp, imp_state_write, exp_state_write, clock, rc) ! !----------------------------------------------------------------------- !*** INITIALIZE THE WRITE GRIDDED COMPONENT. !----------------------------------------------------------------------- ! type(esmf_GridComp) :: wrt_comp type(ESMF_State) :: imp_state_write, exp_state_write type(esmf_Clock) :: clock integer,intent(out) :: rc ! !*** LOCAL VARIABLES type(ESMF_Info) :: info logical, allocatable :: is_moving(:) type(ESMF_StateItem_Flag), allocatable :: itemTypeList(:) character(len=ESMF_MAXSTR),allocatable :: itemNameList(:) integer :: i, j, bundleCount, fieldCount type(ESMF_FieldBundle) :: mirrorFB type(ESMF_Field), allocatable :: fieldList(:) type(ESMF_Grid) :: grid type(ESMF_DistGrid) :: acceptorDG, newAcceptorDG ! ! !----------------------------------------------------------------------- !*********************************************************************** !----------------------------------------------------------------------- ! rc = ESMF_SUCCESS ! call ESMF_InfoGetFromHost(imp_state_write, info=info, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_InfoGetAlloc(info, key="is_moving", values=is_moving, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_StateGet(imp_state_write, itemCount=bundleCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return allocate(itemNameList(bundleCount), itemTypeList(bundleCount)) call ESMF_StateGet(imp_state_write, itemNameList=itemNameList, & itemTypeList=itemTypeList, & !itemorderflag=ESMF_ITEMORDER_ADDORDER, & rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return do i=1, bundleCount if (itemTypeList(i) == ESMF_STATEITEM_FIELDBUNDLE) then if (index(trim(itemNameList(i)), "mirror_")==1) then ! this is a 'mirror_' FieldBundle -> GridTransfer acceptor side call ESMF_StateGet(imp_state_write, itemName=trim(itemNameList(i)), fieldbundle=mirrorFB, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! access the grid that is passed in from the provider side call ESMF_FieldBundleGet(mirrorFB, grid=grid, fieldCount=fieldCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! access the acceptor DistGrid call ESMF_GridGet(grid, distgrid=acceptorDG, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! rebalance the acceptor DistGrid across the local PETs newAcceptorDG = ESMF_DistGridCreate(acceptorDG, balanceflag=.true., rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! create a new Grid on the rebalanced DistGrid grid = ESMF_GridCreate(newAcceptorDG, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! point all of the acceptor fields to the new acceptor Grid allocate(fieldList(fieldCount)) call ESMF_FieldBundleGet(mirrorFB, fieldList=fieldList, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return do j=1, fieldCount call ESMF_FieldEmptySet(fieldList(j), grid=grid, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return enddo ! clean-up deallocate(fieldList) endif else ! anything but a FieldBundle in the state is unexpected here call ESMF_LogSetError(ESMF_RC_ARG_BAD, msg="Only FieldBundles supported in fcstState.", line=__LINE__, file=__FILE__) return endif enddo !----------------------------------------------------------------------- ! end subroutine wrt_initialize_p2 ! !----------------------------------------------------------------------- !####################################################################### !----------------------------------------------------------------------- ! subroutine wrt_initialize_p3(wrt_comp, imp_state_write, exp_state_write, clock, rc) ! !----------------------------------------------------------------------- !*** INITIALIZE THE WRITE GRIDDED COMPONENT. !----------------------------------------------------------------------- ! type(esmf_GridComp) :: wrt_comp type(ESMF_State) :: imp_state_write, exp_state_write type(esmf_Clock) :: clock integer,intent(out) :: rc !*** LOCAL VARIABLES type(ESMF_Info) :: info logical, allocatable :: is_moving(:) type(ESMF_StateItem_Flag), allocatable :: itemTypeList(:) character(len=ESMF_MAXSTR),allocatable :: itemNameList(:) integer :: i, j, bundleCount, fieldCount type(ESMF_FieldBundle) :: mirrorFB type(ESMF_Field), allocatable :: fieldList(:) type(ESMF_TypeKind_Flag) :: typekind integer :: tk integer, allocatable :: gridToFieldMap(:) integer, allocatable :: ungriddedLBound(:) integer, allocatable :: ungriddedUBound(:) ! ! !----------------------------------------------------------------------- !*********************************************************************** !----------------------------------------------------------------------- ! rc = ESMF_SUCCESS ! call ESMF_InfoGetFromHost(imp_state_write, info=info, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_InfoGetAlloc(info, key="is_moving", values=is_moving, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_StateGet(imp_state_write, itemCount=bundleCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return allocate(itemNameList(bundleCount), itemTypeList(bundleCount)) call ESMF_StateGet(imp_state_write, itemNameList=itemNameList, & itemTypeList=itemTypeList, & !itemorderflag=ESMF_ITEMORDER_ADDORDER, & rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return do i=1, bundleCount if (itemTypeList(i) == ESMF_STATEITEM_FIELDBUNDLE) then if (index(trim(itemNameList(i)), "mirror_")==1) then ! this is a 'mirror_' FieldBundle -> GridTransfer acceptor side call ESMF_StateGet(imp_state_write, itemName=trim(itemNameList(i)), fieldbundle=mirrorFB, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! finish creating all the mirror Fields call ESMF_FieldBundleGet(mirrorFB, fieldCount=fieldCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return allocate(fieldList(fieldCount)) call ESMF_FieldBundleGet(mirrorFB, fieldList=fieldList, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return do j=1, fieldCount ! first access information stored on the field needed for completion call ESMF_InfoGetFromHost(fieldList(j), info=info, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_InfoGet(info, key="typekind", value=tk, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return typekind = tk ! convert integer into TypeKind_Flag call ESMF_InfoGetAlloc(info, key="gridToFieldMap", values=gridToFieldMap, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_InfoGetAlloc(info, key="ungriddedLBound", values=ungriddedLBound, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_InfoGetAlloc(info, key="ungriddedUBound", values=ungriddedUBound, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! now complete the field creation call ESMF_FieldEmptyComplete(fieldList(j), typekind=typekind, gridToFieldMap=gridToFieldMap, & ungriddedLBound=ungriddedLBound, ungriddedUBound=ungriddedUBound, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return enddo ! clean-up deallocate(fieldList) endif else ! anything but a FieldBundle in the state is unexpected here call ESMF_LogSetError(ESMF_RC_ARG_BAD, msg="Only FieldBundles supported in fcstState.", line=__LINE__, file=__FILE__) return endif enddo !----------------------------------------------------------------------- ! end subroutine wrt_initialize_p3 ! !----------------------------------------------------------------------- !####################################################################### !----------------------------------------------------------------------- ! subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) ! !----------------------------------------------------------------------- !*** the run step for the write gridded component. !----------------------------------------------------------------------- ! type(ESMF_GridComp) :: wrt_comp type(ESMF_State) :: imp_state_write, exp_state_write type(ESMF_Clock) :: clock integer,intent(out) :: rc ! !----------------------------------------------------------------------- !*** local variables ! TYPE(ESMF_VM) :: VM type(ESMF_FieldBundle) :: file_bundle, mirror_bundle type(ESMF_StateItem_Flag) :: itemType type(ESMF_Time) :: currtime type(ESMF_TimeInterval) :: io_currtimediff type(ESMF_Grid) :: fbgrid, wrtGrid type(ESMF_State),save :: stateGridFB type(optimizeT), save :: optimize(4) type(ESMF_GridComp), save, allocatable :: compsGridFB(:) type(ESMF_RouteHandle) :: rh type(ESMF_RegridMethod_Flag) :: regridmethod integer :: srcTermProcessing ! type(write_wrap) :: wrap type(wrt_internal_state),pointer :: wrt_int_state ! integer :: i,j,n,mype,nolog, grid_id, localPet ! integer :: nf_hours,nf_seconds,nf_minutes real(ESMF_KIND_R8) :: nfhour ! integer :: nbdl, date(6), ndig, nnnn integer :: step=1 ! logical :: opened logical :: lmask_fields ! character(esmf_maxstr) :: filename,compname, traceString character(40) :: cfhour, cform character(20) :: time_iso ! type(ESMF_Grid) :: grid type(ESMF_Info) :: info real(ESMF_KIND_R8), allocatable :: values(:) character(160) :: msgString type(ESMF_Field), allocatable :: fieldList(:) type(ESMF_Array) :: coordArray(2) type(ESMF_DistGrid) :: coordDG type(ESMF_DELayout) :: coordDL integer :: fieldCount, deCount, rootPet integer :: minIndexPTile(2,1), maxIndexPTile(2,1), centerIndex(2) integer, allocatable :: minIndexPDe(:,:), maxIndexPDe(:,:), petMap(:) real(ESMF_KIND_R8), pointer :: farrayPtr(:,:) real(ESMF_KIND_R8) :: centerCoord(2) integer :: ii, jj real(ESMF_KIND_R8), pointer :: lonPtr(:,:), latPtr(:,:) real(ESMF_KIND_R8) :: rot_lon, rot_lat real(ESMF_KIND_R8) :: geo_lon, geo_lat real(ESMF_KIND_R8), parameter :: rtod=180.0/pi real(kind=8) :: MPI_Wtime real(kind=8) :: tbeg real(kind=8) :: wbeg,wend logical :: use_parallel_netcdf logical :: lprnt ! !----------------------------------------------------------------------- !*********************************************************************** !----------------------------------------------------------------------- ! tbeg = MPI_Wtime() rc = esmf_success ! !----------------------------------------------------------------------- !*** get the current write grid comp name, and internal state ! call ESMF_GridCompGet(wrt_comp, name=compname, localPet=localPet, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! Provide log message indicating which wrtComp is active call ESMF_LogWrite("Write component activated: "//trim(compname), & ESMF_LOGMSG_INFO, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! access the internal state call ESMF_GridCompGetInternalState(wrt_Comp, wrap, rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return wrt_int_state => wrap%write_int_state call ESMF_VMGetCurrent(VM,rc=RC) mype = wrt_int_state%mype lprnt = mype == lead_write_task ! print *,'in wrt run, mype=',mype,'lead_write_task=',lead_write_task ! !----------------------------------------------------------------------- !*** get current time and elapsed forecast time call ESMF_ClockGet(clock=CLOCK, currTime=CURRTIME, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_TimeGet(time=currTime,yy=date(1),mm=date(2),dd=date(3),h=date(4), & m=date(5),s=date(6),rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return wrt_int_state%fdate(7) = 1 wrt_int_state%fdate(1:6) = date(1:6) write(time_iso,'(I4,"-",I2.2,"-",I2.2,"T",I2.2,":",I2.2,":",I2.2,"Z")') date(1:6) call ESMF_TimeGet(time=wrt_int_state%IO_BASETIME,yy=date(1),mm=date(2),dd=date(3),h=date(4), & m=date(5),s=date(6),rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return io_currtimediff = currtime - wrt_int_state%IO_BASETIME call ESMF_TimeIntervalGet(timeinterval=io_currtimediff & ,h_r8=nfhour,h=nf_hours,m=nf_minutes,s=nf_seconds,rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if (nf_hours < 0) return if (nsout > 0 .or. lflname_fulltime) then ndig = max(log10(nf_hours+0.5)+1., 3.) write(cform, '("(I",I1,".",I1,",A1,I2.2,A1,I2.2)")') ndig, ndig write(cfhour, cform) nf_hours,'-',nf_minutes,'-',nf_seconds else ndig = max(log10(nf_hours+0.5)+1., 3.) write(cform, '("(I",I1,".",I1,")")') ndig, ndig write(cfhour, cform) nf_hours endif ! if(lprnt) print *,'in wrt run, nfhour=',nfhour,' cfhour=',trim(cfhour) ! !----------------------------------------------------------------------- !*** loop on the "output_" FieldBundles, i.e. files that need to write out !----------------------------------------------------------------------- do i=1, FBCount call ESMF_StateGet(imp_state_write, itemName="output_"//trim(fcstItemNameList(i)), & fieldbundle=file_bundle, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! see whether a "mirror_" FieldBundle exists, i.e. dealing with moving domain that needs updated Regrid() here. call ESMF_StateGet(imp_state_write, itemName="mirror_"//trim(fcstItemNameList(i)), & itemType=itemType, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if (itemType == ESMF_STATEITEM_FIELDBUNDLE) then ! Regrid() for a moving domain call ESMF_LogWrite("Regrid() for moving domain: mirror_"//trim(fcstItemNameList(i)), ESMF_LOGMSG_INFO, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_StateGet(imp_state_write, itemName="mirror_"//trim(fcstItemNameList(i)), & fieldbundle=mirror_bundle, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! Find the centerCoord of the moving domain call ESMF_FieldBundleGet(mirror_bundle, fieldCount=fieldCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return allocate(fieldList(fieldCount)) call ESMF_FieldBundleGet(mirror_bundle, fieldList=fieldList, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_FieldGet(fieldList(1), grid=grid, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return deallocate(fieldList) call ESMF_GridGetCoord(grid, coordDim=1, array=coordArray(1), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_GridGetCoord(grid, coordDim=2, array=coordArray(2), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_ArrayGet(coordArray(1), distgrid=coordDG, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_DistGridGet(coordDG, deCount=deCount, minIndexPTile=minIndexPTile, maxIndexPTile=maxIndexPTile, & delayout=coordDL, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return allocate(petMap(deCount),minIndexPDe(2,deCount), maxIndexPDe(2,deCount)) call ESMF_DELayoutGet(coordDL, petMap=petMap, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_DistGridGet(coordDG, minIndexPDe=minIndexPDe, maxIndexPDe=maxIndexPDe, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return centerIndex(1) = (maxIndexPTile(1,1)-minIndexPTile(1,1)+1)/2 centerIndex(2) = (maxIndexPTile(2,1)-minIndexPTile(2,1)+1)/2 ! write(msgString,*) "Determined centerIndex: ", centerIndex ! call ESMF_LogWrite(trim(msgString), ESMF_LOGMSG_DEBUG, rc=rc) ! if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return do n=1, deCount if (minIndexPDe(1,n)<=centerIndex(1) .and. centerIndex(1)<=maxIndexPDe(1,n) .and. & minIndexPDe(2,n)<=centerIndex(2) .and. centerIndex(2)<=maxIndexPDe(2,n)) then ! found the DE that holds the center coordinate rootPet = petMap(n) if (localPet == rootPet) then ! center DE is on local PET -> fill centerCoord locally call ESMF_ArrayGet(coordArray(1), farrayPtr=farrayPtr, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return centerCoord(1) = farrayPtr(centerIndex(1)-minIndexPDe(1,n)+1,centerIndex(2)-minIndexPDe(2,n)+1) call ESMF_ArrayGet(coordArray(2), farrayPtr=farrayPtr, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return centerCoord(2) = farrayPtr(centerIndex(1)-minIndexPDe(1,n)+1,centerIndex(2)-minIndexPDe(2,n)+1) ! write(msgString,*) "Found centerCoord: ", centerCoord ! call ESMF_LogWrite(trim(msgString), ESMF_LOGMSG_DEBUG, rc=rc) ! if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return endif exit endif enddo deallocate(petMap,minIndexPDe,maxIndexPDe) call ESMF_VMBroadcast(vm, centerCoord, count=2, rootPet=rootPet, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return write(msgString,*) "All PETs know centerCoord in radians: ", centerCoord call ESMF_LogWrite(trim(msgString), ESMF_LOGMSG_DEBUG, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! determine regridmethod if (index(fcstItemNameList(i),"_bilinear") >0 ) then traceString = "-bilinear" regridmethod = ESMF_REGRIDMETHOD_BILINEAR else if (index(fcstItemNameList(i),"_patch") >0) then traceString = "-patch" regridmethod = ESMF_REGRIDMETHOD_PATCH else if (index(fcstItemNameList(i),"_nearest_stod") >0) then traceString = "-nearest_stod" regridmethod = ESMF_REGRIDMETHOD_NEAREST_STOD else if (index(fcstItemNameList(i),"_nearest_dtos") >0) then traceString = "-nearest_dtos" regridmethod = ESMF_REGRIDMETHOD_NEAREST_DTOS else if (index(fcstItemNameList(i),"_conserve") >0) then traceString = "-conserve" regridmethod = ESMF_REGRIDMETHOD_CONSERVE else call ESMF_LogSetError(ESMF_RC_ARG_BAD, & msg="Unable to determine regrid method.", & line=__LINE__, file=__FILE__, rcToReturn=rc) return endif srcTermProcessing = 1 ! have this fixed for bit-for-bit reproducibility ! RegridStore() ! update output grid coordinates based of fcstgrid center lat/lon call ESMF_FieldBundleGet(file_bundle, grid=grid, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_GridGetCoord(grid, coordDim=1, farrayPtr=lonPtr, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_GridGetCoord(grid, coordDim=2, farrayPtr=latPtr, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeGet(mirror_bundle, convention="NetCDF", purpose="FV3", & name="grid_id", value=grid_id, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if (trim(output_grid(grid_id)) == 'regional_latlon_moving' .or. & trim(output_grid(grid_id)) == 'rotated_latlon_moving') then n = grid_id cen_lon(n) = centerCoord(1)*rtod cen_lat(n) = centerCoord(2)*rtod if (cen_lon(n) > 180.0) cen_lon(n) = cen_lon(n) - 360.0 cen_lon(n) = NINT(cen_lon(n)*1000.0)/1000.0 cen_lat(n) = NINT(cen_lat(n)*1000.0)/1000.0 endif if (trim(output_grid(grid_id)) == 'regional_latlon_moving') then lon1(n) = cen_lon(n) - 0.5 * (imo(n)-1) * dlon(n) lat1(n) = cen_lat(n) - 0.5 * (jmo(n)-1) * dlat(n) lon2(n) = cen_lon(n) + 0.5 * (imo(n)-1) * dlon(n) lat2(n) = cen_lat(n) + 0.5 * (jmo(n)-1) * dlat(n) do jj=lbound(lonPtr,2),ubound(lonPtr,2) do ii=lbound(lonPtr,1),ubound(lonPtr,1) lonPtr(ii,jj) = lon1(n) + (lon2(n)-lon1(n))/(imo(n)-1) * (ii-1) latPtr(ii,jj) = lat1(n) + (lat2(n)-lat1(n))/(jmo(n)-1) * (jj-1) wrt_int_state%out_grid_info(n)%latPtr(ii,jj) = latPtr(ii,jj) wrt_int_state%out_grid_info(n)%lonPtr(ii,jj) = lonPtr(ii,jj) enddo enddo else if (trim(output_grid(grid_id)) == 'rotated_latlon_moving') then lon1(n) = - 0.5 * (imo(n)-1) * dlon(n) lat1(n) = - 0.5 * (jmo(n)-1) * dlat(n) lon2(n) = 0.5 * (imo(n)-1) * dlon(n) lat2(n) = 0.5 * (jmo(n)-1) * dlat(n) do jj=lbound(lonPtr,2),ubound(lonPtr,2) do ii=lbound(lonPtr,1),ubound(lonPtr,1) rot_lon = lon1(n) + (lon2(n)-lon1(n))/(imo(n)-1) * (ii-1) rot_lat = lat1(n) + (lat2(n)-lat1(n))/(jmo(n)-1) * (jj-1) call rtll(rot_lon, rot_lat, geo_lon, geo_lat, dble(cen_lon(n)), dble(cen_lat(n))) if (geo_lon < 0.0) geo_lon = geo_lon + 360.0 lonPtr(ii,jj) = geo_lon latPtr(ii,jj) = geo_lat wrt_int_state%out_grid_info(n)%latPtr(ii,jj) = latPtr(ii,jj) wrt_int_state%out_grid_info(n)%lonPtr(ii,jj) = lonPtr(ii,jj) enddo enddo endif call ESMF_TraceRegionEnter("ESMF_FieldBundleRegridStore()"//trim(traceString), rc=rc) call ESMF_FieldBundleRegridStore(mirror_bundle, file_bundle, & regridMethod=regridmethod, routehandle=rh, & unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, & srcTermProcessing=srcTermProcessing, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_TraceRegionExit("ESMF_FieldBundleRegridStore()"//trim(traceString), rc=rc) ! Regrid() call ESMF_TraceRegionEnter("ESMF_FieldBundleRegrid()"//trim(traceString), rc=rc) call ESMF_FieldBundleRegrid(mirror_bundle, file_bundle, & routehandle=rh, termorderflag=(/ESMF_TERMORDER_SRCSEQ/), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_TraceRegionExit("ESMF_FieldBundleRegrid()"//trim(traceString), rc=rc) ! RegridRelease() call ESMF_FieldBundleRegridRelease(routehandle=rh, noGarbage=.true., rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! done call ESMF_LogWrite("Done Regrid() for moving domain: mirror_"//trim(fcstItemNameList(i)), ESMF_LOGMSG_INFO, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return endif !recover fields from cartesian vector and sfc pressure call recover_fields(file_bundle,rc) enddo ! !----------------------------------------------------------------------- !*** do post !----------------------------------------------------------------------- lmask_fields = .false. if( wrt_int_state%write_dopost ) then #ifdef INLINE_POST wbeg = MPI_Wtime() do n=1,ngrids if (trim(output_grid(n)) == 'regional_latlon' .or. & trim(output_grid(n)) == 'regional_latlon_moving' .or. & trim(output_grid(n)) == 'rotated_latlon' .or. & trim(output_grid(n)) == 'rotated_latlon_moving' .or. & trim(output_grid(n)) == 'lambert_conformal') then !mask fields according to sfc pressure do nbdl=1, wrt_int_state%FBCount call mask_fields(wrt_int_state%wrtFB(nbdl),rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return enddo lmask_fields = .true. endif call post_run_fv3(wrt_int_state, n, mype, wrt_mpi_comm, lead_write_task, & itasks, jtasks, nf_hours, nf_minutes, nf_seconds) enddo wend = MPI_Wtime() if (lprnt) then write(*,'(A,F10.5,A,I4.2,A,I2.2)')' actual inline post Time is ',wend-wbeg & ,' at Fcst ',nf_hours,':',nf_minutes endif #else rc = ESMF_RC_NOT_IMPL print *,'inline post not available on this machine' if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__)) return #endif endif ! !----------------------------------------------------------------------- ! ** now loop through output field bundle !----------------------------------------------------------------------- if ( wrt_int_state%output_history ) then file_loop_all: do nbdl=1, wrt_int_state%FBCount ! ! get grid_id call ESMF_AttributeGet(wrt_int_state%wrtFB(nbdl), convention="NetCDF", purpose="FV3", & name="grid_id", value=grid_id, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! update lon1/2 and lat1/2 for regional_latlon_moving if (trim(output_grid(grid_id)) == 'regional_latlon_moving') then call ESMF_AttributeSet(wrt_int_state%wrtFB(nbdl), convention="NetCDF", purpose="FV3", & name="lon1", value=lon1(grid_id), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeSet(wrt_int_state%wrtFB(nbdl), convention="NetCDF", purpose="FV3", & name="lat1", value=lat1(grid_id), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeSet(wrt_int_state%wrtFB(nbdl), convention="NetCDF", purpose="FV3", & name="lon2", value=lon2(grid_id), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeSet(wrt_int_state%wrtFB(nbdl), convention="NetCDF", purpose="FV3", & name="lat2", value=lat2(grid_id), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return endif ! update cen_lon/cen_lat, lon1/2 and lat1/2 for rotated_latlon_moving if (trim(output_grid(grid_id)) == 'rotated_latlon_moving') then call ESMF_AttributeSet(wrt_int_state%wrtFB(nbdl), convention="NetCDF", purpose="FV3", & name="cen_lon", value=cen_lon(grid_id), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeSet(wrt_int_state%wrtFB(nbdl), convention="NetCDF", purpose="FV3", & name="cen_lat", value=cen_lat(grid_id), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeSet(wrt_int_state%wrtFB(nbdl), convention="NetCDF", purpose="FV3", & name="lon1", value=lon1(grid_id), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeSet(wrt_int_state%wrtFB(nbdl), convention="NetCDF", purpose="FV3", & name="lat1", value=lat1(grid_id), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeSet(wrt_int_state%wrtFB(nbdl), convention="NetCDF", purpose="FV3", & name="lon2", value=lon2(grid_id), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeSet(wrt_int_state%wrtFB(nbdl), convention="NetCDF", purpose="FV3", & name="lat2", value=lat2(grid_id), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return endif if(step == 1) then file_bundle = wrt_int_state%wrtFB(nbdl) endif ! FIXME map nbdl to [1:num_files], only used for output_file nnnn = mod(nbdl-1, num_files) + 1 ! set default chunksizes for netcdf output ! (use MPI decomposition size). ! if chunksize parameter set to negative value, ! netcdf library default is used. if (output_file(nnnn)(1:6) == 'netcdf') then if (ichunk2d(grid_id) == 0) then if( wrt_int_state%mype == 0 ) & ichunk2d(grid_id) = wrt_int_state%out_grid_info(grid_id)%i_end - wrt_int_state%out_grid_info(grid_id)%i_start + 1 call mpi_bcast(ichunk2d(grid_id),1,mpi_integer,0,wrt_mpi_comm,rc) endif if (jchunk2d(grid_id) == 0) then if( wrt_int_state%mype == 0 ) & jchunk2d(grid_id) = wrt_int_state%out_grid_info(grid_id)%j_end - wrt_int_state%out_grid_info(grid_id)%j_start + 1 call mpi_bcast(jchunk2d(grid_id),1,mpi_integer,0,wrt_mpi_comm,rc) endif if (ichunk3d(grid_id) == 0) then if( wrt_int_state%mype == 0 ) & ichunk3d(grid_id) = wrt_int_state%out_grid_info(grid_id)%i_end - wrt_int_state%out_grid_info(grid_id)%i_start + 1 call mpi_bcast(ichunk3d(grid_id),1,mpi_integer,0,wrt_mpi_comm,rc) endif if (jchunk3d(grid_id) == 0) then if( wrt_int_state%mype == 0 ) & jchunk3d(grid_id) = wrt_int_state%out_grid_info(grid_id)%j_end - wrt_int_state%out_grid_info(grid_id)%j_start + 1 call mpi_bcast(jchunk3d(grid_id),1,mpi_integer,0,wrt_mpi_comm,rc) endif if (kchunk3d(grid_id) == 0 .and. nbdl == 1) then if( wrt_int_state%mype == 0 ) then call ESMF_FieldBundleGet(wrt_int_state%wrtFB(nbdl), grid=wrtGrid) call ESMF_AttributeGet(wrtGrid, convention="NetCDF", purpose="FV3", & attnestflag=ESMF_ATTNEST_OFF, name='pfull', & itemCount=kchunk3d(grid_id), rc=rc) endif call mpi_bcast(kchunk3d(grid_id),1,mpi_integer,0,wrt_mpi_comm,rc) endif if (wrt_int_state%mype == 0) then print *,'ichunk2d,jchunk2d',ichunk2d(grid_id),jchunk2d(grid_id) print *,'ichunk3d,jchunk3d,kchunk3d',ichunk3d(grid_id),jchunk3d(grid_id),kchunk3d(grid_id) endif endif filename = trim(wrt_int_state%wrtFB_names(nbdl))//'f'//trim(cfhour)//'.nc' if(mype == lead_write_task) print *,'in wrt run,filename= ',nbdl,trim(filename) ! ! set the time Attribute on the grid to carry it into the lower levels call ESMF_FieldBundleGet(file_bundle, grid=fbgrid, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeSet(fbgrid, convention="NetCDF", purpose="FV3", & name="time", value=nfhour, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeSet(fbgrid, convention="NetCDF", purpose="FV3", & name="time_iso", value=trim(time_iso), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return !*** write out grid bundle: ! Provide log message indicating which wrtComp is active call ESMF_LogWrite("before Write component before gridFB ", ESMF_LOGMSG_INFO, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if (trim(output_file(nnnn)) == 'netcdf') then use_parallel_netcdf = .false. else if (trim(output_file(nnnn)) == 'netcdf_parallel') then use_parallel_netcdf = .true. else call ESMF_LogWrite("wrt_run: Unknown output_file",ESMF_LOGMSG_ERROR,rc=RC) call ESMF_Finalize(endflag=ESMF_END_ABORT) endif if (trim(output_grid(grid_id)) == 'cubed_sphere_grid') then wbeg = MPI_Wtime() if (trim(output_file(nnnn)) == 'netcdf_parallel') then call write_netcdf(wrt_int_state%wrtFB(nbdl),trim(filename), & .true., wrt_mpi_comm,wrt_int_state%mype, & grid_id,rc) else call ESMFproto_FieldBundleWrite(gridFB, filename=trim(filename), & convention="NetCDF", purpose="FV3", & status=ESMF_FILESTATUS_REPLACE, & state=stateGridFB, comps=compsGridFB,rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMFproto_FieldBundleWrite(wrt_int_state%wrtFB(nbdl), & filename=trim(filename), convention="NetCDF", & purpose="FV3", status=ESMF_FILESTATUS_OLD, & timeslice=step, state=optimize(nbdl)%state, & comps=optimize(nbdl)%comps, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return end if wend = MPI_Wtime() if (lprnt) then write(*,'(A15,A,F10.5,A,I4.2,A,I2.2,1X,A)')trim(output_file(nnnn)),' Write Time is ',wend-wbeg & ,' at Fcst ',NF_HOURS,':',NF_MINUTES endif else if (trim(output_grid(grid_id)) == 'gaussian_grid' .or. & trim(output_grid(grid_id)) == 'global_latlon') then wbeg = MPI_Wtime() call write_netcdf(wrt_int_state%wrtFB(nbdl),trim(filename), & use_parallel_netcdf, wrt_mpi_comm,wrt_int_state%mype, & grid_id,rc) wend = MPI_Wtime() if (lprnt) then write(*,'(A15,A,F10.5,A,I4.2,A,I2.2,1X,A)')trim(output_file(nnnn)),' Write Time is ',wend-wbeg & ,' at Fcst ',NF_HOURS,':',NF_MINUTES endif else if (trim(output_grid(grid_id)) == 'regional_latlon' .or. & trim(output_grid(grid_id)) == 'regional_latlon_moving' .or. & trim(output_grid(grid_id)) == 'rotated_latlon' .or. & trim(output_grid(grid_id)) == 'rotated_latlon_moving' .or. & trim(output_grid(grid_id)) == 'lambert_conformal') then !mask fields according to sfc pressure if( .not. lmask_fields ) then wbeg = MPI_Wtime() call mask_fields(file_bundle,rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return wend = MPI_Wtime() if (lprnt) then write(*,'(A,F10.5,A,I4.2,A,I2.2)')' mask_fields time is ',wend-wbeg endif endif if (nbits(grid_id) /= 0) then call ESMF_LogWrite("wrt_run: lossy compression is not supported for regional grids",ESMF_LOGMSG_ERROR,rc=RC) call ESMF_Finalize(endflag=ESMF_END_ABORT) end if wbeg = MPI_Wtime() call write_netcdf(wrt_int_state%wrtFB(nbdl),trim(filename), & use_parallel_netcdf, wrt_mpi_comm,wrt_int_state%mype, & grid_id,rc) wend = MPI_Wtime() if (lprnt) then write(*,'(A15,A,F10.5,A,I4.2,A,I2.2,1X,A)')trim(output_file(nnnn)),' Write Time is ',wend-wbeg & ,' at Fcst ',NF_HOURS,':',NF_MINUTES endif else ! unknown output_grid call ESMF_LogWrite("wrt_run: Unknown output_grid",ESMF_LOGMSG_ERROR,rc=RC) call ESMF_Finalize(endflag=ESMF_END_ABORT) endif enddo file_loop_all ! end output history endif ! !** write out log file ! if (mype == lead_write_task) then open(newunit=nolog,file='logf'//trim(cfhour),form='FORMATTED') write(nolog,100)nfhour,idate(1:6) 100 format(' completed fv3gfs fhour=',f10.3,2x,6(i4,2x)) close(nolog) endif ! !----------------------------------------------------------------------- ! call ESMF_VMBarrier(VM, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! write_run_tim = MPI_Wtime() - tbeg ! IF (lprnt) THEN WRITE(*,'(A,F10.5,A,I4.2,A,I2.2)')' total Write Time is ',write_run_tim & ,' at Fcst ',NF_HOURS,':',NF_MINUTES ENDIF ! !----------------------------------------------------------------------- ! END SUBROUTINE wrt_run ! !----------------------------------------------------------------------- !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& !----------------------------------------------------------------------- ! subroutine wrt_finalize(wrt_comp, imp_state_write, exp_state_write, clock, rc) ! !----------------------------------------------------------------------- !*** finalize the write gridded component. !----------------------------------------------------------------------- ! !*** HISTORY ! Feb 2017: J. Wang - deallocate for fv3 ! !----------------------------------------------------------------------- ! type(ESMF_GridComp) :: wrt_comp type(ESMF_State) :: imp_state_write, exp_state_write type(ESMF_Clock) :: clock integer,intent(out) :: rc ! !*** local variables ! integer :: stat type(write_wrap) :: wrap ! !----------------------------------------------------------------------- !*********************************************************************** !----------------------------------------------------------------------- ! rc=ESMF_SUCCESS ! !----------------------------------------------------------------------- !*** retrieve the write component's esmf internal state(used later for !*** post finalization) !----------------------------------------------------------------------- ! call ESMF_GridCompGetInternalState(wrt_comp, wrap, rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return deallocate(wrap%write_int_state,stat=stat) if (ESMF_LogFoundDeallocError(statusToCheck=stat, & msg="Deallocation of internal state memory failed.", & line=__LINE__, file=__FILE__)) return ! !----------------------------------------------------------------------- ! end subroutine wrt_finalize ! !----------------------------------------------------------------------- ! subroutine recover_fields(file_bundle,rc) type(ESMF_FieldBundle), intent(in) :: file_bundle integer, intent(out), optional :: rc ! real, parameter :: rdgas = 287.04, grav = 9.80 real, parameter :: stndrd_atmos_ps = 101325. real, parameter :: stndrd_atmos_lapse = 0.0065 integer i,j,k,ifld,fieldCount,nstt,nend,fieldDimCount,gridDimCount integer istart,iend,jstart,jend,kstart,kend logical uPresent, vPresent type(ESMF_Grid) fieldGrid type(ESMF_Field) ufield, vfield type(ESMF_TypeKind_Flag) typekind character(100) fieldName,uwindname,vwindname type(ESMF_Field), allocatable :: fcstField(:) real(ESMF_KIND_R8), dimension(:,:), pointer :: lon, lat real(ESMF_KIND_R8), dimension(:,:), pointer :: lonloc, latloc real(ESMF_KIND_R4), dimension(:,:), pointer :: pressfc real(ESMF_KIND_R4), dimension(:,:), pointer :: uwind2dr4,vwind2dr4 real(ESMF_KIND_R4), dimension(:,:,:), pointer :: uwind3dr4,vwind3dr4 real(ESMF_KIND_R4), dimension(:,:,:), pointer :: cart3dPtr2dr4 real(ESMF_KIND_R4), dimension(:,:,:,:), pointer :: cart3dPtr3dr4 real(ESMF_KIND_R8) :: coslon, sinlon, sinlat ! ! get filed count call ESMF_FieldBundleGet(file_bundle, fieldCount=fieldCount, & grid=fieldGrid, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! call ESMF_LogWrite("call recover field on wrt comp",ESMF_LOGMSG_INFO,rc=RC) call ESMF_GridGet(fieldgrid, dimCount=gridDimCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_LogWrite("call recover field get coord 1",ESMF_LOGMSG_INFO,rc=RC) call ESMF_GridGetCoord(fieldgrid, coordDim=1, farrayPtr=lon, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return allocate(lonloc(lbound(lon,1):ubound(lon,1),lbound(lon,2):ubound(lon,2))) istart = lbound(lon,1) iend = ubound(lon,1) jstart = lbound(lon,2) jend = ubound(lon,2) !$omp parallel do default(none) shared(lon,lonloc,jstart,jend,istart,iend) & !$omp private(i,j) do j=jstart,jend do i=istart,iend lonloc(i,j) = lon(i,j) * pi/180. enddo enddo call ESMF_LogWrite("call recover field get coord 2",ESMF_LOGMSG_INFO,rc=RC) call ESMF_GridGetCoord(fieldgrid, coordDim=2, farrayPtr=lat, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return allocate(latloc(lbound(lat,1):ubound(lat,1),lbound(lat,2):ubound(lat,2))) istart = lbound(lat,1) iend = ubound(lat,1) jstart = lbound(lat,2) jend = ubound(lat,2) !$omp parallel do default(none) shared(lat,latloc,jstart,jend,istart,iend) & !$omp private(i,j) do j=jstart,jend do i=istart,iend latloc(i,j) = lat(i,j) * pi/180.d0 enddo enddo ! allocate(fcstField(fieldCount)) call ESMF_LogWrite("call recover field get fcstField",ESMF_LOGMSG_INFO,rc=RC) call ESMF_FieldBundleGet(file_bundle, fieldList=fcstField, itemorderflag=ESMF_ITEMORDER_ADDORDER, rc=rc) ! do ifld=1,fieldCount call ESMF_LogWrite("call recover field get fieldname, type dimcount",ESMF_LOGMSG_INFO,rc=RC) call ESMF_FieldGet(fcstField(ifld),name=fieldName,typekind=typekind,dimCount=fieldDimCount, rc=rc) ! convert back wind if(index(trim(fieldName),"vector")>0) then nstt = index(trim(fieldName),"wind")+4 nend = index(trim(fieldName),"vector")-1 if( nend>nstt ) then uwindname = 'u'//fieldName(nstt+1:nend) vwindname = 'v'//fieldName(nstt+1:nend) else uwindname = 'ugrd' vwindname = 'vgrd' endif ! print *,'in get 3D vector wind, uwindname=',trim(uwindname),' v=', trim(vwindname),' fieldname=',trim(fieldname) ! get u , v wind call ESMF_LogWrite("call recover field get u, v field",ESMF_LOGMSG_INFO,rc=RC) call ESMF_FieldBundleGet(file_bundle,trim(uwindname),field=ufield,isPresent=uPresent,rc=rc) call ESMF_FieldBundleGet(file_bundle,trim(vwindname),field=vfield,isPresent=vPresent,rc=rc) if(.not. uPresent .or. .not.vPresent) then rc=990 print *,' ERROR ,the local wind is not present! rc=', rc exit endif ! get field data if ( typekind == ESMF_TYPEKIND_R4 ) then if( fieldDimCount > gridDimCount+1 ) then call ESMF_LogWrite("call recover field get 3d card wind farray",ESMF_LOGMSG_INFO,rc=RC) call ESMF_FieldGet(fcstField(ifld), localDe=0, farrayPtr=cart3dPtr3dr4, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if( ubound(cart3dPtr3dr4,1)-lbound(cart3dPtr3dr4,1)+1/=3) then rc=991 print *,'ERROR, 3D the vector dimension /= 3, rc=',rc exit endif iend = ubound(cart3dPtr3dr4,2) istart = lbound(cart3dPtr3dr4,2) iend = ubound(cart3dPtr3dr4,2) jstart = lbound(cart3dPtr3dr4,3) jend = ubound(cart3dPtr3dr4,3) kstart = lbound(cart3dPtr3dr4,4) kend = ubound(cart3dPtr3dr4,4) call ESMF_FieldGet(ufield, localDe=0, farrayPtr=uwind3dr4,rc=rc) call ESMF_FieldGet(vfield, localDe=0, farrayPtr=vwind3dr4,rc=rc) ! update u , v wind !$omp parallel do default(shared) private(i,j,k,coslon,sinlon,sinlat) do k=kstart,kend !$omp parallel do default(none) shared(uwind3dr4,vwind3dr4,lonloc,latloc,cart3dPtr3dr4,jstart,jend,istart,iend,k) & !$omp private(i,j,coslon,sinlon,sinlat) do j=jstart, jend do i=istart, iend coslon = cos(lonloc(i,j)) sinlon = sin(lonloc(i,j)) sinlat = sin(latloc(i,j)) uwind3dr4(i,j,k) = cart3dPtr3dr4(1,i,j,k) * coslon & + cart3dPtr3dr4(2,i,j,k) * sinlon vwind3dr4(i,j,k) =-cart3dPtr3dr4(1,i,j,k) * sinlat*sinlon & + cart3dPtr3dr4(2,i,j,k) * sinlat*coslon & + cart3dPtr3dr4(3,i,j,k) * cos(latloc(i,j)) enddo enddo enddo else call ESMF_FieldGet(fcstField(ifld), localDe=0, farrayPtr=cart3dPtr2dr4, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if( ubound(cart3dPtr2dr4,1)-lbound(cart3dPtr2dr4,1)+1 /= 3) then rc=991 write(0,*) 'ERROR, 2D the vector dimension /= 3, rc=',rc call ESMF_Finalize(endflag=ESMF_END_ABORT) endif istart = lbound(cart3dPtr2dr4,2) iend = ubound(cart3dPtr2dr4,2) jstart = lbound(cart3dPtr2dr4,3) jend = ubound(cart3dPtr2dr4,3) call ESMF_FieldGet(ufield, localDe=0, farrayPtr=uwind2dr4,rc=rc) call ESMF_FieldGet(vfield, localDe=0, farrayPtr=vwind2dr4,rc=rc) ! update u , v wind !$omp parallel do default(none) shared(uwind2dr4,vwind2dr4,lonloc,latloc,cart3dPtr2dr4,jstart,jend,istart,iend) & !$omp private(i,j,k,coslon,sinlon,sinlat) do j=jstart, jend do i=istart, iend coslon = cos(lonloc(i,j)) sinlon = sin(lonloc(i,j)) sinlat = sin(latloc(i,j)) uwind2dr4(i,j) = cart3dPtr2dr4(1,i,j) * coslon & + cart3dPtr2dr4(2,i,j) * sinlon vwind2dr4(i,j) =-cart3dPtr2dr4(1,i,j) * sinlat*sinlon & + cart3dPtr2dr4(2,i,j) * sinlat*coslon & + cart3dPtr2dr4(3,i,j) * cos(latloc(i,j)) enddo enddo endif endif ! print *,'in 3DCartesian2wind, uwindname=', trim(uwindname),'uPresent =',uPresent, 'vPresent=',vPresent,'fieldDimCount=', & ! fieldDimCount,'gridDimCount=',gridDimCount ! convert back surface pressure else if(index(trim(fieldName),"pressfc")>0) then call ESMF_FieldGet(fcstField(ifld),localDe=0, farrayPtr=pressfc, rc=rc) istart = lbound(pressfc,1) iend = ubound(pressfc,1) jstart = lbound(pressfc,2) jend = ubound(pressfc,2) !$omp parallel do default(none) shared(pressfc,jstart,jend,istart,iend) private(i,j) do j=jstart, jend do i=istart, iend pressfc(i,j) = pressfc(i,j)**(grav/(rdgas*stndrd_atmos_lapse))*stndrd_atmos_ps enddo enddo endif enddo ! deallocate(fcstField) rc = 0 end subroutine recover_fields ! !----------------------------------------------------------------------- ! subroutine mask_fields(file_bundle,rc) type(ESMF_FieldBundle), intent(in) :: file_bundle integer, intent(out), optional :: rc ! integer i,j,k,ifld,fieldCount,fieldDimCount,gridDimCount integer istart,iend,jstart,jend,kstart,kend type(ESMF_Grid) fieldGrid type(ESMF_TypeKind_Flag) typekind type(ESMF_TypeKind_Flag) attTypeKind character(len=ESMF_MAXSTR) fieldName type(ESMF_Field), allocatable :: fcstField(:) real(ESMF_KIND_R4), dimension(:,:), pointer :: var2dPtr2dr4 real(ESMF_KIND_R4), dimension(:,:,:), pointer :: var3dPtr3dr4 real(ESMF_KIND_R4), dimension(:,:,:), pointer :: vect3dPtr2dr4 real(ESMF_KIND_R4), dimension(:,:,:,:), pointer :: vect4dPtr3dr4 real(ESMF_KIND_R4), dimension(:,:), allocatable :: maskwrt logical :: mvispresent=.false. real(ESMF_KIND_R4) :: missing_value_r4=-1.e+10 real(ESMF_KIND_R8) :: missing_value_r8=9.99e20 character(len=ESMF_MAXSTR) :: msg call ESMF_LogWrite("call mask field on wrt comp",ESMF_LOGMSG_INFO,rc=RC) ! get fieldCount call ESMF_FieldBundleGet(file_bundle, fieldCount=fieldCount, & grid=fieldGrid, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__)) return ! bail out ! get gridDimCount call ESMF_GridGet(fieldgrid, dimCount=gridDimCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__)) return ! bail out allocate(fcstField(fieldCount)) call ESMF_LogWrite("call mask field get fcstField",ESMF_LOGMSG_INFO,rc=RC) call ESMF_FieldBundleGet(file_bundle, fieldList=fcstField, itemorderflag=ESMF_ITEMORDER_ADDORDER, rc=rc) ! generate the maskwrt according to surface pressure do ifld=1,fieldCount !call ESMF_LogWrite("call mask field get fieldname, type dimcount",ESMF_LOGMSG_INFO,rc=RC) call ESMF_FieldGet(fcstField(ifld),name=fieldName,typekind=typekind,dimCount=fieldDimCount, rc=rc) !write(msg,*) 'fieldName,typekind,fieldDimCount=',trim(fieldName),typekind,fieldDimCount !call ESMF_LogWrite("call mask field: "//trim(msg),ESMF_LOGMSG_INFO,rc=RC) if (.not. allocated(maskwrt)) then if ( typekind == ESMF_TYPEKIND_R4 .and. fieldDimCount == gridDimCount) then call ESMF_FieldGet(fcstField(ifld),localDe=0, farrayPtr=var2dPtr2dr4, rc=rc) istart = lbound(var2dPtr2dr4,1) iend = ubound(var2dPtr2dr4,1) jstart = lbound(var2dPtr2dr4,2) jend = ubound(var2dPtr2dr4,2) allocate(maskwrt(istart:iend,jstart:jend)) maskwrt(istart:iend,jstart:jend)=1.0 endif endif if(index(trim(fieldName),"pressfc")>0) then call ESMF_FieldGet(fcstField(ifld),localDe=0, farrayPtr=var2dPtr2dr4, rc=rc) istart = lbound(var2dPtr2dr4,1) iend = ubound(var2dPtr2dr4,1) jstart = lbound(var2dPtr2dr4,2) jend = ubound(var2dPtr2dr4,2) if (.not. allocated(maskwrt)) then allocate(maskwrt(istart:iend,jstart:jend)) maskwrt(istart:iend,jstart:jend)=1.0 endif !$omp parallel do default(shared) private(i,j) do j=jstart, jend do i=istart, iend if(abs(var2dPtr2dr4(i,j)-0.) < 1.0e-6) maskwrt(i,j)=0. enddo enddo call ESMF_LogWrite("call mask field pressfc found, maskwrt generated",ESMF_LOGMSG_INFO,rc=RC) exit endif enddo ! loop to mask all fields according to maskwrt do ifld=1,fieldCount !call ESMF_LogWrite("call mask field get fieldname, type dimcount",ESMF_LOGMSG_INFO,rc=RC) call ESMF_FieldGet(fcstField(ifld),name=fieldName,typekind=typekind,dimCount=fieldDimCount, rc=rc) !write(msg,*) 'fieldName,typekind,fieldDimCount=',trim(fieldName),typekind,fieldDimCount !call ESMF_LogWrite("call mask field: "//trim(msg),ESMF_LOGMSG_INFO,rc=RC) ! For vector fields if(index(trim(fieldName),"vector")>0) then ! Only work on ESMF_TYPEKIND_R4 fields for now if ( typekind == ESMF_TYPEKIND_R4 ) then ! 3-d vector fields with 4-d arrays if( fieldDimCount > gridDimCount+1 ) then !call ESMF_LogWrite("call mask field get vector 3d farray",ESMF_LOGMSG_INFO,rc=RC) call ESMF_FieldGet(fcstField(ifld), localDe=0, farrayPtr=vect4dPtr3dr4, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__)) return ! bail out if( ubound(vect4dPtr3dr4,1)-lbound(vect4dPtr3dr4,1)+1/=3 ) then rc=991 write(0,*) 'ERROR, 3D the vector dimension /= 3, rc=',rc call ESMF_Finalize(endflag=ESMF_END_ABORT) endif ! Get the _FillValue from the field attribute if exists call ESMF_AttributeGet(fcstField(ifld), convention="NetCDF", purpose="FV3", & name="_FillValue", typekind=attTypeKind, isPresent=mvispresent, rc=rc) !write(msg,*) 'fieldName,attTypeKind,isPresent=',trim(fieldName),attTypeKind,mvispresent !call ESMF_LogWrite("call mask field: "//trim(msg),ESMF_LOGMSG_INFO,rc=RC) if ( mvispresent ) then if (attTypeKind==ESMF_TYPEKIND_R4) then call ESMF_AttributeGet(fcstField(ifld), convention="NetCDF", purpose="FV3", & name="_FillValue", value=missing_value_r4, isPresent=mvispresent, rc=rc) !write(msg,*) 'fieldName,_FillValue,isPresent=',trim(fieldName),missing_value_r4,mvispresent !call ESMF_LogWrite("call mask field: "//trim(msg),ESMF_LOGMSG_INFO,rc=RC) else if (attTypeKind==ESMF_TYPEKIND_R8) then call ESMF_AttributeGet(fcstField(ifld), convention="NetCDF", purpose="FV3", & name="_FillValue", value=missing_value_r8, isPresent=mvispresent, rc=rc) !write(msg,*) 'fieldName,_FillValue,isPresent=',trim(fieldName),missing_value_r8,mvispresent !call ESMF_LogWrite("call mask field: "//trim(msg),ESMF_LOGMSG_INFO,rc=RC) endif istart = lbound(vect4dPtr3dr4,2) iend = ubound(vect4dPtr3dr4,2) jstart = lbound(vect4dPtr3dr4,3) jend = ubound(vect4dPtr3dr4,3) kstart = lbound(vect4dPtr3dr4,4) kend = ubound(vect4dPtr3dr4,4) !$omp parallel do default(shared) private(i,j,k) do k=kstart,kend do j=jstart, jend do i=istart, iend if (maskwrt(i,j)<1.0 .and. attTypeKind==ESMF_TYPEKIND_R4) vect4dPtr3dr4(:,i,j,k)=missing_value_r4 if (maskwrt(i,j)<1.0 .and. attTypeKind==ESMF_TYPEKIND_R8) vect4dPtr3dr4(:,i,j,k)=missing_value_r8 enddo enddo enddo endif !mvispresent ! 2-d vector fields with 3-d arrays else call ESMF_FieldGet(fcstField(ifld), localDe=0, farrayPtr=vect3dPtr2dr4, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__)) return ! bail out if( ubound(vect3dPtr2dr4,1)-lbound(vect3dPtr2dr4,1)+1 /= 3 ) then rc=991 write(0,*) 'ERROR, 2D the vector dimension /= 3, rc=',rc call ESMF_Finalize(endflag=ESMF_END_ABORT) endif ! Get the _FillValue from the field attribute if exists call ESMF_AttributeGet(fcstField(ifld), convention="NetCDF", purpose="FV3", & name="_FillValue", typekind=attTypeKind, isPresent=mvispresent, rc=rc) !write(msg,*) 'fieldName,attTypeKind,isPresent=',trim(fieldName),attTypeKind,mvispresent !call ESMF_LogWrite("call mask field: "//trim(msg),ESMF_LOGMSG_INFO,rc=RC) if ( mvispresent ) then if (attTypeKind==ESMF_TYPEKIND_R4) then call ESMF_AttributeGet(fcstField(ifld), convention="NetCDF", purpose="FV3", & name="_FillValue", value=missing_value_r4, isPresent=mvispresent, rc=rc) !write(msg,*) 'fieldName,_FillValue,isPresent=',trim(fieldName),missing_value_r4,mvispresent !call ESMF_LogWrite("call mask field: "//trim(msg),ESMF_LOGMSG_INFO,rc=RC) else if (attTypeKind==ESMF_TYPEKIND_R8) then call ESMF_AttributeGet(fcstField(ifld), convention="NetCDF", purpose="FV3", & name="_FillValue", value=missing_value_r8, isPresent=mvispresent, rc=rc) !write(msg,*) 'fieldName,_FillValue,isPresent=',trim(fieldName),missing_value_r8,mvispresent !call ESMF_LogWrite("call mask field: "//trim(msg),ESMF_LOGMSG_INFO,rc=RC) endif istart = lbound(vect3dPtr2dr4,2) iend = ubound(vect3dPtr2dr4,2) jstart = lbound(vect3dPtr2dr4,3) jend = ubound(vect3dPtr2dr4,3) !$omp parallel do default(shared) private(i,j) do j=jstart, jend do i=istart, iend if (maskwrt(i,j)<1.0 .and. attTypeKind==ESMF_TYPEKIND_R4) vect3dPtr2dr4(:,i,j)=missing_value_r4 if (maskwrt(i,j)<1.0 .and. attTypeKind==ESMF_TYPEKIND_R8) vect3dPtr2dr4(:,i,j)=missing_value_r8 enddo enddo endif !mvispresent endif endif ! For non-vector fields else ! Only work on ESMF_TYPEKIND_R4 fields for now if ( typekind == ESMF_TYPEKIND_R4 ) then ! 2-d fields if(fieldDimCount == gridDimCount) then call ESMF_FieldGet(fcstField(ifld),localDe=0, farrayPtr=var2dPtr2dr4, rc=rc) ! Get the _FillValue from the field attribute if exists call ESMF_AttributeGet(fcstField(ifld), convention="NetCDF", purpose="FV3", & name="_FillValue", typekind=attTypeKind, isPresent=mvispresent, rc=rc) !write(msg,*) 'fieldName,attTypeKind,isPresent=',trim(fieldName),attTypeKind,mvispresent !call ESMF_LogWrite("call mask field: "//trim(msg),ESMF_LOGMSG_INFO,rc=RC) if ( mvispresent ) then if (attTypeKind==ESMF_TYPEKIND_R4) then call ESMF_AttributeGet(fcstField(ifld), convention="NetCDF", purpose="FV3", & name="_FillValue", value=missing_value_r4, isPresent=mvispresent, rc=rc) !write(msg,*) 'fieldName,_FillValue,isPresent=',trim(fieldName),missing_value_r4,mvispresent !call ESMF_LogWrite("call mask field: "//trim(msg),ESMF_LOGMSG_INFO,rc=RC) else if (attTypeKind==ESMF_TYPEKIND_R8) then call ESMF_AttributeGet(fcstField(ifld), convention="NetCDF", purpose="FV3", & name="_FillValue", value=missing_value_r8, isPresent=mvispresent, rc=rc) !write(msg,*) 'fieldName,_FillValue,isPresent=',trim(fieldName),missing_value_r8,mvispresent !call ESMF_LogWrite("call mask field: "//trim(msg),ESMF_LOGMSG_INFO,rc=RC) endif istart = lbound(var2dPtr2dr4,1) iend = ubound(var2dPtr2dr4,1) jstart = lbound(var2dPtr2dr4,2) jend = ubound(var2dPtr2dr4,2) !$omp parallel do default(shared) private(i,j) do j=jstart, jend do i=istart, iend if (maskwrt(i,j)<1.0 .and. attTypeKind==ESMF_TYPEKIND_R4) var2dPtr2dr4(i,j)=missing_value_r4 if (maskwrt(i,j)<1.0 .and. attTypeKind==ESMF_TYPEKIND_R8) var2dPtr2dr4(i,j)=missing_value_r8 enddo enddo endif !mvispresent ! 3-d fields else if(fieldDimCount == gridDimCount+1) then call ESMF_FieldGet(fcstField(ifld),localDe=0, farrayPtr=var3dPtr3dr4, rc=rc) ! Get the _FillValue from the field attribute if exists call ESMF_AttributeGet(fcstField(ifld), convention="NetCDF", purpose="FV3", & name="_FillValue", typekind=attTypeKind, isPresent=mvispresent, rc=rc) !write(msg,*) 'fieldName,attTypeKind,isPresent=',trim(fieldName),attTypeKind,mvispresent !call ESMF_LogWrite("call mask field: "//trim(msg),ESMF_LOGMSG_INFO,rc=RC) if ( mvispresent ) then if (attTypeKind==ESMF_TYPEKIND_R4) then call ESMF_AttributeGet(fcstField(ifld), convention="NetCDF", purpose="FV3", & name="_FillValue", value=missing_value_r4, isPresent=mvispresent, rc=rc) !write(msg,*) 'fieldName,_FillValue,isPresent=',trim(fieldName),missing_value_r4,mvispresent !call ESMF_LogWrite("call mask field: "//trim(msg),ESMF_LOGMSG_INFO,rc=RC) else if (attTypeKind==ESMF_TYPEKIND_R8) then call ESMF_AttributeGet(fcstField(ifld), convention="NetCDF", purpose="FV3", & name="_FillValue", value=missing_value_r8, isPresent=mvispresent, rc=rc) !write(msg,*) 'fieldName,_FillValue,isPresent=',trim(fieldName),missing_value_r8,mvispresent !call ESMF_LogWrite("call mask field: "//trim(msg),ESMF_LOGMSG_INFO,rc=RC) endif istart = lbound(var3dPtr3dr4,1) iend = ubound(var3dPtr3dr4,1) jstart = lbound(var3dPtr3dr4,2) jend = ubound(var3dPtr3dr4,2) kstart = lbound(var3dPtr3dr4,3) kend = ubound(var3dPtr3dr4,3) !$omp parallel do default(shared) private(i,j,k) do k=kstart,kend do j=jstart, jend do i=istart, iend if (maskwrt(i,j)<1.0 .and. attTypeKind==ESMF_TYPEKIND_R4) var3dPtr3dr4(i,j,k)=missing_value_r4 if (maskwrt(i,j)<1.0 .and. attTypeKind==ESMF_TYPEKIND_R8) var3dPtr3dr4(i,j,k)=missing_value_r8 enddo enddo enddo endif !mvispresent endif endif endif enddo ! deallocate(maskwrt) deallocate(fcstField) rc = 0 end subroutine mask_fields ! !----------------------------------------------------------------------- ! subroutine ESMFproto_FieldBundleWrite(fieldbundle, fileName, & convention, purpose, status, timeslice, state, comps, rc) type(ESMF_FieldBundle), intent(in) :: fieldbundle character(*), intent(in) :: fileName character(*), intent(in), optional :: convention character(*), intent(in), optional :: purpose type(ESMF_FileStatus_Flag), intent(in), optional :: status integer, intent(in), optional :: timeslice type(ESMF_State), intent(inout), optional :: state type(ESMF_GridComp), allocatable, intent(inout), optional :: comps(:) integer, intent(out), optional :: rc ! Prototype multi-tile implementation for FieldBundleWrite(). ! Produces as many output files as there are tiles. The naming of the ! output files is such that the string in the fileName argument is used ! as the basis. If fileName ends with ".nc", then this suffix is replaced ! by ".tileN.nc", where "N" is the tile number. If fileName does not ! end in ".nc", then ".tileN.nc" will simply be appended. ! ! Restrictions: ! - All Fields in the FieldBundle must have the same tileCount integer :: i, j, ind integer :: fieldCount, tileCount, itemCount type(ESMF_Field), allocatable :: fieldList(:), tileFieldList(:) type(ESMF_Grid) :: grid type(ESMF_Array) :: array type(ESMF_DistGrid) :: distgrid type(ESMF_DELayout) :: delayout type(ESMF_FieldBundle) :: wrtTileFB type(ESMF_FieldBundle), allocatable :: wrtTileFBList(:) character(len=80), allocatable :: itemNameList(:) logical :: stateIsEmpty character(len=80) :: tileFileName character(len=80) :: statusStr integer, allocatable :: petList(:) character(1024) :: msgString type(ESMF_State), allocatable :: ioState(:) integer :: timesliceOpt integer :: urc if (present(rc)) rc = ESMF_SUCCESS call ESMF_VMLogMemInfo("Entering ESMFproto_FieldBundleWrite",rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! query number of fields in fieldbundle call ESMF_FieldBundleGet(fieldbundle, fieldCount=fieldCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! early successful exit if there are no fields present if (fieldCount==0) return ! obtain list of fields in the fieldbundle allocate(fieldList(fieldCount), tileFieldList(fieldCount)) call ESMF_FieldBundleGet(fieldbundle, fieldList=fieldList, & itemorderflag=ESMF_ITEMORDER_ADDORDER, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! determine tileCount by looking at first field call ESMF_FieldGet(fieldList(1), array=array, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_ArrayGet(array, tileCount=tileCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! deal with optional state argument stateIsEmpty = .true. if (present(state)) then if (.not.ESMF_StateIsCreated(state, rc=rc)) then state = ESMF_StateCreate(rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return endif call ESMF_StateGet(state, itemCount=itemCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if (itemCount /= 0) then stateIsEmpty = .false. if (itemCount /= tileCount) then call ESMF_LogSetError(ESMF_RC_ARG_BAD, & msg="Number of items in state must match number of tiles.", & line=__LINE__, file=__FILE__, rcToReturn=rc) return endif allocate(itemNameList(itemCount),wrtTileFBList(itemCount)) call ESMF_StateGet(state, itemNameList=itemNameList, & itemorderflag=ESMF_ITEMORDER_ADDORDER, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return do i=1, itemCount call ESMF_StateGet(state, itemName=itemNameList(i), & fieldbundle=wrtTileFBList(i), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return enddo endif endif ! loop over all the tiles and construct a tile specific fieldbundle call ESMF_LogWrite("In ESMFproto_FieldBundleWrite() before tileCount-loop",& ESMF_LOGMSG_INFO, rc=rc) !TODO: remove this once comps is hidden within state if (present(comps)) then allocate(ioState(tileCount)) if (.not.allocated(comps)) then ! first-write allocate(comps(tileCount)) do i=1, tileCount call ESMF_LogWrite("In ESMFproto_FieldBundleWrite() before "// & "ESMFproto_FieldMakeSingleTile() w/ petList", & ESMF_LOGMSG_INFO, rc=rc) do j=1, fieldCount ! access only tile specific part of field call ESMFproto_FieldMakeSingleTile(fieldList(j), tile=i, & tileField=tileFieldList(j), petList=petList, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return enddo ! write(msgString, *) petList ! call ESMF_LogWrite("In ESMFproto_FieldBundleWrite() after "// & ! "ESMFproto_FieldMakeSingleTile(), petList:"//trim(msgString), & ! ESMF_LOGMSG_INFO, rc=rc) ! create component to handle this tile I/O call ESMF_LogWrite("In ESMFproto_FieldBundleWrite() before "// & "tile-component creation", ESMF_LOGMSG_INFO, rc=rc) comps(i) = ESMF_GridCompCreate(petList=petList, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! convention call ESMF_AttributeSet(comps(i), name="convention", & value=convention, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! purpose call ESMF_AttributeSet(comps(i), name="purpose", & value=purpose, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! timeslice timesliceOpt = -1 ! init if (present(timeslice)) timesliceOpt = timeslice call ESMF_AttributeSet(comps(i), name="timeslice", & value=timesliceOpt, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! status statusStr = "ESMF_FILESTATUS_UNKNOWN" ! default if (present(status)) then if (status==ESMF_FILESTATUS_UNKNOWN) then statusStr="ESMF_FILESTATUS_UNKNOWN" ! default else if (status==ESMF_FILESTATUS_NEW) then statusStr="ESMF_FILESTATUS_NEW" ! default else if (status==ESMF_FILESTATUS_OLD) then statusStr="ESMF_FILESTATUS_OLD" ! default else if (status==ESMF_FILESTATUS_REPLACE) then statusStr="ESMF_FILESTATUS_REPLACE" ! default endif endif call ESMF_AttributeSet(comps(i), name="status", & value=statusStr, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_GridCompSetServices(comps(i), ioCompSS, userRc=urc, rc=rc) if (ESMF_LogFoundError(rcToCheck=urc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_LogWrite("In ESMFproto_FieldBundleWrite() after "// & "tile-component creation", ESMF_LOGMSG_INFO, rc=rc) enddo endif endif do i=1, tileCount if (stateIsEmpty) then ! loop over all the fields and add tile specific part to fieldbundle call ESMF_LogWrite("In ESMFproto_FieldBundleWrite() before "// & "ESMFproto_FieldMakeSingleTile()", ESMF_LOGMSG_INFO, rc=rc) do j=1, fieldCount ! access only tile specific part of field call ESMFproto_FieldMakeSingleTile(fieldList(j), tile=i, & tileField=tileFieldList(j), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return enddo call ESMF_LogWrite("In ESMFproto_FieldBundleWrite() after "// & "ESMFproto_FieldMakeSingleTile()", ESMF_LOGMSG_INFO, rc=rc) ! create tile specific fieldbundle wrtTileFB = ESMF_FieldBundleCreate(fieldList=tileFieldList, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! ensure global attributes on the fieldbundle are passed on by reference call ESMF_AttributeCopy(fieldbundle, wrtTileFB, & attcopy=ESMF_ATTCOPY_REFERENCE, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! store this fieldbundle in state if present if (present(state)) then call ESMF_StateAdd(state, fieldbundleList=(/wrtTileFB/), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return endif else ! state brought in existing fieldbundles wrtTileFB = wrtTileFBList(i) endif ! write out the tile specific fieldbundle if(tileCount>1) then ind=min(index(trim(fileName),".nc",.true.)-1,len_trim(fileName)) write(tileFileName, "(A,A,I1,A)") fileName(1:ind), ".tile", i, ".nc" else tileFileName=trim(fileName) endif if (present(comps)) then ioState(i) = ESMF_StateCreate(rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_StateAdd(ioState(i), (/wrtTileFB/), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeSet(comps(i), name="tileFileName", value=tileFileName, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_GridCompRun(comps(i), importState=ioState(i), userRc=urc, rc=rc) if (ESMF_LogFoundError(rcToCheck=urc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_StateDestroy(ioState(i), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return else call ESMF_LogWrite("In ESMFproto_FieldBundleWrite() before "// & "ESMF_FieldBundleWrite(): "//trim(tileFileName), & ESMF_LOGMSG_INFO, rc=rc) call ESMF_FieldBundleWrite(fieldbundle=wrtTileFB, fileName=tileFileName, & convention=convention, purpose=purpose, status=status, & timeslice=timeslice, overwrite=.true., rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_LogWrite("In ESMFproto_FieldBundleWrite() after "// & "ESMF_FieldBundleWrite()", ESMF_LOGMSG_INFO, rc=rc) endif if (.not.present(state)) then ! local garbage collection of fields do j=1, fieldCount call ESMF_FieldGet(tileFieldList(j), array=array, grid=grid, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_FieldDestroy(tileFieldList(j), noGarbage=.true., rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_GridDestroy(grid, noGarbage=.true., rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_ArrayGet(array, distgrid=distgrid, delayout=delayout, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_ArrayDestroy(array, noGarbage=.true., rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_DistGridDestroy(distgrid, noGarbage=.true., rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_DELayoutDestroy(delayout, noGarbage=.true., rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return enddo ! destroy tile specific fieldbundle call ESMF_FieldBundleDestroy(wrtTileFB, noGarbage=.true., rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return endif enddo if (present(comps)) then deallocate(ioState) endif call ESMF_LogWrite("In ESMFproto_FieldBundleWrite() after tileCount-loop",& ESMF_LOGMSG_INFO, rc=rc) ! deallocate temporary lists deallocate(fieldList, tileFieldList) call ESMF_VMLogMemInfo("Exiting ESMFproto_FieldBundleWrite",rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return end subroutine ESMFproto_FieldBundleWrite !----------------------------------------------------------------------------- subroutine ioCompSS(comp, rc) type(ESMF_GridComp) :: comp integer, intent(out) :: rc rc = ESMF_SUCCESS call ESMF_GridCompSetEntryPoint(comp, ESMF_METHOD_RUN, & userRoutine=ioCompRun, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return end subroutine !----------------------------------------------------------------------------- subroutine ioCompRun(comp, importState, exportState, clock, rc) use netcdf type(ESMF_GridComp) :: comp type(ESMF_State) :: importState, exportState type(ESMF_Clock) :: clock integer, intent(out) :: rc type(ESMF_FieldBundle) :: wrtTileFB character(len=80) :: tileFileName character(len=80) :: convention character(len=80) :: purpose integer :: timeslice character(len=80) :: statusStr type(ESMF_FileStatus_Flag) :: status character(len=80) :: itemNameList(1) integer :: localPet, i, j, k, ind type(ESMF_Grid) :: grid real(ESMF_KIND_R4), allocatable :: valueListr4(:) real(ESMF_KIND_R8), allocatable :: valueListr8(:) integer :: valueCount, fieldCount, udimCount character(80), allocatable :: udimList(:) integer :: ncerr, ncid, dimid, varid type(ESMF_Field), allocatable :: fieldList(:) type(ESMF_Field) :: field logical :: isPresent real(ESMF_KIND_R8) :: time integer :: itemCount, attCount character(len=80), allocatable :: attNameList(:) character(len=80) :: attName type(ESMF_TypeKind_Flag) :: typekind character(len=80) :: valueS integer :: valueI4 real(ESMF_KIND_R4) :: valueR4 real(ESMF_KIND_R8) :: valueR8 logical :: thereAreVerticals integer :: ch_dimid, timeiso_varid character(len=ESMF_MAXSTR) :: time_iso rc = ESMF_SUCCESS ! Access the FieldBundle call ESMF_StateGet(importState, itemNameList=itemNameList, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_StateGet(importState, itemName=itemNameList(1), & fieldBundle=wrtTileFB, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! Access attributes on the component and use as parameters for Write() call ESMF_AttributeGet(comp, name="tileFileName", value=tileFileName, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeGet(comp, name="convention", value=convention, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeGet(comp, name="purpose", value=purpose, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeGet(comp, name="timeslice", value=timeslice, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeGet(comp, name="status", value=statusStr, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if (trim(statusStr) == "ESMF_FILESTATUS_UNKNOWN") then status = ESMF_FILESTATUS_UNKNOWN else if (trim(statusStr) == "ESMF_FILESTATUS_NEW") then status = ESMF_FILESTATUS_NEW else if (trim(statusStr) == "ESMF_FILESTATUS_OLD") then status=ESMF_FILESTATUS_OLD else if (trim(statusStr) == "ESMF_FILESTATUS_REPLACE") then status = ESMF_FILESTATUS_REPLACE endif call ESMF_LogWrite("In ioCompRun() before writing to: "// & trim(tileFileName), ESMF_LOGMSG_INFO, rc=rc) if (status == ESMF_FILESTATUS_OLD) then ! This writes the vectical coordinates and the time dimension into the ! file. Doing this before the large data sets are written, assuming that ! the first time coming into ioCompRun() with this tileFileName, only ! the grid info is written. Second time in, with ESMF_FILESTATUS_OLD, ! the large data sets are written. That is when vertical and time info ! is also written. ! Hoping for better performance because file exists, but is still small. call ESMF_GridCompGet(comp, localPet=localPet, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if (localPet==0) then ! do this work only on the root pet call ESMF_FieldBundleGet(wrtTileFB, grid=grid, fieldCount=fieldCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return allocate(fieldList(fieldCount)) call ESMF_FieldBundleGet(wrtTileFB, fieldList=fieldList, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! open this tile's NetCDF file ncerr = nf90_open(tileFileName, NF90_WRITE, ncid=ncid) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return ! loop over all the fields in the bundle and handle their vectical dims thereAreVerticals = .false. do i=1, fieldCount field = fieldList(i) call ESMF_AttributeGetAttPack(field, convention="NetCDF", purpose="FV3", & isPresent=isPresent, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if (.not.isPresent) cycle ! field does not have the AttPack call ESMF_AttributeGet(field, convention="NetCDF", purpose="FV3", & name="ESMF:ungridded_dim_labels", isPresent=isPresent, & itemCount=udimCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if (udimCount==0 .or. .not.isPresent) cycle ! nothing there to do thereAreVerticals = .true. allocate(udimList(udimCount)) call ESMF_AttributeGet(field, convention="NetCDF", purpose="FV3", & name="ESMF:ungridded_dim_labels", valueList=udimList, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! loop over all ungridded dimension labels do k=1, udimCount call write_out_ungridded_dim_atts(dimLabel=trim(udimList(k)), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return enddo deallocate(udimList) enddo ! fieldCount deallocate(fieldList) if (thereAreVerticals) then ! see if the vertical_dim_labels attribute exists on the grid, and ! if so access it and write out vecticals accordingly call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & name="vertical_dim_labels", isPresent=isPresent, & itemCount=udimCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if (isPresent .and. (udimCount>0) ) then allocate(udimList(udimCount)) call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & name="vertical_dim_labels", valueList=udimList, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! loop over all ungridded dimension labels do k=1, udimCount call write_out_ungridded_dim_atts(dimLabel=trim(udimList(k)), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return enddo deallocate(udimList) endif endif ! inquire if NetCDF file already contains the "time" variable ncerr = nf90_inq_varid(ncid, "time", varid=varid) if (ncerr /= NF90_NOERR) then ! the variable does not exist in the NetCDF file yet -> add it ! access the "time" attribute on the grid call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & name="time", value=time, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & name="time_iso", value=time_iso, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ncerr = nf90_redef(ncid=ncid) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return ncerr = nf90_inq_dimid(ncid, "time", dimid=dimid) if (ncerr /= NF90_NOERR) then ! "time" dimension does not yet exist, define as unlimited dim ncerr = nf90_def_dim(ncid, "time", NF90_UNLIMITED, dimid=dimid) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return endif ncerr = nf90_def_var(ncid, "time", NF90_DOUBLE, & dimids=(/dimid/), varid=varid) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return ncerr = nf90_def_dim(ncid, "nchars", 20, ch_dimid) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return ncerr = nf90_def_var(ncid, "time_iso", NF90_CHAR, [ch_dimid,dimid], timeiso_varid) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return ncerr = nf90_put_att(ncid, timeiso_varid, "long_name", "valid time") if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return ncerr = nf90_put_att(ncid, timeiso_varid, "description", "ISO 8601 datetime string") if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return ncerr = nf90_put_att(ncid, timeiso_varid, "_Encoding", "UTF-8") if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return ncerr = nf90_enddef(ncid=ncid) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return ncerr = nf90_put_var(ncid, varid, values=time) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return ncerr = nf90_put_var(ncid, timeiso_varid, values=[trim(time_iso)]) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return ! loop over all the grid attributes that start with "time:", and ! put them on the "time" variable in the NetCDF file call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & name="TimeAttributes", itemCount=itemCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if (itemCount > 0) then ncerr = nf90_redef(ncid=ncid) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return allocate(attNameList(itemCount)) call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & name="TimeAttributes", valueList=attNameList, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return do i=1, itemCount attName = attNameList(i) call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & name=trim(attNameList(i)), typekind=typekind, rc=rc) ! print *,'in esmf call, att name=',trim(attNameList(i)) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if (typekind==ESMF_TYPEKIND_CHARACTER) then call ESMF_AttributeGet(grid, & convention="NetCDF", purpose="FV3", & name=trim(attNameList(i)), value=valueS, rc=rc) ! print *,'in esmf call, att string value=',trim(valueS) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ncerr = nf90_put_att(ncid, varid, & trim(attName(6:len(attName))), values=valueS) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return else if (typekind==ESMF_TYPEKIND_I4) then call ESMF_AttributeGet(grid, & convention="NetCDF", purpose="FV3", & name=trim(attNameList(i)), value=valueI4, rc=rc) ! print *,'in esmf call, att I4 value=',valueR8 if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ncerr = nf90_put_att(ncid, varid, & trim(attName(6:len(attName))), values=valueI4) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return else if (typekind==ESMF_TYPEKIND_R4) then call ESMF_AttributeGet(grid, & convention="NetCDF", purpose="FV3", & name=trim(attNameList(i)), value=valueR4, rc=rc) ! print *,'in esmf call, att r4 value=',valueR8 if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ncerr = nf90_put_att(ncid, varid, & trim(attName(6:len(attName))), values=valueR4) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return else if (typekind==ESMF_TYPEKIND_R8) then call ESMF_AttributeGet(grid, & convention="NetCDF", purpose="FV3", & name=trim(attNameList(i)), value=valueR8, rc=rc) ! print *,'in esmf call, att r8 value=',valueR8 if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ncerr = nf90_put_att(ncid, varid, & trim(attName(6:len(attName))), values=valueR8) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return endif enddo deallocate(attNameList) ncerr = nf90_enddef(ncid=ncid) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return endif endif ! close the NetCDF file ncerr = nf90_close(ncid=ncid) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return endif call ESMF_LogWrite("In ioCompRun() after "// & "writing vectical and time dimensions.", ESMF_LOGMSG_INFO, rc=rc) endif !TODO: remove this block once the ESMF_FieldBundleWrite() below allows to !TODO: specify the NetCDF 64-bit-offset file format. if (status==ESMF_FILESTATUS_REPLACE) then ! First time in with this filename (therefore 'replace'), create the ! file with 64bit-offset format in order to accommodate larger data ! volume. call ESMF_GridCompGet(comp, localPet=localPet, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if (localPet==0) then ! only single PET to deal with NetCDF ncerr = nf90_create(tileFileName, & cmode=or(nf90_clobber,nf90_64bit_offset), ncid=ncid) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return ncerr = nf90_close(ncid=ncid) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return endif status = ESMF_FILESTATUS_OLD ! switch status to 'OLD' to not overwrite call ESMF_LogWrite("In ioCompRun() after creating the NetCDF file", & ESMF_LOGMSG_INFO, rc=rc) endif if (timeslice==-1) then call ESMF_FieldBundleWrite(fieldbundle=wrtTileFB, fileName=tileFileName, & convention=convention, purpose=purpose, status=status, & overwrite=.true., rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return else call ESMF_FieldBundleWrite(fieldbundle=wrtTileFB, fileName=tileFileName, & convention=convention, purpose=purpose, status=status, & timeslice=timeslice, overwrite=.true., rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return endif call ESMF_LogWrite("In ioCompRun() after "// & "ESMF_FieldBundleWrite()", ESMF_LOGMSG_INFO, rc=rc) contains subroutine write_out_ungridded_dim_atts(dimLabel, rc) character(len=*) :: dimLabel integer, intent(out) :: rc ! inquire if NetCDF file already contains this ungridded dimension ncerr = nf90_inq_varid(ncid, trim(dimLabel), varid=varid) if (ncerr == NF90_NOERR) return ! the variable does not exist in the NetCDF file yet -> add it ! access the undistributed dimension attribute on the grid call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & name=trim(dimLabel), itemCount=valueCount, typekind=typekind, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if( typekind == ESMF_TYPEKIND_R4 ) then allocate(valueListr4(valueCount)) call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & name=trim(dimLabel), valueList=valueListr4, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return else if ( typekind == ESMF_TYPEKIND_R8) then allocate(valueListr8(valueCount)) call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & name=trim(dimLabel), valueList=valueListr8, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return endif ! now add it to the NetCDF file ncerr = nf90_redef(ncid=ncid) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return ncerr = nf90_inq_dimid(ncid, trim(dimLabel), dimid=dimid) if (ncerr /= NF90_NOERR) then ! dimension does not yet exist, and must be defined ncerr = nf90_def_dim(ncid, trim(dimLabel), valueCount, dimid=dimid) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return endif if( typekind == ESMF_TYPEKIND_R4 ) then ncerr = nf90_def_var(ncid, trim(dimLabel), NF90_FLOAT, & dimids=(/dimid/), varid=varid) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return ncerr = nf90_enddef(ncid=ncid) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return ncerr = nf90_put_var(ncid, varid, values=valueListr4) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return deallocate(valueListr4) else if(typekind == ESMF_TYPEKIND_R8) then ncerr = nf90_def_var(ncid, trim(dimLabel), NF90_DOUBLE, & dimids=(/dimid/), varid=varid) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return ncerr = nf90_enddef(ncid=ncid) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return ncerr = nf90_put_var(ncid, varid, values=valueListr8) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return deallocate(valueListr8) endif ! add attributes to this vertical variable call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & attnestflag=ESMF_ATTNEST_OFF, count=attCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if (attCount>0) then ncerr = nf90_redef(ncid=ncid) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return endif ! loop over all the attributes do j=1, attCount call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & attnestflag=ESMF_ATTNEST_OFF, attributeIndex=j, & name=attName, typekind=typekind, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! test for name starting with trim(dimLabel)":" if (index(trim(attName), trim(dimLabel)//":") == 1) then ind = len(trim(dimLabel)//":") ! found a matching attributes if (typekind == ESMF_TYPEKIND_CHARACTER) then call ESMF_AttributeGet(grid, & convention="NetCDF", purpose="FV3", & name=trim(attName), value=valueS, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ncerr = nf90_put_att(ncid, varid, & trim(attName(ind+1:len(attName))), values=valueS) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return else if (typekind == ESMF_TYPEKIND_I4) then call ESMF_AttributeGet(grid, & convention="NetCDF", purpose="FV3", & name=trim(attName), value=valueI4, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ncerr = nf90_put_att(ncid, varid, & trim(attName(ind+1:len(attName))), values=valueI4) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return else if (typekind == ESMF_TYPEKIND_R4) then call ESMF_AttributeGet(grid, & convention="NetCDF", purpose="FV3", & name=trim(attName), value=valueR4, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ncerr = nf90_put_att(ncid, varid, & trim(attName(ind+1:len(attName))), values=valueR4) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return else if (typekind == ESMF_TYPEKIND_R8) then call ESMF_AttributeGet(grid, & convention="NetCDF", purpose="FV3", & name=trim(attName), value=valueR8, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ncerr = nf90_put_att(ncid, varid, & trim(attName(ind+1:len(attName))), values=valueR8) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return endif endif enddo if (attCount>0) then ncerr = nf90_enddef(ncid=ncid) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return endif end subroutine end subroutine ioCompRun !----------------------------------------------------------------------------- subroutine ESMFproto_FieldMakeSingleTile(field, tile, tileField, petList, rc) type(ESMF_Field), intent(in) :: field integer, intent(in) :: tile type(ESMF_Field), intent(out) :: tileField integer, allocatable, intent(inout), optional :: petList(:) integer, intent(out), optional :: rc ! Take in a field on a multi-tile grid and return a field that only ! references a single tile. ! This routine only works with references, no data copies are being ! made. The single tile field that is returned points to the original ! field allocation. ! The original field passed in remains valid. type(ESMF_TypeKind_Flag) :: typekind type(ESMF_Index_Flag) :: indexflag type(ESMF_Grid) :: grid, tileGrid type(ESMF_Array) :: array type(ESMF_DistGrid) :: distgrid type(ESMF_DELayout) :: delayout character(40) :: fieldName integer :: fieldDimCount, gridDimCount integer :: undistDims integer :: localDeCount, deCount, tileCount integer, allocatable :: gridToFieldMap(:) integer, allocatable :: ungriddedLBound(:) integer, allocatable :: ungriddedUBound(:) integer, allocatable :: localDeToDeMap(:), deToTileMap(:) type(ESMF_LocalArray), allocatable :: lArrayList(:) integer :: i integer :: tileDeCount, tileLocalDeCount integer, allocatable :: distgridToArrayMap(:) integer, allocatable :: undistLBound(:), undistUBound(:) integer, allocatable :: petMap(:), tilePetMap(:) integer, allocatable :: minIndexPDe(:,:) integer, allocatable :: maxIndexPDe(:,:) integer, allocatable :: minIndexPTile(:,:) integer, allocatable :: maxIndexPTile(:,:) integer, allocatable :: deBlockList(:,:,:) character(800) :: msg if (present(rc)) rc = ESMF_SUCCESS ! access information from the incoming field call ESMF_FieldGet(field, array=array, typekind=typekind, & dimCount=fieldDimCount, name=fieldName, grid=grid, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! access information from the associated grid call ESMF_GridGet(grid, dimCount=gridDimCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return #if 0 write(msg,*) "fieldDimCount=",fieldDimCount,"gridDimCount=",gridDimCount call ESMF_LogWrite(msg, ESMF_LOGMSG_INFO, rc=rc) #endif ! access list type information from the incoming field allocate(gridToFieldMap(gridDimCount)) undistDims = fieldDimCount-gridDimCount if (undistDims < 0) undistDims = 0 ! this supports replicated dimensions allocate(ungriddedLBound(undistDims)) allocate(ungriddedUBound(undistDims)) call ESMF_FieldGet(field, gridToFieldMap=gridToFieldMap, & ungriddedLBound=ungriddedLBound, ungriddedUBound=ungriddedUBound, & rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! access information from associated array call ESMF_ArrayGet(array, distgrid=distgrid, delayout=delayout, & indexflag=indexflag, localDeCount=localDeCount, deCount=deCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! access list type information from associated array allocate(localDeToDeMap(localDeCount), deToTileMap(deCount)) allocate(distgridToArrayMap(gridDimCount)) allocate(undistLBound(undistDims)) allocate(undistUBound(undistDims)) call ESMF_ArrayGet(array, tileCount=tileCount, & localDeToDeMap=localDeToDeMap, deToTileMap=deToTileMap, & distgridToArrayMap=distgridToArrayMap, & undistLBound=undistLBound, undistUBound=undistUBound, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! access list type information from associated distgrid allocate(minIndexPDe(gridDimCount,deCount)) allocate(maxIndexPDe(gridDimCount,deCount)) allocate(minIndexPTile(gridDimCount,tileCount)) allocate(maxIndexPTile(gridDimCount,tileCount)) call ESMF_DistGridGet(distgrid, & minIndexPDe=minIndexPDe, maxIndexPDe=maxIndexPDe, & minIndexPTile=minIndexPTile, maxIndexPTile=maxIndexPTile, & rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! access list type information from associated delayout allocate(petMap(deCount)) call ESMF_DELayoutGet(delayout, petMap=petMap, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! construct data structures selecting specific tile allocate(lArrayList(localDeCount)) tileLocalDeCount = 0 do i=1, localDeCount if (deToTileMap(localDeToDeMap(i)+1) == tile) then ! localDe is on tile tileLocalDeCount = tileLocalDeCount + 1 call ESMF_ArrayGet(array, localDe=i-1, & localarray=lArrayList(tileLocalDeCount), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return endif enddo allocate(tilePetMap(deCount)) allocate(deBlockList(gridDimCount,2,deCount)) tileDeCount=0 do i=1, deCount if (deToTileMap(i) == tile) then ! DE is on tile tileDeCount = tileDeCount + 1 tilePetMap(tileDeCount) = petMap(i) deBlockList(:,1,tileDeCount) = minIndexPDe(:,i) deBlockList(:,2,tileDeCount) = maxIndexPDe(:,i) endif enddo if (present(petList)) then if (.not.allocated(petList)) then allocate(petList(tileDeCount)) else if (size(petList)/=tileDeCount) then deallocate(petList) allocate(petList(tileDeCount)) endif petList(:) = tilePetMap(1:tileDeCount) endif ! create DELayout and DistGrid that only contain the single tile delayout = ESMF_DELayoutCreate(tilePetMap(1:tileDeCount), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return distgrid = ESMF_DistGridCreate(minIndex=minIndexPTile(:,tile), & maxIndex=maxIndexPTile(:,tile), delayout=delayout, & deBlockList=deBlockList(:,:,1:tileDeCount), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! create an Array that only holds tile specific allocations if (tileLocalDeCount>0) then array = ESMF_ArrayCreate(distgrid, lArrayList(1:tileLocalDeCount), & indexflag=indexflag, & distgridToArrayMap=distgridToArrayMap, undistLBound=undistLBound, & undistUBound=undistUBound, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return else array = ESMF_ArrayCreate(distgrid, typekind, indexflag=indexflag, & distgridToArrayMap=distgridToArrayMap, undistLBound=undistLBound, & undistUBound=undistUBound, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return endif ! create a grid on the new distgrid tileGrid = ESMF_GridCreate(distgrid, indexflag=indexflag, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! alias the Attributes on grid level call ESMF_AttributeCopy(grid, tileGrid, attcopy=ESMF_ATTCOPY_REFERENCE, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! create the tile specific field from the array tileField = ESMF_FieldCreate(tileGrid, array=array, name=fieldName, & gridToFieldMap=gridToFieldMap, ungriddedLBound=ungriddedLBound, & ungriddedUBound=ungriddedUBound, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! alias the Attributes on field level call ESMF_AttributeCopy(field, tileField, attcopy=ESMF_ATTCOPY_REFERENCE, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! local garbage collection deallocate(localDeToDeMap, deToTileMap) deallocate(petMap, tilePetMap) deallocate(minIndexPDe, maxIndexPDe) deallocate(minIndexPTile, maxIndexPTile) deallocate(gridToFieldMap, ungriddedLBound, ungriddedUBound) deallocate(deBlockList) end subroutine ESMFproto_FieldMakeSingleTile ! !----------------------------------------------------------------------- subroutine splat4(idrt,jmax,aslat) implicit none integer,intent(in) :: idrt,jmax real(4),intent(out) :: aslat(jmax) ! integer,parameter :: KD=SELECTED_REAL_KIND(15,45) real(kind=KD) :: pk(jmax/2),pkm1(jmax/2),pkm2(jmax/2) real(kind=KD) :: aslatd(jmax/2),sp,spmax,eps=10.d0*epsilon(sp) integer,PARAMETER :: JZ=50 real(8) bz(jz) data bz / 2.4048255577d0, 5.5200781103d0, & 8.6537279129d0, 11.7915344391d0, 14.9309177086d0, 18.0710639679d0, & 21.2116366299d0, 24.3524715308d0, 27.4934791320d0, 30.6346064684d0, & 33.7758202136d0, 36.9170983537d0, 40.0584257646d0, 43.1997917132d0, & 46.3411883717d0, 49.4826098974d0, 52.6240518411d0, 55.7655107550d0, & 58.9069839261d0, 62.0484691902d0, 65.1899648002d0, 68.3314693299d0, & 71.4729816036d0, 74.6145006437d0, 77.7560256304d0, 80.8975558711d0, & 84.0390907769d0, 87.1806298436d0, 90.3221726372d0, 93.4637187819d0, & 96.6052679510d0, 99.7468198587d0, 102.888374254d0, 106.029930916d0, & 109.171489649d0, 112.313050280d0, 115.454612653d0, 118.596176630d0, & 121.737742088d0, 124.879308913d0, 128.020877005d0, 131.162446275d0, & 134.304016638d0, 137.445588020d0, 140.587160352d0, 143.728733573d0, & 146.870307625d0, 150.011882457d0, 153.153458019d0, 156.295034268d0 / real(8) :: dlt integer :: jhe,jho ! real(8),parameter :: PI=3.14159265358979d0,C=(1.d0-(2.d0/PI)**2)*0.25d0 real(8),parameter :: C=(1.d0-(2.d0/PI)**2)*0.25d0 real(8) r integer jh,n,j ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! GAUSSIAN LATITUDES IF(IDRT.EQ.4) THEN JH=JMAX/2 JHE=(JMAX+1)/2 R=1.d0/SQRT((JMAX+0.5d0)**2+C) DO J=1,MIN(JH,JZ) ASLATD(J)=COS(BZ(J)*R) ENDDO DO J=JZ+1,JH ASLATD(J)=COS((BZ(JZ)+(J-JZ)*PI)*R) ENDDO SPMAX=1.d0 DO WHILE(SPMAX.GT.EPS) SPMAX=0.d0 DO J=1,JH PKM1(J)=1.d0 PK(J)=ASLATD(J) ENDDO DO N=2,JMAX DO J=1,JH PKM2(J)=PKM1(J) PKM1(J)=PK(J) PK(J)=((2*N-1)*ASLATD(J)*PKM1(J)-(N-1)*PKM2(J))/N ENDDO ENDDO DO J=1,JH SP=PK(J)*(1.d0-ASLATD(J)**2)/(JMAX*(PKM1(J)-ASLATD(J)*PK(J))) ASLATD(J)=ASLATD(J)-SP SPMAX=MAX(SPMAX,ABS(SP)) ENDDO ENDDO ! DO J=1,JH ASLAT(J)=ASLATD(J) ASLAT(JMAX+1-J)=-ASLAT(J) ENDDO IF(JHE.GT.JH) THEN ASLAT(JHE)=0.d0 ENDIF ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! EQUALLY-SPACED LATITUDES INCLUDING POLES ELSEIF(IDRT.EQ.0) THEN JH=JMAX/2 JHE=(JMAX+1)/2 JHO=JHE-1 DLT=PI/(JMAX-1) ASLAT(1)=1.d0 DO J=2,JH ASLAT(J)=COS((J-1)*DLT) ENDDO ! DO J=1,JH ASLAT(JMAX+1-J)=-ASLAT(J) ENDDO IF(JHE.GT.JH) THEN ASLAT(JHE)=0.d0 ENDIF ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! EQUALLY-SPACED LATITUDES EXCLUDING POLES ELSEIF(IDRT.EQ.256) THEN JH=JMAX/2 JHE=(JMAX+1)/2 JHO=JHE DLT=PI/JMAX ASLAT(1)=1.d0 DO J=1,JH ASLAT(J)=COS((J-0.5)*DLT) ENDDO DO J=1,JH ASLAT(JMAX+1-J)=-ASLAT(J) ENDDO IF(JHE.GT.JH) THEN ASLAT(JHE)=0.d0 ENDIF ENDIF ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end subroutine splat4 !---------------------------------------------------------------------- subroutine splat8(idrt,jmax,aslat) !$$$ implicit none integer,intent(in) :: idrt,jmax real(8),intent(out) :: aslat(jmax) ! integer,parameter :: KD=SELECTED_REAL_KIND(15,45) real(kind=KD) :: pk(jmax/2),pkm1(jmax/2),pkm2(jmax/2) real(kind=KD) :: aslatd(jmax/2),sp,spmax,eps=10.d0*epsilon(sp) integer,parameter :: jz=50 real(8) bz(jz) data bz / 2.4048255577d0, 5.5200781103d0, & 8.6537279129d0, 11.7915344391d0, 14.9309177086d0, 18.0710639679d0, & 21.2116366299d0, 24.3524715308d0, 27.4934791320d0, 30.6346064684d0, & 33.7758202136d0, 36.9170983537d0, 40.0584257646d0, 43.1997917132d0, & 46.3411883717d0, 49.4826098974d0, 52.6240518411d0, 55.7655107550d0, & 58.9069839261d0, 62.0484691902d0, 65.1899648002d0, 68.3314693299d0, & 71.4729816036d0, 74.6145006437d0, 77.7560256304d0, 80.8975558711d0, & 84.0390907769d0, 87.1806298436d0, 90.3221726372d0, 93.4637187819d0, & 96.6052679510d0, 99.7468198587d0, 102.888374254d0, 106.029930916d0, & 109.171489649d0, 112.313050280d0, 115.454612653d0, 118.596176630d0, & 121.737742088d0, 124.879308913d0, 128.020877005d0, 131.162446275d0, & 134.304016638d0, 137.445588020d0, 140.587160352d0, 143.728733573d0, & 146.870307625d0, 150.011882457d0, 153.153458019d0, 156.295034268d0 / real(8) :: dlt integer(4) :: jhe,jho ! real(8),parameter :: PI=3.14159265358979d0,C=(1.d0-(2.d0/PI)**2)*0.25d0 real(8),parameter :: C=(1.d0-(2.d0/PI)**2)*0.25d0 real(8) r integer jh,n,j ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! GAUSSIAN LATITUDES IF(IDRT.EQ.4) THEN JH=JMAX/2 JHE=(JMAX+1)/2 R=1.d0/SQRT((JMAX+0.5d0)**2+C) DO J=1,MIN(JH,JZ) ASLATD(J)=COS(BZ(J)*R) ENDDO DO J=JZ+1,JH ASLATD(J)=COS((BZ(JZ)+(J-JZ)*PI)*R) ENDDO SPMAX=1.d0 DO WHILE(SPMAX.GT.EPS) SPMAX=0.d0 DO J=1,JH PKM1(J)=1.d0 PK(J)=ASLATD(J) ENDDO DO N=2,JMAX DO J=1,JH PKM2(J)=PKM1(J) PKM1(J)=PK(J) PK(J)=((2*N-1)*ASLATD(J)*PKM1(J)-(N-1)*PKM2(J))/N ENDDO ENDDO DO J=1,JH SP=PK(J)*(1.d0-ASLATD(J)**2)/(JMAX*(PKM1(J)-ASLATD(J)*PK(J))) ASLATD(J)=ASLATD(J)-SP SPMAX=MAX(SPMAX,ABS(SP)) ENDDO ENDDO ! DO J=1,JH ASLAT(J)=ASLATD(J) ASLAT(JMAX+1-J)=-ASLAT(J) ENDDO IF(JHE.GT.JH) THEN ASLAT(JHE)=0.d0 ENDIF ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! EQUALLY-SPACED LATITUDES INCLUDING POLES ELSEIF(IDRT.EQ.0) THEN JH=JMAX/2 JHE=(JMAX+1)/2 JHO=JHE-1 DLT=PI/(JMAX-1) ASLAT(1)=1.d0 DO J=2,JH ASLAT(J)=COS((J-1)*DLT) ENDDO DO J=1,JH ASLAT(JMAX+1-J)=-ASLAT(J) ENDDO IF(JHE.GT.JH) THEN ASLAT(JHE)=0.d0 ENDIF ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! EQUALLY-SPACED LATITUDES EXCLUDING POLES ELSEIF(IDRT.EQ.256) THEN JH=JMAX/2 JHE=(JMAX+1)/2 JHO=JHE DLT=PI/JMAX ASLAT(1)=1.d0 DO J=1,JH ASLAT(J)=COS((J-0.5d0)*DLT) ENDDO ! DO J=1,JH ASLAT(JMAX+1-J)=-ASLAT(J) ENDDO IF(JHE.GT.JH) THEN ASLAT(JHE)=0.d0 ENDIF ENDIF ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end subroutine splat8 ! ! subroutine rtll(tlmd,tphd,almd,aphd,tlm0d,tph0d) !------------------------------------------------------------------------------- real(ESMF_KIND_R8), intent(in) :: tlmd, tphd real(ESMF_KIND_R8), intent(out) :: almd, aphd real(ESMF_KIND_R8), intent(in) :: tph0d, tlm0d !------------------------------------------------------------------------------- ! real(ESMF_KIND_R8), parameter :: pi=3.14159265358979323846 real(ESMF_KIND_R8), parameter :: dtr=pi/180.0 ! real(ESMF_KIND_R8) :: tph0, ctph0, stph0, tlm, tph, stph, ctph, ctlm, stlm, aph, cph real(ESMF_KIND_R8) :: xx, yy !------------------------------------------------------------------------------- ! tph0=tph0d*dtr ctph0=cos(tph0) stph0=sin(tph0) ! tlm=tlmd*dtr tph=tphd*dtr stph=sin(tph) ctph=cos(tph) ctlm=cos(tlm) stlm=sin(tlm) ! xx=stph0*ctph*ctlm+ctph0*stph xx=max(xx,-1.0) xx=min(xx, 1.0) aph=asin(xx) cph=cos(aph) ! xx=(ctph0*ctph*ctlm-stph0*stph)/cph xx=max(xx,-1.0) xx=min(xx, 1.0) xx=acos(xx)/dtr yy=ctph*stlm/cph xx=sign(xx,yy) almd=tlm0d+xx aphd=aph/dtr ! if (almd > 180.0) then almd=almd-360.0 end if if (almd < -180.0) then almd=almd+360.0 end if ! return ! end subroutine rtll ! !----------------------------------------------------------------------- ! subroutine lambert(stlat1,stlat2,c_lat,c_lon,glon,glat,x,y,inv) !------------------------------------------------------------------------------- real(ESMF_KIND_R8), intent(in) :: stlat1,stlat2,c_lat,c_lon real(ESMF_KIND_R8), intent(inout) :: glon, glat real(ESMF_KIND_R8), intent(inout) :: x, y integer, intent(in) :: inv !------------------------------------------------------------------------------- ! real(ESMF_KIND_R8), parameter :: pi = 3.14159265358979323846 real(ESMF_KIND_R8), parameter :: dtor = pi/180.0 real(ESMF_KIND_R8), parameter :: rtod = 180.0/pi real(ESMF_KIND_R8), parameter :: a = 6371200.0 !------------------------------------------------------------------------------- ! inv == 1 (glon,glat) ---> (x,y) ! inv == -1 (x,y) ---> (glon,glat) real(ESMF_KIND_R8) :: xp, yp, en, de, rho, rho0, rho2, dlon, theta, dr2 real(ESMF_KIND_R8) :: h = 1.0 ! For reference see: ! John P. Snyder (1987), Map projections: A working manual (pp. 104-110) ! https://doi.org/10.3133/pp1395 if (stlat1 == stlat2) then en = sin(stlat1*dtor) else en = log(cos(stlat1*dtor)/cos(stlat2*dtor)) / & log(tan((45+0.5*stlat2)*dtor)/tan((45+0.5*stlat1)*dtor)) ! (15-3) endif h = sign(1.0_ESMF_KIND_R8,en) de = a*(cos(stlat1*dtor)*tan((45+0.5*stlat1)*dtor)**en)/en ! (15-2) rho0 = de/(tan((45+0.5*c_lat)*dtor)**en) ! (15-1a) if (inv == 1) then ! FORWARD TRANSFORMATION rho = de/(tan((45+0.5*glat)*dtor)**en) ! (15-1) dlon = modulo(glon-c_lon+180.0+3600.0,360.0)-180.0 theta = en*dlon*dtor ! (14-4) x = rho*sin(theta) ! (14-1) y = rho0-rho*cos(theta) ! (14-2) else if (inv == -1) then ! INVERSE TRANSFORMATION xp = h*x; yp = h*(rho0-y) theta = atan2(xp,yp) ! (14-11) glon = c_lon+(theta/en)*rtod ! (14-9) glon = modulo(glon+180.0+3600.0,360.0)-180.0 rho2 = xp*xp+yp*yp ! (14-10) if (rho2 == 0.0) then glat = h*90.0 else glat = 2.0*atan((de*de/rho2)**(1.0/(2.0*en)))*rtod-90.0 ! (15-5) endif else write (unit=*,fmt=*) " lambert: unknown inv argument" return end if return end subroutine lambert ! !----------------------------------------------------------------------- ! subroutine get_outfile(nfl, filename, outfile_name,noutfile) integer, intent(in) :: nfl character(*), intent(in) :: filename(:,:) character(*), intent(inout) :: outfile_name(:) integer, intent(inout) :: noutfile integer :: i,j,n logical :: found ! noutfile = 0 do i=1,nfl loopj: do j=1, 2000 if( trim(filename(j,i)) == '') exit loopj if( trim(filename(j,i)) == 'none') cycle found = .false. loopn: do n=1, noutfile if(trim(filename(j,i)) == trim(outfile_name(n))) then found = .true. exit loopn endif enddo loopn if (.not.found) then noutfile = noutfile + 1 outfile_name(noutfile) = trim(filename(j,i)) ! print *,'in get outfile,noutfile=', noutfile,' outfile=',trim(filename(j,i)) endif enddo loopj ! enddo end subroutine get_outfile ! !----------------------------------------------------------------------- !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& !----------------------------------------------------------------------- ! end module module_wrt_grid_comp ! !-----------------------------------------------------------------------