source: LMDZ4/branches/pre_V3/libf/phylmd/clouds_gno.F @ 4082

Last change on this file since 4082 was 559, checked in by lmdzadmin, 20 years ago

Initialisations diverses YM
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.0 KB
Line 
1!
2! $Header$
3!
4C
5C================================================================================
6C
7      SUBROUTINE CLOUDS_GNO(klon,ND,R,RS,QSUB,PTCONV,RATQSC,CLDF)
8      IMPLICIT NONE
9C     
10C--------------------------------------------------------------------------------
11C
12C Inputs:
13C
14C  ND----------: Number of vertical levels
15C  R--------ND-: Domain-averaged mixing ratio of total water
16C  RS-------ND-: Mean saturation humidity mixing ratio within the gridbox
17C  QSUB-----ND-: Mixing ratio of condensed water within clouds associated
18C                with SUBGRID-SCALE condensation processes (here, it is
19C                predicted by the convection scheme)
20C Outputs:
21C
22C  PTCONV-----ND-: Point convectif = TRUE
23C  RATQSC-----ND-: Largeur normalisee de la distribution
24C  CLDF-----ND-: Fraction nuageuse
25C
26C--------------------------------------------------------------------------------
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
35c -- parameters controlling the iteration:
36c --    nmax    : maximum nb of iterations (hopefully never reached)
37c --    epsilon : accuracy of the numerical resolution
38c --    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
57c lconv = true si le calcul a converge (entre autre si qsub < min_q)
58       LOGICAL lconv(klon)
59
60cym
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
81C
82C ***          There is no subgrid-scale condensation;        ***
83C ***   the scheme becomes equivalent to an "all-or-nothing"  ***
84C ***             large-scale condensation scheme.            ***
85C
86
87C
88C ***     Some condensation is produced at the subgrid-scale       ***
89C ***                                                              ***
90C ***       PDF = generalized log-normal distribution (GNO)        ***
91C ***   (k<0 because a lower bound is considered for the PDF)      ***
92C ***                                                              ***
93C ***  -> Determine x (the parameter k of the GNO PDF) such        ***
94C ***  that the contribution of subgrid-scale processes to         ***
95C ***  the in-cloud water content is equal to QSUB(K)              ***
96C ***  (equations (13), (14), (15) + Appendix B of the paper)      ***
97C ***                                                              ***
98C ***    Here, an iterative method is used for this purpose        ***
99C ***    (other numerical methods might be more efficient)         ***
100C ***                                                              ***
101C ***          NB: the "error function" is called ERF              ***
102C ***                 (ERF in double precision)                   ***
103C
104
105c  On commence par eliminer les cas pour lesquels on n'a pas
106c  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
115c   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
124c --  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
146c----------------------------------------------------------------------
147c   Debut des nmax iterations pour trouver la solution.
148c----------------------------------------------------------------------
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
163c -- use asymptotic expression of erf for u and v large:
164c ( -> 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
169c                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
179c -- erfv -> 1.0, use an asymptotic expression of erfv for v large:
180
181             erfcu(i) = 1.0-ERF(u(i))
182c  !!! 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
195c -- 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
204c              print*,'ATTENTION !!! xx(',i,') =', xx(i)
205c           print*,'ATTENTION !!! klon,ND,R,RS,QSUB,PTCONV,RATQSC,CLDF',
206c     .klon,ND,R(i,k),RS(i,k),QSUB(i,k),PTCONV(i,k),RATQSC(i,k),
207c     .CLDF(i,k)
208c              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
219c -- test numerical convergence:
220
221c         print*,'avant test ',i,k,lconv(i),u(i),v(i)
222          if ( ABS(dist(i)/beta(i)) .LT. epsilon ) then
223c           print*,'v-u **2',(v(i)-u(i))**2
224c           print*,'exp v-u **2',exp((v(i)-u(i))**2)
225            ptconv(i,K) = .TRUE.
226            lconv(i)=.true.
227c  borne pour l'exponentielle
228            ratqsc(i,k)=min(2.*(v(i)-u(i))**2,20.)
229            ratqsc(i,k)=sqrt(exp(ratqsc(i,k))-1.)
230            CLDF(i,K) = 0.5 * block(i)
231          else
232            xx(i) = xx(i) - dist(i)/fprime(i)
233          endif
234c         print*,'apres test ',i,k,lconv(i)
235
236        endif ! lconv
237                                    enddo       ! vector
238
239c----------------------------------------------------------------------
240c   Fin des nmax iterations pour trouver la solution.
241        ENDDO ! n
242c----------------------------------------------------------------------
243
244500   CONTINUE  ! K
245
246       RETURN
247       END
248
249
250
Note: See TracBrowser for help on using the repository browser.