source: LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_alp.F90 @ 5151

Last change on this file since 5151 was 5144, checked in by abarral, 7 weeks ago

Put YOMCST.h into modules
File size: 16.2 KB
Line 
1MODULE lmdz_thermcell_alp
2  ! $Id: thermcell_main.F90 2351 2015-08-25 15:14:59Z emillour $
3
4CONTAINS
5
6  SUBROUTINE thermcell_alp(ngrid, nlay, ptimestep  &                         ! in
7          , pplay, pplev  &                                        ! in
8          , fm0, entr0, lmax  &                                     ! in
9          , pbl_tke, pctsrf, omega, airephy &                        ! in
10          , zw2, fraca &                                           ! in
11          , pcon, rhobarz, wth3, wmax_sec, lalim, fm, alim_star, zmax &  ! in
12
13          , zcong, ale_bl, alp_bl, lalim_conv, wght_th &                    ! out
14          , zlcl, fraca0, w0, w_conv, therm_tke_max0, env_tke_max0 &   ! out
15          , n2, s2, strig, ale_bl_stat &                                   ! out
16          , therm_tke_max, env_tke_max &                           ! out
17          , alp_bl_det, alp_bl_fluct_m, alp_bl_fluct_tke &          ! out
18          , alp_bl_conv, alp_bl_stat &                             ! out
19          )
20
21    USE indice_sol_mod
22    USE lmdz_thermcell_main, ONLY: thermcell_tke_transport
23    USE lmdz_alpale
24    USE lmdz_yoethf
25    USE lmdz_fcttre, ONLY: foeew, foede, qsats, qsatl, dqsats, dqsatl, thermcep
26    USE lmdz_yomcst
27
28    IMPLICIT NONE
29
30    !=======================================================================
31
32    !   Auteurs: Catherine Rio
33    !   Modifications :
34    !   Nicolas Rochetin et Jean-Yves Grandpeix
35    !         pour la fermeture stochastique. 2012
36    !   Frédéric Hourdin :
37    !         netoyage informatique. 2022
38
39    !=======================================================================
40    !-----------------------------------------------------------------------
41    !   arguments:
42    !   ----------
43
44    !------Entrees
45    INTEGER, INTENT(IN) :: ngrid, nlay
46    REAL, INTENT(IN) :: ptimestep
47    REAL, INTENT(IN) :: pplay(ngrid, nlay), pplev(ngrid, nlay + 1)
48    INTEGER, INTENT(IN), DIMENSION(ngrid) :: lmax, lalim
49    REAL, INTENT(IN), DIMENSION(ngrid) :: zmax, zcong
50    REAL, INTENT(IN), DIMENSION(ngrid, nlay + 1) :: zw2
51    REAL, INTENT(IN), DIMENSION(ngrid, nlay + 1) :: fraca
52    REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: wth3
53    REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: rhobarz
54    REAL, INTENT(IN), DIMENSION(ngrid) :: wmax_sec
55    REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: entr0
56    REAL, INTENT(IN), DIMENSION(ngrid, nlay + 1) :: fm0, fm
57    REAL, INTENT(IN), DIMENSION(ngrid) :: pcon
58    REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: alim_star
59    REAL, INTENT(IN), DIMENSION(ngrid, nlay + 1, nbsrf) :: pbl_tke
60    REAL, INTENT(IN), DIMENSION(ngrid, nbsrf) :: pctsrf
61    REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: omega
62    REAL, INTENT(IN), DIMENSION(ngrid) :: airephy
63    !------Sorties
64    REAL, INTENT(OUT), DIMENSION(ngrid) :: ale_bl, alp_bl
65    REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: wght_th
66    INTEGER, INTENT(OUT), DIMENSION(ngrid) :: lalim_conv
67    REAL, INTENT(OUT), DIMENSION(ngrid) :: zlcl, fraca0, w0, w_conv
68    REAL, INTENT(OUT), DIMENSION(ngrid) :: therm_tke_max0, env_tke_max0, n2, s2, ale_bl_stat, strig
69    REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: therm_tke_max, env_tke_max
70    REAL, INTENT(OUT), DIMENSION(ngrid) :: alp_bl_det, alp_bl_fluct_m, alp_bl_fluct_tke
71    REAL, INTENT(OUT), DIMENSION(ngrid) :: alp_bl_conv, alp_bl_stat
72
73    !=============================================================================================
74    !------Local
75    !=============================================================================================
76
77    REAL susqr2pi, reuler
78    INTEGER ig, k, l
79    INTEGER nsrf
80    REAL rhobarz0(ngrid)                    ! Densité au LCL
81    LOGICAL ok_lcl(ngrid)                   ! Existence du LCL des thermiques
82    INTEGER klcl(ngrid)                     ! Niveau du LCL
83    REAL interp(ngrid)                      ! Coef d'interpolation pour le LCL
84    !--Triggering
85    REAL, parameter :: su_cst = 4e4              ! Surface unite: celle d'un updraft élémentaire
86    REAL, parameter :: hcoef = 1             ! Coefficient directeur pour le calcul de s2
87    REAL, parameter :: hmincoef = 0.3        ! Coefficient directeur pour l'ordonnée à l'origine pour le calcul de s2
88    REAL, parameter :: eps1 = 0.3            ! Fraction de surface occupée par la population 1 : eps1=n1*s1/(fraca0*Sd)
89    REAL, DIMENSION(ngrid) :: hmin         ! Ordonnée à l'origine pour le calcul de s2
90    REAL, DIMENSION(ngrid) :: zmax_moy     ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
91    REAL, parameter :: zmax_moy_coef = 0.33
92    REAL, DIMENSION(ngrid) :: depth        ! Epaisseur moyenne du cumulus
93    REAL, DIMENSION(ngrid) :: zcong_moy
94    REAL, DIMENSION(ngrid) :: w_max                 ! Vitesse max statistique
95    REAL, DIMENSION(ngrid) :: s_max(ngrid)
96    !--Closure
97    REAL, DIMENSION(ngrid, nlay) :: pbl_tke_max       ! Profil de TKE moyenne
98    REAL, DIMENSION(ngrid) :: pbl_tke_max0           ! TKE moyenne au LCL
99    REAL, DIMENSION(ngrid, nlay) :: w_ls              ! Vitesse verticale grande échelle (m/s)
100    REAL, parameter :: coef_m = 1.            ! On considère un rendement pour alp_bl_fluct_m
101    REAL, parameter :: coef_tke = 1.          ! On considère un rendement pour alp_bl_fluct_tke
102    REAL :: zdp
103    REAL, DIMENSION(ngrid) :: alp_int, dp_int
104    REAL, DIMENSION(ngrid) :: fm_tot
105
106    !------------------------------------------------------------
107    !  Initialize output arrays related to stochastic triggering
108    !------------------------------------------------------------
109    DO ig = 1, ngrid
110      zlcl(ig) = 0.
111      fraca0(ig) = 0.
112      w0(ig) = 0.
113      w_conv(ig) = 0.
114      therm_tke_max0(ig) = 0.
115      env_tke_max0(ig) = 0.
116      n2(ig) = 0.
117      s2(ig) = 0.
118      ale_bl_stat(ig) = 0.
119      strig(ig) = 0.
120      alp_bl_det(ig) = 0.
121      alp_bl_fluct_m(ig) = 0.
122      alp_bl_fluct_tke(ig) = 0.
123      alp_bl_conv(ig) = 0.
124      alp_bl_stat(ig) = 0.
125    ENDDO
126    DO l = 1, nlay
127      DO ig = 1, ngrid
128        therm_tke_max(ig, l) = 0.
129        env_tke_max(ig, l) = 0.
130      ENDDO
131    ENDDO
132
133    !------------Test sur le LCL des thermiques
134    do ig = 1, ngrid
135      ok_lcl(ig) = .FALSE.
136      IF ((pcon(ig) > pplay(ig, nlay - 1)) .AND. (pcon(ig) < pplay(ig, 1))) ok_lcl(ig) = .TRUE.
137    enddo
138
139    !------------Localisation des niveaux entourant le LCL et du coef d'interpolation
140    do l = 1, nlay - 1
141      do ig = 1, ngrid
142        IF (ok_lcl(ig)) THEN
143          !ATTENTION,zw2 calcule en pplev
144          !          if ((pplay(ig,l) .ge. pcon(ig)) .AND. (pplay(ig,l+1) .le. pcon(ig))) THEN
145          !          klcl(ig)=l
146          !          interp(ig)=(pcon(ig)-pplay(ig,klcl(ig)))/(pplay(ig,klcl(ig)+1)-pplay(ig,klcl(ig)))
147          !          endif
148          IF ((pplev(ig, l) >= pcon(ig)) .AND. (pplev(ig, l + 1) <= pcon(ig))) THEN
149            klcl(ig) = l
150            interp(ig) = (pcon(ig) - pplev(ig, klcl(ig))) / (pplev(ig, klcl(ig) + 1) - pplev(ig, klcl(ig)))
151          endif
152        endif
153      enddo
154    enddo
155
156    do ig = 1, ngrid
157      !CR:REHABILITATION ZMAX CONTINU
158      IF (ok_lcl(ig)) THEN
159        rhobarz0(ig) = rhobarz(ig, klcl(ig)) + (rhobarz(ig, klcl(ig) + 1) &
160                - rhobarz(ig, klcl(ig))) * interp(ig)
161        zlcl(ig) = (pplev(ig, 1) - pcon(ig)) / (rhobarz0(ig) * RG)
162        zlcl(ig) = min(zlcl(ig), zmax(ig))   ! Si zlcl > zmax alors on pose zlcl = zmax
163      else
164        rhobarz0(ig) = 0.
165        zlcl(ig) = zmax(ig)
166      endif
167    enddo
168    !!jyg fin
169
170    !------------Calcul des propriétés du thermique au LCL
171    IF ((iflag_trig_bl>=1) .OR. (iflag_clos_bl>=1)) THEN
172
173      !-----Initialisation de la TKE moyenne
174      do l = 1, nlay
175        do ig = 1, ngrid
176          pbl_tke_max(ig, l) = 0.
177        enddo
178      enddo
179
180      !-----Calcul de la TKE moyenne
181      do nsrf = 1, nbsrf
182        do l = 1, nlay
183          do ig = 1, ngrid
184            pbl_tke_max(ig, l) = pctsrf(ig, nsrf) * pbl_tke(ig, l, nsrf) + pbl_tke_max(ig, l)
185          enddo
186        enddo
187      enddo
188
189      !-----Initialisations des TKE dans et hors des thermiques
190      do l = 1, nlay
191        do ig = 1, ngrid
192          therm_tke_max(ig, l) = pbl_tke_max(ig, l)
193          env_tke_max(ig, l) = pbl_tke_max(ig, l)
194        enddo
195      enddo
196
197      !-----Calcul de la TKE transportée par les thermiques : therm_tke_max
198      CALL thermcell_tke_transport(ngrid, nlay, ptimestep, fm0, entr0, &  ! in
199              rg, pplev, therm_tke_max)                               ! out
200      !   print *,' thermcell_tke_transport -> '   !!jyg
201
202      !-----Calcul des profils verticaux de TKE hors thermiques : env_tke_max, et de la vitesse verticale grande échelle : W_ls
203      do l = 1, nlay
204        do ig = 1, ngrid
205          pbl_tke_max(ig, l) = fraca(ig, l) * therm_tke_max(ig, l) + (1. - fraca(ig, l)) * env_tke_max(ig, l)         !  Recalcul de TKE moyenne après transport de TKE_TH
206          env_tke_max(ig, l) = (pbl_tke_max(ig, l) - fraca(ig, l) * therm_tke_max(ig, l)) / (1. - fraca(ig, l))       !  Recalcul de TKE dans  l'environnement après transport de TKE_TH
207          w_ls(ig, l) = -1. * omega(ig, l) / (RG * rhobarz(ig, l))                                                !  Vitesse verticale de grande échelle
208        enddo
209      enddo
210      !    print *,' apres w_ls = '   !!jyg
211
212      do ig = 1, ngrid
213        IF (ok_lcl(ig)) THEN
214          fraca0(ig) = fraca(ig, klcl(ig)) + (fraca(ig, klcl(ig) + 1) &
215                  - fraca(ig, klcl(ig))) * interp(ig)
216          w0(ig) = zw2(ig, klcl(ig)) + (zw2(ig, klcl(ig) + 1) &
217                  - zw2(ig, klcl(ig))) * interp(ig)
218          w_conv(ig) = w_ls(ig, klcl(ig)) + (w_ls(ig, klcl(ig) + 1) &
219                  - w_ls(ig, klcl(ig))) * interp(ig)
220          therm_tke_max0(ig) = therm_tke_max(ig, klcl(ig)) &
221                  + (therm_tke_max(ig, klcl(ig) + 1) - therm_tke_max(ig, klcl(ig))) * interp(ig)
222          env_tke_max0(ig) = env_tke_max(ig, klcl(ig)) + (env_tke_max(ig, klcl(ig) + 1) &
223                  - env_tke_max(ig, klcl(ig))) * interp(ig)
224          pbl_tke_max0(ig) = pbl_tke_max(ig, klcl(ig)) + (pbl_tke_max(ig, klcl(ig) + 1) &
225                  - pbl_tke_max(ig, klcl(ig))) * interp(ig)
226          IF (therm_tke_max0(ig)>=20.) therm_tke_max0(ig) = 20.
227          IF (env_tke_max0(ig)>=20.) env_tke_max0(ig) = 20.
228          IF (pbl_tke_max0(ig)>=20.) pbl_tke_max0(ig) = 20.
229        else
230          fraca0(ig) = 0.
231          w0(ig) = 0.
232          !!jyg le 27/04/2012
233          !!     zlcl(ig)=0.
234          !!
235        endif
236      enddo
237
238    ENDIF ! IF ( (iflag_trig_bl.ge.1) .OR. (iflag_clos_bl.ge.1) )
239    !  print *,'ENDIF  ( (iflag_trig_bl.ge.1) .OR. (iflag_clos_bl.ge.1) ) '    !!jyg
240
241    !------------Triggering------------------
242    IF (iflag_trig_bl>=1) THEN
243
244      !-----Initialisations
245      depth(:) = 0.
246      n2(:) = 0.
247      s2(:) = 100. ! some low value, arbitrary
248      s_max(:) = 0.
249
250
251
252      !-----Epaisseur du nuage (depth) et détermination de la queue du spectre de panaches (n2,s2) et du panache le plus gros (s_max)
253      do ig = 1, ngrid
254        zmax_moy(ig) = zlcl(ig) + zmax_moy_coef * (zmax(ig) - zlcl(ig))
255        depth(ig) = zmax_moy(ig) - zlcl(ig)
256        hmin(ig) = hmincoef * zlcl(ig)
257        IF (depth(ig)>=10.) THEN
258          s2(ig) = (hcoef * depth(ig) + hmin(ig))**2
259          n2(ig) = (1. - eps1) * fraca0(ig) * airephy(ig) / s2(ig)
260          !!
261          !!jyg le 27/04/2012
262          !!       s_max(ig)=s2(ig)*log(n2(ig))
263          !!       if (n2(ig) .lt. 1) s_max(ig)=0.
264          s_max(ig) = s2(ig) * log(max(n2(ig), 1.))
265          !!fin jyg
266        else
267          n2(ig) = 0.
268          s_max(ig) = 0.
269        endif
270      enddo
271      !   print *,'avant Calcul de Wmax '    !!jyg
272
273      !CR: calcul de strig
274      IF (iflag_strig==0) THEN
275        strig(:) = s_trig
276      ELSE IF (iflag_strig==1) THEN
277        do ig = 1, ngrid
278          !         zcong_moy(ig)=zlcl(ig)+zmax_moy_coef*(zcong(ig)-zlcl(ig))
279          !         strig(ig)=(hcoef*(zcong_moy(ig)-zlcl(ig))+hmin(ig))**2
280          strig(ig) = (zcong(ig) - zlcl(ig))**2
281        enddo
282      ELSE IF (iflag_strig==2) THEN
283        do ig = 1, ngrid
284          IF (h_trig>zlcl(ig)) THEN
285            strig(ig) = (h_trig - zlcl(ig))**2
286          else
287            strig(ig) = s_trig
288          endif
289        enddo
290      endif
291
292      susqr2pi = su_cst * sqrt(2. * Rpi)
293      reuler = exp(1.)
294      do ig = 1, ngrid
295        IF ((depth(ig)>=10.) .AND. (s_max(ig)>susqr2pi * reuler)) THEN
296          w_max(ig) = w0(ig) * (1. + sqrt(2. * log(s_max(ig) / susqr2pi) - log(2. * log(s_max(ig) / susqr2pi))))
297          ale_bl_stat(ig) = 0.5 * w_max(ig)**2
298        else
299          w_max(ig) = 0.
300          ale_bl_stat(ig) = 0.
301        endif
302      enddo
303
304    ENDIF ! iflag_trig_bl
305    !  print *,'ENDIF  iflag_trig_bl'    !!jyg
306
307    !------------Closure------------------
308
309    IF (iflag_clos_bl>=2) THEN
310
311      !-----Calcul de ALP_BL_STAT
312      do ig = 1, ngrid
313        alp_bl_det(ig) = 0.5 * coef_m * rhobarz0(ig) * (w0(ig)**3) * fraca0(ig) * (1. - 2. * fraca0(ig)) / ((1. - fraca0(ig))**2)
314        alp_bl_fluct_m(ig) = 1.5 * rhobarz0(ig) * fraca0(ig) * (w_conv(ig) + coef_m * w0(ig)) * &
315                (w0(ig)**2)
316        alp_bl_fluct_tke(ig) = 3. * coef_m * rhobarz0(ig) * w0(ig) * fraca0(ig) * (therm_tke_max0(ig) - env_tke_max0(ig)) &
317                + 3. * rhobarz0(ig) * w_conv(ig) * pbl_tke_max0(ig)
318        IF (iflag_clos_bl>=2) THEN
319          alp_bl_conv(ig) = 1.5 * coef_m * rhobarz0(ig) * fraca0(ig) * (fraca0(ig) / (1. - fraca0(ig))) * w_conv(ig) * &
320                  (w0(ig)**2)
321        else
322          alp_bl_conv(ig) = 0.
323        endif
324        alp_bl_stat(ig) = alp_bl_det(ig) + alp_bl_fluct_m(ig) + alp_bl_fluct_tke(ig) + alp_bl_conv(ig)
325      enddo
326
327      !-----Sécurité ALP infinie
328      do ig = 1, ngrid
329        IF (fraca0(ig)>0.98) alp_bl_stat(ig) = 2.
330      enddo
331
332    ENDIF ! (iflag_clos_bl.ge.2)
333
334    !!! fin nrlmd le 10/04/2012
335
336    !      PRINT*,'avant calcul ale et alp'
337    !calcul de ALE et ALP pour la convection
338    alp_bl(:) = 0.
339    ale_bl(:) = 0.
340    !          PRINT*,'ALE,ALP ,l,zw2(ig,l),ale_bl(ig),alp_bl(ig)'
341    do l = 1, nlay
342      do ig = 1, ngrid
343        alp_bl(ig) = max(alp_bl(ig), 0.5 * rhobarz(ig, l) * wth3(ig, l))
344        ale_bl(ig) = max(ale_bl(ig), 0.5 * zw2(ig, l)**2)
345        !          PRINT*,'ALE,ALP',l,zw2(ig,l),ale_bl(ig),alp_bl(ig)
346      enddo
347    enddo
348
349    ! ale sec (max de wmax/2 sous la zone d'inhibition) dans
350    ! le cas iflag_trig_bl=3
351    IF (iflag_trig_bl==3) ale_bl(:) = 0.5 * wmax_sec(:)**2
352
353    !test:calcul de la ponderation des couches pour KE
354    !initialisations
355
356    fm_tot(:) = 0.
357    wght_th(:, :) = 1.
358    lalim_conv(:) = lalim(:)
359
360    do k = 1, nlay
361      do ig = 1, ngrid
362        IF (k<=lalim_conv(ig)) fm_tot(ig) = fm_tot(ig) + fm(ig, k)
363      enddo
364    enddo
365
366    ! assez bizarre car, si on est dans la couche d'alim et que alim_star et
367    ! plus petit que 1.e-10, on prend wght_th=1.
368    do k = 1, nlay
369      do ig = 1, ngrid
370        IF (k<=lalim_conv(ig).AND.alim_star(ig, k)>1.e-10) THEN
371          wght_th(ig, k) = alim_star(ig, k)
372        endif
373      enddo
374    enddo
375
376    !      PRINT*,'apres wght_th'
377    !test pour prolonger la convection
378    do ig = 1, ngrid
379      !v1d  if ((alim_star(ig,1).lt.1.e-10).AND.(therm)) THEN
380      IF ((alim_star(ig, 1)<1.e-10)) THEN
381        lalim_conv(ig) = 1
382        wght_th(ig, 1) = 1.
383        !      PRINT*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
384      endif
385    enddo
386
387    !------------------------------------------------------------------------
388    ! Modif CR/FH 20110310 : alp integree sur la verticale.
389    ! Integrale verticale de ALP.
390    ! wth3 etant aux niveaux inter-couches, on utilise d play comme masse des
391    ! couches
392    !------------------------------------------------------------------------
393
394    alp_int(:) = 0.
395    dp_int(:) = 0.
396    do l = 2, nlay
397      do ig = 1, ngrid
398        IF(l<=lmax(ig)) THEN
399          zdp = pplay(ig, l - 1) - pplay(ig, l)
400          alp_int(ig) = alp_int(ig) + 0.5 * rhobarz(ig, l) * wth3(ig, l) * zdp
401          dp_int(ig) = dp_int(ig) + zdp
402        endif
403      enddo
404    enddo
405
406    IF (iflag_coupl>=3 .AND. iflag_coupl<=5) THEN
407      do ig = 1, ngrid
408        !valeur integree de alp_bl * 0.5:
409        IF (dp_int(ig)>0.) THEN
410          alp_bl(ig) = alp_int(ig) / dp_int(ig)
411        endif
412      enddo!
413    endif
414
415
416    ! Facteur multiplicatif sur alp_bl
417    alp_bl(:) = alp_bl_k * alp_bl(:)
418
419    !------------------------------------------------------------------------
420
421    RETURN
422  END
423END MODULE lmdz_thermcell_alp
Note: See TracBrowser for help on using the repository browser.