diff --git a/atmos_model.F90 b/atmos_model.F90 index 987608232e..0c7d273377 100644 --- a/atmos_model.F90 +++ b/atmos_model.F90 @@ -1724,7 +1724,6 @@ subroutine assign_importdata(rc) if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) > zero) then IPD_Data(nb)%Coupling%tseain_cpl(ix) = datar8(i,j) IPD_Data(nb)%Sfcprop%tsfco(ix) = datar8(i,j) -! IPD_Data(nb)%Sfcprop%tsfc(ix) = datar8(i,j) endif enddo enddo @@ -1747,13 +1746,17 @@ subroutine assign_importdata(rc) IPD_Data(nb)%Coupling%ficein_cpl(ix) = zero IPD_Data(nb)%Coupling%slimskin_cpl(ix) = IPD_Data(nb)%Sfcprop%slmsk(ix) 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) = max(zero, min(datar8(i,j),one)) - IPD_Data(nb)%Sfcprop%slmsk(ix) = 2. !slmsk=2 crashes in gcycle on partial land points + IPD_Data(nb)%Coupling%ficein_cpl(ix) = max(zero, min(one, datar8(i,j)/IPD_Data(nb)%Sfcprop%oceanfrac(ix))) !LHS: ice frac wrt water area + if (IPD_Data(nb)%Coupling%ficein_cpl(ix) > one-epsln) IPD_Data(nb)%Coupling%ficein_cpl(ix)=one + if (IPD_Data(nb)%Coupling%ficein_cpl(ix) >= IPD_control%min_seaice) then + if (abs(one-IPD_Data(nb)%Sfcprop%oceanfrac(ix)) < epsln) IPD_Data(nb)%Sfcprop%slmsk(ix) = 2. !slmsk=2 crashes in gcycle on partial land points IPD_Data(nb)%Coupling%slimskin_cpl(ix) = 4. - elseif (abs(one-IPD_Data(nb)%Sfcprop%oceanfrac(ix)) < epsln) then - IPD_Data(nb)%Sfcprop%slmsk(ix) = zero - IPD_Data(nb)%Coupling%slimskin_cpl(ix) = zero + else + IPD_Data(nb)%Coupling%ficein_cpl(ix) = zero + if (abs(one-IPD_Data(nb)%Sfcprop%oceanfrac(ix)) < epsln) then + IPD_Data(nb)%Sfcprop%slmsk(ix) = zero + IPD_Data(nb)%Coupling%slimskin_cpl(ix) = zero + end if endif endif enddo @@ -1924,7 +1927,7 @@ subroutine assign_importdata(rc) ix = Atm_block%ixp(i,j) if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) > zero) then !if it is ocean or ice get surface temperature from mediator - if(IPD_Data(nb)%Coupling%ficein_cpl(ix) >= IPD_control%min_seaice*IPD_Data(nb)%Sfcprop%oceanfrac(ix)) then + if(IPD_Data(nb)%Coupling%ficein_cpl(ix) >= IPD_control%min_seaice) then IPD_Data(nb)%Sfcprop%tisfc(ix) = IPD_Data(nb)%Coupling%tisfcin_cpl(ix) IPD_Data(nb)%Sfcprop%fice(ix) = IPD_Data(nb)%Coupling%ficein_cpl(ix) IPD_Data(nb)%Sfcprop%hice(ix) = IPD_Data(nb)%Coupling%hicein_cpl(ix) diff --git a/ccpp/build_ccpp.sh b/ccpp/build_ccpp.sh index f6bec6ff92..b5b8c53d73 100755 --- a/ccpp/build_ccpp.sh +++ b/ccpp/build_ccpp.sh @@ -5,7 +5,7 @@ set -eu # List of valid/tested machines VALID_MACHINES=( wcoss_cray wcoss_dell_p3 gaea.intel jet.intel \ - hera.intel \ + hera.intel hera.gnu \ cheyenne.intel cheyenne.intel-impi cheyenne.gnu cheyenne.pgi endeavor.intel \ stampede.intel supermuc_phase2.intel macosx.gnu \ linux.intel linux.gnu linux.pgi ) @@ -128,6 +128,12 @@ fi if [[ "${MAKE_OPT}" == *"STATIC=Y"* ]]; then CCPP_CMAKE_FLAGS="${CCPP_CMAKE_FLAGS} -DSTATIC=ON" else + # hera.gnu uses the NCEPLIBS-external/NCEPLIBS umbrella build libraries, + # which cannot be linked dynamically at this point (missing -fPIC flag) + if [[ "${MACHINE_ID}" == "hera.gnu" ]]; then + echo "Dynamic CCPP build not supported on hera.gnu at this time." + exit 1 + fi # Dynamic builds require linking the NCEPlibs, provide path to them CCPP_CMAKE_FLAGS="${CCPP_CMAKE_FLAGS} -DSTATIC=OFF -DBACIO_LIB4=${BACIO_LIB4} -DSP_LIBd=${SP_LIBd} -DW3NCO_LIBd=${W3NCO_LIBd}" fi diff --git a/ccpp/config/ccpp_prebuild_config.py b/ccpp/config/ccpp_prebuild_config.py index 7dd9997509..6258c4b7ad 100755 --- a/ccpp/config/ccpp_prebuild_config.py +++ b/ccpp/config/ccpp_prebuild_config.py @@ -200,6 +200,7 @@ 'FV3/ccpp/physics/physics/cu_gf_driver.F90' : [ 'slow_physics' ], 'FV3/ccpp/physics/physics/cu_gf_driver_post.F90' : [ 'slow_physics' ], 'FV3/ccpp/physics/physics/moninedmf.f' : [ 'slow_physics' ], + 'FV3/ccpp/physics/physics/moninedmf_hafs.f' : [ 'slow_physics' ], 'FV3/ccpp/physics/physics/moninshoc.f' : [ 'slow_physics' ], 'FV3/ccpp/physics/physics/satmedmfvdif.F' : [ 'slow_physics' ], 'FV3/ccpp/physics/physics/satmedmfvdifq.F' : [ 'slow_physics' ], diff --git a/ccpp/framework b/ccpp/framework index e772109863..d32b965b11 160000 --- a/ccpp/framework +++ b/ccpp/framework @@ -1 +1 @@ -Subproject commit e7721098639ee73c2a69ee0e8423e8905549e240 +Subproject commit d32b965b11882a42d9db522dc13823b7720b63aa diff --git a/ccpp/physics b/ccpp/physics index e7909b4f34..3d45390dcd 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit e7909b4f34336742cab7d4ee687ae467f8fd646c +Subproject commit 3d45390dcd118c3cb0c6fb008e2d21608bbf0648 diff --git a/ccpp/set_compilers.sh b/ccpp/set_compilers.sh index 0f591d6c08..5d797a161d 100755 --- a/ccpp/set_compilers.sh +++ b/ccpp/set_compilers.sh @@ -39,6 +39,13 @@ case "$MACHINE_ID" in export F77=mpiifort export F90=mpiifort ;; + hera.gnu) + export CC=mpicc + export CXX=mpicxx + export FC=mpif90 + export F77=mpif77 + export F90=mpif90 + ;; cheyenne.intel) export CC=mpicc export CXX=mpicxx diff --git a/gfsphysics/GFS_layer/GFS_physics_driver.F90 b/gfsphysics/GFS_layer/GFS_physics_driver.F90 index cf8a1527c1..c2bf254acc 100644 --- a/gfsphysics/GFS_layer/GFS_physics_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_physics_driver.F90 @@ -1069,7 +1069,7 @@ subroutine GFS_physics_driver & ! --- ... xw: transfer ice thickness & concentration from global to local variables !## CCPP ## global to local variable transfer not necessary for these two zice(i) = Sfcprop%hice(i) - fice(i) = Sfcprop%fice(i) ! ice fraction of lake/ocean wrt whole cell + fice(i) = Sfcprop%fice(i) !*## CCPP ##* !## CCPP ##* GFS_surface_composites.F90/GFS_surface_composites_pre_run tice(i) = Sfcprop%tisfc(i) @@ -1126,35 +1126,31 @@ subroutine GFS_physics_driver & !*## CCPP ## !## CCPP ##* GFS_surface_composites.F90/GFS_surface_composites_pre - if (Model%frac_grid) then ! here Sfcprop%fice is fraction of the whole grid that is ice + if (Model%frac_grid) then do i = 1, IM frland(i) = Sfcprop%landfrac(i) if (frland(i) > zero) dry(i) = .true. - tem = one - frland(i) - if (tem > epsln) then + if (frland(i) < one) then if (flag_cice(i)) then - if (fice(i) >= Model%min_seaice*tem) then + if (fice(i) >= Model%min_seaice) then icy(i) = .true. else fice(i) = zero endif else - if (fice(i) >= Model%min_lakeice*tem) then + if (fice(i) >= Model%min_lakeice) then icy(i) = .true. - fice(i) = fice(i)/tem ! fice is fraction of ocean/lake else fice(i) = zero endif endif + if (fice(i) < one) then + wet(i)=.true. !some open ocean/lake water exists + if (.not. Model%cplflx) Sfcprop%tsfco(i) = max(Sfcprop%tsfco(i), Sfcprop%tisfc(i), tgice) + + end if else fice(i) = zero - endif - ! ocean/lake area that is not frozen - if (tem-fice(i) > epsln) then - wet(i) = .true. ! there is some open water! - if (.not. Model%cplflx) Sfcprop%tsfco(i) = max(Sfcprop%tsfco(i), Sfcprop%tisfc(i), tgice) -! if (icy(i)) Sfcprop%tsfco(i) = max(Sfcprop%tsfco(i), tgice) -! if (icy(i)) Sfcprop%tsfco(i) = max(Sfcprop%tisfc(i), tgice) endif enddo else @@ -1524,15 +1520,15 @@ subroutine GFS_physics_driver & if (Model%frac_grid) then do i=1,im - tem = one - Sfcprop%fice(i) - frland(i) + tem = (one - frland(i)) * fice(i) ! tem = ice fraction wrt whole cell if (flag_cice(i)) then adjsfculw(i) = adjsfculw3(i,1) * frland(i) & - + Coupling%ulwsfcin_cpl(i) * Sfcprop%fice(i) & - + adjsfculw3(i,3) * tem + + Coupling%ulwsfcin_cpl(i) * tem & + + adjsfculw3(i,3) * (one - frland(i) - tem) else adjsfculw(i) = adjsfculw3(i,1) * frland(i) & - + adjsfculw3(i,2) * Sfcprop%fice(i) & - + adjsfculw3(i,3) * tem + + adjsfculw3(i,2) * tem & + + adjsfculw3(i,3) * (one - frland(i) - tem) endif enddo else @@ -1540,16 +1536,16 @@ subroutine GFS_physics_driver & if (dry(i)) then ! all land adjsfculw(i) = adjsfculw3(i,1) elseif (icy(i)) then ! ice (and water) - tem = one - Sfcprop%fice(i) + tem = one - fice(i) if (flag_cice(i)) then if (wet(i) .and. adjsfculw3(i,3) /= huge) then - adjsfculw(i) = Coupling%ulwsfcin_cpl(i)*Sfcprop%fice(i) + adjsfculw3(i,3)*tem + adjsfculw(i) = Coupling%ulwsfcin_cpl(i)*fice(i) + adjsfculw3(i,3)*tem else adjsfculw(i) = Coupling%ulwsfcin_cpl(i) endif else if (wet(i) .and. adjsfculw3(i,3) /= huge) then - adjsfculw(i) = adjsfculw3(i,2)*Sfcprop%fice(i) + adjsfculw3(i,3)*tem + adjsfculw(i) = adjsfculw3(i,2)*fice(i) + adjsfculw3(i,3)*tem else adjsfculw(i) = adjsfculw3(i,2) endif @@ -2002,9 +1998,9 @@ subroutine GFS_physics_driver & do i=1, im ! ! Three-way composites (fields from sfc_diff) - txl = Sfcprop%landfrac(i) - txi = Sfcprop%fice(i) ! here Sfcprop%fice is grid fraction that is ice - txo = one - txl - txi + txl = frland(i) + txi = fice(i)*(one - frland(i)) ! txi = ice fraction wrt whole cell + txo = max(zero, one - txl - txi) Sfcprop%zorl(i) = txl*zorl3(i,1) + txi*zorl3(i,2) + txo*zorl3(i,3) cd(i) = txl*cd3(i,1) + txi*cd3(i,2) + txo*cd3(i,3) cdq(i) = txl*cdq3(i,1) + txi*cdq3(i,2) + txo*cdq3(i,3) @@ -2062,8 +2058,7 @@ subroutine GFS_physics_driver & if (.not. flag_cice(i)) then if (islmsk(i) == 2) then ! return updated lake ice thickness & concentration to global array Sfcprop%hice(i) = zice(i) -! Sfcprop%fice(i) = fice(i) * Sfcprop%lakefrac(i) ! fice is fraction of lake area that is frozen - Sfcprop%fice(i) = fice(i) * (one-Sfcprop%landfrac(i)) ! fice is fraction of wet area that is frozen + Sfcprop%fice(i) = fice(i) Sfcprop%tisfc(i) = tice(i) else ! this would be over open ocean or land (no ice fraction) Sfcprop%hice(i) = zero @@ -2121,7 +2116,7 @@ subroutine GFS_physics_driver & Sfcprop%zorlo(i) = zorl3(i,3) if (flag_cice(i) .and. wet(i)) then ! this was already done for lake ice in sfc_sice - txi = Sfcprop%fice(i) + txi = fice(i) txo = one - txi evap(i) = txi * evap3(i,2) + txo * evap3(i,3) hflx(i) = txi * hflx3(i,2) + txo * hflx3(i,3) @@ -2843,34 +2838,31 @@ subroutine GFS_physics_driver & if (Model%cplflx) then do i=1,im if (Sfcprop%oceanfrac(i) > zero) then ! Ocean only, NO LAKES - if (fice(i) == Sfcprop%oceanfrac(i)) then ! use results from CICE + if (Sfcprop%fice(i) > one - epsln) then ! no open water, thus use results from CICE Coupling%dusfci_cpl(i) = dusfc_cice(i) Coupling%dvsfci_cpl(i) = dvsfc_cice(i) Coupling%dtsfci_cpl(i) = dtsfc_cice(i) Coupling%dqsfci_cpl(i) = dqsfc_cice(i) - - elseif (wet(i)) then ! use stress_ocean from sfc_diff for opw component at mixed point - if (icy(i) .or. dry(i)) then - tem1 = max(Diag%q1(i), 1.e-8) - rho = Statein%prsl(i,1) / (con_rd*Diag%t1(i)*(one+con_fvirt*tem1)) - if (wind(i) > zero) then - tem = - rho * stress3(i,3) / wind(i) - Coupling%dusfci_cpl(i) = tem * Statein%ugrs(i,1) ! U-momentum flux - Coupling%dvsfci_cpl(i) = tem * Statein%vgrs(i,1) ! V-momentum flux - else - Coupling%dusfci_cpl(i) = zero - Coupling%dvsfci_cpl(i) = zero - endif - Coupling%dtsfci_cpl(i) = con_cp * rho * hflx3(i,3) ! sensible heat flux over open ocean - Coupling%dqsfci_cpl(i) = con_hvap * rho * evap3(i,3) ! latent heat flux over open ocean + elseif (icy(i) .or. dry(i)) then ! use stress_ocean from sfc_diff for opw component at mixed point + tem1 = max(Diag%q1(i), 1.e-8) + rho = Statein%prsl(i,1) / (con_rd*Diag%t1(i)*(one+con_fvirt*tem1)) + if (wind(i) > zero) then + tem = - rho * stress3(i,3) / wind(i) + Coupling%dusfci_cpl(i) = tem * Statein%ugrs(i,1) ! U-momentum flux + Coupling%dvsfci_cpl(i) = tem * Statein%vgrs(i,1) ! V-momentum flux + else + Coupling%dusfci_cpl(i) = zero + Coupling%dvsfci_cpl(i) = zero + endif + Coupling%dtsfci_cpl(i) = con_cp * rho * hflx3(i,3) ! sensible heat flux over open ocean + Coupling%dqsfci_cpl(i) = con_hvap * rho * evap3(i,3) ! latent heat flux over open ocean ! if (lprnt .and. i == ipr) write(0,*)' hflx33=',hflx3(i,3),' evap33=',evap3(i,3), & ! ' con_cp=',con_cp,' rho=',rho,' con_hvap=',con_hvap - else ! use results from PBL scheme for 100% open ocean - Coupling%dusfci_cpl(i) = dusfc1(i) - Coupling%dvsfci_cpl(i) = dvsfc1(i) - Coupling%dtsfci_cpl(i) = dtsfc1(i) - Coupling%dqsfci_cpl(i) = dqsfc1(i) - endif + else ! use results from PBL scheme for 100% open ocean + Coupling%dusfci_cpl(i) = dusfc1(i) + Coupling%dvsfci_cpl(i) = dvsfc1(i) + Coupling%dtsfci_cpl(i) = dtsfc1(i) + Coupling%dqsfci_cpl(i) = dqsfc1(i) endif Coupling%dusfc_cpl (i) = Coupling%dusfc_cpl(i) + Coupling%dusfci_cpl(i) * dtf diff --git a/gfsphysics/GFS_layer/GFS_typedefs.F90 b/gfsphysics/GFS_layer/GFS_typedefs.F90 index ea56d63a4c..8250bea53c 100644 --- a/gfsphysics/GFS_layer/GFS_typedefs.F90 +++ b/gfsphysics/GFS_layer/GFS_typedefs.F90 @@ -709,7 +709,8 @@ module GFS_typedefs !--- Thompson's microphysical parameters logical :: ltaerosol !< flag for aerosol version - logical :: lradar !< flag for radar reflectivity + logical :: lradar !< flag for radar reflectivity + real(kind=kind_phys) :: nsradar_reset !< seconds between resetting radar reflectivity calculation real(kind=kind_phys) :: ttendlim !< temperature tendency limiter per time step in K/s !--- GFDL microphysical paramters @@ -1785,6 +1786,7 @@ module GFS_typedefs real (kind=kind_phys), pointer :: qss_ice(:) => null() !< real (kind=kind_phys), pointer :: qss_land(:) => null() !< real (kind=kind_phys), pointer :: qss_ocean(:) => null() !< + logical :: radar_reset !< real (kind=kind_phys) :: raddt !< real (kind=kind_phys), pointer :: rainmp(:) => null() !< real (kind=kind_phys), pointer :: raincd(:) => null() !< @@ -2781,11 +2783,11 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & real(kind=kind_phys) :: evpco = 2.0d-5 !< [in] coeff for evaporation of largescale rain real(kind=kind_phys) :: wminco(2) = (/1.0d-5,1.0d-5/) !< [in] water and ice minimum threshold for Zhao !---Max hourly - real(kind=kind_phys) :: avg_max_length = 3600. !< reset value in seconds for max hourly. + real(kind=kind_phys) :: avg_max_length = 3600. !< reset value in seconds for max hourly !--- Ferrier-Aligo microphysical parameters #ifdef CCPP - real(kind=kind_phys) :: rhgrd = 0.98 !< fer_hires microphysics only - logical :: spec_adv = .true. !< Individual cloud species advected + real(kind=kind_phys) :: rhgrd = 0.98 !< fer_hires microphysics only + logical :: spec_adv = .true. !< Individual cloud species advected #endif !--- M-G microphysical parameters integer :: fprcp = 0 !< no prognostic rain and snow (MG) @@ -2823,6 +2825,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & !--- Thompson microphysical parameters logical :: ltaerosol = .false. !< flag for aerosol version logical :: lradar = .false. !< flag for radar reflectivity + real(kind=kind_phys) :: nsradar_reset = -999.0 !< seconds between resetting radar reflectivity calculation, set to <0 for every time step real(kind=kind_phys) :: ttendlim = -999.0 !< temperature tendency limiter, set to <0 to deactivate !--- GFDL microphysical parameters @@ -3119,7 +3122,8 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & mg_do_graupel, mg_do_hail, mg_nccons, mg_nicons, mg_ngcons, & mg_ncnst, mg_ninst, mg_ngnst, sed_supersat, do_sb_physics, & mg_alf, mg_qcmin, mg_do_ice_gmao, mg_do_liq_liu, & - ltaerosol, lradar, lrefres, ttendlim, lgfdlmprad, & + ltaerosol, lradar, nsradar_reset, lrefres, ttendlim, & + lgfdlmprad, & !--- max hourly avg_max_length, & !--- land/surface model control @@ -3408,6 +3412,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & !--- Thompson MP parameters Model%ltaerosol = ltaerosol Model%lradar = lradar + Model%nsradar_reset = nsradar_reset Model%ttendlim = ttendlim !--- F-A MP parameters #ifdef CCPP @@ -4201,7 +4206,10 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & if (Model%me == Model%master) print *,' Using Thompson double moment', & ' microphysics',' ltaerosol = ',Model%ltaerosol, & ' ttendlim =',Model%ttendlim, & - ' lradar =',Model%lradar,Model%num_p3d,Model%num_p2d + ' lradar =',Model%lradar, & + ' nsradar_reset =',Model%nsradar_reset, & + ' num_p3d =',Model%num_p3d, & + ' num_p2d =',Model%num_p2d else if (Model%imp_physics == Model%imp_physics_mg) then ! Morrison-Gettelman Microphysics Model%npdf3d = 0 @@ -4473,6 +4481,7 @@ subroutine control_print(Model) print *, ' Thompson microphysical parameters' print *, ' ltaerosol : ', Model%ltaerosol print *, ' lradar : ', Model%lradar + print *, ' nsradar_reset : ', Model%nsradar_reset print *, ' lrefres : ', Model%lrefres print *, ' ttendlim : ', Model%ttendlim print *, ' ' @@ -6518,6 +6527,12 @@ subroutine interstitial_phys_reset (Interstitial, Model) ! ! Set flag for resetting maximum hourly output fields Interstitial%reset = mod(Model%kdt-1, nint(Model%avg_max_length/Model%dtp)) == 0 + ! Set flag for resetting radar reflectivity calculation + if (Model%nsradar_reset<0) then + Interstitial%radar_reset = .true. + else + Interstitial%radar_reset = mod(Model%kdt-1, nint(Model%nsradar_reset/Model%dtp)) == 0 + end if ! end subroutine interstitial_phys_reset @@ -6705,6 +6720,7 @@ subroutine interstitial_print(Interstitial, Model, mpirank, omprank, blkno) write (0,*) 'sum(Interstitial%qss_ice ) = ', sum(Interstitial%qss_ice ) write (0,*) 'sum(Interstitial%qss_land ) = ', sum(Interstitial%qss_land ) write (0,*) 'sum(Interstitial%qss_ocean ) = ', sum(Interstitial%qss_ocean ) + write (0,*) 'Interstitial%radar_reset = ', Interstitial%radar_reset write (0,*) 'Interstitial%raddt = ', Interstitial%raddt write (0,*) 'sum(Interstitial%raincd ) = ', sum(Interstitial%raincd ) write (0,*) 'sum(Interstitial%raincs ) = ', sum(Interstitial%raincs ) diff --git a/gfsphysics/GFS_layer/GFS_typedefs.meta b/gfsphysics/GFS_layer/GFS_typedefs.meta index 48d26266b8..adb24a3fc5 100644 --- a/gfsphysics/GFS_layer/GFS_typedefs.meta +++ b/gfsphysics/GFS_layer/GFS_typedefs.meta @@ -7565,6 +7565,12 @@ dimensions = (horizontal_dimension) type = real kind = kind_phys +[radar_reset] + standard_name = flag_for_resetting_radar_reflectivity_calculation + long_name = flag for resetting radar reflectivity calculation + units = flag + dimensions = () + type = logical [raddt] standard_name = time_step_for_radiation long_name = radiation time step diff --git a/io/FV3GFS_io.F90 b/io/FV3GFS_io.F90 index 25735d7279..652955816b 100644 --- a/io/FV3GFS_io.F90 +++ b/io/FV3GFS_io.F90 @@ -29,6 +29,7 @@ module FV3GFS_io_mod use diag_data_mod, only: output_fields, max_output_fields use diag_util_mod, only: find_input_field use constants_mod, only: grav, rdgas + use physcons, only: con_tice !saltwater freezing temp (K) ! !--- GFS physics modules !#ifndef CCPP @@ -104,6 +105,7 @@ module FV3GFS_io_mod real, parameter :: missing_value = 9.99e20 real, parameter:: stndrd_atmos_ps = 101325. real, parameter:: stndrd_atmos_lapse = 0.0065 + real, parameter:: drythresh = 1.e-4 !--- miscellaneous other variables logical :: use_wrtgridcomp_output = .FALSE. @@ -935,6 +937,31 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) Sfcprop(nb)%tsfcl(ix) = sfc_var2(i,j,33) !--- sfcl (temp on land portion of a cell) Sfcprop(nb)%zorll(ix) = sfc_var2(i,j,34) !--- zorll (zorl on land portion of a cell) end if + + if(Model%frac_grid) then ! obtain slmsk from landfrac +!! next 5 lines are temporary till lake model is available + if (Sfcprop(nb)%lakefrac(ix) > 0.0) then + Sfcprop(nb)%lakefrac(ix) = nint(Sfcprop(nb)%lakefrac(ix)) + Sfcprop(nb)%landfrac(ix) = 1.-Sfcprop(nb)%lakefrac(ix) + if (Sfcprop(nb)%lakefrac(ix) == 0) Sfcprop(nb)%fice(ix)=0. + end if + Sfcprop(nb)%slmsk(ix) = ceiling(Sfcprop(nb)%landfrac(ix)) + if (Sfcprop(nb)%fice(ix) > 0. .and. Sfcprop(nb)%landfrac(ix)==0.) Sfcprop(nb)%slmsk(ix) = 2 ! land dominates ice if co-exist + else ! obtain landfrac from slmsk + if (Sfcprop(nb)%slmsk(ix) > 1.9) then + Sfcprop(nb)%landfrac(ix) = 0.0 + else + Sfcprop(nb)%landfrac(ix) = Sfcprop(nb)%slmsk(ix) + endif + end if + + if (Sfcprop(nb)%lakefrac(ix) > 0.0) then + Sfcprop(nb)%oceanfrac(ix) = 0.0 ! lake & ocean don't coexist in a cell + if (Sfcprop(nb)%fice(ix) < Model%min_lakeice) Sfcprop(nb)%fice(ix) = 0. + else + Sfcprop(nb)%oceanfrac(ix) = 1.0 - Sfcprop(nb)%landfrac(ix) + if (Sfcprop(nb)%fice(ix) < Model%min_seaice) Sfcprop(nb)%fice(ix) = 0. + endif ! !--- NSSTM variables if ((Model%nstf_name(1) > 0) .and. (Model%nstf_name(2) == 1)) then @@ -1103,7 +1130,7 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) do nb = 1, Atm_block%nblks do ix = 1, Atm_block%blksz(nb) Sfcprop(nb)%sncovr(ix) = 0.0 - if (Sfcprop(nb)%slmsk(ix) > 0.001) then + if (Sfcprop(nb)%landfrac(ix) >= drythresh .or. Sfcprop(nb)%fice(ix) >= Model%min_seaice) then vegtyp = Sfcprop(nb)%vtype(ix) if (vegtyp == 0) vegtyp = 7 rsnow = 0.001*Sfcprop(nb)%weasd(ix)/snupx(vegtyp) @@ -1116,21 +1143,43 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) enddo enddo endif + + if(Model%cplflx .or. Model%frac_grid) then + if (nint(sfc_var2(1,1,33)) == -9999) then + if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing tsfcl') + do nb = 1, Atm_block%nblks + do ix = 1, Atm_block%blksz(nb) + Sfcprop(nb)%tsfcl(ix) = Sfcprop(nb)%tsfco(ix) !--- compute tsfcl from existing variables + enddo + enddo + endif + + if (nint(sfc_var2(1,1,34)) == -9999) then + if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing zorll') + do nb = 1, Atm_block%nblks + do ix = 1, Atm_block%blksz(nb) + Sfcprop(nb)%zorll(ix) = Sfcprop(nb)%zorlo(ix) !--- compute zorll from existing variables + enddo + enddo + endif + endif + !#endif if(Model%frac_grid) then ! 3-way composite do nb = 1, Atm_block%nblks do ix = 1, Atm_block%blksz(nb) - tem = 1.0 - Sfcprop(nb)%landfrac(ix) - Sfcprop(nb)%fice(ix) + Sfcprop(nb)%tsfco(ix) = max(con_tice, Sfcprop(nb)%tsfco(ix)) + tem = (1.-Sfcprop(nb)%landfrac(ix)) * Sfcprop(nb)%fice(ix) ! tem = ice fraction wrt whole cell Sfcprop(nb)%zorl(ix) = Sfcprop(nb)%zorll(ix) * Sfcprop(nb)%landfrac(ix) & - + Sfcprop(nb)%zorll(ix) * Sfcprop(nb)%fice(ix) & !zorl ice = zorl land - + Sfcprop(nb)%zorlo(ix) * tem + + Sfcprop(nb)%zorll(ix) * tem & !zorl ice = zorl land + + Sfcprop(nb)%zorlo(ix) * (1.-Sfcprop(nb)%landfrac(ix)-tem) Sfcprop(nb)%tsfc(ix) = Sfcprop(nb)%tsfcl(ix) * Sfcprop(nb)%landfrac(ix) & - + Sfcprop(nb)%tisfc(ix) * Sfcprop(nb)%fice(ix) & - + Sfcprop(nb)%tsfco(ix) * tem + + Sfcprop(nb)%tisfc(ix) * tem & + + Sfcprop(nb)%tsfco(ix) * (1.-Sfcprop(nb)%landfrac(ix)-tem) enddo enddo - else ! in this case ice fraction is fraction of water fraction + else do nb = 1, Atm_block%nblks do ix = 1, Atm_block%blksz(nb) !--- specify tsfcl/zorll from existing variable tsfco/zorlo @@ -1138,33 +1187,10 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) Sfcprop(nb)%zorll(ix) = Sfcprop(nb)%zorlo(ix) Sfcprop(nb)%zorl(ix) = Sfcprop(nb)%zorlo(ix) Sfcprop(nb)%tsfc(ix) = Sfcprop(nb)%tsfco(ix) - if (abs(1.0-Sfcprop(nb)%slmsk(ix)) < 0.1) then - Sfcprop(nb)%landfrac(ix) = 1.0 ! land - Sfcprop(nb)%lakefrac(ix) = 0.0 - else - Sfcprop(nb)%landfrac(ix) = 0.0 ! water - if (Sfcprop(nb)%lakefrac(ix) > 0.0 .or. & - (Sfcprop(nb)%oro_uf(ix) > Model%min_lake_height .and. .not. Model%ignore_lake) ) then - Sfcprop(nb)%lakefrac(ix) = 1.0 ! lake - else - Sfcprop(nb)%lakefrac(ix) = 0.0 ! ocean - endif - endif enddo enddo endif ! if (Model%frac_grid) - do nb = 1, Atm_block%nblks - do ix = 1, Atm_block%blksz(nb) - if (Sfcprop(nb)%lakefrac(ix) > 0.0) then - Sfcprop(nb)%oceanfrac(ix) = 0.0 ! lake & ocean don't coexist in a cell - else - Sfcprop(nb)%oceanfrac(ix) = 1.0 - Sfcprop(nb)%landfrac(ix) !LHS:ocean frac [0:1] - endif - - enddo - enddo - if (Model%lsm == Model%lsm_noahmp) then if (nint(sfc_var2(1,1,nvar_s2m+19)) == -66666) then if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver:: - Cold start Noah MP ') @@ -1209,7 +1235,7 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) Sfcprop(nb)%smoiseq(ix, 1:4) = missing_value Sfcprop(nb)%zsnsoxy(ix, -2:4) = missing_value - if (Sfcprop(nb)%slmsk(ix) > 0.01) then + if (Sfcprop(nb)%landfrac(ix) >= drythresh) then Sfcprop(nb)%tvxy(ix) = Sfcprop(nb)%tsfcl(ix) Sfcprop(nb)%tgxy(ix) = Sfcprop(nb)%tsfcl(ix)