source: LMDZ6/trunk/libf/phylmd/conlmd.f90

Last change on this file was 5274, checked in by abarral, 51 minutes ago

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