Index: LMDZ4/trunk/libf/dyn3d/adaptdt.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/adaptdt.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/adaptdt.F	(revision 1403)
@@ -1,8 +1,9 @@
 !
-! $Header$
+! $Id$
 !
       subroutine adaptdt(nadv,dtbon,n,pbaru,
      c                   masse)
 
+      USE control_mod
       IMPLICIT NONE
 
@@ -16,5 +17,4 @@
 #include "logic.h"
 #include "temps.h"
-#include "control.h"
 #include "ener.h"
 #include "description.h"
Index: LMDZ4/trunk/libf/dyn3d/advtrac.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/advtrac.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/advtrac.F	(revision 1403)
@@ -16,4 +16,6 @@
 c
       USE infotrac
+      USE control_mod
+ 
 
       IMPLICIT NONE
@@ -27,5 +29,4 @@
 #include "logic.h"
 #include "temps.h"
-#include "control.h"
 #include "ener.h"
 #include "description.h"
@@ -121,5 +122,5 @@
 
       ! ... Flux de masse diaganostiques traceurs
-      flxw = wg / FLOAT(iapp_tracvl)
+      flxw = wg / REAL(iapp_tracvl)
 
 c  test sur l'eventuelle creation de valeurs negatives de la masse
Index: LMDZ4/trunk/libf/dyn3d/bilan_dyn.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/bilan_dyn.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/bilan_dyn.F	(revision 1403)
@@ -423,5 +423,5 @@
          Q_cum(:,:,:,iQ)=Q_cum(:,:,:,iQ)/masse_cum(:,:,:)
       enddo
-      zz=1./float(ncum)
+      zz=1./REAL(ncum)
       ps_cum=ps_cum*zz
       masse_cum=masse_cum*zz
Index: LMDZ4/trunk/libf/dyn3d/caladvtrac.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/caladvtrac.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/caladvtrac.F	(revision 1403)
@@ -1,4 +1,4 @@
 !
-! $Header$
+! $Id$
 !
 c
@@ -9,4 +9,6 @@
 c
       USE infotrac
+      USE control_mod
+ 
       IMPLICIT NONE
 c
@@ -24,5 +26,4 @@
 #include "paramet.h"
 #include "comconst.h"
-#include "control.h"
 
 c   Arguments:
Index: LMDZ4/trunk/libf/dyn3d/calfis.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/calfis.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/calfis.F	(revision 1403)
@@ -31,4 +31,6 @@
 c   .........
       USE infotrac
+      USE control_mod
+ 
 
       IMPLICIT NONE
@@ -96,5 +98,5 @@
 #include "comvert.h"
 #include "comgeom2.h"
-#include "control.h"
+#include "iniprint.h"
 
 c    Arguments :
@@ -149,4 +151,10 @@
       REAL zdpsrf(ngridmx)
 c
+      REAL zdufic(ngridmx,llm),zdvfic(ngridmx,llm)
+      REAL zdtfic(ngridmx,llm),zdqfic(ngridmx,llm,nqtot)
+      REAL jH_cur_split,zdt_split
+      LOGICAL debut_split,lafin_split
+      INTEGER isplit
+
       REAL zsin(iim),zcos(iim),z1(iim)
       REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
@@ -181,8 +189,9 @@
         debut = .TRUE.
         IF (ngridmx.NE.2+(jjm-1)*iim) THEN
-         PRINT*,'STOP dans calfis'
-         PRINT*,'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
-         PRINT*,'  ngridmx  jjm   iim   '
-         PRINT*,ngridmx,jjm,iim
+         write(lunout,*) 'STOP dans calfis'
+         write(lunout,*)
+     &   'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
+         write(lunout,*) '  ngridmx  jjm   iim   '
+         write(lunout,*) ngridmx,jjm,iim
          STOP
         ENDIF
@@ -308,7 +317,7 @@
       CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
       DO l=1,llm
-	 DO ig=1,ngridmx
-	   zphi(ig,l)=zphi(ig,l)-zphis(ig)
-	 ENDDO
+         DO ig=1,ngridmx
+           zphi(ig,l)=zphi(ig,l)-zphis(ig)
+         ENDDO
       ENDDO
 
@@ -408,5 +417,5 @@
             z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
             z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
-	 ENDDO
+         ENDDO
 
          DO i=1,iim
@@ -415,5 +424,5 @@
             zsin(i)    = SIN(rlonv(i))*z1(i)
             zsinbis(i) = SIN(rlonv(i))*z1bis(i)
-	 ENDDO
+         ENDDO
 
          zufi(ngridmx,l)  = SSUM(iim,zcos,1)/pi
@@ -443,11 +452,25 @@
       if (planet_type=="earth") then
 #ifdef CPP_EARTH
-      CALL physiq (ngridmx,
+
+      write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
+      zdt_split=dtphys/nsplit_phys
+      zdufic(:,:)=0.
+      zdvfic(:,:)=0.
+      zdtfic(:,:)=0.
+      zdqfic(:,:,:)=0.
+
+      do isplit=1,nsplit_phys
+
+         jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
+         debut_split=debut.and.isplit==1
+         lafin_split=lafin.and.isplit==nsplit_phys
+
+         CALL physiq (ngridmx,
      .             llm,
-     .             debut,
-     .             lafin,
+     .             debut_split,
+     .             lafin_split,
      .             jD_cur,
-     .             jH_cur,
-     .             dtphys,
+     .             jH_cur_split,
+     .             zdt_split,
      .             zplev,
      .             zplay,
@@ -469,4 +492,21 @@
      .             pducov,
      .             PVteta)
+
+         zufi(:,:)=zufi(:,:)+zdufi(:,:)*zdt_split
+         zvfi(:,:)=zvfi(:,:)+zdvfi(:,:)*zdt_split
+         ztfi(:,:)=ztfi(:,:)+zdtfi(:,:)*zdt_split
+         zqfi(:,:,:)=zqfi(:,:,:)+zdqfi(:,:,:)*zdt_split
+
+         zdufic(:,:)=zdufic(:,:)+zdufi(:,:)
+         zdvfic(:,:)=zdvfic(:,:)+zdvfi(:,:)
+         zdtfic(:,:)=zdtfic(:,:)+zdtfi(:,:)
+         zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:)
+
+      enddo
+      zdufi(:,:)=zdufic(:,:)/nsplit_phys
+      zdvfi(:,:)=zdvfic(:,:)/nsplit_phys
+      zdtfi(:,:)=zdtfic(:,:)/nsplit_phys
+      zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys
+
 #endif
       endif !of if (planet_type=="earth")
Index: LMDZ4/trunk/libf/dyn3d/ce0l.F90
===================================================================
--- LMDZ4/trunk/libf/dyn3d/ce0l.F90	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/ce0l.F90	(revision 1403)
@@ -14,4 +14,5 @@
 !     masque is created in etat0, passed to limit to ensure consistancy.
 !-------------------------------------------------------------------------------
+  USE control_mod
 #ifdef CPP_EARTH
 ! This prog. is designed to work for Earth
@@ -36,5 +37,4 @@
 #include "indicesol.h"
 #include "iniprint.h"
-#include "control.h"
 #include "temps.h"
 #include "logic.h"
Index: LMDZ4/trunk/libf/dyn3d/com_io_dyn.h
===================================================================
--- LMDZ4/trunk/libf/dyn3d/com_io_dyn.h	(revision 1340)
+++ 	(revision )
@@ -1,5 +1,0 @@
-!
-! $Header$
-!
-      integer histid, histvid, histaveid
-      common/com_io_dyn/histid, histvid, histaveid
Index: LMDZ4/trunk/libf/dyn3d/com_io_dyn_mod.F90
===================================================================
--- LMDZ4/trunk/libf/dyn3d/com_io_dyn_mod.F90	(revision 1403)
+++ LMDZ4/trunk/libf/dyn3d/com_io_dyn_mod.F90	(revision 1403)
@@ -0,0 +1,31 @@
+!
+! $Id $
+!
+module com_io_dyn_mod
+
+  implicit none 
+
+! Names of various files for outputs (in the dynamics)
+  ! to store instantaneous values:
+  character(len=18),parameter :: dynhist_file="dyn_hist.nc" ! on scalar grid
+  character(len=18),parameter :: dynhistv_file="dyn_histv.nc" ! on v grid
+  character(len=18),parameter :: dynhistu_file="dyn_histu.nc" ! on u grid
+
+  ! to store averaged values:
+  character(len=18),parameter :: dynhistave_file="dyn_hist_ave.nc"
+  character(len=18),parameter :: dynhistvave_file="dyn_histv_ave.nc"
+  character(len=18),parameter :: dynhistuave_file="dyn_histu_ave.nc"
+  
+! Ids of various files for outputs (in the dynamics)
+
+  ! instantaneous (these are set by inithist.F)
+  integer :: histid
+  integer :: histvid
+  integer :: histuid
+  
+  ! averages (these are set by initdynav.F)
+  integer :: histaveid
+  integer :: histvaveid
+  integer :: histuaveid
+  
+end module com_io_dyn_mod
Index: LMDZ4/trunk/libf/dyn3d/conf_gcm.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/conf_gcm.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/conf_gcm.F	(revision 1403)
@@ -6,4 +6,5 @@
       SUBROUTINE conf_gcm( tapedef, etatinit, clesphy0 )
 c
+      USE control_mod
 #ifdef CPP_IOIPSL
       use IOIPSL
@@ -34,5 +35,4 @@
 #include "dimensions.h"
 #include "paramet.h"
-#include "control.h"
 #include "logic.h"
 #include "serre.h"
@@ -162,4 +162,12 @@
        day_step = 240 
        CALL getin('day_step',day_step)
+
+!Config  Key  = nsplit_phys
+!Config  Desc = nombre de pas par jour
+!Config  Def  = 1 
+!Config  Help = nombre de pas par jour (multiple de iperiod) (
+!Config          ici pour  dt = 1 min ) 
+       nsplit_phys = 1 
+       CALL getin('nsplit_phys',nsplit_phys)
 
 !Config  Key  = iperiod
@@ -573,5 +581,4 @@
       CALL getin('config_inca',config_inca)
 
-
 !Config  Key  = ok_dynzon 
 !Config  Desc = calcul et sortie des transports 
@@ -581,4 +588,21 @@
       ok_dynzon = .FALSE. 
       CALL getin('ok_dynzon',ok_dynzon) 
+
+!Config  Key  = ok_dyn_ins
+!Config  Desc = sorties instantanees dans la dynamique
+!Config  Def  = n 
+!Config  Help = 
+!Config          
+      ok_dyn_ins = .FALSE. 
+      CALL getin('ok_dyn_ins',ok_dyn_ins) 
+
+!Config  Key  = ok_dyn_ave
+!Config  Desc = sorties moyennes dans la dynamique
+!Config  Def  = n 
+!Config  Help = 
+!Config          
+      ok_dyn_ave = .FALSE. 
+      CALL getin('ok_dyn_ave',ok_dyn_ave) 
+
 
       write(lunout,*)' #########################################'
@@ -620,4 +644,6 @@
       write(lunout,*)' config_inca = ', config_inca
       write(lunout,*)' ok_dynzon = ', ok_dynzon 
+      write(lunout,*)' ok_dyn_ins = ', ok_dyn_ins 
+      write(lunout,*)' ok_dyn_ave = ', ok_dyn_ave 
 
       RETURN
@@ -746,10 +772,26 @@
 
 !Config  Key  = ok_dynzon 
-!Config  Desc = calcul et sortie des transports 
+!Config  Desc = sortie des transports zonaux dans la dynamique
 !Config  Def  = n 
-!Config  Help = Permet de mettre en route le calcul des transports 
+!Config  Help =  
 !Config          
        ok_dynzon = .FALSE.
        CALL getin('ok_dynzon',ok_dynzon) 
+
+!Config  Key  = ok_dyn_ins
+!Config  Desc = sorties instantanees dans la dynamique
+!Config  Def  = n 
+!Config  Help = 
+!Config          
+      ok_dyn_ins = .FALSE. 
+      CALL getin('ok_dyn_ins',ok_dyn_ins) 
+
+!Config  Key  = ok_dyn_ave
+!Config  Desc = sorties moyennes dans la dynamique
+!Config  Def  = n 
+!Config  Help = 
+!Config          
+      ok_dyn_ave = .FALSE. 
+      CALL getin('ok_dyn_ave',ok_dyn_ave) 
 
 !Config key = ok_strato
@@ -824,4 +866,6 @@
       write(lunout,*)' config_inca = ', config_inca
       write(lunout,*)' ok_dynzon = ', ok_dynzon
+      write(lunout,*)' ok_dyn_ins = ', ok_dyn_ins 
+      write(lunout,*)' ok_dyn_ave = ', ok_dyn_ave 
       write(lunout,*)' ok_strato = ', ok_strato
       write(lunout,*)' ok_gradsfile = ', ok_gradsfile
Index: LMDZ4/trunk/libf/dyn3d/control.h
===================================================================
--- LMDZ4/trunk/libf/dyn3d/control.h	(revision 1340)
+++ 	(revision )
@@ -1,29 +1,0 @@
-!
-! $Header$
-!
-!
-!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
-!                 veillez à n'utiliser que des ! pour les commentaires
-!                 et à bien positionner les & des lignes de continuation 
-!                 (les placer en colonne 6 et en colonne 73)
-!
-!-----------------------------------------------------------------------
-! INCLUDE 'control.h'
-
-      COMMON/control/nday,day_step,                                     &
-     &              iperiod,iapp_tracvl,iconser,iecri,idissip,iphysiq , &
-     &              periodav,iecrimoy,dayref,anneeref,                  &
-     &              raz_date,offline,ip_ebil_dyn,config_inca,           &
-     &              planet_type,output_grads_dyn,ok_dynzon
-
-      INTEGER   nday,day_step,iperiod,iapp_tracvl,iconser,iecri,        &
-     &          idissip,iphysiq,iecrimoy,dayref,anneeref, raz_date      &
-     &          ,ip_ebil_dyn
-      REAL periodav
-      LOGICAL offline
-      CHARACTER (len=4) :: config_inca
-      CHARACTER(len=10) :: planet_type ! planet type ('earth','mars',...)
-      LOGICAL :: output_grads_dyn ! output dynamics diagnostics in
-                                  ! binary grads file 'dyn.dat' (y/n)
-      LOGICAL :: ok_dynzon
-!-----------------------------------------------------------------------
Index: LMDZ4/trunk/libf/dyn3d/control_mod.F90
===================================================================
--- LMDZ4/trunk/libf/dyn3d/control_mod.F90	(revision 1403)
+++ LMDZ4/trunk/libf/dyn3d/control_mod.F90	(revision 1403)
@@ -0,0 +1,27 @@
+!
+! $Id $
+!
+
+MODULE control_mod
+
+! LF 01/2010
+! Remplacement du fichier et common control
+
+  IMPLICIT NONE
+
+  REAL    :: periodav
+  INTEGER :: nday,day_step,iperiod,iapp_tracvl,nsplit_phys
+  INTEGER :: iconser,iecri,idissip,iphysiq,iecrimoy
+  INTEGER :: dayref,anneeref, raz_date, ip_ebil_dyn
+  LOGICAL :: offline
+  CHARACTER (len=4)  :: config_inca
+  CHARACTER (len=10) :: planet_type ! planet type ('earth','mars',...)
+  LOGICAL output_grads_dyn ! output dynamics diagnostics in
+                           ! binary grads file 'dyn.dat' (y/n)
+  LOGICAL ok_dynzon  ! output zonal transports in dynzon.nc file
+  LOGICAL ok_dyn_ins ! output instantaneous values of fields
+                     ! in the dynamics in NetCDF files dyn_hist*nc
+  LOGICAL ok_dyn_ave ! output averaged values of fields in the dynamics
+                     ! in NetCDF files dyn_hist*ave.nc
+
+END MODULE
Index: LMDZ4/trunk/libf/dyn3d/defrun.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/defrun.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/defrun.F	(revision 1403)
@@ -1,4 +1,4 @@
 !
-! $Header$
+! $Id$
 !
 c
@@ -6,4 +6,6 @@
       SUBROUTINE defrun( tapedef, etatinit, clesphy0 )
 c
+      USE control_mod
+ 
       IMPLICIT NONE
 c-----------------------------------------------------------------------
@@ -28,5 +30,4 @@
 #include "dimensions.h"
 #include "paramet.h"
-#include "control.h"
 #include "logic.h"
 #include "serre.h"
@@ -239,6 +240,6 @@
        clesphy0(i) = 0.
       ENDDO
-                          clesphy0(1) = FLOAT( iflag_con )
-                          clesphy0(2) = FLOAT( nbapp_rad )
+                          clesphy0(1) = REAL( iflag_con )
+                          clesphy0(2) = REAL( nbapp_rad )
 
        IF( cycle_diurne  ) clesphy0(3) =  1.
Index: LMDZ4/trunk/libf/dyn3d/disvert.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/disvert.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/disvert.F	(revision 1403)
@@ -111,5 +111,5 @@
       snorm  = 0.
       DO l = 1, llm
-         x = 2.*asin(1.) * (FLOAT(l)-0.5) / float(llm+1)
+         x = 2.*asin(1.) * (REAL(l)-0.5) / REAL(llm+1)
 
          IF (ok_strato) THEN
@@ -135,9 +135,9 @@
 
       DO l=1,llm
-        nivsigs(l) = FLOAT(l)
+        nivsigs(l) = REAL(l)
       ENDDO
 
       DO l=1,llmp1
-        nivsig(l)= FLOAT(l)
+        nivsig(l)= REAL(l)
       ENDDO
 
Index: LMDZ4/trunk/libf/dyn3d/dynetat0.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/dynetat0.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/dynetat0.F	(revision 1403)
@@ -1,4 +1,4 @@
 !
-! $Header$
+! $Id $
 !
       SUBROUTINE dynetat0(fichnom,vcov,ucov,
@@ -34,4 +34,5 @@
 #include "serre.h"
 #include "logic.h"
+#include "iniprint.h"
 
 c   Arguments:
@@ -58,6 +59,6 @@
       ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
       IF (ierr.NE.NF_NOERR) THEN
-        write(6,*)' Pb d''ouverture du fichier start.nc'
-        write(6,*)' ierr = ', ierr
+        write(lunout,*)'dynetat0: Pb d''ouverture du fichier start.nc'
+        write(lunout,*)' ierr = ', ierr
         CALL ABORT
       ENDIF
@@ -66,5 +67,5 @@
       ierr = NF_INQ_VARID (nid, "controle", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Le champ <controle> est absent"
+         write(lunout,*)"dynetat0: Le champ <controle> est absent"
          CALL abort
       ENDIF
@@ -75,5 +76,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Lecture echoue pour <controle>"
+         write(lunout,*)"dynetat0: Lecture echoue pour <controle>"
          CALL abort
       ENDIF
@@ -121,5 +122,6 @@
 c
 c
-      PRINT*,'rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa
+      write(lunout,*)'dynetat0: rad,omeg,g,cpp,kappa',
+     &               rad,omeg,g,cpp,kappa
 
       IF(   im.ne.iim           )  THEN
@@ -136,5 +138,5 @@
       ierr = NF_INQ_VARID (nid, "rlonu", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Le champ <rlonu> est absent"
+         write(lunout,*)"dynetat0: Le champ <rlonu> est absent"
          CALL abort
       ENDIF
@@ -145,5 +147,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Lecture echouee pour <rlonu>"
+         write(lunout,*)"dynetat0: Lecture echouee pour <rlonu>"
          CALL abort
       ENDIF
@@ -151,5 +153,5 @@
       ierr = NF_INQ_VARID (nid, "rlatu", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Le champ <rlatu> est absent"
+         write(lunout,*)"dynetat0: Le champ <rlatu> est absent"
          CALL abort
       ENDIF
@@ -160,5 +162,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Lecture echouee pour <rlatu>"
+         write(lunout,*)"dynetat0: Lecture echouee pour <rlatu>"
          CALL abort
       ENDIF
@@ -166,5 +168,5 @@
       ierr = NF_INQ_VARID (nid, "rlonv", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Le champ <rlonv> est absent"
+         write(lunout,*)"dynetat0: Le champ <rlonv> est absent"
          CALL abort
       ENDIF
@@ -175,5 +177,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Lecture echouee pour <rlonv>"
+         write(lunout,*)"dynetat0: Lecture echouee pour <rlonv>"
          CALL abort
       ENDIF
@@ -181,5 +183,5 @@
       ierr = NF_INQ_VARID (nid, "rlatv", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Le champ <rlatv> est absent"
+         write(lunout,*)"dynetat0: Le champ <rlatv> est absent"
          CALL abort
       ENDIF
@@ -190,5 +192,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Lecture echouee pour rlatv"
+         write(lunout,*)"dynetat0: Lecture echouee pour rlatv"
          CALL abort
       ENDIF
@@ -196,5 +198,5 @@
       ierr = NF_INQ_VARID (nid, "cu", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Le champ <cu> est absent"
+         write(lunout,*)"dynetat0: Le champ <cu> est absent"
          CALL abort
       ENDIF
@@ -205,5 +207,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Lecture echouee pour <cu>"
+         write(lunout,*)"dynetat0: Lecture echouee pour <cu>"
          CALL abort
       ENDIF
@@ -211,5 +213,5 @@
       ierr = NF_INQ_VARID (nid, "cv", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Le champ <cv> est absent"
+         write(lunout,*)"dynetat0: Le champ <cv> est absent"
          CALL abort
       ENDIF
@@ -220,5 +222,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Lecture echouee pour <cv>"
+         write(lunout,*)"dynetat0: Lecture echouee pour <cv>"
          CALL abort
       ENDIF
@@ -226,5 +228,5 @@
       ierr = NF_INQ_VARID (nid, "aire", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Le champ <aire> est absent"
+         write(lunout,*)"dynetat0: Le champ <aire> est absent"
          CALL abort
       ENDIF
@@ -235,5 +237,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Lecture echouee pour <aire>"
+         write(lunout,*)"dynetat0: Lecture echouee pour <aire>"
          CALL abort
       ENDIF
@@ -241,5 +243,5 @@
       ierr = NF_INQ_VARID (nid, "phisinit", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Le champ <phisinit> est absent"
+         write(lunout,*)"dynetat0: Le champ <phisinit> est absent"
          CALL abort
       ENDIF
@@ -250,5 +252,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Lecture echouee pour <phisinit>"
+         write(lunout,*)"dynetat0: Lecture echouee pour <phisinit>"
          CALL abort
       ENDIF
@@ -256,5 +258,5 @@
       ierr = NF_INQ_VARID (nid, "temps", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Le champ <temps> est absent"
+         write(lunout,*)"dynetat0: Le champ <temps> est absent"
          CALL abort
       ENDIF
@@ -265,5 +267,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Lecture echouee <temps>"
+         write(lunout,*)"dynetat0: Lecture echouee <temps>"
          CALL abort
       ENDIF
@@ -271,5 +273,5 @@
       ierr = NF_INQ_VARID (nid, "ucov", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Le champ <ucov> est absent"
+         write(lunout,*)"dynetat0: Le champ <ucov> est absent"
          CALL abort
       ENDIF
@@ -280,5 +282,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Lecture echouee pour <ucov>"
+         write(lunout,*)"dynetat0: Lecture echouee pour <ucov>"
          CALL abort
       ENDIF
@@ -286,5 +288,5 @@
       ierr = NF_INQ_VARID (nid, "vcov", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Le champ <vcov> est absent"
+         write(lunout,*)"dynetat0: Le champ <vcov> est absent"
          CALL abort
       ENDIF
@@ -295,5 +297,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Lecture echouee pour <vcov>"
+         write(lunout,*)"dynetat0: Lecture echouee pour <vcov>"
          CALL abort
       ENDIF
@@ -301,5 +303,5 @@
       ierr = NF_INQ_VARID (nid, "teta", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Le champ <teta> est absent"
+         write(lunout,*)"dynetat0: Le champ <teta> est absent"
          CALL abort
       ENDIF
@@ -310,5 +312,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Lecture echouee pour <teta>"
+         write(lunout,*)"dynetat0: Lecture echouee pour <teta>"
          CALL abort
       ENDIF
@@ -319,6 +321,7 @@
         ierr =  NF_INQ_VARID (nid, tname(iq), nvarid)
         IF (ierr .NE. NF_NOERR) THEN
-           PRINT*, "dynetat0: Le champ <"//tname(iq)//"> est absent"
-           PRINT*, "          Il est donc initialise a zero"
+           write(lunout,*)"dynetat0: Le champ <"//tname(iq)//
+     &                    "> est absent"
+           write(lunout,*)"          Il est donc initialise a zero"
            q(:,:,iq)=0.
         ELSE
@@ -329,6 +332,6 @@
 #endif
           IF (ierr .NE. NF_NOERR) THEN
-             PRINT*, "dynetat0: Lecture echouee pour "//tname(iq)
-             CALL abort
+            write(lunout,*)"dynetat0: Lecture echouee pour "//tname(iq)
+            CALL abort
           ENDIF
         ENDIF
@@ -338,5 +341,5 @@
       ierr = NF_INQ_VARID (nid, "masse", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Le champ <masse> est absent"
+         write(lunout,*)"dynetat0: Le champ <masse> est absent"
          CALL abort
       ENDIF
@@ -347,5 +350,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Lecture echouee pour <masse>"
+         write(lunout,*)"dynetat0: Lecture echouee pour <masse>"
          CALL abort
       ENDIF
@@ -353,5 +356,5 @@
       ierr = NF_INQ_VARID (nid, "ps", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Le champ <ps> est absent"
+         write(lunout,*)"dynetat0: Le champ <ps> est absent"
          CALL abort
       ENDIF
@@ -362,5 +365,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "dynetat0: Lecture echouee pour <ps>"
+         write(lunout,*)"dynetat0: Lecture echouee pour <ps>"
          CALL abort
       ENDIF
Index: LMDZ4/trunk/libf/dyn3d/dynredem.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/dynredem.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/dynredem.F	(revision 1403)
@@ -8,4 +8,5 @@
 #endif
       USE infotrac
+ 
       IMPLICIT NONE
 c=======================================================================
@@ -25,4 +26,5 @@
 #include "description.h"
 #include "serre.h"
+#include "iniprint.h"
 
 c   Arguments:
@@ -72,9 +74,9 @@
        tab_cntrl(l) = 0.
       ENDDO
-       tab_cntrl(1)  = FLOAT(iim)
-       tab_cntrl(2)  = FLOAT(jjm)
-       tab_cntrl(3)  = FLOAT(llm)
-       tab_cntrl(4)  = FLOAT(day_ref)
-       tab_cntrl(5)  = FLOAT(annee_ref)
+       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
@@ -116,6 +118,6 @@
       ENDIF
 
-       tab_cntrl(30) = FLOAT(iday_end)
-       tab_cntrl(31) = FLOAT(itau_dyn + itaufin)
+       tab_cntrl(30) = REAL(iday_end)
+       tab_cntrl(31) = REAL(itau_dyn + itaufin)
 c
 c    .........................................................
@@ -125,6 +127,7 @@
       ierr = NF_CREATE(fichnom, NF_CLOBBER, nid)
       IF (ierr.NE.NF_NOERR) THEN
-         WRITE(6,*)" Pb d ouverture du fichier "//fichnom
-         WRITE(6,*)' ierr = ', ierr
+         write(lunout,*)"dynredem0: Pb d ouverture du fichier "
+     &                  //trim(fichnom)
+         write(lunout,*)' ierr = ', ierr
          CALL ABORT
       ENDIF
@@ -508,7 +511,8 @@
       ierr = NF_CLOSE(nid) ! fermer le fichier
 
-      PRINT*,'iim,jjm,llm,iday_end',iim,jjm,llm,iday_end
-      PRINT*,'rad,omeg,g,cpp,kappa',
-     ,        rad,omeg,g,cpp,kappa
+      write(lunout,*)'dynredem0: iim,jjm,llm,iday_end',
+     &               iim,jjm,llm,iday_end
+      write(lunout,*)'dynredem0: rad,omeg,g,cpp,kappa',
+     &        rad,omeg,g,cpp,kappa
 
       RETURN
@@ -517,4 +521,6 @@
      .                     vcov,ucov,teta,q,masse,ps)
       USE infotrac
+      USE control_mod
+ 
       IMPLICIT NONE
 c=================================================================
@@ -528,5 +534,6 @@
 #include "comgeom.h"
 #include "temps.h"
-#include "control.h"
+#include "iniprint.h"
+
 
       INTEGER l
@@ -555,5 +562,5 @@
       ierr = NF_OPEN(fichnom, NF_WRITE, nid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "Pb. d ouverture "//fichnom
+         write(lunout,*)"dynredem1: Pb. d ouverture "//trim(fichnom)
          CALL abort
       ENDIF
@@ -564,5 +571,5 @@
       ierr = NF_INQ_VARID(nid, "temps", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         print *, NF_STRERROR(ierr)
+         write(lunout,*) NF_STRERROR(ierr)
          abort_message='Variable temps n est pas definie'
          CALL abort_gcm(modname,abort_message,ierr)
@@ -573,5 +580,5 @@
       ierr = NF_PUT_VAR1_REAL (nid,nvarid,nb,time)
 #endif
-      PRINT*, "Enregistrement pour ", nb, time
+      write(lunout,*) "dynredem1: Enregistrement pour ", nb, time
 
 c
@@ -589,5 +596,5 @@
       ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
 #endif
-       tab_cntrl(31) = FLOAT(itau_dyn + itaufin)
+       tab_cntrl(31) = REAL(itau_dyn + itaufin)
 #ifdef NC_DOUBLE
       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tab_cntrl)
@@ -600,6 +607,7 @@
       ierr = NF_INQ_VARID(nid, "ucov", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "Variable ucov n est pas definie"
-         CALL abort
+         abort_message="Variable ucov n est pas definie"
+         ierr=1
+         CALL abort_gcm(modname,abort_message,ierr)
       ENDIF
 #ifdef NC_DOUBLE
@@ -611,6 +619,7 @@
       ierr = NF_INQ_VARID(nid, "vcov", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "Variable vcov n est pas definie"
-         CALL abort
+         abort_message="Variable vcov n est pas definie"
+         ierr=1
+         CALL abort_gcm(modname,abort_message,ierr)
       ENDIF
 #ifdef NC_DOUBLE
@@ -622,6 +631,7 @@
       ierr = NF_INQ_VARID(nid, "teta", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "Variable teta n est pas definie"
-         CALL abort
+         abort_message="Variable teta n est pas definie"
+         ierr=1
+         CALL abort_gcm(modname,abort_message,ierr)
       ENDIF
 #ifdef NC_DOUBLE
@@ -635,6 +645,7 @@
          ierr_file = NF_OPEN ("start_trac.nc", NF_NOWRITE,nid_trac)
          IF (ierr_file .NE.NF_NOERR) THEN
-            write(6,*)' Pb d''ouverture du fichier start_trac.nc'
-            write(6,*)' ierr = ', ierr_file 
+            write(lunout,*)'dynredem1: Pb d''ouverture du fichier',
+     &                     ' start_trac.nc'
+            write(lunout,*)' ierr = ', ierr_file 
          ENDIF
       END IF
@@ -646,6 +657,7 @@
             ierr = NF_INQ_VARID(nid, tname(iq), nvarid)
             IF (ierr .NE. NF_NOERR) THEN
-               PRINT*, "Variable  tname(iq) n est pas definie"
-               CALL abort
+               abort_message="Variable  tname(iq) n est pas definie"
+               ierr=1
+               CALL abort_gcm(modname,abort_message,ierr)
             ENDIF
 #ifdef NC_DOUBLE
@@ -659,9 +671,12 @@
              ierr = NF_INQ_VARID (nid_trac, tname(iq), nvarid_trac)
              IF (ierr .NE. NF_NOERR) THEN
-                PRINT*, tname(iq),"est absent de start_trac.nc"
+                write(lunout,*) "dynredem1: ",trim(tname(iq)),
+     &                          " est absent de start_trac.nc"
                 ierr = NF_INQ_VARID(nid, tname(iq), nvarid)
                 IF (ierr .NE. NF_NOERR) THEN
-                   PRINT*, "Variable ", tname(iq)," n est pas definie"
-                   CALL abort
+                   abort_message="dynredem1: Variable "//
+     &                     trim(tname(iq))//" n est pas definie"
+                   ierr=1
+                   CALL abort_gcm(modname,abort_message,ierr)
                 ENDIF
 #ifdef NC_DOUBLE
@@ -672,5 +687,6 @@
                 
              ELSE
-                PRINT*, tname(iq), "est present dans start_trac.nc"
+                write(lunout,*) "dynredem1: ",trim(tname(iq)),
+     &              " est present dans start_trac.nc"
 #ifdef NC_DOUBLE
                ierr = NF_GET_VAR_DOUBLE(nid_trac, nvarid_trac, trac_tmp)
@@ -679,11 +695,15 @@
 #endif
                 IF (ierr .NE. NF_NOERR) THEN
-                   PRINT*, "Lecture echouee pour", tname(iq)
-                   CALL abort
+                   abort_message="dynredem1: Lecture echouee pour"//
+     &                    trim(tname(iq))
+                   ierr=1
+                   CALL abort_gcm(modname,abort_message,ierr)
                 ENDIF
                 ierr = NF_INQ_VARID(nid, tname(iq), nvarid)
                 IF (ierr .NE. NF_NOERR) THEN
-                   PRINT*, "Variable ", tname(iq)," n est pas definie"
-                   CALL abort
+                   abort_message="dynredem1: Variable "//
+     &                trim(tname(iq))//" n est pas definie"
+                   ierr=1
+                   CALL abort_gcm(modname,abort_message,ierr)
                 ENDIF
 #ifdef NC_DOUBLE
@@ -699,6 +719,8 @@
              ierr = NF_INQ_VARID(nid, tname(iq), nvarid)
              IF (ierr .NE. NF_NOERR) THEN
-                PRINT*, "Variable  tname(iq) n est pas definie"
-                CALL abort
+                abort_message="dynredem1: Variable "//
+     &                trim(tname(iq))//" n est pas definie"
+                   ierr=1
+                   CALL abort_gcm(modname,abort_message,ierr)
              ENDIF
 #ifdef NC_DOUBLE
@@ -715,6 +737,7 @@
       ierr = NF_INQ_VARID(nid, "masse", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "Variable masse n est pas definie"
-         CALL abort
+         abort_message="dynredem1: Variable masse n est pas definie"
+         ierr=1
+         CALL abort_gcm(modname,abort_message,ierr)
       ENDIF
 #ifdef NC_DOUBLE
@@ -726,6 +749,7 @@
       ierr = NF_INQ_VARID(nid, "ps", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "Variable ps n est pas definie"
-         CALL abort
+         abort_message="dynredem1: Variable ps n est pas definie"
+         ierr=1
+         CALL abort_gcm(modname,abort_message,ierr)
       ENDIF
 #ifdef NC_DOUBLE
Index: LMDZ4/trunk/libf/dyn3d/etat0_netcdf.F90
===================================================================
--- LMDZ4/trunk/libf/dyn3d/etat0_netcdf.F90	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/etat0_netcdf.F90	(revision 1403)
@@ -11,4 +11,5 @@
 ! Note: This routine is designed to work for Earth
 !-------------------------------------------------------------------------------
+  USE control_mod
 #ifdef CPP_EARTH
   USE startvar
@@ -72,5 +73,4 @@
 
 #include "comdissnew.h"
-#include "control.h"
 #include "serre.h"
 #include "clesphys.h"
@@ -103,4 +103,5 @@
   REAL    :: tau_thermals, solarlong0,  seuil_inversion
   INTEGER :: read_climoz ! read ozone climatology
+  REAL    :: alp_offset
 !  Allowed values are 0, 1 and 2
 !     0: do not read an ozone climatology
@@ -132,5 +133,6 @@
                    iflag_thermals,nsplit_thermals,tau_thermals,         &
                    iflag_thermals_ed,iflag_thermals_optflux,            &
-                   iflag_coupl,iflag_clos,iflag_wake, read_climoz )
+                   iflag_coupl,iflag_clos,iflag_wake, read_climoz,      &
+                   alp_offset )
 
 ! co2_ppm0 : initial value of atmospheric CO2 from .def file (co2_ppm value)
Index: LMDZ4/trunk/libf/dyn3d/exner_hyb.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/exner_hyb.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/exner_hyb.F	(revision 1403)
@@ -1,4 +1,4 @@
 !
-! $Header$
+! $Id $
 !
       SUBROUTINE  exner_hyb ( ngrid, ps, p,alpha,beta, pks, pk, pkf )
@@ -51,5 +51,33 @@
       REAL SSUM
 c
+
+      if (llm.eq.1) then
+        ! Specific behaviour for Shallow Water (1 vertical layer) case
       
+        ! Sanity checks
+        if (kappa.ne.1) then
+          call abort_gcm("exner_hyb",
+     &    "kappa!=1 , but running in Shallow Water mode!!",42)
+        endif
+        if (cpp.ne.r) then
+        call abort_gcm("exner_hyb",
+     &    "cpp!=r , but running in Shallow Water mode!!",42)
+        endif
+        
+        ! Compute pks(:),pk(:),pkf(:)
+        
+        DO   ij  = 1, ngrid
+          pks(ij) = (cpp/preff) * ps(ij) 
+          pk(ij,1) = .5*pks(ij)
+        ENDDO
+        
+        CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
+        CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 
+        
+        ! our work is done, exit routine
+        return
+      endif ! of if (llm.eq.1)
+
+     
       unpl2k    = 1.+ 2.* kappa
 c
Index: LMDZ4/trunk/libf/dyn3d/extrapol.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/extrapol.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/extrapol.F	(revision 1403)
@@ -1,4 +1,4 @@
 !
-! $Header$
+! $Id$
 !
 C
@@ -158,5 +158,5 @@
                jlat = jy(k)
                pwork(i,j) = pwork(i,j)
-     $                      + pfild(ilon,jlat) * zmask(k)/FLOAT(inbor)
+     $                      + pfild(ilon,jlat) * zmask(k)/REAL(inbor)
             ENDDO
          ENDIF
Index: LMDZ4/trunk/libf/dyn3d/fluxstokenc.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/fluxstokenc.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/fluxstokenc.F	(revision 1403)
@@ -4,6 +4,6 @@
       SUBROUTINE fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
      . time_step,itau )
-#ifdef CPP_EARTH
-! This routine is designed to work for Earth and with ioipsl
+#ifdef CPP_IOIPSL
+! This routine is designed to work with ioipsl
 
        USE IOIPSL
@@ -114,10 +114,10 @@
       DO l=1,llm
          DO ij = 1,ip1jmp1
-            pbaruc(ij,l) = pbaruc(ij,l)/float(istdyn)
-            tetac(ij,l) = tetac(ij,l)/float(istdyn)
-            phic(ij,l) = phic(ij,l)/float(istdyn)
+            pbaruc(ij,l) = pbaruc(ij,l)/REAL(istdyn)
+            tetac(ij,l) = tetac(ij,l)/REAL(istdyn)
+            phic(ij,l) = phic(ij,l)/REAL(istdyn)
          ENDDO
          DO ij = 1,ip1jm
-            pbarvc(ij,l) = pbarvc(ij,l)/float(istdyn)
+            pbarvc(ij,l) = pbarvc(ij,l)/REAL(istdyn)
          ENDDO
       ENDDO
@@ -141,5 +141,5 @@
 
          iadvtr=0
-	Print*,'ITAU auqel on stoke les fluxmasses',itau
+	write(lunout,*)'ITAU auquel on stoke les fluxmasses',itau
 	
 	call histwrite(fluxid, 'masse', itau, massem,
@@ -167,7 +167,7 @@
 #else
       write(lunout,*)
-     & 'fluxstokenc: Needs Earth physics (and ioipsl) to function'
+     & 'fluxstokenc: Needs IOIPSL to function'
 #endif
-! of #ifdef CPP_EARTH
+! of #ifdef CPP_IOIPSL
       RETURN
       END
Index: LMDZ4/trunk/libf/dyn3d/friction.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/friction.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/friction.F	(revision 1403)
@@ -1,7 +1,10 @@
 !
-! $Header$
+! $Id$
 !
 c=======================================================================
       SUBROUTINE friction(ucov,vcov,pdt)
+
+      USE control_mod
+ 
       IMPLICIT NONE
 
@@ -21,5 +24,4 @@
 #include "paramet.h"
 #include "comgeom2.h"
-#include "control.h"
 #include "comconst.h"
 
Index: LMDZ4/trunk/libf/dyn3d/fxhyp.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/fxhyp.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/fxhyp.F	(revision 1403)
@@ -1,4 +1,4 @@
 !
-! $Header$
+! $Id$
 !
 c
@@ -89,5 +89,5 @@
 
        DO i = 0, nmax2 
-        xtild(i) = - pi + FLOAT(i) * depi /nmax2
+        xtild(i) = - pi + REAL(i) * depi /nmax2
        ENDDO
 
@@ -235,5 +235,5 @@
       DO 1500 i = ii1, ii2
 
-      xlon2 = - pi + (FLOAT(i) + xuv - decalx) * depi / FLOAT(iim) 
+      xlon2 = - pi + (REAL(i) + xuv - decalx) * depi / REAL(iim) 
 
       Xfi    = xlon2
@@ -280,5 +280,5 @@
 550   CONTINUE
 
-       xxprim(i) = depi/ ( FLOAT(iim) * Xprimin )
+       xxprim(i) = depi/ ( REAL(iim) * Xprimin )
        xvrai(i)  =  xi + xzoom
 
Index: LMDZ4/trunk/libf/dyn3d/fxy.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/fxy.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/fxy.F	(revision 1403)
@@ -1,4 +1,4 @@
 !
-! $Header$
+! $Id$
 !
       SUBROUTINE fxy (rlatu,yprimu,rlatv,yprimv,rlatu1,yprimu1,
@@ -32,6 +32,6 @@
 c
        DO j = 1, jjm + 1 
-          rlatu(j) = fy    ( FLOAT( j )        )
-         yprimu(j) = fyprim( FLOAT( j )        )
+          rlatu(j) = fy    ( REAL( j )        )
+         yprimu(j) = fyprim( REAL( j )        )
        ENDDO
 
@@ -39,11 +39,11 @@
        DO j = 1, jjm
 
-         rlatv(j)  = fy    ( FLOAT( j ) + 0.5  )
-         rlatu1(j) = fy    ( FLOAT( j ) + 0.25 ) 
-         rlatu2(j) = fy    ( FLOAT( j ) + 0.75 ) 
+         rlatv(j)  = fy    ( REAL( j ) + 0.5  )
+         rlatu1(j) = fy    ( REAL( j ) + 0.25 ) 
+         rlatu2(j) = fy    ( REAL( j ) + 0.75 ) 
 
-        yprimv(j)  = fyprim( FLOAT( j ) + 0.5  ) 
-        yprimu1(j) = fyprim( FLOAT( j ) + 0.25 )
-        yprimu2(j) = fyprim( FLOAT( j ) + 0.75 )
+        yprimv(j)  = fyprim( REAL( j ) + 0.5  ) 
+        yprimu1(j) = fyprim( REAL( j ) + 0.25 )
+        yprimu2(j) = fyprim( REAL( j ) + 0.75 )
 
        ENDDO
@@ -53,13 +53,13 @@
 c
        DO i = 1, iim + 1
-           rlonv(i)     = fx    (   FLOAT( i )          )
-           rlonu(i)     = fx    (   FLOAT( i ) + 0.5    )
-        rlonm025(i)     = fx    (   FLOAT( i ) - 0.25  )
-        rlonp025(i)     = fx    (   FLOAT( i ) + 0.25  )
+           rlonv(i)     = fx    (   REAL( i )          )
+           rlonu(i)     = fx    (   REAL( i ) + 0.5    )
+        rlonm025(i)     = fx    (   REAL( i ) - 0.25  )
+        rlonp025(i)     = fx    (   REAL( i ) + 0.25  )
 
-         xprimv  (i)    = fxprim (  FLOAT( i )          )
-         xprimu  (i)    = fxprim (  FLOAT( i ) + 0.5    )
-        xprimm025(i)    = fxprim (  FLOAT( i ) - 0.25   )
-        xprimp025(i)    = fxprim (  FLOAT( i ) + 0.25   )
+         xprimv  (i)    = fxprim (  REAL( i )          )
+         xprimu  (i)    = fxprim (  REAL( i ) + 0.5    )
+        xprimm025(i)    = fxprim (  REAL( i ) - 0.25   )
+        xprimp025(i)    = fxprim (  REAL( i ) + 0.25   )
        ENDDO
 
Index: LMDZ4/trunk/libf/dyn3d/fxysinus.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/fxysinus.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/fxysinus.F	(revision 1403)
@@ -1,4 +1,4 @@
 !
-! $Header$
+! $Id$
 !
       SUBROUTINE fxysinus (rlatu,yprimu,rlatv,yprimv,rlatu1,yprimu1,
@@ -32,6 +32,6 @@
 c
        DO j = 1, jjm + 1 
-          rlatu(j) = fy    ( FLOAT( j )        )
-         yprimu(j) = fyprim( FLOAT( j )        )
+          rlatu(j) = fy    ( REAL( j )        )
+         yprimu(j) = fyprim( REAL( j )        )
        ENDDO
 
@@ -39,11 +39,11 @@
        DO j = 1, jjm
 
-         rlatv(j)  = fy    ( FLOAT( j ) + 0.5  )
-         rlatu1(j) = fy    ( FLOAT( j ) + 0.25 ) 
-         rlatu2(j) = fy    ( FLOAT( j ) + 0.75 ) 
+         rlatv(j)  = fy    ( REAL( j ) + 0.5  )
+         rlatu1(j) = fy    ( REAL( j ) + 0.25 ) 
+         rlatu2(j) = fy    ( REAL( j ) + 0.75 ) 
 
-        yprimv(j)  = fyprim( FLOAT( j ) + 0.5  ) 
-        yprimu1(j) = fyprim( FLOAT( j ) + 0.25 )
-        yprimu2(j) = fyprim( FLOAT( j ) + 0.75 )
+        yprimv(j)  = fyprim( REAL( j ) + 0.5  ) 
+        yprimu1(j) = fyprim( REAL( j ) + 0.25 )
+        yprimu2(j) = fyprim( REAL( j ) + 0.75 )
 
        ENDDO
@@ -53,13 +53,13 @@
 c
        DO i = 1, iim + 1
-           rlonv(i)     = fx    (   FLOAT( i )          )
-           rlonu(i)     = fx    (   FLOAT( i ) + 0.5    )
-        rlonm025(i)     = fx    (   FLOAT( i ) - 0.25  )
-        rlonp025(i)     = fx    (   FLOAT( i ) + 0.25  )
+           rlonv(i)     = fx    (   REAL( i )          )
+           rlonu(i)     = fx    (   REAL( i ) + 0.5    )
+        rlonm025(i)     = fx    (   REAL( i ) - 0.25  )
+        rlonp025(i)     = fx    (   REAL( i ) + 0.25  )
 
-         xprimv  (i)    = fxprim (  FLOAT( i )          )
-         xprimu  (i)    = fxprim (  FLOAT( i ) + 0.5    )
-        xprimm025(i)    = fxprim (  FLOAT( i ) - 0.25   )
-        xprimp025(i)    = fxprim (  FLOAT( i ) + 0.25   )
+         xprimv  (i)    = fxprim (  REAL( i )          )
+         xprimu  (i)    = fxprim (  REAL( i ) + 0.5    )
+        xprimm025(i)    = fxprim (  REAL( i ) - 0.25   )
+        xprimp025(i)    = fxprim (  REAL( i ) + 0.25   )
        ENDDO
 
Index: LMDZ4/trunk/libf/dyn3d/fyhyp.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/fyhyp.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/fyhyp.F	(revision 1403)
@@ -1,4 +1,4 @@
 !
-! $Header$
+! $Id$
 !
 c
@@ -75,5 +75,5 @@
        depi     = 2. * pi
        pis2     = pi/2.
-       pisjm    = pi/ FLOAT(jjm)
+       pisjm    = pi/ REAL(jjm)
        epsilon  = 1.e-3
        y0       =  yzoomdeg * pi/180. 
@@ -94,5 +94,5 @@
 
        DO i = 0, nmax2 
-        yt(i) = - pis2  + FLOAT(i)* pi /nmax2
+        yt(i) = - pis2  + REAL(i)* pi /nmax2
        ENDDO
 
@@ -210,5 +210,5 @@
        DO 1500 j =  1,jlat
         yo1   = 0.
-        ylon2 =  - pis2 + pisjm * ( FLOAT(j)  + yuv  -1.)  
+        ylon2 =  - pis2 + pisjm * ( REAL(j)  + yuv  -1.)  
         yfi    = ylon2
 c
Index: LMDZ4/trunk/libf/dyn3d/gcm.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/gcm.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/gcm.F	(revision 1403)
@@ -15,4 +15,5 @@
       USE filtreg_mod
       USE infotrac
+      USE control_mod
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -68,13 +69,15 @@
 #include "logic.h"
 #include "temps.h"
-#include "control.h"
+!!!!!!!!!!!#include "control.h"
 #include "ener.h"
 #include "description.h"
 #include "serre.h"
-#include "com_io_dyn.h"
+!#include "com_io_dyn.h"
 #include "iniprint.h"
 #include "tracstoke.h"
+#ifdef INCA
+! Only INCA needs these informations (from the Earth's physics)
 #include "indicesol.h"
-
+#endif
       INTEGER         longcles
       PARAMETER     ( longcles = 20 )
@@ -181,5 +184,5 @@
       if (planet_type.eq."earth") then
 #ifdef CPP_EARTH
-      CALL Init_Phys_lmdz(iim,jjp1,llm,1,(jjm-1)*iim+2)
+      CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
       call InitComgeomphy
 #endif
@@ -241,16 +244,24 @@
       if (read_start) then
       ! we still need to run iniacademic to initialize some
-      ! constants & fields, if we run the 'newtonian' case:
-        if (iflag_phys.eq.2) then
+      ! constants & fields, if we run the 'newtonian' or 'SW' cases:
+        if (iflag_phys.ne.1) then
           CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
         endif
-!#ifdef CPP_IOIPSL
+
         if (planet_type.eq."earth") then
 #ifdef CPP_EARTH
 ! Load an Earth-format start file
          CALL dynetat0("start.nc",vcov,ucov,
-     .              teta,q,masse,ps,phis, time_0)
+     &              teta,q,masse,ps,phis, time_0)
+#else
+        ! SW model also has Earth-format start files
+        ! (but can be used without the CPP_EARTH directive)
+          if (iflag_phys.eq.0) then
+            CALL dynetat0("start.nc",vcov,ucov,
+     &              teta,q,masse,ps,phis, time_0)
+          endif
 #endif
         endif ! of if (planet_type.eq."earth")
+        
 c       write(73,*) 'ucov',ucov
 c       write(74,*) 'vcov',vcov
@@ -294,5 +305,5 @@
       ENDIF
 
-      zdtvr    = daysec/FLOAT(day_step)
+      zdtvr    = daysec/REAL(day_step)
         IF(dtvr.NE.zdtvr) THEN
          WRITE(lunout,*)
@@ -303,5 +314,14 @@
 C on remet le calendrier à zero si demande
 c
-      if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
+      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'
+      ELSE IF (annee_ref .ne. anneeref .or. day_ref .ne. dayref) THEN
         write(lunout,*)
      .  'GCM: Attention les dates initiales lues dans le fichier'
@@ -309,22 +329,33 @@
      .  ' 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
-        if (raz_date .ne. 1) then
-          write(lunout,*)
-     .    'GCM: On garde les dates du fichier restart'
-        else
-          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'
-        endif
-      ELSE
-        raz_date = 0
-      endif
+        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
+        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
+        write(lunout,*)' Pas de remise a zero'
+      ENDIF
+
+c      if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
+c        write(lunout,*)
+c     .  'GCM: Attention les dates initiales lues dans le fichier'
+c        write(lunout,*)
+c     .  ' restart ne correspondent pas a celles lues dans '
+c        write(lunout,*)' gcm.def'
+c        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
+c        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
+c        if (raz_date .ne. 1) then
+c          write(lunout,*)
+c     .    'GCM: On garde les dates du fichier restart'
+c        else
+c          annee_ref = anneeref
+c          day_ref = dayref
+c          day_ini = dayref
+c          itau_dyn = 0
+c          itau_phy = 0
+c          time_0 = 0.
+c          write(lunout,*)
+c     .   'GCM: On reinitialise a la date lue dans gcm.def'
+c        endif
+c      ELSE
+c        raz_date = 0
+c      endif
 
 #ifdef CPP_IOIPSL
@@ -355,19 +386,23 @@
       nbetatmoy = nday / periodav + 1
 
+      if (iflag_phys.eq.1) then
+      ! these initialisations have already been done (via iniacademic)
+      ! if running in SW or Newtonian mode
 c-----------------------------------------------------------------------
 c   Initialisation des constantes dynamiques :
 c   ------------------------------------------
-      dtvr = zdtvr
-      CALL iniconst
+        dtvr = zdtvr
+        CALL iniconst
 
 c-----------------------------------------------------------------------
 c   Initialisation de la geometrie :
 c   --------------------------------
-      CALL inigeom
+        CALL inigeom
 
 c-----------------------------------------------------------------------
 c   Initialisation du filtre :
 c   --------------------------
-      CALL inifilr
+        CALL inifilr
+      endif ! of if (iflag_phys.eq.1)
 c
 c-----------------------------------------------------------------------
@@ -405,5 +440,5 @@
          if (planet_type.eq."earth") then
 #ifdef CPP_EARTH
-         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys ,
+         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys ,
      ,                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp     )
 #endif
@@ -440,21 +475,22 @@
 
 #ifdef CPP_IOIPSL
-      if ( 1.eq.1) then
       time_step = zdtvr
-      t_ops = iecri * daysec
-      t_wrt = iecri * daysec
-!      CALL inithist(dynhist_file,day_ref,annee_ref,time_step,
-!    .              t_ops, t_wrt, histid, histvid)
-
-!     IF (ok_dynzon) THEN 
-!        t_ops = iperiod * time_step
-!        t_wrt = periodav * daysec
-!        CALL initdynav(dynhistave_file,day_ref,annee_ref,time_step,
-!    .        t_ops, t_wrt, histaveid)
-!     END IF
+      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)
+      END IF
       dtav = iperiod*dtvr/daysec
-      endif
-
-
 #endif
 ! #endif of #ifdef CPP_IOIPSL
Index: LMDZ4/trunk/libf/dyn3d/grid_atob.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/grid_atob.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/grid_atob.F	(revision 1403)
@@ -1,4 +1,4 @@
 !
-! $Header$
+! $Id$
 !
       SUBROUTINE grille_m(imdep, jmdep, xdata, ydata, entree,
@@ -717,14 +717,14 @@
 c Calculs intermediares:
 c
-      xtmp(1) = -180.0 + 360.0/FLOAT(imtmp) / 2.0
+      xtmp(1) = -180.0 + 360.0/REAL(imtmp) / 2.0
       DO i = 2, imtmp
-         xtmp(i) = xtmp(i-1) + 360.0/FLOAT(imtmp)
+         xtmp(i) = xtmp(i-1) + 360.0/REAL(imtmp)
       ENDDO
       DO i = 1, imtmp
          xtmp(i) = xtmp(i) /180.0 * 4.0*ATAN(1.0)
       ENDDO
-      ytmp(1) = -90.0 + 180.0/FLOAT(jmtmp) / 2.0
+      ytmp(1) = -90.0 + 180.0/REAL(jmtmp) / 2.0
       DO j = 2, jmtmp
-         ytmp(j) = ytmp(j-1) + 180.0/FLOAT(jmtmp)
+         ytmp(j) = ytmp(j-1) + 180.0/REAL(jmtmp)
       ENDDO
       DO j = 1, jmtmp
Index: LMDZ4/trunk/libf/dyn3d/grid_noro.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/grid_noro.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/grid_noro.F	(revision 1403)
@@ -1,4 +1,4 @@
 !
-! $Header$
+! $Id$
 !
 c
@@ -93,5 +93,5 @@
       xpi=acos(-1.)
       rad    = 6 371 229.
-      zdeltay=2.*xpi/float(jusn)*rad
+      zdeltay=2.*xpi/REAL(jusn)*rad
 c
 c utilise-t'on un masque lu?
@@ -215,6 +215,6 @@
 c  SUMMATION OVER GRIDPOINT AREA
 c 
-      zleny=xpi/float(jusn)*rad
-      xincr=xpi/2./float(jusn)
+      zleny=xpi/REAL(jusn)*rad
+      xincr=xpi/2./REAL(jusn)
        DO ii = 1, imar+1
        DO jj = 1, jmar
@@ -468,5 +468,5 @@
       DO IS=-1,1
         DO JS=-1,1
-          WEIGHTpb(IS,JS)=1./FLOAT((1+IS**2)*(1+JS**2))
+          WEIGHTpb(IS,JS)=1./REAL((1+IS**2)*(1+JS**2))
           SUM=SUM+WEIGHTpb(IS,JS)
         ENDDO
Index: LMDZ4/trunk/libf/dyn3d/grilles_gcm_netcdf.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/grilles_gcm_netcdf.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/grilles_gcm_netcdf.F	(revision 1403)
@@ -1,4 +1,4 @@
 !
-! $Header$
+! $Id$
 !
 c
@@ -218,13 +218,13 @@
       open (20,file='grille.dat',form='unformatted',access='direct'
      s      ,recl=4*ip1jmp1)
-      write(20,rec=1) ((float(mod(i,2)+mod(j,2)),i=1,iip1),j=1,jjp1)
-      write(20,rec=2) ((float(mod(i,2)*mod(j,2)),i=1,iip1),j=1,jjp1)
+      write(20,rec=1) ((REAL(mod(i,2)+mod(j,2)),i=1,iip1),j=1,jjp1)
+      write(20,rec=2) ((REAL(mod(i,2)*mod(j,2)),i=1,iip1),j=1,jjp1)
       do j=2,jjm
          dlat1(j)=180.*(rlatv(j)-rlatv(j-1))/pi
-c        dlat2(j)=180.*fyprim(float(j))/pi
+c        dlat2(j)=180.*fyprim(REAL(j))/pi
       enddo
       do i=2,iip1
          dlon1(i)=180.*(rlonu(i)-rlonu(i-1))/pi
-c        dlon2(i)=180.*fxprim(float(i))/pi
+c        dlon2(i)=180.*fxprim(REAL(i))/pi
       enddo
       do j=2,jjm
Index: LMDZ4/trunk/libf/dyn3d/guide_mod.F90
===================================================================
--- LMDZ4/trunk/libf/dyn3d/guide_mod.F90	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/guide_mod.F90	(revision 1403)
@@ -62,4 +62,6 @@
   SUBROUTINE guide_init
 
+    USE control_mod
+
     IMPLICIT NONE
   
@@ -67,5 +69,4 @@
     INCLUDE "paramet.h"
     INCLUDE "netcdf.inc"
-    INCLUDE "control.h"
 
     INTEGER                :: error,ncidpl,rid,rcod
@@ -269,4 +270,6 @@
 !=======================================================================
   SUBROUTINE guide_main(itau,ucov,vcov,teta,q,masse,ps)
+
+    USE control_mod
  
     IMPLICIT NONE
@@ -274,5 +277,4 @@
     INCLUDE "dimensions.h"
     INCLUDE "paramet.h"
-    INCLUDE "control.h"
     INCLUDE "comconst.h"
     INCLUDE "comvert.h"
@@ -354,7 +356,7 @@
       dday_step=real(day_step)
       IF (iguide_read.LT.0) THEN
-          tau=ditau/dday_step/FLOAT(iguide_read)
+          tau=ditau/dday_step/REAL(iguide_read)
       ELSE
-          tau=FLOAT(iguide_read)*ditau/dday_step
+          tau=REAL(iguide_read)*ditau/dday_step
       ENDIF
       reste=tau-AINT(tau)
@@ -541,5 +543,5 @@
             ENDDO
         ENDDO 
-        fieldm(:,l)=fieldm(:,l)/FLOAT(imax(typ)-imin(typ)+1)
+        fieldm(:,l)=fieldm(:,l)/REAL(imax(typ)-imin(typ)+1)
     ! Compute forcing
         DO j=1,hsize
Index: LMDZ4/trunk/libf/dyn3d/infotrac.F90
===================================================================
--- LMDZ4/trunk/libf/dyn3d/infotrac.F90	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/infotrac.F90	(revision 1403)
@@ -31,4 +31,7 @@
 
   SUBROUTINE infotrac_init
+
+    USE control_mod
+ 
     IMPLICIT NONE
 !=======================================================================
@@ -49,5 +52,4 @@
 
     INCLUDE "dimensions.h"
-    INCLUDE "control.h"
     INCLUDE "iniprint.h"
 
@@ -217,5 +219,5 @@
           new_iq=new_iq+10 ! 9 tracers added
        ELSE
-          WRITE(lunout,*) 'This choice of advection schema is not available'
+          WRITE(lunout,*) 'This choice of advection schema is not available',iq,hadv(iq),vadv(iq)
           CALL abort_gcm('infotrac_init','Bad choice of advection schema - 1',1)
        END IF
@@ -258,5 +260,5 @@
           iadv(new_iq)=11
        ELSE
-          WRITE(lunout,*)'This choice of advection schema is not available'
+          WRITE(lunout,*)'This choice of advection schema is not available',iq,hadv(iq),vadv(iq)
           CALL abort_gcm('infotrac_init','Bad choice of advection schema - 2',1)
        END IF
Index: LMDZ4/trunk/libf/dyn3d/iniacademic.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/iniacademic.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/iniacademic.F	(revision 1403)
@@ -1,4 +1,4 @@
 !
-! $Header$
+! $Id$
 !
 c
@@ -8,4 +8,6 @@
       USE filtreg_mod
       USE infotrac, ONLY : nqtot
+      USE control_mod
+ 
 
 c%W%    %G%
@@ -44,6 +46,6 @@
 #include "ener.h"
 #include "temps.h"
-#include "control.h"
 #include "iniprint.h"
+#include "logic.h"
 
 c   Arguments:
@@ -84,5 +86,5 @@
         time_0=0.
         day_ref=0
-	annee_ref=0
+        annee_ref=0
 
         im         = iim
@@ -93,5 +95,5 @@
         g      = 9.8
         daysec = 86400.
-        dtvr    = daysec/FLOAT(day_step)
+        dtvr    = daysec/REAL(day_step)
         zdtvr=dtvr
         kappa  = 0.2857143
@@ -105,13 +107,31 @@
         ang0       = 0.
 
+        if (llm.eq.1) then
+          ! specific initializations for the shallow water case
+          kappa=1
+        endif
+        
         CALL iniconst
         CALL inigeom
         CALL inifilr
 
-        ps=0.
-        phis=0.
+        if (llm.eq.1) then
+          ! initialize fields for the shallow water case, if required
+          if (.not.read_start) then
+            phis(:)=0.
+            q(:,:,1)=1.e-10
+            q(:,:,2)=1.e-15
+            q(:,:,3:nqtot)=0.
+            CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
+          endif
+        endif
+
+        if (iflag_phys.eq.2) then
+          ! initializations for the academic case
+          ps(:)=1.e5
+          phis(:)=0.
 c---------------------------------------------------------------------
 
-        taurappel=10.*daysec
+          taurappel=10.*daysec
 
 c---------------------------------------------------------------------
@@ -119,63 +139,61 @@
 c   --------------------------------------
 
-        DO l=1,llm
-         zsig=ap(l)/preff+bp(l)
-         if (zsig.gt.0.3) then
-           lsup=l
-           tetarappell=1./8.*(-log(zsig)-.5)
-           DO j=1,jjp1
+          DO l=1,llm
+            zsig=ap(l)/preff+bp(l)
+            if (zsig.gt.0.3) then
+             lsup=l
+             tetarappell=1./8.*(-log(zsig)-.5)
+             DO j=1,jjp1
              ddsin=sin(rlatu(j))-sin(pi/20.)
              tetajl(j,l)=300.*(1+1./18.*(1.-3.*ddsin*ddsin)+tetarappell)
-           ENDDO
-          else
+             ENDDO
+            else
 c   Choix isotherme au-dessus de 300 mbar
-           do j=1,jjp1
-             tetajl(j,l)=tetajl(j,lsup)*(0.3/zsig)**kappa
-           enddo
-          endif ! of if (zsig.gt.0.3)
-        ENDDO ! of DO l=1,llm
-
-        do l=1,llm
-           do j=1,jjp1
+             do j=1,jjp1
+               tetajl(j,l)=tetajl(j,lsup)*(0.3/zsig)**kappa
+             enddo
+            endif ! of if (zsig.gt.0.3)
+          ENDDO ! of DO l=1,llm
+
+          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
+            enddo
+          enddo
 
 c       call dump2d(jjp1,llm,tetajl,'TEQ   ')
 
-        ps=1.e5
-        phis=0.
-        CALL pression ( ip1jmp1, ap, bp, ps, p       )
-        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-        CALL massdair(p,masse)
+          CALL pression ( ip1jmp1, ap, bp, ps, p       )
+          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+          CALL massdair(p,masse)
 
 c  intialisation du vent et de la temperature
-        teta(:,:)=tetarappel(:,:)
-        CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
-        call ugeostr(phi,ucov)
-        vcov=0.
-        q(:,:,1   )=1.e-10
-        q(:,:,2   )=1.e-15
-        q(:,:,3:nqtot)=0.
+          teta(:,:)=tetarappel(:,:)
+          CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
+          call ugeostr(phi,ucov)
+          vcov=0.
+          q(:,:,1   )=1.e-10
+          q(:,:,2   )=1.e-15
+          q(:,:,3:nqtot)=0.
 
 
 c   perturbation aleatoire sur la temperature
-        idum  = -1
-        zz = ran1(idum)
-        idum  = 0
-        do l=1,llm
-           do ij=iip2,ip1jm
+          idum  = -1
+          zz = ran1(idum)
+          idum  = 0
+          do l=1,llm
+            do ij=iip2,ip1jm
               teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
-           enddo
-        enddo
-
-        do l=1,llm
-           do ij=1,ip1jmp1,iip1
+            enddo
+          enddo
+
+          do l=1,llm
+            do ij=1,ip1jmp1,iip1
               teta(ij+iim,l)=teta(ij,l)
-           enddo
-        enddo
+            enddo
+          enddo
 
 
@@ -187,9 +205,10 @@
 
 c   initialisation d'un traceur sur une colonne
-        j=jjp1*3/4
-        i=iip1/2
-        ij=(j-1)*iip1+i
-        q(ij,:,3)=1.
-      
+          j=jjp1*3/4
+          i=iip1/2
+          ij=(j-1)*iip1+i
+          q(ij,:,3)=1.
+        endif ! of if (iflag_phys.eq.2)
+        
       else
         write(lunout,*)"iniacademic: planet types other than earth",
Index: LMDZ4/trunk/libf/dyn3d/iniconst.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/iniconst.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/iniconst.F	(revision 1403)
@@ -1,6 +1,8 @@
 !
-! $Header$
+! $Id$
 !
       SUBROUTINE iniconst
+
+      USE control_mod
 
       IMPLICIT NONE
@@ -16,6 +18,6 @@
 #include "comconst.h"
 #include "temps.h"
-#include "control.h"
 #include "comvert.h"
+#include "iniprint.h"
 
 
@@ -47,5 +49,5 @@
       r       = cpp * kappa
 
-      PRINT*,' R  CP  Kappa ',  r , cpp,  kappa
+      write(lunout,*)'iniconst: R  CP  Kappa ',  r , cpp,  kappa
 c
 c-----------------------------------------------------------------------
Index: LMDZ4/trunk/libf/dyn3d/inidissip.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/inidissip.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/inidissip.F	(revision 1403)
@@ -11,4 +11,6 @@
 c   -------------
 
+      USE control_mod
+
       IMPLICIT NONE
 #include "dimensions.h"
@@ -17,5 +19,4 @@
 #include "comconst.h"
 #include "comvert.h"
-#include "control.h"
 #include "logic.h"
 
@@ -165,5 +166,5 @@
 
 c     IF(.NOT.lstardis) THEN
-         fact    = rad*24./float(jjm)
+         fact    = rad*24./REAL(jjm)
          fact    = fact*fact
          PRINT*,'coef u ', fact/cdivu, 1./cdivu
Index: LMDZ4/trunk/libf/dyn3d/inigeom.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/inigeom.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/inigeom.F	(revision 1403)
@@ -1,4 +1,4 @@
 !
-! $Header$
+! $Id$
 !
 c
@@ -168,15 +168,15 @@
 c
       IF( nitergdiv.NE.2 ) THEN
-        gamdi_gdiv = coefdis/ ( float(nitergdiv) -2. )
+        gamdi_gdiv = coefdis/ ( REAL(nitergdiv) -2. )
       ELSE
         gamdi_gdiv = 0.
       ENDIF
       IF( nitergrot.NE.2 ) THEN
-        gamdi_grot = coefdis/ ( float(nitergrot) -2. )
+        gamdi_grot = coefdis/ ( REAL(nitergrot) -2. )
       ELSE
         gamdi_grot = 0.
       ENDIF
       IF( niterh.NE.2 ) THEN
-        gamdi_h = coefdis/ ( float(niterh) -2. )
+        gamdi_h = coefdis/ ( REAL(niterh) -2. )
       ELSE
         gamdi_h = 0.
@@ -381,6 +381,6 @@
        yprp               = yprimu2(j-1)
        rlatp              = rlatu2 (j-1)
-ccc       yprp             = fyprim( FLOAT(j) - 0.25 )
-ccc       rlatp            = fy    ( FLOAT(j) - 0.25 )
+ccc       yprp             = fyprim( REAL(j) - 0.25 )
+ccc       rlatp            = fy    ( REAL(j) - 0.25 )
 c
       coslatp             = COS( rlatp )
@@ -416,8 +416,8 @@
         rlatm    = rlatu1 (  j  )
         yprm     = yprimu1(  j  )
-cc         rlatp    = fy    ( FLOAT(j) - 0.25 )
-cc         yprp     = fyprim( FLOAT(j) - 0.25 )
-cc         rlatm    = fy    ( FLOAT(j) + 0.25 )
-cc         yprm     = fyprim( FLOAT(j) + 0.25 )
+cc         rlatp    = fy    ( REAL(j) - 0.25 )
+cc         yprp     = fyprim( REAL(j) - 0.25 )
+cc         rlatm    = fy    ( REAL(j) + 0.25 )
+cc         yprm     = fyprim( REAL(j) + 0.25 )
 
          coslatm  = COS( rlatm )
Index: LMDZ4/trunk/libf/dyn3d/integrd.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/integrd.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/integrd.F	(revision 1403)
@@ -5,4 +5,6 @@
      $  (  nq,vcovm1,ucovm1,tetam1,psm1,massem1,
      $     dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis,finvmaold )
+
+      USE control_mod
 
       IMPLICIT NONE
@@ -32,5 +34,4 @@
 #include "temps.h"
 #include "serre.h"
-#include "control.h"
 
 c   Arguments:
Index: LMDZ4/trunk/libf/dyn3d/interp_horiz.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/interp_horiz.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/interp_horiz.F	(revision 1403)
@@ -1,4 +1,4 @@
 c
-c $Header$
+c $Id$
 c
       subroutine interp_horiz (varo,varn,imo,jmo,imn,jmn,lm,
@@ -101,6 +101,6 @@
            end do 
            do ii =1, imn+1
-             varn(ii,1,l) = totn/float(imn+1)
-             varn(ii,jmn+1,l) = tots/float(imn+1)
+             varn(ii,1,l) = totn/REAL(imn+1)
+             varn(ii,jmn+1,l) = tots/REAL(imn+1)
            end do 
        end do
Index: LMDZ4/trunk/libf/dyn3d/interpre.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/interpre.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/interpre.F	(revision 1403)
@@ -1,8 +1,10 @@
 !
-! $Header$
+! $Id$
 !
        subroutine interpre(q,qppm,w,fluxwppm,masse,
      s            apppm,bpppm,massebx,masseby,pbaru,pbarv,
      s            unatppm,vnatppm,psppm)
+
+      USE control_mod
 
        implicit none
@@ -17,5 +19,4 @@
 #include "logic.h"
 #include "temps.h"
-#include "control.h"
 #include "ener.h"
 #include "description.h"
Index: LMDZ4/trunk/libf/dyn3d/juldate.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/juldate.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/juldate.F	(revision 1403)
@@ -1,4 +1,4 @@
 !
-! $Header$
+! $Id$
 !
 	subroutine juldate(ian,imoi,ijou,oh,om,os,tjd,tjdsec)
Index: LMDZ4/trunk/libf/dyn3d/leapfrog.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/leapfrog.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/leapfrog.F	(revision 1403)
@@ -15,4 +15,5 @@
       USE guide_mod, ONLY : guide_main
       USE write_field
+      USE control_mod
       IMPLICIT NONE
 
@@ -56,9 +57,8 @@
 #include "logic.h"
 #include "temps.h"
-#include "control.h"
 #include "ener.h"
 #include "description.h"
 #include "serre.h"
-#include "com_io_dyn.h"
+!#include "com_io_dyn.h"
 #include "iniprint.h"
 #include "academic.h"
@@ -197,10 +197,10 @@
 
       itau = 0
-c$$$      iday = day_ini+itau/day_step
-c$$$      time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
-c$$$         IF(time.GT.1.) THEN
-c$$$          time = time-1.
-c$$$          iday = iday+1
-c$$$         ENDIF
+c      iday = day_ini+itau/day_step
+c      time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
+c         IF(time.GT.1.) THEN
+c          time = time-1.
+c          iday = iday+1
+c         ENDIF
 
 
@@ -276,4 +276,5 @@
 
       IF( purmats ) THEN
+      ! Purely Matsuno time stepping
          IF( MOD(itau,iconser) .EQ.0.AND.  forward    ) conser = .TRUE.
          IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
@@ -281,8 +282,15 @@
      s          .and. iflag_phys.EQ.1                 ) apphys = .TRUE.
       ELSE
+      ! Leapfrog/Matsuno time stepping 
          IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
          IF( MOD(itau+1,idissip)  .EQ. 0              ) apdiss = .TRUE.
          IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE.
       END IF
+
+! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
+!          supress dissipation step
+      if (llm.eq.1) then
+        apdiss=.false.
+      endif
 
 c-----------------------------------------------------------------------
@@ -522,10 +530,10 @@
             IF(forward. OR. leapf) THEN
               itau= itau + 1
-c$$$              iday= day_ini+itau/day_step
-c$$$              time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
-c$$$                IF(time.GT.1.) THEN
-c$$$                  time = time-1.
-c$$$                  iday = iday+1
-c$$$                ENDIF
+c              iday= day_ini+itau/day_step
+c              time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
+c                IF(time.GT.1.) THEN
+c                  time = time-1.
+c                  iday = iday+1
+c                ENDIF
             ENDIF
 
@@ -559,12 +567,16 @@
                IF (ok_dynzon) THEN
 #ifdef CPP_IOIPSL
-!                  CALL writedynav(histaveid, itau,vcov ,
-!     ,                 ucov,teta,pk,phi,q,masse,ps,phis)
-                  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
                END IF
-
-            ENDIF
+               IF (ok_dyn_ave) THEN
+#ifdef CPP_IOIPSL
+                 CALL writedynav(itau,vcov,
+     &                 ucov,teta,pk,phi,q,masse,ps,phis)
+#endif
+               ENDIF
+
+            ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
 
 c-----------------------------------------------------------------------
@@ -572,7 +584,7 @@
 c   ------------------------------
 
-            IF( MOD(itau,iecri         ).EQ.0) THEN
-c           IF( MOD(itau,iecri*day_step).EQ.0) THEN
-
+            IF( MOD(itau,iecri).EQ.0) THEN
+             ! Ehouarn: output only during LF or Backward Matsuno
+	     if (leapf.or.(.not.leapf.and.(.not.forward))) then
               nbetat = nbetatdem
               CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
@@ -583,12 +595,19 @@
               enddo
 #ifdef CPP_IOIPSL
-c             CALL writehist(histid,histvid,itau,vcov, 
-c     &                      ucov,teta,phi,q,masse,ps,phis)
+              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)
 #endif
 ! For some Grads outputs of fields
-             if (output_grads_dyn) then
+              if (output_grads_dyn) then
 #include "write_grads_dyn.h"
-             endif
-
+              endif
+             endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
             ENDIF ! of IF(MOD(itau,iecri).EQ.0)
 
@@ -645,11 +664,11 @@
 
              itau =  itau + 1
-c$$$             iday = day_ini+itau/day_step
-c$$$             time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
-c$$$
-c$$$                  IF(time.GT.1.) THEN
-c$$$                   time = time-1.
-c$$$                   iday = iday+1
-c$$$                  ENDIF
+c             iday = day_ini+itau/day_step
+c             time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
+c
+c                  IF(time.GT.1.) THEN
+c                   time = time-1.
+c                   iday = iday+1
+c                  ENDIF
 
                forward =  .FALSE.
@@ -660,5 +679,5 @@
                GO TO 2
 
-            ELSE ! of IF(forward)
+            ELSE ! of IF(forward) i.e. backward step
 
               IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
@@ -671,10 +690,14 @@
                IF (ok_dynzon) THEN 
 #ifdef CPP_IOIPSL
-!                  CALL writedynav(histaveid, itau,vcov ,
-!     ,                 ucov,teta,pk,phi,q,masse,ps,phis)
-                  CALL bilan_dyn (2,dtvr*iperiod,dtvr*day_step*periodav,
-     ,                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
-#endif
-               END IF
+                 CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav,
+     &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
+#endif
+               ENDIF
+               IF (ok_dyn_ave) THEN
+#ifdef CPP_IOIPSL
+                 CALL writedynav(itau,vcov,
+     &                 ucov,teta,pk,phi,q,masse,ps,phis)
+#endif
+               ENDIF
 
               ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
@@ -690,6 +713,9 @@
                 enddo
 #ifdef CPP_IOIPSL
-c               CALL writehist( histid, histvid, itau,vcov , 
-c    &                           ucov,teta,phi,q,masse,ps,phis)
+              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)
 #endif
 ! For some Grads outputs
Index: LMDZ4/trunk/libf/dyn3d/limit_netcdf.F90
===================================================================
--- LMDZ4/trunk/libf/dyn3d/limit_netcdf.F90	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/limit_netcdf.F90	(revision 1403)
@@ -20,4 +20,5 @@
 !  *    12/2009: D. Cugnet   (f77->f90, calendars, files from coupled runs)
 !-------------------------------------------------------------------------------
+  USE control_mod
 #ifdef CPP_EARTH
   USE dimphy
@@ -27,5 +28,5 @@
                    NF90_DEF_DIM, NF90_DEF_VAR, NF90_PUT_VAR, NF90_PUT_ATT,     &
                    NF90_NOERR,   NF90_NOWRITE, NF90_DOUBLE,  NF90_GLOBAL,      &
-		   NF90_CLOBBER, NF90_ENDDEF,  NF90_UNLIMITED
+		   NF90_CLOBBER, NF90_ENDDEF,  NF90_UNLIMITED, NF90_FLOAT
   USE inter_barxy_m, only: inter_barxy
 #endif
@@ -45,5 +46,4 @@
 !-------------------------------------------------------------------------------
 ! Local variables:
-#include "control.h"
 #include "logic.h"
 #include "comvert.h"
@@ -293,9 +293,9 @@
   USE dimphy, ONLY : klon
   USE phys_state_var_mod, ONLY : pctsrf
+  USE control_mod
   IMPLICIT NONE
 #include "dimensions.h"
 #include "paramet.h"
 #include "comgeom2.h"
-#include "control.h"
 #include "indicesol.h"
 #include "iniprint.h"
Index: LMDZ4/trunk/libf/dyn3d/ord_coord.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/ord_coord.F	(revision 1340)
+++ 	(revision )
@@ -1,95 +1,0 @@
-!
-! $Header$
-!
-       SUBROUTINE ord_coord ( nmax, xi, xo, decrois )
-
-c   .... Auteur :  P. Le Van  .... 
-c
-c   ... Reordonne eventuellement les coordonnees de la grille donnees ...
-c
-       IMPLICIT NONE
-
-c    .....  Arguments  en entree  .....
-
-       INTEGER nmax
-       REAL xi(nmax)
-
-c    .....  Arguments  en sortie  .....
-c
-       REAL xo(nmax+1) 
-       LOGICAL decrois
-
-c    .... Variables  locales  ....
-
-       REAL xscr(nmax)
-       INTEGER i,ii
-       REAL pi, degres, chmin, chmax, mult
-c
-
-       pi     = 2.*ASIN(1.)
-       degres = 180./pi
-       decrois = .FALSE.
-
-       DO i = 1, nmax 
-        xo(i) = xi(i)
-       ENDDO
-     
-       mult = 1.
-       IF( xo(1).GT.xo(nmax) ) mult = -1.
-
-       CALL minmax(nmax,xo(1),chmin,chmax)
-
-       IF(chmax.LT.6.5 )  THEN
-          DO i = 1,nmax
-           xo(i) = xo(i) * degres 
-          ENDDO
-       ENDIF
-
-       IF( ABS( xo( 1  ) + mult* 90. ). LT .0.001. OR .
-     ,     ABS( xo(nmax) - mult* 90. ). LT .0.001      )  THEN
-        PRINT *,' Reverifier les valeurs de  xidat  pour les donnees .'
-        PRINT *,' Elles doivent correspondre aux interfaces et non aux',
-     , 'ordonnees des donnees,egales a -90. et 90.deg aux 2 extremites '
-         CALL ABORT
-       ENDIF
-
-       IF( xo(1).GT.xo(nmax) )   THEN
-          DO i = 1, nmax 
-            xscr(i) = xo(i)
-          ENDDO
-          DO i = 1, nmax 
-            xo(i+1) = xscr(i)
-          ENDDO
-            xo (   1    ) =   90.
-       ELSE
-            xo ( nmax +1) =   90.
-       ENDIF
-
-       IF ( xo(2).LT.xo(1) ) decrois =.TRUE.
-
-       DO i = 3, nmax 
-
-        IF(decrois.AND.xo(i).GT.xo(i-1) ) THEN
-         PRINT 1
-         PRINT 2,(xo(ii),ii=1,nmax)
-         CALL ABORT 
-        ENDIF
-        IF(.NOT.decrois.AND.xo(i).LT.xo(i-1) ) THEN
-         PRINT 1
-         PRINT 2,(xo(ii),ii=1,nmax)
-         CALL ABORT
-        ENDIF
-
-       ENDDO
-        
-       IF( decrois )  THEN
-c         CALL sort(nmax+1,xo(1))
-        CALL sort(nmax+1,xo)
-       ENDIF
-
-1      FORMAT(5x,' Incoherence dans les valeurs des latitudes de la ',
-     ,  'grille du modele ')
-2      FORMAT(1x,8f8.2)
-
-       RETURN
-       END
Index: LMDZ4/trunk/libf/dyn3d/ord_coordm.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/ord_coordm.F	(revision 1340)
+++ 	(revision )
@@ -1,110 +1,0 @@
-!
-! $Header$
-!
-       SUBROUTINE ord_coordm ( nmax, xi, xo, jjm, jmods, decrois )
-
-c   ....  Auteur :  P. Le Van  .... 
-
-c   ... Reordonne eventuellement les coordonnees de la grille modele ...
-c
-       IMPLICIT NONE
-
-c    .....  Arguments  en entree  .....
-
-       INTEGER nmax,jjm
-       REAL xi(nmax)
-
-c    .....  Arguments  en sortie  .....
-c
-       REAL xo(nmax+1) 
-       LOGICAL decrois
-       INTEGER jmods
-
-c    .... Variables  locales  ....
-
-       REAL xscr(nmax)
-       INTEGER i
-       REAL pi, degres, chmin, chmax,mult
-c
-       DO i = 1, nmax 
-        xo(i) = xi(i)
-       ENDDO
-
-       mult = 1.
-       IF( xo(1).GT.xo(nmax) )  mult = - 1.
-       IF( nmax.EQ.jjm    ) jmods = nmax +1
-       IF( nmax.EQ.jjm +1 ) jmods = nmax -1 
-     
-       pi     = 2.*ASIN(1.)
-       degres = 180./pi
-       decrois = .FALSE.
-
-       CALL minmax(nmax,xo(1),chmin,chmax)
-
-       IF(chmax.LT.6.5 )  THEN
-          DO i = 1,nmax
-           xo(i) = xo(i) * degres 
-          ENDDO
-       ENDIF
-
-       IF( nmax.EQ.jjm )   THEN
-         IF( xo(1).GT.xo(nmax) )   THEN
-           DO i = 1, nmax 
-            xscr(i) = xo(i)
-           ENDDO
-           DO i = 1, nmax 
-            xo(i+1) = xscr(i)
-           ENDDO
-            xo (   1    ) =   90.
-         ELSE
-            xo ( nmax+1 ) =   90.
-         ENDIF
-       ELSE
-          IF( nmax.NE.jjm +1 )   THEN
-             PRINT *,'  Dans la routine ord_coordm , l argument nmax '
-             PRINT *,'  n est pas egal a jjm ni a jjm +1 . Corriger !'
-             CALL ABORT
-          ELSE
-            IF( ABS( xo(1)+ mult * 90.).GT.0.01 )  THEN
-              PRINT *,' Avec nmax =',nmax,'on devrait avoir des',
-     ,    ' ordonnees = 90. deg pour j=1 ou jjm+1 ! '
-             CALL ABORT
-            ELSE
-               IF( xo(1).LT.xo(nmax) )  THEN
-                 DO i = 1, nmax
-                  xscr(i) = xo(i)
-                 ENDDO
-                 DO i = 1, nmax -1
-                  xo(i) = xscr(i+1)
-                 ENDDO
-               ENDIF
-            ENDIF
-          ENDIF
-       ENDIF
-
-       IF ( xo(2).LT.xo(1) ) decrois =.TRUE.
-
-       DO i = 3, nmax 
-
-        IF(decrois.AND.xo(i).GT.xo(i-1) ) THEN
-         PRINT 1
-         CALL ABORT 
-        ENDIF
-        IF(.NOT.decrois.AND.xo(i).LT.xo(i-1) ) THEN
-         PRINT 1
-         CALL ABORT
-        ENDIF
-
-       ENDDO
-        
-       IF( decrois )  THEN
-         CALL sort(jmods,xo(1))
-       ENDIF
-
-
-1      FORMAT(5x,' Incoherence dans les valeurs des latitudes de la ',
-     ,  'grille du modele ')
-2      FORMAT(1x,8f8.2)
-
-       RETURN
-       END
Index: LMDZ4/trunk/libf/dyn3d/ppm3d.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/ppm3d.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/ppm3d.F	(revision 1403)
@@ -1,4 +1,4 @@
 !
-! $Header$
+! $Id$
 !
 
@@ -345,6 +345,6 @@
 C
       PI = 4. * ATAN(1.)
-      DL = 2.*PI / float(IMR)
-      DP =    PI / float(JMR)
+      DL = 2.*PI / REAL(IMR)
+      DP =    PI / REAL(JMR)
 C
       if(IGD.eq.0) then
@@ -388,5 +388,5 @@
       ZTC  = acos(CR1) * (180./PI)
 C
-      JS0 = float(JMR)*(90.-ZTC)/180. + 2
+      JS0 = REAL(JMR)*(90.-ZTC)/180. + 2
       JS0 = max(JS0, J1+1)
       IML = min(6*JS0/(J1-1)+2, 4*IMR/5)
@@ -628,5 +628,5 @@
 C Contribution from the N-S advection
       do i=1,imr*(j2-j1+1)
-      JT = float(J1) - VA(i,j1)
+      JT = REAL(J1) - VA(i,j1)
       wk1(i,j1,2) = VA(i,j1) * (q(i,jt,k,IC) - q(i,jt+1,k,IC))
       enddo
@@ -949,5 +949,5 @@
       IF(IORD.eq.1 .or. j.eq.j1. or. j.eq.j2) THEN
       DO 1406 i=1,IMR
-      iu = float(i) - uc(i,j)
+      iu = REAL(i) - uc(i,j)
 1406  fx1(i) = qtmp(iu)
       ELSE
@@ -957,5 +957,5 @@
       if(IORD.eq.2 .or. j.le.j1vl .or. j.ge.j2vl) then
       DO 1408 i=1,IMR
-      iu = float(i) - uc(i,j)
+      iu = REAL(i) - uc(i,j)
 1408  fx1(i) = qtmp(iu) + DC(iu)*(sign(1.,uc(i,j))-uc(i,j))
       else
@@ -1111,5 +1111,5 @@
       if(JORD.eq.1) then
       DO 1000 i=1,len
-      JT = float(J1) - VC(i,J1)
+      JT = REAL(J1) - VC(i,J1)
 1000  fx(i,j1) = p(i,JT)
       else
@@ -1123,5 +1123,5 @@
       else
       DO 1200 i=1,len
-      JT = float(J1) - VC(i,J1)
+      JT = REAL(J1) - VC(i,J1)
 1200  fx(i,j1) = p(i,JT) + (sign(1.,VC(i,j1))-VC(i,j1))*DC2(i,JT)
       endif
@@ -1358,5 +1358,5 @@
 	do j=j1-1,j2+1
       do i=1,imr
-      JP = float(j)-VA(i,j)
+      JP = REAL(j)-VA(i,j)
       ady(i,j) = VA(i,j)*(wk(i,jp)-wk(i,jp+1))
       enddo
@@ -1582,5 +1582,5 @@
       JMR = JNP-1
       do 55 j=2,JNP
-        ph5  =  -0.5*PI + (FLOAT(J-1)-0.5)*DP
+        ph5  =  -0.5*PI + (REAL(J-1)-0.5)*DP
 55      cose(j) = cos(ph5)
 C
@@ -1834,5 +1834,5 @@
 C
 c      if(first) then
-      DP = 4.*ATAN(1.)/float(JNP-1)
+      DP = 4.*ATAN(1.)/REAL(JNP-1)
       CAP1 = IMR*(1.-COS((j1-1.5)*DP))/DP
 c      first = .false.
@@ -1889,5 +1889,5 @@
 C Check Poles.
       if(q(1,1).lt.0.) then
-      dq = q(1,1)*cap1/float(IMR)*acosp(j1)
+      dq = q(1,1)*cap1/REAL(IMR)*acosp(j1)
       do i=1,imr
       q(i,1) = 0.
@@ -1898,5 +1898,5 @@
 C
       if(q(1,JNP).lt.0.) then
-      dq = q(1,JNP)*cap1/float(IMR)*acosp(j2)
+      dq = q(1,JNP)*cap1/REAL(IMR)*acosp(j2)
       do i=1,imr
       q(i,JNP) = 0.
Index: LMDZ4/trunk/libf/dyn3d/ran1.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/ran1.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/ran1.F	(revision 1403)
@@ -1,4 +1,4 @@
 !
-! $Header$
+! $Id$
 !
       FUNCTION RAN1(IDUM)
@@ -20,5 +20,5 @@
           IX1=MOD(IA1*IX1+IC1,M1)
           IX2=MOD(IA2*IX2+IC2,M2)
-          R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
+          R(J)=(REAL(IX1)+REAL(IX2)*RM2)*RM1
 11      CONTINUE
         IDUM=1
@@ -30,5 +30,5 @@
       IF(J.GT.97.OR.J.LT.1)PAUSE
       RAN1=R(J)
-      R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
+      R(J)=(REAL(IX1)+REAL(IX2)*RM2)*RM1
       RETURN
       END
Index: LMDZ4/trunk/libf/dyn3d/sortvarc.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/sortvarc.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/sortvarc.F	(revision 1403)
@@ -1,4 +1,4 @@
 !
-! $Header$
+! $Id$
 !
       SUBROUTINE sortvarc
@@ -59,5 +59,5 @@
 
        dtvrs1j   = dtvr/daysec
-       rjour     = FLOAT( INT( itau * dtvrs1j ))
+       rjour     = REAL( INT( itau * dtvrs1j ))
        heure     = ( itau*dtvrs1j-rjour ) * 24.
        imjmp1    = iim * jjp1
@@ -129,7 +129,7 @@
       ang   = SSUM(     llm,  angl, 1 )
 
-c      rday = FLOAT(INT ( day_ini + time ))
+c      rday = REAL(INT ( day_ini + time ))
 c
-       rday = FLOAT(INT(time-jD_ref-jH_ref))
+       rday = REAL(INT(time-jD_ref-jH_ref))
       IF(ptot0.eq.0.)  THEN
          PRINT 3500, itau, rday, heure,time
Index: LMDZ4/trunk/libf/dyn3d/sortvarc0.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/sortvarc0.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/sortvarc0.F	(revision 1403)
@@ -1,4 +1,4 @@
 !
-! $Header$
+! $Id$
 !
       SUBROUTINE sortvarc0
@@ -60,5 +60,5 @@
 
        dtvrs1j   = dtvr/daysec
-       rjour     = FLOAT( INT( itau * dtvrs1j ))
+       rjour     = REAL( INT( itau * dtvrs1j ))
        heure     = ( itau*dtvrs1j-rjour ) * 24.
        imjmp1    = iim * jjp1
@@ -130,5 +130,5 @@
       ang0   = SSUM(     llm,  angl, 1 )
 
-      rday = FLOAT(INT (time ))
+      rday = REAL(INT (time ))
 c
       PRINT 3500, itau, rday, heure, time
Index: LMDZ4/trunk/libf/dyn3d/sw_case_williamson91_6.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/sw_case_williamson91_6.F	(revision 1403)
+++ LMDZ4/trunk/libf/dyn3d/sw_case_williamson91_6.F	(revision 1403)
@@ -0,0 +1,140 @@
+!
+! $Id $
+!
+      SUBROUTINE sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
+
+c=======================================================================
+c
+c   Author:    Thomas Dubos      original: 26/01/2010
+c   -------
+c
+c   Subject:
+c   ------
+c   Realise le cas-test 6 de Williamson et al. (1991) : onde de Rossby-Haurwitz
+c
+c   Method:
+c   --------
+c
+c   Interface:
+c   ----------
+c
+c      Input:
+c      ------
+c
+c      Output:
+c      -------
+c
+c=======================================================================
+      IMPLICIT NONE
+c-----------------------------------------------------------------------
+c   Declararations:
+c   ---------------
+
+#include "dimensions.h"
+#include "paramet.h"
+#include "comvert.h"
+#include "comconst.h"
+#include "comgeom.h"
+#include "iniprint.h"
+
+c   Arguments:
+c   ----------
+
+c   variables dynamiques
+      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
+      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
+      REAL ps(ip1jmp1)                       ! pression  au sol
+      REAL masse(ip1jmp1,llm)                ! masse d'air
+      REAL phis(ip1jmp1)                     ! geopotentiel au sol
+
+c   Local:
+c   ------
+
+      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 pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
+      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
+
+      REAL :: sinth,costh,costh2, Ath,Bth,Cth, lon,dps
+      INTEGER i,j,ij
+
+      REAL, PARAMETER    :: rho=1 ! masse volumique de l'air (arbitraire)
+      REAL, PARAMETER    :: K    = 7.848e-6  ! K = \omega
+      REAL, PARAMETER    :: gh0  = 9.80616 * 8e3 
+      INTEGER, PARAMETER :: R0=4, R1=R0+1, R2=R0+2         ! mode 4
+c NB : rad = 6371220 dans W91 (6371229 dans LMDZ)
+c      omeg = 7.292e-5 dans W91 (7.2722e-5 dans LMDZ)
+ 
+      IF(0==0) THEN
+c Williamson et al. (1991) : onde de Rossby-Haurwitz
+         teta = preff/rho/cpp
+c geopotentiel (pression de surface)
+         do j=1,jjp1
+            costh2 = cos(rlatu(j))**2
+            Ath = (R0+1)*(costh2**2) + (2*R0*R0-R0-2)*costh2 - 2*R0*R0
+            Ath = .25*(K**2)*(costh2**(R0-1))*Ath
+            Ath = .5*K*(2*omeg+K)*costh2 + Ath 
+            Bth = (R1*R1+1)-R1*R1*costh2
+            Bth = 2*(omeg+K)*K/(R1*R2) * (costh2**(R0/2))*Bth
+            Cth = R1*costh2 - R2
+            Cth = .25*K*K*(costh2**R0)*Cth
+            do i=1,iip1
+               ij=(j-1)*iip1+i
+               lon = rlonv(i)
+               dps = Ath + Bth*cos(R0*lon) + Cth*cos(2*R0*lon)
+               ps(ij) = rho*(gh0 + (rad**2)*dps)
+            enddo
+         enddo
+         write(lunout,*) 'W91 ps', MAXVAL(ps), MINVAL(ps)
+c vitesse zonale ucov
+         do j=1,jjp1
+            costh  = cos(rlatu(j))
+            costh2 = costh**2
+            Ath = rad*K*costh
+            Bth = R0*(1-costh2)-costh2
+            Bth = rad*K*Bth*(costh**(R0-1))
+            do i=1,iip1
+               ij=(j-1)*iip1+i
+               lon = rlonu(i)
+               ucov(ij,1) = (Ath + Bth*cos(R0*lon))
+            enddo
+         enddo
+         write(lunout,*) 'W91 u', MAXVAL(ucov(:,1)), MINVAL(ucov(:,1))
+         ucov(:,1)=ucov(:,1)*cu
+c vitesse meridienne vcov
+         do j=1,jjm
+            sinth  = sin(rlatv(j))
+            costh  = cos(rlatv(j))
+            Ath = -rad*K*R0*sinth*(costh**(R0-1))
+            do i=1,iip1
+               ij=(j-1)*iip1+i
+               lon = rlonv(i)
+               vcov(ij,1) = Ath*sin(R0*lon)
+            enddo
+         enddo
+         write(lunout,*) 'W91 v', MAXVAL(vcov(:,1)), MINVAL(vcov(:,1))
+         vcov(:,1)=vcov(:,1)*cv
+         
+c         ucov=0
+c         vcov=0
+      ELSE
+c test non-tournant, onde se propageant en latitude
+         do j=1,jjp1
+            do i=1,iip1
+               ij=(j-1)*iip1+i
+               ps(ij) = 1e5*(1 + .1*exp(-100*(1+sin(rlatu(j)))**2) )
+            enddo
+         enddo
+         
+c     rho = preff/(cpp*teta)
+         teta = .01*preff/cpp   ! rho = 100 ; phi = ps/rho = 1e3 ; c=30 m/s = 2600 km/j = 23 degres / j
+         ucov=0.
+         vcov=0.
+      END IF      
+      
+      CALL pression ( ip1jmp1, ap, bp, ps, p       )
+      CALL massdair(p,masse)
+
+      END
+c-----------------------------------------------------------------------
Index: LMDZ4/trunk/libf/dyn3d/tourabs.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/tourabs.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/tourabs.F	(revision 1403)
@@ -57,7 +57,7 @@
         ELSE
          rot( ij,l ) = (vcov(ij+1,l)/cv(ij+1)-vcov(ij,l)/cv(ij))/
-     $                 (2.*pi*RAD*cos(rlatv(j)))*float(iim)
+     $                 (2.*pi*RAD*cos(rlatv(j)))*REAL(iim)
      $                +(ucov(ij+iip1,l)/cu(ij+iip1)-ucov(ij,l)/cu(ij))/
-     $                 (pi*RAD)*(float(jjm)-1.)
+     $                 (pi*RAD)*(REAL(jjm)-1.)
 c
         ENDIF
Index: LMDZ4/trunk/libf/dyn3d/traceurpole.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/traceurpole.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/traceurpole.F	(revision 1403)
@@ -1,6 +1,8 @@
 !
-! $Header$
+! $Id$
 !
           subroutine traceurpole(q,masse)
+
+      USE control_mod
 
           implicit none
@@ -15,5 +17,4 @@
 #include "logic.h"
 #include "temps.h"
-#include "control.h"
 #include "ener.h"
 #include "description.h"
Index: LMDZ4/trunk/libf/dyn3d/ugeostr.F
===================================================================
--- LMDZ4/trunk/libf/dyn3d/ugeostr.F	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/ugeostr.F	(revision 1403)
@@ -40,5 +40,5 @@
             DO i=1,iim
                u(i,j,l)=fact*(phi(i,j+1,l)-phi(i,j,l))
-               um(j,l)=um(j,l)+u(i,j,l)/float(iim)
+               um(j,l)=um(j,l)+u(i,j,l)/REAL(iim)
             ENDDO
          ENDDO
Index: LMDZ4/trunk/libf/dyn3d/write_paramLMDZ_dyn.h
===================================================================
--- LMDZ4/trunk/libf/dyn3d/write_paramLMDZ_dyn.h	(revision 1340)
+++ LMDZ4/trunk/libf/dyn3d/write_paramLMDZ_dyn.h	(revision 1403)
@@ -7,41 +7,41 @@
       itau_w=itau_dyn+itau
 c
-      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(prt_level) 
+      zx_tmp_2d(1:iip1,1:jjp1)=REAL(prt_level) 
       CALL histwrite(nid_ctesGCM, "prt_level", itau_w,
      .               zx_tmp_2d,iip1*jjp1,ndex2d)
 c
-      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(dayref)
+      zx_tmp_2d(1:iip1,1:jjp1)=REAL(dayref)
       CALL histwrite(nid_ctesGCM, "dayref", itau_w,
      .               zx_tmp_2d,iip1*jjp1,ndex2d)
 c
-      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(anneeref)
+      zx_tmp_2d(1:iip1,1:jjp1)=REAL(anneeref)
       CALL histwrite(nid_ctesGCM, "anneeref", itau_w,
      .               zx_tmp_2d,iip1*jjp1,ndex2d)
 c
-      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(raz_date)
+      zx_tmp_2d(1:iip1,1:jjp1)=REAL(raz_date)
       CALL histwrite(nid_ctesGCM, "raz_date", itau_w,
      .               zx_tmp_2d,iip1*jjp1,ndex2d)
 c
-      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(nday)
+      zx_tmp_2d(1:iip1,1:jjp1)=REAL(nday)
       CALL histwrite(nid_ctesGCM, "nday", itau_w,
      .               zx_tmp_2d,iip1*jjp1,ndex2d)
 c
-      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(day_step)
+      zx_tmp_2d(1:iip1,1:jjp1)=REAL(day_step)
       CALL histwrite(nid_ctesGCM, "day_step", itau_w,
      .               zx_tmp_2d,iip1*jjp1,ndex2d)
 c
-      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(iperiod)
+      zx_tmp_2d(1:iip1,1:jjp1)=REAL(iperiod)
       CALL histwrite(nid_ctesGCM, "iperiod", itau_w,
      .               zx_tmp_2d,iip1*jjp1,ndex2d)
 c
-      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(iapp_tracvl)
+      zx_tmp_2d(1:iip1,1:jjp1)=REAL(iapp_tracvl)
       CALL histwrite(nid_ctesGCM, "iapp_tracvl", itau_w,
      .               zx_tmp_2d,iip1*jjp1,ndex2d)
 c
-      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(iconser)
+      zx_tmp_2d(1:iip1,1:jjp1)=REAL(iconser)
       CALL histwrite(nid_ctesGCM, "iconser", itau_w,
      .               zx_tmp_2d,iip1*jjp1,ndex2d)
 c
-      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(iecri)
+      zx_tmp_2d(1:iip1,1:jjp1)=REAL(iecri)
       CALL histwrite(nid_ctesGCM, "iecri", itau_w,
      .               zx_tmp_2d,iip1*jjp1,ndex2d)
@@ -51,5 +51,5 @@
      .               zx_tmp_2d,iip1*jjp1,ndex2d)
 c
-      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(idissip)
+      zx_tmp_2d(1:iip1,1:jjp1)=REAL(idissip)
       CALL histwrite(nid_ctesGCM, "idissip", itau_w,
      .               zx_tmp_2d,iip1*jjp1,ndex2d)
@@ -63,13 +63,13 @@
      .               zx_tmp_2d,iip1*jjp1,ndex2d)
 c
-      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(nitergdiv)
+      zx_tmp_2d(1:iip1,1:jjp1)=REAL(nitergdiv)
       CALL histwrite(nid_ctesGCM, "nitergdiv", itau_w,
      .               zx_tmp_2d,iip1*jjp1,ndex2d)
 c
-      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(nitergrot)
+      zx_tmp_2d(1:iip1,1:jjp1)=REAL(nitergrot)
       CALL histwrite(nid_ctesGCM, "nitergrot", itau_w,
      .               zx_tmp_2d,iip1*jjp1,ndex2d)
 c
-      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(niterh) 
+      zx_tmp_2d(1:iip1,1:jjp1)=REAL(niterh) 
       CALL histwrite(nid_ctesGCM, "niterh", itau_w,
      .               zx_tmp_2d,iip1*jjp1,ndex2d)
@@ -118,9 +118,9 @@
      .               zx_tmp_2d,iip1*jjp1,ndex2d)
 c
-      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(iflag_phys)
+      zx_tmp_2d(1:iip1,1:jjp1)=REAL(iflag_phys)
       CALL histwrite(nid_ctesGCM, "iflag_phys", itau_w,
      .               zx_tmp_2d,iip1*jjp1,ndex2d)
 c
-      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(iphysiq)
+      zx_tmp_2d(1:iip1,1:jjp1)=REAL(iphysiq)
       CALL histwrite(nid_ctesGCM, "iphysiq", itau_w,
      .               zx_tmp_2d,iip1*jjp1,ndex2d)
