Index: LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/1DUTILS.h
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/1DUTILS.h	(revision 3605)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/1DUTILS.h	(revision 3798)
@@ -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
@@ -304,4 +322,6 @@
        rugos = 0.0001
        CALL getin('rugos',rugos)
+! FH/2020/04/08/confinement: Pour le nouveau format standard, la rugosite s'appelle z0
+       CALL getin('z0',rugos)
 
 !Config  Key  = rugosh
@@ -357,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
@@ -518,4 +547,5 @@
        CALL getin('forc_ustar',forc_ustar)
        IF (forc_ustar .EQ. 1) ok_prescr_ust=.true.
+
 
 !Config  Key  = nudging_u
@@ -1246,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/branches/Ocean_skin/libf/phylmd/dyn1d/1D_decl_cases.h
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/1D_decl_cases.h	(revision 3605)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/1D_decl_cases.h	(revision 3798)
@@ -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
@@ -286,5 +287,5 @@
         real temp_nudg_mod_cas(llm),qv_nudg_mod_cas(llm),v_nudg_mod_cas(llm),u_nudg_mod_cas(llm)
         real u_mod_cas(llm),v_mod_cas(llm)
-        real omega_mod_cas(llm)
+        real omega_mod_cas(llm),tke_mod_cas(llm+1)
         real ht_mod_cas(llm),vt_mod_cas(llm),dt_mod_cas(llm),dtrad_mod_cas(llm)
         real hth_mod_cas(llm),vth_mod_cas(llm),dth_mod_cas(llm)
Index: LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/1D_interp_cases.h
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/1D_interp_cases.h	(revision 3605)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/1D_interp_cases.h	(revision 3798)
@@ -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=',     &
@@ -13,8 +13,8 @@
      &       ,u_cas,v_cas,ug_cas,vg_cas                                                     &
      &       ,temp_nudg_cas,qv_nudg_cas,u_nudg_cas,v_nudg_cas                               &
-     &       ,vitw_cas,omega_cas,du_cas,hu_cas,vu_cas                                       &
+     &       ,vitw_cas,omega_cas,tke_cas,du_cas,hu_cas,vu_cas                                       &
      &       ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas                           &
      &       ,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas,lat_cas,sens_cas,ustar_cas       &
-     &       ,uw_cas,vw_cas,q1_cas,q2_cas,tke_cas                                           &
+     &       ,uw_cas,vw_cas,q1_cas,q2_cas,tkes_cas                                           &
 !
      &       ,ts_prof_cas,ps_prof_cas,plev_prof_cas,t_prof_cas,theta_prof_cas,thv_prof_cas  &
@@ -22,22 +22,22 @@
      &       ,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas                                 &
      &       ,temp_nudg_prof_cas,qv_nudg_prof_cas,u_nudg_prof_cas,v_nudg_prof_cas           &
-     &       ,vitw_prof_cas,omega_prof_cas                                                  &
+     &       ,vitw_prof_cas,omega_prof_cas,tke_prof_cas                                     &
      &       ,du_prof_cas,hu_prof_cas,vu_prof_cas                                           &
      &       ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas       &
      &       ,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas                            &
      &       ,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 
+     &       ,sens_prof_cas,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas,tkes_prof_cas)
+! EV tg instead of ts_cur
+             tg = ts_prof_cas 
 !            psurf=plev_prof_cas(1)
              psurf=ps_prof_cas
 
 ! vertical interpolation:
-      CALL interp2_case_vertical_std(play,nlev_cas,plev_prof_cas                                              &
+      CALL interp2_case_vertical_std(play,plev,nlev_cas,plev_prof_cas                                              &
      &         ,t_prof_cas,theta_prof_cas,thv_prof_cas,thl_prof_cas                                       &
      &         ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas                                 &
      &         ,ug_prof_cas,vg_prof_cas                                                                   &
      &         ,temp_nudg_prof_cas,qv_nudg_prof_cas,u_nudg_prof_cas,v_nudg_prof_cas                       &
-     &         ,vitw_prof_cas,omega_prof_cas                                      &
+     &         ,vitw_prof_cas,omega_prof_cas,tke_prof_cas                                                 &
      &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas                   &
      &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas    &
@@ -47,5 +47,5 @@
      &         ,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas                                                 &
      &         ,temp_nudg_mod_cas,qv_nudg_mod_cas,u_nudg_mod_cas,v_nudg_mod_cas                           &
-     &         ,w_mod_cas,omega_mod_cas                                                                   &
+     &         ,w_mod_cas,omega_mod_cas,tke_mod_cas                                                       &
      &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas                         &
      &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas           &
@@ -109,4 +109,6 @@
       do l = 1, llm
 ! Modif w_mod_cas -> omega_mod_cas (MM+MPL 20170309)
+	print*, l, llm
+	print*, play(l), temp(l)
        omega(l) = -w_mod_cas(l)*play(l)*rg/(rd*temp(l))
       enddo
Index: LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/1D_read_forc_cases.h
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/1D_read_forc_cases.h	(revision 3605)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/1D_read_forc_cases.h	(revision 3798)
@@ -27,8 +27,8 @@
      &       ,u_cas,v_cas,ug_cas,vg_cas                                                     &
      &       ,temp_nudg_cas,qv_nudg_cas,u_nudg_cas,v_nudg_cas                               &
-     &       ,vitw_cas,omega_cas,du_cas,hu_cas,vu_cas                                       &
+     &       ,vitw_cas,omega_cas,tke_cas,du_cas,hu_cas,vu_cas                                       &
      &       ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas                           &
      &       ,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas,lat_cas,sens_cas,ustar_cas       &
-     &       ,uw_cas,vw_cas,q1_cas,q2_cas,tke_cas                                           &
+     &       ,uw_cas,vw_cas,q1_cas,q2_cas,tkes_cas                                           &
 !
      &       ,ts_prof_cas,ps_prof_cas,plev_prof_cas,t_prof_cas,theta_prof_cas,thv_prof_cas  &
@@ -36,10 +36,10 @@
      &       ,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas                                 &
      &       ,temp_nudg_prof_cas,qv_nudg_prof_cas,u_nudg_prof_cas,v_nudg_prof_cas           &
-     &       ,vitw_prof_cas,omega_prof_cas                                                  &
+     &       ,vitw_prof_cas,omega_prof_cas,tke_prof_cas                                    &
      &       ,du_prof_cas,hu_prof_cas,vu_prof_cas                                           &
      &       ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas       &
      &       ,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas                            &
      &       ,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)
+     &       ,sens_prof_cas,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas,tkes_prof_cas)
 
       do l = 1, nlev_cas
@@ -49,5 +49,5 @@
 ! vertical interpolation using interpolation routine:
 !      write(*,*)'avant interp vert', t_prof
-      CALL interp2_case_vertical_std(play,nlev_cas,plev_prof_cas                                              &
+      CALL interp2_case_vertical_std(play,plev,nlev_cas,plev_prof_cas                                              &
      &         ,t_prof_cas,theta_prof_cas,thv_prof_cas,thl_prof_cas                                          &
      &         ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas                                 &
@@ -55,5 +55,5 @@
      &       ,temp_nudg_prof_cas,qv_nudg_prof_cas,u_nudg_prof_cas,v_nudg_prof_cas                         &
 
-     &         ,vitw_prof_cas,omega_prof_cas                                                              &
+     &         ,vitw_prof_cas,omega_prof_cas,tke_prof_cas                                                 &
      &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas                   &
      &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas    &
@@ -63,5 +63,5 @@
      &         ,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas                                                 &
      &         ,temp_nudg_mod_cas,qv_nudg_mod_cas,u_nudg_mod_cas,v_nudg_mod_cas                           &
-     &         ,w_mod_cas,omega_mod_cas                                                                   &
+     &         ,w_mod_cas,omega_mod_cas,tke_mod_cas                                                       &
      &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas                         &
      &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas           &
@@ -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)
@@ -95,8 +98,11 @@
        d_u_adv(l) = du_mod_cas(l)+hu_mod_cas(l)+vu_mod_cas(l)
        d_v_adv(l) = dv_mod_cas(l)+hv_mod_cas(l)+vv_mod_cas(l)
+      enddo
 
-!print*,'d_t_adv ',d_t_adv(1:20)*86400
+! Etienne pour initialisation de TKE
 
-      enddo     
+       do l=1,llm+1
+       pbl_tke(:,l,:)=tke_mod_cas(l)
+       enddo     
 
 ! Faut-il multiplier par -1 ? (MPL 20160713)
@@ -108,6 +114,6 @@
        IF (ok_prescr_ust) THEN
        ust=ustar_prof_cas
-       print *,'ust=',ust
        ENDIF
 
+
       endif !forcing_SCM
Index: LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/mod_1D_cases_read2.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/mod_1D_cases_read2.F90	(revision 3605)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/mod_1D_cases_read2.F90	(revision 3798)
@@ -316,5 +316,5 @@
 
 !**********************************************************************************************
-SUBROUTINE read_SCM_cas
+SUBROUTINE old_read_SCM_cas
       implicit none
 
@@ -457,5 +457,5 @@
 
         print*,'Allocations OK'
-        call read_SCM (nid,nlev_cas,nt_cas,                                                                     &
+        call old_read_SCM (nid,nlev_cas,nt_cas,                                                                     &
      &     ap_cas,bp_cas,z_cas,plev_cas,zh_cas,plevh_cas,t_cas,th_cas,thv_cas,thl_cas,qv_cas,                    &
      &     ql_cas,qi_cas,rh_cas,rv_cas,u_cas,v_cas,vitw_cas,omega_cas,ug_cas,vg_cas,du_cas,hu_cas,vu_cas,        &
@@ -470,5 +470,5 @@
 
 
-END SUBROUTINE read_SCM_cas
+END SUBROUTINE old_read_SCM_cas
 
 
@@ -846,5 +846,5 @@
 
 !======================================================================
-      subroutine read_SCM(nid,nlevel,ntime,                                       &
+      subroutine old_read_SCM(nid,nlevel,ntime,                                       &
      &     ap,bp,zz,pp,zzh,pph,temp,theta,thv,thl,qv,ql,qi,rh,rv,u,v,vitw,omega,ug,vg,&
      &     du,hu,vu,dv,hv,vv,dt,ht,vt,dq,hq,vq,                                    &
@@ -1073,5 +1073,5 @@
 
          return 
-         end subroutine read_SCM
+         end subroutine old_read_SCM
 !======================================================================
 
Index: LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/mod_1D_cases_read_std.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/mod_1D_cases_read_std.F90	(revision 3605)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/mod_1D_cases_read_std.F90	(revision 3798)
@@ -18,5 +18,5 @@
         real, allocatable::  t_cas(:,:),q_cas(:,:),qv_cas(:,:),ql_cas(:,:),qi_cas(:,:),rh_cas(:,:)
         real, allocatable::  th_cas(:,:),thv_cas(:,:),thl_cas(:,:),rv_cas(:,:)
-        real, allocatable::  u_cas(:,:),v_cas(:,:),vitw_cas(:,:),omega_cas(:,:)
+        real, allocatable::  u_cas(:,:),v_cas(:,:),vitw_cas(:,:),omega_cas(:,:),tke_cas(:,:)
 
 !forcing
@@ -30,5 +30,5 @@
         real, allocatable::  temp_nudg_cas(:,:),qv_nudg_cas(:,:),u_nudg_cas(:,:),v_nudg_cas(:,:)
         real, allocatable::  lat_cas(:),sens_cas(:),ts_cas(:),ps_cas(:),ustar_cas(:)
-        real, allocatable::  uw_cas(:,:),vw_cas(:,:),q1_cas(:,:),q2_cas(:,:),tke_cas(:)
+        real, allocatable::  uw_cas(:,:),vw_cas(:,:),q1_cas(:,:),q2_cas(:,:),tkes_cas(:)
 
 !champs interpoles
@@ -48,4 +48,5 @@
         real, allocatable::  vitw_prof_cas(:)
         real, allocatable::  omega_prof_cas(:)
+        real, allocatable::  tke_prof_cas(:)
         real, allocatable::  ug_prof_cas(:)
         real, allocatable::  vg_prof_cas(:)
@@ -73,5 +74,5 @@
 
 
-        real lat_prof_cas,sens_prof_cas,ts_prof_cas,ps_prof_cas,ustar_prof_cas,tke_prof_cas
+        real lat_prof_cas,sens_prof_cas,ts_prof_cas,ps_prof_cas,ustar_prof_cas,tkes_prof_cas
         real o3_cas,orog_cas,albedo_cas,emiss_cas,t_skin_cas,q_skin_cas,mom_rough,heat_rough,rugos_cas,sand_cas,clay_cas
       
@@ -92,5 +93,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
@@ -123,5 +124,5 @@
       ierr=NF_INQ_DIMLEN(nid,rid,nlev_cas)
       print*,'OK3 read2: nid,rid,nlev_cas',nid,rid,nlev_cas
-      IF ( .NOT. ( nlev_cas > 10 .AND. nlev_cas < 1000 )) THEN
+      IF ( .NOT. ( nlev_cas > 10 .AND. nlev_cas < 200000 )) THEN
               print*,'Valeur de nlev_cas peu probable'
               STOP
@@ -168,5 +169,5 @@
         allocate(th_cas(nlev_cas,nt_cas),thl_cas(nlev_cas,nt_cas),thv_cas(nlev_cas,nt_cas),rv_cas(nlev_cas,nt_cas))
         allocate(u_cas(nlev_cas,nt_cas),v_cas(nlev_cas,nt_cas),vitw_cas(nlev_cas,nt_cas),omega_cas(nlev_cas,nt_cas))
-
+        allocate(tke_cas(nlev_cas,nt_cas))
 !forcing
         allocate(ht_cas(nlev_cas,nt_cas),vt_cas(nlev_cas,nt_cas),dt_cas(nlev_cas,nt_cas),dtrad_cas(nlev_cas,nt_cas))
@@ -179,5 +180,5 @@
         allocate(temp_nudg_cas(nlev_cas,nt_cas),qv_nudg_cas(nlev_cas,nt_cas))
         allocate(u_nudg_cas(nlev_cas,nt_cas),v_nudg_cas(nlev_cas,nt_cas))
-        allocate(lat_cas(nt_cas),sens_cas(nt_cas),ts_cas(nt_cas),ps_cas(nt_cas),ustar_cas(nt_cas),tke_cas(nt_cas))
+        allocate(lat_cas(nt_cas),sens_cas(nt_cas),ts_cas(nt_cas),ps_cas(nt_cas),ustar_cas(nt_cas),tkes_cas(nt_cas))
         allocate(uw_cas(nlev_cas,nt_cas),vw_cas(nlev_cas,nt_cas),q1_cas(nlev_cas,nt_cas),q2_cas(nlev_cas,nt_cas))
 
@@ -200,4 +201,5 @@
         allocate(vitw_prof_cas(nlev_cas))
         allocate(omega_prof_cas(nlev_cas))
+        allocate(tke_prof_cas(nlev_cas))
         allocate(ug_prof_cas(nlev_cas))
         allocate(vg_prof_cas(nlev_cas))
@@ -228,9 +230,9 @@
         CALL read_SCM (nid,nlev_cas,nt_cas,                                                                     &
      &     ap_cas,bp_cas,z_cas,plev_cas,zh_cas,plevh_cas,t_cas,th_cas,thv_cas,thl_cas,qv_cas,                   &
-     &     ql_cas,qi_cas,rh_cas,rv_cas,u_cas,v_cas,vitw_cas,omega_cas,ug_cas,vg_cas,                            &
+     &     ql_cas,qi_cas,rh_cas,rv_cas,u_cas,v_cas,vitw_cas,omega_cas,tke_cas,ug_cas,vg_cas,                            &
      &     temp_nudg_cas,qv_nudg_cas,u_nudg_cas,v_nudg_cas,                                                     &
      &     du_cas,hu_cas,vu_cas,                                                                                &
      &     dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas,               &
-     &     dr_cas,hr_cas,vr_cas,dtrad_cas,sens_cas,lat_cas,ts_cas,ps_cas,ustar_cas,tke_cas,                      &
+     &     dr_cas,hr_cas,vr_cas,dtrad_cas,sens_cas,lat_cas,ts_cas,ps_cas,ustar_cas,tkes_cas,                      &
      &     uw_cas,vw_cas,q1_cas,q2_cas,orog_cas,albedo_cas,emiss_cas,t_skin_cas,q_skin_cas,mom_rough,heat_rough, &
      &     o3_cas,rugos_cas,clay_cas,sand_cas)
@@ -254,5 +256,5 @@
         deallocate(t_cas,q_cas,qv_cas,ql_cas,qi_cas,rh_cas)
         deallocate(th_cas,thl_cas,thv_cas,rv_cas)
-        deallocate(u_cas,v_cas,vitw_cas,omega_cas)
+        deallocate(u_cas,v_cas,vitw_cas,omega_cas,tke_cas)
         
 !forcing
@@ -265,5 +267,5 @@
         deallocate(ug_cas)
         deallocate(vg_cas)
-        deallocate(lat_cas,sens_cas,ts_cas,ps_cas,ustar_cas,tke_cas,uw_cas,vw_cas,q1_cas,q2_cas)
+        deallocate(lat_cas,sens_cas,ts_cas,ps_cas,ustar_cas,tkes_cas,uw_cas,vw_cas,q1_cas,q2_cas)
 
 !champs interpoles
@@ -283,4 +285,5 @@
         deallocate(vitw_prof_cas)
         deallocate(omega_prof_cas)
+        deallocate(tke_prof_cas)
         deallocate(ug_prof_cas)
         deallocate(vg_prof_cas)
@@ -312,8 +315,8 @@
 !=====================================================================
       SUBROUTINE read_SCM(nid,nlevel,ntime,                                       &
-     &     ap,bp,zz,pp,zzh,pph,temp,theta,thv,thl,qv,ql,qi,rh,rv,u,v,vitw,omega,ug,vg,&
+     &     ap,bp,zz,pp,zzh,pph,temp,theta,thv,thl,qv,ql,qi,rh,rv,u,v,vitw,omega,tke,ug,vg,&
      &     temp_nudg,qv_nudg,u_nudg,v_nudg,                                        &
      &     du,hu,vu,dv,hv,vv,dt,ht,vt,dq,hq,vq,                                    &
-     &     dth,hth,vth,dr,hr,vr,dtrad,sens,flat,ts,ps,ustar,tke,uw,vw,q1,q2,       &
+     &     dth,hth,vth,dr,hr,vr,dtrad,sens,flat,ts,ps,ustar,tkes,uw,vw,q1,q2,       &
      &     orog_cas,albedo_cas,emiss_cas,t_skin_cas,q_skin_cas,mom_rough,          &
      &     heat_rough,o3_cas,rugos_cas,clay_cas,sand_cas)
@@ -334,8 +337,8 @@
       real temp(nlevel,ntime),qv(nlevel,ntime),ql(nlevel,ntime),qi(nlevel,ntime),rh(nlevel,ntime)
       real theta(nlevel,ntime),thv(nlevel,ntime),thl(nlevel,ntime),rv(nlevel,ntime)
-      real u(nlevel,ntime),v(nlevel,ntime),tke(nlevel,ntime)
+      real u(nlevel,ntime),v(nlevel,ntime),tkes(ntime)
       real temp_nudg(nlevel,ntime),qv_nudg(nlevel,ntime),u_nudg(nlevel,ntime),v_nudg(nlevel,ntime)
       real ug(nlevel,ntime),vg(nlevel,ntime)
-      real vitw(nlevel,ntime),omega(nlevel,ntime)
+      real vitw(nlevel,ntime),omega(nlevel,ntime),tke(nlevel,ntime)
       real du(nlevel,ntime),hu(nlevel,ntime),vu(nlevel,ntime)
       real dv(nlevel,ntime),hv(nlevel,ntime),vv(nlevel,ntime)
@@ -365,22 +368,22 @@
      &'temp','qv','ql','qi','u','v','tke','pressure',& ! #5-#12
      ! coordonnees pression + temps #42
-     &'w','omega','ug','vg','uadv','uadvh','uadvv','vadv','vadvh','vadvv','tadv','tadvh','tadvv',& !  #13 - #25
-     &'qvadv','qvadvh','qvadvv','thadv','thadvh','thadvv','thladvh',                             & ! #26 - #33
-     & 'radv','radvh','radvv','radcool','q1','q2','ustress','vstress',                           & ! #34 - #41
-     & 'rh','temp_nudg','qv_nudg','u_nudg','v_nudg',&
-     &'height_f','pressure_forc','tempt','theta','thv','thl','qvt','qlt','qit','rv','ut','vt','tket',&
+     &'w','omega','ug','vg','uadv','uadvh','uadvv','vadv','vadvh','vadvv','temp_adv','tadvh','tadvv',& !  #13 - #25
+     &'qv_adv','qvadvh','qvadvv','thadv','thadvh','thadvv','thladvh',                             & ! #26 - #32
+     & 'radv','radvh','radvv','radcool','q1','q2','ustress','vstress',                           & ! #33 - #40
+     & 'rh','temp_nudging','qv_nudging','u_nudging','v_nudging',                                       & ! #41-45
+     &'height_f','pressure_forc','tempt','theta','thv','thl','qvt','qlt','qit','rv','ut','vt',   & ! #46-58
      ! coordonnees temps #12
-     &'sfc_sens_flx','sfc_lat_flx','ts','ps','ustar',&
+     &'tkes','sfc_sens_flx','sfc_lat_flx','ts','ps','ustar',&
      &'orog','albedo','emiss','t_skin','q_skin','mom_rough','heat_rough',&
      ! scalaires #4
      &'o3','rugos','clay','sand'/
 
-      do i=1,nbvar3d
-        missing_var(i)=0.
-      enddo
-
+!-----------------------------------------------------------------------
+! Checking availability of variable #i in the cas.nc file
+!     missing_var=1 if the variable is missing
 !-----------------------------------------------------------------------
 
        do i=1,nbvar3d
+         missing_var(i)=0.
          ierr=NF_INQ_VARID(nid,name_var(i),var3didin(i)) 
          if(ierr/=NF_NOERR) then
@@ -391,11 +394,12 @@
 
 !-----------------------------------------------------------------------
-! Activation de quelques cles en fonction des variables disponibles
+! Activating keys depending on the presence of specific variables in cas.nc
 !-----------------------------------------------------------------------
 if ( 1 == 1 ) THEN
-            if ( name_var(i) == 'temp_nudg' .and. nint(nudging_t)==0) stop 'Nudging inconsistency temp'
-            if ( name_var(i) == 'qv_nudg' .and. nint(nudging_qv)==0) stop 'Nudging inconsistency qv'
-            if ( name_var(i) == 'u_nudg' .and. nint(nudging_u)==0) stop 'Nudging inconsistency u'
-            if ( name_var(i) == 'v_nudg' .and. nint(nudging_u)==0) stop 'Nudging inconsistency v'
+! A MODIFIER: il faudrait dire nudging_temp mais faut le declarer dans compar1d.h etc...        
+!           if ( name_var(i) == 'temp_nudging' .and. nint(nudging_t)==0) stop 'Nudging inconsistency temp'
+            if ( name_var(i) == 'qv_nudging' .and. nint(nudging_qv)==0) stop 'Nudging inconsistency qv'
+            if ( name_var(i) == 'u_nudging' .and. nint(nudging_u)==0) stop 'Nudging inconsistency u'
+            if ( name_var(i) == 'v_nudging' .and. nint(nudging_u)==0) stop 'Nudging inconsistency v'
     ELSE
              print*,'GUIDAGE : CONSISTENCY CHECK DEACTIVATED FOR TESTS of SANDU/REF'
@@ -403,5 +407,7 @@
 
 !-----------------------------------------------------------------------
-           if(i.LE.4) then     ! Lecture des coord pression en (nlevelp1,lat,lon)
+! Reading variables 1D (N+1) vertical variables (nlevelp1,lat,lon)
+!-----------------------------------------------------------------------
+           if(i.LE.4) then
 #ifdef NC_DOUBLE
            ierr = NF_GET_VAR_DOUBLE(nid,var3didin(i),apbp)
@@ -414,6 +420,9 @@
               stop "getvarup"
            endif
-!-----------------------------------------------------------------------
-           else if(i.gt.4.and.i.LE.12) then   ! Lecture des variables en (time,nlevel,lat,lon)
+
+!-----------------------------------------------------------------------
+!  Reading 1D (N) vertical varialbes    (nlevel,lat,lon)   
+!-----------------------------------------------------------------------
+           else if(i.gt.4.and.i.LE.12) then  
 #ifdef NC_DOUBLE
            ierr = NF_GET_VAR_DOUBLE(nid,var3didin(i),resul1)
@@ -427,6 +436,10 @@
            endif
          print*,'Lecture de la variable #i ',i,name_var(i),minval(resul1),maxval(resul1)
-!-----------------------------------------------------------------------
-           else if(i.gt.12.and.i.LE.54) then   ! Lecture des variables en (time,nlevel,lat,lon)
+
+!-----------------------------------------------------------------------
+!  Reading 2D tim-vertical variables  (time,nlevel,lat,lon)
+!  TBD : seems to be the same as above.
+!-----------------------------------------------------------------------
+           else if(i.gt.12.and.i.LE.57) then
 #ifdef NC_DOUBLE
            ierr = NF_GET_VAR_DOUBLE(nid,var3didin(i),resul)
@@ -439,8 +452,10 @@
               stop "getvarup"
            endif
-
          print*,'Lecture de la variable #i ',i,name_var(i),minval(resul),maxval(resul)
-!-----------------------------------------------------------------------
-           else if (i.gt.54.and.i.LE.65) then   ! Lecture des variables en (time,lat,lon)
+
+!-----------------------------------------------------------------------
+!  Reading 1D time variables (time,lat,lon)
+!-----------------------------------------------------------------------
+           else if (i.gt.57.and.i.LE.63) then
 #ifdef NC_DOUBLE
            ierr = NF_GET_VAR_DOUBLE(nid,var3didin(i),resul2)
@@ -454,6 +469,9 @@
            endif
          print*,'Lecture de la variable #i  ',i,name_var(i),minval(resul2),maxval(resul2)
-!-----------------------------------------------------------------------
-           else     ! Lecture des constantes (lat,lon)
+
+!-----------------------------------------------------------------------
+! Reading scalar variables (lat,lon)
+!-----------------------------------------------------------------------
+           else
 #ifdef NC_DOUBLE
            ierr = NF_GET_VAR_DOUBLE(nid,var3didin(i),resul3)
@@ -469,4 +487,7 @@
            endif
          endif
+
+!-----------------------------------------------------------------------
+! Attributing variables
 !-----------------------------------------------------------------------
          select case(i)
@@ -528,6 +549,6 @@
            case(56) ; u=resul
            case(57) ; v=resul
-           case(58) ; tke=resul
-           case(59) ; sens=resul2   ! donnees indexees en time
+           case(58) ; tkes=resul2   ! donnees indexees en time
+           case(59) ; sens=resul2
            case(60) ; flat=resul2
            case(61) ; ts=resul2
@@ -577,4 +598,6 @@
 
 !**********************************************************************************************
+
+!**********************************************************************************************
         SUBROUTINE interp_case_time_std(day,day1,annee_ref                           &
 !    &         ,year_cas,day_cas,nt_cas,pdt_forc,nlev_cas                         &
@@ -583,9 +606,9 @@
      &         ,qv_cas,ql_cas,qi_cas,u_cas,v_cas                                  &
      &         ,ug_cas,vg_cas,temp_nudg_cas,qv_nudg_cas,u_nudg_cas,v_nudg_cas     &
-     &         ,vitw_cas,omega_cas,du_cas,hu_cas,vu_cas             &
+     &         ,vitw_cas,omega_cas,tke_cas,du_cas,hu_cas,vu_cas             &
      &         ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas               &
      &         ,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas                      &
      &         ,lat_cas,sens_cas,ustar_cas                                        &
-     &         ,uw_cas,vw_cas,q1_cas,q2_cas,tke_cas                               &
+     &         ,uw_cas,vw_cas,q1_cas,q2_cas,tkes_cas                               &
 !
      &         ,ts_prof_cas,ps_prof_cas,plev_prof_cas,t_prof_cas,theta_prof_cas               &
@@ -593,13 +616,17 @@
      &         ,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas                     &
      &         ,temp_nudg_prof_cas,qv_nudg_prof_cas,u_nudg_prof_cas,v_nudg_prof_cas     &
-     &         ,vitw_prof_cas,omega_prof_cas,du_prof_cas,hu_prof_cas,vu_prof_cas  &
+     &         ,vitw_prof_cas,omega_prof_cas,tke_prof_cas,du_prof_cas,hu_prof_cas,vu_prof_cas  &
      &         ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas                   &
      &         ,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas                &
      &         ,hq_prof_cas,vq_prof_cas,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)
-          
-
-        implicit none
+     &         ,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas,tkes_prof_cas)
+
+
+
+
+
+
+       implicit none
 
 !---------------------------------------------------------------------------------------
@@ -621,5 +648,6 @@
         real ts_cas(nt_cas),ps_cas(nt_cas)
         real plev_cas(nlev_cas,nt_cas)
-        real t_cas(nlev_cas,nt_cas),theta_cas(nlev_cas,nt_cas),thv_cas(nlev_cas,nt_cas),thl_cas(nlev_cas,nt_cas)
+        real t_cas(nlev_cas,nt_cas),theta_cas(nlev_cas,nt_cas)
+        real thv_cas(nlev_cas,nt_cas), thl_cas(nlev_cas,nt_cas)
         real qv_cas(nlev_cas,nt_cas),ql_cas(nlev_cas,nt_cas),qi_cas(nlev_cas,nt_cas)
         real u_cas(nlev_cas,nt_cas),v_cas(nlev_cas,nt_cas)
@@ -628,5 +656,5 @@
         real u_nudg_cas(nlev_cas,nt_cas),v_nudg_cas(nlev_cas,nt_cas)
 
-        real vitw_cas(nlev_cas,nt_cas),omega_cas(nlev_cas,nt_cas)
+        real vitw_cas(nlev_cas,nt_cas),omega_cas(nlev_cas,nt_cas),tke_cas(nlev_cas,nt_cas)
         real du_cas(nlev_cas,nt_cas),hu_cas(nlev_cas,nt_cas),vu_cas(nlev_cas,nt_cas)
         real dv_cas(nlev_cas,nt_cas),hv_cas(nlev_cas,nt_cas),vv_cas(nlev_cas,nt_cas)
@@ -635,5 +663,5 @@
         real dtrad_cas(nlev_cas,nt_cas)
         real dq_cas(nlev_cas,nt_cas),hq_cas(nlev_cas,nt_cas),vq_cas(nlev_cas,nt_cas)
-        real lat_cas(nt_cas),sens_cas(nt_cas),tke_cas(nt_cas)
+        real lat_cas(nt_cas),sens_cas(nt_cas),tkes_cas(nt_cas)
         real ustar_cas(nt_cas),uw_cas(nlev_cas,nt_cas),vw_cas(nlev_cas,nt_cas)
         real q1_cas(nlev_cas,nt_cas),q2_cas(nlev_cas,nt_cas)
@@ -648,5 +676,5 @@
         real u_nudg_prof_cas(nlev_cas),v_nudg_prof_cas(nlev_cas)
 
-        real vitw_prof_cas(nlev_cas),omega_prof_cas(nlev_cas)
+        real vitw_prof_cas(nlev_cas),omega_prof_cas(nlev_cas),tke_prof_cas(nlev_cas)
         real du_prof_cas(nlev_cas),hu_prof_cas(nlev_cas),vu_prof_cas(nlev_cas)
         real dv_prof_cas(nlev_cas),hv_prof_cas(nlev_cas),vv_prof_cas(nlev_cas)
@@ -655,5 +683,5 @@
         real dtrad_prof_cas(nlev_cas)
         real dq_prof_cas(nlev_cas),hq_prof_cas(nlev_cas),vq_prof_cas(nlev_cas)
-        real lat_prof_cas,sens_prof_cas,tke_prof_cas,ts_prof_cas,ps_prof_cas,ustar_prof_cas
+        real lat_prof_cas,sens_prof_cas,tkes_prof_cas,ts_prof_cas,ps_prof_cas,ustar_prof_cas
         real uw_prof_cas(nlev_cas),vw_prof_cas(nlev_cas),q1_prof_cas(nlev_cas),q2_prof_cas(nlev_cas)
 ! local:
@@ -739,6 +767,6 @@
        sens_prof_cas = sens_cas(it_cas2)                                 &
      &          -frac*(sens_cas(it_cas2)-sens_cas(it_cas1))
-       tke_prof_cas = tke_cas(it_cas2)                                   &
-     &          -frac*(tke_cas(it_cas2)-tke_cas(it_cas1))
+       tkes_prof_cas = tkes_cas(it_cas2)                                   &
+     &          -frac*(tkes_cas(it_cas2)-tkes_cas(it_cas1))
        ts_prof_cas = ts_cas(it_cas2)                                     &
      &          -frac*(ts_cas(it_cas2)-ts_cas(it_cas1))
@@ -786,4 +814,6 @@
         omega_prof_cas(k) = omega_cas(k,it_cas2)                         &
      &          -frac*(omega_cas(k,it_cas2)-omega_cas(k,it_cas1))
+        tke_prof_cas(k) = tke_cas(k,it_cas2)                         &
+     &          -frac*(tke_cas(k,it_cas2)-tke_cas(k,it_cas1))
         du_prof_cas(k) = du_cas(k,it_cas2)                               &
      &          -frac*(du_cas(k,it_cas2)-du_cas(k,it_cas1))
@@ -833,10 +863,10 @@
 !**********************************************************************************************
 !=====================================================================
-       SUBROUTINE interp2_case_vertical_std(play,nlev_cas,plev_prof_cas                                &
+       SUBROUTINE interp2_case_vertical_std(play,plev,nlev_cas,plev_prof_cas                           &
      &         ,t_prof_cas,th_prof_cas,thv_prof_cas,thl_prof_cas                                       &
      &         ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas                              &
      &         ,ug_prof_cas,vg_prof_cas                                                                &
      &         ,temp_nudg_prof_cas,qv_nudg_prof_cas,u_nudg_prof_cas,v_nudg_prof_cas                    &
-     &         ,vitw_prof_cas,omega_prof_cas                                                           &
+     &         ,vitw_prof_cas,omega_prof_cas,tke_prof_cas                                              &
      &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas                &
      &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas &
@@ -847,5 +877,5 @@
      &         ,ug_mod_cas,vg_mod_cas                                                                  &
      &         ,temp_nudg_mod_cas,qv_nudg_mod_cas,u_nudg_mod_cas,v_nudg_mod_cas                        &
-     &         ,w_mod_cas,omega_mod_cas                                                                &
+     &         ,w_mod_cas,omega_mod_cas,tke_mod_cas                                                    &
      &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas                      &
      &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas        &
@@ -870,9 +900,9 @@
 !       real hq_prof(nlevmax),vq_prof(nlevmax)
  
-       real play(llm), plev_prof_cas(nlev_cas) 
+       real play(llm), plev(llm+1), plev_prof_cas(nlev_cas) 
        real t_prof_cas(nlev_cas),th_prof_cas(nlev_cas),thv_prof_cas(nlev_cas),thl_prof_cas(nlev_cas)
        real qv_prof_cas(nlev_cas),ql_prof_cas(nlev_cas),qi_prof_cas(nlev_cas)
        real u_prof_cas(nlev_cas),v_prof_cas(nlev_cas)
-       real ug_prof_cas(nlev_cas),vg_prof_cas(nlev_cas), vitw_prof_cas(nlev_cas),omega_prof_cas(nlev_cas)
+       real ug_prof_cas(nlev_cas),vg_prof_cas(nlev_cas), vitw_prof_cas(nlev_cas),omega_prof_cas(nlev_cas),tke_prof_cas(nlev_cas)
        real temp_nudg_prof_cas(nlev_cas),qv_nudg_prof_cas(nlev_cas)
        real u_nudg_prof_cas(nlev_cas),v_nudg_prof_cas(nlev_cas)
@@ -887,5 +917,5 @@
        real qv_mod_cas(llm),ql_mod_cas(llm),qi_mod_cas(llm)
        real u_mod_cas(llm),v_mod_cas(llm)
-       real ug_mod_cas(llm),vg_mod_cas(llm), w_mod_cas(llm),omega_mod_cas(llm)
+       real ug_mod_cas(llm),vg_mod_cas(llm), w_mod_cas(llm),omega_mod_cas(llm),tke_mod_cas(llm+1)
        real temp_nudg_mod_cas(llm),qv_nudg_mod_cas(llm)
        real u_nudg_mod_cas(llm),v_nudg_mod_cas(llm)
@@ -899,10 +929,7 @@
        real frac,frac1,frac2,fact
  
-!       do l = 1, llm
-!       print *,'debut interp2, play=',l,play(l)
-!       enddo
-!      do l = 1, nlev_cas
-!      print *,'debut interp2, plev_prof_cas=',l,play(l),plev_prof_cas(l)
-!      enddo
+
+
+! for variables defined at the middle of layers
 
        do l = 1, llm
@@ -932,5 +959,8 @@
          endif
 
+
+
          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))
@@ -1029,8 +1059,8 @@
          ug_mod_cas(l)= ug_prof_cas(nlev_cas)                          !jyg
          vg_mod_cas(l)= vg_prof_cas(nlev_cas)                          !jyg
-         temp_nudg_mod_cas(l)= temp_nudg_prof_cas(nlev_cas)                          !jyg
-         qv_nudg_mod_cas(l)= qv_nudg_prof_cas(nlev_cas)                          !jyg
-         u_nudg_mod_cas(l)= u_nudg_prof_cas(nlev_cas)                          !jyg
-         v_nudg_mod_cas(l)= v_nudg_prof_cas(nlev_cas)                          !jyg
+         temp_nudg_mod_cas(l)= temp_nudg_prof_cas(nlev_cas)            !jyg
+         qv_nudg_mod_cas(l)= qv_nudg_prof_cas(nlev_cas)                !jyg
+         u_nudg_mod_cas(l)= u_nudg_prof_cas(nlev_cas)                  !jyg
+         v_nudg_mod_cas(l)= v_nudg_prof_cas(nlev_cas)                  !jyg
          thv_mod_cas(l)= thv_prof_cas(nlev_cas)                        !jyg
          w_mod_cas(l)= 0.0                                             !jyg
@@ -1057,6 +1087,53 @@
        enddo ! l
 
+! for variables defined at layer interfaces (EV):
+
+
+       do l = 1, llm+1
+
+        if (plev(l).ge.plev_prof_cas(nlev_cas)) then
+
+         mxcalc=l
+         k1=0
+         k2=0
+
+         if (plev(l).le.plev_prof_cas(1)) then
+
+         do k = 1, nlev_cas-1
+          if (plev(l).le.plev_prof_cas(k).and. plev(l).gt.plev_prof_cas(k+1)) then
+            k1=k
+            k2=k+1
+          endif
+         enddo
+
+         if (k1.eq.0 .or. k2.eq.0) then
+          write(*,*) 'PB! k1, k2 = ',k1,k2
+          write(*,*) 'l,plev(l) = ',l,plev(l)/100
+         do k = 1, nlev_cas-1
+          write(*,*) 'k,plev_prof_cas(k) = ',k,plev_prof_cas(k)/100
+         enddo
+         endif
+
+         frac = (plev_prof_cas(k2)-plev(l))/(plev_prof_cas(k2)-plev_prof_cas(k1))
+         tke_mod_cas(l)= tke_prof_cas(k2) - frac*(tke_prof_cas(k2)-tke_prof_cas(k1))
+         else !play>plev_prof_cas(1)
+         k1=1
+         k2=2
+         tke_mod_cas(l)= frac1*tke_prof_cas(k1) - frac2*tke_prof_cas(k2)
+
+         endif ! plev.le.plev_prof_cas(1)
+
+        else ! above max altitude of forcing file
+
+         tke_mod_cas(l)=0.0
+
+        endif ! plev
+
+       enddo ! l
+
+
+
           return
-          end
+          end SUBROUTINE interp2_case_vertical_std
 !***************************************************************************** 
 
Index: LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/old_1D_decl_cases.h
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/old_1D_decl_cases.h	(revision 3605)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/old_1D_decl_cases.h	(revision 3798)
@@ -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/branches/Ocean_skin/libf/phylmd/dyn1d/old_1D_interp_cases.h
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/old_1D_interp_cases.h	(revision 3605)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/old_1D_interp_cases.h	(revision 3798)
@@ -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/branches/Ocean_skin/libf/phylmd/dyn1d/old_1D_read_forc_cases.h
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/old_1D_read_forc_cases.h	(revision 3605)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/old_1D_read_forc_cases.h	(revision 3798)
@@ -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
@@ -1015,6 +1017,6 @@
       if (forcing_SCM) then
 
-         write(*,*),'avant call read_SCM'
-         call read_SCM_cas
+         write(*,*),'avant call old_read_SCM_cas'
+         call old_read_SCM_cas
          write(*,*) 'Forcing read' 
 
@@ -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/branches/Ocean_skin/libf/phylmd/dyn1d/old_lmdz1d.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/old_lmdz1d.F90	(revision 3605)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/old_lmdz1d.F90	(revision 3798)
@@ -632,5 +632,7 @@
 !          (phys_state_var_init is called again in physiq)
       read_climoz = 0
-!
+      nsw=6          ! EV et LF: sinon, falb_dir et falb_dif ne peuvent etre alloues
+
+
       call phys_state_var_init(read_climoz)
 
@@ -728,5 +730,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 +795,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 +846,4 @@
      end if
 
-
         print*,'nat_surf,pctsrf(1,is_oce),pctsrf(1,is_ter)',nat_surf         &
      &        ,pctsrf(1,is_oce),pctsrf(1,is_ter)
@@ -848,7 +852,4 @@
         zpic = zpicinp
         ftsol=tsurf
-        nsw=6 ! on met le nb de bandes SW=6, pour initialiser
-              ! 6 albedo, mais on peut quand meme tourner avec
-              ! moins. Seules les 2 ou 4 premiers seront lus
         falb_dir=albedo
         falb_dif=albedo
@@ -913,8 +914,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/branches/Ocean_skin/libf/phylmd/dyn1d/scm.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/scm.F90	(revision 3605)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/scm.F90	(revision 3798)
@@ -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
@@ -190,6 +166,6 @@
       real :: dt_cooling(llm),d_t_adv(llm),d_t_nudge(llm)
       real :: d_u_nudge(llm),d_v_nudge(llm)
-      real :: du_adv(llm),dv_adv(llm)
-      real :: du_age(llm),dv_age(llm)
+!      real :: d_u_adv(llm),d_v_adv(llm)
+      real :: d_u_age(llm),d_v_age(llm)
       real :: alpha
       real :: ttt
@@ -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/
@@ -273,8 +249,9 @@
       d_u_nudge(:)=0.
       d_v_nudge(:)=0.
-      du_adv(:)=0.
-      dv_adv(:)=0.
-      du_age(:)=0.
-      dv_age(:)=0.
+      d_u_adv(:)=0.
+      d_v_adv(:)=0.
+      d_u_age(:)=0.
+      d_v_age(:)=0.
+      
       
 ! Initialization of Common turb_forcing
@@ -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)
@@ -433,5 +412,6 @@
 !          (phys_state_var_init is called again in physiq)
       read_climoz = 0
-!
+      nsw=6
+
       call phys_state_var_init(read_climoz)
 
@@ -446,10 +426,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 +443,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 +466,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 +497,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 +518,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 +537,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, &
@@ -650,7 +625,4 @@
         zpic = zpicinp
         ftsol=tsurf
-        nsw=6 ! on met le nb de bandes SW=6, pour initialiser
-              ! 6 albedo, mais on peut quand meme tourner avec
-              ! moins. Seules les 2 ou 4 premiers seront lus
         falb_dir=albedo
         falb_dif=albedo
@@ -664,16 +636,9 @@
         prw_ancien = 0.
 !jyg<
-!!        pbl_tke(:,:,:)=1.e-8
-        pbl_tke(:,:,:)=0.
-        pbl_tke(:,2,:)=1.e-2
-        PRINT *, ' pbl_tke dans lmdz1d '
-        if (prt_level .ge. 5) then 
-         DO nsrf = 1,4
-           PRINT *,'pbl_tke(1,:,',nsrf,') ',pbl_tke(1,:,nsrf)
-         ENDDO
-        end if
-
+! Etienne: comment those lines since now the TKE is inialized in 1D_read_forc_cases
+!!      pbl_tke(:,:,:)=1.e-8
+!        pbl_tke(:,:,:)=0.
+!        pbl_tke(:,2,:)=1.e-2
 !>jyg
-
         rain_fall=0.
         snow_fall=0.
@@ -715,8 +680,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 +703,5 @@
 ! to be set at some arbitratry convenient values.
 !------------------------------------------------------------------------
-!Al1 =============== restart option ==========================
+!Al1 =============== restart option ======================================
         if (.not.restart) then
           iflag_pbl = 5
@@ -803,5 +768,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 +774,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 +785,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 +798,4 @@
          print*,'PAS DE TEMPS ',timestep
        endif
-!Al1 demande de restartphy.nc
        if (it.eq.it_end) lastcall=.True.
 
@@ -840,19 +805,17 @@
 
 #include "1D_interp_cases.h"
-   ! Vertical advection
-!  call lstendH(llm,nqtot,omega,d_t_vert_adv,d_q_vert_adv,q,temp,u,v,play)
-!  print*,'B d_t_adv ',d_t_adv(1:20)*86400
-!  print*,'B d_t_vert_adv ',d_t_vert_adv(1:20)*86400
-!  print*,'B dt omega ',omega
 
 !---------------------------------------------------------------------
 !  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
+
 
 !---------------------------------------------------------------------
@@ -872,9 +835,13 @@
       teta=temp*(pzero/play)**rkappa
       do l=2,llm-1
+        ! vertical tendencies computed as d X / d t = -W  d X / d z
         d_u_vert_adv(l)=-w_adv(l)*(u(l+1)-u(l-1))/(z_adv(l+1)-z_adv(l-1))
         d_v_vert_adv(l)=-w_adv(l)*(v(l+1)-v(l-1))/(z_adv(l+1)-z_adv(l-1))
-        d_t_vert_adv(l)=-(w_adv(l)*(teta(l+1)-teta(l-1))/(z_adv(l+1)-z_adv(l-1)))/(pzero/play(l))**rkappa
+        ! d theta / dt = -W d theta / d z, transformed into d temp / d t dividing by (pzero/play(l))**rkappa
+        d_t_vert_adv(l)=-w_adv(l)*(teta(l+1)-teta(l-1))/(z_adv(l+1)-z_adv(l-1)) / (pzero/play(l))**rkappa
         d_q_vert_adv(l,1)=-w_adv(l)*(q(l+1,1)-q(l-1,1))/(z_adv(l+1)-z_adv(l-1))
       enddo
+      d_u_adv(:)=d_u_adv(:)+d_u_vert_adv(:)
+      d_v_adv(:)=d_v_adv(:)+d_v_vert_adv(:)
       d_t_adv(:)=d_t_adv(:)+d_t_vert_adv(:)
       d_q_adv(:,1)=d_q_adv(:,1)+d_q_vert_adv(:,1)
@@ -938,90 +905,32 @@
        fcoriolis=2.*sin(rpi*xlat/180.)*romega
 
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!       if (forcing_radconv .or. forcing_fire) then
-!         fcoriolis=0.0
-!         dt_cooling=0.0
-!         d_t_adv=0.0
-!         d_q_adv=0.0
-!       endif
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!      if (forcing_toga .or. forcing_GCSSold .or. forcing_twpice            &
-!    &    .or.forcing_amma .or. forcing_type.eq.101) then
-!        fcoriolis=0.0 ; ug=0. ; vg=0.
-!      endif
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!      if(forcing_rico) then
-!         dt_cooling=0.
-!      endif
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!CRio:Attention modif sp??cifique cas de Caroline
-!     if (forcing_type==-1) then
-!        fcoriolis=0.
-!       
-!on calcule dt_cooling
-!       do l=1,llm
-!       if (play(l).ge.20000.) then
-!           dt_cooling(l)=-1.5/86400.
-!       elseif ((play(l).ge.10000.).and.((play(l).lt.20000.))) then
-!           dt_cooling(l)=-1.5/86400.*(play(l)-10000.)/(10000.)-1./86400.*(20000.-play(l))/10000.*(temp(l)-200.)
-!       else
-!           dt_cooling(l)=-1.*(temp(l)-200.)/86400.
-!       endif
-!       enddo
-!
-!     endif     
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!     if (forcing_sandu) then
-!        ug(1:llm)=u_mod(1:llm)
-!        vg(1:llm)=v_mod(1:llm)
-!     endif
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
       IF (prt_level >= 5) print*, 'fcoriolis, xlat,mxcalc ', &
                                    fcoriolis, xlat,mxcalc
 
-!       print *,'u-ug=',u-ug
-
-!!!!!!!!!!!!!!!!!!!!!!!!
-! Geostrophic wind
-! Le calcul ci dessous est insuffisamment precis
-!      du_age(1:mxcalc)=fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
-!      dv_age(1:mxcalc)=-fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
-!!!!!!!!!!!!!!!!!!!!!!!!
+!---------------------------------------------------------------------
+! Geostrophic forcing 
+!---------------------------------------------------------------------
+
+      IF ( forc_geo == 0 ) THEN
+              d_u_age(1:mxcalc)=0.
+              d_v_age(1:mxcalc)=0.
+      ELSE
        sfdt = sin(0.5*fcoriolis*timestep)
        cfdt = cos(0.5*fcoriolis*timestep)
-!       print *,'fcoriolis,sfdt,cfdt,timestep',fcoriolis,sfdt,cfdt,timestep
-!
-        du_age(1:mxcalc)= -2.*sfdt/timestep*                                &
+
+        d_u_age(1:mxcalc)= -2.*sfdt/timestep*                                &
      &          (sfdt*(u(1:mxcalc)-ug(1:mxcalc)) -                          &
      &           cfdt*(v(1:mxcalc)-vg(1:mxcalc))  )
 !!     : fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
 !
-       dv_age(1:mxcalc)= -2.*sfdt/timestep*                                 &
+       d_v_age(1:mxcalc)= -2.*sfdt/timestep*                                 &
      &          (cfdt*(u(1:mxcalc)-ug(1:mxcalc)) +                           &
      &           sfdt*(v(1:mxcalc)-vg(1:mxcalc))  )
 !!     : -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
-!
-!!!!!!!!!!!!!!!!!!!!!!!!
+      ENDIF
+!
+!---------------------------------------------------------------------
 !  Nudging
-!!!!!!!!!!!!!!!!!!!!!!!!
+!---------------------------------------------------------------------
       d_t_nudge(:) = 0.
       d_u_nudge(:) = 0.
@@ -1039,4 +948,7 @@
       ENDDO
 
+!---------------------------------------------------------------------
+!  Optional outputs
+!---------------------------------------------------------------------
 #ifdef OUTPUT_PHYS_SCM
       CALL iophys_ecrit('w_adv',klev,'w_adv','K/day',w_adv)
@@ -1056,62 +968,13 @@
 #endif
 
-!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!       if (forcing_fire) THEN
-!               print*,'Enlever cette section rapidement'
-!               stop
-!               
-!
-!!let ww=if ( alt le 1100 ) then alt*-0.00001 else 0
-!!let wt=if ( alt le 1100 ) then min( -3.75e-5 , -7.5e-8*alt)  else 0
-!!let wq=if ( alt le 1100 ) then max( 1.5e-8 , 3e-11*alt)  else 0
-!           d_t_adv=0.
-!           d_q_adv=0.
-!           teta=temp*(pzero/play)**rkappa
-!           d_t_adv=0.
-!           d_q_adv=0.
-!           do l=2,llm-1
-!              if (zlay(l)<=1100) then
-!                  wwww=-0.00001*zlay(l)
-!                  d_t_adv(l)=-wwww*(teta(l)-teta(l+1))/(zlay(l)-zlay(l+1)) /(pzero/play(l))**rkappa
-!                  d_q_adv(l,1:2)=-wwww*(q(l,1:2)-q(l+1,1:2))/(zlay(l)-zlay(l+1))
-!                  d_t_adv(l)=d_t_adv(l)+min(-3.75e-5 , -7.5e-8*zlay(l))
-!                  d_q_adv(l,1)=d_q_adv(l,1)+max( 1.5e-8 , 3e-11*zlay(l))
-!              endif
-!           enddo
-!
-!        endif
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!         call  writefield_phy('dv_age' ,dv_age,llm)
-!         call  writefield_phy('du_age' ,du_age,llm)
-!         call  writefield_phy('du_phys' ,du_phys,llm)
-!         call  writefield_phy('u_tend' ,u,llm)
-!         call  writefield_phy('u_g' ,ug,llm)
-!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!! Increment state variables
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     IF (flag_inhib_forcing == 0) then ! if tendency of forcings should be added
 
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! pour les cas sandu et astex, on reclacule u,v,q,temp et teta dans 1D_nudge_sandu_astex.h
-! au dessus de 700hpa, on relaxe vers les profils initiaux
-!     if (forcing_sandu .OR. forcing_astex) then
-!#include "1D_nudge_sandu_astex.h"
-!      else
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         u(1:mxcalc)=u(1:mxcalc) + timestep*(                                &
      &              du_phys(1:mxcalc)                                       &
-     &             +du_age(1:mxcalc)+du_adv(1:mxcalc)                       &
+     &             +d_u_age(1:mxcalc)+d_u_adv(1:mxcalc)                       &
      &             +d_u_nudge(1:mxcalc) )            
         v(1:mxcalc)=v(1:mxcalc) + timestep*(                                 &
      &              dv_phys(1:mxcalc)                                       &
-     &             +dv_age(1:mxcalc)+dv_adv(1:mxcalc)                       &
+     &             +d_v_age(1:mxcalc)+d_v_adv(1:mxcalc)                       &
      &             +d_v_nudge(1:mxcalc) )
         q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*(                              &
@@ -1125,6 +988,6 @@
      &              temp(1),dt_phys(1),d_t_adv(1),dt_cooling(1)
            print* ,'dv_phys=',dv_phys
-           print* ,'dv_age=',dv_age
-           print* ,'dv_adv=',dv_adv
+           print* ,'d_v_age=',d_v_age
+           print* ,'d_v_adv=',d_v_adv
            print* ,'d_v_nudge=',d_v_nudge
            print*, v
@@ -1134,17 +997,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
@@ -1184,7 +1045,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)
@@ -1192,14 +1051,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,                            &
