source: LMDZ5/branches/testing/libf/phylmd/conlmd.F90 @ 2056

Last change on this file since 2056 was 1999, checked in by Laurent Fairhead, 11 years ago

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