Skip to content

Commit ad9d897

Browse files
grantfirlgrantfirl
committed
Merge branch 'ufs/dev' into feature/tendency_cleanup
2 parents 6be9617 + eee1edf commit ad9d897

3 files changed

Lines changed: 102 additions & 66 deletions

File tree

‎physics/CONV/SAMF/samfdeepcnv.f‎

Lines changed: 40 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -218,10 +218,9 @@ subroutine samfdeepcnv_run (im,km,nn,first_time_step,restart, &
218218
!
219219
! parameters for prognostic sigma closure
220220
real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km),
221-
& omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im),qadv(im,km),
222-
& tentr(im,km)
223-
real(kind=kind_phys) gravinv,invdelt,sigmind,sigminm,sigmins
224-
parameter(sigmind=0.01,sigmins=0.03,sigminm=0.01)
221+
& omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im),qadv(im,km)
222+
real(kind=kind_phys) gravinv,invdelt,sigmind,sigminm,sigmins,
223+
& wc_min, wc_eff
225224
logical flag_shallow, flag_mid
226225
c physical parameters
227226
! parameter(grav=grav,asolfac=0.958)
@@ -363,6 +362,21 @@ subroutine samfdeepcnv_run (im,km,nn,first_time_step,restart, &
363362
!************************************************************************
364363
!
365364
!
365+
! - Initialize parameters related to prognostic closure
366+
if (progsigma) then
367+
if (progomega) then
368+
sigmind = 0.03
369+
sigmins = 0.03
370+
sigminm = 0.03
371+
wc_min = 0.2
372+
else
373+
sigmind = 0.01
374+
sigmins = 0.03
375+
sigminm = 0.03
376+
wc_min = 0.2
377+
endif
378+
endif
379+
366380
km1 = km - 1
367381
!> - Initialize column-integrated and other single-value-per-column variable arrays.
368382
c
@@ -1153,7 +1167,6 @@ subroutine samfdeepcnv_run (im,km,nn,first_time_step,restart, &
11531167
do i=1,im
11541168
if(cnvflg(i) .and.
11551169
& (k > kbcon(i) .and. k < kmax(i))) then
1156-
tentr(i,k)=xlamue(i,k)*fent1(i,k)
11571170
tem = cxlamet(i) * frh(i,k) * fent2(i,k)
11581171
xlamue(i,k) = xlamue(i,k)*fent1(i,k) + tem
11591172
tem1 = cxlamdt(i) * frh(i,k)
@@ -1804,11 +1817,11 @@ subroutine samfdeepcnv_run (im,km,nn,first_time_step,restart, &
18041817
if (progomega) then
18051818
call progomega_calc(first_time_step,restart,im,km,
18061819
& kbcon1,ktcon,omegain,delt,del,zi,cnvflg,omegaout,
1807-
& grav,buo,drag,wush,tentr,bb1,bb2)
1820+
& grav,buo,drag,wush,bb1,bb2)
18081821
do k = 1, km
18091822
do i = 1, im
18101823
if (cnvflg(i)) then
1811-
if(k > kbcon1(i) .and. k < ktcon(i)) then
1824+
if(k >= kbcon1(i) .and. k < ktcon(i)) then
18121825
omega_u(i,k)=omegaout(i,k)
18131826
omega_u(i,k)=MAX(omega_u(i,k),-80.)
18141827
! Convert to m/s for use in convective time-scale:
@@ -1894,7 +1907,7 @@ subroutine samfdeepcnv_run (im,km,nn,first_time_step,restart, &
18941907
do k = 2, km1
18951908
do i = 1, im
18961909
if (cnvflg(i)) then
1897-
if(k > kbcon1(i) .and. k < ktcon(i)) then
1910+
if(k >= kbcon1(i) .and. k < ktcon(i)) then
18981911
dp = 1000. * del(i,k)
18991912
tem = 0.5 * (omega_u(i,k) + omega_u(i,k-1))
19001913
omegac(i) = omegac(i) + tem * dp
@@ -2923,27 +2936,25 @@ subroutine samfdeepcnv_run (im,km,nn,first_time_step,restart, &
29232936
! compute convective turn-over time
29242937
!
29252938
!> - Following Bechtold et al. (2008) \cite bechtold_et_al_2008, the convective adjustment time (dtconv) is set to be proportional to the convective turnover time, which is computed using the mean updraft velocity (wc) and the cloud depth. It is also proportional to the grid size (gdx).
2926-
if(hwrf_samfdeep) then
2927-
do i= 1, im
2928-
if(cnvflg(i)) then
2929-
tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i))
2930-
dtconv(i) = tem / wc(i)
2931-
dtconv(i) = max(dtconv(i),dtmin)
2932-
dtconv(i) = min(dtconv(i),dtmax)
2933-
endif
2934-
enddo
2935-
else
2936-
do i= 1, im
2937-
if(cnvflg(i)) then
2938-
tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i))
2939-
dtconv(i) = tem / wc(i)
2940-
tfac = 1. + gdx(i) / 75000.
2941-
dtconv(i) = tfac * dtconv(i)
2942-
dtconv(i) = max(dtconv(i),dtmin)
2943-
dtconv(i) = min(dtconv(i),dtmax)
2944-
endif
2945-
enddo
2946-
endif
2939+
do i = 1, im
2940+
if (cnvflg(i)) then
2941+
tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i))
2942+
if (progomega) then
2943+
wc_eff = max(wc(i), wc_min)
2944+
dtconv(i) = tem / wc_eff
2945+
else
2946+
dtconv(i) = tem / wc(i)
2947+
endif
2948+
!grid spacing scaling (disabled for HWRF SAMF deep)
2949+
if (.not. hwrf_samfdeep) then
2950+
tfac = 1. + gdx(i) / 75000.
2951+
dtconv(i) = tfac * dtconv(i)
2952+
endif
2953+
!bounds
2954+
dtconv(i) = max(dtconv(i), dtmin)
2955+
dtconv(i) = min(dtconv(i), dtmax)
2956+
endif
2957+
enddo
29472958
!
29482959
!> - Calculate advective time scale (tauadv) using a mean cloud layer wind speed.
29492960
do i= 1, im

‎physics/CONV/SAMF/samfshalcnv.f‎

Lines changed: 43 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -174,7 +174,7 @@ subroutine samfshalcnv_run(im,km,nn,itc,ntc,cliq,cp,cvap, &
174174
& omegac(im),zeta(im,km),dbyo1(im,km),
175175
& sigmab(im),qadv(im,km)
176176
real(kind=kind_phys) gravinv,dxcrtas,invdelt,sigmind,sigmins,
177-
& sigminm
177+
& sigminm,wc_min,wc_eff
178178
logical flag_shallow,flag_mid
179179
c physical parameters
180180
! parameter(g=grav,asolfac=0.89)
@@ -208,8 +208,6 @@ subroutine samfshalcnv_run(im,km,nn,itc,ntc,cliq,cp,cvap, &
208208
! parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5)
209209
parameter(betaw=.03,dxcrtc0=9.e3)
210210
parameter(h1=0.33333333)
211-
! progsigma
212-
parameter(dxcrtas=500.e3,sigmind=0.01,sigmins=0.03,sigminm=0.01)
213211
c local variables and arrays
214212
real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km),
215213
& uo(im,km), vo(im,km), qeso(im,km),
@@ -296,8 +294,10 @@ subroutine samfshalcnv_run(im,km,nn,itc,ntc,cliq,cp,cvap, &
296294

297295
if (progsigma) then
298296
dxcrt=10.e3
297+
dxcrtas=500.e3
299298
else
300299
dxcrt=15.e3
300+
dxcrtas=500.e3
301301
endif
302302

303303
c-----------------------------------------------------------------------
@@ -315,7 +315,21 @@ subroutine samfshalcnv_run(im,km,nn,itc,ntc,cliq,cp,cvap, &
315315
prsl = prslp * 0.001
316316
del = delp * 0.001
317317
!************************************************************************
318-
!
318+
! - Initialize parameters related to prognostic closure
319+
if (progsigma) then
320+
if (progomega) then
321+
sigmind = 0.03
322+
sigmins = 0.03
323+
sigminm = 0.03
324+
wc_min = 0.2
325+
else
326+
sigmind = 0.01
327+
sigmins = 0.03
328+
sigminm = 0.03
329+
wc_min = 0.2
330+
endif
331+
endif
332+
!
319333
km1 = km - 1
320334
c
321335
c initialize arrays
@@ -1540,11 +1554,11 @@ subroutine samfshalcnv_run(im,km,nn,itc,ntc,cliq,cp,cvap, &
15401554
if (progomega) then
15411555
call progomega_calc(first_time_step,restart,im,km,
15421556
& kbcon1,ktcon,omegain,delt,del,zi,cnvflg,omegaout,
1543-
& grav,buo,drag,wush,xlamue,bb1,bb2)
1557+
& grav,buo,drag,wush,bb1,bb2)
15441558
do k = 1, km
15451559
do i = 1, im
15461560
if (cnvflg(i)) then
1547-
if(k > kbcon1(i) .and. k < ktcon(i)) then
1561+
if(k >= kbcon1(i) .and. k < ktcon(i)) then
15481562
omega_u(i,k)=omegaout(i,k)
15491563
omega_u(i,k)=MAX(omega_u(i,k),-80.)
15501564
! Convert to m/s for use in convective time-scale:
@@ -1631,7 +1645,7 @@ subroutine samfshalcnv_run(im,km,nn,itc,ntc,cliq,cp,cvap, &
16311645
do k = 2, km1
16321646
do i = 1, im
16331647
if (cnvflg(i)) then
1634-
if(k > kbcon1(i) .and. k < ktcon(i)) then
1648+
if(k >= kbcon1(i) .and. k < ktcon(i)) then
16351649
dp = 1000. * del(i,k)
16361650
tem = 0.5 * (omega_u(i,k) + omega_u(i,k-1))
16371651
omegac(i) = omegac(i) + tem * dp
@@ -1656,7 +1670,7 @@ subroutine samfshalcnv_run(im,km,nn,itc,ntc,cliq,cp,cvap, &
16561670
do k = 2, km1
16571671
do i = 1, im
16581672
if (cnvflg(i)) then
1659-
if(k > kbcon1(i) .and. k < ktcon(i)) then
1673+
if(k >= kbcon1(i) .and. k < ktcon(i)) then
16601674
if(omega_u(i,k) .ne. 0.)then
16611675
zeta(i,k)=eta(i,k)*(omegac(i)/omega_u(i,k))
16621676
else
@@ -1975,21 +1989,28 @@ subroutine samfshalcnv_run(im,km,nn,itc,ntc,cliq,cp,cvap, &
19751989
! compute convective turn-over time
19761990
!
19771991
!> - Following Bechtold et al. (2008) \cite bechtold_et_al_2008, calculate the convective turnover time using the mean updraft velocity (wc) and the cloud depth. It is also proportional to the grid size (gdx).
1978-
do i= 1, im
1979-
if(cnvflg(i)) then
1980-
tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i))
1981-
dtconv(i) = tem / wc(i)
1982-
if (.not.hwrf_samfshal) then
1983-
tfac = 1. + gdx(i) / 75000.
1984-
dtconv(i) = tfac * dtconv(i)
1985-
endif
1986-
dtconv(i) = max(dtconv(i),dtmin)
1987-
dtconv(i) = max(dtconv(i),dt2)
1988-
dtconv(i) = min(dtconv(i),dtmax)
1989-
endif
1992+
do i = 1, im
1993+
if (cnvflg(i)) then
1994+
tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i))
1995+
if (progomega) then
1996+
wc_eff = max(wc(i), wc_min)
1997+
dtconv(i) = tem / wc_eff
1998+
else
1999+
dtconv(i) = tem / wc(i)
2000+
endif
2001+
! - grid spacing scaling (disabled for HWRF shallow option)
2002+
if (.not. hwrf_samfshal) then
2003+
tfac = 1. + gdx(i) / 75000.
2004+
dtconv(i) = tfac * dtconv(i)
2005+
endif
2006+
! - limits
2007+
dtconv(i) = max(dtconv(i), dtmin)
2008+
dtconv(i) = max(dtconv(i),dt2)
2009+
dtconv(i) = min(dtconv(i), dtmax)
2010+
endif
19902011
enddo
1991-
!
1992-
!> - Calculate advective time scale (tauadv) using a mean cloud layer wind speed.
2012+
!
2013+
! > - Calculate advective time scale (tauadv) using a mean cloud layer wind speed.
19932014
do i= 1, im
19942015
if(cnvflg(i)) then
19952016
sumx(i) = 0.

‎physics/CONV/progomega_calc.f90‎

Lines changed: 19 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@ module progomega
2020
!!\section gen_progomega progomega_calc General Algorithm
2121

2222
subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegain,delt,del, &
23-
zi,cnvflg,omegaout,grav,buo,drag,wush,tentr,bb1,bb2)
23+
zi,cnvflg,omegaout,grav,buo,drag,wush,bb1,bb2)
2424

2525
use machine, only : kind_phys
2626
use funcphys, only : fpvs
@@ -30,20 +30,17 @@ subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegai
3030
integer, intent(in) :: kbcon1(im),ktcon(im)
3131
real(kind=kind_phys), intent(in) :: delt,grav,bb1,bb2
3232
real(kind=kind_phys), intent(in) :: omegain(im,km), del(im,km),zi(im,km)
33-
real(kind=kind_phys), intent(in) :: drag(im,km),buo(im,km),wush(im,km),tentr(im,km)
33+
real(kind=kind_phys), intent(in) :: drag(im,km),buo(im,km),wush(im,km)
3434
real(kind=kind_phys), intent(inout) :: omegaout(im,km)
3535
logical, intent(in) :: cnvflg(im),first_time_step,flag_restart
3636
real(kind=kind_phys) :: termA(im,km),termB(im,km),termC(im,km),omega(im,km)
3737
real(kind=kind_phys) :: RHS(im,km),Kd(im,km)
38-
real(kind=kind_phys) :: dp,dz,entrn,Kdn,discr,wush_pa,lbb1,lbb2,lbb3
38+
real(kind=kind_phys) :: dp,dz,discr,wush_pa,lbb1,lbb2,lbb3
3939
integer :: i,k
4040

41-
entrn = 0.8E-4 !0.5E-4 !m^-1
42-
Kdn = 0.5E-4 !2.9E-4 !m^-1
43-
lbb1 = 0.5 !1.0
44-
lbb2 = 3.2 !3.0
45-
lbb3 = 0.5 !0.5
46-
41+
lbb1 = 1.5
42+
lbb2 = 0.6
43+
lbb3 = 1.2
4744

4845
!Initialization 2D
4946
do k = 1,km
@@ -56,6 +53,16 @@ subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegai
5653
enddo
5754
enddo
5855

56+
do k = 1,km
57+
do i = 1,im
58+
if(cnvflg(i))then
59+
if(omega(i,k) < 1.0E-5) then
60+
omega(i,k) = 0.
61+
endif
62+
endif
63+
enddo
64+
enddo
65+
5966
if(first_time_step .and. .not. flag_restart)then
6067
do k = 1,km
6168
do i = 1,im
@@ -76,10 +83,7 @@ subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegai
7683
do k = 2, km
7784
do i = 1, im
7885
if (cnvflg(i)) then
79-
if (k > kbcon1(i) .and. k < ktcon(i)) then
80-
81-
! Aerodynamic drag parameter
82-
Kd(i,k) = (tentr(i,k)/entrn)*Kdn
86+
if (k >= kbcon1(i) .and. k < ktcon(i)) then
8387

8488
! Scale by dp/dz to have equation in Pa/s
8589
!(dp/dz > 0)
@@ -91,8 +95,8 @@ subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegai
9195
!termC - Adds buoyancy forcing
9296

9397
!Coefficients for the quadratic equation
94-
termA(i,k) = delt * ((lbb1 * drag(i,k) * (dp/dz)) + (Kd(i,k) * (dp/dz)))
95-
termB(i,k) = -1.0 - delt * lbb3 * wush(i,k) * dp/dz
98+
termA(i,k) = delt * ((lbb1 * drag(i,k) * (dp/dz)))
99+
termB(i,k) = 1.0 - delt * lbb3 * wush(i,k) * dp/dz
96100
termC(i,k) = omega(i,k) - delt * lbb2 * buo(i,k) * (dp/dz) &
97101
- delt * omega(i,k) * (omega(i,k-1) - omega(i,k)) / dp
98102
!Compute the discriminant

0 commit comments

Comments
 (0)