source: LMDZ6/trunk/libf/phylmd/StratAer/sulfate_aer_mod.F90 @ 4744

Last change on this file since 4744 was 3663, checked in by oboucher, 5 years ago

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