module_hetn2o5.F
References to this file elsewhere.
1 MODULE module_hetn2o5
2
3
4 !temporary:
5 ! from dentener and crutzen for 5S, 160E on 100mb levels startin at 1000mb
6 REAL, PARAMETER, DIMENSION(10) :: &
7 kheti = (/4.2e-5, 3.8e-5, 1.e-5,3.e-6,4.e-6, 2.e-6,3.e-6,5.e-6,1.3e-5,1.8e-5/)
8
9 !kheta: on aerosol
10 REAL, SAVE, DIMENSION(120) :: kheta
11
12 LOGICAL :: do_aerosol=.TRUE.
13
14
15 ! sticking coefficients for cloud water and cloud ice
16
17 REAL , PARAMETER, PRIVATE :: gammacldw = 0.05, &
18 gammacldi = 0.03 !cms check !
19
20 REAL , PARAMETER, PRIVATE :: rhograul = 400. ,&
21 rhoice = 917.
22
23 REAL, PARAMETER, PRIVATE :: pi = 3.1415926
24
25
26 ! D_g: diffusivity of species in gas phase in m^2/s
27 REAL, PARAMETER, PRIVATE :: D_g = 0.1E-6
28 !in Match this is 1.e-5 (Schwartz ?!)
29
30
31
32
33 ! abar_c: mass mean radius for cloud drops in m
34 REAL, PARAMETER, PRIVATE :: abar_c = 10.E-6
35
36
37
38 ! RSTAR2 : universal gas constant in J/(kmol K)
39 REAL, PARAMETER, PRIVATE :: RSTAR2 = 8314.
40
41
42 ! parameters from Lin scheme
43 REAL , PARAMETER, PRIVATE :: xnor = 8.0e6
44 REAL , PARAMETER, PRIVATE :: xnos = 3.0e6
45 REAL , PARAMETER, PRIVATE :: xnog = 4.0e6
46
47
48
49 CONTAINS
50 SUBROUTINE hetn2o5calc(hetn2o5, rho, T, QC, QR, QI, QS, QG, &
51 rhowater, rhosnow, M_n2o5, &
52 P_QI,P_QS,P_QG, &
53 ids,ide, jds,jde, kds,kde, &
54 ims,ime, jms,jme, kms,kme, &
55 its,ite, jts,jte, kts,kte )
56
57
58
59 IMPLICIT NONE
60
61
62
63 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
64 ims,ime, jms,jme, kms,kme , &
65 its,ite, jts,jte, kts,kte , &
66 P_QI,P_QS,P_QG
67
68
69 REAL, DIMENSION( its:ite , kts:kte , jts:jte), &
70 INTENT(INOUT ) :: hetn2o5
71
72
73
74 REAL, INTENT(IN ) :: rhosnow, &
75 rhowater, &
76 M_n2o5
77
78
79 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
80 INTENT(IN ) :: &
81 QC, &
82 QR, &
83 QI, &
84 QS, &
85 QG, &
86 rho, &
87 T
88
89
90 !local vars
91
92 REAL :: orho, tmp1
93
94
95 ! mass mean radii of the different hydrometeors
96 REAL :: &
97 abar_r, abar_i, abar_s, abar_g
98
99
100 ! nu_th: thermal velocity of species in m/s
101 REAL :: nu_th
102
103 ! tau_Dg: gas diffusion timescale
104 ! tau_i : timescale for mass transfer across interface of hydrometeor
105 REAL :: tau_i, tau_Dg
106
107
108 ! L : water contents in cloud drops, rain drops, hail stones,...
109 ! in (cm^3 H_2O / cm^3 air)
110 REAL :: L
111
112
113 ! loss rates:
114 REAL :: kc, kr, ki, ks, kg
115
116
117 INTEGER :: i,j,k
118
119
120 !----
121
122 j_loop: DO j = jts, jte
123 i_loop: DO i = its, ite
124
125 DO k = kts, kte
126 orho = 1./rho(i,k,j)
127
128
129
130 kc = 0.
131 kr = 0.
132 ki = 0.
133 ks = 0.
134 kg = 0.
135
136
137
138 !! abar_c
139
140
141 IF(QR(i,k,j) .GT. 1.e-12) THEN
142 tmp1=sqrt(pi*rhowater*xnor*orho/QR(i,k,j))
143 !!xlambdar1Dr(k)=sqrt(tmp1)
144 abar_r = 2./MAX(sqrt(tmp1), 1.E-8 ) !if abar is large, kt becomes small
145 ENDIF
146
147
148 !!$ IF(QI(i,k,j) .GT. 1.e-8) THEN
149 !!$ tmp1=sqrt(pi*rhoice*xnoi*orho/QI(i,k,j))
150 !!$ !!xlambdar1Di(k)=sqrt(tmp1)
151 !!$ abar_i = 2./MAX(sqrt(tmp1), 1.E-8 )
152 !!$ ENDIF
153
154
155 abar_i = abar_c
156
157 IF(QS(i,k,j) .GT. 1.e-12) THEN
158 tmp1=sqrt(pi*rhosnow*xnos*orho/QS(i,k,j))
159 !!xlambdar1Ds(k)=sqrt(tmp1)
160 abar_s = 2./MAX(sqrt(tmp1), 1.E-8 )
161 ENDIF
162
163
164
165 IF(QG(i,k,j) .GT. 1.e-12) THEN
166 tmp1=sqrt(pi*rhograul*xnog*orho/QG(i,k,j))
167 !!xlambdar1Dg(k)=sqrt(tmp1)
168 abar_g = 2./MAX(sqrt(tmp1), 1.E-8 )
169 ENDIF
170
171
172
173
174 ! calculate thermal velocity of species
175 nu_th = SQRT(8*RSTAR2*T(i,k,j)/(pi*M_n2o5))
176
177
178
179 ! calculate timescales
180 !cloud droplets
181 IF(QC(i,k,j) .GT. 1.e-12) THEN
182 tau_i = (4.*abar_c) / (3.* nu_th * gammacldw)
183 tau_Dg = (abar_c**2) / (3.* D_g)
184
185 L = (rho(i,k,j) / rhowater) * QC(i,k,j)
186 kc = L/(tau_i + tau_Dg)
187 ENDIF
188
189 !rain
190 IF(QR(i,k,j) .GT. 1.e-12) THEN
191 tau_i = (4.*abar_r) / (3.* nu_th * gammacldw)
192 tau_Dg = (abar_r**2) / (3.* D_g)
193
194 L = (rho(i,k,j) / rhowater) * QR(i,k,j)
195 kr = L/(tau_i + tau_Dg)
196 ENDIF
197
198 !cloud ice
199 IF(QI(i,k,j) .GT. 1.e-12) THEN
200 tau_i = (4.*abar_i) / (3.* nu_th * gammacldi)
201 tau_Dg = (abar_i**2) / (3.* D_g)
202
203 L = (rho(i,k,j) / rhowater) * QI(i,k,j)
204 ki = L/(tau_i + tau_Dg)
205 ENDIF
206
207 !snow
208 IF(QS(i,k,j).GT. 1.e-12) THEN
209 tau_i = (4.*abar_s) / (3.* nu_th * gammacldi)
210 tau_Dg = (abar_s**2) / (3.* D_g)
211
212 L = (rho(i,k,j) / rhowater) * QS(i,k,j)
213 ks = L/(tau_i + tau_Dg)
214 ENDIF
215
216 !graupel
217 IF(QG(i,k,j).GT. 1.e-12) THEN
218 tau_i = (4.*abar_g) / (3.* nu_th * gammacldi)
219 tau_Dg = (abar_g**2) / (3.* D_g)
220
221 L = (rho(i,k,j) / rhowater) * QG(i,k,j)
222 kg = L/(tau_i + tau_Dg)
223 ENDIF
224
225
226
227 !!HERE THE VENTILATION COEFF SHOULD BE INCLUDED AND IT PROBABLY DOES NOT HAPPEN ON ICE!!
228
229 !! hetn2o5(i,k,j) = kc + kr + ki + ks + kg
230
231
232
233
234
235 hetn2o5(i,k,j) = kc + kr + kheta(k) !!+ kg
236 !hetn2o5(i,k,j) = 0. !n2o5_test
237
238
239 ENDDO
240
241 ENDDO i_loop
242 ENDDO j_loop
243
244 END SUBROUTINE hetn2o5calc
245
246
247
248
249
250 !-----------------------------------------------------------------
251
252 SUBROUTINE hetn2o5_ini( pb, pp, z, &
253 ids, ide, jds, jde, kds, kde, &
254 ims, ime, jms, jme, kms, kme, &
255 its, ite, jts, jte, kts, kte)
256
257 #ifdef DM_PARALLEL
258 USE module_dm
259 #endif
260
261 IMPLICIT NONE
262
263
264 INTEGER, INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
265 ims, ime, jms, jme, kms, kme, &
266 its, ite, jts, jte, kts, kte
267
268
269 REAL, DIMENSION( ims:ime, kms:kme, jms:jme), &
270 INTENT(IN ) :: z, pb, pp
271
272 REAL, DIMENSION(10) :: pin !moguntia
273
274
275 LOGICAL , EXTERNAL :: wrf_dm_on_monitor
276
277
278 !local
279
280 REAL, DIMENSION(kms:kme) :: zloc, ploc
281
282 INTEGER :: k, kk, l_low, l_up
283 REAL :: m, b, dp
284
285
286 #define DM_BCAST_MACRO(A) CALL wrf_dm_bcast_bytes ( A , size ( A ) * RWORDSIZE )
287 IF (.NOT. do_aerosol) RETURN
288
289 !currently not really necessary
290 dm_on_monitor: IF ( wrf_dm_on_monitor() ) THEN
291
292 kheta(:) = 0.
293
294
295
296 pin(kms)=1000.
297 DO k=2,10
298 pin(k)=pin(k-1)-100.
299 ENDDO
300
301
302 DO k=kms,kme
303
304 zloc(k) = z(its,k,jts)
305 ploc(k) = pb(its,k,jts) + pp(its,k,jts)
306
307 IF ( zloc(k) .GT. 1.e6 .OR. zloc(k) .LT. 0. ) THEN
308 CALL wrf_error_fatal ("STOP: hetn2o5_ini")
309 ENDIF
310
311 ENDDO
312
313 DO k=kms,kme
314 IF ( ploc(k) .GT. 1.e5 ) THEN
315 kheta(k) = kheti(1)
316 ELSE
317
318 l_low = 11.-CEILING(ploc(k)/1.e4)
319 l_up=l_low+1
320
321
322 dp=ploc(k)/100. - pin(l_up)
323 m=(kheti(l_low)-kheti(l_up))/100.
324 kheta(k)=kheti(l_up) + m*dp
325
326
327
328 !print *,l_low
329 !print *, ploc(k), pin(l_low), pin(l_up)
330 !print *, " ", kheti(l_low),kheti(l_up), kheta(k)
331
332
333
334 END IF
335
336
337 ENDDO
338
339
340
341
342 ENDIF dm_on_monitor
343
344 DM_BCAST_MACRO( kheta )
345
346
347
348 END SUBROUTINE hetn2o5_ini
349
350
351 END MODULE module_hetn2o5
352
353
354
355
356
357