Index: /LMDZ5/trunk/libf/dyn3d/comconst.h
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/comconst.h	(revision 1670)
+++ /LMDZ5/trunk/libf/dyn3d/comconst.h	(revision 1671)
@@ -21,6 +21,7 @@
       REAL dtdiss ! (s) time step for the dissipation
       REAL rad ! (m) radius of the planet
-      REAL r ! Gas constant R=8.31 J.K-1.mol-1
-      REAL cpp   ! Cp
+      REAL r ! Reduced Gas constant r=R/mu
+             ! with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol)
+      REAL cpp   ! Specific heat Cp (J.kg-1.K-1)
       REAL kappa ! kappa=R/Cp 
       REAL cotot
Index: /LMDZ5/trunk/libf/dyn3d/gcm.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/gcm.F	(revision 1670)
+++ /LMDZ5/trunk/libf/dyn3d/gcm.F	(revision 1671)
@@ -433,6 +433,7 @@
 ! Physics:
 #ifdef CPP_PHYS
-         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys ,
-     ,                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp     )
+         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys,
+     &                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp,
+     &                iflag_phys)
 #endif
          call_iniphys=.false.
@@ -457,12 +458,4 @@
  301  FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
  302  FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
-#endif
-
-#ifdef CPP_PHYS
-! Create start file (startphy.nc) and boundary conditions (limit.nc) 
-! for the Earth verstion
-       if (iflag_phys>=100) then
-          call iniaqua(ngridmx,latfi,lonfi,iflag_phys)
-       endif
 #endif
 
Index: /LMDZ5/trunk/libf/dyn3dpar/comconst.h
===================================================================
--- /LMDZ5/trunk/libf/dyn3dpar/comconst.h	(revision 1670)
+++ /LMDZ5/trunk/libf/dyn3dpar/comconst.h	(revision 1671)
@@ -21,6 +21,7 @@
       REAL dtdiss ! (s) time step for the dissipation
       REAL rad ! (m) radius of the planet
-      REAL r ! Gas constant R=8.31 J.K-1.mol-1
-      REAL cpp   ! Cp
+      REAL r ! Reduced Gas constant r=R/mu
+             ! with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol)
+      REAL cpp   ! Specific heat Cp (J.kg-1.K-1)
       REAL kappa ! kappa=R/Cp 
       REAL cotot
Index: /LMDZ5/trunk/libf/dyn3dpar/gcm.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3dpar/gcm.F	(revision 1670)
+++ /LMDZ5/trunk/libf/dyn3dpar/gcm.F	(revision 1671)
@@ -423,5 +423,5 @@
 c   Initialisation de la physique :
 c   -------------------------------
-      IF (call_iniphys.and.iflag_phys.eq.1) THEN
+      IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100)) THEN
          latfi(1)=rlatu(1)
          lonfi(1)=0.
@@ -446,9 +446,10 @@
 ! Physics:
 #ifdef CPP_PHYS
-         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys ,
-     ,                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp     )
+         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys,
+     &                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp,
+     &                iflag_phys)
 #endif
          call_iniphys=.false.
-      ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
+      ENDIF ! of IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100))
 
 
@@ -481,12 +482,4 @@
  301  FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
  302  FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
-#endif
-
-#ifdef CPP_PHYS
-! Create start file (startphy.nc) and boundary conditions (limit.nc) 
-! for the Earth verstion
-       if (iflag_phys>=100) then
-          call iniaqua(ngridmx,latfi,lonfi,iflag_phys)
-       endif
 #endif
 
Index: /LMDZ5/trunk/libf/dyn3dpar/leapfrog_p.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3dpar/leapfrog_p.F	(revision 1670)
+++ /LMDZ5/trunk/libf/dyn3dpar/leapfrog_p.F	(revision 1671)
@@ -139,4 +139,5 @@
       REAL :: secondes
 
+      logical :: physic
       LOGICAL first,callinigrads
 
@@ -208,4 +209,6 @@
 
       itau = 0
+      physic=.true.
+      if (iflag_phys==0.or.iflag_phys==2) physic=.false.
 !      iday = day_ini+itau/day_step
 !      time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
@@ -364,5 +367,5 @@
      s        apdiss = .TRUE.
          IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward 
-     s          .and. iflag_phys.EQ.1                 ) apphys = .TRUE.
+     s          .and. physic                        ) apphys = .TRUE.
       ELSE
       ! Leapfrog/Matsuno time stepping 
@@ -370,5 +373,5 @@
          IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward )
      s        apdiss = .TRUE.
-         IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE.
+         IF( MOD(itau+1,iphysiq).EQ.0.AND.physic) apphys=.TRUE.
       END IF
 
Index: /LMDZ5/trunk/libf/phydev/comcstphy.F90
===================================================================
--- /LMDZ5/trunk/libf/phydev/comcstphy.F90	(revision 1671)
+++ /LMDZ5/trunk/libf/phydev/comcstphy.F90	(revision 1671)
@@ -0,0 +1,8 @@
+module comcstphy
+  
+  real :: rradius ! planet radius (m)
+  real :: rr      ! recuced gas constant: R/molar mass of atm
+  real :: rg      ! gravity
+  real :: rcpp    ! specific heat of the atmosphere
+
+end module comcstphy
Index: DZ5/trunk/libf/phydev/comcstphy.h
===================================================================
--- /LMDZ5/trunk/libf/phydev/comcstphy.h	(revision 1670)
+++ 	(revision )
@@ -1,2 +1,0 @@
-      real rradius,rr,rg,rcpp
-      common/comcstphy/rradius,rr,rg,rcpp
Index: /LMDZ5/trunk/libf/phydev/iniphysiq.F
===================================================================
--- /LMDZ5/trunk/libf/phydev/iniphysiq.F	(revision 1670)
+++ /LMDZ5/trunk/libf/phydev/iniphysiq.F	(revision 1671)
@@ -2,15 +2,15 @@
 ! $Id: iniphysiq.F 1403 2010-07-01 09:02:53Z fairhead $
 !
-c
-c
       SUBROUTINE iniphysiq(ngrid,nlayer,
      $           punjours,
      $           pdayref,ptimestep,
      $           plat,plon,parea,pcu,pcv,
-     $           prad,pg,pr,pcpp)
-      USE dimphy
-      USE mod_grid_phy_lmdz
-      USE mod_phys_lmdz_para
-      USE comgeomphy
+     $           prad,pg,pr,pcpp,iflag_phys)
+      USE dimphy, only : klev
+      USE mod_grid_phy_lmdz, only : klon_glo
+      USE mod_phys_lmdz_para, only : klon_omp,klon_omp_begin,
+     &                               klon_omp_end,klon_mpi_begin
+      USE comgeomphy, only : airephy,cuphy,cvphy,rlond,rlatd
+      USE comcstphy, only : rradius,rg,rr,rcpp
 
       IMPLICIT NONE
@@ -18,18 +18,7 @@
 c=======================================================================
 c
-c   subject:
-c   --------
+c   Initialisation of the physical constants and some positional and 
+c   geometrical arrays for the physics
 c
-c   Initialisation for the physical parametrisations of the LMD 
-c   martian atmospheric general circulation modele.
-c
-c   author: Frederic Hourdin 15 / 10 /93
-c   -------
-c
-c   arguments:
-c   ----------
-c
-c   input:
-c   ------
 c
 c    ngrid                 Size of the horizontal grid.
@@ -37,36 +26,38 @@
 c    nlayer                Number of vertical layers.
 c    pdayref               Day of reference for the simulation
-c    firstcall             True at the first call
-c    lastcall              True at the last call
-c    pday                  Number of days counted from the North. Spring
-c                          equinoxe.
 c
 c=======================================================================
-c
-c-----------------------------------------------------------------------
-c   declarations:
-c   -------------
+ 
  
 cym#include "dimensions.h"
 cym#include "dimphy.h"
 cym#include "comgeomphy.h"
-#include "comcstphy.h"
-      REAL prad,pg,pr,pcpp,punjours
- 
-      INTEGER ngrid,nlayer
-      REAL plat(ngrid),plon(ngrid),parea(klon_glo)
-      REAL pcu(klon_glo),pcv(klon_glo)
-      INTEGER pdayref
+#include "iniprint.h"
+
+      REAL,INTENT(IN) :: prad ! radius of the planet (m)
+      REAL,INTENT(IN) :: pg ! gravitational acceleration (m/s2)
+      REAL,INTENT(IN) :: pr ! ! reduced gas constant R/mu
+      REAL,INTENT(IN) :: pcpp ! specific heat Cp
+      REAL,INTENT(IN) :: punjours ! length (in s) of a standard day
+      INTEGER,INTENT(IN) :: ngrid ! number of horizontal grid points in the physics
+      INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
+      REAL,INTENT(IN) :: plat(ngrid) ! latitudes of the physics grid
+      REAL,INTENT(IN) :: plon(ngrid) ! longitudes of the physics grid
+      REAL,INTENT(IN) :: parea(klon_glo) ! area (m2)
+      REAL,INTENT(IN) :: pcu(klon_glo) ! cu coeff. (u_covariant = cu * u)
+      REAL,INTENT(IN) :: pcv(klon_glo) ! cv coeff. (v_covariant = cv * v)
+      INTEGER,INTENT(IN) :: pdayref ! reference day of for the simulation
+      REAL,INTENT(IN) :: ptimestep !physics time step (s)
+      INTEGER,INTENT(IN) :: iflag_phys ! type of physics to be called
+
       INTEGER :: ibegin,iend,offset
- 
-      REAL ptimestep
       CHARACTER (LEN=20) :: modname='iniphysiq'
       CHARACTER (LEN=80) :: abort_message
  
       IF (nlayer.NE.klev) THEN
-         PRINT*,'STOP in inifis'
-         PRINT*,'Probleme de dimensions :'
-         PRINT*,'nlayer     = ',nlayer
-         PRINT*,'klev   = ',klev
+         write(lunout,*) 'STOP in ',trim(modname)
+         write(lunout,*) 'Problem with dimensions :'
+         write(lunout,*) 'nlayer     = ',nlayer
+         write(lunout,*) 'klev   = ',klev
          abort_message = ''
          CALL abort_gcm (modname,abort_message,1)
@@ -74,13 +65,14 @@
 
       IF (ngrid.NE.klon_glo) THEN
-         PRINT*,'STOP in inifis'
-         PRINT*,'Probleme de dimensions :'
-         PRINT*,'ngrid     = ',ngrid
-         PRINT*,'klon   = ',klon_glo
+         write(lunout,*) 'STOP in ',trim(modname)
+         write(lunout,*) 'Problem with dimensions :'
+         write(lunout,*) 'ngrid     = ',ngrid
+         write(lunout,*) 'klon   = ',klon_glo
          abort_message = ''
          CALL abort_gcm (modname,abort_message,1)
       ENDIF
-c$OMP PARALLEL PRIVATE(ibegin,iend) 
-c$OMP+         SHARED(parea,pcu,pcv,plon,plat)
+
+!$OMP PARALLEL PRIVATE(ibegin,iend) 
+!$OMP+         SHARED(parea,pcu,pcv,plon,plat)
       
       offset=klon_mpi_begin-1
@@ -92,6 +84,5 @@
       rlatd(1:klon_omp)=plat(offset+klon_omp_begin:offset+klon_omp_end)
 
-!     call suphel
-!     prad,pg,pr,pcpp
+! copy some fundamental parameters to physics
       rradius=prad
       rg=pg
@@ -99,16 +90,21 @@
       rcpp=pcpp
 
-      
+!$OMP END PARALLEL
 
-c$OMP END PARALLEL
+!      print*,'ATTENTION !!! TRAVAILLER SUR INIPHYSIQ'
+!      print*,'CONTROLE DES LATITUDES, LONGITUDES, PARAMETRES ...'
 
-      print*,'ATTENTION !!! TRAVAILLER SUR INIPHYSIQ'
-      print*,'CONTROLE DES LATITUDES, LONGITUDES, PARAMETRES ...'
+! Additional initializations for aquaplanets 
+!$OMP PARALLEL
+      if (iflag_phys>=100) then
+        call iniaqua(klon_omp,rlatd,rlond,iflag_phys)
+      endif
+!$OMP END PARALLEL
 
-      RETURN
-9999  CONTINUE
-      abort_message ='Cette version demande les fichier rnatur.dat
-     & et surf.def'
-      CALL abort_gcm (modname,abort_message,1)
+!      RETURN
+!9999  CONTINUE
+!      abort_message ='Cette version demande les fichier rnatur.dat
+!     & et surf.def'
+!      CALL abort_gcm (modname,abort_message,1)
 
       END
Index: /LMDZ5/trunk/libf/phydev/iostart.F90
===================================================================
--- /LMDZ5/trunk/libf/phydev/iostart.F90	(revision 1671)
+++ /LMDZ5/trunk/libf/phydev/iostart.F90	(revision 1671)
@@ -0,0 +1,1 @@
+link ../phylmd/iostart.F90
Index: /LMDZ5/trunk/libf/phydev/phyaqua.F
===================================================================
--- /LMDZ5/trunk/libf/phydev/phyaqua.F	(revision 1670)
+++ /LMDZ5/trunk/libf/phydev/phyaqua.F	(revision 1671)
@@ -1,26 +1,41 @@
-! Routines complementaires pour la physique planetaire.
-
+!
+! $Id: $
+!
 
       subroutine iniaqua(nlon,latfi,lonfi,iflag_phys)
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!  Creation d'un etat initial et de conditions aux limites 
-!  (resp startphy.nc et limit.nc) pour des configurations idealisees 
-! du modele LMDZ dans sa version terrestre.
-!  iflag_phys est un parametre qui controle
-!  iflag_phys = N  
-!    de 100 a 199 : aqua planetes avec SST forcees
-!                 N-100 determine le type de SSTs
-!    de 200 a 299 : terra planetes avec Ts calcule
-!        
+!  Create an initial state (startphy.nc) for the physics
+!  Usefull for idealised cases (e.g. aquaplanets or testcases)
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
+      use phys_state_var_mod, only : rlat,rlon,
+     &                               phys_state_var_init
+      use mod_phys_lmdz_para, only : klon_omp
+      use comgeomphy, only : rlond,rlatd
+      implicit none
+      
+      integer,intent(in) :: nlon,iflag_phys
+      real,intent(in) :: lonfi(nlon),latfi(nlon)
 
-      integer nlon,iflag_phys
-cIM ajout latfi, lonfi
-      REAL, DIMENSION (nlon) :: lonfi, latfi
+! local variables
+      real :: pi
+
+! initializations:
+      pi=2.*asin(1.)
+
+      call phys_state_var_init()
+
+      rlat(1:klon_omp)=rlatd(1:klon_omp)*180./pi
+      rlon(1:klon_omp)=rlond(1:klon_omp)*180./pi
 
 
-      return
+! Here you could create an initial condition for the physics
+! ...
+! ... fill in the fields...
+! ...
+! ... and create a "startphy.nc" file
+!      CALL phyredem ("startphy.nc")
+
       end
 
Index: /LMDZ5/trunk/libf/phydev/phyredem.F
===================================================================
--- /LMDZ5/trunk/libf/phydev/phyredem.F	(revision 1671)
+++ /LMDZ5/trunk/libf/phydev/phyredem.F	(revision 1671)
@@ -0,0 +1,40 @@
+!
+! $Id: $
+!
+      SUBROUTINE phyredem (fichnom)
+
+      use iostart, only : open_restartphy,close_restartphy,
+     &                    put_var,put_field
+      use phys_state_var_mod, only : rlon,rlat
+
+      implicit none
+
+      character(len=*),intent(in) :: fichnom
+
+      integer,parameter :: tab_cntrl_len=100
+      real :: tab_cntrl(tab_cntrl_len)
+
+! open file
+
+      CALL open_restartphy(fichnom)
+
+! tab_cntrl() contains run parameters
+
+      tab_cntrl(:)=0.0
+ 
+
+      CALL put_var("controle","Parametres de controle",tab_cntrl)
+
+! coordinates
+
+      CALL put_field("longitude",
+     .               "Longitudes de la grille physique",rlon)
+     
+      CALL put_field("latitude","Latitudes de la grille physique",rlat)
+
+! close file
+
+      CALL close_restartphy
+!$OMP BARRIER
+
+      END
Index: /LMDZ5/trunk/libf/phydev/phys_state_var_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phydev/phys_state_var_mod.F90	(revision 1671)
+++ /LMDZ5/trunk/libf/phydev/phys_state_var_mod.F90	(revision 1671)
@@ -0,0 +1,33 @@
+!
+! $Id:  $
+!
+MODULE phys_state_var_mod
+!======================================================================
+! Variables saved in startphy.nc
+!======================================================================
+
+USE dimphy, only : klon
+ 
+
+REAL, ALLOCATABLE, SAVE :: rlat(:), rlon(:)
+!$OMP THREADPRIVATE(rlat,rlon)
+
+CONTAINS
+
+!======================================================================
+  SUBROUTINE phys_state_var_init()
+  use dimphy, only : klon
+
+  ALLOCATE(rlat(klon),rlon(klon))
+
+  END SUBROUTINE phys_state_var_init
+
+!======================================================================
+  SUBROUTINE phys_state_var_end
+  use dimphy, only : klon
+
+  deallocate(rlat,rlon)
+
+  END SUBROUTINE phys_state_var_end
+
+END MODULE phys_state_var_mod
Index: /LMDZ5/trunk/libf/phydev/physiq.F90
===================================================================
--- /LMDZ5/trunk/libf/phydev/physiq.F90	(revision 1670)
+++ /LMDZ5/trunk/libf/phydev/physiq.F90	(revision 1671)
@@ -11,7 +11,8 @@
      &            , PVteta)
 
-      USE dimphy
-      USE infotrac
-      USE comgeomphy
+      USE dimphy, only : klon,klev
+      USE infotrac, only : nqtot
+      USE comgeomphy, only : rlatd
+      USE comcstphy, only : rg
 
       IMPLICIT none
@@ -50,5 +51,5 @@
 !======================================================================
 #include "dimensions.h"
-#include "comcstphy.h"
+!#include "comcstphy.h"
 
       integer jjmp1
@@ -101,24 +102,41 @@
     PARAMETER    ( longcles = 20 )
 
-real temp_newton(klon,klev)
+real :: temp_newton(klon,klev)
 integer k
 logical, save :: first=.true.
+!$OMP THREADPRIVATE(first)
 
       REAL clesphy0( longcles      )
 
-d_u=0.
-d_v=0.
-d_t=0.
-d_qx=0.
-d_ps=0.
+! initializations
+if (first) then
+! ...
 
-     d_u(:,1)=-u(:,1)/86400.
-     do k=1,klev
-        temp_newton(:,k)=280.+cos(rlatd(:))*40.-pphi(:,k)/rg*6.e-3
-        d_t(:,k)=(temp_newton(:,k)-t(:,k))/1.e5
-     enddo
+  first=.false.
+endif
+
+! set all tendencies to zero
+d_u(:,:)=0.
+d_v(:,:)=0.
+d_t(:,:)=0.
+d_qx(:,:,:)=0.
+d_ps(:)=0.
+
+! compute tendencies to return to the dynamics:
+! "friction" on the first layer
+d_u(:,1)=-u(:,1)/86400.
+! newtonian rlaxation towards temp_newton()
+do k=1,klev
+  temp_newton(:,k)=280.+cos(rlatd(:))*40.-pphi(:,k)/rg*6.e-3
+  d_t(:,k)=(temp_newton(:,k)-t(:,k))/1.e5
+enddo
 
 
-      print*,'COUCOU PHYDEV'
-      return
-      end
+print*,'COUCOU PHYDEV'
+
+! if lastcall, then it is time to write "restartphy.nc" file
+if (lafin) then
+  call phyredem("restartphy.nc")
+endif
+
+end
Index: /LMDZ5/trunk/libf/phylmd/iniphysiq.F
===================================================================
--- /LMDZ5/trunk/libf/phylmd/iniphysiq.F	(revision 1670)
+++ /LMDZ5/trunk/libf/phylmd/iniphysiq.F	(revision 1671)
@@ -8,9 +8,10 @@
      $           pdayref,ptimestep,
      $           plat,plon,parea,pcu,pcv,
-     $           prad,pg,pr,pcpp)
-      USE dimphy
-      USE mod_grid_phy_lmdz
-      USE mod_phys_lmdz_para
-      USE comgeomphy
+     $           prad,pg,pr,pcpp,iflag_phys)
+      USE dimphy, only : klev
+      USE mod_grid_phy_lmdz, only : klon_glo
+      USE mod_phys_lmdz_para, only : klon_omp,klon_omp_begin,
+     &                               klon_omp_end,klon_mpi_begin
+      USE comgeomphy, only : airephy,cuphy,cvphy,rlond,rlatd
 
       IMPLICIT NONE
@@ -18,18 +19,7 @@
 c=======================================================================
 c
-c   subject:
-c   --------
+c   Initialisation of the physical constants and some positional and 
+c   geometrical arrays for the physics
 c
-c   Initialisation for the physical parametrisations of the LMD 
-c   martian atmospheric general circulation modele.
-c
-c   author: Frederic Hourdin 15 / 10 /93
-c   -------
-c
-c   arguments:
-c   ----------
-c
-c   input:
-c   ------
 c
 c    ngrid                 Size of the horizontal grid.
@@ -37,14 +27,6 @@
 c    nlayer                Number of vertical layers.
 c    pdayref               Day of reference for the simulation
-c    firstcall             True at the first call
-c    lastcall              True at the last call
-c    pday                  Number of days counted from the North. Spring
-c                          equinoxe.
 c
 c=======================================================================
-c
-c-----------------------------------------------------------------------
-c   declarations:
-c   -------------
  
 cym#include "dimensions.h"
@@ -52,21 +34,31 @@
 cym#include "comgeomphy.h"
 #include "YOMCST.h"
-      REAL prad,pg,pr,pcpp,punjours
- 
-      INTEGER ngrid,nlayer
-      REAL plat(ngrid),plon(ngrid),parea(klon_glo)
-      REAL pcu(klon_glo),pcv(klon_glo)
-      INTEGER pdayref
-      INTEGER :: ibegin,iend,offset
- 
-      REAL ptimestep
+#include "iniprint.h"
+
+      REAL,INTENT(IN) :: prad ! radius of the planet (m)
+      REAL,INTENT(IN) :: pg ! gravitational acceleration (m/s2)
+      REAL,INTENT(IN) :: pr ! ! reduced gas constant R/mu
+      REAL,INTENT(IN) :: pcpp ! specific heat Cp
+      REAL,INTENT(IN) :: punjours ! length (in s) of a standard day
+      INTEGER,INTENT(IN) :: ngrid ! number of horizontal grid points in the physics
+      INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
+      REAL,INTENT(IN) :: plat(ngrid) ! latitudes of the physics grid
+      REAL,INTENT(IN) :: plon(ngrid) ! longitudes of the physics grid
+      REAL,INTENT(IN) :: parea(klon_glo) ! area (m2)
+      REAL,INTENT(IN) :: pcu(klon_glo) ! cu coeff. (u_covariant = cu * u)
+      REAL,INTENT(IN) :: pcv(klon_glo) ! cv coeff. (v_covariant = cv * v)
+      INTEGER,INTENT(IN) :: pdayref ! reference day of for the simulation
+      REAL,INTENT(IN) :: ptimestep !physics time step (s)
+      INTEGER,INTENT(IN) :: iflag_phys ! type of physics to be called
+
+      INTEGER :: ibegin,iend,offset 
       CHARACTER (LEN=20) :: modname='iniphysiq'
       CHARACTER (LEN=80) :: abort_message
  
       IF (nlayer.NE.klev) THEN
-         PRINT*,'STOP in inifis'
-         PRINT*,'Probleme de dimensions :'
-         PRINT*,'nlayer     = ',nlayer
-         PRINT*,'klev   = ',klev
+         write(lunout,*) 'STOP in ',trim(modname)
+         write(lunout,*) 'Problem with dimensions :'
+         write(lunout,*) 'nlayer     = ',nlayer
+         write(lunout,*) 'klev   = ',klev
          abort_message = ''
          CALL abort_gcm (modname,abort_message,1)
@@ -74,13 +66,14 @@
 
       IF (ngrid.NE.klon_glo) THEN
-         PRINT*,'STOP in inifis'
-         PRINT*,'Probleme de dimensions :'
-         PRINT*,'ngrid     = ',ngrid
-         PRINT*,'klon   = ',klon_glo
+         write(lunout,*) 'STOP in ',trim(modname)
+         write(lunout,*) 'Problem with dimensions :'
+         write(lunout,*) 'ngrid     = ',ngrid
+         write(lunout,*) 'klon   = ',klon_glo
          abort_message = ''
          CALL abort_gcm (modname,abort_message,1)
       ENDIF
-c$OMP PARALLEL PRIVATE(ibegin,iend) 
-c$OMP+         SHARED(parea,pcu,pcv,plon,plat)
+
+!$OMP PARALLEL PRIVATE(ibegin,iend) 
+!$OMP+         SHARED(parea,pcu,pcv,plon,plat)
       
       offset=klon_mpi_begin-1
@@ -92,16 +85,76 @@
       rlatd(1:klon_omp)=plat(offset+klon_omp_begin:offset+klon_omp_end)
 
+      ! suphel => initialize some physical constants (orbital parameters,
+      !           geoid, gravity, thermodynamical constants, etc.) in the
+      !           physics
       call suphel
+      
+!$OMP END PARALLEL
 
-c$OMP END PARALLEL
+      ! check that physical constants set in 'suphel' are coherent
+      ! with values set in the dynamics:
+      if (RDAY.ne.punjours) then
+        write(lunout,*) "iniphysiq: length of day discrepancy!!!"
+        write(lunout,*) "  in the dynamics punjours=",punjours
+        write(lunout,*) "   but in the physics RDAY=",RDAY
+        if (abs(RDAY-punjours).gt.0.01) then
+          ! stop here if the relative difference is more than 1%
+          abort_message = 'length of day discrepancy'
+          CALL abort_gcm (modname,abort_message,1)
+        endif
+      endif
+      if (RG.ne.pg) then
+        write(lunout,*) "iniphysiq: gravity discrepancy !!!"
+        write(lunout,*) "     in the dynamics pg=",pg
+        write(lunout,*) "  but in the physics RG=",RG
+        if (abs(RG-pg).gt.0.01) then
+          ! stop here if the relative difference is more than 1%
+          abort_message = 'gravity discrepancy'
+          CALL abort_gcm (modname,abort_message,1)
+        endif
+      endif
+      if (RA.ne.prad) then
+        write(lunout,*) "iniphysiq: planet radius discrepancy !!!"
+        write(lunout,*) "   in the dynamics prad=",prad
+        write(lunout,*) "  but in the physics RA=",RA
+        if (abs(RA-prad).gt.0.01) then
+          ! stop here if the relative difference is more than 1%
+          abort_message = 'planet radius discrepancy'
+          CALL abort_gcm (modname,abort_message,1)
+        endif
+      endif
+      if (RD.ne.pr) then
+        write(lunout,*)"iniphysiq: reduced gas constant discrepancy !!!"
+        write(lunout,*)"     in the dynamics pr=",pr
+        write(lunout,*)"  but in the physics RD=",RD
+        if (abs(RD-pr).gt.0.01) then
+          ! stop here if the relative difference is more than 1%
+          abort_message = 'reduced gas constant discrepancy'
+          CALL abort_gcm (modname,abort_message,1)
+        endif
+      endif
+      if (RCPD.ne.pcpp) then
+        write(lunout,*)"iniphysiq: specific heat discrepancy !!!"
+        write(lunout,*)"     in the dynamics pcpp=",pcpp
+        write(lunout,*)"  but in the physics RCPD=",RCPD
+        if (abs(RCPD-pcpp).gt.0.01) then
+          ! stop here if the relative difference is more than 1%
+          abort_message = 'specific heat discrepancy'
+          CALL abort_gcm (modname,abort_message,1)
+        endif
+      endif
 
-      print*,'ATTENTION !!! TRAVAILLER SUR INIPHYSIQ'
-      print*,'CONTROLE DES LATITUDES, LONGITUDES, PARAMETRES ...'
+! Additional initializations for aquaplanets
+!$OMP PARALLEL
+      if (iflag_phys>=100) then
+        call iniaqua(klon_omp,rlatd,rlond,iflag_phys)
+      endif
+!$OMP END PARALLEL
 
-      RETURN
-9999  CONTINUE
-      abort_message ='Cette version demande les fichier rnatur.dat
-     & et surf.def'
-      CALL abort_gcm (modname,abort_message,1)
+!      RETURN
+!9999  CONTINUE
+!      abort_message ='Cette version demande les fichier rnatur.dat
+!     & et surf.def'
+!      CALL abort_gcm (modname,abort_message,1)
 
       END
Index: /LMDZ5/trunk/libf/phylmd/phyaqua.F
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phyaqua.F	(revision 1670)
+++ /LMDZ5/trunk/libf/phylmd/phyaqua.F	(revision 1671)
@@ -16,12 +16,11 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-      use comgeomphy
-      use dimphy
+      use comgeomphy, only : rlatd,rlond
+      use dimphy, only : klon
       use surface_data, only : type_ocean,ok_veget
       use pbl_surface_mod, only : pbl_surface_init
       USE fonte_neige_mod, only : fonte_neige_init
       use phys_state_var_mod
-      use control_mod
-
+      use control_mod, only : dayref,nday,iphysiq
 
       USE IOIPSL 
@@ -35,8 +34,10 @@
 #include "dimsoil.h"
 #include "indicesol.h"
-
-      integer nlon,iflag_phys
+#include "temps.h"
+
+      integer,intent(in) :: nlon,iflag_phys
 cIM ajout latfi, lonfi
-      REAL, DIMENSION (nlon) :: lonfi, latfi
+      real,intent(in) :: lonfi(nlon),latfi(nlon)
+
       INTEGER type_profil,type_aqua
 
@@ -71,6 +72,4 @@
 !      integer demih_pas
 
-      integer day_ini
-
       CHARACTER*80 ans,file_forctl, file_fordat, file_start
       character*100 file,var
@@ -88,5 +87,5 @@
       REAL phy_flic(nlon,360)
 
-      integer, save::  read_climoz ! read ozone climatology
+      integer, save::  read_climoz=0 ! read ozone climatology
 
 
@@ -131,7 +130,10 @@
       type_aqua=iflag_phys/100
       type_profil=iflag_phys-type_aqua*100
-      print*,'type_aqua, type_profil',type_aqua, type_profil
-
-      if (klon.ne.nlon) stop'probleme de dimensions dans iniaqua'
+      print*,'iniaqua:type_aqua, type_profil',type_aqua, type_profil
+
+      if (klon.ne.nlon) then
+        write(*,*)"iniaqua: klon=",klon," nlon=",nlon
+        stop'probleme de dimensions dans iniaqua'
+      endif
       call phys_state_var_init(read_climoz)
 
@@ -154,4 +156,5 @@
 
          day_ini=dayref
+         day_end=day_ini+nday
          airefi=1.
          zcufi=1.
@@ -171,5 +174,5 @@
       radsol=0.
       qsol_f=10.
-      CALL getin('albedo',albedo)
+!      CALL getin('albedo',albedo) ! albedo is set below, depending on type_aqua
       alb_ocean=.true.
       CALL getin('alb_ocean',alb_ocean)
@@ -180,4 +183,7 @@
       qsol(:)    = qsol_f
       rugsrel = 0.0    ! (rugsrel = rugoro)
+      rugoro = 0.0
+      u_ancien = 0.0
+      v_ancien = 0.0
       agesno  = 50.0
 ! Relief plat
@@ -308,5 +314,5 @@
      .     evap, frugs, agesno, tsoil)
 
-        print*,'avant phyredem dans iniaqua'
+        print*,'iniaqua: before phyredem'
 
       falb1=albedo
@@ -329,5 +335,5 @@
       CALL phyredem ("startphy.nc")
 
-        print*,'apres phyredem'
+        print*,'iniaqua: after phyredem'
       call phys_state_var_end
 
@@ -450,24 +456,32 @@
       RETURN
       END
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
       subroutine writelim
      s   (klon,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice,
      s    phy_fter,phy_foce,phy_flic,phy_fsic)
 c
+      use mod_phys_lmdz_para, only: is_mpi_root,is_omp_root
+      use mod_grid_phy_lmdz, only : klon_glo
+      use mod_phys_lmdz_transfert_para, only : gather
 !#include "dimensions.h"
 !#include "dimphy.h"
 #include "netcdf.inc"
  
-      integer klon
-      REAL phy_nat(klon,360)
-      REAL phy_alb(klon,360)
-      REAL phy_sst(klon,360)
-      REAL phy_bil(klon,360)
-      REAL phy_rug(klon,360)
-      REAL phy_ice(klon,360)
-      REAL phy_fter(klon,360)
-      REAL phy_foce(klon,360)
-      REAL phy_flic(klon,360)
-      REAL phy_fsic(klon,360)
- 
+      integer,intent(in) :: klon
+      real,intent(in) :: phy_nat(klon,360)
+      real,intent(in) :: phy_alb(klon,360)
+      real,intent(in) :: phy_sst(klon,360)
+      real,intent(in) :: phy_bil(klon,360)
+      real,intent(in) :: phy_rug(klon,360)
+      real,intent(in) :: phy_ice(klon,360)
+      real,intent(in) :: phy_fter(klon,360)
+      real,intent(in) :: phy_foce(klon,360)
+      real,intent(in) :: phy_flic(klon,360)
+      real,intent(in) :: phy_fsic(klon,360)
+
+      real :: phy_glo(klon_glo,360) ! temporary variable, to store phy_***(:)
+                                    ! on the whole physics grid 
       INTEGER ierr
       INTEGER dimfirst(3)
@@ -480,93 +494,194 @@
       INTEGER id_FTER,id_FOCE,id_FSIC,id_FLIC
  
-      PRINT*, 'Ecriture du fichier limit'
-c
-      ierr = NF_CREATE ("limit.nc", NF_CLOBBER, nid)
-c
-      ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 30,
+      if (is_mpi_root.and.is_omp_root) then
+      
+        PRINT*, 'writelim: Ecriture du fichier limit'
+c
+        ierr = NF_CREATE ("limit.nc", NF_CLOBBER, nid)
+c
+        ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 30,
      .                       "Fichier conditions aux limites")
-      ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
-      ierr = NF_DEF_DIM (nid, "time", NF_UNLIMITED, ntim)
-c
-      dims(1) = ndim
-      dims(2) = ntim
+!!        ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
+        ierr = NF_DEF_DIM (nid, "points_physiques", klon_glo, ndim)
+        ierr = NF_DEF_DIM (nid, "time", NF_UNLIMITED, ntim)
+c
+        dims(1) = ndim
+        dims(2) = ntim
 c
 ccc      ierr = NF_DEF_VAR (nid, "TEMPS", NF_DOUBLE, 1,ntim, id_tim)
-      ierr = NF_DEF_VAR (nid, "TEMPS", NF_FLOAT, 1,ntim, id_tim)
-      ierr = NF_PUT_ATT_TEXT (nid, id_tim, "title", 17,
+        ierr = NF_DEF_VAR (nid, "TEMPS", NF_FLOAT, 1,ntim, id_tim)
+        ierr = NF_PUT_ATT_TEXT (nid, id_tim, "title", 17,
      .                        "Jour dans l annee")
 ccc      ierr = NF_DEF_VAR (nid, "NAT", NF_DOUBLE, 2,dims, id_NAT)
-      ierr = NF_DEF_VAR (nid, "NAT", NF_FLOAT, 2,dims, id_NAT)
-      ierr = NF_PUT_ATT_TEXT (nid, id_NAT, "title", 23,
+        ierr = NF_DEF_VAR (nid, "NAT", NF_FLOAT, 2,dims, id_NAT)
+        ierr = NF_PUT_ATT_TEXT (nid, id_NAT, "title", 23,
      .                        "Nature du sol (0,1,2,3)")
 ccc      ierr = NF_DEF_VAR (nid, "SST", NF_DOUBLE, 2,dims, id_SST)
-      ierr = NF_DEF_VAR (nid, "SST", NF_FLOAT, 2,dims, id_SST)
-      ierr = NF_PUT_ATT_TEXT (nid, id_SST, "title", 35,
+        ierr = NF_DEF_VAR (nid, "SST", NF_FLOAT, 2,dims, id_SST)
+        ierr = NF_PUT_ATT_TEXT (nid, id_SST, "title", 35,
      .                        "Temperature superficielle de la mer")
 ccc      ierr = NF_DEF_VAR (nid, "BILS", NF_DOUBLE, 2,dims, id_BILS)
-      ierr = NF_DEF_VAR (nid, "BILS", NF_FLOAT, 2,dims, id_BILS)
-      ierr = NF_PUT_ATT_TEXT (nid, id_BILS, "title", 32,
+        ierr = NF_DEF_VAR (nid, "BILS", NF_FLOAT, 2,dims, id_BILS)
+        ierr = NF_PUT_ATT_TEXT (nid, id_BILS, "title", 32,
      .                        "Reference flux de chaleur au sol")
 ccc      ierr = NF_DEF_VAR (nid, "ALB", NF_DOUBLE, 2,dims, id_ALB)
-      ierr = NF_DEF_VAR (nid, "ALB", NF_FLOAT, 2,dims, id_ALB)
-      ierr = NF_PUT_ATT_TEXT (nid, id_ALB, "title", 19,
+        ierr = NF_DEF_VAR (nid, "ALB", NF_FLOAT, 2,dims, id_ALB)
+        ierr = NF_PUT_ATT_TEXT (nid, id_ALB, "title", 19,
      .                        "Albedo a la surface")
 ccc      ierr = NF_DEF_VAR (nid, "RUG", NF_DOUBLE, 2,dims, id_RUG)
-      ierr = NF_DEF_VAR (nid, "RUG", NF_FLOAT, 2,dims, id_RUG)
-      ierr = NF_PUT_ATT_TEXT (nid, id_RUG, "title", 8,
+        ierr = NF_DEF_VAR (nid, "RUG", NF_FLOAT, 2,dims, id_RUG)
+        ierr = NF_PUT_ATT_TEXT (nid, id_RUG, "title", 8,
      .                        "Rugosite")
 
-      ierr = NF_DEF_VAR (nid, "FTER", NF_FLOAT, 2,dims, id_FTER)
-      ierr = NF_PUT_ATT_TEXT (nid, id_FTER, "title", 8,"Frac. Terre")
-      ierr = NF_DEF_VAR (nid, "FOCE", NF_FLOAT, 2,dims, id_FOCE)
-      ierr = NF_PUT_ATT_TEXT (nid, id_FOCE, "title", 8,"Frac. Terre")
-      ierr = NF_DEF_VAR (nid, "FSIC", NF_FLOAT, 2,dims, id_FSIC)
-      ierr = NF_PUT_ATT_TEXT (nid, id_FSIC, "title", 8,"Frac. Terre")
-      ierr = NF_DEF_VAR (nid, "FLIC", NF_FLOAT, 2,dims, id_FLIC)
-      ierr = NF_PUT_ATT_TEXT (nid, id_FLIC, "title", 8,"Frac. Terre")
-c
-      ierr = NF_ENDDEF(nid)
-c
-      DO k = 1, 360
-c
-      debut(1) = 1
-      debut(2) = k
-      epais(1) = klon
-      epais(2) = 1
-c
-      print*,'Instant ',k
-#ifdef NC_DOUBLE
-      print*,'NC DOUBLE'
-      ierr = NF_PUT_VAR1_DOUBLE (nid,id_tim,k,DBLE(k))
-      ierr = NF_PUT_VARA_DOUBLE (nid,id_NAT,debut,epais,phy_nat(1,k))
-      ierr = NF_PUT_VARA_DOUBLE (nid,id_SST,debut,epais,phy_sst(1,k))
-      ierr = NF_PUT_VARA_DOUBLE (nid,id_BILS,debut,epais,phy_bil(1,k))
-      ierr = NF_PUT_VARA_DOUBLE (nid,id_ALB,debut,epais,phy_alb(1,k))
-      ierr = NF_PUT_VARA_DOUBLE (nid,id_RUG,debut,epais,phy_rug(1,k))
-      ierr = NF_PUT_VARA_DOUBLE (nid,id_FTER,debut,epais,phy_fter(1,k))
-      ierr = NF_PUT_VARA_DOUBLE (nid,id_FOCE,debut,epais,phy_foce(1,k))
-      ierr = NF_PUT_VARA_DOUBLE (nid,id_FSIC,debut,epais,phy_fsic(1,k))
-      ierr = NF_PUT_VARA_DOUBLE (nid,id_FLIC,debut,epais,phy_flic(1,k))
-#else
-      print*,'NC PAS DOUBLE'
-      ierr = NF_PUT_VAR1_REAL (nid,id_tim,k,FLOAT(k))
-      ierr = NF_PUT_VARA_REAL (nid,id_NAT,debut,epais,phy_nat(1,k))
-      ierr = NF_PUT_VARA_REAL (nid,id_SST,debut,epais,phy_sst(1,k))
-      ierr = NF_PUT_VARA_REAL (nid,id_BILS,debut,epais,phy_bil(1,k))
-      ierr = NF_PUT_VARA_REAL (nid,id_ALB,debut,epais,phy_alb(1,k))
-      ierr = NF_PUT_VARA_REAL (nid,id_RUG,debut,epais,phy_rug(1,k))
-      ierr = NF_PUT_VARA_REAL (nid,id_FTER,debut,epais,phy_fter(1,k))
-      ierr = NF_PUT_VARA_REAL (nid,id_FOCE,debut,epais,phy_foce(1,k))
-      ierr = NF_PUT_VARA_REAL (nid,id_FSIC,debut,epais,phy_fsic(1,k))
-      ierr = NF_PUT_VARA_REAL (nid,id_FLIC,debut,epais,phy_flic(1,k))
-
-#endif
-c
-      ENDDO
-c
-      ierr = NF_CLOSE(nid)
-c
-      return
+        ierr = NF_DEF_VAR (nid, "FTER", NF_FLOAT, 2,dims, id_FTER)
+        ierr = NF_PUT_ATT_TEXT (nid, id_FTER, "title", 8,"Frac. Terre")
+        ierr = NF_DEF_VAR (nid, "FOCE", NF_FLOAT, 2,dims, id_FOCE)
+        ierr = NF_PUT_ATT_TEXT (nid, id_FOCE, "title", 8,"Frac. Terre")
+        ierr = NF_DEF_VAR (nid, "FSIC", NF_FLOAT, 2,dims, id_FSIC)
+        ierr = NF_PUT_ATT_TEXT (nid, id_FSIC, "title", 8,"Frac. Terre")
+        ierr = NF_DEF_VAR (nid, "FLIC", NF_FLOAT, 2,dims, id_FLIC)
+        ierr = NF_PUT_ATT_TEXT (nid, id_FLIC, "title", 8,"Frac. Terre")
+c
+        ierr = NF_ENDDEF(nid)
+c
+
+! write the 'times'
+        do k=1,360
+#ifdef NC_DOUBLE
+          ierr = NF_PUT_VAR1_DOUBLE (nid,id_tim,k,DBLE(k))
+#else
+          ierr = NF_PUT_VAR1_REAL (nid,id_tim,k,FLOAT(k))
+#endif
+        enddo
+
+      endif ! of if (is_mpi_root.and.is_omp_root)
+
+! write the fields, after having collected them on master
+
+      call gather(phy_nat,phy_glo)
+      if (is_mpi_root.and.is_omp_root) then
+#ifdef NC_DOUBLE
+        ierr=NF_PUT_VAR_DOUBLE(nid,id_NAT,phy_glo)
+#else
+        ierr=NF_PUT_VAR_REAL(nid,id_NAT,phy_glo)
+#endif
+        if(ierr.ne.NF_NOERR) then
+          write(*,*) "writelim error with phy_nat"
+          write(*,*) NF_STRERROR(ierr)
+        endif
+      endif
+
+      call gather(phy_sst,phy_glo)
+      if (is_mpi_root.and.is_omp_root) then
+#ifdef NC_DOUBLE
+        ierr=NF_PUT_VAR_DOUBLE(nid,id_SST,phy_glo)
+#else
+        ierr=NF_PUT_VAR_REAL(nid,id_SST,phy_glo)
+#endif
+        if(ierr.ne.NF_NOERR) then
+          write(*,*) "writelim error with phy_sst"
+          write(*,*) NF_STRERROR(ierr)
+        endif
+      endif
+
+      call gather(phy_bil,phy_glo)
+      if (is_mpi_root.and.is_omp_root) then
+#ifdef NC_DOUBLE
+        ierr=NF_PUT_VAR_DOUBLE(nid,id_BILS,phy_glo)
+#else
+        ierr=NF_PUT_VAR_REAL(nid,id_BILS,phy_glo)
+#endif
+        if(ierr.ne.NF_NOERR) then
+          write(*,*) "writelim error with phy_bil"
+          write(*,*) NF_STRERROR(ierr)
+        endif
+      endif
+
+      call gather(phy_alb,phy_glo)
+      if (is_mpi_root.and.is_omp_root) then
+#ifdef NC_DOUBLE
+        ierr=NF_PUT_VAR_DOUBLE(nid,id_ALB,phy_glo)
+#else
+        ierr=NF_PUT_VAR_REAL(nid,id_ALB,phy_glo)
+#endif
+        if(ierr.ne.NF_NOERR) then
+          write(*,*) "writelim error with phy_alb"
+          write(*,*) NF_STRERROR(ierr)
+        endif
+      endif
+
+      call gather(phy_rug,phy_glo)
+      if (is_mpi_root.and.is_omp_root) then
+#ifdef NC_DOUBLE
+        ierr=NF_PUT_VAR_DOUBLE(nid,id_RUG,phy_glo)
+#else
+        ierr=NF_PUT_VAR_REAL(nid,id_RUG,phy_glo)
+#endif
+        if(ierr.ne.NF_NOERR) then
+          write(*,*) "writelim error with phy_rug"
+          write(*,*) NF_STRERROR(ierr)
+        endif
+      endif
+
+      call gather(phy_fter,phy_glo)
+      if (is_mpi_root.and.is_omp_root) then
+#ifdef NC_DOUBLE
+        ierr=NF_PUT_VAR_DOUBLE(nid,id_FTER,phy_glo)
+#else
+        ierr=NF_PUT_VAR_REAL(nid,id_FTER,phy_glo)
+#endif
+        if(ierr.ne.NF_NOERR) then
+          write(*,*) "writelim error with phy_fter"
+          write(*,*) NF_STRERROR(ierr)
+        endif
+      endif
+
+      call gather(phy_foce,phy_glo)
+      if (is_mpi_root.and.is_omp_root) then
+#ifdef NC_DOUBLE
+        ierr=NF_PUT_VAR_DOUBLE(nid,id_FOCE,phy_glo)
+#else
+        ierr=NF_PUT_VAR_REAL(nid,id_FOCE,phy_glo)
+#endif
+        if(ierr.ne.NF_NOERR) then
+          write(*,*) "writelim error with phy_foce"
+          write(*,*) NF_STRERROR(ierr)
+        endif
+      endif
+
+      call gather(phy_fsic,phy_glo)
+      if (is_mpi_root.and.is_omp_root) then
+#ifdef NC_DOUBLE
+        ierr=NF_PUT_VAR_DOUBLE(nid,id_FSIC,phy_glo)
+#else
+        ierr=NF_PUT_VAR_REAL(nid,id_FSIC,phy_glo)
+#endif
+        if(ierr.ne.NF_NOERR) then
+          write(*,*) "writelim error with phy_fsic"
+          write(*,*) NF_STRERROR(ierr)
+        endif
+      endif
+
+      call gather(phy_flic,phy_glo)
+      if (is_mpi_root.and.is_omp_root) then
+#ifdef NC_DOUBLE
+        ierr=NF_PUT_VAR_DOUBLE(nid,id_FLIC,phy_glo)
+#else
+        ierr=NF_PUT_VAR_REAL(nid,id_FLIC,phy_glo)
+#endif
+        if(ierr.ne.NF_NOERR) then
+          write(*,*) "writelim error with phy_flic"
+          write(*,*) NF_STRERROR(ierr)
+        endif
+      endif
+
+!  close file:
+      if (is_mpi_root.and.is_omp_root) then
+        ierr = NF_CLOSE(nid)
+      endif
+
       end
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
       SUBROUTINE profil_sst(nlon, rlatd, type_profil, phy_sst)
