Skip to content

Commit

Permalink
Merge pull request NCAR#10 from benwgreen/thsfc_loc2
Browse files Browse the repository at this point in the history 
Changes necessary to remove GSD_SURFACE_FLUXES_BUGFIX from ccpp-physics
  • Loading branch information
@climbfuji
climbfuji authored Apr 26, 2021
2 parents 2ed6a7c + b6d0ede commit 986211f
Show file tree
Hide file tree
Showing 10 changed files with 164 additions and 72 deletions.
24 changes: 14 additions & 10 deletions physics/GFS_surface_composites.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -373,7 +373,7 @@ end subroutine GFS_surface_composites_post_finalize
!! \htmlinclude GFS_surface_composites_post_run.html
!!
subroutine GFS_surface_composites_post_run ( &
im, kice, km, cplflx, cplwav2atm, frac_grid, flag_cice, islmsk, dry, wet, icy, wind, t1, q1, prsl1, &
im, kice, km, cplflx, cplwav2atm, frac_grid, flag_cice, thsfc_loc, islmsk, dry, wet, icy, wind, t1, q1, prsl1, &
rd, rvrdm1, landfrac, lakefrac, oceanfrac, zorl, zorlo, zorll, zorli, &
cd, cd_wat, cd_lnd, cd_ice, cdq, cdq_wat, cdq_lnd, cdq_ice, rb, rb_wat, rb_lnd, rb_ice, stress, stress_wat, stress_lnd, &
stress_ice, ffmm, ffmm_wat, ffmm_lnd, ffmm_ice, ffhh, ffhh_wat, ffhh_lnd, ffhh_ice, uustar, uustar_wat, uustar_lnd, &
Expand Down Expand Up @@ -410,6 +410,7 @@ subroutine GFS_surface_composites_post_run (
real(kind=kind_phys), dimension(im, km), intent(inout) :: stc

! Additional data needed for calling "stability"
logical, intent(in ) :: thsfc_loc
real(kind=kind_phys), intent(in ) :: grav
real(kind=kind_phys), dimension(:), intent(in ) :: prslki, z1, ztmax_wat, ztmax_lnd, ztmax_ice

Expand All @@ -420,7 +421,7 @@ subroutine GFS_surface_composites_post_run (
integer :: i, k
real(kind=kind_phys) :: txl, txi, txo, wfrac, q0, rho
! For calling "stability"
real(kind=kind_phys) :: tsurf, virtfac, thv1, tvs, z0max, ztmax
real(kind=kind_phys) :: tsurf, virtfac, tv1, thv1, tvs, z0max, ztmax

! Initialize CCPP error handling variables
errmsg = ''
Expand Down Expand Up @@ -475,24 +476,27 @@ subroutine GFS_surface_composites_post_run (

q0 = max( q1(i), qmin )
virtfac = one + rvrdm1 * q0
#ifdef GSD_SURFACE_FLUXES_BUGFIX
thv1 = t1(i) / prslk1(i) * virtfac ! Theta-v at lowest level
tvs = half * (tsfc(i)+tsurf)/prsik1(i) * virtfac
#else
thv1 = t1(i) * prslki(i) * virtfac ! Theta-v at lowest level
tvs = half * (tsfc(i)+tsurf) * virtfac
#endif

tv1 = t1(i) * virtfac ! Virtual temperature in middle of lowest layer
if(thsfc_loc) then ! Use local potential temperature
thv1 = t1(i) * prslki(i) * virtfac ! Theta-v at lowest level
tvs = half * (tsfc(i)+tsurf) * virtfac
else ! Use potential temperature referenced to 1000 hPa
thv1 = t1(i) / prslk1(i) * virtfac ! Theta-v at lowest level
tvs = half * (tsfc(i)+tsurf)/prsik1(i) * virtfac
endif

zorl(i) = exp(txl*log(zorll(i)) + txi*log(zorli(i)) + txo*log(zorlo(i)))
z0max = 0.01_kind_phys * zorl(i)
ztmax = exp(txl*log(ztmax_lnd(i)) + txi*log(ztmax_ice(i)) + txo*log(ztmax_wat(i)))

call stability(z1(i), snowd(i), thv1, wind(i), z0max, ztmax, tvs, grav, & ! inputs
tv1, thsfc_loc, & ! inputs
rb(i), ffmm(i), ffhh(i), fm10(i), fh2(i), cd(i), cdq(i), & ! outputs
stress(i), uustar(i))

! BWG, 2021/02/25: cmm=cd*wind, chh=cdq*wind, so use composite cd, cdq
rho = prsl1(i) / (rd*t1(i)*(one + rvrdm1*q0))
rho = prsl1(i) / (rd*t1(i)*virtfac)
cmm(i) = cd(i)*wind(i) !txl*cmm_lnd(i) + txi*cmm_ice(i) + txo*cmm_wat(i)
chh(i) = rho*cdq(i)*wind(i) !txl*chh_lnd(i) + txi*chh_ice(i) + txo*chh_wat(i)

Expand Down
8 changes: 8 additions & 0 deletions physics/GFS_surface_composites.meta
View file
Original file line number Diff line number Diff line change
Expand Up @@ -908,6 +908,14 @@
type = logical
intent = in
optional = F
[thsfc_loc]
standard_name = flag_for_reference_pressure_theta
long_name = flag for reference pressure in theta calculation
units = flag
dimensions = ()
type = logical
intent = in
optional = F
[islmsk]
standard_name = sea_land_ice_mask
long_name = sea/land/ice mask (=0/1/2)
Expand Down
14 changes: 8 additions & 6 deletions physics/sfc_diag.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ end subroutine sfc_diag_finalize
!! @{
subroutine sfc_diag_run &
& (im,grav,cp,eps,epsm1,ps,u1,v1,t1,q1,prslki, &
& evap,fm,fh,fm10,fh2,tskin,qsurf, &
& evap,fm,fh,fm10,fh2,tskin,qsurf,thsfc_loc, &
& f10m,u10m,v10m,t2m,q2m,errmsg,errflg &
& )
!
Expand All @@ -32,6 +32,7 @@ subroutine sfc_diag_run &
implicit none
!
integer, intent(in) :: im
logical, intent(in) :: thsfc_loc ! Flag for reference pot. temp.
real(kind=kind_phys), intent(in) :: grav,cp,eps,epsm1
real(kind=kind_phys), dimension(im), intent(in) :: &
& ps, u1, v1, t1, q1, tskin, &
Expand Down Expand Up @@ -74,11 +75,12 @@ subroutine sfc_diag_run &
! t2m(i) = t2m(i) * sig2k
wrk = 1.0 - fhi

#ifdef GSD_SURFACE_FLUXES_BUGFIX
t2m(i) = tskin(i)*wrk + t1(i)*fhi - (grav+grav)/cp
#else
t2m(i) = tskin(i)*wrk + t1(i)*prslki(i)*fhi - (grav+grav)/cp
#endif

if(thsfc_loc) then ! Use local potential temperature
t2m(i) = tskin(i)*wrk + t1(i)*prslki(i)*fhi - (grav+grav)/cp
else ! Use potential temperature referenced to 1000 hPa
t2m(i) = tskin(i)*wrk + t1(i)*fhi - (grav+grav)/cp
endif

if(evap(i) >= 0.) then ! for evaporation>0, use inferred qsurf to deduce q2m
q2m(i) = qsurf(i)*wrk + max(qmin,q1(i))*fhi
Expand Down
8 changes: 8 additions & 0 deletions physics/sfc_diag.meta
View file
Original file line number Diff line number Diff line change
Expand Up @@ -168,6 +168,14 @@
kind = kind_phys
intent = in
optional = F
[thsfc_loc]
standard_name = flag_for_reference_pressure_theta
long_name = flag for reference pressure in theta calculation
units = flag
dimensions = ()
type = logical
intent = in
optional = F
[f10m]
standard_name = ratio_of_wind_at_lowest_model_layer_and_wind_at_10m
long_name = ratio of fm10 and fm
Expand Down
71 changes: 52 additions & 19 deletions physics/sfc_diff.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,6 +69,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in)
& flag_iter,redrag, & !intent(in)
& u10m,v10m,sfc_z0_type, & !hafs,z0 type !intent(in)
& wet,dry,icy, & !intent(in)
& thsfc_loc, & !intent(in)
& tskin_wat, tskin_lnd, tskin_ice, & !intent(in)
& tsurf_wat, tsurf_lnd, tsurf_ice, & !intent(in)
& snwdph_wat,snwdph_lnd,snwdph_ice, & !intent(in)
Expand Down Expand Up @@ -97,6 +98,8 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in)
logical, intent(in) :: redrag ! reduced drag coeff. flag for high wind over sea (j.han)
logical, dimension(im), intent(in) :: flag_iter, wet, dry, icy
logical, intent(in) :: thsfc_loc ! Flag for reference pressure in theta calculation
real(kind=kind_phys), dimension(im), intent(in) :: u10m,v10m
real(kind=kind_phys), intent(in) :: rvrdm1, eps, epsm1, grav
real(kind=kind_phys), dimension(im), intent(in) :: &
Expand Down Expand Up @@ -133,6 +136,9 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in)
real(kind=kind_phys) :: rat, thv1, restar, wind10m,
& czilc, tem1, tem2, virtfac
!
real(kind=kind_phys) :: tv1
real(kind=kind_phys) :: tvs, z0, z0max
!
real(kind=kind_phys), parameter ::
Expand Down Expand Up @@ -178,18 +184,26 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in)
ztmax_wat(i) = 1. ! log(1) = 0
virtfac = one + rvrdm1 * max(q1(i),qmin)
thv1 = t1(i) * prslki(i) * virtfac
tv1 = t1(i) * virtfac ! Virtual temperature in middle of lowest layer
if(thsfc_loc) then ! Use local potential temperature
thv1 = t1(i) * prslki(i) * virtfac
else ! Use potential temperature reference to 1000 hPa
thv1 = t1(i) / prslk1(i) * virtfac
endif
! compute stability dependent exchange coefficients
! this portion of the code is presently suppressed
!
if (dry(i)) then ! Some land
#ifdef GSD_SURFACE_FLUXES_BUGFIX
tvs = half * (tsurf_lnd(i)+tskin_lnd(i))/prsik1(i)
& * virtfac
#else
tvs = half * (tsurf_lnd(i)+tskin_lnd(i)) * virtfac
#endif
if(thsfc_loc) then ! Use local potential temperature
tvs = half * (tsurf_lnd(i)+tskin_lnd(i)) * virtfac
else ! Use potential temperature referenced to 1000 hPa
tvs = half * (tsurf_lnd(i)+tskin_lnd(i))/prsik1(i)
& * virtfac
endif
z0max = max(zmin, min(0.01_kp * z0rl_lnd(i), z1(i)))
!** xubin's new z0 over land
tem1 = one - shdmax(i)
Expand Down Expand Up @@ -253,14 +267,21 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in)
call stability
! --- inputs:
& (z1(i), snwdph_lnd(i), thv1, wind(i),
& z0max, ztmax_lnd(i), tvs, grav,
& z0max, ztmax_lnd(i), tvs, grav, tv1, thsfc_loc,
! --- outputs:
& rb_lnd(i), fm_lnd(i), fh_lnd(i), fm10_lnd(i), fh2_lnd(i),
& cm_lnd(i), ch_lnd(i), stress_lnd(i), ustar_lnd(i))
endif ! Dry points
if (icy(i)) then ! Some ice
tvs = half * (tsurf_ice(i)+tskin_ice(i)) * virtfac
if(thsfc_loc) then ! Use local potential temperature
tvs = half * (tsurf_ice(i)+tskin_ice(i)) * virtfac
else ! Use potential temperature referenced to 1000 hPa
tvs = half * (tsurf_ice(i)+tskin_ice(i))/prsik1(i)
& * virtfac
endif
z0max = max(zmin, min(0.01_kp * z0rl_ice(i), z1(i)))
!** xubin's new z0 over land and sea ice
tem1 = one - shdmax(i)
Expand Down Expand Up @@ -288,7 +309,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in)
call stability
! --- inputs:
& (z1(i), snwdph_ice(i), thv1, wind(i),
& z0max, ztmax_ice(i), tvs, grav,
& z0max, ztmax_ice(i), tvs, grav, tv1, thsfc_loc,
! --- outputs:
& rb_ice(i), fm_ice(i), fh_ice(i), fm10_ice(i), fh2_ice(i),
& cm_ice(i), ch_ice(i), stress_ice(i), ustar_ice(i))
Expand All @@ -298,7 +319,14 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in)
! the stuff now put into "stability"
if (wet(i)) then ! Some open ocean
tvs = half * (tsurf_wat(i)+tskin_wat(i)) * virtfac
if(thsfc_loc) then ! Use local potential temperature
tvs = half * (tsurf_wat(i)+tskin_wat(i)) * virtfac
else
tvs = half * (tsurf_wat(i)+tskin_wat(i))/prsik1(i)
& * virtfac
endif
z0 = 0.01_kp * z0rl_wat(i)
z0max = max(zmin, min(z0,z1(i)))
ustar_wat(i) = sqrt(grav * z0 / charnock)
Expand Down Expand Up @@ -332,7 +360,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in)
call stability
! --- inputs:
& (z1(i), snwdph_wat(i), thv1, wind(i),
& z0max, ztmax_wat(i), tvs, grav,
& z0max, ztmax_wat(i), tvs, grav, tv1, thsfc_loc,
! --- outputs:
& rb_wat(i), fm_wat(i), fh_wat(i), fm10_wat(i), fh2_wat(i),
& cm_wat(i), ch_wat(i), stress_wat(i), ustar_wat(i))
Expand Down Expand Up @@ -392,6 +420,7 @@ end subroutine sfc_diff_run
subroutine stability &
! --- inputs:
& ( z1, snwdph, thv1, wind, z0max, ztmax, tvs, grav, &
& tv1, thsfc_loc, &
! --- outputs:
& rb, fm, fh, fm10, fh2, cm, ch, stress, ustar)
!-----
Expand All @@ -400,6 +429,8 @@ subroutine stability &
! --- inputs:
real(kind=kind_phys), intent(in) :: &
& z1, snwdph, thv1, wind, z0max, ztmax, tvs, grav
real(kind=kind_phys), intent(in) :: tv1
logical, intent(in) :: thsfc_loc

! --- outputs:
real(kind=kind_phys), intent(out) :: &
Expand Down Expand Up @@ -435,13 +466,15 @@ subroutine stability &
dtv = thv1 - tvs
adtv = max(abs(dtv),0.001_kp)
dtv = sign(1.,dtv) * adtv
#ifdef GSD_SURFACE_FLUXES_BUGFIX
rb = max(-5000.0_kp, grav * dtv * z1
& / (thv1 * wind * wind))
#else
rb = max(-5000.0_kp, (grav+grav) * dtv * z1
& / ((thv1 + tvs) * wind * wind))
#endif

if(thsfc_loc) then ! Use local potential temperature
rb = max(-5000.0_kp, (grav+grav) * dtv * z1
& / ((thv1 + tvs) * wind * wind))
else ! Use potential temperature referenced to 1000 hPa
rb = max(-5000.0_kp, grav * dtv * z1
& / (tv1 * wind * wind))
endif

tem1 = one / z0max
tem2 = one / ztmax
fm = log((z0max+z1) * tem1)
Expand Down
8 changes: 8 additions & 0 deletions physics/sfc_diff.meta
View file
Original file line number Diff line number Diff line change
Expand Up @@ -250,6 +250,14 @@
type = logical
intent = in
optional = F
[thsfc_loc]
standard_name = flag_for_reference_pressure_theta
long_name = flag for reference pressure in theta calculation
units = flag
dimensions = ()
type = logical
intent = in
optional = F
[tskin_wat]
standard_name = surface_skin_temperature_over_water_interstitial
long_name = surface skin temperature over water (temporary use as interstitial)
Expand Down
41 changes: 24 additions & 17 deletions physics/sfc_nst.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ subroutine sfc_nst_run &
& sinlat, stress, &
& sfcemis, dlwflx, sfcnsw, rain, timestep, kdt, solhr,xcosz, &
& wind, flag_iter, flag_guess, nstf_name1, nstf_name4, &
& nstf_name5, lprnt, ipr, &
& nstf_name5, lprnt, ipr, thsfc_loc, &
& tskin, tsurf, xt, xs, xu, xv, xz, zm, xtts, xzts, dt_cool, & ! --- input/output:
& z_c, c_0, c_d, w_0, w_d, d_conv, ifd, qrain, &
& qsurf, gflux, cmm, chh, evap, hflx, ep, errmsg, errflg & ! --- outputs:
Expand All @@ -50,7 +50,7 @@ subroutine sfc_nst_run &
! prsl1, prslki, wet, use_flake, xlon, sinlat, stress, !
! sfcemis, dlwflx, sfcnsw, rain, timestep, kdt,solhr,xcosz, !
! wind, flag_iter, flag_guess, nstf_name1, nstf_name4, !
! nstf_name5, lprnt, ipr, !
! nstf_name5, lprnt, ipr, thsfc_loc, !
! input/outputs: !
! tskin, tsurf, xt, xs, xu, xv, xz, zm, xtts, xzts, dt_cool, !
! z_c, c_0, c_d, w_0, w_d, d_conv, ifd, qrain, !
Expand Down Expand Up @@ -123,6 +123,7 @@ subroutine sfc_nst_run &
! nstf_name5 : zsea2 in mm 1 !
! lprnt - logical, control flag for check print out 1 !
! ipr - integer, grid index for check print out 1 !
! thsfc_loc- logical, flag for reference pressure in theta 1 !
! !
! input/outputs:
! li added for oceanic components
Expand Down Expand Up @@ -199,6 +200,7 @@ subroutine sfc_nst_run &
& use_flake
! &, icy
logical, intent(in) :: lprnt
logical, intent(in) :: thsfc_loc

! --- input/outputs:
! control variables of dtl system (5+2) and sl (2) and coefficients for d(tz)/d(ts) calculation
Expand Down Expand Up @@ -297,11 +299,13 @@ subroutine sfc_nst_run &
wndmag(i) = sqrt(u1(i)*u1(i) + v1(i)*v1(i))

q0(i) = max(q1(i), 1.0e-8_kp)
#ifdef GSD_SURFACE_FLUXES_BUGFIX
theta1(i) = t1(i) / prslk1(i) ! potential temperature at the middle of lowest model layer
#else
theta1(i) = t1(i) * prslki(i)
#endif

if(thsfc_loc) then ! Use local potential temperature
theta1(i) = t1(i) * prslki(i)
else ! Use potential temperature referenced to 1000 hPa
theta1(i) = t1(i) / prslk1(i) ! potential temperature at the middle of lowest model layer
endif

tv1(i) = t1(i) * (one + rvrdm1*q0(i))
rho_a(i) = prsl1(i) / (rd*tv1(i))
qss(i) = fpvs(tsurf(i)) ! pa
Expand All @@ -322,11 +326,12 @@ subroutine sfc_nst_run &
! at previous time step
evap(i) = elocp * rch(i) * (qss(i) - q0(i))
qsurf(i) = qss(i)
#ifdef GSD_SURFACE_FLUXES_BUGFIX
hflx(i) = rch(i) * (tsurf(i)/prsik1(i) - theta1(i))
#else
hflx(i) = rch(i) * (tsurf(i) - theta1(i))
#endif

if(thsfc_loc) then ! Use local potential temperature
hflx(i) = rch(i) * (tsurf(i) - theta1(i))
else ! Use potential temperature referenced to 1000 hPa
hflx(i) = rch(i) * (tsurf(i)/prsik1(i) - theta1(i))
endif

! if (lprnt .and. i == ipr) print *,' tskin=',tskin(i),' theta1=',
! & theta1(i),' hflx=',hflx(i),' t1=',t1(i),'prslki=',prslki(i)
Expand Down Expand Up @@ -621,11 +626,13 @@ subroutine sfc_nst_run &
qss(i) = eps*qss(i) / (ps(i) + epsm1*qss(i))
qsurf(i) = qss(i)
evap(i) = elocp*rch(i) * (qss(i) - q0(i))
#ifdef GSD_SURFACE_FLUXES_BUGFIX
hflx(i) = rch(i) * (tskin(i)/prsik1(i) - theta1(i))
#else
hflx(i) = rch(i) * (tskin(i) - theta1(i))
#endif
if(thsfc_loc) then ! Use local potential temperature
hflx(i) = rch(i) * (tskin(i) - theta1(i))
else ! Use potential temperature referenced to 1000 hPa
hflx(i) = rch(i) * (tskin(i)/prsik1(i) - theta1(i))
endif
endif
enddo
endif ! if ( nstf_name1 > 1 ) then
Expand Down
Loading

0 comments on commit 986211f

Please sign in to comment.