| cplchm | gfs_control_type | logical flag for chemistry collection | .false.
diff --git a/physics/drag_suite.F90 b/physics/drag_suite.F90
index 56902c631..eb371adb1 100644
--- a/physics/drag_suite.F90
+++ b/physics/drag_suite.F90
@@ -485,7 +485,7 @@ subroutine drag_suite_run( &
varmax_fd = 150., &
beta_ss = 0.1, &
beta_fd = 0.2
- real(kind=kind_phys) :: var_temp
+ real(kind=kind_phys) :: var_temp, var_temp2
! added Beljaars orographic form drag
real(kind=kind_phys), dimension(im,km) :: utendform,vtendform
@@ -1060,7 +1060,9 @@ subroutine drag_suite_run( &
!tauwavex0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*u1(i,3)
var_temp = MIN(varss(i),varmax_ss) + &
MAX(0.,beta_ss*(varss(i)-varmax_ss))
- tauwavex0=0.5*XNBV*xlinv(i)*(2.*var_temp)**2*ro(i,kvar)*u1(i,kvar)
+ ! Note: This is a semi-implicit treatment of the time differencing
+ var_temp2 = 0.5*XNBV*xlinv(i)*(2.*var_temp)**2*ro(i,kvar) ! this is greater than zero
+ tauwavex0=-var_temp2*u1(i,kvar)/(1.+var_temp2*deltim)
tauwavex0=tauwavex0*ss_taper
else
tauwavex0=0.
@@ -1073,7 +1075,8 @@ subroutine drag_suite_run( &
!tauwavey0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*v1(i,3)
var_temp = MIN(varss(i),varmax_ss) + &
MAX(0.,beta_ss*(varss(i)-varmax_ss))
- tauwavey0=0.5*XNBV*xlinv(i)*(2.*var_temp)**2*ro(i,kvar)*v1(i,kvar)
+ ! Note: This is a semi-implicit treatment of the time differencing
+ tauwavey0=-var_temp2*v1(i,kvar)/(1.+var_temp2*deltim)
tauwavey0=tauwavey0*ss_taper
else
tauwavey0=0.
@@ -1154,10 +1157,12 @@ subroutine drag_suite_run( &
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
- utendform(i,k)=-0.0759*wsp*u1(i,k)* &
- EXP(-(zl(i,k)/H_efold)**1.5)*a2*zl(i,k)**(-1.2)*ss_taper
- vtendform(i,k)=-0.0759*wsp*v1(i,k)* &
- EXP(-(zl(i,k)/H_efold)**1.5)*a2*zl(i,k)**(-1.2)*ss_taper
+ var_temp = 0.0759*EXP(-(zl(i,k)/H_efold)**1.5)*a2* &
+ zl(i,k)**(-1.2)*ss_taper ! this is greater than zero
+ ! Note: This is a semi-implicit treatment of the time differencing
+ ! per Beljaars et al. (2004, QJRMS)
+ utendform(i,k) = - var_temp*wsp*u1(i,k)/(1. + var_temp*deltim*wsp)
+ vtendform(i,k) = - var_temp*wsp*v1(i,k)/(1. + var_temp*deltim*wsp)
!IF(zl(i,k) > 4000.) exit
ENDDO
ENDIF
diff --git a/physics/gfdl_cloud_microphys.F90 b/physics/gfdl_cloud_microphys.F90
index fcaaa9b94..1ccedb956 100644
--- a/physics/gfdl_cloud_microphys.F90
+++ b/physics/gfdl_cloud_microphys.F90
@@ -113,7 +113,7 @@ end subroutine gfdl_cloud_microphys_finalize
!! \htmlinclude gfdl_cloud_microphys_run.html
!!
subroutine gfdl_cloud_microphys_run( &
- levs, im, con_g, con_fvirt, con_rd, frland, garea, &
+ levs, im, con_g, con_fvirt, con_rd, frland, garea, islmsk, &
gq0, gq0_ntcw, gq0_ntrw, gq0_ntiw, gq0_ntsw, gq0_ntgl, gq0_ntclamt, &
gt0, gu0, gv0, vvl, prsl, phii, del, &
rain0, ice0, snow0, graupel0, prcp0, sr, &
@@ -136,6 +136,7 @@ subroutine gfdl_cloud_microphys_run( &
integer, intent(in ) :: levs, im
real(kind=kind_phys), intent(in ) :: con_g, con_fvirt, con_rd
real(kind=kind_phys), intent(in ), dimension(1:im) :: frland, garea
+ integer, intent(in ), dimension(1:im) :: islmsk
real(kind=kind_phys), intent(inout), dimension(1:im,1:levs) :: gq0, gq0_ntcw, gq0_ntrw, gq0_ntiw, &
gq0_ntsw, gq0_ntgl, gq0_ntclamt
real(kind=kind_phys), intent(inout), dimension(1:im,1:levs) :: gt0, gu0, gv0
@@ -298,9 +299,11 @@ subroutine gfdl_cloud_microphys_run( &
enddo
enddo
call cloud_diagnosis (1, im, 1, levs, den(1:im,1:levs), &
+ del(1:im,1:levs), islmsk(1:im), &
gq0_ntcw(1:im,1:levs), gq0_ntiw(1:im,1:levs), &
- gq0_ntrw(1:im,1:levs), gq0_ntsw(1:im,1:levs), &
- gq0_ntgl(1:im,1:levs), gt0(1:im,1:levs), &
+ gq0_ntrw(1:im,1:levs), &
+ gq0_ntsw(1:im,1:levs) + gq0_ntgl(1:im,1:levs), &
+ gq0_ntgl(1:im,1:levs)*0.0, gt0(1:im,1:levs), &
rew(1:im,1:levs), rei(1:im,1:levs), rer(1:im,1:levs),&
res(1:im,1:levs), reg(1:im,1:levs))
deallocate(den)
diff --git a/physics/gfdl_cloud_microphys.meta b/physics/gfdl_cloud_microphys.meta
index c2ce3f8f5..7f31637bf 100644
--- a/physics/gfdl_cloud_microphys.meta
+++ b/physics/gfdl_cloud_microphys.meta
@@ -180,6 +180,14 @@
kind = kind_phys
intent = in
optional = F
+[islmsk]
+ standard_name = sea_land_ice_mask
+ long_name = sea/land/ice mask (=0/1/2)
+ units = flag
+ dimensions = (horizontal_dimension)
+ type = integer
+ intent = in
+ optional = F
[gq0]
standard_name = water_vapor_specific_humidity_updated_by_physics
long_name = water vapor specific humidity updated by physics
diff --git a/physics/gfdl_fv_sat_adj.F90 b/physics/gfdl_fv_sat_adj.F90
index 14b3975f3..f5c84cd99 100644
--- a/physics/gfdl_fv_sat_adj.F90
+++ b/physics/gfdl_fv_sat_adj.F90
@@ -49,7 +49,7 @@ module fv_sat_adj
! | gfdl_cloud_microphys_mod |
! ql_gen, qi_gen, qi0_max, ql_mlt, ql0_max, qi_lim, qs_mlt,
! tau_r2g, tau_smlt, tau_i2s, tau_v2l, tau_l2v, tau_imlt, tau_l2r,
-! rad_rain, rad_snow, rad_graupel, dw_ocean, dw_land |
+! rad_rain, rad_snow, rad_graupel, dw_ocean, dw_land, tintqs
!
!
! DH* TODO - MAKE THIS INPUT ARGUMENTS *DH
@@ -64,8 +64,7 @@ module fv_sat_adj
use gfdl_cloud_microphys_mod, only: ql_gen, qi_gen, qi0_max, ql_mlt, ql0_max, qi_lim, qs_mlt
use gfdl_cloud_microphys_mod, only: icloud_f, sat_adj0, t_sub, cld_min
use gfdl_cloud_microphys_mod, only: tau_r2g, tau_smlt, tau_i2s, tau_v2l, tau_l2v, tau_imlt, tau_l2r
- use gfdl_cloud_microphys_mod, only: rad_rain, rad_snow, rad_graupel, dw_ocean, dw_land
-
+ use gfdl_cloud_microphys_mod, only: rad_rain, rad_snow, rad_graupel, dw_ocean, dw_land, tintqs
#ifdef MULTI_GASES
use ccpp_multi_gases_mod, only: multi_gases_init, &
multi_gases_finalize, &
@@ -1030,9 +1029,13 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te,
do i = is, ie
+ if(tintqs) then
+ tin = pt1(i)
+ else
tin = pt1 (i) - (lcp2 (i) * q_cond (i) + icp2 (i) * q_sol (i)) ! minimum temperature
! tin = pt1 (i) - ((lv00 + d0_vap * pt1 (i)) * q_cond (i) + &
! (li00 + dc_ice * pt1 (i)) * q_sol (i)) / (mc_air (i) + qpz (i) * c_vap)
+ endif
! -----------------------------------------------------------------------
! determine saturated specific humidity
@@ -1075,14 +1078,14 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te,
! icloud_f = 2: binary cloud scheme (0 / 1)
! -----------------------------------------------------------------------
- if (rh > 0.75 .and. qpz (i) > 1.e-6) then
+ if (rh > 0.75 .and. qpz (i) > 1.e-8) then
dq = hvar (i) * qpz (i)
q_plus = qpz (i) + dq
q_minus = qpz (i) - dq
if (icloud_f == 2) then
if (qpz (i) > qstar (i)) then
qa (i, j) = 1.
- elseif (qstar (i) < q_plus .and. q_cond (i) > 1.e-6) then
+ elseif (qstar (i) < q_plus .and. q_cond (i) > 1.e-8) then
qa (i, j) = ((q_plus - qstar (i)) / dq) ** 2
qa (i, j) = min (1., qa (i, j))
else
@@ -1102,7 +1105,7 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te,
qa (i, j) = 0.
endif
! impose minimum cloudiness if substantial q_cond (i) exist
- if (q_cond (i) > 1.e-6) then
+ if (q_cond (i) > 1.e-8) then
qa (i, j) = max (cld_min, qa (i, j))
endif
qa (i, j) = min (1., qa (i, j))
diff --git a/physics/gwdc.f b/physics/gwdc.f
index 80898c47b..9909a3100 100644
--- a/physics/gwdc.f
+++ b/physics/gwdc.f
@@ -22,7 +22,7 @@ end subroutine gwdc_pre_init
subroutine gwdc_pre_run ( &
& im, cgwf, dx, work1, work2, dlength, cldf, &
& levs, kbot, ktop, dtp, gt0, gt0_init, del, cumabs, &
- & errmsg, errflg )
+ & do_cnvgwd, errmsg, errflg )
use machine, only : kind_phys
implicit none
@@ -38,6 +38,7 @@ subroutine gwdc_pre_run ( &
real(kind=kind_phys), intent(out) :: &
& dlength(:), cldf(:), cumabs(:)
+ logical, intent(in) :: do_cnvgwd
character(len=*), intent(out) :: errmsg
integer, intent(out) :: errflg
@@ -49,6 +50,14 @@ subroutine gwdc_pre_run ( &
errmsg = ''
errflg = 0
+ ! DH*
+ if (.not. do_cnvgwd) then
+ write(0,*) "ERROR: , GWDC_PRE CALLED BUT DO_CNVGWD FALSE"
+ call sleep(5)
+ stop
+ end if
+ ! *DH
+
do i = 1, im
tem1 = dx(i)
tem2 = tem1
diff --git a/physics/gwdc.meta b/physics/gwdc.meta
index b87529aec..2151cc5f7 100644
--- a/physics/gwdc.meta
+++ b/physics/gwdc.meta
@@ -137,6 +137,14 @@
kind = kind_phys
intent = out
optional = F
+[do_cnvgwd]
+ standard_name = flag_for_convective_gravity_wave_drag
+ long_name = flag for convective gravity wave drag (gwd)
+ units = flag
+ dimensions = ()
+ type = logical
+ intent = in
+ optional = F
[errmsg]
standard_name = ccpp_error_message
long_name = error message for error handling in CCPP
diff --git a/physics/gwdps.f b/physics/gwdps.f
index 366a8b974..0ea2c8754 100644
--- a/physics/gwdps.f
+++ b/physics/gwdps.f
@@ -299,12 +299,8 @@ subroutine gwdps_run( &
! Interface variables
integer, intent(in) :: im, ix, km, imx, kdt, ipr, me
integer, intent(in) :: KPBL(IM) ! Index for the PBL top layer!
- ! DH* adding intent(in) information for the following variables
- ! changes the results on Theia/Intel - skip for bit-for-bit results *DH
-! real(kind=kind_phys), intent(in) :: &
-! & deltim, G, CP, RD, RV, cdmbgwd(2)
- real(kind=kind_phys) deltim, G, CP, RD, RV, cdmbgwd(2)
- ! *DH
+ real(kind=kind_phys), intent(in) :: &
+ & deltim, G, CP, RD, RV, cdmbgwd(4)
real(kind=kind_phys), intent(inout) :: &
& A(IX,KM), B(IX,KM), C(IX,KM)
real(kind=kind_phys), intent(in) :: &
@@ -382,7 +378,8 @@ subroutine gwdps_run( &
real(kind=kind_phys) TAUB(IM), XN(IM), YN(IM), UBAR(IM) &
&, VBAR(IM), ULOW(IM), OA(IM), CLX(IM) &
&, ROLL(IM), ULOI(IM) &
- &, DTFAC(IM), XLINV(IM), DELKS(IM), DELKS1(IM)
+ &, DTFAC(IM), XLINV(IM), DELKS(IM)
+! &, DTFAC(IM), XLINV(IM), DELKS(IM), DELKS1(IM)
!
real(kind=kind_phys) BNV2(IM,KM), TAUP(IM,KM+1), ri_n(IM,KM) &
&, TAUD(IM,KM), RO(IM,KM), VTK(IM,KM) &
@@ -392,7 +389,8 @@ subroutine gwdps_run( &
! real(kind=kind_phys) VELKO(KM-1)
integer kref(IM), kint(im), iwk(im), ipt(im)
! for lm mtn blocking
- integer kreflm(IM), iwklm(im)
+ integer iwklm(im)
+! integer kreflm(IM), iwklm(im)
integer idxzb(im), ktrial, klevm1
!
real(kind=kind_phys) gor, gocp, fv, gr2, bnv, fr &
@@ -470,7 +468,7 @@ subroutine gwdps_run( &
do i=1,npt
iwklm(i) = 2
IDXZB(i) = 0
- kreflm(i) = 0
+! kreflm(i) = 0
enddo
! if (lprnt)
! & print *,' in gwdps_lm.f npt,IM,IX,IY,km,me=',npt,IM,IX,IY,km,me
@@ -552,14 +550,14 @@ subroutine gwdps_run( &
!
DO I = 1, npt
J = ipt(i)
- DELKS(I) = 1.0 / (PRSI(J,1) - PRSI(J,iwklm(i)))
- DELKS1(I) = 1.0 / (PRSL(J,1) - PRSL(J,iwklm(i)))
- UBAR (I) = 0.0
- VBAR (I) = 0.0
- ROLL (I) = 0.0
- PE (I) = 0.0
- EK (I) = 0.0
- BNV2bar(I) = (PRSL(J,1)-PRSL(J,2)) * DELKS1(I) * BNV2LM(I,1)
+ DELKS(I) = 1.0 / (PRSI(J,1) - PRSI(J,iwklm(i)))
+! DELKS1(I) = 1.0 / (PRSI(J,1) - PRSL(J,iwklm(i)))
+ UBAR (I) = 0.0
+ VBAR (I) = 0.0
+ ROLL (I) = 0.0
+ PE (I) = 0.0
+ EK (I) = 0.0
+ BNV2bar(I) = (PRSI(J,1)-PRSL(J,1)) * DELKS(I) * BNV2LM(I,1)
ENDDO
! --- find the dividing stream line height
@@ -567,13 +565,13 @@ subroutine gwdps_run( &
! --- iwklm(i) is the k-index of mtn elvmax elevation
!> - Find the dividing streamline height starting from the level above
!! the maximum mountain height and processing downward.
- DO Ktrial = KMLL, 1, -1
- DO I = 1, npt
- IF ( Ktrial < iwklm(I) .and. kreflm(I) == 0 ) then
- kreflm(I) = Ktrial
- ENDIF
- ENDDO
- ENDDO
+! DO Ktrial = KMLL, 1, -1
+! DO I = 1, npt
+! IF ( Ktrial < iwklm(I) .and. kreflm(I) == 0 ) then
+! kreflm(I) = Ktrial
+! ENDIF
+! ENDDO
+! ENDDO
! print *,' in gwdps_lm.f 4 npt=',npt,kreflm(npt),me
!
! --- in the layer kreflm(I) to 1 find PE (which needs N, ELVMAX)
@@ -582,13 +580,17 @@ subroutine gwdps_run( &
! --- is the vert ave of quantities from the surface to mtn top.
!
DO I = 1, npt
- DO K = 1, Kreflm(I)
+ DO K = 1, iwklm(i)-1
J = ipt(i)
RDELKS = DEL(J,K) * DELKS(I)
UBAR(I) = UBAR(I) + RDELKS * U1(J,K) ! trial Mean U below
VBAR(I) = VBAR(I) + RDELKS * V1(J,K) ! trial Mean V below
ROLL(I) = ROLL(I) + RDELKS * RO(I,K) ! trial Mean RO below
- RDELKS = (PRSL(J,K)-PRSL(J,K+1)) * DELKS1(I)
+ if (k < iwklm(I)-1) then
+ RDELKS = (PRSL(J,K)-PRSL(J,K+1)) * DELKS(I)
+ else
+ RDELKS = (PRSL(J,K)-PRSI(J,K+1)) * DELKS(I)
+ endif
BNV2bar(I) = BNV2bar(I) + BNV2lm(I,K) * RDELKS
! --- these vert ave are for diags, testing and GWD to follow (*j*).
ENDDO
@@ -862,14 +864,14 @@ subroutine gwdps_run( &
J = ipt(i)
kref(I) = MAX(IWK(I), KPBL(J)+1 ) ! reference level
DELKS(I) = 1.0 / (PRSI(J,1) - PRSI(J,kref(I)))
- DELKS1(I) = 1.0 / (PRSL(J,1) - PRSL(J,kref(I)))
+! DELKS1(I) = 1.0 / (PRSI(J,1) - PRSL(J,kref(I)))
UBAR (I) = 0.0
VBAR (I) = 0.0
ROLL (I) = 0.0
KBPS = MAX(KBPS, kref(I))
KMPS = MIN(KMPS, kref(I))
!
- BNV2bar(I) = (PRSL(J,1)-PRSL(J,2)) * DELKS1(I) * BNV2(I,1)
+ BNV2bar(I) = (PRSI(J,1)-PRSL(J,1)) * DELKS(I) * BNV2(I,1)
ENDDO
! print *,' in gwdps_lm.f GWD:15 =',KBPS,KMPS
KBPSP1 = KBPS + 1
@@ -883,7 +885,11 @@ subroutine gwdps_run( &
VBAR(I) = VBAR(I) + RDELKS * V1(J,K) ! Mean V below kref
!
ROLL(I) = ROLL(I) + RDELKS * RO(I,K) ! Mean RO below kref
- RDELKS = (PRSL(J,K)-PRSL(J,K+1)) * DELKS1(I)
+ if (k < kref(i)-1) then
+ RDELKS = (PRSL(J,K)-PRSL(J,K+1)) * DELKS(I)
+ else
+ RDELKS = (PRSL(J,K)-PRSI(J,K+1)) * DELKS(I)
+ endif
BNV2bar(I) = BNV2bar(I) + BNV2(I,K) * RDELKS
ENDIF
ENDDO
diff --git a/physics/gwdps.meta b/physics/gwdps.meta
index 97b6abae3..0a141b208 100644
--- a/physics/gwdps.meta
+++ b/physics/gwdps.meta
@@ -318,7 +318,7 @@
standard_name = multiplication_factors_for_mountain_blocking_and_orographic_gravity_wave_drag
long_name = multiplic. factors for (1) mountain blocking drag coeff. and (2) ref. level orographic gravity wave drag
units = none
- dimensions = (2)
+ dimensions = (4)
type = real
kind = kind_phys
intent = in
diff --git a/physics/machine.F b/physics/machine.F
index ea6198c33..896b665da 100644
--- a/physics/machine.F
+++ b/physics/machine.F
@@ -1,10 +1,8 @@
module machine
-#if 0
!! \section arg_table_machine
!! \htmlinclude machine.html
!!
-#endif
implicit none
diff --git a/physics/module_gfdl_cloud_microphys.F90 b/physics/module_gfdl_cloud_microphys.F90
index ac3795566..2f6e5ec1a 100644
--- a/physics/module_gfdl_cloud_microphys.F90
+++ b/physics/module_gfdl_cloud_microphys.F90
@@ -1,6 +1,9 @@
!> \file gfdl_cloud_microphys.F90
-!! This file contains the column GFDL cloud microphysics ( Chen and Lin (2013)
-!! \cite chen_and_lin_2013 ).
+!! This file contains the full GFDL cloud microphysics (Chen and Lin (2013)
+!! \cite chen_and_lin_2013 and Zhou et al. 2019 \cite zhou2019toward).
+!! The module is paired with 'gfdl_fv_sat_adj', which performs the "fast"
+!! processes
+!>author Shian-Jiann Lin, Linjiong Zhou
!***********************************************************************
!* GNU Lesser General Public License
!*
@@ -285,6 +288,18 @@ module gfdl_cloud_microphys_mod
real :: log_10, tice0, t_wfr
+ integer :: reiflag = 1
+ ! 1: Heymsfield and Mcfarquhar, 1996
+ ! 2: Wyser, 1998
+
+ logical :: tintqs = .false. !< use temperature in the saturation mixing in PDF
+
+ real :: rewmin = 5.0, rewmax = 10.0
+ real :: reimin = 10.0, reimax = 150.0
+ real :: rermin = 10.0, rermax = 10000.0
+ real :: resmin = 150.0, resmax = 10000.0
+ real :: regmin = 300.0, regmax = 10000.0
+
! -----------------------------------------------------------------------
! namelist
! -----------------------------------------------------------------------
@@ -299,7 +314,9 @@ module gfdl_cloud_microphys_mod
tau_i2s, tau_l2r, qi_lim, ql_gen, c_paut, c_psaci, c_pgacs, &
z_slope_liq, z_slope_ice, prog_ccn, c_cracw, alin, clin, tice, &
rad_snow, rad_graupel, rad_rain, cld_min, use_ppm, mono_prof, &
- do_sedi_heat, sedi_transport, do_sedi_w, de_ice, icloud_f, irain_f, mp_print
+ do_sedi_heat, sedi_transport, do_sedi_w, de_ice, icloud_f, irain_f, &
+ mp_print, reiflag, rewmin, rewmax, reimin, reimax, rermin, rermax, &
+ resmin, resmax, regmin, regmax, tintqs
public &
mp_time, t_min, t_sub, tau_r2g, tau_smlt, tau_g2r, dw_land, dw_ocean, &
@@ -311,7 +328,9 @@ module gfdl_cloud_microphys_mod
tau_i2s, tau_l2r, qi_lim, ql_gen, c_paut, c_psaci, c_pgacs, &
z_slope_liq, z_slope_ice, prog_ccn, c_cracw, alin, clin, tice, &
rad_snow, rad_graupel, rad_rain, cld_min, use_ppm, mono_prof, &
- do_sedi_heat, sedi_transport, do_sedi_w, de_ice, icloud_f, irain_f, mp_print
+ do_sedi_heat, sedi_transport, do_sedi_w, de_ice, icloud_f, irain_f, &
+ mp_print, reiflag, rewmin, rewmax, reimin, reimax, rermin, rermax, &
+ resmin, resmax, regmin, regmax, tintqs
contains
@@ -3301,7 +3320,7 @@ subroutine fall_speed (ktop, kbot, den, qs, qi, qg, ql, tk, vts, vti, vtg)
else
tc (k) = tk (k) - tice
vti (k) = (3. + log10 (qi (k) * den (k))) * (tc (k) * (aa * tc (k) + bb) + cc) + dd * tc (k) + ee
- vti (k) = vi0 * exp (log_10 * vti (k))
+ vti (k) = vi0 * exp (log_10 * vti (k)) * 0.8
vti (k) = min (vi_max, max (vf_min, vti (k)))
endif
enddo
@@ -4683,127 +4702,141 @@ end subroutine interpolate_z
!> \ingroup mod_gfdl_cloud_mp
!! The subroutine 'cloud_diagnosis' diagnoses the radius of cloud
!! species.
-subroutine cloud_diagnosis (is, ie, js, je, den, qw, qi, qr, qs, qg, t, &
+!>author Linjiong Zhoum, Shian-Jiann Lin
+! =======================================================================
+subroutine cloud_diagnosis (is, ie, ks, ke, den, delp, lsm, qmw, qmi, qmr, qms, qmg, t, &
rew, rei, rer, res, reg)
-! qcw, qci, qcr, qcs, qcg, rew, rei, rer, res, reg)
implicit none
- integer, intent (in) :: is, ie, js, je
+ integer, intent (in) :: is, ie, ks, ke
+ integer, intent (in), dimension (is:ie) :: lsm ! land sea mask, 0: ocean, 1: land, 2: sea ice
- real, intent (in), dimension (is:ie, js:je) :: den, t
- real, intent (in), dimension (is:ie, js:je) :: qw, qi, qr, qs, qg ! units: kg / kg
+ real, intent (in), dimension (is:ie, ks:ke) :: den, delp, t
+ real, intent (in), dimension (is:ie, ks:ke) :: qmw, qmi, qmr, qms, qmg !< units: kg / kg
-! real, intent (out), dimension (is:ie, js:je) :: qcw, qci, qcr, qcs, qcg ! units: kg / m^3
- real, dimension (is:ie, js:je) :: qcw, qci, qcr, qcs, qcg ! units: kg / m^3
- real, intent (out), dimension (is:ie, js:je) :: rew, rei, rer, res, reg ! units: micron
+ real, intent (out), dimension (is:ie, ks:ke) :: rew, rei, rer, res, reg !< units: micron
- integer :: i, j
+ real, dimension (is:ie, ks:ke) :: qcw, qci, qcr, qcs, qcg !< units: g / m^2
+
+ integer :: i, k
real :: lambdar, lambdas, lambdag
+ real :: dpg, rei_fac, mask, ccn, bw
+ real, parameter :: rho_0 = 50.e-3
real :: rhow = 1.0e3, rhor = 1.0e3, rhos = 1.0e2, rhog = 4.0e2
real :: n0r = 8.0e6, n0s = 3.0e6, n0g = 4.0e6
real :: alphar = 0.8, alphas = 0.25, alphag = 0.5
real :: gammar = 17.837789, gammas = 8.2850630, gammag = 11.631769
-! real :: qmin = 1.0e-5, ccn = 1.0e8, beta = 1.22
-! real :: qmin = 5.0e-6, ccn = 1.0e8, beta = 1.22
- real :: qmin = 9.0e-6, ccn = 1.0e8, beta = 1.22
-! real :: qmin = 1.0e-6, ccn = 1.0e8, beta = 1.22
-! real :: qmin = 1.0e-8, ccn = 1.0e8, beta = 1.22
-! real :: qmin = 1.0e-12, ccn = 1.0e8, beta = 1.22
-
- ! real :: rewmin = 1.0, rewmax = 25.0
- ! real :: reimin = 10.0, reimax = 300.0
- ! real :: rermin = 25.0, rermax = 225.0
- ! real :: resmin = 300, resmax = 1000.0
- ! real :: regmin = 1000.0, regmax = 1.0e5
- real :: rewmin = 5.0, rewmax = 10.0
- real :: reimin = 10.0, reimax = 150.0
-! real :: rermin = 0.0, rermax = 10000.0
-! real :: resmin = 0.0, resmax = 10000.0
-! real :: regmin = 0.0, regmax = 10000.0
- real :: rermin = 50.0, rermax = 10000.0
- real :: resmin = 100.0, resmax = 10000.0
- real :: regmin = 300.0, regmax = 10000.0
+ real :: qmin = 1.0e-12, beta = 1.22
- do j = js, je
+ do k = ks, ke
do i = is, ie
+
+ dpg = abs (delp (i, k)) / grav
+ mask = min (max (real(lsm (i)), 0.0), 2.0)
! -----------------------------------------------------------------------
- ! cloud water (martin et al., 1994)
+ ! cloud water (Martin et al., 1994)
! -----------------------------------------------------------------------
- if (qw (i, j) .gt. qmin) then
- qcw (i, j) = den (i, j) * qw (i, j)
- rew (i, j) = exp (1.0 / 3.0 * log ((3 * qcw (i, j)) / (4 * pi * rhow * ccn))) * 1.0e6
- rew (i, j) = max (rewmin, min (rewmax, rew (i, j)))
+ ccn = 0.80 * (- 1.15e-3 * (ccn_o ** 2) + 0.963 * ccn_o + 5.30) * abs (mask - 1.0) + &
+ 0.67 * (- 2.10e-4 * (ccn_l ** 2) + 0.568 * ccn_l - 27.9) * (1.0 - abs (mask - 1.0))
+
+ if (qmw (i, k) .gt. qmin) then
+ qcw (i, k) = dpg * qmw (i, k) * 1.0e3
+ rew (i, k) = exp (1.0 / 3.0 * log ((3.0 * den (i, k) * qmw (i, k)) / (4.0 * pi * rhow * ccn))) * 1.0e4
+ rew (i, k) = max (rewmin, min (rewmax, rew (i, k)))
else
- qcw (i, j) = 0.0
- rew (i, j) = rewmin
+ qcw (i, k) = 0.0
+ rew (i, k) = rewmin
endif
+
+ if (reiflag .eq. 1) then
! -----------------------------------------------------------------------
- ! cloud ice (heymsfield and mcfarquhar, 1996)
+ ! cloud ice (Heymsfield and Mcfarquhar, 1996)
! -----------------------------------------------------------------------
- if (qi (i, j) .gt. qmin) then
- qci (i, j) = den (i, j) * qi (i, j)
- if (t (i, j) - tice .lt. - 50) then
- rei (i, j) = beta / 9.917 * exp ((1 - 0.891) * log (1.0e3 * qci (i, j))) * 1.0e3
- elseif (t (i, j) - tice .lt. - 40) then
- rei (i, j) = beta / 9.337 * exp ((1 - 0.920) * log (1.0e3 * qci (i, j))) * 1.0e3
- elseif (t (i, j) - tice .lt. - 30) then
- rei (i, j) = beta / 9.208 * exp ((1 - 0.945) * log (1.0e3 * qci (i, j))) * 1.0e3
+ if (qmi (i, k) .gt. qmin) then
+ qci (i, k) = dpg * qmi (i, k) * 1.0e3
+ rei_fac = log (1.0e3 * qmi (i, k) * den (i, k))
+ if (t (i, k) - tice .lt. - 50) then
+ rei (i, k) = beta / 9.917 * exp (0.109 * rei_fac) * 1.0e3
+ elseif (t (i, k) - tice .lt. - 40) then
+ rei (i, k) = beta / 9.337 * exp (0.080 * rei_fac) * 1.0e3
+ elseif (t (i, k) - tice .lt. - 30) then
+ rei (i, k) = beta / 9.208 * exp (0.055 * rei_fac) * 1.0e3
else
- rei (i, j) = beta / 9.387 * exp ((1 - 0.969) * log (1.0e3 * qci (i, j))) * 1.0e3
+ rei (i, k) = beta / 9.387 * exp (0.031 * rei_fac) * 1.0e3
endif
- rei (i, j) = max (reimin, min (reimax, rei (i, j)))
+ rei (i, k) = max (reimin, min (reimax, rei (i, k)))
else
- qci (i, j) = 0.0
- rei (i, j) = reimin
+ qci (i, k) = 0.0
+ rei (i, k) = reimin
endif
+ endif
+
+ if (reiflag .eq. 2) then
+
! -----------------------------------------------------------------------
- ! rain (lin et al., 1983)
+ ! cloud ice (Wyser, 1998)
! -----------------------------------------------------------------------
- if (qr (i, j) .gt. qmin) then
- qcr (i, j) = den (i, j) * qr (i, j)
- lambdar = exp (0.25 * log (pi * rhor * n0r / qcr (i, j)))
- rer (i, j) = 0.5 * exp (log (gammar / 6) / alphar) / lambdar * 1.0e6
- rer (i, j) = max (rermin, min (rermax, rer (i, j)))
+ if (qmi (i, k) .gt. qmin) then
+ qci (i, k) = dpg * qmi (i, k) * 1.0e3
+ bw = - 2. + 1.e-3 * log10 (den (i, k) * qmi (i, k) / rho_0) * max (0.0, tice - t (i, k)) ** 1.5
+ rei (i, k) = 377.4 + bw * (203.3 + bw * (37.91 + 2.3696 * bw))
+ rei (i, k) = max (reimin, min (reimax, rei (i, k)))
else
- qcr (i, j) = 0.0
- rer (i, j) = rermin
+ qci (i, k) = 0.0
+ rei (i, k) = reimin
+ endif
+
endif
! -----------------------------------------------------------------------
- ! snow (lin et al., 1983)
+ ! rain (Lin et al., 1983)
! -----------------------------------------------------------------------
- if (qs (i, j) .gt. qmin) then
- qcs (i, j) = den (i, j) * qs (i, j)
- lambdas = exp (0.25 * log (pi * rhos * n0s / qcs (i, j)))
- res (i, j) = 0.5 * exp (log (gammas / 6) / alphas) / lambdas * 1.0e6
- res (i, j) = max (resmin, min (resmax, res (i, j)))
+ if (qmr (i, k) .gt. qmin) then
+ qcr (i, k) = dpg * qmr (i, k) * 1.0e3
+ lambdar = exp (0.25 * log (pi * rhor * n0r / qmr (i, k) / den (i, k)))
+ rer (i, k) = 0.5 * exp (log (gammar / 6) / alphar) / lambdar * 1.0e6
+ rer (i, k) = max (rermin, min (rermax, rer (i, k)))
else
- qcs (i, j) = 0.0
- res (i, j) = resmin
+ qcr (i, k) = 0.0
+ rer (i, k) = rermin
endif
! -----------------------------------------------------------------------
- ! graupel (lin et al., 1983)
+ ! snow (Lin et al., 1983)
! -----------------------------------------------------------------------
- if (qg (i, j) .gt. qmin) then
- qcg (i, j) = den (i, j) * qg (i, j)
- lambdag = exp (0.25 * log (pi * rhog * n0g / qcg (i, j)))
- reg (i, j) = 0.5 * exp (log (gammag / 6) / alphag) / lambdag * 1.0e6
- reg (i, j) = max (regmin, min (regmax, reg (i, j)))
+ if (qms (i, k) .gt. qmin) then
+ qcs (i, k) = dpg * qms (i, k) * 1.0e3
+ lambdas = exp (0.25 * log (pi * rhos * n0s / qms (i, k) / den (i, k)))
+ res (i, k) = 0.5 * exp (log (gammas / 6) / alphas) / lambdas * 1.0e6
+ res (i, k) = max (resmin, min (resmax, res (i, k)))
+ else
+ qcs (i, k) = 0.0
+ res (i, k) = resmin
+ endif
+
+ ! -----------------------------------------------------------------------
+ ! graupel (Lin et al., 1983)
+ ! -----------------------------------------------------------------------
+
+ if (qmg (i, k) .gt. qmin) then
+ qcg (i, k) = dpg * qmg (i, k) * 1.0e3
+ lambdag = exp (0.25 * log (pi * rhog * n0g / qmg (i, k) / den (i, k)))
+ reg (i, k) = 0.5 * exp (log (gammag / 6) / alphag) / lambdag * 1.0e6
+ reg (i, k) = max (regmin, min (regmax, reg (i, k)))
else
- qcg (i, j) = 0.0
- reg (i, j) = regmin
+ qcg (i, k) = 0.0
+ reg (i, k) = regmin
endif
enddo
diff --git a/physics/module_nst_water_prop.f90 b/physics/module_nst_water_prop.f90
index 20c4dff88..3f3916396 100644
--- a/physics/module_nst_water_prop.f90
+++ b/physics/module_nst_water_prop.f90
@@ -657,7 +657,8 @@ subroutine get_dtzm_point(xt,xz,dt_cool,zc,z1,z2,dtm)
end subroutine get_dtzm_point
!>\ingroup waterprop
- subroutine get_dtzm_2d(xt,xz,dt_cool,zc,wet,icy,z1,z2,nx,ny,dtm)
+ subroutine get_dtzm_2d(xt,xz,dt_cool,zc,wet,z1,z2,nx,ny,dtm)
+!subroutine get_dtzm_2d(xt,xz,dt_cool,zc,wet,icy,z1,z2,nx,ny,dtm)
! ===================================================================== !
! !
! description: get dtm = mean of dT(z) (z1 - z2) with NSST dT(z) !
@@ -695,7 +696,8 @@ subroutine get_dtzm_2d(xt,xz,dt_cool,zc,wet,icy,z1,z2,nx,ny,dtm)
integer, intent(in) :: nx,ny
real (kind=kind_phys), dimension(nx,ny), intent(in) :: xt,xz,dt_cool,zc
- logical, dimension(nx,ny), intent(in) :: wet,icy
+ logical, dimension(nx,ny), intent(in) :: wet
+! logical, dimension(nx,ny), intent(in) :: wet,icy
real (kind=kind_phys), intent(in) :: z1,z2
real (kind=kind_phys), dimension(nx,ny), intent(out) :: dtm
! Local variables
@@ -712,7 +714,8 @@ subroutine get_dtzm_2d(xt,xz,dt_cool,zc,wet,icy,z1,z2,nx,ny,dtm)
!
dtw(i,j) = 0.0
dtc(i,j) = 0.0
- if ( wet(i,j) .and. .not.icy(i,j) ) then
+! if ( wet(i,j) .and. .not.icy(i,j) ) then
+ if ( wet(i,j) ) then
!
! get the mean warming in the range of z=z1 to z=z2
!
@@ -746,16 +749,18 @@ subroutine get_dtzm_2d(xt,xz,dt_cool,zc,wet,icy,z1,z2,nx,ny,dtm)
endif
endif
endif
- endif ! if wet(i,j) .and. .not.icy(i,j)
+ endif ! if ( wet(i,j) .and. .not.icy(i,j) ) then
enddo
enddo
!
! get the mean T departure from Tf in the range of z=z1 to z=z2
+! DH* NEED NTHREADS HERE! TODO
!$omp parallel do private(j,i)
do j = 1, ny
do i= 1, nx
- if ( wet(i,j) .and. .not.icy(i,j)) then
+! if ( wet(i,j) .and. .not.icy(i,j)) then
+ if ( wet(i,j) ) then
dtm(i,j) = dtw(i,j) - dtc(i,j)
endif
enddo
diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90
index 4f1f7dbad..af7a8362e 100755
--- a/physics/module_sf_noahmplsm.f90
+++ b/physics/module_sf_noahmplsm.f90
@@ -286,6 +286,7 @@ subroutine noahmp_sflx (parameters, &
qc , soldn , lwdn , & ! in : forcing
prcpconv, prcpnonc, prcpshcv, prcpsnow, prcpgrpl, prcphail, & ! in : forcing
tbot , co2air , o2air , foln , ficeold , zlvl , & ! in : forcing
+ lheatstrg , & ! in : canopy heat storage
albold , sneqvo , & ! in/out :
stc , sh2o , smc , tah , eah , fwet , & ! in/out :
canliq , canice , tv , tg , qsfc , qsnow , & ! in/out :
@@ -293,9 +294,9 @@ subroutine noahmp_sflx (parameters, &
zwt , wa , wt , wslake , lfmass , rtmass , & ! in/out :
stmass , wood , stblcp , fastcp , lai , sai , & ! in/out :
cm , ch , tauss , & ! in/out :
- smcwtd ,deeprech , rech , & ! in/out :
+ smcwtd ,deeprech , rech , cpfac , & ! in/out :
z0wrf , &
- fsa , fsr , fira , fsh , ssoil , fcev , & ! out :
+ fsa , fsr , fira , fshx , ssoil , fcev , & ! out :
fgev , fctr , ecan , etran , edir , trad , & ! out :
tgb , tgv , t2mv , t2mb , q2v , q2b , & ! out :
runsrf , runsub , apar , psn , sav , sag , & ! out :
@@ -336,6 +337,7 @@ subroutine noahmp_sflx (parameters, &
real , intent(in) :: lwdn !downward longwave radiation (w/m2)
real , intent(in) :: sfcprs !pressure (pa)
real , intent(inout) :: zlvl !reference height (m)
+ logical , intent(in) :: lheatstrg ! flag for canopy heat storage parameterization
real , intent(in) :: cosz !cosine solar zenith angle [0-1]
real , intent(in) :: tbot !bottom condition for soil temp. [k]
real , intent(in) :: foln !foliage nitrogen (%) [1-saturated]
@@ -394,13 +396,14 @@ subroutine noahmp_sflx (parameters, &
real, intent(inout) :: smcwtd !soil water content between bottom of the soil and water table [m3/m3]
real, intent(inout) :: deeprech !recharge to or from the water table when deep [m]
real, intent(inout) :: rech !recharge to or from the water table when shallow [m] (diagnostic)
+ real, intent(inout) :: cpfac ! heat capacity enhancement factor due to heat storage
! output
real , intent(out) :: z0wrf !combined z0 sent to coupled model
real , intent(out) :: fsa !total absorbed solar radiation (w/m2)
real , intent(out) :: fsr !total reflected solar radiation (w/m2)
real , intent(out) :: fira !total net lw rad (w/m2) [+ to atm]
- real , intent(out) :: fsh !total sensible heat (w/m2) [+ to atm]
+ real , intent(out) :: fshx !total sensible heat (w/m2) [+ to atm]
real , intent(out) :: fcev !canopy evap heat (w/m2) [+ to atm]
real , intent(out) :: fgev !ground evap heat (w/m2) [+ to atm]
real , intent(out) :: fctr !transpiration heat (w/m2) [+ to atm]
@@ -450,6 +453,7 @@ subroutine noahmp_sflx (parameters, &
real :: taux !wind stress: e-w (n/m2)
real :: tauy !wind stress: n-s (n/m2)
real :: rhoair !density air (kg/m3)
+ real :: fsh !total sensible heat (w/m2) [+ to atm]
! real, dimension( 1: 5) :: vocflx !voc fluxes [ug c m-2 h-1]
real, dimension(-nsnow+1:nsoil) :: dzsnso !snow/soil layer thickness [m]
real :: thair !potential temperature (k)
@@ -640,6 +644,7 @@ subroutine noahmp_sflx (parameters, &
call energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in
isnow ,dt ,rhoair ,sfcprs ,qair , & !in
sfctmp ,thair ,lwdn ,uu ,vv ,zlvl , & !in
+ lheatstrg , & !in
co2air ,o2air ,solad ,solai ,cosz ,igs , & !in
eair ,tbot ,zsnso ,zsoil , & !in
elai ,esai ,fwet ,foln , & !in
@@ -648,16 +653,16 @@ subroutine noahmp_sflx (parameters, &
z0wrf , &
imelt ,snicev ,snliqv ,epore ,t2m ,fsno , & !out
sav ,sag ,qmelt ,fsa ,fsr ,taux , & !out
- tauy ,fira ,fsh ,fcev ,fgev ,fctr , & !out
+ tauy ,fira ,fsh ,fshx ,fcev ,fgev ,fctr , & !out
trad ,psn ,apar ,ssoil ,btrani ,btran , & !out
ponding,ts ,latheav , latheag , frozen_canopy,frozen_ground, & !out
tv ,tg ,stc ,snowh ,eah ,tah , & !inout
sneqvo ,sneqv ,sh2o ,smc ,snice ,snliq , & !inout
albold ,cm ,ch ,dx ,dz8w ,q2 , & !inout
#ifdef CCPP
- tauss ,errmsg ,errflg , & !inout
+ tauss ,cpfac ,errmsg ,errflg , & !inout
#else
- tauss , & !inout
+ tauss ,cpfac , & !inout
#endif
!jref:start
qc ,qsfc ,psfc , & !in
@@ -709,7 +714,7 @@ subroutine noahmp_sflx (parameters, &
! water and energy balance check
- call error (parameters,swdown ,fsa ,fsr ,fira ,fsh ,fcev , & !in
+ call error (parameters,swdown ,fsa ,fsr ,fira ,fsh ,fcev , & !in
fgev ,fctr ,ssoil ,beg_wb ,canliq ,canice , & !in
sneqv ,wa ,smc ,dzsnso ,prcp ,ecan , & !in
etran ,edir ,runsrf ,runsub ,dt ,nsoil , & !in
@@ -1413,6 +1418,7 @@ end subroutine error
subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in
isnow ,dt ,rhoair ,sfcprs ,qair , & !in
sfctmp ,thair ,lwdn ,uu ,vv ,zref , & !in
+ lheatstrg , & !in
co2air ,o2air ,solad ,solai ,cosz ,igs , & !in
eair ,tbot ,zsnso ,zsoil , & !in
elai ,esai ,fwet ,foln , & !in
@@ -1421,16 +1427,16 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in
z0wrf , &
imelt ,snicev ,snliqv ,epore ,t2m ,fsno , & !out
sav ,sag ,qmelt ,fsa ,fsr ,taux , & !out
- tauy ,fira ,fsh ,fcev ,fgev ,fctr , & !out
+ tauy ,fira ,fsh ,fshx ,fcev ,fgev ,fctr , & !out
trad ,psn ,apar ,ssoil ,btrani ,btran , & !out
ponding,ts ,latheav , latheag , frozen_canopy,frozen_ground, & !out
tv ,tg ,stc ,snowh ,eah ,tah , & !inout
sneqvo ,sneqv ,sh2o ,smc ,snice ,snliq , & !inout
albold ,cm ,ch ,dx ,dz8w ,q2 , & !inout
#ifdef CCPP
- tauss ,errmsg ,errflg, & !inout
+ tauss ,cpfac ,errmsg ,errflg, & !inout
#else
- tauss , & !inout
+ tauss ,cpfac , & !inout
#endif
!jref:start
qc ,qsfc ,psfc , & !in
@@ -1512,6 +1518,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in
real , intent(in) :: igs !growing season index (0=off, 1=on)
real , intent(in) :: zref !reference height (m)
+ logical , intent(in) :: lheatstrg ! flag for canopy heat storage parameterization
real , intent(in) :: tbot !bottom condition for soil temp. (k)
real , dimension(-nsnow+1:nsoil), intent(in) :: zsnso !layer-bottom depth from snow surf [m]
real , dimension( 1:nsoil), intent(in) :: zsoil !layer-bottom depth from soil surf [m]
@@ -1546,6 +1553,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in
real , intent(out) :: tauy !wind stress: n-s (n/m2)
real , intent(out) :: fira !total net lw. rad (w/m2) [+ to atm]
real , intent(out) :: fsh !total sensible heat (w/m2) [+ to atm]
+ real , intent(out) :: fshx !total sensible heat (w/m2) [+ to atm]
real , intent(out) :: fcev !canopy evaporation (w/m2) [+ to atm]
real , intent(out) :: fgev !ground evaporation (w/m2) [+ to atm]
real , intent(out) :: fctr !transpiration (w/m2) [+ to atm]
@@ -1592,6 +1600,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in
real , intent(inout) :: tah !canopy air temperature (k)
real , intent(inout) :: albold !snow albedo at last time step(class type)
real , intent(inout) :: tauss !non-dimensional snow age
+ real , intent(inout) :: cpfac !heat capacity enhancement factor due to heat storage
real , intent(inout) :: cm !momentum drag coefficient
real , intent(inout) :: ch !sensible heat exchange coefficient
real , intent(inout) :: q1
@@ -1693,6 +1702,11 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in
real, parameter :: mpe = 1.e-6
real, parameter :: psiwlt = -150. !metric potential for wilting point (m)
real, parameter :: z0 = 0.01 ! bare-soil roughness length (m) (i.e., under the canopy)
+!
+! parameters for heat storage parametrization
+!
+ real, parameter :: z0min = 0.2 !minimum roughness length for heat storage
+ real, parameter :: z0max = 1.0 !maximum roughness length for heat storage
! ---------------------------------------------------------------------------------------------------
! initialize fluxes from veg. fraction
@@ -1758,6 +1772,13 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in
z0m = z0mg
zpd = zpdg
end if
+!
+! compute heat capacity enhancement factor as a function of z0m to mimic heat storage
+!
+ if (lheatstrg .and. (.not. parameters%urban_flag) ) then
+ cpfac = (z0m - z0min) / (z0max - z0min)
+ cpfac = 1. + min(max(cpfac, 0.0), 1.0)
+ endif
zlvl = max(zpd,parameters%hvt) + zref
if(zpdg >= zlvl) zlvl = zpdg + zref
@@ -1862,7 +1883,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in
latheav = hsub
frozen_canopy = .true.
end if
- gammav = cpair*sfcprs/(0.622*latheav)
+ gammav = cpair*cpfac*sfcprs/(0.622*latheav)
if (tg .gt. tfrz) then
latheag = hvap
@@ -1871,14 +1892,14 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in
latheag = hsub
frozen_ground = .true.
end if
- gammag = cpair*sfcprs/(0.622*latheag)
+ gammag = cpair*cpfac*sfcprs/(0.622*latheag)
! if (sfctmp .gt. tfrz) then
! lathea = hvap
! else
! lathea = hsub
! end if
-! gamma = cpair*sfcprs/(0.622*lathea)
+! gamma = cpair*cpfac*sfcprs/(0.622*lathea)
! surface temperatures of the ground and canopy and energy fluxes
@@ -1891,9 +1912,9 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in
call vege_flux (parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & !in
dt ,sav ,sag ,lwdn ,ur , & !in
uu ,vv ,sfctmp ,thair ,qair , & !in
- eair ,rhoair ,snowh ,vai ,gammav ,gammag , & !in
+ eair ,rhoair ,snowh ,vai ,gammav ,gammag , & !in
fwet ,laisun ,laisha ,cwp ,dzsnso , & !in
- zlvl ,zpd ,z0m ,fveg , & !in
+ zlvl ,cpfac ,zpd ,z0m ,fveg , & !in
z0mg ,emv ,emg ,canliq ,fsno, & !in
canice ,stc ,df ,rssun ,rssha , & !in
rsurf ,latheav ,latheag ,parsun ,parsha ,igs , & !in
@@ -1923,7 +1944,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in
call bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & !in
lwdn ,ur ,uu ,vv ,sfctmp , & !in
thair ,qair ,eair ,rhoair ,snowh , & !in
- dzsnso ,zlvl ,zpdg ,z0mg ,fsno, & !in
+ dzsnso ,zlvl ,zpdg ,z0mg ,fsno, & !in
emg ,stc ,df ,rsurf ,latheag , & !in
gammag ,rhsur ,iloc ,jloc ,q2 ,pahb , & !in
#ifdef CCPP
@@ -1949,6 +1970,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in
tauy = fveg * tauyv + (1.0 - fveg) * tauyb
fira = fveg * irg + (1.0 - fveg) * irb + irc
fsh = fveg * shg + (1.0 - fveg) * shb + shc
+ fshx = fveg * shg/cpfac + (1.0 - fveg) * shb + shc/cpfac
fgev = fveg * evg + (1.0 - fveg) * evb
ssoil = fveg * ghv + (1.0 - fveg) * ghb
fcev = evc
@@ -1967,6 +1989,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in
tauy = tauyb
fira = irb
fsh = shb
+ fshx = shb
fgev = evb
ssoil = ghb
tg = tgb
@@ -3260,7 +3283,8 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , &
uu ,vv ,sfctmp ,thair ,qair , & !in
eair ,rhoair ,snowh ,vai ,gammav ,gammag, & !in
fwet ,laisun ,laisha ,cwp ,dzsnso , & !in
- zlvl ,zpd ,z0m ,fveg , & !in
+ zlvl ,cpfac , & !in
+ zpd ,z0m ,fveg , & !in
z0mg ,emv ,emg ,canliq ,fsno, & !in
canice ,stc ,df ,rssun ,rssha , & !in
rsurf ,latheav ,latheag ,parsun ,parsha ,igs , & !in
@@ -3320,6 +3344,8 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , &
real, intent(in) :: laisun !sunlit leaf area index, one-sided (m2/m2)
real, intent(in) :: laisha !shaded leaf area index, one-sided (m2/m2)
real, intent(in) :: zlvl !reference height (m)
+ real, intent(in) :: cpfac !heat capacity enhancement factor due to heat storage
+
real, intent(in) :: zpd !zero plane displacement (m)
real, intent(in) :: z0m !roughness length, momentum (m)
real, intent(in) :: z0mg !roughness length, momentum, ground (m)
@@ -3449,6 +3475,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , &
real :: kh !turbulent transfer coefficient, sensible heat, (m2/s)
real :: h !temporary sensible heat flux (w/m2)
real :: hg !temporary sensible heat flux (w/m2)
+
real :: moz !monin-obukhov stability parameter
real :: mozg !monin-obukhov stability parameter
real :: mozold !monin-obukhov stability parameter from prior iteration
@@ -3578,6 +3605,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , &
air = -emv*(1.+(1.-emv)*(1.-emg))*lwdn - emv*emg*sb*tg**4
cir = (2.-emv*(1.-emg))*emv*sb
+
! ---------------------------------------------------------------------------------------------
loop1: do iter = 1, niterc ! begin stability iteration
@@ -3674,7 +3702,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , &
cond = cah + cvh + cgh
ata = (sfctmp*cah + tg*cgh) / cond
bta = cvh/cond
- csh = (1.-bta)*rhoair*cpair*cvh
+ csh = (1.-bta)*rhoair*cpair*cpfac*cvh
! prepare for latent heat flux above veg.
@@ -3685,8 +3713,8 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , &
cond = caw + cew + ctw + cgw
aea = (eair*caw + estg*cgw) / cond
bea = (cew+ctw)/cond
- cev = (1.-bea)*cew*rhoair*cpair/gammav ! barlage: change to vegetation v3.6
- ctr = (1.-bea)*ctw*rhoair*cpair/gammav
+ cev = (1.-bea)*cew*rhoair*cpair*cpfac/gammav ! barlage: change to vegetation v3.6
+ ctr = (1.-bea)*ctw*rhoair*cpair*cpfac/gammav
! evaluate surface fluxes with current temperature and solve for dts
@@ -3694,9 +3722,9 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , &
eah = aea + bea*estv ! canopy air e
irc = fveg*(air + cir*tv**4)
- shc = fveg*rhoair*cpair*cvh * ( tv-tah)
- evc = fveg*rhoair*cpair*cew * (estv-eah) / gammav ! barlage: change to v in v3.6
- tr = fveg*rhoair*cpair*ctw * (estv-eah) / gammav
+ shc = fveg*rhoair*cpair*cpfac*cvh * ( tv-tah)
+ evc = fveg*rhoair*cpair*cpfac*cew * (estv-eah) / gammav ! barlage: change to v in v3.6
+ tr = fveg*rhoair*cpair*cpfac*ctw * (estv-eah) / gammav
if (tv > tfrz) then
evc = min(canliq*latheav/dt,evc) ! barlage: add if block for canice in v3.6
else
@@ -3736,8 +3764,8 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , &
air = - emg*(1.-emv)*lwdn - emg*emv*sb*tv**4
cir = emg*sb
- csh = rhoair*cpair/rahg
- cev = rhoair*cpair / (gammag*(rawg+rsurf)) ! barlage: change to ground v3.6
+ csh = rhoair*cpair*cpfac/rahg
+ cev = rhoair*cpair*cpfac / (gammag*(rawg+rsurf)) ! barlage: change to ground v3.6
cgh = 2.*df(isnow+1)/dzsnso(isnow+1)
! write(*,*)'inside tg=',tg,'stc(1)=',stc(1)
@@ -3792,10 +3820,10 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , &
! consistent vegetation air temperature and vapor pressure since tg is not consistent with the tah/eah
! calculation.
-! tah = sfctmp + (shg+shc)/(rhoair*cpair*cah)
-! tah = sfctmp + (shg*fveg+shc)/(rhoair*cpair*cah) ! ground flux need fveg
-! eah = eair + (evc+fveg*(tr+evg))/(rhoair*caw*cpair/gammag )
-! qfx = (qsfc-qair)*rhoair*caw !*cpair/gammag
+! tah = sfctmp + (shg+shc)/(rhoair*cpair*cpfac*cah)
+! tah = sfctmp + (shg*fveg+shc)/(rhoair*cpair*cpfac*cah) ! ground flux need fveg
+! eah = eair + (evc+fveg*(tr+evg))/(rhoair*caw*cpair*cpfac/gammag )
+! qfx = (qsfc-qair)*rhoair*cpfac*caw !*cpair/gammag
! 2m temperature over vegetation ( corrected for low cq2v values )
if (opt_sfc == 1 .or. opt_sfc == 2) then
@@ -3808,7 +3836,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , &
! q2v = (eah*0.622/(sfcprs - 0.378*eah))
q2v = qsfc
else
- t2mv = tah - (shg+shc/fveg)/(rhoair*cpair) * 1./cah2
+ t2mv = tah - (shg+shc/fveg)/(rhoair*cpair*cpfac) * 1./cah2
! q2v = (eah*0.622/(sfcprs - 0.378*eah))- qfx/(rhoair*fv)* 1./vkc * log((2.+z0h)/z0h)
q2v = qsfc - ((evc+tr)/fveg+evg)/(latheav*rhoair) * 1./cq2v
endif
diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90
index 3b2da9c3e..812229f98 100644
--- a/physics/mp_thompson.F90
+++ b/physics/mp_thompson.F90
@@ -395,7 +395,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, &
graupel = max(0.0, delta_graupel_mp/1000.0_kind_phys)
ice = max(0.0, delta_ice_mp/1000.0_kind_phys)
snow = max(0.0, delta_snow_mp/1000.0_kind_phys)
- rain = max(0.0, delta_rain_mp - (delta_graupel_mp + delta_ice_mp + delta_snow_mp)/1000.0_kind_phys)
+ rain = max(0.0, (delta_rain_mp - (delta_graupel_mp + delta_ice_mp + delta_snow_mp))/1000.0_kind_phys)
end subroutine mp_thompson_run
!>@}
diff --git a/physics/rrtmg_lw_pre.F90 b/physics/rrtmg_lw_pre.F90
index 783d65e90..5f128a79a 100644
--- a/physics/rrtmg_lw_pre.F90
+++ b/physics/rrtmg_lw_pre.F90
@@ -43,7 +43,7 @@ subroutine rrtmg_lw_pre_run (Model, Grid, Sfcprop, Radtend, im, tsfg, tsfa, errm
!! emissivity for LW radiation.
call setemis (Grid%xlon, Grid%xlat, Sfcprop%slmsk, & ! --- inputs
Sfcprop%snowd, Sfcprop%sncovr, Sfcprop%zorl, &
- tsfg, tsfa, Sfcprop%hprim, IM, &
+ tsfg, tsfa, Sfcprop%hprime(:,1), IM, &
Radtend%semis) ! --- outputs
endif
diff --git a/physics/rrtmg_sw_pre.F90 b/physics/rrtmg_sw_pre.F90
index de994ba79..8eeb16430 100644
--- a/physics/rrtmg_sw_pre.F90
+++ b/physics/rrtmg_sw_pre.F90
@@ -66,13 +66,13 @@ subroutine rrtmg_sw_pre_run (Model, Grid, Sfcprop, Radtend, im, &
!> - Call module_radiation_surface::setalb() to setup surface albedo.
!! for SW radiation.
- call setalb (Sfcprop%slmsk, Sfcprop%snowd, Sfcprop%sncovr,& ! --- inputs:
- Sfcprop%snoalb, Sfcprop%zorl, Radtend%coszen,&
- tsfg, tsfa, Sfcprop%hprim, Sfcprop%alvsf, &
- Sfcprop%alnsf, Sfcprop%alvwf, Sfcprop%alnwf, &
- Sfcprop%facsf, Sfcprop%facwf, Sfcprop%fice, &
- Sfcprop%tisfc, IM, &
- alb1d, Model%pertalb, & ! mg, sfc-perts
+ call setalb (Sfcprop%slmsk, Sfcprop%snowd, Sfcprop%sncovr, & ! --- inputs:
+ Sfcprop%snoalb, Sfcprop%zorl, Radtend%coszen, &
+ tsfg, tsfa, Sfcprop%hprime(:,1), Sfcprop%alvsf, &
+ Sfcprop%alnsf, Sfcprop%alvwf, Sfcprop%alnwf, &
+ Sfcprop%facsf, Sfcprop%facwf, Sfcprop%fice, &
+ Sfcprop%tisfc, IM, &
+ alb1d, Model%pertalb, & ! mg, sfc-perts
sfcalb) ! --- outputs
!> -# Approximate mean surface albedo from vis- and nir- diffuse values.
diff --git a/physics/set_soilveg.f b/physics/set_soilveg.f
index 60a6395b8..efef0f24b 100644
--- a/physics/set_soilveg.f
+++ b/physics/set_soilveg.f
@@ -136,8 +136,9 @@ subroutine set_soilveg(me,isot,ivet,nlunit)
! ----------------------------------------------------------------------
defined_veg=20
- NROOT_DATA =(/4,4,4,4,4,3,3,3,3,3,3,3,1,3,2,
- & 3,0,3,3,2,0,0,0,0,0,0,0,0,0,0/)
+ NROOT_DATA =(/4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 1, 3, 2,
+ & 3, 1, 3, 3, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/)
+! & 3, 0, 3, 3, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/) ! Moorthi
! ----------------------------------------------------------------------
! VEGETATION CLASS-RELATED ARRAYS
! ----------------------------------------------------------------------
diff --git a/physics/sfc_cice.f b/physics/sfc_cice.f
index 562d00bee..0a1a49c77 100644
--- a/physics/sfc_cice.f
+++ b/physics/sfc_cice.f
@@ -35,15 +35,18 @@ end subroutine sfc_cice_finalize
!! @{
-!! use physcons, only : hvap => con_hvap, cp => con_cp, &
+!! use physcons, only : hvap => con_hvap, cp => con_cp, &
!! & rvrdm1 => con_fvirt, rd => con_rd
!
!-----------------------------------
subroutine sfc_cice_run &
- & ( im, cplflx, cplchm, hvap, cp, rvrdm1, rd, & ! --- inputs:
- & u1, v1, t1, q1, cm, ch, prsl1, prslki, &
- & flag_cice, ddvel, flag_iter, dqsfc, dtsfc, &
- & qsurf, cmm, chh, evap, hflx, & ! --- outputs:
+! --- inputs:
+ & ( im, cplflx, cplchm, hvap, cp, rvrdm1, rd, &
+ & t1, q1, cm, ch, prsl1, &
+ & wind, flag_cice, flag_iter, dqsfc, dtsfc, &
+ & dusfc, dvsfc, &
+! --- outputs:
+ & qsurf, cmm, chh, evap, hflx, stress, &
& errmsg, errflg
& )
@@ -55,40 +58,42 @@ subroutine sfc_cice_run &
! !
! call sfc_cice !
! inputs: !
-! ( im, u1, v1, t1, q1, cm, ch, prsl1, prslki, !
-! islimsk, ddvel, flag_iter, dqsfc, dtsfc, !
+! ( im, cplflx, cplchm, hvap, cp, rvrdm1, rd, !
+! t1, q1, cm, ch, prsl1, !
+! wind, flag_cice, flag_iter, dqsfc, dtsfc, !
+! dusfc, dvsfc, !
! outputs: !
-! qsurf, cmm, chh, evap, hflx) !
+! qsurf, cmm, chh, evap, hflx, stress) !
! !
! ==================== defination of variables ==================== !
! !
! inputs:
! im, - integer, horiz dimension
-! u1, v1 - real, u/v component of surface layer wind
+!! u1, v1 - real, u/v component of surface layer wind
! t1 - real, surface layer mean temperature ( k )
! q1 - real, surface layer mean specific humidity
! cm - real, surface exchange coeff for momentum (m/s)
! ch - real, surface exchange coeff heat & moisture(m/s)
! prsl1 - real, surface layer mean pressure
-! prslki - real, ?
-! islimsk - integer, sea/land/ice mask
-! ddvel - real, ?
+! wind - real, wind speed (m/s)
! flag_iter- logical
! dqsfc - real, latent heat flux
! dtsfc - real, sensible heat flux
+! dusfc - real, zonal momentum stress
+! dvsfc - real, meridional momentum stress
! outputs:
! qsurf - real, specific humidity at sfc
! cmm - real, ?
! chh - real, ?
! evap - real, evaperation from latent heat
! hflx - real, sensible heat
+! stress - real, surface stress
! ==================== end of description ===================== !
!
!
use machine , only : kind_phys
implicit none
-
real (kind=kind_phys), intent(in) :: hvap, cp, rvrdm1, rd
! --- inputs:
@@ -96,24 +101,22 @@ subroutine sfc_cice_run &
logical, intent(in) :: cplflx
logical, intent(in) :: cplchm
- real (kind=kind_phys), dimension(im), intent(in) :: u1, v1, &
- & t1, q1, cm, ch, prsl1, prslki, ddvel, dqsfc, dtsfc
+! real (kind=kind_phys), dimension(im), intent(in) :: u1, v1, &
+ real (kind=kind_phys), dimension(im), intent(in) :: &
+ & t1, q1, cm, ch, prsl1, wind, dqsfc, dtsfc, dusfc, dvsfc
- logical, dimension(im), intent(in) :: flag_cice
-
- logical, intent(in) :: flag_iter(im)
+ logical, intent(in) :: flag_cice(im), flag_iter(im)
! --- outputs:
real (kind=kind_phys), dimension(im), intent(out) :: qsurf, &
- & cmm, chh, evap, hflx
+ & cmm, chh, evap, hflx, stress
!
character(len=*), intent(out) :: errmsg
integer, intent(out) :: errflg
! --- locals:
- real (kind=kind_phys), dimension(im) :: q0, rch, rho, tv1, wind
- real (kind=kind_phys) :: tem
+ real (kind=kind_phys) :: rho, tem
real(kind=kind_phys) :: cpinv, hvapi, elocp
@@ -134,22 +137,17 @@ subroutine sfc_cice_run &
do i = 1, im
if (flag_cice(i) .and. flag_iter(i)) then
- wind(i) = sqrt(u1(i)*u1(i) + v1(i)*v1(i)) &
- & + max(0.0, min(ddvel(i), 30.0))
- wind(i) = max(wind(i), 1.0)
-
- q0(i) = max(q1(i), 1.0e-8)
- tv1(i) = t1(i) * (1.0 + rvrdm1*q0(i))
- rho(i) = prsl1(i) / (rd*tv1(i))
+ rho = prsl1(i) &
+ & / (rd * t1(i) * (1.0 + rvrdm1*max(q1(i), 1.0e-8)))
- cmm(i) = cm(i) * wind(i)
- chh(i) = rho(i) * ch(i) * wind(i)
- rch(i) = chh(i) * cp
+ cmm(i) = wind(i) * cm(i)
+ chh(i) = wind(i) * ch(i) * rho
- qsurf(i) = q1(i) + dqsfc(i) / (elocp*rch(i))
- tem = 1.0 / rho(i)
- hflx(i) = dtsfc(i) * tem * cpinv
- evap(i) = dqsfc(i) * tem * hvapi
+ qsurf(i) = q1(i) + dqsfc(i) / (hvap*chh(i))
+ tem = 1.0 / rho
+ hflx(i) = dtsfc(i) * tem * cpinv
+ evap(i) = dqsfc(i) * tem * hvapi
+ stress(i) = sqrt(dusfc(i)*dusfc(i) + dvsfc(i)*dvsfc(i)) * tem
endif
enddo
diff --git a/physics/sfc_cice.meta b/physics/sfc_cice.meta
index 756c760a4..48aa1f4c8 100644
--- a/physics/sfc_cice.meta
+++ b/physics/sfc_cice.meta
@@ -61,24 +61,6 @@
kind = kind_phys
intent = in
optional = F
-[u1]
- standard_name = x_wind_at_lowest_model_layer
- long_name = u component of surface layer wind
- units = m s-1
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = in
- optional = F
-[v1]
- standard_name = y_wind_at_lowest_model_layer
- long_name = v component of surface layer wind
- units = m s-1
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = in
- optional = F
[t1]
standard_name = air_temperature_at_lowest_model_layer
long_name = surface layer mean temperature
@@ -124,10 +106,10 @@
kind = kind_phys
intent = in
optional = F
-[prslki]
- standard_name = ratio_of_exner_function_between_midlayer_and_interface_at_lowest_model_layer
- long_name = Exner function ratio bt midlayer and interface at 1st layer
- units = ratio
+[wind]
+ standard_name = wind_speed_at_lowest_model_layer
+ long_name = wind speed at lowest model level
+ units = m s-1
dimensions = (horizontal_dimension)
type = real
kind = kind_phys
@@ -141,15 +123,6 @@
type = logical
intent = in
optional = F
-[ddvel]
- standard_name = surface_wind_enhancement_due_to_convection
- long_name = wind enhancement due to convection
- units = m s-1
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = in
- optional = F
[flag_iter]
standard_name = flag_for_iteration
long_name = flag for iteration
@@ -176,6 +149,24 @@
kind = kind_phys
intent = in
optional = F
+[dusfc]
+ standard_name = surface_x_momentum_flux_for_coupling_interstitial
+ long_name = sfc x momentum flux for coupling interstitial
+ units = Pa
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[dvsfc]
+ standard_name = surface_y_momentum_flux_for_coupling_interstitial
+ long_name = sfc y momentum flux for coupling interstitial
+ units = Pa
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
[qsurf]
standard_name = surface_specific_humidity_over_ice
long_name = surface air saturation specific humidity over ice
@@ -221,6 +212,15 @@
kind = kind_phys
intent = inout
optional = F
+[stress]
+ standard_name = surface_wind_stress_over_ice
+ long_name = surface wind stress over ice
+ units = m2 s-2
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = inout
+ optional = F
[errmsg]
standard_name = ccpp_error_message
long_name = error message for error handling in CCPP
diff --git a/physics/sfc_diag_post.F90 b/physics/sfc_diag_post.F90
index 6f5561cc0..767e98db5 100644
--- a/physics/sfc_diag_post.F90
+++ b/physics/sfc_diag_post.F90
@@ -41,7 +41,7 @@ subroutine sfc_diag_post_run (im, lsm, lsm_noahmp, dry, lssav, dtf, con_eps, con
! Initialize CCPP error handling variables
errmsg = ''
errflg = 0
-
+
if (lsm == lsm_noahmp) then
do i=1,im
if(dry(i)) then
@@ -50,7 +50,7 @@ subroutine sfc_diag_post_run (im, lsm, lsm_noahmp, dry, lssav, dtf, con_eps, con
endif
enddo
endif
-
+
if (lssav) then
do i=1,im
tmpmax(i) = max(tmpmax(i),t2m(i))
diff --git a/physics/sfc_diff.f b/physics/sfc_diff.f
index 5ada7288c..4cbf94245 100644
--- a/physics/sfc_diff.f
+++ b/physics/sfc_diff.f
@@ -61,8 +61,8 @@ end subroutine sfc_diff_finalize
!! - Calculate the exchange coefficients:\f$cm\f$, \f$ch\f$, and \f$stress\f$ as inputs of other \a sfc schemes.
!!
subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in)
- & ps,u1,v1,t1,q1,z1, & !intent(in)
- & prsl1,prslki,prsik1,prslk1,ddvel, & !intent(in)
+ & ps,t1,q1,z1,wind, & !intent(in)
+ & prsl1,prslki,prsik1,prslk1, & !intent(in)
& sigmaf,vegtype,shdmax,ivegsrc, & !intent(in)
& z0pert,ztpert, & ! mg, sfc-perts !intent(in)
& flag_iter,redrag, & !intent(in)
@@ -81,27 +81,23 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in)
& fh_ocn, fh_lnd, fh_ice, & !intent(inout)
& fm10_ocn, fm10_lnd, fm10_ice, & !intent(inout)
& fh2_ocn, fh2_lnd, fh2_ice, & !intent(inout)
- & wind , & !intent(inout)
& errmsg, errflg) !intent(out)
!
-! DH* 20190718: prslki can be removed if GSD_SURFACE_FLUXES_BUGFIX is adopted
- use funcphys, only : fpvs
-
implicit none
!
integer, intent(in) :: im, ivegsrc
integer, intent(in) :: sfc_z0_type ! option for calculating surface roughness length over ocean
- integer, dimension(im), intent(in) :: vegtype
+ integer, dimension(im), intent(in) :: vegtype
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 ! added by s.lu
+ logical, dimension(im), intent(in) :: flag_iter, wet, dry, icy
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) :: &
- & ps,u1,v1,t1,q1,z1,prsl1,prslki,prsik1,prslk1, &
- & ddvel, sigmaf,shdmax, &
+ & ps,t1,q1,z1,prsl1,prslki,prsik1,prslk1, &
+ & wind,sigmaf,shdmax, &
& z0pert,ztpert ! mg, sfc-perts
real(kind=kind_phys), dimension(im), intent(in) :: &
& tskin_ocn, tskin_lnd, tskin_ice, &
@@ -118,24 +114,18 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in)
& fm_ocn, fm_lnd, fm_ice, &
& fh_ocn, fh_lnd, fh_ice, &
& fm10_ocn, fm10_lnd, fm10_ice, &
- & fh2_ocn, fh2_lnd, fh2_ice, &
- & wind
+ & fh2_ocn, fh2_lnd, fh2_ice
character(len=*), intent(out) :: errmsg
integer, intent(out) :: errflg
!
! locals
!
- real(kind=kind_phys), dimension(im) :: wind10m
-
integer i
!
- real(kind=kind_phys) :: qs1, rat, thv1, restar,
- & czilc, tem1, tem2
+ real(kind=kind_phys) :: rat, thv1, restar, wind10m,
+ & czilc, tem1, tem2, virtfac
- real(kind=kind_phys) :: tvs_ocn, tvs_lnd, tvs_ice, &
- & z0_ocn, z0_lnd, z0_ice, &
- & z0max_ocn,z0max_lnd,z0max_ice, &
- & ztmax_ocn,ztmax_lnd,ztmax_ice
+ real(kind=kind_phys) :: tvs, z0, z0max, ztmax
!
real(kind=kind_phys), parameter ::
& charnock=.014, z0s_max=.317e-2 &! a limiting value at high winds over sea
@@ -170,73 +160,21 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in)
! write(0,*)'in sfc_diff, sfc_z0_type=',sfc_z0_type
do i=1,im
-
- ztmax_ocn = 0.0 ; ztmax_lnd = 0.0 ; ztmax_ice = 0.0
-
- wind10m(i) = max(sqrt( u10m(i)*u10m(i) + v10m(i)*v10m(i)),
- & 1.0)
-
if(flag_iter(i)) then
- wind(i) = max(sqrt(u1(i)*u1(i) + v1(i)*v1(i))
- & + max(0.0, min(ddvel(i), 30.0)), 1.0)
- tem1 = 1.0 + rvrdm1 * max(q1(i),1.e-8)
-#ifdef GSD_SURFACE_FLUXES_BUGFIX
- thv1 = t1(i) / prslk1(i) * tem1
- tvs_lnd = 0.5 * (tsurf_lnd(i)+tskin_lnd(i))/prsik1(i) * tem1
- tvs_ice = 0.5 * (tsurf_ice(i)+tskin_ice(i))/prsik1(i) * tem1
- tvs_ocn = 0.5 * (tsurf_ocn(i)+tskin_ocn(i))/prsik1(i) * tem1
-#else
- thv1 = t1(i) * prslki(i) * tem1
- tvs_lnd = 0.5 * (tsurf_lnd(i)+tskin_lnd(i)) * tem1
- tvs_ice = 0.5 * (tsurf_ice(i)+tskin_ice(i)) * tem1
- tvs_ocn = 0.5 * (tsurf_ocn(i)+tskin_ocn(i)) * tem1
-#endif
- qs1 = fpvs(t1(i))
- qs1 = max(1.0e-8, eps * qs1 / (prsl1(i) + epsm1 * qs1))
-
- z0_lnd = 0.01 * z0rl_lnd(i)
- z0max_lnd = max(1.0e-6, min(z0_lnd,z1(i)))
- z0_ice = 0.01 * z0rl_ice(i)
- z0max_ice = max(1.0e-6, min(z0_ice,z1(i)))
- z0_ocn = 0.01 * z0rl_ocn(i)
- z0max_ocn = max(1.0e-6, min(z0_ocn,z1(i)))
+ virtfac = 1.0 + rvrdm1 * max(q1(i),1.e-8)
+ thv1 = t1(i) * prslki(i) * virtfac
! compute stability dependent exchange coefficients
! this portion of the code is presently suppressed
!
-
- if (wet(i)) then ! some open ocean
- ustar_ocn(i) = sqrt(grav * z0_ocn / charnock)
-
-!** test xubin's new z0
-
-! ztmax = z0max
-
- restar = max(ustar_ocn(i)*z0max_ocn*visi, 0.000001)
-
-! restar = log(restar)
-! restar = min(restar,5.)
-! restar = max(restar,-5.)
-! rat = aa1 + (bb1 + cc1*restar) * restar
-! rat = rat / (1. + (bb2 + cc2*restar) * restar))
-! rat taken from zeng, zhao and dickinson 1997
-
- rat = min(7.0, 2.67 * sqrt(sqrt(restar)) - 2.57)
- ztmax_ocn = z0max_ocn * exp(-rat)
-
- if (sfc_z0_type == 6) then
- call znot_t_v6(wind10m(i),ztmax_ocn) ! 10-m wind,m/s, ztmax(m)
- else if (sfc_z0_type == 7) then
- call znot_t_v7(wind10m(i),ztmax_ocn) ! 10-m wind,m/s, ztmax(m)
- else if (sfc_z0_type .ne. 0) then
- write(0,*)'no option for sfc_z0_type=',sfc_z0_type
- stop
- endif
-
- endif ! Open ocean
-
- if (dry(i) .or. icy(i)) then ! over land or sea ice
-!** xubin's new z0 over land and sea ice
+ if (dry(i)) then ! Some land
+#ifdef GSD_SURFACE_FLUXES_BUGFIX
+ tvs = 0.5 * (tsurf_lnd(i)+tskin_lnd(i))/prsik1(i) * virtfac
+#else
+ tvs = 0.5 * (tsurf_lnd(i)+tskin_lnd(i)) * virtfac
+#endif
+ z0max = max(1.0e-6, min(0.01 * z0rl_lnd(i), z1(i)))
+!** xubin's new z0 over land
tem1 = 1.0 - shdmax(i)
tem2 = tem1 * tem1
tem1 = 1.0 - tem2
@@ -244,134 +182,175 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in)
if( ivegsrc == 1 ) then
if (vegtype(i) == 10) then
- z0max_lnd = exp( tem2*log01 + tem1*log07 )
+ z0max = exp( tem2*log01 + tem1*log07 )
elseif (vegtype(i) == 6) then
- z0max_lnd = exp( tem2*log01 + tem1*log05 )
+ z0max = exp( tem2*log01 + tem1*log05 )
elseif (vegtype(i) == 7) then
! z0max = exp( tem2*log01 + tem1*log01 )
- z0max_lnd = 0.01
+ z0max = 0.01
elseif (vegtype(i) == 16) then
! z0max = exp( tem2*log01 + tem1*log01 )
- z0max_lnd = 0.01
+ z0max = 0.01
else
- z0max_lnd = exp( tem2*log01 + tem1*log(z0max_lnd) )
+ z0max = exp( tem2*log01 + tem1*log(z0max) )
endif
elseif (ivegsrc == 2 ) then
- if (vegtype(i) == 7) then
- z0max_lnd = exp( tem2*log01 + tem1*log07 )
- elseif (vegtype(i) == 8) then
- z0max_lnd = exp( tem2*log01 + tem1*log05 )
- elseif (vegtype(i) == 9) then
-! z0max = exp( tem2*log01 + tem1*log01 )
- z0max_lnd = 0.01
- elseif (vegtype(i) == 11) then
-! z0max = exp( tem2*log01 + tem1*log01 )
- z0max_lnd = 0.01
- else
- z0max_lnd = exp( tem2*log01 + tem1*log(z0max_lnd) )
- endif
-
- endif ! over land or sea ice
-
- z0max_ice = z0max_lnd
+ if (vegtype(i) == 7) then
+ z0max = exp( tem2*log01 + tem1*log07 )
+ elseif (vegtype(i) == 8) then
+ z0max = exp( tem2*log01 + tem1*log05 )
+ elseif (vegtype(i) == 9) then
+! z0max = exp( tem2*log01 + tem1*log01 )
+ z0max = 0.01
+ elseif (vegtype(i) == 11) then
+! z0max = exp( tem2*log01 + tem1*log01 )
+ z0max = 0.01
+ else
+ z0max = exp( tem2*log01 + tem1*log(z0max) )
+ endif
+ endif
! mg, sfc-perts: add surface perturbations to z0max over land
- if (dry(i) .and. z0pert(i) /= 0.0 ) then
- z0max_lnd = z0max_lnd * (10.**z0pert(i))
+ if (z0pert(i) /= 0.0 ) then
+ z0max = z0max * (10.**z0pert(i))
endif
- z0max_lnd = max(z0max_lnd,1.0e-6)
- z0max_ice = max(z0max_ice,1.0e-6)
+ z0max = max(z0max, 1.0e-6)
! czilc = 10.0 ** (- (0.40/0.07) * z0) ! fei's canopy height dependance of czil
czilc = 0.8
- tem1 = 1.0 - sigmaf(i)
- ztmax_lnd = z0max_lnd*exp( - tem1*tem1
+ tem1 = 1.0 - sigmaf(i)
+ ztmax = z0max*exp( - tem1*tem1
& * czilc*ca*sqrt(ustar_lnd(i)*(0.01/1.5e-05)))
- ztmax_ice = z0max_ice*exp( - tem1*tem1
- & * czilc*ca*sqrt(ustar_ice(i)*(0.01/1.5e-05)))
! mg, sfc-perts: add surface perturbations to ztmax/z0max ratio over land
- if (dry(i) .and. ztpert(i) /= 0.0) then
- ztmax_lnd = ztmax_lnd * (10.**ztpert(i))
+ if (ztpert(i) /= 0.0) then
+ ztmax = ztmax * (10.**ztpert(i))
endif
+ ztmax = max(ztmax, 1.0e-6)
+!
+ call stability
+! --- inputs:
+ & (z1(i), snwdph_lnd(i), thv1, wind(i),
+ & z0max, ztmax, 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))
+ endif ! Dry points
+ if (icy(i)) then ! Some ice
+ tvs = 0.5 * (tsurf_ice(i)+tskin_ice(i)) * virtfac
+ z0max = max(1.0e-6, min(0.01 * z0rl_ice(i), z1(i)))
+!** xubin's new z0 over land and sea ice
+ tem1 = 1.0 - shdmax(i)
+ tem2 = tem1 * tem1
+ tem1 = 1.0 - tem2
- endif ! end of if(sfctype flags) then
+ if( ivegsrc == 1 ) then
- ztmax_ocn = max(ztmax_ocn,1.0e-6)
- ztmax_lnd = max(ztmax_lnd,1.0e-6)
- ztmax_ice = max(ztmax_ice,1.0e-6)
+ z0max = exp( tem2*log01 + tem1*log(z0max) )
+ elseif (ivegsrc == 2 ) then
+ z0max = exp( tem2*log01 + tem1*log(z0max) )
+ endif
-! BWG begin "stability" block, 2019-03-23
- if (wet(i)) then ! Some open ocean
- call stability
-! --- inputs:
- & (z1(i),snwdph_ocn(i),thv1,wind(i),
- & z0max_ocn,ztmax_ocn,tvs_ocn,grav,
-! --- outputs:
- & rb_ocn(i),fm_ocn(i),fh_ocn(i),fm10_ocn(i),fh2_ocn(i),
- & cm_ocn(i),ch_ocn(i),stress_ocn(i),ustar_ocn(i))
- endif ! Open ocean points
+ z0max = max(z0max, 1.0e-6)
- if (dry(i)) then ! Some land
- call stability
-! --- inputs:
- & (z1(i),snwdph_lnd(i),thv1,wind(i),
- & z0max_lnd,ztmax_lnd,tvs_lnd,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))
- endif ! Dry points
+! czilc = 10.0 ** (- (0.40/0.07) * z0) ! fei's canopy height
+! dependance of czil
+ czilc = 0.8
- if (icy(i)) then ! Some ice
- call stability
+ tem1 = 1.0 - sigmaf(i)
+ ztmax = z0max*exp( - tem1*tem1
+ & * czilc*ca*sqrt(ustar_ice(i)*(0.01/1.5e-05)))
+ ztmax = max(ztmax, 1.0e-6)
+!
+ call stability
! --- inputs:
- & (z1(i),snwdph_ice(i),thv1,wind(i),
- & z0max_ice,ztmax_ice,tvs_ice,grav,
+ & (z1(i), snwdph_ice(i), thv1, wind(i),
+ & z0max, ztmax, 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))
+ & 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))
endif ! Icy points
! BWG: Everything from here to end of subroutine was after
! the stuff now put into "stability"
+ if (wet(i)) then ! Some open ocean
+ tvs = 0.5 * (tsurf_ocn(i)+tskin_ocn(i)) * virtfac
+ z0 = 0.01 * z0rl_ocn(i)
+ z0max = max(1.0e-6, min(z0,z1(i)))
+ ustar_ocn(i) = sqrt(grav * z0 / charnock)
+ wind10m = sqrt(u10m(i)*u10m(i)+v10m(i)*v10m(i))
+
+!** test xubin's new z0
+
+! ztmax = z0max
+
+ restar = max(ustar_ocn(i)*z0max*visi, 0.000001)
+
+! restar = log(restar)
+! restar = min(restar,5.)
+! restar = max(restar,-5.)
+! rat = aa1 + (bb1 + cc1*restar) * restar
+! rat = rat / (1. + (bb2 + cc2*restar) * restar))
+! rat taken from zeng, zhao and dickinson 1997
+
+ rat = min(7.0, 2.67 * sqrt(sqrt(restar)) - 2.57)
+ ztmax = max(z0max * exp(-rat), 1.0e-6)
+!
+ if (sfc_z0_type == 6) then
+ call znot_t_v6(wind10m, ztmax) ! 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)
+ else if (sfc_z0_type /= 0) then
+ write(0,*)'no option for sfc_z0_type=',sfc_z0_type
+ stop
+ endif
+!
+ call stability
+! --- inputs:
+ & (z1(i), snwdph_ocn(i), thv1, wind(i),
+ & z0max, ztmax, tvs, grav,
+! --- outputs:
+ & rb_ocn(i), fm_ocn(i), fh_ocn(i), fm10_ocn(i), fh2_ocn(i),
+ & cm_ocn(i), ch_ocn(i), stress_ocn(i), ustar_ocn(i))
!
! update z0 over ocean
!
- if (wet(i)) then
- z0_ocn = (charnock / grav) * ustar_ocn(i) * ustar_ocn(i)
+ if (sfc_z0_type == 0) then
+ z0 = (charnock / grav) * ustar_ocn(i) * ustar_ocn(i)
! mbek -- toga-coare flux algorithm
-! z0 = (charnock / grav) * ustar(i)*ustar(i) + arnu/ustar(i)
+! z0 = (charnock / grav) * ustar(i)*ustar(i) + arnu/ustar(i)
! new implementation of z0
-! cc = ustar(i) * z0 / rnu
-! pp = cc / (1. + cc)
-! ff = grav * arnu / (charnock * ustar(i) ** 3)
-! z0 = arnu / (ustar(i) * ff ** pp)
+! cc = ustar(i) * z0 / rnu
+! pp = cc / (1. + cc)
+! ff = grav * arnu / (charnock * ustar(i) ** 3)
+! z0 = arnu / (ustar(i) * ff ** pp)
+
+ if (redrag) then
+ z0rl_ocn(i) = 100.0 * max(min(z0, z0s_max), 1.e-7)
+ else
+ z0rl_ocn(i) = 100.0 * max(min(z0,.1), 1.e-7)
+ endif
- if (redrag) then
- z0rl_ocn(i) = 100.0 * max(min(z0_ocn, z0s_max), 1.e-7)
+ elseif (sfc_z0_type == 6) then ! wang
+ call znot_m_v6(wind10m, z0) ! wind, m/s, z0, m
+ z0rl_ocn(i) = 100.0 * z0 ! cm
+ elseif (sfc_z0_type == 7) then ! wang
+ call znot_m_v7(wind10m, z0) ! wind, m/s, z0, m
+ z0rl_ocn(i) = 100.0 * z0 ! cm
else
- z0rl_ocn(i) = 100.0 * max(min(z0_ocn,.1), 1.e-7)
+ z0rl_ocn(i) = 1.0e-4
endif
- if (sfc_z0_type == 6) then ! wang
- call znot_m_v6(wind10m(i),z0_ocn) ! wind, m/s, z0, m
- z0rl_ocn(i) = 100.0 * z0_ocn ! cm
- endif !wang
- if (sfc_z0_type == 7) then ! wang
- call znot_m_v7(wind10m(i),z0_ocn) ! wind, m/s, z0, m
- z0rl_ocn(i) = 100.0 * z0_ocn ! cm
- endif !wang
-
-
endif ! end of if(open ocean)
+!
endif ! end of if(flagiter) loop
enddo
@@ -382,8 +361,11 @@ end subroutine sfc_diff_run
!----------------------------------------
!>\ingroup GFS_diff_main
subroutine stability &
- & ( z1, snwdph, thv1, wind, z0max, ztmax, tvs, grav, & ! --- inputs:
- & rb, fm, fh, fm10, fh2, cm, ch, stress, ustar) ! --- outputs:
+! --- inputs:
+ & ( z1, snwdph, thv1, wind, z0max, ztmax, tvs, grav, &
+! --- outputs:
+ & rb, fm, fh, fm10, fh2, cm, ch, stress, ustar)
+!-----
! --- inputs:
real(kind=kind_phys), intent(in) :: &
@@ -431,10 +413,10 @@ subroutine stability &
#endif
tem1 = 1.0 / z0max
tem2 = 1.0 / ztmax
- fm = log((z0max+z1) * tem1)
- fh = log((ztmax+z1) * tem2)
- fm10 = log((z0max+10.) * tem1)
- fh2 = log((ztmax+2.) * tem2)
+ fm = log((z0max+z1) * tem1)
+ fh = log((ztmax+z1) * tem2)
+ fm10 = log((z0max+10.) * tem1)
+ fh2 = log((ztmax+2.) * tem2)
hlinf = rb * fm * fm / fh
hlinf = min(max(hlinf,ztmin1),ztmax1)
!
@@ -543,8 +525,9 @@ end subroutine stability
!! add fitted z0,zt curves for hurricane application (used in HWRF/HMON)
!! Weiguo Wang, 2019-0425
- SUBROUTINE znot_m_v6(uref,znotm)
- IMPLICIT NONE
+ SUBROUTINE znot_m_v6(uref, znotm)
+ use machine , only : kind_phys
+ IMPLICIT NONE
! Calculate areodynamical roughness over water with input 10-m wind
! For low-to-moderate winds, try to match the Cd-U10 relationship from COARE V3.5 (Edson et al. 2013)
! For high winds, try to fit available observational data
@@ -555,53 +538,42 @@ SUBROUTINE znot_m_v6(uref,znotm)
! znotm(meter): areodynamical roughness scale over water
!
- REAL, INTENT(IN) :: uref
- REAL, INTENT(OUT):: znotm
- REAL :: p13, p12, p11, p10
- REAL :: p25, p24, p23, p22, p21, p20
- REAL :: p35, p34, p33, p32, p31, p30
- REAL :: p40
-
- p13 = -1.296521881682694e-02
- p12 = 2.855780863283819e-01
- p11 = -1.597898515251717e+00
- p10 = -8.396975715683501e+00
+ REAL(kind=kind_phys), INTENT(IN) :: uref
+ REAL(kind=kind_phys), INTENT(OUT):: znotm
+ real(kind=kind_phys), parameter :: p13 = -1.296521881682694e-02,
+ & p12 = 2.855780863283819e-01, p11 = -1.597898515251717e+00,
+ & p10 = -8.396975715683501e+00,
- p25 = 3.790846746036765e-10
- p24 = 3.281964357650687e-09
- p23 = 1.962282433562894e-07
- p22 = -1.240239171056262e-06
- p21 = 1.739759082358234e-07
- p20 = 2.147264020369413e-05
+ & p25 = 3.790846746036765e-10, p24 = 3.281964357650687e-09,
+ & p23 = 1.962282433562894e-07, p22 = -1.240239171056262e-06,
+ & p21 = 1.739759082358234e-07, p20 = 2.147264020369413e-05,
- p35 = 1.840430200185075e-07
- p34 = -2.793849676757154e-05
- p33 = 1.735308193700643e-03
- p32 = -6.139315534216305e-02
- p31 = 1.255457892775006e+00
- p30 = -1.663993561652530e+01
+ & p35 = 1.840430200185075e-07, p34 = -2.793849676757154e-05,
+ & p33 = 1.735308193700643e-03, p32 = -6.139315534216305e-02,
+ & p31 = 1.255457892775006e+00, p30 = -1.663993561652530e+01,
- p40 = 4.579369142033410e-04
+ & p40 = 4.579369142033410e-04
+
if (uref >= 0.0 .and. uref <= 6.5 ) then
- znotm = exp( p10 + p11*uref + p12*uref**2 +
- & p13*uref**3)
+ znotm = exp(p10 + uref * (p11 + uref * (p12 + uref*p13)))
elseif (uref > 6.5 .and. uref <= 15.7) then
- znotm = p25*uref**5 + p24*uref**4 + p23*uref**3 +
- & p22*uref**2 + p21*uref + p20
+ znotm = p20 + uref * (p21 + uref * (p22 + uref * (p23
+ & + uref * (p24 + uref * p25))))
elseif (uref > 15.7 .and. uref <= 53.0) then
- znotm = exp( p35*uref**5 + p34*uref**4 +
- & p33*uref**3 + p32*uref**2 + p31*uref + p30 )
+ znotm = exp( p30 + uref * (p31 + uref * (p32 + uref * (p33
+ & + uref * (p34 + uref * p35)))))
elseif ( uref > 53.0) then
znotm = p40
else
print*, 'Wrong input uref value:',uref
endif
- END SUBROUTINE znot_m_v6
+ END SUBROUTINE znot_m_v6
- SUBROUTINE znot_t_v6(uref,znott)
- IMPLICIT NONE
+ SUBROUTINE znot_t_v6(uref, znott)
+ use machine , only : kind_phys
+ IMPLICIT NONE
! Calculate scalar roughness over water with input 10-m wind
! For low-to-moderate winds, try to match the Ck-U10 relationship from COARE algorithm
! For high winds, try to retain the Ck-U10 relationship of FY2015 HWRF
@@ -612,85 +584,61 @@ SUBROUTINE znot_t_v6(uref,znott)
! znott(meter): scalar roughness scale over water
!
- REAL, INTENT(IN) :: uref
- REAL, INTENT(OUT):: znott
-
- REAL :: p00
- REAL :: p15, p14, p13, p12, p11, p10
- REAL :: p25, p24, p23, p22, p21, p20
- REAL :: p35, p34, p33, p32, p31, p30
- REAL :: p45, p44, p43, p42, p41, p40
- REAL :: p56, p55, p54, p53, p52, p51, p50
- REAL :: p60
-
- p00 = 1.100000000000000e-04
-
- p15 = -9.144581627678278e-10
- p14 = 7.020346616456421e-08
- p13 = -2.155602086883837e-06
- p12 = 3.333848806567684e-05
- p11 = -2.628501274963990e-04
- p10 = 8.634221567969181e-04
-
- p25 = -8.654513012535990e-12
- p24 = 1.232380050058077e-09
- p23 = -6.837922749505057e-08
- p22 = 1.871407733439947e-06
- p21 = -2.552246987137160e-05
- p20 = 1.428968311457630e-04
-
- p35 = 3.207515102100162e-12
- p34 = -2.945761895342535e-10
- p33 = 8.788972147364181e-09
- p32 = -3.814457439412957e-08
- p31 = -2.448983648874671e-06
- p30 = 3.436721779020359e-05
-
- p45 = -3.530687797132211e-11
- p44 = 3.939867958963747e-09
- p43 = -1.227668406985956e-08
- p42 = -1.367469811838390e-05
- p41 = 5.988240863928883e-04
- p40 = -7.746288511324971e-03
-
- p56 = -1.187982453329086e-13
- p55 = 4.801984186231693e-11
- p54 = -8.049200462388188e-09
- p53 = 7.169872601310186e-07
- p52 = -3.581694433758150e-05
- p51 = 9.503919224192534e-04
- p50 = -1.036679430885215e-02
-
- p60 = 4.751256171799112e-05
-
- if (uref >= 0.0 .and. uref < 5.9 ) then
+ REAL(kind=kind_phys), INTENT(IN) :: uref
+ REAL(kind=kind_phys), INTENT(OUT):: znott
+ real(kind=kind_phys), parameter :: p00 = 1.100000000000000e-04,
+ & p15 = -9.144581627678278e-10, p14 = 7.020346616456421e-08,
+ & p13 = -2.155602086883837e-06, p12 = 3.333848806567684e-05,
+ & p11 = -2.628501274963990e-04, p10 = 8.634221567969181e-04,
+
+ & p25 = -8.654513012535990e-12, p24 = 1.232380050058077e-09,
+ & p23 = -6.837922749505057e-08, p22 = 1.871407733439947e-06,
+ & p21 = -2.552246987137160e-05, p20 = 1.428968311457630e-04,
+
+ & p35 = 3.207515102100162e-12, p34 = -2.945761895342535e-10,
+ & p33 = 8.788972147364181e-09, p32 = -3.814457439412957e-08,
+ & p31 = -2.448983648874671e-06, p30 = 3.436721779020359e-05,
+
+ & p45 = -3.530687797132211e-11, p44 = 3.939867958963747e-09,
+ & p43 = -1.227668406985956e-08, p42 = -1.367469811838390e-05,
+ & p41 = 5.988240863928883e-04, p40 = -7.746288511324971e-03,
+
+ & p56 = -1.187982453329086e-13, p55 = 4.801984186231693e-11,
+ & p54 = -8.049200462388188e-09, p53 = 7.169872601310186e-07,
+ & p52 = -3.581694433758150e-05, p51 = 9.503919224192534e-04,
+ & p50 = -1.036679430885215e-02,
+
+ & p60 = 4.751256171799112e-05
+
+ if (uref >= 0.0 .and. uref < 5.9 ) then
znott = p00
- elseif (uref >= 5.9 .and. uref <= 15.4) then
- znott = p15*uref**5 + p14*uref**4 + p13*uref**3
- & + p12*uref**2 + p11*uref + p10
- elseif (uref > 15.4 .and. uref <= 21.6) then
- znott = p25*uref**5 + p24*uref**4 + p23*uref**3
- & + p22*uref**2 + p21*uref + p20
- elseif (uref > 21.6 .and. uref <= 42.2) then
- znott = p35*uref**5 + p34*uref**4 + p33*uref**3
- & + p32*uref**2 + p31*uref + p30
- elseif ( uref > 42.2 .and. uref <= 53.3) then
- znott = p45*uref**5 + p44*uref**4 + p43*uref**3
- & + p42*uref**2 + p41*uref + p40
- elseif ( uref > 53.3 .and. uref <= 80.0) then
- znott = p56*uref**6 + p55*uref**5 + p54*uref**4
- & + p53*uref**3 + p52*uref**2 + p51*uref + p50
- elseif ( uref > 80.0) then
+ elseif (uref >= 5.9 .and. uref <= 15.4) then
+ znott = p10 + uref * (p11 + uref * (p12 + uref * (p13
+ & + uref * (p14 + uref * p15))))
+ elseif (uref > 15.4 .and. uref <= 21.6) then
+ znott = p20 + uref * (p21 + uref * (p22 + uref * (p23
+ & + uref * (p24 + uref * p25))))
+ elseif (uref > 21.6 .and. uref <= 42.2) then
+ znott = p30 + uref * (p31 + uref * (p32 + uref * (p33
+ & + uref * (p34 + uref * p35))))
+ elseif ( uref > 42.2 .and. uref <= 53.3) then
+ znott = p40 + uref * (p41 + uref * (p42 + uref * (p43
+ & + uref * (p44 + uref * p45))))
+ elseif ( uref > 53.3 .and. uref <= 80.0) then
+ znott = p50 + uref * (p51 + uref * (p52 + uref * (p53
+ & + uref * (p54 + uref * (p55 + uref * p56)))))
+ elseif ( uref > 80.0) then
znott = p60
- else
+ else
print*, 'Wrong input uref value:',uref
- endif
+ endif
- END SUBROUTINE znot_t_v6
+ END SUBROUTINE znot_t_v6
- SUBROUTINE znot_m_v7(uref,znotm)
- IMPLICIT NONE
+ SUBROUTINE znot_m_v7(uref, znotm)
+ use machine , only : kind_phys
+ IMPLICIT NONE
! Calculate areodynamical roughness over water with input 10-m wind
! For low-to-moderate winds, try to match the Cd-U10 relationship from COARE V3.5 (Edson et al. 2013)
! For high winds, try to fit available observational data
@@ -702,52 +650,41 @@ SUBROUTINE znot_m_v7(uref,znotm)
! znotm(meter): areodynamical roughness scale over water
!
- REAL, INTENT(IN) :: uref
- REAL, INTENT(OUT):: znotm
- REAL :: p13, p12, p11, p10
- REAL :: p25, p24, p23, p22, p21, p20
- REAL :: p35, p34, p33, p32, p31, p30
- REAL :: p40
-
- p13 = -1.296521881682694e-02
- p12 = 2.855780863283819e-01
- p11 = -1.597898515251717e+00
- p10 = -8.396975715683501e+00
+ REAL(kind=kind_phys), INTENT(IN) :: uref
+ REAL(kind=kind_phys), INTENT(OUT):: znotm
- p25 = 3.790846746036765e-10
- p24 = 3.281964357650687e-09
- p23 = 1.962282433562894e-07
- p22 = -1.240239171056262e-06
- p21 = 1.739759082358234e-07
- p20 = 2.147264020369413e-05
+ real(kind=kind_phys), parameter :: p13 = -1.296521881682694e-02,
+ & p12 = 2.855780863283819e-01, p11 = -1.597898515251717e+00,
+ & p10 = -8.396975715683501e+00,
+ & p25 = 3.790846746036765e-10, p24 = 3.281964357650687e-09,
+ & p23 = 1.962282433562894e-07, p22 = -1.240239171056262e-06,
+ & p21 = 1.739759082358234e-07, p20 = 2.147264020369413e-05,
- p35 = 1.897534489606422e-07
- p34 = -3.019495980684978e-05
- p33 = 1.931392924987349e-03
- p32 = -6.797293095862357e-02
- p31 = 1.346757797103756e+00
- p30 = -1.707846930193362e+01
+ & p35 = 1.897534489606422e-07, p34 = -3.019495980684978e-05,
+ & p33 = 1.931392924987349e-03, p32 = -6.797293095862357e-02,
+ & p31 = 1.346757797103756e+00, p30 = -1.707846930193362e+01,
- p40 = 3.371427455376717e-04
+ & p40 = 3.371427455376717e-04
- if (uref >= 0.0 .and. uref <= 6.5 ) then
- znotm = exp( p10 + p11*uref + p12*uref**2 + p13*uref**3)
- elseif (uref > 6.5 .and. uref <= 15.7) then
- znotm = p25*uref**5 + p24*uref**4 + p23*uref**3 +
- & p22*uref**2 + p21*uref + p20
- elseif (uref > 15.7 .and. uref <= 53.0) then
- znotm = exp( p35*uref**5 + p34*uref**4 + p33*uref**3
- & + p32*uref**2 + p31*uref + p30 )
- elseif ( uref > 53.0) then
+ if (uref >= 0.0 .and. uref <= 6.5 ) then
+ znotm = exp( p10 + uref * (p11 + uref * (p12 + uref * p13)))
+ elseif (uref > 6.5 .and. uref <= 15.7) then
+ znotm = p20 + uref * (p21 + uref * (p22 + uref * (p23
+ & + uref * (p24 + uref * p25))))
+ elseif (uref > 15.7 .and. uref <= 53.0) then
+ znotm = exp( p30 + uref * (p31 + uref * (p32 + uref * (p33
+ & + uref * (p34 + uref * p35)))))
+ elseif ( uref > 53.0) then
znotm = p40
- else
+ else
print*, 'Wrong input uref value:',uref
- endif
+ endif
END SUBROUTINE znot_m_v7
- SUBROUTINE znot_t_v7(uref,znott)
- IMPLICIT NONE
+ SUBROUTINE znot_t_v7(uref, znott)
+ use machine , only : kind_phys
+ IMPLICIT NONE
! Calculate scalar roughness over water with input 10-m wind
! For low-to-moderate winds, try to match the Ck-U10 relationship from COARE algorithm
! For high winds, try to retain the Ck-U10 relationship of FY2015 HWRF
@@ -759,79 +696,54 @@ SUBROUTINE znot_t_v7(uref,znott)
! znott(meter): scalar roughness scale over water
!
- REAL, INTENT(IN) :: uref
- REAL, INTENT(OUT):: znott
-
- REAL :: p00
- REAL :: p15, p14, p13, p12, p11, p10
- REAL :: p25, p24, p23, p22, p21, p20
- REAL :: p35, p34, p33, p32, p31, p30
- REAL :: p45, p44, p43, p42, p41, p40
- REAL :: p56, p55, p54, p53, p52, p51, p50
- REAL :: p60
-
- p00 = 1.100000000000000e-04
-
- p15 = -9.193764479895316e-10
- p14 = 7.052217518653943e-08
- p13 = -2.163419217747114e-06
- p12 = 3.342963077911962e-05
- p11 = -2.633566691328004e-04
- p10 = 8.644979973037803e-04
-
- p25 = -9.402722450219142e-12
- p24 = 1.325396583616614e-09
- p23 = -7.299148051141852e-08
- p22 = 1.982901461144764e-06
- p21 = -2.680293455916390e-05
- p20 = 1.484341646128200e-04
-
- p35 = 7.921446674311864e-12
- p34 = -1.019028029546602e-09
- p33 = 5.251986927351103e-08
- p32 = -1.337841892062716e-06
- p31 = 1.659454106237737e-05
- p30 = -7.558911792344770e-05
-
- p45 = -2.694370426850801e-10
- p44 = 5.817362913967911e-08
- p43 = -5.000813324746342e-06
- p42 = 2.143803523428029e-04
- p41 = -4.588070983722060e-03
- p40 = 3.924356617245624e-02
-
- p56 = -1.663918773476178e-13
- p55 = 6.724854483077447e-11
- p54 = -1.127030176632823e-08
- p53 = 1.003683177025925e-06
- p52 = -5.012618091180904e-05
- p51 = 1.329762020689302e-03
- p50 = -1.450062148367566e-02
-
- p60 = 6.840803042788488e-05
+ REAL(kind=kind_phys), INTENT(IN) :: uref
+ REAL(kind=kind_phys), INTENT(OUT):: znott
+
+ real(kind=kind_phys), parameter :: p00 = 1.100000000000000e-04,
+
+ & p15 = -9.193764479895316e-10, p14 = 7.052217518653943e-08,
+ & p13 = -2.163419217747114e-06, p12 = 3.342963077911962e-05,
+ & p11 = -2.633566691328004e-04, p10 = 8.644979973037803e-04,
+
+ & p25 = -9.402722450219142e-12, p24 = 1.325396583616614e-09,
+ & p23 = -7.299148051141852e-08, p22 = 1.982901461144764e-06,
+ & p21 = -2.680293455916390e-05, p20 = 1.484341646128200e-04,
+
+ & p35 = 7.921446674311864e-12, p34 = -1.019028029546602e-09,
+ & p33 = 5.251986927351103e-08, p32 = -1.337841892062716e-06,
+ & p31 = 1.659454106237737e-05, p30 = -7.558911792344770e-05,
+
+ & p45 = -2.694370426850801e-10, p44 = 5.817362913967911e-08,
+ & p43 = -5.000813324746342e-06, p42 = 2.143803523428029e-04,
+ & p41 = -4.588070983722060e-03, p40 = 3.924356617245624e-02,
+
+ & p56 = -1.663918773476178e-13, p55 = 6.724854483077447e-11,
+ & p54 = -1.127030176632823e-08, p53 = 1.003683177025925e-06,
+ & p52 = -5.012618091180904e-05, p51 = 1.329762020689302e-03,
+ & p50 = -1.450062148367566e-02, p60 = 6.840803042788488e-05
if (uref >= 0.0 .and. uref < 5.9 ) then
- znott = p00
- elseif (uref >= 5.9 .and. uref <= 15.4) then
- znott = p15*uref**5 + p14*uref**4 + p13*uref**3 +
- & p12*uref**2 + p11*uref + p10
- elseif (uref > 15.4 .and. uref <= 21.6) then
- znott = p25*uref**5 + p24*uref**4 + p23*uref**3 +
- & p22*uref**2 + p21*uref + p20
- elseif (uref > 21.6 .and. uref <= 42.6) then
- znott = p35*uref**5 + p34*uref**4 + p33*uref**3 +
- & p32*uref**2 + p31*uref + p30
- elseif ( uref > 42.6 .and. uref <= 53.0) then
- znott = p45*uref**5 + p44*uref**4 + p43*uref**3 +
- & p42*uref**2 + p41*uref + p40
- elseif ( uref > 53.0 .and. uref <= 80.0) then
- znott = p56*uref**6 + p55*uref**5 + p54*uref**4 +
- & p53*uref**3 + p52*uref**2 + p51*uref + p50
- elseif ( uref > 80.0) then
+ znott = p00
+ elseif (uref >= 5.9 .and. uref <= 15.4) then
+ znott = p10 + uref * (p11 + uref * (p12 + uref * (p13
+ & + uref * (p14 + uref * p15))))
+ elseif (uref > 15.4 .and. uref <= 21.6) then
+ znott = p20 + uref * (p21 + uref * (p22 + uref * (p23
+ & + uref * (p24 + uref * p25))))
+ elseif (uref > 21.6 .and. uref <= 42.6) then
+ znott = p30 + uref * (p31 + uref * (p32 + uref * (p33
+ & + uref * (p34 + uref * p35))))
+ elseif ( uref > 42.6 .and. uref <= 53.0) then
+ znott = p40 + uref * (p41 + uref * (p42 + uref * (p43
+ & + uref * (p44 + uref * p45))))
+ elseif ( uref > 53.0 .and. uref <= 80.0) then
+ znott = p50 + uref * (p51 + uref * (p52 + uref * (p53
+ & + uref * (p54 + uref * (p55 + uref * p56)))))
+ elseif ( uref > 80.0) then
znott = p60
else
print*, 'Wrong input uref value:',uref
- endif
+ endif
END SUBROUTINE znot_t_v7
diff --git a/physics/sfc_diff.meta b/physics/sfc_diff.meta
index de8acc72a..232b0050f 100644
--- a/physics/sfc_diff.meta
+++ b/physics/sfc_diff.meta
@@ -54,24 +54,6 @@
kind = kind_phys
intent = in
optional = F
-[u1]
- standard_name = x_wind_at_lowest_model_layer
- long_name = x component of 1st model layer wind
- units = m s-1
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = in
- optional = F
-[v1]
- standard_name = y_wind_at_lowest_model_layer
- long_name = y component of 1st model layer wind
- units = m s-1
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = in
- optional = F
[t1]
standard_name = air_temperature_at_lowest_model_layer
long_name = 1st model layer air temperature
@@ -99,6 +81,15 @@
kind = kind_phys
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_dimension)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
[prsl1]
standard_name = air_pressure_at_lowest_model_layer
long_name = Model layer 1 mean pressure
@@ -135,15 +126,6 @@
kind = kind_phys
intent = in
optional = F
-[ddvel]
- standard_name = surface_wind_enhancement_due_to_convection
- long_name = surface wind enhancement due to convection
- units = m s-1
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = in
- optional = F
[sigmaf]
standard_name = bounded_vegetation_area_fraction
long_name = areal fractional cover of green vegetation bounded on the bottom
@@ -613,15 +595,6 @@
kind = kind_phys
intent = inout
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_dimension)
- type = real
- kind = kind_phys
- intent = inout
- optional = F
[errmsg]
standard_name = ccpp_error_message
long_name = error message for error handling in CCPP
diff --git a/physics/sfc_drv.f b/physics/sfc_drv.f
index 4e27c07f1..75afaa6ff 100644
--- a/physics/sfc_drv.f
+++ b/physics/sfc_drv.f
@@ -62,9 +62,9 @@ end subroutine lsm_noah_finalize
! !
! call sfc_drv !
! --- inputs: !
-! ( im, km, ps, u1, v1, t1, q1, soiltyp, vegtype, sigmaf, !
+! ( im, km, ps, t1, q1, soiltyp, vegtype, sigmaf, !
! sfcemis, dlwflx, dswsfc, snet, delt, tg3, cm, ch, !
-! prsl1, prslki, zf, land, ddvel, slopetyp, !
+! prsl1, prslki, zf, land, wind, slopetyp, !
! shdmin, shdmax, snoalb, sfalb, flag_iter, flag_guess, !
! lheatstrg, isot, ivegsrc, !
! --- in/outs: !
@@ -94,7 +94,6 @@ end subroutine lsm_noah_finalize
! im - integer, horiz dimention and num of used pts 1 !
! km - integer, vertical soil layer dimension 1 !
! ps - real, surface pressure (pa) im !
-! u1, v1 - real, u/v component of surface layer wind im !
! t1 - real, surface layer mean temperature (k) im !
! q1 - real, surface layer mean specific humidity im !
! soiltyp - integer, soil type (integer index) im !
@@ -112,7 +111,7 @@ end subroutine lsm_noah_finalize
! prslki - real, dimensionless exner function at layer 1 im !
! zf - real, height of bottom layer (m) im !
! land - logical, = T if a point with any land im !
-! ddvel - real, im !
+! wind - real, wind speed (m/s) im !
! slopetyp - integer, class of sfc slope (integer index) im !
! shdmin - real, min fractional coverage of green veg im !
! shdmax - real, max fractnl cover of green veg (not used) im !
@@ -171,10 +170,10 @@ end subroutine lsm_noah_finalize
!> \section general_noah_drv GFS sfc_drv General Algorithm
!> @{
subroutine lsm_noah_run &
- & ( im, km, grav, cp, hvap, rd, eps, epsm1, rvrdm1, ps, u1, & ! --- inputs:
- & v1, t1, q1, soiltyp, vegtype, sigmaf, &
+ & ( im, km, grav, cp, hvap, rd, eps, epsm1, rvrdm1, ps, & ! --- inputs:
+ & t1, q1, soiltyp, vegtype, sigmaf, &
& sfcemis, dlwflx, dswsfc, snet, delt, tg3, cm, ch, &
- & prsl1, prslki, zf, land, ddvel, slopetyp, &
+ & prsl1, prslki, zf, land, wind, slopetyp, &
& shdmin, shdmax, snoalb, sfalb, flag_iter, flag_guess, &
& lheatstrg, isot, ivegsrc, &
& bexppert, xlaipert, vegfpert,pertvegf, & ! sfc perts, mgehne
@@ -212,9 +211,9 @@ subroutine lsm_noah_run &
integer, dimension(im), intent(in) :: soiltyp, vegtype, slopetyp
- real (kind=kind_phys), dimension(im), intent(in) :: ps, u1, v1, &
+ real (kind=kind_phys), dimension(im), intent(in) :: ps, &
& t1, q1, sigmaf, sfcemis, dlwflx, dswsfc, snet, tg3, cm, &
- & ch, prsl1, prslki, ddvel, shdmin, shdmax, &
+ & ch, prsl1, prslki, wind, shdmin, shdmax, &
& snoalb, sfalb, zf, &
& bexppert, xlaipert, vegfpert
@@ -242,7 +241,7 @@ subroutine lsm_noah_run &
! --- locals:
real (kind=kind_phys), dimension(im) :: rch, rho, &
- & q0, qs1, theta1, wind, weasd_old, snwdph_old, &
+ & q0, qs1, theta1, weasd_old, snwdph_old, &
& tprcp_old, srflag_old, tskin_old, canopy_old
real (kind=kind_phys), dimension(km) :: et, sldpth, stsoil, &
@@ -319,9 +318,6 @@ subroutine lsm_noah_run &
do i = 1, im
if (flag_iter(i) .and. land(i)) then
- wind(i) = max(sqrt( u1(i)*u1(i) + v1(i)*v1(i) ) &
- & + max(0.0, min(ddvel(i), 30.0)), 1.0)
-
q0(i) = max(q1(i), 1.e-8) !* q1=specific humidity at level 1 (kg/kg)
theta1(i) = t1(i) * prslki(i) !* adiabatic temp at level 1 (k)
diff --git a/physics/sfc_drv.meta b/physics/sfc_drv.meta
index f628c6c27..7728ee375 100644
--- a/physics/sfc_drv.meta
+++ b/physics/sfc_drv.meta
@@ -165,24 +165,6 @@
kind = kind_phys
intent = in
optional = F
-[u1]
- standard_name = x_wind_at_lowest_model_layer
- long_name = x component of 1st model layer wind
- units = m s-1
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = in
- optional = F
-[v1]
- standard_name = y_wind_at_lowest_model_layer
- long_name = y component of 1st model layer wind
- units = m s-1
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = in
- optional = F
[t1]
standard_name = air_temperature_at_lowest_model_layer
long_name = 1st model layer air temperature
@@ -227,8 +209,8 @@
intent = in
optional = F
[sfcemis]
- standard_name = surface_longwave_emissivity
- long_name = surface longwave emissivity
+ standard_name = surface_longwave_emissivity_over_land_interstitial
+ long_name = surface lw emissivity in fraction over land (temporary use as interstitial)
units = frac
dimensions = (horizontal_dimension)
type = real
@@ -236,8 +218,8 @@
intent = in
optional = F
[dlwflx]
- standard_name = surface_downwelling_longwave_flux_absorbed_by_ground
- long_name = total sky surface downward longwave flux absorbed by the ground
+ standard_name = surface_downwelling_longwave_flux_absorbed_by_ground_over_land
+ long_name = total sky surface downward longwave flux absorbed by the ground over land
units = W m-2
dimensions = (horizontal_dimension)
type = real
@@ -333,9 +315,9 @@
type = logical
intent = in
optional = F
-[ddvel]
- standard_name = surface_wind_enhancement_due_to_convection
- long_name = surface wind enhancement due to convection
+[wind]
+ standard_name = wind_speed_at_lowest_model_layer
+ long_name = wind speed at lowest model level
units = m s-1
dimensions = (horizontal_dimension)
type = real
diff --git a/physics/sfc_drv_ruc.F90 b/physics/sfc_drv_ruc.F90
index 64e4d4597..a16cfc334 100644
--- a/physics/sfc_drv_ruc.F90
+++ b/physics/sfc_drv_ruc.F90
@@ -143,7 +143,7 @@ subroutine lsm_ruc_run & ! inputs
& sfcemis, dlwflx, dswsfc, snet, delt, tg3, cm, ch, &
& prsl1, zf, ddvel, shdmin, shdmax, alvwf, alnwf, &
& snoalb, sfalb, flag_iter, flag_guess, isot, ivegsrc, fice, &
- & smc, stc, slc, lsm_ruc, lsm, land, &
+ & smc, stc, slc, lsm_ruc, lsm, land, islimsk, &
& imp_physics, imp_physics_gfdl, imp_physics_thompson, &
& smcwlt2, smcref2, wspd, do_mynnsfclay, &
& con_cp, con_rv, con_rd, con_g, con_pi, con_hvap, con_fvirt,& ! constants
@@ -184,6 +184,7 @@ subroutine lsm_ruc_run & ! inputs
con_hvap, con_fvirt
logical, dimension(im), intent(in) :: flag_iter, flag_guess, land
+ integer, dimension(im), intent(in) :: islimsk ! sea/land/ice mask (=0/1/2)
logical, intent(in) :: do_mynnsfclay
! --- in/out:
@@ -384,7 +385,7 @@ subroutine lsm_ruc_run & ! inputs
!> - Set flag for land and ice points.
!- 10may19 - ice points are turned off.
flag(i) = land(i)
- if (land(i) .and. (vegtype(i)==iswater .or. vegtype(i)==isice)) then
+ if (land(i) .and. (vegtype(i)==iswater .or. (vegtype(i)==isice.and.islimsk(i)==2))) then
!write(errmsg,'(a,i0,a,i0)') 'Logic error in sfc_drv_ruc_run: for i=', i, &
! ', land(i) is true but vegtype(i) is water or ice: ', vegtype(i)
!errflg = 1
@@ -897,7 +898,7 @@ subroutine lsm_ruc_run & ! inputs
sfcdew(i) = dew(i,j)
qsurf(i) = qsfc(i,j)
sncovr1(i) = sncovr(i,j)
- stm(i) = soilm(i,j) * 1000.0 ! unit conversion (from m to kg m-2)
+ stm(i) = soilm(i,j)
tsurf(i) = soilt(i,j)
tice(i) = tsurf(i)
diff --git a/physics/sfc_drv_ruc.meta b/physics/sfc_drv_ruc.meta
index 8d06e4785..8128a03dd 100644
--- a/physics/sfc_drv_ruc.meta
+++ b/physics/sfc_drv_ruc.meta
@@ -278,6 +278,14 @@
type = logical
intent = in
optional = F
+[islimsk]
+ standard_name = sea_land_ice_mask
+ long_name = sea/land/ice mask (=0/1/2)
+ units = flag
+ dimensions = (horizontal_dimension)
+ type = integer
+ intent = in
+ optional = F
[rainnc]
standard_name = lwe_thickness_of_explicit_rainfall_amount_from_previous_timestep
long_name = explicit rainfall from previous timestep
diff --git a/physics/sfc_noahmp_drv.f b/physics/sfc_noahmp_drv.f
index a089e84d0..ab9f2af0d 100755
--- a/physics/sfc_noahmp_drv.f
+++ b/physics/sfc_noahmp_drv.f
@@ -41,14 +41,19 @@ end subroutine noahmpdrv_finalize
!> \section arg_table_noahmpdrv_run Argument Table
!! \htmlinclude noahmpdrv_run.html
!!
+! !
+! lheatstrg- logical, flag for canopy heat storage 1 !
+! parameterization !
+! !
!-----------------------------------
subroutine noahmpdrv_run &
!...................................
! --- inputs:
& ( im, km, itime, ps, u1, v1, t1, q1, soiltyp, vegtype, &
& sigmaf, sfcemis, dlwflx, dswsfc, snet, delt, tg3, cm, ch, &
- & prsl1, prslki, zf, dry, ddvel, slopetyp, &
+ & prsl1, prslki, zf, dry, wind, slopetyp, &
& shdmin, shdmax, snoalb, sfalb, flag_iter, flag_guess, &
+ & lheatstrg, &
& idveg, iopt_crs, iopt_btr, iopt_run, iopt_sfc, iopt_frz, &
& iopt_inf, iopt_rad, iopt_alb, iopt_snf, iopt_tbot, &
& iopt_stc, xlatin, xcoszin, iyrlen, julian, &
@@ -118,7 +123,7 @@ subroutine noahmpdrv_run &
real (kind=kind_phys), dimension(im), intent(in) :: ps, u1, v1, &
& t1, q1, sigmaf, sfcemis, dlwflx, dswsfc, snet, tg3, cm, &
- & ch, prsl1, prslki, ddvel, shdmin, shdmax, &
+ & ch, prsl1, prslki, wind, shdmin, shdmax, &
& snoalb, sfalb, zf, &
& rainn_mp,rainc_mp,snow_mp,graupel_mp,ice_mp
@@ -136,7 +141,9 @@ subroutine noahmpdrv_run &
real (kind=kind_phys), intent(in) :: delt
logical, dimension(im), intent(in) :: flag_iter, flag_guess
-
+
+ logical, intent(in) :: lheatstrg
+
real (kind=kind_phys), intent(in) :: con_hvap, con_cp, con_jcal, &
& rhoh2o, con_eps, con_epsm1, con_fvirt, &
& con_rd, con_hfus
@@ -178,7 +185,7 @@ subroutine noahmpdrv_run &
! --- locals:
real (kind=kind_phys), dimension(im) :: rch, rho, &
- & q0, qs1, theta1, tv1, wind, weasd_old, snwdph_old, &
+ & q0, qs1, theta1, tv1, weasd_old, snwdph_old, &
& tprcp_old, srflag_old, tskin_old, canopy_old
real (kind=kind_phys), dimension(km) :: et,stsoil,smsoil, slsoil
@@ -236,6 +243,8 @@ subroutine noahmpdrv_run &
& irb,tr,evc,chleaf,chuc,chv2,chb2, &
& fpice,pahv,pahg,pahb,pah,co2pp,o2pp,ch2b
+ real (kind=kind_phys) :: cpfac
+
integer :: i, k, ice, stype, vtype ,slope,nroot,couple
logical :: flag(im)
logical :: snowng,frzgra
@@ -358,10 +367,6 @@ subroutine noahmpdrv_run &
do i = 1, im
if (flag_iter(i) .and. flag(i)) then
- wind(i) = sqrt( u1(i)*u1(i) + v1(i)*v1(i) ) &
- & + max(0.0, min(ddvel(i), 30.0))
- wind(i) = max(wind(i), 1.0)
-
q0(i) = max(q1(i), 1.e-8) !* q1=specific humidity at level 1 (kg/kg)
theta1(i) = t1(i) * prslki(i) !* adiabatic temp at level 1 (k)
@@ -628,6 +633,10 @@ subroutine noahmpdrv_run &
call noahmp_options(idveg ,iopt_crs,iopt_btr,iopt_run,iopt_sfc, &
& iopt_frz,iopt_inf,iopt_rad,iopt_alb,iopt_snf,iopt_tbot,iopt_stc)
+!
+! initialize heat capacity enhancement factor for heat storage parameterization
+!
+ cpfac = 1.0
if ( vtype == isice_table ) then
@@ -716,6 +725,7 @@ subroutine noahmpdrv_run &
& qc , swdn , lwdn ,& ! in : forcing
& pconv , pnonc , pshcv , psnow , pgrpl , phail ,& ! in : forcing
& tbot , co2pp , o2pp , foln , ficeold , zlvl ,& ! in : forcing
+ & lheatstrg ,& ! in : canopy heat storage
& alboldx , sneqvox ,& ! in/out :
& tsnsox , slsoil , smsoil , tahx , eahx , fwetx ,& ! in/out :
& canliqx , canicex , tvx , tgx , qsfc1d , qsnowx ,& ! in/out :
@@ -723,7 +733,7 @@ subroutine noahmpdrv_run &
& zwtx , wax , wtx , wslakex , lfmassx , rtmassx,& ! in/out :
& stmassx , woodx , stblcpx , fastcpx , xlaix ,xsaix ,& ! in/out :
& cmx , chx , taussx ,& ! in/out :
- & smcwtdx ,deeprechx, rechx ,& ! in/out :
+ & smcwtdx ,deeprechx, rechx , cpfac ,& ! in/out :
& z0wrf ,& ! out
& fsa , fsr , fira , fsh , ssoil , fcev ,& ! out :
& fgev , fctr , ecan , etran , edir , trad ,& ! out :
@@ -864,7 +874,7 @@ subroutine noahmpdrv_run &
! ssoil = -1.0 *ssoil
call penman (sfctmp,sfcprs,chx,t2v,th2,prcp,fdown,ssoil, &
- & q2,q2sat,etp,snowng,frzgra,ffrozp,dqsdt2,emissi,fsno)
+ & cpfac,q2,q2sat,etp,snowng,frzgra,ffrozp,dqsdt2,emissi,fsno)
ep(i) = etp
@@ -1126,7 +1136,7 @@ end subroutine transfer_mp_parameters
subroutine penman (sfctmp,sfcprs,ch,t2v,th2,prcp,fdown,ssoil, &
- & q2,q2sat,etp,snowng,frzgra,ffrozp, &
+ & cpfac,q2,q2sat,etp,snowng,frzgra,ffrozp, &
& dqsdt2,emissi_in,sncovr)
! etp is calcuated right after ssoil
@@ -1141,11 +1151,12 @@ subroutine penman (sfctmp,sfcprs,ch,t2v,th2,prcp,fdown,ssoil, &
implicit none
logical, intent(in) :: snowng, frzgra
real, intent(in) :: ch, dqsdt2,fdown,prcp,ffrozp, &
- & q2, q2sat,ssoil, sfcprs, sfctmp, &
+ & q2, q2sat,ssoil,cpfac, sfcprs, sfctmp, &
& t2v, th2,emissi_in,sncovr
real, intent(out) :: etp
real :: epsca,flx2,rch,rr,t24
real :: a, delta, fnet,rad,rho,emissi,elcp1,lvs
+ real :: elcpx
real, parameter :: elcp = 2.4888e+3, lsubc = 2.501000e+6,cp = 1004.6
real, parameter :: lsubs = 2.83e+6, rd = 287.05, cph2o = 4.1855e+3
@@ -1159,11 +1170,12 @@ subroutine penman (sfctmp,sfcprs,ch,t2v,th2,prcp,fdown,ssoil, &
! prepare partial quantities for penman equation.
! ----------------------------------------------------------------------
emissi=emissi_in
-! elcp1 = (1.0-sncovr)*elcp + sncovr*elcp*lsubs/lsubc
+ elcpx = elcp / cpfac
+! elcp1 = (1.0-sncovr)*elcpx + sncovr*elcpx*lsubs/lsubc
lvs = (1.0-sncovr)*lsubc + sncovr*lsubs
flx2 = 0.0
- delta = elcp * dqsdt2
+ delta = elcpx * dqsdt2
! delta = elcp1 * dqsdt2
t24 = sfctmp * sfctmp * sfctmp * sfctmp
rr = t24 * 6.48e-8 / (sfcprs * ch) + 1.0
@@ -1174,7 +1186,7 @@ subroutine penman (sfctmp,sfcprs,ch,t2v,th2,prcp,fdown,ssoil, &
! adjust the partial sums / products with the latent heat
! effects caused by falling precipitation.
! ----------------------------------------------------------------------
- rch = rho * cp * ch
+ rch = rho * cp * cpfac * ch
if (.not. snowng) then
if (prcp > 0.0) rr = rr + cph2o * prcp / rch
else
@@ -1197,7 +1209,7 @@ subroutine penman (sfctmp,sfcprs,ch,t2v,th2,prcp,fdown,ssoil, &
! ----------------------------------------------------------------------
end if
rad = fnet / rch + th2- sfctmp
- a = elcp * (q2sat - q2)
+ a = elcpx * (q2sat - q2)
! a = elcp1 * (q2sat - q2)
epsca = (a * rr + rad * delta) / (delta + rr)
etp = epsca * rch / lsubc
diff --git a/physics/sfc_noahmp_drv.meta b/physics/sfc_noahmp_drv.meta
index 9baa85082..066bc1e87 100644
--- a/physics/sfc_noahmp_drv.meta
+++ b/physics/sfc_noahmp_drv.meta
@@ -150,22 +150,22 @@
intent= in
optional = F
[sfcemis]
- standard_name = surface_longwave_emissivity
- long_name = surface lw emissivity in fraction
+ standard_name = surface_longwave_emissivity_over_land_interstitial
+ long_name = surface lw emissivity in fraction over land (temporary use as interstitial)
units = frac
dimensions = (horizontal_dimension)
type = real
kind = kind_phys
- intent= in
+ intent = in
optional = F
[dlwflx]
- standard_name = surface_downwelling_longwave_flux_absorbed_by_ground
- long_name = total sky surface downward longwave flux absorbed by the ground
+ standard_name = surface_downwelling_longwave_flux_absorbed_by_ground_over_land
+ long_name = total sky surface downward longwave flux absorbed by the ground over land
units = W m-2
dimensions = (horizontal_dimension)
type = real
kind = kind_phys
- intent= in
+ intent = in
optional = F
[dswsfc]
standard_name = surface_downwelling_shortwave_flux
@@ -256,9 +256,9 @@
type = logical
intent = in
optional = F
-[ddvel]
- standard_name = surface_wind_enhancement_due_to_convection
- long_name = surface wind enhancement due to convection
+[wind]
+ standard_name = wind_speed_at_lowest_model_layer
+ long_name = wind speed at lowest model level
units = m s-1
dimensions = (horizontal_dimension)
type = real
@@ -325,6 +325,14 @@
type = logical
intent = in
optional = F
+[lheatstrg]
+ standard_name = flag_for_canopy_heat_storage
+ long_name = flag for canopy heat storage parameterization
+ units = flag
+ dimensions = ()
+ type = logical
+ intent = in
+ optional = F
[idveg]
standard_name = flag_for_dynamic_vegetation_option
long_name = choice for dynamic vegetation option (see noahmp module for definition)
diff --git a/physics/sfc_nst.f b/physics/sfc_nst.f
index da9b8c87c..ed43a719d 100644
--- a/physics/sfc_nst.f
+++ b/physics/sfc_nst.f
@@ -29,19 +29,21 @@ end subroutine sfc_nst_finalize
!! \section NSST_general_algorithm GFS Near-Surface Sea Temperature Scheme General Algorithm
!> @{
subroutine sfc_nst_run &
+! --- inputs:
& ( im, hvap, cp, hfus, jcal, eps, epsm1, rvrdm1, rd, rhw0, &
& pi, sbc, ps, u1, v1, t1, q1, tref, cm, ch, &
- & prsl1, prslki, prsik1, prslk1, wet, icy, xlon, sinlat, &
+ & prsl1, prslki, prsik1, prslk1, wet, xlon, sinlat, &
& stress, &
& sfcemis, dlwflx, sfcnsw, rain, timestep, kdt, solhr,xcosz, &
- & ddvel, flag_iter, flag_guess, nstf_name1, nstf_name4, &
- & nstf_name5, lprnt, ipr, & ! inputs from here and above
+ & wind, flag_iter, flag_guess, nstf_name1, nstf_name4, &
+ & nstf_name5, lprnt, ipr, &
+! --- input/output:
& 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, & ! in/outs from here and above
- & qsurf, gflux, cmm, chh, evap, hflx, ep, errmsg, errflg & ! outputs
+ & z_c, c_0, c_d, w_0, w_d, d_conv, ifd, qrain, &
+! --- outputs:
+ & qsurf, gflux, cmm, chh, evap, hflx, ep, errmsg, errflg &
& )
-
-! DH* 20190718: prslki can be removed if GSD_SURFACE_FLUXES_BUGFIX is adopted
+!
! ===================================================================== !
! description: !
! !
@@ -51,10 +53,9 @@ subroutine sfc_nst_run &
! call sfc_nst !
! inputs: !
! ( im, ps, u1, v1, t1, q1, tref, cm, ch, !
-! prsl1, prslki, prsik1, prslk1, iwet, iice, xlon, sinlat, !
-! stress, !
+! prsl1, prslki, wet, xlon, sinlat, stress, !
! sfcemis, dlwflx, sfcnsw, rain, timestep, kdt,solhr,xcosz, !
-! ddvel, flag_iter, flag_guess, nstf_name1, nstf_name4, !
+! wind, flag_iter, flag_guess, nstf_name1, nstf_name4, !
! nstf_name5, lprnt, ipr, !
! input/outputs: !
! tskin, tsurf, xt, xs, xu, xv, xz, zm, xtts, xzts, dt_cool, !
@@ -106,17 +107,12 @@ subroutine sfc_nst_run &
! sfcemis - real, sfc lw emissivity (fraction) im !
! dlwflx - real, total sky sfc downward lw flux (w/m**2) im !
! sfcnsw - real, total sky sfc netsw flx into ocean (w/m**2) im !
-! DH*
-! The actual unit of rain passed in is m ! see below line 438, qrain(i) = ...
-! where 1000*rain in the nominator converts m to kg m^2; there is still a
-! time unit 's' missing. Need to double-check what is going on.
-! *DH
! rain - real, rainfall rate (kg/m**2/s) im !
! timestep - real, timestep interval (second) 1 !
! kdt - integer, time step counter 1 !
! solhr - real, fcst hour at the end of prev time step 1 !
! xcosz - real, consine of solar zenith angle 1 !
-! ddvel - real, wind enhancement due to convection (m/s) im !
+! wind - real, wind speed (m/s) im !
! flag_iter- logical, execution or not im !
! when iter = 1, flag_iter = .true. for all grids im !
! when iter = 2, flag_iter = .true. when wind < 2 im !
@@ -197,12 +193,12 @@ subroutine sfc_nst_run &
real (kind=kind_phys), dimension(im), intent(in) :: ps, u1, v1, &
& t1, q1, tref, cm, ch, prsl1, prslki, prsik1, prslk1, &
& xlon,xcosz, &
- & sinlat, stress, sfcemis, dlwflx, sfcnsw, rain, ddvel
+ & sinlat, stress, sfcemis, dlwflx, sfcnsw, rain, wind
real (kind=kind_phys), intent(in) :: timestep
real (kind=kind_phys), intent(in) :: solhr
- logical, dimension(im), intent(in) :: flag_iter, flag_guess, wet, &
- & icy
+ logical, dimension(im), intent(in) :: flag_iter, flag_guess, wet
+! &, icy
logical, intent(in) :: lprnt
! --- input/outputs:
@@ -224,7 +220,7 @@ subroutine sfc_nst_run &
integer :: k,i
!
real (kind=kind_phys), dimension(im) :: q0, qss, rch,
- & rho_a, theta1, tv1, wind, wndmag
+ & rho_a, theta1, tv1, wndmag
real(kind=kind_phys) elocp,tem,cpinv,hvapi
!
@@ -265,13 +261,15 @@ subroutine sfc_nst_run &
! flag for open water and where the iteration is on
!
do i = 1, im
- flag(i) = wet(i) .and. .not.icy(i) .and. flag_iter(i)
+! flag(i) = wet(i) .and. .not.icy(i) .and. flag_iter(i)
+ flag(i) = wet(i) .and. flag_iter(i)
enddo
!
! save nst-related prognostic fields for guess run
!
do i=1, im
- if(wet(i) .and. .not.icy(i) .and. flag_guess(i)) then
+! if(wet(i) .and. .not.icy(i) .and. flag_guess(i)) then
+ if(wet(i) .and. flag_guess(i)) then
xt_old(i) = xt(i)
xs_old(i) = xs(i)
xu_old(i) = xu(i)
@@ -298,8 +296,6 @@ subroutine sfc_nst_run &
nswsfc(i) = sfcnsw(i) ! net solar radiation at the air-sea surface (positive downward)
wndmag(i) = sqrt(u1(i)*u1(i) + v1(i)*v1(i))
- wind(i) = wndmag(i) + max( 0.0, min( ddvel(i), 30.0 ) )
- wind(i) = max( wind(i), 1.0 )
q0(i) = max(q1(i), 1.0e-8)
#ifdef GSD_SURFACE_FLUXES_BUGFIX
@@ -588,8 +584,9 @@ subroutine sfc_nst_run &
! restore nst-related prognostic fields for guess run
do i=1, im
- if(wet(i) .and. .not.icy(i)) then
- if(flag_guess(i)) then ! when it is guess of
+! if (wet(i) .and. .not.icy(i)) then
+ if (wet(i)) then
+ if (flag_guess(i)) then ! when it is guess of
xt(i) = xt_old(i)
xs(i) = xs_old(i)
xu(i) = xu_old(i)
@@ -609,9 +606,9 @@ subroutine sfc_nst_run &
!
if ( nstf_name1 > 1 ) then
tskin(i) = tsurf(i)
- endif ! if nstf_name1 > 1
- endif ! if flag_guess(i)
- endif ! if wet(i) .and. .not.icy(i)
+ endif ! if nstf_name1 > 1 then
+ endif ! if flag_guess(i) then
+ endif ! if wet(i) .and. .not.icy(i) then
enddo
! if (lprnt .and. i == ipr) print *,' beg xz8=',xz(i)
@@ -678,11 +675,8 @@ end subroutine sfc_nst_pre_finalize
!> \section NSST_general_pre_algorithm General Algorithm
!! @{
subroutine sfc_nst_pre_run
- & (im, rlapse, icy, wet, zorl_ocn, zorl_ice, cd_ocn, cd_ice,
- & cdq_ocn, cdq_ice, rb_ocn, rb_ice, stress_ocn, stress_ice,
- & ffmm_ocn, ffmm_ice, ffhh_ocn, ffhh_ice, uustar_ocn,
- & uustar_ice, fm10_ocn, fm10_ice, fh2_ocn, fh2_ice, oro,
- & oro_uf, tsfc_ocn, tsurf_ocn, tseal, errmsg, errflg)
+ & (im, wet, tsfc_ocn, tsurf_ocn, tseal, xt, xz, dt_cool,
+ & z_c, tref, cplflx, errmsg, errflg)
use machine , only : kind_phys
@@ -690,16 +684,14 @@ subroutine sfc_nst_pre_run
! --- inputs:
integer, intent(in) :: im
- logical, dimension(im), intent(in) :: icy, wet
- real (kind=kind_phys), intent(in) :: rlapse
- real (kind=kind_phys), dimension(im), intent(in) :: zorl_ice,
- & cd_ice, cdq_ice, rb_ice, stress_ice, ffmm_ice, ffhh_ice,
- & uustar_ice, fm10_ice, fh2_ice, oro, oro_uf, tsfc_ocn
+ logical, dimension(im), intent(in) :: wet
+ real (kind=kind_phys), dimension(im), intent(in) ::
+ & tsfc_ocn, xt, xz, dt_cool, z_c
+ logical, intent(in) :: cplflx
! --- input/outputs:
- real (kind=kind_phys), dimension(im), intent(inout) :: tsurf_ocn,
- & zorl_ocn, cd_ocn, cdq_ocn, rb_ocn, stress_ocn, ffmm_ocn,
- & ffhh_ocn, uustar_ocn, fm10_ocn, fh2_ocn, tseal
+ real (kind=kind_phys), dimension(im), intent(inout) ::
+ & tsurf_ocn, tseal, tref
! --- outputs:
character(len=*), intent(out) :: errmsg
@@ -707,20 +699,48 @@ subroutine sfc_nst_pre_run
! --- locals
integer :: i
- real(kind=kind_phys) :: tem
+ real(kind=kind_phys), parameter :: zero = 0.0d0,
+ & one = 1.0d0,
+ & half = 0.5d0,
+ & omz1 = 10.0d0
+ real(kind=kind_phys) :: tem1, tem2, dt_warm
! Initialize CCPP error handling variables
errmsg = ''
errflg = 0
do i=1,im
- if (wet(i) .and. .not. icy(i)) then
- tem = (oro(i)-oro_uf(i)) * rlapse
- tseal(i) = tsfc_ocn(i) + tem
- tsurf_ocn(i) = tsurf_ocn(i) + tem
+ if (wet(i)) then
+! tem = (oro(i)-oro_uf(i)) * rlapse
+ ! DH* 20190927 simplyfing this code because tem is zero
+ !tem = zero
+ !tseal(i) = tsfc_ocn(i) + tem
+ tseal(i) = tsfc_ocn(i)
+ !tsurf_ocn(i) = tsurf_ocn(i) + tem
+ ! *DH
endif
enddo
+ if (cplflx) then
+ tem1 = half / omz1
+ do i=1,im
+ if (wet(i)) then
+ tem2 = one / xz(i)
+ dt_warm = (xt(i)+xt(i)) * tem2
+ if ( xz(i) > omz1) then
+ tref(i) = tseal(i) - (one-half*omz1*tem2) * dt_warm &
+ & + z_c(i)*dt_cool(i)*tem1
+ else
+ tref(i) = tseal(i) - (xz(i)*dt_warm &
+ & - z_c(i)*dt_cool(i))*tem1
+ endif
+ tseal(i) = tref(i) + dt_warm - dt_cool(i)
+! - (Sfcprop%oro(i)-Sfcprop%oro_uf(i))*rlapse
+ tsurf_ocn(i) = tseal(i)
+ endif
+ enddo
+ endif
+
return
end subroutine sfc_nst_pre_run
!! @}
@@ -799,11 +819,11 @@ subroutine sfc_nst_post_run &
! & ' dt_cool=',dt_cool(ipr),' dt_warm=',2.0*xt(ipr)/xz(ipr),
! & ' kdt=',kdt
- do i = 1, im
- if (wet(i) .and. .not. icy(i)) then
- tsurf_ocn(i) = tsurf_ocn(i) - (oro(i)-oro_uf(i)) * rlapse
- endif
- enddo
+! do i = 1, im
+! if (wet(i) .and. .not. icy(i)) then
+! tsurf_ocn(i) = tsurf_ocn(i) - (oro(i)-oro_uf(i)) * rlapse
+! endif
+! enddo
! --- ... run nsst model ... ---
@@ -812,12 +832,15 @@ subroutine sfc_nst_post_run &
zsea1 = 0.001*real(nstf_name4)
zsea2 = 0.001*real(nstf_name5)
call get_dtzm_2d (xt, xz, dt_cool, &
- & z_c, wet, icy, zsea1, zsea2, &
+ & z_c, wet, zsea1, zsea2, &
& im, 1, dtzm)
do i = 1, im
- if ( wet(i) .and. .not. icy(i) ) then
- tsfc_ocn(i) = max(271.2,tref(i) + dtzm(i)) - &
- & (oro(i)-oro_uf(i))*rlapse
+! if (wet(i) .and. .not.icy(i)) then
+! if (wet(i) .and. (Model%frac_grid .or. .not. icy(i))) then
+ if (wet(i)) then
+ tsfc_ocn(i) = max(271.2, tref(i) + dtzm(i))
+! tsfc_ocn(i) = max(271.2, tref(i) + dtzm(i)) - &
+! (oro(i)-oro_uf(i))*rlapse
endif
enddo
endif
diff --git a/physics/sfc_nst.meta b/physics/sfc_nst.meta
index 73b585c71..d74f68c0e 100644
--- a/physics/sfc_nst.meta
+++ b/physics/sfc_nst.meta
@@ -234,14 +234,6 @@
type = logical
intent = in
optional = F
-[icy]
- standard_name = flag_nonzero_sea_ice_surface_fraction
- long_name = flag indicating presence of some sea ice surface area fraction
- units = flag
- dimensions = (horizontal_dimension)
- type = logical
- intent = in
- optional = F
[xlon]
standard_name = longitude
long_name = longitude
@@ -270,8 +262,8 @@
intent = in
optional = F
[sfcemis]
- standard_name = surface_longwave_emissivity
- long_name = surface longwave emissivity
+ standard_name = surface_longwave_emissivity_over_ocean_interstitial
+ long_name = surface lw emissivity in fraction over ocean (temporary use as interstitial)
units = frac
dimensions = (horizontal_dimension)
type = real
@@ -279,8 +271,8 @@
intent = in
optional = F
[dlwflx]
- standard_name = surface_downwelling_longwave_flux_absorbed_by_ground
- long_name = total sky sfc downward lw flux absorbed by the ocean
+ standard_name = surface_downwelling_longwave_flux_absorbed_by_ground_over_ocean
+ long_name = total sky surface downward longwave flux absorbed by the ground over ocean
units = W m-2
dimensions = (horizontal_dimension)
type = real
@@ -340,9 +332,9 @@
kind = kind_phys
intent = in
optional = F
-[ddvel]
- standard_name = surface_wind_enhancement_due_to_convection
- long_name = wind enhancement due to convection
+[wind]
+ standard_name = wind_speed_at_lowest_model_layer
+ long_name = wind speed at lowest model level
units = m s-1
dimensions = (horizontal_dimension)
type = real
@@ -679,23 +671,6 @@
type = integer
intent = in
optional = F
-[rlapse]
- standard_name = air_temperature_lapse_rate_constant
- long_name = environmental air temperature lapse rate constant
- units = K m-1
- dimensions = ()
- type = real
- kind = kind_phys
- intent = in
- optional = F
-[icy]
- standard_name = flag_nonzero_sea_ice_surface_fraction
- long_name = flag indicating presence of some sea ice surface area fraction
- units = flag
- dimensions = (horizontal_dimension)
- type = logical
- intent = in
- optional = F
[wet]
standard_name = flag_nonzero_wet_surface_fraction
long_name = flag indicating presence of some ocean or lake surface area fraction
@@ -704,230 +679,85 @@
type = logical
intent = in
optional = F
-[zorl_ocn]
- standard_name = surface_roughness_length_over_ocean_interstitial
- long_name = surface roughness length over ocean (temporary use as interstitial)
- units = cm
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = inout
- optional = F
-[zorl_ice]
- standard_name = surface_roughness_length_over_ice_interstitial
- long_name = surface roughness length over ice (temporary use as interstitial)
- units = cm
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = in
- optional = F
-[cd_ocn]
- standard_name = surface_drag_coefficient_for_momentum_in_air_over_ocean
- long_name = surface exchange coeff for momentum over ocean
- units = none
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = inout
- optional = F
-[cd_ice]
- standard_name = surface_drag_coefficient_for_momentum_in_air_over_ice
- long_name = surface exchange coeff for momentum over ice
- units = none
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = in
- optional = F
-[cdq_ocn]
- standard_name = surface_drag_coefficient_for_heat_and_moisture_in_air_over_ocean
- long_name = surface exchange coeff heat & moisture over ocean
- units = none
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = inout
- optional = F
-[cdq_ice]
- standard_name = surface_drag_coefficient_for_heat_and_moisture_in_air_over_ice
- long_name = surface exchange coeff heat & moisture over ice
- units = none
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = in
- optional = F
-[rb_ocn]
- standard_name = bulk_richardson_number_at_lowest_model_level_over_ocean
- long_name = bulk Richardson number at the surface over ocean
- units = none
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = inout
- optional = F
-[rb_ice]
- standard_name = bulk_richardson_number_at_lowest_model_level_over_ice
- long_name = bulk Richardson number at the surface over ice
- units = none
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = in
- optional = F
-[stress_ocn]
- standard_name = surface_wind_stress_over_ocean
- long_name = surface wind stress over ocean
- units = m2 s-2
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = inout
- optional = F
-[stress_ice]
- standard_name = surface_wind_stress_over_ice
- long_name = surface wind stress over ice
- units = m2 s-2
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = in
- optional = F
-[ffmm_ocn]
- standard_name = Monin_Obukhov_similarity_function_for_momentum_over_ocean
- long_name = Monin-Obukhov similarity function for momentum over ocean
- units = none
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = inout
- optional = F
-[ffmm_ice]
- standard_name = Monin_Obukhov_similarity_function_for_momentum_over_ice
- long_name = Monin-Obukhov similarity function for momentum over ice
- units = none
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = in
- optional = F
-[ffhh_ocn]
- standard_name = Monin_Obukhov_similarity_function_for_heat_over_ocean
- long_name = Monin-Obukhov similarity function for heat over ocean
- units = none
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = inout
- optional = F
-[ffhh_ice]
- standard_name = Monin_Obukhov_similarity_function_for_heat_over_ice
- long_name = Monin-Obukhov similarity function for heat over ice
- units = none
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = in
- optional = F
-[uustar_ocn]
- standard_name = surface_friction_velocity_over_ocean
- long_name = surface friction velocity over ocean
- units = m s-1
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = inout
- optional = F
-[uustar_ice]
- standard_name = surface_friction_velocity_over_ice
- long_name = surface friction velocity over ice
- units = m s-1
+[tsfc_ocn]
+ standard_name = surface_skin_temperature_over_ocean_interstitial
+ long_name = surface skin temperature over ocean (temporary use as interstitial)
+ units = K
dimensions = (horizontal_dimension)
type = real
kind = kind_phys
intent = in
optional = F
-[fm10_ocn]
- standard_name = Monin_Obukhov_similarity_function_for_momentum_at_10m_over_ocean
- long_name = Monin-Obukhov similarity parameter for momentum at 10m over ocean
- units = none
+[tsurf_ocn]
+ standard_name = surface_skin_temperature_after_iteration_over_ocean
+ long_name = surface skin temperature after iteration over ocean
+ units = K
dimensions = (horizontal_dimension)
type = real
kind = kind_phys
intent = inout
optional = F
-[fm10_ice]
- standard_name = Monin_Obukhov_similarity_function_for_momentum_at_10m_over_ice
- long_name = Monin-Obukhov similarity parameter for momentum at 10m over ice
- units = none
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = in
- optional = F
-[fh2_ocn]
- standard_name = Monin_Obukhov_similarity_function_for_heat_at_2m_over_ocean
- long_name = Monin-Obukhov similarity parameter for heat at 2m over ocean
- units = none
+[tseal]
+ standard_name = surface_skin_temperature_for_nsst
+ long_name = ocean surface skin temperature
+ units = K
dimensions = (horizontal_dimension)
type = real
kind = kind_phys
intent = inout
optional = F
-[fh2_ice]
- standard_name = Monin_Obukhov_similarity_function_for_heat_at_2m_over_ice
- long_name = Monin-Obukhov similarity parameter for heat at 2m over ice
- units = none
+[xt]
+ standard_name = diurnal_thermocline_layer_heat_content
+ long_name = heat content in diurnal thermocline layer
+ units = K m
dimensions = (horizontal_dimension)
type = real
kind = kind_phys
intent = in
optional = F
-[oro]
- standard_name = orography
- long_name = orography
+[xz]
+ standard_name = diurnal_thermocline_layer_thickness
+ long_name = diurnal thermocline layer thickness
units = m
dimensions = (horizontal_dimension)
type = real
kind = kind_phys
intent = in
optional = F
-[oro_uf]
- standard_name = orography_unfiltered
- long_name = unfiltered orographyo
- units = m
+[dt_cool]
+ standard_name = sub_layer_cooling_amount
+ long_name = sub-layer cooling amount
+ units = K
dimensions = (horizontal_dimension)
type = real
kind = kind_phys
intent = in
optional = F
-[tsfc_ocn]
- standard_name = surface_skin_temperature_over_ocean_interstitial
- long_name = surface skin temperature over ocean (temporary use as interstitial)
- units = K
+[z_c]
+ standard_name = sub_layer_cooling_thickness
+ long_name = sub-layer cooling thickness
+ units = m
dimensions = (horizontal_dimension)
type = real
kind = kind_phys
intent = in
optional = F
-[tsurf_ocn]
- standard_name = surface_skin_temperature_after_iteration_over_ocean
- long_name = surface skin temperature after iteration over ocean
+[tref]
+ standard_name = sea_surface_reference_temperature
+ long_name = reference/foundation temperature
units = K
dimensions = (horizontal_dimension)
type = real
kind = kind_phys
intent = inout
optional = F
-[tseal]
- standard_name = surface_skin_temperature_for_nsst
- long_name = ocean surface skin temperature
- units = K
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = inout
+[cplflx]
+ standard_name = flag_for_flux_coupling
+ long_name = flag controlling cplflx collection (default off)
+ units = flag
+ dimensions = ()
+ type = logical
+ intent = in
optional = F
[errmsg]
standard_name = ccpp_error_message
diff --git a/physics/sfc_ocean.F b/physics/sfc_ocean.F
index 625e8e5f0..9635f30b8 100644
--- a/physics/sfc_ocean.F
+++ b/physics/sfc_ocean.F
@@ -23,8 +23,8 @@ end subroutine sfc_ocean_finalize
subroutine sfc_ocean_run &
!...................................
! --- inputs:
- & ( im, cp, rd, eps, epsm1, hvap, rvrdm1, ps, u1, v1, t1, q1, &
- & tskin, cm, ch, prsl1, prslki, wet, ddvel, &
+ & ( im, cp, rd, eps, epsm1, hvap, rvrdm1, ps, t1, q1, &
+ & tskin, cm, ch, prsl1, prslki, wet, wind, &
& flag_iter, &
! --- outputs:
& qsurf, cmm, chh, gflux, evap, hflx, ep, &
@@ -38,8 +38,9 @@ subroutine sfc_ocean_run &
! !
! call sfc_ocean !
! inputs: !
-! ( im, ps, u1, v1, t1, q1, tskin, cm, ch, !
-! prsl1, prslki, wet, ddvel, flag_iter, !
+! ( im, ps, t1, q1, tskin, cm, ch, !
+!! ( im, ps, u1, v1, t1, q1, tskin, cm, ch, !
+! prsl1, prslki, wet, wind, flag_iter, !
! outputs: !
! qsurf, cmm, chh, gflux, evap, hflx, ep ) !
! !
@@ -62,7 +63,6 @@ subroutine sfc_ocean_run &
! inputs: size !
! im - integer, horizontal dimension 1 !
! ps - real, surface pressure im !
-! u1, v1 - real, u/v component of surface layer wind im !
! t1 - real, surface layer mean temperature ( k ) im !
! q1 - real, surface layer mean specific humidity im !
! tskin - real, ground surface skin temperature ( k ) im !
@@ -71,7 +71,7 @@ subroutine sfc_ocean_run &
! prsl1 - real, surface layer mean pressure im !
! prslki - real, im !
! wet - logical, =T if any ocean/lak, =F otherwise im !
-! ddvel - real, wind enhancement due to convection (m/s) im !
+! wind - real, wind speed (m/s) im !
! flag_iter- logical, im !
! !
! outputs: !
@@ -95,8 +95,8 @@ subroutine sfc_ocean_run &
real (kind=kind_phys), intent(in) :: cp, rd, eps, epsm1, hvap, &
& rvrdm1
- real (kind=kind_phys), dimension(im), intent(in) :: ps, u1, v1, &
- & t1, q1, tskin, cm, ch, prsl1, prslki, ddvel
+ real (kind=kind_phys), dimension(im), intent(in) :: ps, &
+ & t1, q1, tskin, cm, ch, prsl1, prslki, wind
logical, dimension(im), intent(in) :: flag_iter, wet
@@ -109,7 +109,7 @@ subroutine sfc_ocean_run &
! --- locals:
- real (kind=kind_phys) :: q0, qss, rch, rho, wind, tem, cpinv, &
+ real (kind=kind_phys) :: q0, qss, rch, rho, tem, cpinv, &
& hvapi, elocp
integer :: i
@@ -134,10 +134,6 @@ subroutine sfc_ocean_run &
! rho is density, qss is sat. hum. at surface
if ( flag(i) ) then
-
- wind = max(sqrt(u1(i)*u1(i) + v1(i)*v1(i)) &
- & + max( 0.0, min( ddvel(i), 30.0 ) ), 1.0)
-
q0 = max( q1(i), 1.0e-8 )
rho = prsl1(i) / (rd*t1(i)*(1.0 + rvrdm1*q0))
@@ -151,9 +147,9 @@ subroutine sfc_ocean_run &
! --- ... rcp = rho cp ch v
- rch = rho * cp * ch(i) * wind
- cmm(i) = cm(i) * wind
- chh(i) = rho * ch(i) * wind
+ rch = rho * cp * ch(i) * wind(i)
+ cmm(i) = cm(i) * wind(i)
+ chh(i) = rho * ch(i) * wind(i)
! --- ... sensible and latent heat flux over open water
diff --git a/physics/sfc_ocean.meta b/physics/sfc_ocean.meta
index 4304e344d..d60c1ce2c 100644
--- a/physics/sfc_ocean.meta
+++ b/physics/sfc_ocean.meta
@@ -82,24 +82,6 @@
kind = kind_phys
intent = in
optional = F
-[u1]
- standard_name = x_wind_at_lowest_model_layer
- long_name = x component of surface layer wind
- units = m s-1
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = in
- optional = F
-[v1]
- standard_name = y_wind_at_lowest_model_layer
- long_name = y component of surface layer wind
- units = m s-1
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = in
- optional = F
[t1]
standard_name = air_temperature_at_lowest_model_layer
long_name = surface layer mean temperature
@@ -171,9 +153,9 @@
type = logical
intent = in
optional = F
-[ddvel]
- standard_name = surface_wind_enhancement_due_to_convection
- long_name = wind enhancement due to convection
+[wind]
+ standard_name = wind_speed_at_lowest_model_layer
+ long_name = wind speed at lowest model level
units = m s-1
dimensions = (horizontal_dimension)
type = real
diff --git a/physics/sfc_sice.f b/physics/sfc_sice.f
index 7c2da2415..9471792fa 100644
--- a/physics/sfc_sice.f
+++ b/physics/sfc_sice.f
@@ -41,17 +41,16 @@ end subroutine sfc_sice_finalize
!> @{
subroutine sfc_sice_run &
& ( im, km, sbc, hvap, tgice, cp, eps, epsm1, rvrdm1, grav, & ! --- inputs:
- & t0c, rd, cimin, ps, u1, v1, t1, q1, delt, &
+ & t0c, rd, ps, t1, q1, delt, &
& sfcemis, dlwflx, sfcnsw, sfcdsw, srflag, &
- & cm, ch, prsl1, prslki, prsik1, prslk1, islimsk, ddvel, &
- & flag_iter, lprnt, ipr, &
+ & cm, ch, prsl1, prslki, prsik1, prslk1, islimsk, wind, &
+ & flag_iter, lprnt, ipr, cimin, &
& hice, fice, tice, weasd, tskin, tprcp, stc, ep, & ! --- input/outputs:
& snwdph, qsurf, snowmt, gflux, cmm, chh, evap, hflx, & !
& cplflx, cplchm, flag_cice, islmsk_cice, &
& errmsg, errflg
& )
-! DH* 20190718: prslki can be removed if GSD_SURFACE_FLUXES_BUGFIX is adopted
! ===================================================================== !
! description: !
! !
@@ -59,9 +58,9 @@ subroutine sfc_sice_run &
! !
! call sfc_sice !
! inputs: !
-! ( im, km, ps, u1, v1, t1, q1, delt, !
+! ( im, km, ps, t1, q1, delt, !
! sfcemis, dlwflx, sfcnsw, sfcdsw, srflag, !
-! cm, ch, prsl1, prslki, prsik1, prslk1, islimsk, ddvel, !
+! cm, ch, prsl1, prslki, prsik1, prslk1, islimsk, wind, !
! flag_iter, !
! input/outputs: !
! hice, fice, tice, weasd, tskin, tprcp, stc, ep, !
@@ -93,7 +92,6 @@ subroutine sfc_sice_run &
! inputs: size !
! im, km - integer, horiz dimension and num of soil layers 1 !
! ps - real, surface pressure im !
-! u1, v1 - real, u/v component of surface layer wind im !
! t1 - real, surface layer mean temperature ( k ) im !
! q1 - real, surface layer mean specific humidity im !
! delt - real, time interval (second) 1 !
@@ -109,7 +107,7 @@ subroutine sfc_sice_run &
! prsik1 - real, im !
! prslk1 - real, im !
! islimsk - integer, sea/land/ice mask (=0/1/2) im !
-! ddvel - real, im !
+! wind - real, im !
! flag_iter- logical, im !
! !
! input/outputs: !
@@ -134,7 +132,7 @@ subroutine sfc_sice_run &
! !
! ===================================================================== !
!
- use machine, only: kind_phys
+ use machine, only : kind_phys
use funcphys, only : fpvs
!
implicit none
@@ -156,15 +154,15 @@ subroutine sfc_sice_run &
logical, intent(in) :: cplchm
real (kind=kind_phys), intent(in) :: sbc, hvap, tgice, cp, eps, &
- & epsm1, grav, rvrdm1, t0c, rd, cimin
+ & epsm1, grav, rvrdm1, t0c, rd
- real (kind=kind_phys), dimension(im), intent(in) :: ps, u1, v1, &
+ real (kind=kind_phys), dimension(im), intent(in) :: ps, &
& t1, q1, sfcemis, dlwflx, sfcnsw, sfcdsw, srflag, cm, ch, &
- & prsl1, prslki, prsik1, prslk1, ddvel
+ & prsl1, prslki, prsik1, prslk1, wind
integer, dimension(im), intent(in) :: islimsk
integer, dimension(im), intent(in) :: islmsk_cice
- real (kind=kind_phys), intent(in) :: delt
+ real (kind=kind_phys), intent(in) :: delt, cimin
logical, dimension(im), intent(in) :: flag_iter, flag_cice
@@ -189,7 +187,7 @@ subroutine sfc_sice_run &
& snowd, theta1
real (kind=kind_phys) :: t12, t14, tem, stsice(im,kmi)
- &, hflxi, hflxw, q0, qs1, wind, qssi, qssw
+ &, hflxi, hflxw, q0, qs1, qssi, qssw
real (kind=kind_phys) :: cpinv, hvapi, elocp
integer :: i, k
@@ -266,9 +264,6 @@ subroutine sfc_sice_run &
! dlwflx has been given a negative sign for downward longwave
! sfcnsw is the net shortwave flux (direction: dn-up)
- wind = max(sqrt(u1(i)*u1(i) + v1(i)*v1(i)) &
- & + max(zero, min(ddvel(i), 30.0d0)), one)
-
q0 = max(q1(i), 1.0e-8)
! tsurf(i) = tskin(i)
#ifdef GSD_SURFACE_FLUXES_BUGFIX
@@ -307,8 +302,8 @@ subroutine sfc_sice_run &
! --- ... rcp = rho cp ch v
- cmm(i) = cm(i) * wind
- chh(i) = rho(i) * ch(i) * wind
+ cmm(i) = cm(i) * wind(i)
+ chh(i) = rho(i) * ch(i) * wind(i)
rch(i) = chh(i) * cp
!> - Calculate sensible and latent heat flux over open water & sea ice.
diff --git a/physics/sfc_sice.meta b/physics/sfc_sice.meta
index 1af043885..c9641ffaa 100644
--- a/physics/sfc_sice.meta
+++ b/physics/sfc_sice.meta
@@ -107,15 +107,6 @@
kind = kind_phys
intent = in
optional = F
-[cimin]
- standard_name = minimum_sea_ice_concentration
- long_name = minimum sea ice concentration
- units = frac
- dimensions = ()
- type = real
- kind = kind_phys
- intent = in
- optional = F
[ps]
standard_name = surface_air_pressure
long_name = surface pressure
@@ -125,24 +116,6 @@
kind = kind_phys
intent = in
optional = F
-[u1]
- standard_name = x_wind_at_lowest_model_layer
- long_name = u component of surface layer wind
- units = m s-1
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = in
- optional = F
-[v1]
- standard_name = y_wind_at_lowest_model_layer
- long_name = v component of surface layer wind
- units = m s-1
- dimensions = (horizontal_dimension)
- type = real
- kind = kind_phys
- intent = in
- optional = F
[t1]
standard_name = air_temperature_at_lowest_model_layer
long_name = surface layer mean temperature
@@ -171,8 +144,8 @@
intent = in
optional = F
[sfcemis]
- standard_name = surface_longwave_emissivity
- long_name = sfc lw emissivity
+ standard_name = surface_longwave_emissivity_over_ice_interstitial
+ long_name = surface lw emissivity in fraction over ice (temporary use as interstitial)
units = frac
dimensions = (horizontal_dimension)
type = real
@@ -180,8 +153,8 @@
intent = in
optional = F
[dlwflx]
- standard_name = surface_downwelling_longwave_flux_absorbed_by_ground
- long_name = total sky surface downward longwave flux absorbed by the ground
+ standard_name = surface_downwelling_longwave_flux_absorbed_by_ground_over_ice
+ long_name = total sky surface downward longwave flux absorbed by the ground over ice
units = W m-2
dimensions = (horizontal_dimension)
type = real
@@ -277,9 +250,9 @@
type = integer
intent = in
optional = F
-[ddvel]
- standard_name = surface_wind_enhancement_due_to_convection
- long_name = wind enhancement due to convection
+[wind]
+ standard_name = wind_speed_at_lowest_model_layer
+ long_name = wind speed at lowest model level
units = m s-1
dimensions = (horizontal_dimension)
type = real
@@ -310,6 +283,15 @@
type = integer
intent = in
optional = F
+[cimin]
+ standard_name = lake_ice_minimum
+ long_name = minimum lake ice value
+ units = ???
+ dimensions = ()
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
[hice]
standard_name = sea_ice_thickness
long_name = sea-ice thickness
diff --git a/physics/sfcsub.F b/physics/sfcsub.F
index 7039884f8..6296e7856 100644
--- a/physics/sfcsub.F
+++ b/physics/sfcsub.F
@@ -6146,17 +6146,24 @@ subroutine setrmsk(kpds5,slmask,igaul,jgaul,wlon,rnlat, &
!
ijmax = imax*jmax
rslmsk = 0.
+! TG3 MODS BEGIN
+ if(kpds5 == kpdtsf .and. imax == 138 .and. jmax == 116
+ & .and. kpds4 == 128) then
+! print*,'turn off setrmsk for tg3'
+ lmask = .false.
+
+ elseif(kpds5 == kpdtsf) then
+! TG3 MODS END
!
! surface temperature
!
- if(kpds5.eq.kpdtsf) then
-! lmask=.false.
+ lmask = .false.
call ga2la(slmask,igaul,jgaul,rslmsk,imax,jmax,wlon,rnlat
&, rlnout, rltout, gaus, blno, blto)
! &, dlon, dlat, gaus, blno, blto)
- crit=0.5
+ crit = 0.5
call rof01(rslmsk,ijmax,'ge',crit)
- lmask=.true.
+ lmask = .true.
!
! bucket soil wetness
!
@@ -6164,16 +6171,16 @@ subroutine setrmsk(kpds5,slmask,igaul,jgaul,wlon,rnlat, &
call ga2la(slmask,igaul,jgaul,rslmsk,imax,jmax,wlon,rnlat
&, rlnout, rltout, gaus, blno, blto)
! &, dlon, dlat, gaus, blno, blto)
- crit=0.5
+ crit = 0.5
call rof01(rslmsk,ijmax,'ge',crit)
- lmask=.true.
+ lmask = .true.
! write(6,*) 'wet rslmsk'
! znnt=1.
! call nntprt(rslmsk,ijmax,znnt)
!
! snow depth
!
- elseif(kpds5.eq.kpdsnd) then
+ elseif(kpds5 == kpdsnd) then
if(kpds4 == 192) then ! use the bitmap
rslmsk = 0.
do j = 1, jmax
@@ -7043,51 +7050,51 @@ subroutine clima(lugb,iy,im,id,ih,fh,len,lsoil, &
!
! get tsf climatology for the begining of the forecast
!
- if (fh .gt. 0.0) then
+ if (fh > 0.0) then
!cbosu
if (me == 0) print*,'bosu fh gt 0'
- iy4=iy
- if(iy.lt.101) iy4=1900+iy4
- fha=0
- ida=0
- jda=0
-! fha(2)=nint(fh)
- ida(1)=iy
- ida(2)=im
- ida(3)=id
- ida(5)=ih
+ iy4 = iy
+ if (iy < 101) iy4 = 1900 + iy4
+ fha = 0
+ ida = 0
+ jda = 0
+! fha(2) = nint(fh)
+ ida(1) = iy
+ ida(2) = im
+ ida(3) = id
+ ida(5) = ih
call w3kind(w3kindreal,w3kindint)
if(w3kindreal == 4) then
- fha4=fha
+ fha4 = fha
call w3movdat(fha4,ida,jda)
else
call w3movdat(fha,ida,jda)
endif
- jy=jda(1)
- jm=jda(2)
- jd=jda(3)
- jh=jda(5)
- if (me .eq. 0) write(6,*) ' forecast jy,jm,jd,jh',
- & jy,jm,jd,jh
+ jy = jda(1)
+ jm = jda(2)
+ jd = jda(3)
+ jh = jda(5)
+ if (me == 0) write(6,*) ' forecast jy,jm,jd,jh',
+ & jy,jm,jd,jh
jdow = 0
jdoy = 0
jday = 0
call w3doxdat(jda,jdow,jdoy,jday)
- rjday=jdoy+jda(5)/24.
- if(rjday.lt.dayhf(1)) rjday=rjday+365.
+ rjday = jdoy + jda(5) / 24.
+ if(rjday < dayhf(1)) rjday = rjday + 365.
!
- if (me .eq. 0) write(6,*) 'forecast jy,jm,jd,jh=',jy,jm,jd,jh
+ if (me == 0) write(6,*) 'forecast jy,jm,jd,jh=',jy,jm,jd,jh
!
! for monthly mean climatology
!
monend = 12
do mm=1,monend
- mmm=mm
- mmp=mm+1
- if(rjday.ge.dayhf(mmm).and.rjday.lt.dayhf(mmp)) then
- mon1=mmm
- mon2=mmp
+ mmm = mm
+ mmp = mm + 1
+ if(rjday >= dayhf(mmm) .and. rjday < dayhf(mmp)) then
+ mon1 = mmm
+ mon2 = mmp
go to 10
endif
enddo
@@ -7095,17 +7102,18 @@ subroutine clima(lugb,iy,im,id,ih,fh,len,lsoil, &
call abort
10 continue
wei1m = (dayhf(mon2)-rjday)/(dayhf(mon2)-dayhf(mon1))
- wei2m = (rjday-dayhf(mon1))/(dayhf(mon2)-dayhf(mon1))
- if(mon2.eq.13) mon2=1
- if (me .eq. 0) print *,'rjday,mon1,mon2,wei1m,wei2m=',
- & rjday,mon1,mon2,wei1m,wei2m
+ wei2m = 1.0 - wei1m
+! wei2m = (rjday-dayhf(mon1))/(dayhf(mon2)-dayhf(mon1))
+ if (mon2 == 13) mon2 = 1
+ if (me == 0) print *,'rjday,mon1,mon2,wei1m,wei2m=',
+ & rjday,mon1,mon2,wei1m,wei2m
!
! read monthly mean climatology of tsf
!
kpd7 = -1
do nn=1,2
mon = mon1
- if (nn .eq. 2) mon = mon2
+ if (nn == 2) mon = mon2
call fixrdc(lugb,fntsfc,kpdtsf,kpd7,mon,slmask,
& tsf(1,nn),len,iret
&, imsk, jmsk, slmskh, gaus,blno, blto
@@ -7122,8 +7130,8 @@ subroutine clima(lugb,iy,im,id,ih,fh,len,lsoil, &
!
! compute current jy,jm,jd,jh of forecast and the day of the year
!
- iy4=iy
- if(iy.lt.101) iy4=1900+iy4
+ iy4 = iy
+ if (iy < 101) iy4=1900+iy4
fha = 0
ida = 0
jda = 0
@@ -7133,8 +7141,8 @@ subroutine clima(lugb,iy,im,id,ih,fh,len,lsoil, &
ida(3) = id
ida(5) = ih
call w3kind(w3kindreal,w3kindint)
- if(w3kindreal==4) then
- fha4=fha
+ if(w3kindreal == 4) then
+ fha4 = fha
call w3movdat(fha4,ida,jda)
else
call w3movdat(fha,ida,jda)
@@ -7149,44 +7157,45 @@ subroutine clima(lugb,iy,im,id,ih,fh,len,lsoil, &
jdoy = 0
jday = 0
call w3doxdat(jda,jdow,jdoy,jday)
- rjday = jdoy+jda(5)/24.
- if(rjday.lt.dayhf(1)) rjday=rjday+365.
+ rjday = jdoy + jda(5) / 24.
+ if(rjday < dayhf(1)) rjday = rjday + 365.
- if (me .eq. 0) write(6,*) ' forecast jy,jm,jd,jh,rjday=',
- & jy,jm,jd,jh,rjday
+ if (me == 0) write(6,*) ' forecast jy,jm,jd,jh,rjday=',
+ & jy,jm,jd,jh,rjday
!
- if (me .eq. 0) write(6,*) 'forecast jy,jm,jd,jh=',jy,jm,jd,jh
+ if (me == 0) write(6,*) 'forecast jy,jm,jd,jh=',jy,jm,jd,jh
!
! for monthly mean climatology
!
monend = 12
do mm=1,monend
- mmm=mm
- mmp=mm+1
- if(rjday.ge.dayhf(mmm).and.rjday.lt.dayhf(mmp)) then
- mon1=mmm
- mon2=mmp
+ mmm = mm
+ mmp = mm + 1
+ if(rjday >= dayhf(mmm) .and. rjday < dayhf(mmp)) then
+ mon1 = mmm
+ mon2 = mmp
go to 20
endif
enddo
print *,'wrong rjday',rjday
call abort
20 continue
- wei1m=(dayhf(mon2)-rjday)/(dayhf(mon2)-dayhf(mon1))
- wei2m=(rjday-dayhf(mon1))/(dayhf(mon2)-dayhf(mon1))
- if(mon2.eq.13) mon2=1
- if (me .eq. 0) print *,'rjday,mon1,mon2,wei1m,wei2m=',
- & rjday,mon1,mon2,wei1m,wei2m
+ wei1m = (dayhf(mon2)-rjday)/(dayhf(mon2)-dayhf(mon1))
+ wei2m = 1.0 - wei1m
+! wei2m = (rjday-dayhf(mon1))/(dayhf(mon2)-dayhf(mon1))
+ if (mon2 == 13) mon2 = 1
+ if (me == 0) print *,'rjday,mon1,mon2,wei1m,wei2m=',
+ & rjday,mon1,mon2,wei1m,wei2m
!
! for seasonal mean climatology
!
monend = 4
is = im/3 + 1
- if (is.eq.5) is = 1
+ if (is == 5) is = 1
do mm=1,monend
mmm = mm*3 - 2
mmp = (mm+1)*3 - 2
- if(rjday.ge.dayhf(mmm).and.rjday.lt.dayhf(mmp)) then
+ if(rjday >= dayhf(mmm) .and. rjday < dayhf(mmp)) then
sea1 = mmm
sea2 = mmp
go to 30
@@ -7196,20 +7205,21 @@ subroutine clima(lugb,iy,im,id,ih,fh,len,lsoil, &
call abort
30 continue
wei1s = (dayhf(sea2)-rjday)/(dayhf(sea2)-dayhf(sea1))
- wei2s = (rjday-dayhf(sea1))/(dayhf(sea2)-dayhf(sea1))
- if(sea2.eq.13) sea2=1
- if (me .eq. 0) print *,'rjday,sea1,sea2,wei1s,wei2s=',
- & rjday,sea1,sea2,wei1s,wei2s
+ wei2s = 1.0 - wei1s
+! wei2s = (rjday-dayhf(sea1))/(dayhf(sea2)-dayhf(sea1))
+ if (sea2 == 13) sea2 = 1
+ if (me == 0) print *,'rjday,sea1,sea2,wei1s,wei2s=',
+ & rjday,sea1,sea2,wei1s,wei2s
!
! for summer and winter values (maximum and minimum).
!
monend = 2
is = im/6 + 1
- if (is.eq.3) is = 1
+ if (is == 3) is = 1
do mm=1,monend
mmm = mm*6 - 5
mmp = (mm+1)*6 - 5
- if(rjday.ge.dayhf(mmm).and.rjday.lt.dayhf(mmp)) then
+ if(rjday >= dayhf(mmm) .and. rjday < dayhf(mmp)) then
hyr1 = mmm
hyr2 = mmp
go to 31
@@ -7219,10 +7229,11 @@ subroutine clima(lugb,iy,im,id,ih,fh,len,lsoil, &
call abort
31 continue
wei1y = (dayhf(hyr2)-rjday)/(dayhf(hyr2)-dayhf(hyr1))
- wei2y = (rjday-dayhf(hyr1))/(dayhf(hyr2)-dayhf(hyr1))
- if(hyr2.eq.13) hyr2=1
- if (me .eq. 0) print *,'rjday,hyr1,hyr2,wei1y,wei2y=',
- & rjday,hyr1,hyr2,wei1y,wei2y
+ wei2y = 1.0 - wei1y
+! wei2y = (rjday-dayhf(hyr1))/(dayhf(hyr2)-dayhf(hyr1))
+ if (hyr2 == 13) hyr2 = 1
+ if (me == 0) print *,'rjday,hyr1,hyr2,wei1y,wei2y=',
+ & rjday,hyr1,hyr2,wei1y,wei2y
!
! start reading in climatology and interpolate to the date
!
@@ -7622,7 +7633,7 @@ subroutine clima(lugb,iy,im,id,ih,fh,len,lsoil, &
!
mon1s = mon1 ; mon2s = mon2 ; sea1s = sea1 ; sea2s = sea2
!
- if (me .eq. 0) print *,' mon1s=',mon1s,' mon2s=',mon2s
+ if (me == 0) print *,' mon1s=',mon1s,' mon2s=',mon2s
&,' sea1s=',sea1s,' sea2s=',sea2s
!
k1 = 1 ; k2 = 2
@@ -7680,11 +7691,11 @@ subroutine clima(lugb,iy,im,id,ih,fh,len,lsoil, &
! seasonal mean climatology
!
isx = sea2/3 + 1
- if (isx .eq. 5) isx = 1
- if(isx.eq.1) kpd9 = 12
- if(isx.eq.2) kpd9 = 3
- if(isx.eq.3) kpd9 = 6
- if(isx.eq.4) kpd9 = 9
+ if (isx == 5) isx = 1
+ if (isx == 1) kpd9 = 12
+ if (isx == 2) kpd9 = 3
+ if (isx == 3) kpd9 = 6
+ if (isx == 4) kpd9 = 9
!
! albedo
! there are four albedo fields in this version:
@@ -7720,7 +7731,7 @@ subroutine clima(lugb,iy,im,id,ih,fh,len,lsoil, &
if (me == 0) print*,'bosu 2nd time in clima for month ',
& mon, k1,k2
if ( index(fnalbc, "tileX.nc") == 0) then ! grib file
- kpd7=-1
+ kpd7 = -1
do k = 1, 4
call fixrdc(lugb,fnalbc,kpdalb(k),kpd7,mon,slmask,
& alb(1,k,nn),len,iret
@@ -7737,7 +7748,7 @@ subroutine clima(lugb,iy,im,id,ih,fh,len,lsoil, &
!
! tsf at the current time t
!
- kpd7=-1
+ kpd7 = -1
call fixrdc(lugb,fntsfc,kpdtsf,kpd7,mon,slmask,
& tsf(1,nn),len,iret
&, imsk, jmsk, slmskh, gaus,blno, blto
@@ -7745,13 +7756,13 @@ subroutine clima(lugb,iy,im,id,ih,fh,len,lsoil, &
!
! soil wetness
!
- if(fnwetc(1:8).ne.' ') then
+ if (fnwetc(1:8).ne.' ') then
kpd7=-1
call fixrdc(lugb,fnwetc,kpdwet,kpd7,mon,slmask,
& wet(1,nn),len,iret
&, imsk, jmsk, slmskh, gaus,blno, blto
&, outlat, outlon, me)
- elseif(fnsmcc(1:8).ne.' ') then
+ elseif (fnsmcc(1:8).ne.' ') then
if (index(fnsmcc,'global_soilmcpc.1x1.grb') /= 0) then ! the old climo data
kpd7=-1
call fixrdc(lugb,fnsmcc,kpdsmc,kpd7,mon,slmask,
@@ -7793,13 +7804,13 @@ subroutine clima(lugb,iy,im,id,ih,fh,len,lsoil, &
!
! sea ice
!
- kpd7=-1
- if(fnacnc(1:8).ne.' ') then
+ kpd7 = -1
+ if (fnacnc(1:8).ne.' ') then
call fixrdc(lugb,fnacnc,kpdacn,kpd7,mon,slmask,
& acn(1,nn),len,iret
&, imsk, jmsk, slmskh, gaus,blno, blto
&, outlat, outlon, me)
- elseif(fnaisc(1:8).ne.' ') then
+ elseif (fnaisc(1:8).ne.' ') then
call fixrdc(lugb,fnaisc,kpdais,kpd7,mon,slmask,
& ais(1,nn),len,iret
&, imsk, jmsk, slmskh, gaus,blno, blto
@@ -7819,7 +7830,7 @@ subroutine clima(lugb,iy,im,id,ih,fh,len,lsoil, &
!
! snow cover
!
- if(fnscvc(1:8).ne.' ') then
+ if (fnscvc(1:8).ne.' ') then
kpd7=-1
call fixrdc(lugb,fnscvc,kpdscv,kpd7,mon,slmask,
& scv(1,nn),len,iret
@@ -7830,7 +7841,7 @@ subroutine clima(lugb,iy,im,id,ih,fh,len,lsoil, &
!
! surface roughness
!
- if(fnzorc(1:3) == 'sib') then
+ if (fnzorc(1:3) == 'sib') then
if (me == 0) then
write(6,*) 'roughness length to be set from sib veg type'
endif
@@ -7848,7 +7859,7 @@ subroutine clima(lugb,iy,im,id,ih,fh,len,lsoil, &
!
! vegetation cover
!
- if(fnvegc(1:8).ne.' ') then
+ if (fnvegc(1:8) .ne. ' ') then
if ( index(fnvegc, "tileX.nc") == 0) then ! grib file
kpd7=-1
call fixrdc(lugb,fnvegc,kpdveg,kpd7,mon,slmask,
@@ -7870,35 +7881,35 @@ subroutine clima(lugb,iy,im,id,ih,fh,len,lsoil, &
! when chosen, set the z0 based on the vegetation type.
! for this option to work, namelist variable fnvetc must be
! set to point at the proper vegetation type file.
- if(fnzorc(1:3) == 'sib') then
- if(fnvetc(1:4) == ' ') then
+ if (fnzorc(1:3) == 'sib') then
+ if (fnvetc(1:4) == ' ') then
if (me==0) write(6,*) "must choose sib veg type climo file"
call abort
endif
zorclm = 0.0
do i=1,len
- ivtyp=nint(vet(i))
+ ivtyp = nint(vet(i))
if (ivtyp >= 1 .and. ivtyp <= 13) then
zorclm(i) = z0_sib(ivtyp)
endif
enddo
elseif(fnzorc(1:4) == 'igbp') then
- if(fnvetc(1:4) == ' ') then
- if (me==0) write(6,*) "must choose igbp veg type climo file"
+ if (fnvetc(1:4) == ' ') then
+ if (me == 0) write(6,*) "must choose igbp veg type climo file"
call abort
endif
zorclm = 0.0
do i=1,len
- ivtyp=nint(vet(i))
+ ivtyp = nint(vet(i))
if (ivtyp >= 1 .and. ivtyp <= 20) then
z0_season(1) = z0_igbp_min(ivtyp)
z0_season(7) = z0_igbp_max(ivtyp)
- if(outlat(i) < 0.0)then
+ if (outlat(i) < 0.0) then
zorclm(i) = wei1y * z0_season(hyr2) +
- & wei2y *z0_season(hyr1)
+ & wei2y * z0_season(hyr1)
else
zorclm(i) = wei1y * z0_season(hyr1) +
- & wei2y *z0_season(hyr2)
+ & wei2y * z0_season(hyr2)
endif
endif
enddo
diff --git a/physics/sflx.f b/physics/sflx.f
index 5c0cf08ce..1654a8872 100644
--- a/physics/sflx.f
+++ b/physics/sflx.f
@@ -337,7 +337,8 @@ subroutine gfssflx &! --- input
& psisat, quartz, rch, refkdt, rr, rgl, rsmax, sndens, &
& sncond, sbeta, sn_new, slope, snup, salp, soilwm, soilww, &
& t1v, t24, t2v, th2v, topt, tsnow, zbot, z0
-
+
+ real (kind=kind_phys) :: shdfac0
real (kind=kind_phys), dimension(nsold) :: rtdis, zsoil
logical :: frzgra, snowng
@@ -368,6 +369,7 @@ subroutine gfssflx &! --- input
! vegetation fraction (shdfac) = 0.
!> - Set ice = -1 and green vegetation fraction (shdfac) = 0 for glacial-ice land.
+ shdfac0 = shdfac
ice = icein
if(ivegsrc == 2) then
@@ -420,12 +422,18 @@ subroutine gfssflx &! --- input
!only igbp type has urban
!urban
if(vegtyp == 13)then
- shdfac=0.05
- rsmin=400.0
- smcmax = 0.45
- smcref = 0.42
- smcwlt = 0.40
- smcdry = 0.40
+! shdfac=0.05
+! rsmin=400.0
+! smcmax = 0.45
+! smcref = 0.42
+! smcwlt = 0.40
+! smcdry = 0.40
+ rsmin=400.0*(1-shdfac0)+40.0*shdfac0 ! gvf
+ shdfac=shdfac0 ! gvf
+ smcmax = 0.45*(1-shdfac0)+smcmax*shdfac0
+ smcref = 0.42*(1-shdfac0)+smcref*shdfac0
+ smcwlt = 0.40*(1-shdfac0)+smcwlt*shdfac0
+ smcdry = 0.40*(1-shdfac0)+smcdry*shdfac0
endif
endif
@@ -662,18 +670,21 @@ subroutine gfssflx &! --- input
! --- outputs:
& df1 &
& )
-!> - For IGBP/urban, \f$df1=3.24\f$.
- if(ivegsrc == 1) then
+! if(ivegsrc == 1) then
!only igbp type has urban
!urban
- if ( vegtyp == 13 ) df1=3.24
- endif
+! if ( vegtyp == 13 ) df1=3.24
+! endif
!> - Add subsurface heat flux reduction effect from the
!! overlying green canopy, adapted from section 2.1.2 of
!! \cite peters-lidard_et_al_1997.
-
- df1 = df1 * exp( sbeta*shdfac )
+!wz only urban for igbp type
+ if(ivegsrc == 1 .and. vegtyp == 13) then
+ df1 = 3.24*(1.-shdfac) + shdfac*df1*exp(sbeta*shdfac)
+ else
+ df1 = df1 * exp( sbeta*shdfac )
+ endif
endif ! end if_ice_block
@@ -1499,18 +1510,22 @@ subroutine nopac
! --- outputs:
& df1 &
& )
- if(ivegsrc == 1) then
+! if(ivegsrc == 1) then
!urban
- if ( vegtyp == 13 ) df1=3.24
- endif
+! if ( vegtyp == 13 ) df1=3.24
+! endif
! --- ... vegetation greenness fraction reduction in subsurface heat
! flux via reduction factor, which is convenient to apply here
! to thermal diffusivity that is later used in hrt to compute
! sub sfc heat flux (see additional comments on veg effect
! sub-sfc heat flx in routine sflx)
-
- df1 = df1 * exp( sbeta*shdfac )
+!wz only urban for igbp type
+ if(ivegsrc == 1 .and. vegtyp == 13) then
+ df1 = 3.24*(1.-shdfac) + shdfac*df1*exp(sbeta*shdfac)
+ else
+ df1 = df1 * exp( sbeta*shdfac )
+ endif
! --- ... compute intermediate terms passed to routine hrt (via routine
! shflx below) for use in computing subsurface heat flux in hrt
@@ -2595,8 +2610,8 @@ subroutine snopac
if (t12 <= tfreez) then
t1 = t12
-! ssoil = df1 * (t1 - stc(1)) / dtot
- ssoil = (t1 - stc (1)) * max(7.0, df1/dtot)
+ ssoil = df1 * (t1 - stc(1)) / dtot
+!wz ssoil = (t1 - stc (1)) * max(7.0, df1/dtot)
sneqv = max(0.0, sneqv-esnow2)
flx3 = 0.0
ex = 0.0
@@ -2729,7 +2744,7 @@ subroutine snopac
! skin temp value as revised by shflx.
zz1 = 1.0
- yy = stc(1) - 0.5*ssoil*zsoil(1)*zz1 / df1
+ yy = stc(1) - 0.5*ssoil*zsoil(1)*zz1 / df1
t11 = t1
! --- ... shflx will calc/update the soil temps. note: the sub-sfc heat flux
@@ -3371,6 +3386,7 @@ subroutine shflx &
! --- inputs:
& ( nsoil, stc, smc, smcmax, zsoil, yy, zz1, tbot, &
& zbot, psisat, dt, bexp, df1, quartz, csoil,vegtyp, &
+ & shdfac, &
! --- input/outputs:
& sh2o, &
! --- outputs:
@@ -4037,6 +4053,7 @@ subroutine hrt &
! --- inputs:
& ( nsoil, stc, smc, smcmax, zsoil, yy, zz1, tbot, &
& zbot, psisat, dt, bexp, df1, quartz, csoil, vegtyp, &
+ & shdfac, &
! --- input/outputs:
& sh2o, &
! --- outputs:
@@ -4090,7 +4107,7 @@ subroutine hrt &
real (kind=kind_phys), intent(in) :: stc(nsoil), smc(nsoil), &
& smcmax, zsoil(nsoil), yy, zz1, tbot, zbot, psisat, dt, &
- & bexp, df1, quartz, csoil
+ & bexp, df1, quartz, csoil, shdfac
! --- input/outputs:
real (kind=kind_phys), intent(inout) :: sh2o(nsoil)
@@ -4116,7 +4133,8 @@ subroutine hrt &
if (ivegsrc == 1)then
!urban
if( vegtyp == 13 ) then
- csoil_loc=3.0e6
+! csoil_loc=3.0e6
+ csoil_loc=3.0e6*(1.-shdfac)+csoil*shdfac ! gvf
endif
endif
@@ -4206,7 +4224,7 @@ subroutine hrt &
call snksrc &
! --- inputs:
& ( nsoil, 1, tavg, smc(1), smcmax, psisat, bexp, dt, &
- & qtot, zsoil, &
+ & qtot, zsoil, shdfac, &
! --- input/outputs:
& sh2o(1), &
! --- outputs:
@@ -4248,9 +4266,13 @@ subroutine hrt &
& df1n &
& )
!urban
- if (ivegsrc == 1)then
- if ( vegtyp == 13 ) df1n = 3.24
- endif
+! if (ivegsrc == 1)then
+! if ( vegtyp == 13 ) df1n = 3.24
+! endif
+!wz only urban for igbp type
+ if(ivegsrc == 1 .and. vegtyp == 13) then
+ df1n = 3.24*(1.-shdfac) + shdfac*df1n
+ endif
! --- ... calc the vertical soil temp gradient thru this layer
@@ -4288,9 +4310,13 @@ subroutine hrt &
& df1n &
& )
!urban
- if (ivegsrc == 1)then
- if ( vegtyp == 13 ) df1n = 3.24
- endif
+! if (ivegsrc == 1)then
+! if ( vegtyp == 13 ) df1n = 3.24
+! endif
+!wz only urban for igbp type
+ if(ivegsrc == 1 .and. vegtyp == 13) then
+ df1n = 3.24*(1.-shdfac) + shdfac*df1n
+ endif
! --- ... calc the vertical soil temp gradient thru bottom layer.
@@ -4344,7 +4370,7 @@ subroutine hrt &
call snksrc &
! --- inputs:
& ( nsoil, k, tavg, smc(k), smcmax, psisat, bexp, dt, &
- & qtot, zsoil, &
+ & qtot, zsoil, shdfac, &
! --- input/outputs:
& sh2o(k), &
! --- outputs:
@@ -4759,7 +4785,7 @@ end subroutine rosr12
subroutine snksrc &
! --- inputs:
& ( nsoil, k, tavg, smc, smcmax, psisat, bexp, dt, &
- & qtot, zsoil, &
+ & qtot, zsoil, shdfac, &
! --- input/outputs:
& sh2o, &
! --- outputs:
@@ -4804,7 +4830,7 @@ subroutine snksrc &
integer, intent(in) :: nsoil, k
real (kind=kind_phys), intent(in) :: tavg, smc, smcmax, psisat, &
- & bexp, dt, qtot, zsoil(nsoil)
+ & bexp, dt, qtot, zsoil(nsoil), shdfac
! --- input/outputs:
real (kind=kind_phys), intent(inout) :: sh2o
@@ -4819,9 +4845,13 @@ subroutine snksrc &
! real (kind=kind_phys) :: frh2o
!urban
- if (ivegsrc == 1)then
- if ( vegtyp == 13 ) df1=3.24
- endif
+! if (ivegsrc == 1)then
+! if ( vegtyp == 13 ) df1=3.24
+! endif
+!wz only urban for igbp type
+ if(ivegsrc == 1 .and. vegtyp == 13) then
+ df1 = 3.24*(1.-shdfac) + shdfac*df1
+ endif
!
!===> ... begin here
!
diff --git a/physics/ugwp_driver_v0.f b/physics/ugwp_driver_v0.F
similarity index 74%
rename from physics/ugwp_driver_v0.f
rename to physics/ugwp_driver_v0.F
index a3ca5f96d..52375dd18 100644
--- a/physics/ugwp_driver_v0.f
+++ b/physics/ugwp_driver_v0.F
@@ -11,65 +11,76 @@ module sso_coorde
end module sso_coorde
!
!
+! Routine cires_ugwp_driver_v0 is replaced with cires_ugwp.F90/cires_ugwp_run in CCPP
+#if 0
subroutine cires_ugwp_driver_v0(me, master,
- & im, levs, nmtvr, dtp, kdt, imx,do_tofd,
+ & im, levs, nmtvr, dtp, kdt, imx, do_ugwp, do_tofd,
& cdmbgwd, xlat, xlatd, sinlat, coslat, spgrid,
& ugrs, vgrs, tgrs, qgrs, prsi, prsl, prslk,
- & phii, phil, del, oro_stat, sgh30, kpbl,
+ & phii, phil, del, hprime, oc, oa4, clx, theta,
+ & gamm, sigma, elvmax, sgh30, kpbl,
& dusfcg, dvsfcg, gw_dudt, gw_dvdt, gw_dtdt, gw_kdis,
& tau_tofd, tau_mtb, tau_ogw, tau_ngw,
- & zmtb, zlwb, zogw, du3dt_mtb,du3dt_ogw, du3dt_tms,rdxzb )
+ & zmtb, zlwb, zogw, du3dt_mtb,du3dt_ogw, du3dt_tms,rdxzb,
+ & rain, ntke, tke, lprnt, ipr)
!-----------------------------------------------------------
-! Part 1 "old-revised" gfs-gwdps_v0
+! Part 1 "old-revised" gfs-gwdps_v0 or "old" gwdps (if do_ugwp=.false.)
! Part 2 non-stationary multi-wave GWs FV3GFS-v0
! Part 3 Dissipative version of UGWP-tendency application
! (similar to WAM-2017)
!-----------------------------------------------------------
- use machine, only: kind_phys
-! use physcons, only: con_cp, con_fvirt, con_g, con_rd,
-! & con_rv, con_rerth, con_pi
+ use machine, only : kind_phys
+ use physcons, only : con_cp, con_g, con_rd, con_rv
- use ugwp_wmsdis_init, only : tamp_mpa
+ use ugwp_wmsdis_init, only : tamp_mpa, ilaunch
use sso_coorde, only : pgwd, pgwd4
implicit none
!input
integer, intent(in) :: me, master
- integer, intent(in) :: im, levs, nmtvr, kdt, imx
+ integer, intent(in) :: im, levs, kdt, imx, nmtvr, ntke, ipr
- real(kind=kind_phys), intent(in) :: dtp, cdmbgwd(2)
- logical :: do_tofd
+ real(kind=kind_phys), intent(in) :: dtp, cdmbgwd(4)
+ logical :: do_ugwp, do_tofd, lprnt
integer, intent(in) :: kpbl(im)
real(kind=kind_phys), intent(in), dimension(im) :: xlat, xlatd
&, sgh30, sinlat, coslat, spgrid ! spgrid = tile-area
+ &, rain
real(kind=kind_phys), intent(in), dimension(im,levs) :: ugrs
&, vgrs, tgrs, qgrs, prsi, prsl, prslk, phii, phil, del
- real(kind=kind_phys), intent(in) :: oro_stat(im,nmtvr)
+! real(kind=kind_phys), intent(in) :: oro_stat(im,nmtvr)
+ real(kind=kind_phys), intent(in), dimension(im) :: hprime, oc
+ &, theta, gamm, sigma, elvmax
+ real(kind=kind_phys), intent(in), dimension(im,4) :: oa4, clx
+ real(kind=kind_phys), intent(in) :: tke(im,levs)
!out
real(kind=kind_phys), dimension(im,levs) :: gw_dudt, gw_dvdt
&, gw_dTdt, gw_kdis
!-----locals + diagnostics output
- real(kind=kind_phys), dimension(im,levs) :: Pdvdt, Pdudt
+ real(kind=kind_phys), dimension(im,levs) :: Pdvdt, Pdudt
&, Pdtdt, Pkdis, ed_dudt, ed_dvdt, ed_dTdt
- real(kind=kind_phys), dimension(im) :: dusfcg, dvsfcg
+ real(kind=kind_phys), dimension(im) :: dusfcg, dvsfcg
+
+ real(kind=kind_phys), dimension(im) :: rdxzb, zmtb,
+ & zlwb, zogw, tau_mtb, tau_ogw, tau_tofd, tau_ngw, turb_fac
+ real(kind=kind_phys), dimension(im,levs) :: du3dt_mtb, du3dt_ogw
+ &, du3dt_tms
+ real(kind=kind_phys), dimension(im) :: tem
- real(kind=kind_phys), dimension(im) :: rdxzb, zmtb,
- & zlwb, zogw, tau_mtb, tau_ogw, tau_tofd, tau_ngw
- real(kind=kind_phys), dimension(im,levs) :: du3dt_mtb, du3dt_ogw
- &, du3dt_tms
! locals
- integer :: i, j, k, ix
+ real(kind=kind_phys) :: rfac, tx1
+ integer :: i, j, k, ix
!
! define hprime, oc, oa4, clx, theta, sigma, gamm, elvmax
!
- real(kind=kind_phys), dimension(im) :: hprime,
- & oc, theta, sigma, gamm, elvmax
- real(kind=kind_phys), dimension(im, 4) :: clx, oa4
+! real(kind=kind_phys), dimension(im) :: hprime,
+! & oc, theta, sigma, gamm, elvmax
+! real(kind=kind_phys), dimension(im, 4) :: clx, oa4
!
! switches that activate impact of OGWs and NGWs along with eddy diffusion
!
@@ -80,87 +91,129 @@ subroutine cires_ugwp_driver_v0(me, master,
!
if (me == master .and. kdt < 2) then
print *
- write(6,*) 'FV3GFS execute ugwp_driver_v0 nmtvr=', nmtvr
+ write(6,*) 'FV3GFS execute ugwp_driver_v0 '
+! write(6,*) 'FV3GFS execute ugwp_driver_v0 nmtvr=', nmtvr
write(6,*) ' COORDE EXPER pogw = ' , pogw
write(6,*) ' COORDE EXPER pgwd = ' , pgwd
write(6,*) ' COORDE EXPER pgwd4 = ', pgwd4
print *
endif
-
-! print *, ' NMTVR in driver ', nmtvr
do i=1,im
- hprime(i) = oro_stat(i,1)
- oc(i) = oro_stat(i,2)
- oa4(i,1) = oro_stat(i,3)
- oa4(i,2) = oro_stat(i,4)
- oa4(i,3) = oro_stat(i,5)
- oa4(i,4) = oro_stat(i,6)
- clx(i,1) = oro_stat(i,7)
- clx(i,2) = oro_stat(i,8)
- clx(i,3) = oro_stat(i,9)
- clx(i,4) = oro_stat(i,10)
- theta(i) = oro_stat(i,11)
- gamm(i) = oro_stat(i,12)
- sigma(i) = oro_stat(i,13)
- elvmax(i) = oro_stat(i,14)
-
- zlwb(i) = 0.
+ zlwb(i) = 0.
enddo
!
! 1) ORO stationary GWs
-!
-! pdvdt(:,:) = 0. ; pdudt(:,:) = 0.
-! pkdis(:,:) = 0. ; pdtdt(:,:) = 0.
-! zlwb(:) = 0.
+! ------------------
- CALL GWDPS_V0(IM, levs, imx, do_tofd,
- & Pdvdt, Pdudt, Pdtdt, Pkdis,
- & ugrs, vgrs, tgrs, qgrs,KPBL, prsi,del,prsl,
- & prslk, phii, phil, DTP,KDT,
- & sgh30, HPRIME,OC,OA4, CLX, THETA,SIGMA,GAMM,ELVMAX,
- & DUSFCg, DVSFCg, xlatd, sinlat, coslat, spgrid,
- & cdmbgwd, me, master, rdxzb,
- & zmtb, zogw, tau_mtb, tau_ogw, tau_tofd,
- & du3dt_mtb, du3dt_ogw, du3dt_tms)
-!
-!
-! non-stationary GW-scheme with GMAO/MERRA GW-forcing
+ if (do_ugwp .and. nmtvr == 14) then ! calling revised old GFS gravity wave drag
+ CALL GWDPS_V0(IM, levs, imx, do_tofd,
+ & Pdvdt, Pdudt, Pdtdt, Pkdis,
+ & ugrs , vgrs, tgrs, qgrs,KPBL, prsi,del,prsl,
+ & prslk, phii, phil, DTP,KDT,
+ & sgh30, HPRIME, OC, OA4, CLX, THETA,
+ & SIGMA, GAMM, ELVMAX,
+ & DUSFCg, DVSFCg, xlatd, sinlat, coslat, spgrid,
+ & cdmbgwd(1:2), me, master, rdxzb,
+ & zmtb, zogw, tau_mtb, tau_ogw, tau_tofd,
+ & du3dt_mtb, du3dt_ogw, du3dt_tms)
+!
+ if (me == master .and. kdt < 2) then
+ print *
+ write(6,*) 'FV3GFS finished gwdps_v0 in ugwp_driver_v0 '
+ print *
+ endif
+ else ! calling old GFS gravity wave drag as is
+ do k=1,levs
+ do i=1,im
+ pdvdt(i,k) = 0.0
+ pdudt(i,k) = 0.0
+ pdtdt(i,k) = 0.0
+ pkdis(i,k) = 0.0
+ enddo
+ enddo
+ if (cdmbgwd(1) > 0.0 .or. cdmbgwd(2) > 0.0) then
+ call gwdps(im, im, im, levs, Pdvdt, Pdudt, Pdtdt &
+ &, ugrs, vgrs, tgrs, qgrs &
+ &, kpbl, prsi, del, prsl, prslk, phii, phil, dtp, kdt&
+ &, hprime, oc, oa4, clx, theta, sigma, gamm &
+ &, elvmax, dusfcg, dvsfcg &
+ &, con_g, con_cp, con_rd, con_rv, imx &
+ &, nmtvr, cdmbgwd(1:2), me, lprnt, ipr, rdxzb)
+ endif
+
+ tau_mtb = 0.0 ; tau_ogw = 0.0 ; tau_tofd = 0.0
+ du3dt_mtb = 0.0 ; du3dt_ogw = 0.0 ; du3dt_tms= 0.0
+ endif
!
- if (me == master .and. kdt < 2) then
- print *
- write(6,*) 'FV3GFS finished gwdps_v0 in ugwp_driver_v0 '
- print *
- endif
+ if (cdmbgwd(3) > 0.0) then
+! 2) non-stationary GWs with GEOS-5/MERRA GW-forcing
+! ----------------------------------------------
!--------
! GMAO GEOS-5/MERRA GW-forcing lat-dep
!--------
- call slat_geos5_tamp(im, tamp_mpa, xlatd, tau_ngw)
+ call slat_geos5_tamp(im, tamp_mpa, xlatd, tau_ngw)
-! call slat_geos5(im, xlatd, tau_ngw)
+! call slat_geos5(im, xlatd, tau_ngw)
!
-! 2) non-stationary GWs with GEOS-5/MERRA GW-forcing
+ if (abs(1.0-cdmbgwd(3)) > 1.0e-6) then
+ if (cdmbgwd(4) > 0.0) then
+ do i=1,im
+ turb_fac(i) = 0.0
+ enddo
+ if (ntke > 0) then
+ do k=1,(levs+levs)/3
+ do i=1,im
+ turb_fac(i) = turb_fac(i) + del(i,k) * tke(i,k)
+ tem(i) = tem(i) + del(i,k)
+ enddo
+ enddo
+ do i=1,im
+ turb_fac(i) = turb_fac(i) / tem(i)
+ enddo
+ endif
+ rfac = 86400000 / dtp
+ do i=1,im
+ tx1 = cdmbgwd(4)*min(10.0, max(turb_fac(i),rain(i)*rfac))
+ tau_ngw(i) = tau_ngw(i) * max(0.1, min(5.0, tx1))
+ enddo
+ endif
+ do i=1,im
+ tau_ngw(i) = tau_ngw(i) * cdmbgwd(3)
+ enddo
+ endif
!
- call fv3_ugwp_solv2_v0(im, levs, dtp,
- & tgrs, ugrs, vgrs, qgrs, prsl, prsi, phil, xlatd,
- & sinlat, coslat, gw_dudt, gw_dvdt, gw_dTdt, gw_kdis,
- & tau_ngw, me, master, kdt )
-
- if (me == master .and. kdt < 2) then
- print *
- write(6,*)'FV3GFS finished fv3_ugwp_v0 in ugwp_driver_v0 '
- write(6,*) ' non-stationary GWs with GMAO/MERRA GW-forcing '
- print *
+ call fv3_ugwp_solv2_v0(im, levs, dtp,
+ & tgrs, ugrs, vgrs, qgrs, prsl, prsi,
+ & phil, xlatd, sinlat, coslat,
+ & gw_dudt, gw_dvdt, gw_dTdt, gw_kdis,
+ & tau_ngw, me, master, kdt)
+
+ if (me == master .and. kdt < 2) then
+ print *
+ write(6,*)'FV3GFS finished fv3_ugwp_v0 in ugwp_driver_v0 '
+ write(6,*) ' non-stationary GWs with GMAO/MERRA GW-forcing '
+ print *
+ endif
+ do k=1,levs
+ do i=1,im
+ gw_dtdt(i,k) = pngw*gw_dtdt(i,k) + pogw*Pdtdt(i,k)
+ gw_dudt(i,k) = pngw*gw_dudt(i,k) + pogw*Pdudt(i,k)
+ gw_dvdt(i,k) = pngw*gw_dvdt(i,k) + pogw*Pdvdt(i,k)
+ gw_kdis(i,k) = pngw*gw_kdis(i,k) + pogw*Pkdis(i,k)
+ enddo
+ enddo
+ else
+ do k=1,levs
+ do i=1,im
+ gw_dtdt(i,k) = Pdtdt(i,k)
+ gw_dudt(i,k) = Pdudt(i,k)
+ gw_dvdt(i,k) = Pdvdt(i,k)
+ gw_kdis(i,k) = Pkdis(i,k)
+ enddo
+ enddo
endif
- do k=1,levs
- do i=1,im
- gw_dtdt(i,k) = pngw*gw_dtdt(i,k) + pogw*Pdtdt(i,k)
- gw_dudt(i,k) = pngw*gw_dudt(i,k) + pogw*Pdudt(i,k)
- gw_dvdt(i,k) = pngw*gw_dvdt(i,k) + pogw*Pdvdt(i,k)
- gw_kdis(i,k) = pngw*gw_kdis(i,k) + pogw*Pkdis(i,k)
- enddo
- enddo
if (pogw == 0.0) then
! zmtb = 0.; zogw =0.
tau_mtb = 0.0 ; tau_ogw = 0.0 ; tau_tofd = 0.0
@@ -176,9 +229,13 @@ subroutine cires_ugwp_driver_v0(me, master,
!
! 3) application of "eddy"-diffusion to "smooth" UGWP-related tendencies
!------------------------------------------------------------------------------
- ed_dudt(:,:) = 0.0 ; ed_dvdt(:,:) = 0.0 ; ed_dtdt(:,:) = 0.0
+ do k=1,levs
+ do i=1,im
+ ed_dudt(i,k) = 0.0 ; ed_dvdt(i,k) = 0.0 ; ed_dtdt(i,k) = 0.0
+ enddo
+ enddo
- call edmix_ugwp_v0(im, levs, dtp,
+ call edmix_ugwp_v0(im, levs, dtp,
& tgrs, ugrs, vgrs, qgrs, del,
& prsl, prsi, phil, prslk,
& gw_dudt, gw_dvdt, gw_dTdt, gw_kdis,
@@ -193,14 +250,15 @@ subroutine cires_ugwp_driver_v0(me, master,
enddo
enddo
- end subroutine cires_ugwp_driver_v0
+ end subroutine cires_ugwp_driver_v0
+#endif
!
!=====================================================================
!
!ugwp-v0 subroutines: GWDPS_V0 and fv3_ugwp_solv2_v0
!
!=====================================================================
- SUBROUTINE GWDPS_V0(IM, levs, imx, do_tofd,
+ SUBROUTINE GWDPS_V0(IM, km, imx, do_tofd,
& Pdvdt, Pdudt, Pdtdt, Pkdis, U1,V1,T1,Q1,KPBL,
& PRSI,DEL,PRSL,PRSLK,PHII, PHIL,DTP,KDT,
& sgh30, HPRIME,OC,OA4,CLX4,THETA,vSIGMA,vGAMMA,ELVMAXD,
@@ -236,20 +294,21 @@ SUBROUTINE GWDPS_V0(IM, levs, imx, do_tofd,
!----------------------------------------
implicit none
character(len=8) :: strsolver='PSS-1986' ! current operational solver or 'WAM-2017'
- integer, intent(in) :: im, levs, imx, kdt
+ integer, intent(in) :: im, km, imx, kdt
integer, intent(in) :: me, master
logical, intent(in) :: do_tofd
- real(kind=kind_phys), parameter :: sigfac =3, sigfacS = 0.5
+ real(kind=kind_phys), parameter :: sigfac = 3, sigfacS = 0.5
real(kind=kind_phys) :: ztopH,zlowH,ph_blk, dz_blk
integer, intent(in) :: KPBL(IM) ! Index for the PBL top layer!
real(kind=kind_phys), intent(in) :: dtp ! time step
real(kind=kind_phys), intent(in) :: cdmbgwd(2)
- real(kind=kind_phys), intent(in), dimension(im,levs) ::
+ real(kind=kind_phys), intent(in), dimension(im,km) ::
& u1, v1, t1, q1,
& del, prsl, prslk, phil
- real(kind=kind_phys), intent(in),dimension(im,levs+1):: prsi, phii
- real(kind=kind_phys), intent(in) ::xlatd(im),sinlat(im),coslat(im)
+ real(kind=kind_phys), intent(in),dimension(im,km+1):: prsi, phii
+ real(kind=kind_phys), intent(in) :: xlatd(im),sinlat(im),
+ & coslat(im)
real(kind=kind_phys), intent(in) :: sparea(im)
real(kind=kind_phys), intent(in) :: OC(IM), OA4(im,4), CLX4(im,4)
@@ -259,7 +318,7 @@ SUBROUTINE GWDPS_V0(IM, levs, imx, do_tofd,
real(kind=kind_phys) :: SIGMA(IM), GAMMA(IM)
!output -phys-tend
- real(kind=kind_phys),dimension(im,levs),intent(out) ::
+ real(kind=kind_phys),dimension(im,km),intent(out) ::
& Pdvdt, Pdudt, Pkdis, Pdtdt
! output - diag-coorde
&, dudt_mtb, dudt_ogw, dudt_tms
@@ -267,18 +326,39 @@ SUBROUTINE GWDPS_V0(IM, levs, imx, do_tofd,
real(kind=kind_phys),dimension(im) :: RDXZB, zmtb, zogw
&, tau_ogw, tau_mtb, tau_tofd
&, dusfc, dvsfc
+!
+!---------------------------------------------------------------------
+! # of permissible sub-grid orography hills for "any" resolution < 25
+! correction for "elliptical" hills based on shilmin-area =sgrid/25
+! 4.*gamma*b_ell*b_ell >= shilmin
+! give us limits on [b_ell & gamma *b_ell] > 5 km =sso_min
+! gamma_min = 1/4*shilmin/sso_min/sso_min
+!23.01.2019: cdmb = 4.*192/768_c192=1 x 0.5
+! 192: cdmbgwd = 0.5, 2.5
+! cleff = 2.5*0.5e-5 * sqrt(192./768.) => Lh_eff = 1004. km
+! 6*dx = 240 km 8*dx = 320. ~ 3-5 more effective
+!---------------------------------------------------------------------
+ real(kind=kind_phys) :: gammin = 0.00999999
+ real(kind=kind_phys), parameter :: nhilmax = 25.
+ real(kind=kind_phys), parameter :: sso_min = 3000.
+ logical, parameter :: do_adjoro = .true.
+!
+ real(kind=kind_phys) :: shilmin, sgrmax, sgrmin
+ real(kind=kind_phys) :: belpmin, dsmin, dsmax
+! real(kind=kind_phys) :: arhills(im) ! not used why do we need?
+ real(kind=kind_phys) :: xlingfs
!
! locals
! mean flow
- real(kind=kind_phys) :: RI_N(IM,levs), BNV2(IM,levs), RO(IM,levs)
- real(kind=kind_phys) :: VTK(IM,levs),VTJ(IM,levs),VELCO(IM,levs)
+ real(kind=kind_phys), dimension(im,km) :: RI_N, BNV2, RO
+ &, VTK, VTJ, VELCO
!mtb
- real(kind=kind_phys) :: OA(IM), CLX(IM) , elvmax(im)
- real(kind=kind_phys) :: wk(IM)
- real(kind=kind_phys), dimension(im) :: PE, EK, UP
+ real(kind=kind_phys), dimension(im) :: OA, CLX , elvmax, wk
+ &, PE, EK, UP
- real(kind=kind_phys) :: DB(IM,levs),ANG(IM,levs),UDS(IM, levs)
+ real(kind=kind_phys), dimension(im,km) :: DB, ANG, UDS
+
real(kind=kind_phys) :: ZLEN, DBTMP, R, PHIANG, DBIM, ZR
real(kind=kind_phys) :: ENG0, ENG1, COSANG2, SINANG2
real(kind=kind_phys) :: bgam, cgam, gam2, rnom, rdem
@@ -287,83 +367,61 @@ SUBROUTINE GWDPS_V0(IM, levs, imx, do_tofd,
! Some constants now in "use ugwp_oro_init" + "use ugwp_common"
!
!==================
- real(kind=kind_phys) :: unew, vnew, zpbl, sigflt
- real(kind=kind_phys), dimension(levs) :: utofd1, vtofd1
- &, epstofd1, krf_tofd1
- &, up1, vp1, zpm
- real(kind=kind_phys) :: zsurf
- real(kind=kind_phys),dimension(im, levs) :: axtms, aytms
+ real(kind=kind_phys) :: unew, vnew, zpbl, sigflt, zsurf
+ real(kind=kind_phys), dimension(km) :: utofd1, vtofd1
+ &, epstofd1, krf_tofd1
+ &, up1, vp1, zpm
+ real(kind=kind_phys),dimension(im, km) :: axtms, aytms
!
! OGW
!
LOGICAL ICRILV(IM)
!
- real(kind=kind_phys) :: XN(IM), YN(IM), UBAR(IM),
- & VBAR(IM), ULOW(IM),
- & ROLL(IM), bnv2bar(im), SCOR(IM),
- & DTFAC(IM), XLINV(IM), DELKS(IM), DELKS1(IM)
+ real(kind=kind_phys), dimension(im) :: XN, YN, UBAR, VBAR, ULOW,
+ & ROLL, bnv2bar, SCOR, DTFAC, XLINV, DELKS, DELKS1
!
- real(kind=kind_phys) :: TAUP(IM,levs+1), TAUD(IM,levs)
+ real(kind=kind_phys) :: TAUP(IM,km+1), TAUD(IM,km)
real(kind=kind_phys) :: taub(im), taulin(im), heff, hsat, hdis
- integer :: kref(IM), idxzb(im), ipt(im), k_mtb,k_zlow
- integer :: kreflm(IM), iwklm(im), iwk(im), izlow(im)
- integer :: ktrial, klevm1
+ integer, dimension(im) :: kref, idxzb, ipt, kreflm,
+ & iwklm, iwk, izlow
!
!check what we need
!
- real(kind=kind_phys) :: bnv, fr, ri_gw ,
- & brvf, tem, tem1, tem2, temc, temv,
- & ti, rdz, dw2, shr2, bvf2,
- & rdelks, efact, coefm, gfobnv,
- & scork, rscor, hd, fro, sira,
- & dtaux, dtauy, pkp1log, pklog
-
- integer :: km, kmm1, kmm2, lcap, lcapp1
- &, npt, kbps, kbpsp1,kbpsm1
- &, kmps, idir, nwd, klcap, kp1, kmpbl, kmll
-!---------------------------------------------------------------------
-! # of permissible sub-grid orography hills for "any" resolution < 25
-! correction for "elliptical" hills based on shilmin-area =sgrid/25
-! 4.*gamma*b_ell*b_ell >= shilmin
-! give us limits on [b_ell & gamma *b_ell] > 5 km =sso_min
-! gamma_min = 1/4*shilmin/sso_min/sso_min
-!23.01.2019: cdmb = 4.*192/768_c192=1 x 0.5
-! 192: cdmbgwd = 0.5, 2.5
-! cleff = 2.5*0.5e-5 * sqrt(192./768.) => Lh_eff = 1004. km
-! 6*dx = 240 km 8*dx = 320. ~ 3-5 more effective
-!---------------------------------------------------------------------
- real(kind=kind_phys) :: gammin = 0.00999999
- real(kind=kind_phys) :: shilmin, sgrmax, sgrmin
- real(kind=kind_phys) :: belpmin, dsmin, dsmax
- real(kind=kind_phys), parameter :: nhilmax = 25.
- real(kind=kind_phys), parameter :: sso_min = 3000.
- real(kind=kind_phys) :: xlingfs
- real(kind=kind_phys) :: arhills(im)
- logical, parameter :: do_adjoro = .true.
-!
- integer :: i, j, k
- real(kind=kind_phys) :: grav2, rcpdt, windik, wdir
+ real(kind=kind_phys) :: bnv, fr, ri_gw
+ &, brvf, tem, tem1, tem2, temc, temv
+ &, ti, rdz, dw2, shr2, bvf2
+ &, rdelks, efact, coefm, gfobnv
+ &, scork, rscor, hd, fro, sira
+ &, dtaux, dtauy, pkp1log, pklog
+ &, grav2, rcpdt, windik, wdir
&, sigmin, dxres,sigres,hdxres
&, cdmb4, mtbridge
&, kxridge, inv_b2eff, zw1, zw2
&, belps, aelps, nhills, selps
-!
+
+ integer :: kmm1, kmm2, lcap, lcapp1
+ &, npt, kbps, kbpsp1,kbpsm1
+ &, kmps, idir, nwd, klcap, kp1, kmpbl, kmll
+ &, k_mtb, k_zlow, ktrial, klevm1, i, j, k
+!
rcpdt = 1.0 / (cpd*dtp)
grav2 = grav + grav
!
! mtb-blocking sigma_min and dxres => cires_initialize
!
sgrmax = maxval(sparea) ; sgrmin = minval(sparea)
- dsmax = sqrt(sgrmax) ; dsmin = sqrt(sgrmin)
+ dsmax = sqrt(sgrmax) ; dsmin = sqrt(sgrmin)
dxres = pi2*arad/float(IMX)
hdxres = 0.5*dxres
- shilmin = sgrmin/nhilmax
+! shilmin = sgrmin/nhilmax ! not used - Moorthi
- gammin = min(sso_min/dsmax, 1.)
+! gammin = min(sso_min/dsmax, 1.) ! Moorthi - with this results are not reproducible
+ gammin = min(sso_min/dxres, 1.) ! Moorthi
- sigmin = 2.*hpmin/dsmax !dxres
+! sigmin = 2.*hpmin/dsmax !dxres ! Moorthi - this will not reproduce
+ sigmin = 2.*hpmin/dxres !dxres
! if (kdt == 1) then
! print *, sgrmax, sgrmin , ' min-max sparea '
@@ -371,10 +429,10 @@ SUBROUTINE GWDPS_V0(IM, levs, imx, do_tofd,
! print *, 'dxres/dsmax ', dxres, dsmax
! print *, ' shilmin gammin ', shilmin, gammin
! endif
-
+
kxridge = float(IMX)/arad * cdmbgwd(2)
-
- if (me == master .and. kdt==1) then
+
+ if (me == master .and. kdt == 1) then
print *, ' gwdps_v0 kxridge ', kxridge
print *, ' gwdps_v0 scale2 ', cdmbgwd(2)
print *, ' gwdps_v0 IMX ', imx
@@ -383,7 +441,7 @@ SUBROUTINE GWDPS_V0(IM, levs, imx, do_tofd,
endif
do i=1,im
- idxzb(:) = 0
+ idxzb(i) = 0
zmtb(i) = 0.0
zogw(i) = 0.0
rdxzb(i) = 0.0
@@ -392,9 +450,13 @@ SUBROUTINE GWDPS_V0(IM, levs, imx, do_tofd,
dusfc(i) = 0.0
dvsfc(i) = 0.0
tau_tofd(i) = 0.0
+!
+ ipt(i) = 0
+ sigma(i) = max(vsigma(i), sigmin)
+ gamma(i) = max(vgamma(i), gammin)
enddo
-
- do k=1,levs
+
+ do k=1,km
do i=1,im
pdvdt(i,k) = 0.0
pdudt(i,k) = 0.0
@@ -408,56 +470,48 @@ SUBROUTINE GWDPS_V0(IM, levs, imx, do_tofd,
! ---- for lm and gwd calculation points
- ipt(:) = 0
npt = 0
- sigma = vsigma
- gamma = vgamma
do i = 1,im
- if ( (elvmaxd(i) >= hminmt)
- & .and. (gamma(i) >= gammin)
- & .and. (hprime(i) >= hpmin) ) then
+ if ( elvmaxd(i) >= hminmt .and. hprime(i) >= hpmin ) then
- npt = npt + 1
- ipt(npt) = i
- arhills(i) = 1.0
-!
- if (gamma(i) < gammin) gamma(i) = gammin
- sigres = max(sigmin, sigma(i))
- if (sigma(i) < sigmin) sigma(i)= sigmin
- dxres = sqrt(sparea(i))
- if (2.*hprime(i)/sigres > dxres) sigres=2.*hprime(i)/dxres
- aelps = min(2.*hprime(i)/sigres, 0.5*dxres)
- if (gamma(i) > 0.0 ) belps=min(aelps/gamma(i),.5*dxres)
+ npt = npt + 1
+ ipt(npt) = i
+! arhills(i) = 1.0
+!
+ sigres = max(sigmin, sigma(i))
+! if (sigma(i) < sigmin) sigma(i)= sigmin
+ dxres = sqrt(sparea(i))
+ if (2.*hprime(i)/sigres > dxres) sigres=2.*hprime(i)/dxres
+ aelps = min(2.*hprime(i)/sigres, 0.5*dxres)
+ if (gamma(i) > 0.0 ) belps = min(aelps/gamma(i),.5*dxres)
!
! small-scale "turbulent" oro-scales < sso_min
!
- if( aelps < sso_min .and. do_adjoro) then
+ if( aelps < sso_min .and. do_adjoro) then
! a, b > sso_min upscale ellipse a/b > 0.1 a>sso_min & h/b=>new_sigm
!
- aelps = sso_min
- if (belps < sso_min ) then
- gamma(i) = 1.0
- belps = aelps*gamma(i)
-
- else
- gamma(i) = min(aelps/belps, 1.0)
- endif
- sigma(i) = 2.*hprime(i)/aelps
- gamma(i) = min(aelps/belps, 1.0)
- endif
+ aelps = sso_min
+ if (belps < sso_min ) then
+ gamma(i) = 1.0
+ belps = aelps*gamma(i)
+ else
+ gamma(i) = min(aelps/belps, 1.0)
+ endif
+ sigma(i) = 2.*hprime(i)/aelps
+ gamma(i) = min(aelps/belps, 1.0)
+ endif
- selps = belps*belps*gamma(i)*4. ! ellipse area of the el-c hill
- nhills = sparea(i)/selps
- if (nhills > nhilmax) nhills = nhilmax
- arhills(i) = max(nhills, 1.0)
+ selps = belps*belps*gamma(i)*4. ! ellipse area of the el-c hill
+ nhills = min(nhilmax, sparea(i)/selps)
+! arhills(i) = max(nhills, 1.0)
!333 format( ' nhil: ', I6, 4(2x, F9.3), 2(2x, E9.3))
! if (kdt==1 )
! & write(6,333) nint(nhills)+1,xlatd(i), hprime(i),aelps*1.e-3,
! & belps*1.e-3, sigma(i),gamma(i)
- endif
+ endif
enddo
IF (npt == 0) then
@@ -473,7 +527,7 @@ SUBROUTINE GWDPS_V0(IM, levs, imx, do_tofd,
kreflm(i) = 0
enddo
- do k=1,levs
+ do k=1,km
do i=1,im
db(i,k) = 0.0
ang(i,k) = 0.0
@@ -481,17 +535,16 @@ SUBROUTINE GWDPS_V0(IM, levs, imx, do_tofd,
enddo
enddo
- km = levs
- KMM1 = levs- 1 ; KMM2 = levs - 2 ; KMLL = kmm1
- LCAP = levs ; LCAPP1 = LCAP + 1
-
+ KMM1 = km - 1 ; KMM2 = km - 2 ; KMLL = kmm1
+ LCAP = km ; LCAPP1 = LCAP + 1
+
DO I = 1, npt
j = ipt(i)
ELVMAX(J) = min (ELVMAXd(J)*0. + sigfac * hprime(j), hncrit)
+ izlow(i) = 1 ! surface-level
ENDDO
!
- izlow(:) =1 ! surface-level
- DO K = 1, levs-1
+ DO K = 1, kmm1
DO I = 1, npt
j = ipt(i)
ztopH = sigfac * hprime(j)
@@ -508,7 +561,7 @@ SUBROUTINE GWDPS_V0(IM, levs, imx, do_tofd,
ENDDO
ENDDO
!
- DO K = 1,levs
+ DO K = 1,km
DO I =1,npt
J = ipt(i)
VTJ(I,K) = T1(J,K) * (1.+FV*Q1(J,K))
@@ -520,7 +573,7 @@ SUBROUTINE GWDPS_V0(IM, levs, imx, do_tofd,
!
! check RI_N or RI_MF computation
!
- DO K = 1,levs-1
+ DO K = 1,kmm1
DO I =1,npt
J = ipt(i)
RDZ = grav / (phil(j,k+1) - phil(j,k))
@@ -541,153 +594,154 @@ SUBROUTINE GWDPS_V0(IM, levs, imx, do_tofd,
!
ENDDO
ENDDO
- K = 1
- DO I = 1, npt
- bnv2(i,k) = bnv2(i,k+1)
- ENDDO
+ K = 1
+ DO I = 1, npt
+ bnv2(i,k) = bnv2(i,k+1)
+ ENDDO
!
! level iwklm =>phil(j,k)/g < sigfac * hprime(j) < phil(j,k+1)/g
!
- DO I = 1, npt
- J = ipt(i)
- k_zlow = izlow(I)
- if (k_zlow == iwklm(i)) k_zlow = 1
- DELKS(I) = 1.0 / (PRSI(J,k_zlow) - PRSI(J,iwklm(i)))
-! DELKS1(I) = 1.0 /(PRSL(J,k_zlow) - PRSL(J,iwklm(i)))
- UBAR (I) = 0.0
- VBAR (I) = 0.0
- ROLL (I) = 0.0
- PE (I) = 0.0
- EK (I) = 0.0
- BNV2bar(I) = 0.0
- ENDDO
+ DO I = 1, npt
+ J = ipt(i)
+ k_zlow = izlow(I)
+ if (k_zlow == iwklm(i)) k_zlow = 1
+ DELKS(I) = 1.0 / (PRSI(J,k_zlow) - PRSI(J,iwklm(i)))
+! DELKS1(I) = 1.0 /(PRSL(J,k_zlow) - PRSL(J,iwklm(i)))
+ UBAR (I) = 0.0
+ VBAR (I) = 0.0
+ ROLL (I) = 0.0
+ PE (I) = 0.0
+ EK (I) = 0.0
+ BNV2bar(I) = 0.0
+ ENDDO
!
- DO I = 1, npt
- k_zlow = izlow(I)
- if (k_zlow == iwklm(i)) k_zlow = 1
- DO K = k_zlow, iwklm(I)-1 ! Kreflm(I)= iwklm(I)-1
- J = ipt(i) ! laye-aver Rho, U, V
- RDELKS = DEL(J,K) * DELKS(I)
- UBAR(I) = UBAR(I) + RDELKS * U1(J,K) ! trial Mean U below
- VBAR(I) = VBAR(I) + RDELKS * V1(J,K) ! trial Mean V below
- ROLL(I) = ROLL(I) + RDELKS * RO(I,K) ! trial Mean RO below
-!
- BNV2bar(I) = BNV2bar(I) + .5*(BNV2(I,K)+BNV2(I,K+1))* RDELKS
- ENDDO
+ DO I = 1, npt
+ k_zlow = izlow(I)
+ if (k_zlow == iwklm(i)) k_zlow = 1
+ DO K = k_zlow, iwklm(I)-1 ! Kreflm(I)= iwklm(I)-1
+ J = ipt(i) ! laye-aver Rho, U, V
+ RDELKS = DEL(J,K) * DELKS(I)
+ UBAR(I) = UBAR(I) + RDELKS * U1(J,K) ! trial Mean U below
+ VBAR(I) = VBAR(I) + RDELKS * V1(J,K) ! trial Mean V below
+ ROLL(I) = ROLL(I) + RDELKS * RO(I,K) ! trial Mean RO below
+!
+ BNV2bar(I) = BNV2bar(I) + .5*(BNV2(I,K)+BNV2(I,K+1))* RDELKS
ENDDO
+ ENDDO
!
- DO I = 1, npt
- J = ipt(i)
+ DO I = 1, npt
+ J = ipt(i)
!
! integrate from Ztoph = sigfac*hprime down to Zblk if exists
! find ph_blk, dz_blk like in LM-97 and IFS
!
- ph_blk =0.
- DO K = iwklm(I), 1, -1
- PHIANG = atan2(V1(J,K),U1(J,K))*RAD_TO_DEG
- ANG(I,K) = ( THETA(J) - PHIANG )
- if ( ANG(I,K) > 90. ) ANG(I,K) = ANG(I,K) - 180.
- if ( ANG(I,K) < -90. ) ANG(I,K) = ANG(I,K) + 180.
- ANG(I,K) = ANG(I,K) * DEG_TO_RAD
- UDS(I,K) =
- & MAX(SQRT(U1(J,K)*U1(J,K) + V1(J,K)*V1(J,K)), velmin)
-!
- IF (IDXZB(I) == 0 ) then
- dz_blk=( PHII(J,K+1) - PHII(J,K) ) *rgrav
- PE(I) = PE(I) + BNV2(I,K) *
- & ( ELVMAX(J) - phil(J,K)*rgrav ) * dz_blk
-
- UP(I) = max(UDS(I,K) * cos(ANG(I,K)), velmin)
- EK(I) = 0.5 * UP(I) * UP(I)
-
- ph_blk = ph_blk + dz_blk*sqrt(BNV2(I,K))/UP(I)
+ ph_blk =0.
+ DO K = iwklm(I), 1, -1
+ PHIANG = atan2(V1(J,K),U1(J,K))*RAD_TO_DEG
+ ANG(I,K) = ( THETA(J) - PHIANG )
+ if ( ANG(I,K) > 90. ) ANG(I,K) = ANG(I,K) - 180.
+ if ( ANG(I,K) < -90. ) ANG(I,K) = ANG(I,K) + 180.
+ ANG(I,K) = ANG(I,K) * DEG_TO_RAD
+ UDS(I,K) =
+ & MAX(SQRT(U1(J,K)*U1(J,K) + V1(J,K)*V1(J,K)), velmin)
+!
+ IF (IDXZB(I) == 0 ) then
+ dz_blk = ( PHII(J,K+1) - PHII(J,K) ) *rgrav
+ PE(I) = PE(I) + BNV2(I,K) *
+ & ( ELVMAX(J) - phil(J,K)*rgrav ) * dz_blk
+
+ UP(I) = max(UDS(I,K) * cos(ANG(I,K)), velmin)
+ EK(I) = 0.5 * UP(I) * UP(I)
+
+ ph_blk = ph_blk + dz_blk*sqrt(BNV2(I,K))/UP(I)
! --- Dividing Stream lime is found when PE =exceeds EK. oper-l GFS
-! IF ( PE(I) >= EK(I) ) THEN
- IF ( ph_blk >= fcrit_gfs ) THEN
- IDXZB(I) = K
- zmtb (J) = PHIL(J, K)*rgrav
- RDXZB(J) = real(k, kind=kind_phys)
- ENDIF
-
+! IF ( PE(I) >= EK(I) ) THEN
+ IF ( ph_blk >= fcrit_gfs ) THEN
+ IDXZB(I) = K
+ zmtb (J) = PHIL(J, K)*rgrav
+ RDXZB(J) = real(k, kind=kind_phys)
ENDIF
- ENDDO
+
+ ENDIF
+ ENDDO
!
! Alternative expression: ZMTB = max(Heff*(1. -Fcrit_gfs/Fr), 0)
! fcrit_gfs/fr
!
- goto 788
-
- BNV = SQRT( BNV2bar(I) )
- heff = 2.*min(HPRIME(J),hpmax)
- zw2 = UBAR(I)*UBAR(I)+VBAR(I)*VBAR(I)
- Ulow(i) = sqrt(max(zw2,dw2min))
- Fr = heff*bnv/Ulow(i)
- ZW1 = max(Heff*(1. -fcrit_gfs/fr), 0.0)
- zw2 = phil(j,2)*rgrav
- if (Fr > fcrit_gfs .and. zw1 > zw2 ) then
- do k=2, levs-1
+ goto 788
+
+ BNV = SQRT( BNV2bar(I) )
+ heff = 2.*min(HPRIME(J),hpmax)
+ zw2 = UBAR(I)*UBAR(I)+VBAR(I)*VBAR(I)
+ Ulow(i) = sqrt(max(zw2,dw2min))
+ Fr = heff*bnv/Ulow(i)
+ ZW1 = max(Heff*(1. -fcrit_gfs/fr), 0.0)
+ zw2 = phil(j,2)*rgrav
+ if (Fr > fcrit_gfs .and. zw1 > zw2 ) then
+ do k=2, kmm1
pkp1log = phil(j,k+1) * rgrav
pklog = phil(j,k) * rgrav
- if (zw1 <= pkp1log .and. zw1 >= pklog) exit
- enddo
+ if (zw1 <= pkp1log .and. zw1 >= pklog) exit
+ enddo
IDXZB(I) = K
zmtb (J) = PHIL(J, K)*rgrav
- else
- zmtb (J) = 0.
- IDXZB(I) = 0
- endif
+ else
+ zmtb (J) = 0.
+ IDXZB(I) = 0
+ endif
788 continue
- ENDDO
+ ENDDO
!
! --- The drag for mtn blocked flow
!
- cdmb4 = 0.25*cdmb
- DO I = 1, npt
- J = ipt(i)
+ cdmb4 = 0.25*cdmb
+ DO I = 1, npt
+ J = ipt(i)
!
- IF ( IDXZB(I) > 0 ) then
+ IF ( IDXZB(I) > 0 ) then
! (4.16)-IFS
- gam2 = gamma(j)*gamma(j)
- BGAM = 1.0 - 0.18*gamma(j) - 0.04*gam2
- CGAM = 0.48*gamma(j) + 0.30*gam2
- DO K = IDXZB(I)-1, 1, -1
+ gam2 = gamma(j)*gamma(j)
+ BGAM = 1.0 - 0.18*gamma(j) - 0.04*gam2
+ CGAM = 0.48*gamma(j) + 0.30*gam2
+ DO K = IDXZB(I)-1, 1, -1
- ZLEN = SQRT( ( PHIL(J,IDXZB(I)) - PHIL(J,K) ) /
- & ( PHIL(J,K ) + Grav * hprime(J) ) )
+ ZLEN = SQRT( ( PHIL(J,IDXZB(I)) - PHIL(J,K) ) /
+ & ( PHIL(J,K ) + Grav * hprime(J) ) )
- COSANG2 = cos(ANG(I,K))*cos(ANG(I,K))
- SINANG2 = 1.0 - COSANG2
+ tem = cos(ANG(I,K))
+ COSANG2 = tem * tem
+ SINANG2 = 1.0 - COSANG2
!
! cos =1 sin =0 => 1/R= gam ZR = 2.-gam
! cos =0 sin =1 => 1/R= 1/gam ZR = 2.- 1/gam
!
- rdem = COSANG2 + GAM2 * SINANG2
- rnom = COSANG2*GAM2 + SINANG2
+ rdem = COSANG2 + GAM2 * SINANG2
+ rnom = COSANG2*GAM2 + SINANG2
!
! metOffice Dec 2010
! correction of H. Wells & A. Zadra for the
! aspect ratio of the hill seen by MF
! (1/R , R-inverse below: 2-R)
- rdem = max(rdem, 1.e-6)
- R = sqrt(rnom/rdem)
- ZR = MAX( 2. - R, 0. )
+ rdem = max(rdem, 1.e-6)
+ R = sqrt(rnom/rdem)
+ ZR = MAX( 2. - R, 0. )
- sigres = max(sigmin, sigma(J))
- if (hprime(J)/sigres > dxres) sigres = hprime(J)/dxres
- mtbridge = ZR * sigres*ZLEN / hprime(J)
+ sigres = max(sigmin, sigma(J))
+ if (hprime(J)/sigres > dxres) sigres = hprime(J)/dxres
+ mtbridge = ZR * sigres*ZLEN / hprime(J)
! (4.15)-IFS
-! DBTMP = CDmb4 * mtbridge *
-! & MAX(cos(ANG(I,K)), gamma(J)*sin(ANG(I,K)))
+! DBTMP = CDmb4 * mtbridge *
+! & MAX(cos(ANG(I,K)), gamma(J)*sin(ANG(I,K)))
! (4.16)-IFS
- DBTMP = CDmb4*mtbridge*(bgam* COSANG2 +cgam* SINANG2)
- DB(I,K)= DBTMP * UDS(I,K)
- ENDDO
+ DBTMP = CDmb4*mtbridge*(bgam* COSANG2 +cgam* SINANG2)
+ DB(I,K)= DBTMP * UDS(I,K)
+ ENDDO
!
- endif
- ENDDO
+ endif
+ ENDDO
!
!.............................
!.............................
@@ -724,15 +778,15 @@ SUBROUTINE GWDPS_V0(IM, levs, imx, do_tofd,
! LEVEL ~0.4-0.5 KM from surface or/and PBL-top
! in UGWP-V1: options to modify as Htop ~ (2-3)*Hprime > Zmtb
! in UGWP-V0 we ensured that : Zogw > Zmtb
-!
+!
KBPS = 1
- KMPS = levs
- K_mtb = 1
+ KMPS = km
+ K_mtb = 1
DO I=1,npt
J = ipt(i)
K_mtb = max(1, idxzb(i))
-
+
kref(I) = MAX(IWK(I), KPBL(J)+1 ) ! reference level PBL or smt-else ????
kref(I) = MAX(kref(i), iwklm(i) ) ! iwklm => sigfac*hprime
@@ -746,11 +800,11 @@ SUBROUTINE GWDPS_V0(IM, levs, imx, do_tofd,
ROLL (I) = 0.0
BNV2bar(I)= 0.0
ENDDO
-!
+!
KBPSP1 = KBPS + 1
KBPSM1 = KBPS - 1
- K_mtb = 1
-!
+ K_mtb = 1
+!
DO I = 1,npt
K_mtb = max(1, idxzb(i))
DO K = k_mtb,KBPS !KBPS = MAX(kref) ;KMPS= MIN(kref)
@@ -765,7 +819,7 @@ SUBROUTINE GWDPS_V0(IM, levs, imx, do_tofd,
ENDDO
ENDDO
!
-! orographic asymmetry parameter (OA), and (CLX)
+! orographic asymmetry parameter (OA), and (CLX)
DO I = 1,npt
J = ipt(i)
wdir = atan2(UBAR(I),VBAR(I)) + pi
@@ -777,13 +831,13 @@ SUBROUTINE GWDPS_V0(IM, levs, imx, do_tofd,
!
DO I = 1,npt
DTFAC(I) = 1.0
- ICRILV(I) = .FALSE. ! INITIALIZE CRITICAL LEVEL CONTROL VECTOR
+ ICRILV(I) = .FALSE. ! INITIALIZE CRITICAL LEVEL CONTROL VECTOR
ULOW(I) = MAX(SQRT(UBAR(I)*UBAR(I)+VBAR(I)*VBAR(I)),velmin)
XN(I) = UBAR(I) / ULOW(I)
- YN(I) = VBAR(I) / ULOW(I)
+ YN(I) = VBAR(I) / ULOW(I)
ENDDO
!
- DO K = 1, levs-1
+ DO K = 1, kmm1
DO I = 1,npt
J = ipt(i)
VELCO(I,K) = 0.5 * ((U1(J,K)+U1(J,K+1))*XN(I)
@@ -935,7 +989,7 @@ SUBROUTINE GWDPS_V0(IM, levs, imx, do_tofd,
!
! zero momentum deposition at the top model layer
!
- taup(1:npt,levs+1) = taup(1:npt,levs)
+ taup(1:npt,km+1) = taup(1:npt,km)
!
! Calculate wave acc-n: - (grav)*d(tau)/d(p) = taud
!
@@ -948,7 +1002,7 @@ SUBROUTINE GWDPS_V0(IM, levs, imx, do_tofd,
!------scale MOMENTUM DEPOSITION AT TOP TO 1/2 VALUE
! it is zero now
! DO I = 1,npt
-! TAUD(I, levs) = TAUD(I,levs) * FACTOP
+! TAUD(I, km) = TAUD(I,km) * FACTOP
! ENDDO
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -985,7 +1039,7 @@ SUBROUTINE GWDPS_V0(IM, levs, imx, do_tofd,
! XLINV(I) = max(kxridge, inv_b2eff) ! 0.5*sigma(j)/heff = 1./Lridge
dtfac(:) = 1.0
- call oro_wam_2017(im, levs, npt, ipt, kref, kdt, me, master,
+ call oro_wam_2017(im, km, npt, ipt, kref, kdt, me, master,
& dtp, dxres, taub, u1, v1, t1, xn, yn, bnv2, ro, prsi,prsL,
& del, sigma, hprime, gamma, theta,
& sinlat, xlatd, taup, taud, pkdis)
@@ -1009,16 +1063,16 @@ SUBROUTINE GWDPS_V0(IM, levs, imx, do_tofd,
sigflt = min(sgh30(j), 0.3*hprime(j)) ! cannot exceed 30% of LS-SSO
zsurf = phii(j,1)*rgrav
- do k=1,levs
+ do k=1,km
zpm(k) = phiL(j,k)*rgrav
up1(k) = u1(j,k)
vp1(k) = v1(j,k)
enddo
- call ugwp_tofd1d(levs, sigflt, elvmaxd(j), zsurf, zpbl,
+ call ugwp_tofd1d(km, sigflt, elvmaxd(j), zsurf, zpbl,
& up1, vp1, zpm, utofd1, vtofd1, epstofd1, krf_tofd1)
- do k=1,levs
+ do k=1,km
axtms(j,k) = utofd1(k)
aytms(j,k) = vtofd1(k)
!
@@ -1028,7 +1082,7 @@ SUBROUTINE GWDPS_V0(IM, levs, imx, do_tofd,
pdudt(J,k) = pdudt(J,k) + axtms(j,k)
enddo
!2018-diag
- tau_tofd(J) = sum( utofd1(1:levs)* del(j,1:levs))
+ tau_tofd(J) = sum( utofd1(1:km)* del(j,1:km))
enddo
ENDIF ! do_tofd
@@ -1098,11 +1152,11 @@ SUBROUTINE GWDPS_V0(IM, levs, imx, do_tofd,
tau_ogw(j) = -rgrav * tau_ogw(j)
tau_tofd(J) = -rgrav * tau_tofd(j)
ENDDO
-
+
RETURN
-!============ debug ------------------------------------------------
+!============ debug ------------------------------------------------
if (kdt <= 2 .and. me == 0) then
print *, 'vgw-oro done gwdps_v0 in ugwp-v0 step-proc ', kdt, me
!
@@ -1128,7 +1182,7 @@ SUBROUTINE GWDPS_V0(IM, levs, imx, do_tofd,
print *, maxval(prsL), minval(prsL), ' prsL '
print *, maxval(RO), minval(RO), ' RO-dens '
print *, maxval(bnv2(1:npt,:)), minval(bnv2(1:npt,:)),' BNV2 '
- print *, maxval(kpbl), minval(kpbl), ' kpbl '
+ print *, maxval(kpbl), minval(kpbl), ' kpbl '
print *, maxval(sgh30), maxval(hprime), maxval(elvmax),'oro-d'
print *
do i =1, npt
@@ -1185,9 +1239,9 @@ end subroutine gwdps_v0
! (c) guidance from high-res runs for GW sources and res-aware tune-ups
!23456
!
-! call gwdrag_wam(1, im, ix, levs, ksrc, dtp,
+! call gwdrag_wam(1, im, ix, km, ksrc, dtp,
! & xlat, gw_dudt, gw_dvdt, taux, tauy)
-! call fv3_ugwp_wms17(kid1, im, ix, levs, ksrc_ifs, dtp,
+! call fv3_ugwp_wms17(kid1, im, ix, km, ksrc_ifs, dtp,
! & adt,adu,adv,prsl,prsi,phil,xlat, gw_dudt, gw_dvdt, gw_dtdt, gw_ked,
! & taux,tauy,grav, amol_i, me, lstep_first )
!
@@ -1196,9 +1250,10 @@ end subroutine gwdps_v0
subroutine fv3_ugwp_solv2_v0(klon, klev, dtime,
- & tm1 , um1, vm1, qm1,
- & prsl, prsi, philg, xlatd, sinlat, coslat,
- & pdudt, pdvdt, pdtdt, dked, tau_ngw, mpi_id, master, kdt)
+ & tm1 , um1, vm1, qm1,
+ & prsl, prsi, philg, xlatd, sinlat, coslat,
+ & pdudt, pdvdt, pdtdt, dked, tau_ngw,
+ & mpi_id, master, kdt)
!
@@ -1218,7 +1273,7 @@ subroutine fv3_ugwp_solv2_v0(klon, klev, dtime,
use ugwp_wmsdis_init, only : hpscale, rhp2, bv2min, gssec
&, v_kxw, v_kxw2, tamp_mpa, zfluxglob
&, maxdudt, gw_eff, dked_min
- &, nslope, ilaunch, zms
+ &, nslope, ilaunch, zmsi
&, zci, zdci, zci4, zci3, zci2
&, zaz_fct, zcosang, zsinang
&, nwav, nazd, zcimin, zcimax
@@ -1226,33 +1281,34 @@ subroutine fv3_ugwp_solv2_v0(klon, klev, dtime,
implicit none
!23456
- integer, intent(in) :: klev ! vertical level
- integer, intent(in) :: klon ! horiz tiles
-
- real ,intent(in) :: dtime ! model time step
- real ,intent(in) :: vm1(klon,klev) ! meridional wind
- real ,intent(in) :: um1(klon,klev) ! zonal wind
- real ,intent(in) :: qm1(klon,klev) ! spec. humidity
- real ,intent(in) :: tm1(klon,klev) ! kin temperature
-
- real ,intent(in) :: prsl(klon,klev) ! mid-layer pressure
- real ,intent(in) :: philg(klon,klev) ! m2/s2-phil => meters !!!!! phil =philg/grav
- real ,intent(in) :: prsi(klon,klev+1) ! prsi interface pressure
- real ,intent(in) :: xlatd(klon) ! lat was in radians, now with xlat_d in degrees
- real ,intent(in) :: sinlat(klon)
- real ,intent(in) :: coslat(klon)
- real ,intent(in) :: tau_ngw(klon)
-
- integer, intent(in):: mpi_id, master, kdt
+ integer, intent(in) :: klev ! vertical level
+ integer, intent(in) :: klon ! horiz tiles
+
+ real, intent(in) :: dtime ! model time step
+ real, intent(in) :: vm1(klon,klev) ! meridional wind
+ real, intent(in) :: um1(klon,klev) ! zonal wind
+ real, intent(in) :: qm1(klon,klev) ! spec. humidity
+ real, intent(in) :: tm1(klon,klev) ! kin temperature
+
+ real, intent(in) :: prsl(klon,klev) ! mid-layer pressure
+ real, intent(in) :: philg(klon,klev) ! m2/s2-phil => meters !!!!! phil =philg/grav
+ real, intent(in) :: prsi(klon,klev+1)! prsi interface pressure
+ real, intent(in) :: xlatd(klon) ! lat was in radians, now with xlat_d in degrees
+ real, intent(in) :: sinlat(klon)
+ real, intent(in) :: coslat(klon)
+ real, intent(in) :: tau_ngw(klon)
+
+ integer, intent(in) :: mpi_id, master, kdt
!
!
! out-gw effects
!
- real ,intent(out) :: pdudt(klon,klev) ! zonal momentum tendency
- real ,intent(out) :: pdvdt(klon,klev) ! meridional momentum tendency
- real ,intent(out) :: pdtdt(klon,klev) ! gw-heating (u*ax+v*ay)/cp
- real ,intent(out) :: dked(klon,klev) ! gw-eddy diffusion
+ real, intent(out) :: pdudt(klon,klev) ! zonal momentum tendency
+ real, intent(out) :: pdvdt(klon,klev) ! meridional momentum tendency
+ real, intent(out) :: pdtdt(klon,klev) ! gw-heating (u*ax+v*ay)/cp
+ real, intent(out) :: dked(klon,klev) ! gw-eddy diffusion
real, parameter :: minvel = 0.5 !
+ real, parameter :: epsln = 1.0d-12 !
!vay-2018
@@ -1278,12 +1334,12 @@ subroutine fv3_ugwp_solv2_v0(klon, klev, dtime,
!23456
real :: zul(klon,nazd) ! velocity in azimuthal direction at launch level
real :: zci_min(klon,nazd)
- real :: zcrt(klon,klev,nazd)
+! real :: zcrt(klon,klev,nazd) ! not used - do we need it? Moorthi
real :: zact(klon, nwav, nazd) ! if =1 then critical level encountered => c-u
- real :: zacc(klon, nwav, nazd)
+! real :: zacc(klon, nwav, nazd) ! not used!
!
real :: zpu(klon,klev, nazd) ! momentum flux
- real :: zdfl(klon,klev, nazd)
+! real :: zdfl(klon,klev, nazd)
real :: zfct(klon,klev)
real :: zfnorm(klon) ! normalisation factor
@@ -1298,7 +1354,7 @@ subroutine fv3_ugwp_solv2_v0(klon, klev, dtime,
real :: vm_zflx_mode, vc_zflx_mode
real :: kzw2, kzw3, kdsat, cdf2, cdf1, wdop2
- real :: zang, znorm, zang1, ztx
+! real :: zang, znorm, zang1, ztx
real :: zu, zcin, zcpeak, zcin4, zbvfl4
real :: zcin2, zbvfl2, zcin3, zbvfl3, zcinc
real :: zatmp, zfluxs, zdep, zfluxsq, zulm, zdft, ze1, ze2
@@ -1306,15 +1362,18 @@ subroutine fv3_ugwp_solv2_v0(klon, klev, dtime,
!
real :: zdelp,zrgpts
real :: zthstd,zrhostd,zbvfstd
- real :: tvc1, tvm1
+ real :: tvc1, tvm1, tem1, tem2, tem3
real :: zhook_handle
+ real :: delpi(klon,ilaunch:klev)
! real :: rcpd, grav2cpd
real, parameter :: rcpdl = cpd/grav ! 1/[g/cp] == cp/g
&, grav2cpd = grav/rcpdl ! g*(g/cp)= g^2/cp
+ &, cpdi = 1.0d0/cpd
- real :: fmode, expdis, fdis
+ real :: expdis, fdis
+! real :: fmode, expdis, fdis
real :: v_kzi, v_kzw, v_cdp, v_wdp, sc, tx1
integer :: j, k, inc, jk, jl, iazi
@@ -1355,8 +1414,8 @@ subroutine fv3_ugwp_solv2_v0(klon, klev, dtime,
do jk=1,klev
do jl=1,klon
zpu(jl,jk,iazi) = 0.0
- zcrt(jl,jk,iazi) = 0.0
- zdfl(jl,jk,iazi) = 0.0
+! zcrt(jl,jk,iazi) = 0.0
+! zdfl(jl,jk,iazi) = 0.0
enddo
enddo
enddo
@@ -1381,7 +1440,8 @@ subroutine fv3_ugwp_solv2_v0(klon, klev, dtime,
! zrhohm1(jl,jk) = prsi(jl,jk)*rdi/zthm1(jl,jk) ! rho = p/(RTv)
zrhohm1(jl,jk) = prsi(jl,jk)*rdi*zthm1 ! rho = p/(RTv)
zdelp = phil(jl,jk)-phil(jl,jk-1) !>0 ...... dz-meters
- v_zmet(jl,jk) = 2.*zdelp
+ v_zmet(jl,jk) = zdelp + zdelp
+ delpi(jl,jk) = grav / (prsi(jl,jk-1) - prsi(jl,jk))
vueff(jl,jk) =
& 2.e-5*exp( (phil(jl,jk)+phil(jl,jk-1))*rhp2)+dked_min
!
@@ -1406,9 +1466,9 @@ subroutine fv3_ugwp_solv2_v0(klon, klev, dtime,
enddo
endif
do jl=1,klon
- tx1 = OMEGA2 * SINLAT(JL) / V_KXW
- C2F2(JL) = tx1 * tx1
- zbvfl(jl) = zbvfhm1(jl,ilaunch)
+ tx1 = OMEGA2 * SINLAT(JL) / V_KXW
+ C2F2(JL) = tx1 * tx1
+ zbvfl(jl) = zbvfhm1(jl,ilaunch)
enddo
!
! define intrinsic velocity (relative to launch level velocity) u(z)-u(zo), and coefficinets
@@ -1461,9 +1521,9 @@ subroutine fv3_ugwp_solv2_v0(klon, klev, dtime,
zcin = zci(inc)
zcin4 = zci4(inc)
do jl=1,klon
- zbvfl4 = zbvfl(jl)*zbvfl(jl)
- zbvfl4 = zbvfl4 * zbvfl4
- zcpeak = zbvfl(jl)/zms
+ zbvfl4 = zbvfl(jl) * zbvfl(jl)
+ zbvfl4 = zbvfl4 * zbvfl4
+ zcpeak = zbvfl(jl) * zmsi
zflux(jl,inc,1) = zfct(jl,ilaunch)*
& zbvfl4*zcin*zcpeak/(zbvfl4*zcpeak+zcin4*zcin)
enddo
@@ -1536,7 +1596,8 @@ subroutine fv3_ugwp_solv2_v0(klon, klev, dtime,
! copy zflux into all other azimuths
! --------------------------------
- zact(:,:,:) = 1.0 ; zacc(:,:,:) = 1.0
+! zact(:,:,:) = 1.0 ; zacc(:,:,:) = 1.0
+ zact(:,:,:) = 1.0
do iazi=2, nazd
do inc=1,nwav
do jl=1,klon
@@ -1549,6 +1610,8 @@ subroutine fv3_ugwp_solv2_v0(klon, klev, dtime,
! azimuth do-loop
! --------------------
do iazi=1, nazd
+
+! write(0,*)' iazi=',iazi,' ilaunch=',ilaunch
! vertical do-loop
! ----------------
do jk=ilaunch, klev-1
@@ -1560,44 +1623,52 @@ subroutine fv3_ugwp_solv2_v0(klon, klev, dtime,
! set zact to zero if critical level encountered
! ----------------------------------------------
do inc=1, nwav
- zcin = zci(inc)
+! zcin = zci(inc)
do jl=1,klon
- zatmp = minvel + sign(minvel,zcin-zci_min(jl,iazi))
- zacc(jl,inc,iazi) = zact(jl,inc,iazi)-zatmp
- zact(jl,inc,iazi) = zatmp
+! zatmp = minvel + sign(minvel,zcin-zci_min(jl,iazi))
+! zacc(jl,inc,iazi) = zact(jl,inc,iazi)-zatmp
+! zact(jl,inc,iazi) = zatmp
+ zact(jl,inc,iazi) = minvel
+ & + sign(minvel,zci(inc)-zci_min(jl,iazi))
enddo
enddo
!
+! zdfl not used! - do we need it? Moorthi
! integrate to get critical-level contribution to mom deposition
! ---------------------------------------------------------------
- do inc=1, nwav
- zcinc = zdci(inc)
- do jl=1,klon
- zdfl(jl,jk,iazi) = zdfl(jl,jk,iazi) +
- & zacc(jl,inc,iazi)*zflux(jl,inc,iazi)*zcinc
- enddo
- enddo
+! do inc=1, nwav
+! zcinc = zdci(inc)
+! do jl=1,klon
+! zdfl(jl,jk,iazi) = zdfl(jl,jk,iazi) +
+! & zacc(jl,inc,iazi)*zflux(jl,inc,iazi)*zcinc
+! enddo
+! enddo
! --------------------------------------------
-! get weighted average of phase speed in layer
+! get weighted average of phase speed in layer zcrt is not used - do we need it? Moorthi
! --------------------------------------------
- do jl=1,klon
- if(zdfl(jl,jk,iazi) > 0.0 ) then
- zatmp = zcrt(jl,jk,iazi)
- do inc=1, nwav
- zatmp = zatmp + zci(inc) *
- & zacc(jl,inc,iazi)*zflux(jl,inc,iazi)*zdci(inc)
- enddo
-!
- zcrt(jl,jk,iazi)=zatmp/zdfl(jl,jk,iazi)
- else
- zcrt(jl,jk,iazi)=zcrt(jl,jk-1,iazi)
- endif
- enddo
+! do jl=1,klon
+! write(0,*)' jk=',jk,' jl=',jl,' iazi=',iazi, zdfl(jl,jk,iazi)
+! if(zdfl(jl,jk,iazi) > epsln ) then
+! zatmp = zcrt(jl,jk,iazi)
+! do inc=1, nwav
+! zatmp = zatmp + zci(inc) *
+! & zacc(jl,inc,iazi)*zflux(jl,inc,iazi)*zdci(inc)
+! enddo
+!
+! zcrt(jl,jk,iazi) = zatmp / zdfl(jl,jk,iazi)
+! else
+! zcrt(jl,jk,iazi) = zcrt(jl,jk-1,iazi)
+! endif
+! enddo
!
do inc=1, nwav
zcin = zci(inc)
- zcinc = 1.0 / zcin
+ if (abs(zcin) > epsln) then
+ zcinc = 1.0 / zcin
+ else
+ zcinc = 1.0
+ endif
do jl=1,klon
!=======================================================================
! saturated limit wfit = kzw*kzw*kt; wfdt = wfit/(kxw*cx)*betat
@@ -1632,18 +1703,18 @@ subroutine fv3_ugwp_solv2_v0(klon, klev, dtime,
v_cdp = 0. ! no effects of reflected waves
endif
- fmode = zflux(jl,inc,iazi)
- fdis = fmode*expdis
+! fmode = zflux(jl,inc,iazi)
+! fdis = fmode*expdis
+ fdis = expdis * zflux(jl,inc,iazi)
!
! saturated flux + wave dissipation - Keddy_gwsat in UGWP-V1
! linsatdis = 1.0 , here: u'^2 ~ linsatdis* [v_cdp*v_cdp]
!
- zfluxs= zfct(jl,jk)*v_cdp*v_cdp*zcinc
+ zfluxs = zfct(jl,jk)*v_cdp*v_cdp*zcinc
!
! zfluxs= zfct(jl,jk)*(zcin-zui(jl,jk,iazi))**2/zcin
! flux_tot - sat.flux
!
-
zdep = zact(jl,inc,iazi)* (fdis-zfluxs)
if(zdep > 0.0 ) then
! subs on sat-limit
@@ -1662,7 +1733,7 @@ subroutine fv3_ugwp_solv2_v0(klon, klev, dtime,
zdfdz_v(:,jk,iazi) = 0.0
do inc=1, nwav
- zcinc=zdci(inc) ! dc-integration
+ zcinc = zdci(inc) ! dc-integration
do jl=1,klon
vc_zflx_mode = zact(jl,inc,iazi)*zflux(jl,inc,iazi)
zpu(jl,jk,iazi) = zpu(jl,jk,iazi) + vc_zflx_mode*zcinc
@@ -1673,8 +1744,9 @@ subroutine fv3_ugwp_solv2_v0(klon, klev, dtime,
! later sum over selected azimuths as "non-negative" scalars)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (jk > ilaunch)then
- zdelp = grav/(prsi(jl,jk-1)-prsi(jl,jk))*
- & abs(zcin-zui(jl,jk,iazi)) *zcinc
+! zdelp = grav/(prsi(jl,jk-1)-prsi(jl,jk))*
+! & abs(zcin-zui(jl,jk,iazi)) *zcinc
+ zdelp = delpi(jl,jk) * abs(zcin-zui(jl,jk,iazi)) *zcinc
vm_zflx_mode = zact(jl,inc,iazi)* zflux_z(jl,inc,jk-1)
if (vc_zflx_mode > vm_zflx_mode)
@@ -1690,7 +1762,7 @@ subroutine fv3_ugwp_solv2_v0(klon, klev, dtime,
! --------------
enddo ! end jk do-loop vertical loop
! ---------------
- enddo ! end nazd do-loop
+ enddo ! end nazd do-loop
! ----------------------------------------------------------------------------
! sum contribution for total zonal and meridional flux +
! energy dissipation
@@ -1703,15 +1775,15 @@ subroutine fv3_ugwp_solv2_v0(klon, klev, dtime,
enddo
enddo
+ tem3 = zaz_fct*cpdi
do iazi=1,nazd
+ tem1 = zaz_fct*zcosang(iazi)
+ tem2 = zaz_fct*zsinang(iazi)
do jk=ilaunch, klev-1
do jl=1,klon
- taux(jl,jk) = taux(jl,jk)
- & + zpu(jl,jk,iazi)*zaz_fct*zcosang(iazi) ! zaz_fct - "azimuth"-norm-n
- tauy(jl,jk) = tauy(jl,jk)
- & + zpu(jl,jk,iazi)*zaz_fct*zsinang(iazi)
- pdtdt(jl,jk) = pdtdt(jl,jk)
- & + zdfdz_v(jl,jk,iazi)*zaz_fct/cpd ! eps_dis =sum( +d(flux_e)/dz) > 0.
+ taux(jl,jk) = taux(jl,jk) + tem1 * zpu(jl,jk,iazi) ! zaz_fct - "azimuth"-norm-n
+ tauy(jl,jk) = tauy(jl,jk) + tem2 * zpu(jl,jk,iazi)
+ pdtdt(jl,jk) = pdtdt(jl,jk) + tem3 * zdfdz_v(jl,jk,iazi) ! eps_dis =sum( +d(flux_e)/dz) > 0.
enddo
enddo
@@ -1723,7 +1795,8 @@ subroutine fv3_ugwp_solv2_v0(klon, klev, dtime,
do jk=ilaunch,klev
do jl=1, klon
- zdelp = grav / (prsi(jl,jk-1)-prsi(jl,jk))
+! zdelp = grav / (prsi(jl,jk-1)-prsi(jl,jk))
+ zdelp = delpi(jl,jk)
ze1 = (taux(jl,jk)-taux(jl,jk-1))*zdelp
ze2 = (tauy(jl,jk)-tauy(jl,jk-1))*zdelp
if (abs(ze1) >= maxdudt ) then
@@ -1737,7 +1810,7 @@ subroutine fv3_ugwp_solv2_v0(klon, klev, dtime,
!
! Cx =0 based Cx=/= 0. above
!
- pdtdt(jl,jk) = (ze1*um1(jl,jk) + ze2*vm1(jl,jk))/cpd
+ pdtdt(jl,jk) = (ze1*um1(jl,jk) + ze2*vm1(jl,jk)) * cpdi
!
dked(jl,jk) = max(dked_min, pdtdt(jl,jk)/zbn2(jl,jk))
! if (dked(jl,jk) < 0) dked(jl,jk) = dked_min
@@ -1776,7 +1849,7 @@ end subroutine fv3_ugwp_solv2_v0
! after tests of OGW (new revision) and NGW with MERRA-2 forcing.
!
!-------------------------------------------------------------------------------
- subroutine edmix_ugwp_v0(im, levs, dtp,
+ subroutine edmix_ugwp_v0(im, levs, dtp,
& t1, u1, v1, q1, del,
& prsl, prsi, phil, prslk,
& pdudt, pdvdt, pdTdt, pkdis,
@@ -1848,7 +1921,7 @@ subroutine edmix_ugwp_v0(im, levs, dtp,
!
real(kind=kind_phys), parameter :: iPr_pt = 0.5, dw2min = 1.e-4
real(kind=kind_phys), parameter :: lturb = 30., sc2 = lturb*lturb
- real(kind=kind_phys),parameter :: ulturb=150.,sc2u=ulturb* ulturb
+ real(kind=kind_phys), parameter :: ulturb=150.,sc2u=ulturb* ulturb
real(kind=kind_phys), parameter :: ric =0.25
real(kind=kind_phys), parameter :: rimin = -10., prmin = 0.25
real(kind=kind_phys), parameter :: prmax = 4.0
@@ -1920,7 +1993,7 @@ subroutine edmix_ugwp_v0(im, levs, dtp,
Km(1:levs) = ksum(1:levs) * rho(1:levs)* rho(1:levs)
do j=1, nstab
- call diff_1d_wtend(levs, dtstab, Fw, Fw1, Km,
+ call diff_1d_wtend(levs, dtstab, Fw, Fw1, levs,
& del(i,:), Sw, Sw1)
Fw = Sw
Fw1 = Sw1
@@ -1950,13 +2023,15 @@ subroutine diff_1d_wtend(levs, dt, F, F1, Km, rdp, rdpm, S, S1)
real(kind=kind_phys) :: S(levs), S1(levs), F(levs), F1(levs)
real(kind=kind_phys) :: Km(levs), rdp(levs), rdpm(levs-1)
integer :: i, k
- real(kind=kind_phys) :: Km1, Kp1, ad, cd, bd
+ real(kind=kind_phys) :: Kp1, ad, cd, bd
+! real(kind=kind_phys) :: km1, Kp1, ad, cd, bd
! S(:) = 0.0 ; S1(:) = 0.0
!
! explicit diffusion solver
!
k = 1
- km1 = 0. ; ad =0.
+! km1 = 0. ; ad =0.
+ ad =0.
kp1 = .5*(Km(k)+Km(k+1))
cd = rdp(1)*rdpm(1)*kp1*dt
bd = 1. - cd - ad
@@ -1981,16 +2056,18 @@ subroutine diff_1d_ptend(levs, dt, F, Km, rdp, rdpm, S)
implicit none
integer :: levs
real(kind=kind_phys) :: dt
- real(kind=kind_phys) :: S(levs), S1(levs), F(levs), F1(levs)
+ real(kind=kind_phys) :: S(levs), S1(levs), F(levs), F1(levs)
real(kind=kind_phys) :: Km(levs), rdp(levs), rdpm(levs-1)
integer :: i, k
- real(kind=kind_phys) :: Km1, Kp1, ad, cd, bd
+ real(kind=kind_phys) :: Kp1, ad, cd, bd
+! real(kind=kind_phys) :: km1, Kp1, ad, cd, bd
!
! explicit "eddy" smoother for tendencies
!
k = 1
- km1 = 0. ; ad =0.
+! km1 = 0. ; ad =0.
+ ad =0.
kp1 = .5*(Km(k)+Km(k+1))
cd = rdp(1)*rdpm(1)*kp1*dt
bd = 1. -(cd +ad)
@@ -2003,6 +2080,6 @@ subroutine diff_1d_ptend(levs, dt, F, Km, rdp, rdpm, S)
bd = 1.-(ad +cd)
S(k) = cd*F(k+1) + ad *F(k-1) + bd *F(k)
enddo
- k =levs
+ k = levs
S(k) = F(k)
end subroutine diff_1d_ptend