source: LMDZ6/branches/Amaury_dev/libf/phylmd/conlmd.F90 @ 5202

Last change on this file since 5202 was 5153, checked in by abarral, 3 months ago

Revert FCTTRE to INCLUDE to assess impact of inlining
  • 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: 66.2 KB
Line 
1! $Header$
2
3SUBROUTINE conlmd(dtime, paprs, pplay, t, q, conv_q, d_t, d_q, rain, snow, &
4        ibas, itop)
5  USE dimphy
6  USE lmdz_yoethf
7  USE lmdz_yomcst
8
9  IMPLICIT NONE
10  ! ======================================================================
11  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
12  ! Objet: Schema de convection utilis'e dans le modele du LMD
13  ! Ajustement humide (Manabe) + Ajustement convectif (Kuo)
14  ! ======================================================================
15
16  ! Arguments:
17
18  REAL dtime ! pas d'integration (s)
19  REAL paprs(klon, klev + 1) ! pression inter-couche (Pa)
20  REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
21  REAL t(klon, klev) ! temperature (K)
22  REAL q(klon, klev) ! humidite specifique (kg/kg)
23  REAL conv_q(klon, klev) ! taux de convergence humidite (g/g/s)
24
25  REAL d_t(klon, klev) ! incrementation temperature
26  REAL d_q(klon, klev) ! incrementation humidite
27  REAL rain(klon) ! pluies (mm/s)
28  REAL snow(klon) ! neige (mm/s)
29  INTEGER ibas(klon) ! niveau du bas
30  INTEGER itop(klon) ! niveau du haut
31
32  LOGICAL usekuo ! utiliser convection profonde (schema Kuo)
33  PARAMETER (usekuo = .TRUE.)
34
35  REAL d_t_bis(klon, klev)
36  REAL d_q_bis(klon, klev)
37  REAL rain_bis(klon)
38  REAL snow_bis(klon)
39  INTEGER ibas_bis(klon)
40  INTEGER itop_bis(klon)
41  REAL d_ql(klon, klev), d_ql_bis(klon, klev)
42  REAL rneb(klon, klev), rneb_bis(klon, klev)
43
44  INTEGER i, k
45  REAL zlvdcp, zlsdcp, zdelta, zz, za, zb
46
47  ! cc      CALL fiajh ! ancienne version de Convection Manabe
48  CALL conman &                    ! nouvelle version de Convection
49          ! Manabe
50          (dtime, paprs, pplay, t, q, d_t, d_q, d_ql, rneb, rain, snow, ibas, itop)
51
52  IF (usekuo) THEN
53    ! cc      CALL fiajc ! ancienne version de Convection Kuo
54    CALL conkuo &                  ! nouvelle version de Convection
55            ! Kuo
56            (dtime, paprs, pplay, t, q, conv_q, d_t_bis, d_q_bis, d_ql_bis, &
57            rneb_bis, rain_bis, snow_bis, ibas_bis, itop_bis)
58    DO k = 1, klev
59      DO i = 1, klon
60        d_t(i, k) = d_t(i, k) + d_t_bis(i, k)
61        d_q(i, k) = d_q(i, k) + d_q_bis(i, k)
62        d_ql(i, k) = d_ql(i, k) + d_ql_bis(i, k)
63      END DO
64    END DO
65    DO i = 1, klon
66      rain(i) = rain(i) + rain_bis(i)
67      snow(i) = snow(i) + snow_bis(i)
68      ibas(i) = min(ibas(i), ibas_bis(i))
69      itop(i) = max(itop(i), itop_bis(i))
70    END DO
71  END IF
72
73  ! L'eau liquide convective est dispersee dans l'air:
74
75  DO k = 1, klev
76    DO i = 1, klon
77      zlvdcp = rlvtt / rcpd / (1.0 + rvtmp2 * q(i, k))
78      zlsdcp = rlstt / rcpd / (1.0 + rvtmp2 * q(i, k))
79      zdelta = max(0., sign(1., rtt - t(i, k)))
80      zz = d_ql(i, k) ! re-evap. de l'eau liquide
81      zb = max(0.0, zz)
82      za = -max(0.0, zz) * (zlvdcp * (1. - zdelta) + zlsdcp * zdelta)
83      d_t(i, k) = d_t(i, k) + za
84      d_q(i, k) = d_q(i, k) + zb
85    END DO
86  END DO
87
88END SUBROUTINE conlmd
89SUBROUTINE conman(dtime, paprs, pplay, t, q, d_t, d_q, d_ql, rneb, rain, &
90        snow, ibas, itop)
91  USE dimphy
92  USE lmdz_yoethf
93
94  USE lmdz_yomcst
95
96  IMPLICIT NONE
97 INCLUDE "FCTTRE.h"
98  ! ======================================================================
99  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19970324
100  ! Objet: ajustement humide convectif avec la possibilite de faire
101  ! l'ajustement sur une fraction de la maille.
102  ! Methode: On impose une distribution uniforme pour la vapeur d'eau
103  ! au sein d'une maille. On applique la procedure d'ajustement
104  ! successivement a la totalite, 75%, 50%, 25% et 5% de la maille
105  ! jusqu'a ce que l'ajustement a lieu. J'espere que ceci augmente
106  ! les activites convectives et corrige le biais "trop froid et sec"
107  ! du modele.
108  ! ======================================================================
109
110  REAL dtime ! pas d'integration (s)
111  REAL t(klon, klev) ! temperature (K)
112  REAL q(klon, klev) ! humidite specifique (kg/kg)
113  REAL paprs(klon, klev + 1) ! pression inter-couche (Pa)
114  REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
115
116  REAL d_t(klon, klev) ! incrementation temperature
117  REAL d_q(klon, klev) ! incrementation humidite
118  REAL d_ql(klon, klev) ! incrementation eau liquide
119  REAL rneb(klon, klev) ! nebulosite
120  REAL rain(klon) ! pluies (mm/s)
121  REAL snow(klon) ! neige (mm/s)
122  INTEGER ibas(klon) ! niveau du bas
123  INTEGER itop(klon) ! niveau du haut
124
125  LOGICAL afaire(klon) ! .TRUE. implique l'ajustement
126  LOGICAL accompli(klon) ! .TRUE. si l'ajustement est effectif
127
128  INTEGER nb ! nombre de sous-fractions a considere
129  PARAMETER (nb = 1)
130  ! cc      PARAMETER (nb=3)
131
132  REAL ratqs ! largeur de la distribution pour vapeur d'eau
133  PARAMETER (ratqs = 0.05)
134
135  REAL w_q(klon, klev)
136  REAL w_d_t(klon, klev), w_d_q(klon, klev), w_d_ql(klon, klev)
137  REAL w_rneb(klon, klev)
138  REAL w_rain(klon), w_snow(klon)
139  INTEGER w_ibas(klon), w_itop(klon)
140  REAL zq1, zq2
141  INTEGER i, k, n
142
143  REAL t_coup
144  PARAMETER (t_coup = 234.0)
145  REAL zdp1, zdp2
146  REAL zqs1, zqs2, zdqs1, zdqs2
147  REAL zgamdz
148  REAL zflo ! flotabilite
149  REAL zsat ! sur-saturation
150  REAL zdelta, zcor, zcvm5
151  LOGICAL imprim
152
153  INTEGER ncpt
154  SAVE ncpt
155  !$OMP THREADPRIVATE(ncpt)
156  REAL frac(nb) ! valeur de la maille fractionnelle
157  SAVE frac
158  !$OMP THREADPRIVATE(frac)
159  INTEGER opt_cld(nb) ! option pour le modele nuageux
160  SAVE opt_cld
161  !$OMP THREADPRIVATE(opt_cld)
162  LOGICAL appel1er
163  SAVE appel1er
164  !$OMP THREADPRIVATE(appel1er)
165
166  ! Fonctions thermodynamiques:
167
168  DATA frac/1.0/
169  DATA opt_cld/4/
170  ! cc      DATA frac    / 1.0, 0.50, 0.25/
171  ! cc      DATA opt_cld / 4,   4,    4/
172
173  DATA appel1er/.TRUE./
174  DATA ncpt/0/
175
176  IF (appel1er) THEN
177    PRINT *, 'conman, nb:', nb
178    PRINT *, 'conman, frac:', frac
179    PRINT *, 'conman, opt_cld:', opt_cld
180    appel1er = .FALSE.
181  END IF
182
183  ! Initialiser les sorties a zero:
184
185  DO k = 1, klev
186    DO i = 1, klon
187      d_t(i, k) = 0.0
188      d_q(i, k) = 0.0
189      d_ql(i, k) = 0.0
190      rneb(i, k) = 0.0
191    END DO
192  END DO
193  DO i = 1, klon
194    ibas(i) = klev
195    itop(i) = 1
196    rain(i) = 0.0
197    snow(i) = 0.0
198  END DO
199
200  ! S'il n'y a pas d'instabilite conditionnelle,
201  ! pas la penne de se fatiguer:
202
203  DO i = 1, klon
204    afaire(i) = .FALSE.
205  END DO
206  DO k = 1, klev - 1
207    DO i = 1, klon
208      IF (thermcep) THEN
209        zdelta = max(0., sign(1., rtt - t(i, k)))
210        zcvm5 = r5les * rlvtt * (1. - zdelta) + zdelta * r5ies * rlstt
211        zcvm5 = zcvm5 / rcpd / (1.0 + rvtmp2 * q(i, k))
212        zqs1 = r2es * foeew(t(i, k), zdelta) / pplay(i, k)
213        zqs1 = min(0.5, zqs1)
214        zcor = 1. / (1. - retv * zqs1)
215        zqs1 = zqs1 * zcor
216        zdqs1 = foede(t(i, k), zdelta, zcvm5, zqs1, zcor)
217
218        zdelta = max(0., sign(1., rtt - t(i, k + 1)))
219        zcvm5 = r5les * rlvtt * (1. - zdelta) + zdelta * r5ies * rlstt
220        zcvm5 = zcvm5 / rcpd / (1.0 + rvtmp2 * q(i, k + 1))
221        zqs2 = r2es * foeew(t(i, k + 1), zdelta) / pplay(i, k + 1)
222        zqs2 = min(0.5, zqs2)
223        zcor = 1. / (1. - retv * zqs2)
224        zqs2 = zqs2 * zcor
225        zdqs2 = foede(t(i, k + 1), zdelta, zcvm5, zqs2, zcor)
226      ELSE
227        IF (t(i, k)<t_coup) THEN
228          zqs1 = qsats(t(i, k)) / pplay(i, k)
229          zdqs1 = dqsats(t(i, k), zqs1)
230
231          zqs2 = qsats(t(i, k + 1)) / pplay(i, k + 1)
232          zdqs2 = dqsats(t(i, k + 1), zqs2)
233        ELSE
234          zqs1 = qsatl(t(i, k)) / pplay(i, k)
235          zdqs1 = dqsatl(t(i, k), zqs1)
236
237          zqs2 = qsatl(t(i, k + 1)) / pplay(i, k + 1)
238          zdqs2 = dqsatl(t(i, k + 1), zqs2)
239        END IF
240      END IF
241      zdp1 = paprs(i, k) - paprs(i, k + 1)
242      zdp2 = paprs(i, k + 1) - paprs(i, k + 2)
243      zgamdz = -(pplay(i, k) - pplay(i, k + 1)) / paprs(i, k + 1) / rcpd * (rd * (t(i, &
244              k) * zdp1 + t(i, k + 1) * zdp2) / (zdp1 + zdp2) + rlvtt * (zqs1 * zdp1 + zqs2 * zdp2) / (zdp1 + &
245              zdp2)) / (1.0 + (zdqs1 * zdp1 + zdqs2 * zdp2) / (zdp1 + zdp2))
246      zflo = t(i, k) + zgamdz - t(i, k + 1)
247      zsat = (q(i, k) - zqs1) * zdp1 + (q(i, k + 1) - zqs2) * zdp2
248      IF (zflo>0.0) afaire(i) = .TRUE.
249      ! erreur         IF (zflo.GT.0.0 .AND. zsat.GT.0.0) afaire(i) = .TRUE.
250    END DO
251  END DO
252
253  imprim = mod(ncpt, 48) == 0
254  DO n = 1, nb
255
256    DO k = 1, klev
257      DO i = 1, klon
258        IF (afaire(i)) THEN
259          zq1 = q(i, k) * (1.0 - ratqs)
260          zq2 = q(i, k) * (1.0 + ratqs)
261          w_q(i, k) = zq2 - frac(n) / 2.0 * (zq2 - zq1)
262        END IF
263      END DO
264    END DO
265
266    CALL conmanv(dtime, paprs, pplay, t, w_q, afaire, opt_cld(n), w_d_t, &
267            w_d_q, w_d_ql, w_rneb, w_rain, w_snow, w_ibas, w_itop, accompli, &
268            imprim)
269    DO k = 1, klev
270      DO i = 1, klon
271        IF (afaire(i) .AND. accompli(i)) THEN
272          d_t(i, k) = w_d_t(i, k) * frac(n)
273          d_q(i, k) = w_d_q(i, k) * frac(n)
274          d_ql(i, k) = w_d_ql(i, k) * frac(n)
275          IF (nint(w_rneb(i, k))==1) rneb(i, k) = frac(n)
276        END IF
277      END DO
278    END DO
279    DO i = 1, klon
280      IF (afaire(i) .AND. accompli(i)) THEN
281        rain(i) = w_rain(i) * frac(n)
282        snow(i) = w_snow(i) * frac(n)
283        ibas(i) = min(ibas(i), w_ibas(i))
284        itop(i) = max(itop(i), w_itop(i))
285      END IF
286    END DO
287    DO i = 1, klon
288      IF (afaire(i) .AND. accompli(i)) afaire(i) = .FALSE.
289    END DO
290
291  END DO
292
293  ncpt = ncpt + 1
294
295END SUBROUTINE conman
296SUBROUTINE conmanv(dtime, paprs, pplay, t, q, afaire, opt_cld, d_t, d_q, &
297        d_ql, rneb, rain, snow, ibas, itop, accompli, imprim)
298  USE dimphy
299  USE lmdz_yoethf
300
301  USE lmdz_yomcst
302
303  IMPLICIT NONE
304 INCLUDE "FCTTRE.h"
305  ! ======================================================================
306  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
307  ! Objet: ajustement humide (convection proposee par Manabe).
308  ! Pour une colonne verticale, il peut avoir plusieurs blocs
309  ! necessitant l'ajustement. ibas est le bas du plus bas bloc
310  ! et itop est le haut du plus haut bloc
311  ! ======================================================================
312
313  ! Arguments:
314
315  REAL dtime ! pas d'integration (s)
316  REAL t(klon, klev) ! temperature (K)
317  REAL q(klon, klev) ! humidite specifique (kg/kg)
318  REAL paprs(klon, klev + 1) ! pression inter-couche (Pa)
319  REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
320  INTEGER opt_cld ! comment traiter l'eau liquide
321  LOGICAL afaire(klon) ! .TRUE. si le point est a faire (Input)
322  LOGICAL imprim ! .T. pour imprimer quelques diagnostiques
323
324  REAL d_t(klon, klev) ! incrementation temperature
325  REAL d_q(klon, klev) ! incrementation humidite
326  REAL d_ql(klon, klev) ! incrementation eau liquide
327  REAL rneb(klon, klev) ! nebulosite
328  REAL rain(klon) ! pluies (mm/s)
329  REAL snow(klon) ! neige (mm/s)
330  INTEGER ibas(klon) ! niveau du bas
331  INTEGER itop(klon) ! niveau du haut
332  LOGICAL accompli(klon) ! .TRUE. si l'ajustement a eu lieu (Output)
333
334  ! Quelques options:
335
336  LOGICAL new_top ! re-calculer sommet quand re-ajustement est fait
337  PARAMETER (new_top = .FALSE.)
338  LOGICAL evap_prec ! evaporation de pluie au-dessous de convection
339  PARAMETER (evap_prec = .TRUE.)
340  REAL coef_eva
341  PARAMETER (coef_eva = 1.0E-05)
342  REAL t_coup
343  PARAMETER (t_coup = 234.0)
344  REAL seuil_vap
345  PARAMETER (seuil_vap = 1.0E-10)
346  LOGICAL old_tau ! implique precip nulle, si vrai.
347  PARAMETER (old_tau = .FALSE.)
348  REAL toliq(klon) ! rapport entre l'eau nuageuse et l'eau precipitante
349  REAL dpmin, tomax !Epaisseur faible, rapport eau liquide plus grande
350  PARAMETER (dpmin = 0.15, tomax = 0.97)
351  REAL dpmax, tomin !Epaisseur grande, rapport eau liquide plus faible
352  PARAMETER (dpmax = 0.30, tomin = 0.05)
353  REAL deep_sig, deep_to ! au dela de deep_sig, utiliser deep_to
354  PARAMETER (deep_sig = 0.50, deep_to = 0.05)
355  LOGICAL exigent ! implique un calcul supplementaire pour Qs
356  PARAMETER (exigent = .FALSE.)
357
358  INTEGER kbase
359  PARAMETER (kbase = 0)
360
361  ! Variables locales:
362
363  INTEGER nexpo
364  INTEGER i, k, k1min, k1max, k2min, k2max, is
365  REAL zgamdz(klon, klev - 1)
366  REAL zt(klon, klev), zq(klon, klev)
367  REAL zqs(klon, klev), zdqs(klon, klev)
368  REAL zqmqsdp(klon, klev)
369  REAL ztnew(klon, klev), zqnew(klon, klev)
370  REAL zcond(klon), zvapo(klon), zrapp(klon)
371  REAL zrfl(klon), zrfln, zqev, zqevt
372  REAL zsat(klon) ! sur-saturation
373  REAL zflo(klon) ! flotabilite
374  REAL za(klon), zb(klon), zc(klon)
375  INTEGER k1(klon), k2(klon)
376  REAL zdelta, zcor, zcvm5
377  REAL delp(klon, klev)
378  LOGICAL possible(klon), todo(klon), etendre(klon)
379  LOGICAL aller(klon), todobis(klon)
380  REAL zalfa
381  INTEGER nbtodo, nbdone
382
383  ! Fonctions thermodynamiques:
384
385  DO k = 1, klev
386    DO i = 1, klon
387      delp(i, k) = paprs(i, k) - paprs(i, k + 1)
388    END DO
389  END DO
390
391  ! Initialiser les sorties a zero
392
393  DO k = 1, klev
394    DO i = 1, klon
395      d_t(i, k) = 0.0
396      d_q(i, k) = 0.0
397      d_ql(i, k) = 0.0
398      rneb(i, k) = 0.0
399    END DO
400  END DO
401  DO i = 1, klon
402    ibas(i) = klev
403    itop(i) = 1
404    rain(i) = 0.0
405    snow(i) = 0.0
406    accompli(i) = .FALSE.
407  END DO
408
409  ! Preparations
410
411  DO k = 1, klev
412    DO i = 1, klon
413      IF (afaire(i)) THEN
414        zt(i, k) = t(i, k)
415        zq(i, k) = q(i, k)
416
417        ! Calculer Qs et L/Cp*dQs/dT
418
419        IF (thermcep) THEN
420          zdelta = max(0., sign(1., rtt - zt(i, k)))
421          zcvm5 = r5les * rlvtt * (1. - zdelta) + zdelta * r5ies * rlstt
422          zcvm5 = zcvm5 / rcpd / (1.0 + rvtmp2 * zq(i, k))
423          zqs(i, k) = r2es * foeew(zt(i, k), zdelta) / pplay(i, k)
424          zqs(i, k) = min(0.5, zqs(i, k))
425          zcor = 1. / (1. - retv * zqs(i, k))
426          zqs(i, k) = zqs(i, k) * zcor
427          zdqs(i, k) = foede(zt(i, k), zdelta, zcvm5, zqs(i, k), zcor)
428        ELSE
429          IF (zt(i, k)<t_coup) THEN
430            zqs(i, k) = qsats(zt(i, k)) / pplay(i, k)
431            zdqs(i, k) = dqsats(zt(i, k), zqs(i, k))
432          ELSE
433            zqs(i, k) = qsatl(zt(i, k)) / pplay(i, k)
434            zdqs(i, k) = dqsatl(zt(i, k), zqs(i, k))
435          END IF
436        END IF
437
438        ! Calculer (q-qs)*dp
439        zqmqsdp(i, k) = (zq(i, k) - zqs(i, k)) * delp(i, k)
440      END IF
441    END DO
442  END DO
443
444  ! -----zgama is the moist convective lapse rate (-dT/dz).
445  ! -----zgamdz(*,k) est la difference minimale autorisee des temperatures
446  ! -----entre deux couches (k et k+1), c.a.d. si T(k+1)-T(k) est inferieur
447  ! -----a zgamdz(*,k), alors ces 2 couches sont instables conditionnellement
448
449  DO k = 1, klev - 1
450    DO i = 1, klon
451      IF (afaire(i)) THEN
452        zgamdz(i, k) = -(pplay(i, k) - pplay(i, k + 1)) / paprs(i, k + 1) / rcpd * (rd * (zt(&
453                i, k) * delp(i, k) + zt(i, k + 1) * delp(i, k + 1)) / (delp(i, k) + delp(i, &
454                k + 1)) + rlvtt * (zqs(i, k) * delp(i, k) + zqs(i, k + 1) * delp(i, k + 1)) / (delp(i, &
455                k) + delp(i, k + 1))) / (1.0 + (zdqs(i, k) * delp(i, k) + zdqs(i, k + 1) * delp(i, &
456                k + 1)) / (delp(i, k) + delp(i, k + 1)))
457      END IF
458    END DO
459  END DO
460
461  ! On cherche la presence simultanee d'instabilite conditionnelle
462  ! et de sur-saturation. Sinon, pas la penne de se fatiguer:
463
464  DO i = 1, klon
465    possible(i) = .FALSE.
466  END DO
467  DO k = 2, klev
468    DO i = 1, klon
469      IF (afaire(i)) THEN
470        zflo(i) = zt(i, k - 1) + zgamdz(i, k - 1) - zt(i, k)
471        zsat(i) = zqmqsdp(i, k) + zqmqsdp(i, k - 1)
472        IF (zflo(i)>0.0 .AND. zsat(i)>0.0) possible(i) = .TRUE.
473      END IF
474    END DO
475  END DO
476
477  DO i = 1, klon
478    IF (possible(i)) THEN
479      k1(i) = kbase
480      k2(i) = k1(i) + 1
481    END IF
482  END DO
483
484  810 CONTINUE ! chercher le bas de la colonne a ajuster
485
486  k2min = klev
487  DO i = 1, klon
488    todo(i) = .FALSE.
489    aller(i) = .TRUE.
490    IF (possible(i)) k2min = min(k2min, k2(i))
491  END DO
492  IF (k2min==klev) GO TO 860
493  DO k = k2min, klev - 1
494    DO i = 1, klon
495      IF (possible(i) .AND. k>=k2(i) .AND. aller(i)) THEN
496        zflo(i) = zt(i, k) + zgamdz(i, k) - zt(i, k + 1)
497        zsat(i) = zqmqsdp(i, k) + zqmqsdp(i, k + 1)
498        IF (zflo(i)>0.0 .AND. zsat(i)>0.0) THEN
499          k1(i) = k
500          k2(i) = k + 1
501          todo(i) = .TRUE.
502          aller(i) = .FALSE.
503        END IF
504      END IF
505    END DO
506  END DO
507  DO i = 1, klon
508    IF (possible(i) .AND. aller(i)) THEN
509      todo(i) = .FALSE.
510      k1(i) = klev
511      k2(i) = klev
512    END IF
513  END DO
514
515  ! CC      DO i = 1, klon
516  ! CC      IF (possible(i)) THEN
517  ! CC  811    k2(i) = k2(i) + 1
518  ! CC         IF (k2(i) .GT. klev) THEN
519  ! CC            todo(i) = .FALSE.
520  ! CC            GOTO 812
521  ! CC         ENDIF
522  ! CC         k = k2(i)
523  ! CC         zflo(i) = zt(i,k-1) + zgamdz(i,k-1) - zt(i,k)
524  ! CC         zsat(i) = zqmqsdp(i,k) + zqmqsdp(i,k-1)
525  ! CC         IF (zflo(i).LE.0.0 .OR. zsat(i).LE.0.0) GOTO 811
526  ! CC         k1(i) = k2(i) - 1
527  ! CC         todo(i) = .TRUE.
528  ! CC      ENDIF
529  ! CC  812 CONTINUE
530  ! CC      ENDDO
531
532  820 CONTINUE ! chercher le haut de la colonne
533
534  k2min = klev
535  DO i = 1, klon
536    aller(i) = .TRUE.
537    IF (todo(i)) k2min = min(k2min, k2(i))
538  END DO
539  IF (k2min<klev) THEN
540    DO k = k2min, klev
541      DO i = 1, klon
542        IF (todo(i) .AND. k>k2(i) .AND. aller(i)) THEN
543          zsat(i) = zsat(i) + zqmqsdp(i, k)
544          zflo(i) = zt(i, k - 1) + zgamdz(i, k - 1) - zt(i, k)
545          IF (zflo(i)<=0.0 .OR. zsat(i)<=0.0) THEN
546            aller(i) = .FALSE.
547          ELSE
548            k2(i) = k
549          END IF
550        END IF
551      END DO
552    END DO
553    ! error      is = 0
554    ! error      DO i = 1, klon
555    ! error      IF(todo(i).AND.aller(i)) THEN
556    ! error         is = is + 1
557    ! error         todo(i) = .FALSE.
558    ! error         k2(i) = klev
559    ! error      ENDIF
560    ! error      ENDDO
561    ! error      IF (is.GT.0) THEN
562    ! error         PRINT*, "Bizard. je pourrais continuer mais j arrete"
563    ! error         CALL abort
564    ! error      ENDIF
565  END IF
566
567  ! CC      DO i = 1, klon
568  ! CC      IF (todo(i)) THEN
569  ! CC  821    CONTINUE
570  ! CC         IF (k2(i) .EQ. klev) GOTO 822
571  ! CC         k = k2(i) + 1
572  ! CC         zsat(i) = zsat(i) + zqmqsdp(i,k)
573  ! CC         zflo(i) = zt(i,k-1) + zgamdz(i,k-1) - zt(i,k)
574  ! CC         IF (zflo(i).LE.0.0 .OR. zsat(i).LE.0.0) GOTO 822
575  ! CC         k2(i) = k
576  ! CC         GOTO 821
577  ! CC      ENDIF
578  ! CC  822 CONTINUE
579  ! CC      ENDDO
580
581  830 CONTINUE ! faire l'ajustement en sachant k1 et k2
582
583  is = 0
584  DO i = 1, klon
585    IF (todo(i)) THEN
586      IF (k2(i)<=k1(i)) is = is + 1
587    END IF
588  END DO
589  IF (is>0) THEN
590    PRINT *, 'Impossible: k1 trop grand ou k2 trop petit'
591    PRINT *, 'is=', is
592    CALL abort
593  END IF
594
595  k1min = klev
596  k1max = 1
597  k2max = 1
598  DO i = 1, klon
599    IF (todo(i)) THEN
600      k1min = min(k1min, k1(i))
601      k1max = max(k1max, k1(i))
602      k2max = max(k2max, k2(i))
603    END IF
604  END DO
605
606  DO i = 1, klon
607    IF (todo(i)) THEN
608      k = k1(i)
609      za(i) = 0.
610      zb(i) = (rcpd * (1. + zdqs(i, k)) * (zt(i, k) - za(i)) - rlvtt * (zqs(i, k) - zq(i, &
611              k))) * delp(i, k)
612      zc(i) = delp(i, k) * rcpd * (1. + zdqs(i, k))
613    END IF
614  END DO
615
616  DO k = k1min, k2max
617    DO i = 1, klon
618      IF (todo(i) .AND. k>=(k1(i) + 1) .AND. k<=k2(i)) THEN
619        za(i) = za(i) + zgamdz(i, k - 1)
620        zb(i) = zb(i) + (rcpd * (1. + zdqs(i, k)) * (zt(i, k) - za(i)) - rlvtt * (zqs(i, &
621                k) - zq(i, k))) * delp(i, k)
622        zc(i) = zc(i) + delp(i, k) * rcpd * (1. + zdqs(i, k))
623      END IF
624    END DO
625  END DO
626
627  DO i = 1, klon
628    IF (todo(i)) THEN
629      k = k1(i)
630      ztnew(i, k) = zb(i) / zc(i)
631      zqnew(i, k) = zqs(i, k) + (ztnew(i, k) - zt(i, k)) * rcpd / rlvtt * zdqs(i, k)
632    END IF
633  END DO
634
635  DO k = k1min, k2max
636    DO i = 1, klon
637      IF (todo(i) .AND. k>=(k1(i) + 1) .AND. k<=k2(i)) THEN
638        ztnew(i, k) = ztnew(i, k - 1) + zgamdz(i, k - 1)
639        zqnew(i, k) = zqs(i, k) + (ztnew(i, k) - zt(i, k)) * rcpd / rlvtt * zdqs(i, k)
640      END IF
641    END DO
642  END DO
643
644  ! Quantite de condensation produite pendant l'ajustement:
645
646  DO i = 1, klon
647    zcond(i) = 0.0
648  END DO
649  DO k = k1min, k2max
650    DO i = 1, klon
651      IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) THEN
652        rneb(i, k) = 1.0
653        zcond(i) = zcond(i) + (zq(i, k) - zqnew(i, k)) * delp(i, k) / rg
654      END IF
655    END DO
656  END DO
657
658  ! Si condensation negative, effort completement perdu:
659
660  DO i = 1, klon
661    IF (todo(i) .AND. zcond(i)<=0.) todo(i) = .FALSE.
662  END DO
663
664  ! L'ajustement a ete accompli, meme les calculs accessoires
665  ! ne sont pas encore faits:
666
667  DO i = 1, klon
668    IF (todo(i)) accompli(i) = .TRUE.
669  END DO
670
671  ! =====
672  ! Une fois que la condensation a lieu, on doit construire un
673  ! "modele nuageux" pour partager la condensation entre l'eau
674  ! liquide nuageuse et la precipitation (leur rapport toliq
675  ! est calcule selon l'epaisseur nuageuse). Je suppose que
676  ! toliq=tomax quand l'epaisseur nuageuse est inferieure a dpmin,
677  ! et que toliq=tomin quand l'epaisseur depasse dpmax (interpolation
678  ! lineaire entre dpmin et dpmax).
679  ! =====
680  DO i = 1, klon
681    IF (todo(i)) THEN
682      toliq(i) = tomax - ((paprs(i, k1(i)) - paprs(i, k2(i) + 1)) / paprs(i, 1) - dpmin) &
683              * (tomax - tomin) / (dpmax - dpmin)
684      toliq(i) = max(tomin, min(tomax, toliq(i)))
685      IF (pplay(i, k2(i)) / paprs(i, 1)<=deep_sig) toliq(i) = deep_to
686      IF (old_tau) toliq(i) = 1.0
687    END IF
688  END DO
689  ! =====
690  ! On doit aussi determiner la distribution verticale de
691  ! l'eau nuageuse. Plusieurs options sont proposees:
692
693  ! (0) La condensation precipite integralement (toliq ne sera
694  ! pas utilise).
695  ! (1) L'eau liquide est distribuee entre k1 et k2 et proportionnelle
696  ! a la vapeur d'eau locale.
697  ! (2) Elle est distribuee entre k1 et k2 avec une valeur constante.
698  ! (3) Elle est seulement distribuee aux couches ou la vapeur d'eau
699  ! est effectivement diminuee pendant le processus d'ajustement.
700  ! (4) Elle est en fonction (lineaire ou exponentielle) de la
701  ! distance (epaisseur en pression) avec le niveau k1 (la couche
702  ! k1 n'aura donc pas d'eau liquide).
703  ! =====
704
705  IF (opt_cld==0) THEN
706
707    DO i = 1, klon
708      IF (todo(i)) zrfl(i) = zcond(i) / dtime
709    END DO
710
711  ELSE IF (opt_cld==1) THEN
712
713    DO i = 1, klon
714      IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de vapeur d'eau
715    END DO
716    DO k = k1min, k2max
717      DO i = 1, klon
718        IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) zvapo(i) = zvapo(i) + &
719                zqnew(i, k) * delp(i, k) / rg
720      END DO
721    END DO
722    DO i = 1, klon
723      IF (todo(i)) THEN
724        zrapp(i) = toliq(i) * zcond(i) / zvapo(i)
725        zrapp(i) = max(0., min(1., zrapp(i)))
726        zrfl(i) = (1.0 - toliq(i)) * zcond(i) / dtime
727      END IF
728    END DO
729    DO k = k1min, k2max
730      DO i = 1, klon
731        IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) THEN
732          d_ql(i, k) = d_ql(i, k) + zrapp(i) * zqnew(i, k)
733        END IF
734      END DO
735    END DO
736
737  ELSE IF (opt_cld==2) THEN
738
739    DO i = 1, klon
740      IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de masse
741    END DO
742    DO k = k1min, k2max
743      DO i = 1, klon
744        IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) zvapo(i) = zvapo(i) + &
745                delp(i, k) / rg
746      END DO
747    END DO
748    DO k = k1min, k2max
749      DO i = 1, klon
750        IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) THEN
751          d_ql(i, k) = d_ql(i, k) + toliq(i) * zcond(i) / zvapo(i)
752        END IF
753      END DO
754    END DO
755    DO i = 1, klon
756      IF (todo(i)) zrfl(i) = (1.0 - toliq(i)) * zcond(i) / dtime
757    END DO
758
759  ELSE IF (opt_cld==3) THEN
760
761    DO i = 1, klon
762      IF (todo(i)) zvapo(i) = 0.0 ! quantite de l'eau strictement condensee
763    END DO
764    DO k = k1min, k2max
765      DO i = 1, klon
766        IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) zvapo(i) = zvapo(i) + &
767                max(0.0, zq(i, k) - zqnew(i, k)) * delp(i, k) / rg
768      END DO
769    END DO
770    DO k = k1min, k2max
771      DO i = 1, klon
772        IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i) .AND. zvapo(i)>0.0) d_ql(i, &
773                k) = d_ql(i, k) + toliq(i) * zcond(i) / zvapo(i) * max(0.0, zq(i, k) - zqnew &
774                (i, k))
775      END DO
776    END DO
777    DO i = 1, klon
778      IF (todo(i)) zrfl(i) = (1.0 - toliq(i)) * zcond(i) / dtime
779    END DO
780
781  ELSE IF (opt_cld==4) THEN
782
783    nexpo = 3
784    ! cc         nexpo = 1 ! distribution lineaire
785
786    DO i = 1, klon
787      IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de masse
788    END DO ! (avec ponderation)
789    DO k = k1min, k2max
790      DO i = 1, klon
791        IF (todo(i) .AND. k>=(k1(i) + 1) .AND. k<=k2(i)) zvapo(i) = zvapo(i) + &
792                delp(i, k) / rg * (pplay(i, k1(i)) - pplay(i, k))**nexpo
793      END DO
794    END DO
795    DO k = k1min, k2max
796      DO i = 1, klon
797        IF (todo(i) .AND. k>=(k1(i) + 1) .AND. k<=k2(i)) d_ql(i, k) = d_ql(i, &
798                k) + toliq(i) * zcond(i) / zvapo(i) * (pplay(i, k1(i)) - pplay(i, k))**nexpo
799      END DO
800    END DO
801    DO i = 1, klon
802      IF (todo(i)) zrfl(i) = (1.0 - toliq(i)) * zcond(i) / dtime
803    END DO
804
805  ELSE ! valeur non-prevue pour opt_cld
806
807    PRINT *, 'opt_cld est faux:', opt_cld
808    CALL abort
809
810  END IF ! fin de opt_cld
811
812  ! L'eau precipitante peut etre evaporee:
813
814  zalfa = 0.05
815  IF (evap_prec .AND. (k1max>=2)) THEN
816    DO k = k1max - 1, 1, -1
817      DO i = 1, klon
818        IF (todo(i) .AND. k<k1(i) .AND. zrfl(i)>0.0) THEN
819          zqev = max(0.0, (zqs(i, k) - zq(i, k)) * zalfa)
820          zqevt = coef_eva * (1.0 - zq(i, k) / zqs(i, k)) * sqrt(zrfl(i)) * delp(i, k) / &
821                  pplay(i, k) * zt(i, k) * rd / rg
822          zqevt = max(0.0, min(zqevt, zrfl(i))) * rg * dtime / delp(i, k)
823          zqev = min(zqev, zqevt)
824          zrfln = zrfl(i) - zqev * (delp(i, k)) / rg / dtime
825          zq(i, k) = zq(i, k) - (zrfln - zrfl(i)) * (rg / (delp(i, k))) * dtime
826          zt(i, k) = zt(i, k) + (zrfln - zrfl(i)) * (rg / (delp(i, &
827                  k))) * dtime * rlvtt / rcpd / (1.0 + rvtmp2 * zq(i, k))
828          zrfl(i) = zrfln
829        END IF
830      END DO
831    END DO
832  END IF
833
834  ! La temperature de la premiere couche determine la pluie ou la neige:
835
836  DO i = 1, klon
837    IF (todo(i)) THEN
838      IF (zt(i, 1)>rtt) THEN
839        rain(i) = rain(i) + zrfl(i)
840      ELSE
841        snow(i) = snow(i) + zrfl(i)
842      END IF
843    END IF
844  END DO
845
846  ! Mise a jour de la temperature et de l'humidite
847
848  DO k = k1min, k2max
849    DO i = 1, klon
850      IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) THEN
851        zt(i, k) = ztnew(i, k)
852        zq(i, k) = zqnew(i, k)
853      END IF
854    END DO
855  END DO
856
857  ! Re-calculer certaines variables pour etendre et re-ajuster la colonne
858
859  IF (exigent) THEN
860    DO k = 1, klev
861      DO i = 1, klon
862        IF (todo(i)) THEN
863          IF (thermcep) THEN
864            zdelta = max(0., sign(1., rtt - zt(i, k)))
865            zcvm5 = r5les * rlvtt * (1. - zdelta) + zdelta * r5ies * rlstt
866            zcvm5 = zcvm5 / rcpd / (1.0 + rvtmp2 * zq(i, k))
867            zqs(i, k) = r2es * foeew(zt(i, k), zdelta) / pplay(i, k)
868            zqs(i, k) = min(0.5, zqs(i, k))
869            zcor = 1. / (1. - retv * zqs(i, k))
870            zqs(i, k) = zqs(i, k) * zcor
871            zdqs(i, k) = foede(zt(i, k), zdelta, zcvm5, zqs(i, k), zcor)
872          ELSE
873            IF (zt(i, k)<t_coup) THEN
874              zqs(i, k) = qsats(zt(i, k)) / pplay(i, k)
875              zdqs(i, k) = dqsats(zt(i, k), zqs(i, k))
876            ELSE
877              zqs(i, k) = qsatl(zt(i, k)) / pplay(i, k)
878              zdqs(i, k) = dqsatl(zt(i, k), zqs(i, k))
879            END IF
880          END IF
881        END IF
882      END DO
883    END DO
884  END IF
885
886  IF (exigent) THEN
887    DO k = 1, klev - 1
888      DO i = 1, klon
889        IF (todo(i)) THEN
890          zgamdz(i, k) = -(pplay(i, k) - pplay(i, k + 1)) / paprs(i, k + 1) / rcpd * (rd * (&
891                  zt(i, k) * delp(i, k) + zt(i, k + 1) * delp(i, k + 1)) / (delp(i, k) + delp(i, &
892                  k + 1)) + rlvtt * (zqs(i, k) * delp(i, k) + zqs(i, k + 1) * delp(i, k + 1)) / (delp(i, &
893                  k) + delp(i, k + 1))) / (1.0 + (zdqs(i, k) * delp(i, k) + zdqs(i, k + 1) * delp(i, &
894                  k + 1)) / (delp(i, k) + delp(i, k + 1)))
895        END IF
896      END DO
897    END DO
898  END IF
899
900  ! Puisque l'humidite a ete modifiee, on re-fait (q-qs)*dp
901
902  DO k = 1, klev
903    DO i = 1, klon
904      IF (todo(i)) THEN
905        zqmqsdp(i, k) = (zq(i, k) - zqs(i, k)) * delp(i, k)
906      END IF
907    END DO
908  END DO
909
910  ! Verifier si l'on peut etendre le bas de la colonne
911
912  DO i = 1, klon
913    etendre(i) = .FALSE.
914  END DO
915
916  k1max = 1
917  DO i = 1, klon
918    IF (todo(i) .AND. k1(i)>(kbase + 1)) THEN
919      k = k1(i)
920      zflo(i) = zt(i, k - 1) + zgamdz(i, k - 1) - zt(i, k)
921      zsat(i) = zqmqsdp(i, k) + zqmqsdp(i, k - 1)
922      ! sc voici l'ancienne ligne:
923      ! sc         IF (zflo(i).LE.0.0 .OR. zsat(i).LE.0.0) THEN
924      ! sc sylvain: il faut RESPECTER les 2 criteres:
925      IF (zflo(i)>0.0 .AND. zsat(i)>0.0) THEN
926        etendre(i) = .TRUE.
927        k1(i) = k1(i) - 1
928        k1max = max(k1max, k1(i))
929        aller(i) = .TRUE.
930      END IF
931    END IF
932  END DO
933
934  IF (k1max>(kbase + 1)) THEN
935    DO k = k1max, kbase + 1, -1
936      DO i = 1, klon
937        IF (etendre(i) .AND. k<k1(i) .AND. aller(i)) THEN
938          zsat(i) = zsat(i) + zqmqsdp(i, k)
939          zflo(i) = zt(i, k) + zgamdz(i, k) - zt(i, k + 1)
940          IF (zsat(i)<=0.0 .OR. zflo(i)<=0.0) THEN
941            aller(i) = .FALSE.
942          ELSE
943            k1(i) = k
944          END IF
945        END IF
946      END DO
947    END DO
948    DO i = 1, klon
949      IF (etendre(i) .AND. aller(i)) THEN
950        k1(i) = 1
951      END IF
952    END DO
953  END IF
954
955  ! CC      DO i = 1, klon
956  ! CC      IF (etendre(i)) THEN
957  ! CC  840    k = k1(i)
958  ! CC         IF (k.GT.1) THEN
959  ! CC            zsat(i) = zsat(i) + zqmqsdp(i,k-1)
960  ! CC            zflo(i) = zt(i,k-1) + zgamdz(i,k-1) - zt(i,k)
961  ! CC            IF (zflo(i).GT.0.0 .AND. zsat(i).GT.0.0) THEN
962  ! CC               k1(i) = k - 1
963  ! CC               GOTO 840
964  ! CC            ENDIF
965  ! CC         ENDIF
966  ! CC      ENDIF
967  ! CC      ENDDO
968
969  DO i = 1, klon
970    todobis(i) = todo(i)
971    todo(i) = .FALSE.
972  END DO
973  is = 0
974  DO i = 1, klon
975    IF (etendre(i)) THEN
976      todo(i) = .TRUE.
977      is = is + 1
978    END IF
979  END DO
980  IF (is>0) THEN
981    IF (new_top) THEN
982      GO TO 820 ! chercher de nouveau le sommet k2
983    ELSE
984      GO TO 830 ! supposer que le sommet est celui deja trouve
985    END IF
986  END IF
987
988  DO i = 1, klon
989    possible(i) = .FALSE.
990  END DO
991  is = 0
992  DO i = 1, klon
993    IF (todobis(i) .AND. k2(i)<klev) THEN
994      is = is + 1
995      possible(i) = .TRUE.
996    END IF
997  END DO
998  IF (is>0) GO TO 810 !on cherche en haut d'autres blocks
999  ! a ajuster a partir du sommet de la colonne precedente
1000
1001  860 CONTINUE ! Calculer les tendances et diagnostiques
1002  ! cc      PRINT*, "Apres 860"
1003
1004  DO k = 1, klev
1005    DO i = 1, klon
1006      IF (accompli(i)) THEN
1007        d_t(i, k) = zt(i, k) - t(i, k)
1008        zq(i, k) = max(zq(i, k), seuil_vap)
1009        d_q(i, k) = zq(i, k) - q(i, k)
1010      END IF
1011    END DO
1012  END DO
1013
1014  DO i = 1, klon
1015    IF (accompli(i)) THEN
1016      DO k = 1, klev
1017        IF (rneb(i, k)>0.0) THEN
1018          ibas(i) = k
1019          GO TO 807
1020        END IF
1021      END DO
1022      807   CONTINUE
1023      DO k = klev, 1, -1
1024        IF (rneb(i, k)>0.0) THEN
1025          itop(i) = k
1026          GO TO 808
1027        END IF
1028      END DO
1029      808   CONTINUE
1030    END IF
1031  END DO
1032
1033  IF (imprim) THEN
1034    nbtodo = 0
1035    nbdone = 0
1036    DO i = 1, klon
1037      IF (afaire(i)) nbtodo = nbtodo + 1
1038      IF (accompli(i)) nbdone = nbdone + 1
1039    END DO
1040    PRINT *, 'nbTodo, nbDone=', nbtodo, nbdone
1041  END IF
1042
1043END SUBROUTINE conmanv
1044SUBROUTINE conkuo(dtime, paprs, pplay, t, q, conv_q, d_t, d_q, d_ql, rneb, &
1045        rain, snow, ibas, itop)
1046  USE dimphy
1047  USE lmdz_yoethf
1048
1049  USE lmdz_yomcst
1050
1051  IMPLICIT NONE
1052 INCLUDE "FCTTRE.h"
1053  ! ======================================================================
1054  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
1055  ! Objet: Schema de convection de type Kuo (1965).
1056  ! Cette version du code peut calculer le niveau de depart
1057  ! N.B. version vectorielle (le 6 oct. 1997)
1058  ! ======================================================================
1059
1060  ! Arguments:
1061
1062  REAL dtime ! intervalle du temps (s)
1063  REAL paprs(klon, klev + 1) ! pression a inter-couche (Pa)
1064  REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
1065  REAL t(klon, klev) ! temperature (K)
1066  REAL q(klon, klev) ! humidite specifique
1067  REAL conv_q(klon, klev) ! taux de convergence humidite (g/g/s)
1068
1069  REAL d_t(klon, klev) ! incrementation temperature
1070  REAL d_q(klon, klev) ! incrementation humidite
1071  REAL d_ql(klon, klev) ! incrementation eau liquide
1072  REAL rneb(klon, klev) ! nebulosite
1073  REAL rain(klon) ! pluies (mm/s)
1074  REAL snow(klon) ! neige (mm/s)
1075  INTEGER itop(klon) ! niveau du sommet
1076  INTEGER ibas(klon) ! niveau du bas
1077
1078  LOGICAL ldcum(klon) ! convection existe
1079  LOGICAL todo(klon)
1080
1081  ! Quelsques options:
1082
1083  LOGICAL calcfcl ! calculer le niveau de convection libre
1084  PARAMETER (calcfcl = .TRUE.)
1085  INTEGER ldepar ! niveau fixe de convection libre
1086  PARAMETER (ldepar = 4)
1087  INTEGER opt_cld ! comment traiter l'eau liquide
1088  PARAMETER (opt_cld = 4) ! valeur possible: 0, 1, 2, 3 ou 4
1089  LOGICAL evap_prec ! evaporation de pluie au-dessous de convection
1090  PARAMETER (evap_prec = .TRUE.)
1091  REAL coef_eva
1092  PARAMETER (coef_eva = 1.0E-05)
1093  LOGICAL new_deh ! nouvelle facon de calculer dH
1094  PARAMETER (new_deh = .FALSE.)
1095  REAL t_coup
1096  PARAMETER (t_coup = 234.0)
1097  LOGICAL old_tau ! implique precipitation nulle
1098  PARAMETER (old_tau = .FALSE.)
1099  REAL toliq(klon) ! rapport entre l'eau nuageuse et l'eau precipitante
1100  REAL dpmin, tomax !Epaisseur faible, rapport eau liquide plus grande
1101  PARAMETER (dpmin = 0.15, tomax = 0.97)
1102  REAL dpmax, tomin !Epaisseur grande, rapport eau liquide plus faible
1103  PARAMETER (dpmax = 0.30, tomin = 0.05)
1104  REAL deep_sig, deep_to ! au dela de deep_sig, utiliser deep_to
1105  PARAMETER (deep_sig = 0.50, deep_to = 0.05)
1106
1107  ! Variables locales:
1108
1109  INTEGER nexpo
1110  LOGICAL nuage(klon)
1111  INTEGER i, k, kbmin, kbmax, khmax
1112  REAL ztotal(klon, klev), zdeh(klon, klev)
1113  REAL zgz(klon, klev)
1114  REAL zqs(klon, klev)
1115  REAL zdqs(klon, klev)
1116  REAL ztemp(klon, klev)
1117  REAL zpres(klon, klev)
1118  REAL zconv(klon) ! convergence d'humidite
1119  REAL zvirt(klon) ! convergence virtuelle d'humidite
1120  REAL zfrac(klon) ! fraction convective
1121  INTEGER kb(klon), kh(klon)
1122
1123  REAL zcond(klon), zvapo(klon), zrapp(klon)
1124  REAL zrfl(klon), zrfln, zqev, zqevt
1125  REAL zdelta, zcvm5, zcor
1126  REAL zvar
1127
1128  LOGICAL appel1er
1129  SAVE appel1er
1130  !$OMP THREADPRIVATE(appel1er)
1131
1132  ! Fonctions thermodynamiques
1133
1134  DATA appel1er/.TRUE./
1135
1136  IF (appel1er) THEN
1137    PRINT *, 'conkuo, calcfcl:', calcfcl
1138    IF (.NOT. calcfcl) PRINT *, 'conkuo, ldepar:', ldepar
1139    PRINT *, 'conkuo, opt_cld:', opt_cld
1140    PRINT *, 'conkuo, evap_prec:', evap_prec
1141    PRINT *, 'conkuo, new_deh:', new_deh
1142    appel1er = .FALSE.
1143  END IF
1144
1145  ! Initialiser les sorties a zero
1146
1147  DO k = 1, klev
1148    DO i = 1, klon
1149      d_q(i, k) = 0.0
1150      d_t(i, k) = 0.0
1151      d_ql(i, k) = 0.0
1152      rneb(i, k) = 0.0
1153    END DO
1154  END DO
1155  DO i = 1, klon
1156    rain(i) = 0.0
1157    snow(i) = 0.0
1158    ibas(i) = 0
1159    itop(i) = 0
1160  END DO
1161
1162  ! Calculer la vapeur d'eau saturante Qs et sa derive L/Cp * dQs/dT
1163
1164  DO k = 1, klev
1165    DO i = 1, klon
1166      IF (thermcep) THEN
1167        zdelta = max(0., sign(1., rtt - t(i, k)))
1168        zcvm5 = r5les * rlvtt * (1. - zdelta) + zdelta * r5ies * rlstt
1169        zcvm5 = zcvm5 / rcpd / (1.0 + rvtmp2 * q(i, k))
1170        zqs(i, k) = r2es * foeew(t(i, k), zdelta) / pplay(i, k)
1171        zqs(i, k) = min(0.5, zqs(i, k))
1172        zcor = 1. / (1. - retv * zqs(i, k))
1173        zqs(i, k) = zqs(i, k) * zcor
1174        zdqs(i, k) = foede(t(i, k), zdelta, zcvm5, zqs(i, k), zcor)
1175      ELSE
1176        IF (t(i, k)<t_coup) THEN
1177          zqs(i, k) = qsats(t(i, k)) / pplay(i, k)
1178          zdqs(i, k) = dqsats(t(i, k), zqs(i, k))
1179        ELSE
1180          zqs(i, k) = qsatl(t(i, k)) / pplay(i, k)
1181          zdqs(i, k) = dqsatl(t(i, k), zqs(i, k))
1182        END IF
1183      END IF
1184    END DO
1185  END DO
1186
1187  ! Calculer gz (energie potentielle)
1188
1189  DO i = 1, klon
1190    zgz(i, 1) = rd * t(i, 1) / (0.5 * (paprs(i, 1) + pplay(i, &
1191            1))) * (paprs(i, 1) - pplay(i, 1))
1192  END DO
1193  DO k = 2, klev
1194    DO i = 1, klon
1195      zgz(i, k) = zgz(i, k - 1) + rd * 0.5 * (t(i, k - 1) + t(i, k)) / paprs(i, k) * (pplay(i &
1196              , k - 1) - pplay(i, k))
1197    END DO
1198  END DO
1199
1200  ! Calculer l'energie statique humide saturee (Cp*T + gz + L*Qs)
1201
1202  DO k = 1, klev
1203    DO i = 1, klon
1204      ztotal(i, k) = rcpd * t(i, k) + rlvtt * zqs(i, k) + zgz(i, k)
1205    END DO
1206  END DO
1207
1208  ! Determiner le niveau de depart et calculer la difference de
1209  ! l'energie statique humide saturee (ztotal) entre la couche
1210  ! de depart et chaque couche au-dessus.
1211
1212  IF (calcfcl) THEN
1213    DO k = 1, klev
1214      DO i = 1, klon
1215        zpres(i, k) = pplay(i, k)
1216        ztemp(i, k) = t(i, k)
1217      END DO
1218    END DO
1219    CALL kuofcl(ztemp, q, zgz, zpres, ldcum, kb)
1220    DO i = 1, klon
1221      IF (ldcum(i)) THEN
1222        k = kb(i)
1223        IF (new_deh) THEN
1224          zdeh(i, k) = ztotal(i, k - 1) - ztotal(i, k)
1225        ELSE
1226          zdeh(i, k) = rcpd * (t(i, k - 1) - t(i, k)) - rd * 0.5 * (t(i, k - 1) + t(i, k)) / &
1227                  paprs(i, k) * (pplay(i, k - 1) - pplay(i, k)) + &
1228                  rlvtt * (zqs(i, k - 1) - zqs(i, k))
1229        END IF
1230        zdeh(i, k) = zdeh(i, k) * 0.5
1231      END IF
1232    END DO
1233    DO k = 1, klev
1234      DO i = 1, klon
1235        IF (ldcum(i) .AND. k>=(kb(i) + 1)) THEN
1236          IF (new_deh) THEN
1237            zdeh(i, k) = zdeh(i, k - 1) + (ztotal(i, k - 1) - ztotal(i, k))
1238          ELSE
1239            zdeh(i, k) = zdeh(i, k - 1) + rcpd * (t(i, k - 1) - t(i, k)) - &
1240                    rd * 0.5 * (t(i, k - 1) + t(i, k)) / paprs(i, k) * &
1241                            (pplay(i, k - 1) - pplay(i, k)) + rlvtt * (zqs(i, k - 1) - zqs(i, k))
1242          END IF
1243        END IF
1244      END DO
1245    END DO
1246  ELSE
1247    DO i = 1, klon
1248      k = ldepar
1249      kb(i) = ldepar
1250      ldcum(i) = .TRUE.
1251      IF (new_deh) THEN
1252        zdeh(i, k) = ztotal(i, k - 1) - ztotal(i, k)
1253      ELSE
1254        zdeh(i, k) = rcpd * (t(i, k - 1) - t(i, k)) - rd * 0.5 * (t(i, k - 1) + t(i, k)) / paprs(&
1255                i, k) * (pplay(i, k - 1) - pplay(i, k)) + rlvtt * (zqs(i, k - 1) - zqs(i, k))
1256      END IF
1257      zdeh(i, k) = zdeh(i, k) * 0.5
1258    END DO
1259    DO k = ldepar + 1, klev
1260      DO i = 1, klon
1261        IF (new_deh) THEN
1262          zdeh(i, k) = zdeh(i, k - 1) + (ztotal(i, k - 1) - ztotal(i, k))
1263        ELSE
1264          zdeh(i, k) = zdeh(i, k - 1) + rcpd * (t(i, k - 1) - t(i, k)) - &
1265                  rd * 0.5 * (t(i, k - 1) + t(i, k)) / paprs(i, k) * (pplay(i, k - 1) - pplay(i, k)) + &
1266                  rlvtt * (zqs(i, k - 1) - zqs(i, k))
1267        END IF
1268      END DO
1269    END DO
1270  END IF
1271
1272  ! -----Chercher le sommet du nuage
1273  ! -----Calculer la convergence de l'humidite (en kg/m**2 a un facteur
1274  ! -----psolpa/RG pres) du bas jusqu'au sommet du nuage.
1275  ! -----Calculer la convergence virtuelle pour que toute la maille
1276  ! -----deviennt nuageuse (du bas jusqu'au sommet du nuage)
1277
1278  DO i = 1, klon
1279    nuage(i) = .TRUE.
1280    zconv(i) = 0.0
1281    zvirt(i) = 0.0
1282    kh(i) = -999
1283  END DO
1284  DO k = 1, klev
1285    DO i = 1, klon
1286      IF (k>=kb(i) .AND. ldcum(i)) THEN
1287        nuage(i) = nuage(i) .AND. zdeh(i, k) > 0.0
1288        IF (nuage(i)) THEN
1289          kh(i) = k
1290          zconv(i) = zconv(i) + conv_q(i, k) * dtime * (paprs(i, k) - paprs(i, k + 1))
1291          zvirt(i) = zvirt(i) + (zdeh(i, k) / rlvtt + zqs(i, k) - q(i, k)) * (paprs(i, k) &
1292                  - paprs(i, k + 1))
1293        END IF
1294      END IF
1295    END DO
1296  END DO
1297
1298  DO i = 1, klon
1299    todo(i) = ldcum(i) .AND. kh(i) > kb(i) .AND. zconv(i) > 0.0
1300  END DO
1301
1302  kbmin = klev
1303  kbmax = 0
1304  khmax = 0
1305  DO i = 1, klon
1306    IF (todo(i)) THEN
1307      kbmin = min(kbmin, kb(i))
1308      kbmax = max(kbmax, kb(i))
1309      khmax = max(khmax, kh(i))
1310    END IF
1311  END DO
1312
1313  ! -----Calculer la surface couverte par le nuage
1314
1315  DO i = 1, klon
1316    IF (todo(i)) THEN
1317      zfrac(i) = max(0.0, min(zconv(i) / zvirt(i), 1.0))
1318    END IF
1319  END DO
1320
1321  ! -----Calculs essentiels:
1322
1323  DO i = 1, klon
1324    IF (todo(i)) THEN
1325      zcond(i) = 0.0
1326    END IF
1327  END DO
1328  DO k = kbmin, khmax
1329    DO i = 1, klon
1330      IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
1331        zvar = zdeh(i, k) / (1. + zdqs(i, k))
1332        d_t(i, k) = zvar * zfrac(i) / rcpd
1333        d_q(i, k) = (zvar * zdqs(i, k) / rlvtt + zqs(i, k) - q(i, k)) * zfrac(i) - &
1334                conv_q(i, k) * dtime
1335        zcond(i) = zcond(i) - d_q(i, k) * (paprs(i, k) - paprs(i, k + 1)) / rg
1336        rneb(i, k) = zfrac(i)
1337      END IF
1338    END DO
1339  END DO
1340
1341  DO i = 1, klon
1342    IF (todo(i) .AND. zcond(i)<0.0) THEN
1343      PRINT *, 'WARNING: cond. negative (Kuo) ', i, kb(i), kh(i), zcond(i)
1344      zcond(i) = 0.0
1345      DO k = kb(i), kh(i)
1346        d_t(i, k) = 0.0
1347        d_q(i, k) = 0.0
1348      END DO
1349      todo(i) = .FALSE. ! effort totalement perdu
1350    END IF
1351  END DO
1352
1353  ! =====
1354  ! Une fois que la condensation a lieu, on doit construire un
1355  ! "modele nuageux" pour partager la condensation entre l'eau
1356  ! liquide nuageuse et la precipitation (leur rapport toliq
1357  ! est calcule selon l'epaisseur nuageuse). Je suppose que
1358  ! toliq=tomax quand l'epaisseur nuageuse est inferieure a dpmin,
1359  ! et que toliq=tomin quand l'epaisseur depasse dpmax (interpolation
1360  ! lineaire entre dpmin et dpmax).
1361  ! =====
1362  DO i = 1, klon
1363    IF (todo(i)) THEN
1364      toliq(i) = tomax - ((paprs(i, kb(i)) - paprs(i, kh(i) + 1)) / paprs(i, 1) - dpmin) &
1365              * (tomax - tomin) / (dpmax - dpmin)
1366      toliq(i) = max(tomin, min(tomax, toliq(i)))
1367      IF (pplay(i, kh(i)) / paprs(i, 1)<=deep_sig) toliq(i) = deep_to
1368      IF (old_tau) toliq(i) = 1.0
1369    END IF
1370  END DO
1371  ! =====
1372  ! On doit aussi determiner la distribution verticale de
1373  ! l'eau nuageuse. Plusieurs options sont proposees:
1374
1375  ! (0) La condensation precipite integralement (toliq ne sera
1376  ! pas utilise).
1377  ! (1) L'eau liquide est distribuee entre k1 et k2 et proportionnelle
1378  ! a la vapeur d'eau locale.
1379  ! (2) Elle est distribuee entre k1 et k2 avec une valeur constante.
1380  ! (3) Elle est seulement distribuee aux couches ou la vapeur d'eau
1381  ! est effectivement diminuee pendant le processus d'ajustement.
1382  ! (4) Elle est en fonction (lineaire ou exponentielle) de la
1383  ! distance (epaisseur en pression) avec le niveau k1 (la couche
1384  ! k1 n'aura donc pas d'eau liquide).
1385  ! =====
1386
1387  IF (opt_cld==0) THEN
1388
1389    DO i = 1, klon
1390      IF (todo(i)) zrfl(i) = zcond(i) / dtime
1391    END DO
1392
1393  ELSE IF (opt_cld==1) THEN
1394
1395    DO i = 1, klon
1396      IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de vapeur d'eau
1397    END DO
1398    DO k = kbmin, khmax
1399      DO i = 1, klon
1400        IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
1401          zvapo(i) = zvapo(i) + (q(i, k) + d_q(i, k)) * (paprs(i, k) - paprs(i, k + 1)) / &
1402                  rg
1403        END IF
1404      END DO
1405    END DO
1406    DO i = 1, klon
1407      IF (todo(i)) THEN
1408        zrapp(i) = toliq(i) * zcond(i) / zvapo(i)
1409        zrapp(i) = max(0., min(1., zrapp(i)))
1410      END IF
1411    END DO
1412    DO k = kbmin, khmax
1413      DO i = 1, klon
1414        IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
1415          d_ql(i, k) = zrapp(i) * (q(i, k) + d_q(i, k))
1416        END IF
1417      END DO
1418    END DO
1419    DO i = 1, klon
1420      IF (todo(i)) THEN
1421        zrfl(i) = (1.0 - toliq(i)) * zcond(i) / dtime
1422      END IF
1423    END DO
1424
1425  ELSE IF (opt_cld==2) THEN
1426
1427    DO i = 1, klon
1428      IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de masse
1429    END DO
1430    DO k = kbmin, khmax
1431      DO i = 1, klon
1432        IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
1433          zvapo(i) = zvapo(i) + (paprs(i, k) - paprs(i, k + 1)) / rg
1434        END IF
1435      END DO
1436    END DO
1437    DO k = kbmin, khmax
1438      DO i = 1, klon
1439        IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
1440          d_ql(i, k) = toliq(i) * zcond(i) / zvapo(i)
1441        END IF
1442      END DO
1443    END DO
1444    DO i = 1, klon
1445      IF (todo(i)) THEN
1446        zrfl(i) = (1.0 - toliq(i)) * zcond(i) / dtime
1447      END IF
1448    END DO
1449
1450  ELSE IF (opt_cld==3) THEN
1451
1452    DO i = 1, klon
1453      IF (todo(i)) THEN
1454        zvapo(i) = 0.0 ! quantite de l'eau strictement condensee
1455      END IF
1456    END DO
1457    DO k = kbmin, khmax
1458      DO i = 1, klon
1459        IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
1460          zvapo(i) = zvapo(i) + max(0.0, -d_q(i, k)) * (paprs(i, k) - paprs(i, k + 1)) &
1461                  / rg
1462        END IF
1463      END DO
1464    END DO
1465    DO k = kbmin, khmax
1466      DO i = 1, klon
1467        IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i) .AND. zvapo(i)>0.0) THEN
1468          d_ql(i, k) = d_ql(i, k) + toliq(i) * zcond(i) / zvapo(i) * max(0.0, -d_q(&
1469                  i, k))
1470        END IF
1471      END DO
1472    END DO
1473    DO i = 1, klon
1474      IF (todo(i)) THEN
1475        zrfl(i) = (1.0 - toliq(i)) * zcond(i) / dtime
1476      END IF
1477    END DO
1478
1479  ELSE IF (opt_cld==4) THEN
1480
1481    nexpo = 3
1482    ! cc         nexpo = 1 ! distribution lineaire
1483
1484    DO i = 1, klon
1485      IF (todo(i)) THEN
1486        zvapo(i) = 0.0 ! quantite integrale de masse (avec ponderation)
1487      END IF
1488    END DO
1489    DO k = kbmin, khmax
1490      DO i = 1, klon
1491        IF (todo(i) .AND. k>=(kb(i) + 1) .AND. k<=kh(i)) THEN
1492          zvapo(i) = zvapo(i) + (paprs(i, k) - paprs(i, k + 1)) / rg * (pplay(i, kb(i)) - &
1493                  pplay(i, k))**nexpo
1494        END IF
1495      END DO
1496    END DO
1497    DO k = kbmin, khmax
1498      DO i = 1, klon
1499        IF (todo(i) .AND. k>=(kb(i) + 1) .AND. k<=kh(i)) THEN
1500          d_ql(i, k) = d_ql(i, k) + toliq(i) * zcond(i) / zvapo(i) * (pplay(i, kb(i) &
1501                  ) - pplay(i, k))**nexpo
1502        END IF
1503      END DO
1504    END DO
1505    DO i = 1, klon
1506      IF (todo(i)) THEN
1507        zrfl(i) = (1.0 - toliq(i)) * zcond(i) / dtime
1508      END IF
1509    END DO
1510
1511  ELSE ! valeur non-prevue pour opt_cld
1512
1513    PRINT *, 'opt_cld est faux:', opt_cld
1514    CALL abort
1515
1516  END IF ! fin de opt_cld
1517
1518  ! L'eau precipitante peut etre re-evaporee:
1519
1520  IF (evap_prec .AND. kbmax>=2) THEN
1521    DO k = kbmax, 1, -1
1522      DO i = 1, klon
1523        IF (todo(i) .AND. k<=(kb(i) - 1) .AND. zrfl(i)>0.0) THEN
1524          zqev = max(0.0, (zqs(i, k) - q(i, k)) * zfrac(i))
1525          zqevt = coef_eva * (1.0 - q(i, k) / zqs(i, k)) * sqrt(zrfl(i)) * &
1526                  (paprs(i, k) - paprs(i, k + 1)) / pplay(i, k) * t(i, k) * rd / rg
1527          zqevt = max(0.0, min(zqevt, zrfl(i))) * rg * dtime / &
1528                  (paprs(i, k) - paprs(i, k + 1))
1529          zqev = min(zqev, zqevt)
1530          zrfln = zrfl(i) - zqev * (paprs(i, k) - paprs(i, k + 1)) / rg / dtime
1531          d_q(i, k) = -(zrfln - zrfl(i)) * (rg / (paprs(i, k) - paprs(i, k + 1))) * dtime
1532          d_t(i, k) = (zrfln - zrfl(i)) * (rg / (paprs(i, k) - paprs(i, &
1533                  k + 1))) * dtime * rlvtt / rcpd
1534          zrfl(i) = zrfln
1535        END IF
1536      END DO
1537    END DO
1538  END IF
1539
1540  ! La temperature de la premiere couche determine la pluie ou la neige:
1541
1542  DO i = 1, klon
1543    IF (todo(i)) THEN
1544      IF (t(i, 1)>rtt) THEN
1545        rain(i) = rain(i) + zrfl(i)
1546      ELSE
1547        snow(i) = snow(i) + zrfl(i)
1548      END IF
1549    END IF
1550  END DO
1551
1552END SUBROUTINE conkuo
1553SUBROUTINE kuofcl(pt, pq, pg, pp, ldcum, kcbot)
1554  USE dimphy
1555  USE lmdz_yoethf
1556  USE lmdz_yomcst
1557
1558  IMPLICIT NONE
1559  ! ======================================================================
1560  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19940927
1561  ! adaptation du code de Tiedtke du ECMWF
1562  ! Objet: calculer le niveau de convection libre
1563  ! (FCL: Free Convection Level)
1564  ! ======================================================================
1565  ! Arguments:
1566  ! pt---input-R- temperature (K)
1567  ! pq---input-R- vapeur d'eau (kg/kg)
1568  ! pg---input-R- geopotentiel (g*z ou z est en metre)
1569  ! pp---input-R- pression (Pa)
1570
1571  ! LDCUM---output-L- Y-t-il la convection
1572  ! kcbot---output-I- Niveau du bas de la convection
1573  ! ======================================================================
1574
1575  REAL pt(klon, klev), pq(klon, klev), pg(klon, klev), pp(klon, klev)
1576  INTEGER kcbot(klon)
1577  LOGICAL ldcum(klon)
1578
1579  REAL ztu(klon, klev), zqu(klon, klev), zlu(klon, klev)
1580  REAL zqold(klon), zbuo
1581  INTEGER is, i, k
1582
1583  ! klab=1: on est sous le nuage convectif
1584  ! klab=2: le bas du nuage convectif
1585  ! klab=0: autres couches
1586  INTEGER klab(klon, klev)
1587
1588  ! quand lflag=.TRUE., on est sous le nuage, il faut donc appliquer
1589  ! le processus d'elevation.
1590  LOGICAL lflag(klon)
1591
1592  DO k = 1, klev
1593    DO i = 1, klon
1594      ztu(i, k) = pt(i, k)
1595      zqu(i, k) = pq(i, k)
1596      zlu(i, k) = 0.0
1597      klab(i, k) = 0
1598    END DO
1599  END DO
1600  ! ----------------------------------------------------------------------
1601  DO i = 1, klon
1602    klab(i, 1) = 1
1603    kcbot(i) = 2
1604    ldcum(i) = .FALSE.
1605  END DO
1606
1607  DO k = 2, klev - 1
1608
1609    is = 0
1610    DO i = 1, klon
1611      IF (klab(i, k - 1)==1) is = is + 1
1612      lflag(i) = .FALSE.
1613      IF (klab(i, k - 1)==1) lflag(i) = .TRUE.
1614    END DO
1615    IF (is==0) GO TO 290
1616
1617    ! on eleve le parcel d'air selon l'adiabatique sec
1618
1619    DO i = 1, klon
1620      IF (lflag(i)) THEN
1621        zqu(i, k) = zqu(i, k - 1)
1622        ztu(i, k) = ztu(i, k - 1) + (pg(i, k - 1) - pg(i, k)) / rcpd
1623        zbuo = ztu(i, k) * (1. + retv * zqu(i, k)) - pt(i, k) * (1. + retv * pq(i, k)) + &
1624                0.5
1625        IF (zbuo>0.) klab(i, k) = 1
1626        zqold(i) = zqu(i, k)
1627      END IF
1628    END DO
1629
1630    ! on calcule la condensation eventuelle
1631
1632    CALL adjtq(pp(1, k), ztu(1, k), zqu(1, k), lflag, 1)
1633
1634    ! s'il y a la condensation et la "buoyancy" force est positive
1635    ! c'est bien le bas de la tour de convection
1636
1637    DO i = 1, klon
1638      IF (lflag(i) .AND. zqu(i, k)/=zqold(i)) THEN
1639        klab(i, k) = 2
1640        zlu(i, k) = zlu(i, k) + zqold(i) - zqu(i, k)
1641        zbuo = ztu(i, k) * (1. + retv * zqu(i, k)) - pt(i, k) * (1. + retv * pq(i, k)) + &
1642                0.5
1643        IF (zbuo>0.) THEN
1644          kcbot(i) = k
1645          ldcum(i) = .TRUE.
1646        END IF
1647      END IF
1648    END DO
1649
1650  290 END DO
1651
1652END SUBROUTINE kuofcl
1653SUBROUTINE adjtq(pp, pt, pq, ldflag, kcall)
1654  USE dimphy
1655  USE lmdz_yoethf
1656
1657  USE lmdz_yomcst
1658
1659  IMPLICIT NONE
1660 INCLUDE "FCTTRE.h"
1661  ! ======================================================================
1662  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19940927
1663  ! adaptation du code de Tiedtke du ECMWF
1664  ! Objet: ajustement entre T et Q
1665  ! ======================================================================
1666  ! Arguments:
1667  ! pp---input-R- pression (Pa)
1668  ! pt---input/output-R- temperature (K)
1669  ! pq---input/output-R- vapeur d'eau (kg/kg)
1670  ! ======================================================================
1671  ! TO PRODUCE T,Q AND L VALUES FOR CLOUD ASCENT
1672
1673  ! NOTE: INPUT PARAMETER KCALL DEFINES CALCULATION AS
1674  ! KCALL=0    ENV. T AND QS IN*CUINI*
1675  ! KCALL=1  CONDENSATION IN UPDRAFTS  (E.G. CUBASE, CUASC)
1676  ! KCALL=2  EVAPORATION IN DOWNDRAFTS (E.G. CUDLFS,CUDDRAF)
1677
1678  REAL pt(klon), pq(klon), pp(klon)
1679  LOGICAL ldflag(klon)
1680  INTEGER kcall
1681
1682  REAL t_coup
1683  PARAMETER (t_coup = 234.0)
1684
1685  REAL zcond(klon), zcond1
1686  REAL zdelta, zcvm5, zldcp, zqsat, zcor, zdqsat
1687  INTEGER is, i
1688
1689  DO i = 1, klon
1690    zcond(i) = 0.0
1691  END DO
1692
1693  DO i = 1, klon
1694    IF (ldflag(i)) THEN
1695      zdelta = max(0., sign(1., rtt - pt(i)))
1696      zldcp = rlvtt * (1. - zdelta) + zdelta * rlstt
1697      zldcp = zldcp / rcpd / (1.0 + rvtmp2 * pq(i))
1698      IF (thermcep) THEN
1699        zcvm5 = r5les * rlvtt * (1. - zdelta) + zdelta * r5ies * rlstt
1700        zcvm5 = zcvm5 / rcpd / (1.0 + rvtmp2 * pq(i))
1701        zqsat = r2es * foeew(pt(i), zdelta) / pp(i)
1702        zqsat = min(0.5, zqsat)
1703        zcor = 1. / (1. - retv * zqsat)
1704        zqsat = zqsat * zcor
1705        zdqsat = foede(pt(i), zdelta, zcvm5, zqsat, zcor)
1706      ELSE
1707        IF (pt(i)<t_coup) THEN
1708          zqsat = qsats(pt(i)) / pp(i)
1709          zdqsat = dqsats(pt(i), zqsat)
1710        ELSE
1711          zqsat = qsatl(pt(i)) / pp(i)
1712          zdqsat = dqsatl(pt(i), zqsat)
1713        END IF
1714      END IF
1715      zcond(i) = (pq(i) - zqsat) / (1. + zdqsat)
1716      IF (kcall==1) zcond(i) = max(zcond(i), 0.)
1717      IF (kcall==2) zcond(i) = min(zcond(i), 0.)
1718      pt(i) = pt(i) + zldcp * zcond(i)
1719      pq(i) = pq(i) - zcond(i)
1720    END IF
1721  END DO
1722
1723  is = 0
1724  DO i = 1, klon
1725    IF (zcond(i)/=0.) is = is + 1
1726  END DO
1727  IF (is==0) GO TO 230
1728
1729  DO i = 1, klon
1730    IF (ldflag(i) .AND. zcond(i)/=0.) THEN
1731      zdelta = max(0., sign(1., rtt - pt(i)))
1732      zldcp = rlvtt * (1. - zdelta) + zdelta * rlstt
1733      zldcp = zldcp / rcpd / (1.0 + rvtmp2 * pq(i))
1734      IF (thermcep) THEN
1735        zcvm5 = r5les * rlvtt * (1. - zdelta) + zdelta * r5ies * rlstt
1736        zcvm5 = zcvm5 / rcpd / (1.0 + rvtmp2 * pq(i))
1737        zqsat = r2es * foeew(pt(i), zdelta) / pp(i)
1738        zqsat = min(0.5, zqsat)
1739        zcor = 1. / (1. - retv * zqsat)
1740        zqsat = zqsat * zcor
1741        zdqsat = foede(pt(i), zdelta, zcvm5, zqsat, zcor)
1742      ELSE
1743        IF (pt(i)<t_coup) THEN
1744          zqsat = qsats(pt(i)) / pp(i)
1745          zdqsat = dqsats(pt(i), zqsat)
1746        ELSE
1747          zqsat = qsatl(pt(i)) / pp(i)
1748          zdqsat = dqsatl(pt(i), zqsat)
1749        END IF
1750      END IF
1751      zcond1 = (pq(i) - zqsat) / (1. + zdqsat)
1752      pt(i) = pt(i) + zldcp * zcond1
1753      pq(i) = pq(i) - zcond1
1754    END IF
1755  END DO
1756
1757  230 CONTINUE
1758
1759END SUBROUTINE adjtq
1760SUBROUTINE fiajh(dtime, paprs, pplay, t, q, d_t, d_q, d_ql, rneb, rain, snow, &
1761        ibas, itop)
1762  USE dimphy
1763  USE lmdz_yoethf
1764
1765  USE lmdz_yomcst
1766
1767  IMPLICIT NONE
1768 INCLUDE "FCTTRE.h"
1769
1770  ! Ajustement humide (Schema de convection de Manabe)
1771  ! Arguments:
1772
1773  REAL dtime ! intervalle du temps (s)
1774  REAL t(klon, klev) ! temperature (K)
1775  REAL q(klon, klev) ! humidite specifique (kg/kg)
1776  REAL paprs(klon, klev + 1) ! pression a inter-couche (Pa)
1777  REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
1778
1779  REAL d_t(klon, klev) ! incrementation pour la temperature
1780  REAL d_q(klon, klev) ! incrementation pour vapeur d'eau
1781  REAL d_ql(klon, klev) ! incrementation pour l'eau liquide
1782  REAL rneb(klon, klev) ! fraction nuageuse
1783
1784  REAL rain(klon) ! variable non utilisee
1785  REAL snow(klon) ! variable non utilisee
1786  INTEGER ibas(klon) ! variable non utilisee
1787  INTEGER itop(klon) ! variable non utilisee
1788
1789  REAL t_coup
1790  PARAMETER (t_coup = 234.0)
1791  REAL seuil_vap
1792  PARAMETER (seuil_vap = 1.0E-10)
1793
1794  ! Variables locales:
1795
1796  INTEGER i, k
1797  INTEGER k1, k1p, k2, k2p
1798  LOGICAL itest(klon)
1799  REAL delta_q(klon, klev)
1800  REAL cp_new_t(klev)
1801  REAL cp_delta_t(klev)
1802  REAL new_qb(klev)
1803  REAL v_cptj(klev), v_cptjk1, v_ssig
1804  REAL v_cptt(klon, klev), v_p, v_t
1805  REAL v_qs(klon, klev), v_qsd(klon, klev)
1806  REAL zq1(klon), zq2(klon)
1807  REAL gamcpdz(klon, 2:klev)
1808  REAL zdp, zdpm
1809
1810  REAL zsat ! sur-saturation
1811  REAL zflo ! flotabilite
1812
1813  REAL local_q(klon, klev), local_t(klon, klev)
1814
1815  REAL zdelta, zcor, zcvm5
1816
1817  DO k = 1, klev
1818    DO i = 1, klon
1819      local_q(i, k) = q(i, k)
1820      local_t(i, k) = t(i, k)
1821      rneb(i, k) = 0.0
1822      d_ql(i, k) = 0.0
1823      d_t(i, k) = 0.0
1824      d_q(i, k) = 0.0
1825    END DO
1826  END DO
1827  DO i = 1, klon
1828    rain(i) = 0.0
1829    snow(i) = 0.0
1830    ibas(i) = 0
1831    itop(i) = 0
1832  END DO
1833
1834  ! Calculer v_qs et v_qsd:
1835
1836  DO k = 1, klev
1837    DO i = 1, klon
1838      v_cptt(i, k) = rcpd * local_t(i, k)
1839      v_t = local_t(i, k)
1840      v_p = pplay(i, k)
1841
1842      IF (thermcep) THEN
1843        zdelta = max(0., sign(1., rtt - v_t))
1844        zcvm5 = r5les * rlvtt * (1. - zdelta) + zdelta * r5ies * rlstt
1845        zcvm5 = zcvm5 / rcpd / (1.0 + rvtmp2 * local_q(i, k))
1846        v_qs(i, k) = r2es * foeew(v_t, zdelta) / v_p
1847        v_qs(i, k) = min(0.5, v_qs(i, k))
1848        zcor = 1. / (1. - retv * v_qs(i, k))
1849        v_qs(i, k) = v_qs(i, k) * zcor
1850        v_qsd(i, k) = foede(v_t, zdelta, zcvm5, v_qs(i, k), zcor)
1851      ELSE
1852        IF (v_t<t_coup) THEN
1853          v_qs(i, k) = qsats(v_t) / v_p
1854          v_qsd(i, k) = dqsats(v_t, v_qs(i, k))
1855        ELSE
1856          v_qs(i, k) = qsatl(v_t) / v_p
1857          v_qsd(i, k) = dqsatl(v_t, v_qs(i, k))
1858        END IF
1859      END IF
1860    END DO
1861  END DO
1862
1863  ! Calculer Gamma * Cp * dz: (gamm est le gradient critique)
1864
1865  DO k = 2, klev
1866    DO i = 1, klon
1867      zdp = paprs(i, k) - paprs(i, k + 1)
1868      zdpm = paprs(i, k - 1) - paprs(i, k)
1869      gamcpdz(i, k) = ((rd / rcpd / (zdpm + zdp) * (v_cptt(i, k - 1) * zdpm + &
1870              v_cptt(i, k) * zdp) + rlvtt / (zdpm + zdp) * (v_qs(i, k - 1) * zdpm + &
1871              v_qs(i, k) * zdp)) * (pplay(i, k - 1) - pplay(i, k)) / paprs(i, k)) / (1.0 + (v_qsd(i, &
1872              k - 1) * zdpm + v_qsd(i, k) * zdp) / (zdpm + zdp))
1873    END DO
1874  END DO
1875
1876  ! ------------------------------------ modification des profils instables
1877  DO i = 1, klon
1878    itest(i) = .FALSE.
1879
1880    k1 = 0
1881    k2 = 1
1882
1883    810 CONTINUE ! chercher k1, le bas de la colonne
1884    k2 = k2 + 1
1885    IF (k2>klev) GO TO 9999
1886    zflo = v_cptt(i, k2 - 1) - v_cptt(i, k2) - gamcpdz(i, k2)
1887    zsat = (local_q(i, k2 - 1) - v_qs(i, k2 - 1)) * (paprs(i, k2 - 1) - paprs(i, k2)) + &
1888            (local_q(i, k2) - v_qs(i, k2)) * (paprs(i, k2) - paprs(i, k2 + 1))
1889    IF (zflo<=0.0 .OR. zsat<=0.0) GO TO 810
1890    k1 = k2 - 1
1891    itest(i) = .TRUE.
1892
1893    820 CONTINUE ! chercher k2, le haut de la colonne
1894    IF (k2==klev) GO TO 821
1895    k2p = k2 + 1
1896    zsat = zsat + (paprs(i, k2p) - paprs(i, k2p + 1)) * (local_q(i, k2p) - v_qs(i, k2p))
1897    zflo = v_cptt(i, k2p - 1) - v_cptt(i, k2p) - gamcpdz(i, k2p)
1898    IF (zflo<=0.0 .OR. zsat<=0.0) GO TO 821
1899    k2 = k2p
1900    GO TO 820
1901    821 CONTINUE
1902
1903    ! ------------------------------------------------------ ajustement local
1904    830 CONTINUE ! ajustement proprement dit
1905    v_cptj(k1) = 0.0
1906    zdp = paprs(i, k1) - paprs(i, k1 + 1)
1907    v_cptjk1 = ((1.0 + v_qsd(i, k1)) * (v_cptt(i, k1) + v_cptj(k1)) + rlvtt * (local_q(i, &
1908            k1) - v_qs(i, k1))) * zdp
1909    v_ssig = zdp * (1.0 + v_qsd(i, k1))
1910
1911    k1p = k1 + 1
1912    DO k = k1p, k2
1913      zdp = paprs(i, k) - paprs(i, k + 1)
1914      v_cptj(k) = v_cptj(k - 1) + gamcpdz(i, k)
1915      v_cptjk1 = v_cptjk1 + zdp * ((1.0 + v_qsd(i, k)) * (v_cptt(i, &
1916              k) + v_cptj(k)) + rlvtt * (local_q(i, k) - v_qs(i, k)))
1917      v_ssig = v_ssig + zdp * (1.0 + v_qsd(i, k))
1918    END DO
1919
1920    DO k = k1, k2
1921      cp_new_t(k) = v_cptjk1 / v_ssig - v_cptj(k)
1922      cp_delta_t(k) = cp_new_t(k) - v_cptt(i, k)
1923      new_qb(k) = v_qs(i, k) + v_qsd(i, k) * cp_delta_t(k) / rlvtt
1924      local_q(i, k) = new_qb(k)
1925      local_t(i, k) = cp_new_t(k) / rcpd
1926    END DO
1927
1928    ! --------------------------------------------------- sondage vers le bas
1929    ! -- on redefinit les variables prognostiques dans
1930    ! -- la colonne qui vient d'etre ajustee
1931
1932    DO k = k1, k2
1933      v_cptt(i, k) = rcpd * local_t(i, k)
1934      v_t = local_t(i, k)
1935      v_p = pplay(i, k)
1936
1937      IF (thermcep) THEN
1938        zdelta = max(0., sign(1., rtt - v_t))
1939        zcvm5 = r5les * rlvtt * (1. - zdelta) + zdelta * r5ies * rlstt
1940        zcvm5 = zcvm5 / rcpd / (1.0 + rvtmp2 * local_q(i, k))
1941        v_qs(i, k) = r2es * foeew(v_t, zdelta) / v_p
1942        v_qs(i, k) = min(0.5, v_qs(i, k))
1943        zcor = 1. / (1. - retv * v_qs(i, k))
1944        v_qs(i, k) = v_qs(i, k) * zcor
1945        v_qsd(i, k) = foede(v_t, zdelta, zcvm5, v_qs(i, k), zcor)
1946      ELSE
1947        IF (v_t<t_coup) THEN
1948          v_qs(i, k) = qsats(v_t) / v_p
1949          v_qsd(i, k) = dqsats(v_t, v_qs(i, k))
1950        ELSE
1951          v_qs(i, k) = qsatl(v_t) / v_p
1952          v_qsd(i, k) = dqsatl(v_t, v_qs(i, k))
1953        END IF
1954      END IF
1955    END DO
1956    DO k = 2, klev
1957      zdpm = paprs(i, k - 1) - paprs(i, k)
1958      zdp = paprs(i, k) - paprs(i, k + 1)
1959      gamcpdz(i, k) = ((rd / rcpd / (zdpm + zdp) * (v_cptt(i, k - 1) * zdpm + &
1960              v_cptt(i, k) * zdp) + rlvtt / (zdpm + zdp) * (v_qs(i, k - 1) * zdpm + &
1961              v_qs(i, k) * zdp)) * (pplay(i, k - 1) - pplay(i, k)) / paprs(i, k)) / (1.0 + (v_qsd(i, &
1962              k - 1) * zdpm + v_qsd(i, k) * zdp) / (zdpm + zdp))
1963    END DO
1964
1965    ! Verifier si l'on peut etendre la colonne vers le bas
1966
1967    IF (k1==1) GO TO 841 ! extension echouee
1968    zflo = v_cptt(i, k1 - 1) - v_cptt(i, k1) - gamcpdz(i, k1)
1969    zsat = (local_q(i, k1 - 1) - v_qs(i, k1 - 1)) * (paprs(i, k1 - 1) - paprs(i, k1)) + &
1970            (local_q(i, k1) - v_qs(i, k1)) * (paprs(i, k1) - paprs(i, k1 + 1))
1971    IF (zflo<=0.0 .OR. zsat<=0.0) GO TO 841 ! extension echouee
1972
1973    840 CONTINUE
1974    k1 = k1 - 1
1975    IF (k1==1) GO TO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
1976    zsat = zsat + (local_q(i, k1 - 1) - v_qs(i, k1 - 1)) * (paprs(i, k1 - 1) - paprs(i, k1))
1977    zflo = v_cptt(i, k1 - 1) - v_cptt(i, k1) - gamcpdz(i, k1)
1978    IF (zflo>0.0 .AND. zsat>0.0) THEN
1979      GO TO 840
1980    ELSE
1981      GO TO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
1982    END IF
1983    841 CONTINUE
1984
1985    GO TO 810 ! chercher d'autres blocks en haut
1986
1987  9999 END DO ! boucle sur tous les points
1988  ! -----------------------------------------------------------------------
1989
1990  ! Determiner la fraction nuageuse (hypothese: la nebulosite a lieu
1991  ! a l'endroit ou la vapeur d'eau est diminuee par l'ajustement):
1992
1993  DO k = 1, klev
1994    DO i = 1, klon
1995      IF (itest(i)) THEN
1996        delta_q(i, k) = local_q(i, k) - q(i, k)
1997        IF (delta_q(i, k)<0.) rneb(i, k) = 1.0
1998      END IF
1999    END DO
2000  END DO
2001
2002  ! Distribuer l'eau condensee en eau liquide nuageuse (hypothese:
2003  ! l'eau liquide est distribuee aux endroits ou la vapeur d'eau
2004  ! diminue et d'une maniere proportionnelle a cet diminution):
2005
2006  DO i = 1, klon
2007    IF (itest(i)) THEN
2008      zq1(i) = 0.0
2009      zq2(i) = 0.0
2010    END IF
2011  END DO
2012  DO k = 1, klev
2013    DO i = 1, klon
2014      IF (itest(i)) THEN
2015        zdp = paprs(i, k) - paprs(i, k + 1)
2016        zq1(i) = zq1(i) - delta_q(i, k) * zdp
2017        zq2(i) = zq2(i) - min(0.0, delta_q(i, k)) * zdp
2018      END IF
2019    END DO
2020  END DO
2021  DO k = 1, klev
2022    DO i = 1, klon
2023      IF (itest(i)) THEN
2024        IF (zq2(i)/=0.0) d_ql(i, k) = -min(0.0, delta_q(i, k)) * zq1(i) / zq2(i)
2025      END IF
2026    END DO
2027  END DO
2028
2029  DO k = 1, klev
2030    DO i = 1, klon
2031      local_q(i, k) = max(local_q(i, k), seuil_vap)
2032    END DO
2033  END DO
2034
2035  DO k = 1, klev
2036    DO i = 1, klon
2037      d_t(i, k) = local_t(i, k) - t(i, k)
2038      d_q(i, k) = local_q(i, k) - q(i, k)
2039    END DO
2040  END DO
2041
2042END SUBROUTINE fiajh
2043SUBROUTINE fiajc(dtime, paprs, pplay, t, q, conv_q, d_t, d_q, d_ql, rneb, &
2044        rain, snow, ibas, itop)
2045  USE dimphy
2046  USE lmdz_yoethf
2047
2048  USE lmdz_yomcst
2049
2050  IMPLICIT NONE
2051 INCLUDE "FCTTRE.h"
2052
2053  ! Options:
2054
2055  INTEGER plb ! niveau de depart pour la convection
2056  PARAMETER (plb = 4)
2057
2058  ! Mystere: cette option n'est pas innocente pour les resultats !
2059  ! Qui peut resoudre ce mystere ? (Z.X.Li mars 1995)
2060  LOGICAL vector ! calcul vectorise
2061  PARAMETER (vector = .FALSE.)
2062
2063  REAL t_coup
2064  PARAMETER (t_coup = 234.0)
2065
2066  ! Arguments:
2067
2068  REAL q(klon, klev) ! humidite specifique (kg/kg)
2069  REAL t(klon, klev) ! temperature (K)
2070  REAL paprs(klon, klev + 1) ! pression a inter-couche (Pa)
2071  REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
2072  REAL dtime ! intervalle du temps (s)
2073  REAL conv_q(klon, klev) ! taux de convergence de l'humidite
2074  REAL rneb(klon, klev) ! fraction nuageuse
2075  REAL d_q(klon, klev) ! incrementaion pour la vapeur d'eau
2076  REAL d_ql(klon, klev) ! incrementation pour l'eau liquide
2077  REAL d_t(klon, klev) ! incrementation pour la temperature
2078  REAL rain(klon) ! variable non-utilisee
2079  REAL snow(klon) ! variable non-utilisee
2080  INTEGER itop(klon) ! variable non-utilisee
2081  INTEGER ibas(klon) ! variable non-utilisee
2082
2083  INTEGER kh(klon), i, k
2084  LOGICAL nuage(klon), test(klon, klev)
2085  REAL zconv(klon), zdeh(klon, klev), zvirt(klon)
2086  REAL zdqs(klon, klev), zqs(klon, klev)
2087  REAL ztt, zvar, zfrac(klon)
2088  REAL zq1(klon), zq2(klon)
2089  REAL zdelta, zcor, zcvm5
2090
2091  ! Initialiser les sorties:
2092
2093  DO k = 1, klev
2094    DO i = 1, klon
2095      rneb(i, k) = 0.0
2096      d_ql(i, k) = 0.0
2097      d_t(i, k) = 0.0
2098      d_q(i, k) = 0.0
2099    END DO
2100  END DO
2101  DO i = 1, klon
2102    itop(i) = 0
2103    ibas(i) = 0
2104    rain(i) = 0.0
2105    snow(i) = 0.0
2106  END DO
2107
2108  ! Calculer Qs et L/Cp * dQs/dT:
2109
2110  DO k = 1, klev
2111    DO i = 1, klon
2112      ztt = t(i, k)
2113      IF (thermcep) THEN
2114        zdelta = max(0., sign(1., rtt - ztt))
2115        zcvm5 = r5les * rlvtt * (1. - zdelta) + zdelta * r5ies * rlstt
2116        zcvm5 = zcvm5 / rcpd / (1.0 + rvtmp2 * q(i, k))
2117        zqs(i, k) = r2es * foeew(ztt, zdelta) / pplay(i, k)
2118        zqs(i, k) = min(0.5, zqs(i, k))
2119        zcor = 1. / (1. - retv * zqs(i, k))
2120        zqs(i, k) = zqs(i, k) * zcor
2121        zdqs(i, k) = foede(ztt, zdelta, zcvm5, zqs(i, k), zcor)
2122      ELSE
2123        IF (ztt<t_coup) THEN
2124          zqs(i, k) = qsats(ztt) / pplay(i, k)
2125          zdqs(i, k) = dqsats(ztt, zqs(i, k))
2126        ELSE
2127          zqs(i, k) = qsatl(ztt) / pplay(i, k)
2128          zdqs(i, k) = dqsatl(ztt, zqs(i, k))
2129        END IF
2130      END IF
2131    END DO
2132  END DO
2133
2134  ! Determiner la difference de l'energie totale saturee:
2135
2136  DO i = 1, klon
2137    k = plb
2138    zdeh(i, k) = rcpd * (t(i, k - 1) - t(i, k)) - rd * 0.5 * (t(i, k - 1) + t(i, k)) / paprs(i, k &
2139            ) * (pplay(i, k - 1) - pplay(i, k)) + rlvtt * (zqs(i, k - 1) - zqs(i, k))
2140    zdeh(i, k) = zdeh(i, k) * 0.5 ! on prend la moitie
2141  END DO
2142  DO k = plb + 1, klev
2143    DO i = 1, klon
2144      zdeh(i, k) = zdeh(i, k - 1) + rcpd * (t(i, k - 1) - t(i, k)) - &
2145              rd * 0.5 * (t(i, k - 1) + t(i, k)) / paprs(i, k) * (pplay(i, k - 1) - pplay(i, k)) + &
2146              rlvtt * (zqs(i, k - 1) - zqs(i, k))
2147    END DO
2148  END DO
2149
2150  ! Determiner le sommet du nuage selon l'instabilite
2151  ! Calculer les convergences d'humidite (reelle et virtuelle)
2152
2153  DO i = 1, klon
2154    nuage(i) = .TRUE.
2155    zconv(i) = 0.0
2156    zvirt(i) = 0.0
2157    kh(i) = -999
2158  END DO
2159  DO k = plb, klev
2160    DO i = 1, klon
2161      nuage(i) = nuage(i) .AND. zdeh(i, k) > 0.0
2162      IF (nuage(i)) THEN
2163        kh(i) = k
2164        zconv(i) = zconv(i) + conv_q(i, k) * dtime * (paprs(i, k) - paprs(i, k + 1))
2165        zvirt(i) = zvirt(i) + (zdeh(i, k) / rlvtt + zqs(i, k) - q(i, k)) * (paprs(i, k) - &
2166                paprs(i, k + 1))
2167      END IF
2168    END DO
2169  END DO
2170
2171  IF (vector) THEN
2172
2173    DO k = plb, klev
2174      DO i = 1, klon
2175        IF (k<=kh(i) .AND. kh(i)>plb .AND. zconv(i)>0.0) THEN
2176          test(i, k) = .TRUE.
2177          zfrac(i) = max(0.0, min(zconv(i) / zvirt(i), 1.0))
2178        ELSE
2179          test(i, k) = .FALSE.
2180        END IF
2181      END DO
2182    END DO
2183
2184    DO k = plb, klev
2185      DO i = 1, klon
2186        IF (test(i, k)) THEN
2187          zvar = zdeh(i, k) / (1.0 + zdqs(i, k))
2188          d_q(i, k) = (zvar * zdqs(i, k) / rlvtt + zqs(i, k) - q(i, k)) * zfrac(i) - &
2189                  conv_q(i, k) * dtime
2190          d_t(i, k) = zvar * zfrac(i) / rcpd
2191        END IF
2192      END DO
2193    END DO
2194
2195    DO i = 1, klon
2196      zq1(i) = 0.0
2197      zq2(i) = 0.0
2198    END DO
2199    DO k = plb, klev
2200      DO i = 1, klon
2201        IF (test(i, k)) THEN
2202          IF (d_q(i, k)<0.0) rneb(i, k) = zfrac(i)
2203          zq1(i) = zq1(i) - d_q(i, k) * (paprs(i, k) - paprs(i, k + 1))
2204          zq2(i) = zq2(i) - min(0.0, d_q(i, k)) * (paprs(i, k) - paprs(i, k + 1))
2205        END IF
2206      END DO
2207    END DO
2208
2209    DO k = plb, klev
2210      DO i = 1, klon
2211        IF (test(i, k)) THEN
2212          IF (zq2(i)/=0.) d_ql(i, k) = -min(0.0, d_q(i, k)) * zq1(i) / zq2(i)
2213        END IF
2214      END DO
2215    END DO
2216
2217  ELSE ! (.NOT. vector)
2218
2219    DO i = 1, klon
2220      IF (kh(i)>plb .AND. zconv(i)>0.0) THEN
2221        ! cc         IF (kh(i).LE.plb) GOTO 999 ! il n'y a pas d'instabilite
2222        ! cc         IF (zconv(i).LE.0.0) GOTO 999 ! convergence insuffisante
2223        zfrac(i) = max(0.0, min(zconv(i) / zvirt(i), 1.0))
2224        DO k = plb, kh(i)
2225          zvar = zdeh(i, k) / (1.0 + zdqs(i, k))
2226          d_q(i, k) = (zvar * zdqs(i, k) / rlvtt + zqs(i, k) - q(i, k)) * zfrac(i) - &
2227                  conv_q(i, k) * dtime
2228          d_t(i, k) = zvar * zfrac(i) / rcpd
2229        END DO
2230
2231        zq1(i) = 0.0
2232        zq2(i) = 0.0
2233        DO k = plb, kh(i)
2234          IF (d_q(i, k)<0.0) rneb(i, k) = zfrac(i)
2235          zq1(i) = zq1(i) - d_q(i, k) * (paprs(i, k) - paprs(i, k + 1))
2236          zq2(i) = zq2(i) - min(0.0, d_q(i, k)) * (paprs(i, k) - paprs(i, k + 1))
2237        END DO
2238        DO k = plb, kh(i)
2239          IF (zq2(i)/=0.) d_ql(i, k) = -min(0.0, d_q(i, k)) * zq1(i) / zq2(i)
2240        END DO
2241      END IF
2242    END DO
2243
2244  END IF ! fin de teste sur vector
2245
2246END SUBROUTINE fiajc
Note: See TracBrowser for help on using the repository browser.