Index: trunk/LMDZ.MARS/README
===================================================================
--- trunk/LMDZ.MARS/README	(revision 267)
+++ trunk/LMDZ.MARS/README	(revision 268)
@@ -850,2 +850,15 @@
 - Formula correction for the Reynolds in vdif_cd : will affect the thermal roughness length. Results obtained with the precedent formula (previous commit) probably underestimate by about 1K the temperature in the mixing layer.
 - Added a new subroutine to libf/phymars called surflayer_interpol.F. This routine can be called as a diagnostic for interpolation of horizontal velocity norm and temperature in the surface layer, in 1D or 3D. The output height is controlled by a parameter z_out, to be set in physiq.F. This might be later added to a .def file. z_out is set to 0. by default, resulting in no call to this subroutine.
+
+== 16/08/2011 == AC
+- Minor modification in physiq.F to pass Ch from vdifc to meso_inc_les
+- Major modification to the formulation of integrals in surflayer_interpol.F
+  Now stable for most cases. Some cases with highly negative Monin Obukhov length
+  remain to be explored.
+- Added gustiness to the Richardson computation in vdif_cd.F. Gustiness factor is for now of beta=1., after
+  several comparisons with LES aerodynamic conductances. May be subject to a minor change (+/- 0.1)
+  in the near future. (almost transparent for the user)
+- Minor modifications relative to variables in vdifc.F.
+- Esthetic change to calltherm_mars.F
+- Changed the definition for HFX computation in the LES in meso_inc/meso_inc_les.F. New definition yields
+  very similar results to old one and follows a strict definition of what HFX should be.
Index: trunk/LMDZ.MARS/libf/phymars/calltherm_mars.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/calltherm_mars.F90	(revision 267)
+++ trunk/LMDZ.MARS/libf/phymars/calltherm_mars.F90	(revision 268)
@@ -87,5 +87,9 @@
 !  ===================
  
-         r_aspect_thermals=2.
+         r_aspect_thermals=2.   ! ultimately conrols the amount of mass going through the thermals
+                                ! decreasing it increases the thermals effect. Tests at gcm resolution
+                                ! shows that too low values destabilize the model
+                                ! when changing this value, one should check that the surface layer model
+                                ! outputs the correct Cd*u and Ch*u through changing the gustiness coefficient beta
          nsplit_thermals=40
          call getin("nsplit_thermals",nsplit_thermals)
Index: trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_les.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_les.F	(revision 267)
+++ trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_les.F	(revision 268)
@@ -1,5 +1,12 @@
          DO ig=1,ngrid
           !! sensible heat flux in W/m2
-          hfx(ig) = zflubid(ig)-capcal(ig)*zdtsdif(ig)
+!          hfx(ig) = zflubid(ig)-capcal(ig)*zdtsdif(ig)
+
+! New SL parametrization, correct formulation for hfx :
+
+          hfx(ig) = (pplay(ig,1)/(r*pt(ig,1)))*cpp
+     &    *sqrt((pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1)))
+     &    *zcdh(ig)*(tsurf(ig)-zh(ig,1))
+
           !! u star in similarity theory in m/s
 !          ust(ig) = 0.4
Index: trunk/LMDZ.MARS/libf/phymars/physiq.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 267)
+++ trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 268)
@@ -327,5 +327,5 @@
       INTEGER lmax_th(ngridmx)
       REAL dtke_th(ngridmx,nlayermx+1)
-      REAL zcdv(ngridmx)
+      REAL zcdv(ngridmx), zcdh(ngridmx)
       REAL Teta_out(ngridmx),u_out(ngridmx)  ! Interpolated teta and u at z_out
       REAL z_out                          ! height of interpolation between z0 and z1
@@ -346,4 +346,6 @@
          call zerophys(ngrid*nlayer,dtrad)
          call zerophys(ngrid,fluxrad)
+
+         wmax_th(:)=0.
 
 c        read startfi 
@@ -700,5 +702,5 @@
      $        zdum1,zdum2,zdh,pdq,zflubid,
      $        zdudif,zdvdif,zdhdif,zdtsdif,q2,
-     &        zdqdif,zdqsdif,wmax_th,zcdv)
+     &        zdqdif,zdqsdif,wmax_th,zcdv,zcdh)
 
 #ifdef MESOSCALE
@@ -1679,4 +1681,40 @@
 ! THERMALS STUFF 1D
 
+         z_out=0.
+         if (calltherm .and. (z_out .gt. 0.)) then
+
+         call surflayer_interpol(ngrid,nlayer,z0,g,zzlay,zu,zv,wmax_th
+     &              ,tsurf,zt(:,:)*(zplay(:,:)/zplev(:,:))**rcp
+     &              ,z_out,Teta_out,u_out,ustar,tstar)
+
+         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)
+         else
+           if((.not. calltherm).and.(z_out .gt. 0.)) then
+            print*, 'WARNING : no interpolation in surface-layer :'
+            print*, 'Outputing surface-layer quantities without thermals
+     & does not make sense'
+           endif
+         endif
+
+! ---
+
+
+
          if(calltherm) then
 
@@ -1693,5 +1731,4 @@
      &                         0,wmax_th)
 
-
          co2col(:)=0.
          if (tracer) then
@@ -1712,4 +1749,6 @@
          call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
      &                  tsurf)
+!         call WRITEDIAGFI(ngrid,"u","u wind","m/s",1,zu)
+!         call WRITEDIAGFI(ngrid,"v","v wind","m/s",1,zv)
 
          call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)
Index: trunk/LMDZ.MARS/libf/phymars/surflayer_interpol.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/surflayer_interpol.F	(revision 267)
+++ trunk/LMDZ.MARS/libf/phymars/surflayer_interpol.F	(revision 268)
@@ -5,5 +5,5 @@
 !
 !   Subject: interpolation of temperature and velocity norm in the surface layer
-!   by recomputation of surface layer quantities (Richardson, Prandtl, u*, teta*)
+!   by recomputation of surface layer quantities (Richardson, Prandtl, Reynolds, u*, teta*)
 !   -------  
 !
@@ -42,6 +42,6 @@
       REAL, INTENT(IN) :: wmax(ngrid)
       REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay)
-      REAL, INTENT(INOUT) :: z_out                      ! output height (in m above surface)
-      REAL, INTENT(OUT) :: Teta_out(ngrid),u_out(ngrid) ! interpolated fields at z_out : potential temperature and norm(uv)
+      REAL, INTENT(INOUT) :: z_out                  ! output height (in m above surface)
+      REAL, INTENT(OUT) :: Teta_out(ngrid),u_out(ngrid)! interpolated fields at z_out : potential temperature and norm(uv)
       REAL, INTENT(OUT) :: ustar(ngrid), tstar(ngrid) ! u* and teta*
 
@@ -79,5 +79,5 @@
 ! For output :
 
-      REAL zu2                 ! Large-scale wind at first layer
+      REAL zu2(ngrid)                  ! Large-scale wind at first layer
       REAL L_mo(ngrid)                ! Monin-Obukhov length
 !-----------------------------------------------------------------------
@@ -88,16 +88,4 @@
       ustar(:)=0.
       reynolds(:)=0.
-
-! Original formulation :
-
-!      DO ig=1,ngrid
-!         z1=1.E+0 + pz(ig,1)/pz0(ig)
-!         zcd0=karman/log(z1)
-!         zcd0=zcd0*zcd0
-!         pcdv(ig)=zcd0
-!         pcdh(ig)=zcd0
-!      ENDDO
-
-!      print*,'old : cd,ch; ',pcdv,pcdh
 
 ! New formulation (AC) :
@@ -141,10 +129,22 @@
 ! Richardson number formulation proposed by D.E. England et al. (1995)
 
+!        IF((pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) .lt. 1.) 
+!     &            .and. (wmax(ig) .gt. 0.)) THEN
+          zu2(ig)=
+!     &  (MAX(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1),wmax(ig)**2))
+     &  ( pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) + wmax(ig)**2)
+!     &  pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1)
+!        ELSE
+!        zu2(ig)=pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1)
+!        ENDIF
+
           rib(ig) = (pg/ph(ig,1))
-     &      *pz(ig,1)*pz0(ig)/sqrt(pz(ig,1)*pz0t)
-!     &      *sqrt(pz(ig,1)*pz0(ig))
+!     &      *pz(ig,1)*pz0(ig)/sqrt(pz(ig,1)*pz0t)
+     &      *sqrt(pz(ig,1)*pz0(ig))
      &      *(((log(pz(ig,1)/pz0(ig)))**2)/(log(pz(ig,1)/pz0t)))
-     &      *(ph(ig,1)-pts(ig))
-     &  /(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1))
+     &      *(ph(ig,1)-pts(ig))/zu2(ig)
+
+!     &  /(MAX(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1),wmax(ig)**2))
+!     &  /( pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) + wmax(ig)**2)
 
          else
@@ -208,5 +208,5 @@
       ENDDO
 
-! Large-scale wind at first layer (with gustiness) and
+! Large-scale wind at first layer (without gustiness) and
 ! u* theta* computation
       DO ig=1,ngrid
@@ -216,8 +216,18 @@
            tstar(ig)=0.
          else
-           zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1)+(0.*wmax(ig))**2
-           ustar(ig)=karman*sqrt(fm(ig)*zu2)/(log(pz(ig,1)/pz0(ig)))
-           tstar(ig)=karman*fh(ig)*(ph(ig,1)-pts(ig))
-     &                   /(log(pz(ig,1)/pz0tcomp(ig))*sqrt(fm(ig)))
+
+!           ustar(ig)=karman*sqrt(fm(ig)*zu2(ig))/(log(pz(ig,1)/pz0(ig)))
+!           tstar(ig)=karman*fh(ig)*(ph(ig,1)-pts(ig))
+!     &                   /(log(pz(ig,1)/pz0tcomp(ig))*sqrt(fm(ig)))
+
+!simpler definition of u* and teta*, equivalent to the formula above :
+
+            ustar(ig)=sqrt(pcdv(ig))*sqrt(zu2(ig))                   
+            tstar(ig)=-pcdh(ig)*(pts(ig)-ph(ig,1))/sqrt(pcdv(ig))
+
+           if (tstar(ig) .lt. -50) then
+             print*, fh(ig),rib(ig),(ph(ig,1)-pts(ig))
+     &               ,log(pz(ig,1)/pz0tcomp(ig)),sqrt(fm(ig))
+           endif 
          endif
       ENDDO
@@ -233,4 +243,5 @@
       ENDDO
 
+      DO ig=1,ngrid
       IF(z_out .ge. pz(ig,1)) THEN
            z_out=1.
@@ -252,5 +263,5 @@
      &                  and computation is resumed'
       ENDIF
-
+      ENDDO
 
       print*, 'interpolation of u and teta at z_out=',z_out
@@ -265,23 +276,54 @@
      &                            *(z_out-pz0tcomp(ig))
         ELSEIF (L_mo(ig) .lt. 0.) THEN
-           u_out(ig)=(ustar(ig)/karman)*((
+
+          IF(L_mo(ig) .gt. -1000.) THEN
+           
+           u_out(ig)=(ustar(ig)/karman)*(
      &  2.*atan((1.-16.*z_out/L_mo(ig))**0.25)
-     &  +log((-1.+(1.-16.*z_out/L_mo(ig))**0.25)
-     &  /(-1.+(1.-16.*z_out/L_mo(ig))**0.25)))-(
-     &  2.*atan((1.-16.*pz0(ig)/L_mo(ig))**0.25)
-     &  +log((-1.+(1.-16.*pz0(ig)/L_mo(ig))**0.25)
-     &  /(-1.+(1.-16.*pz0(ig)/L_mo(ig))**0.25)))
+     & -2.*atan((1.-16.*pz0(ig)/L_mo(ig))**0.25)
+     & -2.*log(1.+(1.-16.*z_out/L_mo(ig))**0.25)
+     & +2.*log(1.+(1.-16.*pz0(ig)/L_mo(ig))**0.25)
+     & -   log(1.+sqrt(1.-16.*z_out/L_mo(ig)))
+     & +   log(1.+sqrt(1.-16.*pz0(ig)/L_mo(ig)))
+     & +   log(z_out/pz0(ig))
      &                                  )
-           Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*((
-     &  log((-1.+sqrt(1.-16.*z_out/L_mo(ig)))
-     &  /(1.+sqrt(1.-16.*z_out/L_mo(ig)))))-(
-     &  log((-1.+sqrt(1.-16.*pz0tcomp(ig)/L_mo(ig)))
-     &  /(1.+sqrt(1.-16.*pz0tcomp(ig)/L_mo(ig)))))
+
+           Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*(
+     &  2.*log(1.+sqrt(1.-16.*pz0tcomp(ig)/L_mo(ig)))
+     & -2.*log(1.+sqrt(1.-16.*z_out/L_mo(ig)))
+     & +   log(z_out/pz0tcomp(ig))
      &                                                        )
+
+          ELSE
+
+! We have to treat the case where L is very large and negative,
+! corresponding to a very nearly stable atmosphere (but not quite !)
+! i.e. teta* <0 and almost zero. We choose the cutoff
+! corresponding to L_mo < -1000 and use a 3rd order taylor expansion. The difference
+! between the two expression for z/L = -1/1000 is -1.54324*10^-9
+! (we do that to avoid using r*4 precision, otherwise, we get -inf values)          
+
+           u_out(ig)=(ustar(ig)/karman)*(
+     &     (4./L_mo(ig))*(z_out-pz0(ig))
+     &   + (20./(L_mo(ig))**2)*(z_out**2-pz0(ig)**2)
+     &   + (160./(L_mo(ig))**3)*(z_out**3-pz0(ig)**3)
+     &   + log(z_out/pz0(ig))
+     &                                   )
+
+           Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*(
+     &     (8./L_mo(ig))*(z_out-pz0tcomp(ig))
+     &   + (48./(L_mo(ig))**2)*(z_out**2-pz0tcomp(ig)**2)
+     &   + (1280./(3.*(L_mo(ig))**3))*(z_out**3-pz0tcomp(ig)**3)
+     &   + log(z_out/pz0tcomp(ig))
+     &                                                           )
+
+          ENDIF
         ELSE
            u_out(ig)=0.
-           Teta_out(ig)=0.
+           Teta_out(ig)=pts(ig)
         ENDIF
       ENDDO
+
+! Usefull diagnostics for the interpolation routine :
 
          call WRITEDIAGFI(ngrid,'L',
Index: trunk/LMDZ.MARS/libf/phymars/vdif_cd.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/vdif_cd.F	(revision 267)
+++ trunk/LMDZ.MARS/libf/phymars/vdif_cd.F	(revision 268)
@@ -1,4 +1,4 @@
       SUBROUTINE vdif_cd(ngrid,nlay,pz0,
-     &           pg,pz,pu,pv,pts,ph,pcdv,pcdh)
+     &           pg,pz,pu,pv,wmax,pts,ph,pcdv,pcdh)
       IMPLICIT NONE
 c=======================================================================
@@ -45,4 +45,5 @@
       REAL, INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
       REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay)
+      REAL, INTENT(IN) :: wmax(ngrid)
       REAL, INTENT(OUT) :: pcdv(ngrid),pcdh(ngrid) ! momentum and heat drag coefficient
 
@@ -62,4 +63,5 @@
 
       REAL rib(ngrid)  ! Bulk Richardson number
+      REAL rig(ngrid)  ! Gradient Richardson number
       REAL fm(ngrid) ! stability function for momentum
       REAL fh(ngrid) ! stability function for heat
@@ -79,4 +81,5 @@
       REAL ite
       REAL residual
+      REAL zu2
 c-----------------------------------------------------------------------
 c   couche de surface:
@@ -137,9 +140,22 @@
 c Richardson number formulation proposed by D.E. England et al. (1995)
 
+!         zu2=MAX(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1),0.25*wmax(ig)**2)
+!         zu2=pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1)
+!         zu2=MAX(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1),wmax(ig)**2)
+         zu2= pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) + (1.*wmax(ig))**2
+
+               ! we add the wmax to simulate
+               ! bulk Ri changes due to subgrid wind feeding the thermals
+
+!          rig(ig) = (pg/ph(ig,1))*((pz(ig,1)-pz0(ig))**2
+!     &         /(pz(ig,1)-pz0t))*(ph(ig,1)-pts(ig))
+!     &         /zu2
+
           rib(ig) = (pg/ph(ig,1))
-     &      *pz(ig,1)*pz0(ig)/sqrt(pz(ig,1)*pz0t)
+!     &      *pz(ig,1)*pz0(ig)/sqrt(pz(ig,1)*pz0t)
+     &      *sqrt(pz(ig,1)*pz0(ig))
      &      *(((log(pz(ig,1)/pz0(ig)))**2)/(log(pz(ig,1)/pz0t)))
      &      *(ph(ig,1)-pts(ig))
-     &  /(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1))
+     &  /zu2
 
          else
@@ -204,7 +220,10 @@
 ! Some useful diagnostics :
 
-!        call WRITEDIAGFI(ngrid,'Ri',
-!     &              'Richardson nb','no units',
+!        call WRITEDIAGFI(ngrid,'RiB',
+!     &              'Bulk Richardson nb','no units',
 !     &                         2,rib)
+!                call WRITEDIAGFI(ngrid,'RiG',
+!     &              'Grad Richardson nb','no units',
+!     &                         2,rig)
 !        call WRITEDIAGFI(ngrid,'Pr',
 !     &              'Prandtl nb','no units',
Index: trunk/LMDZ.MARS/libf/phymars/vdifc.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 267)
+++ trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 268)
@@ -5,5 +5,5 @@
      $                pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf,
      $                pdudif,pdvdif,pdhdif,pdtsrf,pq2,
-     $                pdqdif,pdqsdif,wmax,zcdv_true)
+     $                pdqdif,pdqsdif,wmax,zcdv_true,zcdh_true)
       IMPLICIT NONE
 
@@ -78,5 +78,5 @@
       REAL zkv(ngridmx,nlayermx+1),zkh(ngridmx,nlayermx+1)
       REAL zcdv(ngridmx),zcdh(ngridmx)
-      REAL zcdv_true(ngridmx),zcdh_true(ngridmx)
+      REAL zcdv_true(ngridmx),zcdh_true(ngridmx)    ! drag coeff are used by the LES to recompute u* and hfx
       REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
       REAL zh(ngridmx,nlayermx)
@@ -232,5 +232,5 @@
 c       ---------------------
 
-      CALL vdif_cd(ngrid,nlay,pz0,g,pzlay,pu,pv,ptsrf,ph
+      CALL vdif_cd(ngrid,nlay,pz0,g,pzlay,pu,pv,wmax,ptsrf,ph
      &             ,zcdv_true,zcdh_true)
 
@@ -238,26 +238,32 @@
 
         zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1)
-     &             +(1.2*wmax(ig))**2        !subgrid gustiness is used to enhance surface flux only, and not u*,t* computations
+     &             +(1.*wmax(ig))**2        ! subgrid gustiness added by quadratic interpolation with a coeff beta, here beta=1. (tuned from
+                                            ! LES comparisons. This parameter is linked to the thermals settings)
         
-        zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
-        zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
-
+        zcdv(ig)=zcdv_true(ig)*sqrt(zu2)     ! 1 / bulk aerodynamic momentum conductance 
+        zcdh(ig)=zcdh_true(ig)*sqrt(zu2)     ! 1 / bulk aerodynamic heat conductance
+                                             ! these are the quantities to be looked at when comparing surface layers of different models
       ENDDO
 
 ! Some usefull diagnostics for the new surface layer parametrization :
 
-!        call WRITEDIAGFI(ngridmx,'Cd',
+!         call WRITEDIAGFI(ngridmx,'pcdv',
 !     &              'momentum drag','no units',
 !     &                         2,zcdv_true)
-!        call WRITEDIAGFI(ngridmx,'Ch',
+!        call WRITEDIAGFI(ngridmx,'pcdh',
 !     &              'heat drag','no units',
 !     &                         2,zcdh_true)
-!        call WRITEDIAGFI(ngridmx,'u*2/u',
-!     &              'friction velocity','m/s',
+!        call WRITEDIAGFI(ngridmx,'u*',
+!    &              'friction velocity','m/s',
+!     &                         2,sqrt(zcdv_true(:))*sqrt(zu2))
+!        call WRITEDIAGFI(ngridmx,'T*',
+!     &              'friction temperature','K',
+!     &          2,-zcdh_true(:)*(ptsrf(:)-ph(:,1))/sqrt(zcdv_true(:)))
+!        call WRITEDIAGFI(ngridmx,'rm-1',
+!     &              'aerodyn momentum conductance','m/s',
 !     &                         2,zcdv)
-!        call WRITEDIAGFI(ngridmx,'u*T*/(T0-T)',
-!     &              'friction temperature','m/s',
+!        call WRITEDIAGFI(ngridmx,'rh-1',
+!     &              'aerodyn heat conductance','m/s',
 !     &                         2,zcdh)
-
 
 c    ** schema de diffusion turbulente dans la couche limite
