Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Changes necessary to remove GSD_SURFACE_FLUXES_BUGFIX from ccpp-physics #10

Merged
merged 2 commits into from
Apr 26, 2021
Merged

Changes necessary to remove GSD_SURFACE_FLUXES_BUGFIX from ccpp-physics #10

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
b51c5b3
Removing GSD_SURFACE_FLUXES_BUGFIX and replacing with flag thsfc_loc
Apr 26, 2021
b6d0ede
Bugfix
Apr 26, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 14 additions & 10 deletions physics/GFS_surface_composites.F90
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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