Index: trunk/LMDZ.MARS/libf/phymars/physiq.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 524)
+++ trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 528)
@@ -340,6 +340,6 @@
 
       REAL lmax_th_out(ngridmx),zmax_th(ngridmx)
-      REAL, SAVE :: wmax_th(ngridmx)
-      REAL hfmax_th(ngridmx)
+      REAL, SAVE :: wstar(ngridmx)
+      REAL, SAVE :: hfmax_th(ngridmx)
       REAL pdu_th(ngridmx,nlayermx),pdv_th(ngridmx,nlayermx)
       REAL pdt_th(ngridmx,nlayermx),pdq_th(ngridmx,nlayermx,nqmx)
@@ -350,5 +350,5 @@
       REAL z_out                          ! height of interpolation between z0 and z1 [meters]
       REAL ustar(ngridmx),tstar(ngridmx)  ! friction velocity and friction potential temp
-      REAL L_mo(ngridmx)
+      REAL L_mo(ngridmx),wstarpbl(ngridmx),vhf(ngridmx),vvv(ngridmx)
       REAL zu2(ngridmx)
 c=======================================================================
@@ -367,5 +367,5 @@
          fluxrad(:)=0
 
-         wmax_th(:)=0.
+         wstar(:)=0.
 
 c        read startfi 
@@ -577,4 +577,6 @@
                  enddo
               endif
+           else
+              zdtnlte(:,:)=0.
            endif
 
@@ -754,6 +756,10 @@
           DO ig=1, ngridmx
              IF (zh(ig,1) .lt. tsurf(ig)) THEN
-               wmax_th(ig)=1.
-             ENDIF        
+               wstar(ig)=1.
+               hfmax_th(ig)=0.2
+             ELSE
+               wstar(ig)=0.
+               hfmax_th(ig)=0.
+             ENDIF      
           ENDDO
         ENDIF
@@ -771,5 +777,5 @@
      $        zdum1,zdum2,zdh,pdq,zflubid,
      $        zdudif,zdvdif,zdhdif,zdtsdif,q2,
-     &        zdqdif,zdqsdif,wmax_th,zcdv,zcdh)
+     &        zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th)
 
 #ifdef MESOSCALE
@@ -832,5 +838,5 @@
      $ pplay,pplev,pphi,zpopsk,
      $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
-     $ dtke_th,hfmax_th,wmax_th)
+     $ dtke_th,hfmax_th,wstar)
  
          DO l=1,nlayer
@@ -861,5 +867,6 @@
         else   !of if calltherm
         lmax_th(:)=0
-        wmax_th(:)=0.
+        wstar(:)=0.
+        hfmax_th(:)=0.
         lmax_th_out(:)=0.
         end if
@@ -913,10 +920,10 @@
 c   -------------------------------------------
 
-#ifdef MESOSCALE
-      !!! get the actual co2 seasonal cap from Titus observations
-      CALL geticecover( ngrid, 180.*zls/pi,
+      IF (tituscap) THEN
+         !!! get the actual co2 seasonal cap from Titus observations
+         CALL geticecover( ngrid, 180.*zls/pi,
      .                  180.*long/pi, 180.*lati/pi, co2ice )
-      co2ice = co2ice * 10000.
-#endif
+         co2ice = co2ice * 10000.
+      ENDIF
 
       IF (callcond) THEN
@@ -964,5 +971,4 @@
 c        ----------------------------------------
          IF (water) THEN
-         
 
            call watercloud(ngrid,nlayer,ptimestep,
@@ -979,5 +985,5 @@
            ENDDO
            endif
-           
+
 ! increment water vapour and ice atmospheric tracers tendencies
            IF (water) THEN
@@ -1170,5 +1176,4 @@
         ENDDO    ! (iq)
 
-
       endif !  of if (tracer) 
 
@@ -1244,5 +1249,5 @@
       ENDIF  ! of IF (tracer.AND.water.AND.(ngridmx.NE.1))
 
-
+c
 c   9.2 Compute soil temperatures and subsurface heat flux:
 c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -1303,5 +1308,5 @@
        DO ig=1,ngridmx
           DO l=1,nlayermx
-              zh(ig,l) = zt(ig,l)*(zplay(ig,l)/zplev(ig,1))**rcp
+              zh(ig,l) = zt(ig,l)*(zplev(ig,1)/zplay(ig,l))**rcp
           ENDDO
        ENDDO
@@ -1595,4 +1600,5 @@
 #include "fill_save.inc"
 #else
+#ifndef MESOINI
 
 c        ==========================================================
@@ -1647,8 +1653,16 @@
              endif
 
-#ifdef MESOINI
+#else
+     !!! this is to ensure correct initialisation of mesoscale model
+        call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
+     &                  tsurf)
+        call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps) 
+        call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
+     &                  co2ice)
+        call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
+        call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
+        call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
         call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
      &                  emis)
-        call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
         call WRITEDIAGFI(ngrid,"tsoil","Soil temperature",
      &                       "K",3,tsoil)
@@ -1833,27 +1847,9 @@
          z_out=0.
          if (calltherm .and. (z_out .gt. 0.)) then
-         call surflayer_interpol(ngrid,nlayer,z0,g,zzlay,zu,zv,wmax_th
-     &              ,tsurf,zh,z_out,Teta_out,u_out,ustar,tstar,L_mo)
-
-         zu2(:)=sqrt(zu(:,1)*zu(:,1)+zv(:,1)*zv(:,1))
-         call WRITEDIAGFI(ngridmx,'sqrt(zu2)',
-     &              'horizontal velocity norm','m/s',
-     &                         2,zu2)
-
-         call WRITEDIAGFI(ngridmx,'Teta_out',
-     &              'potential temperature at z_out','K',
-     &                         2,Teta_out)
-         call WRITEDIAGFI(ngridmx,'u_out',
-     &              'horizontal velocity norm at z_out','m/s',
-     &                         2,u_out)
-         call WRITEDIAGFI(ngridmx,'u*',
-     &              'friction velocity','m/s',
-     &                         2,ustar)
-         call WRITEDIAGFI(ngridmx,'teta*',
-     &              'friction potential temperature','K',
-     &                         2,tstar)
-         call WRITEDIAGFI(ngrid,'L',
-     &              'Monin Obukhov length','m',
-     &                         2,L_mo)
+
+         call pbl_parameters(ngrid,nlayer,z0,
+     & g,zzlay,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out,
+     & Teta_out,u_out,ustar,tstar,wstarpbl,L_mo,vhf,vvv)
+
          else
            if((.not. calltherm).and.(z_out .gt. 0.)) then
@@ -1892,7 +1888,7 @@
      &              'maximum TH heat flux','K.m/s',
      &                         2,hfmax_th)
-        call WRITEDIAGFI(ngridmx,'wmax_th',
+        call WRITEDIAGFI(ngridmx,'wstar',
      &              'maximum TH vertical velocity','m/s',
-     &                         2,wmax_th)
+     &                         2,wstar)
 
          endif
@@ -1944,24 +1940,9 @@
          z_out=0.
          if (calltherm .and. (z_out .gt. 0.)) then
-         call surflayer_interpol(ngrid,nlayer,z0,g,zzlay,zu,zv,wmax_th
-     &              ,tsurf,zh,z_out,Teta_out,u_out,ustar,tstar,L_mo)
-
-         zu2(:)=sqrt(zu(:,1)*zu(:,1)+zv(:,1)*zv(:,1))
-         call WRITEDIAGFI(ngridmx,'sqrt(zu2)',
-     &              'horizontal velocity norm','m/s',
-     &                         0,zu2)
-
-         call WRITEDIAGFI(ngridmx,'Teta_out',
-     &              'potential temperature at z_out','K',
-     &                         0,Teta_out)
-         call WRITEDIAGFI(ngridmx,'u_out',
-     &              'horizontal velocity norm at z_out','m/s',
-     &                         0,u_out)
-         call WRITEDIAGFI(ngridmx,'u*',
-     &              'friction velocity','m/s',
-     &                         0,ustar)
-         call WRITEDIAGFI(ngridmx,'teta*',
-     &              'friction potential temperature','K',
-     &                         0,tstar)
+
+         call pbl_parameters(ngrid,nlayer,z0,
+     & g,zzlay,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out,
+     & Teta_out,u_out,ustar,tstar,wstarpbl,L_mo,vhf,vvv)
+
          else
            if((.not. calltherm).and.(z_out .gt. 0.)) then
@@ -1977,10 +1958,13 @@
      &              'hauteur du thermique','point',
      &                         0,lmax_th_out)
+        call WRITEDIAGFI(ngridmx,'zmax_th',
+     &              'hauteur du thermique','m',
+     &                         0,zmax_th)
         call WRITEDIAGFI(ngridmx,'hfmax_th',
      &              'maximum TH heat flux','K.m/s',
      &                         0,hfmax_th)
-        call WRITEDIAGFI(ngridmx,'wmax_th',
+        call WRITEDIAGFI(ngridmx,'wstar',
      &              'maximum TH vertical velocity','m/s',
-     &                         0,wmax_th)
+     &                         0,wstar)
 
          co2col(:)=0.
Index: trunk/LMDZ.MARS/libf/phymars/surfini.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/surfini.F	(revision 524)
+++ trunk/LMDZ.MARS/libf/phymars/surfini.F	(revision 528)
@@ -64,7 +64,5 @@
         
 #ifdef MESOSCALE
-
       do ig=1,ngridmx
-
          !write(*,*) "all qsurf to zero. dirty."
          do iq=1,nqmx
@@ -83,7 +81,6 @@
                  watercaptag(ig)  = .false.
                  dryness(ig)      = 1.
-         endif 
-         
-      endif
+         endif  
+      enddo
 #else
        ! print*,'ngridmx',ngridmx
