Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 543)
+++ /trunk/LMDZ.MARS/README	(revision 544)
@@ -1414,2 +1414,14 @@
 >> Code about 10% faster [tests carried out with mesoscale summer Tharsis run]
 >> Check with a GCM run that results similar
+
+== 27/02/12 == AC
+>> Continuation of thermals setting, comparisons with mesoscale results on Case C
+>> Added possibility to call gcm (or 1d) with constant prescribed sensible heat flux, in the spirit of direct comparisons with LES
+   ... This is directly comparable to the variable tke_heat_flux in namelist.input
+   ... Requires the use of yamada4.F from terrestrial GCM (mixes more, seem more numerically stable)
+   ... Usually requires high timesteps (>1000) to avoids crashes. Best approach is to compare
+   height of first model level z1 and teta_1 between LES and 1D, and increase the timesteps until results 
+   between the two models are comparable (might require a slitghly different tke_heat_flux between the two models
+   due to difference in vertical diffusion schemes and subgrid effects)
+   ... Syntax for use is to add "tke_heat_flux = ..." in callphys.def
+>> Corrected some stuff with tke transport in thermals (which is anyway deactivated but one day maybe...)
Index: /trunk/LMDZ.MARS/libf/phymars/callkeys.h
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/callkeys.h	(revision 543)
+++ /trunk/LMDZ.MARS/libf/phymars/callkeys.h	(revision 544)
@@ -18,5 +18,6 @@
      &   ,dustbin,nqchem_min,nltemodel,nircorr
      
-      COMMON/callkeys_r/topdustref,solarcondate,semi,alphan
+      COMMON/callkeys_r/topdustref,solarcondate,semi,alphan,            &
+     &   tke_heat_flux
      
       LOGICAL callrad,calldifv,calladj,callcond,callsoil,               &
@@ -36,4 +37,5 @@
       real alphan
       real solarcondate
+      real tke_heat_flux
 
       integer iddist
Index: /trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90	(revision 543)
+++ /trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90	(revision 544)
@@ -158,12 +158,11 @@
        endif
 
-!       dtke_thermals=.false. !! default setting
-!       !call getin("dtke_thermals",dtke_thermals)
-!
-!       IF(dtke_thermals) THEN
-!          DO l=1,nlayermx
-!              q2_therm(:,l)=0.5*(q2(:,l)+q2(:,l+1))
-!          ENDDO
-!       ENDIF
+       dtke_thermals=.false. !! default setting
+       call getin("dtke_thermals",dtke_thermals)
+       IF(dtke_thermals) THEN
+          DO l=1,nlayermx
+              q2_therm(:,l)=0.5*(q2(:,l)+q2(:,l+1))
+          ENDDO
+       ENDIF
 
 ! **********************************************************************
@@ -338,4 +337,11 @@
       endif
 
+      if (dtke_thermals) then
+      modname='tke'
+      call thermcell_dqup(ngridmx,nlayermx,ptimestep     &
+     &     ,fm_therm,entr_therm,detr_therm,  &
+     &    masse,q2_therm,dq2_therm,modname,zdz)
+      endif
+
       DO ig=1,ngridmx
          hfmax(ig)=MAXVAL(heatFlux(ig,:)+heatFlux_down(ig,:))
@@ -369,12 +375,12 @@
            endif
 
-!           IF(dtke_thermals) THEN
-!              DO l=2,nlayermx
-!                 pbl_dtke(:,l)=0.5*(dq2_therm(:,l-1)+dq2_therm(:,l))
-!              ENDDO
-!  
-!              pbl_dtke(:,1)=0.5*dq2_therm(:,1)
-!              pbl_dtke(:,nlayermx+1)=0.
-!           ENDIF
+           IF(dtke_thermals) THEN
+              DO l=2,nlayermx
+                 pbl_dtke(:,l)=0.5*(dq2_therm(:,l-1)+dq2_therm(:,l))
+              ENDDO
+  
+              pbl_dtke(:,1)=0.5*dq2_therm(:,1)
+              pbl_dtke(:,nlayermx+1)=0.
+           ENDIF
 
 
Index: /trunk/LMDZ.MARS/libf/phymars/inifis.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/inifis.F	(revision 543)
+++ /trunk/LMDZ.MARS/libf/phymars/inifis.F	(revision 544)
@@ -203,4 +203,10 @@
          write(*,*) " topdustref = ",topdustref
 
+         write(*,*) "Prescribed surface thermal flux (H/(rho*cp),K m/s)"
+         tke_heat_flux=0. ! default value
+         call getin("tke_heat_flux",tke_heat_flux)
+         write(*,*) " tke_heat_flux = ",tke_heat_flux
+         write(*,*) " 0 means the usual schemes are computing"
+
          write(*,*) "call radiative transfer ?"
          callrad=.true. ! default value
Index: /trunk/LMDZ.MARS/libf/phymars/physiq.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 543)
+++ /trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 544)
@@ -11,4 +11,5 @@
 #endif
      $            )
+
 
       IMPLICIT NONE
@@ -338,5 +339,5 @@
 
 c Variables for PBL
-
+      REAL zz1(ngridmx)
       REAL lmax_th_out(ngridmx),zmax_th(ngridmx)
       REAL, SAVE :: wstar(ngridmx)
@@ -769,4 +770,9 @@
 #endif
 
+         IF (tke_heat_flux .ne. 0.) THEN
+             zz1(:)=(pt(:,1)+pdt(:,1)*ptimestep)*(r/g)*
+     &                 (-alog(pplay(:,1)/pplev(:,1)))
+             pdt(:,1)=pdt(:,1) + (tke_heat_flux/zz1(:))*zpopsk(:,1)
+         ENDIF
 
 c        Calling vdif (Martian version WITH CO2 condensation)
@@ -1997,4 +2003,11 @@
          call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
          call WRITEDIAGFI(ngrid,"rho","rho","kg.m-3",1,rho)
+!         call WRITEDIAGFI(ngrid,"dtrad","rad. heat. rate",              &
+!     &              "K.s-1",1,dtrad/zpopsk)
+!        call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate',
+!     &                   'w.m-2',1,zdtsw/zpopsk)
+!        call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
+!     &                   'w.m-2',1,zdtlw/zpopsk)
+
 ! or output in diagfi.nc (for testphys1d)
          call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps)
Index: /trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90	(revision 543)
+++ /trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90	(revision 544)
@@ -113,5 +113,5 @@
       REAL zdz,zbuoy(ngridmx,nlayermx),zw2m
       LOGICAL activecell(ngridmx),activetmp(ngridmx)
-      REAL a1,b1,ae,be,ad,bd,fdfu,b1inv,a1inv,omega
+      REAL a1,b1,ae,be,ad,bd,fdfu,b1inv,a1inv,omega,adalim
       INTEGER tic
 
@@ -381,4 +381,5 @@
       b1inv=b1
       omega=0.
+      adalim=0.
 
 ! One good config for 34/35 levels
@@ -389,21 +390,37 @@
 ! Best configuration for 222 levels:
 
-      omega=0.06
-      b1=0.
-      a1=1.
-      a1inv=0.25*a1
-      b1inv=0.0002
+!      omega=0.06
+!      b1=0.
+!      a1=1.
+!      a1inv=0.25*a1
+!      b1inv=0.0002
+!!
+!!
+!!      ae=0.9*ae
 
 ! Best config for norad 222 levels:
-
-!       omega=0.06
-!       a1=1.
-!       b1=0.
-!       a1inv=a1
-!       be=1.1*be
-!       ad = 0.0004
-!       b1inv=0.00035
+! with yamada4 and alim at sqrt(zlev)
+
+       omega=0.06
+       a1=1.
+       b1=0.
+       a1inv=a1
+       be=1.1*be
+       ad = 0.0004
+       b1inv=0.00035
+       adalim=0.
+
+      b1inv=0.00025
 
 ! Trying stuff :
+
+!      omega=0.04
+!!      b1=0.
+!      a1=1.
+!      a1inv=a1
+!      b1inv=0.0005689
+!!      be=1.1*be
+!!      ae=0.96*ae
+
 
 !       omega=0.06
@@ -431,7 +448,7 @@
             if (ztv(ig,1)>=(ztv(ig,2))) then
                alim_star(ig,1)=MAX((ztv(ig,1)-ztv(ig,2)),0.)  &
-     &                       *sqrt(zlev(ig,2))
+!     &                       *sqrt(zlev(ig,2))
 !     &                       /sqrt(zlev(ig,2))
-!      &                       *zlev(ig,2)
+      &                       *zlev(ig,2)
                lalim(ig)=2
                alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,1)
@@ -444,6 +461,7 @@
            if (ztv(ig,l)>(ztv(ig,l+1)+0.) .and. ztv(ig,1)>=ztv(ig,l) .and. (alim_star(ig,l-1) .ne. 0.)) then ! .and. (zlev(ig,l+1) .lt. 1000.)) then
                alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
-     &                       *sqrt(zlev(ig,l+1))
-!      &                       *zlev(ig,2)
+!     &                       *sqrt(zlev(ig,l+1))
+!     &                       /sqrt(zlev(ig,l+1))
+      &                       *zlev(ig,l+1)
                 lalim(ig)=l+1
                alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
@@ -562,5 +580,8 @@
           if(zbuoy(ig,l) .gt. 0.) then
              if(l .lt. lalim(ig)) then
-                detr_star(ig,l)=0.
+
+!                detr_star(ig,l)=0.
+                 detr_star(ig,l) = f_star(ig,l)*zdz*              &
+            &  adalim
              else
 
@@ -628,5 +649,5 @@
 !---------------------------------------------------------------------------
 
-      DO tic=0,0  ! internal convergence loop
+      DO tic=0,5  ! internal convergence loop
       activetmp(:)=activecell(:) .and. f_star(:,l+1)>1.e-10
       do ig=1,ngridmx
@@ -674,5 +695,9 @@
           if(zbuoy(ig,l) .gt. 0.) then
              if(l .lt. lalim(ig)) then
-                detr_star(ig,l)=0.
+
+!                detr_star(ig,l)=0.
+                 detr_star(ig,l) = f_star(ig,l)*zdz*              &
+            &  adalim
+
              else
                  detr_star(ig,l) = f_star(ig,l)*zdz*              &
Index: /trunk/LMDZ.MARS/libf/phymars/vdifc.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 543)
+++ /trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 544)
@@ -10,5 +10,4 @@
 #endif
      &                )
-       USE ioipsl_getincom
       IMPLICIT NONE
 
@@ -369,4 +368,5 @@
 c    ** schema de diffusion turbulente dans la couche limite
 c       ---------------------------------------------------- 
+       IF (tke_heat_flux .eq. 0.) THEN
 
        CALL vdif_kc(ptimestep,g,pzlev,pzlay
@@ -374,8 +374,12 @@
      &              ,pq2,zkv,zkh,zq)
 
-!      pt(:,:)=ph(:,:)*ppopsk(:,:)
-!      CALL yamada4(ngrid,nlay,ptimestep,g,r,pplev,pt
-!     s   ,pzlev,pzlay,pu,pv,ph,zcdv_true,pq2,zkv,zkh,zkq,ust
-!     s   ,8)
+      ELSE
+
+      pt(:,:)=zh(:,:)*ppopsk(:,:)
+      CALL yamada4(ngrid,nlay,ptimestep,g,r,pplev,pt
+     s   ,pzlev,pzlay,pu,pv,ph,zcdv_true,pq2,zkv,zkh,zkq,ust
+     s   ,8)
+
+      ENDIF
 
       if ((doubleq).and.(ngrid.eq.1)) then
@@ -523,7 +527,11 @@
 
       z4st=4.*5.67e-8*ptimestep
+      IF (tke_heat_flux .eq. 0.) THEN
       DO ig=1,ngrid
          zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig)
       ENDDO
+      ELSE
+         zdplanck(:)=0.
+      ENDIF
 
 ! calcul de zc et zd pour la couche top en prenant en compte le terme
Index: /trunk/LMDZ.MARS/libf/phymars/yamada4.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/yamada4.F	(revision 544)
+++ /trunk/LMDZ.MARS/libf/phymars/yamada4.F	(revision 544)
@@ -0,0 +1,648 @@
+      SUBROUTINE yamada4(ngrid,nlay,dt,g,rconst,plev,temp
+     s   ,zlev,zlay,u,v,teta,cd,q2,km,kn,kq,ustar
+     s   ,iflag_pbl)
+      IMPLICIT NONE
+c.......................................................................
+#include "dimensions.h"
+#include "dimphys.h"
+#include "tracer.h"
+#include "callkeys.h"
+c.......................................................................
+c
+c
+c dt : pas de temps
+c g  : g
+c zlev : altitude a chaque niveau (interface inferieure de la couche
+c        de meme indice)
+c zlay : altitude au centre de chaque couche
+c u,v : vitesse au centre de chaque couche
+c       (en entree : la valeur au debut du pas de temps)
+c teta : temperature potentielle au centre de chaque couche
+c        (en entree : la valeur au debut du pas de temps)
+c cd : cdrag
+c      (en entree : la valeur au debut du pas de temps)
+c q2 : $q^2$ au bas de chaque couche
+c      (en entree : la valeur au debut du pas de temps)
+c      (en sortie : la valeur a la fin du pas de temps)
+c km : diffusivite turbulente de quantite de mouvement (au bas de chaque
+c      couche)
+c      (en sortie : la valeur a la fin du pas de temps)
+c kn : diffusivite turbulente des scalaires (au bas de chaque couche)
+c      (en sortie : la valeur a la fin du pas de temps)
+c
+c  iflag_pbl doit valoir entre 6 et 9
+c      l=6, on prend  systematiquement une longueur d equilibre
+c    iflag_pbl=6 : MY 2.0
+c    iflag_pbl=7 : MY 2.0.Fournier
+c    iflag_pbl=8 : MY 2.5
+c    iflag_pbl>=9 : MY 2.5 avec diffusion verticale
+
+c.......................................................................
+      REAL dt,g,rconst
+      real plev(ngrid,nlay+1),temp(ngrid,nlay)
+      real ustar(ngrid)
+      real kmin,qmin,pblhmin(ngrid),coriol(ngrid)
+      REAL zlev(ngrid,nlay+1)
+      REAL zlay(ngrid,nlay)
+      REAL u(ngrid,nlay)
+      REAL v(ngrid,nlay)
+      REAL teta(ngrid,nlay)
+      REAL cd(ngrid)
+      REAL q2(ngrid,nlay+1),qpre
+      REAL unsdz(ngrid,nlay)
+      REAL unsdzdec(ngrid,nlay+1)
+
+      REAL km(ngrid,nlay+1)
+      REAL kmpre(ngrid,nlay+1),tmp2
+      REAL mpre(ngrid,nlay+1)
+      REAL kn(ngrid,nlay+1)
+      REAL kq(ngrid,nlay+1)
+      real ff(ngrid,nlay+1),delta(ngrid,nlay+1)
+      real aa(ngrid,nlay+1),aa0,aa1
+      integer iflag_pbl,ngrid,nlay,nlev
+
+
+      logical first
+      integer ipas
+      save first,ipas
+cFH/IM     data first,ipas/.true.,0/
+      data first,ipas/.false.,0/
+
+      integer ig,k
+
+
+      real ri,zrif,zalpha,zsm,zsn
+      real rif(ngrid,nlay+1),sm(ngrid,nlay+1),alpha(ngrid,nlay)
+
+      real m2(ngrid,nlay+1),dz(ngrid,nlay+1),zq,n2(ngrid,nlay+1)
+      real dtetadz(ngrid,nlay+1)
+      real m2cstat,mcstat,kmcstat
+      real l(ngrid,nlay+1)
+      real,allocatable,save :: l0(:)
+      real sq(ngrid),sqz(ngrid),zz(ngrid,nlay+1)
+      integer iter
+
+      real ric,rifc,b1,kap
+      save ric,rifc,b1,kap
+      data ric,rifc,b1,kap/0.195,0.191,16.6,0.4/
+      real frif,falpha,fsm
+      real fl,zzz,zl0,zq2,zn2
+
+! MARS
+      REAL,SAVE :: q2min,q2max,knmin,kmmin
+      DATA q2min,q2max,knmin,kmmin/1.E-10,1.E+2,1.E-5,1.E-5/
+
+! mars q2min 0.001 mais trop fort !
+!      PARAMETER (gnmin=-10.E+0)
+!      PARAMETER (gnmax=0.0233E+0)
+!      PARAMETER (a1=0.92E+0)
+!      PARAMETER (a2=0.74E+0)
+!      PARAMETER (b2=10.1E+0)
+!      PARAMETER (c1=0.08E+0)
+
+      real rino(ngrid,nlay+1),smyam(ngrid,nlay),styam(ngrid,nlay)
+     s  ,lyam(ngrid,nlay),knyam(ngrid,nlay)
+     s  ,w2yam(ngrid,nlay),t2yam(ngrid,nlay)
+      logical,save :: firstcall=.true.
+      frif(ri)=0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))
+      falpha(ri)=1.318*(0.2231-ri)/(0.2341-ri)
+      fsm(ri)=1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))
+      fl(zzz,zl0,zq2,zn2)=
+     s     max(min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
+     s     ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10))) ,1.)
+
+
+      nlev=nlay+1
+      
+      if (firstcall) then
+	allocate(l0(ngrid))
+	firstcall=.false.
+      endif
+
+
+      if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.10)) then
+           stop'probleme de coherence dans appel a MY'
+      endif
+
+      ipas=ipas+1
+      if (0.eq.1.and.first) then
+      do ig=1,1000
+         ri=(ig-800.)/500.
+         if (ri.lt.ric) then
+            zrif=frif(ri)
+         else
+            zrif=rifc
+         endif
+         if(zrif.lt.0.16) then
+            zalpha=falpha(zrif)
+            zsm=fsm(zrif)
+         else
+            zalpha=1.12
+            zsm=0.085
+         endif
+c     print*,ri,rif,zalpha,zsm
+      enddo
+      endif
+
+c.......................................................................
+c  les increments verticaux
+c.......................................................................
+c
+c!!!!! allerte !!!!!c
+c!!!!! zlev n'est pas declare a nlay !!!!!c
+c!!!!! ---->
+                                                      DO ig=1,ngrid
+            zlev(ig,nlev)=zlay(ig,nlay)
+     &             +( zlay(ig,nlay) - zlev(ig,nlev-1) )
+                                                      ENDDO
+c!!!!! <----
+c!!!!! allerte !!!!!c
+c
+      DO k=1,nlay
+                                                      DO ig=1,ngrid
+        unsdz(ig,k)=1.E+0/(zlev(ig,k+1)-zlev(ig,k))
+                                                      ENDDO
+      ENDDO
+                                                      DO ig=1,ngrid
+      unsdzdec(ig,1)=1.E+0/(zlay(ig,1)-zlev(ig,1))
+                                                      ENDDO
+      DO k=2,nlay
+                                                      DO ig=1,ngrid
+        unsdzdec(ig,k)=1.E+0/(zlay(ig,k)-zlay(ig,k-1))
+                                                     ENDDO
+      ENDDO
+                                                      DO ig=1,ngrid
+      unsdzdec(ig,nlay+1)=1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay))
+                                                     ENDDO
+c
+c.......................................................................
+
+      do k=2,nlay
+                                                          do ig=1,ngrid
+         dz(ig,k)=zlay(ig,k)-zlay(ig,k-1)
+         m2(ig,k)=((u(ig,k)-u(ig,k-1))**2+(v(ig,k)-v(ig,k-1))**2)
+     s             /(dz(ig,k)*dz(ig,k))
+         dtetadz(ig,k)=(teta(ig,k)-teta(ig,k-1))/dz(ig,k)
+         n2(ig,k)=g*2.*dtetadz(ig,k)/(teta(ig,k-1)+teta(ig,k))
+c        n2(ig,k)=0.
+         ri=n2(ig,k)/max(m2(ig,k),1.e-10)
+         if (ri.lt.ric) then
+            rif(ig,k)=frif(ri)
+         else
+            rif(ig,k)=rifc
+         endif
+         if(rif(ig,k).lt.0.16) then
+            alpha(ig,k)=falpha(rif(ig,k))
+            sm(ig,k)=fsm(rif(ig,k))
+         else
+            alpha(ig,k)=1.12
+            sm(ig,k)=0.085
+         endif
+         zz(ig,k)=b1*m2(ig,k)*(1.-rif(ig,k))*sm(ig,k)
+c     print*,'RIF L=',k,rif(ig,k),ri*alpha(ig,k)
+
+
+                                                          enddo
+      enddo
+
+
+c====================================================================
+c   Au premier appel, on determine l et q2 de facon iterative.
+c iterration pour determiner la longueur de melange
+
+
+      if (first.or.iflag_pbl.eq.6) then
+                                                          do ig=1,ngrid
+!      l0(ig)=10. !EARTH
+      l0(ig)=160.
+                                                          enddo
+      do k=2,nlay-1
+                                                          do ig=1,ngrid
+        l(ig,k)=l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
+                                                          enddo
+      enddo
+
+      do iter=1,10
+                                                          do ig=1,ngrid
+         sq(ig)=1.e-10
+         sqz(ig)=1.e-10
+                                                          enddo
+         do k=2,nlay-1
+                                                          do ig=1,ngrid
+           q2(ig,k)=l(ig,k)**2*zz(ig,k)
+           l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
+           zq=sqrt(q2(ig,k))
+           sqz(ig)=sqz(ig)+zq*zlev(ig,k)*(zlay(ig,k)-zlay(ig,k-1))
+           sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
+                                                          enddo
+         enddo
+                                                          do ig=1,ngrid
+         l0(ig)=0.2*sqz(ig)/sq(ig)
+c        l0(ig)=30.
+                                                          enddo
+c      print*,'ITER=',iter,'  L0=',l0
+
+      enddo
+
+c     print*,'Fin de l initialisation de q2 et l0'
+
+      endif ! first
+
+c====================================================================
+c  Calcul de la longueur de melange.
+c====================================================================
+
+c   Mise a jour de l0
+                                                          do ig=1,ngrid
+      sq(ig)=1.e-10
+      sqz(ig)=1.e-10
+                                                          enddo
+      do k=2,nlay-1
+                                                          do ig=1,ngrid
+        zq=sqrt(q2(ig,k))
+        sqz(ig)=sqz(ig)+zq*zlev(ig,k)*(zlay(ig,k)-zlay(ig,k-1))
+        sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
+                                                          enddo
+      enddo
+                                                          do ig=1,ngrid
+      l0(ig)=0.2*sqz(ig)/sq(ig)
+c        l0(ig)=30.
+                                                          enddo
+c      print*,'ITER=',iter,'  L0=',l0
+c   calcul de l(z)
+      do k=2,nlay
+                                                          do ig=1,ngrid
+         l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
+         if(first) then
+           q2(ig,k)=l(ig,k)**2*zz(ig,k)
+         endif
+                                                          enddo
+      enddo
+
+c====================================================================
+c   Yamada 2.0
+c====================================================================
+      if (iflag_pbl.eq.6) then
+
+      do k=2,nlay
+                                                          do ig=1,ngrid
+         q2(ig,k)=l(ig,k)**2*zz(ig,k)
+                                                          enddo
+      enddo
+
+
+      else if (iflag_pbl.eq.7) then
+c====================================================================
+c   Yamada 2.Fournier
+c====================================================================
+
+c  Calcul de l,  km, au pas precedent
+      do k=2,nlay
+                                                          do ig=1,ngrid
+c        print*,'SMML=',sm(ig,k),l(ig,k)
+         delta(ig,k)=q2(ig,k)/(l(ig,k)**2*sm(ig,k))
+         kmpre(ig,k)=l(ig,k)*sqrt(q2(ig,k))*sm(ig,k)
+         mpre(ig,k)=sqrt(m2(ig,k))
+c        print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
+                                                          enddo
+      enddo
+
+      do k=2,nlay-1
+                                                          do ig=1,ngrid
+        m2cstat=max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1,1.e-12)
+        mcstat=sqrt(m2cstat)
+
+c        print*,'M2 L=',k,mpre(ig,k),mcstat
+c
+c  -----{puis on ecrit la valeur de q qui annule l equation de m
+c        supposee en q3}
+c
+        IF (k.eq.2) THEN
+          kmcstat=1.E+0 / mcstat
+     &    *( unsdz(ig,k)*kmpre(ig,k+1)
+     &                        *mpre(ig,k+1)
+     &      +unsdz(ig,k-1)
+     &              *cd(ig)
+     &              *( sqrt(u(ig,3)**2+v(ig,3)**2)
+     &                -mcstat/unsdzdec(ig,k)
+     &                -mpre(ig,k+1)/unsdzdec(ig,k+1) )**2)
+     &      /( unsdz(ig,k)+unsdz(ig,k-1) )
+        ELSE
+          kmcstat=1.E+0 / mcstat
+     &    *( unsdz(ig,k)*kmpre(ig,k+1)
+     &                        *mpre(ig,k+1)
+     &      +unsdz(ig,k-1)*kmpre(ig,k-1)
+     &                          *mpre(ig,k-1) )
+     &      /( unsdz(ig,k)+unsdz(ig,k-1) )
+        ENDIF
+c       print*,'T2 L=',k,tmp2
+        tmp2=kmcstat
+     &      /( sm(ig,k)/q2(ig,k) )
+     &      /l(ig,k)
+!       q2(ig,k)=max(tmp2,1.E-10)**(2./3.)
+! MARS
+        q2(ig,k)=max(tmp2,q2min)**(2./3.)
+c       print*,'Q2 L=',k,q2(ig,k)
+c
+                                                          enddo
+      enddo
+
+      else if (iflag_pbl.ge.8) then
+c====================================================================
+c   Yamada 2.5 a la Didi
+c====================================================================
+
+
+c  Calcul de l,  km, au pas precedent
+      do k=2,nlay
+                                                          do ig=1,ngrid
+c        print*,'SMML=',sm(ig,k),l(ig,k)
+         delta(ig,k)=q2(ig,k)/(l(ig,k)**2*sm(ig,k))
+         if (delta(ig,k).lt.1.e-20) then
+c     print*,'ATTENTION   L=',k,'   Delta=',delta(ig,k)
+            delta(ig,k)=1.e-20
+         endif
+         km(ig,k)=l(ig,k)*sqrt(q2(ig,k))*sm(ig,k)
+         aa0=
+     s   (m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)
+         aa1=
+     s   (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
+c abder      print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)
+         aa(ig,k)=aa1*dt/(delta(ig,k)*l(ig,k))
+c     print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
+         qpre=sqrt(q2(ig,k))
+         if (iflag_pbl.eq.8 ) then
+            if (aa(ig,k).gt.0.) then
+               q2(ig,k)=(qpre+aa(ig,k)*qpre*qpre)**2
+            else
+               q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2
+            endif
+         else ! iflag_pbl=9
+            if (aa(ig,k)*qpre.gt.0.9) then
+               q2(ig,k)=(qpre*10.)**2
+            else
+               q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2
+            endif
+         endif
+!         q2(ig,k)=min(max(q2(ig,k),1.e-10),1.e4)
+! MARS
+         q2(ig,k)=min(max(q2(ig,k),q2min),q2max)
+c     print*,'Q2 L=',k,q2(ig,k),qpre*qpre
+                                                          enddo
+      enddo
+
+!MARS
+!      q2(:,nlay+1)=q2min
+! OU ALORS MARS
+      q2(:,nlay+1)=q2(:,nlay)
+      endif ! Fin du cas 8
+
+c     print*,'OK8'
+
+c====================================================================
+c   Calcul des coefficients de m�ange
+c====================================================================
+      do k=2,nlay
+c     print*,'k=',k
+                                                          do ig=1,ngrid
+cabde      print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)
+         zq=sqrt(q2(ig,k))
+         km(ig,k)=l(ig,k)*zq*sm(ig,k)
+         kn(ig,k)=km(ig,k)*alpha(ig,k)
+! MARS
+!         km(ig,k)=max(l(ig,k)*zq*sm(ig,k),kmmin)
+!         kn(ig,k)=max(km(ig,k)*alpha(ig,k),knmin)
+         kq(ig,k)=l(ig,k)*zq*0.2
+c       print*,'KML=',km(ig,k),kn(ig,k)
+                                                          enddo
+      enddo
+! MARS
+!      km(:,nlay+1)=kmmin
+!      kn(:,nlay+1)=knmin
+!      kq(:,nlay+1)=knmin
+! OU ALORS MARS
+      km(:,nlay+1)=km(:,nlay)
+      kn(:,nlay+1)=kn(:,nlay)
+      kq(:,nlay+1)=kq(:,nlay)
+
+! Transport diffusif vertical de la TKE.
+      if (iflag_pbl.ge.9) then
+!       print*,'YAMADA VDIF'
+        q2(:,1)=q2(:,2)
+        call vdif_q2(dt,g,rconst,ngrid,nlay,plev,temp,kq,q2)
+      endif
+
+c   Traitement des cas noctrunes avec l'introduction d'une longueur
+c   minimale.
+
+c====================================================================
+c   Traitement particulier pour les cas tres stables.
+c   D apres Holtslag Boville.
+
+                                                          do ig=1,ngrid
+      coriol(ig)=1.e-4
+      pblhmin(ig)=0.07*ustar(ig)/max(abs(coriol(ig)),2.546e-5)
+                                                          enddo
+
+!      print*,'pblhmin ',pblhmin
+CTest a remettre 21 11 02
+c test abd 13 05 02      if(0.eq.1) then
+      if(1.eq.1) then
+      do k=2,nlay
+         do ig=1,ngrid
+            if (teta(ig,2).gt.teta(ig,1)) then
+               qmin=ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
+               kmin=kap*zlev(ig,k)*qmin
+            else
+               kmin=-1. ! kmin n'est utilise que pour les SL stables.
+            endif 
+            if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
+c               print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
+c     s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
+               kn(ig,k)=kmin
+               km(ig,k)=kmin
+               kq(ig,k)=kmin
+c   la longueur de melange est suposee etre l= kap z
+c   K=l q Sm d ou q2=(K/l Sm)**2
+               q2(ig,k)=(qmin/sm(ig,k))**2
+            endif
+         enddo
+      enddo
+      endif
+
+c   Diagnostique pour stokage
+
+      if(1.eq.0)then
+      rino=rif
+      smyam(1:ngrid,1)=0.
+      styam(1:ngrid,1)=0.
+      lyam(1:ngrid,1)=0.
+      knyam(1:ngrid,1)=0.
+      w2yam(1:ngrid,1)=0.
+      t2yam(1:ngrid,1)=0.
+
+      smyam(1:ngrid,2:nlay)=sm(1:ngrid,2:nlay)
+      styam(1:ngrid,2:nlay)=sm(1:ngrid,2:nlay)*alpha(1:ngrid,2:nlay)
+      lyam(1:ngrid,2:nlay)=l(1:ngrid,2:nlay)
+      knyam(1:ngrid,2:nlay)=kn(1:ngrid,2:nlay)
+
+c   Estimations de w'2 et T'2 d apres Abdela et McFarlane
+
+      w2yam(1:ngrid,2:nlay)=q2(1:ngrid,2:nlay)*0.24
+     s    +lyam(1:ngrid,2:nlay)*5.17*kn(1:ngrid,2:nlay)
+     s    *n2(1:ngrid,2:nlay)/sqrt(q2(1:ngrid,2:nlay))
+
+      t2yam(1:ngrid,2:nlay)=9.1*kn(1:ngrid,2:nlay)
+     s    *dtetadz(1:ngrid,2:nlay)**2
+     s    /sqrt(q2(1:ngrid,2:nlay))*lyam(1:ngrid,2:nlay)
+      endif
+
+c     print*,'OKFIN'
+      first=.false.
+      return
+      end
+      SUBROUTINE vdif_q2(timestep,gravity,rconst,ngrid,nlay,plev,temp,
+     &  kmy,q2)
+      IMPLICIT NONE
+c.......................................................................
+#include "dimensions.h"
+#include "dimphys.h"
+#include "tracer.h"
+#include "callkeys.h"
+c.......................................................................
+c
+c dt : pas de temps
+
+      real plev(ngrid,nlay+1)
+      real temp(ngrid,nlay)
+      real timestep
+      real gravity,rconst
+      real kstar(ngrid,nlay+1),zz
+      real kmy(ngrid,nlay+1)
+      real q2(ngrid,nlay+1)
+      real deltap(ngrid,nlay+1)
+      real denom(ngrid,nlay+1),alpha(ngrid,nlay+1),beta(ngrid,nlay+1)
+      integer ngrid,nlay
+
+      integer i,k
+
+! 	print*,'RD=',rconst
+      do k=1,nlay
+         do i=1,ngrid
+c test
+!       print*,'i,k',i,k
+! 	print*,'temp(i,k)=',temp(i,k)
+! 	print*,'(plev(i,k)-plev(i,k+1))=',plev(i,k),plev(i,k+1)
+            zz=(plev(i,k)+plev(i,k+1))*gravity/(rconst*temp(i,k))
+            kstar(i,k)=0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz
+     s      /(plev(i,k)-plev(i,k+1))*timestep
+         enddo
+      enddo
+
+      do k=2,nlay
+         do i=1,ngrid
+            deltap(i,k)=0.5*(plev(i,k-1)-plev(i,k+1))
+         enddo
+      enddo
+      do i=1,ngrid
+         deltap(i,1)=0.5*(plev(i,1)-plev(i,2))
+         deltap(i,nlay+1)=0.5*(plev(i,nlay)-plev(i,nlay+1))
+         denom(i,nlay+1)=deltap(i,nlay+1)+kstar(i,nlay)
+         alpha(i,nlay+1)=deltap(i,nlay+1)*q2(i,nlay+1)/denom(i,nlay+1)
+         beta(i,nlay+1)=kstar(i,nlay)/denom(i,nlay+1)
+      enddo
+
+      do k=nlay,2,-1
+         do i=1,ngrid
+            denom(i,k)=deltap(i,k)+(1.-beta(i,k+1))*
+     s      kstar(i,k)+kstar(i,k-1)
+c   correction d un bug 10 01 2001
+            alpha(i,k)=(q2(i,k)*deltap(i,k)
+     s      +kstar(i,k)*alpha(i,k+1))/denom(i,k)
+            beta(i,k)=kstar(i,k-1)/denom(i,k)
+         enddo
+      enddo
+
+c  Si on recalcule q2(1)
+      if(1.eq.0) then
+      do i=1,ngrid
+         denom(i,1)=deltap(i,1)+(1-beta(i,2))*kstar(i,1)
+         q2(i,1)=(q2(i,1)*deltap(i,1)
+     s      +kstar(i,1)*alpha(i,2))/denom(i,1)
+      enddo
+      endif
+c   sinon, on peut sauter cette boucle...
+
+      do k=2,nlay+1
+         do i=1,ngrid
+            q2(i,k)=alpha(i,k)+beta(i,k)*q2(i,k-1)
+         enddo
+      enddo
+
+      return
+      end
+      SUBROUTINE vdif_q2e(timestep,gravity,rconst,ngrid,nlay,
+     &   plev,temp,kmy,q2)
+      IMPLICIT NONE
+c.......................................................................
+#include "dimensions.h"
+#include "dimphys.h"
+#include "tracer.h"
+#include "callkeys.h"
+c.......................................................................
+c
+c dt : pas de temps
+
+      real plev(ngrid,nlay+1)
+      real temp(ngrid,nlay)
+      real timestep
+      real gravity,rconst
+      real kstar(ngrid,nlay+1),zz
+      real kmy(ngrid,nlay+1)
+      real q2(ngrid,nlay+1)
+      real deltap(ngrid,nlay+1)
+      real denom(ngrid,nlay+1),alpha(ngrid,nlay+1),beta(ngrid,nlay+1)
+      integer ngrid,nlay
+
+      integer i,k
+
+      do k=1,nlay
+         do i=1,ngrid
+            zz=(plev(i,k)+plev(i,k+1))*gravity/(rconst*temp(i,k))
+            kstar(i,k)=0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz
+     s      /(plev(i,k)-plev(i,k+1))*timestep
+         enddo
+      enddo
+
+      do k=2,nlay
+         do i=1,ngrid
+            deltap(i,k)=0.5*(plev(i,k-1)-plev(i,k+1))
+         enddo
+      enddo
+      do i=1,ngrid
+         deltap(i,1)=0.5*(plev(i,1)-plev(i,2))
+         deltap(i,nlay+1)=0.5*(plev(i,nlay)-plev(i,nlay+1))
+      enddo
+
+      do k=nlay,2,-1
+         do i=1,ngrid
+            q2(i,k)=q2(i,k)+
+     s      ( kstar(i,k)*(q2(i,k+1)-q2(i,k))
+     s       -kstar(i,k-1)*(q2(i,k)-q2(i,k-1)) )
+     s      /deltap(i,k)
+         enddo
+      enddo
+
+      do i=1,ngrid
+         q2(i,1)=q2(i,1)+
+     s   ( kstar(i,1)*(q2(i,2)-q2(i,1))
+     s                                      )
+     s   /deltap(i,1)
+         q2(i,nlay+1)=q2(i,nlay+1)+
+     s   ( 
+     s    -kstar(i,nlay)*(q2(i,nlay+1)-q2(i,nlay)) )
+     s   /deltap(i,nlay+1)
+      enddo
+
+      return
+      end
