source: LMDZ5/trunk/libf/phylmd/thermcell_alp.F90 @ 2384

Last change on this file since 2384 was 2384, checked in by fhourdin, 9 years ago

Séparation de la partie de calcul de ALP dans thermcell_main
Spliting thermcell_main.F90 in 2, isolating the ALP computation.

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