thermcell_alp.F90

Go to the documentation of this file.

! $Id: thermcell_main.F90 2351 2015-08-25 15:14:59Z emillour $
!
      SUBROUTINE thermcell_alp(ngrid,nlay,ptimestep  &
     &                  ,pplay,pplev  &
     &                  ,fm0,entr0,lmax  &
     &                  ,ale_bl,alp_bl,lalim_conv,wght_th &
     &                  ,zw2,fraca &
!!! ncessaire en plus
     &                  ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &
!!! nrlmd le 10/04/2012
     &                  ,pbl_tke,pctsrf,omega,airephy &
     &                  ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
     &                  ,n2,s2,ale_bl_stat &
     &                  ,therm_tke_max,env_tke_max &
     &                  ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
     &                  ,alp_bl_conv,alp_bl_stat &
!!! fin nrlmd le 10/04/2012
     &)

      USE dimphy
      USE indice_sol_mod
      IMPLICIT NONE

!=======================================================================
!   Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
!   Version du 09.02.07
!   Calcul du transport vertical dans la couche limite en presence
!   de "thermiques" explicitement representes avec processus nuageux
!
!   Reecriture a partir d'un listing papier a Habas, le 14/02/00
!
!   le thermique est suppose homogene et dissipe par melange avec
!   son environnement. la longueur l_mix controle l'efficacite du
!   melange
!
!   Le calcul du transport des differentes especes se fait en prenant
!   en compte:
!     1. un flux de masse montant
!     2. un flux de masse descendant
!     3. un entrainement
!     4. un detrainement
!
! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
!    Introduction of an implicit computation of vertical advection in
!    the environment of thermal plumes in thermcell_dq
!    impl =     0 : explicit, 1 : implicit, -1 : old version
!    controled by iflag_thermals =
!       15, 16 run with impl=-1 : numerical convergence with NPv3
!       17, 18 run with impl=1  : more stable
!    15 and 17 correspond to the activation of the stratocumulus "bidouille"
!
!=======================================================================
!-----------------------------------------------------------------------
!   declarations:
!   -------------

#include "YOMCST.h"
#include "YOETHF.h"
#include "FCTTRE.h"
#include "thermcell.h"

!   arguments:
!   ----------

!IM 140508

      INTEGER ngrid,nlay
      real ptimestep
      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)

!   local:
!   ------


      REAL susqr2pi, reuler

      INTEGER ig,k,l
      INTEGER lmax(klon),lalim(klon)
      real zmax(klon),zw2(klon,klev+1)

!on garde le zmax du pas de temps precedent


      real fraca(klon,klev+1)
      real wth3(klon,klev)
! FH probleme de dimensionnement avec l'allocation dynamique
!     common/comtherm/thetath2,wth2
      real rhobarz(klon,klev)

      real wmax_sec(klon)
      real fm0(klon,klev+1),entr0(klon,klev)
      real fm(klon,klev+1)

!niveau de condensation
      real pcon(klon)

      real alim_star(klon,klev)

!!! nrlmd le 10/04/2012

!------Entrées
      real pbl_tke(klon,klev+1,nbsrf)
      real pctsrf(klon,nbsrf)
      real omega(klon,klev)
      real airephy(klon)
!------Sorties
      real zlcl(klon),fraca0(klon),w0(klon),w_conv(klon)
      real therm_tke_max0(klon),env_tke_max0(klon)
      real n2(klon),s2(klon)
      real ale_bl_stat(klon)
      real therm_tke_max(klon,klev),env_tke_max(klon,klev)
      real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon)
!------Local
      integer nsrf
      real rhobarz0(klon)                    ! Densité au LCL
      logical ok_lcl(klon)                   ! Existence du LCL des thermiques
      integer klcl(klon)                     ! Niveau du LCL
      real interp(klon)                      ! Coef d'interpolation pour le LCL
!--Triggering
      real Su                                ! Surface unité: celle d'un updraft élémentaire
      parameter(su=4e4)
      real hcoef                             ! Coefficient directeur pour le calcul de s2
      parameter(hcoef=1)
      real hmincoef                          ! Coefficient directeur pour l'ordonnée à l'origine pour le calcul de s2 
      parameter(hmincoef=0.3)
      real eps1                              ! Fraction de surface occupée par la population 1 : eps1=n1*s1/(fraca0*Sd)
      parameter(eps1=0.3)
      real hmin(ngrid)                       ! Ordonnée à l'origine pour le calcul de s2
      real zmax_moy(ngrid)                   ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
      real zmax_moy_coef
      parameter(zmax_moy_coef=0.33)
      real depth(klon)                       ! Epaisseur moyenne du cumulus
      real w_max(klon)                       ! Vitesse max statistique 
      real s_max(klon)
!--Closure
      real pbl_tke_max(klon,klev)            ! Profil de TKE moyenne 
      real pbl_tke_max0(klon)                ! TKE moyenne au LCL
      real w_ls(klon,klev)                   ! Vitesse verticale grande échelle (m/s)
      real coef_m                            ! On considère un rendement pour alp_bl_fluct_m
      parameter(coef_m=1.)
      real coef_tke                          ! On considère un rendement pour alp_bl_fluct_tke
      parameter(coef_tke=1.)

!!! fin nrlmd le 10/04/2012

!
      !nouvelles variables pour la convection
      real ale_bl(klon)
      real alp_bl(klon)
      real alp_int(klon),dp_int(klon),zdp
      real fm_tot(klon)
      real wght_th(klon,klev)
      integer lalim_conv(klon)
!v1d     logical therm
!v1d     save therm


!------------------------------------------------------------
!  Initialize output arrays related to stochastic triggering
!------------------------------------------------------------
  DO ig = 1,klon
     zlcl(ig) = 0.
     fraca0(ig) = 0.
     w0(ig) = 0.
     w_conv(ig) = 0.
     therm_tke_max0(ig) = 0.
     env_tke_max0(ig) = 0.
     n2(ig) = 0.
     s2(ig) = 0.
     ale_bl_stat(ig) = 0.
     alp_bl_det(ig) = 0.
     alp_bl_fluct_m(ig) = 0.
     alp_bl_fluct_tke(ig) = 0.
     alp_bl_conv(ig) = 0.
     alp_bl_stat(ig) = 0.
  ENDDO
  DO l = 1,klev
    DO ig = 1,klon
     therm_tke_max(ig,l) = 0.
     env_tke_max(ig,l) = 0.
    ENDDO
  ENDDO
!------------------------------------------------------------


!------------Test sur le LCL des thermiques
    do ig=1,ngrid
      ok_lcl(ig)=.false.
      if ( (pcon(ig) .gt. pplay(ig,klev-1)) .and. (pcon(ig) .lt. pplay(ig,1)) ) ok_lcl(ig)=.true.
    enddo

!------------Localisation des niveaux entourant le LCL et du coef d'interpolation
    do l=1,nlay-1
      do ig=1,ngrid
        if (ok_lcl(ig)) then 
!ATTENTION,zw2 calcule en pplev
!          if ((pplay(ig,l) .ge. pcon(ig)) .and. (pplay(ig,l+1) .le. pcon(ig))) then
!          klcl(ig)=l
!          interp(ig)=(pcon(ig)-pplay(ig,klcl(ig)))/(pplay(ig,klcl(ig)+1)-pplay(ig,klcl(ig)))
!          endif
          if ((pplev(ig,l) .ge. pcon(ig)) .and. (pplev(ig,l+1) .le. pcon(ig))) then
          klcl(ig)=l
          interp(ig)=(pcon(ig)-pplev(ig,klcl(ig)))/(pplev(ig,klcl(ig)+1)-pplev(ig,klcl(ig)))
          endif
        endif
      enddo
    enddo

!------------Hauteur des thermiques
!!jyg le 27/04/2012
!!    do ig =1,ngrid
!!    rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
!! &               -rhobarz(ig,klcl(ig)))*interp(ig)
!!    zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
!!      if ( (.not.ok_lcl(ig)) .or. (zlcl(ig).gt.zmax(ig)) ) zlcl(ig)=zmax(ig) ! Si zclc > zmax alors on pose zlcl = zmax
!!    enddo
    do ig =1,ngrid
!CR:REHABILITATION ZMAX CONTINU
     if (ok_lcl(ig)) then 
      rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
 &               -rhobarz(ig,klcl(ig)))*interp(ig)
      zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*rg)
      zlcl(ig)=min(zlcl(ig),zmax(ig))   ! Si zlcl > zmax alors on pose zlcl = zmax
     else
      rhobarz0(ig)=0.
      zlcl(ig)=zmax(ig)
     endif
    enddo
!!jyg fin

!------------Calcul des propriétés du thermique au LCL 
  IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) THEN 

  !-----Initialisation de la TKE moyenne 
   do l=1,nlay
    do ig=1,ngrid
     pbl_tke_max(ig,l)=0.
    enddo
   enddo

!-----Calcul de la TKE moyenne 
   do nsrf=1,nbsrf
    do l=1,nlay
     do ig=1,ngrid
     pbl_tke_max(ig,l)=pctsrf(ig,nsrf)*pbl_tke(ig,l,nsrf)+pbl_tke_max(ig,l)
     enddo
    enddo
   enddo

!-----Initialisations des TKE dans et hors des thermiques 
   do l=1,nlay
    do ig=1,ngrid
    therm_tke_max(ig,l)=pbl_tke_max(ig,l)
    env_tke_max(ig,l)=pbl_tke_max(ig,l)
    enddo
   enddo

!-----Calcul de la TKE transportée par les thermiques : therm_tke_max
   call thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0,  &
  &           rg,pplev,therm_tke_max)
!   print *,' thermcell_tke_transport -> '   !!jyg

!-----Calcul des profils verticaux de TKE hors thermiques : env_tke_max, et de la vitesse verticale grande échelle : W_ls
   do l=1,nlay
    do ig=1,ngrid
     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
     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
     w_ls(ig,l)=-1.*omega(ig,l)/(rg*rhobarz(ig,l))                                                !  Vitesse verticale de grande échelle
    enddo
   enddo
!    print *,' apres w_ls = '   !!jyg

  do ig=1,ngrid
   if (ok_lcl(ig)) then
     fraca0(ig)=fraca(ig,klcl(ig))+(fraca(ig,klcl(ig)+1) &
 &             -fraca(ig,klcl(ig)))*interp(ig)
     w0(ig)=zw2(ig,klcl(ig))+(zw2(ig,klcl(ig)+1) &
 &         -zw2(ig,klcl(ig)))*interp(ig)
     w_conv(ig)=w_ls(ig,klcl(ig))+(w_ls(ig,klcl(ig)+1) &
 &             -w_ls(ig,klcl(ig)))*interp(ig)
     therm_tke_max0(ig)=therm_tke_max(ig,klcl(ig)) &
 &                     +(therm_tke_max(ig,klcl(ig)+1)-therm_tke_max(ig,klcl(ig)))*interp(ig)
     env_tke_max0(ig)=env_tke_max(ig,klcl(ig))+(env_tke_max(ig,klcl(ig)+1) &
 &                   -env_tke_max(ig,klcl(ig)))*interp(ig)
     pbl_tke_max0(ig)=pbl_tke_max(ig,klcl(ig))+(pbl_tke_max(ig,klcl(ig)+1) &
 &                   -pbl_tke_max(ig,klcl(ig)))*interp(ig)
     if (therm_tke_max0(ig).ge.20.) therm_tke_max0(ig)=20.
     if (env_tke_max0(ig).ge.20.) env_tke_max0(ig)=20.
     if (pbl_tke_max0(ig).ge.20.) pbl_tke_max0(ig)=20.
   else 
     fraca0(ig)=0.
     w0(ig)=0.
!!jyg le 27/04/2012
!!     zlcl(ig)=0.
!!
   endif
  enddo

  ENDIF ! IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) )
!  print *,'ENDIF  ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) '    !!jyg

!------------Triggering------------------
  IF (iflag_trig_bl.ge.1) THEN 

!-----Initialisations
   depth(:)=0.
   n2(:)=0.
   s2(:)=100. ! some low value, arbitrary
   s_max(:)=0.

!-----Epaisseur du nuage (depth) et détermination de la queue du spectre de panaches (n2,s2) et du panache le plus gros (s_max)
   do ig=1,ngrid
     zmax_moy(ig)=zlcl(ig)+zmax_moy_coef*(zmax(ig)-zlcl(ig))
     depth(ig)=zmax_moy(ig)-zlcl(ig)
     hmin(ig)=hmincoef*zlcl(ig)
     if (depth(ig).ge.10.) then 
       s2(ig)=(hcoef*depth(ig)+hmin(ig))**2
       n2(ig)=(1.-eps1)*fraca0(ig)*airephy(ig)/s2(ig)
!!
!!jyg le 27/04/2012
!!       s_max(ig)=s2(ig)*log(n2(ig))
!!       if (n2(ig) .lt. 1) s_max(ig)=0.
       s_max(ig)=s2(ig)*log(max(n2(ig),1.))
!!fin jyg
     else
       n2(ig)=0.
       s_max(ig)=0.
     endif
   enddo
!   print *,'avant Calcul de Wmax '    !!jyg

!-----Calcul de Wmax et ALE_BL_STAT associée
!!jyg le 30/04/2012
!!   do ig=1,ngrid
!!     if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.1.) ) then
!!     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))))
!!     ale_bl_stat(ig)=0.5*w_max(ig)**2
!!     else
!!     w_max(ig)=0.
!!     ale_bl_stat(ig)=0.
!!     endif
!!   enddo
   susqr2pi=su*sqrt(2.*rpi)
   reuler=exp(1.)
   do ig=1,ngrid
     if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.susqr2pi*reuler) ) then
      w_max(ig)=w0(ig)*(1.+sqrt(2.*log(s_max(ig)/susqr2pi)-log(2.*log(s_max(ig)/susqr2pi))))
      ale_bl_stat(ig)=0.5*w_max(ig)**2
     else
      w_max(ig)=0.
      ale_bl_stat(ig)=0.
     endif
   enddo

  ENDIF ! iflag_trig_bl
!  print *,'ENDIF  iflag_trig_bl'    !!jyg

!------------Closure------------------

  IF (iflag_clos_bl.ge.2) THEN 

!-----Calcul de ALP_BL_STAT
  do ig=1,ngrid
  alp_bl_det(ig)=0.5*coef_m*rhobarz0(ig)*(w0(ig)**3)*fraca0(ig)*(1.-2.*fraca0(ig))/((1.-fraca0(ig))**2)
  alp_bl_fluct_m(ig)=1.5*rhobarz0(ig)*fraca0(ig)*(w_conv(ig)+coef_m*w0(ig))* &
 &                   (w0(ig)**2)
  alp_bl_fluct_tke(ig)=3.*coef_m*rhobarz0(ig)*w0(ig)*fraca0(ig)*(therm_tke_max0(ig)-env_tke_max0(ig)) &
 &                    +3.*rhobarz0(ig)*w_conv(ig)*pbl_tke_max0(ig)
    if (iflag_clos_bl.ge.2) then 
    alp_bl_conv(ig)=1.5*coef_m*rhobarz0(ig)*fraca0(ig)*(fraca0(ig)/(1.-fraca0(ig)))*w_conv(ig)* &
 &                   (w0(ig)**2)
    else
    alp_bl_conv(ig)=0.
    endif
  alp_bl_stat(ig)=alp_bl_det(ig)+alp_bl_fluct_m(ig)+alp_bl_fluct_tke(ig)+alp_bl_conv(ig)
  enddo

!-----Sécurité ALP infinie
  do ig=1,ngrid
   if (fraca0(ig).gt.0.98) alp_bl_stat(ig)=2.
  enddo

  ENDIF ! (iflag_clos_bl.ge.2)

!!! fin nrlmd le 10/04/2012

!      print*,'avant calcul ale et alp' 
!calcul de ALE et ALP pour la convection
      alp_bl(:)=0.
      ale_bl(:)=0.
!          print*,'ALE,ALP ,l,zw2(ig,l),ale_bl(ig),alp_bl(ig)'
      do l=1,nlay
      do ig=1,ngrid
           alp_bl(ig)=max(alp_bl(ig),0.5*rhobarz(ig,l)*wth3(ig,l) )
           ale_bl(ig)=max(ale_bl(ig),0.5*zw2(ig,l)**2)
!          print*,'ALE,ALP',l,zw2(ig,l),ale_bl(ig),alp_bl(ig)
      enddo
      enddo

! ale sec (max de wmax/2 sous la zone d'inhibition) dans
! le cas iflag_trig_bl=3
      IF (iflag_trig_bl==3) ale_bl(:)=0.5*wmax_sec(:)**2

!test:calcul de la ponderation des couches pour KE
!initialisations

      fm_tot(:)=0.
      wght_th(:,:)=1.
      lalim_conv(:)=lalim(:)

      do k=1,klev
         do ig=1,ngrid
            if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k)
         enddo
      enddo

! assez bizarre car, si on est dans la couche d'alim et que alim_star et
! plus petit que 1.e-10, on prend wght_th=1.
      do k=1,klev
         do ig=1,ngrid
            if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then
               wght_th(ig,k)=alim_star(ig,k)
            endif
         enddo
      enddo

!      print*,'apres wght_th'
!test pour prolonger la convection
      do ig=1,ngrid
!v1d  if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
      if ((alim_star(ig,1).lt.1.e-10)) then
      lalim_conv(ig)=1
      wght_th(ig,1)=1.
!      print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
      endif
      enddo

!------------------------------------------------------------------------
! Modif CR/FH 20110310 : alp integree sur la verticale.
! Integrale verticale de ALP.
! wth3 etant aux niveaux inter-couches, on utilise d play comme masse des
! couches
!------------------------------------------------------------------------

      alp_int(:)=0.
      dp_int(:)=0.
      do l=2,nlay
        do ig=1,ngrid
           if(l.LE.lmax(ig)) THEN
           zdp=pplay(ig,l-1)-pplay(ig,l)
           alp_int(ig)=alp_int(ig)+0.5*rhobarz(ig,l)*wth3(ig,l)*zdp
           dp_int(ig)=dp_int(ig)+zdp
           endif
        enddo
      enddo

      if (iflag_coupl>=3 .and. iflag_coupl<=5) then
      do ig=1,ngrid
!valeur integree de alp_bl * 0.5:
        if (dp_int(ig)>0.) then
        alp_bl(ig)=alp_int(ig)/dp_int(ig)
        endif
      enddo!
      endif


! Facteur multiplicatif sur alp_bl
      alp_bl(:)=alp_bl_k*alp_bl(:)

!------------------------------------------------------------------------



      return
      end