Context Navigation

sulfate_aer_mod.F90 @ 4996

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

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

Download in other formats:

Original Format