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

Last change on this file was 5160, checked in by abarral, 3 months ago

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