From 7ffe6471c20404091fbbf8f321fbb9ee84a4f36d Mon Sep 17 00:00:00 2001 From: Dusan Jovic <48258889+DusanJovic-NOAA@users.noreply.github.com> Date: Fri, 24 Jan 2020 17:47:02 -0500 Subject: [PATCH] GFSv16 netcdf post ficein cpl (#48) * fv3atm issue #37: fix the real(8) lat/lon in netcdf file * fv3atm #35: Reducing background vertical diffusivities in the inversion layers * fv3atm #24: bug in gfsphysics/physics/moninedmf_hafs.f * fv3atm #18: Optimize netcdf write component and bugfix for post and samfdeepcnv.f * set (0-1) bounds for ficein_cpl * remove cache_size due to lower netcdf verion 4.5.1 on mars * Change ice falling to 0.9 in gfsphysics/physics/gfdl_cloud_microphys.F90 --- .gitmodules | 2 + atmos_model.F90 | 2 +- ccpp/framework | 2 +- ccpp/physics | 2 +- gfsphysics/physics/gfdl_cloud_microphys.F90 | 2 +- .../physics/module_sf_noahmp_glacier.f90 | 0 gfsphysics/physics/module_sf_noahmplsm.f90 | 0 gfsphysics/physics/module_wrf_utl.f90 | 0 gfsphysics/physics/moninedmf_hafs.f | 6 +- gfsphysics/physics/noahmp_tables.f90 | 0 gfsphysics/physics/samfdeepcnv.f | 28 ++--- gfsphysics/physics/satmedmfvdifq.f | 11 +- gfsphysics/physics/sfc_noahmp_drv.f | 0 io/FV3GFS_io.F90 | 2 +- io/module_write_nemsio.F90 | 2 +- io/module_write_netcdf.F90 | 115 ++++++++++-------- io/module_wrt_grid_comp.F90 | 64 ++++++---- io/post_gfs.F90 | 14 ++- 18 files changed, 148 insertions(+), 104 deletions(-) mode change 100755 => 100644 gfsphysics/physics/module_sf_noahmp_glacier.f90 mode change 100755 => 100644 gfsphysics/physics/module_sf_noahmplsm.f90 mode change 100755 => 100644 gfsphysics/physics/module_wrf_utl.f90 mode change 100755 => 100644 gfsphysics/physics/noahmp_tables.f90 mode change 100755 => 100644 gfsphysics/physics/sfc_noahmp_drv.f diff --git a/.gitmodules b/.gitmodules index 949d298df..d253f6966 100644 --- a/.gitmodules +++ b/.gitmodules @@ -5,6 +5,8 @@ [submodule "ccpp/framework"] path = ccpp/framework url = https://github.com/NCAR/ccpp-framework + branch = master [submodule "ccpp/physics"] path = ccpp/physics url = https://github.com/NCAR/ccpp-physics + branch = master diff --git a/atmos_model.F90 b/atmos_model.F90 index 34e1a6c99..9a73e0151 100644 --- a/atmos_model.F90 +++ b/atmos_model.F90 @@ -1720,7 +1720,7 @@ subroutine assign_importdata(rc) IPD_Data(nb)%Coupling%ficein_cpl(ix) = zero if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) > zero) then if (datar8(i,j) >= IPD_control%min_seaice*IPD_Data(nb)%Sfcprop%oceanfrac(ix)) then - IPD_Data(nb)%Coupling%ficein_cpl(ix) = datar8(i,j) + IPD_Data(nb)%Coupling%ficein_cpl(ix) = max(zero, min(datar8(i,j),one)) ! if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) == one) IPD_Data(nb)%Sfcprop%slmsk(ix) = 2. !slmsk=2 crashes in gcycle on partial land points IPD_Data(nb)%Sfcprop%slmsk(ix) = 2. !slmsk=2 crashes in gcycle on partial land points IPD_Data(nb)%Coupling%slimskin_cpl(ix) = 4. diff --git a/ccpp/framework b/ccpp/framework index 7ab419eee..e77210986 160000 --- a/ccpp/framework +++ b/ccpp/framework @@ -1 +1 @@ -Subproject commit 7ab419eeebe133e706d9825d14c5bdc5d190e60d +Subproject commit e7721098639ee73c2a69ee0e8423e8905549e240 diff --git a/ccpp/physics b/ccpp/physics index 39981891a..01ed01fb0 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit 39981891a5b84e0f08acdf4d9e72d31c0c07aa82 +Subproject commit 01ed01fb0b3112e96eb619e0339d88fb0201982f diff --git a/gfsphysics/physics/gfdl_cloud_microphys.F90 b/gfsphysics/physics/gfdl_cloud_microphys.F90 index ba4c814d6..f01486db0 100644 --- a/gfsphysics/physics/gfdl_cloud_microphys.F90 +++ b/gfsphysics/physics/gfdl_cloud_microphys.F90 @@ -3266,7 +3266,7 @@ subroutine fall_speed (ktop, kbot, den, qs, qi, qg, ql, tk, vts, vti, vtg) else tc (k) = tk (k) - tice vti (k) = (3. + log10 (qi (k) * den (k))) * (tc (k) * (aa * tc (k) + bb) + cc) + dd * tc (k) + ee - vti (k) = vi0 * exp (log_10 * vti (k)) * 0.8 + vti (k) = vi0 * exp (log_10 * vti (k)) * 0.9 vti (k) = min (vi_max, max (vf_min, vti (k))) endif enddo diff --git a/gfsphysics/physics/module_sf_noahmp_glacier.f90 b/gfsphysics/physics/module_sf_noahmp_glacier.f90 old mode 100755 new mode 100644 diff --git a/gfsphysics/physics/module_sf_noahmplsm.f90 b/gfsphysics/physics/module_sf_noahmplsm.f90 old mode 100755 new mode 100644 diff --git a/gfsphysics/physics/module_wrf_utl.f90 b/gfsphysics/physics/module_wrf_utl.f90 old mode 100755 new mode 100644 diff --git a/gfsphysics/physics/moninedmf_hafs.f b/gfsphysics/physics/moninedmf_hafs.f index 42de545f2..a035ad3d6 100644 --- a/gfsphysics/physics/moninedmf_hafs.f +++ b/gfsphysics/physics/moninedmf_hafs.f @@ -1360,7 +1360,11 @@ subroutine moninedmf_hafs(ix,im,km,ntrac,ntcw,dv,du,tau,rtg, & tem = 0.5 * (diss(i,k-1)+diss(i,k)) tem = max(tem, 0.) ttend = tem / cp - tau(i,k) = tau(i,k) + 0.5*ttend + if (alpha .gt. 0.0) then + tau(i,k) = tau(i,k) + 0.5*ttend + else + tau(i,k) = tau(i,k) + 0.7*ttend ! in HWRF/HMON, use 0.7 + endif enddo enddo ! diff --git a/gfsphysics/physics/noahmp_tables.f90 b/gfsphysics/physics/noahmp_tables.f90 old mode 100755 new mode 100644 diff --git a/gfsphysics/physics/samfdeepcnv.f b/gfsphysics/physics/samfdeepcnv.f index 76204ebb4..d29410d2f 100644 --- a/gfsphysics/physics/samfdeepcnv.f +++ b/gfsphysics/physics/samfdeepcnv.f @@ -1547,22 +1547,22 @@ subroutine samfdeepcnv(im,ix,km,delt,itc,ntc,ntk,ntr,delp, enddo enddo do i = 1, im - betamn = betas - if(islimsk(i) == 1) betamn = betal - if(ntk > 0) then - betamx = betamn + dbeta - if(tkemean(i) > tkemx) then - beta = betamn - else if(tkemean(i) < tkemn) then - beta = betamx + if(cnvflg(i)) then + betamn = betas + if(islimsk(i) == 1) betamn = betal + if(ntk > 0) then + betamx = betamn + dbeta + if(tkemean(i) > tkemx) then + beta = betamn + else if(tkemean(i) < tkemn) then + beta = betamx + else + tem = (betamx - betamn) * (tkemean(i) - tkemn) + beta = betamx - tem / dtke + endif else - tem = (betamx - betamn) * (tkemean(i) - tkemn) - beta = betamx - tem / dtke + beta = betamn endif - else - beta = betamn - endif - if(cnvflg(i)) then dz = (sumx(i)+zi(i,1))/float(kbcon(i)) tem = 1./float(kbcon(i)) xlamd(i) = (1.-beta**tem)/dz diff --git a/gfsphysics/physics/satmedmfvdifq.f b/gfsphysics/physics/satmedmfvdifq.f index 11c047fd0..1cc0bbe89 100644 --- a/gfsphysics/physics/satmedmfvdifq.f +++ b/gfsphysics/physics/satmedmfvdifq.f @@ -148,7 +148,8 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke, & epsi, beta, chx, cqx, & rdt, rdz, qmin, qlmin, & rimin, rbcr, rbint, tdzmin, - & rlmn, rlmn1, rlmx, elmx, + & rlmn, rlmn1, rlmn2, + & rlmx, elmx, & ttend, utend, vtend, qtend, & zfac, zfmin, vk, spdk2, & tkmin, tkminx, xkzinv, xkgdx, @@ -172,7 +173,8 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke, parameter(gamcrt=3.,gamcrq=0.,sfcfrac=0.1) parameter(vk=0.4,rimin=-100.) parameter(rbcr=0.25,zolcru=-0.02,tdzmin=1.e-3) - parameter(rlmn=30.,rlmn1=5.,rlmx=300.,elmx=300.) + parameter(rlmn=30.,rlmn1=5.,rlmn2=10.) + parameter(rlmx=300.,elmx=300.) parameter(prmin=0.25,prmax=4.0) parameter(pr0=1.0,prtke=1.0,prscu=0.67) parameter(f0=1.e-4,crbmin=0.15,crbmax=0.35) @@ -698,8 +700,9 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke, ! if(tem1 > 1.e-5) then tem1 = tvx(i,k+1)-tvx(i,k) if(tem1 > 0.) then - xkzo(i,k) = min(xkzo(i,k),xkzinv) - xkzmo(i,k) = min(xkzmo(i,k),xkzinv) + xkzo(i,k) = min(xkzo(i,k), xkzinv) + xkzmo(i,k) = min(xkzmo(i,k), xkzinv) + rlmnz(i,k) = min(rlmnz(i,k), rlmn2) endif enddo enddo diff --git a/gfsphysics/physics/sfc_noahmp_drv.f b/gfsphysics/physics/sfc_noahmp_drv.f old mode 100755 new mode 100644 diff --git a/io/FV3GFS_io.F90 b/io/FV3GFS_io.F90 index 1ac6eae7b..972d43fde 100644 --- a/io/FV3GFS_io.F90 +++ b/io/FV3GFS_io.F90 @@ -285,7 +285,7 @@ subroutine FV3GFS_IPD_checksum (Model, IPD_Data, Atm_block) temp2d(i,j,55) = IPD_Data(nb)%Coupling%visbmui(ix) temp2d(i,j,56) = IPD_Data(nb)%Coupling%visdfui(ix) temp2d(i,j,57) = IPD_Data(nb)%Coupling%sfcdsw(ix) - temp2d(i,j,59) = IPD_Data(nb)%Coupling%sfcnsw(ix) + temp2d(i,j,58) = IPD_Data(nb)%Coupling%sfcnsw(ix) temp2d(i,j,59) = IPD_Data(nb)%Coupling%sfcdlw(ix) temp2d(i,j,60) = IPD_Data(nb)%Grid%xlon(ix) temp2d(i,j,61) = IPD_Data(nb)%Grid%xlat(ix) diff --git a/io/module_write_nemsio.F90 b/io/module_write_nemsio.F90 index 3afd66789..e51f64a52 100644 --- a/io/module_write_nemsio.F90 +++ b/io/module_write_nemsio.F90 @@ -51,7 +51,7 @@ subroutine nemsio_first_call(fieldbundle, imo, jmo, & integer, intent(in) :: wrt_mype, wrt_ntasks, wrt_mpi_comm integer, intent(in) :: wrt_nbdl, mybdl integer, intent(in) :: inidate(7) - real, intent(in) :: lat(:), lon(:) + real(8), intent(in) :: lat(:), lon(:) integer, optional,intent(out) :: rc !** local vars diff --git a/io/module_write_netcdf.F90 b/io/module_write_netcdf.F90 index 1fce3d8b9..fc59c75c9 100644 --- a/io/module_write_netcdf.F90 +++ b/io/module_write_netcdf.F90 @@ -42,7 +42,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc real(4), dimension(:,:,:), allocatable :: arrayr4_3d,arrayr4_3d_save real(8), dimension(:,:,:), allocatable :: arrayr8_3d - real(8) rad2dg,x(im),y(jm) + real(8) x(im),y(jm) integer :: fieldCount, fieldDimCount, gridDimCount integer, dimension(:), allocatable :: ungriddedLBound, ungriddedUBound @@ -56,7 +56,8 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc character(len=ESMF_MAXSTR) :: attName, fldName integer :: varival - real(4) :: varr4val, scale_fact, compress_err, offset, dataMin, dataMax + real(4) :: varr4val, scale_fact, offset, dataMin, dataMax + real(4), allocatable, dimension(:) :: compress_err real(8) :: varr8val character(len=ESMF_MAXSTR) :: varcval @@ -71,10 +72,10 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc ! !! ! - rad2dg = 45./atan(1.0) call ESMF_FieldBundleGet(fieldbundle, fieldCount=fieldCount, rc=rc); ESMF_ERR_RETURN(rc) + allocate(compress_err(fieldCount)); compress_err=-999. allocate(fldlev(fieldCount)) ; fldlev = 0 allocate(fcstField(fieldCount)) allocate(varids(fieldCount)) @@ -117,13 +118,13 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc if (mype==0) then if (ideflate == 0) then - ncerr = nf90_create(trim(filename), cmode=IOR(NF90_CLOBBER,NF90_64BIT_OFFSET), & + ncerr = nf90_create(trim(filename), cmode=IOR(IOR(NF90_CLOBBER,NF90_64BIT_OFFSET),NF90_SHARE), & ncid=ncid); NC_ERR_STOP(ncerr) ncerr = nf90_set_fill(ncid, NF90_NOFILL, oldMode); NC_ERR_STOP(ncerr) else ncerr = nf90_create(trim(filename), cmode=IOR(IOR(NF90_CLOBBER,NF90_NETCDF4),NF90_CLASSIC_MODEL), & ncid=ncid); NC_ERR_STOP(ncerr) - ! if compression on use HDF5 format with default _FillValue + ncerr = nf90_set_fill(ncid, NF90_NOFILL, oldMode); NC_ERR_STOP(ncerr) endif ! define dimensions @@ -156,28 +157,32 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc ! define variables if (fldlev(i) == 1) then if (typekind == ESMF_TYPEKIND_R4) then - ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & - (/im_dimid,jm_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr) if (ideflate > 0) then - ! shuffle filter on for lossless compression - ncerr = nf90_def_var_deflate(ncid, varids(i), 1, 1, ideflate) - NC_ERR_STOP(ncerr) + ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & + (/im_dimid,jm_dimid,time_dimid/), varids(i), & + shuffle=.true.,deflate_level=ideflate, & + chunksizes=(/im,jm,1/)); NC_ERR_STOP(ncerr) + else + ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & + (/im_dimid,jm_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr) endif else if (typekind == ESMF_TYPEKIND_R8) then ncerr = nf90_def_var(ncid, trim(fldName), NF90_DOUBLE, & - (/im_dimid,jm_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr) + (/im_dimid,jm_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr) else write(0,*)'Unsupported typekind ', typekind stop end if else if (fldlev(i) > 1) then if (typekind == ESMF_TYPEKIND_R4) then - ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & - (/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr) if (ideflate > 0) then - ! shuffle filter off since lossy compression used - ncerr = nf90_def_var_deflate(ncid, varids(i), 0, 1, ideflate) - NC_ERR_STOP(ncerr) + ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & + (/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i), & + shuffle=.false.,deflate_level=ideflate, & + chunksizes=(/im,jm,1,1/)); NC_ERR_STOP(ncerr) + else + ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & + (/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr) endif else if (typekind == ESMF_TYPEKIND_R8) then ncerr = nf90_def_var(ncid, trim(fldName), NF90_DOUBLE, & @@ -219,8 +224,8 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc call ESMF_AttributeGet(fcstField(i), convention="NetCDF", purpose="FV3", & name=trim(attName), value=varr8val, & rc=rc); ESMF_ERR_RETURN(rc) - if (trim(attName) /= '_FillValue' .or. ideflate == 0) then - ! FIXME: _FillValue must be cast to var type when using NF90_NETCDF4 + if (trim(attName) /= '_FillValue' ) then + ! FIXME: _FillValue must be cast to var type for recent versions of netcdf ncerr = nf90_put_att(ncid, varids(i), trim(attName), varr8val); NC_ERR_STOP(ncerr) endif @@ -236,6 +241,25 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc end do ! i=1,fieldCount +! write grid_xt, grid_yt attributes + if (trim(output_grid) == 'gaussian_grid' .or. & + trim(output_grid) == 'regional_latlon') then + ncerr = nf90_put_att(ncid, im_varid, "long_name", "T-cell longitude"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, im_varid, "units", "degrees_E"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, jm_varid, "long_name", "T-cell latiitude"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees_N"); NC_ERR_STOP(ncerr) + else if (trim(output_grid) == 'rotated_latlon') then + ncerr = nf90_put_att(ncid, im_varid, "long_name", "rotated T-cell longiitude"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, im_varid, "units", "degrees"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, jm_varid, "long_name", "rotated T-cell latiitude"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees"); NC_ERR_STOP(ncerr) + else if (trim(output_grid) == 'lambert_conformal') then + ncerr = nf90_put_att(ncid, im_varid, "long_name", "x-coordinate of projection"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, im_varid, "units", "meters"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, jm_varid, "long_name", "y-coordinate of projection"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, jm_varid, "units", "meters"); NC_ERR_STOP(ncerr) + endif + ncerr = nf90_enddef(ncid); NC_ERR_STOP(ncerr) end if @@ -247,31 +271,19 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc if (mype==0) then if (trim(output_grid) == 'gaussian_grid' .or. & trim(output_grid) == 'regional_latlon') then - ncerr = nf90_put_var(ncid, im_varid, values=rad2dg*arrayr8(:,1) ); NC_ERR_STOP(ncerr) - ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, im_varid, "long_name", "T-cell longitude"); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, im_varid, "units", "degrees_E"); NC_ERR_STOP(ncerr) - ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_var(ncid, im_varid, values=arrayr8(:,1) ); NC_ERR_STOP(ncerr) else if (trim(output_grid) == 'rotated_latlon') then do i=1,im x(i) = lon1 + (lon2-lon1)/(im-1) * (i-1) enddo ncerr = nf90_put_var(ncid, im_varid, values=x ); NC_ERR_STOP(ncerr) - ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, im_varid, "long_name", "rotated T-cell longiitude"); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, im_varid, "units", "degrees"); NC_ERR_STOP(ncerr) - ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) else if (trim(output_grid) == 'lambert_conformal') then do i=1,im x(i) = dx * (i-1) enddo ncerr = nf90_put_var(ncid, im_varid, values=x ); NC_ERR_STOP(ncerr) - ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, im_varid, "long_name", "x-coordinate of projection"); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, im_varid, "units", "meters"); NC_ERR_STOP(ncerr) - ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) endif - ncerr = nf90_put_var(ncid, lon_varid, values=rad2dg*arrayr8 ); NC_ERR_STOP(ncerr) + ncerr = nf90_put_var(ncid, lon_varid, values=arrayr8 ); NC_ERR_STOP(ncerr) endif call ESMF_GridGetCoord(wrtGrid, coordDim=2, array=array, rc=rc); ESMF_ERR_RETURN(rc) @@ -279,31 +291,19 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc if (mype==0) then if (trim(output_grid) == 'gaussian_grid' .or. & trim(output_grid) == 'regional_latlon') then - ncerr = nf90_put_var(ncid, jm_varid, values=rad2dg*arrayr8(1,:) ); NC_ERR_STOP(ncerr) - ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, jm_varid, "long_name", "T-cell latiitude"); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees_N"); NC_ERR_STOP(ncerr) - ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_var(ncid, jm_varid, values=arrayr8(1,:) ); NC_ERR_STOP(ncerr) else if (trim(output_grid) == 'rotated_latlon') then do j=1,jm y(j) = lat1 + (lat2-lat1)/(jm-1) * (j-1) enddo ncerr = nf90_put_var(ncid, jm_varid, values=y ); NC_ERR_STOP(ncerr) - ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, jm_varid, "long_name", "rotated T-cell latiitude"); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees"); NC_ERR_STOP(ncerr) - ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) else if (trim(output_grid) == 'lambert_conformal') then do j=1,jm y(j) = dy * (j-1) enddo ncerr = nf90_put_var(ncid, jm_varid, values=y ); NC_ERR_STOP(ncerr) - ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, jm_varid, "long_name", "y-coordinate of projection"); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, jm_varid, "units", "meters"); NC_ERR_STOP(ncerr) - ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) endif - ncerr = nf90_put_var(ncid, lat_varid, values=rad2dg*arrayr8 ); NC_ERR_STOP(ncerr) + ncerr = nf90_put_var(ncid, lat_varid, values=arrayr8 ); NC_ERR_STOP(ncerr) endif do i=1, fieldCount @@ -344,11 +344,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc arrayr4_3d = quantized(arrayr4_3d_save, nbits, dataMin, dataMax) ! compute max abs compression error, save as a variable ! attribute. - compress_err = maxval(abs(arrayr4_3d_save-arrayr4_3d)) - ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, varids(i), 'max_abs_compression_error', compress_err); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, varids(i), 'nbits', nbits); NC_ERR_STOP(ncerr) - ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) + compress_err(i) = maxval(abs(arrayr4_3d_save-arrayr4_3d)) endif ncerr = nf90_put_var(ncid, varids(i), values=arrayr4_3d, start=(/1,1,1/),count=(/im,jm,lm,1/) ); NC_ERR_STOP(ncerr) end if @@ -363,6 +359,17 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc end do + if (ideflate > 0 .and. nbits > 0 .and. mype == 0) then + ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) + do i=1, fieldCount + if (compress_err(i) > 0) then + ncerr = nf90_put_att(ncid, varids(i), 'max_abs_compression_error', compress_err(i)); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, varids(i), 'nbits', nbits); NC_ERR_STOP(ncerr) + endif + enddo + ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) + endif + deallocate(arrayr4) deallocate(arrayr8) deallocate(arrayr4_3d,arrayr4_3d_save) @@ -484,9 +491,9 @@ subroutine get_grid_attr(grid, prefix, ncid, varid, rc) else if (typekind==ESMF_TYPEKIND_R8) then call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & name=trim(attName), value=varr8val, rc=rc); ESMF_ERR_RETURN(rc) - if (trim(attName) /= '_FillValue' .or. ideflate == 0) then - ! FIXME: _FillValue must be cast to var type when using - ! NF90_NETCDF4. Until this is fixed, using netCDF default _FillValue. + if (trim(attName) /= '_FillValue') then + ! FIXME: _FillValue must be cast to var type for recent versions + ! of netcdf ncerr = nf90_put_att(ncid, varid, trim(attName(ind+1:len(attName))), varr8val); NC_ERR_STOP(ncerr) endif diff --git a/io/module_wrt_grid_comp.F90 b/io/module_wrt_grid_comp.F90 index aad3991d0..84769eaea 100644 --- a/io/module_wrt_grid_comp.F90 +++ b/io/module_wrt_grid_comp.F90 @@ -183,7 +183,7 @@ subroutine wrt_initialize(wrt_comp, imp_state_write, exp_state_write, clock, rc) character(128) :: FBlist_outfilename(100), outfile_name character(128),dimension(:,:), allocatable :: outfilename real(8), dimension(:), allocatable :: slat - real, dimension(:), allocatable :: lat, lon, axesdata + 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 @@ -358,19 +358,20 @@ subroutine wrt_initialize(wrt_comp, imp_state_write, exp_state_write, clock, rc) wrt_int_state%latstart = lat(1) wrt_int_state%latlast = lat(jmo) do j=1,imo - lon(j) = 360./real(imo) *real(j-1) + 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./real(imo) * (i-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 -! print *,'aft wrtgrd, lon=',lonPtr(lbound(lonPtr,1),lbound(lonPtr,2)), & +! 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) @@ -1622,13 +1623,14 @@ subroutine recover_fields(file_bundle,rc) 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 lon, lat + save lonloc, latloc real(ESMF_KIND_R8) :: coslon, sinlon, sinlat ! ! get filed count @@ -1648,9 +1650,18 @@ subroutine recover_fields(file_bundle,rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return - lon = lon * pi/180. -! print *,'in 3DCartesian2wind, lon dim=',lbound(lon,1),ubound(lon,1),lbound(lon,2),ubound(lon,2), & -! 'lon=',lon(lbound(lon,1),lbound(lon,2)), lon(ubound(lon,1),ubound(lon,2)) + 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) @@ -1658,9 +1669,18 @@ subroutine recover_fields(file_bundle,rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return - lat = lat * pi/180. -! print *,'in 3DCartesian2wind, lat dim=',lbound(lat,1),ubound(lat,1),lbound(lat,2),ubound(lat,2), & -! 'lat=',lat(lbound(lon,1),lbound(lon,2)), lat(ubound(lon,1),ubound(lon,2)) + 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 ! @@ -1718,18 +1738,18 @@ subroutine recover_fields(file_bundle,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,lon,lat,cart3dPtr3dr4,jstart,jend,istart,iend,k) & -!!$omp private(i,j,coslon,sinlon,sinlat) +!$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(lon(i,j)) - sinlon = sin(lon(i,j)) - sinlat = sin(lat(i,j)) + 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(lat(i,j)) + + cart3dPtr3dr4(3,i,j,k) * cos(latloc(i,j)) enddo enddo enddo @@ -1749,18 +1769,18 @@ subroutine recover_fields(file_bundle,rc) 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,lon,lat,cart3dPtr2dr4,jstart,jend,istart,iend) & +!$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(lon(i,j)) - sinlon = sin(lon(i,j)) - sinlat = sin(lat(i,j)) + 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(lat(i,j)) + + cart3dPtr2dr4(3,i,j) * cos(latloc(i,j)) enddo enddo endif diff --git a/io/post_gfs.F90 b/io/post_gfs.F90 index 98d4fef50..daa6551bd 100644 --- a/io/post_gfs.F90 +++ b/io/post_gfs.F90 @@ -1146,7 +1146,7 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile, & do i=ista, iend if (arrayr42d(i,j) /= spval) then !set range within (0,1) - sr(i,j) = min(1.,max(0.,sr(i,j))) + sr(i,j) = min(1.,max(0.,arrayr42d(i,j))) else sr(i,j) = spval endif @@ -2295,7 +2295,6 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile, & do j=jsta,jend do i=1,im omga(i,j,l) = (-1.) * wh(i,j,l) * dpres(i,j,l)/zint(i,j,l) - pmid(i,j,l) = rgas*dpres(i,j,l) * t(i,j,l)*(q(i,j,l)*fv+1.0)/grav/zint(i,j,l) zint(i,j,l) = zint(i,j,l) + zint(i,j,l+1) enddo enddo @@ -2318,6 +2317,15 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile, & enddo end do +!compute pmid from averaged two layer pint + do l=lm,1,-1 + do j=jsta,jend + do i=1,im + pmid(i,j,l) = 0.5*(pint(i,j,l)+pint(i,j,l+1)) + enddo + enddo + enddo + !$omp parallel do private(i,j) do j=jsta,jend do i=1,im @@ -2337,7 +2345,7 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile, & end do end do -! compute zmid ??? where is definition of alpint(1) +! compute zmid do l=lm,1,-1 !$omp parallel do private(i,j) do j=jsta,jend