From bee8c37d57110146689aaa9af93e659af306c839 Mon Sep 17 00:00:00 2001 From: Kun Gao Date: Thu, 21 Mar 2024 22:08:44 -0400 Subject: [PATCH 1/4] update SHiELD_physics to use sfc fluxes from coupler --- gsmphys/sfc_diff_gfdl_coupled.f | 365 +++++++++++++++++++++++++++++++ gsmphys/sfc_ocean_gfdl_coupled.f | 152 +++++++++++++ 2 files changed, 517 insertions(+) create mode 100644 gsmphys/sfc_diff_gfdl_coupled.f create mode 100644 gsmphys/sfc_ocean_gfdl_coupled.f diff --git a/gsmphys/sfc_diff_gfdl_coupled.f b/gsmphys/sfc_diff_gfdl_coupled.f new file mode 100644 index 00000000..fab6f752 --- /dev/null +++ b/gsmphys/sfc_diff_gfdl_coupled.f @@ -0,0 +1,365 @@ + subroutine sfc_diff_gfdl_coupled(im,ps,u1,v1,t1,q1,z1, + & snwdph,tskin,z0rl,ztrl,cm,ch,rb, + & prsl1,prslki,islimsk, + & stress,fm,fh, + & ustar,wind,ddvel,fm10,fh2, + & sigmaf,vegtype,shdmax,ivegsrc, + & tsurf,flag_iter) !,redrag, +! & z0s_max) +! & do_z0_moon, do_z0_hwrf15, do_z0_hwrf17, +! & do_z0_hwrf17_hwonly, wind_th_hwrf) + + use machine , only : kind_phys + use funcphys, only : fpvs + use physcons, grav => con_g, cp => con_cp + &, rvrdm1 => con_fvirt, rd => con_rd + &, eps => con_eps, epsm1 => con_epsm1 + + implicit none + +! --- input/output + + integer im, ivegsrc + + real(kind=kind_phys), dimension(im)::ps, u1, v1, t1, q1, z1 + &, tskin, z0rl, ztrl, cm, ch, rb + &, prsl1, prslki, stress + &, fm, fh, ustar, wind, ddvel + &, fm10, fh2, sigmaf, shdmax + &, tsurf, snwdph + integer, dimension(im) ::vegtype, islimsk + + logical flag_iter(im) +! logical redrag +! logical do_z0_moon, do_z0_hwrf15, do_z0_hwrf17 ! kgao +! &, do_z0_hwrf17_hwonly ! kgao + +! --- local + + integer i +! + real(kind=kind_phys) aa, aa0, bb, bb0, dtv, adtv,qs1, + & hl1, hl12, pm, ph, pm10, ph2, rat, + & thv1, tvs, z1i, z0, zt, z0max, ztmax, + & fms, fhs, hl0, hl0inf, hlinf, + & hl110, hlt, hltinf, olinf, + & restar, czilc, tem1, tem2, + & u10m, v10m, ws10m, ws10m_moon, !kgao + & z0_1, zt_1, fm1, fh1, ustar_1, ztmax_1 !kgao +! + +! real(kind=kind_phys),intent(in ) :: z0s_max, wind_th_hwrf ! kgao + + real(kind=kind_phys), parameter :: + & charnock=.014, ca=.4 + &, vis=1.4e-5, rnu=1.51e-5, visi=1.0/vis + &, log01=log(0.01), log05=log(0.05), log07=log(0.07) + &, ztmin1=-999.0 + +!================================================ +! Main program starts here +!================================================ + + do i=1,im + + if(flag_iter(i)) then + +! --- get variables at model lowest layer and surface (water/ice/land) + + wind(i) = max(sqrt(u1(i)*u1(i) + v1(i)*v1(i)) + & + max(0.0, min(ddvel(i), 30.0)), 1.0) + tem1 = 1.0 + rvrdm1 * max(q1(i),1.e-8) + thv1 = t1(i) * prslki(i) * tem1 + tvs = 0.5 * (tsurf(i)+tskin(i)) * tem1 + qs1 = fpvs(t1(i)) + qs1 = max(1.0e-8, eps * qs1 / (prsl1(i) + epsm1 * qs1)) + + !(sea/land/ice mask =0/1/2) + if(islimsk(i) == 1 .or. islimsk(i) == 2) then ! over land or sea ice + +!================================================ +! if over land or sea ice: +! step 1 - get z0/zt +! step 2 - call similarity +!================================================ + +! --- get surface roughness for momentum (z0) + + z0 = 0.01 * z0rl(i) + z0max = max(1.0e-6, min(z0,z1(i))) + + !xubin's new z0 over land and sea ice + tem1 = 1.0 - shdmax(i) ! shdmax is max vegetation area fraction + tem2 = tem1 * tem1 + tem1 = 1.0 - tem2 + + if( ivegsrc == 1 ) then + + if (vegtype(i) == 10) then + z0max = exp( tem2*log01 + tem1*log07 ) + elseif (vegtype(i) == 6) then + z0max = exp( tem2*log01 + tem1*log05 ) + elseif (vegtype(i) == 7) then +! z0max = exp( tem2*log01 + tem1*log01 ) + z0max = 0.01 + elseif (vegtype(i) == 16) then +! z0max = exp( tem2*log01 + tem1*log01 ) + z0max = 0.01 + else + z0max = exp( tem2*log01 + tem1*log(z0max) ) + endif + + elseif (ivegsrc == 2 ) then + + if (vegtype(i) == 7) then + z0max = exp( tem2*log01 + tem1*log07 ) + elseif (vegtype(i) == 8) then + z0max = exp( tem2*log01 + tem1*log05 ) + elseif (vegtype(i) == 9) then +! z0max = exp( tem2*log01 + tem1*log01 ) + z0max = 0.01 + elseif (vegtype(i) == 11) then +! z0max = exp( tem2*log01 + tem1*log01 ) + z0max = 0.01 + else + z0max = exp( tem2*log01 + tem1*log(z0max) ) + endif + + z0max = max(z0max,1.0e-6) + + endif + +! --- get surface roughness for heat (zt) + +! czilc = 10.0 ** (- (0.40/0.07) * z0) ! let czilc depend on canopy height + czilc = 0.8 + + tem1 = 1.0 - sigmaf(i) + ztmax = z0max*exp( - tem1*tem1 + & * czilc*ca*sqrt(ustar(i)*(0.01/1.5e-05))) + + ztmax = max(ztmax,1.0e-6) + +! --- call similarity + + call monin_obukhov_similarity + & (z1(i), snwdph(i), thv1, wind(i), z0max, ztmax, tvs, + & rb(i), fm(i), fh(i), fm10(i), fh2(i), + & cm(i), ch(i), stress(i), ustar(i)) + + elseif (islimsk(i) == 0) then + +!================================================ +! if over water +! redesigned by Kun Gao for coupling with MOM6 +! now do not update z0rl and ustar over ocean points +!================================================ + + ! --- get z0/zt + z0 = 0.01 * z0rl(i) ! kgao: from coupler + zt = 0.01 * ztrl(i) + + z0max = max(1.0e-6, min(z0,z1(i))) + ztmax = max(zt,1.0e-6) + + ! --- call similarity + call monin_obukhov_similarity + & (z1(i), snwdph(i), thv1, wind(i), z0max, ztmax, tvs, + & rb(i), fm(i), fh(i), fm10(i), fh2(i), + & cm(i), ch(i), tem1, tem2) !stress(i), ustar(i)) + + ! kgao: use ustar from coupler to update stress + stress(i) = ustar(i) * ustar(i) + + ! kgao: get zt as in old sfc_diff.f + ! this should be removed once using ztrl from coupler + ustar_1 = sqrt(grav * z0max / charnock) + restar = max(ustar_1*z0max*visi, 0.000001) + rat = min(7.0, 2.67 * sqrt(sqrt(restar)) - 2.57) + zt = z0max * exp(-rat) ! zeng, zhao and dickinson 1997 (eq 25) + + ! kgao: do not update z0 + !z0rl(i) = 100.0 * z0max + ztrl(i) = 100.0 * ztmax + + endif ! end of if(islimsk) loop + endif ! end of if(flagiter) loop + enddo ! end of do i=1,im loop + + return + end subroutine + +! ======================================================================= + + subroutine monin_obukhov_similarity + & ( z1, snwdph, thv1, wind, z0max, ztmax, tvs, + & rb, fm, fh, fm10, fh2, cm, ch, stress, ustar) + +! --- input +! z1 - lowest model level height +! snwdph - surface snow thickness +! wind - wind speed at lowest model layer +! thv1 - virtual potential temp at lowest model layer +! tvs - surface temp +! z0max - surface roughness length for momentum +! ztmax - surface roughness length for heat +! +! --- output +! rb - a bulk richardson number +! fm, fh - similarity function defined at lowest model layer +! fm10, fh2 - similarity function defined at 10m (for momentum) and 2m (for heat) +! cm, ch - surface exchange coefficients for momentum and heat +! stress - surface wind stress +! ustar - surface frictional velocity + + use machine , only : kind_phys + use physcons, grav => con_g + +! --- inputs: + real(kind=kind_phys), intent(in) :: + & z1, snwdph, thv1, wind, z0max, ztmax, tvs + +! --- outputs: + real(kind=kind_phys), intent(out) :: + & rb, fm, fh, fm10, fh2, cm, ch, stress, ustar + +! --- locals: + + real(kind=kind_phys), parameter :: alpha=5., a0=-3.975 + &, a1=12.32, alpha4=4.0*alpha + &, b1=-7.755, b2=6.041, alpha2=alpha+alpha, beta=1.0 + &, a0p=-7.941, a1p=24.75, b1p=-8.705, b2p=7.899 + &, ztmin1=-999.0, ca=.4 + + real(kind=kind_phys) aa, aa0, bb, bb0, dtv, adtv, + & hl1, hl12, pm, ph, pm10, ph2, + & z1i, + & fms, fhs, hl0, hl0inf, hlinf, + & hl110, hlt, hltinf, olinf, + & tem1, tem2, ztmax1 + + z1i = 1.0 / z1 + + tem1 = z0max/z1 + if (abs(1.0-tem1) > 1.0e-6) then + ztmax1 = - beta*log(tem1)/(alpha2*(1.-tem1)) + else + ztmax1 = 99.0 + endif + if( z0max < 0.05 .and. snwdph < 10.0 ) ztmax1 = 99.0 + +! +! compute stability indices (rb and hlinf) +! + dtv = thv1 - tvs + adtv = max(abs(dtv),0.001) + dtv = sign(1.,dtv) * adtv + rb = max(-5000.0, (grav+grav) * dtv * z1 + & / ((thv1 + tvs) * wind * wind)) + tem1 = 1.0 / z0max + tem2 = 1.0 / ztmax + fm = log((z0max+z1) * tem1) + fh = log((ztmax+z1) * tem2) + fm10 = log((z0max+10.) * tem1) + fh2 = log((ztmax+2.) * tem2) + hlinf = rb * fm * fm / fh + hlinf = min(max(hlinf,ztmin1),ztmax1) +! +! stable case +! + if (dtv >= 0.0) then + hl1 = hlinf + if(hlinf > .25) then + tem1 = hlinf * z1i + hl0inf = z0max * tem1 + hltinf = ztmax * tem1 + aa = sqrt(1. + alpha4 * hlinf) + aa0 = sqrt(1. + alpha4 * hl0inf) + bb = aa + bb0 = sqrt(1. + alpha4 * hltinf) + pm = aa0 - aa + log( (aa + 1.)/(aa0 + 1.) ) + ph = bb0 - bb + log( (bb + 1.)/(bb0 + 1.) ) + fms = fm - pm + fhs = fh - ph + hl1 = fms * fms * rb / fhs + hl1 = min(max(hl1, ztmin1), ztmax1) + endif +! +! second iteration +! + tem1 = hl1 * z1i + hl0 = z0max * tem1 + hlt = ztmax * tem1 + aa = sqrt(1. + alpha4 * hl1) + aa0 = sqrt(1. + alpha4 * hl0) + bb = aa + bb0 = sqrt(1. + alpha4 * hlt) + pm = aa0 - aa + log( (1.0+aa)/(1.0+aa0) ) + ph = bb0 - bb + log( (1.0+bb)/(1.0+bb0) ) + hl110 = hl1 * 10. * z1i + hl110 = min(max(hl110, ztmin1), ztmax1) + aa = sqrt(1. + alpha4 * hl110) + pm10 = aa0 - aa + log( (1.0+aa)/(1.0+aa0) ) + hl12 = (hl1+hl1) * z1i + hl12 = min(max(hl12,ztmin1),ztmax1) +! aa = sqrt(1. + alpha4 * hl12) + bb = sqrt(1. + alpha4 * hl12) + ph2 = bb0 - bb + log( (1.0+bb)/(1.0+bb0) ) +! +! unstable case - check for unphysical obukhov length +! + else ! dtv < 0 case + olinf = z1 / hlinf + tem1 = 50.0 * z0max + if(abs(olinf) <= tem1) then + hlinf = -z1 / tem1 + hlinf = min(max(hlinf,ztmin1),ztmax1) + endif +! +! get pm and ph +! + if (hlinf >= -0.5) then + hl1 = hlinf + pm = (a0 + a1*hl1) * hl1 / (1.+ (b1+b2*hl1) *hl1) + ph = (a0p + a1p*hl1) * hl1 / (1.+ (b1p+b2p*hl1)*hl1) + hl110 = hl1 * 10. * z1i + hl110 = min(max(hl110, ztmin1), ztmax1) + pm10 = (a0 + a1*hl110) * hl110 / (1.+(b1+b2*hl110)*hl110) + hl12 = (hl1+hl1) * z1i + hl12 = min(max(hl12, ztmin1), ztmax1) + ph2 = (a0p + a1p*hl12) * hl12 / (1.+(b1p+b2p*hl12)*hl12) + else ! hlinf < 0.05 + hl1 = -hlinf + tem1 = 1.0 / sqrt(hl1) + pm = log(hl1) + 2. * sqrt(tem1) - .8776 + ph = log(hl1) + .5 * tem1 + 1.386 +! pm = log(hl1) + 2.0 * hl1 ** (-.25) - .8776 +! ph = log(hl1) + 0.5 * hl1 ** (-.5) + 1.386 + hl110 = hl1 * 10. * z1i + hl110 = min(max(hl110, ztmin1), ztmax1) + pm10 = log(hl110) + 2.0 / sqrt(sqrt(hl110)) - .8776 +! pm10 = log(hl110) + 2. * hl110 ** (-.25) - .8776 + hl12 = (hl1+hl1) * z1i + hl12 = min(max(hl12, ztmin1), ztmax1) + ph2 = log(hl12) + 0.5 / sqrt(hl12) + 1.386 +! ph2 = log(hl12) + .5 * hl12 ** (-.5) + 1.386 + endif + + endif ! end of if (dtv >= 0 ) then loop +! +! finish the exchange coefficient computation to provide fm and fh +! + fm = fm - pm + fh = fh - ph + fm10 = fm10 - pm10 + fh2 = fh2 - ph2 + cm = ca * ca / (fm * fm) + ch = ca * ca / (fm * fh) + tem1 = 0.00001/z1 + cm = max(cm, tem1) + ch = max(ch, tem1) + stress = cm * wind * wind + ustar = sqrt(stress) + + return + end subroutine monin_obukhov_similarity diff --git a/gsmphys/sfc_ocean_gfdl_coupled.f b/gsmphys/sfc_ocean_gfdl_coupled.f new file mode 100644 index 00000000..c46883a3 --- /dev/null +++ b/gsmphys/sfc_ocean_gfdl_coupled.f @@ -0,0 +1,152 @@ +!----------------------------------- + subroutine sfc_ocean_gfdl_coupled & +!................................... +! --- inputs: + & ( im, ps, u1, v1, t1, q1, tskin, cm, ch, & + & prsl1, prslki, islimsk, ddvel, flag_iter, & + & shflx, lhflx, ! kgao: shflx and lhflx from coupler +! --- outputs: + & qsurf, cmm, chh, gflux, evap, hflx, ep & + & ) + +! ===================================================================== ! +! description: ! +! ! +! usage: ! +! ! +! call sfc_ocean ! +! inputs: ! +! ( im, ps, u1, v1, t1, q1, tskin, cm, ch, ! +! prsl1, prslki, islimsk, ddvel, flag_iter, ! +! outputs: ! +! qsurf, cmm, chh, gflux, evap, hflx, ep ) ! +! ! +! ! +! subprograms/functions called: fpvs ! +! ! +! ! +! program history log: ! +! 2005 -- created from the original progtm to account for ! +! ocean only ! +! oct 2006 -- h. wei added cmm and chh to the output ! +! apr 2009 -- y.-t. hou modified to match the modified gbphys.f ! +! reformatted the code and added program documentation ! +! sep 2009 -- s. moorthi removed rcl and made pa as pressure unit ! +! and furthur reformatted the code ! +! ! +! ! +! ==================== defination of variables ==================== ! +! ! +! inputs: size ! +! im - integer, horizontal dimension 1 ! +! ps - real, surface pressure im ! +! u1, v1 - real, u/v component of surface layer wind im ! +! t1 - real, surface layer mean temperature ( k ) im ! +! q1 - real, surface layer mean specific humidity im ! +! tskin - real, ground surface skin temperature ( k ) im ! +! cm - real, surface exchange coeff for momentum (m/s) im ! +! ch - real, surface exchange coeff heat & moisture(m/s) im ! +! prsl1 - real, surface layer mean pressure im ! +! prslki - real, im ! +! islimsk - integer, sea/land/ice mask (=0/1/2) im ! +! ddvel - real, wind enhancement due to convection (m/s) im ! +! flag_iter- logical, im ! +! ! +! outputs: ! +! qsurf - real, specific humidity at sfc im ! +! cmm - real, im ! +! chh - real, im ! +! gflux - real, ground heat flux (zero for ocean) im ! +! evap - real, evaporation from latent heat flux im ! +! hflx - real, sensible heat flux im ! +! ep - real, potential evaporation im ! +! ! +! ===================================================================== ! +! + use machine , only : kind_phys + use funcphys, only : fpvs + use physcons, only : cp => con_cp, rd => con_rd, eps => con_eps, & + & epsm1 => con_epsm1, hvap => con_hvap, & + & rvrdm1 => con_fvirt +! + implicit none +! +! --- constant parameters: + real (kind=kind_phys), parameter :: cpinv = 1.0/cp & + &, hvapi = 1.0/hvap & + &, elocp = hvap/cp + +! --- inputs: + integer, intent(in) :: im + + real (kind=kind_phys), dimension(im), intent(in) :: ps, u1, v1, & + & t1, q1, tskin, cm, ch, prsl1, prslki, ddvel + integer, dimension(im), intent(in):: islimsk + + logical, intent(in) :: flag_iter(im) + +! --- outputs: + real (kind=kind_phys), dimension(im), intent(out) :: qsurf, & + & cmm, chh, gflux, evap, hflx, ep, + & shflx, lhflx ! kgao + +! --- locals: + + real (kind=kind_phys) :: q0, qss, rch, rho, wind, tem + + integer :: i + + logical :: flag(im) +! +!===> ... begin here +! +! --- ... flag for open water + do i = 1, im + flag(i) = ( islimsk(i) == 0 .and. flag_iter(i) ) + +! --- ... initialize variables. all units are supposedly m.k.s. unless specified +! ps is in pascals, wind is wind speed, +! rho is density, qss is sat. hum. at surface + + if ( flag(i) ) then + + wind = max(sqrt(u1(i)*u1(i) + v1(i)*v1(i)) & + & + max( 0.0, min( ddvel(i), 30.0 ) ), 1.0) + + q0 = max( q1(i), 1.0e-8 ) + rho = prsl1(i) / (rd*t1(i)*(1.0 + rvrdm1*q0)) + + qss = fpvs( tskin(i) ) + qss = eps*qss / (ps(i) + epsm1*qss) + + evap(i) = 0.0 + hflx(i) = 0.0 + ep(i) = 0.0 + gflux(i) = 0.0 + +! --- ... rcp = rho cp ch v + + rch = rho * cp * ch(i) * wind + cmm(i) = cm(i) * wind + chh(i) = rho * ch(i) * wind + +! --- ... sensible and latent heat flux over open water + + !hflx(i) = rch * (tskin(i) - t1(i) * prslki(i)) + !evap(i) = elocp*rch * (qss - q0) + qsurf(i) = qss + tem = 1.0 / rho + !hflx(i) = hflx(i) * tem * cpinv + !evap(i) = evap(i) * tem * hvapi + + !kgao: get hflx and evap from data provided by coupler + hflx(i) = shflx(i) * tem * cpinv + evap(i) = lhflx(i) * tem + + endif + enddo +! + return +!................................... + end subroutine !sfc_ocean +!----------------------------------- From a54267cf4a397468118659143296bc319f0b8c20 Mon Sep 17 00:00:00 2001 From: Kun Gao Date: Thu, 21 Mar 2024 22:09:43 -0400 Subject: [PATCH 2/4] update SHiELD_physics to use sfc fluxes from coupler - part 2 --- FV3GFS/FV3GFS_io.F90 | 50 +++++++++++++++++++++++++++++++- GFS_layer/GFS_physics_driver.F90 | 24 ++++++++++----- GFS_layer/GFS_typedefs.F90 | 16 ++++++++-- 3 files changed, 80 insertions(+), 10 deletions(-) diff --git a/FV3GFS/FV3GFS_io.F90 b/FV3GFS/FV3GFS_io.F90 index f907173c..0ccf8523 100644 --- a/FV3GFS/FV3GFS_io.F90 +++ b/FV3GFS/FV3GFS_io.F90 @@ -4394,6 +4394,30 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Cldprop, & Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%v10m(:) enddo + ! KGao + idx = idx + 1 + Diag(idx)%axes = 2 + Diag(idx)%name = 'hflx' + Diag(idx)%desc = 'surface temp flux' + Diag(idx)%unit = 'K*m/s' + Diag(idx)%mod_name = 'gfs_phys' + allocate (Diag(idx)%data(nblks)) + do nb = 1,nblks + Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%hflx(:) + enddo + + ! KGao + idx = idx + 1 + Diag(idx)%axes = 2 + Diag(idx)%name = 'evap' + Diag(idx)%desc = 'surface moisture flux' + Diag(idx)%unit = 'm/s*kg/kg' + Diag(idx)%mod_name = 'gfs_phys' + allocate (Diag(idx)%data(nblks)) + do nb = 1,nblks + Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%evap(:) + enddo + idx = idx + 1 Diag(idx)%axes = 2 Diag(idx)%name = 'dpt2m' @@ -6438,6 +6462,30 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Cldprop, & Diag(idx)%data(nb)%var2 => Sfcprop(nb)%uustar(:) enddo + ! KGao + idx = idx + 1 + Diag(idx)%axes = 2 + Diag(idx)%name = 'shflx' + Diag(idx)%desc = 'shflx from coupler' + Diag(idx)%unit = 'XXX' + Diag(idx)%mod_name = 'gfs_sfc' + allocate (Diag(idx)%data(nblks)) + do nb = 1,nblks + Diag(idx)%data(nb)%var2 => Sfcprop(nb)%shflx(:) + enddo + + ! KGao + idx = idx + 1 + Diag(idx)%axes = 2 + Diag(idx)%name = 'lhflx' + Diag(idx)%desc = 'lhflx from coupler' + Diag(idx)%unit = 'XXX' + Diag(idx)%mod_name = 'gfs_sfc' + allocate (Diag(idx)%data(nblks)) + do nb = 1,nblks + Diag(idx)%data(nb)%var2 => Sfcprop(nb)%lhflx(:) + enddo + idx = idx + 1 Diag(idx)%axes = 2 Diag(idx)%name = 'slope' @@ -6668,7 +6716,7 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Cldprop, & Diag(idx)%desc = 'surface roughness [m]' Diag(idx)%unit = 'm' Diag(idx)%mod_name = 'gfs_sfc' - Diag(idx)%cnvfac = cn_one/cn_100 + Diag(idx)%cnvfac = cn_one/cn_100 allocate (Diag(idx)%data(nblks)) do nb = 1,nblks Diag(idx)%data(nb)%var2 => Sfcprop(nb)%zorl(:) diff --git a/GFS_layer/GFS_physics_driver.F90 b/GFS_layer/GFS_physics_driver.F90 index 4f4013ba..01e06a95 100644 --- a/GFS_layer/GFS_physics_driver.F90 +++ b/GFS_layer/GFS_physics_driver.F90 @@ -1136,18 +1136,20 @@ subroutine GFS_physics_driver & !!$ endif if (Model%sfc_gfdl) then -! a new and more flexible version of sfc_diff by kgao - call sfc_diff_gfdl(im,Statein%pgr, Statein%ugrs, Statein%vgrs,& +! a version of sfc_diff from coupling with MOM6 by kgao +! Sfcprop%uustar,Sfcprop%zorl not updated over ocean points +! needs to consider ztrl too + call sfc_diff_gfdl_coupled(im,Statein%pgr, Statein%ugrs, Statein%vgrs,& Statein%tgrs, Statein%qgrs, Diag%zlvl, Sfcprop%snowd, & Sfcprop%tsfc, Sfcprop%zorl, Sfcprop%ztrl, cd, & cdq, rb, Statein%prsl(1,1), work3, islmsk, stress, & Sfcprop%ffmm, Sfcprop%ffhh, Sfcprop%uustar, & wind, Tbd%phy_f2d(1,Model%num_p2d), fm10, fh2, & sigmaf, vegtype, Sfcprop%shdmax, Model%ivegsrc, & - tsurf, flag_iter, Model%redrag, Model%z0s_max, & - Model%do_z0_moon, Model%do_z0_hwrf15, & - Model%do_z0_hwrf17, Model%do_z0_hwrf17_hwonly, & - Model%wind_th_hwrf) + tsurf, flag_iter) !, Model%redrag, Model%z0s_max, & + !Model%do_z0_moon, Model%do_z0_hwrf15, & + !Model%do_z0_hwrf17, Model%do_z0_hwrf17_hwonly, & + !Model%wind_th_hwrf) else ! GFS original sfc_diff modified by kgao call sfc_diff (im,Statein%pgr, Statein%ugrs, Statein%vgrs,& @@ -1235,11 +1237,15 @@ subroutine GFS_physics_driver & ! --- ... surface energy balance over ocean - call sfc_ocean & + ! kgao: this version gets hflx and evap over ocean points + ! based on shflx and lhflx from coupler + call sfc_ocean_gfdl_coupled & ! --- inputs: (im, Statein%pgr, Statein%ugrs, Statein%vgrs, Statein%tgrs, & Statein%qgrs, Sfcprop%tsfc, cd, cdq, Statein%prsl(1,1), & work3, islmsk, Tbd%phy_f2d(1,Model%num_p2d), flag_iter, & + ! kgao: shflx and lhflx from coupler + Sfcprop%shflx, Sfcprop%lhflx, & ! --- outputs: qss, Diag%cmm, Diag%chh, gflx, evap, hflx, ep1d) @@ -1399,6 +1405,10 @@ subroutine GFS_physics_driver & Diag%u1(:) = Statein%ugrs(:,1) Diag%v1(:) = Statein%vgrs(:,1) Sfcprop%qsfc(:) = qss(:) + + ! KGao + Diag%hflx(:) = hflx(:) + Diag%evap(:) = evap(:) ! --- ... update near surface fields diff --git a/GFS_layer/GFS_typedefs.F90 b/GFS_layer/GFS_typedefs.F90 index a478234b..ba73d869 100644 --- a/GFS_layer/GFS_typedefs.F90 +++ b/GFS_layer/GFS_typedefs.F90 @@ -224,6 +224,8 @@ module GFS_typedefs real (kind=kind_phys), pointer :: uustar (:) => null() !< boundary layer parameter real (kind=kind_phys), pointer :: oro (:) => null() !< orography real (kind=kind_phys), pointer :: oro_uf (:) => null() !< unfiltered orography + real (kind=kind_phys), pointer :: shflx (:) => null() !< sen heat flux (kgao) + real (kind=kind_phys), pointer :: lhflx (:) => null() !< latent heat flux (kgao) !--- IN/out MYJ scheme real (kind=kind_phys), pointer :: QZ0 (:) => null() !< vapor mixing ratio at z=z0 @@ -1219,8 +1221,10 @@ module GFS_typedefs real (kind=kind_phys), pointer :: totgrpb(:) => null() !< accumulated graupel precipitation in bucket (kg/m2) ! Output - only in physics - real (kind=kind_phys), pointer :: u10m (:) => null() !< 10 meater u/v wind speed - real (kind=kind_phys), pointer :: v10m (:) => null() !< 10 meater u/v wind speed + real (kind=kind_phys), pointer :: u10m (:) => null() !< 10 meter u/v wind speed + real (kind=kind_phys), pointer :: v10m (:) => null() !< 10 meter u/v wind speed + real (kind=kind_phys), pointer :: hflx (:) => null() !< sfc temp flux + real (kind=kind_phys), pointer :: evap (:) => null() !< sfc moisture flux real (kind=kind_phys), pointer :: dpt2m (:) => null() !< 2 meter dew point temperature real (kind=kind_phys), pointer :: zlvl (:) => null() !< layer 1 height (m) real (kind=kind_phys), pointer :: psurf (:) => null() !< surface pressure (Pa) @@ -1561,6 +1565,8 @@ subroutine sfcprop_create (Sfcprop, IM, Model) allocate (Sfcprop%uustar (IM)) allocate (Sfcprop%oro (IM)) allocate (Sfcprop%oro_uf (IM)) + allocate (Sfcprop%shflx (IM)) + allocate (Sfcprop%lhflx (IM)) Sfcprop%slope = clear_val Sfcprop%shdmin = clear_val @@ -1573,6 +1579,8 @@ subroutine sfcprop_create (Sfcprop, IM, Model) Sfcprop%uustar = clear_val Sfcprop%oro = clear_val Sfcprop%oro_uf = clear_val + Sfcprop%shflx = clear_val + Sfcprop%lhflx = clear_val if (Model%myj_pbl) then allocate (Sfcprop%QZ0 (IM)) @@ -3770,6 +3778,8 @@ subroutine diag_create (Diag, IM, Model) allocate (Diag%totgrpb (IM)) allocate (Diag%u10m (IM)) allocate (Diag%v10m (IM)) + allocate (Diag%hflx (IM)) + allocate (Diag%evap (IM)) allocate (Diag%dpt2m (IM)) allocate (Diag%zlvl (IM)) allocate (Diag%psurf (IM)) @@ -4037,6 +4047,8 @@ subroutine diag_phys_zero (Diag, Model, linit, iauwindow_center) !--- Out Diag%u10m = zero Diag%v10m = zero + Diag%hflx = zero + Diag%evap = zero Diag%dpt2m = zero Diag%zlvl = zero Diag%psurf = zero From 4cfb96a6bcb2d9a90c833444c820399aa30e8f09 Mon Sep 17 00:00:00 2001 From: Kun Gao Date: Mon, 25 Mar 2024 23:10:37 -0400 Subject: [PATCH 3/4] update SHiELD_physics to use sfc fluxes from coupler - part 3 --- FV3GFS/FV3GFS_io.F90 | 12 + GFS_layer/GFS_physics_driver.F90 | 33 +- GFS_layer/GFS_typedefs.F90 | 6 +- gsmphys/sfc_diff_coupled.f | 180 +++++++++ gsmphys/sfc_diff_gfdl_coupled.f | 365 ------------------ ...ean_gfdl_coupled.f => sfc_ocean_coupled.f} | 2 +- 6 files changed, 228 insertions(+), 370 deletions(-) create mode 100644 gsmphys/sfc_diff_coupled.f delete mode 100644 gsmphys/sfc_diff_gfdl_coupled.f rename gsmphys/{sfc_ocean_gfdl_coupled.f => sfc_ocean_coupled.f} (99%) diff --git a/FV3GFS/FV3GFS_io.F90 b/FV3GFS/FV3GFS_io.F90 index 0ccf8523..29a415a4 100644 --- a/FV3GFS/FV3GFS_io.F90 +++ b/FV3GFS/FV3GFS_io.F90 @@ -6722,6 +6722,18 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Cldprop, & Diag(idx)%data(nb)%var2 => Sfcprop(nb)%zorl(:) enddo + idx = idx + 1 + Diag(idx)%axes = 2 + Diag(idx)%name = 'ZTRLsfc' + Diag(idx)%desc = 'surface roughness for heat [m]' + Diag(idx)%unit = 'm' + Diag(idx)%mod_name = 'gfs_sfc' + Diag(idx)%cnvfac = cn_one/cn_100 + allocate (Diag(idx)%data(nblks)) + do nb = 1,nblks + Diag(idx)%data(nb)%var2 => Sfcprop(nb)%ztrl(:) + enddo + idx = idx + 1 Diag(idx)%axes = 2 Diag(idx)%name = 'VFRACsfc' diff --git a/GFS_layer/GFS_physics_driver.F90 b/GFS_layer/GFS_physics_driver.F90 index 01e06a95..d6cc3cfb 100644 --- a/GFS_layer/GFS_physics_driver.F90 +++ b/GFS_layer/GFS_physics_driver.F90 @@ -1135,11 +1135,12 @@ subroutine GFS_physics_driver & !else !!$ endif - if (Model%sfc_gfdl) then + ! kgao - need a logic to ensure sfc_coupld is true when coupled + if (Model%sfc_coupled) then ! a version of sfc_diff from coupling with MOM6 by kgao ! Sfcprop%uustar,Sfcprop%zorl not updated over ocean points ! needs to consider ztrl too - call sfc_diff_gfdl_coupled(im,Statein%pgr, Statein%ugrs, Statein%vgrs,& + call sfc_diff_coupled(im,Statein%pgr, Statein%ugrs, Statein%vgrs,& Statein%tgrs, Statein%qgrs, Diag%zlvl, Sfcprop%snowd, & Sfcprop%tsfc, Sfcprop%zorl, Sfcprop%ztrl, cd, & cdq, rb, Statein%prsl(1,1), work3, islmsk, stress, & @@ -1150,6 +1151,21 @@ subroutine GFS_physics_driver & !Model%do_z0_moon, Model%do_z0_hwrf15, & !Model%do_z0_hwrf17, Model%do_z0_hwrf17_hwonly, & !Model%wind_th_hwrf) + + else if (Model%sfc_gfdl) then +! a new and more flexible version of sfc_diff by kgao + call sfc_diff_gfdl(im,Statein%pgr, Statein%ugrs, Statein%vgrs,& + Statein%tgrs, Statein%qgrs, Diag%zlvl, Sfcprop%snowd, & + Sfcprop%tsfc, Sfcprop%zorl, Sfcprop%ztrl, cd, & + cdq, rb, Statein%prsl(1,1), work3, islmsk, stress, & + Sfcprop%ffmm, Sfcprop%ffhh, Sfcprop%uustar, & + wind, Tbd%phy_f2d(1,Model%num_p2d), fm10, fh2, & + sigmaf, vegtype, Sfcprop%shdmax, Model%ivegsrc, & + tsurf, flag_iter, Model%redrag, Model%z0s_max, & + Model%do_z0_moon, Model%do_z0_hwrf15, & + Model%do_z0_hwrf17, Model%do_z0_hwrf17_hwonly, & + Model%wind_th_hwrf) + else ! GFS original sfc_diff modified by kgao call sfc_diff (im,Statein%pgr, Statein%ugrs, Statein%vgrs,& @@ -1237,9 +1253,10 @@ subroutine GFS_physics_driver & ! --- ... surface energy balance over ocean + if (Model%sfc_coupled) then ! kgao: this version gets hflx and evap over ocean points ! based on shflx and lhflx from coupler - call sfc_ocean_gfdl_coupled & + call sfc_ocean_coupled & ! --- inputs: (im, Statein%pgr, Statein%ugrs, Statein%vgrs, Statein%tgrs, & Statein%qgrs, Sfcprop%tsfc, cd, cdq, Statein%prsl(1,1), & @@ -1249,6 +1266,16 @@ subroutine GFS_physics_driver & ! --- outputs: qss, Diag%cmm, Diag%chh, gflx, evap, hflx, ep1d) + else + call sfc_ocean & +! --- inputs: + (im, Statein%pgr, Statein%ugrs, Statein%vgrs, Statein%tgrs, & + Statein%qgrs, Sfcprop%tsfc, cd, cdq, Statein%prsl(1,1), & + work3, islmsk, Tbd%phy_f2d(1,Model%num_p2d), flag_iter, & +! --- outputs: + qss, Diag%cmm, Diag%chh, gflx, evap, hflx, ep1d) + endif + endif ! if ( nstf_name(1) > 0 ) then ! if (lprnt) write(0,*)' sfalb=',sfalb(ipr),' ipr=',ipr & diff --git a/GFS_layer/GFS_typedefs.F90 b/GFS_layer/GFS_typedefs.F90 index ba73d869..224c6f77 100644 --- a/GFS_layer/GFS_typedefs.F90 +++ b/GFS_layer/GFS_typedefs.F90 @@ -623,6 +623,7 @@ module GFS_typedefs logical :: shcnvcw !< flag for shallow convective cloud logical :: redrag !< flag for reduced drag coeff. over sea logical :: sfc_gfdl !< flag for using updated sfc layer scheme + logical :: sfc_coupled !< flag for using sfc layer scheme designed for coupling real(kind=kind_phys) :: z0s_max !< a limiting value for z0 under high winds logical :: do_z0_moon !< flag for using z0 scheme in Moon et al. 2007 (kgao) logical :: do_z0_hwrf15 !< flag for using z0 scheme in 2015 HWRF (kgao) @@ -2191,6 +2192,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & logical :: shcnvcw = .false. !< flag for shallow convective cloud logical :: redrag = .false. !< flag for reduced drag coeff. over sea logical :: sfc_gfdl = .false. !< flag for using new sfc layer scheme by kgao at GFDL + logical :: sfc_coupled = .false. !< flag for using sfc layer scheme designed for coupling real(kind=kind_phys) :: z0s_max = .317e-2 !< a limiting value for z0 under high winds logical :: do_z0_moon = .false. !< flag for using z0 scheme in Moon et al. 2007 logical :: do_z0_hwrf15 = .false. !< flag for using z0 scheme in 2015 HWRF @@ -2421,7 +2423,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & ras, trans_trac, old_monin, cnvgwd, mstrat, moist_adj, & cscnv, cal_pre, do_aw, do_shoc, shocaftcnv, shoc_cld, & h2o_phys, pdfcld, shcnvcw, redrag, sfc_gfdl, z0s_max, & - do_z0_moon, do_z0_hwrf15, do_z0_hwrf17, & + sfc_coupled, do_z0_moon, do_z0_hwrf15, do_z0_hwrf17, & do_z0_hwrf17_hwonly, wind_th_hwrf, & hybedmf, dspheat, lheatstrg, hour_canopy, afac_canopy, & cnvcld, no_pbl, xkzm_lim, xkzm_fac, xkgdx, & @@ -2645,6 +2647,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%shcnvcw = shcnvcw Model%redrag = redrag Model%sfc_gfdl = sfc_gfdl + Model%sfc_coupled = sfc_coupled Model%z0s_max = z0s_max Model%do_z0_moon = do_z0_moon Model%do_z0_hwrf15 = do_z0_hwrf15 @@ -3337,6 +3340,7 @@ subroutine control_print(Model) print *, ' shcnvcw : ', Model%shcnvcw print *, ' redrag : ', Model%redrag print *, ' sfc_gfdl : ', Model%sfc_gfdl + print *, ' sfc_coupled : ', Model%sfc_coupled print *, ' z0s_max : ', Model%z0s_max print *, ' do_z0_moon : ', Model%do_z0_moon print *, ' do_z0_hwrf15 : ', Model%do_z0_hwrf15 diff --git a/gsmphys/sfc_diff_coupled.f b/gsmphys/sfc_diff_coupled.f new file mode 100644 index 00000000..dfe08bcb --- /dev/null +++ b/gsmphys/sfc_diff_coupled.f @@ -0,0 +1,180 @@ + subroutine sfc_diff_coupled(im,ps,u1,v1,t1,q1,z1, + & snwdph,tskin,z0rl,ztrl,cm,ch,rb, + & prsl1,prslki,islimsk, + & stress,fm,fh, + & ustar,wind,ddvel,fm10,fh2, + & sigmaf,vegtype,shdmax,ivegsrc, + & tsurf,flag_iter) !,redrag, +! & z0s_max) +! & do_z0_moon, do_z0_hwrf15, do_z0_hwrf17, +! & do_z0_hwrf17_hwonly, wind_th_hwrf) + + use machine , only : kind_phys + use funcphys, only : fpvs + use physcons, grav => con_g, cp => con_cp + &, rvrdm1 => con_fvirt, rd => con_rd + &, eps => con_eps, epsm1 => con_epsm1 + + implicit none + +! --- input/output + + integer im, ivegsrc + + real(kind=kind_phys), dimension(im)::ps, u1, v1, t1, q1, z1 + &, tskin, z0rl, ztrl, cm, ch, rb + &, prsl1, prslki, stress + &, fm, fh, ustar, wind, ddvel + &, fm10, fh2, sigmaf, shdmax + &, tsurf, snwdph + integer, dimension(im) ::vegtype, islimsk + + logical flag_iter(im) +! logical redrag +! logical do_z0_moon, do_z0_hwrf15, do_z0_hwrf17 ! kgao +! &, do_z0_hwrf17_hwonly ! kgao + +! --- local + + integer i +! + real(kind=kind_phys) aa, aa0, bb, bb0, dtv, adtv,qs1, + & hl1, hl12, pm, ph, pm10, ph2, rat, + & thv1, tvs, z1i, z0, zt, z0max, ztmax, + & fms, fhs, hl0, hl0inf, hlinf, + & hl110, hlt, hltinf, olinf, + & restar, czilc, tem1, tem2, + & u10m, v10m, ws10m, ws10m_moon, !kgao + & z0_1, zt_1, fm1, fh1, ustar_1, ztmax_1 !kgao +! + +! real(kind=kind_phys),intent(in ) :: z0s_max, wind_th_hwrf ! kgao + + real(kind=kind_phys), parameter :: + & charnock=.014, ca=.4 + &, vis=1.4e-5, rnu=1.51e-5, visi=1.0/vis + &, log01=log(0.01), log05=log(0.05), log07=log(0.07) + &, ztmin1=-999.0 + +!================================================ +! Main program starts here +!================================================ + + do i=1,im + + if(flag_iter(i)) then + +! --- get variables at model lowest layer and surface (water/ice/land) + + wind(i) = max(sqrt(u1(i)*u1(i) + v1(i)*v1(i)) + & + max(0.0, min(ddvel(i), 30.0)), 1.0) + tem1 = 1.0 + rvrdm1 * max(q1(i),1.e-8) + thv1 = t1(i) * prslki(i) * tem1 + tvs = 0.5 * (tsurf(i)+tskin(i)) * tem1 + qs1 = fpvs(t1(i)) + qs1 = max(1.0e-8, eps * qs1 / (prsl1(i) + epsm1 * qs1)) + + !(sea/land/ice mask =0/1/2) + if(islimsk(i) == 1 .or. islimsk(i) == 2) then ! over land or sea ice + +!================================================ +! if over land or sea ice: +! step 1 - get z0/zt +! step 2 - call similarity +!================================================ + +! --- get surface roughness for momentum (z0) + + z0 = 0.01 * z0rl(i) + z0max = max(1.0e-6, min(z0,z1(i))) + + !xubin's new z0 over land and sea ice + tem1 = 1.0 - shdmax(i) ! shdmax is max vegetation area fraction + tem2 = tem1 * tem1 + tem1 = 1.0 - tem2 + + if( ivegsrc == 1 ) then + + if (vegtype(i) == 10) then + z0max = exp( tem2*log01 + tem1*log07 ) + elseif (vegtype(i) == 6) then + z0max = exp( tem2*log01 + tem1*log05 ) + elseif (vegtype(i) == 7) then +! z0max = exp( tem2*log01 + tem1*log01 ) + z0max = 0.01 + elseif (vegtype(i) == 16) then +! z0max = exp( tem2*log01 + tem1*log01 ) + z0max = 0.01 + else + z0max = exp( tem2*log01 + tem1*log(z0max) ) + endif + + elseif (ivegsrc == 2 ) then + + if (vegtype(i) == 7) then + z0max = exp( tem2*log01 + tem1*log07 ) + elseif (vegtype(i) == 8) then + z0max = exp( tem2*log01 + tem1*log05 ) + elseif (vegtype(i) == 9) then +! z0max = exp( tem2*log01 + tem1*log01 ) + z0max = 0.01 + elseif (vegtype(i) == 11) then +! z0max = exp( tem2*log01 + tem1*log01 ) + z0max = 0.01 + else + z0max = exp( tem2*log01 + tem1*log(z0max) ) + endif + + z0max = max(z0max,1.0e-6) + + endif + +! --- get surface roughness for heat (zt) + +! czilc = 10.0 ** (- (0.40/0.07) * z0) ! let czilc depend on canopy height + czilc = 0.8 + + tem1 = 1.0 - sigmaf(i) + ztmax = z0max*exp( - tem1*tem1 + & * czilc*ca*sqrt(ustar(i)*(0.01/1.5e-05))) + + ztmax = max(ztmax,1.0e-6) + +! --- call similarity + + call monin_obukhov_similarity + & (z1(i), snwdph(i), thv1, wind(i), z0max, ztmax, tvs, + & rb(i), fm(i), fh(i), fm10(i), fh2(i), + & cm(i), ch(i), stress(i), ustar(i)) + + elseif (islimsk(i) == 0) then + +!================================================ +! if over water +! - redesigned by Kun Gao for coupling with MOM6 +! - not updating z0, zt and ustar over ocean +!================================================ + + ! --- z0/zt from coupler + z0 = 0.01 * z0rl(i) + zt = 0.01 * ztrl(i) + + z0max = max(1.0e-6, min(z0,z1(i))) + ztmax = max(zt,1.0e-6) + + ! --- call similarity + ! kgao: use z0/zt to do sfc diagnosis + call monin_obukhov_similarity + & (z1(i), snwdph(i), thv1, wind(i), z0max, ztmax, tvs, + & rb(i), fm(i), fh(i), fm10(i), fh2(i), + & cm(i), ch(i), tem1, tem2) !stress(i), ustar(i)) + + ! kgao: use ustar from coupler to get stress + stress(i) = ustar(i) * ustar(i) + + endif ! end of if(islimsk) loop + endif ! end of if(flagiter) loop + enddo ! end of do i=1,im loop + + return + end subroutine diff --git a/gsmphys/sfc_diff_gfdl_coupled.f b/gsmphys/sfc_diff_gfdl_coupled.f deleted file mode 100644 index fab6f752..00000000 --- a/gsmphys/sfc_diff_gfdl_coupled.f +++ /dev/null @@ -1,365 +0,0 @@ - subroutine sfc_diff_gfdl_coupled(im,ps,u1,v1,t1,q1,z1, - & snwdph,tskin,z0rl,ztrl,cm,ch,rb, - & prsl1,prslki,islimsk, - & stress,fm,fh, - & ustar,wind,ddvel,fm10,fh2, - & sigmaf,vegtype,shdmax,ivegsrc, - & tsurf,flag_iter) !,redrag, -! & z0s_max) -! & do_z0_moon, do_z0_hwrf15, do_z0_hwrf17, -! & do_z0_hwrf17_hwonly, wind_th_hwrf) - - use machine , only : kind_phys - use funcphys, only : fpvs - use physcons, grav => con_g, cp => con_cp - &, rvrdm1 => con_fvirt, rd => con_rd - &, eps => con_eps, epsm1 => con_epsm1 - - implicit none - -! --- input/output - - integer im, ivegsrc - - real(kind=kind_phys), dimension(im)::ps, u1, v1, t1, q1, z1 - &, tskin, z0rl, ztrl, cm, ch, rb - &, prsl1, prslki, stress - &, fm, fh, ustar, wind, ddvel - &, fm10, fh2, sigmaf, shdmax - &, tsurf, snwdph - integer, dimension(im) ::vegtype, islimsk - - logical flag_iter(im) -! logical redrag -! logical do_z0_moon, do_z0_hwrf15, do_z0_hwrf17 ! kgao -! &, do_z0_hwrf17_hwonly ! kgao - -! --- local - - integer i -! - real(kind=kind_phys) aa, aa0, bb, bb0, dtv, adtv,qs1, - & hl1, hl12, pm, ph, pm10, ph2, rat, - & thv1, tvs, z1i, z0, zt, z0max, ztmax, - & fms, fhs, hl0, hl0inf, hlinf, - & hl110, hlt, hltinf, olinf, - & restar, czilc, tem1, tem2, - & u10m, v10m, ws10m, ws10m_moon, !kgao - & z0_1, zt_1, fm1, fh1, ustar_1, ztmax_1 !kgao -! - -! real(kind=kind_phys),intent(in ) :: z0s_max, wind_th_hwrf ! kgao - - real(kind=kind_phys), parameter :: - & charnock=.014, ca=.4 - &, vis=1.4e-5, rnu=1.51e-5, visi=1.0/vis - &, log01=log(0.01), log05=log(0.05), log07=log(0.07) - &, ztmin1=-999.0 - -!================================================ -! Main program starts here -!================================================ - - do i=1,im - - if(flag_iter(i)) then - -! --- get variables at model lowest layer and surface (water/ice/land) - - wind(i) = max(sqrt(u1(i)*u1(i) + v1(i)*v1(i)) - & + max(0.0, min(ddvel(i), 30.0)), 1.0) - tem1 = 1.0 + rvrdm1 * max(q1(i),1.e-8) - thv1 = t1(i) * prslki(i) * tem1 - tvs = 0.5 * (tsurf(i)+tskin(i)) * tem1 - qs1 = fpvs(t1(i)) - qs1 = max(1.0e-8, eps * qs1 / (prsl1(i) + epsm1 * qs1)) - - !(sea/land/ice mask =0/1/2) - if(islimsk(i) == 1 .or. islimsk(i) == 2) then ! over land or sea ice - -!================================================ -! if over land or sea ice: -! step 1 - get z0/zt -! step 2 - call similarity -!================================================ - -! --- get surface roughness for momentum (z0) - - z0 = 0.01 * z0rl(i) - z0max = max(1.0e-6, min(z0,z1(i))) - - !xubin's new z0 over land and sea ice - tem1 = 1.0 - shdmax(i) ! shdmax is max vegetation area fraction - tem2 = tem1 * tem1 - tem1 = 1.0 - tem2 - - if( ivegsrc == 1 ) then - - if (vegtype(i) == 10) then - z0max = exp( tem2*log01 + tem1*log07 ) - elseif (vegtype(i) == 6) then - z0max = exp( tem2*log01 + tem1*log05 ) - elseif (vegtype(i) == 7) then -! z0max = exp( tem2*log01 + tem1*log01 ) - z0max = 0.01 - elseif (vegtype(i) == 16) then -! z0max = exp( tem2*log01 + tem1*log01 ) - z0max = 0.01 - else - z0max = exp( tem2*log01 + tem1*log(z0max) ) - endif - - elseif (ivegsrc == 2 ) then - - if (vegtype(i) == 7) then - z0max = exp( tem2*log01 + tem1*log07 ) - elseif (vegtype(i) == 8) then - z0max = exp( tem2*log01 + tem1*log05 ) - elseif (vegtype(i) == 9) then -! z0max = exp( tem2*log01 + tem1*log01 ) - z0max = 0.01 - elseif (vegtype(i) == 11) then -! z0max = exp( tem2*log01 + tem1*log01 ) - z0max = 0.01 - else - z0max = exp( tem2*log01 + tem1*log(z0max) ) - endif - - z0max = max(z0max,1.0e-6) - - endif - -! --- get surface roughness for heat (zt) - -! czilc = 10.0 ** (- (0.40/0.07) * z0) ! let czilc depend on canopy height - czilc = 0.8 - - tem1 = 1.0 - sigmaf(i) - ztmax = z0max*exp( - tem1*tem1 - & * czilc*ca*sqrt(ustar(i)*(0.01/1.5e-05))) - - ztmax = max(ztmax,1.0e-6) - -! --- call similarity - - call monin_obukhov_similarity - & (z1(i), snwdph(i), thv1, wind(i), z0max, ztmax, tvs, - & rb(i), fm(i), fh(i), fm10(i), fh2(i), - & cm(i), ch(i), stress(i), ustar(i)) - - elseif (islimsk(i) == 0) then - -!================================================ -! if over water -! redesigned by Kun Gao for coupling with MOM6 -! now do not update z0rl and ustar over ocean points -!================================================ - - ! --- get z0/zt - z0 = 0.01 * z0rl(i) ! kgao: from coupler - zt = 0.01 * ztrl(i) - - z0max = max(1.0e-6, min(z0,z1(i))) - ztmax = max(zt,1.0e-6) - - ! --- call similarity - call monin_obukhov_similarity - & (z1(i), snwdph(i), thv1, wind(i), z0max, ztmax, tvs, - & rb(i), fm(i), fh(i), fm10(i), fh2(i), - & cm(i), ch(i), tem1, tem2) !stress(i), ustar(i)) - - ! kgao: use ustar from coupler to update stress - stress(i) = ustar(i) * ustar(i) - - ! kgao: get zt as in old sfc_diff.f - ! this should be removed once using ztrl from coupler - ustar_1 = sqrt(grav * z0max / charnock) - restar = max(ustar_1*z0max*visi, 0.000001) - rat = min(7.0, 2.67 * sqrt(sqrt(restar)) - 2.57) - zt = z0max * exp(-rat) ! zeng, zhao and dickinson 1997 (eq 25) - - ! kgao: do not update z0 - !z0rl(i) = 100.0 * z0max - ztrl(i) = 100.0 * ztmax - - endif ! end of if(islimsk) loop - endif ! end of if(flagiter) loop - enddo ! end of do i=1,im loop - - return - end subroutine - -! ======================================================================= - - subroutine monin_obukhov_similarity - & ( z1, snwdph, thv1, wind, z0max, ztmax, tvs, - & rb, fm, fh, fm10, fh2, cm, ch, stress, ustar) - -! --- input -! z1 - lowest model level height -! snwdph - surface snow thickness -! wind - wind speed at lowest model layer -! thv1 - virtual potential temp at lowest model layer -! tvs - surface temp -! z0max - surface roughness length for momentum -! ztmax - surface roughness length for heat -! -! --- output -! rb - a bulk richardson number -! fm, fh - similarity function defined at lowest model layer -! fm10, fh2 - similarity function defined at 10m (for momentum) and 2m (for heat) -! cm, ch - surface exchange coefficients for momentum and heat -! stress - surface wind stress -! ustar - surface frictional velocity - - use machine , only : kind_phys - use physcons, grav => con_g - -! --- inputs: - real(kind=kind_phys), intent(in) :: - & z1, snwdph, thv1, wind, z0max, ztmax, tvs - -! --- outputs: - real(kind=kind_phys), intent(out) :: - & rb, fm, fh, fm10, fh2, cm, ch, stress, ustar - -! --- locals: - - real(kind=kind_phys), parameter :: alpha=5., a0=-3.975 - &, a1=12.32, alpha4=4.0*alpha - &, b1=-7.755, b2=6.041, alpha2=alpha+alpha, beta=1.0 - &, a0p=-7.941, a1p=24.75, b1p=-8.705, b2p=7.899 - &, ztmin1=-999.0, ca=.4 - - real(kind=kind_phys) aa, aa0, bb, bb0, dtv, adtv, - & hl1, hl12, pm, ph, pm10, ph2, - & z1i, - & fms, fhs, hl0, hl0inf, hlinf, - & hl110, hlt, hltinf, olinf, - & tem1, tem2, ztmax1 - - z1i = 1.0 / z1 - - tem1 = z0max/z1 - if (abs(1.0-tem1) > 1.0e-6) then - ztmax1 = - beta*log(tem1)/(alpha2*(1.-tem1)) - else - ztmax1 = 99.0 - endif - if( z0max < 0.05 .and. snwdph < 10.0 ) ztmax1 = 99.0 - -! -! compute stability indices (rb and hlinf) -! - dtv = thv1 - tvs - adtv = max(abs(dtv),0.001) - dtv = sign(1.,dtv) * adtv - rb = max(-5000.0, (grav+grav) * dtv * z1 - & / ((thv1 + tvs) * wind * wind)) - tem1 = 1.0 / z0max - tem2 = 1.0 / ztmax - fm = log((z0max+z1) * tem1) - fh = log((ztmax+z1) * tem2) - fm10 = log((z0max+10.) * tem1) - fh2 = log((ztmax+2.) * tem2) - hlinf = rb * fm * fm / fh - hlinf = min(max(hlinf,ztmin1),ztmax1) -! -! stable case -! - if (dtv >= 0.0) then - hl1 = hlinf - if(hlinf > .25) then - tem1 = hlinf * z1i - hl0inf = z0max * tem1 - hltinf = ztmax * tem1 - aa = sqrt(1. + alpha4 * hlinf) - aa0 = sqrt(1. + alpha4 * hl0inf) - bb = aa - bb0 = sqrt(1. + alpha4 * hltinf) - pm = aa0 - aa + log( (aa + 1.)/(aa0 + 1.) ) - ph = bb0 - bb + log( (bb + 1.)/(bb0 + 1.) ) - fms = fm - pm - fhs = fh - ph - hl1 = fms * fms * rb / fhs - hl1 = min(max(hl1, ztmin1), ztmax1) - endif -! -! second iteration -! - tem1 = hl1 * z1i - hl0 = z0max * tem1 - hlt = ztmax * tem1 - aa = sqrt(1. + alpha4 * hl1) - aa0 = sqrt(1. + alpha4 * hl0) - bb = aa - bb0 = sqrt(1. + alpha4 * hlt) - pm = aa0 - aa + log( (1.0+aa)/(1.0+aa0) ) - ph = bb0 - bb + log( (1.0+bb)/(1.0+bb0) ) - hl110 = hl1 * 10. * z1i - hl110 = min(max(hl110, ztmin1), ztmax1) - aa = sqrt(1. + alpha4 * hl110) - pm10 = aa0 - aa + log( (1.0+aa)/(1.0+aa0) ) - hl12 = (hl1+hl1) * z1i - hl12 = min(max(hl12,ztmin1),ztmax1) -! aa = sqrt(1. + alpha4 * hl12) - bb = sqrt(1. + alpha4 * hl12) - ph2 = bb0 - bb + log( (1.0+bb)/(1.0+bb0) ) -! -! unstable case - check for unphysical obukhov length -! - else ! dtv < 0 case - olinf = z1 / hlinf - tem1 = 50.0 * z0max - if(abs(olinf) <= tem1) then - hlinf = -z1 / tem1 - hlinf = min(max(hlinf,ztmin1),ztmax1) - endif -! -! get pm and ph -! - if (hlinf >= -0.5) then - hl1 = hlinf - pm = (a0 + a1*hl1) * hl1 / (1.+ (b1+b2*hl1) *hl1) - ph = (a0p + a1p*hl1) * hl1 / (1.+ (b1p+b2p*hl1)*hl1) - hl110 = hl1 * 10. * z1i - hl110 = min(max(hl110, ztmin1), ztmax1) - pm10 = (a0 + a1*hl110) * hl110 / (1.+(b1+b2*hl110)*hl110) - hl12 = (hl1+hl1) * z1i - hl12 = min(max(hl12, ztmin1), ztmax1) - ph2 = (a0p + a1p*hl12) * hl12 / (1.+(b1p+b2p*hl12)*hl12) - else ! hlinf < 0.05 - hl1 = -hlinf - tem1 = 1.0 / sqrt(hl1) - pm = log(hl1) + 2. * sqrt(tem1) - .8776 - ph = log(hl1) + .5 * tem1 + 1.386 -! pm = log(hl1) + 2.0 * hl1 ** (-.25) - .8776 -! ph = log(hl1) + 0.5 * hl1 ** (-.5) + 1.386 - hl110 = hl1 * 10. * z1i - hl110 = min(max(hl110, ztmin1), ztmax1) - pm10 = log(hl110) + 2.0 / sqrt(sqrt(hl110)) - .8776 -! pm10 = log(hl110) + 2. * hl110 ** (-.25) - .8776 - hl12 = (hl1+hl1) * z1i - hl12 = min(max(hl12, ztmin1), ztmax1) - ph2 = log(hl12) + 0.5 / sqrt(hl12) + 1.386 -! ph2 = log(hl12) + .5 * hl12 ** (-.5) + 1.386 - endif - - endif ! end of if (dtv >= 0 ) then loop -! -! finish the exchange coefficient computation to provide fm and fh -! - fm = fm - pm - fh = fh - ph - fm10 = fm10 - pm10 - fh2 = fh2 - ph2 - cm = ca * ca / (fm * fm) - ch = ca * ca / (fm * fh) - tem1 = 0.00001/z1 - cm = max(cm, tem1) - ch = max(ch, tem1) - stress = cm * wind * wind - ustar = sqrt(stress) - - return - end subroutine monin_obukhov_similarity diff --git a/gsmphys/sfc_ocean_gfdl_coupled.f b/gsmphys/sfc_ocean_coupled.f similarity index 99% rename from gsmphys/sfc_ocean_gfdl_coupled.f rename to gsmphys/sfc_ocean_coupled.f index c46883a3..1cf8d5f8 100644 --- a/gsmphys/sfc_ocean_gfdl_coupled.f +++ b/gsmphys/sfc_ocean_coupled.f @@ -1,5 +1,5 @@ !----------------------------------- - subroutine sfc_ocean_gfdl_coupled & + subroutine sfc_ocean_coupled & !................................... ! --- inputs: & ( im, ps, u1, v1, t1, q1, tskin, cm, ch, & From 20019a92c8a0cf2aa9ecf7dd71c638bb54413ef2 Mon Sep 17 00:00:00 2001 From: Kun Gao Date: Mon, 1 Apr 2024 22:10:17 -0400 Subject: [PATCH 4/4] update comments --- GFS_layer/GFS_physics_driver.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/GFS_layer/GFS_physics_driver.F90 b/GFS_layer/GFS_physics_driver.F90 index d6cc3cfb..887ad770 100644 --- a/GFS_layer/GFS_physics_driver.F90 +++ b/GFS_layer/GFS_physics_driver.F90 @@ -1135,11 +1135,10 @@ subroutine GFS_physics_driver & !else !!$ endif - ! kgao - need a logic to ensure sfc_coupld is true when coupled + ! kgao - need a logic to ensure sfc_coupled is true when coupled with MOM6 if (Model%sfc_coupled) then ! a version of sfc_diff from coupling with MOM6 by kgao -! Sfcprop%uustar,Sfcprop%zorl not updated over ocean points -! needs to consider ztrl too +! Sfcprop%uustar,Sfcprop%zorl,Sfcprop%ztrl are not updated over ocean points call sfc_diff_coupled(im,Statein%pgr, Statein%ugrs, Statein%vgrs,& Statein%tgrs, Statein%qgrs, Diag%zlvl, Sfcprop%snowd, & Sfcprop%tsfc, Sfcprop%zorl, Sfcprop%ztrl, cd, & @@ -1254,7 +1253,8 @@ subroutine GFS_physics_driver & ! --- ... surface energy balance over ocean if (Model%sfc_coupled) then - ! kgao: this version gets hflx and evap over ocean points + ! kgao: this version is for coupling with MOM6, which + ! gets hflx and evap over ocean points ! based on shflx and lhflx from coupler call sfc_ocean_coupled & ! --- inputs: