Skip to content

Commit

Permalink
update GWD code for HR4
Browse files Browse the repository at this point in the history
  • Loading branch information
@Qingfu-Liu
Qingfu-Liu committed May 31, 2024
1 parent 622aa41 commit 6b78332
Show file tree
Hide file tree
Showing 6 changed files with 78 additions and 42 deletions.
61 changes: 35 additions & 26 deletions physics/GWD/drag_suite.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -219,8 +219,8 @@ subroutine drag_suite_run( &
& dusfc_ms,dvsfc_ms,dusfc_bl,dvsfc_bl, &
& dusfc_ss,dvsfc_ss,dusfc_fd,dvsfc_fd, &
& slmsk,br1,hpbl, &
& g, cp, rd, rv, fv, pi, imx, cdmbgwd, me, master, &
& lprnt, ipr, rdxzb, dx, gwd_opt, &
& g, cp, rd, rv, fv, pi, imx, cdmbgwd, alpha_fd, &
& me, master, lprnt, ipr, rdxzb, dx, gwd_opt, &
& do_gsl_drag_ls_bl, do_gsl_drag_ss, do_gsl_drag_tofd, &
& dtend, dtidx, index_of_process_orographic_gwd, &
& index_of_temperature, index_of_x_wind, &
Expand Down Expand Up @@ -327,18 +327,18 @@ subroutine drag_suite_run( &
integer, intent(in) :: gwd_opt
logical, intent(in) :: lprnt
integer, intent(in) :: KPBL(:)
real(kind=kind_phys), intent(in) :: deltim, G, CP, RD, RV, cdmbgwd(:)
real(kind=kind_phys), intent(in) :: deltim, G, CP, RD, RV, &
& cdmbgwd(:), alpha_fd
real(kind=kind_phys), intent(inout), optional :: dtend(:,:,:)
logical, intent(in) :: ldiag3d
integer, intent(in) :: dtidx(:,:)
integer, intent(in) :: index_of_temperature, &
integer, intent(in) :: dtidx(:,:), index_of_temperature, &
& index_of_process_orographic_gwd, index_of_x_wind, index_of_y_wind

integer :: kpblmax
integer, parameter :: ims=1, kms=1, its=1, kts=1
real(kind=kind_phys), intent(in) :: fv, pi
real(kind=kind_phys) :: rcl, cdmb
real(kind=kind_phys) :: g_inv
real(kind=kind_phys) :: g_inv, rd_inv

real(kind=kind_phys), intent(inout) :: &
& dudt(:,:),dvdt(:,:), &
Expand Down Expand Up @@ -444,6 +444,7 @@ subroutine drag_suite_run( &
real(kind=kind_phys), dimension(im,km) :: utendform,vtendform
real(kind=kind_phys) :: a1,a2,wsp
real(kind=kind_phys) :: H_efold
real(kind=kind_phys), parameter :: coeff_fd = 6.325e-3

! critical richardson number for wave breaking : ! larger drag with larger value
real(kind=kind_phys), parameter :: ric = 0.25
Expand Down Expand Up @@ -708,11 +709,12 @@ subroutine drag_suite_run( &
taufb(1:im,1:km+1) = 0.0
komax(1:im) = 0
!
rd_inv = 1./rd
do k = kts,km
do i = its,im
vtj(i,k) = t1(i,k) * (1.+fv*q1(i,k))
vtk(i,k) = vtj(i,k) / prslk(i,k)
ro(i,k) = 1./rd * prsl(i,k) / vtj(i,k) ! density kg/m**3
ro(i,k) = rd_inv * prsl(i,k) / vtj(i,k) ! density kg/m**3
enddo
enddo
!
Expand Down Expand Up @@ -1363,8 +1365,10 @@ subroutine drag_suite_run( &
H_efold = 1500.
DO k=kts,km
wsp=SQRT(uwnd1(i,k)**2 + vwnd1(i,k)**2)
! alpha*beta*Cmd*Ccorr*2.109 = 12.*1.*0.005*0.6*2.109 = 0.0759
var_temp = 0.0759*EXP(-(zl(i,k)/H_efold)**1.5)*a2* &
! Note: In Beljaars et al. (2004):
! alpha_fd*beta*Cmd*Ccorr*2.109 = 12.*1.*0.005*0.6*2.109 = 0.0759
! lump beta*Cmd*Ccorr*2.109 into 1.*0.005*0.6*2.109 = coeff_fd ~ 6.325e-3_kind_phys
var_temp = alpha_fd*coeff_fd*EXP(-(zl(i,k)/H_efold)**1.5)*a2* &
zl(i,k)**(-1.2)*ss_taper(i) ! this is greater than zero
! Note: This is a semi-implicit treatment of the time differencing
! per Beljaars et al. (2004, QJRMS)
Expand Down Expand Up @@ -1427,8 +1431,8 @@ subroutine drag_suite_psl( &
& dusfc_ls,dvsfc_ls,dusfc_bl,dvsfc_bl, &
& dusfc_ss,dvsfc_ss,dusfc_fd,dvsfc_fd, &
& slmsk,br1,hpbl,vtype, &
& g, cp, rd, rv, fv, pi, imx, cdmbgwd, me, master, &
& lprnt, ipr, rdxzb, dx, gwd_opt, &
& g, cp, rd, rv, fv, pi, imx, cdmbgwd, alpha_fd, &
& me, master, lprnt, ipr, rdxzb, dx, gwd_opt, &
& do_gsl_drag_ls_bl, do_gsl_drag_ss, do_gsl_drag_tofd, &
& psl_gwd_dx_factor, &
& dtend, dtidx, index_of_process_orographic_gwd, &
Expand Down Expand Up @@ -1537,7 +1541,8 @@ subroutine drag_suite_psl( &
integer, intent(in) :: gwd_opt
logical, intent(in) :: lprnt
integer, intent(in) :: KPBL(:)
real(kind=kind_phys), intent(in) :: deltim, G, CP, RD, RV, cdmbgwd(:)
real(kind=kind_phys), intent(in) :: deltim, G, CP, RD, RV, &
& cdmbgwd(:), alpha_fd
real(kind=kind_phys), intent(inout) :: dtend(:,:,:)
logical, intent(in) :: ldiag3d
integer, intent(in) :: dtidx(:,:), index_of_temperature, &
Expand All @@ -1547,7 +1552,7 @@ subroutine drag_suite_psl( &
integer, parameter :: ims=1, kms=1, its=1, kts=1
real(kind=kind_phys), intent(in) :: fv, pi
real(kind=kind_phys) :: rcl, cdmb
real(kind=kind_phys) :: g_inv
real(kind=kind_phys) :: g_inv, g_cp

real(kind=kind_phys), intent(inout) :: &
& dudt(:,:),dvdt(:,:), &
Expand Down Expand Up @@ -1649,6 +1654,7 @@ subroutine drag_suite_psl( &
real(kind=kind_phys), dimension(im,km) :: utendform,vtendform
real(kind=kind_phys) :: a1,a2,wsp
real(kind=kind_phys) :: H_efold
real(kind=kind_phys), parameter :: coeff_fd = 6.325e-3

! multification factor of standard deviation : ! larger drag with larger value
!!! real(kind=kind_phys), parameter :: psl_gwd_dx_factor = 6.0
Expand Down Expand Up @@ -1789,18 +1795,18 @@ subroutine drag_suite_psl( &
enddo

!--- calculate scale-aware tapering factors
do i=1,im
if ( dx(i) .ge. dxmax_ls ) then
ls_taper(i) = 1.
else
if ( dx(i) .le. dxmin_ls) then
ls_taper(i) = 0.
do i=1,im
if ( dx(i) .ge. dxmax_ls ) then
ls_taper(i) = 1.
else
ls_taper(i) = 0.5 * ( SIN(pi*(dx(i)-0.5*(dxmax_ls+dxmin_ls))/ &
if ( dx(i) .le. dxmin_ls) then
ls_taper(i) = 0.
else
ls_taper(i) = 0.5 * ( SIN(pi*(dx(i)-0.5*(dxmax_ls+dxmin_ls))/ &
(dxmax_ls-dxmin_ls)) + 1. )
endif
endif
endif
enddo
enddo

! Remove ss_tapering
ss_taper(:) = 1.
Expand Down Expand Up @@ -2072,6 +2078,7 @@ subroutine drag_suite_psl( &
IF ( (do_gsl_drag_ls_bl).and. &
((gwd_opt_ls .EQ. 1).or.(gwd_opt_bl .EQ. 1)) ) then

g_cp = g/cp
do i=its,im

if ( ls_taper(i).GT.1.E-02 ) then
Expand All @@ -2086,7 +2093,7 @@ subroutine drag_suite_psl( &
tem2 = v1(i,k) - v1(i,k+1)
dw2 = rcl*(tem1*tem1 + tem2*tem2)
shr2 = max(dw2,dw2min) * rdz * rdz
bvf2 = g*(g/cp+rdz*(vtj(i,k+1)-vtj(i,k))) * ti
bvf2 = g*(g_cp+rdz*(vtj(i,k+1)-vtj(i,k))) * ti
usqj(i,k) = max(bvf2/shr2,rimin)
bnv2(i,k) = 2.0*g*rdz*(vtk(i,k+1)-vtk(i,k))/(vtk(i,k+1)+vtk(i,k))
bnv2(i,k) = max( bnv2(i,k), bnv2min )
Expand Down Expand Up @@ -2341,8 +2348,10 @@ subroutine drag_suite_psl( &
H_efold = 1500.
DO k=kts,km
wsp=SQRT(u1(i,k)**2 + v1(i,k)**2)
! alpha*beta*Cmd*Ccorr*2.109 = 12.*1.*0.005*0.6*2.109 = 0.0759
var_temp = 0.0759*EXP(-(zl(i,k)/H_efold)**1.5)*a2* &
! Note: In Beljaars et al. (2004):
! alpha_fd*beta*Cmd*Ccorr*2.109 = 12.*1.*0.005*0.6*2.109 = 0.0759
! lump beta*Cmd*Ccorr*2.109 into 1.*0.005*0.6*2.109 = coeff_fd ~ 6.325e-3_kind_phys
var_temp = alpha_fd*coeff_fd*EXP(-(zl(i,k)/H_efold)**1.5)*a2* &
zl(i,k)**(-1.2)*ss_taper(i) ! this is greater than zero
! Note: This is a semi-implicit treatment of the time differencing
! per Beljaars et al. (2004, QJRMS)
Expand Down Expand Up @@ -2465,7 +2474,7 @@ subroutine drag_suite_psl( &
do k = km, kpblmin, -1
if(flag(i).and. k.le.kbmax(i)) then
pe = pe + bnv2(i,k)*(zl(i,kbmax(i))-zl(i,k))* &
del(i,k)/g/ro(i,k)
del(i,k)*g_inv/ro(i,k)
ke = 0.5*((rcs*u1(i,k))**2.+(rcs*v1(i,k))**2.)
!
!---------- apply flow-blocking drag when pe >= ke
Expand Down
10 changes: 9 additions & 1 deletion physics/GWD/drag_suite.meta
View file
Original file line number Diff line number Diff line change
Expand Up @@ -528,6 +528,14 @@
type = real
kind = kind_phys
intent = in
[alpha_fd]
standard_name = alpha_coefficient_for_turbulent_orographic_form_drag
long_name = alpha coefficient for Beljaars et al turbulent orographic form drag
units = none
dimensions = ()
type = real
kind = kind_phys
intent = in
[me]
standard_name = mpi_rank
long_name = rank of the current MPI task
Expand Down Expand Up @@ -603,7 +611,7 @@
[psl_gwd_dx_factor]
standard_name = effective_grid_spacing_of_psl_gwd_suite
long_name = multiplication of grid spacing
units = count
units = 1
dimensions = ()
type = real
kind = kind_phys
Expand Down
13 changes: 7 additions & 6 deletions physics/GWD/ugwpv1_gsldrag.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -309,7 +309,7 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ak, bk, ntrac, lonr, dtp,
do_gwd_opt_psl, psl_gwd_dx_factor, &
do_ugwp_v1, do_ugwp_v1_orog_only, &
do_ugwp_v1_w_gsldrag, gwd_opt, do_tofd, ldiag_ugwp, ugwp_seq_update, &
cdmbgwd, jdat, nmtvr, hprime, oc, theta, sigma, gamma, &
cdmbgwd, alpha_fd, jdat, nmtvr, hprime, oc, theta, sigma, gamma, &
elvmax, clx, oa4, varss,oc1ss,oa4ss,ol4ss, dx, xlat, xlat_d, sinlat, coslat, &
area, rain, br1, hpbl,vtype, kpbl, slmsk, &
ugrs, vgrs, tgrs, q1, prsi, prsl, prslk, phii, phil, del, tau_amf, &
Expand Down Expand Up @@ -375,7 +375,7 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ak, bk, ntrac, lonr, dtp,
! SSO parameters and variables
integer, intent(in) :: gwd_opt !gwd_opt and nmtvr are "redundant" controls
integer, intent(in) :: nmtvr
real(kind=kind_phys), intent(in) :: cdmbgwd(:) ! for gsl_drag
real(kind=kind_phys), intent(in) :: cdmbgwd(:), alpha_fd ! for gsl_drag

real(kind=kind_phys), intent(in), dimension(:) :: hprime, oc, theta, sigma, gamma

Expand Down Expand Up @@ -404,8 +404,7 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ak, bk, ntrac, lonr, dtp,
integer, intent(in), dimension(:) :: vtype

real(kind=kind_phys), intent(in), dimension(:) :: rain
real(kind=kind_phys), intent(in), dimension(:) :: br1, slmsk
real(kind=kind_phys), intent(in), dimension(:) :: hpbl
real(kind=kind_phys), intent(in), dimension(:) :: br1, hpbl, slmsk
!
! moved to GFS_phys_time_vary
! real(kind=kind_phys), intent(in), dimension(:) :: ddy_j1tau, ddy_j2tau
Expand Down Expand Up @@ -563,7 +562,8 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ak, bk, ntrac, lonr, dtp,
du_osscol, dv_osscol, du_ofdcol, dv_ofdcol, &
slmsk,br1,hpbl,vtype,con_g,con_cp,con_rd,con_rv, &
con_fv, con_pi, lonr, &
cdmbgwd(1:2),me,master,lprnt,ipr,rdxzb,dx,gwd_opt, &
cdmbgwd(1:2),alpha_fd,me,master, &
lprnt,ipr,rdxzb,dx,gwd_opt, &
do_gsl_drag_ls_bl,do_gsl_drag_ss,do_gsl_drag_tofd, &
psl_gwd_dx_factor, &
dtend, dtidx, index_of_process_orographic_gwd, &
Expand All @@ -583,7 +583,8 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ak, bk, ntrac, lonr, dtp,
du_osscol, dv_osscol, du_ofdcol, dv_ofdcol, &
slmsk,br1,hpbl,con_g,con_cp,con_rd,con_rv, &
con_fv, con_pi, lonr, &
cdmbgwd(1:2),me,master,lprnt,ipr,rdxzb,dx,gwd_opt, &
cdmbgwd(1:2),alpha_fd,me,master, &
lprnt,ipr,rdxzb,dx,gwd_opt, &
do_gsl_drag_ls_bl,do_gsl_drag_ss,do_gsl_drag_tofd, &
dtend, dtidx, index_of_process_orographic_gwd, &
index_of_temperature, index_of_x_wind, &
Expand Down
10 changes: 9 additions & 1 deletion physics/GWD/ugwpv1_gsldrag.meta
View file
Original file line number Diff line number Diff line change
Expand Up @@ -412,7 +412,7 @@
[psl_gwd_dx_factor]
standard_name = effective_grid_spacing_of_psl_gwd_suite
long_name = multiplication of grid spacing
units = count
units = 1
dimensions = ()
type = real
kind = kind_phys
Expand Down Expand Up @@ -474,6 +474,14 @@
type = real
kind = kind_phys
intent = in
[alpha_fd]
standard_name = alpha_coefficient_for_turbulent_orographic_form_drag
long_name = alpha coefficient for Beljaars et al turbulent orographic form drag
units = none
dimensions = ()
type = real
kind = kind_phys
intent = in
[jdat]
standard_name = date_and_time_of_forecast_in_united_states_order
long_name = current forecast date and time
Expand Down
16 changes: 9 additions & 7 deletions physics/GWD/unified_ugwp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -250,7 +250,7 @@ subroutine unified_ugwp_run(me, master, im, levs, ak,bk, ntrac, dtp, fhzero, kdt
dvsfc_ss,dusfc_fd,dvsfc_fd,dtaux2d_ms,dtauy2d_ms,dtaux2d_bl,dtauy2d_bl, &
dtaux2d_ss,dtauy2d_ss,dtaux2d_fd,dtauy2d_fd,dudt_ngw,dvdt_ngw,dtdt_ngw, &
br1,hpbl,vtype,slmsk, do_tofd, ldiag_ugwp, ugwp_seq_update, &
cdmbgwd, jdat, xlat, xlat_d, sinlat, coslat, area, &
cdmbgwd, alpha_fd, jdat, xlat, xlat_d, sinlat, coslat, area, &
ugrs, vgrs, tgrs, q1, prsi, prsl, prslk, phii, phil, &
del, kpbl, dusfcg, dvsfcg, gw_dudt, gw_dvdt, gw_dtdt, gw_kdis, &
tau_tofd, tau_mtb, tau_ogw, tau_ngw, zmtb, zlwb, zogw, &
Expand Down Expand Up @@ -290,7 +290,7 @@ subroutine unified_ugwp_run(me, master, im, levs, ak,bk, ntrac, dtp, fhzero, kdt
real(kind=kind_phys), intent(in), dimension(:,:) :: del, ugrs, vgrs, tgrs, prsl, prslk, phil
real(kind=kind_phys), intent(in), dimension(:,:) :: prsi, phii
real(kind=kind_phys), intent(in), dimension(:,:) :: q1
real(kind=kind_phys), intent(in) :: dtp, fhzero, cdmbgwd(:)
real(kind=kind_phys), intent(in) :: dtp, fhzero, cdmbgwd(:), alpha_fd
integer, intent(in) :: jdat(:)
logical, intent(in) :: do_tofd, ldiag_ugwp, ugwp_seq_update

Expand Down Expand Up @@ -495,7 +495,7 @@ subroutine unified_ugwp_run(me, master, im, levs, ak,bk, ntrac, dtp, fhzero, kdt

if ( do_gsl_drag_ls_bl.or.do_gsl_drag_ss.or.do_gsl_drag_tofd ) then
!
if (do_gwd_opt_psl) then
if (do_gwd_opt_psl) then
call drag_suite_psl(im,levs,dvdt,dudt,dtdt,uwnd1,vwnd1, &
tgrs,q1,kpbl,prsi,del,prsl,prslk,phii,phil,dtp, &
kdt,hprime,oc,oa4,clx,varss,oc1ss,oa4ss, &
Expand All @@ -506,14 +506,15 @@ subroutine unified_ugwp_run(me, master, im, levs, ak,bk, ntrac, dtp, fhzero, kdt
dusfc_ss,dvsfc_ss,dusfc_fd,dvsfc_fd, &
slmsk,br1,hpbl,vtype,con_g,con_cp,con_rd,con_rv, &
con_fvirt,con_pi,lonr, &
cdmbgwd,me,master,lprnt,ipr,rdxzb,dx,gwd_opt, &
cdmbgwd,alpha_fd,me,master, &
lprnt,ipr,rdxzb,dx,gwd_opt, &
do_gsl_drag_ls_bl,do_gsl_drag_ss,do_gsl_drag_tofd, &
psl_gwd_dx_factor, &
dtend, dtidx, index_of_process_orographic_gwd, &
index_of_temperature, index_of_x_wind, &
index_of_y_wind, ldiag3d, ldiag_ugwp, &
ugwp_seq_update, spp_wts_gwd, spp_gwd, errmsg, errflg)
else
else
call drag_suite_run(im,levs,dvdt,dudt,dtdt,uwnd1,vwnd1, &
tgrs,q1,kpbl,prsi,del,prsl,prslk,phii,phil,dtp, &
kdt,hprime,oc,oa4,clx,varss,oc1ss,oa4ss, &
Expand All @@ -524,13 +525,14 @@ subroutine unified_ugwp_run(me, master, im, levs, ak,bk, ntrac, dtp, fhzero, kdt
dusfc_ss,dvsfc_ss,dusfc_fd,dvsfc_fd, &
slmsk,br1,hpbl,con_g,con_cp,con_rd,con_rv, &
con_fvirt,con_pi,lonr, &
cdmbgwd,me,master,lprnt,ipr,rdxzb,dx,gwd_opt, &
cdmbgwd,alpha_fd,me,master, &
lprnt,ipr,rdxzb,dx,gwd_opt, &
do_gsl_drag_ls_bl,do_gsl_drag_ss,do_gsl_drag_tofd, &
dtend, dtidx, index_of_process_orographic_gwd, &
index_of_temperature, index_of_x_wind, &
index_of_y_wind, ldiag3d, ldiag_ugwp, &
ugwp_seq_update, spp_wts_gwd, spp_gwd, errmsg, errflg)
endif
endif
!
! put zeros due to xy GSL-drag style: dtaux2d_bl,dtauy2d_bl,dtaux2d_ss.......dusfc_ms,dvsfc_ms
!
Expand Down
10 changes: 9 additions & 1 deletion physics/GWD/unified_ugwp.meta
View file
Original file line number Diff line number Diff line change
Expand Up @@ -715,6 +715,14 @@
type = real
kind = kind_phys
intent = in
[alpha_fd]
standard_name = alpha_coefficient_for_turbulent_orographic_form_drag
long_name = alpha coefficient for Beljaars et al turbulent orographic form drag
units = none
dimensions = ()
type = real
kind = kind_phys
intent = in
[jdat]
standard_name = date_and_time_of_forecast_in_united_states_order
long_name = current forecast date and time
Expand Down Expand Up @@ -1262,7 +1270,7 @@
[psl_gwd_dx_factor]
standard_name = effective_grid_spacing_of_psl_gwd_suite
long_name = multiplication of grid spacing
units = count
units = 1
dimensions = ()
type = real
kind = kind_phys
Expand Down

0 comments on commit 6b78332

Please sign in to comment.