Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
… into HEAD
  • Loading branch information
climbfuji committed Jul 24, 2024
2 parents c866747 + 2a50ccc commit 9e736da
Show file tree
Hide file tree
Showing 9 changed files with 94 additions and 87 deletions.
28 changes: 20 additions & 8 deletions physics/CONV/SAMF/samfdeepcnv.f
Original file line number Diff line number Diff line change
Expand Up @@ -215,7 +215,8 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
!
! parameters for prognostic sigma closure
real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km),
& omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im),qadv(im,km)
& omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im),qadv(im,km),
& sigmaoutx(im)
real(kind=kind_phys) gravinv,invdelt,sigmind,sigminm,sigmins
parameter(sigmind=0.01,sigmins=0.03,sigminm=0.01)
logical flag_shallow, flag_mid
Expand Down Expand Up @@ -3423,17 +3424,28 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
endif
enddo
c
c convective cloud water
!
if(progsigma)then
do i = 1, im
sigmaoutx(i)=max(sigmaout(i,1),0.0)
sigmaoutx(i)=min(sigmaoutx(i),1.0)
enddo
endif
c
!> - Calculate convective cloud water.
do k = 1, km
do i = 1, im
if (cnvflg(i) .and. rn(i) > 0.) then
if (k >= kbcon(i) .and. k < ktcon(i)) then
cnvw(i,k) = cnvwt(i,k) * xmb(i) * dt2
do i = 1, im
if (cnvflg(i) .and. rn(i) > 0.) then
if (k >= kbcon(i) .and. k < ktcon(i)) then
cnvw(i,k) = cnvwt(i,k) * xmb(i) * dt2
if(progsigma)then
cnvw(i,k) = cnvw(i,k) * sigmaoutx(i)
else
cnvw(i,k) = cnvw(i,k) * sigmagfm(i)
endif
endif
endif
endif
enddo
enddo
enddo
c
c convective cloud cover
Expand Down
33 changes: 21 additions & 12 deletions physics/CONV/SAMF/samfshalcnv.f
Original file line number Diff line number Diff line change
Expand Up @@ -162,7 +162,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
! parameters for prognostic sigma closure
real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km),
& omegac(im),zeta(im,km),dbyo1(im,km),
& sigmab(im),qadv(im,km)
& sigmab(im),qadv(im,km),sigmaoutx(im)
real(kind=kind_phys) gravinv,dxcrtas,invdelt,sigmind,sigmins,
& sigminm
logical flag_shallow,flag_mid
Expand Down Expand Up @@ -2397,20 +2397,29 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
endif
enddo
c
c convective cloud water
c
!> - Calculate shallow convective cloud water.
if(progsigma)then
do i = 1, im
sigmaoutx(i)=max(sigmaout(i,1),0.0)
sigmaoutx(i)=min(sigmaoutx(i),1.0)
enddo
endif

c convective cloud water
do k = 1, km
do i = 1, im
if (cnvflg(i)) then
if (k >= kbcon(i) .and. k < ktcon(i)) then
cnvw(i,k) = cnvwt(i,k) * xmb(i) * dt2
do i = 1, im
if (cnvflg(i)) then
if (k >= kbcon(i) .and. k < ktcon(i)) then
cnvw(i,k) = cnvwt(i,k) * xmb(i) * dt2
if (progsigma) then
cnvw(i,k) = cnvw(i,k) * sigmaoutx(i)
else
cnvw(i,k) = cnvw(i,k) * sigmagfm(i)
endif
endif
endif
endif
enddo
enddo
enddo

c
c
c convective cloud cover
c
!> - Calculate convective cloud cover, which is used when pdf-based cloud fraction is used (i.e., pdfcld=.true.).
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -593,8 +593,10 @@ subroutine GFS_phys_time_vary_init (

isnow = nint(snowxy(ix))+1 ! snowxy <=0.0, dzsno >= 0.0

! using stc and tgxy to linearly interpolate the snow temp for each layer

do is = isnow,0
tsnoxy(ix,is) = tgxy(ix)
tsnoxy(ix,is) = tgxy(ix) + (( sum(dzsno(isnow:is)) -0.5*dzsno(is) )/snd)*(stc(ix,1)-tgxy(ix))
snliqxy(ix,is) = zero
snicexy(ix,is) = one * dzsno(is) * weasd(ix)/snd
enddo
Expand Down
2 changes: 1 addition & 1 deletion physics/PBL/SATMEDMF/satmedmfvdifq.F
Original file line number Diff line number Diff line change
Expand Up @@ -275,7 +275,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
parameter(elmfac=1.0,elefac=1.0,cql=100.)
parameter(dw2min=1.e-4,dkmax=1000.,xkgdx=1000.)
parameter(qlcr=3.5e-5,zstblmax=2500.)
parameter(xkinv1=0.4,xkinv2=0.3)
parameter(xkinv1=0.15,xkinv2=0.3)
parameter(h1=0.33333333,hcrinv=250.)
parameter(vegflo=0.1,vegfup=1.0,z0lo=0.1,z0up=1.0)
parameter(vc0=1.0,zc0=1.0)
Expand Down
14 changes: 1 addition & 13 deletions physics/SFC_Layer/UFS/sfc_diag_post.F90
Original file line number Diff line number Diff line change
Expand Up @@ -14,15 +14,14 @@ module sfc_diag_post
!!
#endif
subroutine sfc_diag_post_run (im, lsm, lsm_noahmp, opt_diag, dry, lssav, dtf, con_eps, con_epsm1, pgr,&
vegtype,t2mmp,q2mp, t2m, q2m, u10m, v10m, tmpmin, tmpmax, spfhmin, spfhmax, &
t2mmp,q2mp, t2m, q2m, u10m, v10m, tmpmin, tmpmax, spfhmin, spfhmax, &
wind10mmax, u10mmax, v10mmax, dpt2m, errmsg, errflg)

use machine, only: kind_phys, kind_dbl_prec

implicit none

integer, intent(in) :: im, lsm, lsm_noahmp,opt_diag
integer, dimension(:), intent(in) :: vegtype ! vegetation type (integer index)
logical, intent(in) :: lssav
real(kind=kind_phys), intent(in) :: dtf, con_eps, con_epsm1
logical , dimension(:), intent(in) :: dry
Expand All @@ -42,17 +41,6 @@ subroutine sfc_diag_post_run (im, lsm, lsm_noahmp, opt_diag, dry, lssav, dtf, co
errflg = 0

if (lsm == lsm_noahmp) then
! over shrublands use opt_diag=2
do i=1, im
if(dry(i)) then
if (vegtype(i) == 6 .or. vegtype(i) == 7 &
.or. vegtype(i) == 16) then
t2m(i) = t2mmp(i)
q2m(i) = q2mp(i)
endif
endif
enddo

if (opt_diag == 2 .or. opt_diag == 3) then
do i=1,im
if(dry(i)) then
Expand Down
7 changes: 0 additions & 7 deletions physics/SFC_Layer/UFS/sfc_diag_post.meta
Original file line number Diff line number Diff line change
Expand Up @@ -82,13 +82,6 @@
type = real
kind = kind_phys
intent = in
[vegtype]
standard_name = vegetation_type_classification
long_name = vegetation type at each grid cell
units = index
dimensions = (horizontal_loop_extent)
type = integer
intent= in
[t2mmp]
standard_name = temperature_at_2m_from_noahmp
long_name = 2 meter temperature from noahmp
Expand Down
50 changes: 24 additions & 26 deletions physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90
Original file line number Diff line number Diff line change
Expand Up @@ -1989,7 +1989,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in

real (kind=kind_phys), parameter :: mpe = 1.e-6
real (kind=kind_phys), parameter :: psiwlt = -150. !metric potential for wilting point (m)
real (kind=kind_phys), parameter :: z0 = 0.002 ! bare-soil roughness length (m) (i.e., under the canopy)
real (kind=kind_phys), parameter :: z0 = 0.015 ! bare-soil roughness length (m) (i.e., under the canopy)

! ---------------------------------------------------------------------------------------------------
! initialize fluxes from veg. fraction
Expand Down Expand Up @@ -2629,10 +2629,10 @@ subroutine csnow (parameters,isnow ,nsnow ,nsoil ,snice ,snliq ,dzsnso
! thermal conductivity of snow

do iz = isnow+1, 0
! tksno(iz) = 3.2217e-6*bdsnoi(iz)**2. ! stieglitz(yen,1965)
! tksno(iz) = 3.2217e-6*bdsnoi(iz)**2. ! stieglitz(yen,1965)
! tksno(iz) = 2e-2+2.5e-6*bdsnoi(iz)*bdsnoi(iz) ! anderson, 1976
! tksno(iz) = 0.35 ! constant
tksno(iz) = 2.576e-6*bdsnoi(iz)**2. + 0.074 ! verseghy (1991)
tksno(iz) = 0.35 ! constant
! tksno(iz) = 2.576e-6*bdsnoi(iz)**2. + 0.074 ! verseghy (1991)
! tksno(iz) = 2.22*(bdsnoi(iz)/1000.)**1.88 ! douvill(yen, 1981)
enddo

Expand Down Expand Up @@ -4056,11 +4056,6 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , &

end if

! prepare for longwave rad.

air = -emv*(1.+(1.-emv)*(1.-emg))*lwdn - emv*emg*sb*tg**4
cir = (2.-emv*(1.-emg))*emv*sb
!
if(opt_sfc == 4) then

gdx = sqrt(garea1)
Expand Down Expand Up @@ -4207,6 +4202,11 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , &
end if
end if

! prepare for longwave rad.

air = -emv*(1.+(1.-emv)*(1.-emg))*lwdn - emv*emg*sb*tg**4
cir = (2.-emv*(1.-emg))*emv*sb

! prepare for sensible heat flux above veg.

cah = 1./rahc
Expand Down Expand Up @@ -4269,7 +4269,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , &

! update vegetation surface temperature
tv = tv + dtv
! tah = ata + bta*tv ! canopy air t; update here for consistency
tah = ata + bta*tv ! canopy air t; update here for consistency

! for computing m-o length in the next iteration
h = rhoair*cpair*(tah - sfctmp) /rahc
Expand All @@ -4282,15 +4282,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , &
qfx = (qsfc-qair)*rhoair*caw
endif


if (liter == 1) then
exit loop1
endif
if (iter >= 5 .and. abs(dtv) <= 0.01 .and. liter == 0) then
liter = 1
endif

end do loop1 ! end stability iteration
! after canopy balance, do the under-canopy ground balance

! under-canopy fluxes and tg

Expand All @@ -4300,8 +4292,6 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , &
cev = rhoair*cpair / (gammag*(rawg+rsurf)) ! barlage: change to ground v3.6
cgh = 2.*df(isnow+1)/dzsnso(isnow+1)

loop2: do iter = 1, niterg

t = tdc(tg)
call esat(t, esatw, esati, dsatw, dsati)
if (t .gt. 0.) then
Expand All @@ -4327,7 +4317,14 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , &
gh = gh + cgh*dtg
tg = tg + dtg

end do loop2
if (liter == 1) then
exit loop1
endif
if (iter >= 5 .and. abs(dtv) <= 0.01 .and. abs(dtg) <= 0.01 .and. liter == 0) then
liter = 1 ! if conditions are met, then do one final loop
endif

end do loop1

! tah = (cah*sfctmp + cvh*tv + cgh*tg)/(cah + cvh + cgh)

Expand Down Expand Up @@ -5824,7 +5821,8 @@ subroutine thermalz0(parameters, fveg, z0m, z0mg, zlvl,

if (opt_trs == z0heqz0m) then

z0m_out = exp(fveg * log(z0m) + (1.0 - fveg) * log(z0mg))
! z0m_out = exp(fveg * log(z0m) + (1.0 - fveg) * log(z0mg))
z0m_out = fveg * z0m + (1.0 - fveg) * z0mg
z0h_out = z0m_out

elseif (opt_trs == chen09) then
Expand All @@ -5841,7 +5839,7 @@ subroutine thermalz0(parameters, fveg, z0m, z0mg, zlvl,
endif

z0h_out = exp( fveg * log(z0m * exp(-czil*0.4*258.2*sqrt(ustarx*z0m))) + &
(1.0 - fveg) * log(max(z0m/exp(kb_sigma_f0),1.0e-6)) )
(1.0 - fveg) * log(max(z0mg/exp(kb_sigma_f0),1.0e-6)) )

elseif (opt_trs == tessel) then

Expand Down Expand Up @@ -5880,15 +5878,15 @@ subroutine thermalz0(parameters, fveg, z0m, z0mg, zlvl,

z0h_out = z0m_out

elseif (opt_trs == chen09 .or. opt_trs == tessel) then
elseif (opt_trs == tessel) then

if (vegtyp <= 5) then
z0h_out = z0m_out
else
z0h_out = z0m_out * 0.01
endif

elseif (opt_trs == blumel99) then
elseif (opt_trs == chen09 .or. opt_trs == blumel99) then

reyn = ustarx*z0m_out/viscosity ! Blumel99 eqn 36c
if (reyn > 2.0) then
Expand Down
11 changes: 7 additions & 4 deletions physics/SFC_Models/Land/Noahmp/noahmpdrv.F90
Original file line number Diff line number Diff line change
Expand Up @@ -663,6 +663,7 @@ subroutine noahmpdrv_run &
real (kind=kind_phys) :: precip_freeze_frac_in ! used for penman calculation

real (kind=kind_phys) :: virtfac1 ! virtual factor
real (kind=kind_phys) :: tflux ! surface flux temp
real (kind=kind_phys) :: tvs1 ! surface virtual temp
real (kind=kind_phys) :: vptemp ! virtual potential temp

Expand Down Expand Up @@ -944,7 +945,8 @@ subroutine noahmpdrv_run &
t2mmp(i) = temperature_bare_2m
q2mp(i) = spec_humidity_bare_2m

tskin(i) = temperature_ground
tskin(i) = temperature_radiative
tflux = temperature_ground
surface_temperature = temperature_ground
vegetation_fraction = vegetation_frac
ch_vegetated = 0.0
Expand Down Expand Up @@ -1038,7 +1040,8 @@ subroutine noahmpdrv_run &
q2mp(i) = spec_humidity_veg_2m * vegetation_fraction + &
spec_humidity_bare_2m * (1-vegetation_fraction)

tskin(i) = surface_temperature
tskin(i) = temperature_radiative
tflux = surface_temperature

endif ! glacial split ends

Expand Down Expand Up @@ -1194,9 +1197,9 @@ subroutine noahmpdrv_run &
endif

if(thsfc_loc) then ! Use local potential temperature
tvs1 = tskin(i) * virtfac1
tvs1 = tflux * virtfac1
else ! Use potential temperature referenced to 1000 hPa
tvs1 = tskin(i)/prsik1(i) * virtfac1
tvs1 = tflux/prsik1(i) * virtfac1
endif

z0_total = max(min(z0_total,forcing_height),1.0e-6)
Expand Down
Loading

0 comments on commit 9e736da

Please sign in to comment.