My Project
 All Classes Files Functions Variables Macros
clouds_gno.F
Go to the documentation of this file.
1 !
2 ! $Header$
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, delta(klon), beta(klon)
50  real zu2,zv2
51  REAL xx(klon), aux(klon), coeff, block
52  REAL dist, fprime, det
53  REAL pi, u, v, erfcu, erfcv
54  REAL xx1, xx2
55  real erf,hsqrtlog_2,v2
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 !cdir arraycomb
61  cldf(1:klon,1:nd)=0.0 ! cym
62  ratqsc(1:klon,1:nd)=0.0
63  ptconv(1:klon,1:nd)=.false.
64 !cdir end arraycomb
65 
66  pi = acos(-1.)
67  sqrtpi=sqrt(pi)
68  sqrt2=sqrt(2.)
69  hsqrtlog_2=0.5*sqrt(log(2.))
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 = rs(i,k)
77  qsat = max(qsat,min_mu)
78  delta(i) = log(mu(i)/qsat)
79 c 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 c 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 = delta(i) + vmax(i)*vmax(i)
127  if (det.LE.0.0) vmax(i) = vmax0 + 1.0
128  det = delta(i) + vmax(i)*vmax(i)
129 
130  if (det.LE.0.) then
131  xx(i) = -0.0001
132  else
133  zx1=-sqrt2*vmax(i)
134  zx2=sqrt(1.0+delta(i)/(vmax(i)*vmax(i)))
135  xx1=zx1*(1.0-zx2)
136  xx2=zx1*(1.0+zx2)
137  xx(i) = 1.01 * xx1
138  if ( xx1 .GE. 0.0 ) xx(i) = 0.5*xx2
139  endif
140  if (delta(i).LT.0.) xx(i) = -hsqrtlog_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 = delta(i)/(xx(i)*sqrt2) + xx(i)/(2.*sqrt2)
156  v = delta(i)/(xx(i)*sqrt2) - xx(i)/(2.*sqrt2)
157  v2 = v*v
158 
159  IF ( v .GT. vmax(i) ) THEN
160 
161  IF ( abs(u) .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 = exp(-v*v) / v / sqrtpi
175  dist = 0.0
176  fprime = 1.0
177 
178  ELSE
179 
180 c -- erfv -> 1.0, use an asymptotic expression of erfv for v large:
181 
182  erfcu = 1.0-erf(u)
183 c !!! ATTENTION : rajout d'un seuil pour l'exponentiel
184  aux(i) = sqrtpi*erfcu*exp(min(v2,100.))
185  coeff = 1.0 - 0.5/(v2) + 0.75/(v2*v2)
186  block = coeff * exp(-v2) / v / sqrtpi
187  dist = v * aux(i) / coeff - beta(i)
188  fprime = 2.0 / xx(i) * (v2)
189  : * ( exp(-delta(i)) - u * aux(i) / coeff )
190  : / coeff
191 
192  ENDIF ! ABS(u)
193 
194  ELSE
195 
196 c -- general case:
197 
198  erfcu = 1.0-erf(u)
199  erfcv = 1.0-erf(v)
200  block = erfcv
201  dist = erfcu / erfcv - beta(i)
202  zu2=u*u
203  zv2=v2
204  if(zu2.gt.20..or. zv2.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=20.
211  zv2=20.
212  fprime = 0.
213  else
214  fprime = 2. /sqrtpi /xx(i) /(erfcv*erfcv)
215  : * ( erfcv*v*exp(-zu2)
216  : - erfcu*u*exp(-zv2) )
217  endif
218  ENDIF ! x
219 
220 c -- test numerical convergence:
221 
222 ! if (beta(i).lt.1.e-10) then
223 ! print*,'avant test ',i,k,lconv(i),u(i),v(i),beta(i)
224 ! stop
225 ! endif
226  if (abs(fprime).lt.1.e-11) then
227 ! print*,'avant test fprime<.e-11 '
228 ! s ,i,k,lconv(i),u(i),v(i),beta(i),fprime(i)
229 ! print*,'klon,ND,R,RS,QSUB',
230 ! s klon,ND,R(i,k),rs(i,k),qsub(i,k)
231  fprime=sign(1.e-11,fprime)
232  endif
233 
234 
235  if ( abs(dist/beta(i)) .LT. epsilon ) then
236 c print*,'v-u **2',(v(i)-u(i))**2
237 c print*,'exp v-u **2',exp((v(i)-u(i))**2)
238  ptconv(i,k) = .true.
239  lconv(i)=.true.
240 c borne pour l'exponentielle
241  ratqsc(i,k)=min(2.*(v-u)*(v-u),20.)
242  ratqsc(i,k)=sqrt(exp(ratqsc(i,k))-1.)
243  cldf(i,k) = 0.5 * block
244  else
245  xx(i) = xx(i) - dist/fprime
246  endif
247 c print*,'apres test ',i,k,lconv(i)
248 
249  endif ! lconv
250  enddo ! vector
251 
252 c----------------------------------------------------------------------
253 c Fin des nmax iterations pour trouver la solution.
254  ENDDO ! n
255 c----------------------------------------------------------------------
256 
257 500 CONTINUE ! K
258 
259  RETURN
260  END
261 
262 
263