1 | !$Id: sethet.F90 112 2009-01-28 16:40:56Z acosce $ |
---|
2 | !! ========================================================================= |
---|
3 | !! INCA - INteraction with Chemistry and Aerosols |
---|
4 | !! |
---|
5 | !! Copyright Laboratoire des Sciences du Climat et de l'Environnement (LSCE) |
---|
6 | !! Unite mixte CEA-CNRS-UVSQ |
---|
7 | !! |
---|
8 | !! Contributors to this INCA subroutine: |
---|
9 | !! |
---|
10 | !! Didier Hauglustaine, LSCE, hauglustaine@cea.fr |
---|
11 | !! Xue-Xi Tie, NCAR |
---|
12 | !! |
---|
13 | !! Anne Cozic, LSCE, anne.cozic@cea.fr |
---|
14 | !! Yann Meurdesoif, LSCE, yann.meurdesoif@cea.fr |
---|
15 | !! |
---|
16 | !! This software is a computer program whose purpose is to simulate the |
---|
17 | !! atmospheric gas phase and aerosol composition. The model is designed to be |
---|
18 | !! used within a transport model or a general circulation model. This version |
---|
19 | !! of INCA was designed to be coupled to the LMDz GCM. LMDz-INCA accounts |
---|
20 | !! for emissions, transport (resolved and sub-grid scale), photochemical |
---|
21 | !! transformations, and scavenging (dry deposition and washout) of chemical |
---|
22 | !! species and aerosols interactively in the GCM. Several versions of the INCA |
---|
23 | !! model are currently used depending on the envisaged applications with the |
---|
24 | !! chemistry-climate model. |
---|
25 | !! |
---|
26 | !! This software is governed by the CeCILL license under French law and |
---|
27 | !! abiding by the rules of distribution of free software. You can use, |
---|
28 | !! modify and/ or redistribute the software under the terms of the CeCILL |
---|
29 | !! license as circulated by CEA, CNRS and INRIA at the following URL |
---|
30 | !! "http://www.cecill.info". |
---|
31 | !! |
---|
32 | !! As a counterpart to the access to the source code and rights to copy, |
---|
33 | !! modify and redistribute granted by the license, users are provided only |
---|
34 | !! with a limited warranty and the software's author, the holder of the |
---|
35 | !! economic rights, and the successive licensors have only limited |
---|
36 | !! liability. |
---|
37 | !! |
---|
38 | !! In this respect, the user's attention is drawn to the risks associated |
---|
39 | !! with loading, using, modifying and/or developing or reproducing the |
---|
40 | !! software by the user in light of its specific status of free software, |
---|
41 | !! that may mean that it is complicated to manipulate, and that also |
---|
42 | !! therefore means that it is reserved for developers and experienced |
---|
43 | !! professionals having in-depth computer knowledge. Users are therefore |
---|
44 | !! encouraged to load and test the software's suitability as regards their |
---|
45 | !! requirements in conditions enabling the security of their systems and/or |
---|
46 | !! data to be ensured and, more generally, to use and operate it in the |
---|
47 | !! same conditions as regards security. |
---|
48 | !! |
---|
49 | !! The fact that you are presently reading this means that you have had |
---|
50 | !! knowledge of the CeCILL license and that you accept its terms. |
---|
51 | !! ========================================================================= |
---|
52 | |
---|
53 | #include <inca_define.h> |
---|
54 | |
---|
55 | |
---|
56 | SUBROUTINE SETHET( & |
---|
57 | het_rates , & |
---|
58 | press , & |
---|
59 | pdel , & |
---|
60 | lat , & |
---|
61 | zmid , & |
---|
62 | tfld , & |
---|
63 | delt , & |
---|
64 | xhnm , & |
---|
65 | flxr , & |
---|
66 | flxs , & |
---|
67 | flxupd , & |
---|
68 | cldtop , & |
---|
69 | cldbot , & |
---|
70 | cldfr , & |
---|
71 | index , & |
---|
72 | qin ) |
---|
73 | !----------------------------------------------------------------------- |
---|
74 | ! ... In-cloud and below-cloud scavenging of soluble species. |
---|
75 | ! Xue-Xi Tie, NCAR, 1998. |
---|
76 | ! Didier Hauglustaine, IPSL, 2001. |
---|
77 | ! Didier Hauglustaine, IPSL, 05-2002. |
---|
78 | !----------------------------------------------------------------------- |
---|
79 | |
---|
80 | USE CHEM_CONS, ONLY : gravit, uma |
---|
81 | USE SPECIES_NAMES |
---|
82 | |
---|
83 | #if defined(AER) || defined(NMHC) |
---|
84 | USE DRYDEP_ARRAYS, ONLY : heff_3D |
---|
85 | !gaf the Boltzmann Constant is included |
---|
86 | #endif |
---|
87 | |
---|
88 | #ifdef NMHC |
---|
89 | USE AIRPLANE_SRC, ONLY : itrop |
---|
90 | #endif |
---|
91 | |
---|
92 | USE INCA_DIM |
---|
93 | USE CHEM_MODS, ONLY : invariants |
---|
94 | USE DRYDEP_PARAMETERS, ONLY : ndep |
---|
95 | USE INPUT_DATA_TABLES, ONLY : spec_map |
---|
96 | USE RATE_INDEX_MOD |
---|
97 | |
---|
98 | IMPLICIT NONE |
---|
99 | |
---|
100 | !----------------------------------------------------------------------- |
---|
101 | ! ... Dummy arguments |
---|
102 | !----------------------------------------------------------------------- |
---|
103 | INTEGER, INTENT(in) :: lat ! latitude index |
---|
104 | INTEGER, INTENT(in) :: INDEX ! index = 1 for stratiform clouds |
---|
105 | ! index = 2 for convective clouds |
---|
106 | INTEGER, INTENT(in) :: cldtop(PLON) ! cloud top level ( 1 ... 19 ) |
---|
107 | !gaf cloud bottom is included: cldbot |
---|
108 | INTEGER, INTENT(in) :: cldbot(PLON) ! cloud bot level ( 1 ... 19 ) |
---|
109 | REAL, INTENT(in) :: delt ! time step ( s ) |
---|
110 | REAL, INTENT(in) :: press(PLON,PLEV) ! midpoint pressure Pa |
---|
111 | REAL, INTENT(in) :: pdel(PLON,PLEV) ! delta pressure (Pa) |
---|
112 | REAL, INTENT(in) :: qin(PLON,PLEV,PCNST) ! xported species (vmr) |
---|
113 | REAL, INTENT(in) :: zmid(PLON,PLEV) ! midpoint geopot |
---|
114 | REAL, INTENT(in) :: tfld(PLON,PLEV) ! temperature |
---|
115 | REAL, INTENT(in) :: xhnm(PLON,PLEV) ! total density (/cm**3) |
---|
116 | REAL, INTENT(inout) :: het_rates(PLON,PLEV,HETCNT) ! rainout loss rates |
---|
117 | REAL, INTENT(inout) :: flxr(PLON,PLEVP) !liquid water flx kgH2O/m2/s |
---|
118 | REAL, INTENT(inout) :: flxs(PLON,PLEVP) !solid water flx kgH2O/m2/s |
---|
119 | REAL, INTENT(in) :: flxupd(PLON,PLEV) !entrainment flx kgAIR/m2/s |
---|
120 | REAL, INTENT(in) :: cldfr(PLON,PLEV) !cloud fraction |
---|
121 | |
---|
122 | !----------------------------------------------------------------------- |
---|
123 | ! ... Local variables |
---|
124 | !----------------------------------------------------------------------- |
---|
125 | REAL, PARAMETER :: RD = 287.04 ! ideal gas constant/molarmass of air J/molkg |
---|
126 | REAL, PARAMETER :: drym = 28.966 |
---|
127 | REAL, PARAMETER :: Rg = 8.205e-2 ! atm cm3/K/M/g |
---|
128 | ! ... The following numbers are criticial and should be calculated |
---|
129 | ! instead of fixed. A size distribution should be used for |
---|
130 | ! rain drops based on the rain intensity (Seinfeld and Pandis, P. 831). |
---|
131 | ! Then the terminal velocity could be calculated as well (Seinfeild |
---|
132 | ! and Pandis, P. 468). See also Roelofs and Lelieveld (1995). |
---|
133 | ! This will be done in a future version. For now we use typical numbers |
---|
134 | ! provided in Brasseur et al. (1988). --DH, 2001 |
---|
135 | |
---|
136 | !DH 11/2011 These variables updated based on Seinfeld and Pandis. |
---|
137 | REAL, PARAMETER :: xrm = 0.100 |
---|
138 | REAL, PARAMETER :: xum = 300. |
---|
139 | |
---|
140 | REAL, PARAMETER :: xvv = 0.146 |
---|
141 | REAL, PARAMETER :: xdg = 0.112 |
---|
142 | |
---|
143 | ! Here as well, we need something better for LWC. |
---|
144 | REAL, DIMENSION(2) :: lwc = (/ .5 , 2. /) |
---|
145 | #if defined(AERONLY) |
---|
146 | INTEGER :: itrop(PLON) ! to exclude Stratosphere wet dep calculations |
---|
147 | #endif |
---|
148 | INTEGER :: i, k, ktrop, kk |
---|
149 | REAL :: cst, alpha, dz |
---|
150 | REAL, SAVE :: xkgm |
---|
151 | INTEGER, SAVE :: index_NH3 |
---|
152 | INTEGER, SAVE :: index_APp1g, index_APp2g, index_ARp1g, index_ARp2g |
---|
153 | !$OMP THREADPRIVATE(xkgm, index_NH3, index_APp1g, index_APp2g, index_ARp1g, index_ARp2g) |
---|
154 | REAL :: all1, all2, stay |
---|
155 | REAL :: xeqca1, xeqca2, xca1, xca2, xdtm |
---|
156 | REAL :: xxx1, xxx2, yhno3, yh2o2 |
---|
157 | REAL, DIMENSION(PLON) :: xk0, xk1, xk2, work1, work2 |
---|
158 | REAL, DIMENSION(PLON) :: hplus_inv |
---|
159 | REAL, DIMENSION(PLEV) :: xgas1, xgas2 |
---|
160 | REAL, DIMENSION(PLON,PLEV) :: delz, xhno3, xh2o2, xliq, wh2o |
---|
161 | REAL, DIMENSION(PLON,PLEV) :: xhenhno3 ! henry constants |
---|
162 | REAL, DIMENSION(PLON,PLEV) :: xhenh2o2 |
---|
163 | |
---|
164 | !----------------------------------------------------------------------- |
---|
165 | ! ... effective Henry's Law Constants |
---|
166 | !----------------------------------------------------------------------- |
---|
167 | ! effective Henry's Law Constants are used for 19 species only |
---|
168 | ! they are (in that order): |
---|
169 | ! HNO3, H2O2, HNO2, HNO4, CH3OOH, CH3OH, CH2O, C2H5OH, CH3CHO, |
---|
170 | ! CH3COOH, CH3COOOH, CH3COCHO, CH3COCH3, C2H5OOH, MVK, MEK, |
---|
171 | ! PAN, ONITR, ONITU |
---|
172 | ! the mapping garanties the correct hand over |
---|
173 | INTEGER, PARAMETER :: n_effhetrxt = 19 |
---|
174 | INTEGER, PARAMETER :: n_hetrxt = HETCNT-1 |
---|
175 | INTEGER, DIMENSION(n_effhetrxt), SAVE :: mapping1, mapping2 |
---|
176 | ! mapping1 permet de retrouver les variables avec constantes d'henry listees ci-dessus |
---|
177 | ! dans la liste ndep des variables pour lesquelles on a calcule heff_3D (mkdvel) |
---|
178 | ! mapping2 permet de mettre les variables heterogene dans l'ordre que l'on veut |
---|
179 | ! dans le fichier inca***.def. Nous ne sommes plus contraint d'avoir hno3 |
---|
180 | ! en premier et pb210 en dernier. Ni d'avoir en premier les especes de la liste ci-dessus |
---|
181 | |
---|
182 | INTEGER, DIMENSION(n_effhetrxt), SAVE :: het_map |
---|
183 | !$OMP THREADPRIVATE(mapping1,mapping2, het_map) |
---|
184 | ! field containing all effective Henry's Law constants |
---|
185 | ! first entry always has to by HNO3 |
---|
186 | ! last slot reserved for Pb210 |
---|
187 | REAL, DIMENSION(PLON,PLEV,HETCNT) :: xhenconst |
---|
188 | |
---|
189 | |
---|
190 | REAL, DIMENSION(PLON,PLEV) :: wrate |
---|
191 | REAL, DIMENSION(PLON,PLEV) :: zrho |
---|
192 | REAL, DIMENSION(PLON,PLEV) :: totmass |
---|
193 | REAL, DIMENSION(PLON,PLEV) :: massupd |
---|
194 | REAL, DIMENSION(PLON,PLEV) :: flxdwn |
---|
195 | REAL, DIMENSION(PLON,PLEV) :: scaveff |
---|
196 | |
---|
197 | LOGICAL, SAVE :: entered = .FALSE. |
---|
198 | !$OMP THREADPRIVATE(entered) |
---|
199 | #ifdef DUSS |
---|
200 | INTEGER, SAVE :: het_HNO3 |
---|
201 | !$OMP THREADPRIVATE(het_HNO3) |
---|
202 | #endif |
---|
203 | |
---|
204 | !-------------------------------------------------- |
---|
205 | ! VERSION AER calcul des constantes d'henry |
---|
206 | ! et het_rate = 0 |
---|
207 | !-------------------------------------------------- |
---|
208 | |
---|
209 | ! changed for AERONLY MS July 07, AER needs to remove SO2, see aeronly section at end |
---|
210 | #if defined(GES) |
---|
211 | het_rates(:,:,:) = 0. |
---|
212 | #else |
---|
213 | |
---|
214 | !----------------------------------------------------------------- |
---|
215 | ! ... PART 0, for input need |
---|
216 | !----------------------------------------------------------------- |
---|
217 | ! wrate= rainwater tendency kg_H2O/kg_air/s |
---|
218 | ! wh2o = rate of gas h2o into rain water #/cm3/s |
---|
219 | ! xliq = liquid rain water content in the drop g_H2O/m3 |
---|
220 | ! delz = delta altitude m then cm |
---|
221 | ! xhno3= initial hno3 concentration #/cm3 |
---|
222 | ! xh2o2= initial h2o2 concentration #/cm3 |
---|
223 | ! |
---|
224 | ! xrm = mean diameter of drop cm |
---|
225 | ! xum = mean terminal velocity cm/s |
---|
226 | ! xvv = kinetic viscosity cm2/s |
---|
227 | ! xdg = mass transport coefficient cm/s |
---|
228 | ! xkgm = mass flux on rain drop |
---|
229 | ! |
---|
230 | ! lwc = liquid water content g_H2O/m3 |
---|
231 | !----------------------------------------------------------------- |
---|
232 | |
---|
233 | IF( .NOT. entered ) THEN |
---|
234 | xkgm = xdg/xrm*2.+xdg/xrm*0.6*SQRT(xrm*xum/xvv)*(xvv/xdg)**(1./3.) |
---|
235 | |
---|
236 | #ifndef DUSS |
---|
237 | #ifdef NMHC |
---|
238 | het_map = (/id_HNO3, id_H2O2, id_HNO2, id_HNO4, id_CH3OOH, & |
---|
239 | id_CH3OH, id_CH2O, id_C2H5OH, id_CH3CHO, id_CH3COOH, & |
---|
240 | id_CH3COOOH, id_CH3COCHO, id_CH3COCH3, id_C2H5OOH, & |
---|
241 | id_MVK, id_MEK, id_PAN, id_ONITR, id_ONITU /) |
---|
242 | |
---|
243 | |
---|
244 | DO k=1,ndep |
---|
245 | DO i=1,n_effhetrxt |
---|
246 | if ( het_map(i) .eq. spec_map(k)) mapping1(i) = k |
---|
247 | ENDDO |
---|
248 | #ifdef AER |
---|
249 | if (spec_map(k) .eq. id_NH3) index_NH3 = k |
---|
250 | if (spec_map(k) .eq. id_APp1g) index_APp1g = k |
---|
251 | if (spec_map(k) .eq. id_APp2g) index_APp2g = k |
---|
252 | if (spec_map(k) .eq. id_ARp1g) index_ARp1g = k |
---|
253 | if (spec_map(k) .eq. id_ARp2g) index_ARp2g = k |
---|
254 | #endif |
---|
255 | ENDDO |
---|
256 | |
---|
257 | do k=1,HETCNT |
---|
258 | do i=1,n_effhetrxt |
---|
259 | if (tracnam(het_map(i)) .eq. hetname(k)) mapping2(i) = k |
---|
260 | enddo |
---|
261 | enddo |
---|
262 | |
---|
263 | |
---|
264 | #else |
---|
265 | DO k=1,ndep |
---|
266 | if (spec_map(k) .eq. id_NH3) index_NH3 = k |
---|
267 | ENDDO |
---|
268 | |
---|
269 | #endif |
---|
270 | #else |
---|
271 | het_HNO3 = het_PB210 |
---|
272 | #endif |
---|
273 | entered = .TRUE. |
---|
274 | END IF |
---|
275 | |
---|
276 | !----------------------------------------------------------------- |
---|
277 | ! ... First, we derive the liquid-solid water tendency based |
---|
278 | ! on GCM variables |
---|
279 | !----------------------------------------------------------------- |
---|
280 | |
---|
281 | zrho = (xhnm*1.e6) * drym * uma |
---|
282 | delz = pdel / zrho / gravit |
---|
283 | |
---|
284 | #if defined(AERONLY) |
---|
285 | DO i = 1, PLON |
---|
286 | itrop(i)=nint(3./4.*PLEVP) |
---|
287 | END DO |
---|
288 | #endif |
---|
289 | |
---|
290 | DO i = 1, PLON |
---|
291 | ktrop = itrop(i) |
---|
292 | flxr(i,ktrop:PLEVP) = 0. |
---|
293 | flxs(i,ktrop:PLEVP) = 0. |
---|
294 | END DO |
---|
295 | |
---|
296 | DO k = 1, PLEV |
---|
297 | wrate(:,k) = flxr(:,k)-flxr(:,k+1)+flxs(:,k)-flxs(:,k+1) |
---|
298 | flxdwn(:,k) = flxr(:,k)+flxs(:,k) |
---|
299 | END DO |
---|
300 | |
---|
301 | ! [kg_h2o/m2/s] [1/m] [m3/kg_air] -> [kg_h2o/kg_air/s] |
---|
302 | wrate = MAX(0., wrate) /delz/zrho |
---|
303 | delz = delz * 1.e2 !delz from m to cm |
---|
304 | |
---|
305 | DO k = 1,PLEV |
---|
306 | wh2o(:,k) = 29.*wrate(1:PLON,k)*xhnm(:,k) / 18. |
---|
307 | xliq(:,k) = wrate(1:PLON,k)*delt*xhnm(:,k)/6.023e23*29.* 1.e6 |
---|
308 | |
---|
309 | #if !defined(AERONLY) |
---|
310 | xhno3(:,k) = qin(:,k,id_hno3) * xhnm(:,k) |
---|
311 | xh2o2(:,k) = qin(:,k,id_h2o2) * xhnm(:,k) |
---|
312 | #else |
---|
313 | #if !defined(DUSS) |
---|
314 | xhno3(:,k) = qin(:,k,id_hno3) * xhnm(:,k) |
---|
315 | xh2o2(:,k) = invariants(:,k,inv_H2O2)* xhnm(:,k)/invariants(:,k,inv_M) |
---|
316 | #else |
---|
317 | xhno3(:,k) = 1.e-9 * xhnm(:,k) |
---|
318 | xh2o2(:,k) = 1.e-10 * xhnm(:,k) |
---|
319 | #endif |
---|
320 | #endif |
---|
321 | END DO |
---|
322 | |
---|
323 | !----------------------------------------------------------------- |
---|
324 | ! ... PART 0b, Temperature dependent Henry Law constants |
---|
325 | ! based on Chameides (1984) and R. Sander |
---|
326 | !----------------------------------------------------------------- |
---|
327 | |
---|
328 | xhenconst = 0. |
---|
329 | |
---|
330 | #ifdef NMHC |
---|
331 | ! hand over of calculated effective HLCs |
---|
332 | DO i = 1, n_effhetrxt |
---|
333 | xhenconst(:,:,mapping2(i)) = heff_3D(:,:,mapping1(i)) |
---|
334 | END DO |
---|
335 | ! assignment of the rest |
---|
336 | xhenconst(:,:,het_C3H7OOH) = xhenconst(:,:,het_C2H5OOH) ! HLC(c3h7ooh) = HLC(c2h5ooh) |
---|
337 | xhenconst(:,:,het_PROPEOOH) = xhenconst(:,:,het_C2H5OOH) ! HLC(propeooh) = HLC(c2h5ooh) |
---|
338 | xhenconst(:,:,het_PROPAOOH) = xhenconst(:,:,het_C2H5OOH) ! HLC(propaooh) = HLC(c2h5ooh) |
---|
339 | xhenconst(:,:,het_PCHO) = xhenconst(:,:,het_CH3CHO) ! HLC(pcho) = HLC(ch3cho) |
---|
340 | xhenconst(:,:,het_MACROOH) = xhenconst(:,:,het_C2H5OOH) ! HLC(macrooh) = HLC(c2h5ooh) |
---|
341 | xhenconst(:,:,het_MEKOOH) = xhenconst(:,:,het_C2H5OOH) ! HLC(mekooh) = HLC(c2h5ooh) |
---|
342 | xhenconst(:,:,het_ALKANOOH) = xhenconst(:,:,het_C2H5OOH) ! HLC(alkanooh) = HLC(c2h5ooh) |
---|
343 | xhenconst(:,:,het_ALKENOOH) = xhenconst(:,:,het_C2H5OOH) ! HLC(alkenooh) = HLC(c2h5ooh) |
---|
344 | xhenconst(:,:,het_AROMOOH) = xhenconst(:,:,het_C2H5OOH) ! HLC(aromooh) = HLC(c2h5ooh) |
---|
345 | xhenconst(:,:,het_XOOH) = xhenconst(:,:,het_C2H5OOH) ! HLC(xooh) = HLC(c2h5ooh) |
---|
346 | |
---|
347 | #ifdef AER |
---|
348 | ! SOA precursors |
---|
349 | xhenconst(:,:,het_APp1g) = heff_3D(:,:,index_APp1g) |
---|
350 | xhenconst(:,:,het_APp2g) = heff_3D(:,:,index_APp2g) |
---|
351 | xhenconst(:,:,het_ARp1g) = heff_3D(:,:,index_ARp1g) |
---|
352 | xhenconst(:,:,het_ARp2g) = heff_3D(:,:,index_ARp2g) |
---|
353 | #endif |
---|
354 | |
---|
355 | #endif |
---|
356 | |
---|
357 | #if !defined(DUSS) && defined(AER) |
---|
358 | DO k =1, PLEV |
---|
359 | work1(:) = 1. / tfld(:PLON,k) - 1. / 298. |
---|
360 | zrho(:,k) = press(:,k)/(tfld(:,k)*RD) |
---|
361 | hplus_inv(:) = 2.*qin(:,k,id_asso4m)*zrho(:,k)*1.e-3/(lwc(index)*drym) |
---|
362 | hplus_inv(:) = 1./MIN(1.e-3,MAX(2.51188e-6,hplus_inv(:))) |
---|
363 | |
---|
364 | #ifdef STRAT |
---|
365 | xhenconst(:,k,het_HCl ) = 2.0e6*EXP(9000.*work1(:)) |
---|
366 | xhenconst(:,k,het_ClONO2) = xhenconst(:,k,het_HNO3) |
---|
367 | xhenconst(:,k,het_HOCl ) = 660.*EXP( 5900.*work1(:)) |
---|
368 | xhenconst(:,k,het_ClNO2 ) = 4.6e-2 |
---|
369 | xhenconst(:,k,het_BrONO2) = xhenconst(:,k,het_HNO3) |
---|
370 | xhenconst(:,k,het_HOBr ) = 6.1e3 |
---|
371 | xhenconst(:,k,het_BrCl ) = 9.4e-1*EXP(5600.*work1(:)) |
---|
372 | #endif |
---|
373 | xhenconst(:,k,het_SO2) = 1.2*EXP(2900.*work1(:)) |
---|
374 | ! Modif Dimitri Caro 14/12/2005 |
---|
375 | xk0(:) = 1.7e-2 *EXP(2090.*work1(:)) |
---|
376 | xk1(:) = 6.0e-8 *EXP(1120.*work1(:)) |
---|
377 | xhenconst(:,k,het_SO2) = & |
---|
378 | xhenconst(:,k,het_SO2)* (1.+xk0(:)*hplus_inv(:)+xk0(:)*xk1(:)*hplus_inv(:)**2) |
---|
379 | |
---|
380 | xhenconst(:,k,het_H2S) = 1.0e-3*EXP(2300*work1(:)) |
---|
381 | xhenconst(:,k,het_H2S) = & |
---|
382 | xhenconst(:,k,het_H2S)*(1.+5.7e-8*hplus_inv(:)+5.7e-8*1.e-13*hplus_inv(:)**2) |
---|
383 | |
---|
384 | xhenconst(:,k,het_DMS) = 4.8e-01*EXP(3100*work1(:)) |
---|
385 | |
---|
386 | xhenconst(:,k,het_DMSO) = 5.0e+04 |
---|
387 | |
---|
388 | ENDDO |
---|
389 | xhenconst(:,:,het_NH3) = heff_3D(:,:,index_NH3) |
---|
390 | |
---|
391 | |
---|
392 | #endif |
---|
393 | |
---|
394 | #ifdef AERONLY |
---|
395 | DO k = 1,PLEV |
---|
396 | work1(:) = 1. / tfld(:PLON,k) - 1. / 298. |
---|
397 | xk0(:) = 2.6e6*EXP( 8700.*work1(:) ) |
---|
398 | xhenhno3(:,k) = xk0(:) / 5.81756e6 * 7.e11 |
---|
399 | xk2(:) = 9.7e4*EXP( 6600.*work1(:) ) |
---|
400 | xhenh2o2(:,k) = xk2(:) / 178695. * 2.e5 |
---|
401 | END DO |
---|
402 | #else |
---|
403 | xhenhno3(:,:) = xhenconst(:,:,het_HNO3) |
---|
404 | xhenh2o2(:,:) = xhenconst(:,:,het_H2O2) |
---|
405 | #endif |
---|
406 | |
---|
407 | SELECT CASE (index) |
---|
408 | CASE(1) |
---|
409 | !----------------------------------------------------------------- |
---|
410 | ! ... PART 1, solve for high henry constant ( HNO3, H2O2) |
---|
411 | ! This includes in-cloud and below-cloud scavenging. |
---|
412 | !----------------------------------------------------------------- |
---|
413 | DO i = 1,PLON |
---|
414 | xgas1(:) = xhno3(i,:) ! xgas will change during |
---|
415 | xgas2(:) = xh2o2(i,:) ! different levels wash |
---|
416 | DO kk = PLEV, 1, -1 |
---|
417 | stay = 1. |
---|
418 | IF( wh2o(i,kk) /= 0. ) THEN ! finding rain cloud |
---|
419 | all1 = 0. ! accumulation to justisfy saturation |
---|
420 | all2 = 0. |
---|
421 | stay = (zmid(i,kk)*1.e5)/(xum*delt) |
---|
422 | stay = MIN( stay,1.0 ) |
---|
423 | |
---|
424 | !---------------------------------------------------------------- |
---|
425 | ! ... Calculate the saturation concentration eqca |
---|
426 | !---------------------------------------------------------------- |
---|
427 | DO k = kk, 1, -1 ! cal washout below cloud |
---|
428 | xeqca1 = & |
---|
429 | xgas1(k)/(xliq(i,kk)*6.023e20+ 1./(xhenhno3(i,k)*1.36e-22*tfld(i,k)))* & |
---|
430 | xliq(i,kk)*6.023e20 |
---|
431 | xeqca2 = & |
---|
432 | xgas2(k)/(xliq(i,kk)*6.023e20+1./(xhenh2o2(i,k)*1.36e-22*tfld(i,k)))* & |
---|
433 | xliq(i,kk)*6.023e20 |
---|
434 | !----------------------------------------------------------------- |
---|
435 | ! ... Calculate ca; inside cloud concentration in #/cm3(air) |
---|
436 | !----------------------------------------------------------------- |
---|
437 | xca1 = 6.*xkgm*xgas1(k)/(xrm*xum)*delz(i,k) & |
---|
438 | * xliq(i,kk) * 1.e-6 !!! * MIN(cldfr(i,k),1.) |
---|
439 | xca2 = 6.*xkgm*xgas2(k)/(xrm*xum)*delz(i,k) & |
---|
440 | * xliq(i,kk) * 1.e-6 !!! * MIN(cldfr(i,k),1.) |
---|
441 | |
---|
442 | !----------------------------------------------------------------- |
---|
443 | ! ... If is not saturated |
---|
444 | ! hno3(gas)_new = hno3(gas)_old - hno3(h2o) |
---|
445 | ! otherwise |
---|
446 | ! hno3(gas)_new = hno3(gas)_old |
---|
447 | !----------------------------------------------------------------- |
---|
448 | all1 = all1 + xca1 |
---|
449 | all2 = all2 + xca2 |
---|
450 | IF( all1 < xeqca1 ) THEN |
---|
451 | xgas1(k) = MAX( xgas1(k) - xca1,0. ) |
---|
452 | END IF |
---|
453 | IF( all2 < xeqca2 ) THEN |
---|
454 | xgas2(k) = MAX( xgas2(k) - xca2,0. ) |
---|
455 | END IF |
---|
456 | END DO |
---|
457 | END IF |
---|
458 | |
---|
459 | !----------------------------------------------------------------- |
---|
460 | ! ... Calculate the lifetime of washout (second) |
---|
461 | ! after all layers washout |
---|
462 | ! the concentration of hno3 is reduced |
---|
463 | ! then the lifetime xtt is calculated by |
---|
464 | ! |
---|
465 | ! xtt = (xhno3(ini) - xgas1(new))/(dt*xhno3(ini)) |
---|
466 | ! where dt = passing time (s) in vertical |
---|
467 | ! path below the cloud |
---|
468 | ! dt = dz(cm)/um(cm/s) |
---|
469 | !----------------------------------------------------------------- |
---|
470 | xdtm = delz(i,kk) / xum ! the traveling time in each dz |
---|
471 | xxx1 = (xhno3(i,kk) - xgas1(kk)) |
---|
472 | xxx2 = (xh2o2(i,kk) - xgas2(kk)) |
---|
473 | IF( xxx1 /= 0. ) THEN ! if no washout lifetime = 1.e29 |
---|
474 | yhno3 = xhno3(i,kk)/xxx1 * xdtm |
---|
475 | ELSE |
---|
476 | yhno3 = 1.e29 |
---|
477 | END IF |
---|
478 | IF( xxx2 /= 0. ) THEN ! if no washout lifetime = 1.e29 |
---|
479 | yh2o2 = xh2o2(i,kk)/xxx2 * xdtm |
---|
480 | ELSE |
---|
481 | yh2o2 = 1.e29 |
---|
482 | END IF |
---|
483 | |
---|
484 | IF (tfld(i,kk) > 258.) THEN |
---|
485 | het_rates(i,kk,het_HNO3) = MAX( 1. / yhno3, 0. ) * stay |
---|
486 | ELSE |
---|
487 | het_rates(i,kk,het_HNO3) = 0. |
---|
488 | ENDIF |
---|
489 | |
---|
490 | END DO |
---|
491 | END DO |
---|
492 | |
---|
493 | CASE(2) |
---|
494 | !----------------------------------------------------------------- |
---|
495 | ! ... PART 1, solve for high henry constant (HNO3). |
---|
496 | ! Convection only. |
---|
497 | ! The scavenging efficiency will be calculated in a |
---|
498 | ! following model version and directly included in |
---|
499 | ! the convective flux calculation according to Mari et al. |
---|
500 | ! Based on Guelle+Balkanski, Liu et al., and Mari et al. |
---|
501 | !----------------------------------------------------------------- |
---|
502 | |
---|
503 | scaveff = 0. |
---|
504 | alpha = 5.e-4 |
---|
505 | totmass = pdel / gravit |
---|
506 | |
---|
507 | DO i = 1, PLON |
---|
508 | het_rates(i,:,het_HNO3) = 0. |
---|
509 | IF (cldbot(i)==0 .OR. cldtop(i)==0 .OR. cldbot(i)==cldtop(i)) CYCLE |
---|
510 | |
---|
511 | DO k = cldbot(i), cldtop(i) |
---|
512 | dz = (zmid(i,k) - zmid(i,cldbot(i)))*1.e3 |
---|
513 | scaveff(i,k) = 1.-exp(-alpha*dz) |
---|
514 | ENDDO |
---|
515 | |
---|
516 | DO k = cldtop(i), cldbot(i)+1, -1 |
---|
517 | scaveff(i,k) = scaveff(i,k) - scaveff(i,k-1) |
---|
518 | ENDDO |
---|
519 | |
---|
520 | DO k = cldbot(i), cldtop(i) |
---|
521 | massupd(i,k) = MIN(flxupd(i,k)*delt,scaveff(i,k)*totmass(i,k)) |
---|
522 | het_rates(i,k,het_HNO3) = MIN(massupd(i,k)/totmass(i,k), scaveff(i,k)) |
---|
523 | ENDDO |
---|
524 | |
---|
525 | ENDDO |
---|
526 | |
---|
527 | |
---|
528 | WHERE (flxdwn .GT. 0.) |
---|
529 | het_rates(:,:,het_HNO3) = MAX (1.e-29,het_rates(:,:,het_HNO3)/delt) |
---|
530 | ELSEWHERE |
---|
531 | het_rates(:,:,het_HNO3) = 0. |
---|
532 | ENDWHERE |
---|
533 | |
---|
534 | END SELECT |
---|
535 | |
---|
536 | !----------------------------------------------------------------- |
---|
537 | ! ... PART 2, solve other low henry constant. |
---|
538 | ! Seinfeld and Pandis, (6.9) |
---|
539 | !----------------------------------------------------------------- |
---|
540 | #ifndef DUSS |
---|
541 | DO k = 1,PLEV |
---|
542 | |
---|
543 | work1(:) = Rg * tfld(:,k) * (lwc(index)*1.e-6) |
---|
544 | |
---|
545 | DO i = 1, HETCNT |
---|
546 | if ((hetname(i) .eq. "HNO3") .or. (hetname(i) .eq. "PB210")) cycle |
---|
547 | work2(:) = work1(:) * xhenconst(:,k,i) |
---|
548 | work2(:) = work2(:) / (1.+work2(:)) |
---|
549 | het_rates(:,k,i) = het_rates(:,k,het_HNO3) * work2(:) |
---|
550 | END DO |
---|
551 | |
---|
552 | END DO |
---|
553 | |
---|
554 | |
---|
555 | ! For Pb210: assume rate as HNO3 |
---|
556 | ! always in last slot |
---|
557 | het_rates(:,:,het_PB210) = het_rates(:,:,het_HNO3) |
---|
558 | |
---|
559 | #if defined(NMHC) && defined(AER) |
---|
560 | ! SO4 NH4 NO3 CSNO3M CINO3M (use the inclfra coefficients from aerosol_mod.F) |
---|
561 | het_rates(:,:,het_CINO3M) = het_rates(:,:,het_HNO3) * 0.80 |
---|
562 | het_rates(:,:,het_CSNO3M) = het_rates(:,:,het_HNO3) * 0.90 |
---|
563 | het_rates(:,:,het_ASNO3M) = het_rates(:,:,het_HNO3) * 0.80 |
---|
564 | het_rates(:,:,het_ASNH4M) = het_rates(:,:,het_HNO3) * 0.90 |
---|
565 | het_rates(:,:,het_ASSO4M) = het_rates(:,:,het_HNO3) * 0.80 |
---|
566 | #endif |
---|
567 | |
---|
568 | |
---|
569 | #endif |
---|
570 | #endif |
---|
571 | END SUBROUTINE SETHET |
---|