Index: LMDZ5/branches/testing/libf/phylmd/calltherm.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/calltherm.F90	(revision 1665)
+++ LMDZ5/branches/testing/libf/phylmd/calltherm.F90	(revision 1669)
@@ -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/branches/testing/libf/phylmd/concvl.F
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/concvl.F	(revision 1665)
+++ LMDZ5/branches/testing/libf/phylmd/concvl.F	(revision 1669)
@@ -248,5 +248,5 @@
          DO i = 1, klon
           cbmf(i) = 0.
-!          plcl(i) = 0.
+          plcl(i) = 0.
           sigd(i) = 0.
          ENDDO
@@ -256,5 +256,4 @@
       plfc(:)  = 0.
       wbeff(:) = 100.
-      plcl(:) = 0.
 
       DO k = 1, klev+1
@@ -369,5 +368,8 @@
      $              cape,cin,tvp,
      $              dd_t,dd_q,Plim1,Plim2,asupmax,supmax0,
-     $              asupmaxmin,lalim_conv)
+     $              asupmaxmin,lalim_conv,
+!AC!
+     $              da,phi)
+!AC!
       endif  
 C------------------------------------------------------------------
@@ -399,4 +401,16 @@
        ENDDO 
        endif
+
+c!AC!
+       if (iflag_con.eq.3) then
+       DO itra = 1,ntra
+        DO k = 1, klev
+         DO i = 1, klon
+            d_tra(i,k,itra) =dtime*d_tra(i,k,itra) 
+         ENDDO
+        ENDDO
+       ENDDO 
+       endif
+c!AC!
 
       DO k = 1, klev
Index: LMDZ5/branches/testing/libf/phylmd/conf_phys.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/conf_phys.F90	(revision 1665)
+++ LMDZ5/branches/testing/libf/phylmd/conf_phys.F90	(revision 1669)
@@ -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
@@ -282,4 +287,5 @@
 !Config Help = Used in physiq.F
 !
+! - flag_aerosol=0 => no aerosol
 ! - flag_aerosol=1 => so4 only (defaut) 
 ! - flag_aerosol=2 => bc  only 
@@ -289,5 +295,5 @@
 ! - flag_aerosol=6 => all aerosol
 
-  flag_aerosol_omp = 1
+  flag_aerosol_omp = 0
   CALL getin('flag_aerosol',flag_aerosol_omp)
 
@@ -1083,4 +1089,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 +1700,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
@@ -1710,4 +1767,7 @@
 ! il n'est utilisable que lors du couplage avec le SO4 seul 
     IF (ok_ade .OR. ok_aie) THEN 
+       IF ( flag_aerosol .EQ. 0 ) THEN
+          CALL abort_gcm('conf_phys','flag_aerosol=0 not compatible avec ok_ade ou ok_aie=.TRUE.',1)
+       END IF
        IF ( .NOT. new_aod .AND.  flag_aerosol .NE. 1) THEN
           CALL abort_gcm('conf_phys','new_aod=.FALSE. not compatible avec flag_aerosol=1',1)
@@ -1839,4 +1899,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/branches/testing/libf/phylmd/cv3_routines.F
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/cv3_routines.F	(revision 1665)
+++ LMDZ5/branches/testing/libf/phylmd/cv3_routines.F	(revision 1669)
@@ -879,16 +879,16 @@
  110  continue
 
-      do 121 j=1,ntra
-ccccc      do 111 k=1,nl+1
-      do 111 k=1,nd
-       nn=0
-      do 101 i=1,len
-      if(iflag1(i).eq.0)then
-       nn=nn+1
-       tra(nn,k,j)=tra1(i,k,j)
-      endif
- 101  continue
- 111  continue
- 121  continue
+!AC!      do 121 j=1,ntra
+!AC!ccccc      do 111 k=1,nl+1
+!AC!      do 111 k=1,nd
+!AC!       nn=0
+!AC!      do 101 i=1,len
+!AC!      if(iflag1(i).eq.0)then
+!AC!       nn=nn+1
+!AC!       tra(nn,k,j)=tra1(i,k,j)
+!AC!      endif
+!AC! 101  continue
+!AC! 111  continue
+!AC! 121  continue
 
       if (nn.ne.ncum) then
@@ -1633,13 +1633,13 @@
       sij(1:ncum,1:nd,1:nd)=0.0
       
-      do k=1,ntra
-       do j=1,nd  ! instead nlp
-        do i=1,nd ! instead nlp
-         do il=1,ncum
-            traent(il,i,j,k)=tra(il,j,k)
-         enddo
-        enddo
-       enddo
-      enddo
+!AC!      do k=1,ntra
+!AC!       do j=1,nd  ! instead nlp
+!AC!        do i=1,nd ! instead nlp
+!AC!         do il=1,ncum
+!AC!            traent(il,i,j,k)=tra(il,j,k)
+!AC!         enddo
+!AC!        enddo
+!AC!       enddo
+!AC!      enddo
       zm(:,:)=0.
 
@@ -1697,15 +1697,15 @@
  710  continue
 
-       do k=1,ntra
-        do j=minorig,nl
-         do il=1,ncum
-          if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
-     :       (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
-            traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
-     :            +(1.-sij(il,i,j))*tra(il,nk(il),k)
-          endif
-         enddo
-        enddo
-       enddo
+!AC!       do k=1,ntra
+!AC!        do j=minorig,nl
+!AC!         do il=1,ncum
+!AC!          if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
+!AC!     :       (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
+!AC!            traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
+!AC!     :            +(1.-sij(il,i,j))*tra(il,nk(il),k)
+!AC!          endif
+!AC!         enddo
+!AC!        enddo
+!AC!       enddo
 
 c
@@ -1730,13 +1730,13 @@
  750  continue
 
-      do j=1,ntra
-       do i=minorig+1,nl
-        do il=1,ncum
-         if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
-          traent(il,i,i,j)=tra(il,nk(il),j)
-         endif
-        enddo
-       enddo
-      enddo
+!AC!      do j=1,ntra
+!AC!       do i=minorig+1,nl
+!AC!        do il=1,ncum
+!AC!         if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
+!AC!          traent(il,i,i,j)=tra(il,nk(il),j)
+!AC!         endif
+!AC!        enddo
+!AC!       enddo
+!AC!      enddo
 
       do 100 j=minorig,nl
@@ -1904,12 +1904,12 @@
       enddo ! il
 
-      do j=1,ntra
-       do il=1,ncum
-        if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
-     :     .and. csum(il,i).lt.m(il,i) ) then
-         traent(il,i,i,j)=tra(il,nk(il),j)
-        endif
-       enddo
-      enddo
+!AC!      do j=1,ntra
+!AC!       do il=1,ncum
+!AC!        if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
+!AC!     :     .and. csum(il,i).lt.m(il,i) ) then
+!AC!         traent(il,i,i,j)=tra(il,nk(il),j)
+!AC!        endif
+!AC!       enddo
+!AC!      enddo
 789   continue
 c      
@@ -2014,11 +2014,11 @@
          enddo
         enddo
-        do k=1,ntra
-         do i=1,nd
-          do il=1,ncum
-           trap(il,i,k)=tra(il,i,k)
-          enddo
-         enddo
-        enddo
+!AC!        do k=1,ntra
+!AC!         do i=1,nd
+!AC!          do il=1,ncum
+!AC!           trap(il,i,k)=tra(il,i,k)
+!AC!          enddo
+!AC!         enddo
+!AC!        enddo
 c
 c   ***  check whether ep(inb)=0, if so, skip precipitating    ***
@@ -2341,21 +2341,21 @@
 c    ***       find tracer concentrations in precipitating downdraft     ***
 c
-      do j=1,ntra
-       do il = 1,ncum
-       if (i.lt.inb(il) .and. lwork(il)) then
-c
-         if(mplus(il))then
-          trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
-     :              +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
-          trap(il,i,j)=trap(il,i,j)/mp(il,i)
-         else ! if (mplus(il))
-          if(mp(il,i+1).gt.1.0e-16)then
-           trap(il,i,j)=trap(il,i+1,j)
-          endif
-         endif ! (mplus(il)) else if (.not.mplus(il))
-c
-        endif ! (i.lt.inb(il) .and. lwork(il))
-       enddo
-      end do
+!AC!      do j=1,ntra
+!AC!       do il = 1,ncum
+!AC!       if (i.lt.inb(il) .and. lwork(il)) then
+!AC!c
+!AC!         if(mplus(il))then
+!AC!          trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
+!AC!     :              +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
+!AC!          trap(il,i,j)=trap(il,i,j)/mp(il,i)
+!AC!         else ! if (mplus(il))
+!AC!          if(mp(il,i+1).gt.1.0e-16)then
+!AC!           trap(il,i,j)=trap(il,i+1,j)
+!AC!          endif
+!AC!         endif ! (mplus(il)) else if (.not.mplus(il))
+!AC!c
+!AC!        endif ! (i.lt.inb(il) .and. lwork(il))
+!AC!       enddo
+!AC!      end do
 
 400   continue
@@ -2484,11 +2484,11 @@
       enddo
 c       print*,'cv3_yield initialisation 2'
-      do j=1,ntra
-       do i=1,nd
-        do il=1,ncum
-          ftra(il,i,j)=0.0
-        enddo
-       enddo
-      enddo
+!AC!      do j=1,ntra
+!AC!       do i=1,nd
+!AC!        do il=1,ncum
+!AC!          ftra(il,i,j)=0.0
+!AC!        enddo
+!AC!       enddo
+!AC!      enddo
 c       print*,'cv3_yield initialisation 3'
       do i=1,nl
@@ -2649,19 +2649,19 @@
 
 
-      do j=1,ntra
-       do il=1,ncum
-        if (iflag(il) .le. 1) then
-        if (cvflag_grav) then
-         ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il)
-     :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
-     :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
-        else
-         ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
-     :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
-     :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
-        endif
-        endif  ! iflag
-       enddo
-      enddo
+!AC!     do j=1,ntra
+!AC!      do il=1,ncum
+!AC!       if (iflag(il) .le. 1) then
+!AC!       if (cvflag_grav) then
+!AC!        ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il)
+!AC!    :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
+!AC!    :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
+!AC!       else
+!AC!        ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
+!AC!    :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
+!AC!    :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
+!AC!       endif
+!AC!       endif  ! iflag
+!AC!      enddo
+!AC!     enddo
 
        do j=2,nl
@@ -2687,21 +2687,21 @@
       enddo
 
-      do k=1,ntra
-       do j=2,nl
-        do il=1,ncum
-         if (j.le.inb(il) .and. iflag(il) .le. 1) then
-
-          if (cvflag_grav) then
-           ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1)
-     :                *(traent(il,j,1,k)-tra(il,1,k))
-          else
-           ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1)
-     :                *(traent(il,j,1,k)-tra(il,1,k))
-          endif
-
-         endif
-        enddo
-       enddo
-      enddo
+!AC!      do k=1,ntra
+!AC!       do j=2,nl
+!AC!        do il=1,ncum
+!AC!         if (j.le.inb(il) .and. iflag(il) .le. 1) then
+!AC!
+!AC!          if (cvflag_grav) then
+!AC!           ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1)
+!AC!     :                *(traent(il,j,1,k)-tra(il,1,k))
+!AC!          else
+!AC!           ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1)
+!AC!     :                *(traent(il,j,1,k)-tra(il,1,k))
+!AC!          endif
+!AC!
+!AC!         endif
+!AC!        enddo
+!AC!       enddo
+!AC!      enddo
 c      print*,'cv3_yield apres ft'
 c
@@ -2865,21 +2865,21 @@
 1350  continue
 
-      do k=1,ntra
-       do il=1,ncum
-        if (i.le.inb(il) .and. iflag(il) .le. 1) then
-         dpinv=1.0/(ph(il,i)-ph(il,i+1))
-         cpinv=1.0/cpn(il,i)
-         if (cvflag_grav) then
-           ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv
-     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
-     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
-         else
-           ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
-     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
-     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
-         endif
-        endif
-       enddo
-      enddo
+!AC!      do k=1,ntra
+!AC!       do il=1,ncum
+!AC!        if (i.le.inb(il) .and. iflag(il) .le. 1) then
+!AC!         dpinv=1.0/(ph(il,i)-ph(il,i+1))
+!AC!         cpinv=1.0/cpn(il,i)
+!AC!         if (cvflag_grav) then
+!AC!           ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv
+!AC!     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
+!AC!     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
+!AC!         else
+!AC!           ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
+!AC!     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
+!AC!     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
+!AC!         endif
+!AC!        endif
+!AC!       enddo
+!AC!      enddo
 
       do 480 k=1,i-1
@@ -2938,21 +2938,21 @@
 480   continue
 
-      do j=1,ntra
-       do k=1,i-1
-        do il=1,ncum
-         if (i.le.inb(il) .and. iflag(il) .le. 1) then
-          dpinv=1.0/(ph(il,i)-ph(il,i+1))
-          cpinv=1.0/cpn(il,i)
-          if (cvflag_grav) then
-           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
-     :        *(traent(il,k,i,j)-tra(il,i,j))
-          else
-           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
-     :        *(traent(il,k,i,j)-tra(il,i,j))
-          endif
-         endif
-        enddo
-       enddo
-      enddo
+!AC!      do j=1,ntra
+!AC!       do k=1,i-1
+!AC!        do il=1,ncum
+!AC!         if (i.le.inb(il) .and. iflag(il) .le. 1) then
+!AC!          dpinv=1.0/(ph(il,i)-ph(il,i+1))
+!AC!          cpinv=1.0/cpn(il,i)
+!AC!          if (cvflag_grav) then
+!AC!           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
+!AC!     :        *(traent(il,k,i,j)-tra(il,i,j))
+!AC!          else
+!AC!           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
+!AC!     :        *(traent(il,k,i,j)-tra(il,i,j))
+!AC!          endif
+!AC!         endif
+!AC!        enddo
+!AC!       enddo
+!AC!      enddo
 
       do 490 k=i,nl+1
@@ -3004,22 +3004,22 @@
 490   continue
 
-      do j=1,ntra
-       do k=i,nl+1
-        do il=1,ncum
-         if (i.le.inb(il) .and. k.le.inb(il)
-     $                .and. iflag(il) .le. 1) then
-          dpinv=1.0/(ph(il,i)-ph(il,i+1))
-          cpinv=1.0/cpn(il,i)
-          if (cvflag_grav) then
-           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
-     :         *(traent(il,k,i,j)-tra(il,i,j))
-          else
-           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
-     :             *(traent(il,k,i,j)-tra(il,i,j))
-          endif
-         endif ! i and k
-        enddo
-       enddo
-      enddo
+!AC!      do j=1,ntra
+!AC!       do k=i,nl+1
+!AC!        do il=1,ncum
+!AC!         if (i.le.inb(il) .and. k.le.inb(il)
+!AC!     $                .and. iflag(il) .le. 1) then
+!AC!          dpinv=1.0/(ph(il,i)-ph(il,i+1))
+!AC!          cpinv=1.0/cpn(il,i)
+!AC!          if (cvflag_grav) then
+!AC!           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
+!AC!     :         *(traent(il,k,i,j)-tra(il,i,j))
+!AC!          else
+!AC!           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
+!AC!     :             *(traent(il,k,i,j)-tra(il,i,j))
+!AC!          endif
+!AC!         endif ! i and k
+!AC!        enddo
+!AC!       enddo
+!AC!      enddo
 
 c sb: interface with the cloud parameterization:          ! cld
@@ -3052,22 +3052,22 @@
       enddo
 
-      do j=1,ntra
-       do il=1,ncum
-        if (i.le.inb(il) .and. iflag(il) .le. 1) then
-         dpinv=1.0/(ph(il,i)-ph(il,i+1))
-         cpinv=1.0/cpn(il,i)
-
-         if (cvflag_grav) then
-          ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv
-     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
-     :     -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
-         else
-          ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
-     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
-     :     -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
-         endif
-        endif ! i
-       enddo
-      enddo
+!AC!      do j=1,ntra
+!AC!       do il=1,ncum
+!AC!        if (i.le.inb(il) .and. iflag(il) .le. 1) then
+!AC!         dpinv=1.0/(ph(il,i)-ph(il,i+1))
+!AC!         cpinv=1.0/cpn(il,i)
+!AC!
+!AC!         if (cvflag_grav) then
+!AC!          ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv
+!AC!     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
+!AC!     :     -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
+!AC!         else
+!AC!          ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
+!AC!     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
+!AC!     :     -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
+!AC!         endif
+!AC!        endif ! i
+!AC!       enddo
+!AC!      enddo
 
 
@@ -3146,27 +3146,27 @@
 503   continue
 
-      do j=1,ntra
-       do il=1,ncum
-        IF (iflag(il) .le. 1) THEN
-	IF (cvflag_grav) then
-        ex=0.01*grav*ment(il,inb(il),inb(il))
-     :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
-     :      /(ph(i l,inb(il))-ph(il,inb(il)+1))
-        ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
-        ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
-     :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
-     :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
-	else
-        ex=0.1*ment(il,inb(il),inb(il))
-     :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
-     :      /(ph(i l,inb(il))-ph(il,inb(il)+1))
-        ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
-        ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
-     :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
-     :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
-        ENDIF   !cvflag grav
-        ENDIF    !iflag
-       enddo
-      enddo
+!AC!      do j=1,ntra
+!AC!       do il=1,ncum
+!AC!        IF (iflag(il) .le. 1) THEN
+!AC!	IF (cvflag_grav) then
+!AC!        ex=0.01*grav*ment(il,inb(il),inb(il))
+!AC!     :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
+!AC!     :      /(ph(i l,inb(il))-ph(il,inb(il)+1))
+!AC!        ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
+!AC!        ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
+!AC!     :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
+!AC!     :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
+!AC!	else
+!AC!        ex=0.1*ment(il,inb(il),inb(il))
+!AC!     :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
+!AC!     :      /(ph(i l,inb(il))-ph(il,inb(il)+1))
+!AC!        ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
+!AC!        ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
+!AC!     :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
+!AC!     :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
+!AC!        ENDIF   !cvflag grav
+!AC!        ENDIF    !iflag
+!AC!       enddo
+!AC!      enddo
 
 c
@@ -3287,13 +3287,14 @@
       ENDDO
       ENDDO
-      DO j = 1,ntra
-      DO i = 1,nl
-       DO il = 1,ncum
-        IF (iflag(il) .le. 1) THEN
-         ftra(il,i,j) = ftra(il,i,j)/alpha_qpos(il)
-        ENDIF
-       ENDDO
-      ENDDO
-      ENDDO
+
+!AC!      DO j = 1,ntra
+!AC!      DO i = 1,nl
+!AC!       DO il = 1,ncum
+!AC!        IF (iflag(il) .le. 1) THEN
+!AC!         ftra(il,i,j) = ftra(il,i,j)/alpha_qpos(il)
+!AC!        ENDIF
+!AC!       ENDDO
+!AC!      ENDDO
+!AC!      ENDDO
 
 c
@@ -3539,4 +3540,29 @@
         end
 
+!AC!
+      SUBROUTINE cv3_tracer(nloc,len,ncum,nd,na,
+     &                        ment,sij,da,phi)
+        implicit none
+c inputs:
+        integer ncum, nd, na, nloc,len
+        real ment(nloc,na,na),sij(nloc,na,na)
+c ouputs:
+        real da(nloc,na),phi(nloc,na,na)
+c local variables:
+        integer i,j,k
+c        
+        da(:,:)=0.
+c
+        do j=1,na
+          do k=1,na
+            do i=1,ncum
+            da(i,j)=da(i,j)+(1.-sij(i,k,j))*ment(i,k,j)
+            phi(i,j,k)=sij(i,k,j)*ment(i,k,j)
+            end do 
+          end do 
+        end do 
+        return
+        end
+!AC!
 
       SUBROUTINE cv3_uncompress(nloc,len,ncum,nd,ntra,idcum
@@ -3609,12 +3635,12 @@
 
 
-        do 2100 j=1,ntra
-c oct3         do 2110 k=1,nl
-         do 2110 k=1,nd ! oct3
-          do 2120 i=1,ncum
-            ftra1(idcum(i),k,j)=ftra(i,k,j)
- 2120     continue
- 2110    continue
- 2100   continue
+!AC!        do 2100 j=1,ntra
+!AC!c oct3         do 2110 k=1,nl
+!AC!         do 2110 k=1,nd ! oct3
+!AC!          do 2120 i=1,ncum
+!AC!            ftra1(idcum(i),k,j)=ftra(i,k,j)
+!AC! 2120     continue
+!AC! 2110    continue
+!AC! 2100   continue
         return
         end
Index: LMDZ5/branches/testing/libf/phylmd/cv3a_compress.F
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/cv3a_compress.F	(revision 1665)
+++ LMDZ5/branches/testing/libf/phylmd/cv3a_compress.F	(revision 1669)
@@ -116,16 +116,16 @@
  110  continue
 
-      do 121 j=1,ntra
-ccccc      do 111 k=1,nl+1
-      do 111 k=1,nd
-       nn=0
-      do 101 i=1,len
-      if(iflag1(i).eq.0)then
-       nn=nn+1
-       tra(nn,k,j)=tra1(i,k,j)
-      endif
- 101  continue
- 111  continue
- 121  continue
+!AC!      do 121 j=1,ntra
+!AC!ccccc      do 111 k=1,nl+1
+!AC!      do 111 k=1,nd
+!AC!       nn=0
+!AC!      do 101 i=1,len
+!AC!      if(iflag1(i).eq.0)then
+!AC!       nn=nn+1
+!AC!       tra(nn,k,j)=tra1(i,k,j)
+!AC!      endif
+!AC! 101  continue
+!AC! 111  continue
+!AC! 121  continue
 
       if (nn.ne.ncum) then
Index: LMDZ5/branches/testing/libf/phylmd/cv3a_uncompress.F
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/cv3a_uncompress.F	(revision 1665)
+++ LMDZ5/branches/testing/libf/phylmd/cv3a_uncompress.F	(revision 1669)
@@ -9,4 +9,7 @@
      :         ,Plim1,Plim2,asupmax,supmax0
      :         ,asupmaxmin
+!AC!
+     :         ,da,phi
+!AC!
      o         ,iflag1,kbas1,ktop1
      :         ,precip1,cbmf1,plcl1,plfc1,wbeff1,sig1,w01,ptop21
@@ -17,5 +20,8 @@
      :         ,ftd1,fqd1
      :         ,Plim11,Plim21,asupmax1,supmax01
-     :         ,asupmaxmin1     )
+     :         ,asupmaxmin1     
+!AC!
+     :         ,da1,phi1  )
+!AC!
 ***************************************************************
 *                                                             *
@@ -50,5 +56,7 @@
       real asupmax(nloc,nd),supmax0(nloc)
       real asupmaxmin(nloc)
-
+!AC!
+      real da(nloc,nd),phi(nloc,nd,nd)
+!AC!
 c outputs:
       integer iflag1(len),kbas1(len),ktop1(len)
@@ -68,4 +76,7 @@
       real asupmax1(len,nd),supmax01(len)
       real asupmaxmin1(len)
+!AC!
+      real da1(nloc,nd),phi1(nloc,nd,nd)
+!AC!
 c
 c local variables:
@@ -111,5 +122,8 @@
             fqd1(idcum(i),k)=fqd(i,k)
             asupmax1(idcum(i),k)=asupmax(i,k)
- 2010     continue
+!AC!
+            da1(idcum(i),k)=da(i,k)
+!AC!
+ 2010    continue
  2020   continue
 
@@ -119,12 +133,23 @@
 
 
-        do 2100 j=1,ntra
-c oct3         do 2110 k=1,nl
-         do 2110 k=1,nd ! oct3
-          do 2120 i=1,ncum
-            ftra1(idcum(i),k,j)=ftra(i,k,j)
- 2120     continue
- 2110    continue
- 2100   continue
+!AC!        do 2100 j=1,ntra
+!AC!c oct3         do 2110 k=1,nl
+!AC!         do 2110 k=1,nd ! oct3
+!AC!          do 2120 i=1,ncum
+!AC!            ftra1(idcum(i),k,j)=ftra(i,k,j)
+!AC! 2120     continue
+!AC! 2110    continue
+!AC! 2100   continue
+
+!AC!
+       do j=1,nd
+         do k=1,nd 
+          do i=1,ncum
+            phi1(idcum(i),k,j)=phi(i,k,j)
+          end do
+         end do
+        end do
+!AC!
+
 c
 c        do 2220 k2=1,nd
Index: LMDZ5/branches/testing/libf/phylmd/cv3p_mixing.F
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/cv3p_mixing.F	(revision 1665)
+++ LMDZ5/branches/testing/libf/phylmd/cv3p_mixing.F	(revision 1669)
@@ -118,9 +118,14 @@
             elij(i,k,j)=0.0
             hent(i,k,j)=0.0
-            ment(i,k,j)=0.0
-            sij(i,k,j)=0.0
+!AC!            ment(i,k,j)=0.0
+!AC!            sij(i,k,j)=0.0
  385      continue
  390    continue
  400  continue
+
+!AC!
+      ment(1:ncum,1:nd,1:nd)=0.0
+      sij(1:ncum,1:nd,1:nd)=0.0
+!AC!
 
       do k=1,ntra
Index: LMDZ5/branches/testing/libf/phylmd/cva_driver.F
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/cva_driver.F	(revision 1665)
+++ LMDZ5/branches/testing/libf/phylmd/cva_driver.F	(revision 1669)
@@ -20,5 +20,8 @@
      &                   ftd1,fqd1,
      &                   Plim11,Plim21,asupmax1,supmax01,asupmaxmin1
-     &                   ,lalim_conv)
+     &                   ,lalim_conv, 
+!AC!
+     &                   da1,phi1)
+!AC!
 ***************************************************************
 *                                                             *
@@ -171,4 +174,8 @@
       real tvp1(len,nd)
 c
+!AC!
+      real da1(len,nd),phi1(len,nd,nd)
+      real da(len,nd),phi(len,nd,nd)
+!AC!
       real ftd1(len,nd)
       real fqd1(len,nd)
@@ -912,4 +919,16 @@
       endif
 
+!AC!
+!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+! --- passive tracers
+!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+      if (iflag_con.eq.3) then
+       CALL cv3_tracer(nloc,len,ncum,nd,nd,
+     :                  ment,sij,da,phi)
+      endif
+
+!AC!
+
 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 ! --- UNCOMPRESS THE FIELDS
@@ -928,4 +947,7 @@
      :          ,Plim1,Plim2,asupmax,supmax0
      :          ,asupmaxmin
+!AC!
+     :          ,da,phi
+!AC!
      o          ,iflag1,kbas1,ktop1
      o          ,precip1,cbmf1,plcl1,plfc1,wbeff1,sig1,w01,ptop21
@@ -936,5 +958,8 @@
      o          ,ftd1,fqd1
      o          ,Plim11,Plim21,asupmax1,supmax01
-     o          ,asupmaxmin1)
+     o          ,asupmaxmin1
+!AC!
+     o          ,da1,phi1)
+!AC!
       endif
 
Index: LMDZ5/branches/testing/libf/phylmd/inistats.F
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/inistats.F	(revision 1665)
+++ LMDZ5/branches/testing/libf/phylmd/inistats.F	(revision 1669)
@@ -18,5 +18,4 @@
       real, dimension(istime) :: lt
       integer :: nvarid
-      real, dimension(llm) :: pseudoalt
 
       write (*,*) 
Index: LMDZ5/branches/testing/libf/phylmd/phys_output_mod.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/phys_output_mod.F90	(revision 1665)
+++ LMDZ5/branches/testing/libf/phylmd/phys_output_mod.F90	(revision 1669)
@@ -33,4 +33,7 @@
   CHARACTER(len=20), dimension(nfiles), private, save   :: type_ecri
   !$OMP THREADPRIVATE(nhorim, nvertm, zoutm,zdtime,type_ecri)
+ ! swaero_diag : flag indicates if it is necessary to do calculation for some aerosol diagnostics
+  logical, save                                :: swaero_diag=.FALSE.
+
 
   !   integer, save                     :: nid_hf3d 
@@ -260,4 +263,24 @@
   type(ctrl_out),save :: o_wape         = ctrl_out((/ 1, 1, 1, 10, 10, 10 /),'wape')
 
+!!! nrlmd le 10/04/2012
+
+!-------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')
+
+!-------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
@@ -365,6 +388,8 @@
 
   type(ctrl_out),save :: o_topswad      = ctrl_out((/ 2, 10, 10, 10, 10, 10 /),'topswad')
+  type(ctrl_out),save :: o_topswad0     = ctrl_out((/ 2, 10, 10, 10, 10, 10 /),'topswad0')
   type(ctrl_out),save :: o_topswai      = ctrl_out((/ 2, 10, 10, 10, 10, 10 /),'topswai')
   type(ctrl_out),save :: o_solswad      = ctrl_out((/ 2, 10, 10, 10, 10, 10 /),'solswad')
+  type(ctrl_out),save :: o_solswad0     = ctrl_out((/ 2, 10, 10, 10, 10, 10 /),'solswad0')
   type(ctrl_out),save :: o_solswai      = ctrl_out((/ 2, 10, 10, 10, 10, 10 /),'solswai')
 
@@ -432,4 +457,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 +520,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 +558,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 +563,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 +647,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 +709,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'
@@ -792,18 +819,22 @@
     DO iff=1,nfiles
 
+       ! Calculate ecrit_files for all files
+       if ( chtimestep(iff).eq.'DefFreq' ) then
+          ! Par defaut ecrit_files = (ecrit_mensuel ecrit_jour ecrit_hf ...)*86400.
+          ecrit_files(iff)=ecrit_files(iff)*86400.
+       else
+          call convers_timesteps(chtimestep(iff),dtime,ecrit_files(iff)) 
+       endif
+       write(lunout,*)'ecrit_files(',iff,')= ',ecrit_files(iff)
+
+       zoutm(iff) = ecrit_files(iff) ! Frequence ou l on ecrit en seconde
+
        IF (clef_files(iff)) THEN
 
-          if ( chtimestep(iff).eq.'DefFreq' ) then
-             ! Par defaut ecrit_files = (ecrit_mensuel ecrit_jour ecrit_hf ...)*86400.
-             ecrit_files(iff)=ecrit_files(iff)*86400.
-          else
-             call convers_timesteps(chtimestep(iff),dtime,ecrit_files(iff)) 
-          endif
-          write(lunout,*)'ecrit_files(',iff,')= ',ecrit_files(iff)
-
-          zoutm(iff) = ecrit_files(iff) ! Frequence ou l on ecrit en seconde
-
           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 !!
@@ -1099,5 +1130,9 @@
                   o_topswad%flag,o_topswad%name, "ADE at TOA", "W/m2")
              CALL histdef2d(iff,clef_stations(iff), &
+                  o_topswad0%flag,o_topswad0%name, "ADE clear-sky at TOA", "W/m2")
+             CALL histdef2d(iff,clef_stations(iff), &
                   o_solswad%flag,o_solswad%name, "ADE at SRF", "W/m2")
+             CALL histdef2d(iff,clef_stations(iff), &
+                  o_solswad0%flag,o_solswad0%name, "ADE clear-sky at SRF", "W/m2")
 
              CALL histdef2d(iff,clef_stations(iff), &
@@ -1216,6 +1251,4 @@
                      o_wbeff%flag,o_wbeff%name, "Conv. updraft velocity at LFC (<100)", "m/s")
              end if
-             CALL histdef2d(iff,clef_stations(iff), &
-                  o_prw%flag,o_prw%name, "Precipitable water", "kg/m2")
              IF (.NOT.clef_stations(iff)) THEN 
                 !
@@ -1253,4 +1286,6 @@
           ENDIF !iflag_con .GE. 3
 
+          CALL histdef2d(iff,clef_stations(iff), &
+               o_prw%flag,o_prw%name, "Precipitable water", "kg/m2")
           CALL histdef2d(iff,clef_stations(iff), &
                o_s_pblh%flag,o_s_pblh%name, "Boundary Layer Height", "m")
@@ -1318,8 +1353,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 +1412,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 +1518,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 +1529,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 +1542,31 @@
                 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
+
+        IF (iflag_trig_bl>=1) THEN
+ 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")
+       ENDIF  !(iflag_trig_bl>=1)
+
+        IF (iflag_clos_bl>=1) THEN
+ 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")
+       ENDIF  !(iflag_clos_bl>=1)
+
+!!! fin nrlmd le 10/04/2012
 
           CALL histdef3d(iff,clef_stations(iff),o_dtlsc%flag,o_dtlsc%name, "Condensation dT", "K/s")
@@ -1519,5 +1581,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 +1593,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 +1608,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")
@@ -1708,4 +1768,15 @@
 
     ENDDO !  iff
+
+    ! Updated write frequencies due to phys_out_filetimesteps. 
+    ! Write frequencies are now in seconds.  
+    ecrit_mth = ecrit_files(1)
+    ecrit_day = ecrit_files(2)
+    ecrit_hf  = ecrit_files(3)
+    ecrit_ins = ecrit_files(4)
+    ecrit_LES = ecrit_files(5)
+    ecrit_ins = ecrit_files(6)
+
+    write(lunout,*)'swaero_diag=',swaero_diag
     write(lunout,*)'Fin phys_output_mod.F90'
   end subroutine phys_output_open
@@ -1756,4 +1827,11 @@
        endif
     endif
+
+    ! Set swaero_diag=true if at least one of the concerned variables are defined 
+    if (nomvar=='topswad' .OR. nomvar=='topswai' .OR. nomvar=='solswad' .OR. nomvar=='solswai' ) THEN
+       if  ( flag_var(iff)<=lev_files(iff) ) then
+          swaero_diag=.TRUE.
+       end if
+    end if
   end subroutine histdef2d
 
Index: LMDZ5/branches/testing/libf/phylmd/phys_output_write.h
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/phys_output_write.h	(revision 1665)
+++ LMDZ5/branches/testing/libf/phylmd/phys_output_write.h	(revision 1669)
@@ -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),
@@ -631,9 +631,4 @@
         end if
 
-        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 (.NOT.clef_stations(iff)) THEN
       IF (o_cape_max%flag(iff)<=lev_files(iff)) THEN
@@ -671,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
@@ -681,4 +676,9 @@
       
       ENDIF !iflag_con .GE. 3
+
+        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
@@ -801,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),
@@ -809,4 +810,5 @@
      $o_alp_bl%name,itau_w,alp_bl)
        ENDIF
+      ENDIF !iflag_coupl>=1
       ENDIF !(iflag_con.GE.3)
 
@@ -823,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),
@@ -869,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),
@@ -897,4 +899,67 @@
       ENDIF !(iflag_con.EQ.3) 
  
+!!! nrlmd le 10/04/2012
+
+        IF (iflag_trig_bl>=1) THEN
+          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
+       ENDIF  !(iflag_trig_bl>=1)
+
+        IF (iflag_clos_bl>=1) THEN
+         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
+       ENDIF  !(iflag_clos_bl>=1)
+
+!!! fin nrlmd le 10/04/2012
+
       IF (type_ocean=='slab ') THEN
       IF ( o_slab_bils%flag(iff)<=lev_files(iff)) 
@@ -1193,4 +1258,10 @@
      $            topswad_aero)
           ENDIF
+          IF (o_topswad0%flag(iff)<=lev_files(iff)) THEN
+             CALL histwrite_phy(nid_files(iff),
+     $clef_stations(iff),
+     $o_topswad0%name,itau_w,
+     $            topswad0_aero)
+          ENDIF
           IF (o_solswad%flag(iff)<=lev_files(iff)) THEN
              CALL histwrite_phy(nid_files(iff),
@@ -1198,4 +1269,10 @@
      $o_solswad%name,itau_w,
      $            solswad_aero)
+          ENDIF
+          IF (o_solswad0%flag(iff)<=lev_files(iff)) THEN
+             CALL histwrite_phy(nid_files(iff),
+     $clef_stations(iff),
+     $o_solswad0%name,itau_w,
+     $            solswad0_aero)
           ENDIF
 
@@ -1410,4 +1487,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
 
@@ -1637,4 +1719,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
@@ -1680,5 +1768,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
@@ -1744,5 +1832,5 @@
        ENDIF
 
-      endif ! iflag_thermals>1
+      endif ! iflag_thermals>=1
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -1791,5 +1879,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
@@ -1820,8 +1908,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
@@ -2051,17 +2135,29 @@
        ENDIF
 
-       IF (o_mcd%flag(iff)<=lev_files(iff)) THEN
-      zx_tmp_fi3d(1:klon,1:klev)=-1 * (dnwd(1:klon,1:klev)+
-     $                                 dnwd0(1:klon,1:klev)) 
-      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
-     $o_mcd%name,itau_w,zx_tmp_fi3d)
-       ENDIF
-
-       IF (o_dmc%flag(iff)<=lev_files(iff)) THEN
-      zx_tmp_fi3d(1:klon,1:klev)=upwd(1:klon,1:klev) +
-     $  dnwd(1:klon,1:klev)+ dnwd0(1:klon,1:klev) 
-      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
-     $o_dmc%name,itau_w,zx_tmp_fi3d)
-       ENDIF
+       if (iflag_con >= 3) then
+          IF (o_mcd%flag(iff)<=lev_files(iff)) THEN
+             zx_tmp_fi3d(1:klon,1:klev)=-1 * (dnwd(1:klon,1:klev)+
+     $            dnwd0(1:klon,1:klev)) 
+             CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     $            o_mcd%name,itau_w,zx_tmp_fi3d)
+          ENDIF
+
+          IF (o_dmc%flag(iff)<=lev_files(iff)) THEN
+             zx_tmp_fi3d(1:klon,1:klev)=upwd(1:klon,1:klev) +
+     $            dnwd(1:klon,1:klev)+ dnwd0(1:klon,1:klev) 
+             CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     $            o_dmc%name,itau_w,zx_tmp_fi3d)
+          ENDIF
+       else if (iflag_con == 2) then
+          IF (o_mcd%flag(iff) <= lev_files(iff)) THEN
+             CALL histwrite_phy(nid_files(iff), clef_stations(iff),
+     $            o_mcd%name, itau_w, pmfd)
+          ENDIF
+
+          IF (o_dmc%flag(iff) <= lev_files(iff)) THEN
+             CALL histwrite_phy(nid_files(iff), clef_stations(iff),
+     $            o_dmc%name, itau_w, pmfu + pmfd)
+          ENDIF
+       end if
 
        IF (o_ref_liq%flag(iff)<=lev_files(iff)) THEN
Index: LMDZ5/branches/testing/libf/phylmd/phys_state_var_mod.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/phys_state_var_mod.F90	(revision 1665)
+++ LMDZ5/branches/testing/libf/phylmd/phys_state_var_mod.F90	(revision 1669)
@@ -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/branches/testing/libf/phylmd/physiq.F
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/physiq.F	(revision 1665)
+++ LMDZ5/branches/testing/libf/phylmd/physiq.F	(revision 1669)
@@ -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):
@@ -1212,7 +1247,8 @@
       LOGICAL, SAVE :: mskocean_beta
 c$OMP THREADPRIVATE(mskocean_beta)
-      REAL, dimension(klon, klev) :: beta       ! facteur sur cldtaurad et cldemirad pour evaluer les retros liees aux CRF
-      REAL, dimension(klon, klev) :: cldtaurad  ! epaisseur optique pour radlwsw,COSP
-      REAL, dimension(klon, klev) :: cldemirad  ! emissivite pour radlwsw,COSP 
+      REAL, dimension(klon, klev) :: beta         ! facteur sur cldtaurad et cldemirad pour evaluer les retros liees aux CRF
+      REAL, dimension(klon, klev) :: cldtaurad    ! epaisseur optique pour radlwsw,COSP
+      REAL, dimension(klon, klev) :: cldtaupirad  ! epaisseur optique pour radlwsw,COSP cas pre-industrial
+      REAL, dimension(klon, klev) :: cldemirad    ! emissivite pour radlwsw,COSP 
       INTEGER :: nbtr_tmp ! Number of tracer inside concvl
       REAL, dimension(klon,klev) :: sh_in ! Specific humidity entering in phytrac
@@ -1354,5 +1390,4 @@
          solswad(:)=0.
 
-         lambda_th(:,:)=0.
          wmax_th(:)=0.
          tau_overturning_th(:)=0.
@@ -1490,8 +1525,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 +1553,5 @@
        print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
       ENDIF
+
 c
       ALLOCATE(tabCFMIP(nCFMIP))
@@ -1624,19 +1662,7 @@
 
 #endif
-
-
-         ecrit_hf = ecrit_hf * un_jour
-cIM
-         IF(ecrit_day.LE.1.) THEN
-          ecrit_day = ecrit_day * un_jour !en secondes
-         ENDIF
-cIM
-         ecrit_mth = ecrit_mth * un_jour
-         ecrit_ins = ecrit_ins * un_jour
          ecrit_reg = ecrit_reg * un_jour
          ecrit_tra = ecrit_tra * un_jour
-         ecrit_LES = ecrit_LES * un_jour
-c
-
+       
 cXXXPB Positionner date0 pour initialisation de ORCHIDEE
       date0 = jD_ref 
@@ -1735,4 +1761,5 @@
 !
       itap   = itap + 1
+c
 !
 ! Update fraction of the sub-surfaces (pctsrf) and 
@@ -2042,58 +2069,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 +2270,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 +2295,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 +2312,7 @@
           endif
        enddo
-
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
+       endif
        do i=1,klon
           if (alp(i)>alp_max) then
@@ -2586,5 +2629,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 +2650,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 +2761,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 +2948,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 +2976,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,
@@ -2962,5 +3100,5 @@
 cjq - introduce the aerosol direct and first indirect radiative forcings
 cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
-      IF (ok_ade.OR.ok_aie) THEN
+      IF (flag_aerosol .gt. 0) THEN
          IF (.NOT. aerosol_couple)
      &        CALL readaerosol_optic(
@@ -3247,6 +3385,7 @@
 cIM betaCRF
 c
-      cldtaurad = cldtau
-      cldemirad = cldemi
+      cldtaurad   = cldtau
+      cldtaupirad = cldtaupi
+      cldemirad   = cldemi
 c
       if(lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND.
@@ -3265,6 +3404,7 @@
          beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
         endif
-        cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
-        cldemirad(i,k) = cldemi(i,k) * beta(i,k)
+        cldtaurad(i,k)   = cldtau(i,k) * beta(i,k)
+        cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
+        cldemirad(i,k)   = cldemi(i,k) * beta(i,k)
        ENDDO
        ENDDO 
@@ -3287,6 +3427,7 @@
           beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
          endif
-        cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
-        cldemirad(i,k) = cldemi(i,k) * beta(i,k)
+        cldtaurad(i,k)   = cldtau(i,k) * beta(i,k)
+        cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
+        cldemirad(i,k)   = cldemi(i,k) * beta(i,k)
         endif
 c
@@ -3337,5 +3478,5 @@
      s        topsw_aero, topsw0_aero,
      s        solsw_aero, solsw0_aero,
-     e        cldtaupi,
+     e        cldtaupirad,
      s        topswai_aero, solswai_aero)
             
@@ -3351,4 +3492,8 @@
        RCFC12 = RCFC12_act
 c
+      IF (prt_level .GE.10) THEN
+       print *,' ->radlwsw, number 1 '
+      ENDIF
+c
          CALL radlwsw
      e        (dist, rmu0, fract, 
@@ -3356,7 +3501,7 @@
      e        t_seri,q_seri,wo,
      e        cldfra, cldemirad, cldtaurad,
-     e        ok_ade, ok_aie,
+     e        ok_ade, ok_aie, flag_aerosol,
      e        tau_aero, piz_aero, cg_aero,
-     e        cldtaupi,new_aod,
+     e        cldtaupirad,new_aod,
      e        zqsat, flwc, fiwc,
      s        heat,heat0,cool,cool0,radsol,albpla,
@@ -3388,4 +3533,8 @@
        RCFC12 = RCFC12_per
 c
+      IF (prt_level .GE.10) THEN
+       print *,' ->radlwsw, number 2 '
+      ENDIF
+c
          CALL radlwsw
      e        (dist, rmu0, fract, 
@@ -3393,5 +3542,5 @@
      e        t_seri,q_seri,wo,
      e        cldfra, cldemi, cldtau,
-     e        ok_ade, ok_aie,
+     e        ok_ade, ok_aie, flag_aerosol,
      e        tau_aero, piz_aero, cg_aero,
      e        cldtaupi,new_aod,
@@ -3479,4 +3628,8 @@
 c Appeler le programme de parametrisation de l'orographie
 c a l'echelle sous-maille:
+c
+      IF (prt_level .GE.10) THEN
+       print *,' call orography ? ', ok_orodr
+      ENDIF
 c
       IF (ok_orodr) THEN
@@ -3569,4 +3722,5 @@
 
        IF (ok_hines) then
+
          CALL hines_gwd(klon,klev,dtime,paprs,pplay,
      i                  rlat,t_seri,u_seri,v_seri,
@@ -3576,6 +3730,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 +3760,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 +4055,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 +4177,5 @@
 c=============================================================
 
-      if (iflag_thermals>1) then
+      if (iflag_thermals>=1) then
       d_t_lscth=0.
       d_t_lscst=0.
Index: LMDZ5/branches/testing/libf/phylmd/radlwsw.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/radlwsw.F90	(revision 1665)
+++ LMDZ5/branches/testing/libf/phylmd/radlwsw.F90	(revision 1669)
@@ -10,5 +10,5 @@
    t,q,wo,&
    cldfra, cldemi, cldtaupd,&
-   ok_ade, ok_aie,&
+   ok_ade, ok_aie, flag_aerosol,&
    tau_aero, piz_aero, cg_aero,&
    cldtaupi, new_aod, &
@@ -56,4 +56,5 @@
   ! ok_ade---input-L- apply the Aerosol Direct Effect or not?
   ! ok_aie---input-L- apply the Aerosol Indirect Effect or not?
+  ! flag_aerosol-input-I- aerosol flag from 0 to 6
   ! tau_ae, piz_ae, cg_ae-input-R- aerosol optical properties (calculated in aeropt.F)
   ! cldtaupi-input-R- epaisseur optique des nuages dans le visible
@@ -119,4 +120,5 @@
 
   LOGICAL, INTENT(in)  :: ok_ade, ok_aie                                 ! switches whether to use aerosol direct (indirect) effects or not
+  INTEGER, INTENT(in)  :: flag_aerosol                                   ! takes value 0 (no aerosol) or 1 to 6 (aerosols)
   REAL,    INTENT(in)  :: cldfra(KLON,KLEV), cldemi(KLON,KLEV), cldtaupd(KLON,KLEV)
   REAL,    INTENT(in)  :: tau_aero(KLON,KLEV,9,2)                        ! aerosol optical properties (see aeropt.F)
@@ -354,9 +356,9 @@
                zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0,&
                ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,&
-               tau_aero(:,:,5,:), piz_aero(:,:,5,:), cg_aero(:,:,5,:),& 
+               tauaero(:,:,5,:), pizaero(:,:,5,:), cgaero(:,:,5,:),& 
                PTAUA, POMEGAA,&
                ztopswadaero,zsolswadaero,&
                ztopswaiaero,zsolswaiaero,& 
-               ok_ade, ok_aie) 
+               ok_ade, ok_aie, flag_aerosol) 
           
        ELSE ! new_aod=T         
@@ -377,6 +379,5 @@
                zsolsw_aero,zsolsw0_aero,&
                ztopswcf_aero,zsolswcf_aero, & 
-               ok_ade, ok_aie) 
-          
+               ok_ade, ok_aie, flag_aerosol) 
        ENDIF
 
Index: LMDZ5/branches/testing/libf/phylmd/regr_lat_time_climoz_m.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/regr_lat_time_climoz_m.F90	(revision 1665)
+++ LMDZ5/branches/testing/libf/phylmd/regr_lat_time_climoz_m.F90	(revision 1669)
@@ -224,5 +224,5 @@
     ! Get the  number of months:
     call nf95_inq_dimid(ncid_in, "time", dimid)
-    call nf95_inquire_dimension(ncid_in, dimid, len=n_month)
+    call nf95_inquire_dimension(ncid_in, dimid, nclen=n_month)
 
     allocate(o3_in(n_lat, n_plev, n_month, read_climoz))
Index: LMDZ5/branches/testing/libf/phylmd/sw_aeroAR4.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/sw_aeroAR4.F90	(revision 1665)
+++ LMDZ5/branches/testing/libf/phylmd/sw_aeroAR4.F90	(revision 1669)
@@ -18,8 +18,8 @@
      PSOLSWAERO,PSOLSW0AERO,&
      PTOPSWCFAERO,PSOLSWCFAERO,&
-     ok_ade, ok_aie )
+     ok_ade, ok_aie, flag_aerosol )
 
   USE dimphy
-
+  USE phys_output_mod, ONLY : swaero_diag
   IMPLICIT NONE
 
@@ -56,7 +56,8 @@
   !     --------------
   !        ORIGINAL : 89-07-14
-  !        95-01-01   J.-J. MORCRETTE  Direct/Diffuse Albedo
-  !        03-11-27   J. QUAAS Introduce aerosol forcings (based on BOUCHER)
-  !        09-04      A. COZIC - C.DEANDREIS Indroduce NAT/BC/POM/DUST/SS aerosol forcing
+  !        1995-01-01  J.-J. MORCRETTE  Direct/Diffuse Albedo
+  !        2003-11-27  J. QUAAS Introduce aerosol forcings (based on BOUCHER)
+  !        2009-04     A. COZIC - C.DEANDREIS Indroduce NAT/BC/POM/DUST/SS aerosol forcing
+  !        2012-09     O. BOUCHER - reorganise aerosol cases with ok_ade, ok_aie, flag_aerosol
   !     ------------------------------------------------------------------
   !
@@ -82,5 +83,5 @@
 
   REAL(KIND=8) PCLDSW(KDLON,KFLEV)    ! CLOUD FRACTION
-  REAL(KIND=8) PTAU(KDLON,2,KFLEV)    ! CLOUD OPTICAL THICKNESS
+  REAL(KIND=8) PTAU(KDLON,2,KFLEV)    ! CLOUD OPTICAL THICKNESS (pre-industrial value)
   REAL(KIND=8) PCG(KDLON,2,KFLEV)     ! ASYMETRY FACTOR
   REAL(KIND=8) POMEGA(KDLON,2,KFLEV)  ! SINGLE SCATTERING ALBEDO
@@ -132,24 +133,25 @@
   !$OMP THREADPRIVATE(initialized)
 
-  !jq-Introduced for aerosol forcings
+  !jq-local flag introduced for aerosol forcings
   REAL(KIND=8), SAVE :: flag_aer
   !$OMP THREADPRIVATE(flag_aer)
 
   LOGICAL ok_ade, ok_aie    ! use aerosol forcings or not?
+  INTEGER flag_aerosol      ! global flag for aerosol 0 (no aerosol) or 1-5 (aerosols)
   REAL(KIND=8) tauaero(kdlon,kflev,9,2)  ! aerosol optical properties
   REAL(KIND=8) pizaero(kdlon,kflev,9,2)  ! (see aeropt.F)
   REAL(KIND=8) cgaero(kdlon,kflev,9,2)   ! -"-
-  REAL(KIND=8) PTAUA(KDLON,2,KFLEV)    ! CLOUD OPTICAL THICKNESS (pre-industrial value)
+  REAL(KIND=8) PTAUA(KDLON,2,KFLEV)    ! CLOUD OPTICAL THICKNESS (present-day value)
   REAL(KIND=8) POMEGAA(KDLON,2,KFLEV)  ! SINGLE SCATTERING ALBEDO
   REAL(KIND=8) PTOPSWADAERO(KDLON)     ! SHORTWAVE FLUX AT T.O.A.(+AEROSOL DIR)
   REAL(KIND=8) PSOLSWADAERO(KDLON)     ! SHORTWAVE FLUX AT SURFACE(+AEROSOL DIR)
-  REAL(KIND=8) PTOPSWAD0AERO(KDLON)     ! SHORTWAVE FLUX AT T.O.A.(+AEROSOL DIR)
-  REAL(KIND=8) PSOLSWAD0AERO(KDLON)     ! SHORTWAVE FLUX AT SURFACE(+AEROSOL DIR)
+  REAL(KIND=8) PTOPSWAD0AERO(KDLON)    ! SHORTWAVE FLUX AT T.O.A.(+AEROSOL DIR)
+  REAL(KIND=8) PSOLSWAD0AERO(KDLON)    ! SHORTWAVE FLUX AT SURFACE(+AEROSOL DIR)
   REAL(KIND=8) PTOPSWAIAERO(KDLON)     ! SHORTWAVE FLUX AT T.O.A.(+AEROSOL IND)
   REAL(KIND=8) PSOLSWAIAERO(KDLON)     ! SHORTWAVE FLUX AT SURFACE(+AEROSOL IND)
-  REAL(KIND=8) PTOPSWAERO(KDLON,9)	 ! SW TOA AS DRF nat & ant 
-  REAL(KIND=8) PTOPSW0AERO(KDLON,9)	 ! SW SRF AS DRF nat & ant 
-  REAL(KIND=8) PSOLSWAERO(KDLON,9)	 ! SW TOA CS DRF nat & ant
-  REAL(KIND=8) PSOLSW0AERO(KDLON,9)	 ! SW SRF CS DRF nat & ant
+  REAL(KIND=8) PTOPSWAERO(KDLON,9)     ! SW TOA AS DRF nat & ant 
+  REAL(KIND=8) PTOPSW0AERO(KDLON,9)    ! SW SRF AS DRF nat & ant 
+  REAL(KIND=8) PSOLSWAERO(KDLON,9)     ! SW TOA CS DRF nat & ant
+  REAL(KIND=8) PSOLSW0AERO(KDLON,9)    ! SW SRF CS DRF nat & ant
   REAL(KIND=8) PTOPSWCFAERO(KDLON,3)   !  SW TOA AS cloudRF nat & ant 
   REAL(KIND=8) PSOLSWCFAERO(KDLON,3)   !  SW SRF AS cloudRF nat & ant 
@@ -179,24 +181,13 @@
 
 ! Key to define the aerosol effect acting on climate
-! 0: aerosol feedback active according to ok_ade, ok_aie  DEFAULT 
-! 1: no feedback , zero aerosol fluxes are used for climate, diagnostics according to ok_ade_ok_aie
-! 2: feedback according to total aerosol direct effect used for climate, diagnostics according to ok_ade, ok_aie
-! 3: feedback according to natural aerosol direct effect used for climate, diagnostics according to ok_ade_ok_aie
-
-  INTEGER,SAVE :: AEROSOLFEEDBACK_ACTIVE = 0
+! OB: AEROSOLFEEDBACK_ACTIVE is now a LOGICAL
+! TRUE: fluxes use natural and/or anthropogenic aerosols according to ok_ade and ok_aie, DEFAULT
+! FALSE: fluxes use no aerosols (case 1)
+
+  LOGICAL,SAVE :: AEROSOLFEEDBACK_ACTIVE = .TRUE.
 !$OMP THREADPRIVATE(AEROSOLFEEDBACK_ACTIVE)  
 
       CHARACTER (LEN=20) :: modname='sw_aeroAR4'
       CHARACTER (LEN=80) :: abort_message
-
-  IF ((.not. ok_ade) .and. (AEROSOLFEEDBACK_ACTIVE .ge. 2)) THEN
-     abort_message ='Error: direct effect is not activated but assumed to be active - see sw_aeroAR4.F90'
-     CALL abort_gcm (modname,abort_message,1)
-  ENDIF
-  AEROSOLFEEDBACK_ACTIVE=MIN(MAX(AEROSOLFEEDBACK_ACTIVE,0),3)
-  IF  (AEROSOLFEEDBACK_ACTIVE .gt. 3) THEN
-     abort_message ='Error: AEROSOLFEEDBACK_ACTIVE options go only until 3'
-     CALL abort_gcm (modname,abort_message,1)
-  ENDIF
 
   IF(.NOT.initialized) THEN
@@ -209,8 +200,25 @@
      ALLOCATE(ZFSUPAI_AERO(KDLON,KFLEV+1))
      ALLOCATE(ZFSDNAI_AERO(KDLON,KFLEV+1))
-     ALLOCATE(ZFSUP_AERO (KDLON,KFLEV+1,9))
-     ALLOCATE(ZFSDN_AERO (KDLON,KFLEV+1,9))
-     ALLOCATE(ZFSUP0_AERO(KDLON,KFLEV+1,9))
-     ALLOCATE(ZFSDN0_AERO(KDLON,KFLEV+1,9))
+!-OB decrease size of these arrays to what is needed
+!                | direct effect 
+!ind effect      | no aerosol   natural  total
+!natural (PTAU)  |   1            3       2     --ZFSUP/ZFSDN
+!total (PTAUA)   |                5       4     --ZFSUP/ZFSDN
+!no cloud        |   1            3       2     --ZFSUP0/ZFSDN0
+! so we need which case when ? 
+! ok_ade and ok_aie = 4-5, 4-2 and 2 
+! ok_ade and not ok_aie = 2-3 and 2 
+! not ok_ade and ok_aie = 5-3 and 5 
+! not ok_ade and not ok_aie = 3 
+! therefore the cases have the folliwng switches 
+! 3 = not ok_ade or not ok_aie 
+! 4 = ok_ade and ok_aie 
+! 2 = ok_ade 
+! 5 = ok_aie 
+     ALLOCATE(ZFSUP_AERO (KDLON,KFLEV+1,5))
+     ALLOCATE(ZFSDN_AERO (KDLON,KFLEV+1,5))
+     ALLOCATE(ZFSUP0_AERO(KDLON,KFLEV+1,3))
+     ALLOCATE(ZFSDN0_AERO(KDLON,KFLEV+1,3))
+! end OB modif
      ZFSUPAD_AERO(:,:)=0.
      ZFSDNAD_AERO(:,:)=0.
@@ -226,5 +234,5 @@
 
   IF (appel1er) THEN
-     WRITE(lunout,*) 'SW calling frequency : ', swpas
+     WRITE(lunout,*)'SW calling frequency : ', swpas
      WRITE(lunout,*) "   In general, it should be 1"
      appel1er = .FALSE.
@@ -241,6 +249,6 @@
      ENDDO
 
-! clear sky is either computed IF no direct effect is asked for, or for extended diag) 
-     IF (( lev_histmth .ge. 4 ) .or. ( .not. ok_ade )) THEN    
+! clear sky with no aerosols at all is computed IF ACTIVEFEEDBACK_ACTIVE is false or for extended diag
+     IF ( swaero_diag .or. .not. AEROSOLFEEDBACK_ACTIVE .OR. flag_aerosol .EQ. 0 ) THEN    
 
      ! clear-sky: zero aerosol effect
@@ -268,8 +276,8 @@
         ENDDO
      ENDDO
-     ENDIF
-
-! cloudy sky is either computed IF no indirect effect is asked for, or for extended diag) 
-     IF (( lev_histmth .ge. 4 ) .or. ( .not. ok_aie )) THEN    
+     ENDIF ! swaero_diag .or. .not. AEROSOLFEEDBACK_ACTIVE
+
+! cloudy sky with no aerosols at all is either computed IF no indirect effect is asked for, or for extended diag
+     IF ( swaero_diag .or. .not. AEROSOLFEEDBACK_ACTIVE .OR. flag_aerosol .EQ. 0 ) THEN    
      ! cloudy-sky: zero aerosol effect
      flag_aer=0.0
@@ -297,10 +305,41 @@
         ENDDO
      ENDDO
-     ENDIF
-
+     ENDIF ! swaero_diag .or. .not. AEROSOLFEEDBACK_ACTIVE
+
+     IF (flag_aerosol .GT. 0 ) THEN
+
+     IF (ok_ade.and.swaero_diag .or. .not. ok_ade) THEN
+
+        ! clear sky direct effect natural aerosol
+        ! CAS AER (3)
+        flag_aer=1.0
+        CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,&
+             PRMU0,PFRAC,PTAVE,PWV,&
+             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
+        INU = 1
+        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),&
+             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             ZFD, ZFU)
+        INU = 2
+        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),&
+             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             PWV, PQS,&
+             ZFDOWN, ZFUP)
+
+        DO JK = 1 , KFLEV+1
+           DO JL = 1, KDLON
+              ZFSUP0_AERO(JL,JK,3) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
+              ZFSDN0_AERO(JL,JK,3) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
+           ENDDO
+        ENDDO
+     ENDIF !--end not swaero_diag or not ok_ade
 
      IF (ok_ade) THEN
 
-        ! clear sky (Anne Cozic 03/07/2007) direct effect of total aerosol
+        ! clear sky direct effect of total aerosol
         ! CAS AER (2)
         flag_aer=1.0
@@ -329,7 +368,8 @@
         ENDDO
 
-! cloudy sky is either computed IF no indirect effect is asked for, or for extended diag) 
-        IF (( lev_histmth .ge. 2 ) .or. (.not. ok_aie)) THEN  
-        ! cloudy-sky aerosol direct effect of total aerosol
+        ! cloudy-sky with natural aerosols for indirect effect 
+        ! but total aerosols for direct effect
+        ! PTAU
+        ! CAS AER (2)
         flag_aer=1.0
         CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
@@ -356,37 +396,13 @@
            ENDDO
         ENDDO
-        ENDIF
-
-! natural aeroosl clear sky is  computed  for extended diag) 
-        IF ( lev_histmth .ge. 4 ) THEN            
-        ! clear sky direct effect natural aerosol
-        flag_aer=1.0
-        CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,&
-             PRMU0,PFRAC,PTAVE,PWV,&
-             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
-        INU = 1
-        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
-             tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),&
-             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
-             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
-             ZFD, ZFU)
-        INU = 2
-        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
-             tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),&
-             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
-             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
-             PWV, PQS,&
-             ZFDOWN, ZFUP)
-
-        DO JK = 1 , KFLEV+1
-           DO JL = 1, KDLON
-              ZFSUP0_AERO(JL,JK,3) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
-              ZFSDN0_AERO(JL,JK,3) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
-           ENDDO
-        ENDDO
-        ENDIF
-
-! cloud sky natural is for extended diagnostics
-        IF ( lev_histmth .ge. 2 ) THEN
+
+     ENDIF !-end ok_ade
+
+     IF ( .not. ok_ade .or. .not. ok_aie ) THEN
+
+        ! cloudy-sky with natural aerosols for indirect effect 
+        ! and natural aerosols for direct effect
+        ! PTAU
+        ! CAS AER (3)
         ! cloudy-sky direct effect natural aerosol
         flag_aer=1.0
@@ -414,11 +430,13 @@
            ENDDO
         ENDDO
-        ENDIF
-
-     ENDIF ! ok_ade
-
-! cloudy sky needs to be computed in all cases IF ok_aie is activated
-     IF (ok_aie) THEN
-        !jq   cloudy-sky + aerosol direct + aerosol indirect of total aerosol
+
+     ENDIF  !--true/false or false/true
+
+     IF (ok_ade .and. ok_aie) THEN
+
+        ! cloudy-sky with total aerosols for indirect effect 
+        ! and total aerosols for direct effect
+        ! PTAUA
+        ! CAS AER (2)
         flag_aer=1.0
         CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
@@ -438,4 +456,5 @@
              PWV, PQS,&
              ZFDOWN, ZFUP)
+
         DO JK = 1 , KFLEV+1
            DO JL = 1, KDLON
@@ -444,11 +463,46 @@
            ENDDO
         ENDDO
+  
+      ENDIF ! ok_ade .and. ok_aie
+
+     IF (ok_aie) THEN
+        ! cloudy-sky with total aerosols for indirect effect 
+        ! and natural aerosols for direct effect
+        ! PTAUA
+        ! CAS AER (3)
+        flag_aer=1.0
+        CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
+             PRMU0,PFRAC,PTAVE,PWV,&
+             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
+        INU = 1
+        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),&
+             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGAA, ZOZ, ZRMU, ZSEC, PTAUA, ZUD,&
+             ZFD, ZFU)
+        INU = 2
+        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),&
+             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGAA, ZOZ, ZRMU, ZSEC, PTAUA, ZUD,&
+             PWV, PQS,&
+             ZFDOWN, ZFUP)
+  
+        DO JK = 1 , KFLEV+1
+           DO JL = 1, KDLON
+              ZFSUP_AERO(JL,JK,5) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
+              ZFSDN_AERO(JL,JK,5) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL) 
+           ENDDO
+        ENDDO
+
      ENDIF ! ok_aie      
 
+     ENDIF !--if flag_aerosol GT 0
+
      itapsw = 0
   ENDIF
   itapsw = itapsw + 1
 
-  IF  ( AEROSOLFEEDBACK_ACTIVE .eq. 0) THEN
+  IF  ( AEROSOLFEEDBACK_ACTIVE .AND. flag_aerosol .GT. 0 ) THEN
   IF ( ok_ade .and. ok_aie  ) THEN
     ZFSUP(:,:) =    ZFSUP_AERO(:,:,4)
@@ -457,4 +511,5 @@
     ZFSDN0(:,:) =   ZFSDN0_AERO(:,:,2)
   ENDIF
+
   IF ( ok_ade .and. (.not. ok_aie) )  THEN
     ZFSUP(:,:) =    ZFSUP_AERO(:,:,2)
@@ -465,11 +520,21 @@
 
   IF ( (.not. ok_ade) .and. ok_aie  )  THEN
-    print*,'Warning: indirect effect in cloudy regions includes direct aerosol effect'
-    ZFSUP(:,:) =    ZFSUP_AERO(:,:,4)
-    ZFSDN(:,:) =    ZFSDN_AERO(:,:,4)
-    ZFSUP0(:,:) =   ZFSUP0_AERO(:,:,1)
-    ZFSDN0(:,:) =   ZFSDN0_AERO(:,:,1)
-  ENDIF
+    ZFSUP(:,:) =    ZFSUP_AERO(:,:,5)
+    ZFSDN(:,:) =    ZFSDN_AERO(:,:,5)
+    ZFSUP0(:,:) =   ZFSUP0_AERO(:,:,3)
+    ZFSDN0(:,:) =   ZFSDN0_AERO(:,:,3)
+  ENDIF
+
   IF ((.not. ok_ade) .and. (.not. ok_aie)) THEN
+    ZFSUP(:,:) =    ZFSUP_AERO(:,:,3)
+    ZFSDN(:,:) =    ZFSDN_AERO(:,:,3)
+    ZFSUP0(:,:) =   ZFSUP0_AERO(:,:,3)
+    ZFSDN0(:,:) =   ZFSDN0_AERO(:,:,3)
+  ENDIF
+
+! MS the following allows to compute the forcing diagostics without
+! letting the aerosol forcing act on the meteorology
+! SEE logic above
+  ELSE 
     ZFSUP(:,:) =    ZFSUP_AERO(:,:,1)
     ZFSDN(:,:) =    ZFSDN_AERO(:,:,1)
@@ -478,15 +543,5 @@
   ENDIF
 
-! MS the following allows to compute the forcing diagostics without
-! letting the aerosol forcing act on the meteorology
-! SEE logic above
-  ELSEIF  ( AEROSOLFEEDBACK_ACTIVE .gt. 0) THEN
-    ZFSUP(:,:) =    ZFSUP_AERO(:,:,AEROSOLFEEDBACK_ACTIVE)
-    ZFSDN(:,:) =    ZFSDN_AERO(:,:,AEROSOLFEEDBACK_ACTIVE)
-    ZFSUP0(:,:) =   ZFSUP0_AERO(:,:,AEROSOLFEEDBACK_ACTIVE)
-    ZFSDN0(:,:) =   ZFSDN0_AERO(:,:,AEROSOLFEEDBACK_ACTIVE)
-  ENDIF
-  
-
+! Now computes heating rates
   DO k = 1, KFLEV
      kpl1 = k+1
@@ -511,24 +566,33 @@
      PTOPSW(i) = ZFSDN(i,KFLEV+1) - ZFSUP(i,KFLEV+1)
 
-
 ! net anthropogenic forcing direct and 1st indirect effect diagnostics
 ! requires a natural aerosol field read and used 
 ! Difference of net fluxes from double call to radiation
 
-
 IF (ok_ade) THEN
 
 ! indices 1: natural; 2 anthropogenic 
+
 ! TOA/SRF all sky natural forcing
      PSOLSWAERO(i,1) = (ZFSDN_AERO(i,1,3) - ZFSUP_AERO(i,1,3))-(ZFSDN_AERO(i,1,1) - ZFSUP_AERO(i,1,1))
      PTOPSWAERO(i,1) = (ZFSDN_AERO(i,KFLEV+1,3) - ZFSUP_AERO(i,KFLEV+1,3))- (ZFSDN_AERO(i,KFLEV+1,1) - ZFSUP_AERO(i,KFLEV+1,1))
 
+! TOA/SRF clear sky natural forcing
+     PSOLSW0AERO(i,1) = (ZFSDN0_AERO(i,1,3) - ZFSUP0_AERO(i,1,3))-(ZFSDN0_AERO(i,1,1) - ZFSUP0_AERO(i,1,1))
+     PTOPSW0AERO(i,1) = (ZFSDN0_AERO(i,KFLEV+1,3) - ZFSUP0_AERO(i,KFLEV+1,3))-(ZFSDN0_AERO(i,KFLEV+1,1) - ZFSUP0_AERO(i,KFLEV+1,1))
+
+   IF (ok_aie) THEN 
+
+! TOA/SRF all sky anthropogenic forcing
+     PSOLSWAERO(i,2) = (ZFSDN_AERO(i,1,4) - ZFSUP_AERO(i,1,4))-(ZFSDN_AERO(i,1,5) - ZFSUP_AERO(i,1,5))
+     PTOPSWAERO(i,2) = (ZFSDN_AERO(i,KFLEV+1,4) - ZFSUP_AERO(i,KFLEV+1,4))- (ZFSDN_AERO(i,KFLEV+1,5) - ZFSUP_AERO(i,KFLEV+1,5))
+
+   ELSE 
+
 ! TOA/SRF all sky anthropogenic forcing
      PSOLSWAERO(i,2) = (ZFSDN_AERO(i,1,2) - ZFSUP_AERO(i,1,2))-(ZFSDN_AERO(i,1,3) - ZFSUP_AERO(i,1,3))
      PTOPSWAERO(i,2) = (ZFSDN_AERO(i,KFLEV+1,2) - ZFSUP_AERO(i,KFLEV+1,2))- (ZFSDN_AERO(i,KFLEV+1,3) - ZFSUP_AERO(i,KFLEV+1,3))
 
-! TOA/SRF clear sky natural forcing
-     PSOLSW0AERO(i,1) = (ZFSDN0_AERO(i,1,3) - ZFSUP0_AERO(i,1,3))-(ZFSDN0_AERO(i,1,1) - ZFSUP0_AERO(i,1,1))
-     PTOPSW0AERO(i,1) = (ZFSDN0_AERO(i,KFLEV+1,3) - ZFSUP0_AERO(i,KFLEV+1,3))-(ZFSDN0_AERO(i,KFLEV+1,1) - ZFSUP0_AERO(i,KFLEV+1,1))
+   ENDIF
 
 ! TOA/SRF clear sky anthropogenic forcing
@@ -536,4 +600,11 @@
      PTOPSW0AERO(i,2) = (ZFSDN0_AERO(i,KFLEV+1,2) - ZFSUP0_AERO(i,KFLEV+1,2))-(ZFSDN0_AERO(i,KFLEV+1,3) - ZFSUP0_AERO(i,KFLEV+1,3))
 
+! direct anthropogenic forcing , as in old LMDzT, however differences of net fluxes
+     PSOLSWADAERO(i) = PSOLSWAERO(i,2)
+     PTOPSWADAERO(i) = PTOPSWAERO(i,2)
+     PSOLSWAD0AERO(i) = PSOLSW0AERO(i,2)
+     PTOPSWAD0AERO(i) = PTOPSW0AERO(i,2)
+
+! OB: these diagnostics may not always work but who need them
 ! Cloud forcing indices 1: natural; 2 anthropogenic; 3: zero aerosol direct effect
 ! Instantaneously computed cloudy sky direct aerosol effect, cloud forcing due to aerosols above clouds
@@ -552,18 +623,17 @@
      PTOPSWCFAERO(i,3) = (ZFSDN_AERO(i,KFLEV+1,1) - ZFSUP_AERO(i,KFLEV+1,1))- (ZFSDN0_AERO(i,KFLEV+1,1) - ZFSUP0_AERO(i,KFLEV+1,1))
 
-! direct anthropogenic forcing , as in old LMDzT, however differences of net fluxes
-     PSOLSWADAERO(i) = PSOLSWAERO(i,2)
-     PTOPSWADAERO(i) = PTOPSWAERO(i,2)
-     PSOLSWAD0AERO(i) = PSOLSW0AERO(i,2)
-     PTOPSWAD0AERO(i) = PTOPSW0AERO(i,2)
-
 ENDIF
 
-
 IF (ok_aie) THEN
+   IF (ok_ade) THEN 
      PSOLSWAIAERO(i) = (ZFSDN_AERO(i,1,4) - ZFSUP_AERO(i,1,4))-(ZFSDN_AERO(i,1,2) - ZFSUP_AERO(i,1,2))
      PTOPSWAIAERO(i) = (ZFSDN_AERO(i,KFLEV+1,4) - ZFSUP_AERO(i,KFLEV+1,4))-(ZFSDN_AERO(i,KFLEV+1,2) - ZFSUP_AERO(i,KFLEV+1,2))
+   ELSE 
+     PSOLSWAIAERO(i) = (ZFSDN_AERO(i,1,5) - ZFSUP_AERO(i,1,5))-(ZFSDN_AERO(i,1,3) - ZFSUP_AERO(i,1,3))
+     PTOPSWAIAERO(i) = (ZFSDN_AERO(i,KFLEV+1,5) - ZFSUP_AERO(i,KFLEV+1,5))-(ZFSDN_AERO(i,KFLEV+1,3) - ZFSUP_AERO(i,KFLEV+1,3))
+   ENDIF
 ENDIF
 
-  ENDDO
+ENDDO
+
 END SUBROUTINE SW_AEROAR4
Index: LMDZ5/branches/testing/libf/phylmd/thermcell.h
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/thermcell.h	(revision 1665)
+++ LMDZ5/branches/testing/libf/phylmd/thermcell.h	(revision 1669)
@@ -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,10 @@
       common/ctherm5/iflag_thermals_ed,iflag_thermals_optflux
 
+!!! 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/)
+!$OMP THREADPRIVATE(/ctherm6/,/ctherm7/,/ctherm8/)
Index: LMDZ5/branches/testing/libf/phylmd/thermcell_main.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/thermcell_main.F90	(revision 1665)
+++ LMDZ5/branches/testing/libf/phylmd/thermcell_main.F90	(revision 1669)
@@ -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
+
Index: LMDZ5/branches/testing/libf/phylmd/wake.F
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/wake.F	(revision 1665)
+++ LMDZ5/branches/testing/libf/phylmd/wake.F	(revision 1669)
@@ -561,7 +561,8 @@
       ENDDO
 
-c      On evite kupper = 1
+c      On evite kupper = 1 et kupper = klev
       DO i=1,klon
         kupper(i) = max(kupper(i),2)
+        kupper(i) = min(kupper(i),klev-1)
       ENDDO
 
