Index: LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/lmdz_1dutils.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/lmdz_1dutils.f90	(revision 5143)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/lmdz_1dutils.f90	(revision 5144)
@@ -1172,6 +1172,7 @@
     !   sans WTG rajouter une advection horizontale
     !----------------------------------------------------------------------
+    USE lmdz_yomcst
+
     IMPLICIT NONE
-    include "YOMCST.h"
     !        argument
     INTEGER llm
@@ -1244,6 +1245,7 @@
     !   sans WTG rajouter une advection horizontale
     !----------------------------------------------------------------------
+    USE lmdz_yomcst
+
     IMPLICIT NONE
-    include "YOMCST.h"
     !        argument
     INTEGER llm, nqtot
@@ -1318,6 +1320,7 @@
     ! ========================================================
     USE dimphy
-    USE lmdz_YOETHF
+    USE lmdz_yoethf
     USE lmdz_fcttre, ONLY: foeew, foede, qsats, qsatl, dqsats, dqsatl, thermcep
+    USE lmdz_yomcst
 
     IMPLICIT NONE
@@ -1337,8 +1340,4 @@
     INTEGER k, i
     REAL zx_qs
-
-    ! Declaration des constantes et des fonctions thermodynamiques
-
-    include "YOMCST.h"
 
     DO k = 1, klev
@@ -1394,6 +1393,7 @@
     ! ========================================================
     USE dimphy
-    USE lmdz_YOETHF
+    USE lmdz_yoethf
     USE lmdz_fcttre, ONLY: foeew, foede, qsats, qsatl, dqsats, dqsatl, thermcep
+    USE lmdz_yomcst
 
     IMPLICIT NONE
@@ -1424,8 +1424,4 @@
     INTEGER k, i
     REAL zx_qs, rh, tnew, d_rh, rhnew
-
-    ! Declaration des constantes et des fonctions thermodynamiques
-
-    include "YOMCST.h"
 
     print *, 'dtime, tau ', dtime, tau
@@ -1545,7 +1541,8 @@
           &, dth_mod_cas, hth_mod_cas, vth_mod_cas, mxcalc)
 
+    USE lmdz_yomcst
+
     IMPLICIT NONE
 
-    include "YOMCST.h"
     include "dimensions.h"
 
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/lmdz_old_1dconv.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/lmdz_old_1dconv.f90	(revision 5143)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/lmdz_old_1dconv.f90	(revision 5144)
@@ -1,4 +1,4 @@
 MODULE lmdz_old_1dconv
-   PRIVATE  ! -- We'd love to put IMPLICIT NONE;  here...
+  PRIVATE  ! -- We'd love to put IMPLICIT NONE;  here...
   PUBLIC get_uvd, copie, get_uvd2, rdgrads, spaces
 CONTAINS
@@ -7,4 +7,6 @@
           &       ht, hq, hw, hu, hv, hthturb, hqturb, &
           &       Ts, imp_fcg, ts_fcg, Tp_fcg, Turb_fcg)
+
+    USE lmdz_yomcst
 
     IMPLICIT NONE
@@ -14,6 +16,4 @@
     ! pouvoir calculer la convergence et le cisaillement dans la physiq
     !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-    INCLUDE "YOMCST.h"
 
     INTEGER klev
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/lmdz_old_lmdz1d.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/lmdz_old_lmdz1d.F90	(revision 5143)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/lmdz_old_lmdz1d.F90	(revision 5144)
@@ -56,7 +56,7 @@
     USE lmdz_fcs_gcssold, ONLY: imp_fcg_gcssold, ts_fcg_gcssold, Tp_fcg_gcssold, Tp_ini_gcssold, xTurb_fcg_gcssold
     USE lmdz_tsoilnudge, ONLY: nudge_tsoil, isoil_nudge, Tsoil_nudge, tau_soil_nudge
+    USE lmdz_yomcst
 
     INCLUDE "dimensions.h"
-    INCLUDE "YOMCST.h"
     INCLUDE "dimsoil.h"
     INCLUDE "compar1d.h"
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/lmdz_scm.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/lmdz_scm.F90	(revision 5143)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/lmdz_scm.F90	(revision 5144)
@@ -49,7 +49,7 @@
     USE lmdz_fcs_gcssold, ONLY: imp_fcg_gcssold, ts_fcg_gcssold, Tp_fcg_gcssold, Tp_ini_gcssold, xTurb_fcg_gcssold
     USE lmdz_tsoilnudge, ONLY: nudge_tsoil, isoil_nudge, Tsoil_nudge, tau_soil_nudge
+    USE lmdz_yomcst
 
     INCLUDE "dimensions.h"
-    INCLUDE "YOMCST.h"
     INCLUDE "dimsoil.h"
     INCLUDE "compar1d.h"
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/mod_1D_cases_read_std.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/mod_1D_cases_read_std.F90	(revision 5143)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/mod_1D_cases_read_std.F90	(revision 5144)
@@ -1,85 +1,82 @@
-
 ! $Id: mod_1D_cases_read.F90 2373 2015-10-13 17:28:01Z jyg $
 
 MODULE mod_1D_cases_read_std
-  USE netcdf, ONLY:nf90_noerr,nf90_inq_varid,nf90_inq_dimid,nf90_inquire_dimension,nf90_open,nf90_nowrite,&
-          nf90_strerror,nf90_get_var
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  USE netcdf, ONLY: nf90_noerr, nf90_inq_varid, nf90_inq_dimid, nf90_inquire_dimension, nf90_open, nf90_nowrite, &
+          nf90_strerror, nf90_get_var
+
+  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !Declarations specifiques au cas standard
   CHARACTER*80 :: fich_cas
-  ! Discr?tisation 
+  ! Discr?tisation
   INTEGER nlev_cas, nt_cas
 
 
   !profils environnementaux
-  REAL, ALLOCATABLE::  plev_cas(:,:),plevh_cas(:)
-  REAL, ALLOCATABLE::  ap_cas(:),bp_cas(:)
-
-  REAL, ALLOCATABLE::  z_cas(:,:),zh_cas(:)
-  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(:,:),tke_cas(:,:)
+  REAL, ALLOCATABLE :: plev_cas(:, :), plevh_cas(:)
+  REAL, ALLOCATABLE :: ap_cas(:), bp_cas(:)
+
+  REAL, ALLOCATABLE :: z_cas(:, :), zh_cas(:)
+  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(:, :), tke_cas(:, :)
 
   !forcing
-  REAL, ALLOCATABLE::  ht_cas(:,:),vt_cas(:,:),dt_cas(:,:),dtrad_cas(:,:)
-  REAL, ALLOCATABLE::  hth_cas(:,:),vth_cas(:,:),dth_cas(:,:)
-  REAL, ALLOCATABLE::  hq_cas(:,:),vq_cas(:,:),dq_cas(:,:)
-  REAL, ALLOCATABLE::  hr_cas(:,:),vr_cas(:,:),dr_cas(:,:)
-  REAL, ALLOCATABLE::  hu_cas(:,:),vu_cas(:,:),du_cas(:,:)
-  REAL, ALLOCATABLE::  hv_cas(:,:),vv_cas(:,:),dv_cas(:,:)
-  REAL, ALLOCATABLE::  ug_cas(:,:),vg_cas(:,:)
-  REAL, ALLOCATABLE::  temp_nudg_cas(:,:),qv_nudg_cas(:,:),u_nudg_cas(:,:),v_nudg_cas(:,:)
-  REAL, ALLOCATABLE::  invtau_temp_nudg_cas(:,:),invtau_qv_nudg_cas(:,:),invtau_u_nudg_cas(:,:),invtau_v_nudg_cas(:,:)
-  REAL, ALLOCATABLE::  lat_cas(:),sens_cas(:),tskin_cas(:),ts_cas(:),ps_cas(:),ustar_cas(:)
-  REAL, ALLOCATABLE::  uw_cas(:,:),vw_cas(:,:),q1_cas(:,:),q2_cas(:,:),tkes_cas(:)
+  REAL, ALLOCATABLE :: ht_cas(:, :), vt_cas(:, :), dt_cas(:, :), dtrad_cas(:, :)
+  REAL, ALLOCATABLE :: hth_cas(:, :), vth_cas(:, :), dth_cas(:, :)
+  REAL, ALLOCATABLE :: hq_cas(:, :), vq_cas(:, :), dq_cas(:, :)
+  REAL, ALLOCATABLE :: hr_cas(:, :), vr_cas(:, :), dr_cas(:, :)
+  REAL, ALLOCATABLE :: hu_cas(:, :), vu_cas(:, :), du_cas(:, :)
+  REAL, ALLOCATABLE :: hv_cas(:, :), vv_cas(:, :), dv_cas(:, :)
+  REAL, ALLOCATABLE :: ug_cas(:, :), vg_cas(:, :)
+  REAL, ALLOCATABLE :: temp_nudg_cas(:, :), qv_nudg_cas(:, :), u_nudg_cas(:, :), v_nudg_cas(:, :)
+  REAL, ALLOCATABLE :: invtau_temp_nudg_cas(:, :), invtau_qv_nudg_cas(:, :), invtau_u_nudg_cas(:, :), invtau_v_nudg_cas(:, :)
+  REAL, ALLOCATABLE :: lat_cas(:), sens_cas(:), tskin_cas(:), ts_cas(:), ps_cas(:), ustar_cas(:)
+  REAL, ALLOCATABLE :: uw_cas(:, :), vw_cas(:, :), q1_cas(:, :), q2_cas(:, :), tkes_cas(:)
 
   !champs interpoles
-  REAL, ALLOCATABLE::  plev_prof_cas(:)
-  REAL, ALLOCATABLE::  t_prof_cas(:)
-  REAL, ALLOCATABLE::  theta_prof_cas(:)
-  REAL, ALLOCATABLE::  thl_prof_cas(:)
-  REAL, ALLOCATABLE::  thv_prof_cas(:)
-  REAL, ALLOCATABLE::  q_prof_cas(:)
-  REAL, ALLOCATABLE::  qv_prof_cas(:)
-  REAL, ALLOCATABLE::  ql_prof_cas(:)
-  REAL, ALLOCATABLE::  qi_prof_cas(:)
-  REAL, ALLOCATABLE::  rh_prof_cas(:)
-  REAL, ALLOCATABLE::  rv_prof_cas(:)
-  REAL, ALLOCATABLE::  u_prof_cas(:)
-  REAL, ALLOCATABLE::  v_prof_cas(:)
-  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(:)
-  REAL, ALLOCATABLE::  temp_nudg_prof_cas(:),qv_nudg_prof_cas(:),u_nudg_prof_cas(:),v_nudg_prof_cas(:)
-  REAL, ALLOCATABLE::  invtau_temp_nudg_prof_cas(:),invtau_qv_nudg_prof_cas(:),invtau_u_nudg_prof_cas(:),invtau_v_nudg_prof_cas(:)
-
-  REAL, ALLOCATABLE::  ht_prof_cas(:)
-  REAL, ALLOCATABLE::  hth_prof_cas(:)
-  REAL, ALLOCATABLE::  hq_prof_cas(:)
-  REAL, ALLOCATABLE::  vt_prof_cas(:)
-  REAL, ALLOCATABLE::  vth_prof_cas(:)
-  REAL, ALLOCATABLE::  vq_prof_cas(:)
-  REAL, ALLOCATABLE::  dt_prof_cas(:)
-  REAL, ALLOCATABLE::  dth_prof_cas(:)
-  REAL, ALLOCATABLE::  dtrad_prof_cas(:)
-  REAL, ALLOCATABLE::  dq_prof_cas(:)
-  REAL, ALLOCATABLE::  hu_prof_cas(:)
-  REAL, ALLOCATABLE::  hv_prof_cas(:)
-  REAL, ALLOCATABLE::  vu_prof_cas(:)
-  REAL, ALLOCATABLE::  vv_prof_cas(:)
-  REAL, ALLOCATABLE::  du_prof_cas(:)
-  REAL, ALLOCATABLE::  dv_prof_cas(:)
-  REAL, ALLOCATABLE::  uw_prof_cas(:)
-  REAL, ALLOCATABLE::  vw_prof_cas(:)
-  REAL, ALLOCATABLE::  q1_prof_cas(:)
-  REAL, ALLOCATABLE::  q2_prof_cas(:)
-
-
-  REAL o3_cas,lat_prof_cas,sens_prof_cas,ts_prof_cas,tskin_prof_cas,ps_prof_cas,ustar_prof_cas,tkes_prof_cas
-  REAL orog_cas,albedo_cas,emiss_cas,q_skin_cas,mom_rough,heat_rough,rugos_cas,sand_cas,clay_cas
-
+  REAL, ALLOCATABLE :: plev_prof_cas(:)
+  REAL, ALLOCATABLE :: t_prof_cas(:)
+  REAL, ALLOCATABLE :: theta_prof_cas(:)
+  REAL, ALLOCATABLE :: thl_prof_cas(:)
+  REAL, ALLOCATABLE :: thv_prof_cas(:)
+  REAL, ALLOCATABLE :: q_prof_cas(:)
+  REAL, ALLOCATABLE :: qv_prof_cas(:)
+  REAL, ALLOCATABLE :: ql_prof_cas(:)
+  REAL, ALLOCATABLE :: qi_prof_cas(:)
+  REAL, ALLOCATABLE :: rh_prof_cas(:)
+  REAL, ALLOCATABLE :: rv_prof_cas(:)
+  REAL, ALLOCATABLE :: u_prof_cas(:)
+  REAL, ALLOCATABLE :: v_prof_cas(:)
+  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(:)
+  REAL, ALLOCATABLE :: temp_nudg_prof_cas(:), qv_nudg_prof_cas(:), u_nudg_prof_cas(:), v_nudg_prof_cas(:)
+  REAL, ALLOCATABLE :: invtau_temp_nudg_prof_cas(:), invtau_qv_nudg_prof_cas(:), invtau_u_nudg_prof_cas(:), invtau_v_nudg_prof_cas(:)
+
+  REAL, ALLOCATABLE :: ht_prof_cas(:)
+  REAL, ALLOCATABLE :: hth_prof_cas(:)
+  REAL, ALLOCATABLE :: hq_prof_cas(:)
+  REAL, ALLOCATABLE :: vt_prof_cas(:)
+  REAL, ALLOCATABLE :: vth_prof_cas(:)
+  REAL, ALLOCATABLE :: vq_prof_cas(:)
+  REAL, ALLOCATABLE :: dt_prof_cas(:)
+  REAL, ALLOCATABLE :: dth_prof_cas(:)
+  REAL, ALLOCATABLE :: dtrad_prof_cas(:)
+  REAL, ALLOCATABLE :: dq_prof_cas(:)
+  REAL, ALLOCATABLE :: hu_prof_cas(:)
+  REAL, ALLOCATABLE :: hv_prof_cas(:)
+  REAL, ALLOCATABLE :: vu_prof_cas(:)
+  REAL, ALLOCATABLE :: vv_prof_cas(:)
+  REAL, ALLOCATABLE :: du_prof_cas(:)
+  REAL, ALLOCATABLE :: dv_prof_cas(:)
+  REAL, ALLOCATABLE :: uw_prof_cas(:)
+  REAL, ALLOCATABLE :: vw_prof_cas(:)
+  REAL, ALLOCATABLE :: q1_prof_cas(:)
+  REAL, ALLOCATABLE :: q2_prof_cas(:)
+
+  REAL o3_cas, lat_prof_cas, sens_prof_cas, ts_prof_cas, tskin_prof_cas, ps_prof_cas, ustar_prof_cas, tkes_prof_cas
+  REAL orog_cas, albedo_cas, emiss_cas, q_skin_cas, mom_rough, heat_rough, rugos_cas, sand_cas, clay_cas
 
 
@@ -93,94 +90,94 @@
     INCLUDE "date_cas.h"
 
-    INTEGER nid,rid,ierr
-    INTEGER ii,jj,timeid
+    INTEGER nid, rid, ierr
+    INTEGER ii, jj, timeid
     REAL, ALLOCATABLE :: time_val(:)
 
-    fich_cas='cas.nc'
-    PRINT*,'fich_cas ',fich_cas
-    ierr = nf90_open(fich_cas,nf90_nowrite,nid)
-    PRINT*,'fich_cas,nf90_nowrite,nid ',fich_cas,nf90_nowrite,nid
+    fich_cas = 'cas.nc'
+    PRINT*, 'fich_cas ', fich_cas
+    ierr = nf90_open(fich_cas, nf90_nowrite, nid)
+    PRINT*, 'fich_cas,nf90_nowrite,nid ', fich_cas, nf90_nowrite, nid
     IF (ierr/=nf90_noerr) THEN
-       WRITE(*,*) 'ERROR: GROS Pb opening forcings nc file '
-       WRITE(*,*) nf90_strerror(ierr)
-       stop ""
+      WRITE(*, *) 'ERROR: GROS Pb opening forcings nc file '
+      WRITE(*, *) nf90_strerror(ierr)
+      stop ""
     endif
     !.......................................................................
-    ierr=nf90_inq_dimid(nid,'lat',rid)
+    ierr = nf90_inq_dimid(nid, 'lat', rid)
     IF (ierr/=nf90_noerr) THEN
-       PRINT*, 'Oh probleme lecture dimension lat'
+      PRINT*, 'Oh probleme lecture dimension lat'
     ENDIF
-    ierr=nf90_inquire_dimension(nid,rid,len=ii)
-    PRINT*,'OK1 read_SCM_cas: nid,rid,lat',nid,rid,ii
+    ierr = nf90_inquire_dimension(nid, rid, len = ii)
+    PRINT*, 'OK1 read_SCM_cas: nid,rid,lat', nid, rid, ii
     !.......................................................................
-    ierr=nf90_inq_dimid(nid,'lon',rid)
+    ierr = nf90_inq_dimid(nid, 'lon', rid)
     IF (ierr/=nf90_noerr) THEN
-       PRINT*, 'Oh probleme lecture dimension lon'
+      PRINT*, 'Oh probleme lecture dimension lon'
     ENDIF
-    ierr=nf90_inquire_dimension(nid,rid,len=jj)
-    PRINT*,'OK2 read_SCM_cas: nid,rid,lat',nid,rid,jj
+    ierr = nf90_inquire_dimension(nid, rid, len = jj)
+    PRINT*, 'OK2 read_SCM_cas: nid,rid,lat', nid, rid, jj
     !.......................................................................
-    ierr=nf90_inq_dimid(nid,'lev',rid)
+    ierr = nf90_inq_dimid(nid, 'lev', rid)
     IF (ierr/=nf90_noerr) THEN
-       PRINT*, 'Oh probleme lecture dimension nlev'
+      PRINT*, 'Oh probleme lecture dimension nlev'
     ENDIF
-    ierr=nf90_inquire_dimension(nid,rid,len=nlev_cas)
-    PRINT*,'OK3 read_SCM_cas: nid,rid,nlev_cas',nid,rid,nlev_cas
-    IF ( .NOT. ( nlev_cas > 10 .AND. nlev_cas < 200000 )) THEN
-       PRINT*,'Valeur de nlev_cas peu probable'
-       STOP
+    ierr = nf90_inquire_dimension(nid, rid, len = nlev_cas)
+    PRINT*, 'OK3 read_SCM_cas: nid,rid,nlev_cas', nid, rid, nlev_cas
+    IF (.NOT. (nlev_cas > 10 .AND. nlev_cas < 200000)) THEN
+      PRINT*, 'Valeur de nlev_cas peu probable'
+      STOP
     ENDIF
     !.......................................................................
-    ierr=nf90_inq_dimid(nid,'time',rid)
-    nt_cas=0
+    ierr = nf90_inq_dimid(nid, 'time', rid)
+    nt_cas = 0
     IF (ierr/=nf90_noerr) THEN
-       stop 'Oh probleme lecture dimension time'
+      stop 'Oh probleme lecture dimension time'
     ENDIF
-    ierr=nf90_inquire_dimension(nid,rid,len=nt_cas)
-    PRINT*,'OK4 read_SCM_cas: nid,rid,nt_cas',nid,rid,nt_cas
+    ierr = nf90_inquire_dimension(nid, rid, len = nt_cas)
+    PRINT*, 'OK4 read_SCM_cas: nid,rid,nt_cas', nid, rid, nt_cas
     ! Lecture de l'axe des temps
-    PRINT*,'LECTURE DU TEMPS'
-    ierr=nf90_inq_varid(nid,'time',timeid)
+    PRINT*, 'LECTURE DU TEMPS'
+    ierr = nf90_inq_varid(nid, 'time', timeid)
     IF(ierr/=nf90_noerr) THEN
-       print *,'Variable time manquante dans cas.nc:'
-       ierr=nf90_noerr
+      print *, 'Variable time manquante dans cas.nc:'
+      ierr = nf90_noerr
     else
-       allocate(time_val(nt_cas))
-       ierr = nf90_get_var(nid,timeid,time_val)
-       IF(ierr/=nf90_noerr) THEN
-          print *,'A Pb a la lecture de time cas.nc: '
-       endif
+      allocate(time_val(nt_cas))
+      ierr = nf90_get_var(nid, timeid, time_val)
+      IF(ierr/=nf90_noerr) THEN
+        print *, 'A Pb a la lecture de time cas.nc: '
+      endif
     endif
     IF (nt_cas>1) THEN
-       pdt_cas=time_val(2)-time_val(1)
+      pdt_cas = time_val(2) - time_val(1)
     ELSE
-       pdt_cas=0.
+      pdt_cas = 0.
     ENDIF
 
 
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     !profils moyens:
-    allocate(plev_cas(nlev_cas,nt_cas),plevh_cas(nlev_cas+1))        
-    allocate(z_cas(nlev_cas,nt_cas),zh_cas(nlev_cas+1))
-    allocate(ap_cas(nlev_cas+1),bp_cas(nlev_cas+1))
-    allocate(t_cas(nlev_cas,nt_cas),q_cas(nlev_cas,nt_cas),qv_cas(nlev_cas,nt_cas),ql_cas(nlev_cas,nt_cas), &
-         qi_cas(nlev_cas,nt_cas),rh_cas(nlev_cas,nt_cas))
-    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))
+    allocate(plev_cas(nlev_cas, nt_cas), plevh_cas(nlev_cas + 1))
+    allocate(z_cas(nlev_cas, nt_cas), zh_cas(nlev_cas + 1))
+    allocate(ap_cas(nlev_cas + 1), bp_cas(nlev_cas + 1))
+    allocate(t_cas(nlev_cas, nt_cas), q_cas(nlev_cas, nt_cas), qv_cas(nlev_cas, nt_cas), ql_cas(nlev_cas, nt_cas), &
+            qi_cas(nlev_cas, nt_cas), rh_cas(nlev_cas, nt_cas))
+    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))
-    allocate(hq_cas(nlev_cas,nt_cas),vq_cas(nlev_cas,nt_cas),dq_cas(nlev_cas,nt_cas))
-    allocate(hth_cas(nlev_cas,nt_cas),vth_cas(nlev_cas,nt_cas),dth_cas(nlev_cas,nt_cas))
-    allocate(hr_cas(nlev_cas,nt_cas),vr_cas(nlev_cas,nt_cas),dr_cas(nlev_cas,nt_cas))
-    allocate(hu_cas(nlev_cas,nt_cas),vu_cas(nlev_cas,nt_cas),du_cas(nlev_cas,nt_cas))
-    allocate(hv_cas(nlev_cas,nt_cas),vv_cas(nlev_cas,nt_cas),dv_cas(nlev_cas,nt_cas))
-    allocate(ug_cas(nlev_cas,nt_cas),vg_cas(nlev_cas,nt_cas))
-    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(invtau_temp_nudg_cas(nlev_cas,nt_cas),invtau_qv_nudg_cas(nlev_cas,nt_cas))
-    allocate(invtau_u_nudg_cas(nlev_cas,nt_cas),invtau_v_nudg_cas(nlev_cas,nt_cas))
-    allocate(lat_cas(nt_cas),sens_cas(nt_cas),ts_cas(nt_cas),tskin_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))
+    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))
+    allocate(hq_cas(nlev_cas, nt_cas), vq_cas(nlev_cas, nt_cas), dq_cas(nlev_cas, nt_cas))
+    allocate(hth_cas(nlev_cas, nt_cas), vth_cas(nlev_cas, nt_cas), dth_cas(nlev_cas, nt_cas))
+    allocate(hr_cas(nlev_cas, nt_cas), vr_cas(nlev_cas, nt_cas), dr_cas(nlev_cas, nt_cas))
+    allocate(hu_cas(nlev_cas, nt_cas), vu_cas(nlev_cas, nt_cas), du_cas(nlev_cas, nt_cas))
+    allocate(hv_cas(nlev_cas, nt_cas), vv_cas(nlev_cas, nt_cas), dv_cas(nlev_cas, nt_cas))
+    allocate(ug_cas(nlev_cas, nt_cas), vg_cas(nlev_cas, nt_cas))
+    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(invtau_temp_nudg_cas(nlev_cas, nt_cas), invtau_qv_nudg_cas(nlev_cas, nt_cas))
+    allocate(invtau_u_nudg_cas(nlev_cas, nt_cas), invtau_v_nudg_cas(nlev_cas, nt_cas))
+    allocate(lat_cas(nt_cas), sens_cas(nt_cas), ts_cas(nt_cas), tskin_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))
 
 
@@ -205,8 +202,8 @@
     allocate(ug_prof_cas(nlev_cas))
     allocate(vg_prof_cas(nlev_cas))
-    allocate(temp_nudg_prof_cas(nlev_cas),qv_nudg_prof_cas(nlev_cas))
-    allocate(u_nudg_prof_cas(nlev_cas),v_nudg_prof_cas(nlev_cas))
-    allocate(invtau_temp_nudg_prof_cas(nlev_cas),invtau_qv_nudg_prof_cas(nlev_cas))
-    allocate(invtau_u_nudg_prof_cas(nlev_cas),invtau_v_nudg_prof_cas(nlev_cas))
+    allocate(temp_nudg_prof_cas(nlev_cas), qv_nudg_prof_cas(nlev_cas))
+    allocate(u_nudg_prof_cas(nlev_cas), v_nudg_prof_cas(nlev_cas))
+    allocate(invtau_temp_nudg_prof_cas(nlev_cas), invtau_qv_nudg_prof_cas(nlev_cas))
+    allocate(invtau_u_nudg_prof_cas(nlev_cas), invtau_v_nudg_prof_cas(nlev_cas))
     allocate(ht_prof_cas(nlev_cas))
     allocate(hth_prof_cas(nlev_cas))
@@ -230,46 +227,45 @@
     allocate(q2_prof_cas(nlev_cas))
 
-    PRINT*,'Allocations OK'
-    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,tke_cas,ug_cas,vg_cas,                            &
-         temp_nudg_cas,qv_nudg_cas,u_nudg_cas,v_nudg_cas,                                                     &
-         invtau_temp_nudg_cas,invtau_qv_nudg_cas,invtau_u_nudg_cas,invtau_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,tskin_cas,ps_cas,ustar_cas,tkes_cas,            &
-         uw_cas,vw_cas,q1_cas,q2_cas,orog_cas,albedo_cas,emiss_cas,q_skin_cas,mom_rough,heat_rough, &
-         o3_cas,rugos_cas,clay_cas,sand_cas)
-    PRINT*,'read_SCM cas OK'
-    do ii=1,nlev_cas
-       PRINT*,'apres read_SCM_cas, plev_cas=',ii,plev_cas(ii,1)
-       !PRINT*,'apres read_SCM, plev_cas=',ii,omega_cas(ii,nt_cas/2+1)
+    PRINT*, 'Allocations OK'
+    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, tke_cas, ug_cas, vg_cas, &
+            temp_nudg_cas, qv_nudg_cas, u_nudg_cas, v_nudg_cas, &
+            invtau_temp_nudg_cas, invtau_qv_nudg_cas, invtau_u_nudg_cas, invtau_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, tskin_cas, ps_cas, ustar_cas, tkes_cas, &
+            uw_cas, vw_cas, q1_cas, q2_cas, orog_cas, albedo_cas, emiss_cas, q_skin_cas, mom_rough, heat_rough, &
+            o3_cas, rugos_cas, clay_cas, sand_cas)
+    PRINT*, 'read_SCM cas OK'
+    do ii = 1, nlev_cas
+      PRINT*, 'apres read_SCM_cas, plev_cas=', ii, plev_cas(ii, 1)
+      !PRINT*,'apres read_SCM, plev_cas=',ii,omega_cas(ii,nt_cas/2+1)
     enddo
 
-
   END SUBROUTINE read_SCM_cas
 
 
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   SUBROUTINE deallocate2_1D_cases
     !profils environnementaux:
-    deallocate(plev_cas,plevh_cas)
-
-    deallocate(z_cas,zh_cas)
-    deallocate(ap_cas,bp_cas)
-    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,tke_cas)
+    deallocate(plev_cas, plevh_cas)
+
+    deallocate(z_cas, zh_cas)
+    deallocate(ap_cas, bp_cas)
+    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, tke_cas)
 
     !forcing
-    deallocate(ht_cas,vt_cas,dt_cas,dtrad_cas)
-    deallocate(hq_cas,vq_cas,dq_cas)
-    deallocate(hth_cas,vth_cas,dth_cas)
-    deallocate(hr_cas,vr_cas,dr_cas)
-    deallocate(hu_cas,vu_cas,du_cas)
-    deallocate(hv_cas,vv_cas,dv_cas)
+    deallocate(ht_cas, vt_cas, dt_cas, dtrad_cas)
+    deallocate(hq_cas, vq_cas, dq_cas)
+    deallocate(hth_cas, vth_cas, dth_cas)
+    deallocate(hr_cas, vr_cas, dr_cas)
+    deallocate(hu_cas, vu_cas, du_cas)
+    deallocate(hv_cas, vv_cas, dv_cas)
     deallocate(ug_cas)
     deallocate(vg_cas)
-    deallocate(lat_cas,sens_cas,tskin_cas,ts_cas,ps_cas,ustar_cas,tkes_cas,uw_cas,vw_cas,q1_cas,q2_cas)
+    deallocate(lat_cas, sens_cas, tskin_cas, ts_cas, ps_cas, ustar_cas, tkes_cas, uw_cas, vw_cas, q1_cas, q2_cas)
 
     !champs interpoles
@@ -292,6 +288,6 @@
     deallocate(ug_prof_cas)
     deallocate(vg_prof_cas)
-    deallocate(temp_nudg_prof_cas,qv_nudg_prof_cas,u_nudg_prof_cas,v_nudg_prof_cas)
-    deallocate(invtau_temp_nudg_prof_cas,invtau_qv_nudg_prof_cas,invtau_u_nudg_prof_cas,invtau_v_nudg_prof_cas)
+    deallocate(temp_nudg_prof_cas, qv_nudg_prof_cas, u_nudg_prof_cas, v_nudg_prof_cas)
+    deallocate(invtau_temp_nudg_prof_cas, invtau_qv_nudg_prof_cas, invtau_u_nudg_prof_cas, invtau_v_nudg_prof_cas)
     deallocate(ht_prof_cas)
     deallocate(hq_prof_cas)
@@ -319,12 +315,12 @@
 
   !=====================================================================
-  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,tke,ug,vg,&
-       temp_nudg,qv_nudg,u_nudg,v_nudg,                                        &
-       invtau_temp_nudg,invtau_qv_nudg,invtau_u_nudg,invtau_v_nudg,             &
-       du,hu,vu,dv,hv,vv,dt,ht,vt,dq,hq,vq,                                    &
-       dth,hth,vth,dr,hr,vr,dtrad,sens,flat,ts,tskin,ps,ustar,tkes,uw,vw,q1,q2,       &
-       orog_cas,albedo_cas,emiss_cas,q_skin_cas,mom_rough,          &
-       heat_rough,o3_cas,rugos_cas,clay_cas,sand_cas)
+  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, tke, ug, vg, &
+          temp_nudg, qv_nudg, u_nudg, v_nudg, &
+          invtau_temp_nudg, invtau_qv_nudg, invtau_u_nudg, invtau_v_nudg, &
+          du, hu, vu, dv, hv, vv, dt, ht, vt, dq, hq, vq, &
+          dth, hth, vth, dr, hr, vr, dtrad, sens, flat, ts, tskin, ps, ustar, tkes, uw, vw, q1, q2, &
+          orog_cas, albedo_cas, emiss_cas, q_skin_cas, mom_rough, &
+          heat_rough, o3_cas, rugos_cas, clay_cas, sand_cas)
 
     !program reading forcing of the case study
@@ -332,37 +328,36 @@
     INCLUDE "compar1d.h"
 
-    INTEGER ntime,nlevel,k,t
-
-    REAL ap(nlevel+1),bp(nlevel+1)
-    REAL zz(nlevel,ntime),zzh(nlevel+1)
-    REAL pp(nlevel,ntime),pph(nlevel+1)
+    INTEGER ntime, nlevel, k, t
+
+    REAL ap(nlevel + 1), bp(nlevel + 1)
+    REAL zz(nlevel, ntime), zzh(nlevel + 1)
+    REAL pp(nlevel, ntime), pph(nlevel + 1)
     !profils initiaux
-    REAL temp0(nlevel),qv0(nlevel),ql0(nlevel),qi0(nlevel),u0(nlevel),v0(nlevel),tke0(nlevel)
+    REAL temp0(nlevel), qv0(nlevel), ql0(nlevel), qi0(nlevel), u0(nlevel), v0(nlevel), tke0(nlevel)
     REAL pp0(nlevel)
-    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),tkes(ntime)
-    REAL temp_nudg(nlevel,ntime),qv_nudg(nlevel,ntime),u_nudg(nlevel,ntime),v_nudg(nlevel,ntime)
-    REAL invtau_temp_nudg(nlevel,ntime),invtau_qv_nudg(nlevel,ntime),invtau_u_nudg(nlevel,ntime),invtau_v_nudg(nlevel,ntime)
-    REAL ug(nlevel,ntime),vg(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)
-    REAL dt(nlevel,ntime),ht(nlevel,ntime),vt(nlevel,ntime)
-    REAL dtrad(nlevel,ntime)
-    REAL dq(nlevel,ntime),hq(nlevel,ntime),vq(nlevel,ntime)
-    REAL dth(nlevel,ntime),hth(nlevel,ntime),vth(nlevel,ntime),hthl(nlevel,ntime)
-    REAL dr(nlevel,ntime),hr(nlevel,ntime),vr(nlevel,ntime)
-    REAL flat(ntime),sens(ntime),ustar(ntime)
-    REAL uw(nlevel,ntime),vw(nlevel,ntime),q1(nlevel,ntime),q2(nlevel,ntime)
-    REAL ts(ntime),tskin(ntime),ps(ntime)
-    REAL orog_cas,albedo_cas,emiss_cas,q_skin_cas,mom_rough,heat_rough,o3_cas,rugos_cas,clay_cas,sand_cas
-    REAL apbp(nlevel+1),resul(nlevel,ntime),resul1(nlevel),resul2(ntime),resul3
-
-
-    INTEGER nid, ierr,ierr1,ierr2,rid,i,int_test
+    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), tkes(ntime)
+    REAL temp_nudg(nlevel, ntime), qv_nudg(nlevel, ntime), u_nudg(nlevel, ntime), v_nudg(nlevel, ntime)
+    REAL invtau_temp_nudg(nlevel, ntime), invtau_qv_nudg(nlevel, ntime), invtau_u_nudg(nlevel, ntime), invtau_v_nudg(nlevel, ntime)
+    REAL ug(nlevel, ntime), vg(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)
+    REAL dt(nlevel, ntime), ht(nlevel, ntime), vt(nlevel, ntime)
+    REAL dtrad(nlevel, ntime)
+    REAL dq(nlevel, ntime), hq(nlevel, ntime), vq(nlevel, ntime)
+    REAL dth(nlevel, ntime), hth(nlevel, ntime), vth(nlevel, ntime), hthl(nlevel, ntime)
+    REAL dr(nlevel, ntime), hr(nlevel, ntime), vr(nlevel, ntime)
+    REAL flat(ntime), sens(ntime), ustar(ntime)
+    REAL uw(nlevel, ntime), vw(nlevel, ntime), q1(nlevel, ntime), q2(nlevel, ntime)
+    REAL ts(ntime), tskin(ntime), ps(ntime)
+    REAL orog_cas, albedo_cas, emiss_cas, q_skin_cas, mom_rough, heat_rough, o3_cas, rugos_cas, clay_cas, sand_cas
+    REAL apbp(nlevel + 1), resul(nlevel, ntime), resul1(nlevel), resul2(ntime), resul3
+
+    INTEGER nid, ierr, ierr1, ierr2, rid, i, int_test
     INTEGER nbvar3d
-    parameter(nbvar3d=78)
-    INTEGER var3didin(nbvar3d),missing_var(nbvar3d)
+    parameter(nbvar3d = 78)
+    INTEGER var3didin(nbvar3d), missing_var(nbvar3d)
     CHARACTER*13 name_var(1:nbvar3d)
 
@@ -385,23 +380,21 @@
     !     &'o3','rugos','clay','sand'/
 
-
-
     data name_var/ &
-                                ! coordonnees pression (n+1 niveaux) #4
-         'coor_par_a','coor_par_b','zf','pressure_h',& ! #1-#4
-                                ! coordonnees pression (n niveaux) #8
-         'ta','qv','ql','qi','ua','va','tke','pa',& ! #5-#12
-                                ! coordonnees pression + temps #46
-         'wa','wap','ug','vg','tnua_adv','tnua_advh','tnua_advv','tnva_adv','tnva_advh','tnva_advv','tnta_adv','tnta_advh',& !  #13 - #25
-         'tnta_advv','tnqv_adv','tnqv_advh','tnqv_advv','thadv','thadvh','thadvv','thladvh',                             & ! #26 - #32
-         'radv','radvh','radvv','tnta_rad','q1','q2','ustress','vstress',                           & ! #33 - #40
-         'rh','ta_nud','qv_nud','ua_nud','va_nud',                                       & ! #41-45
-         'zh_forc','pa_forc','tat','thetat','thetavt','thetalt','qvt','qlt','qit','rvt','uat','vat',   & ! #46-57
-         'nudging_constant_ta', 'nudging_constant_qv', 'nudging_constant_ua', 'nudging_constant_va',           & ! # 58-61
-                                ! coordonnees temps #12
-         'tkes','hfss','hfls','ts_forc','tskin','ps_forc','ustar', &                     ! 62-68
-                                  ! scalaires 
-         'orog','albedo','emiss','q_skin','z0','z0h',       &                    ! 69-74
-         'O3','rugos','clay','sand'/                                                      ! 75-78
+            ! coordonnees pression (n+1 niveaux) #4
+            'coor_par_a', 'coor_par_b', 'zf', 'pressure_h', & ! #1-#4
+            ! coordonnees pression (n niveaux) #8
+            'ta', 'qv', 'ql', 'qi', 'ua', 'va', 'tke', 'pa', & ! #5-#12
+            ! coordonnees pression + temps #46
+            'wa', 'wap', 'ug', 'vg', 'tnua_adv', 'tnua_advh', 'tnua_advv', 'tnva_adv', 'tnva_advh', 'tnva_advv', 'tnta_adv', 'tnta_advh', & !  #13 - #25
+            'tnta_advv', 'tnqv_adv', 'tnqv_advh', 'tnqv_advv', 'thadv', 'thadvh', 'thadvv', 'thladvh', & ! #26 - #32
+            'radv', 'radvh', 'radvv', 'tnta_rad', 'q1', 'q2', 'ustress', 'vstress', & ! #33 - #40
+            'rh', 'ta_nud', 'qv_nud', 'ua_nud', 'va_nud', & ! #41-45
+            'zh_forc', 'pa_forc', 'tat', 'thetat', 'thetavt', 'thetalt', 'qvt', 'qlt', 'qit', 'rvt', 'uat', 'vat', & ! #46-57
+            'nudging_constant_ta', 'nudging_constant_qv', 'nudging_constant_ua', 'nudging_constant_va', & ! # 58-61
+            ! coordonnees temps #12
+            'tkes', 'hfss', 'hfls', 'ts_forc', 'tskin', 'ps_forc', 'ustar', &                     ! 62-68
+            ! scalaires
+            'orog', 'albedo', 'emiss', 'q_skin', 'z0', 'z0h', &                    ! 69-74
+            'O3', 'rugos', 'clay', 'sand'/                                                      ! 75-78
 
 
@@ -411,13 +404,12 @@
     !-----------------------------------------------------------------------
 
-
-    ierr=nf90_inq_varid(nid,'ta',int_test)
+    ierr = nf90_inq_varid(nid, 'ta', int_test)
     IF(ierr/=nf90_noerr) THEN
-       PRINT*, '++++++++++++++++++++++++++++++'
-       PRINT*, 'variable ta missing in cas.nc '
-       PRINT*, 'You are probably using an obsolete version of the 1D cases'
-       PRINT*, 'please dowload the last version of the 1D archive from https://lmdz.lmd.jussieu.fr/pub/'
-       PRINT*, '++++++++++++++++++++++++++++++'
-       CALL abort_gcm ('mod_1D_cases_read_std','bad version of 1D directory',0)
+      PRINT*, '++++++++++++++++++++++++++++++'
+      PRINT*, 'variable ta missing in cas.nc '
+      PRINT*, 'You are probably using an obsolete version of the 1D cases'
+      PRINT*, 'please dowload the last version of the 1D archive from https://lmdz.lmd.jussieu.fr/pub/'
+      PRINT*, '++++++++++++++++++++++++++++++'
+      CALL abort_gcm ('mod_1D_cases_read_std', 'bad version of 1D directory', 0)
     endif
 
@@ -427,196 +419,195 @@
     !-----------------------------------------------------------------------
 
-    do i=1,nbvar3d
-       missing_var(i)=0.
-       ierr=nf90_inq_varid(nid,name_var(i),var3didin(i))
-       PRINT*, 'name_var(i)', name_var(i), var3didin(i)
-       IF(ierr/=nf90_noerr) THEN
-          print *,'Variable manquante dans cas.nc:',i,name_var(i)
-          ierr=nf90_noerr
-          missing_var(i)=1
-       else
+    do i = 1, nbvar3d
+      missing_var(i) = 0.
+      ierr = nf90_inq_varid(nid, name_var(i), var3didin(i))
+      PRINT*, 'name_var(i)', name_var(i), var3didin(i)
+      IF(ierr/=nf90_noerr) THEN
+        print *, 'Variable manquante dans cas.nc:', i, name_var(i)
+        ierr = nf90_noerr
+        missing_var(i) = 1
+      else
+
+        !-----------------------------------------------------------------------
+        ! Activating keys depending on the presence of specific variables in cas.nc
+        !-----------------------------------------------------------------------
+        IF (1 == 1) THEN
+          ! 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_nud' .AND. nint(nudging_qv)==0) stop 'Nudging inconsistency qv'
+          IF (name_var(i) == 'ua_nud' .AND. nint(nudging_u)==0) stop 'Nudging inconsistency u'
+          IF (name_var(i) == 'va_nud' .AND. nint(nudging_v)==0) stop 'Nudging inconsistency v'
+        ELSE
+          PRINT*, 'GUIDAGE : CONSISTENCY CHECK DEACTIVATED FOR TESTS of SANDU/REF'
+        ENDIF
+
+        !-----------------------------------------------------------------------
+        ! Reading variables 1D (N+1) vertical variables (nlevelp1,lat,lon)
+        !-----------------------------------------------------------------------
+        IF(i<=4) THEN
+          ierr = nf90_get_var(nid, var3didin(i), apbp)
+          print *, 'read_SCM(apbp), on a lu ', i, name_var(i)
+          IF(ierr/=nf90_noerr) THEN
+            print *, 'B Pb a la lecture de cas.nc: ', name_var(i)
+            stop "getvarup"
+          endif
 
           !-----------------------------------------------------------------------
-          ! Activating keys depending on the presence of specific variables in cas.nc
+          !  Reading 1D (N) vertical varialbes    (nlevel,lat,lon)
           !-----------------------------------------------------------------------
-          IF ( 1 == 1 ) THEN
-             ! 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_nud' .AND. nint(nudging_qv)==0) stop 'Nudging inconsistency qv'
-             IF ( name_var(i) == 'ua_nud' .AND. nint(nudging_u)==0) stop 'Nudging inconsistency u'
-             IF ( name_var(i) == 'va_nud' .AND. nint(nudging_v)==0) stop 'Nudging inconsistency v'
-          ELSE
-             PRINT*,'GUIDAGE : CONSISTENCY CHECK DEACTIVATED FOR TESTS of SANDU/REF'
-          ENDIF
+        else IF(i>4.AND.i<=12) THEN
+          ierr = nf90_get_var(nid, var3didin(i), resul1)
+          print *, 'read_SCM(resul1), on a lu ', i, name_var(i)
+          IF(ierr/=nf90_noerr) THEN
+            print *, 'C Pb a la lecture de cas.nc: ', name_var(i)
+            stop "getvarup"
+          endif
+          PRINT*, 'Lecture de la variable #i ', i, name_var(i), minval(resul1), maxval(resul1)
 
           !-----------------------------------------------------------------------
-          ! Reading variables 1D (N+1) vertical variables (nlevelp1,lat,lon)
+          !  Reading 2D tim-vertical variables  (time,nlevel,lat,lon)
+          !  TBD : seems to be the same as above.
           !-----------------------------------------------------------------------
-          IF(i<=4) THEN
-             ierr = nf90_get_var(nid,var3didin(i),apbp)
-             print *,'read_SCM(apbp), on a lu ',i,name_var(i)
-             IF(ierr/=nf90_noerr) THEN
-                print *,'B Pb a la lecture de cas.nc: ',name_var(i)
-                stop "getvarup"
-             endif
-
-             !-----------------------------------------------------------------------
-             !  Reading 1D (N) vertical varialbes    (nlevel,lat,lon)   
-             !-----------------------------------------------------------------------
-          else IF(i>4.AND.i<=12) THEN
-             ierr = nf90_get_var(nid,var3didin(i),resul1)
-             print *,'read_SCM(resul1), on a lu ',i,name_var(i)
-             IF(ierr/=nf90_noerr) THEN
-                print *,'C Pb a la lecture de cas.nc: ',name_var(i)
-                stop "getvarup"
-             endif
-             PRINT*,'Lecture de la variable #i ',i,name_var(i),minval(resul1),maxval(resul1)
-
-             !-----------------------------------------------------------------------
-             !  Reading 2D tim-vertical variables  (time,nlevel,lat,lon)
-             !  TBD : seems to be the same as above.
-             !-----------------------------------------------------------------------
-          else IF(i>12.AND.i<=61) THEN
-             ierr = nf90_get_var(nid,var3didin(i),resul)
-             print *,'read_SCM(resul), on a lu ',i,name_var(i)
-             IF(ierr/=nf90_noerr) THEN
-                print *,'D Pb a la lecture de cas.nc: ',name_var(i)
-                stop "getvarup"
-             endif
-             PRINT*,'Lecture de la variable #i ',i,name_var(i),minval(resul),maxval(resul)
-
-             !-----------------------------------------------------------------------
-             !  Reading 1D time variables (time,lat,lon)
-             !-----------------------------------------------------------------------
-          ELSE IF (i>62.AND.i<=75) THEN
-             ierr = nf90_get_var(nid,var3didin(i),resul2)
-             print *,'read_SCM(resul2), on a lu ',i,name_var(i)
-             IF(ierr/=nf90_noerr) THEN
-                print *,'E Pb a la lecture de cas.nc: ',name_var(i)
-                stop "getvarup"
-             endif
-             PRINT*,'Lecture de la variable #i  ',i,name_var(i),minval(resul2),maxval(resul2)
-
-             !-----------------------------------------------------------------------
-             ! Reading scalar variables (lat,lon)
-             !-----------------------------------------------------------------------
-          else
-             ierr = nf90_get_var(nid,var3didin(i),resul3)
-             print *,'read_SCM(resul3), on a lu ',i,name_var(i)
-             IF(ierr/=nf90_noerr) THEN
-                print *,'F Pb a la lecture de cas.nc: ',name_var(i)
-                stop "getvarup"
-             endif
-             PRINT*,'Lecture de la variable #i ',i,name_var(i),resul3
+        else IF(i>12.AND.i<=61) THEN
+          ierr = nf90_get_var(nid, var3didin(i), resul)
+          print *, 'read_SCM(resul), on a lu ', i, name_var(i)
+          IF(ierr/=nf90_noerr) THEN
+            print *, 'D Pb a la lecture de cas.nc: ', name_var(i)
+            stop "getvarup"
           endif
-       endif
-
-       !-----------------------------------------------------------------------
-       ! Attributing variables
-       !-----------------------------------------------------------------------
-       select case(i)
-          !case(1) ; ap=apbp       ! donnees indexees en nlevel+1
-          ! case(2) ; bp=apbp
-       case(3) ; zzh=apbp
-       case(4) ; pph=apbp
-       case(5) ; temp0=resul1    ! donnees initiales
-       case(6) ; qv0=resul1
-       case(7) ; ql0=resul1
-       case(8) ; qi0=resul1
-       case(9) ; u0=resul1
-       case(10) ; v0=resul1
-       case(11) ; tke0=resul1
-       case(12) ; pp0=resul1
-       case(13) ; vitw=resul    ! donnees indexees en nlevel,time
-       case(14) ; omega=resul
-       case(15) ; ug=resul
-       case(16) ; vg=resul
-       case(17) ; du=resul
-       case(18) ; hu=resul
-       case(19) ; vu=resul
-       case(20) ; dv=resul
-       case(21) ; hv=resul
-       case(22) ; vv=resul
-       case(23) ; dt=resul
-       case(24) ; ht=resul
-       case(25) ; vt=resul
-       case(26) ; dq=resul
-       case(27) ; hq=resul
-       case(28) ; vq=resul
-       case(29) ; dth=resul
-       case(30) ; hth=resul
-       case(31) ; vth=resul
-       case(32) ; hthl=resul
-       case(33) ; dr=resul
-       case(34) ; hr=resul
-       case(35) ; vr=resul
-       case(36) ; dtrad=resul
-       case(37) ; q1=resul
-       case(38) ; q2=resul
-       case(39) ; uw=resul
-       case(40) ; vw=resul
-       case(41) ; rh=resul
-       case(42) ; temp_nudg=resul
-       case(43) ; qv_nudg=resul
-       case(44) ; u_nudg=resul
-       case(45) ; v_nudg=resul
-       case(46) ; zz=resul      ! donnees en time,nlevel pour profil initial
-       case(47) ; pp=resul
-       case(48) ; temp=resul
-       case(49) ; theta=resul
-       case(50) ; thv=resul
-       case(51) ; thl=resul
-       case(52) ; qv=resul
-       case(53) ; ql=resul
-       case(54) ; qi=resul
-       case(55) ; rv=resul
-       case(56) ; u=resul
-       case(57) ; v=resul
-       case(58) ; invtau_temp_nudg=resul
-       case(59) ; invtau_qv_nudg=resul
-       case(60) ; invtau_u_nudg=resul
-       case(61) ; invtau_v_nudg=resul
-       case(62) ; tkes=resul2   ! donnees indexees en time
-       case(63) ; sens=resul2
-       case(64) ; flat=resul2
-       case(65) ; ts=resul2
-       case(66) ; tskin=resul2        
-       case(67) ; ps=resul2
-       case(68) ; ustar=resul2
-       case(69) ; orog_cas=resul3      ! constantes
-       case(70) ; albedo_cas=resul3
-       case(71) ; emiss_cas=resul3
-       case(72) ; q_skin_cas=resul3
-       case(73) ; mom_rough=resul3
-       case(74) ; heat_rough=resul3
-       case(75) ; o3_cas=resul3        
-       case(76) ; rugos_cas=resul3
-       case(77) ; clay_cas=resul3
-       case(78) ; sand_cas=resul3
-       end select
-       resul=0.
-       resul1=0.
-       resul2=0.
-       resul3=0.
+          PRINT*, 'Lecture de la variable #i ', i, name_var(i), minval(resul), maxval(resul)
+
+          !-----------------------------------------------------------------------
+          !  Reading 1D time variables (time,lat,lon)
+          !-----------------------------------------------------------------------
+        ELSE IF (i>62.AND.i<=75) THEN
+          ierr = nf90_get_var(nid, var3didin(i), resul2)
+          print *, 'read_SCM(resul2), on a lu ', i, name_var(i)
+          IF(ierr/=nf90_noerr) THEN
+            print *, 'E Pb a la lecture de cas.nc: ', name_var(i)
+            stop "getvarup"
+          endif
+          PRINT*, 'Lecture de la variable #i  ', i, name_var(i), minval(resul2), maxval(resul2)
+
+          !-----------------------------------------------------------------------
+          ! Reading scalar variables (lat,lon)
+          !-----------------------------------------------------------------------
+        else
+          ierr = nf90_get_var(nid, var3didin(i), resul3)
+          print *, 'read_SCM(resul3), on a lu ', i, name_var(i)
+          IF(ierr/=nf90_noerr) THEN
+            print *, 'F Pb a la lecture de cas.nc: ', name_var(i)
+            stop "getvarup"
+          endif
+          PRINT*, 'Lecture de la variable #i ', i, name_var(i), resul3
+        endif
+      endif
+
+      !-----------------------------------------------------------------------
+      ! Attributing variables
+      !-----------------------------------------------------------------------
+      select case(i)
+        !case(1) ; ap=apbp       ! donnees indexees en nlevel+1
+        ! case(2) ; bp=apbp
+      case(3) ; zzh = apbp
+      case(4) ; pph = apbp
+      case(5) ; temp0 = resul1    ! donnees initiales
+      case(6) ; qv0 = resul1
+      case(7) ; ql0 = resul1
+      case(8) ; qi0 = resul1
+      case(9) ; u0 = resul1
+      case(10) ; v0 = resul1
+      case(11) ; tke0 = resul1
+      case(12) ; pp0 = resul1
+      case(13) ; vitw = resul    ! donnees indexees en nlevel,time
+      case(14) ; omega = resul
+      case(15) ; ug = resul
+      case(16) ; vg = resul
+      case(17) ; du = resul
+      case(18) ; hu = resul
+      case(19) ; vu = resul
+      case(20) ; dv = resul
+      case(21) ; hv = resul
+      case(22) ; vv = resul
+      case(23) ; dt = resul
+      case(24) ; ht = resul
+      case(25) ; vt = resul
+      case(26) ; dq = resul
+      case(27) ; hq = resul
+      case(28) ; vq = resul
+      case(29) ; dth = resul
+      case(30) ; hth = resul
+      case(31) ; vth = resul
+      case(32) ; hthl = resul
+      case(33) ; dr = resul
+      case(34) ; hr = resul
+      case(35) ; vr = resul
+      case(36) ; dtrad = resul
+      case(37) ; q1 = resul
+      case(38) ; q2 = resul
+      case(39) ; uw = resul
+      case(40) ; vw = resul
+      case(41) ; rh = resul
+      case(42) ; temp_nudg = resul
+      case(43) ; qv_nudg = resul
+      case(44) ; u_nudg = resul
+      case(45) ; v_nudg = resul
+      case(46) ; zz = resul      ! donnees en time,nlevel pour profil initial
+      case(47) ; pp = resul
+      case(48) ; temp = resul
+      case(49) ; theta = resul
+      case(50) ; thv = resul
+      case(51) ; thl = resul
+      case(52) ; qv = resul
+      case(53) ; ql = resul
+      case(54) ; qi = resul
+      case(55) ; rv = resul
+      case(56) ; u = resul
+      case(57) ; v = resul
+      case(58) ; invtau_temp_nudg = resul
+      case(59) ; invtau_qv_nudg = resul
+      case(60) ; invtau_u_nudg = resul
+      case(61) ; invtau_v_nudg = resul
+      case(62) ; tkes = resul2   ! donnees indexees en time
+      case(63) ; sens = resul2
+      case(64) ; flat = resul2
+      case(65) ; ts = resul2
+      case(66) ; tskin = resul2
+      case(67) ; ps = resul2
+      case(68) ; ustar = resul2
+      case(69) ; orog_cas = resul3      ! constantes
+      case(70) ; albedo_cas = resul3
+      case(71) ; emiss_cas = resul3
+      case(72) ; q_skin_cas = resul3
+      case(73) ; mom_rough = resul3
+      case(74) ; heat_rough = resul3
+      case(75) ; o3_cas = resul3
+      case(76) ; rugos_cas = resul3
+      case(77) ; clay_cas = resul3
+      case(78) ; sand_cas = resul3
+      end select
+      resul = 0.
+      resul1 = 0.
+      resul2 = 0.
+      resul3 = 0.
     enddo
-    PRINT*,'Lecture de la variable APRES ,sens ',minval(sens),maxval(sens)
-    PRINT*,'Lecture de la variable APRES ,flat ',minval(flat),maxval(flat)
+    PRINT*, 'Lecture de la variable APRES ,sens ', minval(sens), maxval(sens)
+    PRINT*, 'Lecture de la variable APRES ,flat ', minval(flat), maxval(flat)
 
     !CR:ATTENTION EN ATTENTE DE REGLER LA QUESTION DU PAS DE TEMPS INITIAL
-    do t=1,ntime
-       do k=1,nlevel
-          temp(k,t)=temp0(k)
-          qv(k,t)=qv0(k)
-          ql(k,t)=ql0(k)
-          qi(k,t)=qi0(k)
-          u(k,t)=u0(k)
-          v(k,t)=v0(k)
-          tke(k,t)=tke0(k)
-       enddo
+    do t = 1, ntime
+      do k = 1, nlevel
+        temp(k, t) = temp0(k)
+        qv(k, t) = qv0(k)
+        ql(k, t) = ql0(k)
+        qi(k, t) = qi0(k)
+        u(k, t) = u0(k)
+        v(k, t) = v0(k)
+        tke(k, t) = tke0(k)
+      enddo
     enddo
-!!!! TRAVAIL : EN FONCTION DES DECISIONS SUR LA SPECIFICATION DE W
-!!!omega=-vitw*pres*rg/(rd*temp)
+    !!!! TRAVAIL : EN FONCTION DES DECISIONS SUR LA SPECIFICATION DE W
+    !!!omega=-vitw*pres*rg/(rd*temp)
     !-----------------------------------------------------------------------
-
 
   END SUBROUTINE read_SCM
@@ -628,33 +619,28 @@
 
   !**********************************************************************************************
-  SUBROUTINE interp_case_time_std(day,day1,annee_ref                           &
-       !    &         ,year_cas,day_cas,nt_cas,pdt_forc,nlev_cas                         &
-       ,nt_cas,nlev_cas                                                   &
-       ,ts_cas,tskin_cas,ps_cas,plev_cas,t_cas,theta_cas,thv_cas,thl_cas            &
-       ,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     &
-       ,invtau_temp_nudg_cas,invtau_qv_nudg_cas,invtau_u_nudg_cas,invtau_v_nudg_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,tkes_cas                               &
-
-       ,ts_prof_cas,tskin_prof_cas,ps_prof_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     &
-       ,invtau_temp_nudg_prof_cas,invtau_qv_nudg_prof_cas,invtau_u_nudg_prof_cas,invtau_v_nudg_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,tkes_prof_cas)
-
-
-
-
-
+  SUBROUTINE interp_case_time_std(day, day1, annee_ref                           &
+          !    &         ,year_cas,day_cas,nt_cas,pdt_forc,nlev_cas                         &
+          , nt_cas, nlev_cas                                                   &
+          , ts_cas, tskin_cas, ps_cas, plev_cas, t_cas, theta_cas, thv_cas, thl_cas            &
+          , 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     &
+          , invtau_temp_nudg_cas, invtau_qv_nudg_cas, invtau_u_nudg_cas, invtau_v_nudg_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, tkes_cas                               &
+
+          , ts_prof_cas, tskin_prof_cas, ps_prof_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     &
+          , invtau_temp_nudg_prof_cas, invtau_qv_nudg_prof_cas, invtau_u_nudg_prof_cas, invtau_v_nudg_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, tkes_prof_cas)
 
     IMPLICIT NONE
@@ -665,5 +651,5 @@
     ! day: current julian day (e.g. 717538.2)
     ! day1: first day of the simulation
-    ! nt_cas: total nb of data in the forcing 
+    ! nt_cas: total nb of data in the forcing
     ! pdt_cas: total time interval (in sec) between 2 forcing data
     !---------------------------------------------------------------------------------------
@@ -674,57 +660,56 @@
     ! inputs:
     INTEGER annee_ref
-    INTEGER nt_cas,nlev_cas
-    REAL day, day1,day_cas
-    REAL ts_cas(nt_cas),tskin_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)
-    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)
-    REAL ug_cas(nlev_cas,nt_cas),vg_cas(nlev_cas,nt_cas)
-    REAL temp_nudg_cas(nlev_cas,nt_cas),qv_nudg_cas(nlev_cas,nt_cas)
-    REAL u_nudg_cas(nlev_cas,nt_cas),v_nudg_cas(nlev_cas,nt_cas)
-
-    REAL invtau_temp_nudg_cas(nlev_cas,nt_cas),invtau_qv_nudg_cas(nlev_cas,nt_cas)
-    REAL invtau_u_nudg_cas(nlev_cas,nt_cas),invtau_v_nudg_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)
-    REAL dt_cas(nlev_cas,nt_cas),ht_cas(nlev_cas,nt_cas),vt_cas(nlev_cas,nt_cas)
-    REAL dth_cas(nlev_cas,nt_cas),hth_cas(nlev_cas,nt_cas),vth_cas(nlev_cas,nt_cas)
-    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),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)
+    INTEGER nt_cas, nlev_cas
+    REAL day, day1, day_cas
+    REAL ts_cas(nt_cas), tskin_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)
+    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)
+    REAL ug_cas(nlev_cas, nt_cas), vg_cas(nlev_cas, nt_cas)
+    REAL temp_nudg_cas(nlev_cas, nt_cas), qv_nudg_cas(nlev_cas, nt_cas)
+    REAL u_nudg_cas(nlev_cas, nt_cas), v_nudg_cas(nlev_cas, nt_cas)
+
+    REAL invtau_temp_nudg_cas(nlev_cas, nt_cas), invtau_qv_nudg_cas(nlev_cas, nt_cas)
+    REAL invtau_u_nudg_cas(nlev_cas, nt_cas), invtau_v_nudg_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)
+    REAL dt_cas(nlev_cas, nt_cas), ht_cas(nlev_cas, nt_cas), vt_cas(nlev_cas, nt_cas)
+    REAL dth_cas(nlev_cas, nt_cas), hth_cas(nlev_cas, nt_cas), vth_cas(nlev_cas, nt_cas)
+    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), 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)
 
     ! outputs:
     REAL plev_prof_cas(nlev_cas)
-    REAL t_prof_cas(nlev_cas),theta_prof_cas(nlev_cas),thl_prof_cas(nlev_cas),thv_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)
-    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)
-
-    REAL invtau_temp_nudg_prof_cas(nlev_cas),invtau_qv_nudg_prof_cas(nlev_cas)
-    REAL invtau_u_nudg_prof_cas(nlev_cas),invtau_v_nudg_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)
-    REAL dt_prof_cas(nlev_cas),ht_prof_cas(nlev_cas),vt_prof_cas(nlev_cas)
-    REAL dth_prof_cas(nlev_cas),hth_prof_cas(nlev_cas),vth_prof_cas(nlev_cas)
+    REAL t_prof_cas(nlev_cas), theta_prof_cas(nlev_cas), thl_prof_cas(nlev_cas), thv_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)
+    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)
+
+    REAL invtau_temp_nudg_prof_cas(nlev_cas), invtau_qv_nudg_prof_cas(nlev_cas)
+    REAL invtau_u_nudg_prof_cas(nlev_cas), invtau_v_nudg_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)
+    REAL dt_prof_cas(nlev_cas), ht_prof_cas(nlev_cas), vt_prof_cas(nlev_cas)
+    REAL dth_prof_cas(nlev_cas), hth_prof_cas(nlev_cas), vth_prof_cas(nlev_cas)
     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,tkes_prof_cas,ts_prof_cas,tskin_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)
+    REAL dq_prof_cas(nlev_cas), hq_prof_cas(nlev_cas), vq_prof_cas(nlev_cas)
+    REAL lat_prof_cas, sens_prof_cas, tkes_prof_cas, ts_prof_cas, tskin_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:
-    INTEGER it_cas1, it_cas2,k
-    REAL timeit,time_cas1,time_cas2,frac
-
-
-    PRINT*,'Check time',day1,day_ju_ini_cas,day_deb+1,pdt_cas
+    INTEGER it_cas1, it_cas2, k
+    REAL timeit, time_cas1, time_cas2, frac
+
+    PRINT*, 'Check time', day1, day_ju_ini_cas, day_deb + 1, pdt_cas
     !       do k=1,nlev_cas
     !       PRINT*,'debut de interp2_case_time, plev_cas=',k,plev_cas(k,1)
@@ -761,10 +746,10 @@
     !        timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
     !       endif
-    timeit=(day-day_ju_ini_cas)*86400
-    print *,'day=',day
-    print *,'day_ju_ini_cas=',day_ju_ini_cas
-    print *,'pdt_cas=',pdt_cas
-    print *,'timeit=',timeit
-    print *,'nt_cas=',nt_cas
+    timeit = (day - day_ju_ini_cas) * 86400
+    print *, 'day=', day
+    print *, 'day_ju_ini_cas=', day_ju_ini_cas
+    print *, 'pdt_cas=', pdt_cas
+    print *, 'timeit=', timeit
+    print *, 'nt_cas=', nt_cas
 
     ! Determine the closest observation times:
@@ -774,134 +759,133 @@
     !       time_cas2=(it_cas2-1)*pdt_cas
 
-    it_cas1=INT(timeit/pdt_cas)+1
+    it_cas1 = INT(timeit / pdt_cas) + 1
     IF (it_cas1 == nt_cas) THEN
-       it_cas2=it_cas1 
+      it_cas2 = it_cas1
     ELSE
-       it_cas2=it_cas1 + 1
+      it_cas2 = it_cas1 + 1
     ENDIF
-    time_cas1=(it_cas1-1)*pdt_cas
-    time_cas2=(it_cas2-1)*pdt_cas
+    time_cas1 = (it_cas1 - 1) * pdt_cas
+    time_cas2 = (it_cas2 - 1) * pdt_cas
     !     print *,'timeit,pdt_cas,nt_cas=',timeit,pdt_cas,nt_cas
     !     print *,'it_cas1,it_cas2,time_cas1,time_cas2=',it_cas1,it_cas2,time_cas1,time_cas2
 
     IF (it_cas1 > nt_cas) THEN
-       WRITE(*,*) 'PB-stop: day, day_ju_ini_cas,it_cas1, it_cas2, timeit: '            &
-            ,day,day_ju_ini_cas,it_cas1,it_cas2,timeit
-       stop
+      WRITE(*, *) 'PB-stop: day, day_ju_ini_cas,it_cas1, it_cas2, timeit: '            &
+              , day, day_ju_ini_cas, it_cas1, it_cas2, timeit
+      stop
     endif
 
     ! time interpolation:
     IF (it_cas1 == it_cas2) THEN
-       frac=0.
+      frac = 0.
     ELSE
-       frac=(time_cas2-timeit)/(time_cas2-time_cas1)
-       frac=max(frac,0.0)
+      frac = (time_cas2 - timeit) / (time_cas2 - time_cas1)
+      frac = max(frac, 0.0)
     ENDIF
 
     lat_prof_cas = lat_cas(it_cas2)                                   &
-         -frac*(lat_cas(it_cas2)-lat_cas(it_cas1)) 
+            - frac * (lat_cas(it_cas2) - lat_cas(it_cas1))
     sens_prof_cas = sens_cas(it_cas2)                                 &
-         -frac*(sens_cas(it_cas2)-sens_cas(it_cas1))
+            - frac * (sens_cas(it_cas2) - sens_cas(it_cas1))
     tkes_prof_cas = tkes_cas(it_cas2)                                   &
-         -frac*(tkes_cas(it_cas2)-tkes_cas(it_cas1))
+            - 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))
+            - frac * (ts_cas(it_cas2) - ts_cas(it_cas1))
     tskin_prof_cas = tskin_cas(it_cas2)                                     &
-         -frac*(tskin_cas(it_cas2)-tskin_cas(it_cas1))
+            - frac * (tskin_cas(it_cas2) - tskin_cas(it_cas1))
     ps_prof_cas = ps_cas(it_cas2)                                     &
-         -frac*(ps_cas(it_cas2)-ps_cas(it_cas1))
+            - frac * (ps_cas(it_cas2) - ps_cas(it_cas1))
     ustar_prof_cas = ustar_cas(it_cas2)                               &
-         -frac*(ustar_cas(it_cas2)-ustar_cas(it_cas1))
-
-    do k=1,nlev_cas
-       plev_prof_cas(k) = plev_cas(k,it_cas2)                           &     
-            -frac*(plev_cas(k,it_cas2)-plev_cas(k,it_cas1))
-       t_prof_cas(k) = t_cas(k,it_cas2)                                 &        
-            -frac*(t_cas(k,it_cas2)-t_cas(k,it_cas1))
-       !print *,'k,frac,plev_cas1,plev_cas2=',k,frac,plev_cas(k,it_cas1),plev_cas(k,it_cas2)
-       theta_prof_cas(k) = theta_cas(k,it_cas2)                         &                      
-            -frac*(theta_cas(k,it_cas2)-theta_cas(k,it_cas1))
-       thv_prof_cas(k) = thv_cas(k,it_cas2)                             &          
-            -frac*(thv_cas(k,it_cas2)-thv_cas(k,it_cas1))
-       thl_prof_cas(k) = thl_cas(k,it_cas2)                             &              
-            -frac*(thl_cas(k,it_cas2)-thl_cas(k,it_cas1))
-       qv_prof_cas(k) = qv_cas(k,it_cas2)                               &
-            -frac*(qv_cas(k,it_cas2)-qv_cas(k,it_cas1))
-       ql_prof_cas(k) = ql_cas(k,it_cas2)                               &
-            -frac*(ql_cas(k,it_cas2)-ql_cas(k,it_cas1))
-       qi_prof_cas(k) = qi_cas(k,it_cas2)                               &
-            -frac*(qi_cas(k,it_cas2)-qi_cas(k,it_cas1))
-       u_prof_cas(k) = u_cas(k,it_cas2)                                 &
-            -frac*(u_cas(k,it_cas2)-u_cas(k,it_cas1))
-       v_prof_cas(k) = v_cas(k,it_cas2)                                 &
-            -frac*(v_cas(k,it_cas2)-v_cas(k,it_cas1))
-       ug_prof_cas(k) = ug_cas(k,it_cas2)                               &
-            -frac*(ug_cas(k,it_cas2)-ug_cas(k,it_cas1))
-       vg_prof_cas(k) = vg_cas(k,it_cas2)                               &
-            -frac*(vg_cas(k,it_cas2)-vg_cas(k,it_cas1))
-       temp_nudg_prof_cas(k) = temp_nudg_cas(k,it_cas2)                    &
-            -frac*(temp_nudg_cas(k,it_cas2)-temp_nudg_cas(k,it_cas1))
-       qv_nudg_prof_cas(k) = qv_nudg_cas(k,it_cas2)                        &
-            -frac*(qv_nudg_cas(k,it_cas2)-qv_nudg_cas(k,it_cas1))
-       u_nudg_prof_cas(k) = u_nudg_cas(k,it_cas2)                          &
-            -frac*(u_nudg_cas(k,it_cas2)-u_nudg_cas(k,it_cas1))
-       v_nudg_prof_cas(k) = v_nudg_cas(k,it_cas2)                          &
-            -frac*(v_nudg_cas(k,it_cas2)-v_nudg_cas(k,it_cas1))
-       invtau_temp_nudg_prof_cas(k) = invtau_temp_nudg_cas(k,it_cas2)                    &
-            -frac*(invtau_temp_nudg_cas(k,it_cas2)-invtau_temp_nudg_cas(k,it_cas1))
-       invtau_qv_nudg_prof_cas(k) = invtau_qv_nudg_cas(k,it_cas2)                        &
-            -frac*(invtau_qv_nudg_cas(k,it_cas2)-invtau_qv_nudg_cas(k,it_cas1))
-       invtau_u_nudg_prof_cas(k) = invtau_u_nudg_cas(k,it_cas2)                          &
-            -frac*(invtau_u_nudg_cas(k,it_cas2)-invtau_u_nudg_cas(k,it_cas1))
-       invtau_v_nudg_prof_cas(k) = invtau_v_nudg_cas(k,it_cas2)                          &
-            -frac*(invtau_v_nudg_cas(k,it_cas2)-invtau_v_nudg_cas(k,it_cas1))
-       vitw_prof_cas(k) = vitw_cas(k,it_cas2)                           &
-            -frac*(vitw_cas(k,it_cas2)-vitw_cas(k,it_cas1))
-       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))
-       hu_prof_cas(k) = hu_cas(k,it_cas2)                               &
-            -frac*(hu_cas(k,it_cas2)-hu_cas(k,it_cas1))
-       vu_prof_cas(k) = vu_cas(k,it_cas2)                               &
-            -frac*(vu_cas(k,it_cas2)-vu_cas(k,it_cas1))
-       dv_prof_cas(k) = dv_cas(k,it_cas2)                               &
-            -frac*(dv_cas(k,it_cas2)-dv_cas(k,it_cas1))
-       hv_prof_cas(k) = hv_cas(k,it_cas2)                               &
-            -frac*(hv_cas(k,it_cas2)-hv_cas(k,it_cas1))
-       vv_prof_cas(k) = vv_cas(k,it_cas2)                               &
-            -frac*(vv_cas(k,it_cas2)-vv_cas(k,it_cas1))
-       dt_prof_cas(k) = dt_cas(k,it_cas2)                               &
-            -frac*(dt_cas(k,it_cas2)-dt_cas(k,it_cas1))
-       ht_prof_cas(k) = ht_cas(k,it_cas2)                               &
-            -frac*(ht_cas(k,it_cas2)-ht_cas(k,it_cas1))
-       vt_prof_cas(k) = vt_cas(k,it_cas2)                               &
-            -frac*(vt_cas(k,it_cas2)-vt_cas(k,it_cas1))
-       dth_prof_cas(k) = dth_cas(k,it_cas2)                             &
-            -frac*(dth_cas(k,it_cas2)-dth_cas(k,it_cas1))
-       hth_prof_cas(k) = hth_cas(k,it_cas2)                             &
-            -frac*(hth_cas(k,it_cas2)-hth_cas(k,it_cas1))
-       vth_prof_cas(k) = vth_cas(k,it_cas2)                             &
-            -frac*(vth_cas(k,it_cas2)-vth_cas(k,it_cas1))
-       dtrad_prof_cas(k) = dtrad_cas(k,it_cas2)                         &
-            -frac*(dtrad_cas(k,it_cas2)-dtrad_cas(k,it_cas1))
-       dq_prof_cas(k) = dq_cas(k,it_cas2)                               &
-            -frac*(dq_cas(k,it_cas2)-dq_cas(k,it_cas1))
-       hq_prof_cas(k) = hq_cas(k,it_cas2)                               &
-            -frac*(hq_cas(k,it_cas2)-hq_cas(k,it_cas1))
-       vq_prof_cas(k) = vq_cas(k,it_cas2)                               &
-            -frac*(vq_cas(k,it_cas2)-vq_cas(k,it_cas1))
-       uw_prof_cas(k) = uw_cas(k,it_cas2)                                &
-            -frac*(uw_cas(k,it_cas2)-uw_cas(k,it_cas1))
-       vw_prof_cas(k) = vw_cas(k,it_cas2)                                &
-            -frac*(vw_cas(k,it_cas2)-vw_cas(k,it_cas1))
-       q1_prof_cas(k) = q1_cas(k,it_cas2)                                &
-            -frac*(q1_cas(k,it_cas2)-q1_cas(k,it_cas1))
-       q2_prof_cas(k) = q2_cas(k,it_cas2)                                &
-            -frac*(q2_cas(k,it_cas2)-q2_cas(k,it_cas1))
+            - frac * (ustar_cas(it_cas2) - ustar_cas(it_cas1))
+
+    do k = 1, nlev_cas
+      plev_prof_cas(k) = plev_cas(k, it_cas2)                           &
+              - frac * (plev_cas(k, it_cas2) - plev_cas(k, it_cas1))
+      t_prof_cas(k) = t_cas(k, it_cas2)                                 &
+              - frac * (t_cas(k, it_cas2) - t_cas(k, it_cas1))
+      !print *,'k,frac,plev_cas1,plev_cas2=',k,frac,plev_cas(k,it_cas1),plev_cas(k,it_cas2)
+      theta_prof_cas(k) = theta_cas(k, it_cas2)                         &
+              - frac * (theta_cas(k, it_cas2) - theta_cas(k, it_cas1))
+      thv_prof_cas(k) = thv_cas(k, it_cas2)                             &
+              - frac * (thv_cas(k, it_cas2) - thv_cas(k, it_cas1))
+      thl_prof_cas(k) = thl_cas(k, it_cas2)                             &
+              - frac * (thl_cas(k, it_cas2) - thl_cas(k, it_cas1))
+      qv_prof_cas(k) = qv_cas(k, it_cas2)                               &
+              - frac * (qv_cas(k, it_cas2) - qv_cas(k, it_cas1))
+      ql_prof_cas(k) = ql_cas(k, it_cas2)                               &
+              - frac * (ql_cas(k, it_cas2) - ql_cas(k, it_cas1))
+      qi_prof_cas(k) = qi_cas(k, it_cas2)                               &
+              - frac * (qi_cas(k, it_cas2) - qi_cas(k, it_cas1))
+      u_prof_cas(k) = u_cas(k, it_cas2)                                 &
+              - frac * (u_cas(k, it_cas2) - u_cas(k, it_cas1))
+      v_prof_cas(k) = v_cas(k, it_cas2)                                 &
+              - frac * (v_cas(k, it_cas2) - v_cas(k, it_cas1))
+      ug_prof_cas(k) = ug_cas(k, it_cas2)                               &
+              - frac * (ug_cas(k, it_cas2) - ug_cas(k, it_cas1))
+      vg_prof_cas(k) = vg_cas(k, it_cas2)                               &
+              - frac * (vg_cas(k, it_cas2) - vg_cas(k, it_cas1))
+      temp_nudg_prof_cas(k) = temp_nudg_cas(k, it_cas2)                    &
+              - frac * (temp_nudg_cas(k, it_cas2) - temp_nudg_cas(k, it_cas1))
+      qv_nudg_prof_cas(k) = qv_nudg_cas(k, it_cas2)                        &
+              - frac * (qv_nudg_cas(k, it_cas2) - qv_nudg_cas(k, it_cas1))
+      u_nudg_prof_cas(k) = u_nudg_cas(k, it_cas2)                          &
+              - frac * (u_nudg_cas(k, it_cas2) - u_nudg_cas(k, it_cas1))
+      v_nudg_prof_cas(k) = v_nudg_cas(k, it_cas2)                          &
+              - frac * (v_nudg_cas(k, it_cas2) - v_nudg_cas(k, it_cas1))
+      invtau_temp_nudg_prof_cas(k) = invtau_temp_nudg_cas(k, it_cas2)                    &
+              - frac * (invtau_temp_nudg_cas(k, it_cas2) - invtau_temp_nudg_cas(k, it_cas1))
+      invtau_qv_nudg_prof_cas(k) = invtau_qv_nudg_cas(k, it_cas2)                        &
+              - frac * (invtau_qv_nudg_cas(k, it_cas2) - invtau_qv_nudg_cas(k, it_cas1))
+      invtau_u_nudg_prof_cas(k) = invtau_u_nudg_cas(k, it_cas2)                          &
+              - frac * (invtau_u_nudg_cas(k, it_cas2) - invtau_u_nudg_cas(k, it_cas1))
+      invtau_v_nudg_prof_cas(k) = invtau_v_nudg_cas(k, it_cas2)                          &
+              - frac * (invtau_v_nudg_cas(k, it_cas2) - invtau_v_nudg_cas(k, it_cas1))
+      vitw_prof_cas(k) = vitw_cas(k, it_cas2)                           &
+              - frac * (vitw_cas(k, it_cas2) - vitw_cas(k, it_cas1))
+      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))
+      hu_prof_cas(k) = hu_cas(k, it_cas2)                               &
+              - frac * (hu_cas(k, it_cas2) - hu_cas(k, it_cas1))
+      vu_prof_cas(k) = vu_cas(k, it_cas2)                               &
+              - frac * (vu_cas(k, it_cas2) - vu_cas(k, it_cas1))
+      dv_prof_cas(k) = dv_cas(k, it_cas2)                               &
+              - frac * (dv_cas(k, it_cas2) - dv_cas(k, it_cas1))
+      hv_prof_cas(k) = hv_cas(k, it_cas2)                               &
+              - frac * (hv_cas(k, it_cas2) - hv_cas(k, it_cas1))
+      vv_prof_cas(k) = vv_cas(k, it_cas2)                               &
+              - frac * (vv_cas(k, it_cas2) - vv_cas(k, it_cas1))
+      dt_prof_cas(k) = dt_cas(k, it_cas2)                               &
+              - frac * (dt_cas(k, it_cas2) - dt_cas(k, it_cas1))
+      ht_prof_cas(k) = ht_cas(k, it_cas2)                               &
+              - frac * (ht_cas(k, it_cas2) - ht_cas(k, it_cas1))
+      vt_prof_cas(k) = vt_cas(k, it_cas2)                               &
+              - frac * (vt_cas(k, it_cas2) - vt_cas(k, it_cas1))
+      dth_prof_cas(k) = dth_cas(k, it_cas2)                             &
+              - frac * (dth_cas(k, it_cas2) - dth_cas(k, it_cas1))
+      hth_prof_cas(k) = hth_cas(k, it_cas2)                             &
+              - frac * (hth_cas(k, it_cas2) - hth_cas(k, it_cas1))
+      vth_prof_cas(k) = vth_cas(k, it_cas2)                             &
+              - frac * (vth_cas(k, it_cas2) - vth_cas(k, it_cas1))
+      dtrad_prof_cas(k) = dtrad_cas(k, it_cas2)                         &
+              - frac * (dtrad_cas(k, it_cas2) - dtrad_cas(k, it_cas1))
+      dq_prof_cas(k) = dq_cas(k, it_cas2)                               &
+              - frac * (dq_cas(k, it_cas2) - dq_cas(k, it_cas1))
+      hq_prof_cas(k) = hq_cas(k, it_cas2)                               &
+              - frac * (hq_cas(k, it_cas2) - hq_cas(k, it_cas1))
+      vq_prof_cas(k) = vq_cas(k, it_cas2)                               &
+              - frac * (vq_cas(k, it_cas2) - vq_cas(k, it_cas1))
+      uw_prof_cas(k) = uw_cas(k, it_cas2)                                &
+              - frac * (uw_cas(k, it_cas2) - uw_cas(k, it_cas1))
+      vw_prof_cas(k) = vw_cas(k, it_cas2)                                &
+              - frac * (vw_cas(k, it_cas2) - vw_cas(k, it_cas1))
+      q1_prof_cas(k) = q1_cas(k, it_cas2)                                &
+              - frac * (q1_cas(k, it_cas2) - q1_cas(k, it_cas1))
+      q2_prof_cas(k) = q2_cas(k, it_cas2)                                &
+              - frac * (q2_cas(k, it_cas2) - q2_cas(k, it_cas1))
     enddo
-
 
   END SUBROUTINE interp_case_time_std
@@ -909,28 +893,29 @@
   !**********************************************************************************************
   !=====================================================================
-  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                    &
-       ,invtau_temp_nudg_prof_cas,invtau_qv_nudg_prof_cas,invtau_u_nudg_prof_cas,invtau_v_nudg_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                                                 &
-
-       ,t_mod_cas,theta_mod_cas,thv_mod_cas,thl_mod_cas                                        &
-       ,qv_mod_cas,ql_mod_cas,qi_mod_cas,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                        &
-       ,invtau_temp_nudg_mod_cas,invtau_qv_nudg_mod_cas,invtau_u_nudg_mod_cas,invtau_v_nudg_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        &
-       ,dth_mod_cas,hth_mod_cas,vth_mod_cas,mxcalc)
+  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                    &
+          , invtau_temp_nudg_prof_cas, invtau_qv_nudg_prof_cas, invtau_u_nudg_prof_cas, invtau_v_nudg_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                                                 &
+
+          , t_mod_cas, theta_mod_cas, thv_mod_cas, thl_mod_cas                                        &
+          , qv_mod_cas, ql_mod_cas, qi_mod_cas, 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                        &
+          , invtau_temp_nudg_mod_cas, invtau_qv_nudg_mod_cas, invtau_u_nudg_mod_cas, invtau_v_nudg_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        &
+          , dth_mod_cas, hth_mod_cas, vth_mod_cas, mxcalc)
+
+    USE lmdz_yomcst
 
     IMPLICIT NONE
 
-    INCLUDE "YOMCST.h"
     INCLUDE "dimensions.h"
 
@@ -940,7 +925,7 @@
 
     INTEGER nlevmax
-    parameter (nlevmax=41)
-    INTEGER nlev_cas,mxcalc
-    !       real play(llm), plev_prof(nlevmax) 
+    parameter (nlevmax = 41)
+    INTEGER nlev_cas, mxcalc
+    !       real play(llm), plev_prof(nlevmax)
     !       real t_prof(nlevmax),q_prof(nlevmax)
     !       real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
@@ -948,36 +933,36 @@
     !       real hq_prof(nlevmax),vq_prof(nlevmax)
 
-    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),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)
-    REAL invtau_temp_nudg_prof_cas(nlev_cas),invtau_qv_nudg_prof_cas(nlev_cas)
-    REAL invtau_u_nudg_prof_cas(nlev_cas),invtau_v_nudg_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)
-    REAL dt_prof_cas(nlev_cas),ht_prof_cas(nlev_cas),vt_prof_cas(nlev_cas),dtrad_prof_cas(nlev_cas)
-    REAL dth_prof_cas(nlev_cas),hth_prof_cas(nlev_cas),vth_prof_cas(nlev_cas)
-    REAL dq_prof_cas(nlev_cas),hq_prof_cas(nlev_cas),vq_prof_cas(nlev_cas)
-
-    REAL t_mod_cas(llm),theta_mod_cas(llm),thv_mod_cas(llm),thl_mod_cas(llm)
-    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),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)
-    REAL invtau_temp_nudg_mod_cas(llm),invtau_qv_nudg_mod_cas(llm)
-    REAL invtau_u_nudg_mod_cas(llm),invtau_v_nudg_mod_cas(llm)
-    REAL du_mod_cas(llm),hu_mod_cas(llm),vu_mod_cas(llm)
-    REAL dv_mod_cas(llm),hv_mod_cas(llm),vv_mod_cas(llm)
-    REAL dt_mod_cas(llm),ht_mod_cas(llm),vt_mod_cas(llm),dtrad_mod_cas(llm)
-    REAL dth_mod_cas(llm),hth_mod_cas(llm),vth_mod_cas(llm)
-    REAL dq_mod_cas(llm),hq_mod_cas(llm),vq_mod_cas(llm)
-
-    INTEGER l,k,k1,k2
-    REAL frac,frac1,frac2,fact
+    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), 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)
+    REAL invtau_temp_nudg_prof_cas(nlev_cas), invtau_qv_nudg_prof_cas(nlev_cas)
+    REAL invtau_u_nudg_prof_cas(nlev_cas), invtau_v_nudg_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)
+    REAL dt_prof_cas(nlev_cas), ht_prof_cas(nlev_cas), vt_prof_cas(nlev_cas), dtrad_prof_cas(nlev_cas)
+    REAL dth_prof_cas(nlev_cas), hth_prof_cas(nlev_cas), vth_prof_cas(nlev_cas)
+    REAL dq_prof_cas(nlev_cas), hq_prof_cas(nlev_cas), vq_prof_cas(nlev_cas)
+
+    REAL t_mod_cas(llm), theta_mod_cas(llm), thv_mod_cas(llm), thl_mod_cas(llm)
+    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), 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)
+    REAL invtau_temp_nudg_mod_cas(llm), invtau_qv_nudg_mod_cas(llm)
+    REAL invtau_u_nudg_mod_cas(llm), invtau_v_nudg_mod_cas(llm)
+    REAL du_mod_cas(llm), hu_mod_cas(llm), vu_mod_cas(llm)
+    REAL dv_mod_cas(llm), hv_mod_cas(llm), vv_mod_cas(llm)
+    REAL dt_mod_cas(llm), ht_mod_cas(llm), vt_mod_cas(llm), dtrad_mod_cas(llm)
+    REAL dth_mod_cas(llm), hth_mod_cas(llm), vth_mod_cas(llm)
+    REAL dq_mod_cas(llm), hq_mod_cas(llm), vq_mod_cas(llm)
+
+    INTEGER l, k, k1, k2
+    REAL frac, frac1, frac2, fact
 
 
@@ -987,170 +972,168 @@
     do l = 1, llm
 
-       IF (play(l)>=plev_prof_cas(nlev_cas)) THEN
-          mxcalc=l
-          !        print *,'debut interp2, mxcalc=',mxcalc
-          k1=0
-          k2=0
-
-          IF (play(l)<=plev_prof_cas(1)) THEN
-             do k = 1, nlev_cas-1
-                IF (play(l)<=plev_prof_cas(k).AND. play(l)>plev_prof_cas(k+1)) THEN
-                   k1=k
-                   k2=k+1
-                endif
-             enddo
-
-             IF (k1==0 .OR. k2==0) THEN
-                WRITE(*,*) 'PB! k1, k2 = ',k1,k2
-                WRITE(*,*) 'l,play(l) = ',l,play(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)-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))
-             IF(theta_mod_cas(l)/=0) t_mod_cas(l)= theta_mod_cas(l)*(play(l)/100000.)**(RD/RCPD)
-             thv_mod_cas(l)= thv_prof_cas(k2) - frac*(thv_prof_cas(k2)-thv_prof_cas(k1))
-             thl_mod_cas(l)= thl_prof_cas(k2) - frac*(thl_prof_cas(k2)-thl_prof_cas(k1))
-             qv_mod_cas(l)= qv_prof_cas(k2) - frac*(qv_prof_cas(k2)-qv_prof_cas(k1))
-             ql_mod_cas(l)= ql_prof_cas(k2) - frac*(ql_prof_cas(k2)-ql_prof_cas(k1))
-             qi_mod_cas(l)= qi_prof_cas(k2) - frac*(qi_prof_cas(k2)-qi_prof_cas(k1))
-             u_mod_cas(l)= u_prof_cas(k2) - frac*(u_prof_cas(k2)-u_prof_cas(k1))
-             v_mod_cas(l)= v_prof_cas(k2) - frac*(v_prof_cas(k2)-v_prof_cas(k1))
-             ug_mod_cas(l)= ug_prof_cas(k2) - frac*(ug_prof_cas(k2)-ug_prof_cas(k1))
-             vg_mod_cas(l)= vg_prof_cas(k2) - frac*(vg_prof_cas(k2)-vg_prof_cas(k1))
-             temp_nudg_mod_cas(l)= temp_nudg_prof_cas(k2) - frac*(temp_nudg_prof_cas(k2)-temp_nudg_prof_cas(k1))
-             qv_nudg_mod_cas(l)= qv_nudg_prof_cas(k2) - frac*(qv_nudg_prof_cas(k2)-qv_nudg_prof_cas(k1))
-             u_nudg_mod_cas(l)= u_nudg_prof_cas(k2) - frac*(u_nudg_prof_cas(k2)-u_nudg_prof_cas(k1))
-             v_nudg_mod_cas(l)= v_nudg_prof_cas(k2) - frac*(v_nudg_prof_cas(k2)-v_nudg_prof_cas(k1))
-
-             invtau_temp_nudg_mod_cas(l)= invtau_temp_nudg_prof_cas(k2) &
-                  - frac*(invtau_temp_nudg_prof_cas(k2)-invtau_temp_nudg_prof_cas(k1))
-             invtau_qv_nudg_mod_cas(l)= invtau_qv_nudg_prof_cas(k2) - frac*(invtau_qv_nudg_prof_cas(k2)-invtau_qv_nudg_prof_cas(k1))
-             invtau_u_nudg_mod_cas(l)= invtau_u_nudg_prof_cas(k2) - frac*(invtau_u_nudg_prof_cas(k2)-invtau_u_nudg_prof_cas(k1))
-             invtau_v_nudg_mod_cas(l)= invtau_v_nudg_prof_cas(k2) - frac*(invtau_v_nudg_prof_cas(k2)-invtau_v_nudg_prof_cas(k1))
-
-             w_mod_cas(l)= vitw_prof_cas(k2) - frac*(vitw_prof_cas(k2)-vitw_prof_cas(k1))
-             omega_mod_cas(l)= omega_prof_cas(k2) - frac*(omega_prof_cas(k2)-omega_prof_cas(k1))
-             du_mod_cas(l)= du_prof_cas(k2) - frac*(du_prof_cas(k2)-du_prof_cas(k1))
-             hu_mod_cas(l)= hu_prof_cas(k2) - frac*(hu_prof_cas(k2)-hu_prof_cas(k1))
-             vu_mod_cas(l)= vu_prof_cas(k2) - frac*(vu_prof_cas(k2)-vu_prof_cas(k1))
-             dv_mod_cas(l)= dv_prof_cas(k2) - frac*(dv_prof_cas(k2)-dv_prof_cas(k1))
-             hv_mod_cas(l)= hv_prof_cas(k2) - frac*(hv_prof_cas(k2)-hv_prof_cas(k1))
-             vv_mod_cas(l)= vv_prof_cas(k2) - frac*(vv_prof_cas(k2)-vv_prof_cas(k1))
-             dt_mod_cas(l)= dt_prof_cas(k2) - frac*(dt_prof_cas(k2)-dt_prof_cas(k1))
-             ht_mod_cas(l)= ht_prof_cas(k2) - frac*(ht_prof_cas(k2)-ht_prof_cas(k1))
-             vt_mod_cas(l)= vt_prof_cas(k2) - frac*(vt_prof_cas(k2)-vt_prof_cas(k1))
-             dth_mod_cas(l)= dth_prof_cas(k2) - frac*(dth_prof_cas(k2)-dth_prof_cas(k1))
-             hth_mod_cas(l)= hth_prof_cas(k2) - frac*(hth_prof_cas(k2)-hth_prof_cas(k1))
-             vth_mod_cas(l)= vth_prof_cas(k2) - frac*(vth_prof_cas(k2)-vth_prof_cas(k1))
-             dq_mod_cas(l)= dq_prof_cas(k2) - frac*(dq_prof_cas(k2)-dq_prof_cas(k1))
-             hq_mod_cas(l)= hq_prof_cas(k2) - frac*(hq_prof_cas(k2)-hq_prof_cas(k1))
-             vq_mod_cas(l)= vq_prof_cas(k2) - frac*(vq_prof_cas(k2)-vq_prof_cas(k1))
-             dtrad_mod_cas(l)= dtrad_prof_cas(k2) - frac*(dtrad_prof_cas(k2)-dtrad_prof_cas(k1))
-
-          else !play>plev_prof_cas(1)
-
-             k1=1
-             k2=2
-             print *,'interp2_vert, k1,k2=',plev_prof_cas(k1),plev_prof_cas(k2)
-             frac1 = (play(l)-plev_prof_cas(k2))/(plev_prof_cas(k1)-plev_prof_cas(k2))
-             frac2 = (play(l)-plev_prof_cas(k1))/(plev_prof_cas(k1)-plev_prof_cas(k2))
-             t_mod_cas(l)= frac1*t_prof_cas(k1) - frac2*t_prof_cas(k2)
-             theta_mod_cas(l)= frac1*th_prof_cas(k1) - frac2*th_prof_cas(k2)
-             IF(theta_mod_cas(l)/=0) t_mod_cas(l)= theta_mod_cas(l)*(play(l)/100000.)**(RD/RCPD)
-             thv_mod_cas(l)= frac1*thv_prof_cas(k1) - frac2*thv_prof_cas(k2)
-             thl_mod_cas(l)= frac1*thl_prof_cas(k1) - frac2*thl_prof_cas(k2)
-             qv_mod_cas(l)= frac1*qv_prof_cas(k1) - frac2*qv_prof_cas(k2)
-             ql_mod_cas(l)= frac1*ql_prof_cas(k1) - frac2*ql_prof_cas(k2)
-             qi_mod_cas(l)= frac1*qi_prof_cas(k1) - frac2*qi_prof_cas(k2)
-             u_mod_cas(l)= frac1*u_prof_cas(k1) - frac2*u_prof_cas(k2)
-             v_mod_cas(l)= frac1*v_prof_cas(k1) - frac2*v_prof_cas(k2)
-             ug_mod_cas(l)= frac1*ug_prof_cas(k1) - frac2*ug_prof_cas(k2)
-             vg_mod_cas(l)= frac1*vg_prof_cas(k1) - frac2*vg_prof_cas(k2)
-             temp_nudg_mod_cas(l)= frac1*temp_nudg_prof_cas(k1) - frac2*temp_nudg_prof_cas(k2)
-             qv_nudg_mod_cas(l)= frac1*qv_nudg_prof_cas(k1) - frac2*qv_nudg_prof_cas(k2)
-             u_nudg_mod_cas(l)= frac1*u_nudg_prof_cas(k1) - frac2*u_nudg_prof_cas(k2)
-             v_nudg_mod_cas(l)= frac1*v_nudg_prof_cas(k1) - frac2*v_nudg_prof_cas(k2)
-
-             invtau_temp_nudg_mod_cas(l)= frac1*invtau_temp_nudg_prof_cas(k1) - frac2*invtau_temp_nudg_prof_cas(k2)
-             invtau_qv_nudg_mod_cas(l)= frac1*invtau_qv_nudg_prof_cas(k1) - frac2*invtau_qv_nudg_prof_cas(k2)
-             invtau_u_nudg_mod_cas(l)= frac1*invtau_u_nudg_prof_cas(k1) - frac2*invtau_u_nudg_prof_cas(k2)
-             invtau_v_nudg_mod_cas(l)= frac1*invtau_v_nudg_prof_cas(k1) - frac2*invtau_v_nudg_prof_cas(k2)
-
-             w_mod_cas(l)= frac1*vitw_prof_cas(k1) - frac2*vitw_prof_cas(k2)
-             omega_mod_cas(l)= frac1*omega_prof_cas(k1) - frac2*omega_prof_cas(k2)
-             du_mod_cas(l)= frac1*du_prof_cas(k1) - frac2*du_prof_cas(k2)
-             hu_mod_cas(l)= frac1*hu_prof_cas(k1) - frac2*hu_prof_cas(k2)
-             vu_mod_cas(l)= frac1*vu_prof_cas(k1) - frac2*vu_prof_cas(k2)
-             dv_mod_cas(l)= frac1*dv_prof_cas(k1) - frac2*dv_prof_cas(k2)
-             hv_mod_cas(l)= frac1*hv_prof_cas(k1) - frac2*hv_prof_cas(k2)
-             vv_mod_cas(l)= frac1*vv_prof_cas(k1) - frac2*vv_prof_cas(k2)
-             dt_mod_cas(l)= frac1*dt_prof_cas(k1) - frac2*dt_prof_cas(k2)
-             ht_mod_cas(l)= frac1*ht_prof_cas(k1) - frac2*ht_prof_cas(k2)
-             vt_mod_cas(l)= frac1*vt_prof_cas(k1) - frac2*vt_prof_cas(k2)
-             dth_mod_cas(l)= frac1*dth_prof_cas(k1) - frac2*dth_prof_cas(k2)
-             hth_mod_cas(l)= frac1*hth_prof_cas(k1) - frac2*hth_prof_cas(k2)
-             vth_mod_cas(l)= frac1*vth_prof_cas(k1) - frac2*vth_prof_cas(k2)
-             dq_mod_cas(l)= frac1*dq_prof_cas(k1) - frac2*dq_prof_cas(k2)
-             hq_mod_cas(l)= frac1*hq_prof_cas(k1) - frac2*hq_prof_cas(k2)
-             vq_mod_cas(l)= frac1*vq_prof_cas(k1) - frac2*vq_prof_cas(k2)
-             dtrad_mod_cas(l)= frac1*dtrad_prof_cas(k1) - frac2*dtrad_prof_cas(k2)
-
-          endif ! play.le.plev_prof_cas(1)
-
-       else ! above max altitude of forcing file
-
-          !jyg
-          fact=20.*(plev_prof_cas(nlev_cas)-play(l))/plev_prof_cas(nlev_cas) !jyg
-          fact = max(fact,0.)                                           !jyg
-          fact = exp(-fact)                                             !jyg
-          t_mod_cas(l)= t_prof_cas(nlev_cas)                            !jyg
-          theta_mod_cas(l)= th_prof_cas(nlev_cas)                       !jyg
-          thv_mod_cas(l)= thv_prof_cas(nlev_cas)                        !jyg
-          thl_mod_cas(l)= thl_prof_cas(nlev_cas)                        !jyg
-          qv_mod_cas(l)= qv_prof_cas(nlev_cas)*fact                     !jyg
-          ql_mod_cas(l)= ql_prof_cas(nlev_cas)*fact                     !jyg
-          qi_mod_cas(l)= qi_prof_cas(nlev_cas)*fact                     !jyg
-          u_mod_cas(l)= u_prof_cas(nlev_cas)*fact                       !jyg
-          v_mod_cas(l)= v_prof_cas(nlev_cas)*fact                       !jyg
-          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
-
-          invtau_temp_nudg_mod_cas(l)= invtau_temp_nudg_prof_cas(nlev_cas)            !jyg
-          invtau_qv_nudg_mod_cas(l)= invtau_qv_nudg_prof_cas(nlev_cas)                !jyg
-          invtau_u_nudg_mod_cas(l)= invtau_u_nudg_prof_cas(nlev_cas)                  !jyg
-          invtau_v_nudg_mod_cas(l)= invtau_v_nudg_prof_cas(nlev_cas)                  !jyg
-
-          thv_mod_cas(l)= thv_prof_cas(nlev_cas)                        !jyg
-          w_mod_cas(l)= 0.0                                             !jyg
-          omega_mod_cas(l)= 0.0                                         !jyg
-          du_mod_cas(l)= du_prof_cas(nlev_cas)*fact 
-          hu_mod_cas(l)= hu_prof_cas(nlev_cas)*fact                     !jyg
-          vu_mod_cas(l)= vu_prof_cas(nlev_cas)*fact                     !jyg
-          dv_mod_cas(l)= dv_prof_cas(nlev_cas)*fact
-          hv_mod_cas(l)= hv_prof_cas(nlev_cas)*fact                     !jyg
-          vv_mod_cas(l)= vv_prof_cas(nlev_cas)*fact                     !jyg
-          dt_mod_cas(l)= dt_prof_cas(nlev_cas) 
-          ht_mod_cas(l)= ht_prof_cas(nlev_cas)                          !jyg
-          vt_mod_cas(l)= vt_prof_cas(nlev_cas)                          !jyg
-          dth_mod_cas(l)= dth_prof_cas(nlev_cas) 
-          hth_mod_cas(l)= hth_prof_cas(nlev_cas)                        !jyg
-          vth_mod_cas(l)= vth_prof_cas(nlev_cas)                        !jyg
-          dq_mod_cas(l)= dq_prof_cas(nlev_cas)*fact 
-          hq_mod_cas(l)= hq_prof_cas(nlev_cas)*fact                     !jyg
-          vq_mod_cas(l)= vq_prof_cas(nlev_cas)*fact                     !jyg
-          dtrad_mod_cas(l)= dtrad_prof_cas(nlev_cas)*fact               !jyg
-
-       endif ! play
+      IF (play(l)>=plev_prof_cas(nlev_cas)) THEN
+        mxcalc = l
+        !        print *,'debut interp2, mxcalc=',mxcalc
+        k1 = 0
+        k2 = 0
+
+        IF (play(l)<=plev_prof_cas(1)) THEN
+          do k = 1, nlev_cas - 1
+            IF (play(l)<=plev_prof_cas(k).AND. play(l)>plev_prof_cas(k + 1)) THEN
+              k1 = k
+              k2 = k + 1
+            endif
+          enddo
+
+          IF (k1==0 .OR. k2==0) THEN
+            WRITE(*, *) 'PB! k1, k2 = ', k1, k2
+            WRITE(*, *) 'l,play(l) = ', l, play(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) - 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))
+          IF(theta_mod_cas(l)/=0) t_mod_cas(l) = theta_mod_cas(l) * (play(l) / 100000.)**(RD / RCPD)
+          thv_mod_cas(l) = thv_prof_cas(k2) - frac * (thv_prof_cas(k2) - thv_prof_cas(k1))
+          thl_mod_cas(l) = thl_prof_cas(k2) - frac * (thl_prof_cas(k2) - thl_prof_cas(k1))
+          qv_mod_cas(l) = qv_prof_cas(k2) - frac * (qv_prof_cas(k2) - qv_prof_cas(k1))
+          ql_mod_cas(l) = ql_prof_cas(k2) - frac * (ql_prof_cas(k2) - ql_prof_cas(k1))
+          qi_mod_cas(l) = qi_prof_cas(k2) - frac * (qi_prof_cas(k2) - qi_prof_cas(k1))
+          u_mod_cas(l) = u_prof_cas(k2) - frac * (u_prof_cas(k2) - u_prof_cas(k1))
+          v_mod_cas(l) = v_prof_cas(k2) - frac * (v_prof_cas(k2) - v_prof_cas(k1))
+          ug_mod_cas(l) = ug_prof_cas(k2) - frac * (ug_prof_cas(k2) - ug_prof_cas(k1))
+          vg_mod_cas(l) = vg_prof_cas(k2) - frac * (vg_prof_cas(k2) - vg_prof_cas(k1))
+          temp_nudg_mod_cas(l) = temp_nudg_prof_cas(k2) - frac * (temp_nudg_prof_cas(k2) - temp_nudg_prof_cas(k1))
+          qv_nudg_mod_cas(l) = qv_nudg_prof_cas(k2) - frac * (qv_nudg_prof_cas(k2) - qv_nudg_prof_cas(k1))
+          u_nudg_mod_cas(l) = u_nudg_prof_cas(k2) - frac * (u_nudg_prof_cas(k2) - u_nudg_prof_cas(k1))
+          v_nudg_mod_cas(l) = v_nudg_prof_cas(k2) - frac * (v_nudg_prof_cas(k2) - v_nudg_prof_cas(k1))
+
+          invtau_temp_nudg_mod_cas(l) = invtau_temp_nudg_prof_cas(k2) &
+                  - frac * (invtau_temp_nudg_prof_cas(k2) - invtau_temp_nudg_prof_cas(k1))
+          invtau_qv_nudg_mod_cas(l) = invtau_qv_nudg_prof_cas(k2) - frac * (invtau_qv_nudg_prof_cas(k2) - invtau_qv_nudg_prof_cas(k1))
+          invtau_u_nudg_mod_cas(l) = invtau_u_nudg_prof_cas(k2) - frac * (invtau_u_nudg_prof_cas(k2) - invtau_u_nudg_prof_cas(k1))
+          invtau_v_nudg_mod_cas(l) = invtau_v_nudg_prof_cas(k2) - frac * (invtau_v_nudg_prof_cas(k2) - invtau_v_nudg_prof_cas(k1))
+
+          w_mod_cas(l) = vitw_prof_cas(k2) - frac * (vitw_prof_cas(k2) - vitw_prof_cas(k1))
+          omega_mod_cas(l) = omega_prof_cas(k2) - frac * (omega_prof_cas(k2) - omega_prof_cas(k1))
+          du_mod_cas(l) = du_prof_cas(k2) - frac * (du_prof_cas(k2) - du_prof_cas(k1))
+          hu_mod_cas(l) = hu_prof_cas(k2) - frac * (hu_prof_cas(k2) - hu_prof_cas(k1))
+          vu_mod_cas(l) = vu_prof_cas(k2) - frac * (vu_prof_cas(k2) - vu_prof_cas(k1))
+          dv_mod_cas(l) = dv_prof_cas(k2) - frac * (dv_prof_cas(k2) - dv_prof_cas(k1))
+          hv_mod_cas(l) = hv_prof_cas(k2) - frac * (hv_prof_cas(k2) - hv_prof_cas(k1))
+          vv_mod_cas(l) = vv_prof_cas(k2) - frac * (vv_prof_cas(k2) - vv_prof_cas(k1))
+          dt_mod_cas(l) = dt_prof_cas(k2) - frac * (dt_prof_cas(k2) - dt_prof_cas(k1))
+          ht_mod_cas(l) = ht_prof_cas(k2) - frac * (ht_prof_cas(k2) - ht_prof_cas(k1))
+          vt_mod_cas(l) = vt_prof_cas(k2) - frac * (vt_prof_cas(k2) - vt_prof_cas(k1))
+          dth_mod_cas(l) = dth_prof_cas(k2) - frac * (dth_prof_cas(k2) - dth_prof_cas(k1))
+          hth_mod_cas(l) = hth_prof_cas(k2) - frac * (hth_prof_cas(k2) - hth_prof_cas(k1))
+          vth_mod_cas(l) = vth_prof_cas(k2) - frac * (vth_prof_cas(k2) - vth_prof_cas(k1))
+          dq_mod_cas(l) = dq_prof_cas(k2) - frac * (dq_prof_cas(k2) - dq_prof_cas(k1))
+          hq_mod_cas(l) = hq_prof_cas(k2) - frac * (hq_prof_cas(k2) - hq_prof_cas(k1))
+          vq_mod_cas(l) = vq_prof_cas(k2) - frac * (vq_prof_cas(k2) - vq_prof_cas(k1))
+          dtrad_mod_cas(l) = dtrad_prof_cas(k2) - frac * (dtrad_prof_cas(k2) - dtrad_prof_cas(k1))
+
+        else !play>plev_prof_cas(1)
+
+          k1 = 1
+          k2 = 2
+          print *, 'interp2_vert, k1,k2=', plev_prof_cas(k1), plev_prof_cas(k2)
+          frac1 = (play(l) - plev_prof_cas(k2)) / (plev_prof_cas(k1) - plev_prof_cas(k2))
+          frac2 = (play(l) - plev_prof_cas(k1)) / (plev_prof_cas(k1) - plev_prof_cas(k2))
+          t_mod_cas(l) = frac1 * t_prof_cas(k1) - frac2 * t_prof_cas(k2)
+          theta_mod_cas(l) = frac1 * th_prof_cas(k1) - frac2 * th_prof_cas(k2)
+          IF(theta_mod_cas(l)/=0) t_mod_cas(l) = theta_mod_cas(l) * (play(l) / 100000.)**(RD / RCPD)
+          thv_mod_cas(l) = frac1 * thv_prof_cas(k1) - frac2 * thv_prof_cas(k2)
+          thl_mod_cas(l) = frac1 * thl_prof_cas(k1) - frac2 * thl_prof_cas(k2)
+          qv_mod_cas(l) = frac1 * qv_prof_cas(k1) - frac2 * qv_prof_cas(k2)
+          ql_mod_cas(l) = frac1 * ql_prof_cas(k1) - frac2 * ql_prof_cas(k2)
+          qi_mod_cas(l) = frac1 * qi_prof_cas(k1) - frac2 * qi_prof_cas(k2)
+          u_mod_cas(l) = frac1 * u_prof_cas(k1) - frac2 * u_prof_cas(k2)
+          v_mod_cas(l) = frac1 * v_prof_cas(k1) - frac2 * v_prof_cas(k2)
+          ug_mod_cas(l) = frac1 * ug_prof_cas(k1) - frac2 * ug_prof_cas(k2)
+          vg_mod_cas(l) = frac1 * vg_prof_cas(k1) - frac2 * vg_prof_cas(k2)
+          temp_nudg_mod_cas(l) = frac1 * temp_nudg_prof_cas(k1) - frac2 * temp_nudg_prof_cas(k2)
+          qv_nudg_mod_cas(l) = frac1 * qv_nudg_prof_cas(k1) - frac2 * qv_nudg_prof_cas(k2)
+          u_nudg_mod_cas(l) = frac1 * u_nudg_prof_cas(k1) - frac2 * u_nudg_prof_cas(k2)
+          v_nudg_mod_cas(l) = frac1 * v_nudg_prof_cas(k1) - frac2 * v_nudg_prof_cas(k2)
+
+          invtau_temp_nudg_mod_cas(l) = frac1 * invtau_temp_nudg_prof_cas(k1) - frac2 * invtau_temp_nudg_prof_cas(k2)
+          invtau_qv_nudg_mod_cas(l) = frac1 * invtau_qv_nudg_prof_cas(k1) - frac2 * invtau_qv_nudg_prof_cas(k2)
+          invtau_u_nudg_mod_cas(l) = frac1 * invtau_u_nudg_prof_cas(k1) - frac2 * invtau_u_nudg_prof_cas(k2)
+          invtau_v_nudg_mod_cas(l) = frac1 * invtau_v_nudg_prof_cas(k1) - frac2 * invtau_v_nudg_prof_cas(k2)
+
+          w_mod_cas(l) = frac1 * vitw_prof_cas(k1) - frac2 * vitw_prof_cas(k2)
+          omega_mod_cas(l) = frac1 * omega_prof_cas(k1) - frac2 * omega_prof_cas(k2)
+          du_mod_cas(l) = frac1 * du_prof_cas(k1) - frac2 * du_prof_cas(k2)
+          hu_mod_cas(l) = frac1 * hu_prof_cas(k1) - frac2 * hu_prof_cas(k2)
+          vu_mod_cas(l) = frac1 * vu_prof_cas(k1) - frac2 * vu_prof_cas(k2)
+          dv_mod_cas(l) = frac1 * dv_prof_cas(k1) - frac2 * dv_prof_cas(k2)
+          hv_mod_cas(l) = frac1 * hv_prof_cas(k1) - frac2 * hv_prof_cas(k2)
+          vv_mod_cas(l) = frac1 * vv_prof_cas(k1) - frac2 * vv_prof_cas(k2)
+          dt_mod_cas(l) = frac1 * dt_prof_cas(k1) - frac2 * dt_prof_cas(k2)
+          ht_mod_cas(l) = frac1 * ht_prof_cas(k1) - frac2 * ht_prof_cas(k2)
+          vt_mod_cas(l) = frac1 * vt_prof_cas(k1) - frac2 * vt_prof_cas(k2)
+          dth_mod_cas(l) = frac1 * dth_prof_cas(k1) - frac2 * dth_prof_cas(k2)
+          hth_mod_cas(l) = frac1 * hth_prof_cas(k1) - frac2 * hth_prof_cas(k2)
+          vth_mod_cas(l) = frac1 * vth_prof_cas(k1) - frac2 * vth_prof_cas(k2)
+          dq_mod_cas(l) = frac1 * dq_prof_cas(k1) - frac2 * dq_prof_cas(k2)
+          hq_mod_cas(l) = frac1 * hq_prof_cas(k1) - frac2 * hq_prof_cas(k2)
+          vq_mod_cas(l) = frac1 * vq_prof_cas(k1) - frac2 * vq_prof_cas(k2)
+          dtrad_mod_cas(l) = frac1 * dtrad_prof_cas(k1) - frac2 * dtrad_prof_cas(k2)
+
+        endif ! play.le.plev_prof_cas(1)
+
+      else ! above max altitude of forcing file
+
+        !jyg
+        fact = 20. * (plev_prof_cas(nlev_cas) - play(l)) / plev_prof_cas(nlev_cas) !jyg
+        fact = max(fact, 0.)                                           !jyg
+        fact = exp(-fact)                                             !jyg
+        t_mod_cas(l) = t_prof_cas(nlev_cas)                            !jyg
+        theta_mod_cas(l) = th_prof_cas(nlev_cas)                       !jyg
+        thv_mod_cas(l) = thv_prof_cas(nlev_cas)                        !jyg
+        thl_mod_cas(l) = thl_prof_cas(nlev_cas)                        !jyg
+        qv_mod_cas(l) = qv_prof_cas(nlev_cas) * fact                     !jyg
+        ql_mod_cas(l) = ql_prof_cas(nlev_cas) * fact                     !jyg
+        qi_mod_cas(l) = qi_prof_cas(nlev_cas) * fact                     !jyg
+        u_mod_cas(l) = u_prof_cas(nlev_cas) * fact                       !jyg
+        v_mod_cas(l) = v_prof_cas(nlev_cas) * fact                       !jyg
+        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
+
+        invtau_temp_nudg_mod_cas(l) = invtau_temp_nudg_prof_cas(nlev_cas)            !jyg
+        invtau_qv_nudg_mod_cas(l) = invtau_qv_nudg_prof_cas(nlev_cas)                !jyg
+        invtau_u_nudg_mod_cas(l) = invtau_u_nudg_prof_cas(nlev_cas)                  !jyg
+        invtau_v_nudg_mod_cas(l) = invtau_v_nudg_prof_cas(nlev_cas)                  !jyg
+
+        thv_mod_cas(l) = thv_prof_cas(nlev_cas)                        !jyg
+        w_mod_cas(l) = 0.0                                             !jyg
+        omega_mod_cas(l) = 0.0                                         !jyg
+        du_mod_cas(l) = du_prof_cas(nlev_cas) * fact
+        hu_mod_cas(l) = hu_prof_cas(nlev_cas) * fact                     !jyg
+        vu_mod_cas(l) = vu_prof_cas(nlev_cas) * fact                     !jyg
+        dv_mod_cas(l) = dv_prof_cas(nlev_cas) * fact
+        hv_mod_cas(l) = hv_prof_cas(nlev_cas) * fact                     !jyg
+        vv_mod_cas(l) = vv_prof_cas(nlev_cas) * fact                     !jyg
+        dt_mod_cas(l) = dt_prof_cas(nlev_cas)
+        ht_mod_cas(l) = ht_prof_cas(nlev_cas)                          !jyg
+        vt_mod_cas(l) = vt_prof_cas(nlev_cas)                          !jyg
+        dth_mod_cas(l) = dth_prof_cas(nlev_cas)
+        hth_mod_cas(l) = hth_prof_cas(nlev_cas)                        !jyg
+        vth_mod_cas(l) = vth_prof_cas(nlev_cas)                        !jyg
+        dq_mod_cas(l) = dq_prof_cas(nlev_cas) * fact
+        hq_mod_cas(l) = hq_prof_cas(nlev_cas) * fact                     !jyg
+        vq_mod_cas(l) = vq_prof_cas(nlev_cas) * fact                     !jyg
+        dtrad_mod_cas(l) = dtrad_prof_cas(nlev_cas) * fact               !jyg
+
+      endif ! play
 
     enddo ! l
@@ -1158,55 +1141,48 @@
     ! for variables defined at layer interfaces (EV):
 
-
-    do l = 1, llm+1
-
-       IF (plev(l)>=plev_prof_cas(nlev_cas)) THEN
-          mxcalc=l
-          k1=0
-          k2=0
-
-          IF (plev(l)<=plev_prof_cas(1)) THEN
-             do k = 1, nlev_cas-1
-                IF (plev(l)<=plev_prof_cas(k).AND. plev(l)>plev_prof_cas(k+1)) THEN
-                   k1=k
-                   k2=k+1
-                endif
-             enddo
-
-             IF (k1==0 .OR. k2==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
-             frac1 = (play(l)-plev_prof_cas(k2))/(plev_prof_cas(k1)-plev_prof_cas(k2))
-             frac2 = (play(l)-plev_prof_cas(k1))/(plev_prof_cas(k1)-plev_prof_cas(k2))
-             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
+    do l = 1, llm + 1
+
+      IF (plev(l)>=plev_prof_cas(nlev_cas)) THEN
+        mxcalc = l
+        k1 = 0
+        k2 = 0
+
+        IF (plev(l)<=plev_prof_cas(1)) THEN
+          do k = 1, nlev_cas - 1
+            IF (plev(l)<=plev_prof_cas(k).AND. plev(l)>plev_prof_cas(k + 1)) THEN
+              k1 = k
+              k2 = k + 1
+            endif
+          enddo
+
+          IF (k1==0 .OR. k2==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
+          frac1 = (play(l) - plev_prof_cas(k2)) / (plev_prof_cas(k1) - plev_prof_cas(k2))
+          frac2 = (play(l) - plev_prof_cas(k1)) / (plev_prof_cas(k1) - plev_prof_cas(k2))
+          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
 
-
-
-
   END SUBROUTINE  interp2_case_vertical_std
-  !***************************************************************************** 
-
-
-
+  !*****************************************************************************
 
 
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/old_1DUTILS_read_interp.h
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/old_1DUTILS_read_interp.h	(revision 5143)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/old_1DUTILS_read_interp.h	(revision 5144)
@@ -2372,7 +2372,8 @@
             nf90_inq_dimid,nf90_inquire_dimension
 
+            USE lmdz_yomcst
+
       IMPLICIT NONE
 
-      INCLUDE "YOMCST.h"
 
       INTEGER ntime,nlevel
@@ -2945,7 +2946,8 @@
 
       SUBROUTINE read_circ(nlev_circ,cf,lwp,iwp,reliq,reice,t,z,p,pm,h2o,o3,sza)
+
+      USE lmdz_yomcst
       
       parameter (ncm_1=49180)
-      INCLUDE "YOMCST.h"
 
       REAL albsfc(ncm_1), albsfc_w(ncm_1)
@@ -3051,5 +3053,5 @@
       SUBROUTINE read_rtmip(nlev_rtmip,play,plev,t,h2o,o3)
       
-      INCLUDE "YOMCST.h"
+      USE lmdz_yomcst
 
       REAL t(nlev_rtmip), pt(nlev_rtmip),pb(nlev_rtmip),h2o(nlev_rtmip), o3(nlev_rtmip)
