source: LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/conlmd.F90 @ 3713

Last change on this file since 3713 was 2346, checked in by Ehouarn Millour, 9 years ago

Physics/dynamics separation:

  • remove all references to dimensions.h from physics. nbp_lon (==iim) , nbp_lat (==jjm+1) and nbp_lev (==llm) from mod_grid_phy_lmdz should be used instead.
  • added module regular_lonlat_mod in phy_common to store information about the global (lon-lat) grid cell boundaries and centers.

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