source: LMDZ6/branches/Amaury_dev/libf/phylmd/clouds_gno.F90 @ 5441

Last change on this file since 5441 was 5117, checked in by abarral, 5 months ago

rename modules properly lmdz_*
move some unused files to obsolete/
(lint) uppercase fortran keywords

  • 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
RevLine 
[1992]1
[1279]2! $Header$
[524]3
4
[1992]5! ================================================================================
[524]6
[1992]7SUBROUTINE clouds_gno(klon, nd, r, rs, qsub, ptconv, ratqsc, cldf)
8  IMPLICIT NONE
[524]9
[1992]10  ! --------------------------------------------------------------------------------
[524]11
[1992]12  ! Inputs:
[524]13
[1992]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:
[524]21
[1992]22  ! PTCONV-----ND-: Point convectif = TRUE
23  ! RATQSC-----ND-: Largeur normalisee de la distribution
24  ! CLDF-----ND-: Fraction nuageuse
[524]25
[1992]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)
[1279]79      delta(i) = log(mu(i)/qsat)
[1992]80      ! enddo ! vector
[524]81
82
[1992]83      ! ***          There is no subgrid-scale condensation;        ***
84      ! ***   the scheme becomes equivalent to an "all-or-nothing"  ***
85      ! ***             large-scale condensation scheme.            ***
[524]86
87
88
[1992]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)                   ***
[524]104
105
[1992]106      ! On commence par eliminer les cas pour lesquels on n'a pas
107      ! suffisamment d'eau nuageuse.
[524]108
[1992]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.
[524]121        vmax(i) = vmax0
122
[1992]123        beta(i) = qsub(i, k)/mu(i) + exp(-min(0.0,delta(i)))
[524]124
[1992]125        ! --  roots of equation v > vmax:
[524]126
[1279]127        det = delta(i) + vmax(i)*vmax(i)
[1992]128        IF (det<=0.0) vmax(i) = vmax0 + 1.0
[1279]129        det = delta(i) + vmax(i)*vmax(i)
[524]130
[1992]131        IF (det<=0.) THEN
[524]132          xx(i) = -0.0001
[1992]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
[524]142
[1992]143      END IF
[524]144
[1992]145    END DO ! vector
[524]146
[1992]147    ! ----------------------------------------------------------------------
148    ! Debut des nmax iterations pour trouver la solution.
149    ! ----------------------------------------------------------------------
[524]150
[1992]151    DO n = 1, nmax
[524]152
[1992]153      DO i = 1, klon ! vector
154        IF (.NOT. lconv(i)) THEN
[524]155
[1279]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
[524]159
[1992]160          IF (v>vmax(i)) THEN
[524]161
[1992]162            IF (abs(u)>vmax(i) .AND. delta(i)<0.) THEN
[524]163
[1992]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
[5103]169                ! PRINT*,'AUX(',i,',',k,')<0',aux(i),delta(i),beta(i)
[1992]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
[524]176
177            ELSE
178
[1992]179              ! -- erfv -> 1.0, use an asymptotic expression of erfv for v
180              ! large:
[524]181
[1992]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
[524]189
[1992]190            END IF ! ABS(u)
191
[524]192          ELSE
193
[1992]194            ! -- general case:
[524]195
[1992]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
[5103]203              ! PRINT*,'ATTENTION !!! xx(',i,') =', xx(i)
204              ! PRINT*,'ATTENTION !!! klon,ND,R,RS,QSUB,PTCONV,RATQSC,CLDF',
[1992]205              ! .klon,ND,R(i,k),RS(i,k),QSUB(i,k),PTCONV(i,k),RATQSC(i,k),
206              ! .CLDF(i,k)
[5103]207              ! PRINT*,'ATTENTION !!! zu2 zv2 =',zu2(i),zv2(i)
[1992]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
[524]216
[1992]217          ! -- test numerical convergence:
[524]218
[5116]219          ! if (beta(i).lt.1.e-10) THEN
[5103]220          ! PRINT*,'avant test ',i,k,lconv(i),u(i),v(i),beta(i)
[1992]221          ! stop
[5117]222          ! END IF
[1992]223          IF (abs(fprime)<1.E-11) THEN
[5103]224            ! PRINT*,'avant test fprime<.e-11 '
[1992]225            ! s        ,i,k,lconv(i),u(i),v(i),beta(i),fprime(i)
[5103]226            ! PRINT*,'klon,ND,R,RS,QSUB',
[1992]227            ! s        klon,ND,R(i,k),rs(i,k),qsub(i,k)
228            fprime = sign(1.E-11, fprime)
229          END IF
[878]230
231
[1992]232          IF (abs(dist/beta(i))<epsilon) THEN
[5103]233            ! PRINT*,'v-u **2',(v(i)-u(i))**2
234            ! PRINT*,'exp v-u **2',exp((v(i)-u(i))**2)
[1992]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
[1279]242            xx(i) = xx(i) - dist/fprime
[1992]243          END IF
[5103]244          ! PRINT*,'apres test ',i,k,lconv(i)
[524]245
[1992]246        END IF ! lconv
247      END DO ! vector
[524]248
[1992]249      ! ----------------------------------------------------------------------
250      ! Fin des nmax iterations pour trouver la solution.
251    END DO ! n
252    ! ----------------------------------------------------------------------
[524]253
254
[1992]255  END DO
256  ! K
[5105]257
[1992]258END SUBROUTINE clouds_gno
[524]259
260
[1992]261
Note: See TracBrowser for help on using the repository browser.