source: LMDZ5/trunk/libf/phylmd/conlmd.F @ 1907

Last change on this file since 1907 was 1907, checked in by lguez, 10 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

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