source: LMDZ6/trunk/libf/phylmd/lmdz_thermcell_alp.f90 @ 5274

Last change on this file since 5274 was 5274, checked in by abarral, 3 hours ago

Replace yomcst.h by existing module

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