| blob | 1e4aa790558a9f24cc7c1e619a24883cc09bfadc |
1 !WRF:ODEL_LAYER:PHYSICS
2 !
4 MODULE module_ctrans_grell
5 USE module_cu_gd
6 USE module_dep_simple
7 USE module_state_description, only:p_co,p_qv,p_so2,p_hno3,p_hno4,p_n2o5,p_nh3,p_h2o2, &
8 p_o3,p_ora1,p_op1,p_paa,p_sulf,p_so4aj,p_nh4aj,p_no3aj, &
9 p_bc1,p_bc2,p_oc1,p_oc2,p_seas_1,p_seas_2, &
10 p_seas_3,p_seas_4,p_dms, &
11 p_facd,p_mepx,p_pacd
12 USE module_state_description, only:p_cvasoaX,p_cvasoa1,p_cvasoa2,p_cvasoa3,p_cvasoa4,&
13 p_cvbsoaX,p_cvbsoa1,p_cvbsoa2,p_cvbsoa3,p_cvbsoa4
15 USE module_state_description, ONLY: mozart_mosaic_4bin_kpp, &
16 mozart_mosaic_4bin_aq_kpp, &
17 mozcart_kpp, t1_mozcart_kpp, &
18 p_hcho, p_c3h6ooh, p_onit, p_mvk, p_macr, &
19 p_etooh, p_prooh, p_acetp, p_mgly, p_mvkooh, &
20 p_onitr, p_isooh, p_ch3oh, p_c2h5oh, &
21 p_glyald, p_hydrald, p_ald, p_isopn, &
22 p_alkooh, p_mekooh, p_tolooh, p_terpooh, &
23 p_xooh, p_ch3cooh, p_hcooh, p_ch3ooh, &
24 p_so4_a01,p_no3_a01,p_smpa_a01,p_smpbb_a01,&
25 p_glysoa_r1_a01,p_glysoa_r2_a01,&
26 p_glysoa_sfc_a01,p_glysoa_nh4_a01,&
27 p_glysoa_oh_a01,&
28 p_asoaX_a01,p_asoa1_a01,p_asoa2_a01,p_asoa3_a01,p_asoa4_a01,&
29 p_bsoaX_a01,p_bsoa1_a01,p_bsoa2_a01,p_bsoa3_a01,p_bsoa4_a01,&
30 p_biog1_c_a01,p_biog1_o_a01,&
31 p_cl_a01,p_co3_a01,p_nh4_a01,p_na_a01,&
32 p_ca_a01,p_oin_a01,p_oc_a01,p_bc_a01,&
33 p_so4_a02,p_no3_a02,p_smpa_a02,p_smpbb_a02,&
34 p_glysoa_r1_a02,p_glysoa_r2_a02,&
35 p_glysoa_sfc_a02,p_glysoa_nh4_a02,&
36 p_glysoa_oh_a02,&
37 p_asoaX_a02,p_asoa1_a02,p_asoa2_a02,p_asoa3_a02,p_asoa4_a02,&
38 p_bsoaX_a02,p_bsoa1_a02,p_bsoa2_a02,p_bsoa3_a02,p_bsoa4_a02,&
39 p_biog1_c_a02,p_biog1_o_a02,&
40 p_cl_a02,p_co3_a02,p_nh4_a02,p_na_a02,&
41 p_ca_a02,p_oin_a02,p_oc_a02,p_bc_a02,&
42 p_so4_a03,p_no3_a03,p_smpa_a03,p_smpbb_a03,&
43 p_glysoa_r1_a03,p_glysoa_r2_a03,&
44 p_glysoa_sfc_a03,p_glysoa_nh4_a03,&
45 p_glysoa_oh_a03,&
46 p_asoaX_a03,p_asoa1_a03,p_asoa2_a03,p_asoa3_a03,p_asoa4_a03,&
47 p_bsoaX_a03,p_bsoa1_a03,p_bsoa2_a03,p_bsoa3_a03,p_bsoa4_a03,&
48 p_biog1_c_a03,p_biog1_o_a03,&
49 p_cl_a03,p_co3_a03,p_nh4_a03,p_na_a03,&
50 p_ca_a03,p_oin_a03,p_oc_a03,p_bc_a03,&
51 p_so4_a04,p_no3_a04,p_smpa_a04,p_smpbb_a04,&
52 p_glysoa_r1_a04,p_glysoa_r2_a04,&
53 p_glysoa_sfc_a04,p_glysoa_nh4_a04,&
54 p_glysoa_oh_a04,&
55 p_asoaX_a04,p_asoa1_a04,p_asoa2_a04,p_asoa3_a04,p_asoa4_a04,&
56 p_bsoaX_a04,p_bsoa1_a04,p_bsoa2_a04,p_bsoa3_a04,p_bsoa4_a04,&
57 p_biog1_c_a04,p_biog1_o_a04,&
58 p_cl_a04,p_co3_a04,p_nh4_a04,p_na_a04,&
59 p_ca_a04,p_oin_a04,p_oc_a04,p_bc_a04
61 IMPLICIT NONE
63 ! REAL, PARAMETER :: qcldwtr_cutoff = 1.0e-6 ! kg/kg
64 REAL, PARAMETER :: qcldwtr_cutoff = 1.0e-6 ! kg/m3
65 !
67 REAL, PARAMETER :: mwdry = 28.966 ! Molecular mass of dry air (g/mol)
68 REAL, PARAMETER :: mwso4 = 96.00 ! Molecular mass of SO4-- (g/mol)
69 REAL, PARAMETER :: mwno3 = 62.0 ! Molecular mass of NO3- (g/mol)
70 REAL, PARAMETER :: mwnh4 = 18.0985 ! Molecular mass of NH4+ (g/mol)
72 REAL, PARAMETER :: mwoa = 250.0 ! Molecular mass of OA (g/mol)
74 INTEGER, allocatable :: HLC_ndx(:)
76 CONTAINS
78 !-------------------------------------------------------------
79 SUBROUTINE GRELLDRVCT(DT,itimestep,DX, &
80 rho_phy,RAINCV,chem, &
81 U,V,t_phy,moist,dz8w,p_phy, &
82 XLV,CP,G,r_v,z,cu_co_ten, &
83 wd_no3,wd_so4, &
84 wd_nh4,wd_oa, &
85 wd_so2, wd_sulf, wd_hno3, wd_nh3, &
86 wd_cvasoa, wd_cvbsoa, wd_asoa, wd_bsoa, &
87 k22_shallow,kbcon_shallow,ktop_shallow,xmb_shallow, &
88 ishallow,num_moist,numgas,num_chem,chemopt,scalaropt, &
89 conv_tr_wetscav,conv_tr_aqchem, &
90 ids,ide, jds,jde, kds,kde, &
91 ims,ime, jms,jme, kms,kme, &
92 its,ite, jts,jte, kts,kte )
93 ! USE module_configure
94 ! USE module_state_description
95 !-------------------------------------------------------------
96 IMPLICIT NONE
97 !-------------------------------------------------------------
98 INTEGER, INTENT(IN ) :: &
99 numgas,chemopt,scalaropt, &
100 ids,ide, jds,jde, kds,kde, &
101 ims,ime, jms,jme, kms,kme, &
102 its,ite, jts,jte, kts,kte, &
103 ishallow,num_chem,num_moist, &
104 conv_tr_wetscav,conv_tr_aqchem
106 INTEGER, INTENT(IN ) :: ITIMESTEP
108 REAL, INTENT(IN ) :: XLV, R_v
109 REAL, INTENT(IN ) :: CP,G
110 REAL, DIMENSION( ims:ime , kms:kme , jms:jme,num_moist ) , &
111 INTENT(IN ) :: moist
112 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
113 INTENT(IN ) :: &
114 U, &
115 V, &
116 t_phy, &
117 z, &
118 p_phy, &
119 dz8w, &
120 rho_phy
121 !
122 ! on output for control only, purely diagnostic
123 !
124 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
125 INTENT(INOUT ) :: &
126 cu_co_ten
128 !
129 REAL, INTENT(IN ) :: DT, DX
130 !
131 REAL, DIMENSION( ims:ime , kms:kme , jms:jme, num_chem ), &
132 INTENT(INOUT) :: &
133 chem
135 REAL, DIMENSION( ims:ime , jms:jme ), &
136 INTENT(IN) :: RAINCV, xmb_shallow
137 INTEGER, DIMENSION( ims:ime , jms:jme ), &
138 INTENT(IN) :: k22_shallow,kbcon_shallow,ktop_shallow
139 !
140 ! Accumulated wet deposition
141 !
142 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: wd_no3,wd_so4
143 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: wd_nh4,wd_oa, &
144 wd_so2, wd_sulf, wd_hno3, wd_nh3
145 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: &
146 wd_cvasoa,wd_cvbsoa,wd_asoa,wd_bsoa
147 ! LOCAL VARS
148 real, dimension (its:ite,kts:kte) :: &
149 OUTT,OUTQ,OUTQC
150 real, dimension (its:ite) :: &
151 pret, ter11
152 ! Auxiliary wet deposition variables
153 REAL, DIMENSION (its:ite,num_chem) :: wetdep_1d
154 REAL, DIMENSION (ims:ime,jms:jme,num_chem) :: wetdep_2d
155 ! Wet deposition over the current time step
156 REAL, DIMENSION( ims:ime , jms:jme ) :: wdi_no3,wdi_so4
157 REAL, DIMENSION( ims:ime , jms:jme ) :: wdi_nh4,wdi_oa, &
158 wdi_so2, wdi_sulf, wdi_hno3, wdi_nh3
159 REAL, DIMENSION( ims:ime , jms:jme ) :: wdi_cvasoa, wdi_cvbsoa, &
160 wdi_asoa, wdi_bsoa
162 ! basic environmental input includes moisture convergence (mconv)
163 ! omega (omeg), windspeed (us,vs), and a flag (aaeq) to turn off
164 ! convection for this call only and at that particular gridpoint
165 !
166 real, dimension (its:ite,kts:kte) :: &
167 T,TN,q,qo,PO,P,US,VS,hstary
168 real, dimension (its:ite,kts:kte,num_chem) :: &
169 tracer,tracert,tracert3
170 real, dimension (its:ite) :: &
171 Z1,PSUR,AAEQ,xmb3
172 integer, dimension (its:ite) :: &
173 ktop,k23,kbcon3,ktop3
175 INTEGER :: ishallow_gf = 0
176 INTEGER :: nv,i,j,k,ICLDCK,ipr,jpr,npr
177 REAL :: tcrit,dp,dq,epsilc
178 INTEGER :: itf,jtf,ktf,iopt
179 epsilc=1.e-30
180 wetdep_1d(:,:) = 0.0
181 wetdep_2d(:,:,:) = 0.0
182 ! return
183 ! ipr=111
184 ! jpr=40
185 ! if(itimestep.lt.34.or.itimestep.gt.36)ipr=0
186 ! if(itimestep.lt.34.or.itimestep.gt.36)jpr=0
187 ! ipr=61
188 ! jpr=60
189 ipr=0
190 jpr=0
191 npr=1
192 if(p_co.gt.1)npr=p_co
193 tcrit=258.
194 iopt=0
195 itf=MIN(ite,ide-1)
196 ktf=MIN(kte,kde-1)
197 jtf=MIN(jte,jde-1)
198 !
199 !
200 123 continue
201 DO 100 J = jts,jtf
202 if(j.eq.jpr)print *,'dt = ',dt
203 DO I=ITS,ITF
204 ktop(i)=0
205 PSUR(I)=p_phy(I,kts,J)*.01
206 TER11(I)=z(i,kts,j)
207 aaeq(i)=0.
209 !
210 ! rainrate is input for this transport/wet-deposition routine
211 !
212 pret(i)=raincv(i,j)/dt
213 if(pret(i).le.0.)aaeq(i)=20.
214 ENDDO
215 if(Ishallow.eq.1)then
216 DO I=ITS,ITF
217 xmb3(i)=xmb_shallow(i,j)
218 ktop3(i)=ktop_shallow(i,j)
219 k23(i)=k22_shallow(i,j)
220 kbcon3(i)=kbcon_shallow(i,j)
221 ENDDO
222 else
223 DO I=ITS,ITF
224 xmb3(i)=0.
225 ktop3(i)=0
226 k23(i)=0
227 kbcon3(i)=0
228 ENDDO
229 endif
230 DO K=kts,ktf
231 DO I=ITS,ITF
232 po(i,k)=p_phy(i,k,j)*.01
233 P(I,K)=PO(i,k)
234 US(I,K) =u(i,k,j)
235 VS(I,K) =v(i,k,j)
236 T(I,K)=t_phy(i,k,j)
237 q(I,K)=moist(i,k,j,p_qv)
238 IF(Q(I,K).LT.1.E-08)Q(I,K)=1.E-08
239 ENDDO
240 ENDDO
241 do nv=2,num_chem
242 DO K=kts,ktf
243 DO I=ITS,ITF
244 tracer(i,k,nv)=max(epsilc,chem(i,k,j,nv))
245 tracert(i,k,nv)=0.
246 tracert3(i,k,nv)=0.
247 ENDDO
248 ENDDO
249 ENDDO
250 DO K=kts,ktf
251 DO I=ITS,ITF
252 cu_co_ten(i,k,j)=0.
253 ! hstary(i,k)=hstar4(nv)*exp(dhr(nv)*(1./t(i,k)-1./298.))
254 if(i.eq.ipr.and.j.eq.jpr)then
255 print *,k,pret(i),tracer(i,k,npr),p(i,k),z(i,k,j)
256 endif
257 ENDDO
258 ENDDO
259 ! ENDDO
260 !
261 !---- CALL NON_RESOLVED CONVECTIVE TRANSPORT
262 !
263 CALL CUP_ct(ktop,k23,kbcon3,ktop3,xmb3, &
264 tracer,j,AAEQ,T,Q,TER11,PRET,P,tracert,tracert3, &
265 hstary,DT,PSUR,US,VS,tcrit,wetdep_1d, &
266 xlv,r_v,cp,g,ipr,jpr,npr,num_chem,chemopt,scalaropt,&
267 conv_tr_wetscav,conv_tr_aqchem, &
268 ishallow,numgas,ids,ide, jds,jde, kds,kde, &
269 ims,ime, jms,jme, kms,kme, &
270 its,ite, jts,jte, kts,kte )
272 do nv=2,num_chem
273 DO I=its,itf
274 if(pret(i).le.0.)then
275 DO K=kts,ktf
276 tracert(i,k,nv)=0.
277 ENDDO
278 endif
279 enddo
280 enddo
281 if(ishallow.eq.1)then
282 CALL neg_check_ct('shallow',pret,ktop3,epsilc,dt,tracer,tracert3,iopt,num_chem, &
283 its,ite,kts,kte,itf,ktf,ipr,jpr,npr,j)
284 do nv=2,num_chem
285 DO I=its,itf
286 DO K=kts,ktf
287 chem(i,k,j,nv)=max(epsilc,chem(i,k,j,nv)+tracert3(i,k,nv)*dt)
288 enddo
289 ENDDO
290 ENDDO
291 endif
292 !
293 ! now deep convection
294 !
295 CALL neg_check_ct('deep',pret,ktop,epsilc,dt,tracer,tracert,iopt,num_chem, &
296 its,ite,kts,kte,itf,ktf,ipr,jpr,npr,j)
297 do nv=2,num_chem
298 DO I=its,itf
299 if(pret(i).gt.0.)then
300 DO K=kts,ktf
301 chem(i,k,j,nv)=max(epsilc,chem(i,k,j,nv)+tracert(i,k,nv)*dt)
302 if(nv.eq.npr)then
303 cu_co_ten(i,k,j)=tracert(i,k,npr)*dt
304 ! if(i.eq.ipr.and.j.eq.jpr)print *,k,chem(i,k,j,nv),cu_co_ten(i,k,j)
305 endif
306 ENDDO
307 else
308 DO K=kts,ktf
309 tracert(i,k,nv)=0.
310 if(nv.eq.npr)cu_co_ten(i,k,j)=0.
311 enddo
312 endif
313 ENDDO
314 ENDDO
316 wetdep_2d(its:ite,j,:) = wetdep_1d(its:ite,:)
318 100 continue
320 if(chemopt > 0) then ! skip for tracers (chemopt=0)
321 ! Calculate the wet deposition over the time step:
323 wdi_no3(:,:) = 0.0
324 wdi_so4(:,:) = 0.0
325 wdi_nh4(:,:) = 0.0
326 wdi_oa(:,:) = 0.0
327 wdi_so2(:,:) = 0.0
328 wdi_sulf(:,:) = 0.0
329 wdi_hno3(:,:) = 0.0
330 wdi_nh3(:,:) = 0.0
332 wdi_cvasoa(:,:) = 0.0
333 wdi_cvbsoa(:,:) = 0.0
334 wdi_asoa(:,:) = 0.0
335 wdi_bsoa(:,:) = 0.0
337 ! We use the indices of the chem array that point to aerosol outside of cloud water,
338 ! because that's what the cumulus scheme operates with.
340 if (chemopt == mozart_mosaic_4bin_kpp .OR. &
341 chemopt == mozart_mosaic_4bin_aq_kpp) then
343 wdi_no3(its:ite,jts:jte) = wdi_no3(its:ite,jts:jte) + &
344 (wetdep_2d(its:ite,jts:jte,p_no3_a01) + &
345 wetdep_2d(its:ite,jts:jte,p_no3_a02) + &
346 wetdep_2d(its:ite,jts:jte,p_no3_a03) + &
347 wetdep_2d(its:ite,jts:jte,p_no3_a04))*dt*0.001/mwno3 ! mmol/m2
349 wdi_so4(its:ite,jts:jte) = wdi_so4(its:ite,jts:jte) + &
350 (wetdep_2d(its:ite,jts:jte,p_so4_a01) + &
351 wetdep_2d(its:ite,jts:jte,p_so4_a02) + &
352 wetdep_2d(its:ite,jts:jte,p_so4_a03) + &
353 wetdep_2d(its:ite,jts:jte,p_so4_a04))*dt*0.001/mwso4 ! mmol/m2
355 wdi_nh4(its:ite,jts:jte) = wdi_nh4(its:ite,jts:jte) + &
356 (wetdep_2d(its:ite,jts:jte,p_nh4_a01) + &
357 wetdep_2d(its:ite,jts:jte,p_nh4_a02) + &
358 wetdep_2d(its:ite,jts:jte,p_nh4_a03) + &
359 wetdep_2d(its:ite,jts:jte,p_nh4_a04))*dt*0.001/mwnh4 ! mmol/m2
361 if (chemopt == mozart_mosaic_4bin_kpp) then
363 wdi_oa(its:ite,jts:jte) = wdi_oa(its:ite,jts:jte) + &
364 (wetdep_2d(its:ite,jts:jte,p_oc_a01) + &
365 wetdep_2d(its:ite,jts:jte,p_oc_a02) + &
366 wetdep_2d(its:ite,jts:jte,p_oc_a03) + &
367 wetdep_2d(its:ite,jts:jte,p_oc_a04) + &
368 wetdep_2d(its:ite,jts:jte,p_smpa_a01) + &
369 wetdep_2d(its:ite,jts:jte,p_smpa_a02) + &
370 wetdep_2d(its:ite,jts:jte,p_smpa_a03) + &
371 wetdep_2d(its:ite,jts:jte,p_smpa_a04) + &
372 wetdep_2d(its:ite,jts:jte,p_smpbb_a01) + &
373 wetdep_2d(its:ite,jts:jte,p_smpbb_a02) + &
374 wetdep_2d(its:ite,jts:jte,p_smpbb_a03) + &
375 wetdep_2d(its:ite,jts:jte,p_smpbb_a04) + &
376 wetdep_2d(its:ite,jts:jte,p_biog1_c_a01) + &
377 wetdep_2d(its:ite,jts:jte,p_biog1_c_a02) + &
378 wetdep_2d(its:ite,jts:jte,p_biog1_c_a03) + &
379 wetdep_2d(its:ite,jts:jte,p_biog1_c_a04) + &
380 wetdep_2d(its:ite,jts:jte,p_biog1_o_a01) + &
381 wetdep_2d(its:ite,jts:jte,p_biog1_o_a02) + &
382 wetdep_2d(its:ite,jts:jte,p_biog1_o_a03) + &
383 wetdep_2d(its:ite,jts:jte,p_biog1_o_a04) + &
384 wetdep_2d(its:ite,jts:jte,p_glysoa_r1_a01) + &
385 wetdep_2d(its:ite,jts:jte,p_glysoa_r1_a02) + &
386 wetdep_2d(its:ite,jts:jte,p_glysoa_r1_a03) + &
387 wetdep_2d(its:ite,jts:jte,p_glysoa_r1_a04) + &
388 wetdep_2d(its:ite,jts:jte,p_glysoa_r2_a01) + &
389 wetdep_2d(its:ite,jts:jte,p_glysoa_r2_a02) + &
390 wetdep_2d(its:ite,jts:jte,p_glysoa_r2_a03) + &
391 wetdep_2d(its:ite,jts:jte,p_glysoa_r2_a04) + &
392 wetdep_2d(its:ite,jts:jte,p_glysoa_oh_a01) + &
393 wetdep_2d(its:ite,jts:jte,p_glysoa_oh_a02) + &
394 wetdep_2d(its:ite,jts:jte,p_glysoa_oh_a03) + &
395 wetdep_2d(its:ite,jts:jte,p_glysoa_oh_a04) + &
396 wetdep_2d(its:ite,jts:jte,p_glysoa_nh4_a01) + &
397 wetdep_2d(its:ite,jts:jte,p_glysoa_nh4_a02) + &
398 wetdep_2d(its:ite,jts:jte,p_glysoa_nh4_a03) + &
399 wetdep_2d(its:ite,jts:jte,p_glysoa_nh4_a04) + &
400 wetdep_2d(its:ite,jts:jte,p_glysoa_sfc_a01) + &
401 wetdep_2d(its:ite,jts:jte,p_glysoa_sfc_a02) + &
402 wetdep_2d(its:ite,jts:jte,p_glysoa_sfc_a03) + &
403 wetdep_2d(its:ite,jts:jte,p_glysoa_sfc_a04))*dt*0.001/mwoa ! mmol/m2
404 endif
406 if (chemopt == mozart_mosaic_4bin_aq_kpp) then
408 wdi_asoa(its:ite,jts:jte) = wdi_asoa(its:ite,jts:jte) + &
409 (wetdep_2d(its:ite,jts:jte,p_asoaX_a01) + &
410 wetdep_2d(its:ite,jts:jte,p_asoaX_a02) + &
411 wetdep_2d(its:ite,jts:jte,p_asoaX_a03) + &
412 wetdep_2d(its:ite,jts:jte,p_asoaX_a04) + &
413 wetdep_2d(its:ite,jts:jte,p_asoa1_a01) + &
414 wetdep_2d(its:ite,jts:jte,p_asoa1_a02) + &
415 wetdep_2d(its:ite,jts:jte,p_asoa1_a03) + &
416 wetdep_2d(its:ite,jts:jte,p_asoa1_a04) + &
417 wetdep_2d(its:ite,jts:jte,p_asoa2_a01) + &
418 wetdep_2d(its:ite,jts:jte,p_asoa2_a02) + &
419 wetdep_2d(its:ite,jts:jte,p_asoa2_a03) + &
420 wetdep_2d(its:ite,jts:jte,p_asoa2_a04) + &
421 wetdep_2d(its:ite,jts:jte,p_asoa3_a01) + &
422 wetdep_2d(its:ite,jts:jte,p_asoa3_a02) + &
423 wetdep_2d(its:ite,jts:jte,p_asoa3_a03) + &
424 wetdep_2d(its:ite,jts:jte,p_asoa3_a04) + &
425 wetdep_2d(its:ite,jts:jte,p_asoa4_a01) + &
426 wetdep_2d(its:ite,jts:jte,p_asoa4_a02) + &
427 wetdep_2d(its:ite,jts:jte,p_asoa4_a03) + &
428 wetdep_2d(its:ite,jts:jte,p_asoa4_a04))*dt*0.001/150.0 ! mmol/m2
430 wdi_bsoa(its:ite,jts:jte) = wdi_bsoa(its:ite,jts:jte) + &
431 (wetdep_2d(its:ite,jts:jte,p_bsoaX_a01) + &
432 wetdep_2d(its:ite,jts:jte,p_bsoaX_a02) + &
433 wetdep_2d(its:ite,jts:jte,p_bsoaX_a03) + &
434 wetdep_2d(its:ite,jts:jte,p_bsoaX_a04) + &
435 wetdep_2d(its:ite,jts:jte,p_bsoa1_a01) + &
436 wetdep_2d(its:ite,jts:jte,p_bsoa1_a02) + &
437 wetdep_2d(its:ite,jts:jte,p_bsoa1_a03) + &
438 wetdep_2d(its:ite,jts:jte,p_bsoa1_a04) + &
439 wetdep_2d(its:ite,jts:jte,p_bsoa2_a01) + &
440 wetdep_2d(its:ite,jts:jte,p_bsoa2_a02) + &
441 wetdep_2d(its:ite,jts:jte,p_bsoa2_a03) + &
442 wetdep_2d(its:ite,jts:jte,p_bsoa2_a04) + &
443 wetdep_2d(its:ite,jts:jte,p_bsoa3_a01) + &
444 wetdep_2d(its:ite,jts:jte,p_bsoa3_a02) + &
445 wetdep_2d(its:ite,jts:jte,p_bsoa3_a03) + &
446 wetdep_2d(its:ite,jts:jte,p_bsoa3_a04) + &
447 wetdep_2d(its:ite,jts:jte,p_bsoa4_a01) + &
448 wetdep_2d(its:ite,jts:jte,p_bsoa4_a02) + &
449 wetdep_2d(its:ite,jts:jte,p_bsoa4_a03) + &
450 wetdep_2d(its:ite,jts:jte,p_bsoa4_a04))*dt*0.001/180.0 ! mmol/m2
452 wdi_cvasoa(its:ite,jts:jte) = wdi_cvasoa(its:ite,jts:jte) + &
453 (wetdep_2d(its:ite,jts:jte,p_cvasoaX) + &
454 wetdep_2d(its:ite,jts:jte,p_cvasoa1) + &
455 wetdep_2d(its:ite,jts:jte,p_cvasoa2) + &
456 wetdep_2d(its:ite,jts:jte,p_cvasoa3) + &
457 wetdep_2d(its:ite,jts:jte,p_cvasoa4))*dt ! mmol/m2
459 wdi_cvbsoa(its:ite,jts:jte) = wdi_cvbsoa(its:ite,jts:jte) + &
460 (wetdep_2d(its:ite,jts:jte,p_cvbsoaX) + &
461 wetdep_2d(its:ite,jts:jte,p_cvbsoa1) + &
462 wetdep_2d(its:ite,jts:jte,p_cvbsoa2) + &
463 wetdep_2d(its:ite,jts:jte,p_cvbsoa3) + &
464 wetdep_2d(its:ite,jts:jte,p_cvbsoa4))*dt ! mmol/m2
466 wdi_oa(its:ite,jts:jte) = wdi_oa(its:ite,jts:jte) + &
467 (wetdep_2d(its:ite,jts:jte,p_oc_a01) + &
468 wetdep_2d(its:ite,jts:jte,p_oc_a02) + &
469 wetdep_2d(its:ite,jts:jte,p_oc_a03) + &
470 wetdep_2d(its:ite,jts:jte,p_oc_a04) + &
471 wetdep_2d(its:ite,jts:jte,p_asoaX_a01) + &
472 wetdep_2d(its:ite,jts:jte,p_asoaX_a02) + &
473 wetdep_2d(its:ite,jts:jte,p_asoaX_a03) + &
474 wetdep_2d(its:ite,jts:jte,p_asoaX_a04) + &
475 wetdep_2d(its:ite,jts:jte,p_asoa1_a01) + &
476 wetdep_2d(its:ite,jts:jte,p_asoa1_a02) + &
477 wetdep_2d(its:ite,jts:jte,p_asoa1_a03) + &
478 wetdep_2d(its:ite,jts:jte,p_asoa1_a04) + &
479 wetdep_2d(its:ite,jts:jte,p_asoa2_a01) + &
480 wetdep_2d(its:ite,jts:jte,p_asoa2_a02) + &
481 wetdep_2d(its:ite,jts:jte,p_asoa2_a03) + &
482 wetdep_2d(its:ite,jts:jte,p_asoa2_a04) + &
483 wetdep_2d(its:ite,jts:jte,p_asoa3_a01) + &
484 wetdep_2d(its:ite,jts:jte,p_asoa3_a02) + &
485 wetdep_2d(its:ite,jts:jte,p_asoa3_a03) + &
486 wetdep_2d(its:ite,jts:jte,p_asoa3_a04) + &
487 wetdep_2d(its:ite,jts:jte,p_asoa4_a01) + &
488 wetdep_2d(its:ite,jts:jte,p_asoa4_a02) + &
489 wetdep_2d(its:ite,jts:jte,p_asoa4_a03) + &
490 wetdep_2d(its:ite,jts:jte,p_asoa4_a04) + &
491 wetdep_2d(its:ite,jts:jte,p_bsoaX_a01) + &
492 wetdep_2d(its:ite,jts:jte,p_bsoaX_a02) + &
493 wetdep_2d(its:ite,jts:jte,p_bsoaX_a03) + &
494 wetdep_2d(its:ite,jts:jte,p_bsoaX_a04) + &
495 wetdep_2d(its:ite,jts:jte,p_bsoa1_a01) + &
496 wetdep_2d(its:ite,jts:jte,p_bsoa1_a02) + &
497 wetdep_2d(its:ite,jts:jte,p_bsoa1_a03) + &
498 wetdep_2d(its:ite,jts:jte,p_bsoa1_a04) + &
499 wetdep_2d(its:ite,jts:jte,p_bsoa2_a01) + &
500 wetdep_2d(its:ite,jts:jte,p_bsoa2_a02) + &
501 wetdep_2d(its:ite,jts:jte,p_bsoa2_a03) + &
502 wetdep_2d(its:ite,jts:jte,p_bsoa2_a04) + &
503 wetdep_2d(its:ite,jts:jte,p_bsoa3_a01) + &
504 wetdep_2d(its:ite,jts:jte,p_bsoa3_a02) + &
505 wetdep_2d(its:ite,jts:jte,p_bsoa3_a03) + &
506 wetdep_2d(its:ite,jts:jte,p_bsoa3_a04) + &
507 wetdep_2d(its:ite,jts:jte,p_bsoa4_a01) + &
508 wetdep_2d(its:ite,jts:jte,p_bsoa4_a02) + &
509 wetdep_2d(its:ite,jts:jte,p_bsoa4_a03) + &
510 wetdep_2d(its:ite,jts:jte,p_bsoa4_a04) + &
511 wetdep_2d(its:ite,jts:jte,p_glysoa_r1_a01) + &
512 wetdep_2d(its:ite,jts:jte,p_glysoa_r1_a02) + &
513 wetdep_2d(its:ite,jts:jte,p_glysoa_r1_a03) + &
514 wetdep_2d(its:ite,jts:jte,p_glysoa_r1_a04) + &
515 wetdep_2d(its:ite,jts:jte,p_glysoa_r2_a01) + &
516 wetdep_2d(its:ite,jts:jte,p_glysoa_r2_a02) + &
517 wetdep_2d(its:ite,jts:jte,p_glysoa_r2_a03) + &
518 wetdep_2d(its:ite,jts:jte,p_glysoa_r2_a04) + &
519 wetdep_2d(its:ite,jts:jte,p_glysoa_oh_a01) + &
520 wetdep_2d(its:ite,jts:jte,p_glysoa_oh_a02) + &
521 wetdep_2d(its:ite,jts:jte,p_glysoa_oh_a03) + &
522 wetdep_2d(its:ite,jts:jte,p_glysoa_oh_a04) + &
523 wetdep_2d(its:ite,jts:jte,p_glysoa_nh4_a01) + &
524 wetdep_2d(its:ite,jts:jte,p_glysoa_nh4_a02) + &
525 wetdep_2d(its:ite,jts:jte,p_glysoa_nh4_a03) + &
526 wetdep_2d(its:ite,jts:jte,p_glysoa_nh4_a04) + &
527 wetdep_2d(its:ite,jts:jte,p_glysoa_sfc_a01) + &
528 wetdep_2d(its:ite,jts:jte,p_glysoa_sfc_a02) + &
529 wetdep_2d(its:ite,jts:jte,p_glysoa_sfc_a03) + &
530 wetdep_2d(its:ite,jts:jte,p_glysoa_sfc_a04))*dt*0.001/mwoa ! mmol/m2
531 endif
533 else
534 if (p_no3aj .gt.1) wdi_no3(its:ite,jts:jte) = wdi_no3(its:ite,jts:jte) + wetdep_2d(its:ite,jts:jte,p_no3aj)*dt*0.001/mwno3 ! mmol/m2
535 if (p_so4aj .gt.1) wdi_so4(its:ite,jts:jte) = wdi_so4(its:ite,jts:jte) + wetdep_2d(its:ite,jts:jte,p_so4aj)*dt*0.001/mwso4 ! mmol/m2
536 endif
538 if (p_hno3 .gt.1) wdi_hno3(its:ite,jts:jte) = wdi_hno3(its:ite,jts:jte) + wetdep_2d(its:ite,jts:jte,p_hno3)*dt ! mmol/m2
539 if (p_hno4 .gt.1) wdi_hno3(its:ite,jts:jte) = wdi_hno3(its:ite,jts:jte) + wetdep_2d(its:ite,jts:jte,p_hno4)*dt ! mmol/m2
541 if (p_sulf .gt.1) wdi_sulf(its:ite,jts:jte) = wdi_sulf(its:ite,jts:jte) + wetdep_2d(its:ite,jts:jte,p_sulf)*dt ! mmol/m2
542 if (p_so2 .gt.1) wdi_so2(its:ite,jts:jte) = wdi_so2(its:ite,jts:jte) + wetdep_2d(its:ite,jts:jte,p_so2)*dt ! mmol/m2
544 if (p_nh3 .gt.1) wdi_nh3(its:ite,jts:jte) = wdi_nh3(its:ite,jts:jte) + wetdep_2d(its:ite,jts:jte,p_nh3)*dt ! mmol/m2
546 ! Update the accumulated wet deposition:
548 wd_no3(its:ite,jts:jte) = wd_no3(its:ite,jts:jte) + wdi_no3(its:ite,jts:jte) ! mmol/m2
549 wd_so4(its:ite,jts:jte) = wd_so4(its:ite,jts:jte) + wdi_so4(its:ite,jts:jte) ! mmol/m2
550 wd_nh4(its:ite,jts:jte) = wd_nh4(its:ite,jts:jte) + wdi_nh4(its:ite,jts:jte) ! mmol/m2
551 wd_oa (its:ite,jts:jte) = wd_oa (its:ite,jts:jte) + wdi_oa (its:ite,jts:jte) ! mmol/m2
552 wd_so2 (its:ite,jts:jte) = wd_so2 (its:ite,jts:jte) + wdi_so2 (its:ite,jts:jte) ! mmol/m2
553 wd_sulf (its:ite,jts:jte) = wd_sulf (its:ite,jts:jte) + wdi_sulf (its:ite,jts:jte) ! mmol/m2
554 wd_hno3 (its:ite,jts:jte) = wd_hno3 (its:ite,jts:jte) + wdi_hno3 (its:ite,jts:jte) ! mmol/m2
555 wd_nh3 (its:ite,jts:jte) = wd_nh3 (its:ite,jts:jte) + wdi_nh3 (its:ite,jts:jte) ! mmol/m2
557 wd_asoa(its:ite,jts:jte) = wd_asoa(its:ite,jts:jte) + wdi_asoa(its:ite,jts:jte) ! mmol/m2
558 wd_bsoa(its:ite,jts:jte) = wd_bsoa(its:ite,jts:jte) + wdi_bsoa(its:ite,jts:jte) ! mmol/m2
559 wd_cvasoa(its:ite,jts:jte) = wd_cvasoa(its:ite,jts:jte) + wdi_cvasoa(its:ite,jts:jte) ! mmol/m2
560 wd_cvbsoa(its:ite,jts:jte) = wd_cvbsoa(its:ite,jts:jte) + wdi_cvbsoa(its:ite,jts:jte) ! mmol/m2
562 endif
564 END SUBROUTINE GRELLDRVCT
567 SUBROUTINE CUP_ct(ktop,k23,kbcon3,ktop3,xmb3,tracer,J,AAEQ,T,Q,Z1, &
568 PRE,P,tracert,tracert3,hstary,DTIME,PSUR,US,VS,TCRIT, &
569 wetdep,xl,rv,cp,g,ipr,jpr,npr,num_chem,chemopt,scalaropt,&
570 conv_tr_wetscav,conv_tr_aqchem, &
571 ishallow,numgas,ids,ide, jds,jde, kds,kde, &
572 ims,ime, jms,jme, kms,kme, &
573 its,ite, jts,jte, kts,kte )
575 IMPLICIT NONE
577 integer &
578 ,intent (in ) :: &
579 num_chem,ids,ide, jds,jde, kds,kde, &
580 ims,ime, jms,jme, kms,kme,scalaropt,conv_tr_wetscav, &
581 conv_tr_aqchem,ishallow, &
582 its,ite, jts,jte, kts,kte,ipr,jpr,npr,chemopt,numgas
583 integer, intent (in ) :: &
584 j
585 !
586 !
587 !
588 !tracert = output temp tendency (per s)
589 ! pre = input precip
590 real, dimension (its:ite,kts:kte,num_chem) &
591 ,intent (inout ) :: &
592 tracert,tracer,tracert3
593 real, dimension (its:ite) &
594 ,intent (inout ) :: &
595 pre
596 integer, dimension (its:ite) &
597 ,intent (inout ) :: &
598 ktop,k23,kbcon3,ktop3
599 integer, dimension (its:ite) :: &
600 kbcon
601 !
602 ! basic environmental input includes moisture convergence (mconv)
603 ! omega (omeg), windspeed (us,vs), and a flag (aaeq) to turn off
604 ! convection for this call only and at that particular gridpoint
605 !
606 real, dimension (its:ite,kts:kte) &
607 ,intent (in ) :: &
608 T,P,US,VS,HSTARY
609 real, dimension (its:ite,kts:kte) &
610 ,intent (inout) :: &
611 Q
612 real, dimension (its:ite) &
613 ,intent (in ) :: &
614 Z1,PSUR,AAEQ,xmb3
617 real &
618 ,intent (in ) :: &
619 dtime,tcrit,xl,cp,rv,g
622 real, dimension (its:ite,1:3) :: &
623 edtc
624 !
625 !
626 !
627 !***************** the following are your basic environmental
628 ! variables. They carry a "_cup" if they are
629 ! on model cloud levels (staggered). They carry
630 ! an "o"-ending (z becomes zo), if they are the forced
631 ! variables. They are preceded by x (z becomes xz)
632 ! to indicate modification by some typ of cloud
633 !
634 ! z = heights of model levels
635 ! q = environmental mixing ratio
636 ! qes = environmental saturation mixing ratio
637 ! t = environmental temp
638 ! p = environmental pressure
639 ! he = environmental moist static energy
640 ! hes = environmental saturation moist static energy
641 ! z_cup = heights of model cloud levels
642 ! q_cup = environmental q on model cloud levels
643 ! qes_cup = saturation q on model cloud levels
644 ! t_cup = temperature (Kelvin) on model cloud levels
645 ! p_cup = environmental pressure
646 ! he_cup = moist static energy on model cloud levels
647 ! hes_cup = saturation moist static energy on model cloud levels
648 ! gamma_cup = gamma on model cloud levels
649 !
650 !
651 ! hcd = moist static energy in downdraft
652 ! zd normalized downdraft mass flux
653 ! dby = buoancy term
654 ! entr = entrainment rate
655 ! zd = downdraft normalized mass flux
656 ! entr= entrainment rate
657 ! hcd = h in model cloud
658 ! bu = buoancy term
659 ! zd = normalized downdraft mass flux
660 ! gamma_cup = gamma on model cloud levels
661 ! mentr_rate = entrainment rate
662 ! qcd = cloud q (including liquid water) after entrainment
663 ! qrch = saturation q in cloud
664 ! pwd = evaporate at that level
665 ! pwev = total normalized integrated evaoprate (I2)
666 ! entr= entrainment rate
667 ! z1 = terrain elevation
668 ! entr = downdraft entrainment rate
669 ! jmin = downdraft originating level
670 ! kdet = level above ground where downdraft start detraining
671 ! psur = surface pressure
672 ! z1 = terrain elevation
673 ! zd = downdraft normalized mass flux
674 ! zu = updraft normalized mass flux
675 ! mbdt = arbitrary numerical parameter
676 ! dtime = dt over which forcing is applied
677 ! kbcon = LFC of parcel from k22
678 ! k22 = updraft originating level
679 ! dby = buoancy term
680 ! ktop = cloud top (output)
681 ! xmb = total base mass flux
682 ! hc = cloud moist static energy
683 ! hkb = moist static energy at originating level
684 ! mentr_rate = entrainment rate
686 real, dimension (its:ite,kts:kte) :: &
687 he,hes,qes,z,pwdper, &
689 qes_cup,q_cup,he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup, &
691 dby,qc,qrcd,pwd,pw,hcd,qcd,dbyd,hc,qrc,zu,zd,clw_all,zd3, &
692 clw_all3,cd3,pw3,zu3, &
694 ! cd = detrainment function for updraft
695 ! cdd = detrainment function for downdraft
697 cd,cdd,cdd3,scr1,DELLAH,DELLAQ,DELLAT,DELLAQC
699 ! edt = epsilon
700 ! edt = epsilon
701 real, dimension (its:ite) :: &
702 edt,HKB,QKB,edt3, &
703 XMB,PWAV,PWEV,BU,cap_max,cap_max_increment
704 real, dimension (its:ite,kts:kte,num_chem) :: &
705 tr3_c,tr3_up,tr3_pw
706 real, dimension (its:ite,num_chem) :: &
707 trkb3
708 real, dimension (its:ite,kts:kte,num_chem) :: &
709 tr_c,tr_up,tr_dd,tr_dd3,tre_cup,tr_pw,tr_pwd
710 real, dimension (its:ite,num_chem) :: &
711 trkb,wetdep
712 integer, dimension (its:ite) :: &
713 kzdown,KDET,K22,KB,JMIN,kstabi,kstabm, & !-lxz
714 ierr,KBMAX,ierr5
716 integer :: &
717 ki,I,K,KK
718 real :: &
719 day,dz,mbdt,entr_rate,radius,entrd_rate,mentr_rate,mentrd_rate, &
720 zcutdown,edtmax,edtmin,depth_min,zkbmax,z_detr,zktop, &
721 dh,cap_maxs,entr_rate3,mentr_rate3
723 integer :: itf,jtf,ktf
725 itf=MIN(ite,ide-1)
726 ktf=MIN(kte,kde-1)
727 jtf=MIN(jte,jde-1)
729 !sms$distribute end
730 day=86400.
731 !
732 !--- specify entrainmentrate and detrainmentrate
733 !
734 radius=12000.
735 !
736 !--- gross entrainment rate (these may be changed later on in the
737 !--- program, depending what your detrainment is!!)
738 !
739 entr_rate=.2/radius
740 entr_rate3=.2/200.
741 !
742 !--- entrainment of mass
743 !
744 mentrd_rate=0.
745 mentr_rate=entr_rate
746 mentr_rate3=entr_rate3
747 !
748 !--- initial detrainmentrates
749 !
750 do k=kts,ktf
751 do i=its,itf
752 cd(i,k)=0.1*entr_rate
753 cd3(i,k)=entr_rate3
754 cdd(i,k)=0.
755 cdd3(i,k)=0.
756 clw_all(i,k)=0.
757 enddo
758 enddo
759 do ki=1,num_chem
760 do k=kts,ktf
761 do i=its,itf
762 tr_dd3(i,k,ki)=0.
763 enddo
764 enddo
765 enddo
766 !
767 !--- max/min allowed value for epsilon (ratio downdraft base mass flux/updraft
768 ! base mass flux
769 !
770 edtmax=.8
771 edtmin=.2
772 !
773 !--- minimum depth (m), clouds must have
774 !
775 depth_min=500.
776 !
777 !--- maximum depth (mb) of capping
778 !--- inversion (larger cap = no convection)
779 !
780 cap_maxs=175.
781 !sms$to_local(grid_dh: <1, mix :size>, <2, mjx :size>) begin
782 DO 7 i=its,itf
783 kbmax(i)=1
784 cap_max_increment(i)=0.
785 edt(i)=0.
786 edt3(i)=0.
787 kstabm(i)=ktf-1
788 IERR(i)=0
789 IERR5(i)=0
790 if(ktop3(i).lt.2)ierr5(i)=20
791 if(xmb3(i).eq.0.)ierr5(i)=21
792 7 CONTINUE
793 do i=its,itf
794 cap_max(i)=cap_maxs
795 do k=kts,kte
796 zd3(i,k)=0.
797 enddo
798 enddo
799 !
800 !--- max height(m) above ground where updraft air can originate
801 !
802 zkbmax=4000.
803 !
804 !--- height(m) above which no downdrafts are allowed to originate
805 !
806 zcutdown=3000.
807 !
808 !--- depth(m) over which downdraft detrains all its mass
809 !
810 z_detr=1250.
811 !
812 mbdt=dtime*4.E-03
813 !
814 !--- calculate moist static energy, heights, qes
815 !
816 call cup_env(z,qes,he,hes,t,q,p,z1, &
817 psur,ierr,tcrit,0,xl,cp, &
818 ids,ide, jds,jde, kds,kde, &
819 ims,ime, jms,jme, kms,kme, &
820 its,ite, jts,jte, kts,kte)
821 !
822 !--- environmental values on cloud levels
823 !
824 call cup_env_clev(t,qes,q,he,hes,z,p,qes_cup,q_cup,he_cup, &
825 hes_cup,z_cup,p_cup,gamma_cup,t_cup,psur, &
826 ierr,z1,xl,rv,cp, &
827 ids,ide, jds,jde, kds,kde, &
828 ims,ime, jms,jme, kms,kme, &
829 its,ite, jts,jte, kts,kte)
830 call cup_env_clev_tr(tracer,tre_cup,num_chem,ierr, &
831 ids,ide, jds,jde, kds,kde, &
832 ims,ime, jms,jme, kms,kme, &
833 its,ite, jts,jte, kts,kte)
834 do i=its,itf
835 if(ierr(i).eq.0)then
836 !
837 do k=kts,ktf-2
838 if(z_cup(i,k).gt.zkbmax+z1(i))then
839 kbmax(i)=k
840 go to 25
841 endif
842 enddo
843 25 continue
844 !
845 !
846 !--- level where detrainment for downdraft starts
847 !
848 do k=kts,ktf
849 if(z_cup(i,k).gt.z_detr+z1(i))then
850 kdet(i)=k
851 go to 26
852 endif
853 enddo
854 26 continue
855 !
856 endif
857 enddo
858 !
859 !
860 !
861 !------- DETERMINE LEVEL WITH HIGHEST MOIST STATIC ENERGY CONTENT - K22
862 !
863 CALL cup_MAXIMI(HE_CUP,3,KBMAX,K22,ierr, &
864 ids,ide, jds,jde, kds,kde, &
865 ims,ime, jms,jme, kms,kme, &
866 its,ite, jts,jte, kts,kte)
867 DO 36 i=its,itf
868 IF(ierr(I).eq.0.)THEN
869 IF(K22(I).GE.KBMAX(i))ierr(i)=2
870 endif
871 36 CONTINUE
872 !
873 ! done with basic stuff that is needed everywhere, from here on do not do
874 ! deep convection where aaeq .ne.0
875 !
876 DO i=its,itf
877 if(aaeq(i).ne.0.)then
878 ierr(i)=20
879 endif
880 enddo
881 !
882 !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
883 !
884 call cup_kbcon(cap_max_increment,1,k22,kbcon,he_cup,hes_cup, &
885 ierr,kbmax,p_cup,cap_max, &
886 ids,ide, jds,jde, kds,kde, &
887 ims,ime, jms,jme, kms,kme, &
888 its,ite, jts,jte, kts,kte)
889 !
890 !--- increase detrainment in stable layers
891 !
892 CALL cup_minimi(HEs_cup,Kbcon,kstabm,kstabi,ierr, &
893 ids,ide, jds,jde, kds,kde, &
894 ims,ime, jms,jme, kms,kme, &
895 its,ite, jts,jte, kts,kte)
896 do i=its,itf
897 IF(ierr(I).eq.0.)THEN
898 if(kstabm(i)-1.gt.kstabi(i))then
899 do k=kstabi(i),kstabm(i)-1
900 cd(i,k)=cd(i,k-1)+1.5*entr_rate
901 if(cd(i,k).gt.10.0*entr_rate)cd(i,k)=10.0*entr_rate
902 enddo
903 ENDIF
904 ENDIF
905 ENDDO
906 !
907 !--- calculate incloud moist static energy
908 !
909 call cup_up_he(k22,hkb,z_cup,cd,mentr_rate,he_cup,hc, &
910 kbcon,ierr,dby,he,hes_cup, &
911 ids,ide, jds,jde, kds,kde, &
912 ims,ime, jms,jme, kms,kme, &
913 its,ite, jts,jte, kts,kte)
915 !--- DETERMINE CLOUD TOP - KTOP
916 !
917 call cup_ktop(1,dby,kbcon,ktop,ierr, &
918 ids,ide, jds,jde, kds,kde, &
919 ims,ime, jms,jme, kms,kme, &
920 its,ite, jts,jte, kts,kte)
921 DO 37 i=its,itf
922 kzdown(i)=0
923 if(ierr(i).eq.0)then
924 zktop=(z_cup(i,ktop(i))-z1(i))*.6
925 zktop=min(zktop+z1(i),zcutdown+z1(i))
926 do k=kts,ktf
927 if(z_cup(i,k).gt.zktop)then
928 kzdown(i)=k
929 go to 37
930 endif
931 enddo
932 endif
933 37 CONTINUE
934 !
935 !--- DOWNDRAFT ORIGINATING LEVEL - JMIN
936 !
937 call cup_minimi(HEs_cup,K22,kzdown,JMIN,ierr, &
938 ids,ide, jds,jde, kds,kde, &
939 ims,ime, jms,jme, kms,kme, &
940 its,ite, jts,jte, kts,kte)
941 DO 100 i=its,ite
942 IF(ierr(I).eq.0.)THEN
943 if(jmin(i).le.3)then
944 ierr(i)=9
945 go to 100
946 endif
947 !
948 !--- check whether it would have buoyancy, if there where
949 !--- no entrainment/detrainment
950 !
951 101 continue
952 if(jmin(i)-1.lt.KDET(I))kdet(i)=jmin(i)-1
953 if(jmin(i).ge.Ktop(I)-1)jmin(i)=ktop(i)-2
954 ki=jmin(i)
955 hcd(i,ki)=hes_cup(i,ki)
956 DZ=Z_cup(i,Ki+1)-Z_cup(i,Ki)
957 dh=dz*(HCD(i,Ki)-hes_cup(i,ki))
958 dh=0.
959 !
960 do k=ki-1,1,-1
961 hcd(i,k)=hes_cup(i,jmin(i))
962 DZ=Z_cup(i,K+1)-Z_cup(i,K)
963 dh=dh+dz*(HCD(i,K)-hes_cup(i,k))
964 if(dh.gt.0.)then
965 jmin(i)=jmin(i)-1
966 if(jmin(i).gt.3)then
967 go to 101
968 else if(jmin(i).le.3)then
969 ierr(i)=9
970 go to 100
971 endif
972 endif
973 enddo
975 IF(JMIN(I).LE.3)then
976 ierr(i)=4
977 endif
979 ENDIF
980 100 continue
981 !
982 ! - Must have at least depth_min m between cloud convective base
983 ! and cloud top.
984 !
985 do i=its,itf
986 IF(ierr(I).eq.0.)THEN
987 if(jmin(i).le.3)then
988 ierr(i)=9
989 endif
990 IF(-z_cup(I,KBCON(I))+z_cup(I,KTOP(I)).LT.depth_min)then
991 ierr(i)=6
992 endif
993 endif
994 enddo
996 !
997 !c--- normalized updraft mass flux profile
998 !
999 call cup_up_nms(zu,z_cup,mentr_rate,cd,kbcon,ktop,ierr,k22, &
1000 ids,ide, jds,jde, kds,kde, &
1001 ims,ime, jms,jme, kms,kme, &
1002 its,ite, jts,jte, kts,kte)
1003 !
1004 !c--- normalized downdraft mass flux profile,also work on bottom detrainment
1005 !--- in this routine
1006 !
1007 call cup_dd_nms(zd,z_cup,cdd,mentrd_rate,jmin,ierr, &
1008 0,kdet,z1, &
1009 ids,ide, jds,jde, kds,kde, &
1010 ims,ime, jms,jme, kms,kme, &
1011 its,ite, jts,jte, kts,kte)
1012 !
1013 !--- downdraft moist static energy
1014 !
1015 call cup_dd_he(hes_cup,zd,hcd,z_cup,cdd,mentrd_rate, &
1016 jmin,ierr,he,dbyd,he_cup, &
1017 ids,ide, jds,jde, kds,kde, &
1018 ims,ime, jms,jme, kms,kme, &
1019 its,ite, jts,jte, kts,kte)
1020 !
1021 !--- calculate moisture properties of downdraft
1022 !
1023 call cup_dd_moisture(zd,hcd,hes_cup,qcd,qes_cup, &
1024 pwd,q_cup,z_cup,cdd,mentrd_rate,jmin,ierr,gamma_cup, &
1025 pwev,bu,qrcd,q,he,t_cup,1,xl, &
1026 ids,ide, jds,jde, kds,kde, &
1027 ims,ime, jms,jme, kms,kme, &
1028 its,ite, jts,jte, kts,kte)
1029 !
1030 !--- calculate moisture properties of updraft
1031 !
1032 call cup_up_moisture(ierr,z_cup,qc,qrc,pw,pwav, &
1033 kbcon,ktop,cd,dby,mentr_rate,clw_all, &
1034 q,GAMMA_cup,zu,qes_cup,k22,q_cup,xl, &
1035 ids,ide, jds,jde, kds,kde, &
1036 ims,ime, jms,jme, kms,kme, &
1037 its,ite, jts,jte, kts,kte)
1038 !
1039 !--- DETERMINE DOWNDRAFT STRENGTH IN TERMS OF WINDSHEAR
1040 !
1041 call cup_dd_edt(ierr,us,vs,z,ktop,kbcon,edt,p,pwav, &
1042 pwev,edtmax,edtmin,3,edtc, &
1043 ids,ide, jds,jde, kds,kde, &
1044 ims,ime, jms,jme, kms,kme, &
1045 its,ite, jts,jte, kts,kte)
1046 do i=its,itf
1047 if(ierr(i).eq.0)then
1048 edt(i)=edtc(i,2)
1049 endif
1050 enddo
1051 !
1052 ! massflux from precip and normalized cloud properties
1053 !
1054 pwdper=0.
1055 do i=its,itf
1057 if(ierr(i).gt.0)pre(i)=0.
1058 if(ierr(i).eq.0)then
1059 xmb(i)=pre(i)/(pwav(i)+edt(i)*pwev(i))
1060 !
1061 !--- percent of that that is evaporated (pwd is negative)
1062 !
1063 if(i.eq.ipr.and.j.eq.jpr)then
1064 print *,'xmb,edt,pwav = ',xmb(i),edt(i),pwav(i)
1065 print *,'k,pwdper(i,k),pw,pwd(i,k)',z1(i)
1066 endif
1067 do k=1,ktop(i)
1068 pwdper(i,k)=-edt(i)*pwd(i,k)/pwav(i)
1069 if(i.eq.ipr.and.j.eq.jpr)then
1070 print *,k,pwdper(i,k),pw(i,k),pwd(i,k)
1071 endif
1072 enddo
1073 endif
1074 enddo
1075 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1076 !
1077 !!!!! NOW WE HAVE EVREYTHING TO CALCULATE TRACER TRANSPORT AND WET DEPOSITION !!!
1078 !
1079 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1080 !
1081 !--- calculate incloud tracer distribution
1082 !
1083 ! if(j.eq.jpr)print *,'calling up_tracer'
1084 call cup_up_tracer(ierr,tcrit,t,pre,z_cup,p,tracer,tre_cup,tr_up,tr_pw, &
1085 tr_c,hstary,pw,clw_all,kbcon,ktop,cd,mentr_rate,zu,k22,&
1086 numgas,chemopt,scalaropt,conv_tr_wetscav,conv_tr_aqchem, &
1087 ids,ide, jds,jde, kds,kde, &
1088 ims,ime, jms,jme, kms,kme, &
1089 its,ite, jts,jte, kts,kte, &
1090 ipr,jpr,j,npr,num_chem,'deep')
1091 ! if(j.eq.jpr)print *,'called up_tracer'
1092 !
1093 ! for shallow convection, only 3 subroutines required
1094 !
1095 if(ishallow.eq.1)then
1096 call cup_up_nms(zu3,z_cup,mentr_rate3,cd3,kbcon3,ktop3,ierr5,k23, &
1097 ids,ide, jds,jde, kds,kde, &
1098 ims,ime, jms,jme, kms,kme, &
1099 its,ite, jts,jte, kts,kte)
1100 call cup_up_tracer(ierr5,tcrit,t,pre,z_cup,p,tracer,tre_cup,tr3_up,tr3_pw, &
1101 tr3_c,hstary,pw3,clw_all3,kbcon3,ktop3,cd3,mentr_rate3,zu3,k23,&
1102 numgas,chemopt,scalaropt,conv_tr_wetscav,conv_tr_aqchem, &
1103 ids,ide, jds,jde, kds,kde, &
1104 ims,ime, jms,jme, kms,kme, &
1105 its,ite, jts,jte, kts,kte, &
1106 ipr,jpr,j,npr,num_chem,'shallow')
1107 call cup_dellas_tr(ierr5,z_cup,p_cup,tr_dd3,edt3,zd3,cdd3, &
1108 tracer,tracert3,j,mentrd_rate,zu3,g,xmb3, &
1109 cd3,tr3_up,ktop3,k23,kbcon3,mentr_rate3,jmin,tre_cup,kdet, &
1110 k23,ipr,jpr,npr,'shallow',num_chem, &
1111 ids,ide, jds,jde, kds,kde, &
1112 ims,ime, jms,jme, kms,kme, &
1113 its,ite, jts,jte, kts,kte )
1114 endif
1116 call cup_dd_tracer(ierr,z_cup,qrcd,tracer,tre_cup,tr_up,tr_dd, &
1117 tr_pw,tr_pwd,jmin,cdd,mentrd_rate,zd,pwdper,wetdep,xmb,k22, &
1118 numgas,num_chem,ids,ide, jds,jde, kds,kde, &
1119 ims,ime, jms,jme, kms,kme, &
1120 its,ite, jts,jte, kts,kte)
1121 if(j.eq.jpr)print *,'called dd_tracer'
1124 !
1125 if(j.eq.jpr)then
1126 i=ipr
1127 print *,'in 250 loop ',edt(ipr),ierr(ipr)
1128 ! if(ierr(i).eq.0.or.ierr(i).eq.3)then
1129 print *,k22(I),kbcon(i),ktop(i),jmin(i)
1130 print *,edt(i)
1131 do k=kts,ktf
1132 print *,k,z(i,k),he(i,k),hes(i,k)
1133 enddo
1134 do k=1,ktop(i)+1
1135 print *,zu(i,k),zd(i,k),pw(i,k),pwd(i,k)
1136 enddo
1137 print *,'tr_up(i,k,6),tr_dd(i,k,6),tr_pw(i,k,6),tr_pwd(i,k,6)'
1138 do k=1,ktop(i)+1
1139 print *,tr_up(i,k,npr),tr_dd(i,k,npr),tr_pw(i,k,npr),tr_pwd(i,k,npr)
1140 enddo
1141 endif
1142 ! endif
1143 !
1144 !--- calculate transport tendencies
1145 !
1146 !--- 1. in bottom layer
1147 !
1148 call cup_dellabot_tr(ipr,jpr,tre_cup,ierr,z_cup,p,tr_dd,edt, &
1149 zd,cdd,tracer,tracert,j,mentrd_rate,z,g,xmb, &
1150 num_chem,ids,ide, jds,jde, kds,kde, &
1151 ims,ime, jms,jme, kms,kme, &
1152 its,ite, jts,jte, kts,kte)
1153 !
1154 !--- 2. everywhere else
1155 !
1157 call cup_dellas_tr(ierr,z_cup,p_cup,tr_dd,edt,zd,cdd, &
1158 tracer,tracert,j,mentrd_rate,zu,g,xmb, &
1159 cd,tr_up,ktop,k22,kbcon,mentr_rate,jmin,tre_cup,kdet, &
1160 k22,ipr,jpr,npr,'deep',num_chem, &
1161 ids,ide, jds,jde, kds,kde, &
1162 ims,ime, jms,jme, kms,kme, &
1163 its,ite, jts,jte, kts,kte )
1164 if(j.eq.jpr)then
1165 i=ipr
1166 do k=kts,ktf
1167 print *,k,tracer(i,k,npr),tracert(i,k,npr)
1168 enddo
1169 endif
1170 !
1171 ! may need more below for wet deposition......
1172 !
1173 !
1174 ! call cup_output_wd ( &
1175 ! ids,ide, jds,jde, kds,kde, &
1176 ! ims,ime, jms,jme, kms,kme, &
1177 ! its,ite, jts,jte, kts,kte)
1179 END SUBROUTINE CUP_CT
1181 SUBROUTINE cup_dellabot_tr(ipr,jpr,tre_cup,ierr,z_cup,p_cup, &
1182 tr_dd,edt,zd,cdd,tracer,tracert,j,mentrd_rate,z,g,xmb, &
1183 num_chem,ids,ide, jds,jde, kds,kde, &
1184 ims,ime, jms,jme, kms,kme, &
1185 its,ite, jts,jte, kts,kte )
1187 IMPLICIT NONE
1189 integer &
1190 ,intent (in ) :: &
1191 num_chem,ids,ide, jds,jde, kds,kde, &
1192 ims,ime, jms,jme, kms,kme, &
1193 its,ite, jts,jte, kts,kte
1194 integer, intent (in ) :: &
1195 j,ipr,jpr
1196 !
1197 ! ierr error value, maybe modified in this routine
1198 !
1199 real, dimension (its:ite,kts:kte,1:num_chem) &
1200 ,intent (out ) :: &
1201 tracert
1202 real, dimension (its:ite,kts:kte,1:num_chem) &
1203 ,intent (in ) :: &
1204 tre_cup,tracer,tr_dd
1205 real, dimension (its:ite,kts:kte) &
1206 ,intent (in ) :: &
1207 z_cup,p_cup,zd,cdd,z
1208 real, dimension (its:ite) &
1209 ,intent (in ) :: &
1210 edt,xmb
1211 real &
1212 ,intent (in ) :: &
1213 g,mentrd_rate
1214 integer, dimension (its:ite) &
1215 ,intent (inout) :: &
1216 ierr
1217 !
1218 ! local variables in this routine
1219 !
1221 integer i
1222 real detdo,detdo1,detdo2,entdo,dp,dz,subin, &
1223 totmas
1224 !
1225 integer :: itf, ktf, nv, npr
1226 npr=24
1227 itf=MIN(ite,ide-1)
1228 ktf=MIN(kte,kde-1)
1229 !
1230 !
1231 if(j.eq.jpr)print *,'in cup dellabot '
1232 tracert=0.
1233 do 100 i=its,itf
1234 if(ierr(i).ne.0)go to 100
1235 dz=z_cup(i,2)-z_cup(i,1)
1236 DP=100.*(p_cup(i,1)-P_cup(i,2))
1237 detdo1=edt(i)*zd(i,2)*CDD(i,1)*DZ
1238 detdo2=edt(i)*zd(i,1)
1239 entdo=edt(i)*zd(i,2)*mentrd_rate*dz
1240 subin=-EDT(I)*zd(i,2)
1241 detdo=detdo1+detdo2-entdo+subin
1242 do nv=2,num_chem
1243 tracert(I,1,nv)=(detdo1*.5*(tr_dd(i,1,nv)+tr_dd(i,2,nv)) &
1244 +detdo2*tr_dd(i,1,nv) &
1245 +subin*tre_cup(i,2,nv) &
1246 -entdo*tracer(i,1,nv))*g/dp*xmb(i)
1247 enddo
1248 if(j.eq.jpr.and.i.eq.ipr)print *,'in cup dellabot ',tracert(I,1,npr), &
1249 detdo1,detdo2,subin,entdo,tr_dd(i,1,npr),tr_dd(i,2,npr),tracer(i,1,npr)
1250 100 CONTINUE
1252 END SUBROUTINE cup_dellabot_tr
1255 SUBROUTINE cup_dellas_tr(ierr,z_cup,p_cup,tr_dd,edt,zd,cdd, &
1256 tracer,tracert,j,mentrd_rate,zu,g,xmb, &
1257 cd,tr_up,ktop,k22,kbcon,mentr_rate,jmin,tre_cup,kdet,kpbl, &
1258 ipr,jpr,npr,name,num_chem, &
1259 ids,ide, jds,jde, kds,kde, &
1260 ims,ime, jms,jme, kms,kme, &
1261 its,ite, jts,jte, kts,kte )
1263 IMPLICIT NONE
1265 integer &
1266 ,intent (in ) :: &
1267 num_chem,ids,ide, jds,jde, kds,kde, &
1268 ims,ime, jms,jme, kms,kme, &
1269 its,ite, jts,jte, kts,kte
1270 integer, intent (in ) :: &
1271 j,ipr,jpr,npr
1272 !
1273 ! ierr error value, maybe modified in this routine
1274 !
1275 real, dimension (its:ite,kts:kte,1:num_chem) &
1276 ,intent (inout ) :: &
1277 tracert
1278 real, dimension (its:ite,kts:kte,1:num_chem) &
1279 ,intent (in ) :: &
1280 tr_up,tr_dd,tre_cup,tracer
1281 real, dimension (its:ite,kts:kte) &
1282 ,intent (in ) :: &
1283 z_cup,p_cup,zd,cdd,cd,zu
1284 real, dimension (its:ite) &
1285 ,intent (in ) :: &
1286 edt,xmb
1287 real &
1288 ,intent (in ) :: &
1289 g,mentrd_rate,mentr_rate
1290 integer, dimension (its:ite) &
1291 ,intent (in ) :: &
1292 kbcon,ktop,k22,jmin,kdet,kpbl
1293 integer, dimension (its:ite) &
1294 ,intent (inout) :: &
1295 ierr
1296 character *(*), intent (in) :: &
1297 name
1298 !
1299 ! local variables in this routine
1300 !
1302 integer i,k,nv
1303 real detdo1,detdo2,entdo,dp,dz,subin,detdo,entup, &
1304 detup,subdown,entdoj,entupk,detupk,totmas
1305 !
1306 integer :: itf, ktf, kstart
1307 ! npr=24
1308 itf=MIN(ite,ide-1)
1309 ktf=MIN(kte,kde-1)
1310 kstart=kts+1
1311 if(name.eq.'shallow')kstart=kts
1312 !
1313 !
1314 i=ipr
1315 if(j.eq.jpr)then
1316 print *,'in dellas kpbl(i),k22(i),kbcon(i),ktop(i),jmin(i)'
1317 print *,kpbl(i),k22(i),kbcon(i),ktop(i),jmin(i)
1318 endif
1319 do nv=2,num_chem
1320 DO K=kstart,kte
1321 do i=its,itf
1322 tracert(i,k,nv)=0.
1323 enddo
1324 enddo
1325 enddo
1326 !
1327 DO 100 k=kts+1,ktf-1
1328 DO 100 i=its,ite
1329 IF(ierr(i).ne.0)GO TO 100
1330 IF(K.Gt.KTOP(I))GO TO 100
1331 !
1332 !--- SPECIFY DETRAINMENT OF DOWNDRAFT, HAS TO BE CONSISTENT
1333 !--- WITH ZD CALCULATIONS IN SOUNDD.
1334 !
1335 DZ=Z_cup(I,K+1)-Z_cup(I,K)
1336 detdo=edt(i)*CDD(i,K)*DZ*ZD(i,k+1)
1337 entdo=edt(i)*mentrd_rate*dz*zd(i,k+1)
1338 subin=zu(i,k+1)-zd(i,k+1)*edt(i)
1339 entup=0.
1340 detup=0.
1341 if(k.ge.kbcon(i).and.k.lt.ktop(i))then
1342 entup=mentr_rate*dz*zu(i,k)
1343 detup=CD(i,K+1)*DZ*ZU(i,k)
1344 endif
1345 subdown=(zu(i,k)-zd(i,k)*edt(i))
1346 entdoj=0.
1347 entupk=0.
1348 detupk=0.
1349 !
1350 if(k.eq.jmin(i))then
1351 entdoj=edt(i)*zd(i,k)
1352 endif
1354 if(k.eq.k22(i)-1)then
1355 entupk=zu(i,kpbl(i))
1356 endif
1358 if(k.gt.kdet(i))then
1359 detdo=0.
1360 endif
1362 if(k.eq.ktop(i)-0)then
1363 detupk=zu(i,ktop(i))
1364 subin=0.
1365 endif
1366 if(k.lt.kbcon(i))then
1367 detup=0.
1368 endif
1369 !C
1370 !C--- CHANGED DUE TO SUBSIDENCE AND ENTRAINMENT
1371 !C
1372 totmas=subin-subdown+detup-entup-entdo+ &
1373 detdo-entupk-entdoj+detupk
1374 if(j.eq.jpr.and.i.eq.ipr)print *,'k,totmas,sui,sud = ',k, &
1375 totmas,subin,subdown
1376 ! if(j.eq.jpr.and.i.eq.ipr)print *,'updr stuff = ',detup,
1377 ! 1 entup,entupk,detupk
1378 ! if(j.eq.jpr.and.i.eq.ipr)print *,'dddr stuff = ',entdo,
1379 ! 1 detdo,entdoj
1380 if(abs(totmas).gt.1.e-6)then
1381 print *,'*********************',i,j,k,totmas,name
1382 print *,kpbl(i),k22(i),kbcon(i),ktop(i)
1383 !c print *,'updr stuff = ',subin,
1384 !c 1 subdown,detup,entup,entupk,detupk
1385 !c print *,'dddr stuff = ',entdo,
1386 !c 1 detdo,entdoj
1387 CALL wrf_error_fatal ( 'cup_dellas_tr: TOTMAS > CRITICAL VALUE')
1388 endif
1389 dp=100.*(p_cup(i,k-1)-p_cup(i,k))
1390 do nv=2,num_chem
1391 ! tracert(i,k,nv)=(subin*tre_cup(i,k+1,nv) &
1392 ! -subdown*tre_cup(i,k,nv) &
1393 tracert(i,k,nv)=(subin*tracer(i,k+1,nv) &
1394 -subdown*tracer(i,k,nv) &
1395 +detup*.5*(tr_up(i,K+1,nv)+tr_up(i,K,nv)) &
1396 +detdo*.5*(tr_dd(i,K+1,nv)+tr_dd(i,K,nv)) &
1397 -entup*tracer(i,k,nv) &
1398 -entdo*tracer(i,k,nv) &
1399 -entupk*tre_cup(i,k22(i),nv) &
1400 -entdoj*tre_cup(i,jmin(i),nv) &
1401 +detupk*tr_up(i,ktop(i),nv) &
1402 )*g/dp*xmb(i)
1403 enddo
1404 if(i.eq.ipr.and.j.eq.jpr)then
1405 print *,k,tracert(i,k,npr),subin*tre_cup(i,k+1,npr),subdown*tre_cup(i,k,npr), &
1406 detdo*.5*(tr_dd(i,K+1,npr)+tr_dd(i,K,npr))
1407 print *,k,detup*.5*(tr_up(i,K+1,npr)+tr_up(i,K,npr)),detupk*tr_up(i,ktop(i),npr), &
1408 entup*tracer(i,k,npr),entdo*tracer(i,k,npr)
1409 print *,k,entupk*tre_cup(i,k,npr),detupk,tr_up(i,ktop(i),npr)
1410 endif
1412 100 CONTINUE
1414 END SUBROUTINE cup_dellas_tr
1415 SUBROUTINE cup_env_clev_tr(tracer,tre_cup,num_chem,ierr, &
1416 ids,ide, jds,jde, kds,kde, &
1417 ims,ime, jms,jme, kms,kme, &
1418 its,ite, jts,jte, kts,kte)
1419 implicit none
1420 integer &
1421 ,intent (in ) :: &
1422 num_chem,ids,ide, jds,jde, kds,kde, &
1423 ims,ime, jms,jme, kms,kme, &
1424 its,ite, jts,jte, kts,kte
1425 integer, dimension (its:ite) &
1426 ,intent (in) :: &
1427 ierr
1429 real, dimension (its:ite,kts:kte,1:num_chem) &
1430 ,intent (in ) :: &
1431 tracer
1432 real, dimension (its:ite,kts:kte,1:num_chem) &
1433 ,intent (out ) :: &
1434 tre_cup
1435 !
1436 ! local variables in this routine
1437 !
1439 integer :: &
1440 i,k,nv,itf,ktf
1441 itf=MIN(ite,ide-1)
1442 ktf=MIN(kte,kde-1)
1443 do nv=2,num_chem
1444 do k=kts+1,ktf
1445 do i=its,ite
1446 if(ierr(i).eq.0)then
1447 tre_cup(i,k,nv)=.5*(tracer(i,k-1,nv)+tracer(i,k,nv))
1448 endif
1449 enddo
1450 enddo
1451 enddo
1452 do nv=2,num_chem
1453 do i=its,ite
1454 if(ierr(i).eq.0)then
1455 tre_cup(i,kts,nv)=tracer(i,kts,nv)
1456 endif
1457 enddo
1458 enddo
1461 END subroutine cup_env_clev_tr
1464 SUBROUTINE cup_up_tracer(ierr,tcrit,t,pre,z_cup,p,tracer,tre_cup,tr_up, &
1465 tr_pw,tr_c,hstary,cupclw,clw_all,kbcon,ktop,cd,mentr_rate,zu,k22, &
1466 numgas,chemopt,scalaropt,conv_tr_wetscav,conv_tr_aqchem, &
1467 ids,ide, jds,jde, kds,kde, &
1468 ims,ime, jms,jme, kms,kme, &
1469 its,ite, jts,jte, kts,kte,ipr,jpr,j,npr,num_chem,name)
1470 ! USE module_configure
1471 !!!TUCCELLA
1472 USE module_state_description, only: RADM2SORG,RADM2SORG_AQ,RACMSORG_AQ,RACMSORG_KPP, &
1473 RADM2SORG_KPP,RACM_ESRLSORG_KPP,RACM_SOA_VBS_KPP, &
1474 RADM2SORG_AQCHEM,RACMSORG_AQCHEM_KPP,RACM_ESRLSORG_AQCHEM_KPP, &
1475 RACM_SOA_VBS_AQCHEM_KPP, &
1476 RACM_SOA_VBS_HET_KPP, &
1477 CB05_SORG_VBS_AQ_KPP
1478 USE module_ctrans_aqchem
1479 USE module_chem_share, only: get_last_gas
1480 USE module_HLawConst, only: HLC
1481 implicit none
1482 ! Aqeuous species pointers INCLUDE File
1484 !...........PARAMETERS and their descriptions:
1486 INTEGER, PARAMETER :: NGAS = 12 ! number of gas-phase species for AQCHEM
1487 INTEGER, PARAMETER :: NAER = 36 ! number of aerosol species for AQCHEM
1488 INTEGER, PARAMETER :: NLIQS = 41 ! number of liquid-phase species in AQCHEM
1490 !...pointers for the AQCHEM array GAS
1492 INTEGER, PARAMETER :: LSO2 = 1 ! Sulfur Dioxide
1493 INTEGER, PARAMETER :: LHNO3 = 2 ! Nitric Acid
1494 INTEGER, PARAMETER :: LN2O5 = 3 ! Dinitrogen Pentoxide
1495 INTEGER, PARAMETER :: LCO2 = 4 ! Carbon Dioxide
1496 INTEGER, PARAMETER :: LNH3 = 5 ! Ammonia
1497 INTEGER, PARAMETER :: LH2O2 = 6 ! Hydrogen Perioxide
1498 INTEGER, PARAMETER :: LO3 = 7 ! Ozone
1499 INTEGER, PARAMETER :: LFOA = 8 ! Formic Acid
1500 INTEGER, PARAMETER :: LMHP = 9 ! Methyl Hydrogen Peroxide
1501 INTEGER, PARAMETER :: LPAA = 10 ! Peroxyacidic Acid
1502 INTEGER, PARAMETER :: LH2SO4 = 11 ! Sulfuric Acid
1503 INTEGER, PARAMETER :: LHCL = 12 ! Hydrogen Chloride
1505 !...pointers for the AQCHEM array AEROSOL
1507 INTEGER, PARAMETER :: LSO4AKN = 1 ! Aitken-mode Sulfate
1508 INTEGER, PARAMETER :: LSO4ACC = 2 ! Accumulation-mode Sulfate
1509 INTEGER, PARAMETER :: LSO4COR = 3 ! Coarse-mode Sulfate
1510 INTEGER, PARAMETER :: LNH4AKN = 4 ! Aitken-mode Ammonium
1511 INTEGER, PARAMETER :: LNH4ACC = 5 ! Accumulation-mode Ammonium
1512 INTEGER, PARAMETER :: LNO3AKN = 6 ! Aitken-mode Nitrate
1513 INTEGER, PARAMETER :: LNO3ACC = 7 ! Accumulation-mode Nitrate
1514 INTEGER, PARAMETER :: LNO3COR = 8 ! Coarse-mode Nitrate
1515 INTEGER, PARAMETER :: LORGAAKN = 9 ! Aitken-mode anthropogenic SOA
1516 INTEGER, PARAMETER :: LORGAACC = 10 ! Accumulation-mode anthropogenic SOA
1517 INTEGER, PARAMETER :: LORGPAKN = 11 ! Aitken-mode primary organic aerosol
1518 INTEGER, PARAMETER :: LORGPACC = 12 ! Accumulation-mode primary organic aerosol
1519 INTEGER, PARAMETER :: LORGBAKN = 13 ! Aitken-mode biogenic SOA
1520 INTEGER, PARAMETER :: LORGBACC = 14 ! Accumulation-mode biogenic SOA
1521 INTEGER, PARAMETER :: LECAKN = 15 ! Aitken-mode elemental carbon
1522 INTEGER, PARAMETER :: LECACC = 16 ! Accumulation-mode elemental carbon
1523 INTEGER, PARAMETER :: LPRIAKN = 17 ! Aitken-mode primary aerosol
1524 INTEGER, PARAMETER :: LPRIACC = 18 ! Accumulation-mode primary aerosol
1525 INTEGER, PARAMETER :: LPRICOR = 19 ! Coarse-mode primary aerosol
1526 INTEGER, PARAMETER :: LNAAKN = 20 ! Aitken-mode Sodium
1527 INTEGER, PARAMETER :: LNAACC = 21 ! Accumulation-mode Sodium
1528 INTEGER, PARAMETER :: LNACOR = 22 ! Coarse-mode Sodium
1529 INTEGER, PARAMETER :: LCLAKN = 23 ! Aitken-mode Chloride ion
1530 INTEGER, PARAMETER :: LCLACC = 24 ! Accumulation-mode Chloride ion
1531 INTEGER, PARAMETER :: LCLCOR = 25 ! Coarse-mode Chloride ion
1532 INTEGER, PARAMETER :: LNUMAKN = 26 ! Aitken-mode number
1533 INTEGER, PARAMETER :: LNUMACC = 27 ! Accumulation-mode number
1534 INTEGER, PARAMETER :: LNUMCOR = 28 ! Coarse-mode number
1535 INTEGER, PARAMETER :: LSRFAKN = 29 ! Aitken-mode surface area
1536 INTEGER, PARAMETER :: LSRFACC = 30 ! Accumulation-mode surface area
1537 INTEGER, PARAMETER :: LNACL = 31 ! Sodium Chloride aerosol for AE3 only {depreciated in AE4}
1538 INTEGER, PARAMETER :: LCACO3 = 32 ! Calcium Carbonate aerosol (place holder)
1539 INTEGER, PARAMETER :: LMGCO3 = 33 ! Magnesium Carbonate aerosol (place holder)
1540 INTEGER, PARAMETER :: LA3FE = 34 ! Iron aerosol (place holder)
1541 INTEGER, PARAMETER :: LB2MN = 35 ! Manganese aerosol (place holder)
1542 INTEGER, PARAMETER :: LK = 36 ! Potassium aerosol (Cl- tracked separately) (place holder)
1544 !
1545 ! on input
1546 !
1548 ! only local wrf dimensions are need as of now in this routine
1550 integer &
1551 ,intent (in ) :: &
1552 ids,ide, jds,jde, kds,kde,scalaropt, &
1553 numgas, conv_tr_wetscav,conv_tr_aqchem, &
1554 num_chem,ims,ime, jms,jme, kms,kme,chemopt, &
1555 its,ite, jts,jte, kts,kte,ipr,jpr,j,npr
1556 real, dimension (its:ite,kts:kte) &
1557 ,intent (in ) :: &
1558 z_cup,cd,zu,p,hstary,t
1559 real, dimension (its:ite,kts:kte) &
1560 ,intent (inout ) :: &
1561 cupclw,clw_all
1562 real, dimension (its:ite,kts:kte,1:num_chem) &
1563 ,intent (inout ) :: &
1564 tr_up,tr_c,tr_pw
1565 real, dimension (its:ite,kts:kte,1:num_chem) &
1566 ,intent (in ) :: &
1567 tre_cup,tracer
1568 real, dimension (its:ite) &
1569 ,intent (in ) :: &
1570 pre
1572 ! entr= entrainment rate
1573 real &
1574 ,intent (in ) :: &
1575 mentr_rate,tcrit
1576 integer, dimension (its:ite) &
1577 ,intent (in ) :: &
1578 kbcon,ktop,k22
1579 ! ierr error value, maybe modified in this routine
1581 integer, dimension (its:ite) &
1582 ,intent (inout) :: &
1583 ierr
1584 character *(*), intent (in) :: &
1585 name
1586 ! local variables in this routine
1587 !
1588 real :: conc_equi,conc_mxr,partialp,taucld
1589 real :: HLCnst1, HLCnst2
1590 real :: HLCnst3, HLCnst4
1591 real :: HLCnst5, HLCnst6
1592 real :: kh, dk1s, dk2s, heff, tfac
1594 integer :: &
1595 iall,i,k,iwd,nv, HLndx
1596 real :: &
1597 trcc,trch,dh,c0,dz,radius,airm,dens
1598 integer :: &
1599 itf,ktf,iaer,igas
1601 real :: &
1602 frac_so4(4), frac_no3(4), frac_nh4(4), tot_so4, tot_nh4, tot_no3
1604 ! Gas/aqueous phase partitioning for wet scavenging/deposition of gas
1605 ! phase and aerosol species:
1606 real aq_gas_ratio
1608 ! I/O for AQCHEM:
1610 real precip ! Precipitation rate (mm/h)
1611 real, dimension (ngas) :: gas,gaswdep
1612 real, dimension (naer) :: aerosol,aerwdep
1613 real, dimension (nliqs) :: liquid
1614 real hpwdep
1615 real alfa0,alfa2,alfa3 ! Aerosol scavenging coefficients for Aitken mode
1617 real, parameter :: hion = 1.e-5
1618 real, parameter :: hion_inv = 1./hion
1619 real, parameter :: HL_t0 = 298.
1621 itf=MIN(ite,ide-1)
1622 ktf=MIN(kte,kde-1)
1624 !
1625 iall=0
1626 c0=.002
1627 iwd=0
1628 !
1629 !--- no precip for small clouds
1630 !
1631 ! if(mentr_rate.gt.0.)then
1632 ! radius=.2/mentr_rate
1633 ! if(radius.lt.900.)c0=0.
1634 if(name.eq.'shallow')c0=0.
1635 ! if(radius.lt.900.)iall=0
1636 ! endif
1637 do nv=2,num_chem
1638 do k=kts,ktf
1639 do i=its,itf
1640 ! Initialize species mass in rain water:
1641 tr_pw(i,k,nv)=0.
1642 ! Initialize species mass in updraft:
1643 if(ierr(i).eq.0)tr_up(i,k,nv)=tre_cup(i,k,nv)
1644 ! Initialize species mass in cloud + rain water:
1645 tr_c(i,k,nv)=0.
1646 enddo
1647 enddo
1648 enddo
1650 do nv=2,num_chem
1651 do i=its,itf
1652 if(ierr(i).eq.0.)then
1653 do k=k22(i),kbcon(i)-1
1654 tr_up(i,k,nv)=tre_cup(i,k22(i),nv)
1655 enddo
1656 endif
1657 enddo
1658 enddo
1660 DO 100 k=kts+1,ktf-1
1661 DO 100 i=its,itf
1663 IF(ierr(i).ne.0)GO TO 100
1664 IF(K.Lt.KBCON(I))GO TO 100
1665 IF(K.Gt.KTOP(I)+1)GO TO 100
1667 DZ=Z_cup(i,K)-Z_cup(i,K-1)
1669 if(cupclw(i,k).le.0.)cupclw(i,k)=0.
1670 if(clw_all(i,k).le.0.)clw_all(i,k)=0.
1672 !
1673 ! Steady state plume equation, for what could be in cloud before anything
1674 ! happens (kg/kg). tr_up would be the concentration if tracers were conserved.
1675 !
1677 do nv=2,num_chem
1678 if(i.eq.ipr.and.j.eq.jpr.and.nv.eq.npr)print *,k,tr_up(i,K-1,nv),tr_up(i,K,nv),tr_pw(i,k-1,nv),clw_all(i,k),cupclw(i,k)
1679 tr_up(i,K,nv)=(tr_up(i,K-1,nv)*(1.-.5*CD(i,K)*DZ)+mentr_rate* &
1680 DZ*tracer(i,K-1,nv))/(1.+mentr_rate*DZ-.5*cd(i,k)*dz)
1681 tr_up(i,k,nv)=max(1.e-16,tr_up(i,K,nv))
1682 enddo
1684 !
1685 ! Aqueous chemistry
1686 !
1687 !!! TUCCELLA
1688 if ((chemopt .EQ. RADM2SORG .OR. chemopt .EQ. RADM2SORG_AQ .OR. chemopt .EQ. RACMSORG_AQ .OR. &
1689 chemopt .EQ. RACMSORG_KPP .OR. chemopt .EQ. RADM2SORG_KPP .OR. chemopt .EQ. RACM_ESRLSORG_KPP .OR. &
1690 chemopt .EQ. RACM_SOA_VBS_KPP .OR. chemopt .EQ. RADM2SORG_AQCHEM .OR. chemopt .EQ. RACMSORG_AQCHEM_KPP .OR. &
1691 chemopt .EQ. CB05_SORG_VBS_AQ_KPP .OR. &
1692 chemopt .EQ. RACM_SOA_VBS_HET_KPP .OR. &
1693 chemopt .EQ. RACM_ESRLSORG_AQCHEM_KPP .OR. chemopt .EQ. RACM_SOA_VBS_AQCHEM_KPP) &
1694 .and. (conv_tr_aqchem == 1)) then
1696 !
1697 ! For MADE/SORGAM derived schemes with aqueous chemistry
1698 !
1700 ! Air mass density
1701 dens = 0.1*p(i,k)/t(i,k)*mwdry/8.314472 ! kg/m3
1703 ! Column air number density:
1704 airm = 1000.0*dens*dz/mwdry ! mol/m2
1706 ! Wet scavenging initialization for AQCHEM
1708 GASWDEP = 0.0
1709 AERWDEP = 0.0
1710 HPWDEP = 0.0
1712 ! We provide a precipitation rate and aerosol scavenging rates of zero,
1713 ! in order to prevent wet scavenging in AQCHEM (it is treated later):
1715 precip = 0.0 ! mm/hr
1717 alfa0 = 0.0
1718 alfa2 = 0.0
1719 alfa3 = 0.0
1721 ! Gas phase concentrations before aqueous phase chemistry
1722 ! (with units conversion ppm -> mol/mol)
1724 gas(:) = 0.0
1726 gas(lco2) = 380.0e-6
1728 gas(lso2) = tr_up(i,k,p_so2)*1.0e-6
1729 gas(lhno3) = tr_up(i,k,p_hno3)*1.0e-6
1730 gas(ln2o5) = tr_up(i,k,p_n2o5)*1.0e-6
1731 gas(lnh3) = tr_up(i,k,p_nh3)*1.0e-6
1732 gas(lh2o2) = tr_up(i,k,p_h2o2)*1.0e-6
1733 gas(lo3) = tr_up(i,k,p_o3)*1.0e-6
1734 gas(lh2so4) = tr_up(i,k,p_sulf)*1.0e-6
1735 if (chemopt==CB05_SORG_VBS_AQ_KPP) then
1736 gas(lfoa) = tr_up(i,k,p_facd)*1.0e-6
1737 gas(lmhp) = tr_up(i,k,p_mepx)*1.0e-6
1738 gas(lpaa) = tr_up(i,k,p_pacd)*1.0e-6
1739 else
1740 gas(lfoa) = tr_up(i,k,p_ora1)*1.0e-6
1741 gas(lmhp) = tr_up(i,k,p_op1)*1.0e-6
1742 gas(lpaa) = tr_up(i,k,p_paa)*1.0e-6
1743 end if
1745 ! Aerosol mass concentrations before aqueous phase chemistry
1746 ! (with units conversion ug/kg -> mol/mol). Although AQCHEM
1747 ! accounts for much of the aerosol compounds in MADE, they are
1748 ! not treated at the moment by AQCHEM, as the mapping between
1749 ! the organic compound groups in MADE and AQCHEM is not obvious.
1751 aerosol(:) = 0.0
1753 ! We assume all accumulation mode particles are activated in cumulus clouds:
1755 aerosol(lso4acc) = tr_up(i,k,p_so4aj)*1.0e-9*mwdry/mwso4
1756 aerosol(lnh4acc) = tr_up(i,k,p_nh4aj)*1.0e-9*mwdry/mwnh4
1757 aerosol(lno3acc) = tr_up(i,k,p_no3aj)*1.0e-9*mwdry/mwno3
1759 ! Cloud lifetime:
1760 taucld = 1800.0
1762 if (clw_all(i,k)*dens .gt. qcldwtr_cutoff) then ! Cloud water > threshold
1763 CALL AQCHEM( &
1764 t(i,k), &
1765 p(i,k)*100., &
1766 taucld, &
1767 precip, &
1768 clw_all(i,k)*dens, &
1769 clw_all(i,k)*dens, &
1770 airm, &
1771 ALFA0, &
1772 ALFA2, &
1773 ALFA3, &
1774 GAS, &
1775 AEROSOL, &
1776 LIQUID, &
1777 GASWDEP, &
1778 AERWDEP, &
1779 HPWDEP)
1780 endif
1782 ! Gas phase concentrations after aqueous phase chemistry
1783 ! (with units conversion mol/mol -> ppm)
1785 tr_up(i,k,p_so2) = gas(lso2)*1.0e6
1786 tr_up(i,k,p_hno3) = gas(lhno3)*1.0e6
1787 tr_up(i,k,p_n2o5) = gas(ln2o5)*1.0e6
1788 tr_up(i,k,p_nh3) = gas(lnh3)*1.0e6
1789 tr_up(i,k,p_h2o2) = gas(lh2o2)*1.0e6
1790 tr_up(i,k,p_o3) = gas(lo3)*1.0e6
1791 tr_up(i,k,p_sulf) = gas(lh2so4)*1.0e6
1792 if (chemopt==CB05_SORG_VBS_AQ_KPP) then
1793 tr_up(i,k,p_facd) = gas(lfoa)*1.0e6
1794 tr_up(i,k,p_mepx) = gas(lmhp)*1.0e6
1795 tr_up(i,k,p_pacd) = gas(lpaa)*1.0e6
1796 else
1797 tr_up(i,k,p_ora1) = gas(lfoa)*1.0e6
1798 tr_up(i,k,p_op1) = gas(lmhp)*1.0e6
1799 tr_up(i,k,p_paa) = gas(lpaa)*1.0e6
1800 end if
1802 ! Aerosol mass concentrations
1803 ! (with units conversion mol/mol -> ug/kg)
1805 tr_up(i,k,p_so4aj) = aerosol(lso4acc)*1.0e9*mwso4/mwdry
1806 tr_up(i,k,p_nh4aj) = aerosol(lnh4acc)*1.0e9*mwnh4/mwdry
1807 tr_up(i,k,p_no3aj) = aerosol(lno3acc)*1.0e9*mwno3/mwdry
1809 else if ((chemopt .EQ. mozart_mosaic_4bin_kpp .OR. &
1810 chemopt .EQ. mozart_mosaic_4bin_aq_kpp) &
1811 .AND. (conv_tr_aqchem == 1)) then
1813 !
1814 ! For MOSAIC 4bin scheme with aqueous chemistry
1815 !
1817 ! Air mass density
1818 dens = 0.1*p(i,k)/t(i,k)*mwdry/8.314472 ! kg/m3
1820 ! Column air number density:
1821 airm = 1000.0*dens*dz/mwdry ! mol/m2
1823 ! Wet scavenging initialization for AQCHEM
1825 GASWDEP = 0.0
1826 AERWDEP = 0.0
1827 HPWDEP = 0.0
1829 ! We provide a precipitation rate and aerosol scavenging rates of zero,
1830 ! in order to prevent wet scavenging in AQCHEM (it is treated later):
1832 precip = 0.0 ! mm/hr
1834 alfa0 = 0.0
1835 alfa2 = 0.0
1836 alfa3 = 0.0
1838 ! Gas phase concentrations before aqueous phase chemistry
1839 ! (with units conversion ppm -> mol/mol)
1841 gas(:) = 0.0
1843 gas(lco2) = 380.0e-6
1845 gas(lso2) = tr_up(i,k,p_so2)*1.0e-6
1846 gas(lhno3) = tr_up(i,k,p_hno3)*1.0e-6
1847 gas(ln2o5) = tr_up(i,k,p_n2o5)*1.0e-6
1848 gas(lnh3) = tr_up(i,k,p_nh3)*1.0e-6
1849 gas(lh2o2) = tr_up(i,k,p_h2o2)*1.0e-6
1850 gas(lo3) = tr_up(i,k,p_o3)*1.0e-6
1851 gas(lfoa) = tr_up(i,k,p_hcooh)*1.0e-6
1852 gas(lmhp) = tr_up(i,k,p_ch3ooh)*1.0e-6
1853 gas(lpaa) = tr_up(i,k,p_paa)*1.0e-6
1854 gas(lh2so4) = tr_up(i,k,p_sulf)*1.0e-6
1856 ! Aerosol mass concentrations before aqueous phase chemistry
1857 ! (with units conversion ug/kg -> mol/mol). Although AQCHEM
1858 ! accounts for much of the aerosol compounds in MADE, they are
1859 ! not treated at the moment by AQCHEM, as the mapping between
1860 ! the organic compound groups in MADE and AQCHEM is not obvious.
1862 aerosol(:) = 0.0
1864 ! We assume all particles in bins 2 - 4 are activated in cumulus clouds:
1866 ! remember size distribution
1867 ! (if none existed before, frac_x is not set, hence distribute equally as default)
1868 frac_so4(:) = 0.25
1869 frac_nh4(:) = 0.25
1870 frac_no3(:) = 0.25
1872 tot_so4 = tr_up(i,k,p_so4_a01)+tr_up(i,k,p_so4_a02)+&
1873 tr_up(i,k,p_so4_a03)+tr_up(i,k,p_so4_a04)
1874 tot_nh4 = tr_up(i,k,p_nh4_a01)+tr_up(i,k,p_nh4_a02)+&
1875 tr_up(i,k,p_nh4_a03)+tr_up(i,k,p_nh4_a04)
1876 tot_no3 = tr_up(i,k,p_no3_a01)+tr_up(i,k,p_no3_a02)+&
1877 tr_up(i,k,p_no3_a03)+tr_up(i,k,p_no3_a04)
1879 if (tot_so4 > 0.0) then
1880 frac_so4(1) = tr_up(i,k,p_so4_a01) / tot_so4
1881 frac_so4(2) = tr_up(i,k,p_so4_a02) / tot_so4
1882 frac_so4(3) = tr_up(i,k,p_so4_a03) / tot_so4
1883 frac_so4(4) = tr_up(i,k,p_so4_a04) / tot_so4
1884 aerosol(lso4acc) = tot_so4 *1.0e-9*mwdry/mwso4
1885 end if
1887 if (tot_nh4 > 0.0) then
1888 frac_nh4(1) = tr_up(i,k,p_nh4_a01) / tot_nh4
1889 frac_nh4(2) = tr_up(i,k,p_nh4_a02) / tot_nh4
1890 frac_nh4(3) = tr_up(i,k,p_nh4_a03) / tot_nh4
1891 frac_nh4(4) = tr_up(i,k,p_nh4_a04) / tot_nh4
1892 aerosol(lnh4acc) = tot_nh4 *1.0e-9*mwdry/mwnh4
1893 end if
1895 if (tot_no3 > 0.0) then
1896 frac_no3(1) = tr_up(i,k,p_no3_a01) / tot_no3
1897 frac_no3(2) = tr_up(i,k,p_no3_a02) / tot_no3
1898 frac_no3(3) = tr_up(i,k,p_no3_a03) / tot_no3
1899 frac_no3(4) = tr_up(i,k,p_no3_a04) / tot_no3
1900 aerosol(lno3acc) = tot_no3 *1.0e-9*mwdry/mwno3
1901 end if
1903 ! Cloud lifetime:
1904 taucld = 1800.0
1906 if (clw_all(i,k)*dens .gt. qcldwtr_cutoff) then ! Cloud water > threshold
1907 CALL AQCHEM( &
1908 t(i,k), &
1909 p(i,k)*100., &
1910 taucld, &
1911 precip, &
1912 clw_all(i,k)*dens, &
1913 clw_all(i,k)*dens, &
1914 airm, &
1915 ALFA0, &
1916 ALFA2, &
1917 ALFA3, &
1918 GAS, &
1919 AEROSOL, &
1920 LIQUID, &
1921 GASWDEP, &
1922 AERWDEP, &
1923 HPWDEP)
1924 endif
1926 ! Gas phase concentrations after aqueous phase chemistry
1927 ! (with units conversion mol/mol -> ppm)
1929 tr_up(i,k,p_so2) = gas(lso2)*1.0e6
1930 tr_up(i,k,p_hno3) = gas(lhno3)*1.0e6
1931 tr_up(i,k,p_n2o5) = gas(ln2o5)*1.0e6
1932 tr_up(i,k,p_nh3) = gas(lnh3)*1.0e6
1933 tr_up(i,k,p_h2o2) = gas(lh2o2)*1.0e6
1934 tr_up(i,k,p_o3) = gas(lo3)*1.0e6
1935 tr_up(i,k,p_hcooh) = gas(lfoa)*1.0e6
1936 tr_up(i,k,p_ch3ooh) = gas(lmhp)*1.0e6
1937 tr_up(i,k,p_paa) = gas(lpaa)*1.0e6
1938 tr_up(i,k,p_sulf) = gas(lh2so4)*1.0e6
1940 ! Aerosol mass concentrations
1941 ! (with units conversion mol/mol -> ug/kg)
1943 tr_up(i,k,p_so4_a01) = aerosol(lso4acc) * frac_so4(1) * 1.0e9*mwso4/mwdry
1944 tr_up(i,k,p_so4_a02) = aerosol(lso4acc) * frac_so4(2) * 1.0e9*mwso4/mwdry
1945 tr_up(i,k,p_so4_a03) = aerosol(lso4acc) * frac_so4(3) * 1.0e9*mwso4/mwdry
1946 tr_up(i,k,p_so4_a04) = aerosol(lso4acc) * frac_so4(4) * 1.0e9*mwso4/mwdry
1948 tr_up(i,k,p_nh4_a01) = aerosol(lnh4acc) * frac_nh4(1) * 1.0e9*mwnh4/mwdry
1949 tr_up(i,k,p_nh4_a02) = aerosol(lnh4acc) * frac_nh4(2) * 1.0e9*mwnh4/mwdry
1950 tr_up(i,k,p_nh4_a03) = aerosol(lnh4acc) * frac_nh4(3) * 1.0e9*mwnh4/mwdry
1951 tr_up(i,k,p_nh4_a04) = aerosol(lnh4acc) * frac_nh4(4) * 1.0e9*mwnh4/mwdry
1953 tr_up(i,k,p_no3_a01) = aerosol(lno3acc) * frac_no3(1) * 1.0e9*mwno3/mwdry
1954 tr_up(i,k,p_no3_a02) = aerosol(lno3acc) * frac_no3(2) * 1.0e9*mwno3/mwdry
1955 tr_up(i,k,p_no3_a03) = aerosol(lno3acc) * frac_no3(3) * 1.0e9*mwno3/mwdry
1956 tr_up(i,k,p_no3_a04) = aerosol(lno3acc) * frac_no3(4) * 1.0e9*mwno3/mwdry
1958 endif
1960 ! wet scavenging option (turn off by setting conv_tr_wetscav = 0)
1961 !
1962 if (conv_tr_wetscav == 1) then
1965 tfac = (HL_t0 - t(i,k))/(t(i,k)*HL_t0)
1966 do nv=2,num_chem
1968 aq_gas_ratio = 0.0
1969 is_moz_chm: if( chemopt == MOZCART_KPP .or. chemopt == T1_MOZCART_KPP .or. &
1970 chemopt == MOZART_MOSAIC_4BIN_KPP .or. chemopt == MOZART_MOSAIC_4BIN_AQ_KPP ) then
1971 HLndx = HLC_ndx(nv)
1972 if( HLndx /= 0 ) then
1973 HLCnst1 = HLC(HLndx)%hcnst(1) ; HLCnst2 = HLC(HLndx)%hcnst(2)
1974 HLCnst3 = HLC(HLndx)%hcnst(3) ; HLCnst4 = HLC(HLndx)%hcnst(4)
1975 HLCnst5 = HLC(HLndx)%hcnst(5) ; HLCnst6 = HLC(HLndx)%hcnst(6)
1976 kh = HLCnst1 * exp( HLCnst2* tfac )
1977 if( HLCnst3 /= 0. ) then
1978 dk1s = HLCnst3 * exp( HLCnst4* tfac )
1979 else
1980 dk1s = 0.
1981 endif
1982 if( HLCnst5 /= 0. ) then
1983 dk2s = HLCnst5 * exp( HLCnst6* tfac )
1984 else
1985 dk2s = 0.
1986 endif
1987 if( nv /= p_nh3 ) then
1988 heff = kh*(1. + dk1s*hion_inv*(1. + dk2s*hion_inv))
1989 else
1990 heff = kh*(1. + dk1s*hion/dk2s)
1991 endif
1992 aq_gas_ratio = moz_aq_frac(t(i,k), clw_all(i,k)*dens, heff )
1993 endif
1994 else is_moz_chm
1996 ! Fraction of gas phase species that partions into the liquid phase:
1998 ! tried to be consistent with values and species in module_mozcart_wetscav.F
1999 if (nv .eq. p_h2o2) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 8.33e+04, 7379.)
2000 if (nv .eq. p_hno3) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 2.6e+06, 8700.)
2001 if (nv .eq. p_hcho) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 6.30e+03, 6425.)
2002 if (nv .eq. p_ch3ooh) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 3.11e+02, 5241.)
2003 if (nv .eq. p_c3h6ooh) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 2.20e+02, 5653.)
2004 if (nv .eq. p_paa) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 8.37e+02, 5308.)
2005 if (nv .eq. p_hno4) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 1.2e+04, 6900.) ! values from henrys-law.org, Regimbal and Mozurkewich, 1997
2006 if (nv .eq. p_onit) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 1.00e+03, 6000.)
2007 if (nv .eq. p_mvk) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 1.7e-03, 0.)
2008 if (nv .eq. p_macr) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 1.70e-03, 0.)
2009 if (nv .eq. p_etooh) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 3.36e+02, 5995.)
2010 if (nv .eq. p_prooh) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 3.36e+02, 5995.)
2011 if (nv .eq. p_acetp) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 3.36e+02, 5995.)
2012 if (nv .eq. p_mgly) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 3.71e+03, 7541.)
2013 if (nv .eq. p_mvkooh) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 2.6e+06, 8700.)
2014 if (nv .eq. p_onitr) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 7.51e+03, 6485.)
2015 if (nv .eq. p_isooh) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 2.6e+06, 8700.)
2016 if (nv .eq. p_ch3oh) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 2.20e+02, 4934.)
2017 if (nv .eq. p_c2h5oh) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 2.00e+02, 6500.)
2018 if (nv .eq. p_glyald) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 4.14e+04, 4630.)
2019 if (nv .eq. p_hydrald) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 7.00e+01, 6000.)
2020 if (nv .eq. p_ald) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 1.14e+01, 6267.)
2021 if (nv .eq. p_isopn) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 1.00e+01, 0.)
2022 if (nv .eq. p_alkooh) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 3.11e+02, 5241.)
2023 if (nv .eq. p_mekooh) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 3.11e+02, 5241.)
2024 if (nv .eq. p_tolooh) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 3.11e+02, 5241.)
2025 if (nv .eq. p_terpooh) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 3.11e+02, 5241.)
2026 if (nv .eq. p_nh3) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 7.40e+01, 3400.)
2027 if (nv .eq. p_xooh) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 90.5, 5607.)
2028 if (nv .eq. p_ch3cooh) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 4.1e3, 6300.)
2029 if (nv .eq. p_so2) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 1.2, 3100.)
2030 if (nv .eq. p_sulf) aq_gas_ratio = aq_frac(p(i,k)*100., t(i,k), clw_all(i,k)*dens, 1e+11, 0.) ! order of magnitude approx. (Gmitro and Vermeulen, 1964)
2031 ENDIF is_moz_chm
2033 ! if (nv.eq.p_so2) aq_gas_ratio = 1.0
2034 ! if (nv.eq.p_sulf) aq_gas_ratio = 1.0
2035 ! if (nv.eq.p_nh3) aq_gas_ratio = 1.0
2036 ! if (nv.eq.p_hno3) aq_gas_ratio = 1.0
2039 if (nv.gt.numgas) aq_gas_ratio = 0.5
2041 if (nv.eq.p_so4aj) aq_gas_ratio = 1.0
2042 if (nv.eq.p_nh4aj) aq_gas_ratio = 1.0
2043 if (nv.eq.p_no3aj) aq_gas_ratio = 1.0
2045 if (nv.eq.p_bc1 .or. nv.eq.p_oc1 .or. nv.eq.p_dms) aq_gas_ratio=0.
2046 if (nv.eq.p_sulf .or. nv.eq.p_seas_1 .or. nv.eq.p_seas_2) aq_gas_ratio=1.
2047 if (nv.eq.p_seas_3 .or. nv.eq.p_seas_4) aq_gas_ratio=1.
2048 if (nv.eq.p_bc2 .or. nv.eq.p_oc2) aq_gas_ratio=0.8
2050 if (aq_gas_ratio > 0.0) then
2051 tr_c(i,k,nv) = aq_gas_ratio*tr_up(i,k,nv) ! Amount of species cloud + rain water
2052 trch = tr_up(i,k,nv)-tr_c(i,k,nv) ! Amount of species remaining in gas phase
2053 trcc = (tr_up(i,k,nv)-trch)/(1.+c0*dz*zu(i,k)) ! Amount of species cloud + rain water
2054 tr_pw(i,k,nv) = c0*dz*trcc*zu(i,k) ! Amount of species in rain water
2055 tr_up(i,k,nv) = trcc+trch ! Total amount of species in updraft = conc in liq water (trcc) + conc in gas phase (trch)
2056 endif
2058 enddo
2060 endif ! parameterization wetscavenging option
2063 100 CONTINUE
2065 END subroutine cup_up_tracer
2067 ! calculates the fraction (0-1) of a soluble gas that should
2068 ! partition into the liquid phase according to instantaneous
2069 ! Henry's law equilibrium
2070 REAL FUNCTION aq_frac(p, T, q, Kh_298, dHoR)
2072 REAL, INTENT(IN) :: p, & ! air pressure (Pa)
2073 T, & ! air temperature (K)
2074 q, & ! total liquid water content (kg/m3)
2075 Kh_298, & ! Henry's law constant (M/atm == (mol_aq/dm3_aq)/atm)
2076 dHoR ! enthalpy of solution (in K, dH/R)
2078 REAL, PARAMETER :: Rgas = 8.31446 ! ideal gas constant (J mol-1 K-1)
2080 ! local variables
2081 REAL :: Kh, tr_air, tr_aq
2083 ! with van't Hoff's equation as temperature dependence
2084 ! and conversion to SI units ( (mol_aq/m3_aq)/Pa )
2085 Kh = Kh_298 * exp ( dHoR * ( 1.0/T - 1.0/298 ) ) * 101.325
2087 ! moles tracer m-3_air
2088 tr_air = 1 / (Rgas * T)
2090 ! moles tracer m-3 (air)
2091 tr_aq = Kh * (q / 1000.0)
2093 aq_frac = min( 1.0, max( 0.0, tr_aq / (tr_aq + tr_air) ) )
2095 END FUNCTION aq_frac
2097 REAL FUNCTION moz_aq_frac(T, q, heff )
2099 REAL, INTENT(IN) :: T, & ! air temperature (K)
2100 q, & ! total liquid water content (kg/m3)
2101 heff ! Henry's law constant (M/atm == (mol_aq/dm3_aq)/atm)
2103 REAL, PARAMETER :: Rgas = 8.31446 ! ideal gas constant (J mol-1 K-1)
2105 ! local variables
2106 REAL :: tr_air, tr_aq
2108 ! moles tracer m-3_air
2109 tr_air = 1. / (Rgas * T)
2111 ! moles tracer m-3 (air)
2112 tr_aq = heff * (q / 1000.0)
2114 moz_aq_frac = min( 1.0, max( 0.0, tr_aq / (tr_aq + tr_air) ) )
2116 END FUNCTION moz_aq_frac
2120 SUBROUTINE cup_dd_tracer(ierr,z_cup,qrcd,tracer,tre_cup,tr_up,tr_dd, &
2121 tr_pw,tr_pwd,jmin,cdd,entr,zd,pwdper,wetdep,xmb,k22, &
2122 numgas,num_chem,ids,ide, jds,jde, kds,kde, &
2123 ims,ime, jms,jme, kms,kme, &
2124 its,ite, jts,jte, kts,kte)
2126 USE module_model_constants, only: mwdry
2128 implicit none
2129 !
2130 ! on input
2131 !
2133 ! only local wrf dimensions are need as of now in this routine
2135 integer &
2136 ,intent (in ) :: &
2137 numgas,num_chem, &
2138 ids,ide, jds,jde, kds,kde, &
2139 ims,ime, jms,jme, kms,kme, &
2140 its,ite, jts,jte, kts,kte
2141 real, dimension (its:ite,kts:kte) &
2142 ,intent (in ) :: &
2143 pwdper,zd,cdd,qrcd,z_cup
2144 real, dimension (its:ite) &
2145 ,intent (in ) :: &
2146 xmb
2147 real, dimension (its:ite,kts:kte,1:num_chem) &
2148 ,intent (inout ) :: &
2149 tr_dd,tr_pwd,tr_up
2150 real, dimension (its:ite,kts:kte,1:num_chem) &
2151 ,intent (in ) :: &
2152 tre_cup,tracer,tr_pw
2153 real, dimension (its:ite,1:num_chem) &
2154 ,intent (out ) :: &
2155 wetdep
2157 ! entr= entrainment rate
2158 real &
2159 ,intent (in ) :: &
2160 entr
2161 integer, dimension (its:ite) &
2162 ,intent (in ) :: &
2163 jmin
2164 ! ierr error value, maybe modified in this routine
2166 integer, dimension (its:ite) &
2167 ,intent (inout) :: &
2168 ierr,k22
2169 ! local variables in this routine
2170 !
2172 integer :: &
2173 iall,i,k,nv,ki,kj
2174 real :: &
2175 dh,c0,dz,radius
2176 integer :: &
2177 itf,ktf
2179 real :: &
2180 evaporate,condensate
2182 itf=MIN(ite,ide-1)
2183 ktf=MIN(kte,kde-1)
2185 ! Initialize the tracer amount that evaporated from rain water:
2186 tr_pwd(its:ite,kts:kte,1:num_chem) = 0.0
2188 ! Initialize wet deposition:
2189 wetdep(its:ite,1:num_chem) = 0.0
2191 ! Calculate wet deposition with re-evaporation (based on wet scavenging in cup_up_tracer);
2192 ! assume no gas takeup by rain during fall for now
2194 do i=its,itf
2196 IF(ierr(I).eq.0)then
2198 ! why shouldn't that work for aerosols as well?
2199 ! do nv=1,numgas ! Only gas phase species evaporate along with rain water
2200 do nv=1,num_chem
2202 do k=ktf,kts,-1 ! Descending loop over all levels
2204 ! Start with initializing the condensate with the tracer amount in rain water on this level:
2205 condensate = tr_pw(i,k,nv)
2207 do kj=k,kts,-1 ! Descending loop over this and lower levels
2208 ! Evaporated tracer amount at the current level:
2209 evaporate = condensate*pwdper(i,kj)
2210 ! Accumulate the evaporated tracer amount:
2211 tr_pwd(i,kj,nv) = tr_pwd(i,kj,nv) + evaporate
2212 ! Remove the evaporated tracer amount from the condensate:
2213 condensate = max(0.0,condensate - evaporate)
2214 enddo
2216 enddo
2218 ! Calculate the wet deposition as the column integral over the
2219 ! tracer amount in rain water - evaporated tracer amount:
2221 do k=kts,ktf
2222 wetdep(i,nv) = wetdep(i,nv) + tr_pw(i,k,nv) - tr_pwd(i,k,nv)
2223 enddo
2225 wetdep(i,nv) = max(0.0,wetdep(i,nv))
2227 enddo
2229 endif
2231 enddo
2233 ! In downdraft, do only transport of tracers
2235 do i=its,itf
2237 IF(ierr(I).eq.0)then
2239 do nv=1,num_chem
2240 tr_dd(i,jmin(i):ktf,nv)=tre_cup(i,jmin(i):ktf,nv) ! Tracer amount in downdraft
2241 enddo
2243 do ki=jmin(i)-1,1,-1
2244 DZ=Z_cup(i,ki+1)-Z_cup(i,ki)
2245 do nv=1,num_chem
2246 tr_dd(i,ki,nv)=(tr_dd(i,ki+1,nv)*(1.-.5*CDD(i,ki)*DZ) &
2247 +entr*DZ*tracer(i,ki,nv) &
2248 )/(1.+entr*DZ-.5*CDD(i,ki)*DZ)
2249 enddo
2251 enddo
2253 endif
2255 enddo
2257 ! Add evaporation from rain water:
2259 do i=its,itf
2260 IF(ierr(I).eq.0)then
2261 do nv=1,num_chem
2262 do k=kts,ktf
2263 tr_dd(i,k,nv) = tr_dd(i,k,nv) + tr_pwd(i,k,nv)
2264 enddo
2265 enddo
2266 endif
2267 enddo
2269 ! To obtain wet deposition, multiply with base mass flux; units are given
2270 ! assuming that the base mass flux units [xmb] = kg(air)/m2/s:
2272 do nv=1,num_chem
2273 if (nv <= numgas) then
2274 do i=its,itf
2275 if(ierr(I).eq.0)then
2276 wetdep(i,nv)=wetdep(i,nv)*xmb(i)/mwdry ! mmol/m2/s for gas phase species
2277 wetdep(i,nv) = max(0.0,wetdep(i,nv))
2278 endif
2279 enddo
2280 else
2281 do i=its,itf
2282 if(ierr(I).eq.0)then
2283 wetdep(i,nv)=wetdep(i,nv)*xmb(i) ! ug/m2/s for aerosol mass, #/m2/s for aerosol number
2284 wetdep(i,nv) = max(0.0,wetdep(i,nv))
2285 endif
2286 enddo
2287 endif
2288 enddo
2290 END subroutine cup_dd_tracer
2293 SUBROUTINE neg_check_ct(name,pret,ktop,epsilc,dt,q,outq,iopt,num_chem, &
2294 its,ite,kts,kte,itf,ktf,ipr,jpr,npr,j)
2296 INTEGER, INTENT(IN ) :: iopt,num_chem,its,ite,kts,kte,itf,ktf,ipr,jpr,npr,j
2298 real, dimension (its:ite,kts:kte,num_chem ) , &
2299 intent(inout ) :: &
2300 q,outq
2301 real, dimension (its:ite ) , &
2302 intent(in ) :: &
2303 pret
2304 integer, dimension (its:ite ) , &
2305 intent(in ) :: &
2306 ktop
2307 real &
2308 ,intent (in ) :: &
2309 dt,epsilc
2310 real :: tracermin,tracermax,thresh,qmem,qmemf,qmem2,qtest,qmem1
2311 character *(*) name
2312 integer :: nv, i, k
2313 !
2314 ! check whether routine produces negative q's. This can happen, since
2315 ! tendencies are calculated based on forced q's. This should have no
2316 ! influence on conservation properties, it scales linear through all
2317 ! tendencies. Use iopt=0 to test for each tracer seperately, iopt=1
2318 ! for a more severe limitation...
2319 !
2320 thresh=epsilc
2321 ! thresh=1.e-30
2322 if(iopt.eq.0)then
2323 do nv=2,num_chem
2324 do 100 i=its,itf
2325 if(pret(i).le.0.and.name.eq.'deep')go to 100
2326 tracermin=q(i,kts,nv)
2327 tracermax=q(i,kts,nv)
2328 do k=kts+1,kte-1
2329 tracermin=min(tracermin,q(i,k,nv))
2330 tracermax=max(tracermax,q(i,k,nv))
2331 enddo
2332 tracermin=max(tracermin,thresh)
2333 qmemf=1.
2334 !
2335 ! first check for minimum restriction
2336 !
2337 do k=kts,ktop(i)
2338 !
2339 ! tracer tendency
2340 !
2341 qmem=outq(i,k,nv)
2342 !
2343 ! only necessary if there is a tendency
2344 !
2345 if(qmem.lt.0.)then
2346 qtest=q(i,k,nv)+outq(i,k,nv)*dt
2347 if(qtest.lt.tracermin)then
2348 !
2349 ! qmem2 would be the maximum allowable tendency
2350 !
2351 qmem1=outq(i,k,nv)
2352 qmem2=(tracermin-q(i,k,nv))/dt
2353 qmemf=min(qmemf,qmem2/qmem1)
2354 if(qmemf.gt.1.)print *,'something wrong in negct_1',qmem2,qmem1
2355 if(i.eq.ipr.and.j.eq.jpr.and.nv.eq.npr)then
2356 print *,k,qtest,qmem2,qmem1,qmemf
2357 endif
2358 qmemf=max(qmemf,0.)
2359 endif
2360 endif
2361 enddo
2362 do k=kts,ktop(i)
2363 outq(i,k,nv)=outq(i,k,nv)*qmemf
2364 enddo
2365 !
2366 ! now check max
2367 !
2368 qmemf=1.
2369 do k=kts,ktop(i)
2370 !
2371 ! tracer tendency
2372 !
2373 qmem=outq(i,k,nv)
2374 !
2375 ! only necessary if there is a tendency
2376 !
2377 if(qmem.gt.0.)then
2378 qtest=q(i,k,nv)+outq(i,k,nv)*dt
2379 if(qtest.gt.tracermax)then
2380 !
2381 ! qmem2 would be the maximum allowable tendency
2382 !
2383 qmem1=outq(i,k,nv)
2384 qmem2=(tracermax-q(i,k,nv))/dt
2385 qmemf=min(qmemf,qmem2/qmem1)
2386 if(qmemf.gt.1.)print *,'something wrong in negct_2',qmem2,qmem1
2387 if(i.eq.ipr.and.j.eq.jpr.and.nv.eq.npr)then
2388 print *,'2',k,qtest,qmem2,qmem1,qmemf
2389 endif
2390 qmemf=max(qmemf,0.)
2391 endif
2392 endif
2393 enddo
2394 do k=kts,ktop(i)
2395 outq(i,k,nv)=outq(i,k,nv)*qmemf
2396 enddo
2397 100 continue
2398 enddo
2399 !
2400 ! ELSE
2401 !
2402 elseif(iopt.eq.1)then
2403 do i=its,itf
2404 qmemf=1.
2405 do k=kts,ktop(i)
2406 do nv=2,num_chem
2407 !
2408 ! tracer tendency
2409 !
2410 qmem=outq(i,k,nv)
2411 !
2412 ! only necessary if tendency is larger than zero
2413 !
2414 if(qmem.lt.0.)then
2415 qtest=q(i,k,nv)+outq(i,k,nv)*dt
2416 if(qtest.lt.thresh)then
2417 !
2418 ! qmem2 would be the maximum allowable tendency
2419 !
2420 qmem1=outq(i,k,nv)
2421 qmem2=(thresh-q(i,k,nv))/dt
2422 qmemf=min(qmemf,qmem2/qmem1)
2423 qmemf=max(0.,qmemf)
2424 endif
2425 endif
2426 enddo
2427 enddo
2428 do nv=2,num_chem
2429 do k=kts,ktop(i)
2430 outq(i,k,nv)=outq(i,k,nv)*qmemf
2431 enddo
2432 enddo
2433 enddo
2434 endif
2436 END SUBROUTINE neg_check_ct
2438 SUBROUTINE conv_tr_wetscav_init( numgas, num_chem )
2440 use module_state_description, only : param_first_scalar
2441 use module_scalar_tables, only : chem_dname_table
2442 use module_chem_utilities, only : UPCASE
2443 use module_HLawConst
2445 integer, intent(in) :: numgas, num_chem
2447 !----------------------------------------------------------------------
2448 ! local variables
2449 !----------------------------------------------------------------------
2450 integer :: m, m1
2451 integer :: astat
2452 character(len=64) :: HL_tbl_name
2453 character(len=64) :: wrf_spc_name
2455 is_allocated : &
2456 if( .not. allocated(HLC_ndx) ) then
2457 !----------------------------------------------------------------------
2458 ! scan HLawConst table for match with chem_opt scheme gas phase species
2459 !----------------------------------------------------------------------
2460 allocate( HLC_ndx(num_chem),stat=astat )
2461 if( astat /= 0 ) then
2462 call wrf_error_fatal("conv_tr_wetscav_init: failed to allocate HLC_ndx")
2463 endif
2464 HLC_ndx(:) = 0
2465 do m = param_first_scalar,numgas
2466 wrf_spc_name = chem_dname_table(1,m)
2467 call upcase( wrf_spc_name )
2468 do m1 = 1,nHLC
2469 HL_tbl_name = HLC(m1)%name
2470 call upcase( HL_tbl_name )
2471 if( trim(HL_tbl_name) == trim(wrf_spc_name) ) then
2472 HLC_ndx(m) = m1
2473 exit
2474 endif
2475 end do
2476 end do
2477 endif is_allocated
2479 END SUBROUTINE conv_tr_wetscav_init
2481 !-------------------------------------------------------
2482 END MODULE module_ctrans_grell