| blob | bf5829047e7a3efbc10254ac400de5b841409f40 |
1 MODULE module_cu_camzm_driver
3 !Roughly based on zm_conv_intr.F90 from CAM
5 USE module_cam_support, only: pcnst, pcols, pver, pverp
7 IMPLICIT NONE
9 PRIVATE !Default to private
10 PUBLIC :: & !Public entities
11 camzm_driver, &
12 zm_conv_init
14 CONTAINS
16 !------------------------------------------------------------------------
17 SUBROUTINE camzm_driver( &
18 ids,ide, jds,jde, kds,kde &
19 ,ims,ime, jms,jme, kms,kme &
20 ,its,ite, jts,jte, kts,kte &
21 ,itimestep, bl_pbl_physics, sf_sfclay_physics &
22 ,th, t_phy, tsk, tke_pbl, ust, qv, qc, qi &
23 ,mavail, kpbl, pblh, xland, z &
24 ,z_at_w, dz8w, ht, p, p8w, pi_phy, psfc &
25 ,u_phy, v_phy, hfx, qfx, cldfra &
26 ,tpert_camuwpbl &
27 ,dx, dt, stepcu, cudt, curr_secs, adapt_step_flag&
28 ,cudtacttime &
29 ,cape_out, mu_out, md_out, zmdt, zmdq &
30 ,rliq_out, dlf_out &
31 ,pconvb, pconvt, cubot, cutop, raincv, pratec &
32 ,rucuten, rvcuten &
33 ,rthcuten, rqvcuten, rqccuten, rqicuten &
34 ,evaptzm, fzsntzm, evsntzm, evapqzm, zmflxprc &
35 ,zmflxsnw, zmntprpd, zmntsnpd, zmeiheat &
36 ,cmfmc, cmfmcdzm, preccdzm, precz &
37 ,zmmtu, zmmtv, zmupgu, zmupgd, zmvpgu, zmvpgd &
38 ,zmicuu, zmicud, zmicvu, zmicvd &
39 ,zmdice, zmdliq &
40 )
41 ! This routine is based on zm_conv_tend in CAM. It handles the mapping of
42 ! variables from the WRF to the CAM framework for the Zhang-McFarlane
43 ! convective parameterization.
44 !
45 ! Author: William.Gustafson@pnl.gov, Nov. 2009
46 ! Last modified: William.Gustafson@pnl.gov, Nov. 2010
47 !------------------------------------------------------------------------
48 USE shr_kind_mod, only: r8 => shr_kind_r8
49 USE physconst, only: cpair, gravit
50 USE module_cu_camzm, only: convtran, momtran, zm_conv_evap, zm_convr
52 ! Subroutine arguments...
53 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
54 ims,ime, jms,jme, kms,kme, &
55 its,ite, jts,jte, kts,kte, &
56 bl_pbl_physics, & !pbl scheme
57 itimestep, & !time step index
58 sf_sfclay_physics, & !surface layer scheme
59 stepcu !number of time steps between Cu calls
61 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: &
62 cldfra, & !cloud fraction
63 dz8w, & !height between interfaces (m)
64 p, & !pressure at mid-level (Pa)
65 p8w, & !pressure at level interface (Pa)
66 pi_phy, & !exner function, (p0/p)^(R/cpair) (none)
67 qv, & !water vapor mixing ratio (kg/kg-dry air)
68 th, & !potential temperature (K)
69 tke_pbl, & !turbulent kinetic energy from PBL (m2/s2)
70 t_phy, & !temperature (K)
71 u_phy, & !zonal wind component on T points (m/s)
72 v_phy, & !meridional wind component on T points (m/s)
73 z, & !height above sea level at mid-level (m)
74 z_at_w !height above sea level at interface (m)
76 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN), OPTIONAL :: &
77 qc, & !cloud droplet mixing ratio (kg/kg-dry air)
78 qi !cloud ice crystal mixing ratio (kg/kg-dry air)
80 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) :: &
81 dlf_out, & !detraining cloud water tendendcy
82 evaptzm, & !temp. tendency - evap/snow prod from ZM (K/s)
83 fzsntzm, & !temp. tendency - rain to snow conversion from ZM (K/s)
84 evsntzm, & !temp. tendency - snow to rain conversion from ZM (K/s)
85 evapqzm, & !spec. humidity tend. - evaporation from ZM (kg/kg/s)
86 zmflxprc, & !flux of precipitation from ZM (kg/m2/s)
87 zmflxsnw, & !flux of snow from ZM (kg/m2/s)
88 zmntprpd, & !net precipitation production from ZM (kg/kg/s)
89 zmntsnpd, & !net snow production from ZM (kg/kg/s)
90 zmeiheat, & !heating by ice and evaporation from ZM (W/kg)
91 cmfmc, & !convective mass flux--m sub c, deep here but ultimately deep+shallow (kg/m2/s)
92 cmfmcdzm, & !convection mass flux from ZM deep (kg/m2/s)
93 md_out, & !output convective downdraft mass flux (kg/m2/s)
94 mu_out, & !output convective updraft mass flux (kg/m2/s)
95 rucuten, & !UNcoupled zonal wind tendency due to Cu scheme (m/s2)
96 rvcuten, & !UNcoupled meridional wind tendency due to Cu scheme (m/s2)
97 rthcuten, & !UNcoupled potential temperature tendendcy due to cu scheme (K/s)
98 rqvcuten, & !UNcoupled water vapor mixing ratio tendency due to Cu scheme (kg/kg/s)
99 zmdt, & !temp. tendency from moist convection (K/s)
100 zmdq, & !spec. humidity tendency from moist convection (kg/kg/s)
101 zmmtu, & !U tendency from ZM convective momentum transport (m/s2)
102 zmmtv, & !V tendency from ZM convective momentum transport (m/s2)
103 zmupgu, & !zonal force from ZM updraft pressure gradient term (m/s2)
104 zmupgd, & !zonal force from ZM downdraft pressure gradient term (m/s2)
105 zmvpgu, & !meridional force from ZM updraft pressure gradient term (m/s2)
106 zmvpgd, & !meridional force from ZM downdraft pressure gradient term (m/s2)
107 zmicuu, & !ZM in-cloud U updrafts (m/s)
108 zmicud, & !ZM in-cloud U downdrafts (m/s)
109 zmicvu, & !ZM in-cloud V updrafts (m/s)
110 zmicvd, & !ZM in-cloud V downdrafts (m/s)
111 zmdice, & !ZM cloud ice tendency (kg/kg/s)
112 zmdliq !ZM cloud liquid tendency (kg/kg/s)
114 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT), OPTIONAL :: &
115 rqccuten, & !UNcoupled cloud droplet mixing ratio tendency due to Cu scheme (kg/kg/s)
116 rqicuten !UNcoupled ice crystal mixing ratio tendency due to Cu scheme (kg/kg/s)
118 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: &
119 hfx, & !upward heat flux at sfc (W/m2)
120 ht, & !terrain height (m)
121 xland, & !land/water mask, 1=land, 2=water
122 mavail, & !soil moisture availability
123 pblh, & !planetary boundary layer height (m)
124 psfc, & !surface pressure (Pa)
125 qfx, & !upward moisture flux at sfc (kg/m2/s)
126 tpert_camuwpbl, & !temperature perturbation from UW CAM PBL
127 tsk, & !skin temperature (K)
128 ust !u* in similarity theory (m/s)
130 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: &
131 cape_out, & !convective available potential energy (J/kg)
132 cubot, & !level number of base of convection
133 cutop, & !level number of top of convection
134 pconvb, & !pressure of base of convection (Pa)
135 pconvt, & !pressure of top of convection (Pa)
136 pratec, & !rain rate returned to WRF for accumulation (mm/s)
137 preccdzm, & !convection precipitation rate from ZM deep (m/s)
138 precz, & !total precipitation rate from ZM (m/s)
139 raincv, & !time-step convective rain amount (mm)
140 rliq_out !vertical integral of reserved cloud water
142 REAL, INTENT(IN) :: &
143 cudt, & !cumulus time step (min)
144 curr_secs, & !current forecast time (s)
145 dt, & !domain time step (s)
146 dx !grid spacing (m)
148 REAL, INTENT (INOUT) :: &
149 cudtacttime !for adaptive time stepping, next to to run scheme (s)
151 INTEGER, DIMENSION( ims:ime, jms:jme), INTENT(IN) :: &
152 kpbl !index of PBL level
154 LOGICAL, INTENT(IN), OPTIONAL :: &
155 adapt_step_flag !using adaptive time steps?
157 ! Local variables...
158 !Variables dimensioned for input to ZM routines
159 REAL(r8), DIMENSION(pcols, kte+1) :: &
160 mcon, & !convective mass flux--m sub c (sub arg out in CAM)
161 pflx, & !scattered precip flux at each level (sub arg out in CAM)
162 pint8, & !pressure at layer interface (Pa)
163 zi8 !height above sea level at mid-level (m)
165 REAL(r8), DIMENSION(pcols, kte, pcnst) :: &
166 qh8 !specific humidity (kg/kg-moist air)
168 REAL(r8), DIMENSION(pcols, kte, 3) :: &
169 cloud, & !holder for cloud water and ice (q in CAM)
170 cloudtnd, & !holder for cloud tendencies (ptend_loc%q in CAM)
171 fracis !fraction of cloud species that are insoluble
173 REAL(r8), DIMENSION(pcols, kte, 2) :: &
174 icwu, & !in-cloud winds in upraft (m/s)
175 icwd, & !in-cloud winds in downdraft (m/s)
176 pguall, & !apparent force from updraft pres. gradient (~units?)
177 pgdall, & !apparent force from downdraft pres. gradient (~units?)
178 winds, & !wind components (m/s)
179 wind_tends !wind component tendencies (m/s2)
181 REAL(r8), DIMENSION(pcols, kte) :: &
182 cld8, & !cloud fraction
183 cme, & !cmf condensation - evaporation (~units?, sub arg out in CAM)
184 dlf, & !scattered version of the detraining cld h2o tendency (~units?)
185 fake_dpdry, & !place holder for dpdry, delta pres. of dry atmos.
186 flxprec, & !evap outfld: convective-scale flux of precip at interfaces (kg/m2/s)
187 flxsnow, & !evap outfld: convective-scale flux of snow at interfaces (kg/m2/s)
188 ntprprd, & !evap outfld: net precip production in layer (kg/kg/s)
189 ntsnprd, & !evap outfld: net snow production in layer (kg/kg/s)
190 pdel8, & !pressure thickness of layer (between interfaces, Pa)
191 pmid8, & !pressure at layer middle (Pa)
192 ql8, & !cloud liquid water (~units?)
193 qi8, & !cloud ice (~units?)
194 t8, & !temperature (K)
195 zm8, & !height above ground at mid-level (m)
196 qtnd, & !specific humidity tendency (kg/kg/s)
197 rprd, & !rain production rate (kg/kg/s, comes from pbuf array in CAM, add to restart?~)
198 stnd, & !heating rate (dry static energy tendency, W/kg)
199 tend_s_snwprd, & !heating rate of snow production (~units?)
200 tend_s_snwevmlt, & !heating rate of evap/melting of snow (~units?)
201 zdu !detraining mass flux (~units? sub arg out in CAM)
203 REAL(r8), DIMENSION(pcols) :: &
204 cape, & !convective available potential energy (J/kg)
205 jctop, & !row of top-of-deep-convection indices passed out (sub arg out in CAM)
206 jcbot, & !row of base of cloud indices passed out (sub arg out in CAM)
207 landfrac, & !land fraction
208 pblh8, & !planetary boundary layer height (m)
209 phis, & !geopotential at terrain height (m2/s2)
210 prec, & !convective-scale precipitation rate from ZM (m/s)
211 rliq, & !reserved liquid (not yet in cldliq) for energy integrals (units? sub arg out in CAM)
212 snow, & !convective-scale snowfall rate from ZM (m/s)
213 tpert !thermal (convective) temperature excess (K)
215 !Variables that were declared at the module level of zm_conv_intr in
216 !CAM. In that context they needed to be held between calls to
217 !zm_conv_tend and zm_conv_tend2 at the chunk level. In WRF, these vars
218 !are setup to hold the whole "memory" dimension, but as a 1-D vector
219 !instead of a 2-D array as is typically done in WRF. This allows us to
220 !leave the CAM routines in tact. For now, this forces us to immediately
221 !call zm_conv_tend2 before leaving this module, but it allows us to
222 !follow the WRF rules. We can deal with generalizing this for handling
223 !tracer convective transport of aerosols later.~
224 REAL(r8), DIMENSION(pcols, kte, (ime-ims+1)*(jme-jms+1)) :: &
225 dp, & !layer pres. thickness between interfaces (mb)
226 du, &
227 ed, &
228 eu, &
229 md, &
230 mu
232 REAL(r8), DIMENSION(pcols, (ime-ims+1)*(jme-jms+1)) :: &
233 dsubcld ! layer pres. thickness between LCL and maxi (mb)
235 INTEGER, DIMENSION(pcols, (ime-ims+1)*(jme-jms+1)) :: &
236 ideep, & ! holds position of gathered points
237 jt, & ! top-level index of deep cumulus convection
238 maxg ! gathered values of maxi
240 INTEGER, DIMENSION((ime-ims+1)*(jme-jms+1)) :: &
241 lengath ! number of gathered points
243 !Other local vars
244 INTEGER :: i, j, k, kflip, n, ncnst
245 INTEGER :: lchnk !chunk identifier, used to map 2-D to 1-D arrays in WRF
246 INTEGER :: ncol !number of atmospheric columns in chunk
247 LOGICAL, DIMENSION(3) :: l_qt !logical switches for constituent tendency presence
248 LOGICAL, DIMENSION(2) :: l_windt !logical switches for wind tendency presence
249 LOGICAL :: run_param , & !flag for handling alternate cumulus call freq.
250 doing_adapt_dt , decided
251 !
252 ! Check to see if this is a convection timestep...
253 !
255 ! Initialization for adaptive time step.
257 doing_adapt_dt = .FALSE.
258 IF ( PRESENT(adapt_step_flag) ) THEN
259 IF ( adapt_step_flag ) THEN
260 doing_adapt_dt = .TRUE.
261 IF ( cudtacttime .EQ. 0. ) THEN
262 cudtacttime = curr_secs + cudt*60.
263 END IF
264 END IF
265 END IF
267 ! Do we run through this scheme or not?
269 ! Test 1: If this is the initial model time, then yes.
270 ! ITIMESTEP=1
271 ! Test 2: If the user asked for the cumulus to be run every time step, then yes.
272 ! CUDT=0 or STEPCU=1
273 ! Test 3: If not adaptive dt, and this is on the requested cumulus frequency, then yes.
274 ! MOD(ITIMESTEP,STEPCU)=0
275 ! Test 4: If using adaptive dt and the current time is past the last requested activate cumulus time, then yes.
276 ! CURR_SECS >= CUDTACTTIME
278 ! If we do run through the scheme, we set the flag run_param to TRUE and we set the decided flag
279 ! to TRUE. The decided flag says that one of these tests was able to say "yes", run the scheme.
280 ! We only proceed to other tests if the previous tests all have left decided as FALSE.
282 ! If we set run_param to TRUE and this is adaptive time stepping, we set the time to the next
283 ! cumulus run.
285 decided = .FALSE.
286 run_param = .FALSE.
287 IF ( ( .NOT. decided ) .AND. &
288 ( itimestep .EQ. 1 ) ) THEN
289 run_param = .TRUE.
290 decided = .TRUE.
291 END IF
293 IF ( ( .NOT. decided ) .AND. &
294 ( ( cudt .EQ. 0. ) .OR. ( stepcu .EQ. 1 ) ) ) THEN
295 run_param = .TRUE.
296 decided = .TRUE.
297 END IF
299 IF ( ( .NOT. decided ) .AND. &
300 ( .NOT. doing_adapt_dt ) .AND. &
301 ( MOD(itimestep,stepcu) .EQ. 0 ) ) THEN
302 run_param = .TRUE.
303 decided = .TRUE.
304 END IF
306 IF ( ( .NOT. decided ) .AND. &
307 ( doing_adapt_dt ) .AND. &
308 ( curr_secs .GE. cudtacttime ) ) THEN
309 run_param = .TRUE.
310 decided = .TRUE.
311 cudtacttime = curr_secs + cudt*60
312 END IF
314 !Leave the subroutine if it is not yet time to call the cumulus param
315 if( .not. run_param ) return
316 !
317 ! Initialize...
318 !
319 ncol = 1 !chunk size in WRF is 1 since we loop over all columns in a tile
321 cape_out(its:ite, jts:jte) = 0.
322 mu_out(its:ite, kts:kte, jts:jte) = 0.
323 md_out(its:ite, kts:kte, jts:jte) = 0.
324 zmdt(its:ite, kts:kte, jts:jte) = 0.
325 !
326 ! Map variables to inputs for zm_convr and call it...
327 ! Loop over the points in the tile and treat them each as a CAM chunk.
328 !
329 do j = jts,jte
330 do i = its,ite
331 lchnk = (j-jts)*(ite-its+1) + (i-its+1) !1-D index location from the 2-D tile
333 !Flip variables on the layer middles
334 do k = kts,kte
335 kflip = kte-k+1
337 cld8(1,kflip) = cldfra(i,k,j)
338 pdel8(1,kflip) = p8w(i,k,j) - p8w(i,k+1,j)
339 pmid8(1,kflip) = p(i,k,j)
340 qh8(1,kflip,1) = max( qv(i,k,j)/(1.+qv(i,k,j)), 1e-30 ) !values of 0 cause a crash in entropy
341 if( present(qc) ) then
342 ql8(1,kflip) = qc(i,k,j)/(1.+qv(i,k,j)) !Convert to moist mix ratio
343 else
344 ql8(1,kflip) = 0.
345 end if
346 if( present(qi) ) then
347 qi8(1,kflip) = qi(i,k,j)/(1.+qv(i,k,j)) !Used in convtran, ditto for conversion
348 else
349 qi8(1,kflip) = 0.
350 end if
351 t8(1,kflip) = t_phy(i,k,j)
352 zm8(1,kflip) = z(i,k,j) - ht(i,j) !Height above the ground at midlevels
353 end do
355 !Flip variables on the layer interfaces
356 do k = kts,kte+1
357 kflip = kte-k+2
359 pint8(1,kflip) = p8w(i,k,j)
360 zi8(1,kflip) = z_at_w(i,k,j) -ht(i,j) !Height above the ground at interfaces
361 end do
363 !Other necessary conversions for input to ZM
364 if( xland(i,j)==2 ) then
365 landfrac(1) = 1. !land, WRF is all or nothing
366 else
367 landfrac(1) = 0. !water
368 end if
369 pblh8(1) = pblh(i,j)
370 phis(1) = ht(i,j)*gravit
372 call get_tpert(bl_pbl_physics, sf_sfclay_physics, dx, &
373 mavail(i,j), kpbl(i,j), pblh(i,j), &
374 dz8w(i,1,j), psfc(i,j), qv(i,1,j), t_phy(i,1,j), &
375 th(i,1,j), tsk(i,j), tke_pbl(i,:,j), ust(i,j), &
376 u_phy(i,1,j), v_phy(i,1,j), hfx(i,j), qfx(i,j), &
377 tpert_camuwpbl(i,j), kte, &
378 tpert(1))
380 !Call the Zhang-McFarlane (1996) convection parameterization
381 !NOTE: The 0.5*dt is correct and is a nuance of CAM typically
382 ! using 2*dt for physics tendencies. Everywhere in zm_convr
383 ! that dt is used, the dt is multiplied by 2 to get back to
384 ! the 2*dt. Everywhere else in the CAM ZM interface the full
385 ! 2*dt is passed into the subroutines. In WRF we use 1*dt
386 ! in place of CAM's 2*dt since the adjustment is made
387 ! elsewhere.
388 call zm_convr( lchnk, ncol, &
389 t8, qh8, prec, jctop, jcbot, &
390 pblh8, zm8, phis, zi8, qtnd, &
391 stnd, pmid8, pint8, pdel8, &
392 0.5_r8*real(dt,r8), mcon, cme, cape, &
393 tpert, dlf, pflx, zdu, rprd, &
394 mu(:,:,lchnk), md(:,:,lchnk),du(:,:,lchnk),eu(:,:,lchnk),ed(:,:,lchnk), &
395 dp(:,:,lchnk), dsubcld(:,lchnk), jt(:,lchnk), maxg(:,lchnk), ideep(:,lchnk), &
396 lengath(lchnk), ql8, rliq, landfrac )
398 !Start mapping CAM output to WRF output variables. We follow the
399 !same order as in CAM's zm_conv_tend to ease maintenance...
400 do k=kts,kte
401 kflip = kte-k+1
402 dlf_out(i,k,j) = dlf(1,kflip)
403 end do
404 cape_out(i,j) = cape(1)
405 rliq_out(i,j) = rliq(1)
407 !Convert mass flux from reported mb/s to kg/m2/s
408 mcon(:ncol,:pver) = mcon(:ncol,:pver) * 100._r8/gravit
410 ! Store upward and downward mass fluxes in un-gathered arrays
411 ! + convert from mb/s to kg/m2/s
412 do n=1,lengath(lchnk)
413 do k=kts,kte
414 kflip = kte-k+1
415 ! ii = ideep(n,lchnk) <--in WRF this is always 1 because chunk size=1
416 mu_out(i,k,j) = mu(n,kflip,lchnk) * 100._r8/gravit
417 md_out(i,k,j) = md(n,kflip,lchnk) * 100._r8/gravit
418 end do
419 end do
421 do k=kts,kte
422 kflip = kte-k+1
423 zmdt(i,k,j) = stnd(1,kflip)/cpair
424 zmdq(i,k,j) = qtnd(1,kflip)
425 end do
427 !Top and bottom pressure of convection
428 pconvt(i,j) = p8w(i,1,j)
429 pconvb(i,j) = p8w(i,1,j)
430 do n = 1,lengath(lchnk)
431 if (maxg(n,lchnk).gt.jt(n,lchnk)) then
432 pconvt(i,j) = pmid8(ideep(n,lchnk),jt(n,lchnk)) ! gathered array (or jctop ungathered)
433 pconvb(i,j) = pmid8(ideep(n,lchnk),maxg(n,lchnk))! gathered array
434 endif
435 end do
436 cutop(i,j) = jctop(1)
437 cubot(i,j) = jcbot(1)
439 !Add tendency from this process to tendencies arrays. Also,
440 !increment the local state arrays to reflect additional tendency
441 !from zm_convr, i.e. the physics update call in CAM. Note that
442 !we are not readjusting the pressure levels to hydrostatic
443 !balance for the new virtual temperature like is done in CAM. We
444 !will let WRF take care of such things later for the total
445 !tendency.
446 do k=kts,kte
447 kflip = kte-k+1
449 !Convert temperature to potential temperature and
450 !specific humidity to water vapor mixing ratio
451 rthcuten(i,k,j) = zmdt(i,k,j)/pi_phy(i,k,j)
452 rqvcuten(i,k,j) = zmdq(i,k,j)/(1._r8 - zmdq(i,k,j))
454 t8(1,kflip) = t8(1,kflip) + zmdt(i,k,j)*real(dt,r8)
455 qh8(1,kflip,1) = qh8(1,kflip,1) + zmdq(i,k,j)*real(dt,r8)
456 end do
458 ! Determine the phase of the precipitation produced and add latent heat of fusion
459 ! Evaporate some of the precip directly into the environment (Sundqvist)
460 ! Allow this to use the updated state1 (t8 and qh8 in WRF) and the fresh ptend_loc type
461 ! heating and specific humidity tendencies produced
462 qtnd = 0._r8 !re-initialize tendencies (i.e. physics_ptend_init(ptend_loc))
463 stnd = 0._r8
464 call zm_conv_evap(ncol, lchnk, &
465 t8, pmid8, pdel8, qh8(:,:,1), &
466 stnd, tend_s_snwprd, tend_s_snwevmlt, qtnd, &
467 rprd, cld8, real(dt,r8), &
468 prec, snow, ntprprd, ntsnprd , flxprec, flxsnow)
470 ! Parse output variables from zm_conv_evap
471 do k=kts,kte
472 kflip = kte-k+1
474 evaptzm(i,k,j) = stnd(1,kflip)/cpair
475 fzsntzm(i,k,j) = tend_s_snwprd(1,kflip)/cpair
476 evsntzm(i,k,j) = tend_s_snwevmlt(1,kflip)/cpair
477 evapqzm(i,k,j) = qtnd(1,kflip)
478 zmflxprc(i,k,j) = flxprec(1,kflip)
479 zmflxsnw(i,k,j) = flxsnow(1,kflip)
480 zmntprpd(i,k,j) = ntprprd(1,kflip)
481 zmntsnpd(i,k,j) = ntsnprd(1,kflip)
482 zmeiheat(i,k,j) = stnd(1,kflip) !Do we really need this and evaptzm?
483 cmfmc(i,k,j) = mcon(1,kflip) !Set to deep value here, shallow added in UW scheme
484 cmfmcdzm(i,k,j) = mcon(1,kflip)
485 preccdzm(i,j) = prec(1) !Rain rate from just deep
486 precz(i,j) = prec(1) !Rain rate for total convection (just deep right now)
487 pratec(i,j) = prec(1)*1e3 !Rain rate used in WRF for accumulation (mm/s)
488 raincv(i,j) = pratec(i,j)*dt !Rain amount for time step returned back to WRF
489 end do
491 !Add tendency from zm_conv_evap to tendencies arrays. Also,
492 !increment the local state arrays to reflect additional tendency
493 !Note that we are not readjusting the pressure levels to hydrostatic
494 !balance for the new virtual temperature like is done in CAM. We
495 !will let WRF take care of such things later for the total
496 !tendency.
497 do k=kts,kte
498 kflip = kte-k+1
500 !Convert temperature to potential temperature and
501 !specific humidity to water vapor mixing ratio
502 rthcuten(i,k,j) = rthcuten(i,k,j) + &
503 evaptzm(i,k,j)/pi_phy(i,k,j)
504 rqvcuten(i,k,j) = rqvcuten(i,k,j) + &
505 evapqzm(i,k,j)/(1. - qv(i,k,j))
507 t8(1,kflip) = t8(1,kflip) + evaptzm(i,k,j)*real(dt,r8)
508 qh8(1,kflip,1) = qh8(1,kflip,1) + evapqzm(i,k,j)*real(dt,r8)
509 end do
511 ! Momentum transport
512 stnd = 0._r8 !Zero relevant tendencies in preparation
513 wind_tends = 0._r8
514 do k=kts,kte
515 kflip = kte-k+1
516 winds(1,k,1) = u_phy(i,kflip,j)
517 winds(1,k,2) = v_phy(i,kflip,j)
518 end do
519 l_windt(1:2) = .true.
521 call momtran (lchnk, ncol, &
522 l_windt, winds, 2, mu(:,:,lchnk), md(:,:,lchnk), &
523 du(:,:,lchnk), eu(:,:,lchnk), ed(:,:,lchnk), dp(:,:,lchnk), dsubcld(:,lchnk), &
524 jt(:,lchnk),maxg(:,lchnk), ideep(:,lchnk), 1, lengath(lchnk), &
525 itimestep, wind_tends, pguall, pgdall, icwu, icwd, real(dt,r8), stnd )
527 !Add tendency from momtran to tendencies arrays. Also,
528 !increment the local state arrays to reflect additional tendency
529 !Note that we are not readjusting the pressure levels to hydrostatic
530 !balance for the new virtual temperature like is done in CAM. We
531 !will let WRF take care of such things later for the total
532 !tendency.
533 do k=kts,kte
534 kflip = kte-k+1
536 !Convert temperature to potential temperature and
537 !specific humidity to water vapor mixing ratio
538 rucuten(i,k,j) = wind_tends(1,kflip,1)
539 rvcuten(i,k,j) = wind_tends(1,kflip,2)
540 rthcuten(i,k,j) = rthcuten(i,k,j) + &
541 stnd(1,kflip)/cpair/pi_phy(i,k,j)
543 t8(1,kflip) = t8(1,kflip) + stnd(1,kflip)/cpair*real(dt,r8)
544 !winds is not used again so do not bother updating them
545 end do
547 !Parse output arrays for momtran
548 do k=kts,kte
549 kflip = kte-k+1
551 zmmtu(i,k,j) = wind_tends(1,kflip,1) !wind tendencies...
552 zmmtv(i,k,j) = wind_tends(1,kflip,2)
554 zmupgu(i,k,j) = pguall(1,kflip,1) !apparent force pres. grads...
555 zmupgd(i,k,j) = pgdall(1,kflip,1)
556 zmvpgu(i,k,j) = pguall(1,kflip,2)
557 zmvpgd(i,k,j) = pgdall(1,kflip,2)
559 zmicuu(i,k,j) = icwu(1,kflip,1) !in-cloud winds...
560 zmicud(i,k,j) = icwd(1,kflip,1)
561 zmicvu(i,k,j) = icwu(1,kflip,2)
562 zmicvd(i,k,j) = icwd(1,kflip,2)
563 end do
565 !Setup for convective transport of cloud water and ice
566 !~We should set this up to use an integer pointer instead of
567 ! hard-coded numbers for each species so that we can easily
568 ! handle other species. Something to come back to later...
569 l_qt(1) = .false. !do not mix water vapor
570 l_qt(2:3) = .true. !do mix cloud water and ice
571 cloudtnd = 0._r8
572 fracis(1,:,1:3) = 0._r8 !all cloud liquid & ice is soluble
573 ncnst = 3 !number of constituents in cloud array (including vapor)
574 fake_dpdry = 0._r8 !delta pres. for dry atmos. (0 since assuming moist mix ratios)
575 do k=kts,kte
576 kflip = kte-k+1
578 cloud(1,kflip,1) = qh8(1,kflip,1) !We can either use moist mix ratios, as is
579 cloud(1,kflip,2) = ql8(1,kflip) !done here, or else use dry mix ratios, send
580 cloud(1,kflip,3) = qi8(1,kflip) !in appropriate dpdry values, and return the
581 !approp. response from cnst_get_type_byind
582 end do
584 call convtran (lchnk, &
585 l_qt, cloud, ncnst, mu(:,:,lchnk), md(:,:,lchnk), &
586 du(:,:,lchnk), eu(:,:,lchnk), ed(:,:,lchnk), dp(:,:,lchnk), dsubcld(:,lchnk), &
587 jt(:,lchnk), maxg(:,lchnk), ideep(:,lchnk), 1, lengath(lchnk), &
588 itimestep, fracis, cloudtnd, fake_dpdry)
590 !Parse output for convtran and increment tendencies
591 do k=kts,kte
592 kflip = kte-k+1
594 zmdice(i,k,j) = cloudtnd(1,kflip,3)
595 zmdliq(i,k,j) = cloudtnd(1,kflip,2)
597 !Convert cloud tendencies from wet to dry mix ratios
598 if( present(rqccuten) ) then
599 rqccuten(i,k,j) = cloudtnd(1,kflip,2)/(1. - qv(i,k,j))
600 end if
601 if( present(rqicuten) ) then
602 rqicuten(i,k,j) = cloudtnd(1,kflip,3)/(1. - qv(i,k,j))
603 end if
604 end do
606 end do !i-loop
607 end do !j-loop
608 END SUBROUTINE camzm_driver
611 !------------------------------------------------------------------------
612 SUBROUTINE zm_conv_init(rucuten, rvcuten, rthcuten, rqvcuten, &
613 rqccuten, rqicuten, &
614 p_qc, p_qi, param_first_scalar, &
615 restart, &
616 ids, ide, jds, jde, kds, kde, &
617 ims, ime, jms, jme, kms, kme, &
618 its, ite, jts, jte, kts, kte )
619 ! Initialization routine for Zhang-McFarlane cumulus parameterization
620 ! from CAM. The routine with this name in CAM is revamped here to give
621 ! the needed functionality within WRF.
622 !
623 ! Author: William.Gustafson@pnl.gov, Nov. 2009
624 !------------------------------------------------------------------------
625 USE module_cam_esinti, only: esinti
626 USE physconst, only: epsilo, latvap, latice, rh2o, cpair, tmelt
627 USE module_bl_camuwpbl_driver, only: vd_register
628 USE module_cu_camzm, only: zm_convi, zmconv_readnl
630 LOGICAL , INTENT(IN) :: restart
631 INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
632 ims, ime, jms, jme, kms, kme, &
633 its, ite, jts, jte, kts, kte, &
634 p_qc, p_qi, param_first_scalar
636 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(OUT) :: &
637 rucuten, &
638 rvcuten, &
639 rthcuten, &
640 rqvcuten, &
641 rqccuten, &
642 rqicuten
644 integer :: i, itf, j, jtf, k, ktf
645 integer :: limcnv
647 jtf = min(jte,jde-1)
648 ktf = min(kte,kde-1)
649 itf = min(ite,ide-1)
650 !
651 ! Initialize module_cam_support variables...
652 ! This could be moved to a master "cam-init" subroutine once multiple CAM
653 ! parameterizations are implemented in WRF. For now, it doesn't hurt to
654 ! have these possibly initialized more than once since they are
655 ! essentially constant.
656 !
657 pver = ktf-kts+1
658 pverp = pver+1
659 !
660 ! Initialize the saturation vapor pressure look-up table...
661 ! This could be moved to a master cam-init subroutine too if it is needed
662 ! by more than one CAM parameterization. In CAM this is called from
663 ! phys_init.
664 !
665 call esinti(epsilo, latvap, latice, rh2o, cpair, tmelt)
667 !~Need code here to set limcnv to max convective level of 40 mb... To
668 ! properly set an average value for the whole domain would involve doing
669 ! a collective operation across the whole domain to get pressure averages
670 ! by level, which is difficult within the WRF framework. So, we will just
671 ! use the model top for now.
672 !
673 ! Technically, limcnv should be calculated for each grid point at each
674 ! time because the levels in WRF do not have a constant pressure, as
675 ! opposed to CAM. But, that would involve making this variable local
676 ! instead of at the module level as is currently done in CAM. For now,
677 ! we will not worry about this because WRF rarely has domains higher than
678 ! 40 mb. There is also little variability between grid points near the
679 ! model top due to the underlying topography so a constant value would
680 ! be OK across the comain. Also, remember that CAM levels are reversed in
681 ! the vertical from WRF. So, setting limcnv to 1 means the top of the
682 ! domain. Note that because of a bug in the parcel_dilute routine, limcnv
683 ! must be >=2 instead of 1. Otherwise, an array out-of-bounds occurs.
684 limcnv = 2
685 !
686 ! Get the ZM namelist variables (hard-wired for now)...
687 !
688 call zmconv_readnl("hard-wired")
689 !
690 !~need to determine if convection should happen in PBL and set
691 ! no_deep_pbl_in accordingly
692 call zm_convi(limcnv, NO_DEEP_PBL_IN=.false.)
694 !
695 ! Set initial values for arrays dependent on restart conditions
696 !
697 if(.not.restart)then
698 do j=jts,jtf
699 do k=kts,ktf
700 do i=its,itf
701 rucuten(i,k,j) = 0.
702 rvcuten(i,k,j) = 0.
703 rthcuten(i,k,j) = 0.
704 rqvcuten(i,k,j) = 0.
705 if( p_qc > param_first_scalar ) rqccuten(i,k,j) = 0.
706 if( p_qi > param_first_scalar ) rqicuten(i,k,j) = 0.
707 enddo
708 enddo
709 enddo
710 end if
712 END SUBROUTINE zm_conv_init
715 !------------------------------------------------------------------------
716 SUBROUTINE get_tpert(bl_pbl_physics, sf_sfclay_physics, dx, &
717 mavail, kpbl, pblh, dzlowest, &
718 psfc, qvlowest, t_phylowest, thlowest, tsk, tke_pbl, ust, &
719 u_phylowest, v_phylowest, hfx, qfx, tpert_camuwpbl, kte, tpert)
720 ! Calculates the temperature perturbation used to trigger convection.
721 ! This should only be used if tpert is not already provided by CAM's PBL
722 ! scheme. Right now, this only works with the Mellor-Yamada boundary
723 ! layer scheme in combination with the MYJ surface scheme. This is due to
724 ! the need of TKE for the vertical velocity perturbation. This scheme has
725 ! not been generalized to handle all schemes that produce TKE at this
726 ! time, and the non-TKE schemes, e.g. YSU, could probably have an
727 ! appropriate tpert deduced but we have not taken the time yet to do it.
728 !
729 ! Author: William.Gustafson@pnl.gov, Nov. 2009
730 ! Last updated: William.Gustafson@pnl.gov, Nov. 2010
731 !------------------------------------------------------------------------
732 USE shr_kind_mod, only: r8 => shr_kind_r8
733 USE module_model_constants, only: cp, ep_1, ep_2, g, r_d, rcp, &
734 svp1, svp2, svp3, svpt0, xlv
735 USE module_state_description, ONLY : SFCLAYSCHEME &
736 ,MYJSFCSCHEME &
737 ,GFSSFCSCHEME &
738 ,SLABSCHEME &
739 ,LSMSCHEME &
740 ,RUCLSMSCHEME &
741 ,MYJPBLSCHEME &
742 ,CAMUWPBLSCHEME
743 !
744 ! Subroutine arguments...
745 !
746 real, dimension(:), intent(in) :: tke_pbl
747 real, intent(in) :: dx, dzlowest, hfx, mavail, pblh, psfc, qvlowest, &
748 tpert_camuwpbl, tsk, t_phylowest, thlowest, ust, u_phylowest, &
749 v_phylowest
750 integer, intent(in) :: bl_pbl_physics, kpbl, kte, sf_sfclay_physics
751 real(r8),intent(inout) :: tpert
752 !
753 ! Local vars...
754 !
755 real, parameter :: fak = 8.5 !Constant in surface temperature excess
756 real, parameter :: tfac = 1. !Ratio of 'tpert' to (w't')/wpert
757 real, parameter :: wfac = 1. !Ratio of 'wpert' to sqrt(tke)
758 real, parameter :: wpertmin = 1.e-6 !Min PBL eddy vertical velocity perturbation
759 real :: ebrk !In CAM, net mean TKE of CL including
760 !entrainment effect (m2/s2). In WRF,
761 !average TKE within the PBL
762 real :: br2, dthvdz, e1, flux, govrth, psfccmb, qfx, qsfc, rhox, thgb, &
763 thv, tskv, tvcon, vconv, vsgd, wpert, wspd, za
764 integer :: k
765 character(len=250) :: msg
766 logical :: UnstableOrNeutral
767 !
768 ! We can get the perturbation values directly from CAMUWPBLSCHEME if
769 ! available. Otherwise, we have to calculate them.
770 !
771 if( bl_pbl_physics==CAMUWPBLSCHEME ) then
772 tpert = tpert_camuwpbl
773 !
774 !...non-CAMUWPBL. Need MYJ SFC & PBL for now until other schemes
775 ! get coded to give perturbations too.
776 ! First, we need to determine if the conditions are stable or unstable.
777 ! We will do this by mimicing the bulk Richardson calculation from
778 ! module_sf_sfclay.F because the MYJ scheme does not provide this info.
779 ! The logic borrowed from the CuP shallow cumulus scheme. Non-MYJ sfc
780 ! scheme code is commented out for now since we also require MYJ PBL
781 ! scheme.
782 !
783 elseif( bl_pbl_physics==MYJPBLSCHEME ) then
785 UnstableOrNeutral = .false.
786 sfclay_case: SELECT CASE (sf_sfclay_physics)
787 CASE (MYJSFCSCHEME)
788 ! The MYJ sfc scheme does not already provide a stability criteria.
789 ! So, we will mimic the bulk Richardson calculation from
790 ! module_sf_sfclay.F.
792 if( pblh <= 0. ) call wrf_error_fatal( &
793 "CAMZMSCHEME needs a PBL height from a PBL scheme.")
795 za = 0.5*dzlowest
797 e1 = svp1*exp(svp2*(tsk-svpt0)/(tsk-svp3))
798 psfccmb=psfc/1000. !converts from Pa to cmb
799 qsfc = ep_2*e1/(psfccmb-e1)
800 thgb = tsk*(100./psfccmb)**rcp
801 tskv = thgb*(1.+ep_1*qsfc*mavail)
802 tvcon = 1.+ep_1*qvlowest
803 thv = thlowest*tvcon
804 dthvdz = (thv-tskv)
806 govrth = g/thlowest
808 rhox = psfc/(r_d*t_phylowest*tvcon)
809 flux = max(hfx/rhox/cp + ep_1*tskv*qfx/rhox,0.)
810 vconv = (g/tsk*pblh*flux)**.33
811 vsgd = 0.32 * (max(dx/5000.-1.,0.))**.33
812 wspd = sqrt(u_phylowest*u_phylowest+v_phylowest*v_phylowest)
813 wspd = sqrt(wspd*wspd+vconv*vconv+vsgd*vsgd)
814 wspd = max(wspd,0.1)
816 !And finally, the bulk Richardson number...
817 br2 = govrth*za*dthvdz/(wspd*wspd)
819 if( br2 <= 0. ) UnstableOrNeutral = .true.
821 CASE DEFAULT
822 call wrf_error_fatal("CAMZMSCHEME requires MYJSFCSCHEME or else CAMUWPBLSCHEME.")
824 END SELECT sfclay_case
825 !
826 ! The perturbation temperature for unstable conditions...
827 ! The calculation follows the one in caleddy inside eddy_diff.F90 from
828 ! CAM.
829 !
830 !Check that we are using the MJY BL scheme since we are hard-wired to
831 !use TKE and u* from this scheme. At some point this dependency should
832 !be broken and a way needs to be found for other schemes to provide
833 !reasonable tpert values too.
834 if( bl_pbl_physics /= MYJPBLSCHEME ) &
835 call wrf_error_fatal("CAMZMSCHEME requires MYJPBLSCHEME or CAMUWPBLSCHEME")
838 !Recalculate rhox in case a non-MYJ scheme is used to get
839 !stability and rhox is not calculated above. Right now, this is
840 !technically reduncant, but cheap.
841 tvcon = 1.+ep_1*qvlowest
842 rhox = psfc/(r_d*t_phylowest*tvcon)
844 if( UnstableOrNeutral ) then
845 !1st, get avg TKE w/n the PBL as a proxy for ebrk variable in CAM
846 ebrk = 0.
847 do k=1,min(kpbl,kte)
848 ebrk = ebrk + tke_pbl(k)
849 end do
850 ebrk = ebrk/real(kpbl)
852 wpert = max( wfac*sqrt(ebrk), wpertmin )
853 tpert = max( abs(hfx/rhox/cp)*tfac/wpert, 0. )
854 !
855 ! Or, the perturbation temperature for stable conditions...
856 !
857 else !Stable conditions
858 tpert = max( hfx/rhox/cp*fak/ust, 0. )
859 end if
861 else
862 call wrf_error_fatal("CAMZMSCHEME requires MYJPBLSCHEME or CAMUWPBLSCHEME")
864 end if !PBL choice
866 END SUBROUTINE get_tpert
868 END MODULE module_cu_camzm_driver