| blob | c25ffb26cb0eeb06739908cb9fa04ed137c92f18 |
1 MODULE module_cu_camuwshcu_driver
2 USE shr_kind_mod, only: r8 => shr_kind_r8
4 ! Roughly based on convect_shallow_tend in convect_shallow.F90 from CAM
5 ! but tailored for the UW shallow cumulus scheme.
7 IMPLICIT NONE
9 PRIVATE !Default to private
10 PUBLIC :: & !Public entities
11 camuwshcu_driver
13 CONTAINS
15 !------------------------------------------------------------------------
16 SUBROUTINE camuwshcu_driver( &
17 ids,ide, jds,jde, kds,kde &
18 ,ims,ime, jms,jme, kms,kme &
19 ,its,ite, jts,jte, kts,kte &
20 ,num_moist, dt &
21 ,p, p8w, pi_phy, z, z_at_w, dz8w &
22 ,t_phy, u_phy, v_phy &
23 ,moist, qv, qc, qi &
24 ,pblh_in, tke_pbl, cldfra, cldfra_old, cldfrash &
25 ,cush_inout, rainsh, pratesh, snowsh, icwmrsh &
26 ,cmfmc, cmfmc2_inout, rprdsh_inout, cbmf_inout &
27 ,cmfsl, cmflq, dlf, evapcsh_inout &
28 ,rliq, rliq2_inout, cubot, cutop &
29 ,rushten, rvshten, rthshten &
30 ,rqvshten, rqcshten, rqrshten &
31 ,rqishten, rqsshten, rqgshten &
32 ,ht &
33 )
34 ! This routine is based on convect_shallow_tend in CAM. It handles the
35 ! mapping of variables from the WRF to the CAM framework for the UW
36 ! shallow convective parameterization.
37 !
38 ! Author: William.Gustafson@pnl.gov, Jan. 2010
39 !------------------------------------------------------------------------
40 USE module_state_description, only: param_first_scalar, &
41 p_qc, p_qr, p_qi, p_qs, p_qg
42 USE module_cam_support, only: pcols, pver
43 USE constituents, only: cnst_get_ind
44 USE physconst, only: cpair, gravit, latvap
45 USE uwshcu, only: compute_uwshcu_inv
46 USE wv_saturation, only: fqsatd
48 ! Subroutine arguments...
49 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
50 ims,ime, jms,jme, kms,kme, &
51 its,ite, jts,jte, kts,kte, &
52 num_moist
54 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ), INTENT(IN) :: &
55 moist !moist tracer array
57 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: &
58 cldfra, & !cloud fraction
59 cldfra_old, & !previous time step cloud fraction
60 dz8w, & !height between layer interface (m)
61 p, & !pressure at mid-level (Pa)
62 p8w, & !pressure at level interface (Pa)
63 pi_phy, & !exner function, (p0/p)^(R/cpair) (none)
64 qv, & !water vapor mixing ratio (kg/kg-dry air)
65 t_phy, & !temperature (K)
66 tke_pbl, & !turbulent kinetic energy from PBL (m2/s2)
67 u_phy, & !zonal wind component on T points (m/s)
68 v_phy, & !meridional wind component on T points (m/s)
69 z, & !height above sea level at mid-level (m)
70 z_at_w !height above sea level at interface (m)
72 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN), OPTIONAL :: &
73 qc, & !cloud droplet mixing ratio (kg/kg-dry air)
74 qi !cloud ice crystal mixing ratio (kg/kg-dry air)
76 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: &
77 pblh_in, & !height of PBL (m)
78 ht !Terrain height (m)
80 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) :: &
81 cldfrash, & !shallow convective cloud fraction
82 cmfmc, & !deep+shalow cloud fraction (already contains deep part from ZM)
83 cmfmc2_inout, & !shallow cloud fraction
84 cmflq, & !convective flux of total water in energy unit (~units)
85 cmfsl, & !convective flux of liquid water static energy (~units)
86 dlf, & !dq/dt due to export of cloud water (input=deep from ZM, output=deep+shallow)
87 evapcsh_inout, & !output array for evaporation of shallow convection precipitation (kg/kg/s)
88 icwmrsh, & !shallow cumulus in-cloud water mixing ratio (kg/m2)
89 rprdsh_inout, & !dq/dt due to deep(~?) & shallow convective rainout (~units?)
90 rushten, & !UNcoupled zonal wind tend from shallow Cu scheme (m/s2)
91 rvshten, & !UNcoupled meridional wind tend from shallow Cu scheme (m/s2)
92 rthshten, & !UNcoupled potential temperature tendendcy from shallow cu scheme (K/s)
93 rqvshten, & !UNcoupled water vapor mixing ratio tend from shallow Cu scheme (kg/kg/s)
94 rqcshten, & !UNcoupled clod droplet mixing ratio tend from shallow Cu scheme (kg/kg/s)
95 rqrshten, & !UNcoupled raindrop mixing ratio tend from shallow Cu scheme (kg/kg/s)
96 rqishten, & !UNcoupled ice crystal mixing ratio tend from shallow Cu scheme (kg/kg/s)
97 rqsshten, & !UNcoupled snow mixing ratio tend from shallow Cu scheme (kg/kg/s)
98 rqgshten !UNcoupled graupel mixing ratio tend from shallow Cu scheme (kg/kg/s)
100 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: &
101 cbmf_inout, & !cloud base mass flux (kg/m2/s)
102 cubot, & !level number of base of convection
103 cutop, & !level number of top of convection
104 cush_inout, & !convective scale height (~units?)
105 pratesh, & !time-step shallow cumulus precip rate at surface (mm/s)
106 rainsh, & !time-step shallow cumulus precip (rain+snow) at surface (mm)
107 rliq, & !vertically-integrated reserved cloud condensate (m/s)
108 rliq2_inout, & !vertically-integrated reserved cloud condensate for shallow (m/s)
109 snowsh !accumulated convective snow rate at surface for shallow Cu (m/s) ~are these the units we should use for WRF?
111 REAL, INTENT(IN) :: &
112 dt !time step (s)
114 ! Local variables...
115 !Variables dimensioned for input to CAM routines
116 REAL(r8), DIMENSION(pcols, kte, 5) :: & ! dimensioned by 5(=ncnst) associated with CAM constituents (cnst_name size)
117 moist8, & !tracer array for CAM routines
118 tnd_tracer !tracer tendency
120 REAL(r8), DIMENSION(pcols, kte+1) :: &
121 pint8, & !pressure at layer interface (Pa)
122 zi8, & !height above the ground at interfaces (m)
123 tke8, & !turbulent kinetic energy at level interfaces (m2/s2)
124 slflx, & !convective liquid water static energy flux (~units?)
125 qtflx !convective total water flux (~units?)
129 REAL(r8), DIMENSION(pcols, kte) :: &
130 cld8, & !cloud fraction
131 cldold8, & !previous time step cloud fraction ~should this be just the convective part?
132 cmfdqs, & !convective snow production (~units?)
133 cmfmc2, & !cloud fraction
134 evapcsh, & !evaporation of shallow convection precipitation >= 0. (kg/kg/s)
135 iccmr_uw, & !in-cloud cumulus LWC+IWC (kg/m2)
136 icwmr_uw, & !in-cloud cumulus LWC (kg/m2)
137 icimr_uw, & !in-cloud cumulus IWC (kg/m2)
138 pdel8, & !pressure difference between layer interfaces (Pa)
139 pdeldry8, & !pressure difference between layer interfaces for dry atm (Pa)
140 pmid8, & !pressure at layer middle (Pa)
141 qc2, & !dq/dt due to export of cloud water
142 qh8, & !specific humidity (kg/kg-moist air)
143 qc8, & !cloud liquid water (~units?)
144 qi8, & !cloud ice (~units?)
145 qhtnd, & !specific humidity tendency (kg/kg/s)
146 qctnd, & !cloud water mixing ratio tendency
147 qitnd, & !cloud ice mixing ratio tendency
148 rprdsh, & !dq/dt due to deep(~?) & shallow convective rainout (~units?)
149 s8, & !dry static energy (J/kg)
150 shfrc, & !shallow cloud fraction
151 stnd, & !heating rate (dry static energy tendency, W/kg)
152 t8, & !temperature (K)
153 u8, & !environment zonal wind (m/s)
154 utnd, & !zonal wind tendency (m/s2)
155 v8, & !environment meridional wind (m/s)
156 vtnd, & !meridional wind tendency (m/s2)
157 zm8 !height between interfaces (m)
159 REAL(r8), DIMENSION(pcols) :: &
160 cbmf, & !cloud base mass flux (kg/m2/s)
161 cnb2, & !bottom level of convective activity
162 cnt2, & !top level of convective activity
163 cush, & !convective scale height (~units?)
164 pblh, & !pblh height (m)
165 precc, & !convective precip (rain+snow) at surface for shallow Cu (m/s)
166 rliq2, & !vertically-integrated reserved cloud condensate for shallow (m/s)
167 snow !convective snow rate at surface (m/s)
169 !Other local vars
170 REAL(r8) :: ztodt !2 delta-t (s)
171 INTEGER :: i, j, k, kflip, m, mp1
172 INTEGER :: cnb, cnt !index of cloud base and top in CAM world (indices decrease with height)
173 INTEGER :: lchnk !chunk identifier, used to map 2-D to 1-D arrays in WRF
174 INTEGER :: ncnst !number of tracers
175 INTEGER :: ncol !number of atmospheric columns in chunk
176 CHARACTER(LEN=250) :: msg
177 !
178 ! Initialize...
179 !
180 ncol = 1 !chunk size in WRF is 1 since we loop over all columns in a tile
181 ! ncnst = num_moist-1 !currently we only handle the moist array for tracers
182 ncnst = 5 ! this is associated with moist8 array
183 ztodt = 2.*dt
184 !
185 ! Map variables to inputs for zm_convr and call it...
186 ! Loop over the points in the tile and treat them each as a CAM chunk.
187 !
188 ij_loops : do j = jts,jte
189 do i = its,ite
190 lchnk = (j-jts)*(ite-its+1) + (i-its+1) !1-D index location from the 2-D tile
192 !Flip variables on the layer interfaces
193 do k = kts,kte+1
194 kflip = kte-k+2
196 pint8(1,kflip) = p8w(i,k,j)
197 zi8(1,kflip) = z_at_w(i,k,j) - ht(i,j) ! height above the ground at interfaces
198 end do
200 !Flip variables on the layer middles
201 do k = kts,kte
202 kflip = kte-k+1
204 cld8(1,kflip) = cldfra(i,k,j)
205 cldold8(1,kflip) = cldfra_old(i,k,j)
206 pdel8(1,kflip) = p8w(i,k,j) - p8w(i,k+1,j)
207 pmid8(1,kflip) = p(i,k,j)
208 qh8(1,kflip) = max( qv(i,k,j)/(1. + qv(i,k,j)), 1e-30 ) !values of 0 cause a crash in entropy
209 if( present(qc) ) then
210 qc8(1,kflip) = qc(i,k,j)/(1. + qv(i,k,j)) !Convert to moist mix ratio
211 else
212 qc8(1,kflip) = 0.
213 end if
214 if( present(qi) ) then
215 qi8(1,kflip) = qi(i,k,j)/(1. + qv(i,k,j)) !Used in convtran, ditto for conversion
216 else
217 qi8(1,kflip) = 0.
218 end if
219 pdeldry8(1,kflip)= pdel8(1,kflip)*(1._r8 - qh8(1,kflip))
220 t8(1,kflip) = t_phy(i,k,j)
221 s8(1,kflip) = cpair*t8(1,kflip) + gravit*(z(i,k,j)-ht(i,j))
222 u8(1,kflip) = u_phy(i,k,j)
223 v8(1,kflip) = v_phy(i,k,j)
224 zm8(1,kflip) = dz8w(i,k,j)
225 end do
227 !TKE in CAM is on the interfaces but in WRF it is on the layer
228 !middle. We will interpolate the TKE from the to the interfaces
229 !and then just use the lowest TKE for the surface and the highest
230 !TKE at the top.
231 tke8(1,kte+1) = tke_pbl(i,1,j) !surface
232 tke8(1,1) = tke_pbl(i,kte,j) !model top interface
233 do k = kts,kte-1
234 kflip = kte-k+1
235 tke8(1,kflip) = 0.5*(tke_pbl(i,k,j) + tke_pbl(i,k+1,j))
236 end do
238 !Flip the tracer array -
239 !shift tracer dimension down one to remove "blank" index and
240 !convert to wet instead of dry mixing ratios.
241 do k = kts,kte
242 kflip = kte-k+1
243 !!$ do m = 1,ncnst
244 !!$ moist8(1,kflip,m) = moist(i,k,j,m+1)/(1. + qv(i,k,j))
245 !!$ end do
247 !~For now, send replicate part of the tracer array and send zeros for the
248 ! rest since CAM treats condensate diagnostically compared to WRF that is
249 ! prognostic. This avoids issues with hard-wired assumptions in the UW
250 ! ShCu scheme. This should be looked at again when more time is available.
251 ! Set to zero for most then overwrite as needed...
252 moist8(1,kflip,1:ncnst) = 0.
254 moist8(1,kflip,1) = qv(i,k,j)/(1. + qv(i,k,j))
256 call cnst_get_ind( 'CLDLIQ', m )
257 moist8(1,kflip,m) = qc(i,k,j)/(1. + qv(i,k,j))
259 call cnst_get_ind( 'CLDICE', m )
260 moist8(1,kflip,m) = qi(i,k,j)/(1. + qv(i,k,j))
262 call cnst_get_ind( 'NUMLIQ', m )
263 moist8(1,kflip,m) = 0.
265 call cnst_get_ind( 'NUMICE', m )
266 moist8(1,kflip,m) = 0.
267 end do
269 !Some remapping to get arrays to pass into the routine
270 pblh(1) = pblh_in(i,j)
271 cush(1) = cush_inout(i,j)
272 !
273 ! Main guts of the routine...
274 ! This is a bit inefficient because we are flippling the arrays and they
275 ! will then get flipped back again by compute_uwshcu_inv. We are doing
276 ! this to preserve the CAM code as much as possible for maintenance.
277 !
278 call compute_uwshcu_inv( &
279 pcols, pver, ncol, ncnst, ztodt, &
280 pint8, zi8, pmid8, zm8, pdel8, &
281 u8, v8, qh8, qc8, qi8, &
282 t8, s8, moist8, &
283 tke8, cld8, cldold8, pblh, cush, &
284 cmfmc2, slflx, qtflx, &
285 qhtnd, qctnd, qitnd, &
286 stnd, utnd, vtnd, tnd_tracer, &
287 rprdsh, cmfdqs, precc, snow, &
288 evapcsh, shfrc, iccmr_UW, icwmr_UW, &
289 icimr_UW, cbmf, qc2, rliq2, &
290 cnt2, cnb2, fqsatd, lchnk, pdeldry8 )
291 !
292 ! Map output into WRF-dimensioned arrays...
293 !
294 cush_inout(i,j) = cush(1)
296 do k = kts,kte
297 kflip = kte-k+1
299 !Add shallow reserved cloud condensate to deep reserved cloud condensate
300 ! dlf (kg/kg/s, qc in CAM), rliq done below
301 dlf(i,k,j) = dlf(i,k,j) + qc2(1,kflip)
303 evapcsh_inout(i,k,j)= evapcsh(1,kflip)
304 icwmrsh(i,k,j) = icwmr_uw(1,kflip)
306 rprdsh(1,kflip) = rprdsh(1,kflip) + cmfdqs(1,kflip)
307 rprdsh_inout(i,k,j) = rprdsh(1,kflip)
308 !Not doing rprdtot for now since not yet used by other CAM routines in WRF
310 !Tendencies of winds, potential temperature, and moisture
311 !fields treated specifically by UW scheme
312 rushten(i,k,j) = utnd(1,kflip)
313 rvshten(i,k,j) = vtnd(1,kflip)
314 rthshten(i,k,j) = stnd(1,kflip)/cpair/pi_phy(i,k,j)
315 rqvshten(i,k,j) = qhtnd(1,kflip)/(1. - qv(i,k,j))
316 if( p_qc >= param_first_scalar ) &
317 rqcshten(i,k,j) = qctnd(1,kflip)/(1. - qv(i,k,j))
318 if( p_qi >= param_first_scalar ) &
319 rqishten(i,k,j) = qitnd(1,kflip)/(1. - qv(i,k,j))
321 !~Turn off tendencies for most condensates since CAM treats them diagnostically.
322 ! !Tendencies of tracers except qv,qc,qi
323 !!~need to make sure qg tendency is propagated through to application
324 ! do m = 4,ncnst
325 ! mp1 = m+1 !shift to p_ value for the tracer
326 ! if( mp1==p_qr ) then
327 ! rqrshten(i,k,j) = tnd_tracer(1,kflip,m)/(1. - qv(i,k,j))
328 ! else if( mp1==p_qs ) then
329 ! rqsshten(i,k,j) = tnd_tracer(1,kflip,m)/(1. - qv(i,k,j))
330 ! else if( mp1==p_qg ) then
331 ! rqgshten(i,k,j) = tnd_tracer(1,kflip,m)/(1. - qv(i,k,j))
332 ! else
333 ! write(msg,'(a,i3)') "WARNING: UW shallow Cu cannot handle tracer ",m
334 ! call wrf_debug(100, msg)
335 ! end if
336 ! end do
338 !Combine shallow and deep cumulus updraft mass flux
339 cmfmc2_inout(i,k,j) = cmfmc2(1,kflip)
340 cmfmc(i,k,j) = cmfmc(i,k,j) + cmfmc2(1,kflip)
342 end do !k-loop to kte
344 do k = kts,kte+1
345 kflip = kte-k+2
347 !Convective fluxes of 'sl' and 'qt' in energy unit
348 cmfsl(i,k,j) = slflx(1,kflip)
349 cmflq(i,k,j) = qtflx(1,kflip)*latvap
350 end do !k-loop to kte+1
352 !Calculate fractional occurance of shallow convection
353 !~Not doing this since it would require adding time averaging ability across output times
355 !Rain rate for shallow convection
356 ! rainsh(i,j) = precc(1)*1e3
357 pratesh(i,j) = precc(1)*1e3/dt !~this will need changing for adaptive time steps and cudt
359 !Get indices of convection top and bottom based on deep+shallow
360 !Note: cnt2 and cnb2 have indices decreasing with height, but
361 ! cutop and cubot have indicies increasing with height
362 kflip = kte - cutop(i,j) + 1
363 cnt = kflip
364 if( cnt2(1) < kflip ) cnt = cnt2(1)
365 cutop(i,j) = kte - cnt + 1
367 kflip = kte - cubot(i,j) + 1
368 cnb = kflip
369 if( cnb2(1) > kflip ) cnb = cnb2(1)
370 cubot(i,j) = kte - cnb + 1
372 !Add shallow reserved cloud condensate to deep reserved cloud condensate
373 !dlf done above, rliq (m/s)
374 rliq2_inout(i,j) = rliq2(1)
375 rliq(i,j) = rliq(i,j) + rliq2(1)
377 end do
378 end do ij_loops
379 END SUBROUTINE camuwshcu_driver
381 END MODULE module_cu_camuwshcu_driver