source: LMDZ6/branches/Ocean_skin/libf/phylmd/thermcell_alp.F90 @ 3627

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

Poursuite du travail sur les thermiques.
Extraction d'une routine calculant l'alimentation laterale à la base.
+ une correction de Jean-Yves sur des initialisations maquantes dans
thermcell_main (mais qui sont maintenant dans thermcell_alp).

File size: 15.4 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!------------------------------------------------------------
159!  Initialize output arrays related to stochastic triggering
160!------------------------------------------------------------
161  DO ig = 1,klon
162     zlcl(ig) = 0.
163     fraca0(ig) = 0.
164     w0(ig) = 0.
165     w_conv(ig) = 0.
166     therm_tke_max0(ig) = 0.
167     env_tke_max0(ig) = 0.
168     n2(ig) = 0.
169     s2(ig) = 0.
170     ale_bl_stat(ig) = 0.
171     alp_bl_det(ig) = 0.
172     alp_bl_fluct_m(ig) = 0.
173     alp_bl_fluct_tke(ig) = 0.
174     alp_bl_conv(ig) = 0.
175     alp_bl_stat(ig) = 0.
176  ENDDO
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
186!------------Test sur le LCL des thermiques
187    do ig=1,ngrid
188      ok_lcl(ig)=.false.
189      if ( (pcon(ig) .gt. pplay(ig,klev-1)) .and. (pcon(ig) .lt. pplay(ig,1)) ) ok_lcl(ig)=.true.
190    enddo
191
192!------------Localisation des niveaux entourant le LCL et du coef d'interpolation
193    do l=1,nlay-1
194      do ig=1,ngrid
195        if (ok_lcl(ig)) then
196!ATTENTION,zw2 calcule en pplev
197!          if ((pplay(ig,l) .ge. pcon(ig)) .and. (pplay(ig,l+1) .le. pcon(ig))) then
198!          klcl(ig)=l
199!          interp(ig)=(pcon(ig)-pplay(ig,klcl(ig)))/(pplay(ig,klcl(ig)+1)-pplay(ig,klcl(ig)))
200!          endif
201          if ((pplev(ig,l) .ge. pcon(ig)) .and. (pplev(ig,l+1) .le. pcon(ig))) then
202          klcl(ig)=l
203          interp(ig)=(pcon(ig)-pplev(ig,klcl(ig)))/(pplev(ig,klcl(ig)+1)-pplev(ig,klcl(ig)))
204          endif
205        endif
206      enddo
207    enddo
208
209!------------Hauteur des thermiques
210!!jyg le 27/04/2012
211!!    do ig =1,ngrid
212!!    rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
213!! &               -rhobarz(ig,klcl(ig)))*interp(ig)
214!!    zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
215!!      if ( (.not.ok_lcl(ig)) .or. (zlcl(ig).gt.zmax(ig)) ) zlcl(ig)=zmax(ig) ! Si zclc > zmax alors on pose zlcl = zmax
216!!    enddo
217    do ig =1,ngrid
218!CR:REHABILITATION ZMAX CONTINU
219     if (ok_lcl(ig)) then
220      rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
221 &               -rhobarz(ig,klcl(ig)))*interp(ig)
222      zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
223      zlcl(ig)=min(zlcl(ig),zmax(ig))   ! Si zlcl > zmax alors on pose zlcl = zmax
224     else
225      rhobarz0(ig)=0.
226      zlcl(ig)=zmax(ig)
227     endif
228    enddo
229!!jyg fin
230
231!------------Calcul des propriétés du thermique au LCL
232  IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) THEN
233
234  !-----Initialisation de la TKE moyenne
235   do l=1,nlay
236    do ig=1,ngrid
237     pbl_tke_max(ig,l)=0.
238    enddo
239   enddo
240
241!-----Calcul de la TKE moyenne
242   do nsrf=1,nbsrf
243    do l=1,nlay
244     do ig=1,ngrid
245     pbl_tke_max(ig,l)=pctsrf(ig,nsrf)*pbl_tke(ig,l,nsrf)+pbl_tke_max(ig,l)
246     enddo
247    enddo
248   enddo
249
250!-----Initialisations des TKE dans et hors des thermiques
251   do l=1,nlay
252    do ig=1,ngrid
253    therm_tke_max(ig,l)=pbl_tke_max(ig,l)
254    env_tke_max(ig,l)=pbl_tke_max(ig,l)
255    enddo
256   enddo
257
258!-----Calcul de la TKE transportée par les thermiques : therm_tke_max
259   call thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0,  &
260  &           rg,pplev,therm_tke_max)
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!    print *,' apres w_ls = '   !!jyg
272
273  do ig=1,ngrid
274   if (ok_lcl(ig)) then
275     fraca0(ig)=fraca(ig,klcl(ig))+(fraca(ig,klcl(ig)+1) &
276 &             -fraca(ig,klcl(ig)))*interp(ig)
277     w0(ig)=zw2(ig,klcl(ig))+(zw2(ig,klcl(ig)+1) &
278 &         -zw2(ig,klcl(ig)))*interp(ig)
279     w_conv(ig)=w_ls(ig,klcl(ig))+(w_ls(ig,klcl(ig)+1) &
280 &             -w_ls(ig,klcl(ig)))*interp(ig)
281     therm_tke_max0(ig)=therm_tke_max(ig,klcl(ig)) &
282 &                     +(therm_tke_max(ig,klcl(ig)+1)-therm_tke_max(ig,klcl(ig)))*interp(ig)
283     env_tke_max0(ig)=env_tke_max(ig,klcl(ig))+(env_tke_max(ig,klcl(ig)+1) &
284 &                   -env_tke_max(ig,klcl(ig)))*interp(ig)
285     pbl_tke_max0(ig)=pbl_tke_max(ig,klcl(ig))+(pbl_tke_max(ig,klcl(ig)+1) &
286 &                   -pbl_tke_max(ig,klcl(ig)))*interp(ig)
287     if (therm_tke_max0(ig).ge.20.) therm_tke_max0(ig)=20.
288     if (env_tke_max0(ig).ge.20.) env_tke_max0(ig)=20.
289     if (pbl_tke_max0(ig).ge.20.) pbl_tke_max0(ig)=20.
290   else
291     fraca0(ig)=0.
292     w0(ig)=0.
293!!jyg le 27/04/2012
294!!     zlcl(ig)=0.
295!!
296   endif
297  enddo
298
299  ENDIF ! IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) )
300!  print *,'ENDIF  ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) '    !!jyg
301
302!------------Triggering------------------
303  IF (iflag_trig_bl.ge.1) THEN
304
305!-----Initialisations
306   depth(:)=0.
307   n2(:)=0.
308   s2(:)=100. ! some low value, arbitrary
309   s_max(:)=0.
310
311!-----Epaisseur du nuage (depth) et détermination de la queue du spectre de panaches (n2,s2) et du panache le plus gros (s_max)
312   do ig=1,ngrid
313     zmax_moy(ig)=zlcl(ig)+zmax_moy_coef*(zmax(ig)-zlcl(ig))
314     depth(ig)=zmax_moy(ig)-zlcl(ig)
315     hmin(ig)=hmincoef*zlcl(ig)
316     if (depth(ig).ge.10.) then
317       s2(ig)=(hcoef*depth(ig)+hmin(ig))**2
318       n2(ig)=(1.-eps1)*fraca0(ig)*airephy(ig)/s2(ig)
319!!
320!!jyg le 27/04/2012
321!!       s_max(ig)=s2(ig)*log(n2(ig))
322!!       if (n2(ig) .lt. 1) s_max(ig)=0.
323       s_max(ig)=s2(ig)*log(max(n2(ig),1.))
324!!fin jyg
325     else
326       n2(ig)=0.
327       s_max(ig)=0.
328     endif
329   enddo
330!   print *,'avant Calcul de Wmax '    !!jyg
331
332!-----Calcul de Wmax et ALE_BL_STAT associée
333!!jyg le 30/04/2012
334!!   do ig=1,ngrid
335!!     if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.1.) ) then
336!!     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))))
337!!     ale_bl_stat(ig)=0.5*w_max(ig)**2
338!!     else
339!!     w_max(ig)=0.
340!!     ale_bl_stat(ig)=0.
341!!     endif
342!!   enddo
343   susqr2pi=su*sqrt(2.*Rpi)
344   reuler=exp(1.)
345   do ig=1,ngrid
346     if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.susqr2pi*reuler) ) then
347      w_max(ig)=w0(ig)*(1.+sqrt(2.*log(s_max(ig)/susqr2pi)-log(2.*log(s_max(ig)/susqr2pi))))
348      ale_bl_stat(ig)=0.5*w_max(ig)**2
349     else
350      w_max(ig)=0.
351      ale_bl_stat(ig)=0.
352     endif
353   enddo
354
355  ENDIF ! iflag_trig_bl
356!  print *,'ENDIF  iflag_trig_bl'    !!jyg
357
358!------------Closure------------------
359
360  IF (iflag_clos_bl.ge.2) THEN
361
362!-----Calcul de ALP_BL_STAT
363  do ig=1,ngrid
364  alp_bl_det(ig)=0.5*coef_m*rhobarz0(ig)*(w0(ig)**3)*fraca0(ig)*(1.-2.*fraca0(ig))/((1.-fraca0(ig))**2)
365  alp_bl_fluct_m(ig)=1.5*rhobarz0(ig)*fraca0(ig)*(w_conv(ig)+coef_m*w0(ig))* &
366 &                   (w0(ig)**2)
367  alp_bl_fluct_tke(ig)=3.*coef_m*rhobarz0(ig)*w0(ig)*fraca0(ig)*(therm_tke_max0(ig)-env_tke_max0(ig)) &
368 &                    +3.*rhobarz0(ig)*w_conv(ig)*pbl_tke_max0(ig)
369    if (iflag_clos_bl.ge.2) then
370    alp_bl_conv(ig)=1.5*coef_m*rhobarz0(ig)*fraca0(ig)*(fraca0(ig)/(1.-fraca0(ig)))*w_conv(ig)* &
371 &                   (w0(ig)**2)
372    else
373    alp_bl_conv(ig)=0.
374    endif
375  alp_bl_stat(ig)=alp_bl_det(ig)+alp_bl_fluct_m(ig)+alp_bl_fluct_tke(ig)+alp_bl_conv(ig)
376  enddo
377
378!-----Sécurité ALP infinie
379  do ig=1,ngrid
380   if (fraca0(ig).gt.0.98) alp_bl_stat(ig)=2.
381  enddo
382
383  ENDIF ! (iflag_clos_bl.ge.2)
384
385!!! fin nrlmd le 10/04/2012
386
387!      print*,'avant calcul ale et alp'
388!calcul de ALE et ALP pour la convection
389      alp_bl(:)=0.
390      ale_bl(:)=0.
391!          print*,'ALE,ALP ,l,zw2(ig,l),ale_bl(ig),alp_bl(ig)'
392      do l=1,nlay
393      do ig=1,ngrid
394           alp_bl(ig)=max(alp_bl(ig),0.5*rhobarz(ig,l)*wth3(ig,l) )
395           ale_bl(ig)=max(ale_bl(ig),0.5*zw2(ig,l)**2)
396!          print*,'ALE,ALP',l,zw2(ig,l),ale_bl(ig),alp_bl(ig)
397      enddo
398      enddo
399
400! ale sec (max de wmax/2 sous la zone d'inhibition) dans
401! le cas iflag_trig_bl=3
402      IF (iflag_trig_bl==3) ale_bl(:)=0.5*wmax_sec(:)**2
403
404!test:calcul de la ponderation des couches pour KE
405!initialisations
406
407      fm_tot(:)=0.
408      wght_th(:,:)=1.
409      lalim_conv(:)=lalim(:)
410
411      do k=1,klev
412         do ig=1,ngrid
413            if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k)
414         enddo
415      enddo
416
417! assez bizarre car, si on est dans la couche d'alim et que alim_star et
418! plus petit que 1.e-10, on prend wght_th=1.
419      do k=1,klev
420         do ig=1,ngrid
421            if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then
422               wght_th(ig,k)=alim_star(ig,k)
423            endif
424         enddo
425      enddo
426
427!      print*,'apres wght_th'
428!test pour prolonger la convection
429      do ig=1,ngrid
430!v1d  if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
431      if ((alim_star(ig,1).lt.1.e-10)) then
432      lalim_conv(ig)=1
433      wght_th(ig,1)=1.
434!      print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
435      endif
436      enddo
437
438!------------------------------------------------------------------------
439! Modif CR/FH 20110310 : alp integree sur la verticale.
440! Integrale verticale de ALP.
441! wth3 etant aux niveaux inter-couches, on utilise d play comme masse des
442! couches
443!------------------------------------------------------------------------
444
445      alp_int(:)=0.
446      dp_int(:)=0.
447      do l=2,nlay
448        do ig=1,ngrid
449           if(l.LE.lmax(ig)) THEN
450           zdp=pplay(ig,l-1)-pplay(ig,l)
451           alp_int(ig)=alp_int(ig)+0.5*rhobarz(ig,l)*wth3(ig,l)*zdp
452           dp_int(ig)=dp_int(ig)+zdp
453           endif
454        enddo
455      enddo
456
457      if (iflag_coupl>=3 .and. iflag_coupl<=5) then
458      do ig=1,ngrid
459!valeur integree de alp_bl * 0.5:
460        if (dp_int(ig)>0.) then
461        alp_bl(ig)=alp_int(ig)/dp_int(ig)
462        endif
463      enddo!
464      endif
465
466
467! Facteur multiplicatif sur alp_bl
468      alp_bl(:)=alp_bl_k*alp_bl(:)
469
470!------------------------------------------------------------------------
471
472
473
474      return
475      end
Note: See TracBrowser for help on using the repository browser.