Changeset 1669 for LMDZ5/branches/testing/libf/phylmd/thermcell_main.F90

Timestamp:

Oct 16, 2012, 2:41:50 PM (12 years ago)

Author:

Laurent Fairhead

Message:

Version testing basée sur la r1668

http://lmdz.lmd.jussieu.fr/utilisateurs/distribution-du-modele/versions-intermediaires

---

Testing release based on r1668

Location:

LMDZ5/branches/testing

Files:

• . (modified) (1 prop)
• libf/phylmd/thermcell_main.F90 (modified) (6 diffs)

LMDZ5/branches/testing/libf/phylmd/thermcell_main.F90

@@ -10 +10 @@
      &                  ,Ale_bl,Alp_bl,lalim_conv,wght_th &
      &                  ,zmax0, f0,zw2,fraca,ztv &
-     &                  ,zpspsk,ztla,zthl)
+     &                  ,zpspsk,ztla,zthl &
+!!! 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
@@ -47 +56 @@
 #include "iniprint.h"
 #include "thermcell.h"
+!!! nrlmd le 10/04/2012
+#include "indicesol.h"
+!!! fin nrlmd le 10/04/2012
 
 !   arguments:
@@ -77 +89 @@
       integer,save :: lev_out=10
 !$OMP THREADPRIVATE(lev_out)
+
+      REAL susqr2pi, Reuler
 
       INTEGER ig,k,l,ll,ierr
@@ -155 +169 @@
        real seuil
       real csc(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
 
 !
@@ -679 +740 @@
       enddo
 !
+
+!!! nrlmd le 10/04/2012
+
+!------------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 
+          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
+        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)
+!!    zmax(ig)=pphi(ig,lmax(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
+     zmax(ig)=pphi(ig,lmax(ig))/rg
+     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(:)=0.
+   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
+       s2(ig)=0.
+       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.1) 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.1)
+
+!!! fin nrlmd le 10/04/2012
+
       if (prt_level.ge.10) then
          ig=igout
@@ -858 +1120 @@
       end
 
+!!! nrlmd le 10/04/2012                          Transport de la TKE par le thermique moyen pour la fermeture en ALP 
+!                                                         On transporte pbl_tke pour donner therm_tke
+!                                          Copie conforme de la subroutine DTKE dans physiq.F écrite par Frederic Hourdin
+      subroutine thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0,  &
+     &           rg,pplev,therm_tke_max)
+      implicit none
+
+#include "iniprint.h"
+!=======================================================================
+!
+!   Calcul du transport verticale dans la couche limite en presence
+!   de "thermiques" explicitement representes
+!   calcul du dq/dt une fois qu'on connait les ascendances
+!
+!=======================================================================
+
+      integer ngrid,nlay,nsrf
+
+      real ptimestep
+      real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)
+      real entr0(ngrid,nlay),rg
+      real therm_tke_max(ngrid,nlay)
+      real detr0(ngrid,nlay)
+
+
+      real masse(ngrid,nlay),fm(ngrid,nlay+1)
+      real entr(ngrid,nlay)
+      real q(ngrid,nlay)
+      integer lev_out                           ! niveau pour les print
+
+      real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
+
+      real zzm
+
+      integer ig,k
+      integer isrf
+
+
+      lev_out=0
+
+
+      if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
+
+!   calcul du detrainement
+      do k=1,nlay
+         detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
+         masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG
+      enddo
+
+
+! Decalage vertical des entrainements et detrainements.
+      masse(:,1)=0.5*masse0(:,1)
+      entr(:,1)=0.5*entr0(:,1)
+      detr(:,1)=0.5*detr0(:,1)
+      fm(:,1)=0.
+      do k=1,nlay-1
+         masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
+         entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
+         detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
+         fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
+      enddo
+      fm(:,nlay+1)=0.
+
+!!! nrlmd le 16/09/2010
+!   calcul de la valeur dans les ascendances
+!       do ig=1,ngrid
+!          qa(ig,1)=q(ig,1)
+!       enddo
+!!! 
+
+!do isrf=1,nsrf
+
+!   q(:,:)=therm_tke(:,:,isrf)
+   q(:,:)=therm_tke_max(:,:)
+!!! nrlmd le 16/09/2010
+      do ig=1,ngrid
+         qa(ig,1)=q(ig,1)
+      enddo
+!!!
+
+    if (1==1) then
+      do k=2,nlay
+         do ig=1,ngrid
+            if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
+     &         1.e-5*masse(ig,k)) then
+         qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
+     &         /(fm(ig,k+1)+detr(ig,k))
+            else
+               qa(ig,k)=q(ig,k)
+            endif
+            if (qa(ig,k).lt.0.) then
+!               print*,'qa<0!!!'
+            endif
+            if (q(ig,k).lt.0.) then
+!               print*,'q<0!!!'
+            endif
+         enddo
+      enddo
+
+! Calcul du flux subsident
+
+      do k=2,nlay
+         do ig=1,ngrid
+            wqd(ig,k)=fm(ig,k)*q(ig,k)
+            if (wqd(ig,k).lt.0.) then
+!               print*,'wqd<0!!!'
+            endif
+         enddo
+      enddo
+      do ig=1,ngrid
+         wqd(ig,1)=0.
+         wqd(ig,nlay+1)=0.
+      enddo
+
+! Calcul des tendances
+      do k=1,nlay
+         do ig=1,ngrid
+            q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)  &
+     &               -wqd(ig,k)+wqd(ig,k+1))  &
+     &               *ptimestep/masse(ig,k)
+         enddo
+      enddo
+
+ endif
+
+   therm_tke_max(:,:)=q(:,:)
+
+      return
+!!! fin nrlmd le 10/04/2012
+     end
+