source: trunk/LMDZ.GENERIC/libf/phystd/thermcell_alp.F90 @ 2099

Last change on this file since 2099 was 2060, checked in by aboissinot, 6 years ago

Add the thermal plume model (cf. Rio et al. 2010) extended to gas giant.
Specific parameters are set in thermcell_mod.

File size: 15.9 KB
Line 
1!
2!
3!
4      SUBROUTINE thermcell_alp(ngrid,nlay,ptimestep  &
5     &                  ,pplay,pplev  &
6     &                  ,fm0,entr0,lmax  &
7     &                  ,ale_bl,alp_bl,lalim_conv,wght_th &
8     &                  ,zw2,fraca &
9!!! ncessaire en plus
10     &                  ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &
11!!! nrlmd le 10/04/2012
12     &                  ,pbl_tke,pctsrf,omega,airephy &
13     &                  ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
14     &                  ,n2,s2,ale_bl_stat &
15     &                  ,therm_tke_max,env_tke_max &
16     &                  ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
17     &                  ,alp_bl_conv,alp_bl_stat )
18!!! fin nrlmd le 10/04/2012
19
20      USE dimphy
21! AB : why use dimphy to get klon, klev while subroutine arguments ngrid, nlay are the same?
22      USE thermcell_mod
23     
24      IMPLICIT NONE
25
26!=======================================================================
27!   Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
28!   Version du 09.02.07
29!   Calcul du transport vertical dans la couche limite en presence
30!   de "thermiques" explicitement representes avec processus nuageux
31!
32!   Reecriture a partir d'un listing papier a Habas, le 14/02/00
33!
34!   le thermique est suppose homogene et dissipe par melange avec
35!   son environnement. la longueur l_mix controle l'efficacite du
36!   melange
37!
38!   Le calcul du transport des differentes especes se fait en prenant
39!   en compte:
40!     1. un flux de masse montant
41!     2. un flux de masse descendant
42!     3. un entrainement
43!     4. un detrainement
44!
45! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
46!    Introduction of an implicit computation of vertical advection in
47!    the environment of thermal plumes in thermcell_dq
48!    impl =     0 : explicit, 1 : implicit, -1 : old version
49!    controled by iflag_thermals =
50!       15, 16 run with impl=-1 : numerical convergence with NPv3
51!       17, 18 run with impl=1  : more stable
52!    15 and 17 correspond to the activation of the stratocumulus "bidouille"
53!
54!=======================================================================
55
56!-----------------------------------------------------------------------
57!   declarations:
58!   -------------
59
60!   arguments:
61!   ----------
62
63!IM 140508
64
65      INTEGER ngrid,nlay
66      real ptimestep
67      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
68
69!   local:
70!   ------
71
72
73      REAL susqr2pi, reuler
74
75      INTEGER ig,k,l
76      INTEGER lmax(klon),lalim(klon)
77      real zmax(klon),zw2(klon,klev+1)
78
79!on garde le zmax du pas de temps precedent
80
81
82      real fraca(klon,klev+1)
83      real wth3(klon,klev)
84! FH probleme de dimensionnement avec l'allocation dynamique
85!     common/comtherm/thetath2,wth2
86      real rhobarz(klon,klev)
87
88      real wmax_sec(klon)
89      real fm0(klon,klev+1),entr0(klon,klev)
90      real fm(klon,klev+1)
91
92!niveau de condensation
93      real pcon(klon)
94
95      real alim_star(klon,klev)
96
97!!! nrlmd le 10/04/2012
98
99!------Entrées
100      real pbl_tke(klon,klev+1,nbsrf)
101      real pctsrf(klon,nbsrf)
102      real omega(klon,klev)
103      real airephy(klon)
104!------Sorties
105      real zlcl(klon),fraca0(klon),w0(klon),w_conv(klon)
106      real therm_tke_max0(klon),env_tke_max0(klon)
107      real n2(klon),s2(klon)
108      real ale_bl_stat(klon)
109      real therm_tke_max(klon,klev),env_tke_max(klon,klev)
110      real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon)
111!------Local
112      integer nsrf
113      real rhobarz0(klon)                    ! Densité au LCL
114      logical ok_lcl(klon)                   ! Existence du LCL des thermiques
115      integer klcl(klon)                     ! Niveau du LCL
116      real interp(klon)                      ! Coef d'interpolation pour le LCL
117!--Triggering
118      real Su                                ! Surface unité: celle d'un updraft élémentaire
119      parameter(Su=4e4)
120      real hcoef                             ! Coefficient directeur pour le calcul de s2
121      parameter(hcoef=1)
122      real hmincoef                          ! Coefficient directeur pour l'ordonnée à l'origine pour le calcul de s2
123      parameter(hmincoef=0.3)
124      real eps1                              ! Fraction de surface occupée par la population 1 : eps1=n1*s1/(fraca0*Sd)
125      parameter(eps1=0.3)
126      real hmin(ngrid)                       ! Ordonnée à l'origine pour le calcul de s2
127      real zmax_moy(ngrid)                   ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
128      real zmax_moy_coef
129      parameter(zmax_moy_coef=0.33)
130      real depth(klon)                       ! Epaisseur moyenne du cumulus
131      real w_max(klon)                       ! Vitesse max statistique
132      real s_max(klon)
133!--Closure
134      real pbl_tke_max(klon,klev)            ! Profil de TKE moyenne
135      real pbl_tke_max0(klon)                ! TKE moyenne au LCL
136      real w_ls(klon,klev)                   ! Vitesse verticale grande échelle (m/s)
137      real coef_m                            ! On considère un rendement pour alp_bl_fluct_m
138      parameter(coef_m=1.)
139      real coef_tke                          ! On considère un rendement pour alp_bl_fluct_tke
140      parameter(coef_tke=1.)
141
142!!! fin nrlmd le 10/04/2012
143
144!
145      !nouvelles variables pour la convection
146      real ale_bl(klon)
147      real alp_bl(klon)
148      real alp_int(klon),dp_int(klon),zdp
149      real fm_tot(klon)
150      real wght_th(klon,klev)
151      integer lalim_conv(klon)
152!v1d     logical therm
153!v1d     save therm
154     
155     
156!------------------------------------------------------------
157!  Initialize output arrays related to stochastic triggering
158!------------------------------------------------------------
159     
160      DO ig = 1,klon
161         zlcl(ig) = 0.
162         fraca0(ig) = 0.
163         w0(ig) = 0.
164         w_conv(ig) = 0.
165         therm_tke_max0(ig) = 0.
166         env_tke_max0(ig) = 0.
167         n2(ig) = 0.
168         s2(ig) = 0.
169         ale_bl_stat(ig) = 0.
170         alp_bl_det(ig) = 0.
171         alp_bl_fluct_m(ig) = 0.
172         alp_bl_fluct_tke(ig) = 0.
173         alp_bl_conv(ig) = 0.
174         alp_bl_stat(ig) = 0.
175      ENDDO
176     
177      DO l = 1,klev
178         DO ig = 1,klon
179            therm_tke_max(ig,l) = 0.
180            env_tke_max(ig,l) = 0.
181         ENDDO
182      ENDDO
183     
184!------------------------------------------------------------
185!------------Test sur le LCL des thermiques
186      do ig=1,ngrid
187      ok_lcl(ig)=.false.
188      if ( (pcon(ig) .gt. pplay(ig,klev-1)) .and. (pcon(ig) .lt. pplay(ig,1)) ) ok_lcl(ig)=.true.
189   enddo
190
191!------------Localisation des niveaux entourant le LCL et du coef d'interpolation
192    do l=1,nlay-1
193      do ig=1,ngrid
194        if (ok_lcl(ig)) then
195!ATTENTION,zw2 calcule en pplev
196!          if ((pplay(ig,l) .ge. pcon(ig)) .and. (pplay(ig,l+1) .le. pcon(ig))) then
197!          klcl(ig)=l
198!          interp(ig)=(pcon(ig)-pplay(ig,klcl(ig)))/(pplay(ig,klcl(ig)+1)-pplay(ig,klcl(ig)))
199!          endif
200          if ((pplev(ig,l) .ge. pcon(ig)) .and. (pplev(ig,l+1) .le. pcon(ig))) then
201          klcl(ig)=l
202          interp(ig)=(pcon(ig)-pplev(ig,klcl(ig)))/(pplev(ig,klcl(ig)+1)-pplev(ig,klcl(ig)))
203          endif
204        endif
205      enddo
206    enddo
207
208!------------Hauteur des thermiques
209!!jyg le 27/04/2012
210!!    do ig =1,ngrid
211!!    rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
212!!    &               -rhobarz(ig,klcl(ig)))*interp(ig)
213!!    zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
214!!      if ( (.not.ok_lcl(ig)) .or. (zlcl(ig).gt.zmax(ig)) ) zlcl(ig)=zmax(ig) ! Si zclc > zmax alors on pose zlcl = zmax
215!!    enddo
216    do ig =1,ngrid
217!CR:REHABILITATION ZMAX CONTINU
218      if (ok_lcl(ig)) then
219         rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
220         &               -rhobarz(ig,klcl(ig)))*interp(ig)
221         zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
222         zlcl(ig)=min(zlcl(ig),zmax(ig))   ! Si zlcl > zmax alors on pose zlcl = zmax
223      else
224         rhobarz0(ig)=0.
225         zlcl(ig)=zmax(ig)
226      endif
227   enddo
228!!jyg fin
229
230!------------Calcul des propriétés du thermique au LCL
231  IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) THEN
232
233!-----Initialisation de la TKE moyenne
234   do l=1,nlay
235      do ig=1,ngrid
236         pbl_tke_max(ig,l)=0.
237      enddo
238   enddo
239
240!-----Calcul de la TKE moyenne
241   do nsrf=1,nbsrf
242      do l=1,nlay
243         do ig=1,ngrid
244            pbl_tke_max(ig,l)=pctsrf(ig,nsrf)*pbl_tke(ig,l,nsrf)+pbl_tke_max(ig,l)
245         enddo
246      enddo
247   enddo
248
249!-----Initialisations des TKE dans et hors des thermiques
250   do l=1,nlay
251      do ig=1,ngrid
252         therm_tke_max(ig,l)=pbl_tke_max(ig,l)
253         env_tke_max(ig,l)=pbl_tke_max(ig,l)
254      enddo
255   enddo
256
257!-----Calcul de la TKE transportée par les thermiques : therm_tke_max
258   call thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0,  &
259   &           rg,pplev,therm_tke_max)
260   
261!   print *,' thermcell_tke_transport -> '   !!jyg
262
263!-----Calcul des profils verticaux de TKE hors thermiques : env_tke_max, et de la vitesse verticale grande échelle : W_ls
264   do l=1,nlay
265      do ig=1,ngrid
266         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
267         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
268         w_ls(ig,l)=-1.*omega(ig,l)/(RG*rhobarz(ig,l))                                                !  Vitesse verticale de grande échelle
269      enddo
270   enddo
271   
272!   print *,' apres w_ls = '   !!jyg
273   
274   do ig=1,ngrid
275      if (ok_lcl(ig)) then
276         fraca0(ig)=fraca(ig,klcl(ig))+(fraca(ig,klcl(ig)+1) &
277         &             -fraca(ig,klcl(ig)))*interp(ig)
278         w0(ig)=zw2(ig,klcl(ig))+(zw2(ig,klcl(ig)+1) &
279         &         -zw2(ig,klcl(ig)))*interp(ig)
280         w_conv(ig)=w_ls(ig,klcl(ig))+(w_ls(ig,klcl(ig)+1) &
281         &             -w_ls(ig,klcl(ig)))*interp(ig)
282         therm_tke_max0(ig)=therm_tke_max(ig,klcl(ig)) &
283         &                     +(therm_tke_max(ig,klcl(ig)+1)-therm_tke_max(ig,klcl(ig)))*interp(ig)
284         env_tke_max0(ig)=env_tke_max(ig,klcl(ig))+(env_tke_max(ig,klcl(ig)+1) &
285         &                   -env_tke_max(ig,klcl(ig)))*interp(ig)
286         pbl_tke_max0(ig)=pbl_tke_max(ig,klcl(ig))+(pbl_tke_max(ig,klcl(ig)+1) &
287         &                   -pbl_tke_max(ig,klcl(ig)))*interp(ig)
288         
289         if (therm_tke_max0(ig).ge.20.) therm_tke_max0(ig)=20.
290         
291         if (env_tke_max0(ig).ge.20.) env_tke_max0(ig)=20.
292         
293         if (pbl_tke_max0(ig).ge.20.) pbl_tke_max0(ig)=20.
294      else
295         fraca0(ig)=0.
296         w0(ig)=0.
297! jyg le 27/04/2012
298!         zlcl(ig)=0.
299      endif
300   enddo
301   
302   ENDIF ! IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) )
303
304!   print *,'ENDIF  ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) '    !!jyg
305
306!------------Triggering------------------
307  IF (iflag_trig_bl.ge.1) THEN
308
309!-----Initialisations
310   depth(:)=0.
311   n2(:)=0.
312   s2(:)=100. ! some low value, arbitrary
313   s_max(:)=0.
314
315!-----Epaisseur du nuage (depth) et détermination de la queue du spectre de panaches (n2,s2) et du panache le plus gros (s_max)
316   do ig=1,ngrid
317      zmax_moy(ig)=zlcl(ig)+zmax_moy_coef*(zmax(ig)-zlcl(ig))
318      depth(ig)=zmax_moy(ig)-zlcl(ig)
319      hmin(ig)=hmincoef*zlcl(ig)
320     
321      if (depth(ig).ge.10.) then
322         s2(ig)=(hcoef*depth(ig)+hmin(ig))**2
323         n2(ig)=(1.-eps1)*fraca0(ig)*airephy(ig)/s2(ig)
324         
325!!jyg le 27/04/2012
326!!        s_max(ig)=s2(ig)*log(n2(ig))
327!!        if (n2(ig) .lt. 1) s_max(ig)=0.
328          s_max(ig)=s2(ig)*log(max(n2(ig),1.))
329!!fin jyg
330      else
331         n2(ig)=0.
332         s_max(ig)=0.
333      endif
334   enddo
335   
336!   print *,'avant Calcul de Wmax '    !!jyg
337
338!-----Calcul de Wmax et ALE_BL_STAT associée
339!!jyg le 30/04/2012
340!!   do ig=1,ngrid
341!!     if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.1.) ) then
342!!     w_max(ig)=w0(ig)*(1.+sqrt(2.*log(s_max(ig)/su)-log(2.*3.14)-log(2.*log(s_max(ig)/su)-log(2.*3.14))))
343!!     ale_bl_stat(ig)=0.5*w_max(ig)**2
344!!     else
345!!     w_max(ig)=0.
346!!     ale_bl_stat(ig)=0.
347!!     endif
348!!   enddo
349   
350   susqr2pi=su*sqrt(2.*Rpi)
351   reuler=exp(1.)
352   
353   do ig=1,ngrid
354      if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.susqr2pi*reuler) ) then
355         w_max(ig)=w0(ig)*(1.+sqrt(2.*log(s_max(ig)/susqr2pi)-log(2.*log(s_max(ig)/susqr2pi))))
356         ale_bl_stat(ig)=0.5*w_max(ig)**2
357      else
358         w_max(ig)=0.
359         ale_bl_stat(ig)=0.
360      endif
361   enddo
362
363   ENDIF ! iflag_trig_bl
364   
365!   print *,'ENDIF  iflag_trig_bl'    !!jyg
366
367!------------Closure------------------
368
369   IF (iflag_clos_bl.ge.2) THEN
370     
371!-----Calcul de ALP_BL_STAT
372      do ig=1,ngrid
373         alp_bl_det(ig)=0.5*coef_m*rhobarz0(ig)*(w0(ig)**3)*fraca0(ig)*(1.-2.*fraca0(ig))/((1.-fraca0(ig))**2)
374         alp_bl_fluct_m(ig)=1.5*rhobarz0(ig)*fraca0(ig)*(w_conv(ig)+coef_m*w0(ig))* &
375         &                   (w0(ig)**2)
376         alp_bl_fluct_tke(ig)=3.*coef_m*rhobarz0(ig)*w0(ig)*fraca0(ig)*(therm_tke_max0(ig)-env_tke_max0(ig)) &
377         &                    +3.*rhobarz0(ig)*w_conv(ig)*pbl_tke_max0(ig)
378         
379         if (iflag_clos_bl.ge.2) then
380            alp_bl_conv(ig)=1.5*coef_m*rhobarz0(ig)*fraca0(ig)*(fraca0(ig)/(1.-fraca0(ig)))*w_conv(ig)* &
381            &                   (w0(ig)**2)
382         else
383            alp_bl_conv(ig)=0.
384         endif
385         
386         alp_bl_stat(ig)=alp_bl_det(ig)+alp_bl_fluct_m(ig)+alp_bl_fluct_tke(ig)+alp_bl_conv(ig)
387      enddo
388     
389!-----Sécurité ALP infinie
390      do ig=1,ngrid
391         if (fraca0(ig).gt.0.98) alp_bl_stat(ig)=2.
392      enddo
393     
394  ENDIF ! (iflag_clos_bl.ge.2)
395
396!!! fin nrlmd le 10/04/2012
397
398!   print*,'avant calcul ale et alp'
399
400!calcul de ALE et ALP pour la convection
401   alp_bl(:)=0.
402   ale_bl(:)=0.
403     
404!   print*,'ALE,ALP ,l,zw2(ig,l),ale_bl(ig),alp_bl(ig)'
405
406   do l=1,nlay
407      do ig=1,ngrid
408         alp_bl(ig)=max(alp_bl(ig),0.5*rhobarz(ig,l)*wth3(ig,l) )
409         ale_bl(ig)=max(ale_bl(ig),0.5*zw2(ig,l)**2)
410!         print*,'ALE,ALP',l,zw2(ig,l),ale_bl(ig),alp_bl(ig)
411      enddo
412   enddo
413
414!   ale sec (max de wmax/2 sous la zone d'inhibition) dans
415!   le cas iflag_trig_bl=3
416   
417   IF (iflag_trig_bl==3) then
418      ale_bl(:)=0.5*wmax_sec(:)**2
419   ENDIF
420
421!test:calcul de la ponderation des couches pour KE
422!initialisations
423
424      fm_tot(:)=0.
425      wght_th(:,:)=1.
426      lalim_conv(:)=lalim(:)
427
428      do k=1,klev
429         do ig=1,ngrid
430            if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k)
431         enddo
432      enddo
433
434! assez bizarre car, si on est dans la couche d'alim et que alim_star et
435! plus petit que 1.e-10, on prend wght_th=1.
436      do k=1,klev
437         do ig=1,ngrid
438            if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then
439               wght_th(ig,k)=alim_star(ig,k)
440            endif
441         enddo
442      enddo
443
444!      print*,'apres wght_th'
445!test pour prolonger la convection
446      do ig=1,ngrid
447!v1d  if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
448      if ((alim_star(ig,1).lt.1.e-10)) then
449      lalim_conv(ig)=1
450      wght_th(ig,1)=1.
451!      print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
452      endif
453      enddo
454
455!------------------------------------------------------------------------
456! Modif CR/FH 20110310 : alp integree sur la verticale.
457! Integrale verticale de ALP.
458! wth3 etant aux niveaux inter-couches, on utilise d play comme masse des
459! couches
460!------------------------------------------------------------------------
461
462      alp_int(:)=0.
463      dp_int(:)=0.
464      do l=2,nlay
465        do ig=1,ngrid
466           if(l.LE.lmax(ig)) THEN
467           zdp=pplay(ig,l-1)-pplay(ig,l)
468           alp_int(ig)=alp_int(ig)+0.5*rhobarz(ig,l)*wth3(ig,l)*zdp
469           dp_int(ig)=dp_int(ig)+zdp
470           endif
471        enddo
472      enddo
473
474      if (iflag_coupl>=3 .and. iflag_coupl<=5) then
475         do ig=1,ngrid
476            ! valeur integree de alp_bl * 0.5:
477            if (dp_int(ig)>0.) then
478               alp_bl(ig)=alp_int(ig)/dp_int(ig)
479            endif
480         enddo
481      endif
482
483
484! Facteur multiplicatif sur alp_bl
485      alp_bl(:)=alp_bl_k*alp_bl(:)
486
487!------------------------------------------------------------------------
488
489
490
491      return
492      end
Note: See TracBrowser for help on using the repository browser.