1 | ! |
---|
2 | ! $Id: clouds_gno.F 1163 2009-05-20 14:11:21Z evignon $ |
---|
3 | ! |
---|
4 | C |
---|
5 | C================================================================================ |
---|
6 | C |
---|
7 | SUBROUTINE CLOUDS_GNO(klon,ND,R,RS,QSUB,PTCONV,RATQSC,CLDF) |
---|
8 | IMPLICIT NONE |
---|
9 | C |
---|
10 | C-------------------------------------------------------------------------------- |
---|
11 | C |
---|
12 | C Inputs: |
---|
13 | C |
---|
14 | C ND----------: Number of vertical levels |
---|
15 | C R--------ND-: Domain-averaged mixing ratio of total water |
---|
16 | C RS-------ND-: Mean saturation humidity mixing ratio within the gridbox |
---|
17 | C QSUB-----ND-: Mixing ratio of condensed water within clouds associated |
---|
18 | C with SUBGRID-SCALE condensation processes (here, it is |
---|
19 | C predicted by the convection scheme) |
---|
20 | C Outputs: |
---|
21 | C |
---|
22 | C PTCONV-----ND-: Point convectif = TRUE |
---|
23 | C RATQSC-----ND-: Largeur normalisee de la distribution |
---|
24 | C CLDF-----ND-: Fraction nuageuse |
---|
25 | C |
---|
26 | C-------------------------------------------------------------------------------- |
---|
27 | |
---|
28 | |
---|
29 | INTEGER klon,ND |
---|
30 | REAL R(klon,ND), RS(klon,ND), QSUB(klon,ND) |
---|
31 | LOGICAL PTCONV(klon,ND) |
---|
32 | REAL RATQSC(klon,ND) |
---|
33 | REAL CLDF(klon,ND) |
---|
34 | |
---|
35 | c -- parameters controlling the iteration: |
---|
36 | c -- nmax : maximum nb of iterations (hopefully never reached) |
---|
37 | c -- epsilon : accuracy of the numerical resolution |
---|
38 | c -- vmax : v-value above which we use an asymptotic expression for ERF(v) |
---|
39 | |
---|
40 | INTEGER nmax |
---|
41 | PARAMETER ( nmax = 10) |
---|
42 | REAL epsilon, vmax0, vmax(klon) |
---|
43 | PARAMETER ( epsilon = 0.02, vmax0 = 2.0 ) |
---|
44 | |
---|
45 | REAL min_mu, min_Q |
---|
46 | PARAMETER ( min_mu = 1.e-12, min_Q=1.e-12 ) |
---|
47 | |
---|
48 | INTEGER i,K, n, m |
---|
49 | REAL mu(klon), qsat(klon), delta(klon), beta(klon) |
---|
50 | real zu2(klon),zv2(klon) |
---|
51 | REAL xx(klon), aux(klon), coeff(klon), block(klon) |
---|
52 | REAL dist(klon), fprime(klon), det(klon) |
---|
53 | REAL pi, u(klon), v(klon), erfcu(klon), erfcv(klon) |
---|
54 | REAL xx1(klon), xx2(klon) |
---|
55 | real erf,kkk |
---|
56 | real sqrtpi,sqrt2,zx1,zx2,exdel |
---|
57 | c lconv = true si le calcul a converge (entre autre si qsub < min_q) |
---|
58 | LOGICAL lconv(klon) |
---|
59 | |
---|
60 | cym |
---|
61 | cldf(:,:)=0.0 |
---|
62 | |
---|
63 | pi = ACOS(-1.) |
---|
64 | sqrtpi=sqrt(pi) |
---|
65 | sqrt2=sqrt(2.) |
---|
66 | |
---|
67 | ptconv=.false. |
---|
68 | ratqsc=0. |
---|
69 | |
---|
70 | |
---|
71 | DO 500 K = 1, ND |
---|
72 | |
---|
73 | do i=1,klon ! vector |
---|
74 | mu(i) = R(i,K) |
---|
75 | mu(i) = MAX(mu(i),min_mu) |
---|
76 | qsat(i) = RS(i,K) |
---|
77 | qsat(i) = MAX(qsat(i),min_mu) |
---|
78 | delta(i) = log(mu(i)/qsat(i)) |
---|
79 | enddo ! vector |
---|
80 | |
---|
81 | C |
---|
82 | C *** There is no subgrid-scale condensation; *** |
---|
83 | C *** the scheme becomes equivalent to an "all-or-nothing" *** |
---|
84 | C *** large-scale condensation scheme. *** |
---|
85 | C |
---|
86 | |
---|
87 | C |
---|
88 | C *** Some condensation is produced at the subgrid-scale *** |
---|
89 | C *** *** |
---|
90 | C *** PDF = generalized log-normal distribution (GNO) *** |
---|
91 | C *** (k<0 because a lower bound is considered for the PDF) *** |
---|
92 | C *** *** |
---|
93 | C *** -> Determine x (the parameter k of the GNO PDF) such *** |
---|
94 | C *** that the contribution of subgrid-scale processes to *** |
---|
95 | C *** the in-cloud water content is equal to QSUB(K) *** |
---|
96 | C *** (equations (13), (14), (15) + Appendix B of the paper) *** |
---|
97 | C *** *** |
---|
98 | C *** Here, an iterative method is used for this purpose *** |
---|
99 | C *** (other numerical methods might be more efficient) *** |
---|
100 | C *** *** |
---|
101 | C *** NB: the "error function" is called ERF *** |
---|
102 | C *** (ERF in double precision) *** |
---|
103 | C |
---|
104 | |
---|
105 | c On commence par eliminer les cas pour lesquels on n'a pas |
---|
106 | c suffisamment d'eau nuageuse. |
---|
107 | |
---|
108 | do i=1,klon ! vector |
---|
109 | |
---|
110 | IF ( QSUB(i,K) .lt. min_Q ) THEN |
---|
111 | ptconv(i,k)=.false. |
---|
112 | ratqsc(i,k)=0. |
---|
113 | lconv(i) = .true. |
---|
114 | |
---|
115 | c Rien on a deja initialise |
---|
116 | |
---|
117 | ELSE |
---|
118 | |
---|
119 | lconv(i) = .FALSE. |
---|
120 | vmax(i) = vmax0 |
---|
121 | |
---|
122 | beta(i) = QSUB(i,K)/mu(i) + EXP( -MIN(0.0,delta(i)) ) |
---|
123 | |
---|
124 | c -- roots of equation v > vmax: |
---|
125 | |
---|
126 | det(i) = delta(i) + vmax(i)**2. |
---|
127 | if (det(i).LE.0.0) vmax(i) = vmax0 + 1.0 |
---|
128 | det(i) = delta(i) + vmax(i)**2. |
---|
129 | |
---|
130 | if (det(i).LE.0.) then |
---|
131 | xx(i) = -0.0001 |
---|
132 | else |
---|
133 | zx1=-sqrt2*vmax(i) |
---|
134 | zx2=SQRT(1.0+delta(i)/(vmax(i)**2.)) |
---|
135 | xx1(i)=zx1*(1.0-zx2) |
---|
136 | xx2(i)=zx1*(1.0+zx2) |
---|
137 | xx(i) = 1.01 * xx1(i) |
---|
138 | if ( xx1(i) .GE. 0.0 ) xx(i) = 0.5*xx2(i) |
---|
139 | endif |
---|
140 | if (delta(i).LT.0.) xx(i) = -0.5*SQRT(log(2.)) |
---|
141 | |
---|
142 | ENDIF |
---|
143 | |
---|
144 | enddo ! vector |
---|
145 | |
---|
146 | c---------------------------------------------------------------------- |
---|
147 | c Debut des nmax iterations pour trouver la solution. |
---|
148 | c---------------------------------------------------------------------- |
---|
149 | |
---|
150 | DO n = 1, nmax |
---|
151 | |
---|
152 | do i=1,klon ! vector |
---|
153 | if (.not.lconv(i)) then |
---|
154 | |
---|
155 | u(i) = delta(i)/(xx(i)*sqrt2) + xx(i)/(2.*sqrt2) |
---|
156 | v(i) = delta(i)/(xx(i)*sqrt2) - xx(i)/(2.*sqrt2) |
---|
157 | |
---|
158 | IF ( v(i) .GT. vmax(i) ) THEN |
---|
159 | |
---|
160 | IF ( ABS(u(i)) .GT. vmax(i) |
---|
161 | : .AND. delta(i) .LT. 0. ) THEN |
---|
162 | |
---|
163 | c -- use asymptotic expression of erf for u and v large: |
---|
164 | c ( -> analytic solution for xx ) |
---|
165 | exdel=beta(i)*EXP(delta(i)) |
---|
166 | aux(i) = 2.0*delta(i)*(1.-exdel) |
---|
167 | : /(1.+exdel) |
---|
168 | if (aux(i).lt.0.) then |
---|
169 | c print*,'AUX(',i,',',k,')<0',aux(i),delta(i),beta(i) |
---|
170 | aux(i)=0. |
---|
171 | endif |
---|
172 | xx(i) = -SQRT(aux(i)) |
---|
173 | block(i) = EXP(-v(i)*v(i)) / v(i) / sqrtpi |
---|
174 | dist(i) = 0.0 |
---|
175 | fprime(i) = 1.0 |
---|
176 | |
---|
177 | ELSE |
---|
178 | |
---|
179 | c -- erfv -> 1.0, use an asymptotic expression of erfv for v large: |
---|
180 | |
---|
181 | erfcu(i) = 1.0-ERF(u(i)) |
---|
182 | c !!! ATTENTION : rajout d'un seuil pour l'exponentiel |
---|
183 | aux(i) = sqrtpi*erfcu(i)*EXP(min(v(i)*v(i),100.)) |
---|
184 | coeff(i) = 1.0 - 1./2./(v(i)**2.) + 3./4./(v(i)**4.) |
---|
185 | block(i) = coeff(i) * EXP(-v(i)*v(i)) / v(i) / sqrtpi |
---|
186 | dist(i) = v(i) * aux(i) / coeff(i) - beta(i) |
---|
187 | fprime(i) = 2.0 / xx(i) * (v(i)**2.) |
---|
188 | : * ( coeff(i)*EXP(-delta(i)) - u(i) * aux(i) ) |
---|
189 | : / coeff(i) / coeff(i) |
---|
190 | |
---|
191 | ENDIF ! ABS(u) |
---|
192 | |
---|
193 | ELSE |
---|
194 | |
---|
195 | c -- general case: |
---|
196 | |
---|
197 | erfcu(i) = 1.0-ERF(u(i)) |
---|
198 | erfcv(i) = 1.0-ERF(v(i)) |
---|
199 | block(i) = erfcv(i) |
---|
200 | dist(i) = erfcu(i) / erfcv(i) - beta(i) |
---|
201 | zu2(i)=u(i)*u(i) |
---|
202 | zv2(i)=v(i)*v(i) |
---|
203 | if(zu2(i).gt.20..or. zv2(i).gt.20.) then |
---|
204 | c print*,'ATTENTION !!! xx(',i,') =', xx(i) |
---|
205 | c print*,'ATTENTION !!! klon,ND,R,RS,QSUB,PTCONV,RATQSC,CLDF', |
---|
206 | c .klon,ND,R(i,k),RS(i,k),QSUB(i,k),PTCONV(i,k),RATQSC(i,k), |
---|
207 | c .CLDF(i,k) |
---|
208 | c print*,'ATTENTION !!! zu2 zv2 =',zu2(i),zv2(i) |
---|
209 | zu2(i)=20. |
---|
210 | zv2(i)=20. |
---|
211 | fprime(i) = 0. |
---|
212 | else |
---|
213 | fprime(i) = 2. /sqrtpi /xx(i) /erfcv(i)**2. |
---|
214 | : * ( erfcv(i)*v(i)*EXP(-zu2(i)) |
---|
215 | : - erfcu(i)*u(i)*EXP(-zv2(i)) ) |
---|
216 | endif |
---|
217 | ENDIF ! x |
---|
218 | |
---|
219 | c -- test numerical convergence: |
---|
220 | |
---|
221 | ! if (beta(i).lt.1.e-10) then |
---|
222 | ! print*,'avant test ',i,k,lconv(i),u(i),v(i),beta(i) |
---|
223 | ! stop |
---|
224 | ! endif |
---|
225 | if (abs(fprime(i)).lt.1.e-11) then |
---|
226 | ! print*,'avant test fprime<.e-11 ' |
---|
227 | ! s ,i,k,lconv(i),u(i),v(i),beta(i),fprime(i) |
---|
228 | ! print*,'klon,ND,R,RS,QSUB', |
---|
229 | ! s klon,ND,R(i,k),rs(i,k),qsub(i,k) |
---|
230 | fprime(i)=sign(1.e-11,fprime(i)) |
---|
231 | endif |
---|
232 | |
---|
233 | |
---|
234 | if ( ABS(dist(i)/beta(i)) .LT. epsilon ) then |
---|
235 | c print*,'v-u **2',(v(i)-u(i))**2 |
---|
236 | c print*,'exp v-u **2',exp((v(i)-u(i))**2) |
---|
237 | ptconv(i,K) = .TRUE. |
---|
238 | lconv(i)=.true. |
---|
239 | c borne pour l'exponentielle |
---|
240 | ratqsc(i,k)=min(2.*(v(i)-u(i))**2,20.) |
---|
241 | ratqsc(i,k)=sqrt(exp(ratqsc(i,k))-1.) |
---|
242 | CLDF(i,K) = 0.5 * block(i) |
---|
243 | else |
---|
244 | xx(i) = xx(i) - dist(i)/fprime(i) |
---|
245 | endif |
---|
246 | c print*,'apres test ',i,k,lconv(i) |
---|
247 | |
---|
248 | endif ! lconv |
---|
249 | enddo ! vector |
---|
250 | |
---|
251 | c---------------------------------------------------------------------- |
---|
252 | c Fin des nmax iterations pour trouver la solution. |
---|
253 | ENDDO ! n |
---|
254 | c---------------------------------------------------------------------- |
---|
255 | |
---|
256 | 500 CONTINUE ! K |
---|
257 | |
---|
258 | RETURN |
---|
259 | END |
---|
260 | |
---|
261 | |
---|
262 | |
---|