source: LMDZ6/trunk/libf/phylmd/thermcell_alp.F90 @ 4283

Last change on this file since 4283 was 4089, checked in by fhourdin, 2 years ago

Reecriture des thermiques

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