source: LMDZ.3.3/branches/rel-1-0-patch/libf/phylmd/conlmd.F @ 2167

Last change on this file since 2167 was 42, checked in by lmdz, 25 years ago

Ajustement de parametres et introduction de modifs de S. Cheinet. L.Li
Lf

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