| blob | 64d64c773959f974a05b722434e989c9727e70b6 |
1 MODULE module_cu_gf_ctrans
2 real, parameter::g=9.81
3 INTEGER, allocatable :: HLC_ndx(:)
4 !++ GF CTRAN ++ lyy 01/2020
5 #if ( WRF_CHEM == 1 )
6 contains
8 SUBROUTINE ctrans_gf(numgas,num_chem,tracer,chemopt,traceropt &
9 ,tracert,conv_tr_wetscav,conv_tr_aqchem &
10 ,po,po_cup,zo_cup &
11 ,zuo,zdo,pwo,pwdo,pwevo,pwavo &
12 ,up_massentro,up_massdetro &
13 ,dd_massentro,dd_massdetro &
14 ,tempco,clw_all &
15 ,ktop,k22,kbcon,jmin &
16 ,xmb,ierr,edto &
17 ,itf,ktf,its,ite,kts,kte &
18 ,ishallow)
20 IMPLICIT NONE
21 integer,intent (in ) :: &
22 itf,ktf,its,ite, kts,kte, &
23 numgas,num_chem,conv_tr_wetscav,conv_tr_aqchem, &
24 chemopt,traceropt,ishallow
25 integer,dimension (its:ite),intent (in ) :: &
26 kbcon,ktop,k22,ierr,jmin
27 real, dimension (its:ite),intent (in ) :: &
28 pwevo,pwavo,xmb,edto
29 real,dimension (its:ite,kts:kte),intent (in ):: &
30 po,po_cup,zo_cup,zuo,zdo,pwo,pwdo,tempco,clw_all, &
31 up_massentro,up_massdetro,dd_massentro,dd_massdetro
32 REAL,DIMENSION(its:ite , kts:kte , num_chem),INTENT(IN):: &
33 tracer
34 REAL,DIMENSION(its:ite , kts:kte , num_chem),INTENT(INOUT):: &
35 tracert
36 !local
37 INTEGER :: nv,i,k
38 real, dimension(its:ite, kts:kte, num_chem) :: &
39 tr_up,tr_dd,tr_pw,tot_up_pw,tre_cup,tr_pwd
40 real::dp
42 !-1) get mass mixing ratios at the cloud levels
43 call cup_env_clev_tr_gf(tracer,tre_cup,num_chem,ierr &
44 ,itf,ktf,its,ite,kts,kte)
45 !-2) determine in-cloud tracer mixing ratios
46 ! 2a) chem - updraft
47 call cup_up_tracer_gf(tracer,tre_cup,num_chem,numgas &
48 ,chemopt,traceropt &
49 ,conv_tr_wetscav,conv_tr_aqchem &
50 ,tr_up,tr_pw,tot_up_pw &
51 ,zo_cup,po,po_cup,clw_all,tempco &
52 ,zuo,up_massentro,up_massdetro &
53 ,k22,kbcon,ktop,ierr &
54 ,itf,ktf,its,ite,kts,kte)
55 ! 2b) chem - downdraft
56 if (ishallow==0) then
57 call cup_dd_tracer_gf(num_chem,tracer,tre_cup &
58 ,tot_up_pw &
59 ,tr_dd,tr_pwd,po_cup,pwdo &
60 ,pwevo,pwavo,edto,zdo &
61 ,dd_massentro,dd_massdetro &
62 ,jmin,ierr &
63 ,itf,ktf,its,ite,kts,kte)
64 endif
65 !-3) determine the vertical transport
66 do i=its,itf
67 if (ierr(i)==0) then
68 do k=kts,ktop(i)
69 dp=100.*(po_cup(i,k)-po_cup(i,k+1))
70 do nv=2,num_chem
71 tracert(i,k,nv)=-(zuo(i,k+1)*(tr_up(i,k+1,nv)-tre_cup(i,k+1,nv)) - &
72 zuo(i,k )*(tr_up(i,k ,nv)-tre_cup(i,k ,nv)))*g/dp*xmb(i) &
73 +(zdo(i,k+1)*(tr_dd(i,k+1,nv)-tre_cup(i,k+1,nv)) - &
74 zdo(i,k )*(tr_dd(i,k ,nv)-tre_cup(i,k ,nv)))*g/dp*edto(i)*xmb(i)
75 enddo ! nv
76 enddo ! k
77 endif ! ierr
78 enddo ! i
80 END SUBROUTINE ctrans_gf
82 SUBROUTINE cup_env_clev_tr_gf(tracer,tre_cup,num_chem,ierr &
83 ,itf,ktf,its,ite,kts,kte)
84 IMPLICIT NONE
85 integer,intent (in ) :: &
86 itf,ktf,its,ite, kts,kte,num_chem
87 integer,dimension (its:ite),intent (in ) :: &
88 ierr
89 REAL,DIMENSION(its:ite, kts:kte, num_chem),INTENT(IN):: &
90 tracer
91 real,dimension(its:ite, kts:kte, num_chem),INTENT(INOUT) :: &
92 tre_cup
93 !local
94 integer::i,k,nv
95 integer,parameter :: clev_opt=2 !-use option 2
97 if (clev_opt == 1) then
98 !original version
99 do i=its,itf
100 if (ierr(i).eq.0) then
101 do nv=2,num_chem
102 do k=kts+1,ktf
103 tre_cup(i,k,nv)=0.5*(tracer(i,k-1,nv)+tracer(i,k,nv))
104 enddo !k
105 tre_cup(i,kts,nv)=tracer(i,kts,nv)
106 enddo !nv
107 endif !ierr
108 enddo !i
109 else
110 ! version 2: tre_cup(k+1/2)=tracer(k) => smoother profiles
111 do i=its,itf
112 if (ierr(i).eq.0) then
113 do nv=2,num_chem
114 do k=kts,ktf
115 tre_cup(i,k,nv)=tracer(i,k,nv)
116 enddo !k
117 enddo !nv
118 endif !ierr
119 enddo !i
120 endif !clev_opt
122 END SUBROUTINE cup_env_clev_tr_gf
124 SUBROUTINE cup_up_tracer_gf(tracer,tre_cup,num_chem,numgas &
125 ,chemopt,traceropt &
126 ,conv_tr_wetscav,conv_tr_aqchem &
127 ,tr_up,tr_pw,tot_up_pw &
128 ,zo_cup,p,po_cup,clw_all,t &
129 ,zuo,up_massentro,up_massdetro &
130 ,k22,kbcon,ktop,ierr &
131 ,itf,ktf,its,ite,kts,kte)
132 IMPLICIT NONE
133 integer,intent (in ) :: &
134 itf,ktf,its,ite, kts,kte, &
135 numgas,num_chem,conv_tr_wetscav,conv_tr_aqchem, &
136 chemopt,traceropt
137 integer,dimension (its:ite),intent (in ) :: &
138 kbcon,ktop,k22,ierr
139 real,dimension (its:ite,kts:kte),intent (in ):: &
140 p,po_cup,zo_cup,zuo,t,clw_all, &
141 up_massentro,up_massdetro
142 REAL,DIMENSION(its:ite , kts:kte , num_chem),INTENT(IN):: &
143 tracer,tre_cup
144 REAL,DIMENSION(its:ite , kts:kte , num_chem),INTENT(OUT):: &
145 tr_up,tr_pw
146 REAL,DIMENSION(its:ite , num_chem),INTENT(OUT):: &
147 tot_up_pw
148 !local
149 INTEGER :: nv,i,k
150 REAL,DIMENSION(its:ite,num_chem):: &
151 tc_b
152 REAL ::XZZ,XZD,XZE,denom,dz,dp
153 ! initialize
154 do i=its,itf
155 if (ierr(i)==0) then
156 do k=kts,ktf
157 do nv=2,num_chem
158 tr_up(i,k,nv)=tre_cup(i,k,nv)
159 tr_pw(i,k,nv)=0.
160 tot_up_pw(i,nv)=0.
161 enddo ! nv
162 enddo ! k
163 endif ! ierr
164 enddo ! i
166 ! below k22
167 do i=its,itf
168 if (ierr(i)==0) then
169 do nv=2,num_chem
170 call get_cloud_bc_chem(kte,tre_cup(i,1:kte,nv),tc_b(i,nv),k22(i))
171 do k=kts,k22(i)
172 tr_up(i,k,nv)=tc_b(i,nv)
173 enddo ! k
174 enddo ! nv
175 endif ! ierr
176 enddo ! i
178 ! above k22
179 DO i=its,itf
180 if (ierr(i)==0) then
181 do k=k22(i)+1,ktop(i)+1
182 dz=zo_cup(i,k)-zo_cup(i,k-1)
183 dp=100.*(po_cup(i,k)-po_cup(i,k+1))
184 !--entr, detr, mass flux--
185 XZZ=zuo(i,k-1)
186 XZD=0.5*up_massdetro(i,k-1)
187 XZE=up_massentro(i,k-1)
188 denom=(XZZ-XZD+XZE)
189 !--transport + mixing
190 do nv=2,num_chem
191 if (denom>0) then
192 tr_up(i,k,nv)=(tr_up(i,k-1,nv)*XZZ-tr_up(i,k-1,nv)*XZD &
193 +tracer(i,k-1,nv)*XZE)/denom
194 else
195 tr_up(i,k,nv)=tr_up(i,k-1,nv)
196 endif ! denom
197 enddo ! nv
198 !--Aqueous Chemistry--
199 if ((conv_tr_aqchem==1).and.(chemopt>0)) then
200 call aqchem_gf(chemopt,num_chem,p(i,k),t(i,k) &
201 ,dz,tr_up(i,k,:),clw_all(i,k) &
202 )
203 endif
204 !--wet scavenging--
205 if ((conv_tr_wetscav==1).and.(chemopt>0)) then
206 do nv=2,num_chem
207 call wetscav(tr_up(i,k,nv),tr_pw(i,k,nv) &
208 ,zuo(i,k),nv,p(i,k),t(i,k),clw_all(i,k),dz &
209 ,chemopt,numgas,num_chem)
210 tot_up_pw(i,nv)=tot_up_pw(i,nv)+tr_pw(i,k,nv)*dp/g
211 enddo ! nv
212 endif ! scav
213 enddo ! k
214 endif ! ierr=0
215 ENDDO ! i
216 END SUBROUTINE cup_up_tracer_gf
218 SUBROUTINE cup_dd_tracer_gf(num_chem,tracer,tre_cup &
219 ,tot_up_pw &
220 ,tr_dd,tr_pwd,po_cup,pwdo &
221 ,pwevo,pwavo,edto,zdo &
222 ,dd_massentro,dd_massdetro &
223 ,jmin,ierr &
224 ,itf,ktf,its,ite,kts,kte)
225 IMPLICIT NONE
226 integer,intent (in ) :: &
227 itf,ktf,its,ite, kts,kte, &
228 num_chem
229 integer,dimension (its:ite),intent (in ) :: &
230 ierr,jmin
231 real, dimension (its:ite),intent (in ) :: &
232 pwevo,pwavo,edto
233 real,dimension (its:ite,kts:kte),intent (in ):: &
234 po_cup,zdo,pwdo,dd_massentro,dd_massdetro
235 REAL,DIMENSION(its:ite , kts:kte , num_chem),INTENT(IN):: &
236 tracer,tre_cup
237 REAL,DIMENSION(its:ite, num_chem),INTENT(IN):: &
238 tot_up_pw
239 REAL,DIMENSION(its:ite , kts:kte , num_chem),INTENT(OUT):: &
240 tr_pwd,tr_dd
241 !local
242 INTEGER :: nv,i,k
243 real:: frac_evap,dp,XZZ,XZD,XZE,denom,pwdper
244 do i=its,itf
245 if (ierr(i)==0) then
246 do k=kts,ktf
247 do nv=2,num_chem
248 tr_dd(i,k,nv)=0.
249 tr_pwd(i,k,nv)=0.
250 enddo ! nv
251 enddo ! k
252 endif ! ierr
253 enddo ! i
255 do i=its,itf
256 if (ierr(i)==0) then
257 !--- fration of the total rain that was evaporated
258 frac_evap = - pwevo(i)/(1.e-16+pwavo(i))
259 !--- scalar concentration in-cloud - downdraft
261 !--- at k=jmim
262 k=jmin(i)
263 pwdper = pwdo(i,k)/(1.e-16+pwevo(i)) *frac_evap ! > 0
264 dp= 100.*(po_cup(i,k)-po_cup(i,k+1))
265 do nv=2,num_chem
266 !--downdrafts will be initiate with a mixture of 50% environmental and in-cloud concentrations
267 tr_dd(i,k,nv)=tre_cup(i,k,nv)
268 tr_pwd(i,k,nv)=-pwdper*tot_up_pw(i,nv)*g/dp
269 tr_dd(i,k,nv)=tr_dd(i,k,nv)-tr_pwd(i,k,nv)
270 enddo ! nv
271 !
272 !--- calculate downdraft mass terms
273 do k=jmin(i)-1,kts,-1
274 XZZ= zdo(i,k+1)
275 XZD= 0.5*dd_massdetro(i,k )
276 XZE= dd_massentro(i,k )
277 denom = (XZZ-XZD+XZE)
278 !-- transport + mixing
279 do nv=2,num_chem
280 if(denom > 0.) then
281 tr_dd(i,k,nv) = (tr_dd(i,k+1,nv)*XZZ-tr_dd(i,k+1,nv)*XZD &
282 +tracer(i,k,nv)*XZE) / denom
283 else
284 tr_dd(i,k,nv) = tr_dd(i,k+1,nv)
285 endif
286 enddo
287 !-- evaporation term
288 dp= 100.*(po_cup(i,k)-po_cup(i,k+1))
289 !-- fraction of evaporated precip per layer
290 pwdper= pwdo(i,k)/(1.e-16+pwevo(i))! > 0
291 !-- fraction of the total precip that was actually evaporated at layer k
292 pwdper= pwdper*frac_evap
293 !-- sanity check
294 pwdper= min(1.,max(pwdper,0.))
295 do nv=2,num_chem
296 !-- amount evaporated by the downdraft from the precipitation
297 tr_pwd(i,k,nv)=-pwdper* tot_up_pw(i,nv)*g/dp
298 ! < 0. => source term for the downdraft tracer concentration
299 !-- final tracer in the downdraft
300 tr_dd(i,k,nv)= tr_dd(i,k,nv)-tr_pwd(i,k,nv) ! observe that -tr_pwd is > 0.
301 enddo ! nv
302 enddo ! k
303 endif ! ierr
304 enddo ! i
306 END SUBROUTINE cup_dd_tracer_gf
309 SUBROUTINE wetscav(tr_up1d,tr_pw1d &
310 ,zu1d,nv,p1d,t1d,clw_all1d,dz &
311 ,chemopt,numgas,num_chem)
312 USE module_HLawConst, only:HLC
313 USE module_state_description, ONLY: mozart_mosaic_4bin_kpp, &
314 mozart_mosaic_4bin_aq_kpp, &
315 mozcart_kpp, t1_mozcart_kpp, &
316 p_nh3,p_h2o2,p_hno3,p_hcho,p_ch3ooh, &
317 p_c3h6ooh,p_paa,p_hno4,p_onit,p_mvk, &
318 p_macr,p_etooh,p_prooh,p_acetp,p_mgly, &
319 p_mvkooh,p_onitr,p_isooh,p_ch3oh,p_c2h5oh, &
320 p_glyald,p_hydrald,p_ald,p_isopn,p_alkooh, &
321 p_mekooh,p_tolooh,p_terpooh,p_nh3,p_xooh, &
322 p_ch3cooh,p_so2,p_sulf,p_so4aj,p_nh4aj, &
323 p_no3aj,p_bc1,p_oc1,p_dms,p_sulf,p_seas_1, &
324 p_seas_2,p_seas_3,p_seas_4,p_bc2,p_oc2, &
325 p_hcooh
326 IMPLICIT NONE
327 integer,intent (in) :: &
328 chemopt,nv,numgas,num_chem
329 real,intent (in ):: &
330 p1d,t1d,clw_all1d,dz,zu1d
331 REAL,INTENT(INOUT):: &
332 tr_up1d,tr_pw1d
333 !local
334 real::tr_c,trch,trcc,c0,dens,tfac, &
335 aq_gas_ratio,kh,dk1s,dk2s, &
336 HLCnst1,HLCnst2,HLCnst3, &
337 HLCnst4,HLCnst5,HLCnst6, &
338 heff
339 integer::HLndx
340 real, parameter :: hion = 1.e-5
341 real, parameter :: hion_inv = 1./hion
342 real, parameter :: HL_t0 = 298.
343 REAL, PARAMETER :: mwdry = 28.966 ! Molecular mass of dry air (g/mol)
344 !++ ice retention (Li et al. 2019 JGR)
345 integer,parameter :: USE_ICE_RETENTION=1
346 real::reteff
347 !--
348 tfac=(HL_t0-t1d)/(t1d*HL_t0)
349 aq_gas_ratio=0.0
350 dens=0.1*p1d/t1d*mwdry/8.314472 !kg/m3
351 !++ ice retention
352 reteff = 0.
353 if( nv == p_h2o2 ) then
354 reteff=.64
355 elseif( nv == p_hno3 ) then
356 reteff=1.
357 elseif( nv == p_hcooh ) then
358 reteff=.64
359 elseif( nv == p_ch3ooh ) then
360 reteff=.02
361 elseif( nv == p_so2 ) then
362 reteff= .02
363 elseif( nv == p_hcooh ) then
364 reteff= .68
365 endif
366 c0=0.004
367 if (t1d < 273.15) c0=c0*exp(0.07*(t1d-273.15))
368 !--
369 !chem moz
370 if( chemopt == MOZCART_KPP .or. chemopt == T1_MOZCART_KPP .or. &
371 chemopt == MOZART_MOSAIC_4BIN_KPP .or. &
372 chemopt == MOZART_MOSAIC_4BIN_AQ_KPP ) then
373 ! This setup is not ideal, as the the sub is just a duplicate version of the
374 ! init that occurs on the /chem side. But, it will work for now.
375 if ( .not. allocated(HLC_ndx) ) then
376 call conv_tr_wetscav_init_phys( numgas, num_chem )
377 endif
378 HLndx = HLC_ndx(nv)
380 if( HLndx /= 0 ) then
381 HLCnst1 = HLC(HLndx)%hcnst(1)
382 HLCnst2 = HLC(HLndx)%hcnst(2)
383 HLCnst3 = HLC(HLndx)%hcnst(3)
384 HLCnst4 = HLC(HLndx)%hcnst(4)
385 HLCnst5 = HLC(HLndx)%hcnst(5)
386 HLCnst6 = HLC(HLndx)%hcnst(6)
387 kh = HLCnst1 * exp( HLCnst2* tfac )
388 if( HLCnst3 /= 0. ) then
389 dk1s = HLCnst3 * exp( HLCnst4* tfac )
390 else
391 dk1s = 0.
392 endif
393 if( HLCnst5 /= 0. ) then
394 dk2s = HLCnst5 * exp( HLCnst6* tfac )
395 else
396 dk2s = 0.
397 endif
398 if( nv /= p_nh3 ) then
399 heff = kh*(1. + dk1s*hion_inv*(1. + dk2s*hion_inv))
400 else
401 heff = kh*(1. + dk1s*hion/dk2s)
402 endif
403 aq_gas_ratio = moz_aq_frac(t1d, clw_all1d*dens, heff)
404 endif ! HLndx
405 else !chem moz
406 ! Fraction of gas phase species that partions into the liquid phase:
407 ! tried to be consistent with values and species in module_mozcart_wetscav.F
408 if (nv .eq. p_h2o2) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 8.33e+04, 7379.)
409 if (nv .eq. p_hno3) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 2.6e+06, 8700.)
410 if (nv .eq. p_hcho) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 6.30e+03, 6425.)
411 if (nv .eq. p_ch3ooh) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 3.11e+02, 5241.)
412 if (nv .eq. p_c3h6ooh) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 2.20e+02, 5653.)
413 if (nv .eq. p_paa) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 8.37e+02, 5308.)
414 if (nv .eq. p_hno4) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 1.2e+04, 6900.) ! values from henrys-law.org, Regimbal and Mozurkewich, 1997
415 if (nv .eq. p_onit) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 1.00e+03, 6000.)
416 if (nv .eq. p_mvk) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 1.7e-03, 0.)
417 if (nv .eq. p_macr) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 1.70e-03, 0.)
418 if (nv .eq. p_etooh) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 3.36e+02, 5995.)
419 if (nv .eq. p_prooh) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 3.36e+02, 5995.)
420 if (nv .eq. p_acetp) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 3.36e+02, 5995.)
421 if (nv .eq. p_mgly) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 3.71e+03, 7541.)
422 if (nv .eq. p_mvkooh) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 2.6e+06, 8700.)
423 if (nv .eq. p_onitr) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 7.51e+03, 6485.)
424 if (nv .eq. p_isooh) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 2.6e+06, 8700.)
425 if (nv .eq. p_ch3oh) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 2.20e+02, 4934.)
426 if (nv .eq. p_c2h5oh) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 2.00e+02, 6500.)
427 if (nv .eq. p_glyald) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 4.14e+04, 4630.)
428 if (nv .eq. p_hydrald) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 7.00e+01, 6000.)
429 if (nv .eq. p_ald) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 1.14e+01, 6267.)
430 if (nv .eq. p_isopn) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 1.00e+01, 0.)
431 if (nv .eq. p_alkooh) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 3.11e+02, 5241.)
432 if (nv .eq. p_mekooh) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 3.11e+02, 5241.)
433 if (nv .eq. p_tolooh) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 3.11e+02, 5241.)
434 if (nv .eq. p_terpooh) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 3.11e+02, 5241.)
435 if (nv .eq. p_nh3) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 7.40e+01, 3400.)
436 if (nv .eq. p_xooh) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 90.5, 5607.)
437 if (nv .eq. p_ch3cooh) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 4.1e3, 6300.)
438 if (nv .eq. p_so2) aq_gas_ratio = aq_frac(p1d*100., t1d, clw_all1d*dens, 1.2, 3100.)
439 endif !chem moz
440 !++ ice retention
441 if ((USE_ICE_RETENTION==1).and.(t1d < 273.15)) aq_gas_ratio=reteff*aq_gas_ratio
442 !--
443 ! if (nv.eq.p_so2) aq_gas_ratio = 1.0
444 ! if (nv.eq.p_sulf) aq_gas_ratio = 1.0
445 ! if (nv.eq.p_nh3) aq_gas_ratio = 1.0
446 ! if (nv.eq.p_hno3) aq_gas_ratio = 1.0
447 if (nv.gt.numgas) aq_gas_ratio = 0.5
448 if (nv.eq.p_so4aj) aq_gas_ratio = 1.0
449 if (nv.eq.p_nh4aj) aq_gas_ratio = 1.0
450 if (nv.eq.p_no3aj) aq_gas_ratio = 1.0
451 if (nv.eq.p_bc1 .or. nv.eq.p_oc1 .or. nv.eq.p_dms) aq_gas_ratio=0.
452 if (nv.eq.p_sulf .or. nv.eq.p_seas_1 .or. nv.eq.p_seas_2) aq_gas_ratio=1.
453 if (nv.eq.p_seas_3 .or. nv.eq.p_seas_4) aq_gas_ratio=1.
454 if (nv.eq.p_bc2 .or. nv.eq.p_oc2) aq_gas_ratio=0.8
456 if (aq_gas_ratio > 0.0) then
457 tr_c = aq_gas_ratio*tr_up1d ! Amount of species cloud + rain water
458 trch = tr_up1d-tr_c ! Amount of species remaining in gas phase
459 trcc = tr_c/(1.+c0*dz*zu1d) ! Amount of species cloud
460 tr_pw1d = c0*dz*trcc*zu1d ! Amount of species in rain water
461 tr_up1d = trcc+trch ! Total amount of species in updraft = conc in liq water (trcc) + conc in gas phase (trch)
462 endif
463 END SUBROUTINE wetscav
465 SUBROUTINE conv_tr_wetscav_init_phys( numgas, num_chem )
467 use module_state_description, only : param_first_scalar
468 use module_scalar_tables, only : chem_dname_table
469 use module_chem_utilities, only : UPCASE
470 use module_HLawConst
472 integer, intent(in) :: numgas, num_chem
474 !----------------------------------------------------------------------
475 ! local variables
476 !----------------------------------------------------------------------
477 integer :: m, m1
478 integer :: astat
479 character(len=64) :: HL_tbl_name
480 character(len=64) :: wrf_spc_name
482 is_allocated : &
483 if( .not. allocated(HLC_ndx) ) then
484 !----------------------------------------------------------------------
485 ! scan HLawConst table for match with chem_opt scheme gas phase species
486 !----------------------------------------------------------------------
487 allocate( HLC_ndx(num_chem),stat=astat )
488 if( astat /= 0 ) then
489 call wrf_error_fatal("conv_tr_wetscav_init: failed to allocate HLC_ndx")
490 endif
491 HLC_ndx(:) = 0
492 do m = param_first_scalar,numgas
493 wrf_spc_name = chem_dname_table(1,m)
494 call upcase( wrf_spc_name )
495 do m1 = 1,nHLC
496 HL_tbl_name = HLC(m1)%name
497 call upcase( HL_tbl_name )
498 if( trim(HL_tbl_name) == trim(wrf_spc_name) ) then
499 HLC_ndx(m) = m1
500 exit
501 endif
502 end do
503 end do
504 endif is_allocated
506 END SUBROUTINE conv_tr_wetscav_init_phys
511 REAL FUNCTION moz_aq_frac(T, q, heff )
512 implicit none
513 REAL, INTENT(IN) :: T, & ! air temperature (K)
514 q, & ! total liquid water content (kg/m3)
515 heff ! Henry's law constant (M/atm == (mol_aq/dm3_aq)/atm)
516 REAL, PARAMETER :: Rgas = 8.31446 ! ideal gas constant (J mol-1 K-1)
518 ! local variables
519 REAL :: tr_air, tr_aq
520 ! moles tracer m-3_air
521 tr_air = 1. / (Rgas * T)
522 ! moles tracer m-3 (air)
523 tr_aq = heff * (q / 1000.0)
524 moz_aq_frac = min( 1.0, max( 0.0, tr_aq / (tr_aq + tr_air) ) )
525 END FUNCTION moz_aq_frac
526 REAL FUNCTION aq_frac(p, T, q, Kh_298, dHoR)
527 REAL, INTENT(IN) :: p, & ! air pressure (Pa)
528 T, & ! air temperature (K)
529 q, & ! total liquid water content (kg/m3)
530 Kh_298, & ! Henry's law constant (M/atm == (mol_aq/dm3_aq)/atm)
531 dHoR ! enthalpy of solution (in K, dH/R)
532 REAL, PARAMETER :: Rgas = 8.31446 ! ideal gas constant (J mol-1 K-1)
533 ! local variables
534 REAL :: Kh, tr_air, tr_aq
535 ! with van't Hoff's equation as temperature dependence
536 ! and conversion to SI units ( (mol_aq/m3_aq)/Pa )
537 Kh = Kh_298 * exp ( dHoR * ( 1.0/T - 1.0/298 ) ) * 101.325
538 ! moles tracer m-3_air
539 tr_air = 1 / (Rgas * T)
540 ! moles tracer m-3 (air)
541 tr_aq = Kh * (q / 1000.0)
542 aq_frac = min( 1.0, max( 0.0, tr_aq / (tr_aq + tr_air) ) )
543 END FUNCTION aq_frac
544 SUBROUTINE aqchem_gf(chemopt,num_chem,p1d,t1d &
545 ,dz,tr_up1d,clw_all1d &
546 )
547 USE module_ctrans_aqchem
548 USE module_state_description, only:RADM2SORG,RADM2SORG_AQ, &
549 RACMSORG_AQ,RACMSORG_KPP,RADM2SORG_KPP, &
550 RACM_ESRLSORG_KPP,RACM_SOA_VBS_KPP, &
551 RADM2SORG_AQCHEM,RACMSORG_AQCHEM_KPP, &
552 CB05_SORG_VBS_AQ_KPP,RACM_SOA_VBS_HET_KPP, &
553 RACM_ESRLSORG_AQCHEM_KPP,RACM_SOA_VBS_AQCHEM_KPP, &
554 mozart_mosaic_4bin_kpp,mozart_mosaic_4bin_aq_kpp, &
555 p_so2,p_hno3,p_n2o5,p_nh3,p_h2o2,p_o3,p_sulf, &
556 p_facd,p_mepx,p_pacd,p_ora1,p_op1,p_paa,p_hcooh, &
557 p_ch3ooh,p_sulf,p_so4aj,p_nh4aj,p_no3aj, &
558 p_so4_a01,p_so4_a02,p_so4_a03,p_so4_a04,p_nh4_a01, &
559 p_nh4_a02,p_nh4_a03,p_nh4_a04,p_no3_a01,p_no3_a02, &
560 p_no3_a03,p_no3_a04
563 implicit none
564 INTEGER,INTENT(IN) :: chemopt,num_chem
565 real,intent(in) ::p1d,t1d,dz,clw_all1d
566 real,dimension(num_chem), intent(inout) :: tr_up1d
567 ! Aqeuous species pointers INCLUDE File
569 !...........PARAMETERS and their descriptions:
571 INTEGER, PARAMETER :: NGAS = 12 ! number of gas-phase species for AQCHEM
572 INTEGER, PARAMETER :: NAER = 36 ! number of aerosol species for AQCHEM
573 INTEGER, PARAMETER :: NLIQS = 41 ! number of liquid-phase species in AQCHEM
575 !...pointers for the AQCHEM array GAS
577 INTEGER, PARAMETER :: LSO2 = 1 ! Sulfur Dioxide
578 INTEGER, PARAMETER :: LHNO3 = 2 ! Nitric Acid
579 INTEGER, PARAMETER :: LN2O5 = 3 ! Dinitrogen Pentoxide
580 INTEGER, PARAMETER :: LCO2 = 4 ! Carbon Dioxide
581 INTEGER, PARAMETER :: LNH3 = 5 ! Ammonia
582 INTEGER, PARAMETER :: LH2O2 = 6 ! Hydrogen Perioxide
583 INTEGER, PARAMETER :: LO3 = 7 ! Ozone
584 INTEGER, PARAMETER :: LFOA = 8 ! Formic Acid
585 INTEGER, PARAMETER :: LMHP = 9 ! Methyl Hydrogen Peroxide
586 INTEGER, PARAMETER :: LPAA = 10 ! Peroxyacidic Acid
587 INTEGER, PARAMETER :: LH2SO4 = 11 ! Sulfuric Acid
588 INTEGER, PARAMETER :: LHCL = 12 ! Hydrogen Chloride
590 !...pointers for the AQCHEM array AEROSOL
592 INTEGER, PARAMETER :: LSO4AKN = 1 ! Aitken-mode Sulfate
593 INTEGER, PARAMETER :: LSO4ACC = 2 ! Accumulation-mode Sulfate
594 INTEGER, PARAMETER :: LSO4COR = 3 ! Coarse-mode Sulfate
595 INTEGER, PARAMETER :: LNH4AKN = 4 ! Aitken-mode Ammonium
596 INTEGER, PARAMETER :: LNH4ACC = 5 ! Accumulation-mode Ammonium
597 INTEGER, PARAMETER :: LNO3AKN = 6 ! Aitken-mode Nitrate
598 INTEGER, PARAMETER :: LNO3ACC = 7 ! Accumulation-mode Nitrate
599 INTEGER, PARAMETER :: LNO3COR = 8 ! Coarse-mode Nitrate
600 INTEGER, PARAMETER :: LORGAAKN = 9 ! Aitken-mode anthropogenic SOA
601 INTEGER, PARAMETER :: LORGAACC = 10 ! Accumulation-mode anthropogenic SOA
602 INTEGER, PARAMETER :: LORGPAKN = 11 ! Aitken-mode primary organic aerosol
603 INTEGER, PARAMETER :: LORGPACC = 12 ! Accumulation-mode primary organic aerosol
604 INTEGER, PARAMETER :: LORGBAKN = 13 ! Aitken-mode biogenic SOA
605 INTEGER, PARAMETER :: LORGBACC = 14 ! Accumulation-mode biogenic SOA
606 INTEGER, PARAMETER :: LECAKN = 15 ! Aitken-mode elemental carbon
607 INTEGER, PARAMETER :: LECACC = 16 ! Accumulation-mode elemental carbon
608 INTEGER, PARAMETER :: LPRIAKN = 17 ! Aitken-mode primary aerosol
609 INTEGER, PARAMETER :: LPRIACC = 18 ! Accumulation-mode primary aerosol
610 INTEGER, PARAMETER :: LPRICOR = 19 ! Coarse-mode primary aerosol
611 INTEGER, PARAMETER :: LNAAKN = 20 ! Aitken-mode Sodium
612 INTEGER, PARAMETER :: LNAACC = 21 ! Accumulation-mode Sodium
613 INTEGER, PARAMETER :: LNACOR = 22 ! Coarse-mode Sodium
614 INTEGER, PARAMETER :: LCLAKN = 23 ! Aitken-mode Chloride ion
615 INTEGER, PARAMETER :: LCLACC = 24 ! Accumulation-mode Chloride ion
616 INTEGER, PARAMETER :: LCLCOR = 25 ! Coarse-mode Chloride ion
617 INTEGER, PARAMETER :: LNUMAKN = 26 ! Aitken-mode number
618 INTEGER, PARAMETER :: LNUMACC = 27 ! Accumulation-mode number
619 INTEGER, PARAMETER :: LNUMCOR = 28 ! Coarse-mode number
620 INTEGER, PARAMETER :: LSRFAKN = 29 ! Aitken-mode surface area
621 INTEGER, PARAMETER :: LSRFACC = 30 ! Accumulation-mode surface area
622 INTEGER, PARAMETER :: LNACL = 31 ! Sodium Chloride aerosol for AE3 only {depreciated in AE4}
623 INTEGER, PARAMETER :: LCACO3 = 32 ! Calcium Carbonate aerosol (place holder)
624 INTEGER, PARAMETER :: LMGCO3 = 33 ! Magnesium Carbonate aerosol (place holder)
625 INTEGER, PARAMETER :: LA3FE = 34 ! Iron aerosol (place holder)
626 INTEGER, PARAMETER :: LB2MN = 35 ! Manganese aerosol (place holder)
627 INTEGER, PARAMETER :: LK = 36 ! Potassium aerosol (Cl- tracked separately) (place holder)
628 real,parameter::mwdry=28.966 ! Molecular mass of dry air (g/mol)
629 REAL, PARAMETER :: mwso4 = 96.00 ! Molecular mass of SO4-- (g/mol)
630 REAL, PARAMETER :: mwno3 = 62.0 ! Molecular mass of NO3- (g/mol)
631 REAL, PARAMETER :: mwnh4 = 18.0985 ! Molecular mass of NH4+ (g/mol)
632 REAL, PARAMETER :: qcldwtr_cutoff = 1.0e-6 ! kg/m3
634 ! I/O for AQCHEM:
636 real precip,dens,airm,taucld ! Precipitation rate (mm/h)
637 real, dimension (ngas) :: gas,gaswdep
638 real, dimension (naer) :: aerosol,aerwdep
639 real, dimension (nliqs) :: liquid
640 real hpwdep
641 real alfa0,alfa2,alfa3 ! Aerosol scavenging coefficients for Aitken mode
642 real :: &
643 frac_so4(4), frac_no3(4), frac_nh4(4), tot_so4, tot_nh4, tot_no3
645 !
646 ! Aqueous chemistry
647 !
648 !!! TUCCELLA
649 if ((chemopt .EQ. RADM2SORG .OR. chemopt .EQ. RADM2SORG_AQ .OR. chemopt .EQ. RACMSORG_AQ .OR. &
650 chemopt .EQ. RACMSORG_KPP .OR. chemopt .EQ. RADM2SORG_KPP .OR. chemopt .EQ. RACM_ESRLSORG_KPP .OR. &
651 chemopt .EQ. RACM_SOA_VBS_KPP .OR. chemopt .EQ. RADM2SORG_AQCHEM .OR. chemopt .EQ. RACMSORG_AQCHEM_KPP .OR. &
652 chemopt .EQ. CB05_SORG_VBS_AQ_KPP .OR. &
653 chemopt .EQ. RACM_SOA_VBS_HET_KPP .OR. &
654 chemopt .EQ. RACM_ESRLSORG_AQCHEM_KPP .OR. chemopt .EQ. RACM_SOA_VBS_AQCHEM_KPP) &
655 ) then
657 !
658 ! For MADE/SORGAM derived schemes with aqueous chemistry
659 !
661 ! Air mass density
662 dens = 0.1*p1d/t1d*mwdry/8.314472 ! kg/m3
665 ! Column air number density:
666 airm = 1000.0*dens*dz/mwdry ! mol/m2
668 ! Wet scavenging initialization for AQCHEM
670 GASWDEP = 0.0
671 AERWDEP = 0.0
672 HPWDEP = 0.0
674 ! We provide a precipitation rate and aerosol scavenging rates of zero,
675 ! in order to prevent wet scavenging in AQCHEM (it is treated later):
677 precip = 0.0 ! mm/hr
679 alfa0 = 0.0
680 alfa2 = 0.0
681 alfa3 = 0.0
683 ! Gas phase concentrations before aqueous phase chemistry
684 ! (with units conversion ppm -> mol/mol)
686 gas(:) = 0.0
688 gas(lco2) = 380.0e-6
690 gas(lso2) = tr_up1d(p_so2)*1.0e-6
691 gas(lhno3) = tr_up1d(p_hno3)*1.0e-6
692 gas(ln2o5) = tr_up1d(p_n2o5)*1.0e-6
693 gas(lnh3) = tr_up1d(p_nh3)*1.0e-6
694 gas(lh2o2) = tr_up1d(p_h2o2)*1.0e-6
695 gas(lo3) = tr_up1d(p_o3)*1.0e-6
696 gas(lh2so4) = tr_up1d(p_sulf)*1.0e-6
697 if (chemopt==CB05_SORG_VBS_AQ_KPP) then
698 gas(lfoa) = tr_up1d(p_facd)*1.0e-6
699 gas(lmhp) = tr_up1d(p_mepx)*1.0e-6
700 gas(lpaa) = tr_up1d(p_pacd)*1.0e-6
701 else
702 gas(lfoa) = tr_up1d(p_ora1)*1.0e-6
703 gas(lmhp) = tr_up1d(p_op1)*1.0e-6
704 gas(lpaa) = tr_up1d(p_paa)*1.0e-6
705 end if
707 ! Aerosol mass concentrations before aqueous phase chemistry
708 ! (with units conversion ug/kg -> mol/mol). Although AQCHEM
709 ! accounts for much of the aerosol compounds in MADE, they are
710 ! not treated at the moment by AQCHEM, as the mapping between
711 ! the organic compound groups in MADE and AQCHEM is not obvious.
713 aerosol(:) = 0.0
715 ! We assume all accumulation mode particles are activated in cumulus clouds:
717 aerosol(lso4acc) = tr_up1d(p_so4aj)*1.0e-9*mwdry/mwso4
718 aerosol(lnh4acc) = tr_up1d(p_nh4aj)*1.0e-9*mwdry/mwnh4
719 aerosol(lno3acc) = tr_up1d(p_no3aj)*1.0e-9*mwdry/mwno3
721 ! Cloud lifetime:
722 taucld = 1800.0
724 if (clw_all1d*dens .gt. qcldwtr_cutoff) then ! Cloud water > threshold
725 CALL AQCHEM( &
726 t1d, &
727 p1d*100., &
728 taucld, &
729 precip, &
730 clw_all1d*dens, &
731 clw_all1d*dens, &
732 airm, &
733 ALFA0, &
734 ALFA2, &
735 ALFA3, &
736 GAS, &
737 AEROSOL, &
738 LIQUID, &
739 GASWDEP, &
740 AERWDEP, &
741 HPWDEP)
742 endif
744 ! Gas phase concentrations after aqueous phase chemistry
745 ! (with units conversion mol/mol -> ppm)
747 tr_up1d(p_so2) = gas(lso2)*1.0e6
748 tr_up1d(p_hno3) = gas(lhno3)*1.0e6
749 tr_up1d(p_n2o5) = gas(ln2o5)*1.0e6
750 tr_up1d(p_nh3) = gas(lnh3)*1.0e6
751 tr_up1d(p_h2o2) = gas(lh2o2)*1.0e6
752 tr_up1d(p_o3) = gas(lo3)*1.0e6
753 tr_up1d(p_sulf) = gas(lh2so4)*1.0e6
754 if (chemopt==CB05_SORG_VBS_AQ_KPP) then
755 tr_up1d(p_facd) = gas(lfoa)*1.0e6
756 tr_up1d(p_mepx) = gas(lmhp)*1.0e6
757 tr_up1d(p_pacd) = gas(lpaa)*1.0e6
758 else
759 tr_up1d(p_ora1) = gas(lfoa)*1.0e6
760 tr_up1d(p_op1) = gas(lmhp)*1.0e6
761 tr_up1d(p_paa) = gas(lpaa)*1.0e6
762 end if
764 ! Aerosol mass concentrations
765 ! (with units conversion mol/mol -> ug/kg)
767 tr_up1d(p_so4aj) = aerosol(lso4acc)*1.0e9*mwso4/mwdry
768 tr_up1d(p_nh4aj) = aerosol(lnh4acc)*1.0e9*mwnh4/mwdry
769 tr_up1d(p_no3aj) = aerosol(lno3acc)*1.0e9*mwno3/mwdry
770 else if ((chemopt .EQ. mozart_mosaic_4bin_kpp .OR. &
771 chemopt .EQ. mozart_mosaic_4bin_aq_kpp) &
772 ) then
774 !
775 ! For MOSAIC 4bin scheme with aqueous chemistry
776 !
778 ! Air mass density
779 dens = 0.1*p1d/t1d*mwdry/8.314472 ! kg/m3
781 ! Column air number density:
782 airm = 1000.0*dens*dz/mwdry ! mol/m2
784 ! Wet scavenging initialization for AQCHEM
786 GASWDEP = 0.0
787 AERWDEP = 0.0
788 HPWDEP = 0.0
790 ! We provide a precipitation rate and aerosol scavenging rates of zero,
791 ! in order to prevent wet scavenging in AQCHEM (it is treated later):
793 precip = 0.0 ! mm/hr
795 alfa0 = 0.0
796 alfa2 = 0.0
797 alfa3 = 0.0
799 ! Gas phase concentrations before aqueous phase chemistry
800 ! (with units conversion ppm -> mol/mol)
802 gas(:) = 0.0
804 gas(lco2) = 380.0e-6
806 gas(lso2) = tr_up1d(p_so2)*1.0e-6
807 gas(lhno3) = tr_up1d(p_hno3)*1.0e-6
808 gas(ln2o5) = tr_up1d(p_n2o5)*1.0e-6
809 gas(lnh3) = tr_up1d(p_nh3)*1.0e-6
810 gas(lh2o2) = tr_up1d(p_h2o2)*1.0e-6
811 gas(lo3) = tr_up1d(p_o3)*1.0e-6
812 gas(lfoa) = tr_up1d(p_hcooh)*1.0e-6
813 gas(lmhp) = tr_up1d(p_ch3ooh)*1.0e-6
814 gas(lpaa) = tr_up1d(p_paa)*1.0e-6
815 gas(lh2so4) = tr_up1d(p_sulf)*1.0e-6
817 ! Aerosol mass concentrations before aqueous phase chemistry
818 ! (with units conversion ug/kg -> mol/mol). Although AQCHEM
819 ! accounts for much of the aerosol compounds in MADE, they are
820 ! not treated at the moment by AQCHEM, as the mapping between
821 ! the organic compound groups in MADE and AQCHEM is not obvious.
823 aerosol(:) = 0.0
825 ! We assume all particles in bins 2 - 4 are activated in cumulus clouds:
827 ! remember size distribution
828 ! (if none existed before, frac_x is not set, hence distribute equally as default)
829 frac_so4(:) = 0.25
830 frac_nh4(:) = 0.25
831 frac_no3(:) = 0.25
833 tot_so4 = tr_up1d(p_so4_a01)+tr_up1d(p_so4_a02)+&
834 tr_up1d(p_so4_a03)+tr_up1d(p_so4_a04)
835 tot_nh4 = tr_up1d(p_nh4_a01)+tr_up1d(p_nh4_a02)+&
836 tr_up1d(p_nh4_a03)+tr_up1d(p_nh4_a04)
837 tot_no3 = tr_up1d(p_no3_a01)+tr_up1d(p_no3_a02)+&
838 tr_up1d(p_no3_a03)+tr_up1d(p_no3_a04)
840 if (tot_so4 > 0.0) then
841 frac_so4(1) = tr_up1d(p_so4_a01) / tot_so4
842 frac_so4(2) = tr_up1d(p_so4_a02) / tot_so4
843 frac_so4(3) = tr_up1d(p_so4_a03) / tot_so4
844 frac_so4(4) = tr_up1d(p_so4_a04) / tot_so4
845 aerosol(lso4acc) = tot_so4 *1.0e-9*mwdry/mwso4
846 end if
848 if (tot_nh4 > 0.0) then
849 frac_nh4(1) = tr_up1d(p_nh4_a01) / tot_nh4
850 frac_nh4(2) = tr_up1d(p_nh4_a02) / tot_nh4
851 frac_nh4(3) = tr_up1d(p_nh4_a03) / tot_nh4
852 frac_nh4(4) = tr_up1d(p_nh4_a04) / tot_nh4
853 aerosol(lnh4acc) = tot_nh4 *1.0e-9*mwdry/mwnh4
854 end if
856 if (tot_no3 > 0.0) then
857 frac_no3(1) = tr_up1d(p_no3_a01) / tot_no3
858 frac_no3(2) = tr_up1d(p_no3_a02) / tot_no3
859 frac_no3(3) = tr_up1d(p_no3_a03) / tot_no3
860 frac_no3(4) = tr_up1d(p_no3_a04) / tot_no3
861 aerosol(lno3acc) = tot_no3 *1.0e-9*mwdry/mwno3
862 end if
864 ! Cloud lifetime:
865 taucld = 1800.0
867 if (clw_all1d*dens .gt. qcldwtr_cutoff) then ! Cloud water > threshold
868 CALL AQCHEM( &
869 t1d, &
870 p1d*100., &
871 taucld, &
872 precip, &
873 clw_all1d*dens, &
874 clw_all1d*dens, &
875 airm, &
876 ALFA0, &
877 ALFA2, &
878 ALFA3, &
879 GAS, &
880 AEROSOL, &
881 LIQUID, &
882 GASWDEP, &
883 AERWDEP, &
884 HPWDEP)
885 endif
887 ! Gas phase concentrations after aqueous phase chemistry
888 ! (with units conversion mol/mol -> ppm)
890 tr_up1d(p_so2) = gas(lso2)*1.0e6
891 tr_up1d(p_hno3) = gas(lhno3)*1.0e6
892 tr_up1d(p_n2o5) = gas(ln2o5)*1.0e6
893 tr_up1d(p_nh3) = gas(lnh3)*1.0e6
894 tr_up1d(p_h2o2) = gas(lh2o2)*1.0e6
895 tr_up1d(p_o3) = gas(lo3)*1.0e6
896 tr_up1d(p_hcooh) = gas(lfoa)*1.0e6
897 tr_up1d(p_ch3ooh) = gas(lmhp)*1.0e6
898 tr_up1d(p_paa) = gas(lpaa)*1.0e6
899 tr_up1d(p_sulf) = gas(lh2so4)*1.0e6
901 ! Aerosol mass concentrations
902 ! (with units conversion mol/mol -> ug/kg)
904 tr_up1d(p_so4_a01) = aerosol(lso4acc) * frac_so4(1) * 1.0e9*mwso4/mwdry
905 tr_up1d(p_so4_a02) = aerosol(lso4acc) * frac_so4(2) * 1.0e9*mwso4/mwdry
906 tr_up1d(p_so4_a03) = aerosol(lso4acc) * frac_so4(3) * 1.0e9*mwso4/mwdry
907 tr_up1d(p_so4_a04) = aerosol(lso4acc) * frac_so4(4) * 1.0e9*mwso4/mwdry
909 tr_up1d(p_nh4_a01) = aerosol(lnh4acc) * frac_nh4(1) * 1.0e9*mwnh4/mwdry
910 tr_up1d(p_nh4_a02) = aerosol(lnh4acc) * frac_nh4(2) * 1.0e9*mwnh4/mwdry
911 tr_up1d(p_nh4_a03) = aerosol(lnh4acc) * frac_nh4(3) * 1.0e9*mwnh4/mwdry
912 tr_up1d(p_nh4_a04) = aerosol(lnh4acc) * frac_nh4(4) * 1.0e9*mwnh4/mwdry
914 tr_up1d(p_no3_a01) = aerosol(lno3acc) * frac_no3(1) * 1.0e9*mwno3/mwdry
915 tr_up1d(p_no3_a02) = aerosol(lno3acc) * frac_no3(2) * 1.0e9*mwno3/mwdry
916 tr_up1d(p_no3_a03) = aerosol(lno3acc) * frac_no3(3) * 1.0e9*mwno3/mwdry
917 tr_up1d(p_no3_a04) = aerosol(lno3acc) * frac_no3(4) * 1.0e9*mwno3/mwdry
918 endif
919 END SUBROUTINE aqchem_gf
921 SUBROUTINE neg_check_chem(ktop,dt,q,outq,iopt,num_chem, &
922 its,ite,kts,kte,itf)
923 implicit none
924 INTEGER,INTENT(IN) :: iopt,num_chem,its,ite,kts,kte,itf
926 real,dimension(its:ite,kts:kte,num_chem), &
927 intent(inout) :: &
928 outq
929 real,dimension(its:ite,kts:kte,num_chem), &
930 intent(in ) :: &
931 q
932 integer,dimension(its:ite), &
933 intent(in ) :: &
934 ktop
935 real,intent(in ) :: &
936 dt
937 real :: tracermin,tracermax,thresh,qmem,qmemf,qmem2,qtest,qmem1
938 integer :: nv, i, k
939 !
940 ! check whether routine produces negative q's. This can happen, since
941 ! tendencies are calculated based on forced q's. This should have no
942 ! influence on conservation properties, it scales linear through all
943 ! tendencies. Use iopt=0 to test for each tracer seperately, iopt=1
944 ! for a more severe limitation...
945 !
946 ! thresh=epsilc
947 thresh=1.e-30
948 if (iopt==0) then
949 do nv=2,num_chem
950 do 100 i=its,itf
951 tracermin=q(i,kts,nv)
952 tracermax=q(i,kts,nv)
953 do k=kts+1,kte-1
954 tracermin=min(tracermin,q(i,k,nv))
955 tracermax=max(tracermax,q(i,k,nv))
956 enddo ! k
957 tracermin=max(tracermin,thresh)
958 qmemf=1.
959 !
960 ! first check for minimum restriction
961 !
962 do k=kts,ktop(i)
963 !
964 ! tracer tendency
965 !
966 qmem=outq(i,k,nv)
967 !
968 ! only necessary if there is a tendency
969 !
970 if(qmem.lt.0.)then
971 qtest=q(i,k,nv)+outq(i,k,nv)*dt
972 if(qtest.lt.tracermin)then
973 !
974 ! qmem2 would be the maximum allowable tendency
975 !
976 qmem1=outq(i,k,nv)
977 qmem2=(tracermin-q(i,k,nv))/dt
978 qmemf=min(qmemf,qmem2/qmem1)
979 if(qmemf.gt.1.)print *,'something wrong in negct_1',qmem2,qmem1
980 qmemf=max(qmemf,0.)
981 endif ! qtest
982 endif ! qmem
983 enddo ! k
984 do k=kts,ktop(i)
985 outq(i,k,nv)=outq(i,k,nv)*qmemf
986 enddo ! k
987 !
988 ! now check max
989 !
990 qmemf=1.
991 do k=kts,ktop(i)
992 !
993 ! tracer tendency
994 !
995 qmem=outq(i,k,nv)
996 !
997 ! only necessary if there is a tendency
998 !
999 if(qmem.gt.0.)then
1000 qtest=q(i,k,nv)+outq(i,k,nv)*dt
1001 if(qtest.gt.tracermax)then
1002 !
1003 ! qmem2 would be the maximum allowable tendency
1004 !
1005 qmem1=outq(i,k,nv)
1006 qmem2=(tracermax-q(i,k,nv))/dt
1007 qmemf=min(qmemf,qmem2/qmem1)
1008 if(qmemf.gt.1.)print *,'something wrong in negct_2',qmem2,qmem1
1009 qmemf=max(qmemf,0.)
1010 endif ! qtest
1011 endif ! qmem
1012 enddo ! k
1013 do k=kts,ktop(i)
1014 outq(i,k,nv)=outq(i,k,nv)*qmemf
1015 enddo ! k
1016 100 continue ! i
1017 enddo ! nv
1018 !
1019 ! ELSE
1020 !
1021 elseif(iopt.eq.1)then
1022 do i=its,itf
1023 qmemf=1.
1024 do k=kts,ktop(i)
1025 do nv=2,num_chem
1026 !
1027 ! tracer tendency
1028 !
1029 qmem=outq(i,k,nv)
1030 !
1031 ! only necessary if tendency is larger than zero
1032 !
1033 if(qmem.lt.0.)then
1034 qtest=q(i,k,nv)+outq(i,k,nv)*dt
1035 if(qtest.lt.thresh)then
1036 !
1037 ! qmem2 would be the maximum allowable tendency
1038 !
1039 qmem1=outq(i,k,nv)
1040 qmem2=(thresh-q(i,k,nv))/dt
1041 qmemf=min(qmemf,qmem2/qmem1)
1042 qmemf=max(0.,qmemf)
1043 endif ! qtest
1044 endif ! qmem
1045 enddo ! nv
1046 enddo ! k
1047 do nv=2,num_chem
1048 do k=kts,ktop(i)
1049 outq(i,k,nv)=outq(i,k,nv)*qmemf
1050 enddo ! k
1051 enddo ! nv
1052 enddo ! i
1053 endif !iopt
1055 END SUBROUTINE neg_check_chem
1056 SUBROUTINE get_cloud_bc_chem(mzp,array,x_aver,k22,add)
1057 implicit none
1058 integer, intent(in) :: mzp,k22
1059 real , intent(in) :: array(mzp)
1060 real , optional , intent(in) :: add
1061 real , intent(out) :: x_aver
1062 integer :: i,local_order_aver,order_aver
1064 !-- dimension of the average
1065 !-- a) to pick the value at k22 level, instead of a average between
1066 !-- k22-order_aver, ..., k22-1, k22 set order_aver=1)
1067 !-- b) to average between 1 and k22 => set order_aver = k22
1068 order_aver = 3 !=> average between k22, k22-1 and k22-2
1070 local_order_aver=min(k22,order_aver)
1072 x_aver=0.
1073 do i = 1,local_order_aver
1074 x_aver = x_aver + array(k22-i+1)
1075 enddo
1076 x_aver = x_aver/float(local_order_aver)
1077 if(present(add)) x_aver = x_aver + add
1079 end SUBROUTINE get_cloud_bc_chem
1080 #endif
1081 !-- GF CTRAN --
1082 END MODULE module_cu_gf_ctrans