source: LMDZ.3.3/branches/rel-LF/libf/phylmd/clouds_gno.F

Last change on this file was 546, checked in by lmdzadmin, 20 years ago

Suite des adaptations au portage sur SGI et VPP pour Prism (Caubel & Demory)
LF

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