Index: LMDZ5/trunk/libf/phylmd/physiq.F
===================================================================
--- LMDZ5/trunk/libf/phylmd/physiq.F	(revision 1628)
+++ LMDZ5/trunk/libf/phylmd/physiq.F	(revision 1638)
@@ -180,5 +180,5 @@
       real facteur,zfratqs1,zfratqs2
 
-      REAL lambda_th(klon,klev),zz,znum,zden
+      REAL zz,znum,zden
       REAL wmax_th(klon)
       REAL zmax_th(klon)
@@ -614,6 +614,6 @@
       REAL dd_t(klon,klev),dd_q(klon,klev)
 
-      real, save :: alp_bl_prescr=0.1
-      real, save :: ale_bl_prescr=4.
+      real, save :: alp_bl_prescr=0.
+      real, save :: ale_bl_prescr=0.
 
       real, save :: ale_max=1000.
@@ -689,4 +689,39 @@
       REAL ztla(klon,klev) 
       REAL zthl(klon,klev)
+
+ccc nrlmd le 10/04/2012
+
+c--------Stochastic Boundary Layer Triggering: ALE_BL--------
+c---Propriétés du thermiques au LCL 
+      real zlcl_th(klon)                                     ! Altitude du LCL calculé continument (pcon dans thermcell_main.F90)
+      real fraca0(klon)                                      ! Fraction des thermiques au LCL
+      real w0(klon)                                          ! Vitesse des thermiques au LCL
+      real w_conv(klon)                                      ! Vitesse verticale de grande échelle au LCL
+      real therm_tke_max0(klon)                              ! TKE dans les thermiques au LCL 
+      real env_tke_max0(klon)                                ! TKE dans l'environnement au LCL 
+
+c---Spectre de thermiques de type 2 au LCL
+      real n2(klon),s2(klon)
+      real ale_bl_stat(klon)
+
+c---Déclenchement stochastique
+      integer :: tau_trig(klon)
+      real proba_notrig(klon)
+      real random_notrig(klon)
+
+c--------Statistical Boundary Layer Closure: ALP_BL--------
+c---Profils de TKE dans et hors du thermique
+      real pbl_tke_input(klon,klev+1,nbsrf)
+      real therm_tke_max(klon,klev)                          ! Profil de TKE dans les thermiques
+      real env_tke_max(klon,klev)                            ! Profil de TKE dans l'environnement
+
+c---Fermeture statistique
+      real alp_bl_det(klon)                                     ! ALP déterministe du thermique unique
+      real alp_bl_fluct_m(klon)                                 ! ALP liée aux fluctuations de flux de masse sous-nuageux
+      real alp_bl_fluct_tke(klon)                               ! ALP liée aux fluctuations d'énergie cinétique sous-nuageuse 
+      real alp_bl_conv(klon)                                    ! ALP liée à grande échelle
+      real alp_bl_stat(klon)                                    ! ALP totale
+
+ccc fin nrlmd le 10/04/2012
 
 c Variables locales pour la couche limite (al1):
@@ -791,4 +826,7 @@
 cIM
       EXTERNAL haut2bas  !variables de haut en bas
+      INTEGER lnblnk1
+      EXTERNAL lnblnk1   !enleve les blancs a la fin d'une variable de type
+                         !caracter
       EXTERNAL ini_undefSTD  !initialise a 0 une variable a 1 niveau de pression
       EXTERNAL undefSTD      !somme les valeurs definies d'1 var a 1 niveau de pression
@@ -1354,5 +1392,4 @@
          solswad(:)=0.
 
-         lambda_th(:,:)=0.
          wmax_th(:)=0.
          tau_overturning_th(:)=0.
@@ -1490,8 +1527,10 @@
 cCR:04.12.07: initialisations poches froides
 c Controle de ALE et ALP pour la fermeture convective (jyg)
-         CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr
+          if (iflag_wake>=1) then
+            CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr
      s                  ,alp_bl_prescr, ale_bl_prescr)
 c 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
 c        print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
+          endif
 
         do i = 1,klon
@@ -1516,4 +1555,5 @@
        print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
       ENDIF
+
 c
       ALLOCATE(tabCFMIP(nCFMIP))
@@ -1735,4 +1775,5 @@
 !
       itap   = itap + 1
+c
 !
 ! Update fraction of the sub-surfaces (pctsrf) and 
@@ -2042,58 +2083,57 @@
 c
 
-      if (iflag_pbl/=0) then 
-
-        CALL pbl_surface( 
-     e       dtime,     date0,     itap,    days_elapsed+1,
-     e       debut,     lafin,
-     e       rlon,      rlat,      rugoro,  rmu0,     
-     e       rain_fall, snow_fall, solsw,   sollw,    
-     e       t_seri,    q_seri,    u_seri,  v_seri,   
-     e       pplay,     paprs,     pctsrf,            
-     +       ftsol,     falb1,     falb2,   u10m,   v10m,
-     s       sollwdown, cdragh,    cdragm,  u1,    v1,
-     s       albsol1,   albsol2,   sens,    evap,  
-     s       zxtsol,    zxfluxlat, zt2m,    qsat2m, 
-     s       d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf,
-     s       coefh,     coefm,     slab_wfbils,                
-     d       qsol,      zq2m,      s_pblh,  s_lcl,
-     d       s_capCL,   s_oliqCL,  s_cteiCL,s_pblT,
-     d       s_therm,   s_trmb1,   s_trmb2, s_trmb3,
-     d       zxrugs,    zu10m,     zv10m,   fder,
-     d       zxqsurf,   rh2m,      zxfluxu, zxfluxv,
-     d       frugs,     agesno,    fsollw,  fsolsw,
-     d       d_ts,      fevap,     fluxlat, t2m,
-     d       wfbils,    wfbilo,    fluxt,   fluxu,  fluxv,
-     -       dsens,     devap,     zxsnow,
-     -       zxfluxt,   zxfluxq,   q2m,     fluxq, pbl_tke )
+      if (iflag_pbl/=0) then
+
+      CALL pbl_surface( 
+     e     dtime,     date0,     itap,    days_elapsed+1,
+     e     debut,     lafin,
+     e     rlon,      rlat,      rugoro,  rmu0,     
+     e     rain_fall, snow_fall, solsw,   sollw,    
+     e     t_seri,    q_seri,    u_seri,  v_seri,   
+     e     pplay,     paprs,     pctsrf,            
+     +     ftsol,     falb1,     falb2,   u10m,   v10m,
+     s     sollwdown, cdragh,    cdragm,  u1,    v1,
+     s     albsol1,   albsol2,   sens,    evap,  
+     s     zxtsol,    zxfluxlat, zt2m,    qsat2m, 
+     s     d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf,
+     s     coefh,     coefm,     slab_wfbils,                
+     d     qsol,      zq2m,      s_pblh,  s_lcl,
+     d     s_capCL,   s_oliqCL,  s_cteiCL,s_pblT,
+     d     s_therm,   s_trmb1,   s_trmb2, s_trmb3,
+     d     zxrugs,    zu10m,     zv10m,   fder,
+     d     zxqsurf,   rh2m,      zxfluxu, zxfluxv,
+     d     frugs,     agesno,    fsollw,  fsolsw,
+     d     d_ts,      fevap,     fluxlat, t2m,
+     d     wfbils,    wfbilo,    fluxt,   fluxu,  fluxv,
+     -     dsens,     devap,     zxsnow,
+     -     zxfluxt,   zxfluxq,   q2m,     fluxq, pbl_tke )
 
 
 !-----------------------------------------------------------------------------------------
 ! ajout des tendances de la diffusion turbulente
-        CALL add_phys_tend(d_u_vdf,d_v_vdf,d_t_vdf,d_q_vdf,dql0,'vdf')
+      CALL add_phys_tend(d_u_vdf,d_v_vdf,d_t_vdf,d_q_vdf,dql0,'vdf')
 !-----------------------------------------------------------------------------------------
 
-        if (mydebug) then
-          call writefield_phy('u_seri',u_seri,llm)
-          call writefield_phy('v_seri',v_seri,llm)
-          call writefield_phy('t_seri',t_seri,llm)
-          call writefield_phy('q_seri',q_seri,llm)
-        endif
-
-
-        IF (ip_ebil_phy.ge.2) THEN 
-          ztit='after surface_main'
-          CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
+      if (mydebug) then
+        call writefield_phy('u_seri',u_seri,llm)
+        call writefield_phy('v_seri',v_seri,llm)
+        call writefield_phy('t_seri',t_seri,llm)
+        call writefield_phy('q_seri',q_seri,llm)
+      endif
+
+
+      IF (ip_ebil_phy.ge.2) THEN 
+        ztit='after surface_main'
+        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
      e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
      s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-          call diagphy(airephy,ztit,ip_ebil_phy
+         call diagphy(airephy,ztit,ip_ebil_phy
      e      , zero_v, zero_v, zero_v, zero_v, sens
      e      , evap  , zero_v, zero_v, ztsol
      e      , d_h_vcol, d_qt, d_ec
      s      , fs_bound, fq_bound )
-        END IF 
+      END IF 
 
       ENDIF
-
 c =================================================================== c
 c   Calcul de Qsat
@@ -2244,10 +2284,23 @@
 cdans le thermique sinon
        if (iflag_coupl.eq.0) then
-          if (debut.and.prt_level.gt.9)WRITE(lunout,*) 'ALE&ALP imposes'
-          Ale_bl(1:klon) = ale_bl_prescr
-          Alp_bl(1:klon) = alp_bl_prescr
+          if (debut.and.prt_level.gt.9)
+     $                     WRITE(lunout,*)'ALE et ALP imposes'
+          do i = 1,klon
+con ne couple que ale
+c           ALE(i) = max(ale_wake(i),Ale_bl(i))
+            ALE(i) = max(ale_wake(i),ale_bl_prescr)
+con ne couple que alp
+c           ALP(i) = alp_wake(i) + Alp_bl(i)
+            ALP(i) = alp_wake(i) + alp_bl_prescr
+          enddo
        else
          IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
-       endif
+!         do i = 1,klon
+!             ALE(i) = max(ale_wake(i),Ale_bl(i))
+! avant        ALP(i) = alp_wake(i) + Alp_bl(i)
+!             ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
+!         write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
+!         write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
+!         enddo
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -2256,7 +2309,11 @@
 ! w si <0
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
        do i = 1,klon
           ALE(i) = max(ale_wake(i),Ale_bl(i))
+ccc nrlmd le 10/04/2012----------Stochastic triggering--------------
+          if (iflag_trig_bl.ge.1) then
+             ALE(i) = max(ale_wake(i),Ale_bl_trig(i))
+          endif
+ccc fin nrlmd le 10/04/2012
           if (alp_offset>=0.) then
             ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
@@ -2269,7 +2326,7 @@
           endif
        enddo
-
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
+       endif
        do i=1,klon
           if (alp(i)>alp_max) then
@@ -2586,5 +2643,16 @@
 
 
-         if (iflag_thermals.gt.1) then
+ccc nrlmd le 10/04/2012
+         DO k=1,klev+1
+           DO i=1,klon
+           pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
+           pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
+           pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
+           pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
+           ENDDO
+         ENDDO
+ccc fin nrlmd le 10/04/2012
+
+         if (iflag_thermals>=1) then
          call calltherm(pdtphys
      s      ,pplay,paprs,pphi,weak_inversion
@@ -2596,5 +2664,85 @@
 con rajoute ale et alp, et les caracteristiques de la couche alim
      s      ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca
-     s      ,ztv,zpspsk,ztla,zthl)
+     s      ,ztv,zpspsk,ztla,zthl
+ccc nrlmd le 10/04/2012
+     e      ,pbl_tke_input,pctsrf,omega,airephy
+     s      ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0
+     s      ,n2,s2,ale_bl_stat
+     s      ,therm_tke_max,env_tke_max
+     s      ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke
+     s      ,alp_bl_conv,alp_bl_stat
+ccc fin nrlmd le 10/04/2012
+     s                 )
+
+ccc nrlmd le 10/04/2012
+c-----------Stochastic triggering-----------
+      if (iflag_trig_bl.ge.1) then
+c
+        IF (prt_level .GE. 10) THEN
+         print *,'cin, ale_bl_stat, alp_bl_stat ',
+     $            cin, ale_bl_stat, alp_bl_stat
+        ENDIF
+
+c----Initialisations
+      do i=1,klon
+      proba_notrig(i)=1.
+      random_notrig(i)=1e6*ale_bl_stat(i)-int(1e6*ale_bl_stat(i))
+        if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then 
+        tau_trig(i)=tau_trig_shallow
+        else
+        tau_trig(i)=tau_trig_deep
+        endif
+      enddo
+c
+        IF (prt_level .GE. 10) THEN
+         print *,'random_notrig, tau_trig ',
+     $            random_notrig, tau_trig
+          print *,'s_trig,s2,n2 ',
+     $             s_trig,s2,n2
+        ENDIF
+
+c----Tirage aléatoire et calcul de ale_bl_trig
+      do i=1,klon
+        if ( (ale_bl_stat(i) .gt. abs(cin(i))+1.e-10) )  then
+        proba_notrig(i)=(1.-exp(-s_trig/s2(i)))**
+     $                  (n2(i)*dtime/tau_trig(i))
+c        print *, 'proba_notrig(i) ',proba_notrig(i)
+          if (random_notrig(i) .ge. proba_notrig(i)) then 
+          ale_bl_trig(i)=ale_bl_stat(i)
+          else
+          ale_bl_trig(i)=0.
+          endif
+        else
+        proba_notrig(i)=1.
+        random_notrig(i)=0.
+        ale_bl_trig(i)=0.
+        endif
+      enddo
+c
+        IF (prt_level .GE. 10) THEN
+         print *,'proba_notrig, ale_bl_trig ',
+     $            proba_notrig, ale_bl_trig
+        ENDIF
+
+      endif !(iflag_trig_bl)
+
+c-----------Statistical closure-----------
+      if (iflag_clos_bl.ge.1) then 
+
+        do i=1,klon
+        alp_bl(i)=alp_bl_stat(i)
+        enddo
+
+      else
+
+      alp_bl_stat(:)=0.
+
+      endif !(iflag_clos_bl)
+
+        IF (prt_level .GE. 10) THEN
+         print *,'ale_bl_trig, alp_bl_stat ',ale_bl_trig, alp_bl_stat
+        ENDIF
+
+ccc fin nrlmd le 10/04/2012
 
 ! ----------------------------------------------------------------------
@@ -2627,5 +2775,5 @@
 c  ==============
 
-! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
+! Dans le cas oÃ¹ on active les thermiques, on fait partir l'ajustement
 ! a partir du sommet des thermiques.
 ! Dans le cas contraire, on demarre au niveau 1.
@@ -2814,5 +2962,5 @@
 ! FH 22/09/2009
 ! La ligne ci-dessous faisait osciller le modele et donnait une solution
-! asymptotique bidon et d\'ependant fortement du pas de temps.
+! assymptotique bidon et dÃ©pendant fortement du pas de temps.
 !        ratqs(:,:)=sqrt(ratqs(:,:)**2+ratqss(:,:)**2)
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -2842,4 +2990,8 @@
 c Appeler le processus de condensation a grande echelle
 c et le processus de precipitation
+c-------------------------------------------------------------------------
+      IF (prt_level .GE.10) THEN
+       print *,' ->fisrtilp '
+      ENDIF
 c-------------------------------------------------------------------------
       CALL fisrtilp(dtime,paprs,pplay,
@@ -3351,4 +3503,8 @@
        RCFC12 = RCFC12_act
 c
+      IF (prt_level .GE.10) THEN
+       print *,' ->radlwsw, number 1 '
+      ENDIF
+c
          CALL radlwsw
      e        (dist, rmu0, fract, 
@@ -3388,4 +3544,8 @@
        RCFC12 = RCFC12_per
 c
+      IF (prt_level .GE.10) THEN
+       print *,' ->radlwsw, number 2 '
+      ENDIF
+c
          CALL radlwsw
      e        (dist, rmu0, fract, 
@@ -3480,4 +3640,8 @@
 c a l'echelle sous-maille:
 c
+      IF (prt_level .GE.10) THEN
+       print *,' call orography ? ', ok_orodr
+      ENDIF
+c
       IF (ok_orodr) THEN
 c
@@ -3569,4 +3733,5 @@
 
        IF (ok_hines) then
+
          CALL hines_gwd(klon,klev,dtime,paprs,pplay,
      i                  rlat,t_seri,u_seri,v_seri,
@@ -3576,6 +3741,10 @@
 c  ajout des tendances
         CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,'hin')
+
       ENDIF
-
+c
+
+c
+cIM cf. FLott BEG
 C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
 
@@ -3602,6 +3771,5 @@
 cIM calcul composantes axiales du moment angulaire et couple des montagnes
 c
-      IF (is_sequential .and. ok_orodr) THEN
-      
+      IF (is_sequential .and. ok_orodr) THEN 
         CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur,
      C                 ra,rg,romega,
@@ -3898,4 +4066,8 @@
 c Convertir les incrementations en tendances
 c
+      IF (prt_level .GE.10) THEN
+        print *,'Convertir les incrementations en tendances '
+      ENDIF
+c
       if (mydebug) then
         call writefield_phy('u_seri',u_seri,llm)
@@ -4016,5 +4188,5 @@
 c=============================================================
 
-      if (iflag_thermals>1) then
+      if (iflag_thermals>=1) then
       d_t_lscth=0.
       d_t_lscst=0.
