source: LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_plume_6A.f90 @ 5157

Last change on this file since 5157 was 5119, checked in by abarral, 4 months ago

enforce PRIVATE by default in several modules, expose PUBLIC as needed
move eigen.f90 to obsolete/
(lint) aslong the way
  • 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: 43.1 KB
Line 
1! $Id: lmdz_thermcell_plume_6A.f90 5119 2024-07-24 16:46:45Z abarral $
2
3
4MODULE lmdz_thermcell_plume_6A
5
6CONTAINS
7
8  SUBROUTINE thermcell_plume_6A(itap, ngrid, nlay, ptimestep, ztv, zthl, po, zl, rhobarz, &
9          zlev, pplev, pphi, zpspsk, alim_star, alim_star_tot, &
10          lalim, f0, detr_star, entr_star, f_star, csc, ztva, &
11          ztla, zqla, zqta, zha, zw2, w_est, ztva_est, zqsatth, lmix, lmix_bis, linter &
12          , lev_out, lunout1, igout)
13    !     &           ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
14    !--------------------------------------------------------------------------
15    !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
16    !--------------------------------------------------------------------------
17
18    USE lmdz_thermcell_ini, ONLY: prt_level, fact_thermals_ed_dz, iflag_thermals_ed, RLvCP, RETV, RG
19    USE lmdz_thermcell_ini, ONLY: fact_epsilon, betalpha, afact, fact_shell
20    USE lmdz_thermcell_ini, ONLY: detr_min, entr_min, detr_q_coef, detr_q_power
21    USE lmdz_thermcell_ini, ONLY: mix0, thermals_flag_alim
22    USE lmdz_thermcell_alim, ONLY: thermcell_alim
23    USE lmdz_thermcell_qsat, ONLY: thermcell_qsat
24
25    IMPLICIT NONE
26
27    INTEGER, INTENT(IN) :: itap, lev_out, lunout1, igout, ngrid, nlay
28    REAL, INTENT(IN) :: ptimestep
29    REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: ztv
30    REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: zthl
31    REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: po
32    REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: zl
33    REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: rhobarz
34    REAL, INTENT(IN), DIMENSION(ngrid, nlay + 1) :: zlev
35    REAL, INTENT(IN), DIMENSION(ngrid, nlay + 1) :: pplev
36    REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: pphi
37    REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: zpspsk
38    REAL, INTENT(IN), DIMENSION(ngrid) :: f0
39
40    INTEGER, INTENT(OUT) :: lalim(ngrid)
41    REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: alim_star
42    REAL, INTENT(OUT), DIMENSION(ngrid) :: alim_star_tot
43    REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: detr_star
44    REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: entr_star
45    REAL, INTENT(OUT), DIMENSION(ngrid, nlay + 1) :: f_star
46    REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: csc
47    REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: ztva
48    REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: ztla
49    REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: zqla
50    REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: zqta
51    REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: zha
52    REAL, INTENT(OUT), DIMENSION(ngrid, nlay + 1) :: zw2
53    REAL, INTENT(OUT), DIMENSION(ngrid, nlay + 1) :: w_est
54    REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: ztva_est
55    REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: zqsatth
56    INTEGER, INTENT(OUT), DIMENSION(ngrid) :: lmix
57    INTEGER, INTENT(OUT), DIMENSION(ngrid) :: lmix_bis
58    REAL, INTENT(OUT), DIMENSION(ngrid) :: linter
59
60    REAL zdw2, zdw2bis
61    REAL zw2modif
62    REAL zw2fact, zw2factbis
63    REAL, DIMENSION(ngrid, nlay) :: zeps
64
65    REAL, DIMENSION(ngrid) :: wmaxa
66
67    INTEGER ig, l, k, lt, it, lm
68    INTEGER nbpb
69
70    REAL, DIMENSION(ngrid, nlay) :: detr
71    REAL, DIMENSION(ngrid, nlay) :: entr
72    REAL, DIMENSION(ngrid, nlay + 1) :: wa_moy
73    REAL, DIMENSION(ngrid, nlay) :: ztv_est
74    REAL, DIMENSION(ngrid) :: ztemp, zqsat
75    REAL, DIMENSION(ngrid, nlay) :: zqla_est
76    REAL, DIMENSION(ngrid, nlay) :: zta_est
77
78    REAL, DIMENSION(ngrid, nlay) :: zbuoy, gamma, zdqt
79    REAL zdz, zalpha, zw2m
80    REAL, DIMENSION(ngrid, nlay) :: zbuoyjam, zdqtjam
81    REAL zbuoybis, zdz2, zdz3, lmel, entrbis, zdzbis
82    REAL, DIMENSION(ngrid) :: d_temp
83    REAL ztv1, ztv2, factinv, zinv, zlmel
84    REAL zlmelup, zlmeldwn, zlt, zltdwn, zltup
85    REAL atv1, atv2, btv1, btv2
86    REAL ztv_est1, ztv_est2
87    REAL zcor, zdelta, zcvm5, qlbef
88    REAL zbetalpha, coefzlmel
89    REAL eps
90    LOGICAL Zsat
91    LOGICAL, DIMENSION(ngrid) :: active, activetmp
92    REAL fact_gamma, fact_gamma2, fact_epsilon2
93    REAL coefc
94    REAL, DIMENSION(ngrid, nlay) :: c2
95
96    IF (ngrid==1) PRINT*, 'THERMCELL PLUME MODIFIE 2014/07/11'
97    Zsat = .FALSE.
98    ! Initialisation
99
100    zbetalpha = betalpha / (1. + betalpha)
101
102
103    ! Initialisations des variables r?elles
104    IF (1==1) THEN
105      ztva(:, :) = ztv(:, :)
106      ztva_est(:, :) = ztva(:, :)
107      ztv_est(:, :) = ztv(:, :)
108      ztla(:, :) = zthl(:, :)
109      zqta(:, :) = po(:, :)
110      zqla(:, :) = 0.
111      zha(:, :) = ztva(:, :)
112    else
113      ztva(:, :) = 0.
114      ztv_est(:, :) = 0.
115      ztva_est(:, :) = 0.
116      ztla(:, :) = 0.
117      zqta(:, :) = 0.
118      zha(:, :) = 0.
119    endif
120
121    zqla_est(:, :) = 0.
122    zqsatth(:, :) = 0.
123    zqla(:, :) = 0.
124    detr_star(:, :) = 0.
125    entr_star(:, :) = 0.
126    alim_star(:, :) = 0.
127    alim_star_tot(:) = 0.
128    csc(:, :) = 0.
129    detr(:, :) = 0.
130    entr(:, :) = 0.
131    zw2(:, :) = 0.
132    zbuoy(:, :) = 0.
133    zbuoyjam(:, :) = 0.
134    gamma(:, :) = 0.
135    zeps(:, :) = 0.
136    w_est(:, :) = 0.
137    f_star(:, :) = 0.
138    wa_moy(:, :) = 0.
139    linter(:) = 1.
140    !     linter(:)=1.
141    ! Initialisation des variables entieres
142    lmix(:) = 1
143    lmix_bis(:) = 2
144    wmaxa(:) = 0.
145
146    ! Initialisation a 0  en cas de sortie dans replay
147    zqsat(:) = 0.
148    zta_est(:, :) = 0.
149    zdqt(:, :) = 0.
150    zdqtjam(:, :) = 0.
151    c2(:, :) = 0.
152
153
154    !-------------------------------------------------------------------------
155    ! On ne considere comme actif que les colonnes dont les deux premieres
156    ! couches sont instables.
157    !-------------------------------------------------------------------------
158
159    active(:) = ztv(:, 1)>ztv(:, 2)
160    d_temp(:) = 0. ! Pour activer un contraste de temperature a la base
161    ! du panache
162    !  Cet appel pourrait être fait avant thermcell_plume dans thermcell_main
163    CALL thermcell_alim(thermals_flag_alim, ngrid, nlay, ztv, d_temp, zlev, alim_star, lalim)
164
165    !------------------------------------------------------------------------------
166    ! Calcul dans la premiere couche
167    ! On decide dans cette version que le thermique n'est actif que si la premiere
168    ! couche est instable.
169    ! Pourrait etre change si on veut que le thermiques puisse se d??clencher
170    ! dans une couche l>1
171    !------------------------------------------------------------------------------
172    do ig = 1, ngrid
173      ! Le panache va prendre au debut les caracteristiques de l'air contenu
174      ! dans cette couche.
175      IF (active(ig)) THEN
176        ztla(ig, 1) = zthl(ig, 1)
177        zqta(ig, 1) = po(ig, 1)
178        zqla(ig, 1) = zl(ig, 1)
179        !cr: attention, prise en compte de f*(1)=1
180        f_star(ig, 2) = alim_star(ig, 1)
181        zw2(ig, 2) = 2. * RG * (ztv(ig, 1) - ztv(ig, 2)) / ztv(ig, 2)  &
182                & * (zlev(ig, 2) - zlev(ig, 1))  &
183                & * 0.4 * pphi(ig, 1) / (pphi(ig, 2) - pphi(ig, 1))
184        w_est(ig, 2) = zw2(ig, 2)
185      endif
186    enddo
187
188    !==============================================================================
189    !boucle de calcul de la vitesse verticale dans le thermique
190    !==============================================================================
191    do l = 2, nlay - 1
192      !==============================================================================
193
194
195      ! On decide si le thermique est encore actif ou non
196      ! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
197      do ig = 1, ngrid
198        active(ig) = active(ig) &
199                &                 .AND. zw2(ig, l)>1.e-10 &
200                &                 .AND. f_star(ig, l) + alim_star(ig, l)>1.e-10
201      enddo
202
203
204
205      !---------------------------------------------------------------------------
206      ! calcul des proprietes thermodynamiques et de la vitesse de la couche l
207      ! sans tenir compte du detrainement et de l'entrainement dans cette
208      ! couche
209      ! C'est a dire qu'on suppose
210      ! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
211      ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
212      ! avant) a l'alimentation pour avoir un calcul plus propre
213      !---------------------------------------------------------------------------
214
215      ztemp(:) = zpspsk(:, l) * ztla(:, l - 1)
216      CALL thermcell_qsat(ngrid, active, pplev(:, l), ztemp, zqta(:, l - 1), zqsat(:))
217      do ig = 1, ngrid
218        !       PRINT*,'active',active(ig),ig,l
219        IF(active(ig)) THEN
220          zqla_est(ig, l) = max(0., zqta(ig, l - 1) - zqsat(ig))
221          ztva_est(ig, l) = ztla(ig, l - 1) * zpspsk(ig, l) + RLvCp * zqla_est(ig, l)
222          zta_est(ig, l) = ztva_est(ig, l)
223          ztva_est(ig, l) = ztva_est(ig, l) / zpspsk(ig, l)
224          ztva_est(ig, l) = ztva_est(ig, l) * (1. + RETV * (zqta(ig, l - 1)  &
225                  - zqla_est(ig, l)) - zqla_est(ig, l))
226
227
228          !Modif AJAM
229
230          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
231          zdz = zlev(ig, l + 1) - zlev(ig, l)
232          lmel = fact_thermals_ed_dz * zlev(ig, l)
233          !        lmel=0.09*zlev(ig,l)
234          zlmel = zlev(ig, l) + lmel
235          zlmelup = zlmel + (zdz / 2)
236          zlmeldwn = zlmel - (zdz / 2)
237
238          lt = l + 1
239          zlt = zlev(ig, lt)
240          zdz3 = zlev(ig, lt + 1) - zlt
241          zltdwn = zlt - zdz3 / 2
242          zltup = zlt + zdz3 / 2
243
244          !=========================================================================
245          ! 3. Calcul de la flotabilite modifie par melange avec l'air au dessus
246          !=========================================================================
247
248          !--------------------------------------------------
249          IF (iflag_thermals_ed<8) THEN
250            !--------------------------------------------------
251            !AJ052014: J'ai remplac?? la boucle do par un do while
252            ! afin de faire moins de calcul dans la boucle
253            !--------------------------------------------------
254            do while (zlmelup>zltup)
255              lt = lt + 1
256              zlt = zlev(ig, lt)
257              zdz3 = zlev(ig, lt + 1) - zlt
258              zltdwn = zlt - zdz3 / 2
259              zltup = zlt + zdz3 / 2
260            enddo
261            !--------------------------------------------------
262            !AJ052014: Si iflag_thermals_ed<8 (par ex 6), alors
263            ! on cherche o?? se trouve l'altitude d'inversion
264            ! en calculant ztv1 (interpolation de la valeur de
265            ! theta au niveau lt en utilisant les niveaux lt-1 et
266            ! lt-2) et ztv2 (interpolation avec les niveaux lt+1
267            ! et lt+2). Si theta r??ellement calcul??e au niveau lt
268            ! comprise entre ztv1 et ztv2, alors il y a inversion
269            ! et on calcule son altitude zinv en supposant que ztv(lt)
270            ! est une combinaison lineaire de ztv1 et ztv2.
271            ! Ensuite, on calcule la flottabilite en comparant
272            ! la temperature de la couche l a celle de l'air situe
273            ! l+lmel plus haut, ce qui necessite de savoir quel fraction
274            ! de cet air est au-dessus ou en-dessous de l'inversion
275            !--------------------------------------------------
276            atv1 = (ztv(ig, lt - 1) - ztv(ig, lt - 2)) / (zlev(ig, lt - 1) - zlev(ig, lt - 2))
277            btv1 = (ztv(ig, lt - 2) * zlev(ig, lt - 1) - ztv(ig, lt - 1) * zlev(ig, lt - 2)) &
278                    / (zlev(ig, lt - 1) - zlev(ig, lt - 2))
279            atv2 = (ztv(ig, lt + 2) - ztv(ig, lt + 1)) / (zlev(ig, lt + 2) - zlev(ig, lt + 1))
280            btv2 = (ztv(ig, lt + 1) * zlev(ig, lt + 2) - ztv(ig, lt + 2) * zlev(ig, lt + 1)) &
281                    / (zlev(ig, lt + 2) - zlev(ig, lt + 1))
282
283            ztv1 = atv1 * zlt + btv1
284            ztv2 = atv2 * zlt + btv2
285
286            IF (ztv(ig, lt)>ztv1.AND.ztv(ig, lt)<ztv2) THEN
287              !--------------------------------------------------
288              !AJ052014: D??calage de zinv qui est entre le haut
289              !          et le bas de la couche lt
290              !--------------------------------------------------
291              factinv = (ztv2 - ztv(ig, lt)) / (ztv2 - ztv1)
292              zinv = zltdwn + zdz3 * factinv
293
294              IF (zlmeldwn>=zinv) THEN
295                ztv_est(ig, l) = atv2 * zlmel + btv2
296                zbuoyjam(ig, l) = fact_shell * RG * (ztva_est(ig, l) - ztv_est(ig, l)) / ztv_est(ig, l) &
297                        + (1. - fact_shell) * zbuoy(ig, l)
298              elseif (zlmelup>=zinv) THEN
299                ztv_est2 = atv2 * 0.5 * (zlmelup + zinv) + btv2
300                ztv_est1 = atv1 * 0.5 * (zinv + zlmeldwn) + btv1
301                ztv_est(ig, l) = ((zlmelup - zinv) / zdz) * ztv_est2 + ((zinv - zlmeldwn) / zdz) * ztv_est1
302
303                zbuoyjam(ig, l) = fact_shell * RG * (((zlmelup - zinv) / zdz) * (ztva_est(ig, l) - &
304                        ztv_est2) / ztv_est2 + ((zinv - zlmeldwn) / zdz) * (ztva_est(ig, l) - &
305                        ztv_est1) / ztv_est1) + (1. - fact_shell) * zbuoy(ig, l)
306
307              else
308                ztv_est(ig, l) = atv1 * zlmel + btv1
309                zbuoyjam(ig, l) = fact_shell * RG * (ztva_est(ig, l) - ztv_est(ig, l)) / ztv_est(ig, l) &
310                        + (1. - fact_shell) * zbuoy(ig, l)
311              endif
312
313            else ! if (ztv(ig,lt).gt.ztv1.AND.ztv(ig,lt).lt.ztv2) THEN
314              IF (zlmeldwn>zltdwn) THEN
315                zbuoyjam(ig, l) = fact_shell * RG * ((ztva_est(ig, l) - &
316                        ztv(ig, lt)) / ztv(ig, lt)) + (1. - fact_shell) * zbuoy(ig, l)
317              else
318                zbuoyjam(ig, l) = fact_shell * RG * (((zlmelup - zltdwn) / zdz) * (ztva_est(ig, l) - &
319                        ztv(ig, lt)) / ztv(ig, lt) + ((zltdwn - zlmeldwn) / zdz) * (ztva_est(ig, l) - &
320                        ztv(ig, lt - 1)) / ztv(ig, lt - 1)) + (1. - fact_shell) * zbuoy(ig, l)
321
322              endif
323
324              !          zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- &
325              !    &          ztv1)/ztv1+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
326              !    &          ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
327              !         zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- &
328              !    &          po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- &
329              !     &          po(ig,lt-1))/po(ig,lt-1))
330            endif ! if (ztv(ig,lt).gt.ztv1.AND.ztv(ig,lt).lt.ztv2) THEN
331          else  !   if (iflag_thermals_ed.lt.8) THEN
332            lt = l + 1
333            zlt = zlev(ig, lt)
334            zdz2 = zlev(ig, lt) - zlev(ig, l)
335
336            do while (lmel>zdz2)
337              lt = lt + 1
338              zlt = zlev(ig, lt)
339              zdz2 = zlt - zlev(ig, l)
340            enddo
341            zdz3 = zlev(ig, lt + 1) - zlt
342            zltdwn = zlev(ig, lt) - zdz3 / 2
343            zlmelup = zlmel + (zdz / 2)
344            coefzlmel = Min(1., (zlmelup - zltdwn) / zdz)
345            zbuoyjam(ig, l) = 1. * RG * (coefzlmel * (ztva_est(ig, l) - &
346                    ztv(ig, lt)) / ztv(ig, lt) + (1. - coefzlmel) * (ztva_est(ig, l) - &
347                    ztv(ig, lt - 1)) / ztv(ig, lt - 1)) + 0. * zbuoy(ig, l)
348          endif !   if (iflag_thermals_ed.lt.8) THEN
349          !------------------------------------------------
350          !AJAM:nouveau calcul de w?
351          !------------------------------------------------
352          zdz = zlev(ig, l + 1) - zlev(ig, l)
353          zdzbis = zlev(ig, l) - zlev(ig, l - 1)
354          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
355
356          zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
357          zw2factbis = fact_epsilon * 2. * zdzbis / (1. + betalpha)
358          zdw2 = afact * zbuoy(ig, l) / fact_epsilon
359          zdw2bis = afact * zbuoy(ig, l - 1) / fact_epsilon
360          !              zdw2bis=0.5*(zdw2+zdw2bis)
361          lm = Max(1, l - 2)
362          !              zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) &
363          !    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
364          !              zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) &
365          !    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
366          !             w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
367          !             w_est(ig,l+1)=(zdz/zdzbis)*Max(0.0001,exp(-zw2fact)* &
368          !    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
369          !    &                     Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
370          !              w_est(ig,l+1)=Max(0.0001,(1-exp(-zw2fact))*zdw2+w_est(ig,l)*exp(-zw2fact))
371
372          !--------------------------------------------------
373          !AJ052014: J'ai remplac? w_est(ig,l) par zw2(ig,l)
374          !--------------------------------------------------
375          IF (iflag_thermals_ed==8) THEN
376            ! Ancienne version
377            !             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
378            !    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
379            !    &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
380
381            w_est(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (w_est(ig, l) - zdw2) + zdw2)
382
383            ! Nouvelle version Arnaud
384          else
385            !             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
386            !    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
387            !    &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
388
389            w_est(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (w_est(ig, l) - zdw2bis) + zdw2)
390
391            !             w_est(ig,l+1)=Max(0.0001,(zdz/(zdzbis+zdz))*(exp(-zw2fact)* &
392            !    &                     (w_est(ig,l)-zdw2bis)+zdw2)+(zdzbis/(zdzbis+zdz))* &
393            !    &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2bis))
394
395
396
397            !            w_est(ig,l+1)=Max(0.0001,(w_est(ig,l)+zdw2bis*zw2fact)*exp(-zw2fact))
398
399            !             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ &
400            !    &                      (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis))
401
402            !             w_est(ig,l+1)=Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
403
404          endif
405
406          IF (iflag_thermals_ed<6) THEN
407            zalpha = f0(ig) * f_star(ig, l) / sqrt(w_est(ig, l + 1)) / rhobarz(ig, l)
408            !              fact_epsilon=0.0005/(zalpha+0.025)**0.5
409            !              fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
410            fact_epsilon = 0.0002 / (zalpha + 0.1)
411            zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
412            zw2factbis = fact_epsilon * 2. * zdzbis / (1. + betalpha)
413            zdw2 = afact * zbuoy(ig, l) / fact_epsilon
414            zdw2bis = afact * zbuoy(ig, l - 1) / fact_epsilon
415            !              w_est(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
416
417            !              w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
418            !    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
419            !    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
420
421            w_est(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (w_est(ig, l) - zdw2bis) + zdw2)
422
423          endif
424          !--------------------------------------------------
425          !AJ052014: J'ai comment? ce if plus n?cessaire puisqu'
426          !on fait max(0.0001,.....)
427          !--------------------------------------------------
428
429          !             if (w_est(ig,l+1).lt.0.) THEN
430          !               w_est(ig,l+1)=zw2(ig,l)
431          !                w_est(ig,l+1)=0.0001
432          !             endif
433
434        endif
435      enddo
436
437
438      !-------------------------------------------------
439      !calcul des taux d'entrainement et de detrainement
440      !-------------------------------------------------
441
442      do ig = 1, ngrid
443        IF (active(ig)) THEN
444          !          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
445          zw2m = w_est(ig, l + 1)
446          !          zw2m=zw2(ig,l)
447          zdz = zlev(ig, l + 1) - zlev(ig, l)
448          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
449          !          zbuoybis=zbuoy(ig,l)+RG*0.1/300.
450          zbuoybis = zbuoy(ig, l)
451          zalpha = f0(ig) * f_star(ig, l) / sqrt(w_est(ig, l + 1)) / rhobarz(ig, l)
452          zdqt(ig, l) = max(zqta(ig, l - 1) - po(ig, l), 0.) / po(ig, l)
453
454
455          !          entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0.,  &
456          !    &     afact*zbuoybis/zw2m - fact_epsilon )
457
458          !          entr_star(ig,l)=MAX(0.,f_star(ig,l)*zdz*zbetalpha*  &
459          !    &     afact*zbuoybis/zw2m - fact_epsilon )
460
461
462
463          !          zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
464
465          !=========================================================================
466          ! 4. Calcul de l'entrainement et du detrainement
467          !=========================================================================
468
469          !          entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0.,  &
470          !    &     afact*zbuoyjam(ig,l)/zw2m - fact_epsilon )
471          !          entrbis=entr_star(ig,l)
472
473          IF (iflag_thermals_ed<6) THEN
474            fact_epsilon = 0.0002 / (zalpha + 0.1)
475          endif
476
477          detr_star(ig, l) = f_star(ig, l) * zdz             &
478                  * (mix0 * 0.1 / (zalpha + 0.001)               &
479                          + MAX(detr_min, -afact * zbetalpha * zbuoyjam(ig, l) / zw2m   &
480                                  + detr_q_coef * (zdqt(ig, l) / zw2m)**detr_q_power))
481
482          !          detr_star(ig,l)=(zdz/zdzbis)*detr_star(ig,l)+ &
483          !    &                          ((zdzbis-zdz)/zdzbis)*detr_star(ig,l-1)
484
485          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
486
487          entr_star(ig, l) = f_star(ig, l) * zdz * (&
488                  mix0 * 0.1 / (zalpha + 0.001)               &
489                          + zbetalpha * MAX(entr_min, &
490                          afact * zbuoyjam(ig, l) / zw2m - fact_epsilon))
491
492
493          !          entr_star(ig,l)=f_star(ig,l)*zdz* (         &
494          !    &       mix0 * 0.1 / (zalpha+0.001)               &
495          !    &     + MAX(entr_min,                   &
496          !    &     zbetalpha*afact*zbuoyjam(ig,l)/zw2m - fact_epsilon +  &
497          !    &     detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power))
498
499
500          !          entr_star(ig,l)=(zdz/zdzbis)*entr_star(ig,l)+ &
501          !    &                          ((zdzbis-zdz)/zdzbis)*entr_star(ig,l-1)
502
503          !          entr_star(ig,l)=Max(0.,f_star(ig,l)*zdz*zbetalpha*  &
504          !    &     afact*zbuoy(ig,l)/zw2m &
505          !    &     - 1.*fact_epsilon)
506
507
508          ! En dessous de lalim, on prend le max de alim_star et entr_star pour
509          ! alim_star et 0 sinon
510          IF (l<lalim(ig)) THEN
511            alim_star(ig, l) = max(alim_star(ig, l), entr_star(ig, l))
512            entr_star(ig, l) = 0.
513          endif
514          !        if (l.lt.lalim(ig).AND.alim_star(ig,l)>alim_star(ig,l-1)) THEN
515          !          alim_star(ig,l)=entrbis
516          !        endif
517
518          !        PRINT*,'alim0',zlev(ig,l),entr_star(ig,l),detr_star(ig,l),zw2m,zbuoy(ig,l),f_star(ig,l)
519          ! Calcul du flux montant normalise
520          f_star(ig, l + 1) = f_star(ig, l) + alim_star(ig, l) + entr_star(ig, l)  &
521                  - detr_star(ig, l)
522
523        endif
524      enddo
525
526
527      !============================================================================
528      ! 5. calcul de la vitesse verticale en melangeant Tl et qt du thermique
529      !===========================================================================
530
531      activetmp(:) = active(:) .AND. f_star(:, l + 1)>1.e-10
532      do ig = 1, ngrid
533        IF (activetmp(ig)) THEN
534          Zsat = .FALSE.
535          ztla(ig, l) = (f_star(ig, l) * ztla(ig, l - 1) + &
536                  (alim_star(ig, l) + entr_star(ig, l)) * zthl(ig, l))  &
537                  / (f_star(ig, l + 1) + detr_star(ig, l))
538          zqta(ig, l) = (f_star(ig, l) * zqta(ig, l - 1) + &
539                  (alim_star(ig, l) + entr_star(ig, l)) * po(ig, l))  &
540                  / (f_star(ig, l + 1) + detr_star(ig, l))
541
542        endif
543      enddo
544
545      ztemp(:) = zpspsk(:, l) * ztla(:, l)
546      CALL thermcell_qsat(ngrid, activetmp, pplev(:, l), ztemp, zqta(:, l), zqsatth(:, l))
547      do ig = 1, ngrid
548        IF (activetmp(ig)) THEN
549          ! on ecrit de maniere conservative (sat ou non)
550          !          T = Tl +Lv/Cp ql
551          zqla(ig, l) = max(0., zqta(ig, l) - zqsatth(ig, l))
552          ztva(ig, l) = ztla(ig, l) * zpspsk(ig, l) + RLvCp * zqla(ig, l)
553          ztva(ig, l) = ztva(ig, l) / zpspsk(ig, l)
554          !on rajoute le calcul de zha pour diagnostiques (temp potentielle)
555          zha(ig, l) = ztva(ig, l)
556          ztva(ig, l) = ztva(ig, l) * (1. + RETV * (zqta(ig, l)  &
557                  - zqla(ig, l)) - zqla(ig, l))
558          zbuoy(ig, l) = RG * (ztva(ig, l) - ztv(ig, l)) / ztv(ig, l)
559          zdz = zlev(ig, l + 1) - zlev(ig, l)
560          zdzbis = zlev(ig, l) - zlev(ig, l - 1)
561          zeps(ig, l) = (entr_star(ig, l) + alim_star(ig, l)) / (f_star(ig, l) * zdz)
562          !!!!!!!          fact_epsilon=0.002
563          zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
564          zw2factbis = fact_epsilon * 2. * zdzbis / (1. + betalpha)
565          zdw2 = afact * zbuoy(ig, l) / (fact_epsilon)
566          zdw2bis = afact * zbuoy(ig, l - 1) / (fact_epsilon)
567          !              zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) &
568          !    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
569          !              lm=Max(1,l-2)
570          !              zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) &
571          !    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
572          !            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
573          !            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ &
574          !     &                   (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis))
575          !            zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
576          !            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
577          !    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
578          !    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
579          IF (iflag_thermals_ed==8) THEN
580            zw2(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (zw2(ig, l) - zdw2) + zdw2)
581          else
582            zw2(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (zw2(ig, l) - zdw2bis) + zdw2)
583          endif
584          !            zw2(ig,l+1)=Max(0.0001,(zdz/(zdz+zdzbis))*(exp(-zw2fact)* &
585          !    &                     (zw2(ig,l)-zdw2)+zdw2bis)+(zdzbis/(zdz+zdzbis))* &
586          !    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2bis))
587
588          IF (iflag_thermals_ed<6) THEN
589            zalpha = f0(ig) * f_star(ig, l) / sqrt(zw2(ig, l + 1)) / rhobarz(ig, l)
590            !           fact_epsilon=0.0005/(zalpha+0.025)**0.5
591            !           fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
592            fact_epsilon = 0.0002 / (zalpha + 0.1)**1
593            zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
594            zw2factbis = fact_epsilon * 2. * zdzbis / (1. + betalpha)
595            zdw2 = afact * zbuoy(ig, l) / (fact_epsilon)
596            zdw2bis = afact * zbuoy(ig, l - 1) / (fact_epsilon)
597
598            !            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
599            !    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
600            !    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
601            !            zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
602            zw2(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (zw2(ig, l) - zdw2bis) + zdw2)
603
604          endif
605
606        endif
607      enddo
608
609      IF (prt_level>=20) PRINT*, 'coucou calcul detr 460: ig, l', ig, l
610
611      !===========================================================================
612      ! 6. initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
613      !===========================================================================
614
615      nbpb = 0
616      do ig = 1, ngrid
617        IF (zw2(ig, l + 1)>0. .AND. zw2(ig, l + 1)<1.e-10) THEN
618          !               stop 'On tombe sur le cas particulier de thermcell_dry'
619          !               PRINT*,'On tombe sur le cas particulier de thermcell_plume'
620          nbpb = nbpb + 1
621          zw2(ig, l + 1) = 0.
622          linter(ig) = l + 1
623        endif
624
625        IF (zw2(ig, l + 1)<0.) THEN
626          linter(ig) = (l * (zw2(ig, l + 1) - zw2(ig, l))  &
627                  - zw2(ig, l)) / (zw2(ig, l + 1) - zw2(ig, l))
628          zw2(ig, l + 1) = 0.
629          !+CR:04/05/12:correction calcul linter pour calcul de zmax continu
630        elseif (f_star(ig, l + 1)<0.) THEN
631          linter(ig) = (l * (f_star(ig, l + 1) - f_star(ig, l))  &
632                  - f_star(ig, l)) / (f_star(ig, l + 1) - f_star(ig, l))
633          zw2(ig, l + 1) = 0.
634          !fin CR:04/05/12
635        endif
636
637        wa_moy(ig, l + 1) = sqrt(zw2(ig, l + 1))
638
639        IF (wa_moy(ig, l + 1)>wmaxa(ig)) THEN
640          !   lmix est le niveau de la couche ou w (wa_moy) est maximum
641          !on rajoute le calcul de lmix_bis
642          IF (zqla(ig, l)<1.e-10) THEN
643            lmix_bis(ig) = l + 1
644          endif
645          lmix(ig) = l + 1
646          wmaxa(ig) = wa_moy(ig, l + 1)
647        endif
648      enddo
649
650      IF (nbpb>0) THEN
651        PRINT*, 'WARNING on tombe ', nbpb, ' x sur un pb pour l=', l, ' dans thermcell_plume'
652      endif
653
654      !=========================================================================
655      ! FIN DE LA BOUCLE VERTICALE
656    enddo
657    !=========================================================================
658
659    !on recalcule alim_star_tot
660    do ig = 1, ngrid
661      alim_star_tot(ig) = 0.
662    enddo
663    do ig = 1, ngrid
664      do l = 1, lalim(ig) - 1
665        alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig, l)
666      enddo
667    enddo
668
669    IF (prt_level>=20) PRINT*, 'coucou calcul detr 470: ig, l', ig, l
670    RETURN
671  end
672
673
674  SUBROUTINE thermcell_plume_5B(itap, ngrid, nlay, ptimestep, ztv, zthl, po, zl, rhobarz, &
675          zlev, pplev, pphi, zpspsk, alim_star, alim_star_tot, &
676          lalim, f0, detr_star, entr_star, f_star, csc, ztva, &
677          ztla, zqla, zqta, zha, zw2, w_est, ztva_est, zqsatth, lmix, lmix_bis, linter &
678          , lev_out, lunout1, igout)
679    !&           ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
680
681    !--------------------------------------------------------------------------
682    !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
683    ! Version conforme a l'article de Rio et al. 2010.
684    ! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
685    !--------------------------------------------------------------------------
686
687    USE lmdz_thermcell_ini, ONLY: prt_level, fact_thermals_ed_dz, iflag_thermals_ed, RLvCP, RETV, RG
688    USE lmdz_thermcell_qsat, ONLY: thermcell_qsat
689    IMPLICIT NONE
690
691    INTEGER itap
692    INTEGER lunout1, igout
693    INTEGER ngrid, nlay
694    REAL ptimestep
695    REAL ztv(ngrid, nlay)
696    REAL zthl(ngrid, nlay)
697    REAL, INTENT(IN) :: po(ngrid, nlay)
698    REAL zl(ngrid, nlay)
699    REAL rhobarz(ngrid, nlay)
700    REAL zlev(ngrid, nlay + 1)
701    REAL pplev(ngrid, nlay + 1)
702    REAL pphi(ngrid, nlay)
703    REAL zpspsk(ngrid, nlay)
704    REAL alim_star(ngrid, nlay)
705    REAL f0(ngrid)
706    INTEGER lalim(ngrid)
707    INTEGER lev_out                           ! niveau pour les print
708    INTEGER nbpb
709
710    REAL alim_star_tot(ngrid)
711
712    REAL ztva(ngrid, nlay)
713    REAL ztla(ngrid, nlay)
714    REAL zqla(ngrid, nlay)
715    REAL zqta(ngrid, nlay)
716    REAL zha(ngrid, nlay)
717
718    REAL detr_star(ngrid, nlay)
719    REAL coefc
720    REAL entr_star(ngrid, nlay)
721    REAL detr(ngrid, nlay)
722    REAL entr(ngrid, nlay)
723
724    REAL csc(ngrid, nlay)
725
726    REAL zw2(ngrid, nlay + 1)
727    REAL w_est(ngrid, nlay + 1)
728    REAL f_star(ngrid, nlay + 1)
729    REAL wa_moy(ngrid, nlay + 1)
730
731    REAL ztva_est(ngrid, nlay)
732    REAL zqla_est(ngrid, nlay)
733    REAL zqsatth(ngrid, nlay)
734    REAL zta_est(ngrid, nlay)
735    REAL zbuoyjam(ngrid, nlay)
736    REAL ztemp(ngrid), zqsat(ngrid)
737    REAL zdw2
738    REAL zw2modif
739    REAL zw2fact
740    REAL zeps(ngrid, nlay)
741
742    REAL linter(ngrid)
743    INTEGER lmix(ngrid)
744    INTEGER lmix_bis(ngrid)
745    REAL    wmaxa(ngrid)
746
747    INTEGER ig, l, k
748
749    REAL zdz, zbuoy(ngrid, nlay), zalpha, gamma(ngrid, nlay), zdqt(ngrid, nlay), zw2m
750    REAL zbuoybis
751    REAL zcor, zdelta, zcvm5, qlbef, zdz2
752    REAL betalpha, zbetalpha
753    REAL eps, afact
754    LOGICAL Zsat
755    LOGICAL active(ngrid), activetmp(ngrid)
756    REAL fact_gamma, fact_epsilon, fact_gamma2, fact_epsilon2
757    REAL c2(ngrid, nlay)
758    Zsat = .FALSE.
759    ! Initialisation
760
761    fact_epsilon = 0.002
762    betalpha = 0.9
763    afact = 2. / 3.
764
765    zbetalpha = betalpha / (1. + betalpha)
766
767
768    ! Initialisations des variables reeles
769    IF (1==1) THEN
770      ztva(:, :) = ztv(:, :)
771      ztva_est(:, :) = ztva(:, :)
772      ztla(:, :) = zthl(:, :)
773      zqta(:, :) = po(:, :)
774      zha(:, :) = ztva(:, :)
775    else
776      ztva(:, :) = 0.
777      ztva_est(:, :) = 0.
778      ztla(:, :) = 0.
779      zqta(:, :) = 0.
780      zha(:, :) = 0.
781    endif
782
783    zqla_est(:, :) = 0.
784    zqsatth(:, :) = 0.
785    zqla(:, :) = 0.
786    detr_star(:, :) = 0.
787    entr_star(:, :) = 0.
788    alim_star(:, :) = 0.
789    alim_star_tot(:) = 0.
790    csc(:, :) = 0.
791    detr(:, :) = 0.
792    entr(:, :) = 0.
793    zw2(:, :) = 0.
794    zbuoy(:, :) = 0.
795    zbuoyjam(:, :) = 0.
796    gamma(:, :) = 0.
797    zeps(:, :) = 0.
798    w_est(:, :) = 0.
799    f_star(:, :) = 0.
800    wa_moy(:, :) = 0.
801    linter(:) = 1.
802    !     linter(:)=1.
803    ! Initialisation des variables entieres
804    lmix(:) = 1
805    lmix_bis(:) = 2
806    wmaxa(:) = 0.
807    lalim(:) = 1
808
809
810    !-------------------------------------------------------------------------
811    ! On ne considere comme actif que les colonnes dont les deux premieres
812    ! couches sont instables.
813    !-------------------------------------------------------------------------
814    active(:) = ztv(:, 1)>ztv(:, 2)
815
816    !-------------------------------------------------------------------------
817    ! Definition de l'alimentation
818    !-------------------------------------------------------------------------
819    do l = 1, nlay - 1
820      do ig = 1, ngrid
821        IF (ztv(ig, l)> ztv(ig, l + 1) .AND. ztv(ig, 1)>=ztv(ig, l)) THEN
822          alim_star(ig, l) = MAX((ztv(ig, l) - ztv(ig, l + 1)), 0.)  &
823                  * sqrt(zlev(ig, l + 1))
824          lalim(ig) = l + 1
825          alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig, l)
826        endif
827      enddo
828    enddo
829    do l = 1, nlay
830      do ig = 1, ngrid
831        IF (alim_star_tot(ig) > 1.e-10) THEN
832          alim_star(ig, l) = alim_star(ig, l) / alim_star_tot(ig)
833        endif
834      enddo
835    enddo
836    alim_star_tot(:) = 1.
837
838
839
840    !------------------------------------------------------------------------------
841    ! Calcul dans la premiere couche
842    ! On decide dans cette version que le thermique n'est actif que si la premiere
843    ! couche est instable.
844    ! Pourrait etre change si on veut que le thermiques puisse se d??clencher
845    ! dans une couche l>1
846    !------------------------------------------------------------------------------
847    do ig = 1, ngrid
848      ! Le panache va prendre au debut les caracteristiques de l'air contenu
849      ! dans cette couche.
850      IF (active(ig)) THEN
851        ztla(ig, 1) = zthl(ig, 1)
852        zqta(ig, 1) = po(ig, 1)
853        zqla(ig, 1) = zl(ig, 1)
854        !cr: attention, prise en compte de f*(1)=1
855        f_star(ig, 2) = alim_star(ig, 1)
856        zw2(ig, 2) = 2. * RG * (ztv(ig, 1) - ztv(ig, 2)) / ztv(ig, 2)  &
857                & * (zlev(ig, 2) - zlev(ig, 1))  &
858                & * 0.4 * pphi(ig, 1) / (pphi(ig, 2) - pphi(ig, 1))
859        w_est(ig, 2) = zw2(ig, 2)
860      endif
861    enddo
862
863    !==============================================================================
864    !boucle de calcul de la vitesse verticale dans le thermique
865    !==============================================================================
866    do l = 2, nlay - 1
867      !==============================================================================
868
869
870      ! On decide si le thermique est encore actif ou non
871      ! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
872      do ig = 1, ngrid
873        active(ig) = active(ig) &
874                &                 .AND. zw2(ig, l)>1.e-10 &
875                &                 .AND. f_star(ig, l) + alim_star(ig, l)>1.e-10
876      enddo
877
878
879
880      !---------------------------------------------------------------------------
881      ! calcul des proprietes thermodynamiques et de la vitesse de la couche l
882      ! sans tenir compte du detrainement et de l'entrainement dans cette
883      ! couche
884      ! C'est a dire qu'on suppose
885      ! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
886      ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
887      ! avant) a l'alimentation pour avoir un calcul plus propre
888      !---------------------------------------------------------------------------
889
890      ztemp(:) = zpspsk(:, l) * ztla(:, l - 1)
891      CALL thermcell_qsat(ngrid, active, pplev(:, l), ztemp, zqta(:, l - 1), zqsat(:))
892
893      do ig = 1, ngrid
894        !       PRINT*,'active',active(ig),ig,l
895        IF(active(ig)) THEN
896          zqla_est(ig, l) = max(0., zqta(ig, l - 1) - zqsat(ig))
897          ztva_est(ig, l) = ztla(ig, l - 1) * zpspsk(ig, l) + RLvCp * zqla_est(ig, l)
898          zta_est(ig, l) = ztva_est(ig, l)
899          ztva_est(ig, l) = ztva_est(ig, l) / zpspsk(ig, l)
900          ztva_est(ig, l) = ztva_est(ig, l) * (1. + RETV * (zqta(ig, l - 1)  &
901                  - zqla_est(ig, l)) - zqla_est(ig, l))
902
903          !------------------------------------------------
904          !AJAM:nouveau calcul de w?
905          !------------------------------------------------
906          zdz = zlev(ig, l + 1) - zlev(ig, l)
907          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
908
909          zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
910          zdw2 = (afact) * zbuoy(ig, l) / (fact_epsilon)
911          w_est(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (w_est(ig, l) - zdw2) + zdw2)
912
913          IF (w_est(ig, l + 1)<0.) THEN
914            w_est(ig, l + 1) = zw2(ig, l)
915          endif
916        endif
917      enddo
918
919
920      !-------------------------------------------------
921      !calcul des taux d'entrainement et de detrainement
922      !-------------------------------------------------
923
924      do ig = 1, ngrid
925        IF (active(ig)) THEN
926          zw2m = max(0.5 * (w_est(ig, l) + w_est(ig, l + 1)), 0.1)
927          zw2m = w_est(ig, l + 1)
928          zdz = zlev(ig, l + 1) - zlev(ig, l)
929          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
930          !          zbuoybis=zbuoy(ig,l)+RG*0.1/300.
931          zbuoybis = zbuoy(ig, l)
932          zalpha = f0(ig) * f_star(ig, l) / sqrt(w_est(ig, l + 1)) / rhobarz(ig, l)
933          zdqt(ig, l) = max(zqta(ig, l - 1) - po(ig, l), 0.) / po(ig, l)
934
935          entr_star(ig, l) = f_star(ig, l) * zdz * zbetalpha * MAX(0., &
936                  afact * zbuoybis / zw2m - fact_epsilon)
937
938          detr_star(ig, l) = f_star(ig, l) * zdz                        &
939                  * MAX(1.e-3, -afact * zbetalpha * zbuoy(ig, l) / zw2m          &
940                          + 0.012 * (zdqt(ig, l) / zw2m)**0.5)
941
942          ! En dessous de lalim, on prend le max de alim_star et entr_star pour
943          ! alim_star et 0 sinon
944          IF (l<lalim(ig)) THEN
945            alim_star(ig, l) = max(alim_star(ig, l), entr_star(ig, l))
946            entr_star(ig, l) = 0.
947          endif
948
949          ! Calcul du flux montant normalise
950          f_star(ig, l + 1) = f_star(ig, l) + alim_star(ig, l) + entr_star(ig, l)  &
951                  - detr_star(ig, l)
952
953        endif
954      enddo
955
956
957      !----------------------------------------------------------------------------
958      !calcul de la vitesse verticale en melangeant Tl et qt du thermique
959      !---------------------------------------------------------------------------
960      activetmp(:) = active(:) .AND. f_star(:, l + 1)>1.e-10
961      do ig = 1, ngrid
962        IF (activetmp(ig)) THEN
963          Zsat = .FALSE.
964          ztla(ig, l) = (f_star(ig, l) * ztla(ig, l - 1) + &
965                  (alim_star(ig, l) + entr_star(ig, l)) * zthl(ig, l))  &
966                  / (f_star(ig, l + 1) + detr_star(ig, l))
967          zqta(ig, l) = (f_star(ig, l) * zqta(ig, l - 1) + &
968                  (alim_star(ig, l) + entr_star(ig, l)) * po(ig, l))  &
969                  / (f_star(ig, l + 1) + detr_star(ig, l))
970
971        endif
972      enddo
973
974      ztemp(:) = zpspsk(:, l) * ztla(:, l)
975      CALL thermcell_qsat(ngrid, activetmp, pplev(:, l), ztemp, zqta(:, l), zqsatth(:, l))
976
977      do ig = 1, ngrid
978        IF (activetmp(ig)) THEN
979          ! on ecrit de maniere conservative (sat ou non)
980          !          T = Tl +Lv/Cp ql
981          zqla(ig, l) = max(0., zqta(ig, l) - zqsatth(ig, l))
982          ztva(ig, l) = ztla(ig, l) * zpspsk(ig, l) + RLvCp * zqla(ig, l)
983          ztva(ig, l) = ztva(ig, l) / zpspsk(ig, l)
984          !on rajoute le calcul de zha pour diagnostiques (temp potentielle)
985          zha(ig, l) = ztva(ig, l)
986          ztva(ig, l) = ztva(ig, l) * (1. + RETV * (zqta(ig, l)  &
987                  - zqla(ig, l)) - zqla(ig, l))
988          zbuoy(ig, l) = RG * (ztva(ig, l) - ztv(ig, l)) / ztv(ig, l)
989          zdz = zlev(ig, l + 1) - zlev(ig, l)
990          zeps(ig, l) = (entr_star(ig, l) + alim_star(ig, l)) / (f_star(ig, l) * zdz)
991
992          zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
993          zdw2 = afact * zbuoy(ig, l) / (fact_epsilon)
994          zw2(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (zw2(ig, l) - zdw2) + zdw2)
995        endif
996      enddo
997
998      IF (prt_level>=20) PRINT*, 'coucou calcul detr 460: ig, l', ig, l
999
1000      !---------------------------------------------------------------------------
1001      !initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
1002      !---------------------------------------------------------------------------
1003
1004      nbpb = 0
1005      do ig = 1, ngrid
1006        IF (zw2(ig, l + 1)>0. .AND. zw2(ig, l + 1)<1.e-10) THEN
1007          !               stop 'On tombe sur le cas particulier de thermcell_dry'
1008          !               PRINT*,'On tombe sur le cas particulier de thermcell_plume'
1009          nbpb = nbpb + 1
1010          zw2(ig, l + 1) = 0.
1011          linter(ig) = l + 1
1012        endif
1013
1014        IF (zw2(ig, l + 1)<0.) THEN
1015          linter(ig) = (l * (zw2(ig, l + 1) - zw2(ig, l))  &
1016                  - zw2(ig, l)) / (zw2(ig, l + 1) - zw2(ig, l))
1017          zw2(ig, l + 1) = 0.
1018        elseif (f_star(ig, l + 1)<0.) THEN
1019          linter(ig) = (l * (f_star(ig, l + 1) - f_star(ig, l))  &
1020                  - f_star(ig, l)) / (f_star(ig, l + 1) - f_star(ig, l))
1021          !           PRINT*,"linter plume", linter(ig)
1022          zw2(ig, l + 1) = 0.
1023        endif
1024
1025        wa_moy(ig, l + 1) = sqrt(zw2(ig, l + 1))
1026
1027        IF (wa_moy(ig, l + 1)>wmaxa(ig)) THEN
1028          !   lmix est le niveau de la couche ou w (wa_moy) est maximum
1029          !on rajoute le calcul de lmix_bis
1030          IF (zqla(ig, l)<1.e-10) THEN
1031            lmix_bis(ig) = l + 1
1032          endif
1033          lmix(ig) = l + 1
1034          wmaxa(ig) = wa_moy(ig, l + 1)
1035        endif
1036      enddo
1037
1038      IF (nbpb>0) THEN
1039        PRINT*, 'WARNING on tombe ', nbpb, ' x sur un pb pour l=', l, ' dans thermcell_plume'
1040      endif
1041
1042      !=========================================================================
1043      ! FIN DE LA BOUCLE VERTICALE
1044    enddo
1045    !=========================================================================
1046
1047    !on recalcule alim_star_tot
1048    do ig = 1, ngrid
1049      alim_star_tot(ig) = 0.
1050    enddo
1051    do ig = 1, ngrid
1052      do l = 1, lalim(ig) - 1
1053        alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig, l)
1054      enddo
1055    enddo
1056
1057    IF (prt_level>=20) PRINT*, 'coucou calcul detr 470: ig, l', ig, l
1058  END
1059END MODULE lmdz_thermcell_plume_6A
Note: See TracBrowser for help on using the repository browser.