source: lmdz_wrf/trunk/WRFV3/phys/module_cu_camuwshcu_driver.F @ 354

Last change on this file since 354 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 18.8 KB
Line 
1MODULE module_cu_camuwshcu_driver
2  USE shr_kind_mod,    only: r8 => shr_kind_r8
3
4! Roughly based on convect_shallow_tend in convect_shallow.F90 from CAM
5! but tailored for the UW shallow cumulus scheme.
6
7  IMPLICIT NONE
8
9  PRIVATE                  !Default to private
10  PUBLIC :: &              !Public entities
11       camuwshcu_driver
12
13CONTAINS
14
15!------------------------------------------------------------------------
16SUBROUTINE 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
47
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
53
54  REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ), INTENT(IN) :: &
55                              moist    !moist tracer array
56
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)
71
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)
75
76  REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: &
77                            pblh_in, & !height of PBL (m)
78                            ht         !Terrain height (m)
79
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)
99
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?
110
111  REAL, INTENT(IN) :: &
112                                 dt    !time step (s)
113
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
119
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?)
126                                                           
127
128
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)
158
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)
168
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
191
192        !Flip variables on the layer interfaces
193        do k = kts,kte+1
194           kflip = kte-k+2
195
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
199
200        !Flip variables on the layer middles
201        do k = kts,kte
202           kflip = kte-k+1
203
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
226
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
237
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
246
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.
253
254           moist8(1,kflip,1) = qv(i,k,j)/(1. + qv(i,k,j))
255
256           call cnst_get_ind( 'CLDLIQ', m )
257           moist8(1,kflip,m) = qc(i,k,j)/(1. + qv(i,k,j))
258
259           call cnst_get_ind( 'CLDICE', m )
260           moist8(1,kflip,m) = qi(i,k,j)/(1. + qv(i,k,j))
261
262           call cnst_get_ind( 'NUMLIQ', m )
263           moist8(1,kflip,m) = 0.
264
265           call cnst_get_ind( 'NUMICE', m )
266           moist8(1,kflip,m) = 0.
267        end do
268
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)
295
296        do k = kts,kte
297           kflip = kte-k+1
298
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)
302
303           evapcsh_inout(i,k,j)= evapcsh(1,kflip)
304           icwmrsh(i,k,j)      = icwmr_uw(1,kflip)
305
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
309
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))
320
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
337
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)
341
342        end do !k-loop to kte
343
344        do k = kts,kte+1
345           kflip = kte-k+2
346
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
351
352        !Calculate fractional occurance of shallow convection
353        !~Not doing this since it would require adding time averaging ability across output times
354
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
358
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
366
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
371
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)
376
377     end do
378  end do ij_loops
379END SUBROUTINE camuwshcu_driver
380
381END MODULE module_cu_camuwshcu_driver
Note: See TracBrowser for help on using the repository browser.