Index: LMDZ5/trunk/libf/phylmd/calltherm.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/calltherm.F90	(revision 1635)
+++ LMDZ5/trunk/libf/phylmd/calltherm.F90	(revision 1638)
@@ -8,5 +8,14 @@
      &      ,fm_therm,entr_therm,detr_therm,zqasc,clwcon0,lmax,ratqscth,  &
      &       ratqsdiff,zqsatth,Ale_bl,Alp_bl,lalim_conv,wght_th, &
-     &       zmax0,f0,zw2,fraca,ztv,zpspsk,ztla,zthl)
+     &       zmax0,f0,zw2,fraca,ztv,zpspsk,ztla,zthl &
+!!! nrlmd le 10/04/2012
+     &      ,pbl_tke,pctsrf,omega,airephy &
+     &      ,zlcl_th,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
@@ -16,4 +25,8 @@
 #include "thermcell.h"
 #include "iniprint.h"
+
+!!! nrlmd le 10/04/2012
+#include "indicesol.h"
+!!! fin nrlmd le 10/04/2012
 
 !IM 140508
@@ -75,4 +88,18 @@
       !on garde le zmax du pas de temps precedent
       real zmax0(klon), f0(klon)
+
+!!! nrlmd le 10/04/2012
+      real pbl_tke(klon,klev+1,nbsrf)
+      real pctsrf(klon,nbsrf)
+      real omega(klon,klev)
+      real airephy(klon)
+      real zlcl_th(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)
+!!! fin nrlmd le 10/04/2012
+
 !********************************************************
 
@@ -220,5 +247,14 @@
      &      ,Ale,Alp,lalim_conv,wght_th &
      &      ,zmax0,f0,zw2,fraca,ztv,zpspsk &
-     &      ,ztla,zthl)
+     &      ,ztla,zthl &
+!!! nrlmd le 10/04/2012
+     &      ,pbl_tke,pctsrf,omega,airephy &
+     &      ,zlcl_th,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
+     &                         )
            if (prt_level.gt.10) write(lunout,*)'Apres thermcell_main OK'
          else
@@ -227,5 +263,20 @@
          endif
 
-       flag_bidouille_stratocu=iflag_thermals.eq.14.or.iflag_thermals.eq.16
+! Attention : les noms sont contre intuitif.
+! flag_bidouille_stratocu est .true. si on ne fait pas de bidouille.
+! Il aurait mieux valu avoir un nobidouille_stratocu
+! Et pour simplifier :
+! nobidouille_stratocu=.not.(iflag_thermals==13.or.iflag_thermals=15)
+! Ce serait bien de changer, mai en prenant le temps de vérifier que ca
+! fait bien ce qu'on croit.
+
+       flag_bidouille_stratocu=iflag_thermals<=12.or.iflag_thermals==14.or.iflag_thermals==16
+
+      if (iflag_thermals<=12) then
+         lmax=1
+         do k=1,klev-1
+            zdetr_therm(:,k)=zentr_therm(:,k)+zfm_therm(:,k)-zfm_therm(:,k+1)
+         enddo
+      endif
 
       fact(:)=0.
@@ -267,5 +318,4 @@
 
        DO i=1,klon
-        if(prt_level.GE.10) print*,'calltherm i Alp_bl Alp Ale_bl Ale',i,Alp_bl(i),Alp(i),Ale_bl(i),Ale(i)
             fm_therm(i,klev+1)=0.
             Ale_bl(i)=Ale_bl(i)+Ale(i)/REAL(nsplit_thermals)
@@ -273,4 +323,5 @@
             Alp_bl(i)=Alp_bl(i)+Alp(i)/REAL(nsplit_thermals)
 !            write(23,*)'ALP CALLTHERM',Alp_bl(i),Alp(i)
+        if(prt_level.GE.10) print*,'calltherm i Alp_bl Alp Ale_bl Ale',i,Alp_bl(i),Alp(i),Ale_bl(i),Ale(i)
        ENDDO
 
Index: LMDZ5/trunk/libf/phylmd/conf_phys.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/conf_phys.F90	(revision 1635)
+++ LMDZ5/trunk/libf/phylmd/conf_phys.F90	(revision 1638)
@@ -110,4 +110,9 @@
   integer,SAVE :: iflag_thermals_omp,nsplit_thermals_omp
   real,save :: tau_thermals_omp,alp_bl_k_omp
+!!! nrlmd le 10/04/2012
+  integer,SAVE :: iflag_trig_bl_omp,iflag_clos_bl_omp
+  integer,SAVE :: tau_trig_shallow_omp,tau_trig_deep_omp
+  real,SAVE    :: s_trig_omp
+!!! fin nrlmd le 10/04/2012
   real :: alp_offset
   REAL, SAVE :: alp_offset_omp
@@ -1083,4 +1088,48 @@
   alp_bl_k_omp = 1.
   call getin('alp_bl_k',alp_bl_k_omp)
+
+!!! nrlmd le 10/04/2012
+
+!Config Key  = iflag_trig_bl
+!Config Desc =  
+!Config Def  = 0
+!Config Help = 
+!
+  iflag_trig_bl_omp = 0
+  call getin('iflag_trig_bl',iflag_trig_bl_omp)
+
+!Config Key  = s_trig_bl
+!Config Desc =  
+!Config Def  = 0
+!Config Help = 
+!
+  s_trig_omp = 2e7
+  call getin('s_trig',s_trig_omp)
+
+!Config Key  = tau_trig_shallow
+!Config Desc =  
+!Config Def  = 0
+!Config Help = 
+!
+  tau_trig_shallow_omp = 600
+  call getin('tau_trig_shallow',tau_trig_shallow_omp)
+
+!Config Key  = tau_trig_deep
+!Config Desc =  
+!Config Def  = 0
+!Config Help = 
+!
+  tau_trig_deep_omp = 1800
+  call getin('tau_trig_deep',tau_trig_deep_omp)
+
+!Config Key  = iflag_clos_bl
+!Config Desc =  
+!Config Def  = 0
+!Config Help = 
+!
+  iflag_clos_bl_omp = 0
+  call getin('iflag_clos_bl',iflag_clos_bl_omp)
+
+!!! fin nrlmd le 10/04/2012
 
 !
@@ -1650,4 +1699,11 @@
     tau_thermals = tau_thermals_omp
     alp_bl_k = alp_bl_k_omp
+!!! nrlmd le 10/04/2012
+    iflag_trig_bl = iflag_trig_bl_omp
+    s_trig = s_trig_omp
+    tau_trig_shallow = tau_trig_shallow_omp
+    tau_trig_deep = tau_trig_deep_omp
+    iflag_clos_bl = iflag_clos_bl_omp
+!!! fin nrlmd le 10/04/2012
     iflag_coupl = iflag_coupl_omp
     iflag_clos = iflag_clos_omp
@@ -1839,4 +1895,11 @@
   write(lunout,*)' iflag_wake = ', iflag_wake
   write(lunout,*)' alp_offset = ', alp_offset
+!!! nrlmd le 10/04/2012
+  write(lunout,*)' iflag_trig_bl = ', iflag_trig_bl
+  write(lunout,*)' s_trig = ', s_trig
+  write(lunout,*)' tau_trig_shallow = ', tau_trig_shallow
+  write(lunout,*)' tau_trig_deep = ', tau_trig_deep
+  write(lunout,*)' iflag_clos_bl = ', iflag_clos_bl
+!!! fin nrlmd le 10/04/2012
 
   write(lunout,*)' lonmin lonmax latmin latmax bilKP_ins =',&
Index: LMDZ5/trunk/libf/phylmd/phys_output_mod.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/phys_output_mod.F90	(revision 1635)
+++ LMDZ5/trunk/libf/phylmd/phys_output_mod.F90	(revision 1638)
@@ -260,4 +260,36 @@
   type(ctrl_out),save :: o_wape         = ctrl_out((/ 1, 1, 1, 10, 10, 10 /),'wape')
 
+!!! nrlmd le 10/04/2012
+
+!-------Propriétés du thermiques au LCL 
+  type(ctrl_out),save :: o_zlcl_th        = ctrl_out((/ 1, 1, 1, 6, 10, 10 /),'zlcl_th')
+  type(ctrl_out),save :: o_fraca0         = ctrl_out((/ 1, 1, 1, 6, 10, 10 /),'fraca0')
+  type(ctrl_out),save :: o_w0             = ctrl_out((/ 1, 1, 1, 6, 10, 10 /),'w0')
+  type(ctrl_out),save :: o_w_conv         = ctrl_out((/ 1, 1, 1, 6, 10, 10 /),'w_conv')
+  type(ctrl_out),save :: o_therm_tke_max0 = ctrl_out((/ 1, 1, 1, 6, 10, 10 /),'therm_tke_max0')
+  type(ctrl_out),save :: o_env_tke_max0   = ctrl_out((/ 1, 1, 1, 6, 10, 10 /),'env_tke_max0')
+
+!-------Spectre de thermiques de type 2 au LCL
+  type(ctrl_out),save :: o_n2                = ctrl_out((/ 1, 1, 1, 6, 10, 10 /),'n2')
+  type(ctrl_out),save :: o_s2                = ctrl_out((/ 1, 1, 1, 6, 10, 10 /),'s2')
+                                                                             
+!-------Déclenchement stochastique                                           
+  type(ctrl_out),save :: o_proba_notrig      = ctrl_out((/ 1, 1, 1, 6, 10, 10 /),'proba_notrig')
+  type(ctrl_out),save :: o_random_notrig     = ctrl_out((/ 1, 1, 1, 6, 10, 10 /),'random_notrig')
+  type(ctrl_out),save :: o_ale_bl_stat       = ctrl_out((/ 1, 1, 1, 6, 10, 10 /),'ale_bl_stat')
+  type(ctrl_out),save :: o_ale_bl_trig       = ctrl_out((/ 1, 1, 1, 6, 10, 10 /),'ale_bl_trig')
+                                                                             
+!-------Profils de TKE dans et hors du thermique                             
+  type(ctrl_out),save :: o_therm_tke_max     = ctrl_out((/ 1, 1, 1, 6, 10, 10 /),'therm_tke_max')
+  type(ctrl_out),save :: o_env_tke_max       = ctrl_out((/ 1, 1, 1, 6, 10, 10 /),'env_tke_max')
+
+!-------Fermeture statistique
+  type(ctrl_out),save :: o_alp_bl_det        = ctrl_out((/ 1, 1, 1, 10, 10, 10 /),'alp_bl_det')
+  type(ctrl_out),save :: o_alp_bl_fluct_m    = ctrl_out((/ 1, 1, 1, 10, 10, 10 /),'alp_bl_fluct_m')
+  type(ctrl_out),save :: o_alp_bl_fluct_tke  = ctrl_out((/ 1, 1, 1, 10, 10, 10 /),'alp_bl_fluct_tke')
+  type(ctrl_out),save :: o_alp_bl_conv       = ctrl_out((/ 1, 1, 1, 10, 10, 10 /),'alp_bl_conv')
+  type(ctrl_out),save :: o_alp_bl_stat       = ctrl_out((/ 1, 1, 1, 10, 10, 10 /),'alp_bl_stat')
+
+!!! fin nrlmd le 10/04/2012
 
   ! Champs interpolles sur des niveaux de pression ??? a faire correctement
@@ -432,4 +464,5 @@
   type(ctrl_out),save :: o_ovap         = ctrl_out((/ 2, 3, 4, 10, 10, 10 /),'ovap')
   type(ctrl_out),save :: o_ovapinit     = ctrl_out((/ 2, 10, 10, 10, 10, 10 /),'ovapinit')
+  type(ctrl_out),save :: o_oliq         = ctrl_out((/ 2, 3, 4, 10, 10, 10 /),'oliq')
   type(ctrl_out),save :: o_wvapp        = ctrl_out((/ 2, 10, 10, 10, 10, 10 /),'wvapp')
   type(ctrl_out),save :: o_geop         = ctrl_out((/ 2, 3, 10, 10, 10, 10 /),'geop')
@@ -494,4 +527,5 @@
   type(ctrl_out),save :: o_dtcon        = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'dtcon')
   type(ctrl_out),save :: o_ducon        = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'ducon')
+  type(ctrl_out),save :: o_dvcon        = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'dvcon')
   type(ctrl_out),save :: o_dqcon        = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'dqcon')
   type(ctrl_out),save :: o_dtwak        = ctrl_out((/ 4, 5, 10, 10, 10, 10 /),'dtwak')
@@ -531,5 +565,4 @@
   type(ctrl_out),save :: o_e_th         = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'e_th')
   type(ctrl_out),save :: o_w_th         = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'w_th')
-  type(ctrl_out),save :: o_lambda_th    = ctrl_out((/ 10, 10, 10, 10, 10, 10 /),'lambda_th')
   type(ctrl_out),save :: o_ftime_th     = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'ftime_th')
   type(ctrl_out),save :: o_q_th         = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'q_th')
@@ -537,5 +570,5 @@
   type(ctrl_out),save :: o_d_th         = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'d_th')
   type(ctrl_out),save :: o_f0_th        = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'f0_th')
-  type(ctrl_out),save :: o_zmax_th      = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'zmax_th')
+  type(ctrl_out),save :: o_zmax_th      = ctrl_out((/ 4,  4,  4,  5, 10, 10 /),'zmax_th')
   type(ctrl_out),save :: o_dqthe        = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'dqthe')
   type(ctrl_out),save :: o_dtajs        = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'dtajs')
@@ -621,4 +654,5 @@
     USE infotrac
     USE ioipsl
+!    USE phys_cal_mod, only : hour
     USE mod_phys_lmdz_para
     USE aero_mod, only : naero_spc,name_aero
@@ -682,9 +716,9 @@
     !                 entre [phys_out_lonmin,phys_out_lonmax] et [phys_out_latmin,phys_out_latmax]
 
-    logical, dimension(nfiles), save  :: phys_out_regfkey       = (/ .false., .false., .false., .false., .false., .false. /)
-    real, dimension(nfiles), save     :: phys_out_lonmin        = (/ -180., -180., -180., -180., -180., -180. /)
-    real, dimension(nfiles), save     :: phys_out_lonmax        = (/ 180., 180., 180., 180., 180., 180. /)
-    real, dimension(nfiles), save     :: phys_out_latmin        = (/ -90., -90., -90., -90., -90., -90. /)
-    real, dimension(nfiles), save     :: phys_out_latmax        = (/ 90., 90., 90., 90., 90., 90. /)
+    logical, dimension(nfiles), save  :: phys_out_regfkey       = (/ .false., .false., .false.,  .false., .false., .false. /)
+    real, dimension(nfiles), save     :: phys_out_lonmin        = (/   -180.,   -180.,   -180.,    -180.,   -180.,   -180. /)
+    real, dimension(nfiles), save     :: phys_out_lonmax        = (/    180.,    180.,    180.,     180.,    180.,    180. /)
+    real, dimension(nfiles), save     :: phys_out_latmin        = (/    -90.,    -90.,    -90.,     -90.,    -90.,    -90. /)
+    real, dimension(nfiles), save     :: phys_out_latmax        = (/     90.,     90.,     90.,     90.,     90.,     90. /)
 
     write(lunout,*) 'Debut phys_output_mod.F90'
@@ -805,5 +839,8 @@
 
           idayref = day_ref
-          CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
+          CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)       
+! correction pour l heure initiale                               !jyg 
+!                                                                !jyg 
+!      CALL ymds2ju(annee_ref, 1, idayref, hour, zjulian)         !jyg
 
 !!!!!!!!!!!!!!!!! Traitement dans le cas ou l'on veut stocker sur un domaine limite !!
@@ -1318,8 +1355,10 @@
           ! Couplage conv-CL
           IF (iflag_con.GE.3) THEN
+             IF (iflag_coupl>=1) THEN
                 CALL histdef2d(iff,clef_stations(iff), &
                      o_ale_bl%flag,o_ale_bl%name, "ALE BL", "m2/s2") 
                 CALL histdef2d(iff,clef_stations(iff), &
                      o_alp_bl%flag,o_alp_bl%name, "ALP BL", "m2/s2") 
+             ENDIF
           ENDIF !(iflag_con.GE.3)
 
@@ -1375,4 +1414,5 @@
           CALL histdef3d(iff,clef_stations(iff),o_theta%flag,o_theta%name, "Potential air temperature", "K" )
           CALL histdef3d(iff,clef_stations(iff),o_ovap%flag,o_ovap%name, "Specific humidity", "kg/kg" )
+          CALL histdef3d(iff,clef_stations(iff),o_oliq%flag,o_oliq%name, "Condensed water", "kg/kg" )
           CALL histdef3d(iff,clef_stations(iff), &
                o_ovapinit%flag,o_ovapinit%name, "Specific humidity (begin of timestep)", "kg/kg" )
@@ -1480,4 +1520,6 @@
                o_ducon%flag,o_ducon%name, "Convection du", "m/s2")
           CALL histdef3d(iff,clef_stations(iff), &
+               o_dvcon%flag,o_dvcon%name, "Convection dv", "m/s2")
+          CALL histdef3d(iff,clef_stations(iff), &
                o_dqcon%flag,o_dqcon%name, "Convection dQ", "(kg/kg)/s")
 
@@ -1489,4 +1531,10 @@
                 CALL histdef2d(iff,clef_stations(iff), &
                      o_alp_wk%flag,o_alp_wk%name, "ALP WK", "m2/s2")
+                CALL histdef2d(iff,clef_stations(iff), &
+                     o_ale%flag,o_ale%name, "ALE", "m2/s2")
+                CALL histdef2d(iff,clef_stations(iff), &
+                     o_alp%flag,o_alp%name, "ALP", "W/m2")
+                CALL histdef2d(iff,clef_stations(iff),o_cin%flag,o_cin%name, "Convective INhibition", "m2/s2")
+                CALL histdef2d(iff,clef_stations(iff),o_wape%flag,o_WAPE%name, "WAPE", "m2/s2")
                 CALL histdef2d(iff,clef_stations(iff),o_wake_h%flag,o_wake_h%name, "wake_h", "-")
                 CALL histdef2d(iff,clef_stations(iff),o_wake_s%flag,o_wake_s%name, "wake_s", "-")
@@ -1496,15 +1544,36 @@
                 CALL histdef3d(iff,clef_stations(iff),o_wake_deltaq%flag,o_wake_deltaq%name, "wake_deltaq", " ")
                 CALL histdef3d(iff,clef_stations(iff),o_wake_omg%flag,o_wake_omg%name, "wake_omg", " ")
-                CALL histdef2d(iff,clef_stations(iff),o_wape%flag,o_WAPE%name, "WAPE", "m2/s2")
              ENDIF
-             CALL histdef2d(iff,clef_stations(iff), &
-                     o_ale%flag,o_ale%name, "ALE", "m2/s2")
-             CALL histdef2d(iff,clef_stations(iff), &
-                     o_alp%flag,o_alp%name, "ALP", "W/m2")
-             CALL histdef2d(iff,clef_stations(iff),o_cin%flag,o_cin%name, "Convective INhibition", "m2/s2")
              CALL histdef3d(iff,clef_stations(iff),o_Vprecip%flag,o_Vprecip%name, "precipitation vertical profile", "-")
              CALL histdef3d(iff,clef_stations(iff),o_ftd%flag,o_ftd%name, "tend temp due aux descentes precip", "-")
              CALL histdef3d(iff,clef_stations(iff),o_fqd%flag,o_fqd%name,"tend vap eau due aux descentes precip", "-")
           ENDIF !(iflag_con.EQ.3)
+
+!!! nrlmd le 10/04/2012
+
+ CALL histdef2d(iff,clef_stations(iff),o_zlcl_th%flag,o_zlcl_th%name, "Altitude du LCL", "m")
+ CALL histdef2d(iff,clef_stations(iff),o_fraca0%flag,o_fraca0%name, "Fraction de Thermique au LCL", "")
+ CALL histdef2d(iff,clef_stations(iff),o_w0%flag,o_w0%name, "Vitesse thermique au LCL", "m/s")
+ CALL histdef2d(iff,clef_stations(iff),o_w_conv%flag,o_w_conv%name, "Vitesse verticale grande échelle au LCL", "m/s")
+ CALL histdef2d(iff,clef_stations(iff),o_therm_tke_max0%flag,o_therm_tke_max0%name, "TKE thermique au LCL", "m2/s2")
+ CALL histdef2d(iff,clef_stations(iff),o_env_tke_max0%flag,o_env_tke_max0%name, "TKE environnement au LCL", "m2/s2")
+
+ CALL histdef2d(iff,clef_stations(iff),o_n2%flag,o_n2%name, "Nombre de panaches de type 2", " ")
+ CALL histdef2d(iff,clef_stations(iff),o_s2%flag,o_s2%name, "Surface moyenne des panaches de type 2", "m2")
+
+ CALL histdef2d(iff,clef_stations(iff),o_proba_notrig%flag,o_proba_notrig%name, "Probabilité de non-déclenchement", " ")
+ CALL histdef2d(iff,clef_stations(iff),o_random_notrig%flag,o_random_notrig%name, "Tirage aléatoire de non-déclenchement", " ")
+ CALL histdef2d(iff,clef_stations(iff),o_ale_bl_trig%flag,o_ale_bl_trig%name, "ALE_BL_STAT + Condition P>Pseuil", "m2/s2")
+ CALL histdef2d(iff,clef_stations(iff),o_ale_bl_stat%flag,o_ale_bl_stat%name, "ALE_BL_STAT", "m2/s2")
+
+ CALL histdef3d(iff,clef_stations(iff),o_therm_tke_max%flag,o_therm_tke_max%name, "TKE thermique", "m2/s2")
+ CALL histdef3d(iff,clef_stations(iff),o_env_tke_max%flag,o_env_tke_max%name, "TKE environnement", "m2/s2")
+ CALL histdef2d(iff,clef_stations(iff),o_alp_bl_det%flag,o_alp_bl_det%name, "ALP_BL_DET", "W/m2")
+ CALL histdef2d(iff,clef_stations(iff),o_alp_bl_fluct_m%flag,o_alp_bl_fluct_m%name, "ALP_BL_FLUCT_M", "W/m2")
+ CALL histdef2d(iff,clef_stations(iff),o_alp_bl_fluct_tke%flag,o_alp_bl_fluct_tke%name, "ALP_BL_FLUCT_TKE", "W/m2")
+ CALL histdef2d(iff,clef_stations(iff),o_alp_bl_conv%flag,o_alp_bl_conv%name, "ALP_BL_CONV", "W/m2")
+ CALL histdef2d(iff,clef_stations(iff),o_alp_bl_stat%flag,o_alp_bl_stat%name, "ALP_BL_STAT", "W/m2")
+
+!!! fin nrlmd le 10/04/2012
 
           CALL histdef3d(iff,clef_stations(iff),o_dtlsc%flag,o_dtlsc%name, "Condensation dT", "K/s")
@@ -1519,5 +1588,5 @@
           CALL histdef3d(iff,clef_stations(iff),o_dtthe%flag,o_dtthe%name, "Thermal dT", "K/s")
 
-          if(iflag_thermals.gt.1) THEN
+          if(iflag_thermals.ge.1) THEN
              CALL histdef3d(iff,clef_stations(iff),o_dqlscth%flag,o_dqlscth%name, "dQ therm.", "(kg/kg)/s")
              CALL histdef3d(iff,clef_stations(iff),o_dqlscst%flag,o_dqlscst%name, "dQ strat.", "(kg/kg)/s")
@@ -1531,6 +1600,4 @@
              CALL histdef3d(iff,clef_stations(iff),o_e_th%flag,o_e_th%name,"Thermal plume entrainment","K/s")
              CALL histdef3d(iff,clef_stations(iff),o_w_th%flag,o_w_th%name,"Thermal plume vertical velocity","m/s")
-             CALL histdef3d(iff,clef_stations(iff), &
-                  o_lambda_th%flag,o_lambda_th%name,"Thermal plume vertical velocity","m/s")
              CALL histdef2d(iff,clef_stations(iff), &
                   o_ftime_th%flag,o_ftime_th%name,"Fraction of time Shallow convection occurs"," ")
@@ -1548,5 +1615,5 @@
              CALL histdef3d(iff,clef_stations(iff), &
                   o_dqthe%flag,o_dqthe%name, "Thermal dQ", "(kg/kg)/s")
-          endif !iflag_thermals.gt.1
+          endif !iflag_thermals.ge.1
           CALL histdef3d(iff,clef_stations(iff), &
                o_dtajs%flag,o_dtajs%name, "Dry adjust. dT", "K/s")
Index: LMDZ5/trunk/libf/phylmd/phys_output_write.h
===================================================================
--- LMDZ5/trunk/libf/phylmd/phys_output_write.h	(revision 1635)
+++ LMDZ5/trunk/libf/phylmd/phys_output_write.h	(revision 1638)
@@ -508,5 +508,5 @@
         ENDIF
 
-       if (iflag_pbl>1 .and. lev_histday.gt.10 ) then
+       if (iflag_pbl>1 .and. lev_files(iff).gt.10 ) then
         IF (o_tke_srf(nsrf)%flag(iff)<=lev_files(iff)) THEN
         CALL histwrite_phy(nid_files(iff),clef_stations(iff),
@@ -666,5 +666,5 @@
 
        IF (o_mc%flag(iff)<=lev_files(iff)) THEN
-        if(iflag_thermals.gt.1)then
+        if(iflag_thermals>=1)then
          zx_tmp_fi3d=dnwd+dnwd0+upwd+fm_therm(:,1:klev)
         else
@@ -677,7 +677,8 @@
       ENDIF !iflag_con .GE. 3
 
-      IF (o_prw%flag(iff) <= lev_files(iff))
-     $     CALL histwrite_phy(nid_files(iff), clef_stations(iff),
-     $     o_prw%name, itau_w, prw)
+        IF (o_prw%flag(iff)<=lev_files(iff)) THEN
+      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     $o_prw%name,itau_w,prw)
+        ENDIF
 
         IF (o_s_pblh%flag(iff)<=lev_files(iff)) THEN
@@ -800,4 +801,5 @@
 ! Couplage convection-couche limite
       IF (iflag_con.GE.3) THEN
+      IF (iflag_coupl>=1) THEN
        IF (o_ale_bl%flag(iff)<=lev_files(iff)) THEN
        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
@@ -808,4 +810,5 @@
      $o_alp_bl%name,itau_w,alp_bl)
        ENDIF
+      ENDIF !iflag_coupl>=1
       ENDIF !(iflag_con.GE.3)
 
@@ -822,4 +825,16 @@
        ENDIF
 
+       IF (o_ale%flag(iff)<=lev_files(iff)) THEN
+       CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     $o_ale%name,itau_w,ale)
+       ENDIF
+       IF (o_alp%flag(iff)<=lev_files(iff)) THEN
+       CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     $o_alp%name,itau_w,alp)
+       ENDIF
+       IF (o_cin%flag(iff)<=lev_files(iff)) THEN
+       CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     $o_cin%name,itau_w,cin)
+       ENDIF
        IF (o_wape%flag(iff)<=lev_files(iff)) THEN
        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
@@ -868,16 +883,4 @@
       ENDIF ! iflag_wake>=1
 
-       IF (o_ale%flag(iff)<=lev_files(iff)) THEN
-       CALL histwrite_phy(nid_files(iff),clef_stations(iff),
-     $o_ale%name,itau_w,ale)
-       ENDIF
-       IF (o_alp%flag(iff)<=lev_files(iff)) THEN
-       CALL histwrite_phy(nid_files(iff),clef_stations(iff),
-     $o_alp%name,itau_w,alp)
-       ENDIF
-       IF (o_cin%flag(iff)<=lev_files(iff)) THEN
-       CALL histwrite_phy(nid_files(iff),clef_stations(iff),
-     $o_cin%name,itau_w,cin)
-       ENDIF
         IF (o_Vprecip%flag(iff)<=lev_files(iff)) THEN
        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
@@ -896,4 +899,103 @@
       ENDIF !(iflag_con.EQ.3) 
  
+!!! nrlmd le 10/04/2012
+
+        IF (o_zlcl_th%flag(iff)<=lev_files(iff)) THEN
+        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     s                     o_zlcl_th%name,itau_w,zlcl_th)
+        ENDIF
+
+        IF (o_fraca0%flag(iff)<=lev_files(iff)) THEN
+        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     s                     o_fraca0%name,itau_w,fraca0)
+        ENDIF
+
+        IF (o_w0%flag(iff)<=lev_files(iff)) THEN
+        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     s                     o_w0%name,itau_w,w0)
+        ENDIF
+
+        IF (o_w_conv%flag(iff)<=lev_files(iff)) THEN
+        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     s                     o_w_conv%name,itau_w,w_conv)
+        ENDIF
+
+        IF (o_therm_tke_max0%flag(iff)<=lev_files(iff)) THEN
+        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     s                     o_therm_tke_max0%name,itau_w,therm_tke_max0)
+        ENDIF
+
+        IF (o_env_tke_max0%flag(iff)<=lev_files(iff)) THEN
+        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     s                     o_env_tke_max0%name,itau_w,env_tke_max0)
+        ENDIF
+
+        IF (o_n2%flag(iff)<=lev_files(iff)) THEN
+      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     s                     o_n2%name,itau_w,n2)
+       ENDIF
+
+       IF (o_s2%flag(iff)<=lev_files(iff)) THEN
+      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     s                     o_s2%name,itau_w,s2)
+       ENDIF
+
+        IF (o_proba_notrig%flag(iff)<=lev_files(iff)) THEN
+      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     s                     o_proba_notrig%name,itau_w,proba_notrig)
+       ENDIF
+
+       IF (o_random_notrig%flag(iff)<=lev_files(iff)) THEN
+      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     s                     o_random_notrig%name,itau_w,random_notrig)
+       ENDIF
+
+       IF (o_ale_bl_stat%flag(iff)<=lev_files(iff)) THEN
+      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     s                     o_ale_bl_stat%name,itau_w,ale_bl_stat)
+       ENDIF
+
+       IF (o_ale_bl_trig%flag(iff)<=lev_files(iff)) THEN
+      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     s                     o_ale_bl_trig%name,itau_w,ale_bl_trig)
+       ENDIF
+
+       IF (o_therm_tke_max%flag(iff)<=lev_files(iff)) THEN
+        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     s                     o_therm_tke_max%name,itau_w,therm_tke_max)
+        ENDIF
+
+        IF (o_env_tke_max%flag(iff)<=lev_files(iff)) THEN
+        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     s                     o_env_tke_max%name,itau_w,env_tke_max)
+        ENDIF
+
+       IF (o_alp_bl_det%flag(iff)<=lev_files(iff)) THEN
+      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     s                     o_alp_bl_det%name,itau_w,alp_bl_det)
+       ENDIF
+
+       IF (o_alp_bl_fluct_m%flag(iff)<=lev_files(iff)) THEN
+      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     s                     o_alp_bl_fluct_m%name,itau_w,alp_bl_fluct_m)
+       ENDIF
+
+       IF (o_alp_bl_fluct_tke%flag(iff)<=lev_files(iff)) THEN
+      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     s                o_alp_bl_fluct_tke%name,itau_w,alp_bl_fluct_tke)
+       ENDIF
+
+       IF (o_alp_bl_conv%flag(iff)<=lev_files(iff)) THEN
+      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     s                     o_alp_bl_conv%name,itau_w,alp_bl_conv)
+       ENDIF
+
+       IF (o_alp_bl_stat%flag(iff)<=lev_files(iff)) THEN
+      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     s                     o_alp_bl_stat%name,itau_w,alp_bl_stat)
+       ENDIF
+
+!!! fin nrlmd le 10/04/2012
+
       IF (type_ocean=='slab ') THEN
       IF ( o_slab_bils%flag(iff)<=lev_files(iff)) 
@@ -1409,4 +1511,9 @@
       CALL histwrite_phy(nid_files(iff),clef_stations(iff),
      $                   o_ovap%name,itau_w,q_seri)
+       ENDIF
+                                                               
+       IF (o_oliq%flag(iff)<=lev_files(iff)) THEN              
+      CALL histwrite_phy(nid_files(iff),clef_stations(iff),    
+     $                   o_oliq%name,itau_w,ql_seri)            
        ENDIF
 
@@ -1636,4 +1743,10 @@
        ENDIF
 
+       IF (o_dvcon%flag(iff)<=lev_files(iff)) THEN
+      zx_tmp_fi3d(1:klon,1:klev)=d_v_con(1:klon,1:klev)/pdtphys
+      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     $o_dvcon%name,itau_w,zx_tmp_fi3d)
+       ENDIF
+
        IF (o_dqcon%flag(iff)<=lev_files(iff)) THEN
       zx_tmp_fi3d(1:klon,1:klev)=d_q_con(1:klon,1:klev)/pdtphys
@@ -1679,5 +1792,5 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 ! Sorties specifiques a la separation thermiques/non thermiques
-       if (iflag_thermals>1) then
+       if (iflag_thermals>=1) then
 
        IF (o_dtlscth%flag(iff)<=lev_files(iff)) THEN
@@ -1743,5 +1856,5 @@
        ENDIF
 
-      endif ! iflag_thermals>1
+      endif ! iflag_thermals>=1
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -1790,5 +1903,5 @@
        ENDIF
 
-       IF (iflag_thermals.gt.1) THEN
+       IF (iflag_thermals>=1) THEN
         IF (o_ftime_th%flag(iff)<=lev_files(iff)) THEN
 ! Pour l instant 0 a y reflichir pour les thermiques
@@ -1819,8 +1932,4 @@
         ENDIF
 
-        IF (o_lambda_th%flag(iff)<=lev_files(iff)) THEN
-        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
-     s                     o_lambda_th%name,itau_w,lambda_th)
-        ENDIF
 
         IF (o_a_th%flag(iff)<=lev_files(iff)) THEN
Index: LMDZ5/trunk/libf/phylmd/phys_state_var_mod.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/phys_state_var_mod.F90	(revision 1635)
+++ LMDZ5/trunk/libf/phylmd/phys_state_var_mod.F90	(revision 1638)
@@ -346,4 +346,9 @@
 !$OMP THREADPRIVATE(ccm)
 
+!!! nrlmd le 10/04/2012
+      REAL,SAVE,ALLOCATABLE :: ale_bl_trig(:)
+!$OMP THREADPRIVATE(ale_bl_trig)
+!!! fin nrlmd le 10/04/2012
+
 CONTAINS
 
@@ -496,4 +501,8 @@
       ALLOCATE(tau_aero(klon,klev,naero_grp,nbands),piz_aero(klon,klev,naero_grp,nbands),cg_aero(klon,klev,naero_grp,nbands))
       ALLOCATE(ccm(klon,klev,nbands))
+
+!!! nrlmd le 10/04/2012
+      ALLOCATE(ale_bl_trig(klon))
+!!! fin nrlmd le 10/04/2012
 
 END SUBROUTINE phys_state_var_init
@@ -603,4 +612,8 @@
       deallocate(ccm)
        
+!!! nrlmd le 10/04/2012
+      deallocate(ale_bl_trig)
+!!! fin nrlmd le 10/04/2012
+
 END SUBROUTINE phys_state_var_end
 
Index: LMDZ5/trunk/libf/phylmd/physiq.F
===================================================================
--- LMDZ5/trunk/libf/phylmd/physiq.F	(revision 1635)
+++ 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.
Index: LMDZ5/trunk/libf/phylmd/thermcell.h
===================================================================
--- LMDZ5/trunk/libf/phylmd/thermcell.h	(revision 1635)
+++ LMDZ5/trunk/libf/phylmd/thermcell.h	(revision 1638)
@@ -1,7 +1,14 @@
       integer            :: iflag_thermals,nsplit_thermals
+
+!!! nrlmd le 10/04/2012
+      integer            :: iflag_trig_bl,iflag_clos_bl
+      integer            :: tau_trig_shallow,tau_trig_deep
+      real               :: s_trig
+!!! fin nrlmd le 10/04/2012
+
       real,parameter     :: r_aspect_thermals=2.,l_mix_thermals=30.
       real               :: alp_bl_k
       real               :: tau_thermals
-      integer,parameter  :: w2di_thermals=1
+      integer,parameter  :: w2di_thermals=0
       integer            :: isplit
 
@@ -14,3 +21,9 @@
       common/ctherm5/iflag_thermals_ed,iflag_thermals_optflux
 
-!$OMP THREADPRIVATE(/ctherm1/,/ctherm2/,/ctherm4/,/ctherm5/)
+!!! nrlmd le 10/04/2012
+      common/ctherm6/iflag_trig_bl,iflag_clos_bl
+      common/ctherm7/tau_trig_shallow,tau_trig_deep
+      common/ctherm8/s_trig
+!!! fin nrlmd le 10/04/2012
+
+!$OMP THREADPRIVATE(/ctherm1/,/ctherm2/,/ctherm4/,/ctherm5/,/ctherm6/,/ctherm7/,,/ctherm6/,/ctherm8/)
Index: LMDZ5/trunk/libf/phylmd/thermcell_main.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/thermcell_main.F90	(revision 1635)
+++ LMDZ5/trunk/libf/phylmd/thermcell_main.F90	(revision 1638)
@@ -10,5 +10,14 @@
      &                  ,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,4 +56,7 @@
 #include "iniprint.h"
 #include "thermcell.h"
+!!! nrlmd le 10/04/2012
+#include "indicesol.h"
+!!! fin nrlmd le 10/04/2012
 
 !   arguments:
@@ -77,4 +89,6 @@
       integer,save :: lev_out=10
 !$OMP THREADPRIVATE(lev_out)
+
+      REAL susqr2pi, Reuler
 
       INTEGER ig,k,l,ll,ierr
@@ -155,4 +169,51 @@
        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,4 +740,205 @@
       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,2 +1120,133 @@
       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
+
