Context Navigation

sulfate_aer_mod.F90 @ 4750

Last change on this file since 4750 was 4750, checked in by dcugnet, 7 months ago
StratAer? : new strat composition (from reprobus) and density (Tabazadeh, 1994) routines.
Property svn:keywords set to `Id`
File size: 20.0 KB

Line
1	MODULE sulfate_aer_mod
2
3	! microphysical routines based on UPMC aerosol model by Slimane Bekki
4	! adapted for stratospheric sulfate aerosol in LMDZ by Christoph Kleinschmitt
5
6	CONTAINS
7
8	!*******************************************************************
9	SUBROUTINE STRACOMP_BIN(sh,t_seri,pplay)
10	!
11	! Aerosol H2SO4 weight fraction as a function of PH2O and temperature
12	! INPUT:
13	! sh: VMR of H2O
14	! t_seri: temperature (K)
15	! pplay: middle layer pression (Pa)
16	!
17	! OUTPUT:
18	! R2SO4: aerosol H2SO4 weight fraction (percent)
19
20	USE dimphy, ONLY : klon,klev ! nb of longitude and altitude bands
21	USE aerophys
22	USE phys_local_var_mod, ONLY: R2SO4
23
24	IMPLICIT NONE
25
26	REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature
27	REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression in the middle of each layer (Pa)
28	REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! specific humidity
29
30	REAL ks(7)
31	REAL t,qh2o,ptot,pw
32	REAL a,b,c,det
33	REAL xsb,msb
34
35	INTEGER ilon,ilev
36	DATA ks/-21.661,2724.2,51.81,-15732.0,47.004,-6969.0,-4.6183/
37
38	!*******************************************************************
39	!*** liquid aerosols process
40	!*******************************************************************
41	! BINARIES LIQUID AEROROLS:
42
43	DO ilon=1,klon
44	DO ilev=1,klev
45
46	t = max(t_seri(ilon,ilev),185.)
47	qh2o=sh(ilon,ilev)/18.*28.9
48	ptot=pplay(ilon,ilev)/100.
49	pw = qh2o*ptot/1013.0
50	pw = min(pw,2.e-3/1013.)
51	pw = max(pw,2.e-5/1013.)
52
53	!*******************************************************************
54	!*** binaries aerosols h2so4/h2o
55	!*******************************************************************
56	a = ks(3) + ks(4)/t
57	b = ks(1) + ks(2)/t
58	c = ks(5) + ks(6)/t + ks(7)*log(t) - log(pw)
59
60	det = b*2 - 4.a*c
61
62	IF (det > 0.) THEN
63	xsb = (-b - sqrt(det))/(2.*a)
64	msb = 55.51*xsb/(1.0 - xsb)
65	ELSE
66	msb = 0.
67	ENDIF
68	R2SO4(ilon,ilev) = 100msb0.098076/(1.0 + msb*0.098076)
69
70	! H2SO4 min dilution: 0.5%
71	R2SO4(ilon,ilev) = max( R2SO4(ilon,ilev), 0.005 )
72	ENDDO
73	ENDDO
74	100 RETURN
75
76	END SUBROUTINE STRACOMP_BIN
77
78	!********************************************************************
79	SUBROUTINE STRACOMP(sh,t_seri,pplay)
80
81	! AEROSOL H2SO4 WEIGHT FRACTION AS A FUNCTION OF PH2O AND TEMPERATURE
82	! ----------------------------------------------------------------
83	! INPUT:
84	! H2O: VMR of H2O
85	! t_seri: temperature (K)
86	! PMB: pressure (mb)
87	! klon: number of latitude bands in the model domain
88	! klev: number of altitude bands in the model domain
89	! for IFS: perhaps add another dimension for longitude
90	!
91	! OUTPUT:
92	! R2SO4: aerosol H2SO4 weight fraction (percent)
93
94	USE dimphy, ONLY : klon,klev
95	USE aerophys
96	USE phys_local_var_mod, ONLY: R2SO4
97
98	IMPLICIT NONE
99
100	REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature
101	REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa)
102	REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique
103
104	REAL PMB(klon,klev), H2O(klon,klev)
105	!
106	! working variables
107	INTEGER I,J,K
108	REAL TP, PH2O, VAL, A, B
109	! local variables to be saved on exit
110	INTEGER INSTEP
111	INTEGER, PARAMETER :: N=16, M=28
112	DATA INSTEP/0/
113	REAL F(N,M)
114	REAL XC(N)
115	REAL YC(M)
116	REAL XC1, XC16, YC1, YC28
117	!
118	SAVE INSTEP,F,XC,YC,XC1,XC16,YC1,YC28
119	!$OMP THREADPRIVATE(INSTEP,F,XC,YC,XC1,XC16,YC1,YC28)
120
121	! convert pplay (in Pa) to PMB (in mb)
122	PMB(:,:)=pplay(:,:)/100.0
123
124	! convert specific humidity sh (in kg/kg) to VMR H2O
125	H2O(:,:)=sh(:,:)*mAIRmol/mH2Omol
126
127	IF(INSTEP.EQ.0) THEN
128
129	INSTEP=1
130	XC(1)=0.01
131	XC(2)=0.1
132	XC(3)=0.5
133	XC(4)=1.0
134	XC(5)=1.5
135	XC(6)=2.0
136	XC(7)=3.0
137	XC(8)=5.0
138	XC(9)=6.0
139	XC(10)=8.0
140	XC(11)=10.0
141	XC(12)=12.0
142	XC(13)=15.0
143	XC(14)=20.0
144	XC(15)=30.0
145	XC(16)=100.0
146	!
147	YC(1)=175.0
148	DO I=2,28
149	YC(I)=YC(I-1)+5.0
150	ENDDO
151
152	! CONVERSION mb IN 1.0E-4mB
153	DO I=1,16
154	XC(I)=XC(I)*1.0E-4
155	ENDDO
156	!
157	XC1=XC(1)+1.E-10
158	XC16=XC(16)-1.E-8
159	YC1=YC(1)+1.E-5
160	YC28=YC(28)-1.E-5
161
162	F(6,4)=43.45
163	F(6,5)=53.96
164	F(6,6)=60.62
165	F(6,7)=65.57
166	F(6,8)=69.42
167	F(6,9)=72.56
168	F(6,10)=75.17
169	F(6,11)=77.38
170	F(6,12)=79.3
171	F(6,13)=80.99
172	F(6,14)=82.5
173	F(6,15)=83.92
174	F(6,16)=85.32
175	F(6,17)=86.79
176	F(6,18)=88.32
177	!
178	! ADD FACTOR BECAUSE THE SLOP IS TOO IMPORTANT
179	! NOT FOR THIS ONE BUT THE REST
180	! LOG DOESN'T WORK
181	A=(F(6,5)-F(6,4))/( (YC(5)-YC(4))*2.0)
182	B=-A*YC(4) + F(6,4)
183	F(6,1)=A*YC(1) + B
184	F(6,2)=A*YC(2) + B
185	F(6,3)=A*YC(3) + B
186	!
187	F(7,4)=37.02
188	F(7,5)=49.46
189	F(7,6)=57.51
190	F(7,7)=63.12
191	F(7,8)=67.42
192	F(7,9)=70.85
193	F(7,10)=73.70
194	F(7,11)=76.09
195	F(7,12)=78.15
196	F(7,13)=79.96
197	F(7,14)=81.56
198	F(7,15)=83.02
199	F(7,16)=84.43
200	F(7,17)=85.85
201	F(7,18)=87.33
202	!
203	A=(F(7,5)-F(7,4))/( (YC(5)-YC(4))*2.0)
204	B=-A*YC(4) + F(7,4)
205	F(7,1)=A*YC(1) + B
206	F(7,2)=A*YC(2) + B
207	F(7,3)=A*YC(3) + B
208	!
209	F(8,4)=25.85
210	F(8,5)=42.26
211	F(8,6)=52.78
212	F(8,7)=59.55
213	F(8,8)=64.55
214	F(8,9)=68.45
215	F(8,10)=71.63
216	F(8,11)=74.29
217	F(8,12)=76.56
218	F(8,13)=78.53
219	F(8,14)=80.27
220	F(8,15)=81.83
221	F(8,16)=83.27
222	F(8,17)=84.67
223	F(8,18)=86.10
224	!
225	A=(F(8,5)-F(8,4))/( (YC(5)-YC(4))*2.5 )
226	B=-A*YC(4) + F(8,4)
227	F(8,1)=A*YC(1) + B
228	F(8,2)=A*YC(2) + B
229	F(8,3)=A*YC(3) + B
230	!
231	F(9,4)=15.38
232	F(9,5)=39.35
233	F(9,6)=50.73
234	F(9,7)=58.11
235	F(9,8)=63.41
236	F(9,9)=67.52
237	F(9,10)=70.83
238	F(9,11)=73.6
239	F(9,12)=75.95
240	F(9,13)=77.98
241	F(9,14)=79.77
242	F(9,15)=81.38
243	F(9,16)=82.84
244	F(9,17)=84.25
245	F(9,18)=85.66
246	!
247	A=(F(9,5)-F(9,4))/( (YC(5)-YC(4))*7.0)
248	B=-A*YC(4) + F(9,4)
249	F(9,1)=A*YC(1) + B
250	F(9,2)=A*YC(2) + B
251	F(9,3)=A*YC(3) + B
252	!
253	F(10,4)=0.0
254	F(10,5)=34.02
255	F(10,6)=46.93
256	F(10,7)=55.61
257	F(10,8)=61.47
258	F(10,9)=65.94
259	F(10,10)=69.49
260	F(10,11)=72.44
261	F(10,12)=74.93
262	F(10,13)=77.08
263	F(10,14)=78.96
264	F(10,15)=80.63
265	F(10,16)=82.15
266	F(10,17)=83.57
267	F(10,18)=84.97
268	!
269	A=(F(10,6)-F(10,5))/( (YC(6)-YC(5))*1.5)
270	B=-A*YC(5) + F(10,5)
271	F(10,1)=A*YC(1) + B
272	F(10,2)=A*YC(2) + B
273	F(10,3)=A*YC(3) + B
274	F(10,4)=A*YC(4) + B
275	!
276	F(11,4)=0.0
277	F(11,5)=29.02
278	F(11,6)=43.69
279	F(11,7)=53.44
280	F(11,8)=59.83
281	F(11,9)=64.62
282	F(11,10)=68.39
283	F(11,11)=71.48
284	F(11,12)=74.10
285	F(11,13)=76.33
286	F(11,14)=78.29
287	F(11,15)=80.02
288	F(11,16)=81.58
289	F(11,17)=83.03
290	F(11,18)=84.44
291	!
292	A=(F(11,6)-F(11,5))/( (YC(6)-YC(5))*2.5 )
293	B=-A*YC(5) + F(11,5)
294	F(11,1)=A*YC(1) + B
295	F(11,2)=A*YC(2) + B
296	F(11,3)=A*YC(3) + B
297	F(11,4)=A*YC(4) + B
298	!
299	F(12,4)=0.0
300	F(12,5)=23.13
301	F(12,6)=40.86
302	F(12,7)=51.44
303	F(12,8)=58.38
304	F(12,9)=63.47
305	F(12,10)=67.43
306	F(12,11)=70.66
307	F(12,12)=73.38
308	F(12,13)=75.70
309	F(12,14)=77.72
310	F(12,15)=79.51
311	F(12,16)=81.11
312	F(12,17)=82.58
313	F(12,18)=83.99
314	!
315	A=(F(12,6)-F(12,5))/( (YC(6)-YC(5))*3.5 )
316	B=-A*YC(5) + F(12,5)
317	F(12,1)=A*YC(1) + B
318	F(12,2)=A*YC(2) + B
319	F(12,3)=A*YC(3) + B
320	F(12,4)=A*YC(4) + B
321	!
322	F(13,4)=0.0
323	F(13,5)=0.0
324	F(13,6)=36.89
325	F(13,7)=48.63
326	F(13,8)=56.46
327	F(13,9)=61.96
328	F(13,10)=66.19
329	F(13,11)=69.6
330	F(13,12)=72.45
331	F(13,13)=74.89
332	F(13,14)=76.99
333	F(13,15)=78.85
334	F(13,16)=80.50
335	F(13,17)=82.02
336	F(13,18)=83.44
337	!
338	A=(F(13,7)-F(13,6))/( (YC(7)-YC(6))*2.0)
339	B=-A*YC(6) + F(13,6)
340	F(13,1)=A*YC(1) + B
341	F(13,2)=A*YC(2) + B
342	F(13,3)=A*YC(3) + B
343	F(13,4)=A*YC(4) + B
344	F(13,5)=A*YC(5) + B
345	!
346	F(14,4)=0.0
347	F(14,5)=0.0
348	F(14,6)=30.82
349	F(14,7)=44.49
350	F(14,8)=53.69
351	F(14,9)=59.83
352	F(14,10)=64.47
353	F(14,11)=68.15
354	F(14,12)=71.19
355	F(14,13)=73.77
356	F(14,14)=76.0
357	F(14,15)=77.95
358	F(14,16)=79.69
359	F(14,17)=81.26
360	F(14,18)=82.72
361	!
362	A=(F(14,7)-F(14,6))/( (YC(7)-YC(6))*2.5 )
363	B=-A*YC(6) + F(14,6)
364	F(14,1)=A*YC(1) + B
365	F(14,2)=A*YC(2) + B
366	F(14,3)=A*YC(3) + B
367	F(14,4)=A*YC(4) + B
368	F(14,5)=A*YC(5) + B
369	!
370	F(15,4)=0.0
371	F(15,5)=0.0
372	F(15,6)=0.0
373	F(15,7)=37.71
374	F(15,8)=48.49
375	F(15,9)=56.40
376	F(15,10)=61.75
377	F(15,11)=65.89
378	F(15,12)=69.25
379	F(15,13)=72.07
380	F(15,14)=74.49
381	F(15,15)=76.59
382	F(15,16)=78.45
383	F(15,17)=80.12
384	F(15,18)=81.64
385	!
386	A=(F(15,8)-F(15,7))/( (YC(8)-YC(7))*1.5)
387	B=-A*YC(7) + F(15,7)
388	F(15,1)=A*YC(1) + B
389	F(15,2)=A*YC(2) + B
390	F(15,3)=A*YC(3) + B
391	F(15,4)=A*YC(4) + B
392	F(15,5)=A*YC(5) + B
393	F(15,6)=A*YC(6) + B
394
395	! SUPPOSE THAT AT GIVEN AND PH2O<2mB,
396	! %H2SO4 = A *LOG(PH2O) +B
397	! XC(1-5) :EXTENSION LEFT (LOW H2O)
398	DO J=1,18
399	A=(F(6,J)-F(7,J))/(LOG(XC(6))-LOG(XC(7)))
400	B=-A*LOG(XC(6)) + F(6,J)
401	DO K=1,5
402	F(K,J)=A*LOG(XC(K)) + B
403	ENDDO
404	ENDDO
405
406	! XC(16) :EXTENSION RIGHT (HIGH H2O)
407	DO J=1,18
408	A=(F(15,J)-F(14,J))/(XC(15)-XC(14))
409	B=-A*XC(15) + F(15,J)
410	F(16,J)=A*XC(16) + B
411	! F(16,2)=1.0
412	ENDDO
413
414	! YC(16-25) :EXTENSION DOWN (HIGH T)
415	DO I=1,16
416	A=(F(I,18)-F(I,17))/(YC(18)-YC(17))
417	B=-A*YC(18) + F(I,18)
418	DO K=19,28
419	F(I,K)=A*YC(K) + B
420	ENDDO
421	ENDDO
422
423	! MANUAL CORRECTIONS
424	DO J=1,10
425	F(1,J)=94.0
426	ENDDO
427
428	DO J=1,6
429	F(2,J)=77.0 +REAL(J)
430	ENDDO
431
432	DO J=1,7
433	F(16,J)=9.0
434	ENDDO
435
436	DO I=1,16
437	DO J=1,28
438	IF (F(I,J).LT.9.0) F(I,J)=30.0
439	IF (F(I,J).GT.99.99) F(I,J)=99.99
440	ENDDO
441	ENDDO
442
443	ENDIF
444
445	DO I=1,klon
446	DO J=1,klev
447	TP=t_seri(I,J)
448	IF (TP.LT.175.1) TP=175.1
449	! Partial pressure of H2O (mb)
450	PH2O =PMB(I,J)*H2O(I,J)
451	IF (PH2O.LT.XC1) THEN
452	R2SO4(I,J)=99.99
453	! PH2O=XC(1)+1.0E-10
454	ELSE
455	IF (PH2O.GT.XC16) PH2O=XC16
456	! SIMPLE LINEAR INTERPOLATIONS
457	CALL FIND(PH2O,TP,XC,YC,F,VAL,N,M)
458	IF (PMB(I,J).GE.10.0.AND.VAL.LT.60.0) VAL=60.0
459	R2SO4(I,J)=VAL
460	ENDIF
461	ENDDO
462	ENDDO
463
464	END SUBROUTINE
465
466	!****************************************************************
467	SUBROUTINE STRAACT(ACTSO4)
468
469	! H2SO4 ACTIVITY (GIAUQUE) AS A FUNCTION OF H2SO4 WP
470	! ----------------------------------------
471	! INPUT:
472	! H2SO4: VMR of H2SO4
473	! klon: number of latitude bands in the model domain
474	! klev: number of altitude bands in the model domain
475	! for IFS: perhaps add another dimension for longitude
476	!
477	! OUTPUT:
478	! ACTSO4: H2SO4 activity (percent)
479
480	USE dimphy, ONLY : klon,klev
481	USE phys_local_var_mod, ONLY: R2SO4
482
483	IMPLICIT NONE
484
485	REAL ACTSO4(klon,klev)
486
487	! Working variables
488	INTEGER NN,I,J,JX,JX1
489	REAL TC,TB,TA,XT
490	PARAMETER (NN=109)
491	REAL XC(NN), X(NN)
492
493	! H2SO4 activity
494	DATA X/ &
495	& 0.0,0.25,0.78,1.437,2.19,3.07,4.03,5.04,6.08 &
496	& ,7.13,8.18,14.33,18.59,28.59,39.17,49.49 &
497	& ,102.4,157.8,215.7,276.9,341.6,409.8,481.5,556.6 &
498	& ,635.5,719.,808.,902.,1000.,1103.,1211.,1322.,1437.,1555. &
499	& ,1677.,1800.,1926.,2054.,2183.,2312.,2442.,2572.,2701.,2829. &
500	& ,2955.,3080.,3203.,3325.,3446.,3564.,3681.,3796.,3910.,4022. &
501	& ,4134.,4351.,4564.,4771.,4974.,5171.,5364.,5551.,5732.,5908. &
502	& ,6079.,6244.,6404.,6559.,6709.,6854.,6994.,7131.,7264.,7393. &
503	& ,7520.,7821.,8105.,8373.,8627.,8867.,9093.,9308.,9511.,9703. &
504	& ,9885.,10060.,10225.,10535.,10819.,11079.,11318.,11537. &
505	& ,11740.,12097.,12407.,12676.,12915.,13126.,13564.,13910. &
506	& ,14191.,14423.,14617.,14786.,10568.,15299.,15491.,15654. &
507	& ,15811./
508	! H2SO4 weight fraction (percent)
509	DATA XC/ &
510	& 100.0,99.982,99.963,99.945,99.927,99.908,99.890,99.872 &
511	& ,99.853,99.835,99.817,99.725,99.634,99.452,99.270 &
512	& ,99.090,98.196,97.319,96.457,95.610,94.777,93.959,93.156 &
513	& ,92.365,91.588,90.824,90.073,89.334,88.607,87.892,87.188 &
514	& ,86.495,85.814,85.143,84.482,83.832,83.191,82.560,81.939 &
515	& ,81.327,80.724,80.130,79.545,78.968,78.399,77.839,77.286 &
516	& ,76.741,76.204,75.675,75.152,74.637,74.129,73.628,73.133 &
517	& ,72.164,71.220,70.300,69.404,68.530,67.678,66.847,66.037 &
518	& ,65.245,64.472,63.718,62.981,62.261,61.557,60.868,60.195 &
519	& ,59.537,58.893,58.263,57.646,56.159,54.747,53.405,52.126 &
520	& ,50.908,49.745,48.634,47.572,46.555,45.580,44.646,43.749 &
521	& ,42.059,40.495,39.043,37.691,36.430,35.251,33.107,31.209 &
522	& ,29.517,27.999,26.629,23.728,21.397,19.482,17.882,16.525 &
523	& ,15.360,13.461,11.980,10.792,9.819,8.932/
524
525	DO I=1,klon
526	DO J=1,klev
527	! HERE LINEAR INTERPOLATIONS
528	XT=R2SO4(I,J)
529	CALL POSACT(XT,XC,NN,JX)
530	JX1=JX+1
531	IF(JX.EQ.0) THEN
532	ACTSO4(I,J)=0.0
533	ELSE IF(JX.GE.NN) THEN
534	ACTSO4(I,J)=15811.0
535	ELSE
536	TC=XT -XC(JX)
537	TB=X(JX1) -X(JX)
538	TA=XC(JX1) -XC(JX)
539	TA=TB/TA
540	ACTSO4(I,J)=X(JX) + TA*TC
541	ENDIF
542	ENDDO
543	ENDDO
544
545	END SUBROUTINE
546
547	!****************************************************************
548	SUBROUTINE DENH2SA_TABA(t_seri)
549
550	! AERSOL DENSITY AS A FUNCTION OF H2SO4 WEIGHT PERCENT AND T
551	! from Tabazadeh et al. (1994) abaques
552	! ---------------------------------------------
553
554	!
555	! INPUT:
556	! R2SO4: aerosol H2SO4 weight fraction (percent)
557	! t_seri: temperature (K)
558	! klon: number of latitude bands in the model domain
559	! klev: number of altitude bands in the model domain
560	! for IFS: perhaps add another dimension for longitude
561	!
562	! OUTPUT:
563	! DENSO4: aerosol mass density (gr/cm3 = aerosol mass/aerosol volume)
564	!
565	USE dimphy, ONLY : klon,klev
566	USE phys_local_var_mod, ONLY: R2SO4, DENSO4
567
568	IMPLICIT NONE
569
570	REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature
571
572	INTEGER i,j
573
574	!----------------------------------------------------------------------
575	! ... Local variables
576	!----------------------------------------------------------------------
577	real, parameter :: a9 = -268.2616e4, a10 = 576.4288e3
578
579	real :: a0, a1, a2, a3, a4, a5, a6, a7 ,a8
580	real :: c1, c2, c3, c4, w
581
582
583	! Loop on model domain (2 dimension for UPMC model; 3 for IFS)
584	DO i=1,klon
585	DO j=1,klev
586	!----------------------------------------------------------------------
587	! ... Temperature variables
588	!----------------------------------------------------------------------
589	c1 = t_seri(I,J)- 273.15
590	c2 = c1**2
591	c3 = c1*c2
592	c4 = c1*c3
593	!----------------------------------------------------------------------
594	! Polynomial Coefficients
595	!----------------------------------------------------------------------
596	a0 = 999.8426 + 334.5402e-4c1 - 569.1304e-5c2
597	a1 = 547.2659 - 530.0445e-2c1 + 118.7671e-4c2 + 599.0008e-6*c3
598	a2 = 526.295e1 + 372.0445e-1c1 + 120.1909e-3c2 - 414.8594e-5c3 + 119.7973e-7c4
599	a3 = -621.3958e2 - 287.7670c1 - 406.4638e-3c2 + 111.9488e-4c3 + 360.7768e-7c4
600	a4 = 409.0293e3 + 127.0854e1c1 + 326.9710e-3c2 - 137.7435e-4c3 - 263.3585e-7c4
601	a5 = -159.6989e4 - 306.2836e1c1 + 136.6499e-3c2 + 637.3031e-5*c3
602	a6 = 385.7411e4 + 408.3717e1c1 - 192.7785e-3c2
603	a7 = -580.8064e4 - 284.4401e1*c1
604	a8 = 530.1976e4 + 809.1053*c1
605	!----------------------------------------------------------------------
606	! ... Summation
607	!----------------------------------------------------------------------
608	! w : H2SO4 Weight fraction
609	w=r2SO4(i,j)*0.01
610	DENSO4(i,j) = 0.001(a0 + w(a1 + w(a2 + w(a3 + w*(a4 + &
611	w(a5 + w(a6 + w(a7 + w(a8 + w(a9 + wa10))))))))))
612	DENSO4(i,j) = max (0.0, DENSO4(i,j) )
613
614	ENDDO
615	ENDDO
616
617	END SUBROUTINE DENH2SA_TABA
618
619	!****************************************************************
620	SUBROUTINE DENH2SA(t_seri)
621
622	! AERSOL DENSITY AS A FUNCTION OF H2SO4 WEIGHT PERCENT AND T
623	! ---------------------------------------------
624	! VERY ROUGH APPROXIMATION (SEE FOR WATER IN HANDBOOK
625	! LINEAR 2% FOR 30 DEGREES with RESPECT TO WATER)
626	!
627	! INPUT:
628	! R2SO4: aerosol H2SO4 weight fraction (percent)
629	! t_seri: temperature (K)
630	! klon: number of latitude bands in the model domain
631	! klev: number of altitude bands in the model domain
632	! for IFS: perhaps add another dimension for longitude
633	!
634	! OUTPUT:
635	! DENSO4: aerosol mass density (gr/cm3 = aerosol mass/aerosol volume)
636	!
637	USE dimphy, ONLY : klon,klev
638	USE phys_local_var_mod, ONLY: R2SO4, DENSO4
639
640	IMPLICIT NONE
641
642	REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature
643
644	INTEGER I,J
645
646	! Loop on model domain (2 dimension for UPMC model; 3 for IFS)
647	DO I=1,klon
648	DO J=1,klev
649	! RO AT 20C
650	DENSO4(I,J)=0.78681252E-5R2SO4(I,J)R2SO4(I,J)+ 0.82185978E-2*R2SO4(I,J)+0.97968381
651	DENSO4(I,J)=DENSO4(I,J)* ( 1.0 - (t_seri(I,J)-293.0)*0.02/30.0 )
652	ENDDO
653	ENDDO
654
655	END SUBROUTINE
656
657	!***********************************************************
658	SUBROUTINE FIND(X,Y,XC,YC,F,VAL,N,M)
659	!
660	! BI-LINEAR INTERPOLATION
661
662	! INPUT:
663	! X: Partial pressure of H2O (mb)
664	! Y: temperature (K)
665	! XC: Table partial pressure of H2O (mb)
666	! YC: Table temperature (K)
667	! F: Table aerosol H2SO4 weight fraction=f(XC,YC) (percent)
668	!
669	! OUTPUT:
670	! VAL: aerosol H2SO4 weight fraction (percent)
671
672	IMPLICIT NONE
673
674	INTEGER N,M
675	REAL X,Y,XC(N),YC(M),F(N,M),VAL
676	!
677	! working variables
678	INTEGER IERX,IERY,JX,JY,JXP1,JYP1
679	REAL SXY,SX1Y,SX1Y1,SXY1,TA,TB,T,UA,UB,U
680
681	IERX=0
682	IERY=0
683	CALL POSITION(XC,X,N,JX,IERX)
684	CALL POSITION(YC,Y,M,JY,IERY)
685
686	IF(JX.EQ.0.OR.IERY.EQ.1) THEN
687	VAL=99.99
688	RETURN
689	ENDIF
690
691	IF(JY.EQ.0.OR.IERX.EQ.1) THEN
692	VAL=9.0
693	RETURN
694	ENDIF
695
696	JXP1=JX+1
697	JYP1=JY+1
698	SXY=F(JX, JY )
699	SX1Y=F(JXP1,JY )
700	SX1Y1=F(JXP1,JYP1)
701	SXY1=F(JX, JYP1)
702
703	! x-slope.
704	TA=X -XC(JX)
705	TB=XC(JXP1)-XC(JX)
706	T=TA/TB
707
708	! y-slope.
709	UA=Y -YC(JY)
710	UB=YC(JYP1)-YC(JY)
711	U=UA/UB
712
713	! Use bilinear interpolation to determine function at point X,Y.
714	VAL=(1.-T)(1.-U)SXY + T(1.0-U)SX1Y + TUSX1Y1 + (1.0-T)USXY1
715
716	IF(VAL.LT.9.0) VAL=9.0
717	IF(VAL.GT.99.99) VAL=99.99
718
719	RETURN
720	END SUBROUTINE
721	!****************************************************************
722	SUBROUTINE POSITION(XC,X,N,JX,IER)
723
724	IMPLICIT NONE
725
726	INTEGER N,JX,IER,I
727	REAL X,XC(N)
728
729	IER=0
730	IF(X.LT.XC(1)) THEN
731	JX=0
732	ELSE
733	DO 10 I=1,N
734	IF (X.LT.XC(I)) GO TO 20
735	10 CONTINUE
736	IER=1
737	20 JX=I-1
738	ENDIF
739
740	RETURN
741	END SUBROUTINE
742	!********************************************************************
743	SUBROUTINE POSACT(XT,X,N,JX)
744
745	! POSITION OF XT IN THE ARRAY X
746	! -----------------------------------------------
747
748	IMPLICIT NONE
749
750	INTEGER N
751	REAL XT,X(N)
752	! Working variables
753	INTEGER JX,I
754
755	IF(XT.GT.X(1)) THEN
756	JX=0
757	ELSE
758	DO 10 I=1,N
759	IF (XT.GT.X(I)) GO TO 20
760	10 CONTINUE
761	20 JX=I
762	ENDIF
763
764	RETURN
765	END SUBROUTINE
766
767	END MODULE sulfate_aer_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format