source: LMDZ6/branches/contrails/libf/phylmd/conlmd.f90 @ 5452

Last change on this file since 5452 was 5285, checked in by abarral, 2 months ago

As discussed internally, remove generic ONLY: ... for new _mod_h modules
  • 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: 63.6 KB
Line 
1
2! $Header$
3
4SUBROUTINE conlmd(dtime, paprs, pplay, t, q, conv_q, d_t, d_q, rain, snow, &
5    ibas, itop)
6  USE dimphy
7  USE yomcst_mod_h
8  USE yoethf_mod_h
9IMPLICIT 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
88  RETURN
89END SUBROUTINE conlmd
90SUBROUTINE conman(dtime, paprs, pplay, t, q, d_t, d_q, d_ql, rneb, rain, &
91    snow, ibas, itop)
92  USE dimphy
93  USE yomcst_mod_h
94  USE yoethf_mod_h
95  IMPLICIT NONE
96  ! ======================================================================
97  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19970324
98  ! Objet: ajustement humide convectif avec la possibilite de faire
99  ! l'ajustement sur une fraction de la maille.
100  ! Methode: On impose une distribution uniforme pour la vapeur d'eau
101  ! au sein d'une maille. On applique la procedure d'ajustement
102  ! successivement a la totalite, 75%, 50%, 25% et 5% de la maille
103  ! jusqu'a ce que l'ajustement a lieu. J'espere que ceci augmente
104  ! les activites convectives et corrige le biais "trop froid et sec"
105  ! du modele.
106  ! ======================================================================
107  REAL dtime ! pas d'integration (s)
108  REAL t(klon, klev) ! temperature (K)
109  REAL q(klon, klev) ! humidite specifique (kg/kg)
110  REAL paprs(klon, klev+1) ! pression inter-couche (Pa)
111  REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
112
113  REAL d_t(klon, klev) ! incrementation temperature
114  REAL d_q(klon, klev) ! incrementation humidite
115  REAL d_ql(klon, klev) ! incrementation eau liquide
116  REAL rneb(klon, klev) ! nebulosite
117  REAL rain(klon) ! pluies (mm/s)
118  REAL snow(klon) ! neige (mm/s)
119  INTEGER ibas(klon) ! niveau du bas
120  INTEGER itop(klon) ! niveau du haut
121
122  LOGICAL afaire(klon) ! .TRUE. implique l'ajustement
123  LOGICAL accompli(klon) ! .TRUE. si l'ajustement est effectif
124
125  INTEGER nb ! nombre de sous-fractions a considere
126  PARAMETER (nb=1)
127  ! cc      PARAMETER (nb=3)
128
129  REAL ratqs ! largeur de la distribution pour vapeur d'eau
130  PARAMETER (ratqs=0.05)
131
132  REAL w_q(klon, klev)
133  REAL w_d_t(klon, klev), w_d_q(klon, klev), w_d_ql(klon, klev)
134  REAL w_rneb(klon, klev)
135  REAL w_rain(klon), w_snow(klon)
136  INTEGER w_ibas(klon), w_itop(klon)
137  REAL zq1, zq2
138  INTEGER i, k, n
139
140  REAL t_coup
141  PARAMETER (t_coup=234.0)
142  REAL zdp1, zdp2
143  REAL zqs1, zqs2, zdqs1, zdqs2
144  REAL zgamdz
145  REAL zflo ! flotabilite
146  REAL zsat ! sur-saturation
147  REAL zdelta, zcor, zcvm5
148  LOGICAL imprim
149
150  INTEGER ncpt
151  SAVE ncpt
152  !$OMP THREADPRIVATE(ncpt)
153  REAL frac(nb) ! valeur de la maille fractionnelle
154  SAVE frac
155  !$OMP THREADPRIVATE(frac)
156  INTEGER opt_cld(nb) ! option pour le modele nuageux
157  SAVE opt_cld
158  !$OMP THREADPRIVATE(opt_cld)
159  LOGICAL appel1er
160  SAVE appel1er
161  !$OMP THREADPRIVATE(appel1er)
162
163  ! Fonctions thermodynamiques:
164
165  include "FCTTRE.h"
166
167  DATA frac/1.0/
168  DATA opt_cld/4/
169  ! cc      DATA frac    / 1.0, 0.50, 0.25/
170  ! cc      DATA opt_cld / 4,   4,    4/
171
172  DATA appel1er/.TRUE./
173  DATA ncpt/0/
174
175  IF (appel1er) THEN
176    PRINT *, 'conman, nb:', nb
177    PRINT *, 'conman, frac:', frac
178    PRINT *, 'conman, opt_cld:', opt_cld
179    appel1er = .FALSE.
180  END IF
181
182  ! Initialiser les sorties a zero:
183
184  DO k = 1, klev
185    DO i = 1, klon
186      d_t(i, k) = 0.0
187      d_q(i, k) = 0.0
188      d_ql(i, k) = 0.0
189      rneb(i, k) = 0.0
190    END DO
191  END DO
192  DO i = 1, klon
193    ibas(i) = klev
194    itop(i) = 1
195    rain(i) = 0.0
196    snow(i) = 0.0
197  END DO
198
199  ! S'il n'y a pas d'instabilite conditionnelle,
200  ! pas la penne de se fatiguer:
201
202  DO i = 1, klon
203    afaire(i) = .FALSE.
204  END DO
205  DO k = 1, klev - 1
206    DO i = 1, klon
207      IF (thermcep) THEN
208        zdelta = max(0., sign(1.,rtt-t(i,k)))
209        zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
210        zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*q(i,k))
211        zqs1 = r2es*foeew(t(i,k), zdelta)/pplay(i, k)
212        zqs1 = min(0.5, zqs1)
213        zcor = 1./(1.-retv*zqs1)
214        zqs1 = zqs1*zcor
215        zdqs1 = foede(t(i,k), zdelta, zcvm5, zqs1, zcor)
216
217        zdelta = max(0., sign(1.,rtt-t(i,k+1)))
218        zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
219        zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*q(i,k+1))
220        zqs2 = r2es*foeew(t(i,k+1), zdelta)/pplay(i, k+1)
221        zqs2 = min(0.5, zqs2)
222        zcor = 1./(1.-retv*zqs2)
223        zqs2 = zqs2*zcor
224        zdqs2 = foede(t(i,k+1), zdelta, zcvm5, zqs2, zcor)
225      ELSE
226        IF (t(i,k)<t_coup) THEN
227          zqs1 = qsats(t(i,k))/pplay(i, k)
228          zdqs1 = dqsats(t(i,k), zqs1)
229
230          zqs2 = qsats(t(i,k+1))/pplay(i, k+1)
231          zdqs2 = dqsats(t(i,k+1), zqs2)
232        ELSE
233          zqs1 = qsatl(t(i,k))/pplay(i, k)
234          zdqs1 = dqsatl(t(i,k), zqs1)
235
236          zqs2 = qsatl(t(i,k+1))/pplay(i, k+1)
237          zdqs2 = dqsatl(t(i,k+1), zqs2)
238        END IF
239      END IF
240      zdp1 = paprs(i, k) - paprs(i, k+1)
241      zdp2 = paprs(i, k+1) - paprs(i, k+2)
242      zgamdz = -(pplay(i,k)-pplay(i,k+1))/paprs(i, k+1)/rcpd*(rd*(t(i, &
243        k)*zdp1+t(i,k+1)*zdp2)/(zdp1+zdp2)+rlvtt*(zqs1*zdp1+zqs2*zdp2)/(zdp1+ &
244        zdp2))/(1.0+(zdqs1*zdp1+zdqs2*zdp2)/(zdp1+zdp2))
245      zflo = t(i, k) + zgamdz - t(i, k+1)
246      zsat = (q(i,k)-zqs1)*zdp1 + (q(i,k+1)-zqs2)*zdp2
247      IF (zflo>0.0) afaire(i) = .TRUE.
248      ! erreur         IF (zflo.GT.0.0 .AND. zsat.GT.0.0) afaire(i) = .TRUE.
249    END DO
250  END DO
251
252  imprim = mod(ncpt, 48) == 0
253  DO n = 1, nb
254
255    DO k = 1, klev
256      DO i = 1, klon
257        IF (afaire(i)) THEN
258          zq1 = q(i, k)*(1.0-ratqs)
259          zq2 = q(i, k)*(1.0+ratqs)
260          w_q(i, k) = zq2 - frac(n)/2.0*(zq2-zq1)
261        END IF
262      END DO
263    END DO
264
265    CALL conmanv(dtime, paprs, pplay, t, w_q, afaire, opt_cld(n), w_d_t, &
266      w_d_q, w_d_ql, w_rneb, w_rain, w_snow, w_ibas, w_itop, accompli, &
267      imprim)
268    DO k = 1, klev
269      DO i = 1, klon
270        IF (afaire(i) .AND. accompli(i)) THEN
271          d_t(i, k) = w_d_t(i, k)*frac(n)
272          d_q(i, k) = w_d_q(i, k)*frac(n)
273          d_ql(i, k) = w_d_ql(i, k)*frac(n)
274          IF (nint(w_rneb(i,k))==1) rneb(i, k) = frac(n)
275        END IF
276      END DO
277    END DO
278    DO i = 1, klon
279      IF (afaire(i) .AND. accompli(i)) THEN
280        rain(i) = w_rain(i)*frac(n)
281        snow(i) = w_snow(i)*frac(n)
282        ibas(i) = min(ibas(i), w_ibas(i))
283        itop(i) = max(itop(i), w_itop(i))
284      END IF
285    END DO
286    DO i = 1, klon
287      IF (afaire(i) .AND. accompli(i)) afaire(i) = .FALSE.
288    END DO
289
290  END DO
291
292  ncpt = ncpt + 1
293
294  RETURN
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 yomcst_mod_h
300  USE yoethf_mod_h
301IMPLICIT NONE
302  ! ======================================================================
303  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
304  ! Objet: ajustement humide (convection proposee par Manabe).
305  ! Pour une colonne verticale, il peut avoir plusieurs blocs
306  ! necessitant l'ajustement. ibas est le bas du plus bas bloc
307  ! et itop est le haut du plus haut bloc
308  ! ======================================================================
309
310
311  ! Arguments:
312
313  REAL dtime ! pas d'integration (s)
314  REAL t(klon, klev) ! temperature (K)
315  REAL q(klon, klev) ! humidite specifique (kg/kg)
316  REAL paprs(klon, klev+1) ! pression inter-couche (Pa)
317  REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
318  INTEGER opt_cld ! comment traiter l'eau liquide
319  LOGICAL afaire(klon) ! .TRUE. si le point est a faire (Input)
320  LOGICAL imprim ! .T. pour imprimer quelques diagnostiques
321
322  REAL d_t(klon, klev) ! incrementation temperature
323  REAL d_q(klon, klev) ! incrementation humidite
324  REAL d_ql(klon, klev) ! incrementation eau liquide
325  REAL rneb(klon, klev) ! nebulosite
326  REAL rain(klon) ! pluies (mm/s)
327  REAL snow(klon) ! neige (mm/s)
328  INTEGER ibas(klon) ! niveau du bas
329  INTEGER itop(klon) ! niveau du haut
330  LOGICAL accompli(klon) ! .TRUE. si l'ajustement a eu lieu (Output)
331
332  ! Quelques options:
333
334  LOGICAL new_top ! re-calculer sommet quand re-ajustement est fait
335  PARAMETER (new_top=.FALSE.)
336  LOGICAL evap_prec ! evaporation de pluie au-dessous de convection
337  PARAMETER (evap_prec=.TRUE.)
338  REAL coef_eva
339  PARAMETER (coef_eva=1.0E-05)
340  REAL t_coup
341  PARAMETER (t_coup=234.0)
342  REAL seuil_vap
343  PARAMETER (seuil_vap=1.0E-10)
344  LOGICAL old_tau ! implique precip nulle, si vrai.
345  PARAMETER (old_tau=.FALSE.)
346  REAL toliq(klon) ! rapport entre l'eau nuageuse et l'eau precipitante
347  REAL dpmin, tomax !Epaisseur faible, rapport eau liquide plus grande
348  PARAMETER (dpmin=0.15, tomax=0.97)
349  REAL dpmax, tomin !Epaisseur grande, rapport eau liquide plus faible
350  PARAMETER (dpmax=0.30, tomin=0.05)
351  REAL deep_sig, deep_to ! au dela de deep_sig, utiliser deep_to
352  PARAMETER (deep_sig=0.50, deep_to=0.05)
353  LOGICAL exigent ! implique un calcul supplementaire pour Qs
354  PARAMETER (exigent=.FALSE.)
355
356  INTEGER kbase
357  PARAMETER (kbase=0)
358
359  ! Variables locales:
360
361  INTEGER nexpo
362  INTEGER i, k, k1min, k1max, k2min, k2max, is
363  REAL zgamdz(klon, klev-1)
364  REAL zt(klon, klev), zq(klon, klev)
365  REAL zqs(klon, klev), zdqs(klon, klev)
366  REAL zqmqsdp(klon, klev)
367  REAL ztnew(klon, klev), zqnew(klon, klev)
368  REAL zcond(klon), zvapo(klon), zrapp(klon)
369  REAL zrfl(klon), zrfln, zqev, zqevt
370  REAL zsat(klon) ! sur-saturation
371  REAL zflo(klon) ! flotabilite
372  REAL za(klon), zb(klon), zc(klon)
373  INTEGER k1(klon), k2(klon)
374  REAL zdelta, zcor, zcvm5
375  REAL delp(klon, klev)
376  LOGICAL possible(klon), todo(klon), etendre(klon)
377  LOGICAL aller(klon), todobis(klon)
378  REAL zalfa
379  INTEGER nbtodo, nbdone
380
381  ! Fonctions thermodynamiques:
382
383  include "FCTTRE.h"
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
484810 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
532820 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
581830 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
1001860 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
1022807   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
1029808   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
1043  RETURN
1044END SUBROUTINE conmanv
1045SUBROUTINE conkuo(dtime, paprs, pplay, t, q, conv_q, d_t, d_q, d_ql, rneb, &
1046    rain, snow, ibas, itop)
1047  USE dimphy
1048  USE yomcst_mod_h
1049  USE yoethf_mod_h
1050IMPLICIT NONE
1051  ! ======================================================================
1052  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
1053  ! Objet: Schema de convection de type Kuo (1965).
1054  ! Cette version du code peut calculer le niveau de depart
1055  ! N.B. version vectorielle (le 6 oct. 1997)
1056  ! ======================================================================
1057
1058
1059  ! Arguments:
1060
1061  REAL dtime ! intervalle du temps (s)
1062  REAL paprs(klon, klev+1) ! pression a inter-couche (Pa)
1063  REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
1064  REAL t(klon, klev) ! temperature (K)
1065  REAL q(klon, klev) ! humidite specifique
1066  REAL conv_q(klon, klev) ! taux de convergence humidite (g/g/s)
1067
1068  REAL d_t(klon, klev) ! incrementation temperature
1069  REAL d_q(klon, klev) ! incrementation humidite
1070  REAL d_ql(klon, klev) ! incrementation eau liquide
1071  REAL rneb(klon, klev) ! nebulosite
1072  REAL rain(klon) ! pluies (mm/s)
1073  REAL snow(klon) ! neige (mm/s)
1074  INTEGER itop(klon) ! niveau du sommet
1075  INTEGER ibas(klon) ! niveau du bas
1076
1077  LOGICAL ldcum(klon) ! convection existe
1078  LOGICAL todo(klon)
1079
1080  ! Quelsques options:
1081
1082  LOGICAL calcfcl ! calculer le niveau de convection libre
1083  PARAMETER (calcfcl=.TRUE.)
1084  INTEGER ldepar ! niveau fixe de convection libre
1085  PARAMETER (ldepar=4)
1086  INTEGER opt_cld ! comment traiter l'eau liquide
1087  PARAMETER (opt_cld=4) ! valeur possible: 0, 1, 2, 3 ou 4
1088  LOGICAL evap_prec ! evaporation de pluie au-dessous de convection
1089  PARAMETER (evap_prec=.TRUE.)
1090  REAL coef_eva
1091  PARAMETER (coef_eva=1.0E-05)
1092  LOGICAL new_deh ! nouvelle facon de calculer dH
1093  PARAMETER (new_deh=.FALSE.)
1094  REAL t_coup
1095  PARAMETER (t_coup=234.0)
1096  LOGICAL old_tau ! implique precipitation nulle
1097  PARAMETER (old_tau=.FALSE.)
1098  REAL toliq(klon) ! rapport entre l'eau nuageuse et l'eau precipitante
1099  REAL dpmin, tomax !Epaisseur faible, rapport eau liquide plus grande
1100  PARAMETER (dpmin=0.15, tomax=0.97)
1101  REAL dpmax, tomin !Epaisseur grande, rapport eau liquide plus faible
1102  PARAMETER (dpmax=0.30, tomin=0.05)
1103  REAL deep_sig, deep_to ! au dela de deep_sig, utiliser deep_to
1104  PARAMETER (deep_sig=0.50, deep_to=0.05)
1105
1106  ! Variables locales:
1107
1108  INTEGER nexpo
1109  LOGICAL nuage(klon)
1110  INTEGER i, k, kbmin, kbmax, khmax
1111  REAL ztotal(klon, klev), zdeh(klon, klev)
1112  REAL zgz(klon, klev)
1113  REAL zqs(klon, klev)
1114  REAL zdqs(klon, klev)
1115  REAL ztemp(klon, klev)
1116  REAL zpres(klon, klev)
1117  REAL zconv(klon) ! convergence d'humidite
1118  REAL zvirt(klon) ! convergence virtuelle d'humidite
1119  REAL zfrac(klon) ! fraction convective
1120  INTEGER kb(klon), kh(klon)
1121
1122  REAL zcond(klon), zvapo(klon), zrapp(klon)
1123  REAL zrfl(klon), zrfln, zqev, zqevt
1124  REAL zdelta, zcvm5, zcor
1125  REAL zvar
1126
1127  LOGICAL appel1er
1128  SAVE appel1er
1129  !$OMP THREADPRIVATE(appel1er)
1130
1131  ! Fonctions thermodynamiques
1132
1133  include "FCTTRE.h"
1134
1135  DATA appel1er/.TRUE./
1136
1137  IF (appel1er) THEN
1138    PRINT *, 'conkuo, calcfcl:', calcfcl
1139    IF (.NOT. calcfcl) PRINT *, 'conkuo, ldepar:', ldepar
1140    PRINT *, 'conkuo, opt_cld:', opt_cld
1141    PRINT *, 'conkuo, evap_prec:', evap_prec
1142    PRINT *, 'conkuo, new_deh:', new_deh
1143    appel1er = .FALSE.
1144  END IF
1145
1146  ! Initialiser les sorties a zero
1147
1148  DO k = 1, klev
1149    DO i = 1, klon
1150      d_q(i, k) = 0.0
1151      d_t(i, k) = 0.0
1152      d_ql(i, k) = 0.0
1153      rneb(i, k) = 0.0
1154    END DO
1155  END DO
1156  DO i = 1, klon
1157    rain(i) = 0.0
1158    snow(i) = 0.0
1159    ibas(i) = 0
1160    itop(i) = 0
1161  END DO
1162
1163  ! Calculer la vapeur d'eau saturante Qs et sa derive L/Cp * dQs/dT
1164
1165  DO k = 1, klev
1166    DO i = 1, klon
1167      IF (thermcep) THEN
1168        zdelta = max(0., sign(1.,rtt-t(i,k)))
1169        zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
1170        zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*q(i,k))
1171        zqs(i, k) = r2es*foeew(t(i,k), zdelta)/pplay(i, k)
1172        zqs(i, k) = min(0.5, zqs(i,k))
1173        zcor = 1./(1.-retv*zqs(i,k))
1174        zqs(i, k) = zqs(i, k)*zcor
1175        zdqs(i, k) = foede(t(i,k), zdelta, zcvm5, zqs(i,k), zcor)
1176      ELSE
1177        IF (t(i,k)<t_coup) THEN
1178          zqs(i, k) = qsats(t(i,k))/pplay(i, k)
1179          zdqs(i, k) = dqsats(t(i,k), zqs(i,k))
1180        ELSE
1181          zqs(i, k) = qsatl(t(i,k))/pplay(i, k)
1182          zdqs(i, k) = dqsatl(t(i,k), zqs(i,k))
1183        END IF
1184      END IF
1185    END DO
1186  END DO
1187
1188  ! Calculer gz (energie potentielle)
1189
1190  DO i = 1, klon
1191    zgz(i, 1) = rd*t(i, 1)/(0.5*(paprs(i,1)+pplay(i, &
1192      1)))*(paprs(i,1)-pplay(i,1))
1193  END DO
1194  DO k = 2, klev
1195    DO i = 1, klon
1196      zgz(i, k) = zgz(i, k-1) + rd*0.5*(t(i,k-1)+t(i,k))/paprs(i, k)*(pplay(i &
1197        ,k-1)-pplay(i,k))
1198    END DO
1199  END DO
1200
1201  ! Calculer l'energie statique humide saturee (Cp*T + gz + L*Qs)
1202
1203  DO k = 1, klev
1204    DO i = 1, klon
1205      ztotal(i, k) = rcpd*t(i, k) + rlvtt*zqs(i, k) + zgz(i, k)
1206    END DO
1207  END DO
1208
1209  ! Determiner le niveau de depart et calculer la difference de
1210  ! l'energie statique humide saturee (ztotal) entre la couche
1211  ! de depart et chaque couche au-dessus.
1212
1213  IF (calcfcl) THEN
1214    DO k = 1, klev
1215      DO i = 1, klon
1216        zpres(i, k) = pplay(i, k)
1217        ztemp(i, k) = t(i, k)
1218      END DO
1219    END DO
1220    CALL kuofcl(ztemp, q, zgz, zpres, ldcum, kb)
1221    DO i = 1, klon
1222      IF (ldcum(i)) THEN
1223        k = kb(i)
1224        IF (new_deh) THEN
1225          zdeh(i, k) = ztotal(i, k-1) - ztotal(i, k)
1226        ELSE
1227          zdeh(i, k) = rcpd*(t(i,k-1)-t(i,k)) - rd*0.5*(t(i,k-1)+t(i,k))/ &
1228            paprs(i, k)*(pplay(i,k-1)-pplay(i,k)) + &
1229            rlvtt*(zqs(i,k-1)-zqs(i,k))
1230        END IF
1231        zdeh(i, k) = zdeh(i, k)*0.5
1232      END IF
1233    END DO
1234    DO k = 1, klev
1235      DO i = 1, klon
1236        IF (ldcum(i) .AND. k>=(kb(i)+1)) THEN
1237          IF (new_deh) THEN
1238            zdeh(i, k) = zdeh(i, k-1) + (ztotal(i,k-1)-ztotal(i,k))
1239          ELSE
1240            zdeh(i, k) = zdeh(i, k-1) + rcpd*(t(i,k-1)-t(i,k)) - &
1241              rd*0.5*(t(i,k-1)+t(i,k))/paprs(i, k)* &
1242              (pplay(i,k-1)-pplay(i,k)) + rlvtt*(zqs(i,k-1)-zqs(i,k))
1243          END IF
1244        END IF
1245      END DO
1246    END DO
1247  ELSE
1248    DO i = 1, klon
1249      k = ldepar
1250      kb(i) = ldepar
1251      ldcum(i) = .TRUE.
1252      IF (new_deh) THEN
1253        zdeh(i, k) = ztotal(i, k-1) - ztotal(i, k)
1254      ELSE
1255        zdeh(i, k) = rcpd*(t(i,k-1)-t(i,k)) - rd*0.5*(t(i,k-1)+t(i,k))/paprs( &
1256          i, k)*(pplay(i,k-1)-pplay(i,k)) + rlvtt*(zqs(i,k-1)-zqs(i,k))
1257      END IF
1258      zdeh(i, k) = zdeh(i, k)*0.5
1259    END DO
1260    DO k = ldepar + 1, klev
1261      DO i = 1, klon
1262        IF (new_deh) THEN
1263          zdeh(i, k) = zdeh(i, k-1) + (ztotal(i,k-1)-ztotal(i,k))
1264        ELSE
1265          zdeh(i, k) = zdeh(i, k-1) + rcpd*(t(i,k-1)-t(i,k)) - &
1266            rd*0.5*(t(i,k-1)+t(i,k))/paprs(i, k)*(pplay(i,k-1)-pplay(i,k)) + &
1267            rlvtt*(zqs(i,k-1)-zqs(i,k))
1268        END IF
1269      END DO
1270    END DO
1271  END IF
1272
1273  ! -----Chercher le sommet du nuage
1274  ! -----Calculer la convergence de l'humidite (en kg/m**2 a un facteur
1275  ! -----psolpa/RG pres) du bas jusqu'au sommet du nuage.
1276  ! -----Calculer la convergence virtuelle pour que toute la maille
1277  ! -----deviennt nuageuse (du bas jusqu'au sommet du nuage)
1278
1279  DO i = 1, klon
1280    nuage(i) = .TRUE.
1281    zconv(i) = 0.0
1282    zvirt(i) = 0.0
1283    kh(i) = -999
1284  END DO
1285  DO k = 1, klev
1286    DO i = 1, klon
1287      IF (k>=kb(i) .AND. ldcum(i)) THEN
1288        nuage(i) = nuage(i) .AND. zdeh(i, k) > 0.0
1289        IF (nuage(i)) THEN
1290          kh(i) = k
1291          zconv(i) = zconv(i) + conv_q(i, k)*dtime*(paprs(i,k)-paprs(i,k+1))
1292          zvirt(i) = zvirt(i) + (zdeh(i,k)/rlvtt+zqs(i,k)-q(i,k))*(paprs(i,k) &
1293            -paprs(i,k+1))
1294        END IF
1295      END IF
1296    END DO
1297  END DO
1298
1299  DO i = 1, klon
1300    todo(i) = ldcum(i) .AND. kh(i) > kb(i) .AND. zconv(i) > 0.0
1301  END DO
1302
1303  kbmin = klev
1304  kbmax = 0
1305  khmax = 0
1306  DO i = 1, klon
1307    IF (todo(i)) THEN
1308      kbmin = min(kbmin, kb(i))
1309      kbmax = max(kbmax, kb(i))
1310      khmax = max(khmax, kh(i))
1311    END IF
1312  END DO
1313
1314  ! -----Calculer la surface couverte par le nuage
1315
1316  DO i = 1, klon
1317    IF (todo(i)) THEN
1318      zfrac(i) = max(0.0, min(zconv(i)/zvirt(i),1.0))
1319    END IF
1320  END DO
1321
1322  ! -----Calculs essentiels:
1323
1324  DO i = 1, klon
1325    IF (todo(i)) THEN
1326      zcond(i) = 0.0
1327    END IF
1328  END DO
1329  DO k = kbmin, khmax
1330    DO i = 1, klon
1331      IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
1332        zvar = zdeh(i, k)/(1.+zdqs(i,k))
1333        d_t(i, k) = zvar*zfrac(i)/rcpd
1334        d_q(i, k) = (zvar*zdqs(i,k)/rlvtt+zqs(i,k)-q(i,k))*zfrac(i) - &
1335          conv_q(i, k)*dtime
1336        zcond(i) = zcond(i) - d_q(i, k)*(paprs(i,k)-paprs(i,k+1))/rg
1337        rneb(i, k) = zfrac(i)
1338      END IF
1339    END DO
1340  END DO
1341
1342  DO i = 1, klon
1343    IF (todo(i) .AND. zcond(i)<0.0) THEN
1344      PRINT *, 'WARNING: cond. negative (Kuo) ', i, kb(i), kh(i), zcond(i)
1345      zcond(i) = 0.0
1346      DO k = kb(i), kh(i)
1347        d_t(i, k) = 0.0
1348        d_q(i, k) = 0.0
1349      END DO
1350      todo(i) = .FALSE. ! effort totalement perdu
1351    END IF
1352  END DO
1353
1354  ! =====
1355  ! Une fois que la condensation a lieu, on doit construire un
1356  ! "modele nuageux" pour partager la condensation entre l'eau
1357  ! liquide nuageuse et la precipitation (leur rapport toliq
1358  ! est calcule selon l'epaisseur nuageuse). Je suppose que
1359  ! toliq=tomax quand l'epaisseur nuageuse est inferieure a dpmin,
1360  ! et que toliq=tomin quand l'epaisseur depasse dpmax (interpolation
1361  ! lineaire entre dpmin et dpmax).
1362  ! =====
1363  DO i = 1, klon
1364    IF (todo(i)) THEN
1365      toliq(i) = tomax - ((paprs(i,kb(i))-paprs(i,kh(i)+1))/paprs(i,1)-dpmin) &
1366        *(tomax-tomin)/(dpmax-dpmin)
1367      toliq(i) = max(tomin, min(tomax,toliq(i)))
1368      IF (pplay(i,kh(i))/paprs(i,1)<=deep_sig) toliq(i) = deep_to
1369      IF (old_tau) toliq(i) = 1.0
1370    END IF
1371  END DO
1372  ! =====
1373  ! On doit aussi determiner la distribution verticale de
1374  ! l'eau nuageuse. Plusieurs options sont proposees:
1375
1376  ! (0) La condensation precipite integralement (toliq ne sera
1377  ! pas utilise).
1378  ! (1) L'eau liquide est distribuee entre k1 et k2 et proportionnelle
1379  ! a la vapeur d'eau locale.
1380  ! (2) Elle est distribuee entre k1 et k2 avec une valeur constante.
1381  ! (3) Elle est seulement distribuee aux couches ou la vapeur d'eau
1382  ! est effectivement diminuee pendant le processus d'ajustement.
1383  ! (4) Elle est en fonction (lineaire ou exponentielle) de la
1384  ! distance (epaisseur en pression) avec le niveau k1 (la couche
1385  ! k1 n'aura donc pas d'eau liquide).
1386  ! =====
1387
1388  IF (opt_cld==0) THEN
1389
1390    DO i = 1, klon
1391      IF (todo(i)) zrfl(i) = zcond(i)/dtime
1392    END DO
1393
1394  ELSE IF (opt_cld==1) THEN
1395
1396    DO i = 1, klon
1397      IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de vapeur d'eau
1398    END DO
1399    DO k = kbmin, khmax
1400      DO i = 1, klon
1401        IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
1402          zvapo(i) = zvapo(i) + (q(i,k)+d_q(i,k))*(paprs(i,k)-paprs(i,k+1))/ &
1403            rg
1404        END IF
1405      END DO
1406    END DO
1407    DO i = 1, klon
1408      IF (todo(i)) THEN
1409        zrapp(i) = toliq(i)*zcond(i)/zvapo(i)
1410        zrapp(i) = max(0., min(1.,zrapp(i)))
1411      END IF
1412    END DO
1413    DO k = kbmin, khmax
1414      DO i = 1, klon
1415        IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
1416          d_ql(i, k) = zrapp(i)*(q(i,k)+d_q(i,k))
1417        END IF
1418      END DO
1419    END DO
1420    DO i = 1, klon
1421      IF (todo(i)) THEN
1422        zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
1423      END IF
1424    END DO
1425
1426  ELSE IF (opt_cld==2) THEN
1427
1428    DO i = 1, klon
1429      IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de masse
1430    END DO
1431    DO k = kbmin, khmax
1432      DO i = 1, klon
1433        IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
1434          zvapo(i) = zvapo(i) + (paprs(i,k)-paprs(i,k+1))/rg
1435        END IF
1436      END DO
1437    END DO
1438    DO k = kbmin, khmax
1439      DO i = 1, klon
1440        IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
1441          d_ql(i, k) = toliq(i)*zcond(i)/zvapo(i)
1442        END IF
1443      END DO
1444    END DO
1445    DO i = 1, klon
1446      IF (todo(i)) THEN
1447        zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
1448      END IF
1449    END DO
1450
1451  ELSE IF (opt_cld==3) THEN
1452
1453    DO i = 1, klon
1454      IF (todo(i)) THEN
1455        zvapo(i) = 0.0 ! quantite de l'eau strictement condensee
1456      END IF
1457    END DO
1458    DO k = kbmin, khmax
1459      DO i = 1, klon
1460        IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
1461          zvapo(i) = zvapo(i) + max(0.0, -d_q(i,k))*(paprs(i,k)-paprs(i,k+1)) &
1462            /rg
1463        END IF
1464      END DO
1465    END DO
1466    DO k = kbmin, khmax
1467      DO i = 1, klon
1468        IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i) .AND. zvapo(i)>0.0) THEN
1469          d_ql(i, k) = d_ql(i, k) + toliq(i)*zcond(i)/zvapo(i)*max(0.0, -d_q( &
1470            i,k))
1471        END IF
1472      END DO
1473    END DO
1474    DO i = 1, klon
1475      IF (todo(i)) THEN
1476        zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
1477      END IF
1478    END DO
1479
1480  ELSE IF (opt_cld==4) THEN
1481
1482    nexpo = 3
1483    ! cc         nexpo = 1 ! distribution lineaire
1484
1485    DO i = 1, klon
1486      IF (todo(i)) THEN
1487        zvapo(i) = 0.0 ! quantite integrale de masse (avec ponderation)
1488      END IF
1489    END DO
1490    DO k = kbmin, khmax
1491      DO i = 1, klon
1492        IF (todo(i) .AND. k>=(kb(i)+1) .AND. k<=kh(i)) THEN
1493          zvapo(i) = zvapo(i) + (paprs(i,k)-paprs(i,k+1))/rg*(pplay(i,kb(i))- &
1494            pplay(i,k))**nexpo
1495        END IF
1496      END DO
1497    END DO
1498    DO k = kbmin, khmax
1499      DO i = 1, klon
1500        IF (todo(i) .AND. k>=(kb(i)+1) .AND. k<=kh(i)) THEN
1501          d_ql(i, k) = d_ql(i, k) + toliq(i)*zcond(i)/zvapo(i)*(pplay(i,kb(i) &
1502            )-pplay(i,k))**nexpo
1503        END IF
1504      END DO
1505    END DO
1506    DO i = 1, klon
1507      IF (todo(i)) THEN
1508        zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
1509      END IF
1510    END DO
1511
1512  ELSE ! valeur non-prevue pour opt_cld
1513
1514    PRINT *, 'opt_cld est faux:', opt_cld
1515    CALL abort
1516
1517  END IF ! fin de opt_cld
1518
1519  ! L'eau precipitante peut etre re-evaporee:
1520
1521  IF (evap_prec .AND. kbmax>=2) THEN
1522    DO k = kbmax, 1, -1
1523      DO i = 1, klon
1524        IF (todo(i) .AND. k<=(kb(i)-1) .AND. zrfl(i)>0.0) THEN
1525          zqev = max(0.0, (zqs(i,k)-q(i,k))*zfrac(i))
1526          zqevt = coef_eva*(1.0-q(i,k)/zqs(i,k))*sqrt(zrfl(i))* &
1527            (paprs(i,k)-paprs(i,k+1))/pplay(i, k)*t(i, k)*rd/rg
1528          zqevt = max(0.0, min(zqevt,zrfl(i)))*rg*dtime/ &
1529            (paprs(i,k)-paprs(i,k+1))
1530          zqev = min(zqev, zqevt)
1531          zrfln = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))/rg/dtime
1532          d_q(i, k) = -(zrfln-zrfl(i))*(rg/(paprs(i,k)-paprs(i,k+1)))*dtime
1533          d_t(i, k) = (zrfln-zrfl(i))*(rg/(paprs(i,k)-paprs(i, &
1534            k+1)))*dtime*rlvtt/rcpd
1535          zrfl(i) = zrfln
1536        END IF
1537      END DO
1538    END DO
1539  END IF
1540
1541  ! La temperature de la premiere couche determine la pluie ou la neige:
1542
1543  DO i = 1, klon
1544    IF (todo(i)) THEN
1545      IF (t(i,1)>rtt) THEN
1546        rain(i) = rain(i) + zrfl(i)
1547      ELSE
1548        snow(i) = snow(i) + zrfl(i)
1549      END IF
1550    END IF
1551  END DO
1552
1553  RETURN
1554END SUBROUTINE conkuo
1555SUBROUTINE kuofcl(pt, pq, pg, pp, ldcum, kcbot)
1556  USE yomcst_mod_h
1557  USE dimphy
1558  USE yoethf_mod_h
1559  IMPLICIT NONE
1560  ! ======================================================================
1561  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19940927
1562  ! adaptation du code de Tiedtke du ECMWF
1563  ! Objet: calculer le niveau de convection libre
1564  ! (FCL: Free Convection Level)
1565  ! ======================================================================
1566  ! Arguments:
1567  ! pt---input-R- temperature (K)
1568  ! pq---input-R- vapeur d'eau (kg/kg)
1569  ! pg---input-R- geopotentiel (g*z ou z est en metre)
1570  ! pp---input-R- pression (Pa)
1571
1572  ! LDCUM---output-L- Y-t-il la convection
1573  ! kcbot---output-I- Niveau du bas de la convection
1574  ! ======================================================================
1575
1576  REAL pt(klon, klev), pq(klon, klev), pg(klon, klev), pp(klon, klev)
1577  INTEGER kcbot(klon)
1578  LOGICAL ldcum(klon)
1579
1580  REAL ztu(klon, klev), zqu(klon, klev), zlu(klon, klev)
1581  REAL zqold(klon), zbuo
1582  INTEGER is, i, k
1583
1584  ! klab=1: on est sous le nuage convectif
1585  ! klab=2: le bas du nuage convectif
1586  ! klab=0: autres couches
1587  INTEGER klab(klon, klev)
1588
1589  ! quand lflag=.true., on est sous le nuage, il faut donc appliquer
1590  ! le processus d'elevation.
1591  LOGICAL lflag(klon)
1592
1593  DO k = 1, klev
1594    DO i = 1, klon
1595      ztu(i, k) = pt(i, k)
1596      zqu(i, k) = pq(i, k)
1597      zlu(i, k) = 0.0
1598      klab(i, k) = 0
1599    END DO
1600  END DO
1601  ! ----------------------------------------------------------------------
1602  DO i = 1, klon
1603    klab(i, 1) = 1
1604    kcbot(i) = 2
1605    ldcum(i) = .FALSE.
1606  END DO
1607
1608  DO k = 2, klev - 1
1609
1610    is = 0
1611    DO i = 1, klon
1612      IF (klab(i,k-1)==1) is = is + 1
1613      lflag(i) = .FALSE.
1614      IF (klab(i,k-1)==1) lflag(i) = .TRUE.
1615    END DO
1616    IF (is==0) GO TO 290
1617
1618    ! on eleve le parcel d'air selon l'adiabatique sec
1619
1620    DO i = 1, klon
1621      IF (lflag(i)) THEN
1622        zqu(i, k) = zqu(i, k-1)
1623        ztu(i, k) = ztu(i, k-1) + (pg(i,k-1)-pg(i,k))/rcpd
1624        zbuo = ztu(i, k)*(1.+retv*zqu(i,k)) - pt(i, k)*(1.+retv*pq(i,k)) + &
1625          0.5
1626        IF (zbuo>0.) klab(i, k) = 1
1627        zqold(i) = zqu(i, k)
1628      END IF
1629    END DO
1630
1631    ! on calcule la condensation eventuelle
1632
1633    CALL adjtq(pp(1,k), ztu(1,k), zqu(1,k), lflag, 1)
1634
1635    ! s'il y a la condensation et la "buoyancy" force est positive
1636    ! c'est bien le bas de la tour de convection
1637
1638    DO i = 1, klon
1639      IF (lflag(i) .AND. zqu(i,k)/=zqold(i)) THEN
1640        klab(i, k) = 2
1641        zlu(i, k) = zlu(i, k) + zqold(i) - zqu(i, k)
1642        zbuo = ztu(i, k)*(1.+retv*zqu(i,k)) - pt(i, k)*(1.+retv*pq(i,k)) + &
1643          0.5
1644        IF (zbuo>0.) THEN
1645          kcbot(i) = k
1646          ldcum(i) = .TRUE.
1647        END IF
1648      END IF
1649    END DO
1650
1651290 END DO
1652
1653  RETURN
1654END SUBROUTINE kuofcl
1655SUBROUTINE adjtq(pp, pt, pq, ldflag, kcall)
1656  USE yomcst_mod_h
1657  USE dimphy
1658  USE yoethf_mod_h
1659  IMPLICIT NONE
1660  ! ======================================================================
1661  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19940927
1662  ! adaptation du code de Tiedtke du ECMWF
1663  ! Objet: ajustement entre T et Q
1664  ! ======================================================================
1665  ! Arguments:
1666  ! pp---input-R- pression (Pa)
1667  ! pt---input/output-R- temperature (K)
1668  ! pq---input/output-R- vapeur d'eau (kg/kg)
1669  ! ======================================================================
1670  ! TO PRODUCE T,Q AND L VALUES FOR CLOUD ASCENT
1671
1672  ! NOTE: INPUT PARAMETER KCALL DEFINES CALCULATION AS
1673  ! KCALL=0    ENV. T AND QS IN*CUINI*
1674  ! KCALL=1  CONDENSATION IN UPDRAFTS  (E.G. CUBASE, CUASC)
1675  ! KCALL=2  EVAPORATION IN DOWNDRAFTS (E.G. CUDLFS,CUDDRAF)
1676
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  include "FCTTRE.h"
1689
1690  DO i = 1, klon
1691    zcond(i) = 0.0
1692  END DO
1693
1694  DO i = 1, klon
1695    IF (ldflag(i)) THEN
1696      zdelta = max(0., sign(1.,rtt-pt(i)))
1697      zldcp = rlvtt*(1.-zdelta) + zdelta*rlstt
1698      zldcp = zldcp/rcpd/(1.0+rvtmp2*pq(i))
1699      IF (thermcep) THEN
1700        zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
1701        zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*pq(i))
1702        zqsat = r2es*foeew(pt(i), zdelta)/pp(i)
1703        zqsat = min(0.5, zqsat)
1704        zcor = 1./(1.-retv*zqsat)
1705        zqsat = zqsat*zcor
1706        zdqsat = foede(pt(i), zdelta, zcvm5, zqsat, zcor)
1707      ELSE
1708        IF (pt(i)<t_coup) THEN
1709          zqsat = qsats(pt(i))/pp(i)
1710          zdqsat = dqsats(pt(i), zqsat)
1711        ELSE
1712          zqsat = qsatl(pt(i))/pp(i)
1713          zdqsat = dqsatl(pt(i), zqsat)
1714        END IF
1715      END IF
1716      zcond(i) = (pq(i)-zqsat)/(1.+zdqsat)
1717      IF (kcall==1) zcond(i) = max(zcond(i), 0.)
1718      IF (kcall==2) zcond(i) = min(zcond(i), 0.)
1719      pt(i) = pt(i) + zldcp*zcond(i)
1720      pq(i) = pq(i) - zcond(i)
1721    END IF
1722  END DO
1723
1724  is = 0
1725  DO i = 1, klon
1726    IF (zcond(i)/=0.) is = is + 1
1727  END DO
1728  IF (is==0) GO TO 230
1729
1730  DO i = 1, klon
1731    IF (ldflag(i) .AND. zcond(i)/=0.) THEN
1732      zdelta = max(0., sign(1.,rtt-pt(i)))
1733      zldcp = rlvtt*(1.-zdelta) + zdelta*rlstt
1734      zldcp = zldcp/rcpd/(1.0+rvtmp2*pq(i))
1735      IF (thermcep) THEN
1736        zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
1737        zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*pq(i))
1738        zqsat = r2es*foeew(pt(i), zdelta)/pp(i)
1739        zqsat = min(0.5, zqsat)
1740        zcor = 1./(1.-retv*zqsat)
1741        zqsat = zqsat*zcor
1742        zdqsat = foede(pt(i), zdelta, zcvm5, zqsat, zcor)
1743      ELSE
1744        IF (pt(i)<t_coup) THEN
1745          zqsat = qsats(pt(i))/pp(i)
1746          zdqsat = dqsats(pt(i), zqsat)
1747        ELSE
1748          zqsat = qsatl(pt(i))/pp(i)
1749          zdqsat = dqsatl(pt(i), zqsat)
1750        END IF
1751      END IF
1752      zcond1 = (pq(i)-zqsat)/(1.+zdqsat)
1753      pt(i) = pt(i) + zldcp*zcond1
1754      pq(i) = pq(i) - zcond1
1755    END IF
1756  END DO
1757
1758230 CONTINUE
1759  RETURN
1760END SUBROUTINE adjtq
1761SUBROUTINE fiajh(dtime, paprs, pplay, t, q, d_t, d_q, d_ql, rneb, rain, snow, &
1762    ibas, itop)
1763  USE dimphy
1764  USE yomcst_mod_h
1765  USE yoethf_mod_h
1766IMPLICIT NONE
1767
1768  ! Ajustement humide (Schema de convection de Manabe)
1769  ! .
1770
1771
1772  ! Arguments:
1773
1774  REAL dtime ! intervalle du temps (s)
1775  REAL t(klon, klev) ! temperature (K)
1776  REAL q(klon, klev) ! humidite specifique (kg/kg)
1777  REAL paprs(klon, klev+1) ! pression a inter-couche (Pa)
1778  REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
1779
1780  REAL d_t(klon, klev) ! incrementation pour la temperature
1781  REAL d_q(klon, klev) ! incrementation pour vapeur d'eau
1782  REAL d_ql(klon, klev) ! incrementation pour l'eau liquide
1783  REAL rneb(klon, klev) ! fraction nuageuse
1784
1785  REAL rain(klon) ! variable non utilisee
1786  REAL snow(klon) ! variable non utilisee
1787  INTEGER ibas(klon) ! variable non utilisee
1788  INTEGER itop(klon) ! variable non utilisee
1789
1790  REAL t_coup
1791  PARAMETER (t_coup=234.0)
1792  REAL seuil_vap
1793  PARAMETER (seuil_vap=1.0E-10)
1794
1795  ! Variables locales:
1796
1797  INTEGER i, k
1798  INTEGER k1, k1p, k2, k2p
1799  LOGICAL itest(klon)
1800  REAL delta_q(klon, klev)
1801  REAL cp_new_t(klev)
1802  REAL cp_delta_t(klev)
1803  REAL new_qb(klev)
1804  REAL v_cptj(klev), v_cptjk1, v_ssig
1805  REAL v_cptt(klon, klev), v_p, v_t
1806  REAL v_qs(klon, klev), v_qsd(klon, klev)
1807  REAL zq1(klon), zq2(klon)
1808  REAL gamcpdz(klon, 2:klev)
1809  REAL zdp, zdpm
1810
1811  REAL zsat ! sur-saturation
1812  REAL zflo ! flotabilite
1813
1814  REAL local_q(klon, klev), local_t(klon, klev)
1815
1816  REAL zdelta, zcor, zcvm5
1817
1818  include "FCTTRE.h"
1819
1820  DO k = 1, klev
1821    DO i = 1, klon
1822      local_q(i, k) = q(i, k)
1823      local_t(i, k) = t(i, k)
1824      rneb(i, k) = 0.0
1825      d_ql(i, k) = 0.0
1826      d_t(i, k) = 0.0
1827      d_q(i, k) = 0.0
1828    END DO
1829  END DO
1830  DO i = 1, klon
1831    rain(i) = 0.0
1832    snow(i) = 0.0
1833    ibas(i) = 0
1834    itop(i) = 0
1835  END DO
1836
1837  ! Calculer v_qs et v_qsd:
1838
1839  DO k = 1, klev
1840    DO i = 1, klon
1841      v_cptt(i, k) = rcpd*local_t(i, k)
1842      v_t = local_t(i, k)
1843      v_p = pplay(i, k)
1844
1845      IF (thermcep) THEN
1846        zdelta = max(0., sign(1.,rtt-v_t))
1847        zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
1848        zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*local_q(i,k))
1849        v_qs(i, k) = r2es*foeew(v_t, zdelta)/v_p
1850        v_qs(i, k) = min(0.5, v_qs(i,k))
1851        zcor = 1./(1.-retv*v_qs(i,k))
1852        v_qs(i, k) = v_qs(i, k)*zcor
1853        v_qsd(i, k) = foede(v_t, zdelta, zcvm5, v_qs(i,k), zcor)
1854      ELSE
1855        IF (v_t<t_coup) THEN
1856          v_qs(i, k) = qsats(v_t)/v_p
1857          v_qsd(i, k) = dqsats(v_t, v_qs(i,k))
1858        ELSE
1859          v_qs(i, k) = qsatl(v_t)/v_p
1860          v_qsd(i, k) = dqsatl(v_t, v_qs(i,k))
1861        END IF
1862      END IF
1863    END DO
1864  END DO
1865
1866  ! Calculer Gamma * Cp * dz: (gamm est le gradient critique)
1867
1868  DO k = 2, klev
1869    DO i = 1, klon
1870      zdp = paprs(i, k) - paprs(i, k+1)
1871      zdpm = paprs(i, k-1) - paprs(i, k)
1872      gamcpdz(i, k) = ((rd/rcpd/(zdpm+zdp)*(v_cptt(i,k-1)*zdpm+ &
1873        v_cptt(i,k)*zdp)+rlvtt/(zdpm+zdp)*(v_qs(i,k-1)*zdpm+ &
1874        v_qs(i,k)*zdp))*(pplay(i,k-1)-pplay(i,k))/paprs(i,k))/(1.0+(v_qsd(i, &
1875        k-1)*zdpm+v_qsd(i,k)*zdp)/(zdpm+zdp))
1876    END DO
1877  END DO
1878
1879  ! ------------------------------------ modification des profils instables
1880  DO i = 1, klon
1881    itest(i) = .FALSE.
1882
1883    k1 = 0
1884    k2 = 1
1885
1886810 CONTINUE ! chercher k1, le bas de la colonne
1887    k2 = k2 + 1
1888    IF (k2>klev) GO TO 9999
1889    zflo = v_cptt(i, k2-1) - v_cptt(i, k2) - gamcpdz(i, k2)
1890    zsat = (local_q(i,k2-1)-v_qs(i,k2-1))*(paprs(i,k2-1)-paprs(i,k2)) + &
1891      (local_q(i,k2)-v_qs(i,k2))*(paprs(i,k2)-paprs(i,k2+1))
1892    IF (zflo<=0.0 .OR. zsat<=0.0) GO TO 810
1893    k1 = k2 - 1
1894    itest(i) = .TRUE.
1895
1896820 CONTINUE ! chercher k2, le haut de la colonne
1897    IF (k2==klev) GO TO 821
1898    k2p = k2 + 1
1899    zsat = zsat + (paprs(i,k2p)-paprs(i,k2p+1))*(local_q(i,k2p)-v_qs(i,k2p))
1900    zflo = v_cptt(i, k2p-1) - v_cptt(i, k2p) - gamcpdz(i, k2p)
1901    IF (zflo<=0.0 .OR. zsat<=0.0) GO TO 821
1902    k2 = k2p
1903    GO TO 820
1904821 CONTINUE
1905
1906    ! ------------------------------------------------------ ajustement local
1907830 CONTINUE ! ajustement proprement dit
1908    v_cptj(k1) = 0.0
1909    zdp = paprs(i, k1) - paprs(i, k1+1)
1910    v_cptjk1 = ((1.0+v_qsd(i,k1))*(v_cptt(i,k1)+v_cptj(k1))+rlvtt*(local_q(i, &
1911      k1)-v_qs(i,k1)))*zdp
1912    v_ssig = zdp*(1.0+v_qsd(i,k1))
1913
1914    k1p = k1 + 1
1915    DO k = k1p, k2
1916      zdp = paprs(i, k) - paprs(i, k+1)
1917      v_cptj(k) = v_cptj(k-1) + gamcpdz(i, k)
1918      v_cptjk1 = v_cptjk1 + zdp*((1.0+v_qsd(i,k))*(v_cptt(i, &
1919        k)+v_cptj(k))+rlvtt*(local_q(i,k)-v_qs(i,k)))
1920      v_ssig = v_ssig + zdp*(1.0+v_qsd(i,k))
1921    END DO
1922
1923    DO k = k1, k2
1924      cp_new_t(k) = v_cptjk1/v_ssig - v_cptj(k)
1925      cp_delta_t(k) = cp_new_t(k) - v_cptt(i, k)
1926      new_qb(k) = v_qs(i, k) + v_qsd(i, k)*cp_delta_t(k)/rlvtt
1927      local_q(i, k) = new_qb(k)
1928      local_t(i, k) = cp_new_t(k)/rcpd
1929    END DO
1930
1931    ! --------------------------------------------------- sondage vers le bas
1932    ! -- on redefinit les variables prognostiques dans
1933    ! -- la colonne qui vient d'etre ajustee
1934
1935    DO k = k1, k2
1936      v_cptt(i, k) = rcpd*local_t(i, k)
1937      v_t = local_t(i, k)
1938      v_p = pplay(i, k)
1939
1940      IF (thermcep) THEN
1941        zdelta = max(0., sign(1.,rtt-v_t))
1942        zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
1943        zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*local_q(i,k))
1944        v_qs(i, k) = r2es*foeew(v_t, zdelta)/v_p
1945        v_qs(i, k) = min(0.5, v_qs(i,k))
1946        zcor = 1./(1.-retv*v_qs(i,k))
1947        v_qs(i, k) = v_qs(i, k)*zcor
1948        v_qsd(i, k) = foede(v_t, zdelta, zcvm5, v_qs(i,k), zcor)
1949      ELSE
1950        IF (v_t<t_coup) THEN
1951          v_qs(i, k) = qsats(v_t)/v_p
1952          v_qsd(i, k) = dqsats(v_t, v_qs(i,k))
1953        ELSE
1954          v_qs(i, k) = qsatl(v_t)/v_p
1955          v_qsd(i, k) = dqsatl(v_t, v_qs(i,k))
1956        END IF
1957      END IF
1958    END DO
1959    DO k = 2, klev
1960      zdpm = paprs(i, k-1) - paprs(i, k)
1961      zdp = paprs(i, k) - paprs(i, k+1)
1962      gamcpdz(i, k) = ((rd/rcpd/(zdpm+zdp)*(v_cptt(i,k-1)*zdpm+ &
1963        v_cptt(i,k)*zdp)+rlvtt/(zdpm+zdp)*(v_qs(i,k-1)*zdpm+ &
1964        v_qs(i,k)*zdp))*(pplay(i,k-1)-pplay(i,k))/paprs(i,k))/(1.0+(v_qsd(i, &
1965        k-1)*zdpm+v_qsd(i,k)*zdp)/(zdpm+zdp))
1966    END DO
1967
1968    ! Verifier si l'on peut etendre la colonne vers le bas
1969
1970    IF (k1==1) GO TO 841 ! extension echouee
1971    zflo = v_cptt(i, k1-1) - v_cptt(i, k1) - gamcpdz(i, k1)
1972    zsat = (local_q(i,k1-1)-v_qs(i,k1-1))*(paprs(i,k1-1)-paprs(i,k1)) + &
1973      (local_q(i,k1)-v_qs(i,k1))*(paprs(i,k1)-paprs(i,k1+1))
1974    IF (zflo<=0.0 .OR. zsat<=0.0) GO TO 841 ! extension echouee
1975
1976840 CONTINUE
1977    k1 = k1 - 1
1978    IF (k1==1) GO TO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
1979    zsat = zsat + (local_q(i,k1-1)-v_qs(i,k1-1))*(paprs(i,k1-1)-paprs(i,k1))
1980    zflo = v_cptt(i, k1-1) - v_cptt(i, k1) - gamcpdz(i, k1)
1981    IF (zflo>0.0 .AND. zsat>0.0) THEN
1982      GO TO 840
1983    ELSE
1984      GO TO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
1985    END IF
1986841 CONTINUE
1987
1988    GO TO 810 ! chercher d'autres blocks en haut
1989
19909999 END DO ! boucle sur tous les points
1991  ! -----------------------------------------------------------------------
1992
1993  ! Determiner la fraction nuageuse (hypothese: la nebulosite a lieu
1994  ! a l'endroit ou la vapeur d'eau est diminuee par l'ajustement):
1995
1996  DO k = 1, klev
1997    DO i = 1, klon
1998      IF (itest(i)) THEN
1999        delta_q(i, k) = local_q(i, k) - q(i, k)
2000        IF (delta_q(i,k)<0.) rneb(i, k) = 1.0
2001      END IF
2002    END DO
2003  END DO
2004
2005  ! Distribuer l'eau condensee en eau liquide nuageuse (hypothese:
2006  ! l'eau liquide est distribuee aux endroits ou la vapeur d'eau
2007  ! diminue et d'une maniere proportionnelle a cet diminution):
2008
2009  DO i = 1, klon
2010    IF (itest(i)) THEN
2011      zq1(i) = 0.0
2012      zq2(i) = 0.0
2013    END IF
2014  END DO
2015  DO k = 1, klev
2016    DO i = 1, klon
2017      IF (itest(i)) THEN
2018        zdp = paprs(i, k) - paprs(i, k+1)
2019        zq1(i) = zq1(i) - delta_q(i, k)*zdp
2020        zq2(i) = zq2(i) - min(0.0, delta_q(i,k))*zdp
2021      END IF
2022    END DO
2023  END DO
2024  DO k = 1, klev
2025    DO i = 1, klon
2026      IF (itest(i)) THEN
2027        IF (zq2(i)/=0.0) d_ql(i, k) = -min(0.0, delta_q(i,k))*zq1(i)/zq2(i)
2028      END IF
2029    END DO
2030  END DO
2031
2032  DO k = 1, klev
2033    DO i = 1, klon
2034      local_q(i, k) = max(local_q(i,k), seuil_vap)
2035    END DO
2036  END DO
2037
2038  DO k = 1, klev
2039    DO i = 1, klon
2040      d_t(i, k) = local_t(i, k) - t(i, k)
2041      d_q(i, k) = local_q(i, k) - q(i, k)
2042    END DO
2043  END DO
2044
2045  RETURN
2046END SUBROUTINE fiajh
2047SUBROUTINE fiajc(dtime, paprs, pplay, t, q, conv_q, d_t, d_q, d_ql, rneb, &
2048    rain, snow, ibas, itop)
2049  USE dimphy
2050  USE yomcst_mod_h
2051  USE yoethf_mod_h
2052IMPLICIT NONE
2053
2054
2055
2056  ! Options:
2057
2058  INTEGER plb ! niveau de depart pour la convection
2059  PARAMETER (plb=4)
2060
2061  ! Mystere: cette option n'est pas innocente pour les resultats !
2062  ! Qui peut resoudre ce mystere ? (Z.X.Li mars 1995)
2063  LOGICAL vector ! calcul vectorise
2064  PARAMETER (vector=.FALSE.)
2065
2066  REAL t_coup
2067  PARAMETER (t_coup=234.0)
2068
2069  ! Arguments:
2070
2071  REAL q(klon, klev) ! humidite specifique (kg/kg)
2072  REAL t(klon, klev) ! temperature (K)
2073  REAL paprs(klon, klev+1) ! pression a inter-couche (Pa)
2074  REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
2075  REAL dtime ! intervalle du temps (s)
2076  REAL conv_q(klon, klev) ! taux de convergence de l'humidite
2077  REAL rneb(klon, klev) ! fraction nuageuse
2078  REAL d_q(klon, klev) ! incrementaion pour la vapeur d'eau
2079  REAL d_ql(klon, klev) ! incrementation pour l'eau liquide
2080  REAL d_t(klon, klev) ! incrementation pour la temperature
2081  REAL rain(klon) ! variable non-utilisee
2082  REAL snow(klon) ! variable non-utilisee
2083  INTEGER itop(klon) ! variable non-utilisee
2084  INTEGER ibas(klon) ! variable non-utilisee
2085
2086  INTEGER kh(klon), i, k
2087  LOGICAL nuage(klon), test(klon, klev)
2088  REAL zconv(klon), zdeh(klon, klev), zvirt(klon)
2089  REAL zdqs(klon, klev), zqs(klon, klev)
2090  REAL ztt, zvar, zfrac(klon)
2091  REAL zq1(klon), zq2(klon)
2092  REAL zdelta, zcor, zcvm5
2093
2094  include "FCTTRE.h"
2095
2096  ! Initialiser les sorties:
2097
2098  DO k = 1, klev
2099    DO i = 1, klon
2100      rneb(i, k) = 0.0
2101      d_ql(i, k) = 0.0
2102      d_t(i, k) = 0.0
2103      d_q(i, k) = 0.0
2104    END DO
2105  END DO
2106  DO i = 1, klon
2107    itop(i) = 0
2108    ibas(i) = 0
2109    rain(i) = 0.0
2110    snow(i) = 0.0
2111  END DO
2112
2113  ! Calculer Qs et L/Cp * dQs/dT:
2114
2115  DO k = 1, klev
2116    DO i = 1, klon
2117      ztt = t(i, k)
2118      IF (thermcep) THEN
2119        zdelta = max(0., sign(1.,rtt-ztt))
2120        zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
2121        zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*q(i,k))
2122        zqs(i, k) = r2es*foeew(ztt, zdelta)/pplay(i, k)
2123        zqs(i, k) = min(0.5, zqs(i,k))
2124        zcor = 1./(1.-retv*zqs(i,k))
2125        zqs(i, k) = zqs(i, k)*zcor
2126        zdqs(i, k) = foede(ztt, zdelta, zcvm5, zqs(i,k), zcor)
2127      ELSE
2128        IF (ztt<t_coup) THEN
2129          zqs(i, k) = qsats(ztt)/pplay(i, k)
2130          zdqs(i, k) = dqsats(ztt, zqs(i,k))
2131        ELSE
2132          zqs(i, k) = qsatl(ztt)/pplay(i, k)
2133          zdqs(i, k) = dqsatl(ztt, zqs(i,k))
2134        END IF
2135      END IF
2136    END DO
2137  END DO
2138
2139  ! Determiner la difference de l'energie totale saturee:
2140
2141  DO i = 1, klon
2142    k = plb
2143    zdeh(i, k) = rcpd*(t(i,k-1)-t(i,k)) - rd*0.5*(t(i,k-1)+t(i,k))/paprs(i, k &
2144      )*(pplay(i,k-1)-pplay(i,k)) + rlvtt*(zqs(i,k-1)-zqs(i,k))
2145    zdeh(i, k) = zdeh(i, k)*0.5 ! on prend la moitie
2146  END DO
2147  DO k = plb + 1, klev
2148    DO i = 1, klon
2149      zdeh(i, k) = zdeh(i, k-1) + rcpd*(t(i,k-1)-t(i,k)) - &
2150        rd*0.5*(t(i,k-1)+t(i,k))/paprs(i, k)*(pplay(i,k-1)-pplay(i,k)) + &
2151        rlvtt*(zqs(i,k-1)-zqs(i,k))
2152    END DO
2153  END DO
2154
2155  ! Determiner le sommet du nuage selon l'instabilite
2156  ! Calculer les convergences d'humidite (reelle et virtuelle)
2157
2158  DO i = 1, klon
2159    nuage(i) = .TRUE.
2160    zconv(i) = 0.0
2161    zvirt(i) = 0.0
2162    kh(i) = -999
2163  END DO
2164  DO k = plb, klev
2165    DO i = 1, klon
2166      nuage(i) = nuage(i) .AND. zdeh(i, k) > 0.0
2167      IF (nuage(i)) THEN
2168        kh(i) = k
2169        zconv(i) = zconv(i) + conv_q(i, k)*dtime*(paprs(i,k)-paprs(i,k+1))
2170        zvirt(i) = zvirt(i) + (zdeh(i,k)/rlvtt+zqs(i,k)-q(i,k))*(paprs(i,k)- &
2171          paprs(i,k+1))
2172      END IF
2173    END DO
2174  END DO
2175
2176  IF (vector) THEN
2177
2178
2179    DO k = plb, klev
2180      DO i = 1, klon
2181        IF (k<=kh(i) .AND. kh(i)>plb .AND. zconv(i)>0.0) THEN
2182          test(i, k) = .TRUE.
2183          zfrac(i) = max(0.0, min(zconv(i)/zvirt(i),1.0))
2184        ELSE
2185          test(i, k) = .FALSE.
2186        END IF
2187      END DO
2188    END DO
2189
2190    DO k = plb, klev
2191      DO i = 1, klon
2192        IF (test(i,k)) THEN
2193          zvar = zdeh(i, k)/(1.0+zdqs(i,k))
2194          d_q(i, k) = (zvar*zdqs(i,k)/rlvtt+zqs(i,k)-q(i,k))*zfrac(i) - &
2195            conv_q(i, k)*dtime
2196          d_t(i, k) = zvar*zfrac(i)/rcpd
2197        END IF
2198      END DO
2199    END DO
2200
2201    DO i = 1, klon
2202      zq1(i) = 0.0
2203      zq2(i) = 0.0
2204    END DO
2205    DO k = plb, klev
2206      DO i = 1, klon
2207        IF (test(i,k)) THEN
2208          IF (d_q(i,k)<0.0) rneb(i, k) = zfrac(i)
2209          zq1(i) = zq1(i) - d_q(i, k)*(paprs(i,k)-paprs(i,k+1))
2210          zq2(i) = zq2(i) - min(0.0, d_q(i,k))*(paprs(i,k)-paprs(i,k+1))
2211        END IF
2212      END DO
2213    END DO
2214
2215    DO k = plb, klev
2216      DO i = 1, klon
2217        IF (test(i,k)) THEN
2218          IF (zq2(i)/=0.) d_ql(i, k) = -min(0.0, d_q(i,k))*zq1(i)/zq2(i)
2219        END IF
2220      END DO
2221    END DO
2222
2223  ELSE ! (.NOT. vector)
2224
2225    DO i = 1, klon
2226      IF (kh(i)>plb .AND. zconv(i)>0.0) THEN
2227        ! cc         IF (kh(i).LE.plb) GOTO 999 ! il n'y a pas d'instabilite
2228        ! cc         IF (zconv(i).LE.0.0) GOTO 999 ! convergence insuffisante
2229        zfrac(i) = max(0.0, min(zconv(i)/zvirt(i),1.0))
2230        DO k = plb, kh(i)
2231          zvar = zdeh(i, k)/(1.0+zdqs(i,k))
2232          d_q(i, k) = (zvar*zdqs(i,k)/rlvtt+zqs(i,k)-q(i,k))*zfrac(i) - &
2233            conv_q(i, k)*dtime
2234          d_t(i, k) = zvar*zfrac(i)/rcpd
2235        END DO
2236
2237        zq1(i) = 0.0
2238        zq2(i) = 0.0
2239        DO k = plb, kh(i)
2240          IF (d_q(i,k)<0.0) rneb(i, k) = zfrac(i)
2241          zq1(i) = zq1(i) - d_q(i, k)*(paprs(i,k)-paprs(i,k+1))
2242          zq2(i) = zq2(i) - min(0.0, d_q(i,k))*(paprs(i,k)-paprs(i,k+1))
2243        END DO
2244        DO k = plb, kh(i)
2245          IF (zq2(i)/=0.) d_ql(i, k) = -min(0.0, d_q(i,k))*zq1(i)/zq2(i)
2246        END DO
2247      END IF
2248    END DO
2249
2250  END IF ! fin de teste sur vector
2251
2252  RETURN
2253END SUBROUTINE fiajc
Note: See TracBrowser for help on using the repository browser.