Changeset 4171 for LMDZ6/branches/LMDZ-ECRAD/libf/phylmd/thermcell_alp.F90

Timestamp:

Jun 17, 2022, 4:24:49 PM (2 years ago)

Author:

lguez

Message:

Sync latest trunk changes to branch LMDZ-ECRAD.

Location:

LMDZ6/branches/LMDZ-ECRAD

Files:

• . (modified) (1 prop)
• libf/phylmd/thermcell_alp.F90 (modified) (8 diffs)

LMDZ6/branches/LMDZ-ECRAD/libf/phylmd/thermcell_alp.F90

@@ -1 +1 @@
 ! $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
+      SUBROUTINE thermcell_alp(ngrid,nlay,ptimestep  &                         ! in
+     &                  ,pplay,pplev  &                                        ! in
+     &                  ,fm0,entr0,lmax  &                                     ! in
+     &                  ,pbl_tke,pctsrf,omega,airephy &                        ! in
+     &                  ,zw2,fraca &                                           ! in
+     &                  ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &  ! in
+!
+     &                  ,ale_bl,alp_bl,lalim_conv,wght_th &                    ! out
+     &                  ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &   ! out
+     &                  ,n2,s2,ale_bl_stat &                                   ! out
+     &                  ,therm_tke_max,env_tke_max &                           ! out
+     &                  ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &          ! out
+     &                  ,alp_bl_conv,alp_bl_stat &                             ! out
      &)
 
-      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"
-!
+!   Auteurs: Catherine Rio
+!   Modifications :
+!   Nicolas Rochetin et Jean-Yves Grandpeix
+!         pour la fermeture stochastique. 2012
+!   Frédéric Hourdin :
+!         netoyage informatique. 2022
+!   
 !=======================================================================
 !-----------------------------------------------------------------------
@@ -58 +36 @@
 #include "YOETHF.h"
 #include "FCTTRE.h"
-#include "thermcell.h"
+#include "alpale.h"
 
 !   arguments:
 !   ----------
 
-!IM 140508
-
-      INTEGER ngrid,nlay
-      real ptimestep
-      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
-
-!   local:
-!   ------
-
+!------Entrees
+      integer, intent(in) :: ngrid,nlay
+      real, intent(in) :: ptimestep
+      real, intent(in) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1)
+      integer, intent(in), dimension(ngrid) ::lmax,lalim
+      real, intent(in), dimension(ngrid) :: zmax
+      real, intent(in), dimension(ngrid,nlay+1) :: zw2
+      real, intent(in), dimension(ngrid,nlay+1) :: fraca
+      real, intent(in), dimension(ngrid,nlay) :: wth3
+      real, intent(in), dimension(ngrid,nlay) :: rhobarz
+      real, intent(in), dimension(ngrid) :: wmax_sec
+      real, intent(in), dimension(ngrid,nlay) :: entr0
+      real, intent(in), dimension(ngrid,nlay+1) :: fm0,fm
+      real, intent(in), dimension(ngrid) :: pcon
+      real, intent(in), dimension(ngrid,nlay) :: alim_star
+      real, intent(in), dimension(ngrid,nlay+1,nbsrf) :: pbl_tke
+      real, intent(in), dimension(ngrid,nbsrf) :: pctsrf
+      real, intent(in), dimension(ngrid,nlay) :: omega
+      real, intent(in), dimension(ngrid) :: airephy
+!------Sorties
+      real, intent(out), dimension(ngrid) :: ale_bl,alp_bl
+      real, intent(out), dimension(ngrid,nlay) :: wght_th
+      integer, intent(out), dimension(ngrid) :: lalim_conv
+      real, intent(out), dimension(ngrid) :: zlcl,fraca0,w0,w_conv
+      real, intent(out), dimension(ngrid) :: therm_tke_max0,env_tke_max0,n2,s2,ale_bl_stat
+      real, intent(out), dimension(ngrid,nlay) :: therm_tke_max,env_tke_max
+      real, intent(out), dimension(ngrid) :: alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke
+      real, intent(out), dimension(ngrid) :: alp_bl_conv,alp_bl_stat
+
+!=============================================================================================
+!------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
+      real rhobarz0(ngrid)                    ! Densité au LCL
+      logical ok_lcl(ngrid)                   ! Existence du LCL des thermiques
+      integer klcl(ngrid)                     ! Niveau du LCL
+      real interp(ngrid)                      ! 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)
+      real, parameter :: su_cst=4e4              ! Surface unite: celle d'un updraft élémentaire
+      real, parameter :: hcoef=1             ! Coefficient directeur pour le calcul de s2
+      real, parameter :: hmincoef=0.3        ! Coefficient directeur pour l'ordonnée à l'origine pour le calcul de s2 
+      real, parameter :: eps1=0.3            ! Fraction de surface occupée par la population 1 : eps1=n1*s1/(fraca0*Sd)
+      real, dimension(ngrid) :: hmin         ! Ordonnée à l'origine pour le calcul de s2
+      real, dimension(ngrid) :: zmax_moy     ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
+      real, parameter :: zmax_moy_coef=0.33
+      real, dimension(ngrid) :: depth        ! Epaisseur moyenne du cumulus
+      real, dimension(ngrid) ::  w_max                 ! Vitesse max statistique 
+      real, dimension(ngrid) ::  s_max(ngrid)
 !--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
-
+      real, dimension(ngrid,nlay) :: pbl_tke_max       ! Profil de TKE moyenne 
+      real, dimension(ngrid) :: pbl_tke_max0           ! TKE moyenne au LCL
+      real, dimension(ngrid,nlay) :: w_ls              ! Vitesse verticale grande échelle (m/s)
+      real, parameter :: coef_m=1.            ! On considère un rendement pour alp_bl_fluct_m
+      real, parameter :: coef_tke=1.          ! On considère un rendement pour alp_bl_fluct_tke
+      real :: zdp
+      real, dimension(ngrid) :: alp_int,dp_int
+      real, dimension(ngrid) :: fm_tot
 
 !------------------------------------------------------------
 !  Initialize output arrays related to stochastic triggering
 !------------------------------------------------------------
-  DO ig = 1,klon
+  DO ig = 1,ngrid
      zlcl(ig) = 0.
      fraca0(ig) = 0.
@@ -175 +121 @@
      alp_bl_stat(ig) = 0.
   ENDDO
-  DO l = 1,klev
-    DO ig = 1,klon
+  DO l = 1,nlay
+    DO ig = 1,ngrid
      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.
+      if ( (pcon(ig) .gt. pplay(ig,nlay-1)) .and. (pcon(ig) .lt. pplay(ig,1)) ) ok_lcl(ig)=.true.
     enddo
 
@@ -207 +151 @@
     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
@@ -257 +193 @@
 
 !-----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)
+   call thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0,  &  ! in
+  &           rg,pplev,therm_tke_max)                               ! out
 !   print *,' thermcell_tke_transport -> '   !!jyg
 
@@ -330 +266 @@
 !   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)
+   susqr2pi=su_cst*sqrt(2.*Rpi)
    reuler=exp(1.)
    do ig=1,ngrid
@@ -409 +334 @@
       lalim_conv(:)=lalim(:)
 
-      do k=1,klev
+      do k=1,nlay
          do ig=1,ngrid
             if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k)
@@ -417 +342 @@
 ! 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 k=1,nlay
          do ig=1,ngrid
             if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then