GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: phylmd/clouds_gno.F90 Lines: 68 71 95.8 %
Date: 2023-06-30 12:56:34 Branches: 41 48 85.4 %

Line Branch Exec Source
1
2
! $Header$
3
4
5
! ================================================================================
6
7
144
SUBROUTINE 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
288
  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
288
  REAL mu(klon), qsat, delta(klon), beta(klon)
51
  REAL zu2, zv2
52
288
  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
144
  LOGICAL lconv(klon)
60
61
  ! cdir arraycomb
62

5588064
  cldf(1:klon, 1:nd) = 0.0 ! cym
63

5588064
  ratqsc(1:klon, 1:nd) = 0.0
64

5588064
  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
5760
  DO k = 1, nd
73
74
5587920
    DO i = 1, klon ! vector
75
5582304
      mu(i) = r(i, k)
76
5582304
      mu(i) = max(mu(i), min_mu)
77
5582304
      qsat = rs(i, k)
78
5582304
      qsat = max(qsat, min_mu)
79
5582304
      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
5587920
      IF (qsub(i,k)<min_q) THEN
112
5333407
        ptconv(i, k) = .FALSE.
113
5333407
        ratqsc(i, k) = 0.
114
5333407
        lconv(i) = .TRUE.
115
116
        ! Rien on a deja initialise
117
118
      ELSE
119
120
248897
        lconv(i) = .FALSE.
121
248897
        vmax(i) = vmax0
122
123
248897
        beta(i) = qsub(i, k)/mu(i) + exp(-min(0.0,delta(i)))
124
125
        ! --  roots of equation v > vmax:
126
127
248897
        det = delta(i) + vmax(i)*vmax(i)
128
248897
        IF (det<=0.0) vmax(i) = vmax0 + 1.0
129
248897
        det = delta(i) + vmax(i)*vmax(i)
130
131
248897
        IF (det<=0.) THEN
132
          xx(i) = -0.0001
133
        ELSE
134
248897
          zx1 = -sqrt2*vmax(i)
135
248897
          zx2 = sqrt(1.0+delta(i)/(vmax(i)*vmax(i)))
136
248897
          xx1 = zx1*(1.0-zx2)
137
248897
          xx2 = zx1*(1.0+zx2)
138
248897
          xx(i) = 1.01*xx1
139
248897
          IF (xx1>=0.0) xx(i) = 0.5*xx2
140
        END IF
141
248897
        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
61920
    DO n = 1, nmax
152
153
55884816
      DO i = 1, klon ! vector
154
55879200
        IF (.NOT. lconv(i)) THEN
155
156
656763
          u = delta(i)/(xx(i)*sqrt2) + xx(i)/(2.*sqrt2)
157
656763
          v = delta(i)/(xx(i)*sqrt2) - xx(i)/(2.*sqrt2)
158
656763
          v2 = v*v
159
160
656763
          IF (v>vmax(i)) THEN
161
162

118954
            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
105721
              exdel = beta(i)*exp(delta(i))
167
105721
              aux(i) = 2.0*delta(i)*(1.-exdel)/(1.+exdel)
168
105721
              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
105721
              xx(i) = -sqrt(aux(i))
173
105721
              block = exp(-v*v)/v/sqrtpi
174
              dist = 0.0
175
105721
              fprime = 1.0
176
177
            ELSE
178
179
              ! -- erfv -> 1.0, use an asymptotic expression of erfv for v
180
              ! large:
181
182
13233
              erfcu = 1.0 - erf(u)
183
              ! !!! ATTENTION : rajout d'un seuil pour l'exponentiel
184
13233
              aux(i) = sqrtpi*erfcu*exp(min(v2,100.))
185
13233
              coeff = 1.0 - 0.5/(v2) + 0.75/(v2*v2)
186
13233
              block = coeff*exp(-v2)/v/sqrtpi
187
13233
              dist = v*aux(i)/coeff - beta(i)
188
13233
              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
537809
            erfcu = 1.0 - erf(u)
197
537809
            erfcv = 1.0 - erf(v)
198
            block = erfcv
199
537809
            dist = erfcu/erfcv - beta(i)
200
537809
            zu2 = u*u
201
            zv2 = v2
202

537809
            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
537809
                (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
656763
          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
656763
          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
248897
            ptconv(i, k) = .TRUE.
236
248897
            lconv(i) = .TRUE.
237
            ! borne pour l'exponentielle
238
248897
            ratqsc(i, k) = min(2.*(v-u)*(v-u), 20.)
239
248897
            ratqsc(i, k) = sqrt(exp(ratqsc(i,k))-1.)
240
248897
            cldf(i, k) = 0.5*block
241
          ELSE
242
407866
            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
144
  RETURN
258
END SUBROUTINE clouds_gno
259
260
261