source: trunk/LMDZ.TITAN/libf/phytitan/calc_ysat.F90 @ 2236

Last change on this file since 2236 was 1966, checked in by jvatant, 7 years ago

Improve computation of saturation profiles that was done in planetary average
-> Now done for every grid point ( enables sensitivity to winter low temps )
-> firstcall "inicondens.F90" replaced by every-step-call-to "calc_ysat.F90"
--JVO
File size: 6.6 KB
Line 
1SUBROUTINE calc_ysat(nlon,nlay,press,temp,ysat)
2
3  !     ==============================================================================
4  !     Purpose
5  !     -------
6  !     Compute saturation profiles (mol/mol) for chemical tracers
7  !     
8  !     Authors
9  !     -------
10  !     S. Lebonnois
11  !     -> inicondens.F in old Titan, with T,P in planetary average
12  !     J. Vatant d'Ollone (2017)
13  !     -> Adapt to new physics, move to F90 and compute every grid point
14  !     ==============================================================================
15
16
17  !-----------------------------------------------------------------------
18  !     Declarations:
19  !     -------------
20
21  USE comchem_h, only: nkim, cnames
22
23  IMPLICIT NONE
24
25  !     Arguments :
26  !     -----------
27  INTEGER, INTENT(IN)                           :: nlon  ! # of grid points in the chunk
28  INTEGER, INTENT(IN)                           :: nlay  ! # of vertical layes
29  REAL, DIMENSION(nlon,nlay),      INTENT(IN)   :: press ! Mid-layers pressure (Pa)
30  REAL, DIMENSION(nlon,nlay),      INTENT(IN)   :: temp  ! Mid-layers temperature (K)
31  REAL, DIMENSION(nlon,nlay,nkim), INTENT(OUT)  :: ysat  ! Saturation profiles (mol/mol)
32
33  !     Local variables :
34  !     -----------------
35  INTEGER                           ::  ic
36  REAL,DIMENSION(:,:), ALLOCATABLE  ::  x
37  !     -------------------------------------------------------------------------
38
39  ALLOCATE(x(nlon,nlay))
40
41  !     By default, ysat=1, meaning we do not condense
42  ysat(:,:,:) = 1.0
43
44  !     Main loop
45
46  do ic=1,nkim
47
48     if(trim(cnames(ic)).eq."CH4") then
49
50        where (temp(:,:).lt.90.65)                                     
51           ysat(:,:,ic) = 10.0**(4.42507e0 - ( ( ( 1165560.7e0 / temp(:,:)                  &
52                -  115352.19e0 ) / temp(:,:) + 4055.6016e0 ) / temp(:,:)                    &
53                + 453.92414e0 ) / temp(:,:) ) / press(:,:) * 1013.25e0
54        elsewhere
55           ysat(:,:,ic) = 10.0**(3.901408e0 - ( ( 154567.02e0 / temp(:,:) - 1598.8512e0 )   &
56                / temp(:,:) + 437.54809e0 ) / temp(:,:)) / press(:,:) * 1013.25e0
57        endwhere
58
59        ! Forcing CH4 to 1.4% minimum               
60        where (ysat(:,:,ic).lt.0.014) ysat(:,:,ic) = 0.014
61
62     else if(trim(cnames(ic)).eq."C2H2") then
63
64        ysat(:,:,ic) = 10.0**(6.09748e0-1644.1e0/temp(:,:)+7.42346e0                        &
65             * alog10(1.0e3/temp(:,:)) ) / press(:,:)*1013.25e0/760.0
66
67     else if(trim(cnames(ic)).eq."C2H4") then
68
69        where (temp(:,:).lt.89.0)
70           ysat(:,:,ic) = 10.0**(1.5477e0 + (1.0e0/temp(:,:) - 0.011e0)                     &
71                * (16537.0e0*(1.0e0/temp(:,:) - 0.011e0) - 1038.1e0))                       &
72                / press(:,:) * 1.01325e0 / 760.0
73        elsewhere (temp(:,:).lt.104.0)
74           ysat(:,:,ic) = 10.0**(8.724e0 - 901.6e0/(temp(:,:) - 2.555e0) )                  &
75                / press(:,:) * 1013.25e0 / 760.0
76        elsewhere (temp(:,:).lt.120.0)
77           ysat(:,:,ic) = 10.0**(50.79e0 - 1703.0e0/temp(:,:) - 17.141e0                    &
78                * alog10(temp(:,:)) ) / press(:,:) * 1013.25e0 / 760.0
79        elsewhere (temp(:,:).lt.155.0)
80           ysat(:,:,ic) = 10.0**(6.74756e0 - 585.0e0/(temp(:,:) - 18.16e0) )                &
81                / press(:,:) * 1013.25e0 / 760.0
82        endwhere
83
84     else if(trim(cnames(ic)).eq."C2H6") then
85
86        where (temp(:,:).lt.90.)
87           ysat(:,:,ic) = 10.0**(10.01e0-1085.0e0/(temp(:,:)-0.561e0) )                     &
88                / press(:,:) * 1013.25e0 / 760.0e0
89        elsewhere
90           ysat(:,:,ic) = 10.0**(5.9366e0 - 1086.17e0/temp(:,:) + 3.83464e0                 &
91                * alog10(1.0e3/temp(:,:)) ) / press(:,:)*1013.25e0/760.0
92        endwhere
93
94     else if((trim(cnames(ic)).eq."CH3CCH") .or. (trim(cnames(ic)).eq."CH2CCH2")) then
95
96        ysat(:,:,ic) = 10.0**(2.8808e0 - 4.5e0*(249.9e0 - temp(:,:))                        &
97             / (1.15e0*temp(:,:) - 37.485e0) ) / press(:,:) * 1013.25e0 / 760.0e0
98
99     else if(trim(cnames(ic)).eq."C3H6")  then
100
101        ysat(:,:,ic) = 10.0**(7.4463e0 - 1028.5654e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
102
103     else if(trim(cnames(ic)).eq."C3H8")  then
104
105        ysat(:,:,ic) = 10.0**(7.217e0 - 994.30251e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
106
107     else if((trim(cnames(ic)).eq."C4H2").or.(trim(cnames(ic)).eq."C4H2s")) then
108
109        ysat(:,:,ic) = 10.0**(96.26781e0 - 4651.872e0/temp(:,:) - 31.68595e0                &
110             * alog10(temp(:,:)) ) / press(:,:) * 1013.25e0 / 760.0e0
111
112     else if(trim(cnames(ic)).eq."C4H4")  then
113
114        ysat(:,:,ic)= 1.0e3 * exp(9.3898e0 - 2203.57/(temp(:,:)-43.15e0) ) / press(:,:)
115
116     else if(trim(cnames(ic)).eq."C4H6")  then
117
118        ysat(:,:,ic)= 10.0**(2.8808e0 - 4.6e0*(262.3e0 - temp(:,:))                         &
119             /(1.15e0*temp(:,:) - 39.345e0) ) / press(:,:) * 1013.25e0 / 760.0e0
120
121     else if(trim(cnames(ic)).eq."C4H10")  then
122
123        ysat(:,:,ic)= 10.0**(8.446e0 - 1461.2e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
124
125     else if(trim(cnames(ic)).eq."C6H2")  then
126
127        ysat(:,:,ic)= 10.0**(4.666e0 - 4956e0/temp(:,:) + 25.845e0 *                        &
128             alog10(1.0e3/temp(:,:)) ) / press(:,:) * 1013.25e0 / 760.0e0
129
130     else if(trim(cnames(ic)).eq."C8H2")  then
131
132        ysat(:,:,ic)= 10.0**(3.95e0 - 6613e0/temp(:,:) + 35.055e0 *                         &
133             alog10(1.0e3/temp(:,:)) ) / press(:,:) * 1013.25e0 / 760.0e0
134
135     else if(trim(cnames(ic)).eq."AC6H6")  then
136
137        x = 1.0e0 - temp(:,:) / 562.2e0
138        ysat(:,:,ic)= 48.9e3 * exp( ( 1.33213 * x**1.5 - 6.98273 * x                        &
139             - x**3 * (2.62863 + 3.33399 * x**3) ) * 562.2e0/temp(:,:) ) / press(:,:)
140
141     else if(trim(cnames(ic)).eq."HCN")  then
142
143        ysat(:,:,ic)= 10.0**(8.6165e0 - 1516.5e0/(temp(:,:) - 26.2e0) ) / press(:,:) * 1013.25e0 / 760.0e0
144
145     else if(trim(cnames(ic)).eq."CH3CN")  then
146
147        ysat(:,:,ic)= 10.0**(8.458e0 - 1911.7e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
148
149     else if(trim(cnames(ic)).eq."C2H3CN")  then
150
151        ysat(:,:,ic)= 10.0**(9.3051e0 - 2782.21/(temp(:,:) - 51.15e0) ) / press(:,:) * 1013.25e0 / 760.0e0
152
153     else if(trim(cnames(ic)).eq."NCCN")  then
154
155        ysat(:,:,ic)=  10.0**(7.454e0 - 1832e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
156
157     else if(trim(cnames(ic)).eq."HC3N")  then
158
159        ysat(:,:,ic)= 10.0**(7.7446e0 - 1453.5609e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
160
161     else if(trim(cnames(ic)).eq."C4N2")  then
162
163        ysat(:,:,ic)= 10.0**(8.269e0 - 2155.0e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
164
165     endif
166  enddo
167
168END SUBROUTINE calc_ysat
Note: See TracBrowser for help on using the repository browser.