Index: LMDZ6/trunk/libf/phylmd/calbeta.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/calbeta.F90	(revision 3779)
+++ LMDZ6/trunk/libf/phylmd/calbeta.F90	(revision 3780)
@@ -7,5 +7,9 @@
   USE dimphy
   USE indice_sol_mod
+
   IMPLICIT none
+
+#include "flux_arp.h"
+
 !======================================================================
 ! Auteur(s): Z.X. Li (LMD/CNRS) (adaptation du GCM au LMD)
@@ -82,4 +86,11 @@
      ENDDO
   ENDIF
+
+  ! EV: when beta is prescribed for 1D cases:
+  IF (knon.EQ.1 .AND. ok_prescr_beta) THEN
+     DO i = 1, knon
+          vbeta(i)=betaevap
+      ENDDO
+  ENDIF
   
 END SUBROUTINE calbeta
Index: LMDZ6/trunk/libf/phylmd/change_srf_frac_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/change_srf_frac_mod.F90	(revision 3779)
+++ LMDZ6/trunk/libf/phylmd/change_srf_frac_mod.F90	(revision 3780)
@@ -183,5 +183,4 @@
            tsurf, alb_dir,alb_dif, ustar, u10m, v10m, pbl_tke)
 
-
     ELSE
        ! No modifcation should be done
Index: LMDZ6/trunk/libf/phylmd/dyn1d/1DUTILS.h
===================================================================
--- LMDZ6/trunk/libf/phylmd/dyn1d/1DUTILS.h	(revision 3779)
+++ LMDZ6/trunk/libf/phylmd/dyn1d/1DUTILS.h	(revision 3780)
@@ -233,4 +233,11 @@
        CALL getin('ok_flux_surf',ok_flux_surf)
 
+!Config  Key  = ok_forc_tsurf
+!Config  Desc = forcage ou non par la Ts
+!Config  Def  = false
+!Config  Help = forcage ou non par la Ts
+       ok_forc_tsurf=.false.
+       CALL getin('ok_forc_tsurf',ok_forc_tsurf)
+
 !Config  Key  = ok_prescr_ust
 !Config  Desc = ustar impose ou non
@@ -239,4 +246,12 @@
        ok_prescr_ust = .false.
        CALL getin('ok_prescr_ust',ok_prescr_ust)
+
+
+!Config  Key  = ok_prescr_beta
+!Config  Desc = betaevap impose ou non
+!Config  Def  = false
+!Config  Help = betaevap impose ou non
+       ok_prescr_beta = .false.
+       CALL getin('ok_prescr_beta',ok_prescr_beta)
 
 !Config  Key  = ok_old_disvert
@@ -280,5 +295,5 @@
 !Config  Desc = surface temperature
 !Config  Def  = 290.
-!Config  Help = not used if type_ts_forcing=1 in lmdz1d.F
+!Config  Help = surface temperature
        tsurf = 290.
        CALL getin('tsurf',tsurf)
@@ -297,4 +312,7 @@
        zsurf = 0.
        CALL getin('zsurf',zsurf)
+! EV pour accord avec format standard       
+       CALL getin('zorog',zsurf)
+
 
 !Config  Key  = rugos
@@ -359,4 +377,13 @@
        qsolinp = 1.
        CALL getin('qsolinp',qsolinp)
+
+
+
+!Config  Key  = betaevap
+!Config  Desc = beta for actual evaporation when prescribed
+!Config  Def  = 1.0
+!Config  Help = 
+       betaevap = 1.
+       CALL getin('betaevap',betaevap)      
 
 !Config  Key  = zpicinp
@@ -520,4 +547,5 @@
        CALL getin('forc_ustar',forc_ustar)
        IF (forc_ustar .EQ. 1) ok_prescr_ust=.true.
+
 
 !Config  Key  = nudging_u
@@ -1248,25 +1276,25 @@
       END
 
-!======================================================================
-       SUBROUTINE read_tsurf1d(knon,sst_out)
-
-! This subroutine specifies the surface temperature to be used in 1D simulations
-
-      USE dimphy, ONLY : klon
-
-      INTEGER, INTENT(IN)                  :: knon     ! nomber of points on compressed grid
-      REAL, DIMENSION(klon), INTENT(OUT)   :: sst_out  ! tsurf used to force the single-column model
-
-       INTEGER :: i
-! COMMON defined in lmdz1d.F:
-       real ts_cur
-       common /sst_forcing/ts_cur
-
-       DO i = 1, knon
-        sst_out(i) = ts_cur
-       ENDDO
-
-      END SUBROUTINE read_tsurf1d
-
+!!======================================================================
+!       SUBROUTINE read_tsurf1d(knon,sst_out)
+!
+!! This subroutine specifies the surface temperature to be used in 1D simulations
+!
+!      USE dimphy, ONLY : klon
+!
+!      INTEGER, INTENT(IN)                  :: knon     ! nomber of points on compressed grid
+!      REAL, DIMENSION(klon), INTENT(OUT)   :: sst_out  ! tsurf used to force the single-column model
+!
+!       INTEGER :: i
+!! COMMON defined in lmdz1d.F:
+!       real ts_cur
+!       common /sst_forcing/ts_cur
+
+!       DO i = 1, knon
+!        sst_out(i) = ts_cur
+!       ENDDO
+!
+!      END SUBROUTINE read_tsurf1d
+!
 !===============================================================
       subroutine advect_vert(llm,w,dt,q,plev)
Index: LMDZ6/trunk/libf/phylmd/dyn1d/1D_decl_cases.h
===================================================================
--- LMDZ6/trunk/libf/phylmd/dyn1d/1D_decl_cases.h	(revision 3779)
+++ LMDZ6/trunk/libf/phylmd/dyn1d/1D_decl_cases.h	(revision 3780)
@@ -34,10 +34,11 @@
         real w_mod(llm), t_mod(llm),q_mod(llm)
         real u_mod(llm),v_mod(llm), ht_mod(llm),vt_mod(llm),ug_mod(llm),vg_mod(llm)
-	real temp_nudg_mod(llm),qv_nudg_mod(llm),u_nudg_mod(llm),v_nudg_mod(llm)
+	    real temp_nudg_mod(llm),qv_nudg_mod(llm),u_nudg_mod(llm),v_nudg_mod(llm)
         real hq_mod(llm),vq_mod(llm),qv_mod(llm),ql_mod(llm),qt_mod(llm)
         real th_mod(llm)
 
-        real ts_cur
-        common /sst_forcing/ts_cur ! also in read_tsurf1d.F
+! EV comment these lines
+!        real ts_cur
+!        common /sst_forcing/ts_cur ! also in read_tsurf1d.F
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 ! Declarations specifiques au cas RICO
Index: LMDZ6/trunk/libf/phylmd/dyn1d/1D_interp_cases.h
===================================================================
--- LMDZ6/trunk/libf/phylmd/dyn1d/1D_interp_cases.h	(revision 3779)
+++ LMDZ6/trunk/libf/phylmd/dyn1d/1D_interp_cases.h	(revision 3780)
@@ -1,4 +1,4 @@
 
-         print*,'FORCING CASE forcing_case2'
+        print*,'FORCING CASE forcing_case2'
 !       print*,                                                             &
 !    & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/pdt_cas=',     &
@@ -28,6 +28,6 @@
      &       ,dth_prof_cas,hth_prof_cas,vth_prof_cas,lat_prof_cas                           &
      &       ,sens_prof_cas,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas,tke_prof_cas)
-
-             ts_cur = ts_prof_cas 
+! EV tg instead of ts_cur
+             tg = ts_prof_cas 
 !            psurf=plev_prof_cas(1)
              psurf=ps_prof_cas
Index: LMDZ6/trunk/libf/phylmd/dyn1d/1D_read_forc_cases.h
===================================================================
--- LMDZ6/trunk/libf/phylmd/dyn1d/1D_read_forc_cases.h	(revision 3779)
+++ LMDZ6/trunk/libf/phylmd/dyn1d/1D_read_forc_cases.h	(revision 3780)
@@ -70,7 +70,10 @@
 
 ! initial and boundary conditions :
-!      tsurf = ts_prof_cas
+!     tsurf = ts_prof_cas
       psurf = ps_prof_cas
-      ts_cur = ts_prof_cas 
+      !EV tg instead of ts_cur
+      tg = ts_prof_cas 
+      print*, 'tg=', tg
+
       do l = 1, llm
        temp(l) = t_mod_cas(l)
@@ -108,6 +111,6 @@
        IF (ok_prescr_ust) THEN
        ust=ustar_prof_cas
-       print *,'ust=',ust
        ENDIF
 
+
       endif !forcing_SCM
Index: LMDZ6/trunk/libf/phylmd/dyn1d/mod_1D_cases_read_std.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/dyn1d/mod_1D_cases_read_std.F90	(revision 3779)
+++ LMDZ6/trunk/libf/phylmd/dyn1d/mod_1D_cases_read_std.F90	(revision 3780)
@@ -92,5 +92,5 @@
       REAL, ALLOCATABLE :: time_val(:)
 
-      print*,'ON EST VRAIMENT LA'
+      print*,'ON EST VRAIMENT DASN MOD_1D_CASES_READ_STD'
       fich_cas='cas.nc'
       print*,'fich_cas ',fich_cas
@@ -924,4 +924,6 @@
 !      enddo
 
+!       print*, 'plev_prof_cas', plev_prof_cas
+!       print*, 'play', play
        do l = 1, llm
 
@@ -951,4 +953,5 @@
 
          frac = (plev_prof_cas(k2)-play(l))/(plev_prof_cas(k2)-plev_prof_cas(k1))
+         
          t_mod_cas(l)= t_prof_cas(k2) - frac*(t_prof_cas(k2)-t_prof_cas(k1))
          theta_mod_cas(l)= th_prof_cas(k2) - frac*(th_prof_cas(k2)-th_prof_cas(k1))
@@ -1075,4 +1078,5 @@
        enddo ! l
 
+
           return
           end SUBROUTINE interp2_case_vertical_std
Index: LMDZ6/trunk/libf/phylmd/dyn1d/old_1D_decl_cases.h
===================================================================
--- LMDZ6/trunk/libf/phylmd/dyn1d/old_1D_decl_cases.h	(revision 3779)
+++ LMDZ6/trunk/libf/phylmd/dyn1d/old_1D_decl_cases.h	(revision 3780)
@@ -37,6 +37,6 @@
         real th_mod(llm)
 
-        real ts_cur
-        common /sst_forcing/ts_cur ! also in read_tsurf1d.F
+        !real ts_cur
+        !common /sst_forcing/ts_cur ! also in read_tsurf1d.F
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 ! Declarations specifiques au cas RICO
Index: LMDZ6/trunk/libf/phylmd/dyn1d/old_1D_interp_cases.h
===================================================================
--- LMDZ6/trunk/libf/phylmd/dyn1d/old_1D_interp_cases.h	(revision 3779)
+++ LMDZ6/trunk/libf/phylmd/dyn1d/old_1D_interp_cases.h	(revision 3780)
@@ -62,6 +62,6 @@
      &             ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof    &
      &             ,ht_prof,vt_prof,hq_prof,vq_prof)
-
-        if (type_ts_forcing.eq.1) ts_cur = ts_prof ! SST used in read_tsurf1d
+! EV: tg instead of ts_cur
+        if (type_ts_forcing.eq.1) tg = ts_prof ! 
 
 ! vertical interpolation:
@@ -113,6 +113,6 @@
 !     print *,'llm l omega_profd',llm,l,omega_profd(l)
 !     enddo
-
-        if (type_ts_forcing.eq.1) ts_cur = tg_prof ! SST used in read_tsurf1d
+! EV tg instead of ts_cur
+        if (type_ts_forcing.eq.1) tg = tg_prof ! SST used 
 
 ! vertical interpolation:
@@ -206,6 +206,6 @@
      &             ,ug_gabls4,vg_gabls4,ht_gabls4,hq_gabls4,tg_gabls4                            &
      &             ,ug_profg,vg_profg,ht_profg,hq_profg,tg_profg)
-
-        if (type_ts_forcing.eq.1) ts_cur = tg_prof ! SST used in read_tsurf1d
+!EV tg instead of ts_cur
+        if (type_ts_forcing.eq.1) tg = tg_prof ! SST used
 
 ! vertical interpolation:
@@ -499,6 +499,6 @@
      &             ,nlev_sandu                                              &
      &             ,ts_sandu,ts_prof)
-
-        if (type_ts_forcing.eq.1) ts_cur = ts_prof ! SST used in read_tsurf1d
+! EV tg instead of ts_cur
+        if (type_ts_forcing.eq.1) tg = ts_prof ! SST used in read_tsurf1d
 
 ! vertical interpolation:
@@ -582,7 +582,6 @@
      &             ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof    &
      &             ,ufa_prof,vfa_prof)
-
-        if (type_ts_forcing.eq.1) ts_cur = ts_prof ! SST used in read_tsurf1d
-
+! EV tg instead of ts_cur
+        if (type_ts_forcing.eq.1) tg = ts_prof ! SST used 
 ! vertical interpolation:
       CALL interp_astex_vertical(play,nlev_astex,plev_profa                 &
@@ -675,6 +674,7 @@
      &       ,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas,lat_prof_cas               &
      &       ,sens_prof_cas,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas)
-
-             ts_cur = ts_prof_cas 
+! EV tg instead of ts_cur
+
+             tg = ts_prof_cas 
              psurf=plev_prof_cas(1)
 
@@ -850,6 +850,7 @@
      &       ,dth_prof_cas,hth_prof_cas,vth_prof_cas,lat_prof_cas                           &
      &       ,sens_prof_cas,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas,tke_prof_cas)
-
-             ts_cur = ts_prof_cas 
+! EV tg instead of ts_cur
+
+             tg = ts_prof_cas 
 !            psurf=plev_prof_cas(1)
              psurf=ps_prof_cas
Index: LMDZ6/trunk/libf/phylmd/dyn1d/old_1D_read_forc_cases.h
===================================================================
--- LMDZ6/trunk/libf/phylmd/dyn1d/old_1D_read_forc_cases.h	(revision 3779)
+++ LMDZ6/trunk/libf/phylmd/dyn1d/old_1D_read_forc_cases.h	(revision 3780)
@@ -875,6 +875,7 @@
 
 ! initial and boundary conditions :
-!      tsurf = ts_prof_cas
-      ts_cur = ts_prof_cas 
+!     tsurf = ts_prof_cas
+! EV tg instead of ts_cur
+      tg= ts_prof_cas 
       psurf=plev_prof_cas(1)
       write(*,*) 'SST initiale: ',tsurf
@@ -965,5 +966,6 @@
 ! initial and boundary conditions :
 !      tsurf = ts_prof_cas
-      ts_cur = ts_prof_cas 
+! EV tg instead of ts_cur
+      tg = ts_prof_cas 
       psurf=plev_prof_cas(1)
       write(*,*) 'SST initiale: ',tsurf
@@ -1063,5 +1065,7 @@
 ! initial and boundary conditions :
 !      tsurf = ts_prof_cas
-      ts_cur = ts_prof_cas 
+! EV tg instead of ts_cur
+
+      tg = ts_prof_cas 
       psurf=plev_prof_cas(1)
       write(*,*) 'SST initiale: ',tsurf
Index: LMDZ6/trunk/libf/phylmd/dyn1d/old_lmdz1d.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/dyn1d/old_lmdz1d.F90	(revision 3779)
+++ LMDZ6/trunk/libf/phylmd/dyn1d/old_lmdz1d.F90	(revision 3780)
@@ -728,5 +728,7 @@
 
 !Al1 pour SST forced, appell?? depuis ocean_forced_noice
-      ts_cur = tsurf ! SST used in read_tsurf1d
+! EV tg instead of ts_cur
+
+      tg = tsurf ! SST used in read_tsurf1d
 !=====================================================================
 ! Initialisation de la physique : 
@@ -791,4 +793,5 @@
 
         fder=0.
+        print *, 'snsrf', snsrf
         snsrf(1,:)=snowmass ! masse de neige des sous surface
         qsurfsrf(1,:)=qsurf ! humidite de l'air des sous surface
@@ -841,5 +844,4 @@
      end if
 
-
         print*,'nat_surf,pctsrf(1,is_oce),pctsrf(1,is_ter)',nat_surf         &
      &        ,pctsrf(1,is_oce),pctsrf(1,is_ter)
@@ -851,6 +853,8 @@
               ! 6 albedo, mais on peut quand meme tourner avec
               ! moins. Seules les 2 ou 4 premiers seront lus
+              print*, 'les albedos sont sont', albedo, falb_dir
         falb_dir=albedo
         falb_dif=albedo
+             print*, falb_dir
         rugoro=rugos 
         t_ancien(1,:)=temp(:)
@@ -913,8 +917,8 @@
         v_ancien(1,:)=v(:)
 
-u10m=0.
-v10m=0.
-ale_wake=0.
-ale_bl_stat=0.
+        u10m=0.
+        v10m=0.
+        ale_wake=0.
+        ale_bl_stat=0.
  
 !------------------------------------------------------------------------
Index: LMDZ6/trunk/libf/phylmd/dyn1d/scm.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/dyn1d/scm.F90	(revision 3779)
+++ LMDZ6/trunk/libf/phylmd/dyn1d/scm.F90	(revision 3780)
@@ -75,28 +75,4 @@
       real :: zcufi    = 1.
       real :: zcvfi    = 1.
-
-!-      real :: nat_surf
-!-      logical :: ok_flux_surf
-!-      real :: fsens
-!-      real :: flat
-!-      real :: tsurf
-!-      real :: rugos
-!-      real :: qsol(1:2)
-!-      real :: qsurf
-!-      real :: psurf
-!-      real :: zsurf
-!-      real :: albedo
-!-
-!-      real :: time     = 0.
-!-      real :: time_ini
-!-      real :: xlat 
-!-      real :: xlon 
-!-      real :: wtsurf 
-!-      real :: wqsurf 
-!-      real :: restart_runoff 
-!-      real :: xagesno 
-!-      real :: qsolinp 
-!-      real :: zpicinp 
-!-
       real :: fnday 
       real :: day, daytime 
@@ -141,6 +117,4 @@
         logical :: forcing_case2   = .false.
         logical :: forcing_SCM   = .false.
-        integer :: type_ts_forcing ! 0 = SST constant; 1 = SST read from a file
-!                                                            (cf read_tsurf1d.F)
 
 !flag forcings
@@ -148,4 +122,6 @@
         logical :: nudge_thermo=.false.
         logical :: cptadvw=.true.
+
+
 !=====================================================================
 ! DECLARATIONS FOR EACH CASE
@@ -248,5 +224,5 @@
 !
       integer :: it_end ! iteration number of the last call
-!Al1
+!Al1,plev,play,phi,phis,presnivs,
       integer ecrit_slab_oc !1=ecrit,-1=lit,0=no file
       data ecrit_slab_oc/-1/
@@ -278,4 +254,5 @@
       d_v_age(:)=0.
       
+      
 ! Initialization of Common turb_forcing
        dtime_frcg = 0.
@@ -290,5 +267,4 @@
 ! OPTIONS OF THE 1D SIMULATION (lmdz1d.def => unicol.def)
 !---------------------------------------------------------------------
-!Al1
         call conf_unicol
 !Al1 moves this gcssold var from common fcg_gcssold to 
@@ -296,5 +272,4 @@
 ! --------------------------------------------------------------------
         close(1)
-!Al1
         write(*,*) 'lmdz1d.def lu => unicol.def'
 
@@ -302,6 +277,7 @@
        year_ini_cas=1997
        ! It is possible that those parameters are run twice.
-
        ! A REVOIR : LIRE PEUT ETRE AN MOIS JOUR DIRECETEMENT
+
+
        call getin('anneeref',year_ini_cas)
        call getin('dayref',day_deb)
@@ -309,9 +285,8 @@
        call getin('time_ini',heure_ini_cas)
 
-        type_ts_forcing = 0
-        IF (nat_surf==0) type_ts_forcing=1 ! SST forcee sur OCEAN
-        print*,'NATURE DE LA SURFACE ',nat_surf
+      print*,'NATURE DE LA SURFACE ',nat_surf
 !
 ! Initialization of the logical switch for nudging
+
      jcode = iflag_nudge
      do i = 1,nudge_max
@@ -319,7 +294,7 @@
        jcode = jcode/10
      enddo
-!---------------------------------------------------------------------
+!-----------------------------------------------------------------------
 !  Definition of the run
-!---------------------------------------------------------------------
+!-----------------------------------------------------------------------
 
       call conf_gcm( 99, .TRUE. )
@@ -343,4 +318,6 @@
       allocate( phy_flic(year_len)) ! Fraction de glace
       phy_flic(:)=0.0
+
+
 !-----------------------------------------------------------------------
 !   Choix du calendrier
@@ -373,4 +350,6 @@
 !      Le numero du jour est dans "day". L heure est traitee separement.
 !      La date complete est dans "daytime" (l'unite est le jour).
+
+
       if (nday>0) then
          fnday=nday
@@ -409,5 +388,5 @@
 ! Initialization of dimensions, geometry and initial state
 !---------------------------------------------------------------------
-!      call init_phys_lmdz(1,1,llm,1,(/1/)) ! job now done via iniphysiq
+!     call init_phys_lmdz(1,1,llm,1,(/1/)) ! job now done via iniphysiq
 !     but we still need to initialize dimphy module (klon,klev,etc.)  here.
       call init_dimphy1D(1,llm)
@@ -446,10 +425,4 @@
 !!! Feedback forcing values for Gateaux differentiation (al1)
 !!!=====================================================================
-!!! Surface Planck forcing bracketing call radiation
-!!      surf_Planck = 0.
-!!      surf_Conv   = 0.
-!!      write(*,*) 'Gateaux-dif Planck,Conv:',surf_Planck,surf_Conv
-!!! a mettre dans le lmdz1d.def ou autre
-!!
 !!
       qsol = qsolinp
@@ -469,8 +442,4 @@
       ENDIF
       print*,'Flux sol ',fsens,flat
-!!      ok_flux_surf=.false.
-!!      fsens=-wtsurf*rcpd*rho(1)
-!!      flat=-wqsurf*rlvtt*rho(1)
-!!!!
 
 ! Vertical discretization and pressure levels at half and mid levels:
@@ -496,5 +465,5 @@
       plev =ap+bp*psurf
       play = 0.5*(plev(1:llm)+plev(2:llm+1))
-      zlay=-rd*300.*log(play/psurf)/rg ! moved after reading profiles
+      zlay=-rd*300.*log(play/psurf)/rg ! moved after reading profiles. 
 
       IF (forcing_type .eq. 59) THEN
@@ -527,8 +496,10 @@
       print*,'mxcalc=',mxcalc
 !     print*,'zlay=',zlay(mxcalc)
-      print*,'play=',play(mxcalc)
-
-!Al1 pour SST forced, appell?? depuis ocean_forced_noice
-      ts_cur = tsurf ! SST used in read_tsurf1d
+!      print*,'play=',play(mxcalc)
+
+!! When surface temperature is forced
+      tg= tsurf ! surface T used in read_tsurf1d
+
+
 !=====================================================================
 ! Initialisation de la physique : 
@@ -546,4 +517,6 @@
 ! airefi,zcufi,zcvfi initialises au debut de ce programme
 ! rday,ra,rg,rd,rcpd declares dans YOMCST.h et calcules dans suphel.F
+
+
       day_step = float(nsplit_phys)*day_step/float(iphysiq)
       write (*,*) 'Time step divided by nsplit_phys (=',nsplit_phys,')'
@@ -563,4 +536,5 @@
      ! e.g. for cell boundaries, which are meaningless in 1D; so pad these 
      ! with '0.' when necessary
+
       call iniphysiq(iim,jjm,llm, &
            1,comm_lmdz, &
@@ -653,4 +627,5 @@
               ! 6 albedo, mais on peut quand meme tourner avec
               ! moins. Seules les 2 ou 4 premiers seront lus
+              print*, 'les albedos sont', albedo
         falb_dir=albedo
         falb_dif=albedo
@@ -664,7 +639,8 @@
         prw_ancien = 0.
 !jyg<
-!!        pbl_tke(:,:,:)=1.e-8
+!!      pbl_tke(:,:,:)=1.e-8
         pbl_tke(:,:,:)=0.
-        pbl_tke(:,2,:)=1.e-2
+! EV: pourquoi????
+!        pbl_tke(:,2,:)=1.e-2
         PRINT *, ' pbl_tke dans lmdz1d '
         if (prt_level .ge. 5) then 
@@ -675,5 +651,4 @@
 
 !>jyg
-
         rain_fall=0.
         snow_fall=0.
@@ -715,8 +690,8 @@
         v_ancien(1,:)=v(:)
  
-u10m=0.
-v10m=0.
-ale_wake=0.
-ale_bl_stat=0.
+        u10m=0.
+        v10m=0.
+        ale_wake=0.
+        ale_bl_stat=0.
 
 !------------------------------------------------------------------------
@@ -738,5 +713,5 @@
 ! to be set at some arbitratry convenient values.
 !------------------------------------------------------------------------
-!Al1 =============== restart option ==========================
+!Al1 =============== restart option ======================================
         if (.not.restart) then
           iflag_pbl = 5
@@ -803,5 +778,5 @@
        print*,'plev,play,phi,phis,presnivs,u,v,temp,q,omega2'
        print*,'temp(1),q(1,1),u(1),v(1),plev(1),phis :'
-       print*,temp(1),q(1,1),u(1),v(1),plev(1),phis
+       print*,temp(1),q(1,1),u(1),v(1),plev(1),phis(1)
 ! raz for safety
        do l=1,llm
@@ -809,5 +784,5 @@
        enddo
       endif
-!Al1 ================  end restart =================================
+!======================  end restart =================================
       IF (ecrit_slab_oc.eq.1) then
          open(97,file='div_slab.dat',STATUS='UNKNOWN')
@@ -820,9 +795,10 @@
        CALL iophys_ini
 #endif
+
+!=====================================================================
 ! START OF THE TEMPORAL LOOP :
 !=====================================================================
            
       it_end = nint(fnday*day_step)
-!test JLD     it_end = 10
       do while(it.le.it_end)
 
@@ -832,5 +808,4 @@
          print*,'PAS DE TEMPS ',timestep
        endif
-!Al1 demande de restartphy.nc
        if (it.eq.it_end) lastcall=.True.
 
@@ -844,10 +819,13 @@
 !  Geopotential :
 !---------------------------------------------------------------------
-
+!        phis(1)=zsurf*RG
+!        phi(1)=phis(1)+RD*temp(1)*(plev(1)-play(1))/(.5*(plev(1)+play(1)))
         phi(1)=RD*temp(1)*(plev(1)-play(1))/(.5*(plev(1)+play(1)))
+
         do l = 1, llm-1
           phi(l+1)=phi(l)+RD*(temp(l)+temp(l+1))*                           &
      &    (play(l)-play(l+1))/(play(l)+play(l+1))
         enddo
+
 
 !---------------------------------------------------------------------
@@ -950,6 +928,5 @@
        sfdt = sin(0.5*fcoriolis*timestep)
        cfdt = cos(0.5*fcoriolis*timestep)
-!       print *,'fcoriolis,sfdt,cfdt,timestep',fcoriolis,sfdt,cfdt,timestep
-!
+
         d_u_age(1:mxcalc)= -2.*sfdt/timestep*                                &
      &          (sfdt*(u(1:mxcalc)-ug(1:mxcalc)) -                          &
@@ -1030,17 +1007,15 @@
         temp(1:mxcalc)=temp(1:mxcalc)+timestep*(                            &
      &              dt_phys(1:mxcalc)                                       &
-     &             +d_t_adv(1:mxcalc)                                      &
-     &             +d_t_nudge(1:mxcalc)                                      &
+     &             +d_t_adv(1:mxcalc)                                       &
+     &             +d_t_nudge(1:mxcalc)                                     &
      &             +dt_cooling(1:mxcalc))  ! Taux de chauffage ou refroid.
 
 
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!=======================================================================
 !! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!     endif  ! forcing_sandu or forcing_astex
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!=======================================================================
 
         teta=temp*(pzero/play)**rkappa
-!
+
 !---------------------------------------------------------------------
 !   Nudge soil temperature if requested
@@ -1080,7 +1055,5 @@
 
 !  incremente day time
-!        print*,'daytime bef',daytime,1./day_step
         daytime = daytime+1./day_step
-!Al1dbg
         day = int(daytime+0.1/day_step)
 !        time = max(daytime-day,0.0)
@@ -1088,14 +1061,12 @@
 !cc        time = real(mod(it,day_step))/day_step
         time = time_ini/24.+real(mod(it,day_step))/day_step
-!        print*,'daytime nxt time',daytime,time
         it=it+1
 
       enddo
 
-!Al1
       if (ecrit_slab_oc.ne.-1) close(97)
 
 !Al1 Call to 1D equivalent of dynredem (an,mois,jour,heure ?)
-! -------------------------------------
+! ---------------------------------------------------------------------------
        call dyn1dredem("restart1dyn.nc",                                    &
      &              plev,play,phi,phis,presnivs,                            &
Index: LMDZ6/trunk/libf/phylmd/flux_arp.h
===================================================================
--- LMDZ6/trunk/libf/phylmd/flux_arp.h	(revision 3779)
+++ LMDZ6/trunk/libf/phylmd/flux_arp.h	(revision 3780)
@@ -1,17 +1,23 @@
 !
 ! $Id: flux_arp.h 2010-08-04 17:02:56Z lahellec $
+! Modif EV, 10/2020
 !
       logical :: ok_flux_surf
       logical :: ok_prescr_ust !for prescribed ustar
+      logical :: ok_prescr_beta
+      logical :: ok_forc_tsurf
+
+
       real :: fsens
       real :: flat
+      real :: betaevap
       real :: ust
       real :: tg
 
-      common /flux_arp/fsens,flat,ust,tg,ok_flux_surf,ok_prescr_ust
+      common /flux_arp/fsens,flat,ust,tg,ok_flux_surf,ok_prescr_ust,ok_prescr_beta,betaevap,ok_forc_tsurf
 
 !$OMP THREADPRIVATE(/flux_arp/)
 
 
+      
 
-
Index: LMDZ6/trunk/libf/phylmd/ocean_forced_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/ocean_forced_mod.F90	(revision 3779)
+++ LMDZ6/trunk/libf/phylmd/ocean_forced_mod.F90	(revision 3780)
@@ -38,5 +38,5 @@
     INCLUDE "YOMCST.h"
     INCLUDE "clesphys.h"
-
+    INCLUDE "flux_arp.h"
 
 ! Input arguments
@@ -76,6 +76,6 @@
     REAL, DIMENSION(klon)       :: u1_lay, v1_lay
     LOGICAL                     :: check=.FALSE.
-    REAL, DIMENSION(klon) :: sens_prec_liq, sens_prec_sol    
-    REAL, DIMENSION(klon) :: lat_prec_liq, lat_prec_sol    
+    REAL, DIMENSION(klon)       :: sens_prec_liq, sens_prec_sol    
+    REAL, DIMENSION(klon)       :: lat_prec_liq, lat_prec_sol    
 
 !****************************************************************************************
@@ -92,5 +92,10 @@
 !!jyg    if (knon.eq.1) then ! single-column model
     if (klon_glo.eq.1) then ! single-column model
-      CALL read_tsurf1d(knon,tsurf_lim) ! new
+      ! EV: now surface Tin flux_arp.h
+      !CALL read_tsurf1d(knon,tsurf_lim) ! new
+       DO i = 1, knon
+        tsurf_lim(i) = tg
+       ENDDO
+
     else ! GCM
       CALL limit_read_sst(knon,knindex,tsurf_lim)
@@ -104,7 +109,12 @@
 !****************************************************************************************
 ! Set some variables for calcul_fluxs
-    cal = 0.
-    beta = 1.
-    dif_grnd = 0.
+    !cal = 0.
+    !beta = 1.
+    !dif_grnd = 0.
+    ! EV: use calbeta to calculate beta
+    
+    CALL calbeta(dtime, is_oce, knon, snow, beta*0., beta, cal, dif_grnd)
+
+
     alb_neig(:) = 0.
     agesno(:) = 0.
@@ -172,4 +182,5 @@
     INCLUDE "YOMCST.h"
     INCLUDE "clesphys.h"
+    INCLUDE "flux_arp.h"
 
 ! Input arguments
@@ -233,6 +244,6 @@
     tsurf_tmp(:) = tsurf_in(:)
 
-! calculate the parameters cal, beta, capsol and dif_grnd
-    CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, capsol, dif_grnd)
+! calculate the parameters cal, beta, capsol and dif_grnd and then recalculate cal
+    CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, cal, dif_grnd)
 
     
@@ -249,5 +260,5 @@
     ENDIF
 
-    beta = 1.0
+    !beta = 1.0
     sens_prec_liq = 0.; sens_prec_sol = 0.; lat_prec_liq = 0.; lat_prec_sol = 0.
 
@@ -304,29 +315,6 @@
   END SUBROUTINE ocean_forced_ice
 
+!
 !************************************************************************
-! 1D case
-!************************************************************************
-  SUBROUTINE read_tsurf1d(knon,sst_out)
-
-! This subroutine specifies the surface temperature to be used in 1D simulations
-
-      USE dimphy, ONLY : klon
-
-      INTEGER, INTENT(IN)                  :: knon     ! nomber of points on compressed grid
-      REAL, DIMENSION(klon), INTENT(OUT)   :: sst_out  ! tsurf used to force the single-column model
-
-       INTEGER :: i
-! COMMON defined in lmdz1d.F:
-       real ts_cur
-       common /sst_forcing/ts_cur
-
-       DO i = 1, knon
-        sst_out(i) = ts_cur
-       ENDDO
-
-      END SUBROUTINE read_tsurf1d
-
-!
-!************************************************************************
 !
 END MODULE ocean_forced_mod
Index: LMDZ6/trunk/libf/phylmd/ocean_slab_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/ocean_slab_mod.F90	(revision 3779)
+++ LMDZ6/trunk/libf/phylmd/ocean_slab_mod.F90	(revision 3780)
@@ -421,7 +421,11 @@
 !    
 !****************************************************************************************
-    cal(:)      = 0. ! infinite thermal inertia
-    beta(:)     = 1. ! wet surface
-    dif_grnd(:) = 0. ! no diffusion into ground
+    !cal(:)      = 0. ! infinite thermal inertia
+    !beta(:)     = 1. ! wet surface
+    !dif_grnd(:) = 0. ! no diffusion into ground
+    ! EV: use calbeta
+    CALL calbeta(dtime, is_oce, knon, snow,qsurf, beta, cal, dif_grnd)
+
+
     
 ! Suppose zero surface speed
@@ -742,6 +746,10 @@
 ! set beta, cal, compute conduction fluxes inside ice/snow 
     slab_bilg(:)=0.
-    dif_grnd(:)=0.
-    beta(:) = 1.
+    !dif_grnd(:)=0.
+    !beta(:) = 1.
+    ! EV: use calbeta to calculate beta and then recalculate properly cal
+    CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, cal, dif_grnd)
+
+
     DO i=1,knon
     ki=knindex(i)
Index: LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 3779)
+++ LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 3780)
@@ -206,5 +206,5 @@
 !jyg<
 !!       zxfluxt,   zxfluxq,   q2m,      flux_q, tke,   &
-       zxfluxt,   zxfluxq,   q2m,      flux_q, tke_x,   &
+       zxfluxt,   zxfluxq,   q2m,      flux_q, tke_x,  &
 !>jyg
 !!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
@@ -504,6 +504,7 @@
     CHARACTER(len=8), DIMENSION(nbsrf), SAVE :: cl_surf
 !$OMP THREADPRIVATE(cl_surf)
-    REAL, SAVE                               :: beta_land         ! beta for wx_dts
-!$OMP THREADPRIVATE(beta_land)
+ ! EV Ne sert plus:
+ !   REAL, SAVE                               :: beta_land         ! beta for wx_dts
+!!$OMP THREADPRIVATE(beta_land)
 
 ! Other local variables
@@ -845,8 +846,10 @@
        ! Initialize ok_flux_surf (for 1D model)
        if (klon_glo>1) ok_flux_surf=.FALSE.
+       if (klon_glo>1) ok_forc_tsurf=.FALSE.
+
 
        ! intialize beta_land
-       beta_land = 0.5
-       call getin_p('beta_land', beta_land)
+       !beta_land = 0.5
+       !call getin_p('beta_land', beta_land)
        
        ! Initilize debug IO
@@ -947,4 +950,5 @@
 !!    tke(:,:,is_ave)=0.
     tke_x(:,:,is_ave)=0.
+
     wake_dltke(:,:,is_ave)=0.
 !>jyg
@@ -978,4 +982,5 @@
 !!    d_t_diss= 0.0 ;d_u = 0.0     ; d_v = 0.0 
     yqsol = 0.0    
+
     ytke=0.
 !FC 
@@ -1388,4 +1393,5 @@
              ytke_w(j,k)      = tke_x(i,k,nsrf)+wake_dltke(i,k,nsrf)
              ywake_dltke(j,k) = wake_dltke(i,k,nsrf)
+            
 !>jyg
           ENDDO
@@ -1462,4 +1468,5 @@
       ENDDO
      ENDIF
+
         IF (prt_level >=10) print *,'clcdrag -> ycdragh ', ycdragh
        ELSE  !(iflag_split .eq.0)
@@ -1545,4 +1552,5 @@
       print *,' args coef_diff_turb: ycdragh ', ycdragh
       print *,' args coef_diff_turb: ytke ', ytke
+
        ENDIF
         CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
@@ -1574,4 +1582,5 @@
       print *,' args coef_diff_turb: ycdragh_x ', ycdragh_x
       print *,' args coef_diff_turb: ytke_x ', ytke_x
+
        ENDIF
         CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
@@ -2020,7 +2029,15 @@
 !
 !****************************************************************************************
-
-!!!
-!!! jyg le 10/04/2013
+!!
+!!!
+!!! jyg le 10/04/2013 et EV 10/2020
+
+        IF (ok_forc_tsurf) THEN
+            DO j=1,knon
+                ytsurf_new(j)=tg
+                y_d_ts(j) = ytsurf_new(j) - yts(j) 
+            ENDDO
+        ENDIF ! ok_forc_tsurf
+
 !!!
         IF (ok_flux_surf) THEN
@@ -2451,4 +2468,5 @@
               tke_x(i,k,nsrf)    = ytke(j,k)
               tke_x(i,k,is_ave) = tke_x(i,k,is_ave) + ytke(j,k)*ypct(j)
+
 !>jyg
            ENDDO
@@ -2464,6 +2482,7 @@
 !!            tke(i,k,is_ave) = tke(i,k,is_ave) + tke(i,k,nsrf)*ypct(j)
             tke_x(i,k,nsrf)   = ytke_x(j,k)
-            tke_x(i,k,is_ave)   = tke_x(i,k,is_ave) + tke_x(i,k,nsrf)*ypct(j)
+            tke_x(i,k,is_ave)   = tke_x(i,k,is_ave) + tke_x(i,k,nsrf)*ypct(j)        
             wake_dltke(i,k,is_ave)   = wake_dltke(i,k,is_ave) + wake_dltke(i,k,nsrf)*ypct(j)
+            
 
 !>jyg
Index: LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 3779)
+++ LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 3780)
@@ -16,7 +16,6 @@
       REAL, SAVE, ALLOCATABLE :: u_seri(:,:), v_seri(:,:)
       !$OMP THREADPRIVATE(u_seri, v_seri)
-      REAL, SAVE, ALLOCATABLE :: l_mixmin(:,:,:), l_mix(:,:,:)
-      !$OMP THREADPRIVATE(l_mixmin, l_mix)
-
+      REAL, SAVE, ALLOCATABLE :: l_mixmin(:,:,:), l_mix(:,:,:), tke_dissip(:,:,:)
+      !$OMP THREADPRIVATE(l_mixmin, l_mix, tke_dissip)
       REAL, SAVE, ALLOCATABLE :: tr_seri(:,:,:)
       !$OMP THREADPRIVATE(tr_seri)
@@ -445,6 +444,6 @@
       REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: ref_liq_pi, ref_ice_pi
 !$OMP THREADPRIVATE(ref_liq_pi, ref_ice_pi)
-      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: zx_rh
-!$OMP THREADPRIVATE(zx_rh)
+      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: zx_rh, zx_rhl, zx_rhi
+!$OMP THREADPRIVATE(zx_rh, zx_rhl, zx_rhi)
       REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: prfl, psfl, fraca
 !$OMP THREADPRIVATE(prfl, psfl, fraca)
@@ -561,6 +560,6 @@
       ALLOCATE(t_seri(klon,klev),q_seri(klon,klev),ql_seri(klon,klev),qs_seri(klon,klev))
       ALLOCATE(u_seri(klon,klev),v_seri(klon,klev))
-      ALLOCATE(l_mixmin(klon,klev+1,nbsrf), l_mix(klon,klev+1,nbsrf))
-      l_mix(:,:,:)=0. ; l_mixmin(:,:,:)=0. ! doit etre initialse car pas toujours remplis
+      ALLOCATE(l_mixmin(klon,klev+1,nbsrf), l_mix(klon,klev+1,nbsrf), tke_dissip(klon,klev+1,nbsrf))
+      l_mix(:,:,:)=0. ; l_mixmin(:,:,:)=0. ; tke_dissip(:,:,:)=0. ! doit etre initialse car pas toujours remplis
 
       ALLOCATE(tr_seri(klon,klev,nbtr))
@@ -780,5 +779,5 @@
       ALLOCATE(ref_liq(klon, klev), ref_ice(klon, klev), theta(klon, klev))
       ALLOCATE(ref_liq_pi(klon, klev), ref_ice_pi(klon, klev))
-      ALLOCATE(zphi(klon, klev), zx_rh(klon, klev))
+      ALLOCATE(zphi(klon, klev), zx_rh(klon, klev), zx_rhl(klon,klev), zx_rhi(klon,klev))
       ALLOCATE(pmfd(klon, klev), pmfu(klon, klev))
 
@@ -879,5 +878,5 @@
       DEALLOCATE(t_seri,q_seri,ql_seri,qs_seri)
       DEALLOCATE(u_seri,v_seri)
-      DEALLOCATE(l_mixmin,l_mix)
+      DEALLOCATE(l_mixmin,l_mix, tke_dissip)
 
       DEALLOCATE(tr_seri)
@@ -1074,5 +1073,5 @@
       DEALLOCATE(ref_liq, ref_ice, theta)
       DEALLOCATE(ref_liq_pi, ref_ice_pi)
-      DEALLOCATE(zphi, zx_rh)
+      DEALLOCATE(zphi, zx_rh, zx_rhl, zx_rhi)
       DEALLOCATE(pmfd, pmfu)
 
Index: LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 3779)
+++ LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 3780)
@@ -1005,9 +1005,10 @@
   TYPE(ctrl_out), SAVE :: o_tke = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'tke ', 'TKE', 'm2/s2', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_tke_dissip = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'tke_dissip ', 'TKE DISSIPATION', 'm2/s3', (/ ('', i=1, 10) /))   
   TYPE(ctrl_out), SAVE :: o_tke_max = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'tke_max', 'TKE max', 'm2/s2',                                  &
       (/ 't_max(X)', 't_max(X)', 't_max(X)', 't_max(X)', 't_max(X)', &
          't_max(X)', 't_max(X)', 't_max(X)', 't_max(X)', 't_max(X)' /))
-
   TYPE(ctrl_out), SAVE, DIMENSION(4) :: o_tke_srf      = (/             &
       ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11, 11/),'tke_ter',       &
@@ -1050,4 +1051,5 @@
       ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/),'l_mix_sic',       &
       "min PBL mixing length "//clnsurf(4),"m", (/ ('', i=1, 10) /)) /)
+
 
   TYPE(ctrl_out), SAVE, DIMENSION(4) :: o_tke_max_srf  = (/                          &
@@ -1434,4 +1436,8 @@
   TYPE(ctrl_out), SAVE :: o_rhum = ctrl_out((/ 2, 5, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'rhum', 'Relative humidity', '-', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_rhl = ctrl_out((/ 2, 6, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'rhl', 'Relative humidity wrt liquid', '%', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_rhi = ctrl_out((/ 2, 6, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'rhi', 'Relative humidity wrt ice', '%', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_ozone = ctrl_out((/ 2, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'ozone', 'Ozone mole fraction', '-', (/ ('', i=1, 10) /))
Index: LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 3779)
+++ LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 3780)
@@ -133,8 +133,8 @@
          o_vitu, o_vitv, o_vitw, o_pres, o_paprs, &
          o_zfull, o_zhalf, o_rneb, o_rnebjn, o_rnebcon, &
-         o_rnebls, o_rneblsvol, o_rhum, o_ozone, o_ozone_light, &
+         o_rnebls, o_rneblsvol, o_rhum, o_rhl, o_rhi, o_ozone, o_ozone_light, &
          o_duphy, o_dtphy, o_dqphy, o_dqphy2d, o_dqlphy, o_dqlphy2d, &
          o_dqsphy, o_dqsphy2d, o_albe_srf, o_z0m_srf, o_z0h_srf, &
-         o_ages_srf, o_snow_srf, o_alb1, o_alb2, o_tke, &
+         o_ages_srf, o_snow_srf, o_alb1, o_alb2, o_tke, o_tke_dissip, &
          o_tke_max, o_kz, o_kz_max, o_clwcon, &
          o_dtdyn, o_dqdyn, o_dqdyn2d, o_dqldyn, o_dqldyn2d, & 
@@ -251,5 +251,5 @@
          zt2m_cor,zq2m_cor,zu10m_cor,zv10m_cor, zrh2m_cor, zqsat2m_cor, &
          t2m_min_mon, t2m_max_mon, evap, &
-         l_mixmin,l_mix, &
+         l_mixmin,l_mix, tke_dissip, &
          zu10m, zv10m, zq2m, zustar, zxqsurf, &
          rain_lsc, rain_num, snow_lsc, bils, sens, fder, &
@@ -297,5 +297,5 @@
          ql_seri, qs_seri, tr_seri, &
          zphi, u_seri, v_seri, omega, cldfra, &
-         rneb, rnebjn, rneblsvol, zx_rh, d_t_dyn,  & 
+         rneb, rnebjn, rneblsvol, zx_rh, zx_rhl, zx_rhi, d_t_dyn,  & 
          d_q_dyn,  d_ql_dyn, d_qs_dyn, &
          d_q_dyn2d,  d_ql_dyn2d, d_qs_dyn2d, &
@@ -1059,4 +1059,6 @@
              CALL histwrite_phy(o_l_mixmin(nsrf),  l_mixmin(:,1:klev,nsrf))
              CALL histwrite_phy(o_tke_max_srf(nsrf),  pbl_tke(:,1:klev,nsrf))
+
+                         
           ENDIF
 !jyg<
@@ -1069,6 +1071,20 @@
 !            ENDIF
 
-
        ENDDO
+       
+                
+        IF (iflag_pbl > 1) THEN
+          zx_tmp_fi3d=0.
+          IF (vars_defined) THEN
+             DO nsrf=1,nbsrf
+                DO k=1,klev
+                   zx_tmp_fi3d(:,k)=zx_tmp_fi3d(:,k) &
+                        +pctsrf(:,nsrf)*tke_dissip(:,k,nsrf)
+                ENDDO
+             ENDDO
+          ENDIF
+          
+          CALL histwrite_phy(o_tke_dissip, zx_tmp_fi3d)    
+       ENDIF
 
        IF (vars_defined) zx_tmp_fi2d(1 : klon) = sens_prec_liq_o(1 : klon, 1)
@@ -1695,4 +1711,7 @@
        CALL histwrite_phy(o_rnebjn, zx_tmp_fi3d)
        CALL histwrite_phy(o_rhum, zx_rh)
+       CALL histwrite_phy(o_rhl, zx_rhl*100.)
+       CALL histwrite_phy(o_rhi, zx_rhi*100.)
+
        
        IF (vars_defined) zx_tmp_fi3d = wo(:, :, 1) * dobson_u * 1e3 / zmasse / rmo3 * rmd
@@ -1753,6 +1772,6 @@
           ENDIF
           CALL histwrite_phy(o_tke, zx_tmp_fi3d)
-
-          CALL histwrite_phy(o_tke_max, zx_tmp_fi3d)
+          CALL histwrite_phy(o_tke_max, zx_tmp_fi3d)  
+
        ENDIF
 
Index: LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90	(revision 3779)
+++ LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90	(revision 3780)
@@ -449,5 +449,4 @@
 
 include "clesphys.h"
-
       IF (is_initialized) RETURN
       is_initialized=.TRUE.
Index: LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 3779)
+++ LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 3780)
@@ -273,5 +273,5 @@
        ref_liq, ref_ice, theta,  &
        ref_liq_pi, ref_ice_pi,  &
-       zphi, zx_rh,  &
+       zphi, zx_rh, zx_rhl, zx_rhi,  &
        pmfd, pmfu,  &
        !
@@ -3738,4 +3738,6 @@
           ENDIF
           zx_rh(i,k) = q_seri(i,k)/zx_qs
+          zx_rhl(i,k) = q_seri(i,k)/(qsatl(zx_t)/pplay(i,k))
+          zx_rhi(i,k) = q_seri(i,k)/(qsats(zx_t)/pplay(i,k))
           zqsat(i,k)=zx_qs
        ENDDO
Index: LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90	(revision 3779)
+++ LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90	(revision 3780)
@@ -241,4 +241,10 @@
 ! 
 !****************************************************************************************
+
+    ! EV: use calbeta
+    CALL calbeta(dtime, is_lic, knon, snow, qsol, beta, cal, dif_grnd)
+
+
+    ! use soil model and recalculate properly cal
     IF (soil_model) THEN 
        CALL soil(dtime, is_lic, knon, snow, tsurf, tsoil, soilcap, soilflux)
@@ -255,6 +261,6 @@
 !
 !****************************************************************************************
-    beta(:) = 1.0
-    dif_grnd(:) = 0.0
+ !   beta(:) = 1.0
+ !   dif_grnd(:) = 0.0
 
 ! Suppose zero surface speed
@@ -291,4 +297,5 @@
 !****************************************************************************************
     CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))  
+
     WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
     zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0)))
@@ -303,5 +310,5 @@
 !IM: KstaTER0.77 & LMD_ARMIP6    
 
-! Attantion: alb1 and alb2 are the same!
+! Attantion: alb1 and alb2 are not the same!
     alb1(1:knon)  = alb_vis_sno_lic
     alb2(1:knon)  = alb_nir_sno_lic
Index: LMDZ6/trunk/libf/phylmd/yamada4.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/yamada4.F90	(revision 3779)
+++ LMDZ6/trunk/libf/phylmd/yamada4.F90	(revision 3780)
@@ -6,5 +6,6 @@
   USE dimphy
   USE ioipsl_getin_p_mod, ONLY : getin_p
-
+  USE phys_local_var_mod, only: tke_dissip
+  
   IMPLICIT NONE
   include "iniprint.h"
@@ -56,8 +57,10 @@
   !             iflag_pbl=11 -> the model starts with NP from start files created by ce0l
   !                          -> the model can run with longer time-steps.
-  !             2016/11/30 (EV etienne.vignon@univ-grenoble-alpes.fr)
+  !             2016/11/30 (EV etienne.vignon@lmd.ipsl.fr)
   !               On met tke (=q2/2) en entr??e plut??t que q2
   !               On corrige l'update de la tke
-  !
+  !             2020/10/01 (EV)
+  !               On ajoute la dissipation de la TKE en diagnostique de sortie
+  !                  
   ! Inpout/Output :
   !==============
@@ -121,5 +124,5 @@
   REAL,SAVE :: viscom,viscoh
   !$OMP THREADPRIVATE( hboville,viscom,viscoh)
-  INTEGER ig, k
+  INTEGER ig, jg, k
   REAL ri, zrif, zalpha, zsm, zsn
   REAL rif(klon, klev+1), sm(klon, klev+1), alpha(klon, klev)
@@ -416,4 +419,9 @@
   ELSE IF (iflag_pbl>=10) THEN
 
+    shear(:,:)=0.
+    buoy(:,:)=0.
+    dissip(:,:)=0.
+    km(:,:)=0.
+        
     IF (yamada4_num>=1) THEN
  
@@ -424,8 +432,8 @@
       shear(ig,k)=km(ig, k)*m2(ig, k)
       buoy(ig,k)=km(ig, k)*m2(ig, k)*(-1.*rif(ig,k))
-      dissip(ig,k)=((sqrt(q2(ig,k)))**3)/(b1*l(ig,k))
+      dissip(ig,k)=min(max(((sqrt(q2(ig,k)))**3)/(b1*l(ig,k)),1.E-12),1.E4)
      ENDDO
     ENDDO
-
+    
     IF (yamada4_num==1) THEN ! Schema du MAR tel quel
        DO k = 2, klev - 1
@@ -478,4 +486,9 @@
         ENDDO
        ENDDO
+       
+      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+      !! Attention, yamada4_num=5 est inexacte car néglige les termes de flottabilité
+      !! en conditions instables
+      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     ELSE IF (yamada4_num==5) THEN ! version modifiee avec integration exacte pour la dissipation
        DO k = 2, klev - 1
@@ -507,13 +520,23 @@
        DO k = 2, klev - 1
          DO ig=1,ngrid
+         !tkeprov=q2(ig,k)/ydeux
+         !tkeprov=tkeprov+max(buoy(ig,k)+shear(ig,k),0.)*dt
+         !disseff=dissip(ig,k)-min(0.,buoy(ig,k)+shear(ig,k)) 
+         !tkeexp=exp(-dt*disseff/tkeprov)
+         !tkeprov= tkeprov*tkeexp
+         !q2(ig,k)=tkeprov*ydeux
+         ! En cas stable, on traite la flotabilite comme la
+         ! dissipation, en supposant que dissipeff/TKE est constant.
+         ! Puis on prend la solution exacte
+         ! 
+         ! With drag and dissipation from high vegetation (EV & FC, 05/10/2020)
+         winds(ig,k)=sqrt(u(ig,k)**2+v(ig,k)**2)
          tkeprov=q2(ig,k)/ydeux
-         tkeprov=tkeprov+max(buoy(ig,k)+shear(ig,k),0.)*dt
-         disseff=dissip(ig,k)-min(0.,buoy(ig,k)+shear(ig,k))
+         tkeprov=tkeprov+max(buoy(ig,k)+shear(ig,k)+drgpro(ig,k)*(winds(ig,k))**3,0.)*dt
+         disseff=dissip(ig,k)-min(0.,buoy(ig,k)+shear(ig,k)+drgpro(ig,k)*(winds(ig,k))**3) + drgpro(ig,k)*tkeprov
          tkeexp=exp(-dt*disseff/tkeprov)
          tkeprov= tkeprov*tkeexp
-         q2(ig,k)=tkeprov*ydeux
-         ! En cas stable, on traite la flotabilite comme la
-         ! dissipation, en supposant que buoy/q2^3 est constant.
-         ! Puis on prend la solution exacte
+         q2(ig,k)=tkeprov*ydeux        
+         
         ENDDO
        ENDDO
@@ -533,5 +556,5 @@
 !     q2(1:ngrid, k) = q2(1:ngrid, k) + dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k))
       q2(1:ngrid, k) = min(max(q2(1:ngrid,k),1.E-10), 1.E4)
-       q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(yun*l(1:ngrid,k)*b1))
+      q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(yun*l(1:ngrid,k)*b1))
 !     q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(2*l(1:ngrid,k)*b1))
       q2(1:ngrid, k) = q2(1:ngrid, k)*q2(1:ngrid, k)
@@ -725,4 +748,17 @@
 
 !============================================================================
+! Diagnostique de la dissipation
+!============================================================================
+
+! Diagnostics
+ tke_dissip(1:ngrid,:,nsrf)=0.
+ DO k=2,klev
+    DO ig=1,ngrid
+       jg=ni(ig)
+       tke_dissip(jg,k,nsrf)=dissip(ig,k)
+    ENDDO
+ ENDDO
+ 
+!=============================================================================
 
   RETURN
@@ -1017,5 +1053,5 @@
 !=====================================================================
 
-
+  l1(1:ngrid,:)=0.
   IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN
 
@@ -1069,5 +1105,5 @@
    l2(1:ngrid,:)=0.0
    l_mixmin(1:ngrid,:,nsrf)=0.
-   l_mix(1:ngrid,:,nsrf)=0.
+   l_mix(1:ngrid,:,nsrf)=1.E-5
 
    IF (nsrf .EQ. 1) THEN
@@ -1135,7 +1171,8 @@
 
 
- DO k=2,klev
+ DO k=1,klev+1
     DO ig=1,ngrid
        lmix(ig,k)=MAX(MAX(l1(ig,k), l2(ig,k)),lmixmin)
+       lmix(ig,k)=MAX(lmix(ig,k),1.E-5)
    ENDDO
  ENDDO
@@ -1143,5 +1180,5 @@
 ! Diagnostics
 
- DO k=2,klev
+ DO k=2,klev+1
     DO ig=1,ngrid
        jg=ni(ig)
