Skip to content
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
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
34 changes: 20 additions & 14 deletions physics/GFS_surface_composites.F90
Original file line number Diff line number Diff line change
Expand Up @@ -351,7 +351,7 @@ end subroutine GFS_surface_composites_post_finalize
!!
#endif
subroutine GFS_surface_composites_post_run ( &
im, kice, km, cplflx, cplwav2atm, frac_grid, flag_cice, islmsk, dry, wet, icy, landfrac, lakefrac, oceanfrac, &
im, kice, km, cplflx, cplwav2atm, frac_grid, flag_cice, islmsk, dry, wet, icy, wind, landfrac, lakefrac, oceanfrac, &
zorl, zorlo, zorll, zorli, zorl_wat, zorl_lnd, zorl_ice, &
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 All @@ -367,7 +367,7 @@ subroutine GFS_surface_composites_post_run (
logical, intent(in) :: cplflx, frac_grid, cplwav2atm
logical, dimension(im), intent(in) :: flag_cice, dry, wet, icy
integer, dimension(im), intent(in) :: islmsk
real(kind=kind_phys), dimension(im), intent(in) :: landfrac, lakefrac, oceanfrac, &
real(kind=kind_phys), dimension(im), intent(in) :: wind, landfrac, lakefrac, oceanfrac, &
zorl_wat, zorl_lnd, zorl_ice, cd_wat, cd_lnd, cd_ice, cdq_wat, cdq_lnd, cdq_ice, rb_wat, rb_lnd, rb_ice, stress_wat, &
stress_lnd, stress_ice, ffmm_wat, ffmm_lnd, ffmm_ice, ffhh_wat, ffhh_lnd, ffhh_ice, uustar_wat, uustar_lnd, uustar_ice, &
fm10_wat, fm10_lnd, fm10_ice, fh2_wat, fh2_lnd, fh2_ice, tsurf_wat, tsurf_lnd, tsurf_ice, cmm_wat, cmm_lnd, cmm_ice, &
Expand Down Expand Up @@ -408,20 +408,24 @@ subroutine GFS_surface_composites_post_run (
txi = cice(i) * wfrac ! txi = ice fraction wrt whole cell
txo = max(zero, wfrac-txi) ! txo = open water fraction

zorl(i) = txl*zorl_lnd(i) + txi*zorl_ice(i) + txo*zorl_wat(i)
cd(i) = txl*cd_lnd(i) + txi*cd_ice(i) + txo*cd_wat(i)
cdq(i) = txl*cdq_lnd(i) + txi*cdq_ice(i) + txo*cdq_wat(i)
rb(i) = txl*rb_lnd(i) + txi*rb_ice(i) + txo*rb_wat(i)
stress(i) = txl*stress_lnd(i) + txi*stress_ice(i) + txo*stress_wat(i)
ffmm(i) = txl*ffmm_lnd(i) + txi*ffmm_ice(i) + txo*ffmm_wat(i)
ffhh(i) = txl*ffhh_lnd(i) + txi*ffhh_ice(i) + txo*ffhh_wat(i)
uustar(i) = txl*uustar_lnd(i) + txi*uustar_ice(i) + txo*uustar_wat(i)
fm10(i) = txl*fm10_lnd(i) + txi*fm10_ice(i) + txo*fm10_wat(i)
fh2(i) = txl*fh2_lnd(i) + txi*fh2_ice(i) + txo*fh2_wat(i)
! BWG zorl(i) = txl*zorl_lnd(i) + txi*zorl_ice(i) + txo*zorl_wat(i)
! BWG cd(i) = txl*cd_lnd(i) + txi*cd_ice(i) + txo*cd_wat(i)
! BWG cdq(i) = txl*cdq_lnd(i) + txi*cdq_ice(i) + txo*cdq_wat(i)
! BWG rb(i) = txl*rb_lnd(i) + txi*rb_ice(i) + txo*rb_wat(i)
! BWG stress(i) = txl*stress_lnd(i) + txi*stress_ice(i) + txo*stress_wat(i)
! BWG ffmm(i) = txl*ffmm_lnd(i) + txi*ffmm_ice(i) + txo*ffmm_wat(i)
! BWG ffhh(i) = txl*ffhh_lnd(i) + txi*ffhh_ice(i) + txo*ffhh_wat(i)
! BWG uustar(i) = txl*uustar_lnd(i) + txi*uustar_ice(i) + txo*uustar_wat(i)
! BWG fm10(i) = txl*fm10_lnd(i) + txi*fm10_ice(i) + txo*fm10_wat(i)
! BWG fh2(i) = txl*fh2_lnd(i) + txi*fh2_ice(i) + txo*fh2_wat(i)

!tsurf(i) = txl*tsurf_lnd(i) + txi*tice(i) + txo*tsurf_wat(i)
!tsurf(i) = txl*tsurf_lnd(i) + txi*tsurf_ice(i) + txo*tsurf_wat(i) ! not used again! Moorthi
cmm(i) = txl*cmm_lnd(i) + txi*cmm_ice(i) + txo*cmm_wat(i)
chh(i) = txl*chh_lnd(i) + txi*chh_ice(i) + txo*chh_wat(i)

! BWG, 2021/02/25: cmm=cd*wind, chh=cdq*wind, so use composite cd, cdq
Copy link
Copy Markdown
Collaborator Author

@benwgreen benwgreen Mar 1, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Bug: chh = rho * cdq * wind (cmm = cd*wind is still correct)

cmm(i) = cd(i)*wind(i) !txl*cmm_lnd(i) + txi*cmm_ice(i) + txo*cmm_wat(i)
chh(i) = cdq(i)*wind(i) !txl*chh_lnd(i) + txi*chh_ice(i) + txo*chh_wat(i)

!gflx(i) = txl*gflx_lnd(i) + txi*gflx_ice(i) + txo*gflx_wat(i)
ep1d(i) = txl*ep1d_lnd(i) + txi*ep1d_ice(i) + txo*ep1d_wat(i)
!weasd(i) = txl*weasd_lnd(i) + txi*weasd_ice(i) + txo*weasd_wat(i)
Expand All @@ -441,6 +445,8 @@ subroutine GFS_surface_composites_post_run (
qss(i) = txl*qss_lnd(i) + txi*qss_ice(i) + txo*qss_wat(i)
gflx(i) = txl*gflx_lnd(i) + txi*gflx_ice(i) + txo*gflx_wat(i)
endif

! BWG, 2021/02/25: Need to change composite skin temperature base on ULW (Fanglin)
tsfc(i) = txl*tsfc_lnd(i) + txi*tice(i) + txo*tsfc_wat(i)
Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@barlage would we need to change composite formula here also, or just simply use your composite straight from sfc_diff?

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@benwgreen for now, using the composite here should be sufficient; I have to admit I'm ignorant on what tsfc is used for after this point, i.e., PBL, radiation, both? @yangfanglin might be able to comment; I will state here that this is where I think we need two surface temperatures, one for fluxes and one for radiation. This does depend on what the PBL/radiation use for a boundary condition (flux vs. temperature) but to be as general as possible, we probably should have two temperatures.

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@barlage After thinking more about @climbfuji's comments on what to do if sfc_diff isn't called, for the immediate term it is probably best to leave this line as-is. Once we come up with a way to better generalize the solution, it could very well be that a value from sfc_diff (via meta) ends up here. We will keep you posted as we work through the problem.


zorll(i) = zorl_lnd(i)
Expand Down
9 changes: 9 additions & 0 deletions physics/GFS_surface_composites.meta
Original file line number Diff line number Diff line change
Expand Up @@ -917,6 +917,15 @@
type = logical
intent = in
optional = F
[wind]
standard_name = wind_speed_at_lowest_model_layer
long_name = wind speed at lowest model level
units = m s-1
dimensions = (horizontal_loop_extent)
type = real
kind = kind_phys
intent = in
optional = F
[landfrac]
standard_name = land_area_fraction
long_name = fraction of horizontal grid area occupied by land
Expand Down
176 changes: 144 additions & 32 deletions physics/sfc_diff.f
Original file line number Diff line number Diff line change
Expand Up @@ -71,17 +71,19 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !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)
& landfrac, cice, & !intent(in) -- for use with frac_grid
& islmsk, frac_grid, & !intent(in) -- for use with frac_grid
& z0rl_wat, z0rl_lnd, z0rl_ice, & !intent(inout)
& z0rl_wav, & !intent(inout)
& ustar_wat, ustar_lnd, ustar_ice, & !intent(inout)
& cm_wat, cm_lnd, cm_ice, & !intent(inout)
& ch_wat, ch_lnd, ch_ice, & !intent(inout)
& rb_wat, rb_lnd, rb_ice, & !intent(inout)
& stress_wat,stress_lnd,stress_ice, & !intent(inout)
& fm_wat, fm_lnd, fm_ice, & !intent(inout)
& fh_wat, fh_lnd, fh_ice, & !intent(inout)
& fm10_wat, fm10_lnd, fm10_ice, & !intent(inout)
& fh2_wat, fh2_lnd, fh2_ice, & !intent(inout)
& z0rl_wav, z0rl_cmp, & !intent(inout)
& ustar_wat, ustar_lnd, ustar_ice, ustar_cmp, & !intent(inout)
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Instead of passing all the composite variables to sfc_diff_run, call stabilty from GFS_surface_composites_post_run for the composite variables?

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We just need to think about what to do if sfc_diff is not used, i.e. for example for MYNNSFC. Maybe call the same stability function for the moment and then work with @joeolson42 to call the appropriate function from MYNNSFC in GFS_surface_composites_post_run?

& cm_wat, cm_lnd, cm_ice, cm_cmp, & !intent(inout)
& ch_wat, ch_lnd, ch_ice, ch_cmp, & !intent(inout)
& rb_wat, rb_lnd, rb_ice, rb_cmp, & !intent(inout)
& stress_wat,stress_lnd,stress_ice,stress_cmp, & !intent(inout)
& fm_wat, fm_lnd, fm_ice, fm_cmp, & !intent(inout)
& fh_wat, fh_lnd, fh_ice, fh_cmp, & !intent(inout)
& fm10_wat, fm10_lnd, fm10_ice, fm10_cmp, & !intent(inout)
& fh2_wat, fh2_lnd, fh2_ice, fh2_cmp, & !intent(inout)
& errmsg, errflg) !intent(out)
!
implicit none
Expand All @@ -107,17 +109,25 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in)
& snwdph_wat,snwdph_lnd,snwdph_ice

real(kind=kind_phys), dimension(im), intent(in) :: z0rl_wav

real(kind=kind_phys), dimension(im), intent(in) :: &
& landfrac, cice

integer, dimension(im), intent(in) :: islmsk ! For compositing

logical, intent(in) :: frac_grid ! For compositing

real(kind=kind_phys), dimension(im), intent(inout) :: &
& z0rl_wat, z0rl_lnd, z0rl_ice, &
& ustar_wat, ustar_lnd, ustar_ice, &
& cm_wat, cm_lnd, cm_ice, &
& ch_wat, ch_lnd, ch_ice, &
& rb_wat, rb_lnd, rb_ice, &
& stress_wat,stress_lnd,stress_ice, &
& fm_wat, fm_lnd, fm_ice, &
& fh_wat, fh_lnd, fh_ice, &
& fm10_wat, fm10_lnd, fm10_ice, &
& fh2_wat, fh2_lnd, fh2_ice
& z0rl_wat, z0rl_lnd, z0rl_ice, z0rl_cmp, &
& ustar_wat, ustar_lnd, ustar_ice, ustar_cmp, &
& cm_wat, cm_lnd, cm_ice, cm_cmp, &
& ch_wat, ch_lnd, ch_ice, ch_cmp, &
& rb_wat, rb_lnd, rb_ice, rb_cmp, &
& stress_wat,stress_lnd,stress_ice,stress_cmp, &
& fm_wat, fm_lnd, fm_ice, fm_cmp, &
& fh_wat, fh_lnd, fh_ice, fh_cmp, &
& fm10_wat, fm10_lnd, fm10_ice, fm10_cmp, &
& fh2_wat, fh2_lnd, fh2_ice, fh2_cmp
character(len=*), intent(out) :: errmsg
integer, intent(out) :: errflg
!
Expand All @@ -128,7 +138,14 @@ 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) :: tvs, z0, z0max, ztmax
real(kind=kind_phys) :: tvs, z0, z0max

real(kind=kind_phys), dimension(im) :: &
& ztmax_wat, ztmax_lnd, ztmax_ice

real(kind=kind_phys) :: txl, txi, txo, wfrac ! For fractional
real(kind=kind_phys) :: snwdph_cmp, ztmax_cmp! For fractional
real(kind=kind_phys) :: tskin_cmp, tsurf_cmp ! For fractional
!
real(kind=kind_phys), parameter ::
& one=1.0_kp, zero=0.0_kp, half=0.5_kp, qmin=1.0e-8_kp
Expand Down Expand Up @@ -166,6 +183,12 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in)

do i=1,im
if(flag_iter(i)) then

! BWG: Need to initialize ztmax arrays
ztmax_lnd(i) = 1. ! log(1) = 0
ztmax_ice(i) = 1. ! log(1) = 0
ztmax_wat(i) = 1. ! log(1) = 0

virtfac = one + rvrdm1 * max(q1(i),qmin)
thv1 = t1(i) * prslki(i) * virtfac

Expand Down Expand Up @@ -229,20 +252,20 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in)
czilc = 0.8_kp

tem1 = 1.0_kp - sigmaf(i)
ztmax = z0max*exp( - tem1*tem1
ztmax_lnd(i) = z0max*exp( - tem1*tem1
& * czilc*ca*sqrt(ustar_lnd(i)*(0.01/1.5e-05)))


! mg, sfc-perts: add surface perturbations to ztmax/z0max ratio over land
if (ztpert(i) /= zero) then
ztmax = ztmax * (10.0_kp**ztpert(i))
ztmax_lnd(i) = ztmax_lnd(i) * (10.0_kp**ztpert(i))
endif
ztmax = max(ztmax, zmin)
ztmax_lnd(i) = max(ztmax_lnd(i), zmin)
!
call stability
! --- inputs:
& (z1(i), snwdph_lnd(i), thv1, wind(i),
& z0max, ztmax, tvs, grav,
& z0max, ztmax_lnd(i), tvs, grav,
! --- 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))
Expand Down Expand Up @@ -270,14 +293,14 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in)
czilc = 0.8_kp

tem1 = 1.0_kp - sigmaf(i)
ztmax = z0max*exp( - tem1*tem1
ztmax_ice(i) = z0max*exp( - tem1*tem1
& * czilc*ca*sqrt(ustar_ice(i)*(0.01/1.5e-05)))
ztmax = max(ztmax, 1.0e-6)
ztmax_ice(i) = max(ztmax_ice(i), 1.0e-6)
!
call stability
! --- inputs:
& (z1(i), snwdph_ice(i), thv1, wind(i),
& z0max, ztmax, tvs, grav,
& z0max, ztmax_ice(i), tvs, grav,
! --- 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 Down Expand Up @@ -307,12 +330,12 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in)
! rat taken from zeng, zhao and dickinson 1997

rat = min(7.0_kp, 2.67_kp * sqrt(sqrt(restar)) - 2.57_kp)
ztmax = max(z0max * exp(-rat), zmin)
ztmax_wat(i) = max(z0max * exp(-rat), zmin)
!
if (sfc_z0_type == 6) then
call znot_t_v6(wind10m, ztmax) ! 10-m wind,m/s, ztmax(m)
call znot_t_v6(wind10m, ztmax_wat(i)) ! 10-m wind,m/s, ztmax(m)
else if (sfc_z0_type == 7) then
call znot_t_v7(wind10m, ztmax) ! 10-m wind,m/s, ztmax(m)
call znot_t_v7(wind10m, ztmax_wat(i)) ! 10-m wind,m/s, ztmax(m)
else if (sfc_z0_type > 0) then
write(0,*)'no option for sfc_z0_type=',sfc_z0_type
stop
Expand All @@ -321,7 +344,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, tvs, grav,
& z0max, ztmax_wat(i), tvs, grav,
! --- 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 @@ -372,6 +395,95 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in)
endif ! end of if(flagiter) loop
enddo

! BWG, 2021/02/23: For fractional grid, get composite values
if (frac_grid) then ! If fractional grid is on...
do i=1,im ! Loop over horizontal
if(flag_iter(i)) then
virtfac = one + rvrdm1 * max(q1(i),qmin)
#ifdef GSD_SURFACE_FLUXES_BUGFIX
thv1 = t1(i) / prslk1(i) * virtfac ! Theta-v at lowest level
#else
thv1 = t1(i) * prslki(i) * virtfac ! Theta-v at lowest level
#endif

! Three-way composites (fields from sfc_diff)
txl = landfrac(i) ! land fraction
wfrac = one - txl ! ocean fraction
txi = cice(i) * wfrac ! txi = ice fraction wrt whole cell
txo = max(zero, wfrac-txi) ! txo = open water fraction

! Composite inputs to "stability" function
snwdph_cmp = txl*snwdph_lnd(i) + txi*snwdph_ice(i)
tsurf_cmp = txl* tsurf_lnd(i) + txi* tsurf_ice(i) &
& + txo* tsurf_wat(i)
tskin_cmp = txl* tskin_lnd(i) + txi* tskin_ice(i) &
& + txo* tskin_wat(i)
#ifdef GSD_SURFACE_FLUXES_BUGFIX
tvs = half * (tsurf_cmp+tskin_cmp)/prsik1(i)
& * virtfac
#else
tvs = half * (tsurf_cmp+tskin_cmp) * virtfac
#endif
z0rl_cmp(i) = txl*log(z0rl_lnd(i)) + txi*log(z0rl_ice(i)) &
& + txo*log(z0rl_wat(i))
z0rl_cmp(i) = exp(z0rl_cmp(i))
z0max = 0.01_kp * z0rl_cmp(i)

ztmax_cmp = txl*log(ztmax_lnd(i))+txi*log(ztmax_ice(i)) &
& + txo*log(ztmax_wat(i))
ztmax_cmp = exp(ztmax_cmp)
!
call stability
! --- inputs:
& (z1(i), snwdph_cmp, thv1, wind(i),
& z0max, ztmax_cmp, tvs, grav,
! --- outputs:
& rb_cmp(i), fm_cmp(i), fh_cmp(i), fm10_cmp(i), fh2_cmp(i),
& cm_cmp(i), ch_cmp(i), stress_cmp(i), ustar_cmp(i))

endif ! end of if(flagiter) loop
enddo ! End of loop over horizontal
else ! If frac_grid is false
do i=1,im ! Loop over horizontal
if(flag_iter(i)) then
if (islmsk(i) == 1) then ! Land
z0rl_cmp(i) = z0rl_lnd(i)
ustar_cmp(i) = ustar_lnd(i)
cm_cmp(i) = cm_lnd(i)
ch_cmp(i) = ch_lnd(i)
rb_cmp(i) = rb_lnd(i)
stress_cmp(i) = stress_lnd(i)
fm_cmp(i) = fm_lnd(i)
fh_cmp(i) = fh_lnd(i)
fm10_cmp(i) = fm10_lnd(i)
fh2_cmp(i) = fh2_lnd(i)
elseif (islmsk(i) == 0) then ! Open water
z0rl_cmp(i) = z0rl_wat(i)
ustar_cmp(i) = ustar_wat(i)
cm_cmp(i) = cm_wat(i)
ch_cmp(i) = ch_wat(i)
rb_cmp(i) = rb_wat(i)
stress_cmp(i) = stress_wat(i)
fm_cmp(i) = fm_wat(i)
fh_cmp(i) = fh_wat(i)
fm10_cmp(i) = fm10_wat(i)
fh2_cmp(i) = fh2_wat(i)
else ! if (islmsk(i) == 2) ! Ice
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is incomplete. There needs to be composite between water and ice even in non-fractional grid case.

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@SMoorthi-emc currently the only composited variables in the non-frac case (Lines 579-601 of original GFS_surface_composites.F90) are: zorl, tsfc, evap, hflx, stress, qss, ep1d. This portion of the proposed code change probably won't survive as-is given Dom's comments, but if it were to survive, then the water/ice composite here would only be done for z0rl and stress?

z0rl_cmp(i) = z0rl_ice(i)
ustar_cmp(i) = ustar_ice(i)
cm_cmp(i) = cm_ice(i)
ch_cmp(i) = ch_ice(i)
rb_cmp(i) = rb_ice(i)
stress_cmp(i) = stress_ice(i)
fm_cmp(i) = fm_ice(i)
fh_cmp(i) = fh_ice(i)
fm10_cmp(i) = fm10_ice(i)
fh2_cmp(i) = fh2_ice(i)
endif
endif ! end of if(flagiter) loop
enddo ! End of loop over horizontal
endif ! End of getting composite values for fractional grid

return
end subroutine sfc_diff_run
!> @}
Expand Down
Loading