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