source: LMDZ6/trunk/libf/phylmd/clouds_gno.F90 @ 3141

Last change on this file since 3141 was 1992, checked in by lguez, 11 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.1 KB
Line 
1
2! $Header$
3
4
5! ================================================================================
6
7SUBROUTINE 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  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, delta(klon), beta(klon)
51  REAL zu2, zv2
52  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  LOGICAL lconv(klon)
60
61  ! cdir arraycomb
62  cldf(1:klon, 1:nd) = 0.0 ! cym
63  ratqsc(1:klon, 1:nd) = 0.0
64  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  DO 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 = rs(i, k)
78      qsat = max(qsat, min_mu)
79      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      IF (qsub(i,k)<min_q) THEN
112        ptconv(i, k) = .FALSE.
113        ratqsc(i, k) = 0.
114        lconv(i) = .TRUE.
115
116        ! 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        ! --  roots of equation v > vmax:
126
127        det = delta(i) + vmax(i)*vmax(i)
128        IF (det<=0.0) vmax(i) = vmax0 + 1.0
129        det = delta(i) + vmax(i)*vmax(i)
130
131        IF (det<=0.) THEN
132          xx(i) = -0.0001
133        ELSE
134          zx1 = -sqrt2*vmax(i)
135          zx2 = sqrt(1.0+delta(i)/(vmax(i)*vmax(i)))
136          xx1 = zx1*(1.0-zx2)
137          xx2 = zx1*(1.0+zx2)
138          xx(i) = 1.01*xx1
139          IF (xx1>=0.0) xx(i) = 0.5*xx2
140        END IF
141        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    DO n = 1, nmax
152
153      DO i = 1, klon ! vector
154        IF (.NOT. lconv(i)) THEN
155
156          u = delta(i)/(xx(i)*sqrt2) + xx(i)/(2.*sqrt2)
157          v = delta(i)/(xx(i)*sqrt2) - xx(i)/(2.*sqrt2)
158          v2 = v*v
159
160          IF (v>vmax(i)) THEN
161
162            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              exdel = beta(i)*exp(delta(i))
167              aux(i) = 2.0*delta(i)*(1.-exdel)/(1.+exdel)
168              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              xx(i) = -sqrt(aux(i))
173              block = exp(-v*v)/v/sqrtpi
174              dist = 0.0
175              fprime = 1.0
176
177            ELSE
178
179              ! -- erfv -> 1.0, use an asymptotic expression of erfv for v
180              ! large:
181
182              erfcu = 1.0 - erf(u)
183              ! !!! 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)*(exp(-delta(i))-u*aux(i)/coeff)/coeff
189
190            END IF ! ABS(u)
191
192          ELSE
193
194            ! -- general case:
195
196            erfcu = 1.0 - erf(u)
197            erfcv = 1.0 - erf(v)
198            block = erfcv
199            dist = erfcu/erfcv - beta(i)
200            zu2 = u*u
201            zv2 = v2
202            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                (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          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          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            ptconv(i, k) = .TRUE.
236            lconv(i) = .TRUE.
237            ! borne pour l'exponentielle
238            ratqsc(i, k) = min(2.*(v-u)*(v-u), 20.)
239            ratqsc(i, k) = sqrt(exp(ratqsc(i,k))-1.)
240            cldf(i, k) = 0.5*block
241          ELSE
242            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  RETURN
258END SUBROUTINE clouds_gno
259
260
261
Note: See TracBrowser for help on using the repository browser.