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