source: LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_plume.f90 @ 5136

Last change on this file since 5136 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
File size: 17.5 KB
Line 
1MODULE lmdz_thermcell_plume
2
3  ! $Id: thermcell_plume.F90 3074 2017-11-15 13:31:44Z fhourdin $
4
5CONTAINS
6
7  SUBROUTINE thermcell_plume(itap, ngrid, nlay, ptimestep, ztv, zthl, po, zl, rhobarz, &
8          zlev, pplev, pphi, zpspsk, alim_star, alim_star_tot, &
9          lalim, f0, detr_star, entr_star, f_star, csc, ztva, &
10          ztla, zqla, zqta, zha, zw2, w_est, ztva_est, zqsatth, lmix, lmix_bis, linter &
11          , lev_out, lunout1, igout)
12    !     &           ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
13    !--------------------------------------------------------------------------
14    ! Auhtors : Catherine Rio, Frédéric Hourdin, Arnaud Jam
15
16    !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
17    !   This versions starts from a cleaning of thermcell_plume_6A (2019/01/20)
18    !   thermcell_plume_6A is activate for flag_thermas_ed < 10
19    !   thermcell_plume_5B for flag_thermas_ed < 20
20    !   thermcell_plume for flag_thermals_ed>= 20
21    !   Various options are controled by the flag_thermals_ed parameter
22    !   = 20 : equivalent to thermcell_plume_6A with flag_thermals_ed=8
23    !   = 21 : the Jam strato-cumulus modif is not activated in detrainment
24    !   = 29 : an other way to compute the modified buoyancy (to be tested)
25    !--------------------------------------------------------------------------
26
27    USE lmdz_thermcell_ini, ONLY: prt_level, fact_thermals_ed_dz, iflag_thermals_ed, RLvCP, RETV, RG
28    USE lmdz_thermcell_ini, ONLY: fact_epsilon, betalpha, afact, fact_shell
29    USE lmdz_thermcell_ini, ONLY: detr_min, entr_min, detr_q_coef, detr_q_power
30    USE lmdz_thermcell_ini, ONLY: mix0, thermals_flag_alim
31    USE lmdz_thermcell_alim, ONLY: thermcell_alim
32    USE lmdz_thermcell_qsat, ONLY: thermcell_qsat
33
34    IMPLICIT NONE
35
36    INTEGER, INTENT(IN) :: itap, lev_out, lunout1, igout, ngrid, nlay
37    REAL, INTENT(IN) :: ptimestep
38    REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: ztv
39    REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: zthl
40    REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: po
41    REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: zl
42    REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: rhobarz
43    REAL, INTENT(IN), DIMENSION(ngrid, nlay + 1) :: zlev
44    REAL, INTENT(IN), DIMENSION(ngrid, nlay + 1) :: pplev
45    REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: pphi
46    REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: zpspsk
47    REAL, INTENT(IN), DIMENSION(ngrid) :: f0
48
49    INTEGER, INTENT(OUT) :: lalim(ngrid)
50    REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: alim_star
51    REAL, INTENT(OUT), DIMENSION(ngrid) :: alim_star_tot
52    REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: detr_star
53    REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: entr_star
54    REAL, INTENT(OUT), DIMENSION(ngrid, nlay + 1) :: f_star
55    REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: csc
56    REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: ztva
57    REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: ztla
58    REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: zqla
59    REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: zqta
60    REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: zha
61    REAL, INTENT(OUT), DIMENSION(ngrid, nlay + 1) :: zw2
62    REAL, INTENT(OUT), DIMENSION(ngrid, nlay + 1) :: w_est
63    REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: ztva_est
64    REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: zqsatth
65    INTEGER, INTENT(OUT), DIMENSION(ngrid) :: lmix(ngrid)
66    INTEGER, INTENT(OUT), DIMENSION(ngrid) :: lmix_bis(ngrid)
67    REAL, INTENT(OUT), DIMENSION(ngrid) :: linter(ngrid)
68
69    REAL, DIMENSION(ngrid, nlay + 1) :: wa_moy
70    REAL, DIMENSION(ngrid, nlay) :: entr, detr
71    REAL, DIMENSION(ngrid, nlay) :: ztv_est
72    REAL, DIMENSION(ngrid, nlay) :: zqla_est
73    REAL, DIMENSION(ngrid, nlay) :: zta_est
74    REAL, DIMENSION(ngrid) :: ztemp, zqsat
75    REAL zdw2, zdw2bis
76    REAL zw2modif
77    REAL zw2fact, zw2factbis
78    REAL, DIMENSION(ngrid, nlay) :: zeps
79
80    REAL, DIMENSION(ngrid) :: wmaxa
81
82    INTEGER ig, l, k, lt, it, lm, nbpb
83
84    REAL, DIMENSION(ngrid, nlay) :: zbuoy, gamma, zdqt
85    REAL zdz, zalpha, zw2m
86    REAL, DIMENSION(ngrid, nlay) :: zbuoyjam, zdqtjam
87    REAL zdz2, zdz3, lmel, entrbis, zdzbis
88    REAL, DIMENSION(ngrid) :: d_temp
89    REAL ztv1, ztv2, factinv, zinv, zlmel
90    REAL zlmelup, zlmeldwn, zlt, zltdwn, zltup
91    REAL atv1, atv2, btv1, btv2
92    REAL ztv_est1, ztv_est2
93    REAL zcor, zdelta, zcvm5, qlbef
94    REAL zbetalpha, coefzlmel
95    REAL eps
96    LOGICAL Zsat
97    LOGICAL, DIMENSION(ngrid) :: active, activetmp
98    REAL fact_gamma, fact_gamma2, fact_epsilon2
99
100    REAL, DIMENSION(ngrid, nlay) :: c2
101
102    IF (ngrid==1) PRINT*, 'THERMCELL PLUME MODIFIE 2014/07/11'
103    Zsat = .FALSE.
104    ! Initialisation
105
106    zbetalpha = betalpha / (1. + betalpha)
107
108
109    ! Initialisations des variables r?elles
110    IF (1==1) THEN
111      ztva(:, :) = ztv(:, :)
112      ztva_est(:, :) = ztva(:, :)
113      ztv_est(:, :) = ztv(:, :)
114      ztla(:, :) = zthl(:, :)
115      zqta(:, :) = po(:, :)
116      zqla(:, :) = 0.
117      zha(:, :) = ztva(:, :)
118    else
119      ztva(:, :) = 0.
120      ztv_est(:, :) = 0.
121      ztva_est(:, :) = 0.
122      ztla(:, :) = 0.
123      zqta(:, :) = 0.
124      zha(:, :) = 0.
125    endif
126
127    zqla_est(:, :) = 0.
128    zqsatth(:, :) = 0.
129    zqla(:, :) = 0.
130    detr_star(:, :) = 0.
131    entr_star(:, :) = 0.
132    alim_star(:, :) = 0.
133    alim_star_tot(:) = 0.
134    csc(:, :) = 0.
135    detr(:, :) = 0.
136    entr(:, :) = 0.
137    zw2(:, :) = 0.
138    zbuoy(:, :) = 0.
139    zbuoyjam(:, :) = 0.
140    gamma(:, :) = 0.
141    zeps(:, :) = 0.
142    w_est(:, :) = 0.
143    f_star(:, :) = 0.
144    wa_moy(:, :) = 0.
145    linter(:) = 1.
146    !     linter(:)=1.
147    ! Initialisation des variables entieres
148    lmix(:) = 1
149    lmix_bis(:) = 2
150    wmaxa(:) = 0.
151
152
153    !-------------------------------------------------------------------------
154    ! On ne considere comme actif que les colonnes dont les deux premieres
155    ! couches sont instables.
156    !-------------------------------------------------------------------------
157
158    active(:) = ztv(:, 1)>ztv(:, 2)
159    d_temp(:) = 0. ! Pour activer un contraste de temperature a la base
160    ! du panache
161    !  Cet appel pourrait être fait avant thermcell_plume dans thermcell_main
162    CALL thermcell_alim(thermals_flag_alim, ngrid, nlay, ztv, d_temp, zlev, alim_star, lalim)
163
164    !------------------------------------------------------------------------------
165    ! Calcul dans la premiere couche
166    ! On decide dans cette version que le thermique n'est actif que si la premiere
167    ! couche est instable.
168    ! Pourrait etre change si on veut que le thermiques puisse se d??clencher
169    ! dans une couche l>1
170    !------------------------------------------------------------------------------
171    do ig = 1, ngrid
172      ! Le panache va prendre au debut les caracteristiques de l'air contenu
173      ! dans cette couche.
174      IF (active(ig)) THEN
175        ztla(ig, 1) = zthl(ig, 1)
176        zqta(ig, 1) = po(ig, 1)
177        zqla(ig, 1) = zl(ig, 1)
178        !cr: attention, prise en compte de f*(1)=1
179        f_star(ig, 2) = alim_star(ig, 1)
180        zw2(ig, 2) = 2. * RG * (ztv(ig, 1) - ztv(ig, 2)) / ztv(ig, 2)  &
181                & * (zlev(ig, 2) - zlev(ig, 1))  &
182                & * 0.4 * pphi(ig, 1) / (pphi(ig, 2) - pphi(ig, 1))
183        w_est(ig, 2) = zw2(ig, 2)
184      endif
185    enddo
186
187    !==============================================================================
188    !boucle de calcul de la vitesse verticale dans le thermique
189    !==============================================================================
190    do l = 2, nlay - 1
191      !==============================================================================
192
193
194      ! On decide si le thermique est encore actif ou non
195      ! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
196      do ig = 1, ngrid
197        active(ig) = active(ig) &
198                &                 .AND. zw2(ig, l)>1.e-10 &
199                &                 .AND. f_star(ig, l) + alim_star(ig, l)>1.e-10
200      enddo
201
202
203
204      !---------------------------------------------------------------------------
205      ! calcul des proprietes thermodynamiques et de la vitesse de la couche l
206      ! sans tenir compte du detrainement et de l'entrainement dans cette
207      ! couche
208      ! C'est a dire qu'on suppose
209      ! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
210      ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
211      ! avant) a l'alimentation pour avoir un calcul plus propre
212      !---------------------------------------------------------------------------
213
214      ztemp(:) = zpspsk(:, l) * ztla(:, l - 1)
215      CALL thermcell_qsat(ngrid, active, pplev(:, l), ztemp, zqta(:, l - 1), zqsat(:))
216      do ig = 1, ngrid
217        !       PRINT*,'active',active(ig),ig,l
218        IF(active(ig)) THEN
219          zqla_est(ig, l) = max(0., zqta(ig, l - 1) - zqsat(ig))
220          ztva_est(ig, l) = ztla(ig, l - 1) * zpspsk(ig, l) + RLvCp * zqla_est(ig, l)
221          zta_est(ig, l) = ztva_est(ig, l)
222          ztva_est(ig, l) = ztva_est(ig, l) / zpspsk(ig, l)
223          ztva_est(ig, l) = ztva_est(ig, l) * (1. + RETV * (zqta(ig, l - 1)  &
224                  - zqla_est(ig, l)) - zqla_est(ig, l))
225
226
227          !Modif AJAM
228
229          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
230          zdz = zlev(ig, l + 1) - zlev(ig, l)
231          lmel = fact_thermals_ed_dz * zlev(ig, l)
232          !        lmel=0.09*zlev(ig,l)
233          zlmel = zlev(ig, l) + lmel
234          zlmelup = zlmel + (zdz / 2)
235          zlmeldwn = zlmel - (zdz / 2)
236
237          lt = l + 1
238          zlt = zlev(ig, lt)
239          zdz3 = zlev(ig, lt + 1) - zlt
240          zltdwn = zlt - zdz3 / 2
241          zltup = zlt + zdz3 / 2
242
243          !=========================================================================
244          ! 3. Calcul de la flotabilite modifie par melange avec l'air au dessus
245          !=========================================================================
246
247          !--------------------------------------------------
248          lt = l + 1
249          zlt = zlev(ig, lt)
250          zdz2 = zlev(ig, lt) - zlev(ig, l)
251
252          do while (lmel>zdz2)
253            lt = lt + 1
254            zlt = zlev(ig, lt)
255            zdz2 = zlt - zlev(ig, l)
256          enddo
257          zdz3 = zlev(ig, lt + 1) - zlt
258          zltdwn = zlev(ig, lt) - zdz3 / 2
259          zlmelup = zlmel + (zdz / 2)
260          coefzlmel = Min(1., (zlmelup - zltdwn) / zdz)
261          zbuoyjam(ig, l) = 1. * RG * (coefzlmel * (ztva_est(ig, l) - &
262                  ztv(ig, lt)) / ztv(ig, lt) + (1. - coefzlmel) * (ztva_est(ig, l) - &
263                  ztv(ig, lt - 1)) / ztv(ig, lt - 1)) + 0. * zbuoy(ig, l)
264
265          !------------------------------------------------
266          !AJAM:nouveau calcul de w?
267          !------------------------------------------------
268          zdz = zlev(ig, l + 1) - zlev(ig, l)
269          zdzbis = zlev(ig, l) - zlev(ig, l - 1)
270          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
271          zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
272          zw2factbis = fact_epsilon * 2. * zdzbis / (1. + betalpha)
273          zdw2 = afact * zbuoy(ig, l) / fact_epsilon
274          zdw2bis = afact * zbuoy(ig, l - 1) / fact_epsilon
275          lm = Max(1, l - 2)
276          w_est(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (w_est(ig, l) - zdw2) + zdw2)
277        endif
278      enddo
279
280
281      !-------------------------------------------------
282      !calcul des taux d'entrainement et de detrainement
283      !-------------------------------------------------
284
285      do ig = 1, ngrid
286        IF (active(ig)) THEN
287          !          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
288          zw2m = w_est(ig, l + 1)
289          zdz = zlev(ig, l + 1) - zlev(ig, l)
290          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
291          zalpha = f0(ig) * f_star(ig, l) / sqrt(w_est(ig, l + 1)) / rhobarz(ig, l)
292          zdqt(ig, l) = max(zqta(ig, l - 1) - po(ig, l), 0.) / po(ig, l)
293
294          !=========================================================================
295          ! 4. Calcul de l'entrainement et du detrainement
296          !=========================================================================
297
298          detr_star(ig, l) = f_star(ig, l) * zdz             &
299                  * (mix0 * 0.1 / (zalpha + 0.001)               &
300                          + MAX(detr_min, -afact * zbetalpha * zbuoyjam(ig, l) / zw2m   &
301                                  + detr_q_coef * (zdqt(ig, l) / zw2m)**detr_q_power))
302
303          IF (iflag_thermals_ed == 20) THEN
304            entr_star(ig, l) = f_star(ig, l) * zdz * (&
305                    mix0 * 0.1 / (zalpha + 0.001)               &
306                            + zbetalpha * MAX(entr_min, &
307                            afact * zbuoyjam(ig, l) / zw2m - fact_epsilon))
308          else
309            entr_star(ig, l) = f_star(ig, l) * zdz * (&
310                    mix0 * 0.1 / (zalpha + 0.001)               &
311                            + zbetalpha * MAX(entr_min, &
312                            afact * zbuoy(ig, l) / zw2m - fact_epsilon))
313          endif
314
315          ! En dessous de lalim, on prend le max de alim_star et entr_star pour
316          ! alim_star et 0 sinon
317          IF (l<lalim(ig)) THEN
318            alim_star(ig, l) = max(alim_star(ig, l), entr_star(ig, l))
319            entr_star(ig, l) = 0.
320          endif
321          f_star(ig, l + 1) = f_star(ig, l) + alim_star(ig, l) + entr_star(ig, l)  &
322                  - detr_star(ig, l)
323
324        endif
325      enddo
326
327
328      !============================================================================
329      ! 5. calcul de la vitesse verticale en melangeant Tl et qt du thermique
330      !===========================================================================
331
332      activetmp(:) = active(:) .AND. f_star(:, l + 1)>1.e-10
333      do ig = 1, ngrid
334        IF (activetmp(ig)) THEN
335          Zsat = .FALSE.
336          ztla(ig, l) = (f_star(ig, l) * ztla(ig, l - 1) + &
337                  (alim_star(ig, l) + entr_star(ig, l)) * zthl(ig, l))  &
338                  / (f_star(ig, l + 1) + detr_star(ig, l))
339          zqta(ig, l) = (f_star(ig, l) * zqta(ig, l - 1) + &
340                  (alim_star(ig, l) + entr_star(ig, l)) * po(ig, l))  &
341                  / (f_star(ig, l + 1) + detr_star(ig, l))
342
343        endif
344      enddo
345
346      ztemp(:) = zpspsk(:, l) * ztla(:, l)
347      CALL thermcell_qsat(ngrid, activetmp, pplev(:, l), ztemp, zqta(:, l), zqsatth(:, l))
348      do ig = 1, ngrid
349        IF (activetmp(ig)) THEN
350          ! on ecrit de maniere conservative (sat ou non)
351          !          T = Tl +Lv/Cp ql
352          zqla(ig, l) = max(0., zqta(ig, l) - zqsatth(ig, l))
353          ztva(ig, l) = ztla(ig, l) * zpspsk(ig, l) + RLvCp * zqla(ig, l)
354          ztva(ig, l) = ztva(ig, l) / zpspsk(ig, l)
355          !on rajoute le calcul de zha pour diagnostiques (temp potentielle)
356          zha(ig, l) = ztva(ig, l)
357          ztva(ig, l) = ztva(ig, l) * (1. + RETV * (zqta(ig, l)  &
358                  - zqla(ig, l)) - zqla(ig, l))
359          zbuoy(ig, l) = RG * (ztva(ig, l) - ztv(ig, l)) / ztv(ig, l)
360          zdz = zlev(ig, l + 1) - zlev(ig, l)
361          zdzbis = zlev(ig, l) - zlev(ig, l - 1)
362          zeps(ig, l) = (entr_star(ig, l) + alim_star(ig, l)) / (f_star(ig, l) * zdz)
363          zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
364          zw2factbis = fact_epsilon * 2. * zdzbis / (1. + betalpha)
365          zdw2 = afact * zbuoy(ig, l) / (fact_epsilon)
366          zdw2bis = afact * zbuoy(ig, l - 1) / (fact_epsilon)
367          zw2(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (zw2(ig, l) - zdw2) + zdw2)
368        endif
369      enddo
370
371      IF (prt_level>=20) PRINT*, 'coucou calcul detr 460: ig, l', ig, l
372
373      !===========================================================================
374      ! 6. initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
375      !===========================================================================
376
377      nbpb = 0
378      do ig = 1, ngrid
379        IF (zw2(ig, l + 1)>0. .AND. zw2(ig, l + 1)<1.e-10) THEN
380          !               stop 'On tombe sur le cas particulier de thermcell_dry'
381          !               PRINT*,'On tombe sur le cas particulier de thermcell_plume'
382          nbpb = nbpb + 1
383          zw2(ig, l + 1) = 0.
384          linter(ig) = l + 1
385        endif
386
387        IF (zw2(ig, l + 1)<0.) THEN
388          linter(ig) = (l * (zw2(ig, l + 1) - zw2(ig, l))  &
389                  - zw2(ig, l)) / (zw2(ig, l + 1) - zw2(ig, l))
390          zw2(ig, l + 1) = 0.
391          !+CR:04/05/12:correction calcul linter pour calcul de zmax continu
392        elseif (f_star(ig, l + 1)<0.) THEN
393          linter(ig) = (l * (f_star(ig, l + 1) - f_star(ig, l))  &
394                  - f_star(ig, l)) / (f_star(ig, l + 1) - f_star(ig, l))
395          zw2(ig, l + 1) = 0.
396          !fin CR:04/05/12
397        endif
398
399        wa_moy(ig, l + 1) = sqrt(zw2(ig, l + 1))
400
401        IF (wa_moy(ig, l + 1)>wmaxa(ig)) THEN
402          !   lmix est le niveau de la couche ou w (wa_moy) est maximum
403          !on rajoute le calcul de lmix_bis
404          IF (zqla(ig, l)<1.e-10) THEN
405            lmix_bis(ig) = l + 1
406          endif
407          lmix(ig) = l + 1
408          wmaxa(ig) = wa_moy(ig, l + 1)
409        endif
410      enddo
411
412      IF (nbpb>0) THEN
413        PRINT*, 'WARNING on tombe ', nbpb, ' x sur un pb pour l=', l, ' dans thermcell_plume'
414      endif
415
416      !=========================================================================
417      ! FIN DE LA BOUCLE VERTICALE
418    enddo
419    !=========================================================================
420
421    !on recalcule alim_star_tot
422    do ig = 1, ngrid
423      alim_star_tot(ig) = 0.
424    enddo
425    do ig = 1, ngrid
426      do l = 1, lalim(ig) - 1
427        alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig, l)
428      enddo
429    enddo
430
431    IF (prt_level>=20) PRINT*, 'coucou calcul detr 470: ig, l', ig, l
432    RETURN
433  END
434END MODULE lmdz_thermcell_plume
Note: See TracBrowser for help on using the repository browser.