!----------------------------------------------------------------------- ! module module_wrt_grid_comp ! !----------------------------------------------------------------------- !*** This module includes the funcationailty of write gridded component. !----------------------------------------------------------------------- !*** At intialization step, write grid is defined. The firecast 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 mirrror field bundle on write grid component. Also the IO_BaseTime !*** is set to the intial 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 !*** writen 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 ! Aug 2017: J. Wang - add nemsio binary output for Gaussian grid ! Mar 2018: S Moorthi - changing cfhour to accomodate up to 99999 hours ! Aug 2019: J. Wang - add inline post ! !--------------------------------------------------------------------------------- ! use fms_io_mod, only: field_exist, read_data use esmf use write_internal_state use module_fv3_io_def, only : num_pes_fcst,lead_wrttask, last_wrttask, & n_group, num_files, app_domain, & filename_base, output_grid, output_file, & imo,jmo,ichunk2d,jchunk2d,write_nemsioflip,& ichunk3d,jchunk3d,kchunk3d, & nsout => nsout_io, & cen_lon, cen_lat, & lon1, lat1, lon2, lat2, dlon, dlat, & stdlat1, stdlat2, dx, dy, iau_offset use module_write_nemsio, only : nemsio_first_call, write_nemsio use module_write_netcdf, only : write_netcdf use physcons, only : pi => con_pi use post_gfs, only : post_run_gfs, post_getattr_gfs use module_write_netcdf_parallel, only : write_netcdf_parallel ! !----------------------------------------------------------------------- ! implicit none ! include 'mpif.h' ! !----------------------------------------------------------------------- private ! !----------------------------------------------------------------------- ! integer,parameter :: filename_maxstr=255 real, parameter :: rdgas=287.04, grav=9.80 real, parameter :: stndrd_atmos_ps = 101325. real, parameter :: stndrd_atmos_lapse = 0.0065 ! 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 :: mytile !<-- the tile number in write task integer,save :: wrt_mpi_comm !<-- the mpi communicator in the write comp integer,save :: idate(7) logical,save :: first_init=.false. logical,save :: first_run=.false. logical,save :: first_getlatlon=.true. logical,save :: first_getmaskwrt=.true. !<-- for mask the output grid of the write comp logical,save :: change_wrtidate=.false. ! !----------------------------------------------------------------------- ! type(wrt_internal_state),pointer :: wrt_int_state ! The internal state pointer. type(ESMF_FieldBundle) :: gridFB integer :: FBcount character(len=80),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, & userRoutine=wrt_initialize, rc=rc) if(rc/=0) write(*,*)'Error: write grid comp, initial' ! call ESMF_GridCompSetEntryPoint(wrt_comp, ESMF_METHOD_RUN, & userRoutine=wrt_run, rc=rc) if(rc/=0) write(*,*)'Error: write grid comp, run' ! call ESMF_GridCompSetEntryPoint(wrt_comp, ESMF_METHOD_FINALIZE, & userRoutine=wrt_finalize, rc=rc) if(rc/=0) write(*,*)'Error: write grid comp, run' end subroutine SetServices ! !----------------------------------------------------------------------- !####################################################################### !----------------------------------------------------------------------- ! subroutine wrt_initialize(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 :: ISTAT, 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, axesname,longname type(ESMF_Config) :: cf type(ESMF_DELayout) :: delayout type(ESMF_Grid) :: wrtGrid, fcstGrid type(ESMF_Array) :: array_work, array type(ESMF_FieldBundle) :: fieldbdl_work type(ESMF_Field) :: field_work, field character(len=80) :: attrValueSList(2) type(ESMF_StateItem_Flag), allocatable :: fcstItemTypeList(:) type(ESMF_FieldBundle) :: fcstFB, wrtFB, fieldbundle type(ESMF_Field), allocatable :: fcstField(:) type(ESMF_TypeKind_Flag) :: typekind character(len=80), allocatable :: fieldnamelist(:) integer :: fieldDimCount, gridDimCount 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 integer :: attCount, axeslen, 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 type(ESMF_TimeInterval) :: IAU_offsetTI type(ESMF_DataCopy_Flag) :: copyflag=ESMF_DATACOPY_REFERENCE ! real(8),parameter :: PI=3.14159265358979d0 character(256) :: gridfile ! logical,save :: first=.true. logical :: lprnt !test integer myattCount real(ESMF_KIND_R8),dimension(:,:), pointer :: glatPtr, glonPtr ! !----------------------------------------------------------------------- !*********************************************************************** !----------------------------------------------------------------------- ! 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 allocate(wrt_int_state%wrtFB(num_files)) ! 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 if( wrt_int_state%write_dopost ) then #ifdef NO_INLINE_POST 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 call esmf_configgetattribute(cf,wrt_int_state%post_nlunit,default=777,label='nlunit:',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%post_namelist,default='itag', & label ='post_namelist:',rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__)) return endif ! !----------------------------------------------------------------------- !*** Create the cubed sphere grid with field on PETs !*** first try: Create cubed sphere grid from file !----------------------------------------------------------------------- ! if ( trim(output_grid) == 'cubed_sphere_grid' ) then mytile = mod(wrt_int_state%mype,ntasks)+1 if ( trim(app_domain) == 'global' ) then do tl=1,6 decomptile(1,tl) = 1 decomptile(2,tl) = jidx 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 = ESMF_GridCreateMosaic(filename="INPUT/"//trim(gridfile), & regDecompPTile=decomptile,tileFilePath="INPUT/", & 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(trim(app_domain) == 'nested') then gridfile='grid.nest02.tile7.nc' else if(trim(app_domain) == 'regional') then gridfile='grid.tile7.halo0.nc' endif 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 wrtGrid = ESMF_GridCreate(filename="INPUT/"//trim(gridfile), & fileformat=ESMF_FILEFORMAT_GRIDSPEC, regDecomp=regDecomp, & delayout=delayout, isSphere=.false., indexflag=ESMF_INDEX_DELOCAL, & rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return print *,'in nested/regional cubed_sphere grid, regDecomp=',regDecomp,' PetMap=',petMap(1),petMap(ntasks), & 'gridfile=',trim(gridfile) deallocate(petMap) endif else if ( trim(output_grid) == 'gaussian_grid') then wrtgrid = ESMF_GridCreate1PeriDim(minIndex=(/1,1/), & maxIndex=(/imo,jmo/), regDecomp=(/1,ntasks/), & indexflag=ESMF_INDEX_GLOBAL, & name='wrt_grid',rc=rc) ! indexflag=ESMF_INDEX_GLOBAL, coordSys=ESMF_COORDSYS_SPH_DEG if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_GridAddCoord(wrtgrid, staggerLoc=ESMF_STAGGERLOC_CENTER, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_GridGetCoord(wrtgrid, coordDim=1, farrayPtr=lonPtr, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_GridGetCoord(wrtgrid, coordDim=2, farrayPtr=latPtr, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! allocate(slat(jmo), lat(jmo), lon(imo)) call splat(4, jmo, slat) if(write_nemsioflip) then do j=1,jmo lat(j) = asin(slat(j)) * radi enddo else do j=1,jmo lat(jmo-j+1) = asin(slat(j)) * radi enddo endif wrt_int_state%latstart = lat(1) wrt_int_state%latlast = lat(jmo) do j=1,imo lon(j) = 360.d0/real(imo,8) *real(j-1,8) enddo wrt_int_state%lonstart = lon(1) wrt_int_state%lonlast = lon(imo) do j=lbound(latPtr,2),ubound(latPtr,2) do i=lbound(lonPtr,1),ubound(lonPtr,1) lonPtr(i,j) = 360.d0/real(imo,8) * real(i-1,8) latPtr(i,j) = lat(j) enddo enddo ! print *,'aft wrtgrd, Gaussian, dimi,i=',lbound(lonPtr,1),ubound(lonPtr,1), & ! ' j=',lbound(lonPtr,2),ubound(lonPtr,2),'imo=',imo,'jmo=',jmo ! if(wrt_int_state%mype==0) print *,'aft wrtgrd, lon=',lonPtr(1:5,1), & ! 'lat=',latPtr(1,1:5),'imo,jmo=',imo,jmo ! lonPtr(lbound(lonPtr,1),ubound(lonPtr,2)),'lat=',latPtr(lbound(lonPtr,1),lbound(lonPtr,2)), & ! latPtr(lbound(lonPtr,1),ubound(lonPtr,2)) wrt_int_state%lat_start = lbound(latPtr,2) wrt_int_state%lat_end = ubound(latPtr,2) wrt_int_state%lon_start = lbound(lonPtr,1) wrt_int_state%lon_end = ubound(lonPtr,1) allocate( wrt_int_state%lat_start_wrtgrp(wrt_int_state%petcount)) allocate( wrt_int_state%lat_end_wrtgrp (wrt_int_state%petcount)) call mpi_allgather(wrt_int_state%lat_start,1,MPI_INTEGER, & wrt_int_state%lat_start_wrtgrp, 1, MPI_INTEGER, wrt_mpi_comm, rc) call mpi_allgather(wrt_int_state%lat_end, 1,MPI_INTEGER, & wrt_int_state%lat_end_wrtgrp, 1, MPI_INTEGER, wrt_mpi_comm, rc) if( lprnt ) print *,'aft wrtgrd, Gaussian, dimj_start=',wrt_int_state%lat_start_wrtgrp, & 'dimj_end=',wrt_int_state%lat_end_wrtgrp, 'wrt_group=',n_group allocate( wrt_int_state%latPtr(wrt_int_state%lon_start:wrt_int_state%lon_end, & wrt_int_state%lat_start:wrt_int_state%lat_end)) allocate( wrt_int_state%lonPtr(wrt_int_state%lon_start:wrt_int_state%lon_end, & wrt_int_state%lat_start:wrt_int_state%lat_end)) do j=wrt_int_state%lat_start,wrt_int_state%lat_end do i=wrt_int_state%lon_start,wrt_int_state%lon_end wrt_int_state%latPtr(i,j) = latPtr(i,j) wrt_int_state%lonPtr(i,j) = lonPtr(i,j) enddo enddo wrt_int_state%im = imo wrt_int_state%jm = jmo wrt_int_state%post_maptype = 4 deallocate(slat) else if ( trim(output_grid) == 'latlon_grid') then wrtgrid = ESMF_GridCreate1PeriDimUfrm(maxIndex=(/imo, jmo/), & minCornerCoord=(/0._ESMF_KIND_R8, -80._ESMF_KIND_R8/), & maxCornerCoord=(/360._ESMF_KIND_R8, 80._ESMF_KIND_R8/), & staggerLocList=(/ESMF_STAGGERLOC_CENTER, ESMF_STAGGERLOC_CORNER/), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! if(wrt_int_state%mype == lead_write_task) print *,'af wrtgrd, latlon,rc=',rc, & ! 'imo=',imo,' jmo=',jmo call ESMF_GridGetCoord(wrtgrid, coordDim=1, farrayPtr=lonPtr, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__)) return ! bail out call ESMF_GridGetCoord(wrtgrid, coordDim=2, farrayPtr=latPtr, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__)) return ! bail out wrt_int_state%im = imo wrt_int_state%jm = jmo wrt_int_state%lat_start = lbound(latPtr,2) wrt_int_state%lat_end = ubound(latPtr,2) wrt_int_state%lon_start = lbound(lonPtr,2) wrt_int_state%lon_end = ubound(lonPtr,2) allocate( wrt_int_state%latPtr(wrt_int_state%lon_start:wrt_int_state%lon_end, & wrt_int_state%lat_start:wrt_int_state%lat_end)) allocate( wrt_int_state%lonPtr(wrt_int_state%lon_start:wrt_int_state%lon_end, & wrt_int_state%lat_start:wrt_int_state%lat_end)) do j=wrt_int_state%lat_start,wrt_int_state%lat_end do i=wrt_int_state%lon_start,wrt_int_state%lon_end wrt_int_state%latPtr(i,j) = latPtr(i,j) wrt_int_state%lonPtr(i,j) = lonPtr(i,j) enddo enddo wrt_int_state%post_maptype = 0 else if ( trim(output_grid) == 'regional_latlon' .or. & trim(output_grid) == 'rotated_latlon' .or. & trim(output_grid) == 'lambert_conformal' ) then wrtgrid = ESMF_GridCreate1PeriDim(minIndex=(/1,1/), & maxIndex=(/imo,jmo/), regDecomp=(/1,ntasks/), & 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, staggerLoc=ESMF_STAGGERLOC_CENTER, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_GridGetCoord(wrtgrid, coordDim=1, farrayPtr=lonPtr, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_GridGetCoord(wrtgrid, coordDim=2, farrayPtr=latPtr, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if ( trim(output_grid) == 'regional_latlon' ) then do j=lbound(lonPtr,2),ubound(lonPtr,2) do i=lbound(lonPtr,1),ubound(lonPtr,1) lonPtr(i,j) = lon1 + (lon2-lon1)/(imo-1) * (i-1) latPtr(i,j) = lat1 + (lat2-lat1)/(jmo-1) * (j-1) enddo enddo else if ( trim(output_grid) == 'rotated_latlon' ) then do j=lbound(lonPtr,2),ubound(lonPtr,2) do i=lbound(lonPtr,1),ubound(lonPtr,1) rot_lon = lon1 + (lon2-lon1)/(imo-1) * (i-1) rot_lat = lat1 + (lat2-lat1)/(jmo-1) * (j-1) call rtll(rot_lon, rot_lat, geo_lon, geo_lat, dble(cen_lon), dble(cen_lat)) if (geo_lon < 0.0) geo_lon = geo_lon + 360.0 lonPtr(i,j) = geo_lon latPtr(i,j) = geo_lat enddo enddo else if ( trim(output_grid) == 'lambert_conformal' ) then lon1_r8 = dble(lon1) lat1_r8 = dble(lat1) call lambert(dble(stdlat1),dble(stdlat2),dble(cen_lat),dble(cen_lon), & 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 * (i-1) y = y1 + dy * (j-1) call lambert(dble(stdlat1),dble(stdlat2),dble(cen_lat),dble(cen_lon), & 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: Unknown output_grid ", trim(output_grid) call ESMF_LogWrite("wrt_initialize: Unknown output_grid "//trim(output_grid),ESMF_LOGMSG_ERROR,rc=RC) call ESMF_Finalize(endflag=ESMF_END_ABORT) endif ! !----------------------------------------------------------------------- !*** 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 ! ! Create field bundle !------------------------------------------------------------------- ! !--- check grid dim count first call ESMF_GridGet(wrtgrid, dimCount=gridDimCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! !--- Look at the incoming FieldBundles in the imp_state_write, and mirror them ! 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, 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 ! create a mirror FieldBundle and add it to importState fieldbundle = 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, (/fieldbundle/), 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, fieldbundle, attcopy=ESMF_ATTCOPY_REFERENCE, 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 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 mirror field CALL ESMF_LogWrite("call field create on wrt comp",ESMF_LOGMSG_INFO,rc=RC) field_work = ESMF_FieldCreate(wrtGrid, 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(i), 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 mirror field to the mirror FieldBundle call ESMF_FieldBundleAdd(fieldbundle, (/field_work/), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! local garbage collection deallocate(gridToFieldMap, ungriddedLBound, ungriddedUBound) enddo ! call ESMF_AttributeCopy(fcstGrid, wrtGrid, & 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 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="mirror_"//trim(fcstItemNameList(n)), & fieldbundle=fcstFB, rc=rc) if( index(trim(fcstItemNameList(n)),trim(FBlist_outfilename(i))) > 0 ) then ! ! copy the mirror 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_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) == '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) == '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, rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="jm", value=jmo, rc=rc) else if (trim(output_grid) == 'latlon_grid') then 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=0, rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="lat1", value=-80, rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="lon2", value=360, rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="lat2", value=80, rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="dlon", value=360.0/(imo-1), rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="dlat", value=160.0/(jmo-1), rc=rc) else if (trim(output_grid) == 'regional_latlon') then 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, rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="lat1", value=lat1, rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="lon2", value=lon2, rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="lat2", value=lat2, rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="dlon", value=dlon, rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="dlat", value=dlat, rc=rc) else if (trim(output_grid) == 'rotated_latlon') 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) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="cen_lon", value=cen_lon, rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="cen_lat", value=cen_lat, rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="lon1", value=lon1, rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="lat1", value=lat1, rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="lon2", value=lon2, rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="lat2", value=lat2, rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="dlon", value=dlon, rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="dlat", value=dlat, rc=rc) else if (trim(output_grid) == '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, rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="cen_lat", value=cen_lat, rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="stdlat1", value=stdlat1, rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="stdlat2", value=stdlat2, rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="nx", value=imo, rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="ny", value=jmo, rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="lat1", value=lat1, rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="lon1", value=lon1, rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="dx", value=dx, rc=rc) call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & name="dy", value=dy, 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 ! add the transfer attributes from importState to grid call ESMF_AttributeAdd(wrtgrid, 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, 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, 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, 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, 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, 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, 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 deallocate(attNameList, attNameList2, typekindList) ! !*** 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 call ESMF_AttributeGet(wrtGrid, convention="NetCDF", purpose="FV3", & name="ESMF:gridded_dim_labels", valueList=attrValueSList, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_GridGetCoord(wrtGrid, coordDim=1, & staggerloc=ESMF_STAGGERLOC_CENTER, array=array, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! write(0,*) 'create gridFB,fieldname=',trim(attrValueSList(1)),trim(attrValueSList(2)), & ! 'lon value=',array(1:5) field = ESMF_FieldCreate(wrtGrid, array, name=trim(attrValueSList(1)), 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, coordDim=2, & staggerloc=ESMF_STAGGERLOC_CENTER, array=array, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! write(0,*) 'create gridFB,fieldname=',trim(attrValueSList(1)),trim(attrValueSList(2)), & ! 'lat value=',array(1:5,1),array(1,1:5) field = ESMF_FieldCreate(wrtGrid, array, name=trim(attrValueSList(2)), 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 ! !----------------------------------------------------------------------- !*** SET THE FIRST HISTORY FILE'S TIME INDEX. !----------------------------------------------------------------------- ! wrt_int_state%NFHOUR = 0 ! !----------------------------------------------------------------------- !*** Initialize for POST !----------------------------------------------------------------------- ! call ESMF_LogWrite("before initialize for POST", ESMF_LOGMSG_INFO, rc=rc) if(trim(output_grid) == 'gaussian_grid' .and. wrt_int_state%write_dopost) then do i= 1, wrt_int_state%FBcount call post_getattr_gfs(wrt_int_state,wrt_int_state%wrtFB(i)) enddo endif ! !----------------------------------------------------------------------- !*** Initialize for nemsio file !----------------------------------------------------------------------- ! call ESMF_LogWrite("before initialize for nemsio file", ESMF_LOGMSG_INFO, rc=rc) do i= 1, wrt_int_state%FBcount if (trim(output_grid) == 'gaussian_grid' .and. trim(output_file(i)) == 'nemsio') then ! if (lprnt) write(0,*) 'in wrt initial, befnemsio_first_call wrt_int_state%FBcount=',wrt_int_state%FBcount call nemsio_first_call(wrt_int_state%wrtFB(i), imo, jmo, & wrt_int_state%mype, ntasks, wrt_mpi_comm, & wrt_int_state%FBcount, i, idate, lat, lon, rc) endif enddo call ESMF_LogWrite("after initialize for nemsio file", ESMF_LOGMSG_INFO, rc=rc) ! !----------------------------------------------------------------------- ! IF(RC /= ESMF_SUCCESS) THEN WRITE(0,*)"FAIL: Write_Initialize." ! ELSE ! WRITE(0,*)"PASS: Write_Initialize." ENDIF ! ! write_init_tim = MPI_Wtime() - btim0 ! !----------------------------------------------------------------------- ! end subroutine wrt_initialize ! !----------------------------------------------------------------------- !####################################################################### !----------------------------------------------------------------------- ! 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 type(ESMF_Time) :: currtime type(ESMF_TypeKind_Flag) :: datatype type(ESMF_Field) :: field_work type(ESMF_Grid) :: grid_work, fbgrid, wrtgrid type(ESMF_Array) :: array_work type(ESMF_State),save :: stateGridFB type(optimizeT), save :: optimize(4) type(ESMF_GridComp), save, allocatable :: compsGridFB(:) ! type(write_wrap) :: wrap type(wrt_internal_state),pointer :: wrt_int_state ! integer,dimension(:),allocatable,save :: ih_int, ih_real ! INTEGER,SAVE :: NPOSN_1,NPOSN_2 ! integer :: i,j,n,mype,nolog ! integer :: nf_hours,nf_seconds, nf_minutes, & nseconds,nseconds_num,nseconds_den ! integer :: id integer :: nbdl, idx, date(6), ndig integer :: step=1 ! REAL :: DEGRAD ! logical :: opened logical,save :: first=.true. logical,save :: file_first=.true. ! character(filename_maxstr) :: filename,compname,bundle_name character(40) :: cfhour, cform character(10) :: stepString character(80) :: attrValueS integer :: attrValueI real :: attrValueR real(ESMF_KIND_R8) :: time ! !----------------------------------------------------------------------- ! real(kind=8) :: wait_time, MPI_Wtime real(kind=8) :: times,times2,etim character(10) :: timeb real(kind=8) :: tbeg,tend real(kind=8) :: wbeg,wend integer fieldcount, dimCount real(kind=ESMF_KIND_R8), dimension(:,:,:), pointer :: datar8 real(kind=ESMF_KIND_R8), dimension(:,:), pointer :: datar82d ! integer myattCount logical lprnt ! !----------------------------------------------------------------------- !*********************************************************************** !----------------------------------------------------------------------- ! tbeg = MPI_Wtime() rc = esmf_success ! !----------------------------------------------------------------------- !*** get the current write grid comp name, id, and internal state ! call ESMF_GridCompGet(wrt_comp, name=compname, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! print *,'in wrt run. compname=',trim(compname),' rc=',rc ! instance id from name read(compname(10:11),"(I2)") id ! 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) wrt_int_state%fdate(7) = 1 wrt_int_state%fdate(1:6) = date(1:6) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! if(mype == lead_write_task) print *,'in wrt run, curr time=',date ! 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 ! print *,'in wrt run, io_baseline time=',date ! wrt_int_state%IO_CURRTIMEDIFF = CURRTIME-wrt_int_state%IO_BASETIME ! call ESMF_TimeIntervalGet(timeinterval=wrt_int_state%IO_CURRTIMEDIFF & ,h =nf_hours & !<-- Hours of elapsed time ,m =nf_minutes & !<-- Minutes of elapsed time ,s =nseconds & !<-- Seconds of elapsed time ,sN =nseconds_num & !<-- Numerator of fractional elapsed seconds ,sD =nseconds_den & !<-- denominator of fractional elapsed seconds ,rc =RC) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! if (lprnt) print *,'in wrt run, nf_hours=',nf_hours,nf_minutes,nseconds, & ! 'nseconds_num=',nseconds_num,nseconds_den,'mype=',mype ! nf_seconds = nf_hours*3600+nf_minuteS*60+nseconds+real(nseconds_num)/real(nseconds_den) wrt_int_state%nfhour = nf_seconds/3600. nf_hours = int(nf_seconds/3600.) if(mype == lead_write_task) print *,'in write grid comp, nf_hours=',nf_hours ! if iau_offset > nf_hours, don't write out anything if (nf_hours < 0) return nf_minutes = int((nf_seconds-nf_hours*3600.)/60.) nseconds = int(nf_seconds-nf_hours*3600.-nf_minutes*60.) ! if (nf_seconds-nf_hours*3600 > 0 .and. nsout > 0) then if (nsout > 0) 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,':',nseconds 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, nf_hours=',nf_hours,nf_minutes,nseconds, & 'nseconds_num=',nseconds_num,nseconds_den,' FBCount=',FBCount,' cfhour=',trim(cfhour) ! if(lprnt) print *,'in wrt run, cfhour=',cfhour, & ! print *,'in wrt run, cfhour=',cfhour, & ! ' nf_seconds=',nf_seconds,wrt_int_state%nfhour ! access the time Attribute which is updated by the driver each time call ESMF_LogWrite("before Write component get time", ESMF_LOGMSG_INFO, rc=rc) call ESMF_AttributeGet(imp_state_write, 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_LogWrite("before Write component af get time", ESMF_LOGMSG_INFO, rc=rc) ! !----------------------------------------------------------------------- !*** loop on the files that need to write out !----------------------------------------------------------------------- do i=1, FBCount call ESMF_LogWrite("before Write component get mirror file bundle", ESMF_LOGMSG_INFO, rc=rc) call ESMF_StateGet(imp_state_write, itemName="mirror_"//trim(fcstItemNameList(i)), & fieldbundle=file_bundle, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_LogWrite("before Write component af get mirror file bundle", ESMF_LOGMSG_INFO, rc=rc) !recover fields from cartesian vector and sfc pressure call recover_fields(file_bundle,rc) enddo ! !----------------------------------------------------------------------- !*** do post !----------------------------------------------------------------------- if( wrt_int_state%write_dopost ) then call post_run_gfs(wrt_int_state, mype, wrt_mpi_comm, lead_write_task, & nf_hours, nf_minutes,nseconds) endif ! !----------------------------------------------------------------------- ! ** now loop through output field bundle !----------------------------------------------------------------------- if ( wrt_int_state%output_history ) then file_loop_all: do nbdl=1, wrt_int_state%FBCount ! if(step == 1) then file_bundle = wrt_int_state%wrtFB(nbdl) endif ! 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(nbdl)(1:6) == 'netcdf') then if (ichunk2d == 0) then if( wrt_int_state%mype == 0 ) & ichunk2d = wrt_int_state%lon_end-wrt_int_state%lon_start+1 call mpi_bcast(ichunk2d,1,mpi_integer,0,wrt_mpi_comm,rc) endif if (jchunk2d == 0) then if( wrt_int_state%mype == 0 ) & jchunk2d = wrt_int_state%lat_end-wrt_int_state%lat_start+1 call mpi_bcast(jchunk2d,1,mpi_integer,0,wrt_mpi_comm,rc) endif if (ichunk3d == 0) then if( wrt_int_state%mype == 0 ) & ichunk3d = wrt_int_state%lon_end-wrt_int_state%lon_start+1 call mpi_bcast(ichunk3d,1,mpi_integer,0,wrt_mpi_comm,rc) endif if (jchunk3d == 0) then if( wrt_int_state%mype == 0 ) & jchunk3d = wrt_int_state%lat_end-wrt_int_state%lat_start+1 call mpi_bcast(jchunk3d,1,mpi_integer,0,wrt_mpi_comm,rc) endif if (kchunk3d == 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, rc=rc) endif call mpi_bcast(kchunk3d,1,mpi_integer,0,wrt_mpi_comm,rc) endif if (wrt_int_state%mype == 0) then print *,'ichunk2d,jchunk2d',ichunk2d,jchunk2d print *,'ichunk3d,jchunk3d,kchunk3d',ichunk3d,jchunk3d,kchunk3d endif endif if ( trim(output_file(nbdl)) == 'nemsio' ) then filename = trim(wrt_int_state%wrtFB_names(nbdl))//'f'//trim(cfhour)//'.nemsio' else filename = trim(wrt_int_state%wrtFB_names(nbdl))//'f'//trim(cfhour)//'.nc' endif ! if(mype == lead_write_task) print *,'in wrt run,filename=',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=real(wrt_int_state%nfhour,ESMF_KIND_R8), 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_grid) == 'cubed_sphere_grid') then wbeg = MPI_Wtime() 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 wend = MPI_Wtime() if (lprnt) then write(*,'(A,F10.5,A,I4.2,A,I2.2)')' actual netcdf Write Time is ',wend-wbeg & ,' at Fcst ',NF_HOURS,':',NF_MINUTES endif else if (trim(output_grid) == 'gaussian_grid') then if (trim(output_file(nbdl)) == 'nemsio') then wbeg = MPI_Wtime() call write_nemsio(file_bundle,trim(filename),nf_hours, nf_minutes, & nseconds, nseconds_num, nseconds_den,nbdl, rc) wend = MPI_Wtime() if (lprnt) then write(*,'(A,F10.5,A,I4.2,A,I2.2)')' nemsio Write Time is ',wend-wbeg & ,' at Fcst ',NF_HOURS,':',NF_MINUTES endif else if (trim(output_file(nbdl)) == 'netcdf') then wbeg = MPI_Wtime() call write_netcdf(file_bundle,wrt_int_state%wrtFB(nbdl),trim(filename), & wrt_mpi_comm,wrt_int_state%mype,imo,jmo,& ichunk2d,jchunk2d,ichunk3d,jchunk3d,kchunk3d,rc) wend = MPI_Wtime() if (lprnt) then write(*,'(A,F10.5,A,I4.2,A,I2.2)')' netcdf Write Time is ',wend-wbeg & ,' at Fcst ',NF_HOURS,':',NF_MINUTES endif else if (trim(output_file(nbdl)) == 'netcdf_parallel') then #ifdef NO_PARALLEL_NETCDF rc = ESMF_RC_NOT_IMPL print *,'netcdf_parallel not available on this machine' if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__)) return #endif wbeg = MPI_Wtime() call write_netcdf_parallel(file_bundle,wrt_int_state%wrtFB(nbdl), & trim(filename), wrt_mpi_comm,wrt_int_state%mype,imo,jmo,& ichunk2d,jchunk2d,ichunk3d,jchunk3d,kchunk3d,rc) wend = MPI_Wtime() if (lprnt) then write(*,'(A,F10.5,A,I4.2,A,I2.2)')' parallel netcdf Write Time is ',wend-wbeg & ,' at Fcst ',NF_HOURS,':',NF_MINUTES endif else if (trim(output_file(nbdl)) == 'netcdf_esmf') then wbeg = MPI_Wtime() 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 wend = MPI_Wtime() if (lprnt) then write(*,'(A,F10.5,A,I4.2,A,I2.2)')' netcdf_esmf Write Time is ',wend-wbeg & ,' at Fcst ',NF_HOURS,':',NF_MINUTES endif else ! unknown output_file call ESMF_LogWrite("wrt_run: Unknown output_file",ESMF_LOGMSG_ERROR,rc=RC) call ESMF_Finalize(endflag=ESMF_END_ABORT) endif else if (trim(output_grid) == 'regional_latlon' .or. & trim(output_grid) == 'rotated_latlon' .or. & trim(output_grid) == 'lambert_conformal') then !mask fields according to sfc pressure !if (mype == lead_write_task) print *,'before mask_fields' wbeg = MPI_Wtime() call ESMF_LogWrite("before mask_fields for wrt field bundle", ESMF_LOGMSG_INFO, rc=rc) !call mask_fields(wrt_int_state%wrtFB(nbdl),rc) call mask_fields(file_bundle,rc) !if (mype == lead_write_task) print *,'after mask_fields' call ESMF_LogWrite("after mask_fields for wrt field bundle", ESMF_LOGMSG_INFO, rc=rc) wend = MPI_Wtime() if (mype == lead_write_task) then write(*,'(A,F10.5,A,I4.2,A,I2.2)')' mask_fields time is ',wend-wbeg endif if (trim(output_file(nbdl)) == 'netcdf') then wbeg = MPI_Wtime() call write_netcdf(file_bundle,wrt_int_state%wrtFB(nbdl),trim(filename), & wrt_mpi_comm,wrt_int_state%mype,imo,jmo,& ichunk2d,jchunk2d,ichunk3d,jchunk3d,kchunk3d,rc) wend = MPI_Wtime() if (mype == lead_write_task) then write(*,'(A,F10.5,A,I4.2,A,I2.2)')' netcdf Write Time is ',wend-wbeg & ,' at Fcst ',NF_HOURS,':',NF_MINUTES endif else ! unknown output_file call ESMF_LogWrite("wrt_run: Unknown output_file",ESMF_LOGMSG_ERROR,rc=RC) call ESMF_Finalize(endflag=ESMF_END_ABORT) 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 do n=701,900 inquire(n,opened=OPENED) if(.not.opened)then nolog = n exit endif enddo ! open(nolog,file='logf'//trim(cfhour),form='FORMATTED') write(nolog,100)wrt_int_state%nfhour,idate(1:6) 100 format(' completed fv3gfs fhour=',f10.3,2x,6(i4,2x)) close(nolog) endif ! ! !----------------------------------------------------------------------- ! call ESMF_VMBarrier(VM, rc=rc) ! 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 ! IF(RC /= ESMF_SUCCESS) THEN WRITE(0,*)"FAIL: WRITE_RUN" ! ELSE ! WRITE(0,*)"PASS: WRITE_RUN" 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) deallocate(wrap%write_int_state,stat=stat) ! if (ESMF_LogFoundDeallocError(statusToCheck=stat, & msg="Deallocation of internal state memory failed.", & line=__LINE__, file=__FILE__)) return ! !----------------------------------------------------------------------- ! IF(RC /= ESMF_SUCCESS)THEN WRITE(0,*)'FAIL: Write_Finalize.' ! ELSE ! WRITE(0,*)'PASS: Write_Finalize.' ENDIF ! !----------------------------------------------------------------------- ! end subroutine wrt_finalize ! !----------------------------------------------------------------------- ! subroutine recover_fields(file_bundle,rc) type(ESMF_FieldBundle), intent(in) :: file_bundle integer, intent(out), optional :: rc ! integer i,j,k,ifld,fieldCount,nstt,nend,fieldDimCount,gridDimCount integer istart,iend,jstart,jend,kstart,kend,km 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), dimension(:,:,:,:), pointer :: cart3dPtr3dr8 save lonloc, latloc 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 if( first_getlatlon ) then 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 first_getlatlon = .false. endif ! 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 print *,'ERROR, 2D the vector dimension /= 3, rc=',rc exit 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,nstt,nend,fieldDimCount,gridDimCount integer istart,iend,jstart,jend,kstart,kend,km 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 save maskwrt 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 if( first_getmaskwrt ) then 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 first_getmaskwrt = .false. endif !first_getmaskwrt ! 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 print *,'ERROR, 3D the vector dimension /= 3, rc=',rc exit 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 print *,'ERROR, 2D the vector dimension /= 3, rc=',rc exit 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(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 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 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_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 ! 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,d1=1.d0 integer :: jhe,jho,j0=0 ! 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,js,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,d1=1.d0 integer(4) :: jhe,jho,j0=0 ! 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,js,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) lat/lon to grid ! inv == -1 (x,y) ---> (glon,glat) grid to lat/lon real(ESMF_KIND_R8) :: en,f,rho,rho0, dlon, theta, xp, yp 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)) ENDIF f=(cos(stlat1*dtor)*tan((45+0.5*stlat1)*dtor)**en)/en rho0=a*f/(tan((45+0.5*c_lat)*dtor)**en) if (inv == 1) then ! FORWARD TRANSFORMATION rho=a*f/(tan((45+0.5*glat)*dtor)**en) dlon=modulo(glon-c_lon+180+3600,360.)-180.D0 theta=en*dlon*dtor x=rho*sin(theta) y=rho0-rho*cos(theta) else if (inv == -1) then ! INVERSE TRANSFORMATION y=rho0-y rho = sqrt(x*x+y*y) theta=atan2(x,y) glon=c_lon+(theta/en)*rtod glon=modulo(glon+180+3600,360.)-180.D0 ! glat=(2.0*atan((a*f/rho)**(1.0/en))-0.5*pi)*rtod glat=(0.5*pi-2.0*atan((rho/(a*f))**(1.0/en)))*rtod 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,idx 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 ! !-----------------------------------------------------------------------