diff --git a/fv3_cap.F90 b/fv3_cap.F90 index e8e482099..d04abdce2 100644 --- a/fv3_cap.F90 +++ b/fv3_cap.F90 @@ -32,7 +32,7 @@ module fv3gfs_cap_mod cplprint_flag,output_1st_tstep_rst, & first_kdt - use module_fv3_io_def, only: num_pes_fcst,write_groups, & + use module_fv3_io_def, only: num_pes_fcst,write_groups,app_domain, & num_files, filename_base, & wrttasks_per_group, n_group, & lead_wrttask, last_wrttask, & @@ -322,6 +322,10 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) label ='write_tasks_per_group:',rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + CALL ESMF_ConfigGetAttribute(config=CF,value=app_domain, default="global", & + label ='app_domain:',rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + if(mype == 0) print *,'af nems config,restart_interval=',restart_interval, & 'quilting=',quilting,'write_groups=',write_groups,wrttasks_per_group, & 'calendar=',trim(calendar),'calendar_type=',calendar_type @@ -692,6 +696,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) isrctermprocessing = 1 call ESMF_FieldBundleRegridStore(fcstFB(j), wrtFB(j,i), & regridMethod=regridmethod, routehandle=routehandle(j,i), & + unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, & srcTermProcessing=isrctermprocessing, rc=rc) ! if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return diff --git a/io/FV3GFS_io.F90 b/io/FV3GFS_io.F90 index e7f7dbd57..4b2d1e426 100644 --- a/io/FV3GFS_io.F90 +++ b/io/FV3GFS_io.F90 @@ -101,8 +101,7 @@ module FV3GFS_io_mod logical :: uwork_set = .false. character(128) :: uwindname integer, parameter, public :: DIAG_SIZE = 500 -! real(kind=kind_phys), parameter :: missing_value = 1.d30 - real(kind=kind_phys), parameter :: missing_value = 9.99e20 + real, parameter :: missing_value = 9.99e20 real, parameter:: stndrd_atmos_ps = 101325. real, parameter:: stndrd_atmos_lapse = 0.0065 diff --git a/io/module_fv3_io_def.F90 b/io/module_fv3_io_def.F90 index d039ccc73..17e9f308d 100644 --- a/io/module_fv3_io_def.F90 +++ b/io/module_fv3_io_def.F90 @@ -15,6 +15,7 @@ module module_fv3_io_def logical :: write_nemsioflip logical :: write_fsyncflag integer :: num_files + character(255) :: app_domain character(255) :: output_grid character(255) :: output_file integer :: imo,jmo diff --git a/io/module_wrt_grid_comp.F90 b/io/module_wrt_grid_comp.F90 index 405af2841..aad3991d0 100644 --- a/io/module_wrt_grid_comp.F90 +++ b/io/module_wrt_grid_comp.F90 @@ -32,7 +32,7 @@ module module_wrt_grid_comp use esmf use write_internal_state use module_fv3_io_def, only : num_pes_fcst,lead_wrttask, last_wrttask, & - n_group, num_files, & + n_group, num_files, app_domain, & filename_base, output_grid, output_file, & imo, jmo, write_nemsioflip, & nsout => nsout_io, & @@ -71,6 +71,7 @@ module module_wrt_grid_comp 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. ! !----------------------------------------------------------------------- @@ -147,10 +148,12 @@ subroutine wrt_initialize(wrt_comp, imp_state_write, exp_state_write, clock, rc) 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 @@ -163,6 +166,7 @@ subroutine wrt_initialize(wrt_comp, imp_state_write, exp_state_write, clock, rc) 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(:) @@ -196,6 +200,7 @@ subroutine wrt_initialize(wrt_comp, imp_state_write, exp_state_write, clock, rc) logical :: lprnt !test integer myattCount + real(ESMF_KIND_R8),dimension(:,:), pointer :: glatPtr, glonPtr ! !----------------------------------------------------------------------- !*********************************************************************** @@ -279,24 +284,47 @@ subroutine wrt_initialize(wrt_comp, imp_state_write, exp_state_write, clock, rc) if ( trim(output_grid) == 'cubed_sphere_grid' ) then mytile = mod(wrt_int_state%mype,ntasks)+1 - 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 ( 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 - 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) + 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 - wrtgrid = ESMF_GridCreateMosaic(filename="INPUT/"//trim(gridfile), & - regDecompPTile=decomptile,tileFilePath="INPUT/", & - staggerlocList=(/ESMF_STAGGERLOC_CENTER, ESMF_STAGGERLOC_CORNER/), & - name='wrt_grid', rc=rc) + ! 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 - 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/), & @@ -1440,6 +1468,19 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) 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) == 'netcdf') then wbeg = MPI_Wtime() @@ -1749,6 +1790,269 @@ subroutine recover_fields(file_bundle,rc) 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, & diff --git a/io/post_gfs.F90 b/io/post_gfs.F90 index 3ea19c555..9f3a0a80f 100644 --- a/io/post_gfs.F90 +++ b/io/post_gfs.F90 @@ -133,6 +133,7 @@ subroutine post_run_gfs(wrt_int_state,mypei,mpicomp,lead_write, & ! ifhr = mynfhr ifmin = mynfmin + if (ifhr == 0 ) ifmin = 0 if(mype==0) print *,'bf set_postvars,ifmin=',ifmin,'ifhr=',ifhr call set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile, & setvar_sfcfile)