Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/abort_gcm.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d/abort_gcm.F90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d/abort_gcm.F90	(revision 5118)
@@ -8,6 +8,5 @@
   !! ug Pour les sorties XIOS
   USE lmdz_wxios
-
-  include "iniprint.h"
+  USE lmdz_iniprint, ONLY: lunout, prt_level
 
   !
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/advtrac.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d/advtrac.f90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d/advtrac.f90	(revision 5118)
@@ -16,4 +16,5 @@
   USE lmdz_description, ONLY: descript
   USE lmdz_libmath, ONLY: minmax
+  USE lmdz_iniprint, ONLY: lunout, prt_level
 
   IMPLICIT NONE
@@ -23,5 +24,4 @@
   include "comdissip.h"
   include "comgeom2.h"
-  include "iniprint.h"
 
   !---------------------------------------------------------------------------
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/bilan_dyn.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d/bilan_dyn.F90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d/bilan_dyn.F90	(revision 5118)
@@ -13,4 +13,5 @@
   USE comvert_mod, ONLY: presnivs
   USE temps_mod, ONLY: annee_ref, day_ref, itau_dyn
+  USE lmdz_iniprint, ONLY: lunout, prt_level
 
   IMPLICIT NONE
@@ -19,5 +20,4 @@
   include "paramet.h"
   include "comgeom2.h"
-  include "iniprint.h"
 
   !====================================================================
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/conf_gcm.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d/conf_gcm.f90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d/conf_gcm.f90	(revision 5118)
@@ -16,4 +16,5 @@
           alphax, alphay, taux, tauy
   USE temps_mod, ONLY: calend, year_len
+  USE lmdz_iniprint, ONLY: lunout, prt_level
 
   IMPLICIT NONE
@@ -35,5 +36,4 @@
   include "paramet.h"
   include "comdissnew.h"
-  include "iniprint.h"
 
   !   local:
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/dynetat0.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d/dynetat0.F90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d/dynetat0.F90	(revision 5118)
@@ -21,4 +21,5 @@
   USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_INCA
   USE lmdz_description, ONLY: descript
+  USE lmdz_iniprint, ONLY: lunout, prt_level
 
   IMPLICIT NONE
@@ -26,5 +27,4 @@
   include "paramet.h"
   include "comgeom2.h"
-  include "iniprint.h"
 !===============================================================================
 ! Arguments:
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/dynredem.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d/dynredem.F90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d/dynredem.F90	(revision 5118)
@@ -1,60 +1,60 @@
-SUBROUTINE dynredem0(fichnom,iday_end,phis)
-
-!-------------------------------------------------------------------------------
-! Write the NetCDF restart file (initialization).
-!-------------------------------------------------------------------------------
+SUBROUTINE dynredem0(fichnom, iday_end, phis)
+
+  !-------------------------------------------------------------------------------
+  ! Write the NetCDF restart file (initialization).
+  !-------------------------------------------------------------------------------
   USE IOIPSL
   USE lmdz_strings, ONLY: maxlen
   USE infotrac, ONLY: nqtot, tracers
-  USE netcdf, ONLY: nf90_create, nf90_def_dim, nf90_inq_varid, nf90_global,    &
-                    nf90_close,  nf90_put_att, nf90_unlimited, nf90_clobber,   &
-                    nf90_64bit_offset
+  USE netcdf, ONLY: nf90_create, nf90_def_dim, nf90_inq_varid, nf90_global, &
+          nf90_close, nf90_put_att, nf90_unlimited, nf90_clobber, &
+          nf90_64bit_offset
   USE dynredem_mod, ONLY: cre_var, put_var1, put_var2, err, modname, fil
-  USE comvert_mod,  ONLY: ap, bp, presnivs, pa, preff, nivsig, nivsigs
+  USE comvert_mod, ONLY: ap, bp, presnivs, pa, preff, nivsig, nivsigs
   USE comconst_mod, ONLY: cpp, daysec, dtvr, g, kappa, omeg, rad
   USE logic_mod, ONLY: fxyhypb, ysinus
-  USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy, &
-                              taux,tauy
+  USE serre_mod, ONLY: clon, clat, grossismx, grossismy, dzoomx, dzoomy, &
+          taux, tauy
   USE temps_mod, ONLY: annee_ref, day_ref, itau_dyn, itaufin, start_time
-  USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
+  USE ener_mod, ONLY: etot0, ptot0, ztot0, stot0, ang0
   USE lmdz_description, ONLY: descript
-  
+  USE lmdz_iniprint, ONLY: lunout, prt_level
+
   IMPLICIT NONE
   include "dimensions.h"
   include "paramet.h"
   include "comgeom2.h"
-  include "iniprint.h"
-!===============================================================================
-! Arguments:
-  CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
-  INTEGER,          INTENT(IN) :: iday_end         !--- 
-  REAL,             INTENT(IN) :: phis(iip1, jjp1) !--- GROUND GEOPOTENTIAL
-!===============================================================================
-! Local variables:
+  !===============================================================================
+  ! Arguments:
+  CHARACTER(LEN = *), INTENT(IN) :: fichnom          !--- FILE NAME
+  INTEGER, INTENT(IN) :: iday_end         !---
+  REAL, INTENT(IN) :: phis(iip1, jjp1) !--- GROUND GEOPOTENTIAL
+  !===============================================================================
+  ! Local variables:
   INTEGER :: iq
-  INTEGER, PARAMETER :: length=100
-  REAL    :: tab_cntrl(length)                     !--- RUN PARAMETERS TABLE
-!   For NetCDF:
-  CHARACTER(LEN=maxlen) :: unites
+  INTEGER, PARAMETER :: length = 100
+  REAL :: tab_cntrl(length)                     !--- RUN PARAMETERS TABLE
+  !   For NetCDF:
+  CHARACTER(LEN = maxlen) :: unites
   INTEGER :: indexID
   INTEGER :: rlonuID, rlonvID, rlatuID, rlatvID
   INTEGER :: sID, sigID, nID, timID
   INTEGER :: yyears0, jjour0, mmois0
-  REAL    :: zjulian, hours
-!===============================================================================
-  modname='dynredem0'; fil=fichnom
+  REAL :: zjulian, hours
+  !===============================================================================
+  modname = 'dynredem0'; fil = fichnom
   CALL ymds2ju(annee_ref, 1, iday_end, 0.0, zjulian)
   CALL ju2ymds(zjulian, yyears0, mmois0, jjour0, hours)
 
-  tab_cntrl(:)  = 0.
-  tab_cntrl(1)  = REAL(iim)
-  tab_cntrl(2)  = REAL(jjm)
-  tab_cntrl(3)  = REAL(llm)
-  tab_cntrl(4)  = REAL(day_ref)
-  tab_cntrl(5)  = REAL(annee_ref)
-  tab_cntrl(6)  = rad
-  tab_cntrl(7)  = omeg
-  tab_cntrl(8)  = g
-  tab_cntrl(9)  = cpp
+  tab_cntrl(:) = 0.
+  tab_cntrl(1) = REAL(iim)
+  tab_cntrl(2) = REAL(jjm)
+  tab_cntrl(3) = REAL(llm)
+  tab_cntrl(4) = REAL(day_ref)
+  tab_cntrl(5) = REAL(annee_ref)
+  tab_cntrl(6) = rad
+  tab_cntrl(7) = omeg
+  tab_cntrl(8) = g
+  tab_cntrl(9) = cpp
   tab_cntrl(10) = kappa
   tab_cntrl(11) = daysec
@@ -68,5 +68,5 @@
   tab_cntrl(19) = preff
 
-!    .....    parameters for zoom    ......   
+  !    .....    parameters for zoom    ......
   tab_cntrl(20) = clon
   tab_cntrl(21) = clat
@@ -74,5 +74,5 @@
   tab_cntrl(23) = grossismy
 
-  IF ( fxyhypb )   THEN
+  IF (fxyhypb)   THEN
     tab_cntrl(24) = 1.
     tab_cntrl(25) = dzoomx
@@ -88,60 +88,60 @@
     tab_cntrl(28) = 0.
     tab_cntrl(29) = 0.
-    IF( ysinus )  tab_cntrl(27) = 1.
+    IF(ysinus)  tab_cntrl(27) = 1.
   END IF
   tab_cntrl(30) = REAL(iday_end)
   tab_cntrl(31) = REAL(itau_dyn + itaufin)
-! start_time: start_time of simulation (not necessarily 0.)
+  ! start_time: start_time of simulation (not necessarily 0.)
   tab_cntrl(32) = start_time
 
-!--- File creation
-  CALL err(nf90_create(fichnom,IOR(nf90_clobber,nf90_64bit_offset),nid))
-
-!--- Some global attributes
-  CALL err(nf90_put_att(nid,nf90_global,"title","Fichier demarrage dynamique"))
-
-!--- Dimensions
-  CALL err(nf90_def_dim(nid,"index", length, indexID))
-  CALL err(nf90_def_dim(nid,"rlonu", iip1,   rlonuID))
-  CALL err(nf90_def_dim(nid,"rlatu", jjp1,   rlatuID))
-  CALL err(nf90_def_dim(nid,"rlonv", iip1,   rlonvID))
-  CALL err(nf90_def_dim(nid,"rlatv", jjm,    rlatvID))
-  CALL err(nf90_def_dim(nid,"sigs",  llm,        sID))
-  CALL err(nf90_def_dim(nid,"sig",   llmp1,    sigID))
-  CALL err(nf90_def_dim(nid,"temps", nf90_unlimited, timID))
-
-!--- Define and save invariant fields
-  CALL put_var1(nid,"controle","Parametres de controle" ,[indexID],tab_cntrl)
-  CALL put_var1(nid,"rlonu"   ,"Longitudes des points U",[rlonuID],rlonu)
-  CALL put_var1(nid,"rlatu"   ,"Latitudes des points U" ,[rlatuID],rlatu)
-  CALL put_var1(nid,"rlonv"   ,"Longitudes des points V",[rlonvID],rlonv)
-  CALL put_var1(nid,"rlatv"   ,"Latitudes des points V" ,[rlatvID],rlatv)
-  CALL put_var1(nid,"nivsigs" ,"Numero naturel des couches s"    ,[sID]  ,nivsigs)
-  CALL put_var1(nid,"nivsig"  ,"Numero naturel des couches sigma",[sigID],nivsig)
-  CALL put_var1(nid,"ap"      ,"Coefficient A pour hybride"      ,[sigID],ap)
-  CALL put_var1(nid,"bp"      ,"Coefficient B pour hybride"      ,[sigID],bp)
-  CALL put_var1(nid,"presnivs",""                                ,[sID]  ,presnivs)
-! covariant <-> contravariant <-> natural conversion coefficients
-  CALL put_var2(nid,"cu","Coefficient de passage pour U",[rlonuID,rlatuID],cu)
-  CALL put_var2(nid,"cv","Coefficient de passage pour V",[rlonvID,rlatvID],cv)
-  CALL put_var2(nid,"aire","Aires de chaque maille"     ,[rlonvID,rlatuID],aire)
-  CALL put_var2(nid,"phisinit","Geopotentiel au sol"    ,[rlonvID,rlatuID],phis)
-
-!--- Define fields saved later
-  WRITE(unites,"('days since ',i4,'-',i2.2,'-',i2.2,' 00:00:00')") &
-               yyears0,mmois0,jjour0
-  CALL cre_var(nid,"temps","Temps de simulation",[timID],unites)
-  CALL cre_var(nid,"ucov" ,"Vitesse U"  ,[rlonuID,rlatuID,sID,timID])
-  CALL cre_var(nid,"vcov" ,"Vitesse V"  ,[rlonvID,rlatvID,sID,timID])
-  CALL cre_var(nid,"teta" ,"Temperature",[rlonvID,rlatuID,sID,timID])
-  DO iq=1,nqtot
-    CALL cre_var(nid,tracers(iq)%name,tracers(iq)%longName,[rlonvID,rlatuID,sID,timID])
+  !--- File creation
+  CALL err(nf90_create(fichnom, IOR(nf90_clobber, nf90_64bit_offset), nid))
+
+  !--- Some global attributes
+  CALL err(nf90_put_att(nid, nf90_global, "title", "Fichier demarrage dynamique"))
+
+  !--- Dimensions
+  CALL err(nf90_def_dim(nid, "index", length, indexID))
+  CALL err(nf90_def_dim(nid, "rlonu", iip1, rlonuID))
+  CALL err(nf90_def_dim(nid, "rlatu", jjp1, rlatuID))
+  CALL err(nf90_def_dim(nid, "rlonv", iip1, rlonvID))
+  CALL err(nf90_def_dim(nid, "rlatv", jjm, rlatvID))
+  CALL err(nf90_def_dim(nid, "sigs", llm, sID))
+  CALL err(nf90_def_dim(nid, "sig", llmp1, sigID))
+  CALL err(nf90_def_dim(nid, "temps", nf90_unlimited, timID))
+
+  !--- Define and save invariant fields
+  CALL put_var1(nid, "controle", "Parametres de controle", [indexID], tab_cntrl)
+  CALL put_var1(nid, "rlonu", "Longitudes des points U", [rlonuID], rlonu)
+  CALL put_var1(nid, "rlatu", "Latitudes des points U", [rlatuID], rlatu)
+  CALL put_var1(nid, "rlonv", "Longitudes des points V", [rlonvID], rlonv)
+  CALL put_var1(nid, "rlatv", "Latitudes des points V", [rlatvID], rlatv)
+  CALL put_var1(nid, "nivsigs", "Numero naturel des couches s", [sID], nivsigs)
+  CALL put_var1(nid, "nivsig", "Numero naturel des couches sigma", [sigID], nivsig)
+  CALL put_var1(nid, "ap", "Coefficient A pour hybride", [sigID], ap)
+  CALL put_var1(nid, "bp", "Coefficient B pour hybride", [sigID], bp)
+  CALL put_var1(nid, "presnivs", "", [sID], presnivs)
+  ! covariant <-> contravariant <-> natural conversion coefficients
+  CALL put_var2(nid, "cu", "Coefficient de passage pour U", [rlonuID, rlatuID], cu)
+  CALL put_var2(nid, "cv", "Coefficient de passage pour V", [rlonvID, rlatvID], cv)
+  CALL put_var2(nid, "aire", "Aires de chaque maille", [rlonvID, rlatuID], aire)
+  CALL put_var2(nid, "phisinit", "Geopotentiel au sol", [rlonvID, rlatuID], phis)
+
+  !--- Define fields saved later
+  WRITE(unites, "('days since ',i4,'-',i2.2,'-',i2.2,' 00:00:00')") &
+          yyears0, mmois0, jjour0
+  CALL cre_var(nid, "temps", "Temps de simulation", [timID], unites)
+  CALL cre_var(nid, "ucov", "Vitesse U", [rlonuID, rlatuID, sID, timID])
+  CALL cre_var(nid, "vcov", "Vitesse V", [rlonvID, rlatvID, sID, timID])
+  CALL cre_var(nid, "teta", "Temperature", [rlonvID, rlatuID, sID, timID])
+  DO iq = 1, nqtot
+    CALL cre_var(nid, tracers(iq)%name, tracers(iq)%longName, [rlonvID, rlatuID, sID, timID])
   END DO
-  CALL cre_var(nid,"masse","Masse d air"    ,[rlonvID,rlatuID,sID,timID])
-  CALL cre_var(nid,"ps"   ,"Pression au sol",[rlonvID,rlatuID    ,timID])
+  CALL cre_var(nid, "masse", "Masse d air", [rlonvID, rlatuID, sID, timID])
+  CALL cre_var(nid, "ps", "Pression au sol", [rlonvID, rlatuID, timID])
   CALL err(nf90_close (nid))
 
-  WRITE(lunout,*)TRIM(modname)//': iim,jjm,llm,iday_end',iim,jjm,llm,iday_end
-  WRITE(lunout,*)TRIM(modname)//': rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa
+  WRITE(lunout, *)TRIM(modname) // ': iim,jjm,llm,iday_end', iim, jjm, llm, iday_end
+  WRITE(lunout, *)TRIM(modname) // ': rad,omeg,g,cpp,kappa', rad, omeg, g, cpp, kappa
 
 END SUBROUTINE dynredem0
@@ -152,91 +152,90 @@
 !-------------------------------------------------------------------------------
 
-SUBROUTINE dynredem1(fichnom,time,vcov,ucov,teta,q,masse,ps)
-
-!-------------------------------------------------------------------------------
-! Purpose: Write the NetCDF restart file (append).
-!-------------------------------------------------------------------------------
+SUBROUTINE dynredem1(fichnom, time, vcov, ucov, teta, q, masse, ps)
+
+  !-------------------------------------------------------------------------------
+  ! Purpose: Write the NetCDF restart file (append).
+  !-------------------------------------------------------------------------------
   USE lmdz_strings, ONLY: maxlen
   USE infotrac, ONLY: nqtot, tracers, type_trac
   USE control_mod
-  USE netcdf,   ONLY: nf90_open,  nf90_nowrite, nf90_get_var, nf90_inq_varid,  &
-                      nf90_close, nf90_write,   nf90_put_var, nf90_noerr
+  USE netcdf, ONLY: nf90_open, nf90_nowrite, nf90_get_var, nf90_inq_varid, &
+          nf90_close, nf90_write, nf90_put_var, nf90_noerr
   USE dynredem_mod, ONLY: dynredem_write_u, dynredem_write_v, dynredem_read_u, &
-                          err, modname, fil, msg
+          err, modname, fil, msg
   USE temps_mod, ONLY: itau_dyn, itaufin
   USE lmdz_description, ONLY: descript
-  
+  USE lmdz_iniprint, ONLY: lunout, prt_level
+
   IMPLICIT NONE
   include "dimensions.h"
   include "paramet.h"
   include "comgeom.h"
-  include "iniprint.h"
-!===============================================================================
-! Arguments:
-  CHARACTER(LEN=*), INTENT(IN) :: fichnom              !-- FILE NAME
-  REAL, INTENT(IN)    ::  time                         !-- TIME
-  REAL, INTENT(IN)    ::  vcov(iip1,jjm, llm)          !-- V COVARIANT WIND
-  REAL, INTENT(IN)    ::  ucov(iip1,jjp1,llm)          !-- U COVARIANT WIND
-  REAL, INTENT(IN)    ::  teta(iip1,jjp1,llm)          !-- POTENTIAL TEMPERATURE
-  REAL, INTENT(INOUT) ::     q(iip1,jjp1,llm,nqtot)    !-- TRACERS
-  REAL, INTENT(IN)    :: masse(iip1,jjp1,llm)          !-- MASS PER CELL
-  REAL, INTENT(IN)    ::    ps(iip1,jjp1)              !-- GROUND PRESSURE
-!===============================================================================
-! Local variables:
+  !===============================================================================
+  ! Arguments:
+  CHARACTER(LEN = *), INTENT(IN) :: fichnom              !-- FILE NAME
+  REAL, INTENT(IN) :: time                         !-- TIME
+  REAL, INTENT(IN) :: vcov(iip1, jjm, llm)          !-- V COVARIANT WIND
+  REAL, INTENT(IN) :: ucov(iip1, jjp1, llm)          !-- U COVARIANT WIND
+  REAL, INTENT(IN) :: teta(iip1, jjp1, llm)          !-- POTENTIAL TEMPERATURE
+  REAL, INTENT(INOUT) :: q(iip1, jjp1, llm, nqtot)    !-- TRACERS
+  REAL, INTENT(IN) :: masse(iip1, jjp1, llm)          !-- MASS PER CELL
+  REAL, INTENT(IN) :: ps(iip1, jjp1)              !-- GROUND PRESSURE
+  !===============================================================================
+  ! Local variables:
   INTEGER :: iq, nid, vID, ierr, nid_trac, vID_trac
-  INTEGER, SAVE :: nb=0
-  INTEGER, PARAMETER :: length=100
-  REAL               :: tab_cntrl(length) ! tableau des parametres du run
-  CHARACTER(LEN=maxlen) :: var, dum
-  LOGICAL            :: lread_inca
-!===============================================================================
-
-  modname='dynredem1'; fil=fichnom
-  CALL err(nf90_open(fil,nf90_write,nid),"open",fil)
-
-!--- Write/extend time coordinate
+  INTEGER, SAVE :: nb = 0
+  INTEGER, PARAMETER :: length = 100
+  REAL :: tab_cntrl(length) ! tableau des parametres du run
+  CHARACTER(LEN = maxlen) :: var, dum
+  LOGICAL :: lread_inca
+  !===============================================================================
+
+  modname = 'dynredem1'; fil = fichnom
+  CALL err(nf90_open(fil, nf90_write, nid), "open", fil)
+
+  !--- Write/extend time coordinate
   nb = nb + 1
-  var="temps"
-  CALL err(nf90_inq_varid(nid,var,vID),"inq",var)
-  CALL err(nf90_put_var(nid,vID,[time]),"put",var)
-  WRITE(lunout,*)TRIM(modname)//": Saving for ", nb, time
-
-!--- Rewrite control table (itaufin undefined in dynredem0)
-  var="controle"
-  CALL err(nf90_inq_varid(nid,var,vID),"inq",var)
-  CALL err(nf90_get_var(nid,vID,tab_cntrl),"get",var)
-  tab_cntrl(31)=DBLE(itau_dyn + itaufin)
-  CALL err(nf90_inq_varid(nid,var,vID),"inq",var)
-  CALL err(nf90_put_var(nid,vID,tab_cntrl),"put",var)
-
-!--- Save fields
-  CALL dynredem_write_u(nid,"ucov" ,ucov ,llm)
-  CALL dynredem_write_v(nid,"vcov" ,vcov ,llm)
-  CALL dynredem_write_u(nid,"teta" ,teta ,llm)
-  CALL dynredem_write_u(nid,"masse",masse,llm)
-  CALL dynredem_write_u(nid,"ps"   ,ps   ,1)
-
-!--- Tracers in file "start_trac.nc" (added by Anne)
-  lread_inca=.FALSE.; fil="start_trac.nc"
-  IF(ANY(type_trac == ['inca','inco'])) INQUIRE(FILE=fil,EXIST=lread_inca)
-  IF(lread_inca) CALL err(nf90_open(fil,nf90_nowrite,nid_trac),"open")
-
-!--- Save tracers
-  DO iq=1,nqtot; var=TRIM(tracers(iq)%name); ierr=-1
-    IF(lread_inca) THEN                  !--- Possibly read from "start_trac.nc"
-      fil="start_trac.nc"
-      ierr=nf90_inq_varid(nid_trac,var,vID_trac)
-      dum='inq'; IF(ierr==nf90_noerr) dum='fnd'
-      WRITE(lunout,*)msg(dum,var)
-
-
-      IF(ierr==nf90_noerr) CALL dynredem_read_u(nid_trac,var,q(:,:,:,iq),llm)
-    END IF
-    fil=fichnom
-    CALL dynredem_write_u(nid,var,q(:,:,:,iq),llm)
+  var = "temps"
+  CALL err(nf90_inq_varid(nid, var, vID), "inq", var)
+  CALL err(nf90_put_var(nid, vID, [time]), "put", var)
+  WRITE(lunout, *)TRIM(modname) // ": Saving for ", nb, time
+
+  !--- Rewrite control table (itaufin undefined in dynredem0)
+  var = "controle"
+  CALL err(nf90_inq_varid(nid, var, vID), "inq", var)
+  CALL err(nf90_get_var(nid, vID, tab_cntrl), "get", var)
+  tab_cntrl(31) = DBLE(itau_dyn + itaufin)
+  CALL err(nf90_inq_varid(nid, var, vID), "inq", var)
+  CALL err(nf90_put_var(nid, vID, tab_cntrl), "put", var)
+
+  !--- Save fields
+  CALL dynredem_write_u(nid, "ucov", ucov, llm)
+  CALL dynredem_write_v(nid, "vcov", vcov, llm)
+  CALL dynredem_write_u(nid, "teta", teta, llm)
+  CALL dynredem_write_u(nid, "masse", masse, llm)
+  CALL dynredem_write_u(nid, "ps", ps, 1)
+
+  !--- Tracers in file "start_trac.nc" (added by Anne)
+  lread_inca = .FALSE.; fil = "start_trac.nc"
+  IF(ANY(type_trac == ['inca', 'inco'])) INQUIRE(FILE = fil, EXIST = lread_inca)
+  IF(lread_inca) CALL err(nf90_open(fil, nf90_nowrite, nid_trac), "open")
+
+  !--- Save tracers
+  DO iq = 1, nqtot; var = TRIM(tracers(iq)%name); ierr = -1
+  IF(lread_inca) THEN                  !--- Possibly read from "start_trac.nc"
+    fil = "start_trac.nc"
+    ierr = nf90_inq_varid(nid_trac, var, vID_trac)
+    dum = 'inq'; IF(ierr==nf90_noerr) dum = 'fnd'
+    WRITE(lunout, *)msg(dum, var)
+
+    IF(ierr==nf90_noerr) CALL dynredem_read_u(nid_trac, var, q(:, :, :, iq), llm)
+  END IF
+  fil = fichnom
+  CALL dynredem_write_u(nid, var, q(:, :, :, iq), llm)
   END DO
-  CALL err(nf90_close(nid),"close")
-  fil="start_trac.nc"
-  IF(lread_inca) CALL err(nf90_close(nid_trac),"close")
+  CALL err(nf90_close(nid), "close")
+  fil = "start_trac.nc"
+  IF(lread_inca) CALL err(nf90_close(nid_trac), "close")
 
 END SUBROUTINE dynredem1
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/fluxstokenc.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d/fluxstokenc.F90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d/fluxstokenc.F90	(revision 5118)
@@ -6,4 +6,5 @@
 
   USE IOIPSL
+  USE lmdz_iniprint, ONLY: lunout, prt_level
   !
   ! Auteur :  F. Hourdin
@@ -18,5 +19,4 @@
   include "comgeom.h"
   include "tracstoke.h"
-  include "iniprint.h"
 
   REAL :: time_step, t_wrt, t_ops
@@ -157,4 +157,3 @@
   ENDIF ! if iadvtr.EQ.istdyn
 
-
 END SUBROUTINE fluxstokenc
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/friction.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d/friction.F90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d/friction.F90	(revision 5118)
@@ -7,4 +7,5 @@
   USE IOIPSL
   USE comconst_mod, ONLY: pi
+  USE lmdz_iniprint, ONLY: lunout, prt_level
   IMPLICIT NONE
 
@@ -24,5 +25,4 @@
   include "paramet.h"
   include "comgeom2.h"
-  include "iniprint.h"
   include "academic.h"
 
@@ -125,5 +125,4 @@
   ENDIF
 
-
 END SUBROUTINE friction
 
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/gcm.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d/gcm.F90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d/gcm.F90	(revision 5118)
@@ -1,3 +1,2 @@
-
 ! $Id$
 
@@ -13,6 +12,6 @@
   USE control_mod
   USE mod_const_mpi, ONLY: COMM_LMDZ
-  USE temps_mod, ONLY: calend,start_time,annee_ref,day_ref, &
-                     itau_dyn,itau_phy,day_ini,jD_ref,jH_ref,day_end
+  USE temps_mod, ONLY: calend, start_time, annee_ref, day_ref, &
+          itau_dyn, itau_phy, day_ini, jD_ref, jH_ref, day_end
   USE comconst_mod, ONLY: cpp, daysec, dtphys, dtvr, g, r, rad
   USE logic_mod, ONLY: ecripar, iflag_phys, read_start
@@ -21,5 +20,5 @@
   USE lmdz_description, ONLY: descript
 
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
   ! A nettoyer. On ne veut qu'une ou deux routines d'interface 
@@ -27,5 +26,6 @@
   ! AB 2024/07/20: remplace CPP key by fortran logical, but ^ still relevant, see later use of iniphys later on
   USE iniphysiq_mod, ONLY: iniphysiq
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  USE lmdz_iniprint, ONLY: lunout, prt_level
 
   IMPLICIT NONE
@@ -65,5 +65,4 @@
   include "comdissnew.h"
   include "comgeom.h"
-  include "iniprint.h"
   include "tracstoke.h"
 
@@ -71,13 +70,13 @@
 
   !   variables dynamiques
-  REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
-  REAL teta(ip1jmp1,llm)                 ! temperature potentielle 
-  REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes
+  REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants
+  REAL teta(ip1jmp1, llm)                 ! temperature potentielle
+  REAL, ALLOCATABLE, DIMENSION(:, :, :) :: q! champs advectes
   REAL ps(ip1jmp1)                       ! pression  au sol
-!  REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
-  REAL masse(ip1jmp1,llm)                ! masse d'air
+  !  REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
+  REAL masse(ip1jmp1, llm)                ! masse d'air
   REAL phis(ip1jmp1)                     ! geopotentiel au sol
-!  REAL phi(ip1jmp1,llm)                  ! geopotentiel
-!  REAL w(ip1jmp1,llm)                    ! vitesse verticale
+  !  REAL phi(ip1jmp1,llm)                  ! geopotentiel
+  !  REAL w(ip1jmp1,llm)                    ! vitesse verticale
 
   ! variables dynamiques intermediaire pour le transport
@@ -89,5 +88,4 @@
 
   LOGICAL lafin
-
 
   REAL time_step, t_wrt, t_ops
@@ -101,14 +99,13 @@
   !     tansformation d'energie cinetique en energie thermique
   !     cree par la dissipation
-!  REAL dhecdt(ip1jmp1,llm)
+  !  REAL dhecdt(ip1jmp1,llm)
   !      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
   !      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
-!  CHARACTER (len=15) :: ztit
+  !  CHARACTER (len=15) :: ztit
   !-jld 
 
-
-  CHARACTER (LEN=80) :: dynhist_file, dynhistave_file
-  CHARACTER (LEN=20) :: modname
-  CHARACTER (LEN=80) :: abort_message
+  CHARACTER (LEN = 80) :: dynhist_file, dynhistave_file
+  CHARACTER (LEN = 20) :: modname
+  CHARACTER (LEN = 80) :: abort_message
   ! locales pour gestion du temps
   INTEGER :: an, mois, jour
@@ -123,5 +120,5 @@
   modname = 'gcm'
   descript = 'Run GCM LMDZ'
-  lafin    = .FALSE.
+  lafin = .FALSE.
   dynhist_file = 'dyn_hist.nc'
   dynhistave_file = 'dyn_hist_ave.nc'
@@ -133,17 +130,17 @@
   !  ---------------------------------------
 
-  CALL conf_gcm( 99, .TRUE.)
+  CALL conf_gcm(99, .TRUE.)
 
   IF (mod(iphysiq, iperiod) /= 0) CALL abort_gcm("conf_gcm", &
-       "iphysiq must be a multiple of iperiod", 1)
-
-  use_filtre_fft=.FALSE.
-  CALL getin('use_filtre_fft',use_filtre_fft)
+          "iphysiq must be a multiple of iperiod", 1)
+
+  use_filtre_fft = .FALSE.
+  CALL getin('use_filtre_fft', use_filtre_fft)
   IF (use_filtre_fft) CALL abort_gcm("gcm", 'FFT filter is not available in ' &
           // 'the sequential version of the dynamics.', 1)
 
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Initialisation de XIOS
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
   IF (using_xios) THEN
@@ -152,13 +149,13 @@
 
 
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! FH 2008/05/02
   ! A nettoyer. On ne veut qu'une ou deux routines d'interface 
   ! dynamique -> physique pour l'initialisation
-!#ifdef CPP_PHYS
-!  CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
-!  !      CALL InitComgeomphy ! now done in iniphysiq
-!#endif
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !#ifdef CPP_PHYS
+  !  CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
+  !  !      CALL InitComgeomphy ! now done in iniphysiq
+  !#endif
+  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !-----------------------------------------------------------------------
   !   Choix du calendrier
@@ -168,15 +165,15 @@
 
   IF (calend == 'earth_360d') THEN
-     CALL ioconf_calendar('360_day')
-     WRITE(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
+    CALL ioconf_calendar('360_day')
+    WRITE(lunout, *)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
   ELSE IF (calend == 'earth_365d') THEN
-     CALL ioconf_calendar('noleap')
-     WRITE(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
+    CALL ioconf_calendar('noleap')
+    WRITE(lunout, *)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
   ELSE IF (calend == 'gregorian') THEN
-     CALL ioconf_calendar('gregorian')
-     WRITE(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile'
+    CALL ioconf_calendar('gregorian')
+    WRITE(lunout, *)'CALENDRIER CHOISI: Terrestre bissextile'
   else
-     abort_message = 'Mauvais choix de calendrier'
-     CALL abort_gcm(modname,abort_message,1)
+    abort_message = 'Mauvais choix de calendrier'
+    CALL abort_gcm(modname, abort_message, 1)
   ENDIF
   !-----------------------------------------------------------------------
@@ -196,5 +193,5 @@
 
   ! Allocation de la tableau q : champs advectes   
-  allocate(q(ip1jmp1,llm,nqtot))
+  allocate(q(ip1jmp1, llm, nqtot))
 
   !-----------------------------------------------------------------------
@@ -204,21 +201,21 @@
   !  lecture du fichier start.nc
   IF (read_start) THEN
-     ! we still need to run iniacademic to initialize some
-     ! constants & fields, if we run the 'newtonian' or 'SW' cases:
-     IF (iflag_phys/=1) THEN
-        CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
-     endif
-
-     !        if (planet_type.EQ."earth") THEN
-     ! Load an Earth-format start file
-     CALL dynetat0("start.nc",vcov,ucov, &
-          teta,q,masse,ps,phis, time_0)
-     !        endif ! of if (planet_type.EQ."earth")
-
-     !       WRITE(73,*) 'ucov',ucov
-     !       WRITE(74,*) 'vcov',vcov
-     !       WRITE(75,*) 'teta',teta
-     !       WRITE(76,*) 'ps',ps
-     !       WRITE(77,*) 'q',q
+    ! we still need to run iniacademic to initialize some
+    ! constants & fields, if we run the 'newtonian' or 'SW' cases:
+    IF (iflag_phys/=1) THEN
+      CALL iniacademic(vcov, ucov, teta, q, masse, ps, phis, time_0)
+    endif
+
+    !        if (planet_type.EQ."earth") THEN
+    ! Load an Earth-format start file
+    CALL dynetat0("start.nc", vcov, ucov, &
+            teta, q, masse, ps, phis, time_0)
+    !        endif ! of if (planet_type.EQ."earth")
+
+    !       WRITE(73,*) 'ucov',ucov
+    !       WRITE(74,*) 'vcov',vcov
+    !       WRITE(75,*) 'teta',teta
+    !       WRITE(76,*) 'ps',ps
+    !       WRITE(77,*) 'q',q
 
   ENDIF ! of if (read_start)
@@ -226,10 +223,10 @@
 
   ! le cas echeant, creation d un etat initial
-  IF (prt_level > 9) WRITE(lunout,*) &
-       'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
+  IF (prt_level > 9) WRITE(lunout, *) &
+          'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
   IF (.NOT.read_start) THEN
-     start_time=0.
-     annee_ref=anneeref
-     CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
+    start_time = 0.
+    annee_ref = anneeref
+    CALL iniacademic(vcov, ucov, teta, q, masse, ps, phis, time_0)
   ENDIF
 
@@ -240,20 +237,20 @@
   !  on recalcule eventuellement le pas de temps
 
-  IF(MOD(day_step,iperiod)/=0) THEN
-     abort_message =  &
-          'Il faut choisir un nb de pas par jour multiple de iperiod'
-     CALL abort_gcm(modname,abort_message,1)
-  ENDIF
-
-  IF(MOD(day_step,iphysiq)/=0) THEN
-     abort_message =  &
-          'Il faut choisir un nb de pas par jour multiple de iphysiq'
-     CALL abort_gcm(modname,abort_message,1)
-  ENDIF
-
-  zdtvr    = daysec/REAL(day_step)
+  IF(MOD(day_step, iperiod)/=0) THEN
+    abort_message = &
+            'Il faut choisir un nb de pas par jour multiple de iperiod'
+    CALL abort_gcm(modname, abort_message, 1)
+  ENDIF
+
+  IF(MOD(day_step, iphysiq)/=0) THEN
+    abort_message = &
+            'Il faut choisir un nb de pas par jour multiple de iphysiq'
+    CALL abort_gcm(modname, abort_message, 1)
+  ENDIF
+
+  zdtvr = daysec / REAL(day_step)
   IF(dtvr/=zdtvr) THEN
-     WRITE(lunout,*) &
-          'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
+    WRITE(lunout, *) &
+            'WARNING!!! changement de pas de temps', dtvr, '>', zdtvr
   ENDIF
 
@@ -261,31 +258,31 @@
 
   IF (start_time /= starttime) THEN
-     WRITE(lunout,*)' GCM: Attention l''heure de depart lue dans le' &
-          ,' fichier restart ne correspond pas a celle lue dans le run.def'
-     IF (raz_date == 1) THEN
-        WRITE(lunout,*)'Je prends l''heure lue dans run.def'
-        start_time = starttime
-     ELSE
-        CALL abort_gcm("gcm", "'Je m''arrete'", 1)
-     ENDIF
+    WRITE(lunout, *)' GCM: Attention l''heure de depart lue dans le' &
+            , ' fichier restart ne correspond pas a celle lue dans le run.def'
+    IF (raz_date == 1) THEN
+      WRITE(lunout, *)'Je prends l''heure lue dans run.def'
+      start_time = starttime
+    ELSE
+      CALL abort_gcm("gcm", "'Je m''arrete'", 1)
+    ENDIF
   ENDIF
   IF (raz_date == 1) THEN
-     annee_ref = anneeref
-     day_ref = dayref
-     day_ini = dayref
-     itau_dyn = 0
-     itau_phy = 0
-     time_0 = 0.
-     WRITE(lunout,*) &
-          'GCM: On reinitialise a la date lue dans gcm.def'
+    annee_ref = anneeref
+    day_ref = dayref
+    day_ini = dayref
+    itau_dyn = 0
+    itau_phy = 0
+    time_0 = 0.
+    WRITE(lunout, *) &
+            'GCM: On reinitialise a la date lue dans gcm.def'
   ELSE IF (annee_ref /= anneeref .OR. day_ref /= dayref) THEN
-     WRITE(lunout,*) &
-          'GCM: Attention les dates initiales lues dans le fichier'
-     WRITE(lunout,*) &
-          ' restart ne correspondent pas a celles lues dans '
-     WRITE(lunout,*)' gcm.def'
-     WRITE(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
-     WRITE(lunout,*)' day_ref=',day_ref," dayref=",dayref
-     WRITE(lunout,*)' Pas de remise a zero'
+    WRITE(lunout, *) &
+            'GCM: Attention les dates initiales lues dans le fichier'
+    WRITE(lunout, *) &
+            ' restart ne correspondent pas a celles lues dans '
+    WRITE(lunout, *)' gcm.def'
+    WRITE(lunout, *)' annee_ref=', annee_ref, " anneeref=", anneeref
+    WRITE(lunout, *)' day_ref=', day_ref, " dayref=", dayref
+    WRITE(lunout, *)' Pas de remise a zero'
   ENDIF
 
@@ -323,30 +320,29 @@
   CALL ioconf_startdate(INT(jD_ref), jH_ref)
 
-  WRITE(lunout,*)'DEBUG'
-  WRITE(lunout,*)'annee_ref, mois, day_ref, heure, jD_ref'
-  WRITE(lunout,*)annee_ref, mois, day_ref, heure, jD_ref
-  CALL ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
-  WRITE(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
-  WRITE(lunout,*)jD_ref+jH_ref,an, mois, jour, heure
-
+  WRITE(lunout, *)'DEBUG'
+  WRITE(lunout, *)'annee_ref, mois, day_ref, heure, jD_ref'
+  WRITE(lunout, *)annee_ref, mois, day_ref, heure, jD_ref
+  CALL ju2ymds(jD_ref + jH_ref, an, mois, jour, heure)
+  WRITE(lunout, *)'jD_ref+jH_ref,an, mois, jour, heure'
+  WRITE(lunout, *)jD_ref + jH_ref, an, mois, jour, heure
 
   IF (iflag_phys==1) THEN
-     ! these initialisations have already been done (via iniacademic)
-     ! if running in SW or Newtonian mode
-     !-----------------------------------------------------------------------
-     !   Initialisation des constantes dynamiques :
-     !   ------------------------------------------
-     dtvr = zdtvr
-     CALL iniconst
-
-     !-----------------------------------------------------------------------
-     !   Initialisation de la geometrie :
-     !   --------------------------------
-     CALL inigeom
-
-     !-----------------------------------------------------------------------
-     !   Initialisation du filtre :
-     !   --------------------------
-     CALL inifilr
+    ! these initialisations have already been done (via iniacademic)
+    ! if running in SW or Newtonian mode
+    !-----------------------------------------------------------------------
+    !   Initialisation des constantes dynamiques :
+    !   ------------------------------------------
+    dtvr = zdtvr
+    CALL iniconst
+
+    !-----------------------------------------------------------------------
+    !   Initialisation de la geometrie :
+    !   --------------------------------
+    CALL inigeom
+
+    !-----------------------------------------------------------------------
+    !   Initialisation du filtre :
+    !   --------------------------
+    CALL inifilr
   ENDIF ! of if (iflag_phys.EQ.1)
 
@@ -355,6 +351,6 @@
   !   ----------------------------------
 
-  CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   , &
-       tetagdiv, tetagrot , tetatemp, vert_prof_dissip)
+  CALL inidissip(lstardis, nitergdiv, nitergrot, niterh, &
+          tetagdiv, tetagrot, tetatemp, vert_prof_dissip)
 
   !  numero de stockage pour les fichiers de redemarrage:
@@ -364,19 +360,18 @@
   !   ------------------------
 
-
   IF (nday>=0) THEN
-     day_end = day_ini + nday
+    day_end = day_ini + nday
   else
-     day_end = day_ini - nday/day_step
-  ENDIF
-  WRITE(lunout,300)day_ini,day_end
-300 FORMAT('1'/,15x,'run du jour',i7,2x,'au jour',i7//)
+    day_end = day_ini - nday / day_step
+  ENDIF
+  WRITE(lunout, 300)day_ini, day_end
+  300 FORMAT('1'/, 15x, 'run du jour', i7, 2x, 'au jour', i7//)
 
   CALL ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
-  write (lunout,301)jour, mois, an
+  write (lunout, 301)jour, mois, an
   CALL ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
-  write (lunout,302)jour, mois, an
-301 FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
-302 FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
+  write (lunout, 302)jour, mois, an
+  301 FORMAT('1'/, 15x, 'run du ', i2, '/', i2, '/', i4)
+  302 FORMAT('1'/, 15x, '    au ', i2, '/', i2, '/', i4)
 
   !-----------------------------------------------------------------------
@@ -387,9 +382,9 @@
     ! Physics:
     IF (CPPKEY_PHYS) THEN
-      CALL iniphysiq(iim,jjm,llm, &
-          (jjm-1)*iim+2,comm_lmdz, &
-          daysec,day_ini,dtphys/nsplit_phys, &
-          rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp, &
-          iflag_phys)
+      CALL iniphysiq(iim, jjm, llm, &
+              (jjm - 1) * iim + 2, comm_lmdz, &
+              daysec, day_ini, dtphys / nsplit_phys, &
+              rlatu, rlatv, rlonu, rlonv, aire, cu, cv, rad, g, r, cpp, &
+              iflag_phys)
     END IF
   ENDIF ! of IF ((iflag_phys==1).OR.(iflag_phys>=100))
@@ -405,26 +400,26 @@
   time_step = zdtvr
   IF (ok_dyn_ins) THEN
-     ! initialize output file for instantaneous outputs
-     ! t_ops = iecri * daysec ! do operations every t_ops
-     t_ops =((1.0*iecri)/day_step) * daysec  
-     t_wrt = daysec ! iecri * daysec ! write output every t_wrt
-     CALL inithist(day_ref,annee_ref,time_step, &
-          t_ops,t_wrt)
-  ENDIF
-
-  IF (ok_dyn_ave) THEN 
-     ! initialize output file for averaged outputs
-     t_ops = iperiod * time_step ! do operations every t_ops
-     t_wrt = periodav * daysec   ! write output every t_wrt
-     CALL initdynav(day_ref,annee_ref,time_step, &
-          t_ops,t_wrt)
+    ! initialize output file for instantaneous outputs
+    ! t_ops = iecri * daysec ! do operations every t_ops
+    t_ops = ((1.0 * iecri) / day_step) * daysec
+    t_wrt = daysec ! iecri * daysec ! write output every t_wrt
+    CALL inithist(day_ref, annee_ref, time_step, &
+            t_ops, t_wrt)
+  ENDIF
+
+  IF (ok_dyn_ave) THEN
+    ! initialize output file for averaged outputs
+    t_ops = iperiod * time_step ! do operations every t_ops
+    t_wrt = periodav * daysec   ! write output every t_wrt
+    CALL initdynav(day_ref, annee_ref, time_step, &
+            t_ops, t_wrt)
   END IF
-  dtav = iperiod*dtvr/daysec
+  dtav = iperiod * dtvr / daysec
 
   !  Choix des frequences de stokage pour le offline
   !      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
   !      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
-  istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
-  istphy=istdyn/iphysiq     
+  istdyn = day_step / 4     ! stockage toutes les 6h=1jour/12
+  istphy = istdyn / iphysiq
 
   !-----------------------------------------------------------------------
@@ -438,6 +433,5 @@
   !       WRITE(78,*) 'q',q
 
-
-  CALL leapfrog(ucov,vcov,teta,ps,masse,phis,q,time_0)
+  CALL leapfrog(ucov, vcov, teta, ps, masse, phis, q, time_0)
 
 END PROGRAM gcm
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/guide_mod.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d/guide_mod.F90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d/guide_mod.F90	(revision 5118)
@@ -1,11 +1,10 @@
-
 ! $Id$
 
 MODULE guide_mod
 
-!=======================================================================
-!   Auteur:  F.Hourdin
-!            F. Codron 01/09
-!=======================================================================
+  !=======================================================================
+  !   Auteur:  F.Hourdin
+  !            F. Codron 01/09
+  !=======================================================================
 
   USE getparam, ONLY: ini_getparam, fin_getparam, getpar
@@ -20,53 +19,53 @@
   IMPLICIT NONE
 
-! ---------------------------------------------
-! Declarations des cles logiques et parametres 
-! ---------------------------------------------
-  INTEGER, PRIVATE, SAVE  :: iguide_read,iguide_int,iguide_sav
-  INTEGER, PRIVATE, SAVE  :: nlevnc, guide_plevs
-  LOGICAL, PRIVATE, SAVE  :: guide_u,guide_v,guide_T,guide_Q,guide_P
-  LOGICAL, PRIVATE, SAVE  :: guide_hr,guide_teta  
-  LOGICAL, PRIVATE, SAVE  :: guide_BL,guide_reg,guide_add,gamma4,guide_zon 
-  LOGICAL, PRIVATE, SAVE  :: invert_p,invert_y,ini_anal
-  LOGICAL, PRIVATE, SAVE  :: guide_2D,guide_sav,guide_modele
-!FC
-  LOGICAL, PRIVATE, SAVE  :: convert_Pa
-  
-  REAL, PRIVATE, SAVE     :: tau_min_u,tau_max_u
-  REAL, PRIVATE, SAVE     :: tau_min_v,tau_max_v
-  REAL, PRIVATE, SAVE     :: tau_min_T,tau_max_T
-  REAL, PRIVATE, SAVE     :: tau_min_Q,tau_max_Q
-  REAL, PRIVATE, SAVE     :: tau_min_P,tau_max_P
-
-  REAL, PRIVATE, SAVE     :: lat_min_g,lat_max_g
-  REAL, PRIVATE, SAVE     :: lon_min_g,lon_max_g
-  REAL, PRIVATE, SAVE     :: tau_lon,tau_lat
-
-  REAL, PRIVATE, SAVE     :: plim_guide_BL
-
-  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: alpha_u,alpha_v 
-  REAL, ALLOCATABLE, DIMENSION(:, :), PRIVATE, SAVE     :: alpha_T,alpha_Q 
-  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: alpha_P,alpha_pcor
-  
-! ---------------------------------------------
-! Variables de guidage
-! ---------------------------------------------
-! Variables des fichiers de guidage
-  REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE   :: unat1,unat2
-  REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE   :: vnat1,vnat2
-  REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE   :: tnat1,tnat2
-  REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE   :: qnat1,qnat2
-  REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE   :: pnat1,pnat2
-  REAL, ALLOCATABLE, DIMENSION(:,:),   PRIVATE, SAVE   :: psnat1,psnat2
-  REAL, ALLOCATABLE, DIMENSION(:),     PRIVATE, SAVE   :: apnc,bpnc
-! Variables aux dimensions du modele
-  REAL, ALLOCATABLE, DIMENSION(:,:),   PRIVATE, SAVE   :: ugui1,ugui2
-  REAL, ALLOCATABLE, DIMENSION(:,:),   PRIVATE, SAVE   :: vgui1,vgui2
-  REAL, ALLOCATABLE, DIMENSION(:,:),   PRIVATE, SAVE   :: tgui1,tgui2
-  REAL, ALLOCATABLE, DIMENSION(:,:),   PRIVATE, SAVE   :: qgui1,qgui2
-  REAL, ALLOCATABLE, DIMENSION(:),   PRIVATE, SAVE   :: psgui1,psgui2
+  ! ---------------------------------------------
+  ! Declarations des cles logiques et parametres
+  ! ---------------------------------------------
+  INTEGER, PRIVATE, SAVE :: iguide_read, iguide_int, iguide_sav
+  INTEGER, PRIVATE, SAVE :: nlevnc, guide_plevs
+  LOGICAL, PRIVATE, SAVE :: guide_u, guide_v, guide_T, guide_Q, guide_P
+  LOGICAL, PRIVATE, SAVE :: guide_hr, guide_teta
+  LOGICAL, PRIVATE, SAVE :: guide_BL, guide_reg, guide_add, gamma4, guide_zon
+  LOGICAL, PRIVATE, SAVE :: invert_p, invert_y, ini_anal
+  LOGICAL, PRIVATE, SAVE :: guide_2D, guide_sav, guide_modele
+  !FC
+  LOGICAL, PRIVATE, SAVE :: convert_Pa
+
+  REAL, PRIVATE, SAVE :: tau_min_u, tau_max_u
+  REAL, PRIVATE, SAVE :: tau_min_v, tau_max_v
+  REAL, PRIVATE, SAVE :: tau_min_T, tau_max_T
+  REAL, PRIVATE, SAVE :: tau_min_Q, tau_max_Q
+  REAL, PRIVATE, SAVE :: tau_min_P, tau_max_P
+
+  REAL, PRIVATE, SAVE :: lat_min_g, lat_max_g
+  REAL, PRIVATE, SAVE :: lon_min_g, lon_max_g
+  REAL, PRIVATE, SAVE :: tau_lon, tau_lat
+
+  REAL, PRIVATE, SAVE :: plim_guide_BL
+
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_u, alpha_v
+  REAL, ALLOCATABLE, DIMENSION(:, :), PRIVATE, SAVE :: alpha_T, alpha_Q
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_P, alpha_pcor
+
+  ! ---------------------------------------------
+  ! Variables de guidage
+  ! ---------------------------------------------
+  ! Variables des fichiers de guidage
+  REAL, ALLOCATABLE, DIMENSION(:, :, :), PRIVATE, SAVE :: unat1, unat2
+  REAL, ALLOCATABLE, DIMENSION(:, :, :), PRIVATE, SAVE :: vnat1, vnat2
+  REAL, ALLOCATABLE, DIMENSION(:, :, :), PRIVATE, SAVE :: tnat1, tnat2
+  REAL, ALLOCATABLE, DIMENSION(:, :, :), PRIVATE, SAVE :: qnat1, qnat2
+  REAL, ALLOCATABLE, DIMENSION(:, :, :), PRIVATE, SAVE :: pnat1, pnat2
+  REAL, ALLOCATABLE, DIMENSION(:, :), PRIVATE, SAVE :: psnat1, psnat2
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: apnc, bpnc
+  ! Variables aux dimensions du modele
+  REAL, ALLOCATABLE, DIMENSION(:, :), PRIVATE, SAVE :: ugui1, ugui2
+  REAL, ALLOCATABLE, DIMENSION(:, :), PRIVATE, SAVE :: vgui1, vgui2
+  REAL, ALLOCATABLE, DIMENSION(:, :), PRIVATE, SAVE :: tgui1, tgui2
+  REAL, ALLOCATABLE, DIMENSION(:, :), PRIVATE, SAVE :: qgui1, qgui2
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: psgui1, psgui2
 
 CONTAINS
-!=======================================================================
+  !=======================================================================
 
   SUBROUTINE guide_init
@@ -77,535 +76,534 @@
 
     IMPLICIT NONE
-  
+
     INCLUDE "dimensions.h"
     INCLUDE "paramet.h"
 
-    INTEGER                :: error,ncidpl,rid,rcod
-    CHARACTER (len = 80)   :: abort_message
-    CHARACTER (len = 20)   :: modname = 'guide_init'
-    CHARACTER (len = 20)   :: namedim
-
-! ---------------------------------------------
-! Lecture des parametres:  
-! ---------------------------------------------
+    INTEGER :: error, ncidpl, rid, rcod
+    CHARACTER (len = 80) :: abort_message
+    CHARACTER (len = 20) :: modname = 'guide_init'
+    CHARACTER (len = 20) :: namedim
+
+    ! ---------------------------------------------
+    ! Lecture des parametres:
+    ! ---------------------------------------------
     CALL ini_getparam("nudging_parameters_out.txt")
-! Variables guidees
-    CALL getpar('guide_u',.TRUE.,guide_u,'guidage de u')
-    CALL getpar('guide_v',.TRUE.,guide_v,'guidage de v')
-    CALL getpar('guide_T',.TRUE.,guide_T,'guidage de T')
-    CALL getpar('guide_P',.TRUE.,guide_P,'guidage de P')
-    CALL getpar('guide_Q',.TRUE.,guide_Q,'guidage de Q')
-    CALL getpar('guide_hr',.TRUE.,guide_hr,'guidage de Q par H.R')
-    CALL getpar('guide_teta',.FALSE.,guide_teta,'guidage de T par Teta')
-
-    CALL getpar('guide_add',.FALSE.,guide_add,'forçage constant?')
-    CALL getpar('guide_zon',.FALSE.,guide_zon,'guidage moy zonale')
+    ! Variables guidees
+    CALL getpar('guide_u', .TRUE., guide_u, 'guidage de u')
+    CALL getpar('guide_v', .TRUE., guide_v, 'guidage de v')
+    CALL getpar('guide_T', .TRUE., guide_T, 'guidage de T')
+    CALL getpar('guide_P', .TRUE., guide_P, 'guidage de P')
+    CALL getpar('guide_Q', .TRUE., guide_Q, 'guidage de Q')
+    CALL getpar('guide_hr', .TRUE., guide_hr, 'guidage de Q par H.R')
+    CALL getpar('guide_teta', .FALSE., guide_teta, 'guidage de T par Teta')
+
+    CALL getpar('guide_add', .FALSE., guide_add, 'forçage constant?')
+    CALL getpar('guide_zon', .FALSE., guide_zon, 'guidage moy zonale')
     IF (guide_zon .AND. abs(grossismx - 1.) > 0.01) &
-         CALL abort_gcm("guide_init", &
-         "zonal nudging requires grid regular in longitude", 1)
-
-!   Constantes de rappel. Unite : fraction de jour
-    CALL getpar('tau_min_u',0.02,tau_min_u,'Cste de rappel min, u')
-    CALL getpar('tau_max_u', 10.,tau_max_u,'Cste de rappel max, u')
-    CALL getpar('tau_min_v',0.02,tau_min_v,'Cste de rappel min, v')
-    CALL getpar('tau_max_v', 10.,tau_max_v,'Cste de rappel max, v')
-    CALL getpar('tau_min_T',0.02,tau_min_T,'Cste de rappel min, T')
-    CALL getpar('tau_max_T', 10.,tau_max_T,'Cste de rappel max, T')
-    CALL getpar('tau_min_Q',0.02,tau_min_Q,'Cste de rappel min, Q')
-    CALL getpar('tau_max_Q', 10.,tau_max_Q,'Cste de rappel max, Q')
-    CALL getpar('tau_min_P',0.02,tau_min_P,'Cste de rappel min, P')
-    CALL getpar('tau_max_P', 10.,tau_max_P,'Cste de rappel max, P')
-    CALL getpar('gamma4',.FALSE.,gamma4,'Zone sans rappel elargie')
-    CALL getpar('guide_BL',.TRUE.,guide_BL,'guidage dans C.Lim')
-    CALL getpar('plim_guide_BL',85000.,plim_guide_BL,'BL top presnivs value')
-
-
-! Sauvegarde du forçage
-    CALL getpar('guide_sav',.FALSE.,guide_sav,'sauvegarde guidage')
-    CALL getpar('iguide_sav',4,iguide_sav,'freq. sauvegarde guidage')
+            CALL abort_gcm("guide_init", &
+                    "zonal nudging requires grid regular in longitude", 1)
+
+    !   Constantes de rappel. Unite : fraction de jour
+    CALL getpar('tau_min_u', 0.02, tau_min_u, 'Cste de rappel min, u')
+    CALL getpar('tau_max_u', 10., tau_max_u, 'Cste de rappel max, u')
+    CALL getpar('tau_min_v', 0.02, tau_min_v, 'Cste de rappel min, v')
+    CALL getpar('tau_max_v', 10., tau_max_v, 'Cste de rappel max, v')
+    CALL getpar('tau_min_T', 0.02, tau_min_T, 'Cste de rappel min, T')
+    CALL getpar('tau_max_T', 10., tau_max_T, 'Cste de rappel max, T')
+    CALL getpar('tau_min_Q', 0.02, tau_min_Q, 'Cste de rappel min, Q')
+    CALL getpar('tau_max_Q', 10., tau_max_Q, 'Cste de rappel max, Q')
+    CALL getpar('tau_min_P', 0.02, tau_min_P, 'Cste de rappel min, P')
+    CALL getpar('tau_max_P', 10., tau_max_P, 'Cste de rappel max, P')
+    CALL getpar('gamma4', .FALSE., gamma4, 'Zone sans rappel elargie')
+    CALL getpar('guide_BL', .TRUE., guide_BL, 'guidage dans C.Lim')
+    CALL getpar('plim_guide_BL', 85000., plim_guide_BL, 'BL top presnivs value')
+
+
+    ! Sauvegarde du forçage
+    CALL getpar('guide_sav', .FALSE., guide_sav, 'sauvegarde guidage')
+    CALL getpar('iguide_sav', 4, iguide_sav, 'freq. sauvegarde guidage')
     ! frequences f>0: fx/jour; f<0: tous les f jours; f=0: 1 seule fois.
     IF (iguide_sav>0) THEN
-       iguide_sav=day_step/iguide_sav
+      iguide_sav = day_step / iguide_sav
     ELSE if (iguide_sav == 0) THEN
-       iguide_sav = huge(0)
+      iguide_sav = huge(0)
     ELSE
-       iguide_sav=day_step*iguide_sav
-    ENDIF
-
-! Guidage regional seulement (sinon constant ou suivant le zoom)
-    CALL getpar('guide_reg',.FALSE.,guide_reg,'guidage regional')
-    CALL getpar('lat_min_g',-90.,lat_min_g,'Latitude mini guidage ')
-    CALL getpar('lat_max_g', 90.,lat_max_g,'Latitude maxi guidage ')
-    CALL getpar('lon_min_g',-180.,lon_min_g,'longitude mini guidage ')
-    CALL getpar('lon_max_g', 180.,lon_max_g,'longitude maxi guidage ')
-    CALL getpar('tau_lat', 5.,tau_lat,'raideur lat guide regional ')
-    CALL getpar('tau_lon', 5.,tau_lon,'raideur lon guide regional ')
-
-! Parametres pour lecture des fichiers
-    CALL getpar('iguide_read',4,iguide_read,'freq. lecture guidage')
-    CALL getpar('iguide_int',4,iguide_int,'freq. interpolation vert')
+      iguide_sav = day_step * iguide_sav
+    ENDIF
+
+    ! Guidage regional seulement (sinon constant ou suivant le zoom)
+    CALL getpar('guide_reg', .FALSE., guide_reg, 'guidage regional')
+    CALL getpar('lat_min_g', -90., lat_min_g, 'Latitude mini guidage ')
+    CALL getpar('lat_max_g', 90., lat_max_g, 'Latitude maxi guidage ')
+    CALL getpar('lon_min_g', -180., lon_min_g, 'longitude mini guidage ')
+    CALL getpar('lon_max_g', 180., lon_max_g, 'longitude maxi guidage ')
+    CALL getpar('tau_lat', 5., tau_lat, 'raideur lat guide regional ')
+    CALL getpar('tau_lon', 5., tau_lon, 'raideur lon guide regional ')
+
+    ! Parametres pour lecture des fichiers
+    CALL getpar('iguide_read', 4, iguide_read, 'freq. lecture guidage')
+    CALL getpar('iguide_int', 4, iguide_int, 'freq. interpolation vert')
     IF (iguide_int==0) THEN
-        iguide_int=1
+      iguide_int = 1
     ELSEIF (iguide_int>0) THEN
-        iguide_int=day_step/iguide_int
+      iguide_int = day_step / iguide_int
     ELSE
-        iguide_int=day_step*iguide_int
-    ENDIF
-    CALL getpar('guide_plevs',0,guide_plevs,'niveaux pression fichiers guidage')
+      iguide_int = day_step * iguide_int
+    ENDIF
+    CALL getpar('guide_plevs', 0, guide_plevs, 'niveaux pression fichiers guidage')
     ! Pour compatibilite avec ancienne version avec guide_modele
-    CALL getpar('guide_modele',.FALSE.,guide_modele,'niveaux pression ap+bp*psol')
+    CALL getpar('guide_modele', .FALSE., guide_modele, 'niveaux pression ap+bp*psol')
     IF (guide_modele) THEN
-        guide_plevs=1
-    ENDIF
-!FC
-    CALL getpar('convert_Pa',.TRUE.,convert_Pa,'Convert Pressure levels in Pa')
+      guide_plevs = 1
+    ENDIF
+    !FC
+    CALL getpar('convert_Pa', .TRUE., convert_Pa, 'Convert Pressure levels in Pa')
     ! Fin raccord
-    CALL getpar('ini_anal',.FALSE.,ini_anal,'Etat initial = analyse')
-    CALL getpar('guide_invertp',.TRUE.,invert_p,'niveaux p inverses')
-    CALL getpar('guide_inverty',.TRUE.,invert_y,'inversion N-S')
-    CALL getpar('guide_2D',.FALSE.,guide_2D,'fichier guidage lat-P')
+    CALL getpar('ini_anal', .FALSE., ini_anal, 'Etat initial = analyse')
+    CALL getpar('guide_invertp', .TRUE., invert_p, 'niveaux p inverses')
+    CALL getpar('guide_inverty', .TRUE., invert_y, 'inversion N-S')
+    CALL getpar('guide_2D', .FALSE., guide_2D, 'fichier guidage lat-P')
 
     CALL fin_getparam
-    
-! ---------------------------------------------
-! Determination du nombre de niveaux verticaux
-! des fichiers guidage
-! ---------------------------------------------
-    ncidpl=-99
+
+    ! ---------------------------------------------
+    ! Determination du nombre de niveaux verticaux
+    ! des fichiers guidage
+    ! ---------------------------------------------
+    ncidpl = -99
     IF (guide_plevs==1) THEN
-       IF (ncidpl==-99) THEN
-          rcod=nf90_open('apbp.nc',nf90_nowrite, ncidpl)
-          IF (rcod/=nf90_noerr) THEN
-             abort_message=' Nudging error -> no file apbp.nc'
-             CALL abort_gcm(modname,abort_message,1)
-          endif
-       endif
+      IF (ncidpl==-99) THEN
+        rcod = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
+        IF (rcod/=nf90_noerr) THEN
+          abort_message = ' Nudging error -> no file apbp.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        endif
+      endif
     elseif (guide_plevs==2) THEN
-       IF (ncidpl==-99) THEN
-          rcod=nf90_open('P.nc',nf90_nowrite,ncidpl)
-          IF (rcod/=nf90_noerr) THEN
-             abort_message=' Nudging error -> no file P.nc'
-             CALL abort_gcm(modname,abort_message,1)
-          endif
-       endif
+      IF (ncidpl==-99) THEN
+        rcod = nf90_open('P.nc', nf90_nowrite, ncidpl)
+        IF (rcod/=nf90_noerr) THEN
+          abort_message = ' Nudging error -> no file P.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        endif
+      endif
 
     elseif (guide_u) THEN
-           IF (ncidpl==-99) THEN
-               rcod=nf90_open('u.nc',nf90_nowrite,ncidpl)
-               IF (rcod/=nf90_noerr) THEN
-                  CALL abort_gcm(modname, &
-                       ' Nudging error -> no file u.nc',1)
-               endif
-           endif
+      IF (ncidpl==-99) THEN
+        rcod = nf90_open('u.nc', nf90_nowrite, ncidpl)
+        IF (rcod/=nf90_noerr) THEN
+          CALL abort_gcm(modname, &
+                  ' Nudging error -> no file u.nc', 1)
+        endif
+      endif
 
     elseif (guide_v) THEN
-           IF (ncidpl==-99) THEN
-               rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
-               IF (rcod/=nf90_noerr) THEN
-                  CALL abort_gcm(modname, &
-                       ' Nudging error -> no file v.nc',1)
-               endif
-           endif
+      IF (ncidpl==-99) THEN
+        rcod = nf90_open('v.nc', nf90_nowrite, ncidpl)
+        IF (rcod/=nf90_noerr) THEN
+          CALL abort_gcm(modname, &
+                  ' Nudging error -> no file v.nc', 1)
+        endif
+      endif
     elseif (guide_T) THEN
-           IF (ncidpl==-99) THEN
-               rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
-               IF (rcod/=nf90_noerr) THEN
-                  CALL abort_gcm(modname, &
-                       ' Nudging error -> no file T.nc',1)
-               endif
-           endif
+      IF (ncidpl==-99) THEN
+        rcod = nf90_open('T.nc', nf90_nowrite, ncidpl)
+        IF (rcod/=nf90_noerr) THEN
+          CALL abort_gcm(modname, &
+                  ' Nudging error -> no file T.nc', 1)
+        endif
+      endif
     elseif (guide_Q) THEN
-           IF (ncidpl==-99) THEN
-               rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
-               IF (rcod/=nf90_noerr) THEN
-                  CALL abort_gcm(modname, &
-                       ' Nudging error -> no file hur.nc',1)
-               endif
-           endif
-
-
-    endif 
-    error=nf90_inq_dimid(ncidpl,'LEVEL',rid)
-    IF (error/=nf90_noerr) error=nf90_inq_dimid(ncidpl,'PRESSURE',rid)
+      IF (ncidpl==-99) THEN
+        rcod = nf90_open('hur.nc', nf90_nowrite, ncidpl)
+        IF (rcod/=nf90_noerr) THEN
+          CALL abort_gcm(modname, &
+                  ' Nudging error -> no file hur.nc', 1)
+        endif
+      endif
+
+    endif
+    error = nf90_inq_dimid(ncidpl, 'LEVEL', rid)
+    IF (error/=nf90_noerr) error = nf90_inq_dimid(ncidpl, 'PRESSURE', rid)
     IF (error/=nf90_noerr) THEN
-        CALL abort_gcm(modname,'Nudging: error reading pressure levels',1)
-    ENDIF
-    error=nf90_inquire_dimension(ncidpl,rid,len=nlevnc)
-    WRITE(*,*)trim(modname)//' : number of vertical levels nlevnc', nlevnc
+      CALL abort_gcm(modname, 'Nudging: error reading pressure levels', 1)
+    ENDIF
+    error = nf90_inquire_dimension(ncidpl, rid, len = nlevnc)
+    WRITE(*, *)trim(modname) // ' : number of vertical levels nlevnc', nlevnc
     rcod = nf90_close(ncidpl)
 
-! ---------------------------------------------
-! Allocation des variables
-! ---------------------------------------------
-    abort_message='nudging allocation error'
+    ! ---------------------------------------------
+    ! Allocation des variables
+    ! ---------------------------------------------
+    abort_message = 'nudging allocation error'
 
     ALLOCATE(apnc(nlevnc), stat = error)
-    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
+    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     ALLOCATE(bpnc(nlevnc), stat = error)
-    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
-    apnc=0.;bpnc=0.
+    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
+    apnc = 0.;bpnc = 0.
 
     ALLOCATE(alpha_pcor(llm), stat = error)
-    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
+    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     ALLOCATE(alpha_u(ip1jmp1), stat = error)
-    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
+    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     ALLOCATE(alpha_v(ip1jm), stat = error)
-    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
+    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     ALLOCATE(alpha_T(iip1, jjp1), stat = error)
-    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
+    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     ALLOCATE(alpha_Q(iip1, jjp1), stat = error)
-    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
+    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     ALLOCATE(alpha_P(ip1jmp1), stat = error)
-    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
-    alpha_u=0.;alpha_v=0;alpha_T=0;alpha_Q=0;alpha_P=0
-    
+    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
+    alpha_u = 0.;alpha_v = 0;alpha_T = 0;alpha_Q = 0;alpha_P = 0
+
     IF (guide_u) THEN
-        ALLOCATE(unat1(iip1,jjp1,nlevnc), stat = error)
-        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
-        ALLOCATE(ugui1(ip1jmp1,llm), stat = error)
-        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
-        ALLOCATE(unat2(iip1,jjp1,nlevnc), stat = error)
-        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
-        ALLOCATE(ugui2(ip1jmp1,llm), stat = error)
-        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
-        unat1=0.;unat2=0.;ugui1=0.;ugui2=0.
+      ALLOCATE(unat1(iip1, jjp1, nlevnc), stat = error)
+      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
+      ALLOCATE(ugui1(ip1jmp1, llm), stat = error)
+      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
+      ALLOCATE(unat2(iip1, jjp1, nlevnc), stat = error)
+      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
+      ALLOCATE(ugui2(ip1jmp1, llm), stat = error)
+      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
+      unat1 = 0.;unat2 = 0.;ugui1 = 0.;ugui2 = 0.
     ENDIF
 
     IF (guide_T) THEN
-        ALLOCATE(tnat1(iip1,jjp1,nlevnc), stat = error)
-        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
-        ALLOCATE(tgui1(ip1jmp1,llm), stat = error)
-        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
-        ALLOCATE(tnat2(iip1,jjp1,nlevnc), stat = error)
-        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
-        ALLOCATE(tgui2(ip1jmp1,llm), stat = error)
-        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
-        tnat1=0.;tnat2=0.;tgui1=0.;tgui2=0.
-    ENDIF
-     
+      ALLOCATE(tnat1(iip1, jjp1, nlevnc), stat = error)
+      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
+      ALLOCATE(tgui1(ip1jmp1, llm), stat = error)
+      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
+      ALLOCATE(tnat2(iip1, jjp1, nlevnc), stat = error)
+      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
+      ALLOCATE(tgui2(ip1jmp1, llm), stat = error)
+      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
+      tnat1 = 0.;tnat2 = 0.;tgui1 = 0.;tgui2 = 0.
+    ENDIF
+
     IF (guide_Q) THEN
-        ALLOCATE(qnat1(iip1,jjp1,nlevnc), stat = error)
-        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
-        ALLOCATE(qgui1(ip1jmp1,llm), stat = error)
-        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
-        ALLOCATE(qnat2(iip1,jjp1,nlevnc), stat = error)
-        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
-        ALLOCATE(qgui2(ip1jmp1,llm), stat = error)
-        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
-        qnat1=0.;qnat2=0.;qgui1=0.;qgui2=0.
+      ALLOCATE(qnat1(iip1, jjp1, nlevnc), stat = error)
+      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
+      ALLOCATE(qgui1(ip1jmp1, llm), stat = error)
+      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
+      ALLOCATE(qnat2(iip1, jjp1, nlevnc), stat = error)
+      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
+      ALLOCATE(qgui2(ip1jmp1, llm), stat = error)
+      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
+      qnat1 = 0.;qnat2 = 0.;qgui1 = 0.;qgui2 = 0.
     ENDIF
 
     IF (guide_v) THEN
-        ALLOCATE(vnat1(iip1,jjm,nlevnc), stat = error)
-        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
-        ALLOCATE(vgui1(ip1jm,llm), stat = error)
-        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
-        ALLOCATE(vnat2(iip1,jjm,nlevnc), stat = error)
-        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
-        ALLOCATE(vgui2(ip1jm,llm), stat = error)
-        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
-        vnat1=0.;vnat2=0.;vgui1=0.;vgui2=0.
+      ALLOCATE(vnat1(iip1, jjm, nlevnc), stat = error)
+      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
+      ALLOCATE(vgui1(ip1jm, llm), stat = error)
+      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
+      ALLOCATE(vnat2(iip1, jjm, nlevnc), stat = error)
+      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
+      ALLOCATE(vgui2(ip1jm, llm), stat = error)
+      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
+      vnat1 = 0.;vnat2 = 0.;vgui1 = 0.;vgui2 = 0.
     ENDIF
 
     IF (guide_plevs==2) THEN
-        ALLOCATE(pnat1(iip1,jjp1,nlevnc), stat = error)
-        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
-        ALLOCATE(pnat2(iip1,jjp1,nlevnc), stat = error)
-        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
-        pnat1=0.;pnat2=0.;
+      ALLOCATE(pnat1(iip1, jjp1, nlevnc), stat = error)
+      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
+      ALLOCATE(pnat2(iip1, jjp1, nlevnc), stat = error)
+      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
+      pnat1 = 0.;pnat2 = 0.;
     ENDIF
 
     IF (guide_P.OR.guide_plevs==1) THEN
-        ALLOCATE(psnat1(iip1,jjp1), stat = error)
-        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
-        ALLOCATE(psnat2(iip1,jjp1), stat = error)
-        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
-        psnat1=0.;psnat2=0.;
+      ALLOCATE(psnat1(iip1, jjp1), stat = error)
+      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
+      ALLOCATE(psnat2(iip1, jjp1), stat = error)
+      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
+      psnat1 = 0.;psnat2 = 0.;
     ENDIF
     IF (guide_P) THEN
-        ALLOCATE(psgui2(ip1jmp1), stat = error)
-        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
-        ALLOCATE(psgui1(ip1jmp1), stat = error)
-        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
-        psgui1=0.;psgui2=0.
-    ENDIF
-
-! ---------------------------------------------
-!   Lecture du premier etat de guidage.
-! ---------------------------------------------
+      ALLOCATE(psgui2(ip1jmp1), stat = error)
+      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
+      ALLOCATE(psgui1(ip1jmp1), stat = error)
+      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
+      psgui1 = 0.;psgui2 = 0.
+    ENDIF
+
+    ! ---------------------------------------------
+    !   Lecture du premier etat de guidage.
+    ! ---------------------------------------------
     IF (guide_2D) THEN
-        CALL guide_read2D(1)
+      CALL guide_read2D(1)
     ELSE
-        CALL guide_read(1)
-    ENDIF
-    IF (guide_v) vnat1=vnat2
-    IF (guide_u) unat1=unat2
-    IF (guide_T) tnat1=tnat2
-    IF (guide_Q) qnat1=qnat2
-    IF (guide_plevs==2) pnat1=pnat2
-    IF (guide_P.OR.guide_plevs==1) psnat1=psnat2
+      CALL guide_read(1)
+    ENDIF
+    IF (guide_v) vnat1 = vnat2
+    IF (guide_u) unat1 = unat2
+    IF (guide_T) tnat1 = tnat2
+    IF (guide_Q) qnat1 = qnat2
+    IF (guide_plevs==2) pnat1 = pnat2
+    IF (guide_P.OR.guide_plevs==1) psnat1 = psnat2
 
   END SUBROUTINE guide_init
 
-!=======================================================================
-  SUBROUTINE guide_main(itau,ucov,vcov,teta,q,masse,ps)
+  !=======================================================================
+  SUBROUTINE guide_main(itau, ucov, vcov, teta, q, masse, ps)
 
     USE exner_hyb_m, ONLY: exner_hyb
     USE exner_milieu_m, ONLY: exner_milieu
     USE control_mod, ONLY: day_step, iperiod
-    USE comconst_mod, ONLY: cpp, dtvr, daysec,kappa
+    USE comconst_mod, ONLY: cpp, dtvr, daysec, kappa
     USE comvert_mod, ONLY: ap, bp, preff, presnivs, pressure_exner
- 
+    USE lmdz_iniprint, ONLY: lunout, prt_level
+
     IMPLICIT NONE
-  
+
     INCLUDE "dimensions.h"
     INCLUDE "paramet.h"
-    INCLUDE "iniprint.h"
 
 
     ! Variables entree
-    INTEGER,                       INTENT(IN)    :: itau !pas de temps
-    REAL, DIMENSION (ip1jmp1,llm), INTENT(INOUT) :: ucov,teta,q,masse
-    REAL, DIMENSION (ip1jm,llm),   INTENT(INOUT) :: vcov
-    REAL, DIMENSION (ip1jmp1),     INTENT(INOUT) :: ps
+    INTEGER, INTENT(IN) :: itau !pas de temps
+    REAL, DIMENSION (ip1jmp1, llm), INTENT(INOUT) :: ucov, teta, q, masse
+    REAL, DIMENSION (ip1jm, llm), INTENT(INOUT) :: vcov
+    REAL, DIMENSION (ip1jmp1), INTENT(INOUT) :: ps
 
     ! Variables locales
-    LOGICAL, SAVE :: first=.TRUE.
-    LOGICAL       :: f_out ! sortie guidage
-    REAL, DIMENSION (ip1jmp1,llm) :: f_add ! var aux: champ de guidage
-    REAL :: pk(ip1jmp1,llm) ! Exner at mid-layers
+    LOGICAL, SAVE :: first = .TRUE.
+    LOGICAL :: f_out ! sortie guidage
+    REAL, DIMENSION (ip1jmp1, llm) :: f_add ! var aux: champ de guidage
+    REAL :: pk(ip1jmp1, llm) ! Exner at mid-layers
     REAL :: pks(ip1jmp1) ! Exner at the surface
     REAL :: unskap ! 1./kappa
-    REAL, DIMENSION (ip1jmp1,llmp1) :: p ! Pressure at inter-layers
+    REAL, DIMENSION (ip1jmp1, llmp1) :: p ! Pressure at inter-layers
     ! Compteurs temps:
-    INTEGER, SAVE :: step_rea,count_no_rea,itau_test ! lecture guidage
-    REAL          :: ditau, dday_step
-    REAL          :: tau,reste ! position entre 2 etats de guidage
-    REAL, SAVE    :: factt ! pas de temps en fraction de jour
-    
-    INTEGER       :: l
-    CHARACTER(LEN=20) :: modname="guide_main"
-    CHARACTER (len = 80)   :: abort_message
-
-
-!-----------------------------------------------------------------------
-! Initialisations au premier passage
-!-----------------------------------------------------------------------
+    INTEGER, SAVE :: step_rea, count_no_rea, itau_test ! lecture guidage
+    REAL :: ditau, dday_step
+    REAL :: tau, reste ! position entre 2 etats de guidage
+    REAL, SAVE :: factt ! pas de temps en fraction de jour
+
+    INTEGER :: l
+    CHARACTER(LEN = 20) :: modname = "guide_main"
+    CHARACTER (len = 80) :: abort_message
+
+
+    !-----------------------------------------------------------------------
+    ! Initialisations au premier passage
+    !-----------------------------------------------------------------------
     IF (first) THEN
-        first=.FALSE.
-        CALL guide_init 
-        itau_test=1001
-        step_rea=1
-        count_no_rea=0
-! Calcul des constantes de rappel
-        factt=dtvr*iperiod/daysec 
-        CALL tau2alpha(3,iip1,jjm ,factt,tau_min_v,tau_max_v,alpha_v)
-        CALL tau2alpha(2,iip1,jjp1,factt,tau_min_u,tau_max_u,alpha_u)
-        CALL tau2alpha(1,iip1,jjp1,factt,tau_min_T,tau_max_T,alpha_T)
-        CALL tau2alpha(1,iip1,jjp1,factt,tau_min_P,tau_max_P,alpha_P)
-        CALL tau2alpha(1,iip1,jjp1,factt,tau_min_Q,tau_max_Q,alpha_Q)
-! correction de rappel dans couche limite
-        IF (guide_BL) THEN
-             alpha_pcor(:)=1.
-        else
-            do l=1,llm
-               alpha_pcor(l)=(1.+tanh(((plim_guide_BL-presnivs(l))/preff)/0.05))/2.
-            enddo
-        endif
-! ini_anal: etat initial egal au guidage        
-        IF (ini_anal) THEN
-            CALL guide_interp(ps,teta)
-            IF (guide_u) ucov=ugui2
-            IF (guide_v) vcov=ugui2
-            IF (guide_T) teta=tgui2
-            IF (guide_Q) q=qgui2
-            IF (guide_P) THEN
-                ps=psgui2
-                CALL pression(ip1jmp1,ap,bp,ps,p)
-                CALL massdair(p,masse)
-            ENDIF
-
-        ENDIF
-! Verification structure guidage
-!        IF (guide_u) THEN
-!            CALL writefield('unat',unat1)
-!            CALL writefield('ucov',RESHAPE(ucov,(/iip1,jjp1,llm/)))
-!        ENDIF
-!        IF (guide_T) THEN
-!            CALL writefield('tnat',tnat1)
-!            CALL writefield('teta',RESHAPE(teta,(/iip1,jjp1,llm/)))
-!        ENDIF
+      first = .FALSE.
+      CALL guide_init
+      itau_test = 1001
+      step_rea = 1
+      count_no_rea = 0
+      ! Calcul des constantes de rappel
+      factt = dtvr * iperiod / daysec
+      CALL tau2alpha(3, iip1, jjm, factt, tau_min_v, tau_max_v, alpha_v)
+      CALL tau2alpha(2, iip1, jjp1, factt, tau_min_u, tau_max_u, alpha_u)
+      CALL tau2alpha(1, iip1, jjp1, factt, tau_min_T, tau_max_T, alpha_T)
+      CALL tau2alpha(1, iip1, jjp1, factt, tau_min_P, tau_max_P, alpha_P)
+      CALL tau2alpha(1, iip1, jjp1, factt, tau_min_Q, tau_max_Q, alpha_Q)
+      ! correction de rappel dans couche limite
+      IF (guide_BL) THEN
+        alpha_pcor(:) = 1.
+      else
+        do l = 1, llm
+          alpha_pcor(l) = (1. + tanh(((plim_guide_BL - presnivs(l)) / preff) / 0.05)) / 2.
+        enddo
+      endif
+      ! ini_anal: etat initial egal au guidage
+      IF (ini_anal) THEN
+        CALL guide_interp(ps, teta)
+        IF (guide_u) ucov = ugui2
+        IF (guide_v) vcov = ugui2
+        IF (guide_T) teta = tgui2
+        IF (guide_Q) q = qgui2
+        IF (guide_P) THEN
+          ps = psgui2
+          CALL pression(ip1jmp1, ap, bp, ps, p)
+          CALL massdair(p, masse)
+        ENDIF
+
+      ENDIF
+      ! Verification structure guidage
+      !        IF (guide_u) THEN
+      !            CALL writefield('unat',unat1)
+      !            CALL writefield('ucov',RESHAPE(ucov,(/iip1,jjp1,llm/)))
+      !        ENDIF
+      !        IF (guide_T) THEN
+      !            CALL writefield('tnat',tnat1)
+      !            CALL writefield('teta',RESHAPE(teta,(/iip1,jjp1,llm/)))
+      !        ENDIF
 
     ENDIF !first
 
-!-----------------------------------------------------------------------
-! Lecture des fichiers de guidage ?
-!-----------------------------------------------------------------------
+    !-----------------------------------------------------------------------
+    ! Lecture des fichiers de guidage ?
+    !-----------------------------------------------------------------------
     IF (iguide_read/=0) THEN
-      ditau=real(itau)
-      dday_step=real(day_step)
+      ditau = real(itau)
+      dday_step = real(day_step)
       IF (iguide_read<0) THEN
-          tau=ditau/dday_step/REAL(iguide_read)
+        tau = ditau / dday_step / REAL(iguide_read)
       ELSE
-          tau=REAL(iguide_read)*ditau/dday_step
-      ENDIF
-      reste=tau-AINT(tau)
+        tau = REAL(iguide_read) * ditau / dday_step
+      ENDIF
+      reste = tau - AINT(tau)
       IF (reste==0.) THEN
-          IF (itau_test==itau) THEN
-            WRITE(lunout,*)trim(modname)//' second pass in advreel at itau=',&
-            itau
-              abort_message='stopped'
-              CALL abort_gcm(modname,abort_message,1) 
+        IF (itau_test==itau) THEN
+          WRITE(lunout, *)trim(modname) // ' second pass in advreel at itau=', &
+                  itau
+          abort_message = 'stopped'
+          CALL abort_gcm(modname, abort_message, 1)
+        ELSE
+          IF (guide_v) vnat1 = vnat2
+          IF (guide_u) unat1 = unat2
+          IF (guide_T) tnat1 = tnat2
+          IF (guide_Q) qnat1 = qnat2
+          IF (guide_plevs==2) pnat1 = pnat2
+          IF (guide_P.OR.guide_plevs==1) psnat1 = psnat2
+          step_rea = step_rea + 1
+          itau_test = itau
+          WRITE(*, *)trim(modname) // ' Reading nudging files, step ', &
+                  step_rea, 'after ', count_no_rea, ' skips'
+          IF (guide_2D) THEN
+            CALL guide_read2D(step_rea)
           ELSE
-              IF (guide_v) vnat1=vnat2
-              IF (guide_u) unat1=unat2
-              IF (guide_T) tnat1=tnat2
-              IF (guide_Q) qnat1=qnat2
-              IF (guide_plevs==2) pnat1=pnat2
-              IF (guide_P.OR.guide_plevs==1) psnat1=psnat2
-              step_rea=step_rea+1
-              itau_test=itau
-              WRITE(*,*)trim(modname)//' Reading nudging files, step ',&
-                     step_rea,'after ',count_no_rea,' skips'
-              IF (guide_2D) THEN
-                  CALL guide_read2D(step_rea)
-              ELSE
-                  CALL guide_read(step_rea)
-              ENDIF
-              count_no_rea=0
+            CALL guide_read(step_rea)
           ENDIF
+          count_no_rea = 0
+        ENDIF
       ELSE
-        count_no_rea=count_no_rea+1
+        count_no_rea = count_no_rea + 1
 
       ENDIF
     ENDIF !iguide_read=0
 
-!-----------------------------------------------------------------------
-! Interpolation et conversion des champs de guidage
-!-----------------------------------------------------------------------
-    IF (MOD(itau,iguide_int)==0) THEN
-        CALL guide_interp(ps,teta)
-    ENDIF
-! Repartition entre 2 etats de guidage
+    !-----------------------------------------------------------------------
+    ! Interpolation et conversion des champs de guidage
+    !-----------------------------------------------------------------------
+    IF (MOD(itau, iguide_int)==0) THEN
+      CALL guide_interp(ps, teta)
+    ENDIF
+    ! Repartition entre 2 etats de guidage
     IF (iguide_read/=0) THEN
-        tau=reste
+      tau = reste
     ELSE
-        tau=1.
-    ENDIF
-
-!-----------------------------------------------------------------------
-!   Ajout des champs de guidage 
-!-----------------------------------------------------------------------
-! Sauvegarde du guidage?
-    f_out=((MOD(itau,iguide_sav)==0).AND.guide_sav)
+      tau = 1.
+    ENDIF
+
+    !-----------------------------------------------------------------------
+    !   Ajout des champs de guidage
+    !-----------------------------------------------------------------------
+    ! Sauvegarde du guidage?
+    f_out = ((MOD(itau, iguide_sav)==0).AND.guide_sav)
     IF (f_out) THEN
       ! compute pressures at layer interfaces
-      CALL pression(ip1jmp1,ap,bp,ps,p)
+      CALL pression(ip1jmp1, ap, bp, ps, p)
       IF (pressure_exner) THEN
-        CALL exner_hyb(ip1jmp1,ps,p,pks,pk)
+        CALL exner_hyb(ip1jmp1, ps, p, pks, pk)
       else
-        CALL exner_milieu(ip1jmp1,ps,p,pks,pk)
-      endif
-      unskap=1./kappa
+        CALL exner_milieu(ip1jmp1, ps, p, pks, pk)
+      endif
+      unskap = 1. / kappa
       ! Now compute pressures at mid-layer
-      do l=1,llm
-        p(:,l)=preff*(pk(:,l)/cpp)**unskap
+      do l = 1, llm
+        p(:, l) = preff * (pk(:, l) / cpp)**unskap
       enddo
-      CALL guide_out("SP",jjp1,llm,p(:,1:llm))
-    ENDIF
-    
+      CALL guide_out("SP", jjp1, llm, p(:, 1:llm))
+    ENDIF
+
     IF (guide_u) THEN
-        IF (guide_add) THEN
-           f_add=(1.-tau)*ugui1+tau*ugui2
-        else
-           f_add=(1.-tau)*ugui1+tau*ugui2-ucov
-        endif 
-        IF (guide_zon) CALL guide_zonave(1,jjp1,llm,f_add)
-        CALL guide_addfield(ip1jmp1,llm,f_add,alpha_u)
-        IF (f_out) CALL guide_out("ua",jjp1,llm,(1.-tau)*ugui1+tau*ugui2)
-        IF (f_out) CALL guide_out("u",jjp1,llm,ucov)
-        IF (f_out) CALL guide_out("ucov",jjp1,llm,f_add/factt)
-        ucov=ucov+f_add
+      IF (guide_add) THEN
+        f_add = (1. - tau) * ugui1 + tau * ugui2
+      else
+        f_add = (1. - tau) * ugui1 + tau * ugui2 - ucov
+      endif
+      IF (guide_zon) CALL guide_zonave(1, jjp1, llm, f_add)
+      CALL guide_addfield(ip1jmp1, llm, f_add, alpha_u)
+      IF (f_out) CALL guide_out("ua", jjp1, llm, (1. - tau) * ugui1 + tau * ugui2)
+      IF (f_out) CALL guide_out("u", jjp1, llm, ucov)
+      IF (f_out) CALL guide_out("ucov", jjp1, llm, f_add / factt)
+      ucov = ucov + f_add
     endif
 
     IF (guide_T) THEN
-        IF (guide_add) THEN
-           f_add=(1.-tau)*tgui1+tau*tgui2
-        else
-           f_add=(1.-tau)*tgui1+tau*tgui2-teta
-        endif 
-        IF (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
-        CALL guide_addfield(ip1jmp1,llm,f_add,alpha_T)
-        IF (f_out) CALL guide_out("teta",jjp1,llm,f_add/factt)
-        teta=teta+f_add
+      IF (guide_add) THEN
+        f_add = (1. - tau) * tgui1 + tau * tgui2
+      else
+        f_add = (1. - tau) * tgui1 + tau * tgui2 - teta
+      endif
+      IF (guide_zon) CALL guide_zonave(2, jjp1, llm, f_add)
+      CALL guide_addfield(ip1jmp1, llm, f_add, alpha_T)
+      IF (f_out) CALL guide_out("teta", jjp1, llm, f_add / factt)
+      teta = teta + f_add
     endif
 
     IF (guide_P) THEN
-        IF (guide_add) THEN
-           f_add(1:ip1jmp1,1)=(1.-tau)*psgui1+tau*psgui2
-        else
-           f_add(1:ip1jmp1,1)=(1.-tau)*psgui1+tau*psgui2-ps
-        endif 
-        IF (guide_zon) CALL guide_zonave(2,jjp1,1,f_add(1:ip1jmp1,1))
-        CALL guide_addfield(ip1jmp1,1,f_add(1:ip1jmp1,1),alpha_P)
-!        IF (f_out) CALL guide_out("ps",jjp1,1,f_add(1:ip1jmp1,1)/factt)
-        ps=ps+f_add(1:ip1jmp1,1)
-        CALL pression(ip1jmp1,ap,bp,ps,p)
-        CALL massdair(p,masse)
+      IF (guide_add) THEN
+        f_add(1:ip1jmp1, 1) = (1. - tau) * psgui1 + tau * psgui2
+      else
+        f_add(1:ip1jmp1, 1) = (1. - tau) * psgui1 + tau * psgui2 - ps
+      endif
+      IF (guide_zon) CALL guide_zonave(2, jjp1, 1, f_add(1:ip1jmp1, 1))
+      CALL guide_addfield(ip1jmp1, 1, f_add(1:ip1jmp1, 1), alpha_P)
+      !        IF (f_out) CALL guide_out("ps",jjp1,1,f_add(1:ip1jmp1,1)/factt)
+      ps = ps + f_add(1:ip1jmp1, 1)
+      CALL pression(ip1jmp1, ap, bp, ps, p)
+      CALL massdair(p, masse)
     endif
 
     IF (guide_Q) THEN
-        IF (guide_add) THEN
-           f_add=(1.-tau)*qgui1+tau*qgui2
-        else
-           f_add=(1.-tau)*qgui1+tau*qgui2-q
-        endif 
-        IF (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
-        CALL guide_addfield(ip1jmp1,llm,f_add,alpha_Q)
-        IF (f_out) CALL guide_out("q",jjp1,llm,f_add/factt)
-        q=q+f_add
+      IF (guide_add) THEN
+        f_add = (1. - tau) * qgui1 + tau * qgui2
+      else
+        f_add = (1. - tau) * qgui1 + tau * qgui2 - q
+      endif
+      IF (guide_zon) CALL guide_zonave(2, jjp1, llm, f_add)
+      CALL guide_addfield(ip1jmp1, llm, f_add, alpha_Q)
+      IF (f_out) CALL guide_out("q", jjp1, llm, f_add / factt)
+      q = q + f_add
     endif
 
     IF (guide_v) THEN
-        IF (guide_add) THEN
-           f_add(1:ip1jm,:)=(1.-tau)*vgui1+tau*vgui2
-        else
-           f_add(1:ip1jm,:)=(1.-tau)*vgui1+tau*vgui2-vcov
-        endif 
-        IF (guide_zon) CALL guide_zonave(2,jjm,llm,f_add(1:ip1jm,:))
-        CALL guide_addfield(ip1jm,llm,f_add(1:ip1jm,:),alpha_v)
-        IF (f_out) CALL guide_out("v",jjm,llm,vcov)
-        IF (f_out) CALL guide_out("va",jjm,llm,(1.-tau)*vgui1+tau*vgui2)
-        IF (f_out) CALL guide_out("vcov",jjm,llm,f_add(1:ip1jm,:)/factt)
-        vcov=vcov+f_add(1:ip1jm,:)
+      IF (guide_add) THEN
+        f_add(1:ip1jm, :) = (1. - tau) * vgui1 + tau * vgui2
+      else
+        f_add(1:ip1jm, :) = (1. - tau) * vgui1 + tau * vgui2 - vcov
+      endif
+      IF (guide_zon) CALL guide_zonave(2, jjm, llm, f_add(1:ip1jm, :))
+      CALL guide_addfield(ip1jm, llm, f_add(1:ip1jm, :), alpha_v)
+      IF (f_out) CALL guide_out("v", jjm, llm, vcov)
+      IF (f_out) CALL guide_out("va", jjm, llm, (1. - tau) * vgui1 + tau * vgui2)
+      IF (f_out) CALL guide_out("vcov", jjm, llm, f_add(1:ip1jm, :) / factt)
+      vcov = vcov + f_add(1:ip1jm, :)
     endif
   END SUBROUTINE guide_main
 
-!=======================================================================
-  SUBROUTINE guide_addfield(hsize,vsize,field,alpha)
-! field1=a*field1+alpha*field2
+  !=======================================================================
+  SUBROUTINE guide_addfield(hsize, vsize, field, alpha)
+    ! field1=a*field1+alpha*field2
 
     IMPLICIT NONE
 
     ! input variables
-    INTEGER,                      INTENT(IN)    :: hsize
-    INTEGER,                      INTENT(IN)    :: vsize
-    REAL, DIMENSION(hsize),       INTENT(IN)    :: alpha 
-    REAL, DIMENSION(hsize,vsize), INTENT(INOUT) :: field
+    INTEGER, INTENT(IN) :: hsize
+    INTEGER, INTENT(IN) :: vsize
+    REAL, DIMENSION(hsize), INTENT(IN) :: alpha
+    REAL, DIMENSION(hsize, vsize), INTENT(INOUT) :: field
 
     ! Local variables
     INTEGER :: l
 
-    do l=1,vsize
-        field(:,l)=alpha*field(:,l)*alpha_pcor(l)
+    do l = 1, vsize
+      field(:, l) = alpha * field(:, l) * alpha_pcor(l)
     enddo
 
   END SUBROUTINE guide_addfield
 
-!=======================================================================
-  SUBROUTINE guide_zonave(typ,hsize,vsize,field)
+  !=======================================================================
+  SUBROUTINE guide_zonave(typ, hsize, vsize, field)
 
     USE comconst_mod, ONLY: pi
-    
+
     IMPLICIT NONE
 
@@ -613,316 +611,66 @@
     INCLUDE "paramet.h"
     INCLUDE "comgeom.h"
-    
+
     ! input/output variables
-    INTEGER,                           INTENT(IN)    :: typ
-    INTEGER,                           INTENT(IN)    :: vsize
-    INTEGER,                           INTENT(IN)    :: hsize
-    REAL, DIMENSION(hsize*iip1,vsize), INTENT(INOUT) :: field
+    INTEGER, INTENT(IN) :: typ
+    INTEGER, INTENT(IN) :: vsize
+    INTEGER, INTENT(IN) :: hsize
+    REAL, DIMENSION(hsize * iip1, vsize), INTENT(INOUT) :: field
 
     ! Local variables
-    LOGICAL, SAVE                :: first=.TRUE.
+    LOGICAL, SAVE :: first = .TRUE.
     INTEGER, DIMENSION (2), SAVE :: imin, imax ! averaging domain
-    INTEGER                      :: i,j,l,ij
-    REAL, DIMENSION (iip1)       :: lond       ! longitude in Deg.
-    REAL, DIMENSION (hsize,vsize):: fieldm     ! zon-averaged field
+    INTEGER :: i, j, l, ij
+    REAL, DIMENSION (iip1) :: lond       ! longitude in Deg.
+    REAL, DIMENSION (hsize, vsize) :: fieldm     ! zon-averaged field
 
     IF (first) THEN
-        first=.FALSE.
-!Compute domain for averaging
-        lond=rlonu*180./pi
-        imin(1)=1;imax(1)=iip1;
-        imin(2)=1;imax(2)=iip1;
-        IF (guide_reg) THEN
-            DO i=1,iim
-                IF (lond(i)<lon_min_g) imin(1)=i
-                IF (lond(i)<=lon_max_g) imax(1)=i
-            ENDDO
-            lond=rlonv*180./pi
-            DO i=1,iim
-                IF (lond(i)<lon_min_g) imin(2)=i
-                IF (lond(i)<=lon_max_g) imax(2)=i
-            ENDDO
-        ENDIF
-    ENDIF
-
-    fieldm=0.
-    DO l=1,vsize
-    ! Compute zonal average
-        DO j=1,hsize
-            DO i=imin(typ),imax(typ)
-                ij=(j-1)*iip1+i
-                fieldm(j,l)=fieldm(j,l)+field(ij,l)
-            ENDDO
-        ENDDO 
-        fieldm(:,l)=fieldm(:,l)/REAL(imax(typ)-imin(typ)+1)
-    ! Compute forcing
-        DO j=1,hsize
-            DO i=1,iip1
-                ij=(j-1)*iip1+i
-                field(ij,l)=fieldm(j,l)
-            ENDDO
+      first = .FALSE.
+      !Compute domain for averaging
+      lond = rlonu * 180. / pi
+      imin(1) = 1;imax(1) = iip1;
+      imin(2) = 1;imax(2) = iip1;
+      IF (guide_reg) THEN
+        DO i = 1, iim
+          IF (lond(i)<lon_min_g) imin(1) = i
+          IF (lond(i)<=lon_max_g) imax(1) = i
         ENDDO
+        lond = rlonv * 180. / pi
+        DO i = 1, iim
+          IF (lond(i)<lon_min_g) imin(2) = i
+          IF (lond(i)<=lon_max_g) imax(2) = i
+        ENDDO
+      ENDIF
+    ENDIF
+
+    fieldm = 0.
+    DO l = 1, vsize
+      ! Compute zonal average
+      DO j = 1, hsize
+        DO i = imin(typ), imax(typ)
+          ij = (j - 1) * iip1 + i
+          fieldm(j, l) = fieldm(j, l) + field(ij, l)
+        ENDDO
+      ENDDO
+      fieldm(:, l) = fieldm(:, l) / REAL(imax(typ) - imin(typ) + 1)
+      ! Compute forcing
+      DO j = 1, hsize
+        DO i = 1, iip1
+          ij = (j - 1) * iip1 + i
+          field(ij, l) = fieldm(j, l)
+        ENDDO
+      ENDDO
     ENDDO
 
   END SUBROUTINE guide_zonave
 
-!=======================================================================
-  SUBROUTINE guide_interp(psi,teta)
-  
-  USE exner_hyb_m, ONLY: exner_hyb
-  USE exner_milieu_m, ONLY: exner_milieu
-  USE comconst_mod, ONLY: kappa, cpp
-  USE comvert_mod, ONLY: preff, pressure_exner, bp, ap
-  USE lmdz_q_sat, ONLY: q_sat
-  IMPLICIT NONE
-
-  include "dimensions.h"
-  include "paramet.h"
-  include "comgeom2.h"
-
-  REAL, DIMENSION (iip1,jjp1),     INTENT(IN) :: psi ! Psol gcm
-  REAL, DIMENSION (iip1,jjp1,llm), INTENT(IN) :: teta ! Temp. Pot. gcm
-
-  LOGICAL, SAVE                      :: first=.TRUE.
-  ! Variables pour niveaux pression:
-  REAL, DIMENSION (iip1,jjp1,nlevnc) :: plnc1,plnc2 !niveaux pression guidage
-  REAL, DIMENSION (iip1,jjp1,llm)    :: plunc,plsnc !niveaux pression modele
-  REAL, DIMENSION (iip1,jjm,llm)     :: plvnc       !niveaux pression modele
-  REAL, DIMENSION (iip1,jjp1,llmp1)  :: p           ! pression intercouches 
-  REAL, DIMENSION (iip1,jjp1,llm)    :: pls, pext   ! var intermediaire
-  REAL, DIMENSION (iip1,jjp1,llm)    :: pbarx 
-  REAL, DIMENSION (iip1,jjm,llm)     :: pbary 
-  ! Variables pour fonction Exner (P milieu couche)
-  REAL, DIMENSION (iip1,jjp1,llm)    :: pk
-  REAL, DIMENSION (iip1,jjp1)        :: pks    
-  REAL                               :: prefkap,unskap
-  ! Pression de vapeur saturante
-  REAL, DIMENSION (ip1jmp1,llm)      :: qsat
-  !Variables intermediaires interpolation
-  REAL, DIMENSION (iip1,jjp1,llm)    :: zu1,zu2 
-  REAL, DIMENSION (iip1,jjm,llm)     :: zv1,zv2
-  
-  INTEGER                            :: i,j,l,ij
-  CHARACTER(LEN=20),PARAMETER :: modname="guide_interp"
-  
-    WRITE(*,*)trim(modname)//': interpolate nudging variables'
-! -----------------------------------------------------------------
-! Calcul des niveaux de pression champs guidage
-! -----------------------------------------------------------------
-IF (guide_modele) THEN
-    do i=1,iip1
-        do j=1,jjp1
-            do l=1,nlevnc
-                plnc2(i,j,l)=apnc(l)+bpnc(l)*psnat2(i,j)
-                plnc1(i,j,l)=apnc(l)+bpnc(l)*psnat1(i,j)
-            enddo
-        enddo
-    enddo
-else
-    do i=1,iip1
-        do j=1,jjp1
-            do l=1,nlevnc
-                plnc2(i,j,l)=apnc(l)
-                plnc1(i,j,l)=apnc(l)
-           enddo
-        enddo
-    enddo
-
-END IF
-    IF (first) THEN
-        first=.FALSE.
-        WRITE(*,*)trim(modname)//' : check vertical level order'
-        WRITE(*,*)trim(modname)//' LMDZ :'
-        do l=1,llm
-          WRITE(*,*)trim(modname)//' PL(',l,')=',(ap(l)+ap(l+1))/2. &
-                  +psi(1,jjp1)*(bp(l)+bp(l+1))/2.
-        enddo
-        WRITE(*,*)trim(modname)//' nudging file :'
-        do l=1,nlevnc
-          WRITE(*,*)trim(modname)//' PL(',l,')=',plnc2(1,1,l)
-        enddo
-        WRITE(*,*)trim(modname)//' invert ordering: invert_p=',invert_p
-        IF (guide_u) THEN
-            do l=1,nlevnc
-              WRITE(*,*)trim(modname)//' U(',l,')=',unat2(1,1,l)
-            enddo
-        endif
-        IF (guide_T) THEN
-            do l=1,nlevnc
-              WRITE(*,*)trim(modname)//' T(',l,')=',tnat2(1,1,l)
-            enddo
-        endif
-    endif
-    
-! -----------------------------------------------------------------
-! Calcul niveaux pression modele 
-! -----------------------------------------------------------------
-    CALL pression( ip1jmp1, ap, bp, psi, p )
-    IF (pressure_exner) THEN
-      CALL exner_hyb(ip1jmp1,psi,p,pks,pk)
-    else
-      CALL exner_milieu(ip1jmp1,psi,p,pks,pk)
-    endif
-!    ....  Calcul de pls , pression au milieu des couches ,en Pascals
-    unskap=1./kappa
-    prefkap =  preff  ** kappa
-    DO l = 1, llm
-        DO j=1,jjp1
-            DO i =1, iip1
-                pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
-            ENDDO
-        ENDDO
-    ENDDO
-
-!   calcul des pressions pour les grilles u et v
-    do l=1,llm
-        do j=1,jjp1
-            do i=1,iip1
-                pext(i,j,l)=pls(i,j,l)*aire(i,j)
-            enddo
-        enddo
-    enddo
-    CALL massbar(pext, pbarx, pbary )
-    do l=1,llm
-        do j=1,jjp1
-            do i=1,iip1
-                plunc(i,j,l)=pbarx(i,j,l)/aireu(i,j)
-                plsnc(i,j,l)=pls(i,j,l)
-            enddo
-        enddo
-    enddo
-    do l=1,llm
-        do j=1,jjm
-            do i=1,iip1
-                plvnc(i,j,l)=pbary(i,j,l)/airev(i,j)
-            enddo
-        enddo
-    enddo
-
-! -----------------------------------------------------------------
-! Interpolation champs guidage sur niveaux modele (+inversion N/S)
-! Conversion en variables gcm (ucov, vcov...)
-! -----------------------------------------------------------------
-    IF (guide_P) THEN
-        do j=1,jjp1
-            do i=1,iim
-                ij=(j-1)*iip1+i
-                psgui1(ij)=psnat1(i,j)
-                psgui2(ij)=psnat2(i,j)
-            enddo
-            psgui1(iip1*j)=psnat1(1,j)
-            psgui2(iip1*j)=psnat2(1,j)
-        enddo
-    endif
-
-    IF (guide_u) THEN
-        CALL pres2lev(unat1,zu1,nlevnc,llm,plnc1,plunc,iip1,jjp1,invert_p)
-        CALL pres2lev(unat2,zu2,nlevnc,llm,plnc2,plunc,iip1,jjp1,invert_p)
-        do l=1,llm
-            do j=1,jjp1
-                do i=1,iim
-                    ij=(j-1)*iip1+i
-                    ugui1(ij,l)=zu1(i,j,l)*cu(i,j)
-                    ugui2(ij,l)=zu2(i,j,l)*cu(i,j)
-                enddo
-                ugui1(j*iip1,l)=ugui1((j-1)*iip1+1,l)    
-                ugui2(j*iip1,l)=ugui2((j-1)*iip1+1,l)    
-            enddo
-            do i=1,iip1
-                ugui1(i,l)=0.
-                ugui1(ip1jm+i,l)=0.
-                ugui2(i,l)=0.
-                ugui2(ip1jm+i,l)=0.
-            enddo
-        enddo
-    ENDIF
-    
-    IF (guide_T) THEN
-        CALL pres2lev(tnat1,zu1,nlevnc,llm,plnc1,plsnc,iip1,jjp1,invert_p)
-        CALL pres2lev(tnat2,zu2,nlevnc,llm,plnc2,plsnc,iip1,jjp1,invert_p)
-        do l=1,llm
-            do j=1,jjp1
-                IF (guide_teta) THEN
-                    do i=1,iim
-                        ij=(j-1)*iip1+i
-                        tgui1(ij,l)=zu1(i,j,l)
-                        tgui2(ij,l)=zu2(i,j,l)
-                    enddo
-                ELSE
-                    do i=1,iim
-                        ij=(j-1)*iip1+i
-                        tgui1(ij,l)=zu1(i,j,l)*cpp/pk(i,j,l)
-                        tgui2(ij,l)=zu2(i,j,l)*cpp/pk(i,j,l)
-                    enddo
-                ENDIF
-                tgui1(j*iip1,l)=tgui1((j-1)*iip1+1,l)    
-                tgui2(j*iip1,l)=tgui2((j-1)*iip1+1,l)    
-            enddo
-            do i=1,iip1
-                tgui1(i,l)=tgui1(1,l)
-                tgui1(ip1jm+i,l)=tgui1(ip1jm+1,l) 
-                tgui2(i,l)=tgui2(1,l)
-                tgui2(ip1jm+i,l)=tgui2(ip1jm+1,l) 
-            enddo
-        enddo
-    ENDIF
-
-    IF (guide_v) THEN
-
-        CALL pres2lev(vnat1,zv1,nlevnc,llm,plnc1(:,:jjm,:),plvnc,iip1,jjm,invert_p)
-        CALL pres2lev(vnat2,zv2,nlevnc,llm,plnc2(:,:jjm,:),plvnc,iip1,jjm,invert_p)
-
-        do l=1,llm
-            do j=1,jjm
-                do i=1,iim
-                    ij=(j-1)*iip1+i
-                    vgui1(ij,l)=zv1(i,j,l)*cv(i,j)
-                    vgui2(ij,l)=zv2(i,j,l)*cv(i,j)
-                enddo
-                vgui1(j*iip1,l)=vgui1((j-1)*iip1+1,l)    
-                vgui2(j*iip1,l)=vgui2((j-1)*iip1+1,l)    
-            enddo
-        enddo
-    ENDIF
-    
-    IF (guide_Q) THEN
-        ! On suppose qu'on a la bonne variable dans le fichier de guidage:
-        ! Hum.Rel si guide_hr, Hum.Spec. sinon.
-        CALL pres2lev(qnat1,zu1,nlevnc,llm,plnc1,plsnc,iip1,jjp1,invert_p)
-        CALL pres2lev(qnat2,zu2,nlevnc,llm,plnc2,plsnc,iip1,jjp1,invert_p)
-        do l=1,llm
-            do j=1,jjp1
-                do i=1,iim
-                    ij=(j-1)*iip1+i
-                    qgui1(ij,l)=zu1(i,j,l)
-                    qgui2(ij,l)=zu2(i,j,l)
-                enddo
-                qgui1(j*iip1,l)=qgui1((j-1)*iip1+1,l)    
-                qgui2(j*iip1,l)=qgui2((j-1)*iip1+1,l)    
-            enddo
-            do i=1,iip1
-                qgui1(i,l)=qgui1(1,l)
-                qgui1(ip1jm+i,l)=qgui1(ip1jm+1,l) 
-                qgui2(i,l)=qgui2(1,l)
-                qgui2(ip1jm+i,l)=qgui2(ip1jm+1,l) 
-            enddo
-        enddo
-        IF (guide_hr) THEN
-            CALL q_sat(iip1*jjp1*llm,teta*pk/cpp,plsnc,qsat)
-            qgui1=qgui1*qsat*0.01 !hum. rel. en %
-            qgui2=qgui2*qsat*0.01 
-        ENDIF
-    ENDIF
-
-  END SUBROUTINE guide_interp
-
-!=======================================================================
-  SUBROUTINE tau2alpha(typ,pim,pjm,factt,taumin,taumax,alpha)
-
-! Calcul des constantes de rappel alpha (=1/tau)
-
-    USE comconst_mod, ONLY: pi
-    USE serre_mod, ONLY: clon, clat, grossismx, grossismy
-    
+  !=======================================================================
+  SUBROUTINE guide_interp(psi, teta)
+
+    USE exner_hyb_m, ONLY: exner_hyb
+    USE exner_milieu_m, ONLY: exner_milieu
+    USE comconst_mod, ONLY: kappa, cpp
+    USE comvert_mod, ONLY: preff, pressure_exner, bp, ap
+    USE lmdz_q_sat, ONLY: q_sat
     IMPLICIT NONE
 
@@ -931,157 +679,406 @@
     include "comgeom2.h"
 
-! input arguments :
+    REAL, DIMENSION (iip1, jjp1), INTENT(IN) :: psi ! Psol gcm
+    REAL, DIMENSION (iip1, jjp1, llm), INTENT(IN) :: teta ! Temp. Pot. gcm
+
+    LOGICAL, SAVE :: first = .TRUE.
+    ! Variables pour niveaux pression:
+    REAL, DIMENSION (iip1, jjp1, nlevnc) :: plnc1, plnc2 !niveaux pression guidage
+    REAL, DIMENSION (iip1, jjp1, llm) :: plunc, plsnc !niveaux pression modele
+    REAL, DIMENSION (iip1, jjm, llm) :: plvnc       !niveaux pression modele
+    REAL, DIMENSION (iip1, jjp1, llmp1) :: p           ! pression intercouches
+    REAL, DIMENSION (iip1, jjp1, llm) :: pls, pext   ! var intermediaire
+    REAL, DIMENSION (iip1, jjp1, llm) :: pbarx
+    REAL, DIMENSION (iip1, jjm, llm) :: pbary
+    ! Variables pour fonction Exner (P milieu couche)
+    REAL, DIMENSION (iip1, jjp1, llm) :: pk
+    REAL, DIMENSION (iip1, jjp1) :: pks
+    REAL :: prefkap, unskap
+    ! Pression de vapeur saturante
+    REAL, DIMENSION (ip1jmp1, llm) :: qsat
+    !Variables intermediaires interpolation
+    REAL, DIMENSION (iip1, jjp1, llm) :: zu1, zu2
+    REAL, DIMENSION (iip1, jjm, llm) :: zv1, zv2
+
+    INTEGER :: i, j, l, ij
+    CHARACTER(LEN = 20), PARAMETER :: modname = "guide_interp"
+
+    WRITE(*, *)trim(modname) // ': interpolate nudging variables'
+    ! -----------------------------------------------------------------
+    ! Calcul des niveaux de pression champs guidage
+    ! -----------------------------------------------------------------
+    IF (guide_modele) THEN
+      do i = 1, iip1
+        do j = 1, jjp1
+          do l = 1, nlevnc
+            plnc2(i, j, l) = apnc(l) + bpnc(l) * psnat2(i, j)
+            plnc1(i, j, l) = apnc(l) + bpnc(l) * psnat1(i, j)
+          enddo
+        enddo
+      enddo
+    else
+      do i = 1, iip1
+        do j = 1, jjp1
+          do l = 1, nlevnc
+            plnc2(i, j, l) = apnc(l)
+            plnc1(i, j, l) = apnc(l)
+          enddo
+        enddo
+      enddo
+
+    END IF
+    IF (first) THEN
+      first = .FALSE.
+      WRITE(*, *)trim(modname) // ' : check vertical level order'
+      WRITE(*, *)trim(modname) // ' LMDZ :'
+      do l = 1, llm
+        WRITE(*, *)trim(modname) // ' PL(', l, ')=', (ap(l) + ap(l + 1)) / 2. &
+                + psi(1, jjp1) * (bp(l) + bp(l + 1)) / 2.
+      enddo
+      WRITE(*, *)trim(modname) // ' nudging file :'
+      do l = 1, nlevnc
+        WRITE(*, *)trim(modname) // ' PL(', l, ')=', plnc2(1, 1, l)
+      enddo
+      WRITE(*, *)trim(modname) // ' invert ordering: invert_p=', invert_p
+      IF (guide_u) THEN
+        do l = 1, nlevnc
+          WRITE(*, *)trim(modname) // ' U(', l, ')=', unat2(1, 1, l)
+        enddo
+      endif
+      IF (guide_T) THEN
+        do l = 1, nlevnc
+          WRITE(*, *)trim(modname) // ' T(', l, ')=', tnat2(1, 1, l)
+        enddo
+      endif
+    endif
+
+    ! -----------------------------------------------------------------
+    ! Calcul niveaux pression modele
+    ! -----------------------------------------------------------------
+    CALL pression(ip1jmp1, ap, bp, psi, p)
+    IF (pressure_exner) THEN
+      CALL exner_hyb(ip1jmp1, psi, p, pks, pk)
+    else
+      CALL exner_milieu(ip1jmp1, psi, p, pks, pk)
+    endif
+    !    ....  Calcul de pls , pression au milieu des couches ,en Pascals
+    unskap = 1. / kappa
+    prefkap = preff  ** kappa
+    DO l = 1, llm
+      DO j = 1, jjp1
+        DO i = 1, iip1
+          pls(i, j, l) = preff * (pk(i, j, l) / cpp) ** unskap
+        ENDDO
+      ENDDO
+    ENDDO
+
+    !   calcul des pressions pour les grilles u et v
+    do l = 1, llm
+      do j = 1, jjp1
+        do i = 1, iip1
+          pext(i, j, l) = pls(i, j, l) * aire(i, j)
+        enddo
+      enddo
+    enddo
+    CALL massbar(pext, pbarx, pbary)
+    do l = 1, llm
+      do j = 1, jjp1
+        do i = 1, iip1
+          plunc(i, j, l) = pbarx(i, j, l) / aireu(i, j)
+          plsnc(i, j, l) = pls(i, j, l)
+        enddo
+      enddo
+    enddo
+    do l = 1, llm
+      do j = 1, jjm
+        do i = 1, iip1
+          plvnc(i, j, l) = pbary(i, j, l) / airev(i, j)
+        enddo
+      enddo
+    enddo
+
+    ! -----------------------------------------------------------------
+    ! Interpolation champs guidage sur niveaux modele (+inversion N/S)
+    ! Conversion en variables gcm (ucov, vcov...)
+    ! -----------------------------------------------------------------
+    IF (guide_P) THEN
+      do j = 1, jjp1
+        do i = 1, iim
+          ij = (j - 1) * iip1 + i
+          psgui1(ij) = psnat1(i, j)
+          psgui2(ij) = psnat2(i, j)
+        enddo
+        psgui1(iip1 * j) = psnat1(1, j)
+        psgui2(iip1 * j) = psnat2(1, j)
+      enddo
+    endif
+
+    IF (guide_u) THEN
+      CALL pres2lev(unat1, zu1, nlevnc, llm, plnc1, plunc, iip1, jjp1, invert_p)
+      CALL pres2lev(unat2, zu2, nlevnc, llm, plnc2, plunc, iip1, jjp1, invert_p)
+      do l = 1, llm
+        do j = 1, jjp1
+          do i = 1, iim
+            ij = (j - 1) * iip1 + i
+            ugui1(ij, l) = zu1(i, j, l) * cu(i, j)
+            ugui2(ij, l) = zu2(i, j, l) * cu(i, j)
+          enddo
+          ugui1(j * iip1, l) = ugui1((j - 1) * iip1 + 1, l)
+          ugui2(j * iip1, l) = ugui2((j - 1) * iip1 + 1, l)
+        enddo
+        do i = 1, iip1
+          ugui1(i, l) = 0.
+          ugui1(ip1jm + i, l) = 0.
+          ugui2(i, l) = 0.
+          ugui2(ip1jm + i, l) = 0.
+        enddo
+      enddo
+    ENDIF
+
+    IF (guide_T) THEN
+      CALL pres2lev(tnat1, zu1, nlevnc, llm, plnc1, plsnc, iip1, jjp1, invert_p)
+      CALL pres2lev(tnat2, zu2, nlevnc, llm, plnc2, plsnc, iip1, jjp1, invert_p)
+      do l = 1, llm
+        do j = 1, jjp1
+          IF (guide_teta) THEN
+            do i = 1, iim
+              ij = (j - 1) * iip1 + i
+              tgui1(ij, l) = zu1(i, j, l)
+              tgui2(ij, l) = zu2(i, j, l)
+            enddo
+          ELSE
+            do i = 1, iim
+              ij = (j - 1) * iip1 + i
+              tgui1(ij, l) = zu1(i, j, l) * cpp / pk(i, j, l)
+              tgui2(ij, l) = zu2(i, j, l) * cpp / pk(i, j, l)
+            enddo
+          ENDIF
+          tgui1(j * iip1, l) = tgui1((j - 1) * iip1 + 1, l)
+          tgui2(j * iip1, l) = tgui2((j - 1) * iip1 + 1, l)
+        enddo
+        do i = 1, iip1
+          tgui1(i, l) = tgui1(1, l)
+          tgui1(ip1jm + i, l) = tgui1(ip1jm + 1, l)
+          tgui2(i, l) = tgui2(1, l)
+          tgui2(ip1jm + i, l) = tgui2(ip1jm + 1, l)
+        enddo
+      enddo
+    ENDIF
+
+    IF (guide_v) THEN
+
+      CALL pres2lev(vnat1, zv1, nlevnc, llm, plnc1(:, :jjm, :), plvnc, iip1, jjm, invert_p)
+      CALL pres2lev(vnat2, zv2, nlevnc, llm, plnc2(:, :jjm, :), plvnc, iip1, jjm, invert_p)
+
+      do l = 1, llm
+        do j = 1, jjm
+          do i = 1, iim
+            ij = (j - 1) * iip1 + i
+            vgui1(ij, l) = zv1(i, j, l) * cv(i, j)
+            vgui2(ij, l) = zv2(i, j, l) * cv(i, j)
+          enddo
+          vgui1(j * iip1, l) = vgui1((j - 1) * iip1 + 1, l)
+          vgui2(j * iip1, l) = vgui2((j - 1) * iip1 + 1, l)
+        enddo
+      enddo
+    ENDIF
+
+    IF (guide_Q) THEN
+      ! On suppose qu'on a la bonne variable dans le fichier de guidage:
+      ! Hum.Rel si guide_hr, Hum.Spec. sinon.
+      CALL pres2lev(qnat1, zu1, nlevnc, llm, plnc1, plsnc, iip1, jjp1, invert_p)
+      CALL pres2lev(qnat2, zu2, nlevnc, llm, plnc2, plsnc, iip1, jjp1, invert_p)
+      do l = 1, llm
+        do j = 1, jjp1
+          do i = 1, iim
+            ij = (j - 1) * iip1 + i
+            qgui1(ij, l) = zu1(i, j, l)
+            qgui2(ij, l) = zu2(i, j, l)
+          enddo
+          qgui1(j * iip1, l) = qgui1((j - 1) * iip1 + 1, l)
+          qgui2(j * iip1, l) = qgui2((j - 1) * iip1 + 1, l)
+        enddo
+        do i = 1, iip1
+          qgui1(i, l) = qgui1(1, l)
+          qgui1(ip1jm + i, l) = qgui1(ip1jm + 1, l)
+          qgui2(i, l) = qgui2(1, l)
+          qgui2(ip1jm + i, l) = qgui2(ip1jm + 1, l)
+        enddo
+      enddo
+      IF (guide_hr) THEN
+        CALL q_sat(iip1 * jjp1 * llm, teta * pk / cpp, plsnc, qsat)
+        qgui1 = qgui1 * qsat * 0.01 !hum. rel. en %
+        qgui2 = qgui2 * qsat * 0.01
+      ENDIF
+    ENDIF
+
+  END SUBROUTINE guide_interp
+
+  !=======================================================================
+  SUBROUTINE tau2alpha(typ, pim, pjm, factt, taumin, taumax, alpha)
+
+    ! Calcul des constantes de rappel alpha (=1/tau)
+
+    USE comconst_mod, ONLY: pi
+    USE serre_mod, ONLY: clon, clat, grossismx, grossismy
+
+    IMPLICIT NONE
+
+    include "dimensions.h"
+    include "paramet.h"
+    include "comgeom2.h"
+
+    ! input arguments :
     INTEGER, INTENT(IN) :: typ    ! u(2),v(3), ou scalaire(1)
-    INTEGER, INTENT(IN) :: pim,pjm ! dimensions en lat, lon
-    REAL, INTENT(IN)    :: factt   ! pas de temps en fraction de jour
-    REAL, INTENT(IN)    :: taumin,taumax
-! output arguments:
-    REAL, DIMENSION(pim,pjm), INTENT(OUT) :: alpha 
-  
-!  local variables:
-    LOGICAL, SAVE               :: first=.TRUE.
-    REAL, SAVE                  :: gamma,dxdy_min,dxdy_max
-    REAL, DIMENSION (iip1,jjp1) :: zdx,zdy
-    REAL, DIMENSION (iip1,jjp1) :: dxdys,dxdyu
-    REAL, DIMENSION (iip1,jjm)  :: dxdyv
+    INTEGER, INTENT(IN) :: pim, pjm ! dimensions en lat, lon
+    REAL, INTENT(IN) :: factt   ! pas de temps en fraction de jour
+    REAL, INTENT(IN) :: taumin, taumax
+    ! output arguments:
+    REAL, DIMENSION(pim, pjm), INTENT(OUT) :: alpha
+
+    !  local variables:
+    LOGICAL, SAVE :: first = .TRUE.
+    REAL, SAVE :: gamma, dxdy_min, dxdy_max
+    REAL, DIMENSION (iip1, jjp1) :: zdx, zdy
+    REAL, DIMENSION (iip1, jjp1) :: dxdys, dxdyu
+    REAL, DIMENSION (iip1, jjm) :: dxdyv
     REAL dxdy_
-    REAL zlat,zlon
-    REAL alphamin,alphamax,xi
-    INTEGER i,j,ilon,ilat
-    CHARACTER(LEN=20),parameter :: modname="tau2alpha"
-    CHARACTER (len = 80)   :: abort_message
-
-
-    alphamin=factt/taumax
-    alphamax=factt/taumin
+    REAL zlat, zlon
+    REAL alphamin, alphamax, xi
+    INTEGER i, j, ilon, ilat
+    CHARACTER(LEN = 20), parameter :: modname = "tau2alpha"
+    CHARACTER (len = 80) :: abort_message
+
+    alphamin = factt / taumax
+    alphamax = factt / taumin
     IF (guide_reg.OR.guide_add) THEN
-        alpha=alphamax
-!-----------------------------------------------------------------------
-! guide_reg: alpha=alpha_min dans region, 0. sinon.
-!-----------------------------------------------------------------------
-        IF (guide_reg) THEN
-            do j=1,pjm
-                do i=1,pim
-                    IF (typ==2) THEN
-                       zlat=rlatu(j)*180./pi
-                       zlon=rlonu(i)*180./pi
-                    elseif (typ==1) THEN
-                       zlat=rlatu(j)*180./pi
-                       zlon=rlonv(i)*180./pi
-                    elseif (typ==3) THEN
-                       zlat=rlatv(j)*180./pi
-                       zlon=rlonv(i)*180./pi
-                    endif
-                    alpha(i,j)=alphamax/16.* &
-                              (1.+tanh((zlat-lat_min_g)/tau_lat))* &
-                              (1.+tanh((lat_max_g-zlat)/tau_lat))* &
-                              (1.+tanh((zlon-lon_min_g)/tau_lon))* &
-                              (1.+tanh((lon_max_g-zlon)/tau_lon))
-                enddo
-            enddo
-        ENDIF
+      alpha = alphamax
+      !-----------------------------------------------------------------------
+      ! guide_reg: alpha=alpha_min dans region, 0. sinon.
+      !-----------------------------------------------------------------------
+      IF (guide_reg) THEN
+        do j = 1, pjm
+          do i = 1, pim
+            IF (typ==2) THEN
+              zlat = rlatu(j) * 180. / pi
+              zlon = rlonu(i) * 180. / pi
+            elseif (typ==1) THEN
+              zlat = rlatu(j) * 180. / pi
+              zlon = rlonv(i) * 180. / pi
+            elseif (typ==3) THEN
+              zlat = rlatv(j) * 180. / pi
+              zlon = rlonv(i) * 180. / pi
+            endif
+            alpha(i, j) = alphamax / 16. * &
+                    (1. + tanh((zlat - lat_min_g) / tau_lat)) * &
+                    (1. + tanh((lat_max_g - zlat) / tau_lat)) * &
+                    (1. + tanh((zlon - lon_min_g) / tau_lon)) * &
+                    (1. + tanh((lon_max_g - zlon) / tau_lon))
+          enddo
+        enddo
+      ENDIF
     ELSE
-!-----------------------------------------------------------------------
-! Sinon, alpha varie entre alpha_min et alpha_max suivant le zoom.
-!-----------------------------------------------------------------------
-!Calcul de l'aire des mailles
-        do j=2,jjm
-            do i=2,iip1
-               zdx(i,j)=0.5*(cu(i-1,j)+cu(i,j))/cos(rlatu(j))
-            enddo
-            zdx(1,j)=zdx(iip1,j)
-        enddo
-        do j=2,jjm
-            do i=1,iip1
-               zdy(i,j)=0.5*(cv(i,j-1)+cv(i,j))
-            enddo
-        enddo
-        do i=1,iip1
-            zdx(i,1)=zdx(i,2)
-            zdx(i,jjp1)=zdx(i,jjm)
-            zdy(i,1)=zdy(i,2)
-            zdy(i,jjp1)=zdy(i,jjm)
-        enddo
-        do j=1,jjp1
-            do i=1,iip1
-               dxdys(i,j)=sqrt(zdx(i,j)*zdx(i,j)+zdy(i,j)*zdy(i,j))
-            enddo
-        enddo
-        IF (typ==2) THEN
-            do j=1,jjp1
-                do i=1,iim
-                   dxdyu(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
-                enddo
-                dxdyu(iip1,j)=dxdyu(1,j)
-            enddo
-        ENDIF
-        IF (typ==3) THEN
-            do j=1,jjm
-                do i=1,iip1
-                   dxdyv(i,j)=0.5*(dxdys(i,j)+dxdys(i,j+1))
-                enddo
-            enddo
-        ENDIF
-! Premier appel: calcul des aires min et max et de gamma.
-        IF (first) THEN 
-            first=.FALSE.
-            ! coordonnees du centre du zoom
-            CALL coordij(clon,clat,ilon,ilat) 
-            ! aire de la maille au centre du zoom
-            dxdy_min=dxdys(ilon,ilat)
-            ! dxdy maximale de la maille
-            dxdy_max=0.
-            do j=1,jjp1
-                do i=1,iip1
-                     dxdy_max=max(dxdy_max,dxdys(i,j))
-                enddo
-            enddo
-            ! Calcul de gamma
-            IF (abs(grossismx-1.)<0.1.OR.abs(grossismy-1.)<0.1) THEN
-              WRITE(*,*)trim(modname)//' ATTENTION modele peu zoome'
-              WRITE(*,*)trim(modname)//' ATTENTION on prend une constante de guidage cste'
-              gamma=0.
+      !-----------------------------------------------------------------------
+      ! Sinon, alpha varie entre alpha_min et alpha_max suivant le zoom.
+      !-----------------------------------------------------------------------
+      !Calcul de l'aire des mailles
+      do j = 2, jjm
+        do i = 2, iip1
+          zdx(i, j) = 0.5 * (cu(i - 1, j) + cu(i, j)) / cos(rlatu(j))
+        enddo
+        zdx(1, j) = zdx(iip1, j)
+      enddo
+      do j = 2, jjm
+        do i = 1, iip1
+          zdy(i, j) = 0.5 * (cv(i, j - 1) + cv(i, j))
+        enddo
+      enddo
+      do i = 1, iip1
+        zdx(i, 1) = zdx(i, 2)
+        zdx(i, jjp1) = zdx(i, jjm)
+        zdy(i, 1) = zdy(i, 2)
+        zdy(i, jjp1) = zdy(i, jjm)
+      enddo
+      do j = 1, jjp1
+        do i = 1, iip1
+          dxdys(i, j) = sqrt(zdx(i, j) * zdx(i, j) + zdy(i, j) * zdy(i, j))
+        enddo
+      enddo
+      IF (typ==2) THEN
+        do j = 1, jjp1
+          do i = 1, iim
+            dxdyu(i, j) = 0.5 * (dxdys(i, j) + dxdys(i + 1, j))
+          enddo
+          dxdyu(iip1, j) = dxdyu(1, j)
+        enddo
+      ENDIF
+      IF (typ==3) THEN
+        do j = 1, jjm
+          do i = 1, iip1
+            dxdyv(i, j) = 0.5 * (dxdys(i, j) + dxdys(i, j + 1))
+          enddo
+        enddo
+      ENDIF
+      ! Premier appel: calcul des aires min et max et de gamma.
+      IF (first) THEN
+        first = .FALSE.
+        ! coordonnees du centre du zoom
+        CALL coordij(clon, clat, ilon, ilat)
+        ! aire de la maille au centre du zoom
+        dxdy_min = dxdys(ilon, ilat)
+        ! dxdy maximale de la maille
+        dxdy_max = 0.
+        do j = 1, jjp1
+          do i = 1, iip1
+            dxdy_max = max(dxdy_max, dxdys(i, j))
+          enddo
+        enddo
+        ! Calcul de gamma
+        IF (abs(grossismx - 1.)<0.1.OR.abs(grossismy - 1.)<0.1) THEN
+          WRITE(*, *)trim(modname) // ' ATTENTION modele peu zoome'
+          WRITE(*, *)trim(modname) // ' ATTENTION on prend une constante de guidage cste'
+          gamma = 0.
+        else
+          gamma = (dxdy_max - 2. * dxdy_min) / (dxdy_max - dxdy_min)
+          WRITE(*, *)trim(modname) // ' gamma=', gamma
+          IF (gamma<1.e-5) THEN
+            WRITE(*, *)trim(modname) // ' gamma =', gamma, '<1e-5'
+            abort_message = 'stopped'
+            CALL abort_gcm(modname, abort_message, 1)
+          endif
+          gamma = log(0.5) / log(gamma)
+          IF (gamma4) THEN
+            gamma = min(gamma, 4.)
+          endif
+          WRITE(*, *)trim(modname) // ' gamma=', gamma
+        endif
+      ENDIF !first
+
+      do j = 1, pjm
+        do i = 1, pim
+          IF (typ==1) THEN
+            dxdy_ = dxdys(i, j)
+            zlat = rlatu(j) * 180. / pi
+          elseif (typ==2) THEN
+            dxdy_ = dxdyu(i, j)
+            zlat = rlatu(j) * 180. / pi
+          elseif (typ==3) THEN
+            dxdy_ = dxdyv(i, j)
+            zlat = rlatv(j) * 180. / pi
+          endif
+          IF (abs(grossismx - 1.)<0.1.OR.abs(grossismy - 1.)<0.1) THEN
+            ! pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
+            alpha(i, j) = alphamin
+          else
+            xi = ((dxdy_max - dxdy_) / (dxdy_max - dxdy_min))**gamma
+            xi = min(xi, 1.)
+            IF(lat_min_g<=zlat .AND. zlat<=lat_max_g) THEN
+              alpha(i, j) = xi * alphamin + (1. - xi) * alphamax
             else
-              gamma=(dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)
-              WRITE(*,*)trim(modname)//' gamma=',gamma
-              IF (gamma<1.e-5) THEN
-                WRITE(*,*)trim(modname)//' gamma =',gamma,'<1e-5'
-                abort_message='stopped'
-                CALL abort_gcm(modname,abort_message,1)
-              endif
-              gamma=log(0.5)/log(gamma)
-              IF (gamma4) THEN
-                gamma=min(gamma,4.)
-              endif
-              WRITE(*,*)trim(modname)//' gamma=',gamma
+              alpha(i, j) = 0.
             endif
-        ENDIF !first
-
-        do j=1,pjm
-            do i=1,pim
-                IF (typ==1) THEN
-                   dxdy_=dxdys(i,j)
-                   zlat=rlatu(j)*180./pi
-                elseif (typ==2) THEN
-                   dxdy_=dxdyu(i,j)
-                   zlat=rlatu(j)*180./pi
-                elseif (typ==3) THEN
-                   dxdy_=dxdyv(i,j)
-                   zlat=rlatv(j)*180./pi
-                endif
-                IF (abs(grossismx-1.)<0.1.OR.abs(grossismy-1.)<0.1) THEN
-                ! pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
-                    alpha(i,j)=alphamin
-                else
-                    xi=((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma
-                    xi=min(xi,1.)
-                    IF(lat_min_g<=zlat .AND. zlat<=lat_max_g) THEN
-                        alpha(i,j)=xi*alphamin+(1.-xi)*alphamax
-                    else
-                        alpha(i,j)=0.
-                    endif
-                endif
-            enddo
-        enddo
+          endif
+        enddo
+      enddo
     ENDIF ! guide_reg
 
@@ -1090,5 +1087,5 @@
   END SUBROUTINE tau2alpha
 
-!=======================================================================
+  !=======================================================================
   SUBROUTINE guide_read(timestep)
     IMPLICIT NONE
@@ -1097,295 +1094,294 @@
     include "paramet.h"
 
-    INTEGER, INTENT(IN)   :: timestep
-
-    LOGICAL, SAVE         :: first=.TRUE.
-! Identification fichiers et variables NetCDF:
-    INTEGER, SAVE         :: ncidu,varidu,ncidv,varidv,ncidp,varidp
-    INTEGER, SAVE         :: ncidQ,varidQ,ncidt,varidt,ncidps,varidps
-    INTEGER               :: ncidpl,varidpl,varidap,varidbp,dimid,lendim
-! Variables auxiliaires NetCDF:
-    INTEGER, DIMENSION(4) :: start,count
-    INTEGER               :: status,rcode
-    CHARACTER (len = 80)   :: abort_message
-    CHARACTER (len = 20)   :: modname = 'guide_read'
-    CHARACTER (len = 20)   :: namedim
-
-! -----------------------------------------------------------------
-! Premier appel: initialisation de la lecture des fichiers
-! -----------------------------------------------------------------
+    INTEGER, INTENT(IN) :: timestep
+
+    LOGICAL, SAVE :: first = .TRUE.
+    ! Identification fichiers et variables NetCDF:
+    INTEGER, SAVE :: ncidu, varidu, ncidv, varidv, ncidp, varidp
+    INTEGER, SAVE :: ncidQ, varidQ, ncidt, varidt, ncidps, varidps
+    INTEGER :: ncidpl, varidpl, varidap, varidbp, dimid, lendim
+    ! Variables auxiliaires NetCDF:
+    INTEGER, DIMENSION(4) :: start, count
+    INTEGER :: status, rcode
+    CHARACTER (len = 80) :: abort_message
+    CHARACTER (len = 20) :: modname = 'guide_read'
+    CHARACTER (len = 20) :: namedim
+
+    ! -----------------------------------------------------------------
+    ! Premier appel: initialisation de la lecture des fichiers
+    ! -----------------------------------------------------------------
     IF (first) THEN
-         ncidpl=-99
-         WRITE(*,*) trim(modname)//': opening nudging files '
-! Niveaux de pression si non constants
-         IF (guide_plevs==1) THEN
-             WRITE(*,*) trim(modname)//' Reading nudging on model levels'
-             rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
-             IF (rcode/=nf90_noerr) THEN
-              abort_message='Nudging: error -> no file apbp.nc'
-              CALL abort_gcm(modname,abort_message,1)
-             ENDIF
-             rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
-             IF (rcode/=nf90_noerr) THEN
-              abort_message='Nudging: error -> no AP variable in file apbp.nc'
-              CALL abort_gcm(modname,abort_message,1)
-             ENDIF
-             rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
-             IF (rcode/=nf90_noerr) THEN
-              abort_message='Nudging: error -> no BP variable in file apbp.nc'
-              CALL abort_gcm(modname,abort_message,1)
-             ENDIF
-             WRITE(*,*) trim(modname)//' ncidpl,varidap',ncidpl,varidap
-         endif
-
-! Pression si guidage sur niveaux P variables
-         IF (guide_plevs==2) THEN
-             rcode = nf90_open('P.nc', nf90_nowrite, ncidp)
-             IF (rcode/=nf90_noerr) THEN
-              abort_message='Nudging: error -> no file P.nc'
-              CALL abort_gcm(modname,abort_message,1)
-             ENDIF
-             rcode = nf90_inq_varid(ncidp, 'PRES', varidp)
-             IF (rcode/=nf90_noerr) THEN
-              abort_message='Nudging: error -> no PRES variable in file P.nc'
-              CALL abort_gcm(modname,abort_message,1)
-             ENDIF
-             WRITE(*,*) trim(modname)//' ncidp,varidp',ncidp,varidp
-             IF (ncidpl==-99) ncidpl=ncidp
-         endif
-
-! Vent zonal
-         IF (guide_u) THEN
-             rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
-             IF (rcode/=nf90_noerr) THEN
-              abort_message='Nudging: error -> no file u.nc'
-              CALL abort_gcm(modname,abort_message,1)
-             ENDIF
-             rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
-             IF (rcode/=nf90_noerr) THEN
-              abort_message='Nudging: error -> no UWND variable in file u.nc'
-              CALL abort_gcm(modname,abort_message,1)
-             ENDIF
-             WRITE(*,*) trim(modname)//' ncidu,varidu',ncidu,varidu
-             IF (ncidpl==-99) ncidpl=ncidu
-
-             status=nf90_inq_dimid(ncidu, "LONU", dimid)
-             status=nf90_inquire_dimension(ncidu,dimid,namedim,lendim)
-             IF (lendim /= iip1) THEN
-                abort_message='dimension LONU different from iip1 in u.nc'
-                CALL abort_gcm(modname,abort_message,1)
-             ENDIF
-
-             status=nf90_inq_dimid(ncidu, "LATU", dimid)
-             status=nf90_inquire_dimension(ncidu,dimid,namedim,lendim)
-             IF (lendim /= jjp1) THEN
-                abort_message='dimension LATU different from jjp1 in u.nc'
-                CALL abort_gcm(modname,abort_message,1)
-             ENDIF
-
-         endif
-
-! Vent meridien
-         IF (guide_v) THEN
-             rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
-             IF (rcode/=nf90_noerr) THEN
-              abort_message='Nudging: error -> no file v.nc'
-              CALL abort_gcm(modname,abort_message,1)
-             ENDIF
-             rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
-             IF (rcode/=nf90_noerr) THEN
-              abort_message='Nudging: error -> no VWND variable in file v.nc'
-              CALL abort_gcm(modname,abort_message,1)
-             ENDIF
-             WRITE(*,*) trim(modname)//' ncidv,varidv',ncidv,varidv
-             IF (ncidpl==-99) ncidpl=ncidv
-             
-             status=nf90_inq_dimid(ncidv, "LONV", dimid)
-             status=nf90_inquire_dimension(ncidv,dimid,namedim,lendim)
-             
-                IF (lendim /= iip1) THEN
-                abort_message='dimension LONV different from iip1 in v.nc'
-                CALL abort_gcm(modname,abort_message,1)
-             ENDIF
-
-
-             status=nf90_inq_dimid(ncidv, "LATV", dimid)
-             status=nf90_inquire_dimension(ncidv,dimid,namedim,lendim)
-             IF (lendim /= jjm) THEN
-                abort_message='dimension LATV different from jjm in v.nc'
-                CALL abort_gcm(modname,abort_message,1)
-             ENDIF
-        
-         endif
-
-! Temperature
-         IF (guide_T) THEN
-             rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
-             IF (rcode/=nf90_noerr) THEN
-              abort_message='Nudging: error -> no file T.nc'
-              CALL abort_gcm(modname,abort_message,1)
-             ENDIF
-             rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
-             IF (rcode/=nf90_noerr) THEN
-              abort_message='Nudging: error -> no AIR variable in file T.nc'
-              CALL abort_gcm(modname,abort_message,1)
-             ENDIF
-             WRITE(*,*) trim(modname)//' ncidT,varidT',ncidt,varidt
-             IF (ncidpl==-99) ncidpl=ncidt
-
-             status=nf90_inq_dimid(ncidt, "LONV", dimid)
-             status=nf90_inquire_dimension(ncidt,dimid,namedim,lendim)
-             IF (lendim /= iip1) THEN
-                abort_message='dimension LONV different from iip1 in T.nc'
-                CALL abort_gcm(modname,abort_message,1)
-             ENDIF
-
-             status=nf90_inq_dimid(ncidt, "LATU", dimid)
-             status=nf90_inquire_dimension(ncidt,dimid,namedim,lendim)
-             IF (lendim /= jjp1) THEN
-                abort_message='dimension LATU different from jjp1 in T.nc'
-                CALL abort_gcm(modname,abort_message,1)
-             ENDIF
-
-         endif
-
-! Humidite
-         IF (guide_Q) THEN
-             rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
-             IF (rcode/=nf90_noerr) THEN
-              abort_message='Nudging: error -> no file hur.nc'
-              CALL abort_gcm(modname,abort_message,1)
-             ENDIF
-             rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
-             IF (rcode/=nf90_noerr) THEN
-              abort_message='Nudging: error -> no RH variable in file hur.nc'
-              CALL abort_gcm(modname,abort_message,1)
-             ENDIF
-             WRITE(*,*) trim(modname)//' ncidQ,varidQ',ncidQ,varidQ
-             IF (ncidpl==-99) ncidpl=ncidQ
-
-             status=nf90_inq_dimid(ncidQ, "LONV", dimid)
-             status=nf90_inquire_dimension(ncidQ,dimid,namedim,lendim)
-             IF (lendim /= iip1) THEN
-                abort_message='dimension LONV different from iip1 in hur.nc'
-                CALL abort_gcm(modname,abort_message,1)
-             ENDIF
-
-             status=nf90_inq_dimid(ncidQ, "LATU", dimid)
-             status=nf90_inquire_dimension(ncidQ,dimid,namedim,lendim)
-             IF (lendim /= jjp1) THEN
-                abort_message='dimension LATU different from jjp1 in hur.nc'
-                CALL abort_gcm(modname,abort_message,1)
-             ENDIF
-
-         endif
-
-! Pression de surface
-         IF ((guide_P).OR.(guide_modele)) THEN
-             rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
-             IF (rcode/=nf90_noerr) THEN
-              abort_message='Nudging: error -> no file ps.nc'
-              CALL abort_gcm(modname,abort_message,1)
-             ENDIF
-             rcode = nf90_inq_varid(ncidps, 'SP', varidps)
-             IF (rcode/=nf90_noerr) THEN
-              abort_message='Nudging: error -> no SP variable in file ps.nc'
-              CALL abort_gcm(modname,abort_message,1)
-             ENDIF
-             WRITE(*,*) trim(modname)//' ncidps,varidps',ncidps,varidps
-         endif
-! Coordonnee verticale
-         IF (guide_plevs==0) THEN
-              rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
-              IF (rcode/=0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
-              WRITE(*,*) trim(modname)//' ncidpl,varidpl',ncidpl,varidpl
-         endif
-! Coefs ap, bp pour calcul de la pression aux differents niveaux
-         IF (guide_plevs==1) THEN
-             status=nf90_get_var(ncidpl,varidap,apnc,[1],[nlevnc])
-             status=nf90_get_var(ncidpl,varidbp,bpnc,[1],[nlevnc])
-         ELSEIF (guide_plevs==0) THEN
-             status=nf90_get_var(ncidpl,varidpl,apnc,[1],[nlevnc])
-!FC Pour les corrections la pression est deja en Pascals on commente la ligne ci-dessous
-             IF(convert_Pa) apnc=apnc*100.! conversion en Pascals
-             bpnc(:)=0.
-         endif
-         first=.FALSE.
-     endif ! (first)
-
-! -----------------------------------------------------------------
-!   lecture des champs u, v, T, Q, ps
-! -----------------------------------------------------------------
-
-!  dimensions pour les champs scalaires et le vent zonal
-     start(1)=1
-     start(2)=1
-     start(3)=1
-     start(4)=timestep
-
-     count(1)=iip1
-     count(2)=jjp1
-     count(3)=nlevnc
-     count(4)=1
-
-! Pression 
-     IF (guide_plevs==2) THEN
-         status=nf90_get_var(ncidp,varidp,pnat2,start,count)
-         IF (invert_y) THEN
-!           PRINT*,"Invertion impossible actuellement"
-!           CALL abort_gcm(modname,abort_message,1)
-           CALL invert_lat(iip1,jjp1,nlevnc,pnat2)
-         ENDIF
-     endif
-
-!  Vent zonal
-     IF (guide_u) THEN
-         status=nf90_get_var(ncidu,varidu,unat2,start,count)
-         IF (invert_y) THEN
-           CALL invert_lat(iip1,jjp1,nlevnc,unat2)
-         ENDIF
-     endif
-
-!  Temperature
-     IF (guide_T) THEN
-         status=nf90_get_var(ncidt,varidt,tnat2,start,count)
-         IF (invert_y) THEN
-           CALL invert_lat(iip1,jjp1,nlevnc,tnat2)
-         ENDIF
-     endif
-
-!  Humidite
-     IF (guide_Q) THEN
-         status=nf90_get_var(ncidQ,varidQ,qnat2,start,count)
-         IF (invert_y) THEN
-           CALL invert_lat(iip1,jjp1,nlevnc,qnat2)
-         ENDIF
-         
-     endif
-
-!  Vent meridien
-     IF (guide_v) THEN
-         count(2)=jjm
-         status=nf90_get_var(ncidv,varidv,vnat2,start,count)
-         IF (invert_y) THEN
-           CALL invert_lat(iip1,jjm,nlevnc,vnat2)
-         ENDIF
-     endif
-
-!  Pression de surface
-     IF ((guide_P).OR.(guide_modele))  THEN
-         start(3)=timestep
-         start(4)=0
-         count(2)=jjp1
-         count(3)=1
-         count(4)=0
-         status=nf90_get_var(ncidps,varidps,psnat2,start,count)
-         IF (invert_y) THEN
-           CALL invert_lat(iip1,jjp1,1,psnat2)
-         ENDIF
-     endif
+      ncidpl = -99
+      WRITE(*, *) trim(modname) // ': opening nudging files '
+      ! Niveaux de pression si non constants
+      IF (guide_plevs==1) THEN
+        WRITE(*, *) trim(modname) // ' Reading nudging on model levels'
+        rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no file apbp.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no AP variable in file apbp.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no BP variable in file apbp.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        WRITE(*, *) trim(modname) // ' ncidpl,varidap', ncidpl, varidap
+      endif
+
+      ! Pression si guidage sur niveaux P variables
+      IF (guide_plevs==2) THEN
+        rcode = nf90_open('P.nc', nf90_nowrite, ncidp)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no file P.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        rcode = nf90_inq_varid(ncidp, 'PRES', varidp)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no PRES variable in file P.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        WRITE(*, *) trim(modname) // ' ncidp,varidp', ncidp, varidp
+        IF (ncidpl==-99) ncidpl = ncidp
+      endif
+
+      ! Vent zonal
+      IF (guide_u) THEN
+        rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no file u.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no UWND variable in file u.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        WRITE(*, *) trim(modname) // ' ncidu,varidu', ncidu, varidu
+        IF (ncidpl==-99) ncidpl = ncidu
+
+        status = nf90_inq_dimid(ncidu, "LONU", dimid)
+        status = nf90_inquire_dimension(ncidu, dimid, namedim, lendim)
+        IF (lendim /= iip1) THEN
+          abort_message = 'dimension LONU different from iip1 in u.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+
+        status = nf90_inq_dimid(ncidu, "LATU", dimid)
+        status = nf90_inquire_dimension(ncidu, dimid, namedim, lendim)
+        IF (lendim /= jjp1) THEN
+          abort_message = 'dimension LATU different from jjp1 in u.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+
+      endif
+
+      ! Vent meridien
+      IF (guide_v) THEN
+        rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no file v.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no VWND variable in file v.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        WRITE(*, *) trim(modname) // ' ncidv,varidv', ncidv, varidv
+        IF (ncidpl==-99) ncidpl = ncidv
+
+        status = nf90_inq_dimid(ncidv, "LONV", dimid)
+        status = nf90_inquire_dimension(ncidv, dimid, namedim, lendim)
+
+        IF (lendim /= iip1) THEN
+          abort_message = 'dimension LONV different from iip1 in v.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+
+        status = nf90_inq_dimid(ncidv, "LATV", dimid)
+        status = nf90_inquire_dimension(ncidv, dimid, namedim, lendim)
+        IF (lendim /= jjm) THEN
+          abort_message = 'dimension LATV different from jjm in v.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+
+      endif
+
+      ! Temperature
+      IF (guide_T) THEN
+        rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no file T.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no AIR variable in file T.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        WRITE(*, *) trim(modname) // ' ncidT,varidT', ncidt, varidt
+        IF (ncidpl==-99) ncidpl = ncidt
+
+        status = nf90_inq_dimid(ncidt, "LONV", dimid)
+        status = nf90_inquire_dimension(ncidt, dimid, namedim, lendim)
+        IF (lendim /= iip1) THEN
+          abort_message = 'dimension LONV different from iip1 in T.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+
+        status = nf90_inq_dimid(ncidt, "LATU", dimid)
+        status = nf90_inquire_dimension(ncidt, dimid, namedim, lendim)
+        IF (lendim /= jjp1) THEN
+          abort_message = 'dimension LATU different from jjp1 in T.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+
+      endif
+
+      ! Humidite
+      IF (guide_Q) THEN
+        rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no file hur.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no RH variable in file hur.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        WRITE(*, *) trim(modname) // ' ncidQ,varidQ', ncidQ, varidQ
+        IF (ncidpl==-99) ncidpl = ncidQ
+
+        status = nf90_inq_dimid(ncidQ, "LONV", dimid)
+        status = nf90_inquire_dimension(ncidQ, dimid, namedim, lendim)
+        IF (lendim /= iip1) THEN
+          abort_message = 'dimension LONV different from iip1 in hur.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+
+        status = nf90_inq_dimid(ncidQ, "LATU", dimid)
+        status = nf90_inquire_dimension(ncidQ, dimid, namedim, lendim)
+        IF (lendim /= jjp1) THEN
+          abort_message = 'dimension LATU different from jjp1 in hur.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+
+      endif
+
+      ! Pression de surface
+      IF ((guide_P).OR.(guide_modele)) THEN
+        rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no file ps.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        rcode = nf90_inq_varid(ncidps, 'SP', varidps)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no SP variable in file ps.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        WRITE(*, *) trim(modname) // ' ncidps,varidps', ncidps, varidps
+      endif
+      ! Coordonnee verticale
+      IF (guide_plevs==0) THEN
+        rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
+        IF (rcode/=0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
+        WRITE(*, *) trim(modname) // ' ncidpl,varidpl', ncidpl, varidpl
+      endif
+      ! Coefs ap, bp pour calcul de la pression aux differents niveaux
+      IF (guide_plevs==1) THEN
+        status = nf90_get_var(ncidpl, varidap, apnc, [1], [nlevnc])
+        status = nf90_get_var(ncidpl, varidbp, bpnc, [1], [nlevnc])
+      ELSEIF (guide_plevs==0) THEN
+        status = nf90_get_var(ncidpl, varidpl, apnc, [1], [nlevnc])
+        !FC Pour les corrections la pression est deja en Pascals on commente la ligne ci-dessous
+        IF(convert_Pa) apnc = apnc * 100.! conversion en Pascals
+        bpnc(:) = 0.
+      endif
+      first = .FALSE.
+    endif ! (first)
+
+    ! -----------------------------------------------------------------
+    !   lecture des champs u, v, T, Q, ps
+    ! -----------------------------------------------------------------
+
+    !  dimensions pour les champs scalaires et le vent zonal
+    start(1) = 1
+    start(2) = 1
+    start(3) = 1
+    start(4) = timestep
+
+    count(1) = iip1
+    count(2) = jjp1
+    count(3) = nlevnc
+    count(4) = 1
+
+    ! Pression
+    IF (guide_plevs==2) THEN
+      status = nf90_get_var(ncidp, varidp, pnat2, start, count)
+      IF (invert_y) THEN
+        !           PRINT*,"Invertion impossible actuellement"
+        !           CALL abort_gcm(modname,abort_message,1)
+        CALL invert_lat(iip1, jjp1, nlevnc, pnat2)
+      ENDIF
+    endif
+
+    !  Vent zonal
+    IF (guide_u) THEN
+      status = nf90_get_var(ncidu, varidu, unat2, start, count)
+      IF (invert_y) THEN
+        CALL invert_lat(iip1, jjp1, nlevnc, unat2)
+      ENDIF
+    endif
+
+    !  Temperature
+    IF (guide_T) THEN
+      status = nf90_get_var(ncidt, varidt, tnat2, start, count)
+      IF (invert_y) THEN
+        CALL invert_lat(iip1, jjp1, nlevnc, tnat2)
+      ENDIF
+    endif
+
+    !  Humidite
+    IF (guide_Q) THEN
+      status = nf90_get_var(ncidQ, varidQ, qnat2, start, count)
+      IF (invert_y) THEN
+        CALL invert_lat(iip1, jjp1, nlevnc, qnat2)
+      ENDIF
+
+    endif
+
+    !  Vent meridien
+    IF (guide_v) THEN
+      count(2) = jjm
+      status = nf90_get_var(ncidv, varidv, vnat2, start, count)
+      IF (invert_y) THEN
+        CALL invert_lat(iip1, jjm, nlevnc, vnat2)
+      ENDIF
+    endif
+
+    !  Pression de surface
+    IF ((guide_P).OR.(guide_modele))  THEN
+      start(3) = timestep
+      start(4) = 0
+      count(2) = jjp1
+      count(3) = 1
+      count(4) = 0
+      status = nf90_get_var(ncidps, varidps, psnat2, start, count)
+      IF (invert_y) THEN
+        CALL invert_lat(iip1, jjp1, 1, psnat2)
+      ENDIF
+    endif
 
   END SUBROUTINE guide_read
 
-!=======================================================================
+  !=======================================================================
   SUBROUTINE guide_read2D(timestep)
     IMPLICIT NONE
@@ -1394,257 +1390,257 @@
     include "paramet.h"
 
-    INTEGER, INTENT(IN)   :: timestep
-
-    LOGICAL, SAVE         :: first=.TRUE.
-! Identification fichiers et variables NetCDF:
-    INTEGER, SAVE         :: ncidu,varidu,ncidv,varidv,ncidp,varidp
-    INTEGER, SAVE         :: ncidQ,varidQ,ncidt,varidt,ncidps,varidps
-    INTEGER               :: ncidpl,varidpl,varidap,varidbp
-! Variables auxiliaires NetCDF:
-    INTEGER, DIMENSION(4) :: start,count
-    INTEGER               :: status,rcode
-! Variables for 3D extension:
-    REAL, DIMENSION (jjp1,llm) :: zu
-    REAL, DIMENSION (jjm,llm)  :: zv
-    INTEGER               :: i
-
-    CHARACTER (len = 80)   :: abort_message
-    CHARACTER (len = 20)   :: modname = 'guide_read2D'
-! -----------------------------------------------------------------
-! Premier appel: initialisation de la lecture des fichiers
-! -----------------------------------------------------------------
+    INTEGER, INTENT(IN) :: timestep
+
+    LOGICAL, SAVE :: first = .TRUE.
+    ! Identification fichiers et variables NetCDF:
+    INTEGER, SAVE :: ncidu, varidu, ncidv, varidv, ncidp, varidp
+    INTEGER, SAVE :: ncidQ, varidQ, ncidt, varidt, ncidps, varidps
+    INTEGER :: ncidpl, varidpl, varidap, varidbp
+    ! Variables auxiliaires NetCDF:
+    INTEGER, DIMENSION(4) :: start, count
+    INTEGER :: status, rcode
+    ! Variables for 3D extension:
+    REAL, DIMENSION (jjp1, llm) :: zu
+    REAL, DIMENSION (jjm, llm) :: zv
+    INTEGER :: i
+
+    CHARACTER (len = 80) :: abort_message
+    CHARACTER (len = 20) :: modname = 'guide_read2D'
+    ! -----------------------------------------------------------------
+    ! Premier appel: initialisation de la lecture des fichiers
+    ! -----------------------------------------------------------------
     IF (first) THEN
-         ncidpl=-99
-         WRITE(*,*)trim(modname)//' : opening nudging files '
-! Ap et Bp si niveaux de pression hybrides
-         IF (guide_plevs==1) THEN
-           WRITE(*,*)trim(modname)//' Reading nudging on model levels'
-           rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
-           IF (rcode/=nf90_noerr) THEN
-             abort_message='Nudging: error -> no file apbp.nc'
-           CALL abort_gcm(modname,abort_message,1)
-           ENDIF
-           rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
-           IF (rcode/=nf90_noerr) THEN
-             abort_message='Nudging: error -> no AP variable in file apbp.nc'
-           CALL abort_gcm(modname,abort_message,1)
-           ENDIF
-           rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
-           IF (rcode/=nf90_noerr) THEN
-             abort_message='Nudging: error -> no BP variable in file apbp.nc'
-             CALL abort_gcm(modname,abort_message,1)
-           ENDIF
-           WRITE(*,*)trim(modname)//'ncidpl,varidap',ncidpl,varidap
-         endif
-! Pression
-         IF (guide_plevs==2) THEN
-           rcode = nf90_open('P.nc', nf90_nowrite, ncidp)
-           IF (rcode/=nf90_noerr) THEN
-             abort_message='Nudging: error -> no file P.nc'
-             CALL abort_gcm(modname,abort_message,1)
-           ENDIF
-           rcode = nf90_inq_varid(ncidp, 'PRES', varidp)
-           IF (rcode/=nf90_noerr) THEN
-             abort_message='Nudging: error -> no PRES variable in file P.nc'
-             CALL abort_gcm(modname,abort_message,1)
-           ENDIF
-           WRITE(*,*)trim(modname)//' ncidp,varidp',ncidp,varidp
-           IF (ncidpl==-99) ncidpl=ncidp
-         endif
-! Vent zonal
-         IF (guide_u) THEN
-           rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
-           IF (rcode/=nf90_noerr) THEN
-             abort_message='Nudging: error -> no file u.nc'
-             CALL abort_gcm(modname,abort_message,1)
-           ENDIF
-           rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
-           IF (rcode/=nf90_noerr) THEN
-             abort_message='Nudging: error -> no UWND variable in file u.nc'
-             CALL abort_gcm(modname,abort_message,1)
-           ENDIF
-           WRITE(*,*)trim(modname)//' ncidu,varidu',ncidu,varidu
-           IF (ncidpl==-99) ncidpl=ncidu
-         endif
-! Vent meridien
-         IF (guide_v) THEN
-           rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
-           IF (rcode/=nf90_noerr) THEN
-             abort_message='Nudging: error -> no file v.nc'
-             CALL abort_gcm(modname,abort_message,1)
-           ENDIF
-           rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
-           IF (rcode/=nf90_noerr) THEN
-             abort_message='Nudging: error -> no VWND variable in file v.nc'
-             CALL abort_gcm(modname,abort_message,1)
-           ENDIF
-           WRITE(*,*)trim(modname)//' ncidv,varidv',ncidv,varidv
-           IF (ncidpl==-99) ncidpl=ncidv
-         endif
-! Temperature
-         IF (guide_T) THEN
-           rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
-           IF (rcode/=nf90_noerr) THEN
-             abort_message='Nudging: error -> no file T.nc'
-             CALL abort_gcm(modname,abort_message,1)
-           ENDIF
-           rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
-           IF (rcode/=nf90_noerr) THEN
-             abort_message='Nudging: error -> no AIR variable in file T.nc'
-             CALL abort_gcm(modname,abort_message,1)
-           ENDIF
-           WRITE(*,*)trim(modname)//' ncidT,varidT',ncidt,varidt
-           IF (ncidpl==-99) ncidpl=ncidt
-         endif
-! Humidite
-         IF (guide_Q) THEN
-           rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
-           IF (rcode/=nf90_noerr) THEN
-             abort_message='Nudging: error -> no file hur.nc'
-             CALL abort_gcm(modname,abort_message,1)
-           ENDIF
-           rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
-           IF (rcode/=nf90_noerr) THEN
-             abort_message='Nudging: error -> no RH,variable in file hur.nc'
-             CALL abort_gcm(modname,abort_message,1)
-           ENDIF
-           WRITE(*,*)trim(modname)//' ncidQ,varidQ',ncidQ,varidQ
-           IF (ncidpl==-99) ncidpl=ncidQ
-         endif
-! Pression de surface
-         IF ((guide_P).OR.(guide_modele)) THEN
-           rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
-           IF (rcode/=nf90_noerr) THEN
-             abort_message='Nudging: error -> no file ps.nc'
-             CALL abort_gcm(modname,abort_message,1)
-           ENDIF
-           rcode = nf90_inq_varid(ncidps, 'SP', varidps)
-           IF (rcode/=nf90_noerr) THEN
-             abort_message='Nudging: error -> no SP variable in file ps.nc'
-             CALL abort_gcm(modname,abort_message,1)
-           ENDIF
-           WRITE(*,*)trim(modname)//' ncidps,varidps',ncidps,varidps
-         endif
-! Coordonnee verticale
-         IF (guide_plevs==0) THEN
-           rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
-           IF (rcode/=0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
-           WRITE(*,*)trim(modname)//' ncidpl,varidpl',ncidpl,varidpl
-         endif
-! Coefs ap, bp pour calcul de la pression aux differents niveaux
-         IF (guide_plevs==1) THEN
-             status=nf90_get_var(ncidpl,varidap,apnc,[1],[nlevnc])
-             status=nf90_get_var(ncidpl,varidbp,bpnc,[1],[nlevnc])
-         elseif (guide_plevs==0) THEN
-             status=nf90_get_var(ncidpl,varidpl,apnc,[1],[nlevnc])
-             apnc=apnc*100.! conversion en Pascals
-             bpnc(:)=0.
-         endif
-         first=.FALSE.
-     endif ! (first)
-
-! -----------------------------------------------------------------
-!   lecture des champs u, v, T, Q, ps
-! -----------------------------------------------------------------
-
-!  dimensions pour les champs scalaires et le vent zonal
-     start(1)=1
-     start(2)=1
-     start(3)=1
-     start(4)=timestep
-
-     count(1)=1
-     count(2)=jjp1
-     count(3)=nlevnc
-     count(4)=1
-
-!  Pression
-     IF (guide_plevs==2) THEN
-         status=nf90_get_var(ncidp,varidp,zu,start,count)
-         DO i=1,iip1
-             pnat2(i,:,:)=zu(:,:)
-         ENDDO
-
-         IF (invert_y) THEN
-!           PRINT*,"Invertion impossible actuellement"
-!           CALL abort_gcm(modname,abort_message,1)
-           CALL invert_lat(iip1,jjp1,nlevnc,pnat2)
-         ENDIF
-     endif
-!  Vent zonal
-     IF (guide_u) THEN
-         status=nf90_get_var(ncidu,varidu,zu,start,count)
-         DO i=1,iip1
-             unat2(i,:,:)=zu(:,:)
-         ENDDO
-
-         IF (invert_y) THEN
-           CALL invert_lat(iip1,jjp1,nlevnc,unat2)
-         ENDIF
-
-     endif
-
-!  Temperature
-     IF (guide_T) THEN
-         status=nf90_get_var(ncidt,varidt,zu,start,count)
-         DO i=1,iip1
-             tnat2(i,:,:)=zu(:,:)
-         ENDDO
-
-         IF (invert_y) THEN
-           CALL invert_lat(iip1,jjp1,nlevnc,tnat2)
-         ENDIF
-
-     endif
-
-!  Humidite
-     IF (guide_Q) THEN
-         status=nf90_get_var(ncidQ,varidQ,zu,start,count)
-         DO i=1,iip1
-             qnat2(i,:,:)=zu(:,:)
-         ENDDO
-
-         IF (invert_y) THEN
-           CALL invert_lat(iip1,jjp1,nlevnc,qnat2)
-         ENDIF
-
-     endif
-
-!  Vent meridien
-     IF (guide_v) THEN
-         count(2)=jjm
-         status=nf90_get_var(ncidv,varidv,zv,start,count)
-         DO i=1,iip1
-             vnat2(i,:,:)=zv(:,:)
-         ENDDO
-
-         IF (invert_y) THEN
-           CALL invert_lat(iip1,jjm,nlevnc,vnat2)
-         ENDIF
-
-     endif
-
-!  Pression de surface
-     IF ((guide_P).OR.(guide_plevs==1))  THEN
-         start(3)=timestep
-         start(4)=0
-         count(2)=jjp1
-         count(3)=1
-         count(4)=0
-         status=nf90_get_var(ncidps,varidps,zu(:,1),start,count)
-         DO i=1,iip1
-             psnat2(i,:)=zu(:,1)
-         ENDDO
-
-         IF (invert_y) THEN
-           CALL invert_lat(iip1,jjp1,1,psnat2)
-         ENDIF
-
-     endif
+      ncidpl = -99
+      WRITE(*, *)trim(modname) // ' : opening nudging files '
+      ! Ap et Bp si niveaux de pression hybrides
+      IF (guide_plevs==1) THEN
+        WRITE(*, *)trim(modname) // ' Reading nudging on model levels'
+        rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no file apbp.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no AP variable in file apbp.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no BP variable in file apbp.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        WRITE(*, *)trim(modname) // 'ncidpl,varidap', ncidpl, varidap
+      endif
+      ! Pression
+      IF (guide_plevs==2) THEN
+        rcode = nf90_open('P.nc', nf90_nowrite, ncidp)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no file P.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        rcode = nf90_inq_varid(ncidp, 'PRES', varidp)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no PRES variable in file P.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        WRITE(*, *)trim(modname) // ' ncidp,varidp', ncidp, varidp
+        IF (ncidpl==-99) ncidpl = ncidp
+      endif
+      ! Vent zonal
+      IF (guide_u) THEN
+        rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no file u.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no UWND variable in file u.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        WRITE(*, *)trim(modname) // ' ncidu,varidu', ncidu, varidu
+        IF (ncidpl==-99) ncidpl = ncidu
+      endif
+      ! Vent meridien
+      IF (guide_v) THEN
+        rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no file v.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no VWND variable in file v.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        WRITE(*, *)trim(modname) // ' ncidv,varidv', ncidv, varidv
+        IF (ncidpl==-99) ncidpl = ncidv
+      endif
+      ! Temperature
+      IF (guide_T) THEN
+        rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no file T.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no AIR variable in file T.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        WRITE(*, *)trim(modname) // ' ncidT,varidT', ncidt, varidt
+        IF (ncidpl==-99) ncidpl = ncidt
+      endif
+      ! Humidite
+      IF (guide_Q) THEN
+        rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no file hur.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no RH,variable in file hur.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        WRITE(*, *)trim(modname) // ' ncidQ,varidQ', ncidQ, varidQ
+        IF (ncidpl==-99) ncidpl = ncidQ
+      endif
+      ! Pression de surface
+      IF ((guide_P).OR.(guide_modele)) THEN
+        rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no file ps.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        rcode = nf90_inq_varid(ncidps, 'SP', varidps)
+        IF (rcode/=nf90_noerr) THEN
+          abort_message = 'Nudging: error -> no SP variable in file ps.nc'
+          CALL abort_gcm(modname, abort_message, 1)
+        ENDIF
+        WRITE(*, *)trim(modname) // ' ncidps,varidps', ncidps, varidps
+      endif
+      ! Coordonnee verticale
+      IF (guide_plevs==0) THEN
+        rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
+        IF (rcode/=0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
+        WRITE(*, *)trim(modname) // ' ncidpl,varidpl', ncidpl, varidpl
+      endif
+      ! Coefs ap, bp pour calcul de la pression aux differents niveaux
+      IF (guide_plevs==1) THEN
+        status = nf90_get_var(ncidpl, varidap, apnc, [1], [nlevnc])
+        status = nf90_get_var(ncidpl, varidbp, bpnc, [1], [nlevnc])
+      elseif (guide_plevs==0) THEN
+        status = nf90_get_var(ncidpl, varidpl, apnc, [1], [nlevnc])
+        apnc = apnc * 100.! conversion en Pascals
+        bpnc(:) = 0.
+      endif
+      first = .FALSE.
+    endif ! (first)
+
+    ! -----------------------------------------------------------------
+    !   lecture des champs u, v, T, Q, ps
+    ! -----------------------------------------------------------------
+
+    !  dimensions pour les champs scalaires et le vent zonal
+    start(1) = 1
+    start(2) = 1
+    start(3) = 1
+    start(4) = timestep
+
+    count(1) = 1
+    count(2) = jjp1
+    count(3) = nlevnc
+    count(4) = 1
+
+    !  Pression
+    IF (guide_plevs==2) THEN
+      status = nf90_get_var(ncidp, varidp, zu, start, count)
+      DO i = 1, iip1
+        pnat2(i, :, :) = zu(:, :)
+      ENDDO
+
+      IF (invert_y) THEN
+        !           PRINT*,"Invertion impossible actuellement"
+        !           CALL abort_gcm(modname,abort_message,1)
+        CALL invert_lat(iip1, jjp1, nlevnc, pnat2)
+      ENDIF
+    endif
+    !  Vent zonal
+    IF (guide_u) THEN
+      status = nf90_get_var(ncidu, varidu, zu, start, count)
+      DO i = 1, iip1
+        unat2(i, :, :) = zu(:, :)
+      ENDDO
+
+      IF (invert_y) THEN
+        CALL invert_lat(iip1, jjp1, nlevnc, unat2)
+      ENDIF
+
+    endif
+
+    !  Temperature
+    IF (guide_T) THEN
+      status = nf90_get_var(ncidt, varidt, zu, start, count)
+      DO i = 1, iip1
+        tnat2(i, :, :) = zu(:, :)
+      ENDDO
+
+      IF (invert_y) THEN
+        CALL invert_lat(iip1, jjp1, nlevnc, tnat2)
+      ENDIF
+
+    endif
+
+    !  Humidite
+    IF (guide_Q) THEN
+      status = nf90_get_var(ncidQ, varidQ, zu, start, count)
+      DO i = 1, iip1
+        qnat2(i, :, :) = zu(:, :)
+      ENDDO
+
+      IF (invert_y) THEN
+        CALL invert_lat(iip1, jjp1, nlevnc, qnat2)
+      ENDIF
+
+    endif
+
+    !  Vent meridien
+    IF (guide_v) THEN
+      count(2) = jjm
+      status = nf90_get_var(ncidv, varidv, zv, start, count)
+      DO i = 1, iip1
+        vnat2(i, :, :) = zv(:, :)
+      ENDDO
+
+      IF (invert_y) THEN
+        CALL invert_lat(iip1, jjm, nlevnc, vnat2)
+      ENDIF
+
+    endif
+
+    !  Pression de surface
+    IF ((guide_P).OR.(guide_plevs==1))  THEN
+      start(3) = timestep
+      start(4) = 0
+      count(2) = jjp1
+      count(3) = 1
+      count(4) = 0
+      status = nf90_get_var(ncidps, varidps, zu(:, 1), start, count)
+      DO i = 1, iip1
+        psnat2(i, :) = zu(:, 1)
+      ENDDO
+
+      IF (invert_y) THEN
+        CALL invert_lat(iip1, jjp1, 1, psnat2)
+      ENDIF
+
+    endif
 
   END SUBROUTINE guide_read2D
-  
-!=======================================================================
-  SUBROUTINE guide_out(varname,hsize,vsize,field)
+
+  !=======================================================================
+  SUBROUTINE guide_out(varname, hsize, vsize, field)
 
     USE comconst_mod, ONLY: pi
@@ -1657,162 +1653,164 @@
     INCLUDE "paramet.h"
     INCLUDE "comgeom2.h"
-    
+
     ! Variables entree
-    CHARACTER*(*), INTENT(IN)                          :: varname
-    INTEGER,   INTENT (IN)                         :: hsize,vsize
-    REAL, DIMENSION (iip1,hsize,vsize), INTENT(IN) :: field
+    CHARACTER*(*), INTENT(IN) :: varname
+    INTEGER, INTENT (IN) :: hsize, vsize
+    REAL, DIMENSION (iip1, hsize, vsize), INTENT(IN) :: field
 
     ! Variables locales
-    INTEGER, SAVE :: timestep=0
+    INTEGER, SAVE :: timestep = 0
     ! Identites fichier netcdf
-    INTEGER       :: nid, id_lonu, id_lonv, id_latu, id_latv, id_tim, id_lev
-    INTEGER       :: vid_lonu,vid_lonv,vid_latu,vid_latv,vid_cu,vid_cv,vid_lev
-    INTEGER       :: vid_au,vid_av, varid_alpha_t, varid_alpha_q
+    INTEGER :: nid, id_lonu, id_lonv, id_latu, id_latv, id_tim, id_lev
+    INTEGER :: vid_lonu, vid_lonv, vid_latu, vid_latv, vid_cu, vid_cv, vid_lev
+    INTEGER :: vid_au, vid_av, varid_alpha_t, varid_alpha_q
     INTEGER, DIMENSION (3) :: dim3
-    INTEGER, DIMENSION (4) :: dim4,count,start
-    INTEGER                :: ierr, varid,l
-    REAL, DIMENSION (iip1,hsize,vsize) :: field2
-    CHARACTER(LEN=20),PARAMETER :: modname="guide_out"
-
-    WRITE(*,*)trim(modname)//': output timestep',timestep,'var ',varname
+    INTEGER, DIMENSION (4) :: dim4, count, start
+    INTEGER :: ierr, varid, l
+    REAL, DIMENSION (iip1, hsize, vsize) :: field2
+    CHARACTER(LEN = 20), PARAMETER :: modname = "guide_out"
+
+    WRITE(*, *)trim(modname) // ': output timestep', timestep, 'var ', varname
     IF (timestep==0) THEN
-! ----------------------------------------------
-! initialisation fichier de sortie
-! ----------------------------------------------
-! Ouverture du fichier
-        ierr=nf90_create("guide_ins.nc",IOR(nf90_clobber,nf90_64bit_offset),nid)
-! Definition des dimensions
-        ierr=nf90_def_dim(nid,"LONU",iip1,id_lonu)
-        ierr=nf90_def_dim(nid,"LONV",iip1,id_lonv)
-        ierr=nf90_def_dim(nid,"LATU",jjp1,id_latu)
-        ierr=nf90_def_dim(nid,"LATV",jjm,id_latv)
-        ierr=nf90_def_dim(nid,"LEVEL",llm,id_lev)
-        ierr=nf90_def_dim(nid,"TIME",nf90_unlimited,id_tim)
-
-! Creation des variables dimensions
-        ierr=nf90_def_var(nid,"LONU",nf90_float,id_lonu,vid_lonu)
-        ierr=nf90_def_var(nid,"LONV",nf90_float,id_lonv,vid_lonv)
-        ierr=nf90_def_var(nid,"LATU",nf90_float,id_latu,vid_latu)
-        ierr=nf90_def_var(nid,"LATV",nf90_float,id_latv,vid_latv)
-        ierr=nf90_def_var(nid,"LEVEL",nf90_float,id_lev,vid_lev)
-        ierr=nf90_def_var(nid,"cu",nf90_float,(/id_lonu,id_latu/),vid_cu)
-        ierr=nf90_def_var(nid,"cv",nf90_float,(/id_lonv,id_latv/),vid_cv)
-        ierr=nf90_def_var(nid,"au",nf90_float,(/id_lonu,id_latu/),vid_au)
-        ierr=nf90_def_var(nid,"av",nf90_float,(/id_lonv,id_latv/),vid_av)
-        CALL nf95_def_var(nid, "alpha_T", nf90_float, (/id_lonv, id_latu/), &
-             varid_alpha_t)
-        CALL nf95_def_var(nid, "alpha_q", nf90_float, (/id_lonv, id_latu/), &
-             varid_alpha_q)
-        
-        ierr=nf90_enddef(nid)
-
-! Enregistrement des variables dimensions
-        ierr = nf90_put_var(nid,vid_lonu,rlonu*180./pi)
-        ierr = nf90_put_var(nid,vid_lonv,rlonv*180./pi)
-        ierr = nf90_put_var(nid,vid_latu,rlatu*180./pi)
-        ierr = nf90_put_var(nid,vid_latv,rlatv*180./pi)
-        ierr = nf90_put_var(nid,vid_lev,presnivs)
-        ierr = nf90_put_var(nid,vid_cu,cu)
-        ierr = nf90_put_var(nid,vid_cv,cv)
-        ierr = nf90_put_var(nid,vid_au,alpha_u)
-        ierr = nf90_put_var(nid,vid_av,alpha_v)
-        CALL nf95_put_var(nid, varid_alpha_t, alpha_t)
-        CALL nf95_put_var(nid, varid_alpha_q, alpha_q)
-! --------------------------------------------------------------------
-! Création des variables sauvegardées
-! --------------------------------------------------------------------
-        ierr = nf90_redef(nid)
-! Pressure (GCM)
-        dim4=(/id_lonv,id_latu,id_lev,id_tim/)
-        ierr = nf90_def_var(nid,"SP",nf90_float,dim4,varid)
-! Surface pressure (guidage)
-        IF (guide_P) THEN
-            dim3=(/id_lonv,id_latu,id_tim/)
-            ierr = nf90_def_var(nid,"ps",nf90_float,dim3,varid)
-        ENDIF
-! Zonal wind
-        IF (guide_u) THEN
-            dim4=(/id_lonu,id_latu,id_lev,id_tim/)
-            ierr = nf90_def_var(nid,"u",nf90_float,dim4,varid)
-            ierr = nf90_def_var(nid,"ua",nf90_float,dim4,varid)
-            ierr = nf90_def_var(nid,"ucov",nf90_float,dim4,varid)
-        ENDIF
-! Merid. wind
-        IF (guide_v) THEN
-            dim4=(/id_lonv,id_latv,id_lev,id_tim/)
-            ierr = nf90_def_var(nid,"v",nf90_float,dim4,varid)
-            ierr = nf90_def_var(nid,"va",nf90_float,dim4,varid)
-            ierr = nf90_def_var(nid,"vcov",nf90_float,dim4,varid)
-        ENDIF
-! Pot. Temperature
-        IF (guide_T) THEN
-            dim4=(/id_lonv,id_latu,id_lev,id_tim/)
-            ierr = nf90_def_var(nid,"teta",nf90_float,dim4,varid)
-        ENDIF
-! Specific Humidity
-        IF (guide_Q) THEN
-            dim4=(/id_lonv,id_latu,id_lev,id_tim/)
-            ierr = nf90_def_var(nid,"q",nf90_float,dim4,varid)
-        ENDIF
-        
-        ierr = nf90_enddef(nid)
-        ierr = nf90_close(nid)
+      ! ----------------------------------------------
+      ! initialisation fichier de sortie
+      ! ----------------------------------------------
+      ! Ouverture du fichier
+      ierr = nf90_create("guide_ins.nc", IOR(nf90_clobber, nf90_64bit_offset), nid)
+      ! Definition des dimensions
+      ierr = nf90_def_dim(nid, "LONU", iip1, id_lonu)
+      ierr = nf90_def_dim(nid, "LONV", iip1, id_lonv)
+      ierr = nf90_def_dim(nid, "LATU", jjp1, id_latu)
+      ierr = nf90_def_dim(nid, "LATV", jjm, id_latv)
+      ierr = nf90_def_dim(nid, "LEVEL", llm, id_lev)
+      ierr = nf90_def_dim(nid, "TIME", nf90_unlimited, id_tim)
+
+      ! Creation des variables dimensions
+      ierr = nf90_def_var(nid, "LONU", nf90_float, id_lonu, vid_lonu)
+      ierr = nf90_def_var(nid, "LONV", nf90_float, id_lonv, vid_lonv)
+      ierr = nf90_def_var(nid, "LATU", nf90_float, id_latu, vid_latu)
+      ierr = nf90_def_var(nid, "LATV", nf90_float, id_latv, vid_latv)
+      ierr = nf90_def_var(nid, "LEVEL", nf90_float, id_lev, vid_lev)
+      ierr = nf90_def_var(nid, "cu", nf90_float, (/id_lonu, id_latu/), vid_cu)
+      ierr = nf90_def_var(nid, "cv", nf90_float, (/id_lonv, id_latv/), vid_cv)
+      ierr = nf90_def_var(nid, "au", nf90_float, (/id_lonu, id_latu/), vid_au)
+      ierr = nf90_def_var(nid, "av", nf90_float, (/id_lonv, id_latv/), vid_av)
+      CALL nf95_def_var(nid, "alpha_T", nf90_float, (/id_lonv, id_latu/), &
+              varid_alpha_t)
+      CALL nf95_def_var(nid, "alpha_q", nf90_float, (/id_lonv, id_latu/), &
+              varid_alpha_q)
+
+      ierr = nf90_enddef(nid)
+
+      ! Enregistrement des variables dimensions
+      ierr = nf90_put_var(nid, vid_lonu, rlonu * 180. / pi)
+      ierr = nf90_put_var(nid, vid_lonv, rlonv * 180. / pi)
+      ierr = nf90_put_var(nid, vid_latu, rlatu * 180. / pi)
+      ierr = nf90_put_var(nid, vid_latv, rlatv * 180. / pi)
+      ierr = nf90_put_var(nid, vid_lev, presnivs)
+      ierr = nf90_put_var(nid, vid_cu, cu)
+      ierr = nf90_put_var(nid, vid_cv, cv)
+      ierr = nf90_put_var(nid, vid_au, alpha_u)
+      ierr = nf90_put_var(nid, vid_av, alpha_v)
+      CALL nf95_put_var(nid, varid_alpha_t, alpha_t)
+      CALL nf95_put_var(nid, varid_alpha_q, alpha_q)
+      ! --------------------------------------------------------------------
+      ! Création des variables sauvegardées
+      ! --------------------------------------------------------------------
+      ierr = nf90_redef(nid)
+      ! Pressure (GCM)
+      dim4 = (/id_lonv, id_latu, id_lev, id_tim/)
+      ierr = nf90_def_var(nid, "SP", nf90_float, dim4, varid)
+      ! Surface pressure (guidage)
+      IF (guide_P) THEN
+        dim3 = (/id_lonv, id_latu, id_tim/)
+        ierr = nf90_def_var(nid, "ps", nf90_float, dim3, varid)
+      ENDIF
+      ! Zonal wind
+      IF (guide_u) THEN
+        dim4 = (/id_lonu, id_latu, id_lev, id_tim/)
+        ierr = nf90_def_var(nid, "u", nf90_float, dim4, varid)
+        ierr = nf90_def_var(nid, "ua", nf90_float, dim4, varid)
+        ierr = nf90_def_var(nid, "ucov", nf90_float, dim4, varid)
+      ENDIF
+      ! Merid. wind
+      IF (guide_v) THEN
+        dim4 = (/id_lonv, id_latv, id_lev, id_tim/)
+        ierr = nf90_def_var(nid, "v", nf90_float, dim4, varid)
+        ierr = nf90_def_var(nid, "va", nf90_float, dim4, varid)
+        ierr = nf90_def_var(nid, "vcov", nf90_float, dim4, varid)
+      ENDIF
+      ! Pot. Temperature
+      IF (guide_T) THEN
+        dim4 = (/id_lonv, id_latu, id_lev, id_tim/)
+        ierr = nf90_def_var(nid, "teta", nf90_float, dim4, varid)
+      ENDIF
+      ! Specific Humidity
+      IF (guide_Q) THEN
+        dim4 = (/id_lonv, id_latu, id_lev, id_tim/)
+        ierr = nf90_def_var(nid, "q", nf90_float, dim4, varid)
+      ENDIF
+
+      ierr = nf90_enddef(nid)
+      ierr = nf90_close(nid)
     ENDIF ! timestep=0
 
-! --------------------------------------------------------------------
-! Enregistrement du champ
-! --------------------------------------------------------------------
-    ierr=nf90_open("guide_ins.nc",nf90_write,nid)
-
-    IF (varname=="SP") timestep=timestep+1
-
-    ierr = nf90_inq_varid(nid,varname,varid)
+    ! --------------------------------------------------------------------
+    ! Enregistrement du champ
+    ! --------------------------------------------------------------------
+    ierr = nf90_open("guide_ins.nc", nf90_write, nid)
+
+    IF (varname=="SP") timestep = timestep + 1
+
+    ierr = nf90_inq_varid(nid, varname, varid)
     SELECT CASE (varname)
-    CASE ("SP","ps")
-        start=(/1,1,1,timestep/)
-        count=(/iip1,jjp1,llm,1/)
-    CASE ("v","va","vcov")
-        start=(/1,1,1,timestep/)
-        count=(/iip1,jjm,llm,1/)
+    CASE ("SP", "ps")
+      start = (/1, 1, 1, timestep/)
+      count = (/iip1, jjp1, llm, 1/)
+    CASE ("v", "va", "vcov")
+      start = (/1, 1, 1, timestep/)
+      count = (/iip1, jjm, llm, 1/)
     CASE DEFAULT
-        start=(/1,1,1,timestep/)
-        count=(/iip1,jjp1,llm,1/)
+      start = (/1, 1, 1, timestep/)
+      count = (/iip1, jjp1, llm, 1/)
     END SELECT
 
     SELECT CASE (varname)
-    CASE("u","ua")
-        DO l=1,llm ; field2(:,2:jjm,l)=field(:,2:jjm,l)/cu(:,2:jjm) ; ENDDO
-        field2(:,1,:)=0. ; field2(:,jjp1,:)=0.
-    CASE("v","va")
-        DO l=1,llm ; field2(:,:,l)=field(:,:,l)/cv(:,:) ; ENDDO
+    CASE("u", "ua")
+      DO l = 1, llm ; field2(:, 2:jjm, l) = field(:, 2:jjm, l) / cu(:, 2:jjm) ;
+      ENDDO
+      field2(:, 1, :) = 0. ; field2(:, jjp1, :) = 0.
+    CASE("v", "va")
+      DO l = 1, llm ; field2(:, :, l) = field(:, :, l) / cv(:, :) ;
+      ENDDO
     CASE DEFAULT
-        field2=field
+      field2 = field
     END SELECT
 
-    ierr = nf90_put_var(nid,varid,field2,start,count)
+    ierr = nf90_put_var(nid, varid, field2, start, count)
     ierr = nf90_close(nid)
 
   END SUBROUTINE guide_out
-    
-  
-!===========================================================================
-  SUBROUTINE correctbid(iim,nl,x)
-    INTEGER iim,nl
-    REAL x(iim+1,nl)
-    INTEGER i,l
+
+
+  !===========================================================================
+  SUBROUTINE correctbid(iim, nl, x)
+    INTEGER iim, nl
+    REAL x(iim + 1, nl)
+    INTEGER i, l
     REAL zz
 
-    do l=1,nl
-        do i=2,iim-1
-            IF(abs(x(i,l))>1.e10) THEN
-               zz=0.5*(x(i-1,l)+x(i+1,l))
-              PRINT*,'correction ',i,l,x(i,l),zz
-               x(i,l)=zz
-            endif
-         enddo
-     enddo
+    do l = 1, nl
+      do i = 2, iim - 1
+        IF(abs(x(i, l))>1.e10) THEN
+          zz = 0.5 * (x(i - 1, l) + x(i + 1, l))
+          PRINT*, 'correction ', i, l, x(i, l), zz
+          x(i, l) = zz
+        endif
+      enddo
+    enddo
 
   END SUBROUTINE  correctbid
 
-!===========================================================================
+  !===========================================================================
 END MODULE guide_mod
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/iniacademic.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d/iniacademic.F90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d/iniacademic.F90	(revision 5118)
@@ -1,10 +1,9 @@
-
 ! $Id$
 
-SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
+SUBROUTINE iniacademic(vcov, ucov, teta, q, masse, ps, phis, time_0)
 
   USE lmdz_filtreg, ONLY: inifilr
-  USE infotrac,    ONLY: nqtot, niso, iqIsoPha, tracers, getKey, isoName
-  USE control_mod, ONLY: day_step,planet_type
+  USE infotrac, ONLY: nqtot, niso, iqIsoPha, tracers, getKey, isoName
+  USE control_mod, ONLY: day_step, planet_type
   USE exner_hyb_m, ONLY: exner_hyb
   USE exner_milieu_m, ONLY: exner_milieu
@@ -15,8 +14,9 @@
   USE comvert_mod, ONLY: ap, bp, preff, pa, presnivs, pressure_exner
   USE temps_mod, ONLY: annee_ref, day_ini, day_ref
-  USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
+  USE ener_mod, ONLY: etot0, ptot0, ztot0, stot0, ang0
   USE lmdz_readTracFiles, ONLY: addPhase
-  USE netcdf, ONLY: nf90_nowrite,nf90_open,nf90_noerr,nf90_inq_varid,nf90_close,nf90_get_var
+  USE netcdf, ONLY: nf90_nowrite, nf90_open, nf90_noerr, nf90_inq_varid, nf90_close, nf90_get_var
   USE lmdz_ran1, ONLY: ran1
+  USE lmdz_iniprint, ONLY: lunout, prt_level
 
   !   Author:    Frederic Hourdin      original: 15/01/93
@@ -33,40 +33,39 @@
   include "comgeom.h"
   include "academic.h"
-  include "iniprint.h"
 
   !   Arguments:
   !   ----------
 
-  REAL,INTENT(OUT) :: time_0
+  REAL, INTENT(OUT) :: time_0
 
   !   fields
-  REAL,INTENT(OUT) :: vcov(ip1jm,llm) ! meridional covariant wind
-  REAL,INTENT(OUT) :: ucov(ip1jmp1,llm) ! zonal covariant wind
-  REAL,INTENT(OUT) :: teta(ip1jmp1,llm) ! potential temperature (K)
-  REAL,INTENT(OUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers (.../kg_of_air)
-  REAL,INTENT(OUT) :: ps(ip1jmp1) ! surface pressure (Pa)
-  REAL,INTENT(OUT) :: masse(ip1jmp1,llm) ! air mass in grid cell (kg)
-  REAL,INTENT(OUT) :: phis(ip1jmp1) ! surface geopotential
+  REAL, INTENT(OUT) :: vcov(ip1jm, llm) ! meridional covariant wind
+  REAL, INTENT(OUT) :: ucov(ip1jmp1, llm) ! zonal covariant wind
+  REAL, INTENT(OUT) :: teta(ip1jmp1, llm) ! potential temperature (K)
+  REAL, INTENT(OUT) :: q(ip1jmp1, llm, nqtot) ! advected tracers (.../kg_of_air)
+  REAL, INTENT(OUT) :: ps(ip1jmp1) ! surface pressure (Pa)
+  REAL, INTENT(OUT) :: masse(ip1jmp1, llm) ! air mass in grid cell (kg)
+  REAL, INTENT(OUT) :: phis(ip1jmp1) ! surface geopotential
 
   !   Local:
   !   ------
 
-  REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
+  REAL p (ip1jmp1, llmp1)               ! pression aux interfac.des couches
   REAL pks(ip1jmp1)                      ! exner au  sol
-  REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
-  REAL phi(ip1jmp1,llm)                  ! geopotentiel
-  REAL ddsin,zsig,tetapv,w_pv  ! variables auxiliaires
+  REAL pk(ip1jmp1, llm)                   ! exner au milieu des couches
+  REAL phi(ip1jmp1, llm)                  ! geopotentiel
+  REAL ddsin, zsig, tetapv, w_pv  ! variables auxiliaires
   REAL tetastrat ! potential temperature in the stratosphere, in K
-  REAL tetajl(jjp1,llm)
-  INTEGER i,j,l,lsup,ij, iq, iName, iPhase, iqParent
-
-  INTEGER :: nid_relief,varid,ierr
-  REAL, DIMENSION(iip1,jjp1) :: relief
-
-  REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
-  REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
+  REAL tetajl(jjp1, llm)
+  INTEGER i, j, l, lsup, ij, iq, iName, iPhase, iqParent
+
+  INTEGER :: nid_relief, varid, ierr
+  REAL, DIMENSION(iip1, jjp1) :: relief
+
+  REAL teta0, ttp, delt_y, delt_z, eps ! Constantes pour profil de T
+  REAL k_f, k_c_a, k_c_s         ! Constantes de rappel
   LOGICAL ok_geost             ! Initialisation vent geost. ou nul
   LOGICAL ok_pv                ! Polar Vortex
-  REAL phi_pv,dphi_pv,gam_pv,tetanoise   ! Constantes pour polar vortex 
+  REAL phi_pv, dphi_pv, gam_pv, tetanoise   ! Constantes pour polar vortex
 
   REAL zz
@@ -74,18 +73,18 @@
 
   REAL zdtvr, tnat, alpha_ideal
-  LOGICAL,PARAMETER :: tnat1=.TRUE.
-  
-  CHARACTER(LEN=*),parameter :: modname="iniacademic"
-  CHARACTER(LEN=80) :: abort_message
+  LOGICAL, PARAMETER :: tnat1 = .TRUE.
+
+  CHARACTER(LEN = *), parameter :: modname = "iniacademic"
+  CHARACTER(LEN = 80) :: abort_message
 
   ! Sanity check: verify that options selected by user are not incompatible
   IF ((iflag_phys==1).AND. .NOT. read_start) THEN
-    WRITE(lunout,*) trim(modname)," error: if read_start is set to ", &
-    " false then iflag_phys should not be 1"
-    WRITE(lunout,*) "You most likely want an aquaplanet initialisation", &
-    " (iflag_phys >= 100)"
-    CALL abort_gcm(modname,"incompatible iflag_phys==1 and read_start==.FALSE.",1)
+    WRITE(lunout, *) trim(modname), " error: if read_start is set to ", &
+            " false then iflag_phys should not be 1"
+    WRITE(lunout, *) "You most likely want an aquaplanet initialisation", &
+            " (iflag_phys >= 100)"
+    CALL abort_gcm(modname, "incompatible iflag_phys==1 and read_start==.FALSE.", 1)
   ENDIF
-  
+
   !-----------------------------------------------------------------------
   ! 1. Initializations for Earth-like case
@@ -95,22 +94,22 @@
   CALL conf_planete
 
-  time_0=0.
-  day_ref=1
-  annee_ref=0
-
-  im         = iim
-  jm         = jjm
-  day_ini    = 1
-  dtvr    = daysec/REAL(day_step)
-  zdtvr=dtvr
-  etot0      = 0.
-  ptot0      = 0.
-  ztot0      = 0.
-  stot0      = 0.
-  ang0       = 0.
+  time_0 = 0.
+  day_ref = 1
+  annee_ref = 0
+
+  im = iim
+  jm = jjm
+  day_ini = 1
+  dtvr = daysec / REAL(day_step)
+  zdtvr = dtvr
+  etot0 = 0.
+  ptot0 = 0.
+  ztot0 = 0.
+  stot0 = 0.
+  ang0 = 0.
 
   IF (llm == 1) THEN
-     ! specific initializations for the shallow water case
-     kappa=1
+    ! specific initializations for the shallow water case
+    kappa = 1
   ENDIF
 
@@ -126,237 +125,237 @@
   IF (.NOT. read_start) THEN
 
-     !------------------------------------------------------------------
-     ! Lecture eventuelle d'un fichier de relief interpollee sur la grille
-     ! du modele.
-     ! On suppose que le fichier relief_in.nc est stoké sur une grille
-     ! iim*jjp1
-     ! Facile a créer à partir de la commande
-     ! cdo remapcon,fichier_output_phys.nc Relief.nc relief_in.nc
-     !------------------------------------------------------------------
-
-     relief=0.
-     ierr = nf90_open ('relief_in.nc', nf90_nowrite,nid_relief)
-     IF (ierr==nf90_noerr) THEN
-         ierr=nf90_inq_varid(nid_relief,'RELIEF',varid)
-         IF (ierr==nf90_noerr) THEN
-              ierr=nf90_get_var(nid_relief,varid,relief(1:iim,1:jjp1))
-              relief(iip1,:)=relief(1,:)
-         else
-              CALL abort_gcm ('iniacademic','variable RELIEF pas la',1)
-         endif
-     endif
-     ierr = nf90_close (nid_relief)
-
-     !------------------------------------------------------------------
-     ! Initialisation du geopotentiel au sol et de la pression
-     !------------------------------------------------------------------
-
-     PRINT*,'relief=',minval(relief),maxval(relief),'g=',g
-     do j=1,jjp1
-        do i=1,iip1
-           phis((j-1)*iip1+i)=g*relief(i,j)
+    !------------------------------------------------------------------
+    ! Lecture eventuelle d'un fichier de relief interpollee sur la grille
+    ! du modele.
+    ! On suppose que le fichier relief_in.nc est stoké sur une grille
+    ! iim*jjp1
+    ! Facile a créer à partir de la commande
+    ! cdo remapcon,fichier_output_phys.nc Relief.nc relief_in.nc
+    !------------------------------------------------------------------
+
+    relief = 0.
+    ierr = nf90_open ('relief_in.nc', nf90_nowrite, nid_relief)
+    IF (ierr==nf90_noerr) THEN
+      ierr = nf90_inq_varid(nid_relief, 'RELIEF', varid)
+      IF (ierr==nf90_noerr) THEN
+        ierr = nf90_get_var(nid_relief, varid, relief(1:iim, 1:jjp1))
+        relief(iip1, :) = relief(1, :)
+      else
+        CALL abort_gcm ('iniacademic', 'variable RELIEF pas la', 1)
+      endif
+    endif
+    ierr = nf90_close (nid_relief)
+
+    !------------------------------------------------------------------
+    ! Initialisation du geopotentiel au sol et de la pression
+    !------------------------------------------------------------------
+
+    PRINT*, 'relief=', minval(relief), maxval(relief), 'g=', g
+    do j = 1, jjp1
+      do i = 1, iip1
+        phis((j - 1) * iip1 + i) = g * relief(i, j)
+      enddo
+    enddo
+    PRINT*, 'phis=', minval(phis), maxval(phis), 'g=', g
+
+    ! ground geopotential
+    !phis(:)=0.
+    ps(:) = preff
+    CALL pression (ip1jmp1, ap, bp, ps, p)
+
+    IF (pressure_exner) THEN
+      CALL exner_hyb(ip1jmp1, ps, p, pks, pk)
+    else
+      CALL exner_milieu(ip1jmp1, ps, p, pks, pk)
+    endif
+    CALL massdair(p, masse)
+  ENDIF
+
+  IF (llm == 1) THEN
+    ! initialize fields for the shallow water case, if required
+    IF (.NOT.read_start) THEN
+      phis(:) = 0.
+      q(:, :, :) = 0
+      CALL sw_case_williamson91_6(vcov, ucov, teta, masse, ps)
+    endif
+  ENDIF
+
+  academic_case: if (iflag_phys >= 2) THEN
+    ! initializations
+
+    ! 1. local parameters
+    ! by convention, winter is in the southern hemisphere
+    ! Geostrophic wind or no wind?
+    ok_geost = .TRUE.
+    CALL getin('ok_geost', ok_geost)
+    ! Constants for Newtonian relaxation and friction
+    k_f = 1.                !friction
+    CALL getin('k_j', k_f)
+    k_f = 1. / (daysec * k_f)
+    k_c_s = 4.  !cooling surface
+    CALL getin('k_c_s', k_c_s)
+    k_c_s = 1. / (daysec * k_c_s)
+    k_c_a = 40. !cooling free atm
+    CALL getin('k_c_a', k_c_a)
+    k_c_a = 1. / (daysec * k_c_a)
+    ! Constants for Teta equilibrium profile
+    teta0 = 315.     ! mean Teta (S.H. 315K)
+    CALL getin('teta0', teta0)
+    ttp = 200.       ! Tropopause temperature (S.H. 200K)
+    CALL getin('ttp', ttp)
+    eps = 0.         ! Deviation to N-S symmetry(~0-20K)
+    CALL getin('eps', eps)
+    delt_y = 60.     ! Merid Temp. Gradient (S.H. 60K)
+    CALL getin('delt_y', delt_y)
+    delt_z = 10.     ! Vertical Gradient (S.H. 10K)
+    CALL getin('delt_z', delt_z)
+    ! Polar vortex
+    ok_pv = .FALSE.
+    CALL getin('ok_pv', ok_pv)
+    phi_pv = -50.            ! Latitude of edge of vortex
+    CALL getin('phi_pv', phi_pv)
+    phi_pv = phi_pv * pi / 180.
+    dphi_pv = 5.             ! Width of the edge
+    CALL getin('dphi_pv', dphi_pv)
+    dphi_pv = dphi_pv * pi / 180.
+    gam_pv = 4.              ! -dT/dz vortex (in K/km)
+    CALL getin('gam_pv', gam_pv)
+    tetanoise = 0.005
+    CALL getin('tetanoise', tetanoise)
+
+    ! 2. Initialize fields towards which to relax
+    ! Friction
+    knewt_g = k_c_a
+    DO l = 1, llm
+      zsig = presnivs(l) / preff
+      knewt_t(l) = (k_c_s - k_c_a) * MAX(0., (zsig - 0.7) / 0.3)
+      kfrict(l) = k_f * MAX(0., (zsig - 0.7) / 0.3)
+    ENDDO
+    DO j = 1, jjp1
+      clat4((j - 1) * iip1 + 1:j * iip1) = cos(rlatu(j))**4
+    ENDDO
+
+    ! Potential temperature
+    DO l = 1, llm
+      zsig = presnivs(l) / preff
+      tetastrat = ttp * zsig**(-kappa)
+      tetapv = tetastrat
+      IF ((ok_pv).AND.(zsig<0.1)) THEN
+        tetapv = tetastrat * (zsig * 10.)**(kappa * cpp * gam_pv / 1000. / g)
+      ENDIF
+      DO j = 1, jjp1
+        ! Troposphere
+        ddsin = sin(rlatu(j))
+        tetajl(j, l) = teta0 - delt_y * ddsin * ddsin + eps * ddsin &
+                - delt_z * (1. - ddsin * ddsin) * log(zsig)
+        IF (planet_type=="giant") THEN
+          tetajl(j, l) = teta0 + (delt_y * &
+                  ((sin(rlatu(j) * 3.14159 * eps + 0.0001))**2)   &
+                  / ((rlatu(j) * 3.14159 * eps + 0.0001)**2))     &
+                  - delt_z * log(zsig)
+        endif
+        ! Profil stratospherique isotherme (+vortex)
+        w_pv = (1. - tanh((rlatu(j) - phi_pv) / dphi_pv)) / 2.
+        tetastrat = tetastrat * (1. - w_pv) + tetapv * w_pv
+        tetajl(j, l) = MAX(tetajl(j, l), tetastrat)
+      ENDDO
+    ENDDO
+
+    !          CALL writefield('theta_eq',tetajl)
+
+    do l = 1, llm
+      do j = 1, jjp1
+        do i = 1, iip1
+          ij = (j - 1) * iip1 + i
+          tetarappel(ij, l) = tetajl(j, l)
         enddo
-     enddo
-     PRINT*,'phis=',minval(phis),maxval(phis),'g=',g
-
-     ! ground geopotential
-     !phis(:)=0.
-     ps(:)=preff
-     CALL pression ( ip1jmp1, ap, bp, ps, p       )
-
-     IF (pressure_exner) THEN
-       CALL exner_hyb( ip1jmp1, ps, p, pks, pk)
-     else
-       CALL exner_milieu(ip1jmp1,ps,p,pks,pk)
-     endif
-     CALL massdair(p,masse)
-  ENDIF
-
-  IF (llm == 1) THEN
-     ! initialize fields for the shallow water case, if required
-     IF (.NOT.read_start) THEN
-        phis(:)=0.
-        q(:,:,:)=0
-        CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
-     endif
-  ENDIF
-
-  academic_case: if (iflag_phys >= 2) THEN
-     ! initializations
-
-     ! 1. local parameters
-     ! by convention, winter is in the southern hemisphere
-     ! Geostrophic wind or no wind?
-     ok_geost=.TRUE.
-     CALL getin('ok_geost',ok_geost)
-     ! Constants for Newtonian relaxation and friction
-     k_f=1.                !friction 
-     CALL getin('k_j',k_f)
-     k_f=1./(daysec*k_f)
-     k_c_s=4.  !cooling surface
-     CALL getin('k_c_s',k_c_s)
-     k_c_s=1./(daysec*k_c_s)
-     k_c_a=40. !cooling free atm
-     CALL getin('k_c_a',k_c_a)
-     k_c_a=1./(daysec*k_c_a)
-     ! Constants for Teta equilibrium profile
-     teta0=315.     ! mean Teta (S.H. 315K)
-     CALL getin('teta0',teta0)
-     ttp=200.       ! Tropopause temperature (S.H. 200K)
-     CALL getin('ttp',ttp)
-     eps=0.         ! Deviation to N-S symmetry(~0-20K)
-     CALL getin('eps',eps)
-     delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
-     CALL getin('delt_y',delt_y)
-     delt_z=10.     ! Vertical Gradient (S.H. 10K)
-     CALL getin('delt_z',delt_z)
-     ! Polar vortex
-     ok_pv=.FALSE.
-     CALL getin('ok_pv',ok_pv)
-     phi_pv=-50.            ! Latitude of edge of vortex
-     CALL getin('phi_pv',phi_pv)
-     phi_pv=phi_pv*pi/180.
-     dphi_pv=5.             ! Width of the edge
-     CALL getin('dphi_pv',dphi_pv)
-     dphi_pv=dphi_pv*pi/180.
-     gam_pv=4.              ! -dT/dz vortex (in K/km)
-     CALL getin('gam_pv',gam_pv)
-     tetanoise=0.005
-     CALL getin('tetanoise',tetanoise)
-
-     ! 2. Initialize fields towards which to relax
-     ! Friction
-     knewt_g=k_c_a
-     DO l=1,llm
-        zsig=presnivs(l)/preff
-        knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
-        kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
-     ENDDO
-     DO j=1,jjp1
-        clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
-     ENDDO
-
-     ! Potential temperature 
-     DO l=1,llm
-        zsig=presnivs(l)/preff
-        tetastrat=ttp*zsig**(-kappa)
-        tetapv=tetastrat
-        IF ((ok_pv).AND.(zsig<0.1)) THEN
-           tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
-        ENDIF
-        DO j=1,jjp1
-           ! Troposphere
-           ddsin=sin(rlatu(j))
-           tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin &
-                -delt_z*(1.-ddsin*ddsin)*log(zsig)
-           IF (planet_type=="giant") THEN
-             tetajl(j,l)=teta0+(delt_y*                   &
-                ((sin(rlatu(j)*3.14159*eps+0.0001))**2)   &
-                / ((rlatu(j)*3.14159*eps+0.0001)**2))     &
-                -delt_z*log(zsig)
-           endif
-           ! Profil stratospherique isotherme (+vortex)
-           w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
-           tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
-           tetajl(j,l)=MAX(tetajl(j,l),tetastrat)  
-        ENDDO
-     ENDDO
-
-     !          CALL writefield('theta_eq',tetajl)
-
-     do l=1,llm
-        do j=1,jjp1
-           do i=1,iip1
-              ij=(j-1)*iip1+i
-              tetarappel(ij,l)=tetajl(j,l)
-           enddo
+      enddo
+    enddo
+
+    ! 3. Initialize fields (if necessary)
+    IF (.NOT. read_start) THEN
+      ! bulk initialization of temperature
+      IF (iflag_phys>10000) THEN
+        ! Particular case to impose a constant temperature T0=0.01*iflag_physx
+        teta(:, :) = 0.01 * iflag_phys / (pk(:, :) / cpp)
+      ELSE
+        teta(:, :) = tetarappel(:, :)
+      ENDIF
+      ! geopotential
+      CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
+
+      DO l = 1, llm
+        PRINT*, 'presnivs,play,l', presnivs(l), (pk(1, l) / cpp)**(1. / kappa) * preff
+        !pks(ij) = (cpp/preff) * ps(ij)
+        !pk(ij,1) = .5*pks(ij)
+        ! pk = cpp * (p/preff)^kappa
+      ENDDO
+
+      ! winds
+      IF (ok_geost) THEN
+        CALL ugeostr(phi, ucov)
+      else
+        ucov(:, :) = 0.
+      endif
+      vcov(:, :) = 0.
+
+      ! bulk initialization of tracers
+      IF (planet_type=="earth") THEN
+        ! Earth: first two tracers will be water
+        do iq = 1, nqtot
+          q(:, :, iq) = 0.
+          IF(tracers(iq)%name == addPhase('H2O', 'g')) q(:, :, iq) = 1.e-10
+          IF(tracers(iq)%name == addPhase('H2O', 'l')) q(:, :, iq) = 1.e-15
+
+          ! CRisi: init des isotopes
+          ! distill de Rayleigh très simplifiée
+          iName = tracers(iq)%iso_iName
+          IF (niso <= 0 .OR. iName <= 0) CYCLE
+          iPhase = tracers(iq)%iso_iPhase
+          iqParent = tracers(iq)%iqParent
+          IF(tracers(iq)%iso_iZone == 0) THEN
+            IF (tnat1) THEN
+              tnat = 1.0
+              alpha_ideal = 1.0
+              WRITE(*, *) 'Attention dans iniacademic: alpha_ideal=1'
+            else
+              IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
+                      CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
+            endif
+            q(:, :, iq) = q(:, :, iqParent) * tnat * (q(:, :, iqParent) / 30.e-3)**(alpha_ideal - 1.)
+          ELSE !IF(tracers(iq)%iso_iZone == 0) THEN
+            IF(tracers(iq)%iso_iZone == 1) THEN
+              ! correction le 14 mai 2024 pour que tous les traceurs soient de la couleur 1.
+              ! Sinon, on va avoir des porblèmes de conservation de masse de traceurs.
+              q(:, :, iq) = q(:, :, iqIsoPha(iName, iPhase))
+            else !IF(tracers(iq)%iso_iZone == 1) THEN
+              q(:, :, iq) = 0.
+            endif !IF(tracers(iq)%iso_iZone == 1) THEN
+          END IF !IF(tracers(iq)%iso_iZone == 0) THEN
         enddo
-     enddo
-
-     ! 3. Initialize fields (if necessary)
-     IF (.NOT. read_start) THEN
-        ! bulk initialization of temperature
-        IF (iflag_phys>10000) THEN
-        ! Particular case to impose a constant temperature T0=0.01*iflag_physx
-           teta(:,:)= 0.01*iflag_phys/(pk(:,:)/cpp)
-        ELSE
-           teta(:,:)=tetarappel(:,:)
-        ENDIF
-        ! geopotential
-        CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
-
-        DO l=1,llm
-          PRINT*,'presnivs,play,l',presnivs(l),(pk(1,l)/cpp)**(1./kappa)*preff
-         !pks(ij) = (cpp/preff) * ps(ij)
-         !pk(ij,1) = .5*pks(ij)
-         ! pk = cpp * (p/preff)^kappa
-        ENDDO
-
-        ! winds
-        IF (ok_geost) THEN
-           CALL ugeostr(phi,ucov)
-        else
-           ucov(:,:)=0.
-        endif
-        vcov(:,:)=0.
-
-        ! bulk initialization of tracers
-        IF (planet_type=="earth") THEN
-           ! Earth: first two tracers will be water
-           do iq=1,nqtot
-              q(:,:,iq)=0.
-              IF(tracers(iq)%name == addPhase('H2O', 'g')) q(:,:,iq)=1.e-10
-              IF(tracers(iq)%name == addPhase('H2O', 'l')) q(:,:,iq)=1.e-15
-
-              ! CRisi: init des isotopes
-              ! distill de Rayleigh très simplifiée
-              iName    = tracers(iq)%iso_iName
-              IF (niso <= 0 .OR. iName <= 0) CYCLE
-              iPhase   = tracers(iq)%iso_iPhase
-              iqParent = tracers(iq)%iqParent
-              IF(tracers(iq)%iso_iZone == 0) THEN
-                 IF (tnat1) THEN
-                         tnat=1.0
-                         alpha_ideal=1.0
-                         WRITE(*,*) 'Attention dans iniacademic: alpha_ideal=1'
-                 else
-                    IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
-                    CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
-                 endif
-                 q(:,:,iq) = q(:,:,iqParent)*tnat*(q(:,:,iqParent)/30.e-3)**(alpha_ideal-1.)
-              ELSE !IF(tracers(iq)%iso_iZone == 0) THEN
-                 IF(tracers(iq)%iso_iZone == 1) THEN
-                    ! correction le 14 mai 2024 pour que tous les traceurs soient de la couleur 1.
-                    ! Sinon, on va avoir des porblèmes de conservation de masse de traceurs.
-                    q(:,:,iq) = q(:,:,iqIsoPha(iName,iPhase))
-                 else !IF(tracers(iq)%iso_iZone == 1) THEN
-                    q(:,:,iq) = 0.
-                 endif !IF(tracers(iq)%iso_iZone == 1) THEN
-              END IF !IF(tracers(iq)%iso_iZone == 0) THEN
-           enddo
-        else
-           q(:,:,:)=0
-        endif ! of if (planet_type=="earth")
-
-        CALL check_isotopes_seq(q,1,ip1jmp1,'iniacademic_loc')
-
-        ! add random perturbation to temperature
-        idum  = -1
-        zz = ran1(idum)
-        idum  = 0
-        do l=1,llm
-           do ij=iip2,ip1jm
-              teta(ij,l)=teta(ij,l)*(1.+tetanoise*ran1(idum))
-           enddo
+      else
+        q(:, :, :) = 0
+      endif ! of if (planet_type=="earth")
+
+      CALL check_isotopes_seq(q, 1, ip1jmp1, 'iniacademic_loc')
+
+      ! add random perturbation to temperature
+      idum = -1
+      zz = ran1(idum)
+      idum = 0
+      do l = 1, llm
+        do ij = iip2, ip1jm
+          teta(ij, l) = teta(ij, l) * (1. + tetanoise * ran1(idum))
         enddo
-
-        ! maintain periodicity in longitude
-        do l=1,llm
-           do ij=1,ip1jmp1,iip1
-              teta(ij+iim,l)=teta(ij,l)
-           enddo
+      enddo
+
+      ! maintain periodicity in longitude
+      do l = 1, llm
+        do ij = 1, ip1jmp1, iip1
+          teta(ij + iim, l) = teta(ij, l)
         enddo
-
-     ENDIF ! of IF (.NOT. read_start)
+      enddo
+
+    ENDIF ! of IF (.NOT. read_start)
   ENDIF academic_case
 
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/integrd.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d/integrd.F90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d/integrd.F90	(revision 5118)
@@ -11,4 +11,5 @@
   USE comvert_mod, ONLY: ap, bp
   USE temps_mod, ONLY: dt
+  USE lmdz_iniprint, ONLY: lunout, prt_level
 
   IMPLICIT NONE
@@ -33,5 +34,4 @@
   include "paramet.h"
   include "comgeom.h"
-  include "iniprint.h"
 
   !   Arguments:
@@ -250,4 +250,3 @@
   END IF
 
-
 END SUBROUTINE integrd
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/leapfrog.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d/leapfrog.F90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d/leapfrog.F90	(revision 5118)
@@ -26,4 +26,5 @@
   USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_PHYS
   USE lmdz_description, ONLY: descript
+  USE lmdz_iniprint, ONLY: lunout, prt_level
 
   IMPLICIT NONE
@@ -64,5 +65,4 @@
   include "comdissnew.h"
   include "comgeom.h"
-  include "iniprint.h"
   include "academic.h"
 
@@ -227,8 +227,8 @@
   !   -----
 
-  jD_cur = jD_ref + day_ini - day_ref +                             &
-  (itau+1)/day_step
-  jH_cur = jH_ref + start_time +                                    &
-          mod(itau+1, day_step)/float(day_step)
+  jD_cur = jD_ref + day_ini - day_ref + &
+          (itau + 1) / day_step
+  jH_cur = jH_ref + start_time + &
+          mod(itau + 1, day_step) / float(day_step)
   jD_cur = jD_cur + int(jH_cur)
   jH_cur = jH_cur - int(jH_cur)
@@ -237,5 +237,5 @@
 
   IF (ok_guide) THEN
-    CALL guide_main(itau,ucov,vcov,teta,q,masse,ps)
+    CALL guide_main(itau, ucov, vcov, teta, q, masse, ps)
   ENDIF
 
@@ -346,7 +346,6 @@
       !maf stokage du flux de masse pour traceurs OFF-LINE
 
-       CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis, &
-             dtvr, itau)
-
+      CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
+              dtvr, itau)
 
     ENDIF ! of IF (offline)
@@ -405,6 +404,6 @@
     !       jH_cur = jH_ref +                                            &
     ! &              (itau * dtvr / daysec - int(itau * dtvr / daysec))
-    jD_cur = jD_ref + day_ini - day_ref +                        &
-    (itau+1)/day_step
+    jD_cur = jD_ref + day_ini - day_ref + &
+            (itau + 1) / day_step
 
     IF (planet_type =="generic") THEN
@@ -413,6 +412,6 @@
     ENDIF
 
-    jH_cur = jH_ref + start_time +                               &
-            mod(itau+1, day_step)/float(day_step)
+    jH_cur = jH_ref + start_time + &
+            mod(itau + 1, day_step) / float(day_step)
     jD_cur = jD_cur + int(jH_cur)
     jH_cur = jH_cur - int(jH_cur)
@@ -650,10 +649,10 @@
 
       IF (ok_dynzon) THEN
-             CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav, &
-                   ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
+        CALL bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &
+                ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)
       END IF
       IF (ok_dyn_ave) THEN
-             CALL writedynav(itau,vcov, &
-                   ucov,teta,pk,phi,q,masse,ps,phis)
+        CALL writedynav(itau, vcov, &
+                ucov, teta, pk, phi, q, masse, ps, phis)
       ENDIF
 
@@ -675,13 +674,13 @@
           vnat(:, l) = vcov(:, l) / cv(:)
         enddo
-          IF (ok_dyn_ins) THEN
-            ! WRITE(lunout,*) "leapfrog: CALL writehist, itau=",itau
-           CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
-            ! CALL WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
-            ! CALL WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
-           ! CALL WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
-           !  CALL WriteField('ps',reshape(ps,(/iip1,jmp1/)))
-           !  CALL WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
-          endif ! of if (ok_dyn_ins)
+        IF (ok_dyn_ins) THEN
+          ! WRITE(lunout,*) "leapfrog: CALL writehist, itau=",itau
+          CALL writehist(itau, vcov, ucov, teta, phi, q, masse, ps, phis)
+          ! CALL WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
+          ! CALL WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
+          ! CALL WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
+          !  CALL WriteField('ps',reshape(ps,(/iip1,jmp1/)))
+          !  CALL WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
+        endif ! of if (ok_dyn_ins)
         ! For some Grads outputs of fields
         IF (output_grads_dyn) THEN
@@ -781,10 +780,10 @@
 
         IF (ok_dynzon) THEN
-             CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav, &
-                   ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
+          CALL bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &
+                  ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)
         ENDIF
         IF (ok_dyn_ave) THEN
-             CALL writedynav(itau,vcov, &
-                   ucov,teta,pk,phi,q,masse,ps,phis)
+          CALL writedynav(itau, vcov, &
+                  ucov, teta, pk, phi, q, masse, ps, phis)
         ENDIF
 
@@ -799,9 +798,9 @@
           vnat(:, l) = vcov(:, l) / cv(:)
         enddo
-          IF (ok_dyn_ins) THEN
-             ! WRITE(lunout,*) "leapfrog: CALL writehist (b)",
-  ! &                        itau,iecri
-            CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
-          endif ! of if (ok_dyn_ins)
+        IF (ok_dyn_ins) THEN
+          ! WRITE(lunout,*) "leapfrog: CALL writehist (b)",
+          ! &                        itau,iecri
+          CALL writehist(itau, vcov, ucov, teta, phi, q, masse, ps, phis)
+        endif ! of if (ok_dyn_ins)
         ! For some Grads outputs
         IF (output_grads_dyn) THEN
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/sw_case_williamson91_6.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d/sw_case_williamson91_6.F90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d/sw_case_williamson91_6.F90	(revision 5118)
@@ -27,4 +27,5 @@
   USE comconst_mod, ONLY: cpp, omeg, rad
   USE comvert_mod, ONLY: ap, bp, preff
+  USE lmdz_iniprint, ONLY: lunout, prt_level
 
   IMPLICIT NONE
@@ -36,5 +37,4 @@
   include "paramet.h"
   include "comgeom.h"
-  include "iniprint.h"
 
   !   Arguments:
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/top_bound.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d/top_bound.F90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d/top_bound.F90	(revision 5118)
@@ -6,4 +6,5 @@
           tau_top_bound
   USE comvert_mod, ONLY: presnivs, preff, scaleheight
+  USE lmdz_iniprint, ONLY: lunout, prt_level
 
   IMPLICIT NONE
@@ -53,5 +54,4 @@
 
   include "comdissipn.h"
-  include "iniprint.h"
 
   !   Arguments:
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/vlsplt.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d/vlsplt.F90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d/vlsplt.F90	(revision 5118)
@@ -107,4 +107,5 @@
   USE infotrac, ONLY: nqtot,tracers, & ! CRisi
         min_qParent,min_qMass,min_ratio ! MVals et CRisi
+  USE lmdz_iniprint, ONLY: lunout, prt_level
 
   ! Auteurs:   P.Le Van, F.Hourdin, F.Forget
@@ -121,5 +122,4 @@
   include "dimensions.h"
   include "paramet.h"
-  include "iniprint.h"
   !
   !
