Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/LIM_SRC_2/limmsh_2.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/LIM_SRC_2/limmsh_2.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/LIM_SRC_2/limmsh_2.F90	(revision 2007)
@@ -66,4 +66,7 @@
       ENDIF
       
+      IF( jphgr_msh == 2 .OR. jphgr_msh == 3 .OR. jphgr_msh == 5 )   &
+          &      CALL ctl_stop(' Coriolis parameter in LIM not set for f- or beta-plane' )
+
       !----------------------------------------------------------                          
       !    Initialization of local and some global (common) variables 
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 2007)
@@ -80,5 +80,4 @@
       INTEGER  ::   ifvt, i1mfr, idfr               ! some switches
       INTEGER  ::   iflt, ial, iadv, ifral, ifrdv
-      INTEGER  ::   ii0, ii1, ij0, ij1  ! temporary integers
       REAL(wp) ::   zrdtir           ! 1. / rdt_ice
       REAL(wp) ::   zqsr  , zqns     ! solar & non solar heat flux
@@ -109,14 +108,17 @@
          sice_r(:,:) = sice
          !
-         IF( cp_cfg == "orca"  .AND. jp_cfg == 2 ) THEN
-            !                                        ! =======================
-            !                                        !  ORCA_R2 configuration
-            !                                        ! =======================
-            ii0 = 145   ;   ii1 = 180        ! Baltic Sea
-            ij0 = 113   ;   ij1 = 130   ;   soce_r(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0
-                                            sice_r(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 2.e0
-         ENDIF
-         !
-      ENDIF
+         IF( cp_cfg == "orca" ) THEN
+           !   ocean/ice salinity in the Baltic sea 
+           DO jj = 1, jpj
+              DO ji = 1, jpi
+                 IF( glamt(ji,jj) >= 14. .AND.  glamt(ji,jj) <= 32. .AND. gphit(ji,jj) >= 54. .AND. gphit(ji,jj) <= 66. ) THEN 
+                   soce_r(ji,jj) = 4.e0 
+                   sice_r(ji,jj) = 2.e0
+                 END IF
+              END DO
+           END DO
+           !
+         END IF
+      END IF
 
       !------------------------------------------!
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90	(revision 2007)
@@ -360,5 +360,6 @@
          END DO
       ENDIF
-      
+
+      CALL lbc_lnk( frld , 'T', 1. )      
       
       ! Select points for lateral accretion (this occurs when heat exchange
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/LIM_SRC_2/limtrp_2.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/LIM_SRC_2/limtrp_2.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/LIM_SRC_2/limtrp_2.F90	(revision 2007)
@@ -4,4 +4,8 @@
    !! LIM 2.0 transport ice model : sea-ice advection/diffusion
    !!======================================================================
+   !! History :  LIM  !  2000-01 (UCL)  Original code
+   !!            2.0  !  2001-05 (G. Madec, R. Hordoir) opa norm
+   !!             -   !  2004-01 (G. Madec, C. Ethe)  F90, mpp
+   !!----------------------------------------------------------------------
 #if defined key_lim2
    !!----------------------------------------------------------------------
@@ -11,27 +15,23 @@
    !!   lim_trp_init_2 : initialization and namelist read
    !!----------------------------------------------------------------------
-   !! * Modules used
-   USE phycst
-   USE dom_oce
+   USE phycst          ! physical constant
+   USE sbc_oce         ! ocean surface boundary condition
+   USE dom_oce         ! ocean domain
    USE in_out_manager  ! I/O manager
-   USE dom_ice_2
-   USE ice_2
-   USE limistate_2
-   USE limadv_2
-   USE limhdf_2
-   USE lbclnk
-   USE lib_mpp
+   USE dom_ice_2       ! LIM-2 domain
+   USE ice_2           ! LIM-2 variables
+   USE limistate_2     ! LIM-2 initial state
+   USE limadv_2        ! LIM-2 advection
+   USE limhdf_2        ! LIM-2 horizontal diffusion
+   USE lbclnk          ! lateral boundary conditions -- MPP exchanges
+   USE lib_mpp         ! MPP library
 
    IMPLICIT NONE
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC lim_trp_2     ! called by sbc_ice_lim_2
-
-   !! * Shared module variables
-   REAL(wp), PUBLIC  ::   &  !:
-      bound  = 0.e0          !: boundary condit. (0.0 no-slip, 1.0 free-slip)
-
-   !! * Module variables
+   PUBLIC   lim_trp_2   ! called by sbc_ice_lim_2
+
+   REAL(wp), PUBLIC  ::   bound  = 0.e0   !: boundary condit. (0.0 no-slip, 1.0 free-slip)
+
    REAL(wp)  ::           &  ! constant values
       epsi06 = 1.e-06  ,  &
@@ -44,7 +44,7 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   !! NEMO/LIM 3.2,  UCL-LOCEAN-IPSL (2010) 
    !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -62,34 +62,16 @@
       !!
       !! ** action :
-      !!
-      !! History :
-      !!   1.0  !  00-01 (LIM)  Original code
-      !!        !  01-05 (G. Madec, R. Hordoir) opa norm
-      !!   2.0  !  04-01 (G. Madec, C. Ethe)  F90, mpp
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt     ! number of iteration
-
-      INTEGER  ::   ji, jj, jk,   &  ! dummy loop indices
-         &          initad           ! number of sub-timestep for the advection
-
-      REAL(wp) ::  &                              
-         zindb  ,  &
-         zacrith, &
-         zindsn , &
-         zindic , &
-         zusvosn, &
-         zusvoic, &
-         zignm  , &
-         zindhe , &
-         zvbord , &
-         zcfl   , &
-         zusnit , &
-         zrtt, ztsn, ztic1, ztic2
-
-      REAL(wp), DIMENSION(jpi,jpj)  ::   &  ! temporary workspace
-         zui_u , zvi_v , zsm   ,         &
-         zs0ice, zs0sn , zs0a  ,         &
-         zs0c0 , zs0c1 , zs0c2 ,         &
-         zs0st
+      !!
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      INTEGER  ::   initad       ! number of sub-timestep for the advection
+      REAL(wp) ::   zindb  , zindsn , zindic, zacrith   ! local scalars
+      REAL(wp) ::   zusvosn, zusvoic, zignm , zindhe    !   -      -
+      REAL(wp) ::   zvbord , zcfl   , zusnit            !   -      -
+      REAL(wp) ::   zrtt   , ztsn   , ztic1 , ztic2     !   -      -
+      REAL(wp), DIMENSION(jpi,jpj)  ::   zui_u , zvi_v , zsm             ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj)  ::   zs0ice, zs0sn , zs0a            !  -      -
+      REAL(wp), DIMENSION(jpi,jpj)  ::   zs0c0 , zs0c1 , zs0c2 , zs0st   !  -      -
       !---------------------------------------------------------------------
 
@@ -105,6 +87,5 @@
          ! ice velocities at ocean U- and V-points (zui_u,zvi_v)
          ! ---------------------------------------
-         ! zvbord factor between 1 and 2 to take into account slip or no-slip boundary conditions.        
-         zvbord = 1.0 + ( 1.0 - bound )
+         zvbord = 1.0 + ( 1.0 - bound )      ! zvbord=2 no-slip, =0 free slip boundary conditions        
          DO jj = 1, jpjm1
             DO ji = 1, jpim1   ! NO vector opt.
@@ -113,7 +94,6 @@
             END DO
          END DO
-         ! Lateral boundary conditions on zui_u, zvi_v
-         CALL lbc_lnk( zui_u, 'U', -1. )
-         CALL lbc_lnk( zvi_v, 'V', -1. )
+         CALL lbc_lnk( zui_u, 'U', -1. )   ;   CALL lbc_lnk( zvi_v, 'V', -1. )         ! Lateral boundary conditions
+
 
          ! CFL test for stability
@@ -122,28 +102,27 @@
          zcfl  = MAX( zcfl, MAXVAL( ABS( zui_u(1:jpim1, :     ) ) * rdt_ice / e1u(1:jpim1, :     ) ) )
          zcfl  = MAX( zcfl, MAXVAL( ABS( zvi_v( :     ,1:jpjm1) ) * rdt_ice / e2v( :     ,1:jpjm1) ) )
-
-         IF (lk_mpp ) CALL mpp_max(zcfl)
-
-         IF ( zcfl > 0.5 .AND. lwp )   WRITE(numout,*) 'lim_trp_2 : violation of cfl criterion the ',nday,'th day, cfl = ',zcfl
+         !
+         IF(lk_mpp)   CALL mpp_max( zcfl )
+         !
+         IF( zcfl > 0.5 .AND. lwp )   WRITE(numout,*) 'lim_trp_2 : violation of cfl criterion the ',nday,'th day, cfl = ', zcfl
 
          ! content of properties
          ! ---------------------
          zs0sn (:,:) =  hsnm(:,:) * area(:,:)                 ! Snow volume.
-         zs0ice(:,:) =  hicm (:,:) * area(:,:)                ! Ice volume.
-         zs0a  (:,:) =  ( 1.0 - frld(:,:) ) * area(:,:)       ! Surface covered by ice.
-         zs0c0 (:,:) =  tbif(:,:,1) / rt0_snow * zs0sn(:,:)   ! Heat content of the snow layer.
+         zs0ice(:,:) =  hicm(:,:) * area(:,:)                 ! Ice volume.
+         zs0a  (:,:) =  ( 1.0 - frld(:,:) )    * area  (:,:)  ! Surface covered by ice.
+         zs0c0 (:,:) =  tbif(:,:,1) / rt0_snow * zs0sn (:,:)  ! Heat content of the snow layer.
          zs0c1 (:,:) =  tbif(:,:,2) / rt0_ice  * zs0ice(:,:)  ! Heat content of the first ice layer.
          zs0c2 (:,:) =  tbif(:,:,3) / rt0_ice  * zs0ice(:,:)  ! Heat content of the second ice layer.
-         zs0st (:,:) =  qstoif(:,:) / xlic     * zs0a(:,:)    ! Heat reservoir for brine pockets.
+         zs0st (:,:) =  qstoif(:,:) / xlic     * zs0a  (:,:)  ! Heat reservoir for brine pockets.
          
  
-         ! Advection 
+         ! Advection (Prather scheme)
          ! ---------
-         ! If ice drift field is too fast, use an appropriate time step for advection.         
-         initad = 1 + INT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) )
-         zusnit = 1.0 / REAL( initad ) 
-         
-         IF ( MOD( nday , 2 ) == 0) THEN
-            DO jk = 1,initad
+         initad = 1 + INT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) )   ! If ice drift field is too fast,          
+         zusnit = 1.0 / REAL( initad )                              ! split the ice time step in two
+         !
+         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0) THEN        !==  odd ice time step:  adv_x then adv_y  ==!
+            DO jk = 1, initad
                CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice )
                CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice )
@@ -161,5 +140,5 @@
                CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0st , sxst , sxxst , syst , syyst , sxyst  )
             END DO
-         ELSE
+         ELSE                                                 !==  even ice time step:  adv_x then adv_y  ==!
             DO jk = 1, initad
                CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice )
@@ -182,4 +161,5 @@
          ! recover the properties from their contents
          ! ------------------------------------------
+!!gm Define in limmsh one for all area = 1 /area  (CPU time saved !)
          zs0ice(:,:) = zs0ice(:,:) / area(:,:)
          zs0sn (:,:) = zs0sn (:,:) / area(:,:)
@@ -205,4 +185,12 @@
             END DO
          END DO
+!!gm more readable coding: (and avoid an error in F90 with sign of zero)
+!        DO jj = 1, jpjm1          ! NB: has not to be defined on jpj line and jpi row
+!           DO ji = 1 , fs_jpim1   ! vector opt.
+!              IF( MIN( zs0a(ji,jj) , zs0a(ji+1,jj) ) == 0.e0 )   pahu(ji,jj) = 0.e0
+!              IF( MIN( zs0a(ji,jj) , zs0a(ji,jj+1) ) == 0.e0 )   pahv(ji,jj) = 0.e0
+!           END DO
+!        END DO
+!!gm end
 
          ! diffusion
@@ -216,8 +204,9 @@
          CALL lim_hdf_2( zs0st  )
 
-         zs0ice(:,:) = MAX( rzero, zs0ice(:,:) * area(:,:) )    !!bug:  est-ce utile
-         zs0sn (:,:) = MAX( rzero, zs0sn (:,:) * area(:,:) )    !!bug:  cf /area  juste apres
-         zs0a  (:,:) = MAX( rzero, zs0a  (:,:) * area(:,:) )    !! suppression des 2 change le resultat...
-         zs0c0 (:,:) = MAX( rzero, zs0c0 (:,:) * area(:,:) )
+!!gm see comment this can be skipped
+         zs0ice(:,:) = MAX( rzero, zs0ice(:,:) * area(:,:) )    !!bug:  useless
+         zs0sn (:,:) = MAX( rzero, zs0sn (:,:) * area(:,:) )    !!bug:  cf /area  just below
+         zs0a  (:,:) = MAX( rzero, zs0a  (:,:) * area(:,:) )    !! caution: the suppression of the 2 changes 
+         zs0c0 (:,:) = MAX( rzero, zs0c0 (:,:) * area(:,:) )    !! the last digit of the results
          zs0c1 (:,:) = MAX( rzero, zs0c1 (:,:) * area(:,:) )
          zs0c2 (:,:) = MAX( rzero, zs0c2 (:,:) * area(:,:) )
@@ -225,14 +214,11 @@
 
 
-         ! -------------------------------------------------------------------!
-         !   Up-dating and limitation of sea ice properties after transport   !
-         ! -------------------------------------------------------------------!
-
-         ! Up-dating and limitation of sea ice properties after transport.
+         !-------------------------------------------------------------------!
+         !   Updating and limitation of sea ice properties after transport   !
+         !-------------------------------------------------------------------!
          DO jj = 1, jpj
-!!!iii      zindhe = REAL( MAX( 0, isign(1, jj - njeq ) ) )              !ibug mpp  !!bugmpp  njeq!
             zindhe = MAX( 0.e0, SIGN( 1.e0, fcor(1,jj) ) )              ! = 0 for SH, =1 for NH
             DO ji = 1, jpi
-
+               !
                ! Recover mean values over the grid squares.
                zs0sn (ji,jj) = MAX( rzero, zs0sn (ji,jj)/area(ji,jj) )
@@ -272,7 +258,7 @@
             END DO
          END DO
-         
+         !
       ENDIF
-      
+      !
    END SUBROUTINE lim_trp_2
 
@@ -284,17 +270,13 @@
       !! ** Purpose :   initialization of ice advection parameters
       !!
-      !! ** Method  : Read the namicetrp namelist and check the parameter 
-      !!       values called at the first timestep (nit000)
+      !! ** Method  :   Read the namicetrp namelist and check the parameter 
+      !!              values called at the first timestep (nit000)
       !!
       !! ** input   :   Namelist namicetrp
-      !!
-      !! history :
-      !!   2.0  !  03-08 (C. Ethe)  Original code
       !!-------------------------------------------------------------------
       NAMELIST/namicetrp/ bound
       !!-------------------------------------------------------------------
-
-      ! Read Namelist namicetrp
-      REWIND ( numnam_ice )
+      !
+      REWIND ( numnam_ice )      ! Read Namelist namicetrp
       READ   ( numnam_ice  , namicetrp )
       IF(lwp) THEN
@@ -304,5 +286,5 @@
          WRITE(numout,*) '   boundary conditions (0. no-slip, 1. free-slip) bound  = ', bound
       ENDIF
-            
+      !
    END SUBROUTINE lim_trp_init_2
 
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/LIM_SRC_2/limwri_2.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/LIM_SRC_2/limwri_2.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/LIM_SRC_2/limwri_2.F90	(revision 2007)
@@ -313,6 +313,6 @@
       CALL histwrite( kid, "iicetemp", kt, sist(:,:) - rt0, jpi*jpj, (/1/) )
       CALL histwrite( kid, "ioceflxb", kt, fbif           , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicevelv", kt, u_ice          , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicevelu", kt, v_ice          , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicevelu", kt, u_ice          , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicevelv", kt, v_ice          , jpi*jpj, (/1/) )
       CALL histwrite( kid, "isstempe", kt, sst_m          , jpi*jpj, (/1/) )
       CALL histwrite( kid, "isssalin", kt, sss_m          , jpi*jpj, (/1/) )
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/LIM_SRC_3/limmsh.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/LIM_SRC_3/limmsh.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/LIM_SRC_3/limmsh.F90	(revision 2007)
@@ -53,4 +53,7 @@
          WRITE(numout,*) '~~~~~~~'
       ENDIF
+
+      IF( jphgr_msh == 2 .OR. jphgr_msh == 3 .OR. jphgr_msh == 5 )   &
+          &      CALL ctl_stop(' Coriolis parameter in LIM not set for f- or beta-plane' )
 
       !                           !==  coriolis factor & Equator position ==!
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 2007)
@@ -204,5 +204,5 @@
          zusnit = 1.0 / REAL( initad ) 
 
-         IF ( MOD( nday , 2 ) == 0) THEN
+         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0) THEN        !==  odd ice time step:  adv_x then adv_y  ==!
             DO jk = 1,initad
                !--- ice open water area
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/NST_SRC/agrif2model.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/NST_SRC/agrif2model.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/NST_SRC/agrif2model.F90	(revision 2007)
@@ -73,4 +73,23 @@
 
    END SUBROUTINE Agrif_clustering_def
+
+   SUBROUTINE Agrif_comm_def(modelcomm)
+
+      !!---------------------------------------------
+      !!   *** ROUTINE Agrif_clustering_def ***
+      !!--------------------------------------------- 
+      Use Agrif_Types
+      Use lib_mpp
+
+      IMPLICIT NONE
+
+      INTEGER :: modelcomm
+
+#if defined key_mpp_mpi
+      modelcomm = mpi_comm_opa
+#endif
+      Return
+
+   END SUBROUTINE Agrif_comm_def
 #else
    SUBROUTINE Agrif2Model
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/NST_SRC/agrif_user.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/NST_SRC/agrif_user.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/NST_SRC/agrif_user.F90	(revision 2007)
@@ -73,4 +73,7 @@
       USE tradmp
 #endif
+#if defined key_obc   ||   defined key_esopa
+      USE obc_par
+#endif
       USE sol_oce
       USE in_out_manager
@@ -105,4 +108,8 @@
       ! no tracer damping on fine grids
       lk_tradmp = .FALSE.
+#endif
+#if defined key_obc || defined key_esopa
+      ! no open boundary on fine grids
+      lk_obc = .FALSE.
 #endif
       ! 1. Declaration of the type of variable which have to be interpolated
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OFF_SRC/DOM/domrea.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OFF_SRC/DOM/domrea.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OFF_SRC/DOM/domrea.F90	(revision 2007)
@@ -215,11 +215,10 @@
 
  
-         DO jk = 1,jpk
-            gdept(:,:,jk) = gdept_0(jk)
-            gdepw(:,:,jk) = gdepw_0(jk)
-         END DO
-         
-
          IF( ln_zps ) THEN   
+            ! Vertical coordinates and scales factors
+            CALL iom_get( inum4, jpdom_unknown, 'gdept_0', gdept_0 ) ! depth
+            CALL iom_get( inum4, jpdom_unknown, 'gdepw_0', gdepw_0 )
+            CALL iom_get( inum4, jpdom_unknown, 'e3t_0'  , e3t_0   )
+            CALL iom_get( inum4, jpdom_unknown, 'e3w_0'  , e3w_0   )
                                       ! z-coordinate - partial steps
             IF( nmsh <= 6 ) THEN                                   !    ! 3D vertical scale factors
@@ -233,5 +232,5 @@
             END IF
 
-            IF( nmsh <= 3 ) THEN                                   !    ! 3D depth
+            IF( iom_varid( inum4, 'gdept', ldstop = .FALSE. ) > 0 ) THEN
               CALL iom_get( inum4, jpdom_data, 'gdept', gdept ) ! scale factors
               CALL iom_get( inum4, jpdom_data, 'gdepw', gdepw )
@@ -240,4 +239,9 @@
               CALL iom_get( inum4, jpdom_data, 'hdepw', hdepw )
          
+              DO jk = 1,jpk
+                gdept(:,:,jk) = gdept_0(jk)
+                gdepw(:,:,jk) = gdepw_0(jk)
+              ENDDO
+
               DO jj = 1, jpj
                 DO ji = 1, jpi
@@ -252,12 +256,8 @@
                 END DO
               END DO
+
             ENDIF
 
          ENDIF
-         ! Vertical coordinates and scales factors
-         CALL iom_get( inum4, jpdom_unknown, 'gdept_0', gdept_0 ) ! depth
-         CALL iom_get( inum4, jpdom_unknown, 'gdepw_0', gdepw_0 )
-         CALL iom_get( inum4, jpdom_unknown, 'e3t_0'  , e3t_0   )
-         CALL iom_get( inum4, jpdom_unknown, 'e3w_0'  , e3w_0   )
 # endif
          IF( ln_zco ) THEN
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OFF_SRC/DOM/phycst.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OFF_SRC/DOM/phycst.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OFF_SRC/DOM/phycst.F90	(revision 2007)
@@ -37,6 +37,5 @@
       rtt      = 273.16_wp  ,  &  !: triple point of temperature (Kelvin)
       rt0      = 273.15_wp  ,  &  !: freezing point of water (Kelvin)
-      rau0     = 1020._wp   ,  &  !: volumic mass of reference (kg/m3)
-      rauw     = 1000._wp   ,  &  !: density of pure water (kg/m3)
+      rau0     = 1035._wp   ,  &  !: volumic mass of reference (kg/m3)
       rcp      =    4.e+3_wp,  &  !: ocean specific heat
       ro0cpr                      !: = 1. / ( rau0 * rcp )
@@ -127,5 +126,4 @@
       ro0cpr = 1. / ( rau0 * rcp )
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '          volumic mass of pure water         rauw   = ', rauw, ' kg/m^3'
       IF(lwp) WRITE(numout,*) '          volumic mass of reference          rau0   = ', rau0, ' kg/m^3'
       IF(lwp) WRITE(numout,*) '          ocean specific heat                rcp    = ', rcp
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OFF_SRC/lib_mpp.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OFF_SRC/lib_mpp.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OFF_SRC/lib_mpp.F90	(revision 2007)
@@ -105,5 +105,5 @@
    !! ========================= !!
 !$AGRIF_DO_NOT_TREAT
-#  include <mpif.h>
+   INCLUDE mpif.h
 !$AGRIF_END_DO_NOT_TREAT
 
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OFF_SRC/trc_oce.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OFF_SRC/trc_oce.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OFF_SRC/trc_oce.F90	(revision 2007)
@@ -39,5 +39,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
-   !! $Id: trc_oce.F90 1423 2009-05-06 16:22:01Z ctlod $ 
+   !! $Id: trc_oce.F90 1834 2010-04-14 11:54:19Z cetlod $ 
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    !!----------------------------------------------------------------------
@@ -126,11 +126,11 @@
       zrgb(1,51) =  3.162   ;   zrgb(2,51) = 0.22703   ;   zrgb(3,51) = 0.16599   ;   zrgb(4,51) = 0.46601
       zrgb(1,52) =  3.548   ;   zrgb(2,52) = 0.24433   ;   zrgb(3,52) = 0.17334   ;   zrgb(4,52) = 0.47313
-      zrgb(1,53) =  3.981   ;   zrgb(2,53) = 0.26301   ;   zrgb(3,53) = 0.18126   ;   zrgb(4,54) = 0.48080
-      zrgb(1,54) =  4.467   ;   zrgb(2,54) = 0.28320   ;   zrgb(3,54) = 0.18981   ;   zrgb(4,55) = 0.48909
-      zrgb(1,55) =  5.012   ;   zrgb(2,55) = 0.30502   ;   zrgb(3,55) = 0.19903   ;   zrgb(4,56) = 0.49803
-      zrgb(1,56) =  5.623   ;   zrgb(2,56) = 0.32858   ;   zrgb(3,56) = 0.20898   ;   zrgb(4,57) = 0.50768
-      zrgb(1,57) =  6.310   ;   zrgb(2,57) = 0.35404   ;   zrgb(3,57) = 0.21971   ;   zrgb(4,58) = 0.51810
-      zrgb(1,58) =  7.079   ;   zrgb(2,58) = 0.38154   ;   zrgb(3,58) = 0.23129   ;   zrgb(4,59) = 0.52934
-      zrgb(1,59) =  7.943   ;   zrgb(2,59) = 0.41125   ;   zrgb(3,59) = 0.24378   ;   zrgb(4,50) = 0.54147
+      zrgb(1,53) =  3.981   ;   zrgb(2,53) = 0.26301   ;   zrgb(3,53) = 0.18126   ;   zrgb(4,53) = 0.48080
+      zrgb(1,54) =  4.467   ;   zrgb(2,54) = 0.28320   ;   zrgb(3,54) = 0.18981   ;   zrgb(4,54) = 0.48909
+      zrgb(1,55) =  5.012   ;   zrgb(2,55) = 0.30502   ;   zrgb(3,55) = 0.19903   ;   zrgb(4,55) = 0.49803
+      zrgb(1,56) =  5.623   ;   zrgb(2,56) = 0.32858   ;   zrgb(3,56) = 0.20898   ;   zrgb(4,56) = 0.50768
+      zrgb(1,57) =  6.310   ;   zrgb(2,57) = 0.35404   ;   zrgb(3,57) = 0.21971   ;   zrgb(4,57) = 0.51810
+      zrgb(1,58) =  7.079   ;   zrgb(2,58) = 0.38154   ;   zrgb(3,58) = 0.23129   ;   zrgb(4,58) = 0.52934
+      zrgb(1,59) =  7.943   ;   zrgb(2,59) = 0.41125   ;   zrgb(3,59) = 0.24378   ;   zrgb(4,59) = 0.54147
       zrgb(1,60) =  8.912   ;   zrgb(2,60) = 0.44336   ;   zrgb(3,60) = 0.25725   ;   zrgb(4,60) = 0.55457
       zrgb(1,61) = 10.000   ;   zrgb(2,61) = 0.47804   ;   zrgb(3,61) = 0.27178   ;   zrgb(4,61) = 0.56870
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90	(revision 2007)
@@ -175,7 +175,8 @@
       thick0(:,:) = 0.e0
       DO jk = 1, jpkm1
-         vol0        = vol0        + SUM( area (:,:) * tmask(:,:,jk) ) * e3t_0(jk) 
-         thick0(:,:) = thick0(:,:) +    tmask_i(:,:) * tmask(:,:,jk)   * e3t_0(jk)
-      END DO
+         vol0        = vol0        + SUM( area (:,:) * tmask(:,:,jk) * fse3t_0(:,:,jk) )
+         thick0(:,:) = thick0(:,:) +    tmask_i(:,:) * tmask(:,:,jk) * fse3t_0(:,:,jk)
+      END DO
+      IF( lk_mpp )   CALL mpp_sum( vol0 )
       
       CALL iom_open ( 'data_1m_salinity_nomask', inum )
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DIA/diadimg.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DIA/diadimg.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DIA/diadimg.F90	(revision 2007)
@@ -10,4 +10,5 @@
    USE dom_oce         ! ocean space and time domain
    USE in_out_manager  ! I/O manager
+   USE daymod          ! calendar
 
    IMPLICIT NONE
@@ -21,5 +22,5 @@
    !!----------------------------------------------------------------------
    !!   OPA 9.0 , LOCEAN-IPSL (2005) 
-   !! $Id$ 
+   !! $Header$ 
    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
    !!----------------------------------------------------------------------
@@ -56,6 +57,7 @@
     INTEGER :: jk, jn           ! dummy loop indices
     INTEGER :: irecl4,             &    ! record length in bytes
-         &       inum,             &    ! logical unit
-         &       irec                   ! current record to be written
+         &       inum,             &    ! logical unit (set to 14)
+         &       irec,             &    ! current record to be written
+         &       irecend                ! record number where nclit... are stored
     REAL(sp)                    :: zdx,zdy,zspval,zwest,ztimm
     REAL(sp)                    :: zsouth
@@ -69,5 +71,5 @@
     !! * Initialisations
 
-    irecl4 = MAX(jpi*jpj*sp , 84+18*sp + (jpk+8)*jpnij*sp  )
+    irecl4 = MAX(jpi*jpj*sp , 84+(18+1+jpk)*sp )
 
     zspval=0.0_sp    ! special values on land
@@ -101,6 +103,7 @@
 
     IF ( ln_dimgnnn  ) THEN
+     irecl4 = MAX(jpi*jpj*sp , 84+(18+jpk)*sp + 8*jpnij*sp  )
        WRITE(clname,'(a,a,i3.3)') TRIM(cd_name),'.',narea
-       CALL ctl_opn( inum, clname, 'REPLACE', 'UNFORMATTED', 'DIRECT', irecl4, numout, lwp )
+       CALL ctl_opn(inum, clname,'UNKNOWN','UNFORMATTED','DIRECT',irecl4,numout,lwp)
        WRITE(inum,REC=1 ) clver, cd_text, irecl4, &
             &     jpi,jpj, klev, 1 , 1 ,            &
@@ -127,16 +130,36 @@
        ENDIF
     ELSE
+       clver='@!03'           ! dimg string identifier
+       ! note that version @!02 is optimized with respect to record length.
+       ! The vertical dep variable is reduced to klev instead of klev*jpnij :
+       !   this is OK for jpnij < 181 (jpk=46)
+       ! for more processors, irecl4 get huge and that's why we switch to '@!03':
+       !  In this case we just add an extra integer to the standard dimg structure,
+       !  which is a record number where the arrays nlci etc... starts (1 per record)
+       
        !! Standard dimgproc (1 file per variable, all procs. write to this file )
        !! * Open file
-       CALL ctl_opn( inum, cd_name, 'REPLACE', 'UNFORMATTED', 'DIRECT', irecl4, numout, lwp )
+       CALL ctl_opn(inum, cd_name,'UNKNOWN','UNFORMATTED','DIRECT',irecl4,numout,lwp)
 
        !! * Write header on record #1
+       irecend=1 + klev*jpnij 
        IF(lwp) WRITE(inum,REC=1 ) clver, cd_text, irecl4, &
-            &     jpi,jpj, klev*jpnij, 1 , 1 ,            &
+            &     jpi,jpj, klev, 1 , 1 ,            &
             &     zwest, zsouth, zdx, zdy, zspval,  &
-            &     (z4dep(1:klev),jn=1,jpnij),       &
+            &     z4dep(1:klev),       &
             &     ztimm,                            &
-            &     narea, jpnij,jpiglo,jpjglo,jpizoom, jpjzoom,    &    ! extension to dimg for mpp output
-            &     nlcit,nlcjt, nldit, nldjt, nleit, nlejt, nimppt, njmppt  !
+            &     narea, jpnij,jpiglo,jpjglo,jpizoom, jpjzoom, irecend 
+       IF (lwp ) THEN
+         WRITE(inum,REC=irecend + 1 ) nlcit
+         WRITE(inum,REC=irecend + 2 ) nlcjt
+         WRITE(inum,REC=irecend + 3 ) nldit
+         WRITE(inum,REC=irecend + 4 ) nldjt
+         WRITE(inum,REC=irecend + 5 ) nleit
+         WRITE(inum,REC=irecend + 6 ) nlejt
+         WRITE(inum,REC=irecend + 7 ) nimppt
+         WRITE(inum,REC=irecend + 8 ) njmppt
+       ENDIF
+      !   &    ! extension to dimg for mpp output
+      !   &     nlcit,nlcjt, nldit, nldjt, nleit, nlejt, nimppt, njmppt  !
 
        !! * Write klev levels
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DIA/dianam.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DIA/dianam.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DIA/dianam.F90	(revision 2007)
@@ -129,7 +129,5 @@
  
       cdfnam = TRIM(cexper)//TRIM(clave)//"_"//TRIM(cldate1)//"_"//TRIM(cldate2)//"_"//TRIM(cdsuff)
-#if defined key_agrif
-      if ( .NOT. Agrif_Root() ) cdfnam = TRIM(Agrif_CFixed())//'_'//TRIM(cdfnam)
-#endif    
+      IF( .NOT. Agrif_Root() ) cdfnam = TRIM(Agrif_CFixed())//'_'//TRIM(cdfnam)
 
    END SUBROUTINE dia_nam
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 2007)
@@ -362,4 +362,95 @@
 #endif
 
+            ! Transports
+            ! T times V on T points (include bolus velocities)
+#if defined key_diaeiv 
+            DO jj = 2, jpj
+               DO ji = 1, jpi
+                  vt(ji,jj,:) = tn(ji,jj,:) * ( vn(ji,jj,:) + vn(ji,jj-1,:) + u_eiv(ji,jj,:) + u_eiv(ji,jj-1,:) )*0.5
+                  vs(ji,jj,:) = sn(ji,jj,:) * ( vn(ji,jj,:) + vn(ji,jj-1,:) + v_eiv(ji,jj,:) + v_eiv(ji,jj-1,:) )*0.5
+               END DO
+            END DO
+#else
+            DO jj = 2, jpj
+               DO ji = 1, jpi
+                  vt(ji,jj,:) = tn(ji,jj,:) * ( vn(ji,jj,:) + vn(ji,jj-1,:) )*0.5
+                  vs(ji,jj,:) = sn(ji,jj,:) * ( vn(ji,jj,:) + vn(ji,jj-1,:) )*0.5
+               END DO
+            END DO
+#endif 
+            CALL lbc_lnk( vs, 'V', -1. )   ;   CALL lbc_lnk( vt, 'V', -1. )
+
+            ht_glo(:) = SUM( ptr_vjk( vt(:,:,:)), 2 )
+            st_glo(:) = SUM( ptr_vjk( vs(:,:,:)), 2 )
+
+            IF ( ln_subbas ) THEN 
+               ht_atl(:) = SUM( ptr_vjk( vt (:,:,:), abasin(:,:)*sbasin(:,:)), 2 )
+               ht_pac(:) = SUM( ptr_vjk( vt (:,:,:), pbasin(:,:)*sbasin(:,:)), 2 )
+               ht_ind(:) = SUM( ptr_vjk( vt (:,:,:), ibasin(:,:)*sbasin(:,:)), 2 )
+               ht_ipc(:) = SUM( ptr_vjk( vt (:,:,:), dbasin(:,:)*sbasin(:,:)), 2 )
+               st_atl(:) = SUM( ptr_vjk( vs (:,:,:), abasin(:,:)*sbasin(:,:)), 2 )
+               st_pac(:) = SUM( ptr_vjk( vs (:,:,:), pbasin(:,:)*sbasin(:,:)), 2 )
+               st_ind(:) = SUM( ptr_vjk( vs (:,:,:), ibasin(:,:)*sbasin(:,:)), 2 )
+               st_ipc(:) = SUM( ptr_vjk( vs (:,:,:), dbasin(:,:)*sbasin(:,:)), 2 )
+            ENDIF
+
+            ! poleward tracer transports: 
+            ! overturning components:
+            IF ( ln_ptrcomp ) THEN 
+               pht_ove_glo(:) = SUM( v_msf_glo(:,:) * tn_jk_glo(:,:), 2 )   ! SUM over jk
+               pst_ove_glo(:) = SUM( v_msf_glo(:,:) * sn_jk_glo(:,:), 2 )  
+               IF ( ln_subbas ) THEN 
+                  pht_ove_atl(:) = SUM( v_msf_atl(:,:) * tn_jk_atl(:,:), 2 )   ! SUM over jk
+                  pst_ove_atl(:) = SUM( v_msf_atl(:,:) * sn_jk_atl(:,:), 2 )  
+                  pht_ove_pac(:) = SUM( v_msf_pac(:,:) * tn_jk_pac(:,:), 2 )   ! SUM over jk
+                  pst_ove_pac(:) = SUM( v_msf_pac(:,:) * sn_jk_pac(:,:), 2 )  
+                  pht_ove_ind(:) = SUM( v_msf_ind(:,:) * tn_jk_ind(:,:), 2 )   ! SUM over jk
+                  pst_ove_ind(:) = SUM( v_msf_ind(:,:) * sn_jk_ind(:,:), 2 )  
+                  pht_ove_ipc(:) = SUM( v_msf_ipc(:,:) * tn_jk_ipc(:,:), 2 )   ! SUM over jk
+                  pst_ove_ipc(:) = SUM( v_msf_ipc(:,:) * sn_jk_ipc(:,:), 2 )  
+               END IF
+            END IF
+
+            ! Bolus component
+#if defined key_diaeiv
+            pht_eiv_glo(:) = SUM( v_msf_eiv_glo(:,:) * tn_jk_glo(:,:), 2 )   ! SUM over jk
+            pst_eiv_glo(:) = SUM( v_msf_eiv_glo(:,:) * sn_jk_glo(:,:), 2 )   ! SUM over jk
+            IF ( ln_subbas ) THEN 
+               pht_eiv_atl(:) = SUM( v_msf_eiv_glo(:,:) * tn_jk_atl(:,:), 2 )   ! SUM over jk
+               pst_eiv_atl(:) = SUM( v_msf_eiv_glo(:,:) * sn_jk_atl(:,:), 2 )   ! SUM over jk
+               pht_eiv_pac(:) = SUM( v_msf_eiv_pac(:,:) * tn_jk_pac(:,:), 2 )   ! SUM over jk
+               pst_eiv_pac(:) = SUM( v_msf_eiv_pac(:,:) * sn_jk_pac(:,:), 2 )   ! SUM over jk
+               pht_eiv_ind(:) = SUM( v_msf_eiv_ind(:,:) * tn_jk_ind(:,:), 2 )   ! SUM over jk
+               pst_eiv_ind(:) = SUM( v_msf_eiv_ind(:,:) * sn_jk_ind(:,:), 2 )   ! SUM over jk
+               pht_eiv_ipc(:) = SUM( v_msf_eiv_ipc(:,:) * tn_jk_ipc(:,:), 2 )   ! SUM over jk
+               pst_eiv_ipc(:) = SUM( v_msf_eiv_ipc(:,:) * sn_jk_ipc(:,:), 2 )   ! SUM over jk
+            ENDIF
+#endif
+
+            ! conversion in PW and G g
+            zpwatt = zpwatt * rau0 * rcp
+            pht_adv(:) = pht_adv(:) * zpwatt  
+            pht_ldf(:) = pht_ldf(:) * zpwatt
+            pst_adv(:) = pst_adv(:) * zggram
+            pst_ldf(:) = pst_ldf(:) * zggram
+            IF ( ln_ptrcomp ) THEN 
+               pht_ove_glo(:) = pht_ove_glo(:) * zpwatt
+               pst_ove_glo(:) = pst_ove_glo(:) * zggram
+            END IF
+#if defined key_diaeiv
+            pht_eiv_glo(:) = pht_eiv_glo(:) * zpwatt
+            pst_eiv_glo(:) = pst_eiv_glo(:) * zggram
+#endif
+            IF( ln_subbas ) THEN
+               ht_atl(:) = ht_atl(:) * zpwatt
+               ht_pac(:) = ht_pac(:) * zpwatt
+               ht_ind(:) = ht_ind(:) * zpwatt
+               ht_ipc(:) = ht_ipc(:) * zpwatt
+               st_atl(:) = st_atl(:) * zggram 
+               st_pac(:) = st_pac(:) * zggram
+               st_ind(:) = st_ind(:) * zggram
+               st_ipc(:) = st_ipc(:) * zggram
+            ENDIF
+
             ! "Meridional" Stream-Function
             DO jk = 2,jpk 
@@ -394,95 +485,4 @@
                v_msf_ind(:,:) = v_msf_ind(:,:) * zsverdrup
                v_msf_ipc(:,:) = v_msf_ipc(:,:) * zsverdrup
-            ENDIF
-
-            ! Transports
-            ! T times V on T points (include bolus velocities)
-#if defined key_diaeiv 
-            DO jj = 2, jpj
-               DO ji = 1, jpi
-                  vt(ji,jj,:) = tn(ji,jj,:) * ( vn(ji,jj,:) + vn(ji,jj-1,:) + u_eiv(ji,jj,:) + u_eiv(ji,jj-1,:) )*0.5
-                  vs(ji,jj,:) = sn(ji,jj,:) * ( vn(ji,jj,:) + vn(ji,jj-1,:) + v_eiv(ji,jj,:) + v_eiv(ji,jj-1,:) )*0.5
-               END DO
-            END DO
-#else
-            DO jj = 2, jpj
-               DO ji = 1, jpi
-                  vt(ji,jj,:) = tn(ji,jj,:) * ( vn(ji,jj,:) + vn(ji,jj-1,:) )*0.5
-                  vs(ji,jj,:) = sn(ji,jj,:) * ( vn(ji,jj,:) + vn(ji,jj-1,:) )*0.5
-               END DO
-            END DO
-#endif 
-            CALL lbc_lnk( vs, 'V', -1. )   ;   CALL lbc_lnk( vt, 'V', -1. )
-
-            ht_glo(:) = SUM( ptr_vjk( vt(:,:,:)), 2 )
-            st_glo(:) = SUM( ptr_vjk( vs(:,:,:)), 2 )
-
-            IF ( ln_subbas ) THEN 
-               ht_atl(:) = SUM( ptr_vjk( vt (:,:,:), abasin(:,:)*sbasin(:,:)), 2 )
-               ht_pac(:) = SUM( ptr_vjk( vt (:,:,:), pbasin(:,:)*sbasin(:,:)), 2 )
-               ht_ind(:) = SUM( ptr_vjk( vt (:,:,:), ibasin(:,:)*sbasin(:,:)), 2 )
-               ht_ipc(:) = SUM( ptr_vjk( vt (:,:,:), dbasin(:,:)*sbasin(:,:)), 2 )
-               st_atl(:) = SUM( ptr_vjk( vs (:,:,:), abasin(:,:)*sbasin(:,:)), 2 )
-               st_pac(:) = SUM( ptr_vjk( vs (:,:,:), pbasin(:,:)*sbasin(:,:)), 2 )
-               st_ind(:) = SUM( ptr_vjk( vs (:,:,:), ibasin(:,:)*sbasin(:,:)), 2 )
-               st_ipc(:) = SUM( ptr_vjk( vs (:,:,:), dbasin(:,:)*sbasin(:,:)), 2 )
-            ENDIF
-
-            ! poleward tracer transports: 
-            ! overturning components:
-            IF ( ln_ptrcomp ) THEN 
-               pht_ove_glo(:) = SUM( v_msf_glo(:,:) * tn_jk_glo(:,:), 2 )   ! SUM over jk
-               pst_ove_glo(:) = SUM( v_msf_glo(:,:) * sn_jk_glo(:,:), 2 )  
-               IF ( ln_subbas ) THEN 
-                  pht_ove_atl(:) = SUM( v_msf_atl(:,:) * tn_jk_atl(:,:), 2 )   ! SUM over jk
-                  pst_ove_atl(:) = SUM( v_msf_atl(:,:) * sn_jk_atl(:,:), 2 )  
-                  pht_ove_pac(:) = SUM( v_msf_pac(:,:) * tn_jk_pac(:,:), 2 )   ! SUM over jk
-                  pst_ove_pac(:) = SUM( v_msf_pac(:,:) * sn_jk_pac(:,:), 2 )  
-                  pht_ove_ind(:) = SUM( v_msf_ind(:,:) * tn_jk_ind(:,:), 2 )   ! SUM over jk
-                  pst_ove_ind(:) = SUM( v_msf_ind(:,:) * sn_jk_ind(:,:), 2 )  
-                  pht_ove_ipc(:) = SUM( v_msf_ipc(:,:) * tn_jk_ipc(:,:), 2 )   ! SUM over jk
-                  pst_ove_ipc(:) = SUM( v_msf_ipc(:,:) * sn_jk_ipc(:,:), 2 )  
-               END IF
-            END IF
-
-            ! Bolus component
-#if defined key_diaeiv
-            pht_eiv_glo(:) = SUM( v_msf_eiv_glo(:,:) * tn_jk_glo(:,:), 2 )   ! SUM over jk
-            pst_eiv_glo(:) = SUM( v_msf_eiv_glo(:,:) * sn_jk_glo(:,:), 2 )   ! SUM over jk
-            IF ( ln_subbas ) THEN 
-               pht_eiv_atl(:) = SUM( v_msf_eiv_glo(:,:) * tn_jk_atl(:,:), 2 )   ! SUM over jk
-               pst_eiv_atl(:) = SUM( v_msf_eiv_glo(:,:) * sn_jk_atl(:,:), 2 )   ! SUM over jk
-               pht_eiv_pac(:) = SUM( v_msf_eiv_pac(:,:) * tn_jk_pac(:,:), 2 )   ! SUM over jk
-               pst_eiv_pac(:) = SUM( v_msf_eiv_pac(:,:) * sn_jk_pac(:,:), 2 )   ! SUM over jk
-               pht_eiv_ind(:) = SUM( v_msf_eiv_ind(:,:) * tn_jk_ind(:,:), 2 )   ! SUM over jk
-               pst_eiv_ind(:) = SUM( v_msf_eiv_ind(:,:) * sn_jk_ind(:,:), 2 )   ! SUM over jk
-               pht_eiv_ipc(:) = SUM( v_msf_eiv_ipc(:,:) * tn_jk_ipc(:,:), 2 )   ! SUM over jk
-               pst_eiv_ipc(:) = SUM( v_msf_eiv_ipc(:,:) * sn_jk_ipc(:,:), 2 )   ! SUM over jk
-            ENDIF
-#endif
-
-            ! conversion in PW and G g
-            zpwatt = zpwatt * rau0 * rcp
-            pht_adv(:) = pht_adv(:) * zpwatt  
-            pht_ldf(:) = pht_ldf(:) * zpwatt
-            pst_adv(:) = pst_adv(:) * zggram
-            pst_ldf(:) = pst_ldf(:) * zggram
-            IF ( ln_ptrcomp ) THEN 
-               pht_ove_glo(:) = pht_ove_glo(:) * zpwatt
-               pst_ove_glo(:) = pst_ove_glo(:) * zggram
-            END IF
-#if defined key_diaeiv
-            pht_eiv_glo(:) = pht_eiv_glo(:) * zpwatt
-            pst_eiv_glo(:) = pst_eiv_glo(:) * zggram
-#endif
-            IF( ln_subbas ) THEN
-               ht_atl(:) = ht_atl(:) * zpwatt
-               ht_pac(:) = ht_pac(:) * zpwatt
-               ht_ind(:) = ht_ind(:) * zpwatt
-               ht_ipc(:) = ht_ipc(:) * zpwatt
-               st_atl(:) = st_atl(:) * zggram 
-               st_pac(:) = st_pac(:) * zggram
-               st_ind(:) = st_ind(:) * zggram
-               st_ipc(:) = st_ipc(:) * zggram
             ENDIF
          ENDIF
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 2007)
@@ -629,7 +629,5 @@
       ! Define name, frequency of output and means
       clname = cdfile_name
-#if defined key_agrif
-      if ( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
-#endif
+      IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
       zdt  = rdt
       zsto = rdt
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 2007)
@@ -219,4 +219,13 @@
 #else
    LOGICAL, PUBLIC, PARAMETER ::   lk_agrif = .FALSE.   !: agrif flag
+
+CONTAINS
+   LOGICAL FUNCTION Agrif_Root()
+      Agrif_Root = .TRUE.
+   END FUNCTION Agrif_Root
+
+   CHARACTER(len=3) FUNCTION Agrif_CFixed()
+     Agrif_CFixed = '0' 
+   END FUNCTION Agrif_CFixed
 #endif
 
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 2007)
@@ -166,21 +166,17 @@
       ENDIF
 
-#if defined key_agrif
       IF( Agrif_Root() ) THEN
-#endif
-      SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
-      CASE (  1 ) 
-         CALL ioconf_calendar('gregorian')
-         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year'
-      CASE (  0 )
-         CALL ioconf_calendar('noleap')
-         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year'
-      CASE ( 30 )
-         CALL ioconf_calendar('360d')
-         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year'
-      END SELECT
-#if defined key_agrif
-      ENDIF
-#endif
+         SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
+         CASE (  1 ) 
+            CALL ioconf_calendar('gregorian')
+            IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year'
+         CASE (  0 )
+            CALL ioconf_calendar('noleap')
+            IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year'
+         CASE ( 30 )
+            CALL ioconf_calendar('360d')
+            IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year'
+         END SELECT
+      ENDIF
 
       REWIND( numnam )             ! Namelist namdom : space & time domain (bathymetry, mesh, timestep)
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 2007)
@@ -270,5 +270,5 @@
          
 #if defined key_agrif && defined key_eel_r6
-         IF (.Not.Agrif_Root()) THEN
+         IF( .NOT. Agrif_Root() ) THEN
            glam0  = Agrif_Parent(glam0) + (Agrif_ix())*Agrif_Parent(ppe1_m) * 1.e-3
            gphi0  = Agrif_Parent(gphi0) + (Agrif_iy())*Agrif_Parent(ppe2_m) * 1.e-3
@@ -465,5 +465,5 @@
          
 #if defined key_agrif && defined key_eel_r6
-         IF (.Not.Agrif_Root()) THEN
+         IF( .NOT. Agrif_Root() ) THEN
            zphi0 = ppgphi0 - FLOAT( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad)
          ENDIF
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 2007)
@@ -62,14 +62,29 @@
       IF( lk_zco )   CALL ctl_stop( 'dom_vvl : key_zco is incompatible with variable volume option key_vvl')
 
-      fsdept(:,:,:) = gdept (:,:,:)
-      fsdepw(:,:,:) = gdepw (:,:,:)
-      fsde3w(:,:,:) = gdep3w(:,:,:)
-      fse3t (:,:,:) = e3t   (:,:,:)
-      fse3u (:,:,:) = e3u   (:,:,:)
-      fse3v (:,:,:) = e3v   (:,:,:)
-      fse3f (:,:,:) = e3f   (:,:,:)
-      fse3w (:,:,:) = e3w   (:,:,:)
-      fse3uw(:,:,:) = e3uw  (:,:,:)
-      fse3vw(:,:,:) = e3vw  (:,:,:)
+      IF( ln_zco) THEN
+         DO jk = 1, jpk
+            gdept(:,:,jk) = gdept_0(jk)
+            gdepw(:,:,jk) = gdepw_0(jk)
+            gdep3w(:,:,jk) = gdepw_0(jk)
+            e3t (:,:,jk) = e3t_0(jk)
+            e3u (:,:,jk) = e3t_0(jk)
+            e3v (:,:,jk) = e3t_0(jk)
+            e3f (:,:,jk) = e3t_0(jk)
+            e3w (:,:,jk) = e3w_0(jk)
+            e3uw(:,:,jk) = e3w_0(jk)
+            e3vw(:,:,jk) = e3w_0(jk)
+         END DO
+      ELSE
+         fsdept(:,:,:) = gdept (:,:,:)
+         fsdepw(:,:,:) = gdepw (:,:,:)
+         fsde3w(:,:,:) = gdep3w(:,:,:)
+         fse3t (:,:,:) = e3t   (:,:,:)
+         fse3u (:,:,:) = e3u   (:,:,:)
+         fse3v (:,:,:) = e3v   (:,:,:)
+         fse3f (:,:,:) = e3f   (:,:,:)
+         fse3w (:,:,:) = e3w   (:,:,:)
+         fse3uw(:,:,:) = e3uw  (:,:,:)
+         fse3vw(:,:,:) = e3vw  (:,:,:)
+      ENDIF
 
       !                                 !==  mu computation  ==!
@@ -139,4 +154,20 @@
       CALL lbc_lnk( sshf_b, 'F', 1. )   ;   CALL lbc_lnk( sshf_n, 'F', 1. )
       !
+         DO jk = 1, jpkm1
+            fsdept(:,:,jk) = fsdept_n(:,:,jk)          ! now local depths stored in fsdep. arrays
+            fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
+            fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
+            !
+            fse3t (:,:,jk) = fse3t_n (:,:,jk)          ! vertical scale factors stored in fse3. arrays
+            fse3u (:,:,jk) = fse3u_n (:,:,jk)
+            fse3v (:,:,jk) = fse3v_n (:,:,jk)
+            fse3f (:,:,jk) = fse3f_n (:,:,jk)
+            fse3w (:,:,jk) = fse3w_n (:,:,jk)
+            fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
+            fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
+         END DO
+
+
+
    END SUBROUTINE dom_vvl
 
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90	(revision 2007)
@@ -45,10 +45,16 @@
       !!      domhgr, domzgr, and dommsk. Note: the file contain depends on
       !!      the vertical coord. used (z-coord, partial steps, s-coord)
-      !!                    nmsh = 1  :   'mesh_mask.nc' file
+      !!            MOD(nmsh, 3) = 1  :   'mesh_mask.nc' file
       !!                         = 2  :   'mesh.nc' and mask.nc' files
-      !!                         = 3  :   'mesh_hgr.nc', 'mesh_zgr.nc' and
+      !!                         = 0  :   'mesh_hgr.nc', 'mesh_zgr.nc' and
       !!                                  'mask.nc' files
       !!      For huge size domain, use option 2 or 3 depending on your 
       !!      vertical coordinate.
+      !!
+      !!      if     nmsh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw]
+      !!      if 3 < nmsh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays 
+      !!                        corresponding to the depth of the bottom points hdep[tw]
+      !!      if 6 < nmsh <= 9: write 2D arrays corresponding to the depth and the
+      !!                        thickness of the bottom points hdep[tw] and e3[tw]_ps
       !!
       !! ** output file : 
@@ -241,5 +247,5 @@
       !                                     !        close the files 
       !                                     ! ============================
-      SELECT CASE ( nmsh )
+      SELECT CASE ( MOD(nmsh, 3) )
       CASE ( 1 )                
          CALL iom_close( inum0 )
@@ -247,5 +253,5 @@
          CALL iom_close( inum1 )
          CALL iom_close( inum2 )
-      CASE ( 3 )
+      CASE ( 0 )
          CALL iom_close( inum2 )
          CALL iom_close( inum3 )
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90	(revision 2007)
@@ -4,9 +4,9 @@
    !!     Definition of of both ocean and ice parameters used in the code
    !!=====================================================================
-   !! History :        !  90-10  (C. Levy - G. Madec)  Original code
-   !!                  !  91-11  (G. Madec)
-   !!                  !  91-12  (M. Imbard)
-   !!             8.5  !  02-08  (G. Madec, C. Ethe)  F90, add ice constants
-   !!             9.0  !  06-08  (G. Madec)  style 
+   !! History :   OPA  !  1990-10  (C. Levy - G. Madec)  Original code
+   !!             8.1  !  1991-11  (G. Madec, M. Imbard)  cosmetic changes
+   !!   NEMO      1.0  !  2002-08  (G. Madec, C. Ethe)  F90, add ice constants
+   !!              -   !  2006-08  (G. Madec)  style 
+   !!             3.2  !  2006-08  (S. Masson, G. Madec)  suppress useless variables + style 
    !!----------------------------------------------------------------------
 
@@ -24,67 +24,64 @@
    REAL(wp), PUBLIC ::   rpi = 3.141592653589793_wp             !: pi
    REAL(wp), PUBLIC ::   rad = 3.141592653589793_wp / 180._wp   !: conversion from degre into radian
-   REAL(wp), PUBLIC ::   rsmall = 0.5 * EPSILON( 1. )           !: smallest real computer value
+   REAL(wp), PUBLIC ::   rsmall = 0.5 * EPSILON( 1.e0 )         !: smallest real computer value
    
-   REAL(wp), PUBLIC ::          & !:
-      rday = 24.*60.*60.  ,     & !: day (s)
-      rsiyea              ,     & !: sideral year (s)
-      rsiday              ,     & !: sideral day (s)
-      raamo =  12._wp     ,     & !: number of months in one year
-      rjjhh =  24._wp     ,     & !: number of hours in one day
-      rhhmm =  60._wp     ,     & !: number of minutes in one hour
-      rmmss =  60._wp     ,     & !: number of seconds in one minute
-!!!   omega = 7.292115083046061e-5_wp ,  &  !: change the last digit!
-      omega               ,    &  !: earth rotation parameter
-      ra    = 6371229._wp ,    &  !: earth radius (meter)
-      grav  = 9.80665_wp          !: gravity (m/s2)
+   REAL(wp), PUBLIC ::   rday = 24.*60.*60.       !: day (s)
+   REAL(wp), PUBLIC ::   rsiyea                   !: sideral year (s)
+   REAL(wp), PUBLIC ::   rsiday                   !: sideral day (s)
+   REAL(wp), PUBLIC ::   raamo =  12._wp          !: number of months in one year
+   REAL(wp), PUBLIC ::   rjjhh =  24._wp          !: number of hours in one day
+   REAL(wp), PUBLIC ::   rhhmm =  60._wp          !: number of minutes in one hour
+   REAL(wp), PUBLIC ::   rmmss =  60._wp          !: number of seconds in one minute
+!! REAL(wp), PUBLIC ::   omega = 7.292115083046061e-5_wp ,  &  !: change the last digit!
+   REAL(wp), PUBLIC ::   omega                    !: earth rotation parameter
+   REAL(wp), PUBLIC ::   ra    = 6371229._wp      !: earth radius (meter)
+   REAL(wp), PUBLIC ::   grav  = 9.80665_wp       !: gravity (m/s2)
    
-   REAL(wp), PUBLIC ::         &  !:
-      rtt      = 273.16_wp  ,  &  !: triple point of temperature (Kelvin)
-      rt0      = 273.15_wp  ,  &  !: freezing point of water (Kelvin)
+   REAL(wp), PUBLIC ::   rtt      = 273.16_wp     !: triple point of temperature (Kelvin)
+   REAL(wp), PUBLIC ::   rt0      = 273.15_wp     !: freezing point of water (Kelvin)
 #if defined key_lim3
-      rt0_snow = 273.16_wp  ,  &  !: melting point of snow  (Kelvin)
-      rt0_ice  = 273.16_wp  ,  &  !: melting point of ice   (Kelvin)
+   REAL(wp), PUBLIC ::   rt0_snow = 273.16_wp     !: melting point of snow  (Kelvin)
+   REAL(wp), PUBLIC ::   rt0_ice  = 273.16_wp     !: melting point of ice   (Kelvin)
 #else
-      rt0_snow = 273.15_wp  ,  &  !: melting point of snow  (Kelvin)
-      rt0_ice  = 273.05_wp  ,  &  !: melting point of ice   (Kelvin)
+   REAL(wp), PUBLIC ::   rt0_snow = 273.15_wp     !: melting point of snow  (Kelvin)
+   REAL(wp), PUBLIC ::   rt0_ice  = 273.05_wp     !: melting point of ice   (Kelvin)
 #endif
-      rau0     = 1035._wp   ,  &  !: volumic mass of reference (kg/m3)
-      rauw     = 1000._wp   ,  &  !: volumic mass of pure water (kg/m3)
-      rcp      =    4.e+3_wp,  &  !: ocean specific heat
-      ro0cpr                      !: = 1. / ( rau0 * rcp )
 
-   REAL(wp), PUBLIC ::            &  !:
+   REAL(wp), PUBLIC ::   rau0     = 1020._wp      !: reference volumic mass (density)  (kg/m3)
+   REAL(wp), PUBLIC ::   rau0r                    !: reference specific volume         (m3/kg)
+   REAL(wp), PUBLIC ::   rcp      =    4.e+3_wp   !: ocean specific heat
+   REAL(wp), PUBLIC ::   ro0cpr                   !: = 1. / ( rau0 * rcp )
+
 #if defined key_lim3
-      rcdsn   =   0.31_wp     ,   &  !: thermal conductivity of snow
-      rcdic   =   2.034396_wp ,   &  !: thermal conductivity of fresh ice
-      cpic    = 2067.0        ,   &
-      ! add the following lines
-      lsub    = 2.834e+6      ,   &  !: pure ice latent heat of sublimation (J.kg-1)
-      lfus    = 0.334e+6      ,   &  !: latent heat of fusion of fresh ice   (J.kg-1)
-      rhoic   = 917._wp       ,   &  !: volumic mass of sea ice (kg/m3)
-      tmut    =   0.054       ,   &  !: decrease of seawater meltpoint with salinity
+   REAL(wp), PUBLIC ::   rcdsn   =   0.31_wp      !: thermal conductivity of snow
+   REAL(wp), PUBLIC ::   rcdic   =   2.034396_wp  !: thermal conductivity of fresh ice
+   REAL(wp), PUBLIC ::   cpic    = 2067.0         !: specific heat of sea ice
+   REAL(wp), PUBLIC ::   lsub    = 2.834e+6       !: pure ice latent heat of sublimation (J.kg-1)
+   REAL(wp), PUBLIC ::   lfus    = 0.334e+6       !: latent heat of fusion of fresh ice   (J.kg-1)
+   REAL(wp), PUBLIC ::   rhoic   = 917._wp        !: volumic mass of sea ice (kg/m3)
+   REAL(wp), PUBLIC ::   tmut    =   0.054        !: decrease of seawater meltpoint with salinity
 #else
-      rcdsn   =   0.22_wp     ,   &  !: conductivity of the snow
-      rcdic   =   2.034396_wp ,   &  !: conductivity of the ice
-      rcpsn   =   6.9069e+5_wp,   &  !: density times specific heat for snow
-      rcpic   =   1.8837e+6_wp,   &  !: volumetric latent heat fusion of sea ice
-      xlsn    = 110.121e+6_wp ,   &  !: volumetric latent heat fusion of snow
-      xlic    = 300.33e+6_wp  ,   &  !: volumetric latent heat fusion of ice
-      xsn     =   2.8e+6      ,   &  !: latent heat of sublimation of snow
-      rhoic   = 900._wp       ,   &  !: volumic mass of sea ice (kg/m3)
+   REAL(wp), PUBLIC ::   rcdsn   =   0.22_wp      !: conductivity of the snow
+   REAL(wp), PUBLIC ::   rcdic   =   2.034396_wp  !: conductivity of the ice
+   REAL(wp), PUBLIC ::   rcpsn   =   6.9069e+5_wp !: density times specific heat for snow
+   REAL(wp), PUBLIC ::   rcpic   =   1.8837e+6_wp !: volumetric latent heat fusion of sea ice
+   REAL(wp), PUBLIC ::   xlsn    = 110.121e+6_wp  !: volumetric latent heat fusion of snow
+   REAL(wp), PUBLIC ::   xlic    = 300.33e+6_wp   !: volumetric latent heat fusion of ice
+   REAL(wp), PUBLIC ::   xsn     =   2.8e+6       !: latent heat of sublimation of snow
+   REAL(wp), PUBLIC ::   rhoic   = 900._wp        !: volumic mass of sea ice (kg/m3)
 #endif
-      rhosn   = 330._wp       ,   &  !: volumic mass of snow (kg/m3)
-      emic    =   0.97_wp     ,   &  !: emissivity of snow or ice
-      sice    =   6.0_wp      ,   &  !: salinity of ice (psu)
-      soce    =  34.7_wp      ,   &  !: salinity of sea (psu)
-      cevap   =   2.5e+6_wp   ,   &  !: latent heat of evaporation (water)
-      srgamma =   0.9_wp      ,   &  !: correction factor for solar radiation (Oberhuber, 1974)
-      vkarmn  =   0.4_wp      ,   &  !: von Karman constant
-      stefan  =   5.67e-8_wp         !: Stefan-Boltzmann constant 
-      !!----------------------------------------------------------------------
-      !!  OPA 9.0 , LOCEAN-IPSL (2005) 
-      !! $Id$ 
-      !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
-      !!----------------------------------------------------------------------
+   REAL(wp), PUBLIC ::   rhosn   = 330._wp        !: volumic mass of snow (kg/m3)
+   REAL(wp), PUBLIC ::   emic    =   0.97_wp      !: emissivity of snow or ice
+   REAL(wp), PUBLIC ::   sice    =   6.0_wp       !: reference salinity of ice (psu)
+   REAL(wp), PUBLIC ::   soce    =  34.7_wp       !: reference salinity of sea (psu)
+   REAL(wp), PUBLIC ::   cevap   =   2.5e+6_wp    !: latent heat of evaporation (water)
+   REAL(wp), PUBLIC ::   srgamma =   0.9_wp       !: correction factor for solar radiation (Oberhuber, 1974)
+   REAL(wp), PUBLIC ::   vkarmn  =   0.4_wp       !: von Karman constant
+   REAL(wp), PUBLIC ::   stefan  =   5.67e-8_wp   !: Stefan-Boltzmann constant 
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
+   !! $Id$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
    
 CONTAINS
@@ -99,11 +96,17 @@
       !!----------------------------------------------------------------------
 
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) ' phy_cst : initialization of ocean parameters and constants'
-      IF(lwp) WRITE(numout,*) ' ~~~~~~~'
+      !                                   ! Define additional parameters
+      rsiyea = 365.25 * rday * 2. * rpi / 6.283076
+      rsiday = rday / ( 1. + rday / rsiyea )
+      omega  = 2. * rpi / rsiday 
 
-      ! Ocean Parameters
-      ! ----------------
-      IF(lwp) THEN
+      rau0r  = 1. /   rau0  
+      ro0cpr = 1. / ( rau0 * rcp )
+
+
+      IF(lwp) THEN                        ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) ' phy_cst : initialization of ocean parameters and constants'
+         WRITE(numout,*) ' ~~~~~~~'
          WRITE(numout,*) '       Domain info'
          WRITE(numout,*) '          dimension of model'
@@ -118,47 +121,30 @@
          WRITE(numout,*) '             jpnij   : ', jpnij
          WRITE(numout,*) '          lateral domain boundary condition type : jperio  = ', jperio
-      ENDIF
-
-      ! Define constants
-      ! ----------------
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '       Constants'
-
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '          mathematical constant                 rpi = ', rpi
-
-      rsiyea = 365.25 * rday * 2. * rpi / 6.283076
-      rsiday = rday / ( 1. + rday / rsiyea )
-      omega  = 2. * rpi / rsiday 
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '          day                                rday   = ', rday,   ' s'
-      IF(lwp) WRITE(numout,*) '          sideral year                       rsiyea = ', rsiyea, ' s'
-      IF(lwp) WRITE(numout,*) '          sideral day                        rsiday = ', rsiday, ' s'
-      IF(lwp) WRITE(numout,*) '          omega                              omega  = ', omega,  ' s-1'
-
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '          nb of months per year               raamo = ', raamo, ' months'
-      IF(lwp) WRITE(numout,*) '          nb of hours per day                 rjjhh = ', rjjhh, ' hours'
-      IF(lwp) WRITE(numout,*) '          nb of minutes per hour              rhhmm = ', rhhmm, ' mn'
-      IF(lwp) WRITE(numout,*) '          nb of seconds per minute            rmmss = ', rmmss, ' s'
-
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '          earth radius                         ra   = ', ra, ' m'
-      IF(lwp) WRITE(numout,*) '          gravity                              grav = ', grav , ' m/s^2'
-
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '          triple point of temperature      rtt      = ', rtt     , ' K'
-      IF(lwp) WRITE(numout,*) '          freezing point of water          rt0      = ', rt0     , ' K'
-      IF(lwp) WRITE(numout,*) '          melting point of snow            rt0_snow = ', rt0_snow, ' K'
-      IF(lwp) WRITE(numout,*) '          melting point of ice             rt0_ice  = ', rt0_ice , ' K'
-
-      ro0cpr = 1. / ( rau0 * rcp )
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '          volumic mass of pure water         rauw   = ', rauw, ' kg/m^3'
-      IF(lwp) WRITE(numout,*) '          volumic mass of reference          rau0   = ', rau0, ' kg/m^3'
-      IF(lwp) WRITE(numout,*) '          ocean specific heat                rcp    = ', rcp
-      IF(lwp) WRITE(numout,*) '                       1. / ( rau0 * rcp ) = ro0cpr = ', ro0cpr
-
-      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) '       Constants'
+         WRITE(numout,*)
+         WRITE(numout,*) '          mathematical constant                 rpi = ', rpi
+         WRITE(numout,*) '          day                                rday   = ', rday,   ' s'
+         WRITE(numout,*) '          sideral year                       rsiyea = ', rsiyea, ' s'
+         WRITE(numout,*) '          sideral day                        rsiday = ', rsiday, ' s'
+         WRITE(numout,*) '          omega                              omega  = ', omega,  ' s-1'
+         WRITE(numout,*)
+         WRITE(numout,*) '          nb of months per year               raamo = ', raamo, ' months'
+         WRITE(numout,*) '          nb of hours per day                 rjjhh = ', rjjhh, ' hours'
+         WRITE(numout,*) '          nb of minutes per hour              rhhmm = ', rhhmm, ' mn'
+         WRITE(numout,*) '          nb of seconds per minute            rmmss = ', rmmss, ' s'
+         WRITE(numout,*)
+         WRITE(numout,*) '          earth radius                         ra   = ', ra, ' m'
+         WRITE(numout,*) '          gravity                              grav = ', grav , ' m/s^2'
+         WRITE(numout,*)
+         WRITE(numout,*) '          triple point of temperature      rtt      = ', rtt     , ' K'
+         WRITE(numout,*) '          freezing point of water          rt0      = ', rt0     , ' K'
+         WRITE(numout,*) '          melting point of snow            rt0_snow = ', rt0_snow, ' K'
+         WRITE(numout,*) '          melting point of ice             rt0_ice  = ', rt0_ice , ' K'
+         WRITE(numout,*)
+         WRITE(numout,*) '          ocean reference volumic mass       rau0   = ', rau0 , ' kg/m^3'
+         WRITE(numout,*) '          ocean reference specific volume    rau0r  = ', rau0r, ' m^3/Kg'
+         WRITE(numout,*) '          ocean specific heat                rcp    = ', rcp
+         WRITE(numout,*) '                       1. / ( rau0 * rcp ) = ro0cpr = ', ro0cpr
          WRITE(numout,*)
          WRITE(numout,*) '          thermal conductivity of the snow          = ', rcdsn   , ' J/s/m/K'
@@ -184,8 +170,6 @@
          WRITE(numout,*) '          von Karman constant                       = ', vkarmn 
          WRITE(numout,*) '          Stefan-Boltzmann constant                 = ', stefan  , ' J/s/m^2/K^4'
-
          WRITE(numout,*)
          WRITE(numout,*) '          conversion: degre ==> radian          rad = ', rad
-
          WRITE(numout,*)
          WRITE(numout,*) '          smallest real computer value       rsmall = ', rsmall
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90	(revision 2007)
@@ -123,33 +123,27 @@
 
 #if defined key_obc
-#if defined key_agrif
-         IF (Agrif_Root() ) THEN
-#endif
-         ! open boundaries (div must be zero behind the open boundary)
-         !  mpp remark: The zeroing of hdivn can probably be extended to 1->jpi/jpj for the correct row/column
-         IF( lp_obc_east  )   hdivn(nie0p1:nie1p1,nje0  :nje1  ,jk) = 0.e0      ! east
-         IF( lp_obc_west  )   hdivn(niw0  :niw1  ,njw0  :njw1  ,jk) = 0.e0      ! west
-         IF( lp_obc_north )   hdivn(nin0  :nin1  ,njn0p1:njn1p1,jk) = 0.e0      ! north
-         IF( lp_obc_south )   hdivn(nis0  :nis1  ,njs0  :njs1  ,jk) = 0.e0      ! south
-#if defined key_agrif
-         ENDIF
-#endif
+         IF( Agrif_Root() ) THEN
+            ! open boundaries (div must be zero behind the open boundary)
+            !  mpp remark: The zeroing of hdivn can probably be extended to 1->jpi/jpj for the correct row/column
+            IF( lp_obc_east  )   hdivn(nie0p1:nie1p1,nje0  :nje1  ,jk) = 0.e0      ! east
+            IF( lp_obc_west  )   hdivn(niw0  :niw1  ,njw0  :njw1  ,jk) = 0.e0      ! west
+            IF( lp_obc_north )   hdivn(nin0  :nin1  ,njn0p1:njn1p1,jk) = 0.e0      ! north
+            IF( lp_obc_south )   hdivn(nis0  :nis1  ,njs0  :njs1  ,jk) = 0.e0      ! south
+         ENDIF
 #endif         
 #if defined key_bdy
          ! unstructured open boundaries (div must be zero behind the open boundary)
          DO jj = 1, jpj
-           DO ji = 1, jpi
-             hdivn(ji,jj,jk)=hdivn(ji,jj,jk)*bdytmask(ji,jj)
-           END DO
+            DO ji = 1, jpi
+               hdivn(ji,jj,jk)=hdivn(ji,jj,jk)*bdytmask(ji,jj)
+            END DO
          END DO
 #endif         
-#if defined key_agrif
-         if ( .NOT. AGRIF_Root() ) then
-           IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn(nlci-1 , :     ,jk) = 0.e0      ! east
-           IF ((nbondi == -1).OR.(nbondi == 2)) hdivn(2      , :     ,jk) = 0.e0      ! west
-           IF ((nbondj ==  1).OR.(nbondj == 2)) hdivn(:      ,nlcj-1 ,jk) = 0.e0      ! north
-           IF ((nbondj == -1).OR.(nbondj == 2)) hdivn(:      ,2      ,jk) = 0.e0      ! south
-         endif
-#endif    
+         IF( .NOT. AGRIF_Root() ) THEN
+            IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn(nlci-1 , :     ,jk) = 0.e0      ! east
+            IF ((nbondi == -1).OR.(nbondi == 2)) hdivn(2      , :     ,jk) = 0.e0      ! west
+            IF ((nbondj ==  1).OR.(nbondj == 2)) hdivn(:      ,nlcj-1 ,jk) = 0.e0      ! north
+            IF ((nbondj == -1).OR.(nbondj == 2)) hdivn(:      ,2      ,jk) = 0.e0      ! south
+         ENDIF
 
          !                                             ! --------
@@ -341,16 +335,12 @@
 
 #if defined key_obc
-#if defined key_agrif
-         IF ( Agrif_Root() ) THEN
-#endif
-         ! open boundaries (div must be zero behind the open boundary)
-         !  mpp remark: The zeroing of hdivn can probably be extended to 1->jpi/jpj for the correct row/column
-         IF( lp_obc_east  )   hdivn(nie0p1:nie1p1,nje0  :nje1  ,jk) = 0.e0      ! east
-         IF( lp_obc_west  )   hdivn(niw0  :niw1  ,njw0  :njw1  ,jk) = 0.e0      ! west
-         IF( lp_obc_north )   hdivn(nin0  :nin1  ,njn0p1:njn1p1,jk) = 0.e0      ! north
-         IF( lp_obc_south )   hdivn(nis0  :nis1  ,njs0  :njs1  ,jk) = 0.e0      ! south
-#if defined key_agrif
-         ENDIF
-#endif
+         IF( Agrif_Root() ) THEN
+            ! open boundaries (div must be zero behind the open boundary)
+            !  mpp remark: The zeroing of hdivn can probably be extended to 1->jpi/jpj for the correct row/column
+            IF( lp_obc_east  )   hdivn(nie0p1:nie1p1,nje0  :nje1  ,jk) = 0.e0      ! east
+            IF( lp_obc_west  )   hdivn(niw0  :niw1  ,njw0  :njw1  ,jk) = 0.e0      ! west
+            IF( lp_obc_north )   hdivn(nin0  :nin1  ,njn0p1:njn1p1,jk) = 0.e0      ! north
+            IF( lp_obc_south )   hdivn(nis0  :nis1  ,njs0  :njs1  ,jk) = 0.e0      ! south
+         ENDIF
 #endif         
 #if defined key_bdy
@@ -362,12 +352,10 @@
          END DO
 #endif        
-#if defined key_agrif
-         if ( .NOT. AGRIF_Root() ) then
+         IF( .NOT. AGRIF_Root() ) THEN
             IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn(nlci-1 , :     ,jk) = 0.e0      ! east
             IF ((nbondi == -1).OR.(nbondi == 2)) hdivn(2      , :     ,jk) = 0.e0      ! west
             IF ((nbondj ==  1).OR.(nbondj == 2)) hdivn(:      ,nlcj-1 ,jk) = 0.e0      ! north
             IF ((nbondj == -1).OR.(nbondj == 2)) hdivn(:      ,2      ,jk) = 0.e0      ! south
-         endif
-#endif    
+         ENDIF
 
          !                                             ! --------
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 2007)
@@ -146,5 +146,5 @@
 # if defined key_obc
       !                                !* OBC open boundaries
-      CALL obc_dyn( kt )
+      IF( lk_obc )   CALL obc_dyn( kt )
       !
       IF ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90	(revision 2007)
@@ -186,6 +186,6 @@
 
 #if defined key_obc
-      CALL obc_dyn( kt )      ! Update velocities on each open boundary with the radiation algorithm
-      CALL obc_vol( kt )      ! Correction of the barotropic componant velocity to control the volume of the system
+      IF( lk_obc )   CALL obc_dyn( kt )   ! Update velocities on each open boundary with the radiation algorithm
+      IF( lk_obc )   CALL obc_vol( kt )   ! Correction of the barotropic componant velocity to control the volume of the system
 #endif
 #if defined key_bdy
@@ -315,6 +315,11 @@
 #if defined key_obc
             ! caution : grad D = 0 along open boundaries
-            spgu(ji,jj) = z2dt * ztdgu * obcumask(ji,jj)
-            spgv(ji,jj) = z2dt * ztdgv * obcvmask(ji,jj)
+            IF( Agrif_Root() ) THEN
+               spgu(ji,jj) = z2dt * ztdgu * obcumask(ji,jj)
+               spgv(ji,jj) = z2dt * ztdgv * obcvmask(ji,jj)
+            ELSE
+               spgu(ji,jj) = z2dt * ztdgu
+               spgv(ji,jj) = z2dt * ztdgv
+            ENDIF
 #elif defined key_bdy
             ! caution : grad D = 0 along open boundaries
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 2007)
@@ -157,12 +157,8 @@
 
 #if defined key_obc
-# if defined key_agrif
       IF ( Agrif_Root() ) THEN 
-# endif
          ssha(:,:) = ssha(:,:) * obctmsk(:,:)
          CALL lbc_lnk( ssha, 'T', 1. )  ! absolutly compulsory !! (jmm)
-# if defined key_agrif
-      ENDIF
-# endif
+      ENDIF
 #endif
 
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 2007)
@@ -43,5 +43,5 @@
    LOGICAL, PUBLIC, PARAMETER ::   lk_iomput = .FALSE.       !: iom_put flag
 #endif
-   PUBLIC iom_init, iom_open, iom_close, iom_setkt, iom_varid, iom_get, iom_gettime, iom_rstput, iom_put
+   PUBLIC iom_init, iom_swap, iom_open, iom_close, iom_setkt, iom_varid, iom_get, iom_gettime, iom_rstput, iom_put
 
    PRIVATE iom_rp0d, iom_rp1d, iom_rp2d, iom_rp3d
@@ -86,5 +86,6 @@
       !!----------------------------------------------------------------------
       ! read the xml file
-      CALL event__parse_xml_file( 'iodef.xml' )   ! <- to get from the nameliste (namrun)...
+      IF( Agrif_Root() ) CALL event__parse_xml_file( 'iodef.xml' )   ! <- to get from the nameliste (namrun)...
+      CALL iom_swap
 
       ! calendar parameters
@@ -119,4 +120,22 @@
 
    END SUBROUTINE iom_init
+
+
+   SUBROUTINE iom_swap
+      !!---------------------------------------------------------------------
+      !!                   ***  SUBROUTINE  iom_swap  ***
+      !!
+      !! ** Purpose :  swap context between different agrif grid for xmlio_server
+      !!---------------------------------------------------------------------
+#if defined key_iomput
+
+     IF( TRIM(Agrif_CFixed()) == '0' ) THEN
+        CALL event__swap_context("nemo")
+     ELSE
+        CALL event__swap_context(TRIM(Agrif_CFixed())//"_nemo")
+     ENDIF
+
+#endif
+   END SUBROUTINE iom_swap
 
 
@@ -164,10 +183,4 @@
       ! if iom_open is called for the first time: initialize iom_file(:)%nfid to 0
       ! (could be done when defining iom_file in f95 but not in f90)
-#if ! defined key_agrif
-      IF( iom_open_init == 0 ) THEN
-         iom_file(:)%nfid = 0
-         iom_open_init = 1
-      ENDIF
-#else
       IF( Agrif_Root() ) THEN
          IF( iom_open_init == 0 ) THEN
@@ -176,5 +189,4 @@
          ENDIF
       ENDIF
-#endif
       ! do we read or write the file?
       IF( PRESENT(ldwrt) ) THEN   ;   llwrt = ldwrt
@@ -199,5 +211,4 @@
       ! =============
       clname   = trim(cdname)
-#if defined key_agrif
       IF ( .NOT. Agrif_Root() .AND. .NOT. lliof ) THEN
          iln    = INDEX(clname,'/') 
@@ -206,5 +217,4 @@
          clname=TRIM(cltmpn)//TRIM(Agrif_CFixed())//'_'//TRIM(clname)
       ENDIF
-#endif    
       ! which suffix should we use?
       SELECT CASE (iolib)
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/OBC/obc_oce.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/OBC/obc_oce.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/OBC/obc_oce.F90	(revision 2007)
@@ -27,5 +27,4 @@
    !
    !                                  !!* Namelist namobc: open boundary condition *
-   INTEGER           ::   nn_nbobc    = 2        !: number of open boundaries ( 1=< nbobc =< 4 ) 
    INTEGER           ::   nn_obcdta   = 0        !:  = 0 use the initial state as obc data
    !                                             !   = 1 read obc data in obcxxx.dta files
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/OBC/obc_par.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/OBC/obc_par.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/OBC/obc_par.F90	(revision 2007)
@@ -25,5 +25,9 @@
    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
    !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_obc = .TRUE.   !: Ocean Boundary Condition flag
+#if ! defined key_agrif
+   LOGICAL, PUBLIC, PARAMETER ::   lk_obc = .TRUE.     !: Ocean Boundary Condition flag
+#else
+   LOGICAL, PUBLIC            ::   lk_obc = .TRUE.     !: Ocean Boundary Condition flag
+#endif
 
 # if defined key_eel_r5
@@ -43,10 +47,10 @@
    !! open boundary parameter
    !!---------------------------------------------------------------------
-   INTEGER, PARAMETER ::     &  !: time dimension of the BCS fields on input
+   INTEGER ::     &  !: time dimension of the BCS fields on input
       jptobc  =	      2 
    !! * EAST open boundary
-   LOGICAL, PARAMETER ::     &  !:
+   LOGICAL ::     &  !:
       lp_obc_east = .FALSE.     !: to active or not the East open boundary
-   INTEGER, PARAMETER ::     &  !:
+   INTEGER ::     & 
       jpieob  = jpiglo-2,    &  !: i-localization of the East open boundary (must be ocean U-point)
       jpjed   =        2,    &  !: j-starting indice of the East open boundary (must be land T-point)
@@ -56,7 +60,7 @@
 
    !! * WEST open boundary
-   LOGICAL, PARAMETER ::     &  !:
+   LOGICAL ::     &  !:
       lp_obc_west = .FALSE.     !: to active or not the West open boundary
-   INTEGER, PARAMETER ::     &  !:
+   INTEGER ::     &  !:
       jpiwob  =	       2,    &  !: i-localization of the West open boundary (must be ocean U-point)
       jpjwd   =	       2,    &  !: j-starting indice of the West open boundary (must be land T-point)
@@ -66,7 +70,7 @@
 
    !! * NORTH open boundary
-   LOGICAL, PARAMETER ::     &  !:
+   LOGICAL ::     &  !:
       lp_obc_north = .FALSE.    !: to active or not the North open boundary
-   INTEGER, PARAMETER ::     &  !:
+   INTEGER ::     &  !:
       jpjnob  = jpjglo-2,    &  !: j-localization of the North open boundary (must be ocean V-point)
       jpind   =        2,    &  !: i-starting indice of the North open boundary (must be land T-point)
@@ -76,7 +80,7 @@
 
    !! * SOUTH open boundary
-   LOGICAL, PARAMETER ::     &  !:
+   LOGICAL ::     &  !:
       lp_obc_south = .FALSE.    !: to active or not the South open boundary
-   INTEGER, PARAMETER ::     &  !:
+   INTEGER ::     &  !:
       jpjsob  =        2,    &  !: j-localization of the South open boundary (must be ocean V-point)
       jpisd   =        2,    &  !: i-starting indice of the South open boundary (must be land T-point)
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/OBC/obc_par_EEL_R5.h90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/OBC/obc_par_EEL_R5.h90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/OBC/obc_par_EEL_R5.h90	(revision 2007)
@@ -15,5 +15,10 @@
    LOGICAL, PARAMETER ::     &  !:
       lp_obc_east = .TRUE.      !: to active or not the East open boundary
-   INTEGER, PARAMETER ::     &  !:
+
+     INTEGER   &
+#if !defined key_agrif
+     , PARAMETER   &
+#endif
+    ::     &
       jpieob  = jpiglo-2,    &  !: i-localization of the East open boundary (must be ocean U-point)
       jpjed   =        2,    &  !: j-starting indice of the East open boundary (must be land T-point)
@@ -25,5 +30,10 @@
    LOGICAL, PARAMETER ::     &  !:
       lp_obc_west = .TRUE.      !: to active or not the West open boundary
-   INTEGER, PARAMETER ::     &
+
+     INTEGER   &
+#if !defined key_agrif
+     , PARAMETER   &
+#endif
+    ::     &
       jpiwob  =	       2,    &  !: i-localization of the West open boundary (must be ocean U-point)
       jpjwd   =	       2,    &  !: j-starting indice of the West open boundary (must be land T-point)
@@ -35,5 +45,10 @@
    LOGICAL, PARAMETER ::     &  !:
       lp_obc_north = .FALSE.    !: to active or not the North open boundary
-   INTEGER, PARAMETER ::     &  !:
+
+     INTEGER   &
+#if !defined key_agrif
+     , PARAMETER   &
+#endif
+    ::     &
       jpjnob  = jpjglo-2,    &  !: j-localization of the North open boundary (must be ocean V-point)
       jpind   =        2,    &  !: i-starting indice of the North open boundary (must be land T-point)
@@ -45,5 +60,10 @@
    LOGICAL, PARAMETER ::     &  !:
       lp_obc_south = .FALSE.    !: to active or not the South open boundary
-   INTEGER, PARAMETER ::     &  !:
+
+     INTEGER   &
+#if !defined key_agrif
+     , PARAMETER   &
+#endif
+    ::     &
       jpjsob  =        2,    &  !: j-localization of the South open boundary (must be ocean V-point)
       jpisd   =        2,    &  !: i-starting indice of the South open boundary (must be land T-point)
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/OBC/obc_par_POMME_R025.h90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/OBC/obc_par_POMME_R025.h90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/OBC/obc_par_POMME_R025.h90	(revision 2007)
@@ -4,21 +4,23 @@
    !! open boundary parameter : POMME configuration
    !!---------------------------------------------------------------------
-     INTEGER, PARAMETER ::     &  !: time dimension of the BCS fields on input
-      jptobc  =	      14
+   INTEGER, PARAMETER ::  jptobc  =	      14   
+   !: time dimension of the BCS fields on input
 
    !! * EAST open boundary
    LOGICAL, PARAMETER ::     &  !:
       lp_obc_east = .TRUE.      !:
-   INTEGER, PARAMETER ::     &  !:
 
+     INTEGER   &
+#if !defined key_agrif
+     , PARAMETER   & 
+#endif
+    ::     & 
       ! * default values *
       !jpieob = jpiglo-2,    &  !: i-localization of the East open boundary (must be ocean U-point)
       !jpjed  =        2,    &  !: j-starting indice of the East open boundary (must be land T-point)
       !jpjef  = jpjglo-1,    &  !: j-ending   indice of the East open boundary (must be land T-point)
-
       jpieob = jpiglo-2,    &  !: i-localization of the East open boundary (must be ocean U-point)
       jpjed  =        1,    &  !: j-starting indice of the East open boundary (must be land T-point)
       jpjef  =   jpjglo,    &  !: j-ending   indice of the East open boundary (must be land T-point)
-
       jpjedp1 =  jpjed+1,    &  !: first ocean point         "                 "
       jpjefm1 =  jpjef-1        !: last  ocean point         "                 "
@@ -27,15 +29,17 @@
    LOGICAL, PARAMETER ::     &  !:
       lp_obc_west = .TRUE.     !: to active or not the West open boundary
-   INTEGER, PARAMETER ::     &  !:
 
+     INTEGER   &
+#if !defined key_agrif
+     , PARAMETER   & 
+#endif
+    ::     & 
       ! * default values *
       !jpiwob  =        2,   &  !: i-localization of the West open boundary (must be ocean U-point)
       !jpjwd   =        2,   &  !: j-starting indice of the West open boundary (must be land T-point)
       !jpjwf   = jpjglo-1,   &  !: j-ending   indice of the West open boundary (must be land T-point)
-
       jpiwob  =        2,   &  !: i-localization of the West open boundary (must be ocean U-point)
       jpjwd   =        1,   &  !: j-starting indice of the West open boundary (must be land T-point)
       jpjwf   =   jpjglo,   &  !: j-ending   indice of the West open boundary (must be land T-point)
-
       jpjwdp1 =  jpjwd+1,    &  !: first ocean point         "                 "
       jpjwfm1 =  jpjwf-1        !: last  ocean point         "                 "
@@ -44,15 +48,17 @@
    LOGICAL, PARAMETER ::     &  !:
       lp_obc_north = .TRUE.     !:
-   INTEGER, PARAMETER ::     &  !:
 
+     INTEGER   &
+#if !defined key_agrif
+     , PARAMETER   & 
+#endif
+    ::     & 
       ! * default values *
       !jpjnob = jpjglo-2,    &  !: j-localization of the North open boundary (must be ocean V-point)
       !jpind  =        2,    &  !: i-starting indice of the North open boundary (must be land T-point)
       !jpinf  = jpiglo-1,    &  !: i-ending   indice of the North open boundary (must be land T-point)
-
       jpjnob = jpjglo-2,    &  !: j-localization of the North open boundary (must be ocean V-point)
       jpind  =        1,    &  !: i-starting indice of the North open boundary (must be land T-point)
       jpinf  =   jpiglo,    &  !: i-ending   indice of the North open boundary (must be land T-point)
-
       jpindp1 =  jpind+1,    &  !: first ocean point         "                 "
       jpinfm1 =  jpinf-1        !: last  ocean point         "                 "
@@ -61,15 +67,17 @@
    LOGICAL, PARAMETER ::     &  !:
       lp_obc_south = .TRUE.     !: INDICE to active or not the South open boundary
-   INTEGER, PARAMETER ::     &  !:
 
+     INTEGER   &
+#if !defined key_agrif
+     , PARAMETER   & 
+#endif
+    ::     & 
       ! * default values *
       !jpjsob =        2,    &  !: j-localization of the South open boundary (must be ocean V-point)
       !jpisd  =        2,    &  !: i-starting indice of the South open boundary (must be land T-point)
       !jpisf  = jpiglo-1,    &  !: i-ending   indice of the South open boundary (must be land T-point)
-
       jpjsob =        2,    &  !: j-localization of the South open boundary (must be ocean V-point)
       jpisd  =        1,    &  !: i-starting indice of the South open boundary (must be land T-point)
       jpisf  =   jpiglo,    &  !: i-ending   indice of the South open boundary (must be land T-point)
-
       jpisdp1 =  jpisd+1,    &  !: first ocean point         "                 "
       jpisfm1 =  jpisf-1        !: last  ocean point         "                 "
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/OBC/obcdta.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/OBC/obcdta.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/OBC/obcdta.F90	(revision 2007)
@@ -469,5 +469,9 @@
        ENDIF
     ELSE
+#if defined key_agrif
+       IF ( ASSOCIATED(ztcobc) ) DEALLOCATE ( ztcobc )
+#else
        IF ( ALLOCATED(ztcobc) ) DEALLOCATE ( ztcobc )
+#endif
        ALLOCATE (ztcobc(itobc))
        DO ji=1,1   ! use a dummy loop to read ztcobc only once
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/OBC/obcfla.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/OBC/obcfla.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/OBC/obcfla.F90	(revision 2007)
@@ -82,10 +82,8 @@
 
       DO ji = nie0, nie1
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               ua_e(ji,jj) = (  ubtfoe(jj) + sqrt( grav*hu(ji,jj) )           &
-                  &            * ( ( sshn_e(ji,jj) + sshn_e(ji+1,jj) ) * 0.5  &
-                  &            - sshfoe(jj) )  ) * uemsk(jj,jk)
-            END DO
+         DO jj = 1, jpj
+            ua_e(ji,jj) = (  ubtfoe(jj) * hur(ji,jj) + sqrt( grav*hur(ji,jj) )   &
+               &            * ( ( sshn_e(ji,jj) + sshn_e(ji+1,jj) ) * 0.5        &
+               &            - sshfoe(jj) )  ) * uemsk(jj,jk)
          END DO
       END DO
@@ -97,5 +95,5 @@
             ssha_e(ji,jj) = ssha_e(ji,jj) * ( 1. - temsk(jj,1) ) &
                &            + temsk(jj,1) * sshfoe(jj)
-            va_e(ji,jj) = vbtfoe(jj) * uemsk(jj,jk)
+            va_e(ji,jj) = vbtfoe(jj) * hvr(ji,jj) * uemsk(jj,jk)
          END DO
       END DO
@@ -120,11 +118,9 @@
 
       DO ji = niw0, niw1
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               ua_e(ji,jj) = ( ubtfow(jj) - sqrt( grav * hu(ji,jj) )          &
-                  &            * ( ( sshn_e(ji,jj) + sshn_e(ji+1,jj) ) * 0.5  &
-                  &                - sshfow(jj) ) ) * uwmsk(jj,jk)
-               va_e(ji,jj) = vbtfow(jj) * uwmsk(jj,jk)
-            END DO
+         DO jj = 1, jpj
+            ua_e(ji,jj) = ( ubtfow(jj) * hur(ji,jj) - sqrt( grav * hur(ji,jj) )   &
+               &            * ( ( sshn_e(ji,jj) + sshn_e(ji+1,jj) ) * 0.5         &
+               &                - sshfow(jj) ) ) * uwmsk(jj,jk)
+            va_e(ji,jj) = vbtfow(jj) * hvr(ji,jj) * uwmsk(jj,jk)
          END DO
          DO jj = 1, jpj
@@ -155,10 +151,8 @@
 
       DO jj = njn0, njn1
-         DO jk = 1, jpkm1
-            DO ji = 1, jpi
-               va_e(ji,jj) = ( vbtfon(ji) + sqrt( grav * hv(ji,jj) )           &
-                  &            * ( ( sshn_e(ji,jj) + sshn_e(ji,jj+1) ) * 0.5   &
-                  &                - sshfon(ji) ) ) * vnmsk(ji,jk)
-            END DO
+         DO ji = 1, jpi
+            va_e(ji,jj) = ( vbtfon(ji) * hvr(ji,jj) + sqrt( grav * hvr(ji,jj) )   &
+               &            * ( ( sshn_e(ji,jj) + sshn_e(ji,jj+1) ) * 0.5         &
+               &                - sshfon(ji) ) ) * vnmsk(ji,jk)
          END DO
       END DO
@@ -170,5 +164,5 @@
             ssha_e(ji,jj) = ssha_e(ji,jj) * ( 1. - tnmsk(ji,1) ) &
                &            + sshfon(ji) * tnmsk(ji,1)
-            ua_e(ji,jj) = ubtfon(ji) * vnmsk(ji,jk)
+            ua_e(ji,jj) = ubtfon(ji) * hur(ji,jj) * vnmsk(ji,jk)
          END DO
       END DO
@@ -193,11 +187,9 @@
 
       DO jj = njs0, njs1
-         DO jk = 1, jpkm1
-            DO ji = 1, jpi
-               va_e(ji,jj) = ( vbtfos(ji) - sqrt( grav * hv(ji,jj) )            &
-                  &            * ( ( sshn_e(ji,jj) + sshn_e(ji,jj+1) ) * 0.5    &
-                  &                - sshfos(ji) ) ) * vsmsk(ji,jk)
-               ua_e(ji,jj) = ubtfos(ji) * vsmsk(ji,jk)
-            END DO
+         DO ji = 1, jpi
+            va_e(ji,jj) = ( vbtfos(ji) * hvr(ji,jj) - sqrt( grav * hvr(ji,jj) )   &
+               &            * ( ( sshn_e(ji,jj) + sshn_e(ji,jj+1) ) * 0.5         &
+               &                - sshfos(ji) ) ) * vsmsk(ji,jk)
+            ua_e(ji,jj) = ubtfos(ji) * hur(ji,jj) * vsmsk(ji,jk)
          END DO
          DO ji = 1, jpi
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/OBC/obcini.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/OBC/obcini.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/OBC/obcini.F90	(revision 2007)
@@ -70,5 +70,4 @@
 
       ! convert DOCTOR namelist name into the OLD names
-      nbobc    = nn_nbobc
       nobc_dta = nn_obcdta
       cffile   = cn_obcdta
@@ -101,5 +100,5 @@
       IF(lwp) WRITE(numout,*) 'obc_init : initialization of open boundaries'
       IF(lwp) WRITE(numout,*) '~~~~~~~~'
-      IF(lwp) WRITE(numout,*) '   Number of open boundaries    nn_nbobc = ', nn_nbobc
+      IF(lwp) WRITE(numout,*) '   Number of open boundaries    nbobc = ', nbobc
       IF(lwp) WRITE(numout,*)
 
@@ -306,4 +305,5 @@
       IF( lp_obc_east ) THEN
          !... (jpjed,jpjefm1),jpieob
+         bmask(nie0p1:nie1p1,nje0:nje1m1) = 0.e0
 
          ! ... initilization to zero
@@ -341,4 +341,5 @@
       IF( lp_obc_north ) THEN
          ! ... jpjnob,(jpind,jpisfm1)
+         bmask(nin0:nin1m1,njn0p1:njn1p1) = 0.e0
 
          ! ... initilization to zero
@@ -440,5 +441,5 @@
             END DO
          END IF
-  
+
          IF( lp_obc_north ) THEN ! ... North open boundary lateral surface
             DO jj = njn0, njn1
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/OBC/obcrst.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/OBC/obcrst.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/OBC/obcrst.F90	(revision 2007)
@@ -96,5 +96,5 @@
          ! -------------
 
-         CALL ctl_opn( inum, 'restart.obc.output', 'REPLACE', 'UNFORMATTED', 'DIRECT', nreclo, numout, lwp )
+         CALL ctl_opn( inum, 'restart.obc.output', 'UNKNOWN', 'UNFORMATTED', 'DIRECT', nreclo, numout, lwp )
  
          ! 1.2 Write header
@@ -322,5 +322,5 @@
       ! 0.1 Open files
       ! ---------------
-      CALL ctl_opn( inum, 'restart.obc', 'REPLACE', 'UNFORMATTED', 'DIRECT', nreclo, numout, lwp )
+      CALL ctl_opn( inum, 'restart.obc', 'UNKNOWN', 'UNFORMATTED', 'DIRECT', nreclo, numout, lwp )
 
       ! 1. Read
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90	(revision 2007)
@@ -184,5 +184,5 @@
                         &                    nday   + 1                   - nmonth_len(nmonth) * COUNT((/llnxtmth/)), llstop )
 
-                     IF( sd(jf)%num == 0 .AND. .NOT. llstop ) THEN    ! next year file does not exist
+                     IF( sd(jf)%num <= 0 .AND. .NOT. llstop ) THEN    ! next year file does not exist
                         CALL ctl_warn('next year/month/day file: '//TRIM(sd(jf)%clname)//     &
                                 &     ' not present -> back to current year/month/day')
@@ -368,5 +368,5 @@
          
          ! if previous year/month/day file does not exist, we switch to the current year/month/day
-         IF( llprev .AND. sdjf%num == 0 ) THEN
+         IF( llprev .AND. sdjf%num <= 0 ) THEN
             CALL ctl_warn( 'previous year/month/day file: '//TRIM(sdjf%clname)//' not present -> back to current year/month/day')
             ! we force to read the first record of the current year/month/day instead of last record of previous year/month/day
@@ -399,5 +399,5 @@
       ENDIF
 
-      IF( sdjf%num == 0 )   CALL fld_clopn( sdjf, nyear, nmonth, nday )   ! make sure current year/month/day file is opened
+      IF( sdjf%num <= 0 )   CALL fld_clopn( sdjf, nyear, nmonth, nday )   ! make sure current year/month/day file is opened
 
       sdjf%nswap_sec = nsec_year + nsec1jan000 - 1   ! force read/update the after data in the following part of fld_read 
@@ -815,6 +815,5 @@
             WRITE(aname,'(a3,i2.2)') 'src',jn
             data_tmp(:,:) = 0
-            CALL iom_get ( inum, jpdom_unknown, aname, data_tmp(1:nlci,1:nlcj), &
-                           kstart=(/nimpp,njmpp/), kcount=(/nlci,nlcj/) )
+            CALL iom_get ( inum, jpdom_data, aname, data_tmp(:,:) )
             data_src(:,:) = INT(data_tmp(:,:))
             ref_wgts(nxt_wgt)%data_jpj(:,:,jn) = 1 + (data_src(:,:)-1) / ref_wgts(nxt_wgt)%ddims(1)
@@ -825,15 +824,14 @@
             aname = ' '
             WRITE(aname,'(a3,i2.2)') 'wgt',jn
-            ref_wgts(nxt_wgt)%data_wgt(1:nlci,1:nlcj,jn) = 0.0
-            CALL iom_get ( inum, jpdom_unknown, aname, ref_wgts(nxt_wgt)%data_wgt(1:nlci,1:nlcj,jn), &
-                           kstart=(/nimpp,njmpp/), kcount=(/nlci,nlcj/) )
+            ref_wgts(nxt_wgt)%data_wgt(:,:,jn) = 0.0
+            CALL iom_get ( inum, jpdom_data, aname, ref_wgts(nxt_wgt)%data_wgt(:,:,jn) )
          END DO
          CALL iom_close (inum)
  
          ! find min and max indices in grid
-         ref_wgts(nxt_wgt)%botleft(1) = MINVAL(ref_wgts(nxt_wgt)%data_jpi(1:nlci,1:nlcj,:))
-         ref_wgts(nxt_wgt)%botleft(2) = MINVAL(ref_wgts(nxt_wgt)%data_jpj(1:nlci,1:nlcj,:))
-         ref_wgts(nxt_wgt)%topright(1) = MAXVAL(ref_wgts(nxt_wgt)%data_jpi(1:nlci,1:nlcj,:))
-         ref_wgts(nxt_wgt)%topright(2) = MAXVAL(ref_wgts(nxt_wgt)%data_jpj(1:nlci,1:nlcj,:))
+         ref_wgts(nxt_wgt)%botleft(1) = MINVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
+         ref_wgts(nxt_wgt)%botleft(2) = MINVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
+         ref_wgts(nxt_wgt)%topright(1) = MAXVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
+         ref_wgts(nxt_wgt)%topright(2) = MAXVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
 
          ! and therefore dimensions of the input box
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90	(revision 2007)
@@ -311,8 +311,8 @@
 
       ! lateral boundary cond.: T-, U-, V-, F-pts, sgn
-      CALL lbc_lnk( gcost, 'T', 1. )   ;   CALL lbc_lnk( gsint, 'T', -1. )
-      CALL lbc_lnk( gcosu, 'U', 1. )   ;   CALL lbc_lnk( gsinu, 'U', -1. )
-      CALL lbc_lnk( gcosv, 'V', 1. )   ;   CALL lbc_lnk( gsinv, 'V', -1. )
-      CALL lbc_lnk( gcosf, 'F', 1. )   ;   CALL lbc_lnk( gsinf, 'F', -1. )
+      CALL lbc_lnk( gcost, 'T', -1. )   ;   CALL lbc_lnk( gsint, 'T', -1. )
+      CALL lbc_lnk( gcosu, 'U', -1. )   ;   CALL lbc_lnk( gsinu, 'U', -1. )
+      CALL lbc_lnk( gcosv, 'V', -1. )   ;   CALL lbc_lnk( gsinv, 'V', -1. )
+      CALL lbc_lnk( gcosf, 'F', -1. )   ;   CALL lbc_lnk( gsinf, 'F', -1. )
 
    END SUBROUTINE angle
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 2007)
@@ -23,4 +23,5 @@
    USE sbc_oce         ! Surface boundary condition: ocean fields
    USE sbc_ice         ! Surface boundary condition: ice fields
+   USE phycst          ! physical constants
 #if defined key_lim3
    USE par_ice         ! ice parameters
@@ -45,5 +46,4 @@
    USE lib_mpp         ! distribued memory computing library
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE phycst, ONLY : xlsn, rhosn, xlic, rhoic
 #if defined key_cpl_carbon_cycle
    USE p4zflx, ONLY : oce_co2
@@ -274,5 +274,6 @@
       srcv(jpr_itz2)%clname = 'O_ITauz2'      ! 3rd   -      -         -     - 
       ! 
-      srcv(jpr_otx1:jpr_itz2)%nsgn = -1                           ! Vectors: change of sign at north fold
+      ! Vectors: change of sign at north fold ONLY if on the local grid
+      IF( TRIM( cn_rcv_tau(3) ) == 'local grid' )   srcv(jpr_otx1:jpr_itz2)%nsgn = -1.
       
       !                                                           ! Set grid and action
@@ -714,7 +715,6 @@
          IF( srcv(jpr_qnsoce)%laction )   qns(:,:) = frcv(:,:,jpr_qnsoce)
          IF( srcv(jpr_qnsmix)%laction )   qns(:,:) = frcv(:,:,jpr_qnsmix)        
-         !   energy for melting solid precipitation over free ocean
-         zcoef = xlsn / rhosn
-         qns(:,:) = qns(:,:) - frcv(:,:,jpr_snow) * zcoef
+         qns(:,:) = qns(:,:) - frcv(:,:,jpr_snow) * lfus              ! add the latent heat of solid precip. melting 
+
          !                                                       ! solar flux over the ocean          (qsr)
          IF( srcv(jpr_qsroce)%laction )   qsr(:,:) = frcv(:,:,jpr_qsroce) 
@@ -1117,9 +1117,6 @@
             &                                                   +          pist(:,:,1)   * zicefr(:,:,1) ) )
       END SELECT
-      !                                                           ! snow melting heat flux ....
-      !   energy for melting solid precipitation over ice-free ocean
-      zcoef = xlsn / rhosn
-      ztmp(:,:) = p_frld(:,:,1) * zsnow(:,:) * zcoef
-      pqns_tot(:,:) = pqns_tot(:,:) - ztmp(:,:)
+      ztmp(:,:) = p_frld(:,:,1) * zsnow(:,:) * lfus               ! add the latent heat of solid precip. melting
+      pqns_tot(:,:) = pqns_tot(:,:) - ztmp(:,:)                   ! over free ocean 
       IF( lk_diaar5 )   CALL iom_put( 'hflx_snow_cea', ztmp + zsnow(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average)
 !!gm
@@ -1130,9 +1127,7 @@
 !!
 !! similar job should be done for snow and precipitation temperature
-      !                                                           ! Iceberg melting heat flux ....
-      !   energy for iceberg melting 
-      IF( srcv(jpr_cal)%laction ) THEN 
-         zcoef = xlic / rhoic
-         ztmp(:,:) = frcv(:,:,jpr_cal) * zcoef
+      !                                     
+      IF( srcv(jpr_cal)%laction ) THEN                            ! Iceberg melting 
+         ztmp(:,:) = frcv(:,:,jpr_cal) * lfus                     ! add the latent heat of iceberg melting 
          pqns_tot(:,:) = pqns_tot(:,:) - ztmp(:,:)
          IF( lk_diaar5 )   CALL iom_put( 'hflx_cal_cea', ztmp + frcv(:,:,jpr_cal) * zcptn(:,:) )   ! heat flux from calving
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 2007)
@@ -65,5 +65,5 @@
       INTEGER  ::   inum                  ! temporary logical unit
       INTEGER  ::   ikty, iyear           ! 
-      REAL(wp) ::   z_emp, z_emp_nsrf       ! temporary scalars
+      REAL(wp) ::   z_emp, z_emp_nsrf, zsum_emp, zsum_erp       ! temporary scalars
       REAL(wp) ::   zsurf_neg, zsurf_pos, zsurf_tospread
       REAL(wp), DIMENSION(jpi,jpj) ::   ztmsk_neg, ztmsk_pos, ztmsk_tospread
@@ -165,4 +165,6 @@
             !
             IF( lk_mpp )   CALL  mpp_sum( z_emp )
+            IF( lk_mpp )   CALL  mpp_sum( zsurf_neg )
+            IF( lk_mpp )   CALL  mpp_sum( zsurf_pos )
             
             IF( z_emp < 0.e0 ) THEN
@@ -177,7 +179,12 @@
 
             ! emp global mean over <0 or >0 erp area
-            z_emp_nsrf = SUM( e1e2_i(:,:) * z_emp ) / ( zsurf_tospread + rsmall )
+            zsum_emp = SUM( e1e2_i(:,:) * z_emp )
+            IF( lk_mpp )   CALL  mpp_sum( zsum_emp )
+            z_emp_nsrf =  zsum_emp / ( zsurf_tospread + rsmall )
             ! weight to respect erp field 2D structure 
-            z_wgt(:,:) = ztmsk_tospread(:,:) * erp(:,:) / ( SUM( ztmsk_tospread(:,:) * erp(:,:) * e1e2_i(:,:) ) + rsmall )
+            zsum_erp = SUM( ztmsk_tospread(:,:) * erp(:,:) * e1e2_i(:,:) )
+            IF( lk_mpp )   CALL  mpp_sum( zsum_erp )
+            z_wgt(:,:) = ztmsk_tospread(:,:) * erp(:,:) / ( zsum_erp + rsmall )
+
             ! final correction term to apply
             zerp_cor(:,:) = -1. * z_emp_nsrf * zsurf_tospread * z_wgt(:,:)
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 2007)
@@ -85,12 +85,10 @@
 !!gm here no overwrite, test all option via namelist change: require more incore memory
 !!gm  IF( lk_sbc_cpl       ) THEN   ;   ln_cpl      = .TRUE.   ;   ELSE   ;   ln_cpl      = .FALSE.   ;   ENDIF
-#if defined key_agrif
+
       IF ( Agrif_Root() ) THEN
-#endif
         IF( lk_lim2 )            nn_ice      = 2
         IF( lk_lim3 )            nn_ice      = 3
-#if defined key_agrif
-      ENDIF
-#endif
+      ENDIF
+      !
       IF( cp_cfg == 'gyre' ) THEN
           ln_ana      = .TRUE.   
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/SOL/solmat.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/SOL/solmat.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/SOL/solmat.F90	(revision 2007)
@@ -99,45 +99,46 @@
       
 #  elif defined key_dynspg_flt && defined key_obc
-
-      DO jj = 2, jpjm1                      ! matrix of free surface elliptic system with open boundaries
-         DO ji = 2, jpim1
-            zcoef = z2dt * z2dt * grav * bmask(ji,jj)
-            !                                    ! south coefficient
-            IF( lp_obc_south .AND. ( jj == njs0p1 ) ) THEN
-               zcoefs = -zcoef * hv(ji,jj-1) * e1v(ji,jj-1)/e2v(ji,jj-1)*(1.-vsmsk(ji,1))
-            ELSE
-               zcoefs = -zcoef * hv(ji,jj-1) * e1v(ji,jj-1)/e2v(ji,jj-1)
-            END IF
-            gcp(ji,jj,1) = zcoefs
-            !
-            !                                    ! west coefficient
-            IF( lp_obc_west  .AND. ( ji == niw0p1 ) ) THEN
-               zcoefw = -zcoef * hu(ji-1,jj) * e2u(ji-1,jj)/e1u(ji-1,jj)*(1.-uwmsk(jj,1))
-            ELSE
-               zcoefw = -zcoef * hu(ji-1,jj) * e2u(ji-1,jj)/e1u(ji-1,jj)
-            END IF
-            gcp(ji,jj,2) = zcoefw
-            !
-            !                                    ! east coefficient
-            IF( lp_obc_east  .AND. ( ji == nie0 ) ) THEN
-               zcoefe = -zcoef * hu(ji,jj) * e2u(ji,jj)/e1u(ji,jj)*(1.-uemsk(jj,1))
-            ELSE
-               zcoefe = -zcoef * hu(ji,jj) * e2u(ji,jj)/e1u(ji,jj)
-            END IF
-            gcp(ji,jj,3) = zcoefe
-            !
-            !                                    ! north coefficient
-            IF( lp_obc_north .AND. ( jj == njn0 ) ) THEN
-               zcoefn = -zcoef * hv(ji,jj) * e1v(ji,jj)/e2v(ji,jj)*(1.-vnmsk(ji,1))
-            ELSE
+      IF( Agrif_Root() ) THEN
+         DO jj = 2, jpjm1                      ! matrix of free surface elliptic system with open boundaries
+            DO ji = 2, jpim1
+               zcoef = z2dt * z2dt * grav * bmask(ji,jj)
+               !                                    ! south coefficient
+               IF( lp_obc_south .AND. ( jj == njs0p1 ) ) THEN
+                  zcoefs = -zcoef * hv(ji,jj-1) * e1v(ji,jj-1)/e2v(ji,jj-1)*(1.-vsmsk(ji,1))
+               ELSE
+                  zcoefs = -zcoef * hv(ji,jj-1) * e1v(ji,jj-1)/e2v(ji,jj-1)
+               END IF
+               gcp(ji,jj,1) = zcoefs
+               !
+               !                                    ! west coefficient
+               IF( lp_obc_west  .AND. ( ji == niw0p1 ) ) THEN
+                  zcoefw = -zcoef * hu(ji-1,jj) * e2u(ji-1,jj)/e1u(ji-1,jj)*(1.-uwmsk(jj,1))
+               ELSE
+                  zcoefw = -zcoef * hu(ji-1,jj) * e2u(ji-1,jj)/e1u(ji-1,jj)
+               END IF
+               gcp(ji,jj,2) = zcoefw
+               !
+               !                                    ! east coefficient
+               IF( lp_obc_east  .AND. ( ji == nie0 ) ) THEN
+                  zcoefe = -zcoef * hu(ji,jj) * e2u(ji,jj)/e1u(ji,jj)*(1.-uemsk(jj,1))
+               ELSE
+                  zcoefe = -zcoef * hu(ji,jj) * e2u(ji,jj)/e1u(ji,jj)
+               END IF
+               gcp(ji,jj,3) = zcoefe
+               !
+               !                                    ! north coefficient
+               IF( lp_obc_north .AND. ( jj == njn0 ) ) THEN
+                  zcoefn = -zcoef * hv(ji,jj) * e1v(ji,jj)/e2v(ji,jj)*(1.-vnmsk(ji,1))
+               ELSE
                zcoefn = -zcoef * hv(ji,jj) * e1v(ji,jj)/e2v(ji,jj)
-            END IF
-            gcp(ji,jj,4) = zcoefn
-            !
-            !                                    ! diagonal coefficient
-            gcdmat(ji,jj) = e1t(ji,jj)*e2t(ji,jj)*bmask(ji,jj)   &
-               &            - zcoefs -zcoefw -zcoefe -zcoefn
+               END IF
+               gcp(ji,jj,4) = zcoefn
+               !
+               !                                    ! diagonal coefficient
+               gcdmat(ji,jj) = e1t(ji,jj)*e2t(ji,jj)*bmask(ji,jj)   &
+                  &            - zcoefs -zcoefw -zcoefe -zcoefn
+            END DO
          END DO
-      END DO
+      ENDIF
 #endif
 
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 2007)
@@ -179,4 +179,9 @@
       END DO
 
+      ! "zonal" mean advective heat and salt transport
+      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
+         pht_adv(:) = ptr_vj( ztv(:,:,:) )
+         pst_adv(:) = ptr_vj( zsv(:,:,:) )
+      ENDIF
 
       ! Save the intermediate i / j / k advective trends for diagnostics
@@ -366,6 +371,6 @@
       ! "zonal" mean advective heat and salt transport
       IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
-         pht_adv(:) = ptr_vj( ztv(:,:,:) )
-         pst_adv(:) = ptr_vj( zsv(:,:,:) )
+         pht_adv(:) = ptr_vj( ztv(:,:,:) ) + pht_adv(:)
+         pst_adv(:) = ptr_vj( zsv(:,:,:) ) + pst_adv(:)
       ENDIF
       !
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 2007)
@@ -38,4 +38,5 @@
    USE agrif_opa_update
    USE agrif_opa_interp
+   USE obc_oce 
 
    IMPLICIT NONE
@@ -101,5 +102,5 @@
       !
 #if defined key_obc
-      CALL obc_tra( kt )               ! OBC open boundaries
+      IF( lk_obc )   CALL obc_tra( kt )  ! OBC open boundaries
 #endif
 #if defined key_bdy
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 2007)
@@ -134,5 +134,7 @@
                zta = ro0cpr * qns(ji,jj) * zse3t &                   ! temperature : heat flux
                 &    - emp(ji,jj) * zsrau * tn(ji,jj,1)  * zse3t     ! & cooling/heating effet of EMP flux
-               zsa = 0.e0                                            ! No salinity concent./dilut. effect
+               zsa = ( emps(ji,jj) - emp(ji,jj) ) &
+                &                 * zsrau * sn(ji,jj,1)  * zse3t     ! concent./dilut. effect due to sea-ice 
+                                                                     ! melt/formation and (possibly) SSS restoration
             ELSE
                zta = ro0cpr * qns(ji,jj) * zse3t     ! temperature : heat flux
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/lib_mpp.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/lib_mpp.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/lib_mpp.F90	(revision 2007)
@@ -103,5 +103,5 @@
    !! ========================= !!
 !$AGRIF_DO_NOT_TREAT
-   INCLUDE 'mpif.h'
+   INCLUDE mpif.h
 !$AGRIF_END_DO_NOT_TREAT
    
@@ -112,5 +112,7 @@
    INTEGER ::   mppsize        ! number of process
    INTEGER ::   mpprank        ! process number  [ 0 - size-1 ]
-   INTEGER ::   mpi_comm_opa   ! opa local communicator
+!$AGRIF_DO_NOT_TREAT
+   INTEGER, PUBLIC ::   mpi_comm_opa   ! opa local communicator
+!$AGRIF_END_DO_NOT_TREAT
 
    ! variables used in case of sea-ice
@@ -191,67 +193,13 @@
       WRITE(ldtxt(6),*) '      size in bytes of exported buffer   nn_buffer   = ', nn_buffer
 
-#if defined key_agrif
-      IF( Agrif_Root() ) THEN
-#endif
-         !!bug RB : should be clean to use Agrif in coupled mode
-#if ! defined key_agrif
-         CALL mpi_initialized ( mpi_was_called, code )
-         IF( code /= MPI_SUCCESS ) THEN
-            WRITE(*, cform_err)
-            WRITE(*, *) 'lib_mpp: Error in routine mpi_initialized'
-            CALL mpi_abort( mpi_comm_world, code, ierr )
-         ENDIF
-
-         IF( PRESENT(localComm) .and. mpi_was_called ) THEN
-            mpi_comm_opa = localComm
-            SELECT CASE ( cn_mpi_send )
-            CASE ( 'S' )                ! Standard mpi send (blocking)
-               WRITE(ldtxt(7),*) '           Standard blocking mpi send (send)'
-            CASE ( 'B' )                ! Buffer mpi send (blocking)
-               WRITE(ldtxt(7),*) '           Buffer blocking mpi send (bsend)'
-               CALL mpi_init_opa( ierr ) 
-            CASE ( 'I' )                ! Immediate mpi send (non-blocking send)
-               WRITE(ldtxt(7),*) '           Immediate non-blocking send (isend)'
-               l_isend = .TRUE.
-            CASE DEFAULT
-               WRITE(ldtxt(7),cform_err)
-               WRITE(ldtxt(8),*) '           bad value for cn_mpi_send = ', cn_mpi_send
-               nstop = nstop + 1
-            END SELECT
-         ELSE IF ( PRESENT(localComm) .and. .not. mpi_was_called ) THEN
-            WRITE(ldtxt(7),*) ' lib_mpp: You cannot provide a local communicator '
-            WRITE(ldtxt(8),*) '          without calling MPI_Init before ! '
-            nstop = nstop + 1
-         ELSE
-#endif
-            SELECT CASE ( cn_mpi_send )
-            CASE ( 'S' )                ! Standard mpi send (blocking)
-               WRITE(ldtxt(7),*) '           Standard blocking mpi send (send)'
-               CALL mpi_init( ierr )
-            CASE ( 'B' )                ! Buffer mpi send (blocking)
-               WRITE(ldtxt(7),*) '           Buffer blocking mpi send (bsend)'
-               CALL mpi_init_opa( ierr )
-            CASE ( 'I' )                ! Immediate mpi send (non-blocking send)
-               WRITE(ldtxt(7),*) '           Immediate non-blocking send (isend)'
-               l_isend = .TRUE.
-               CALL mpi_init( ierr )
-            CASE DEFAULT
-               WRITE(ldtxt(7),cform_err)
-               WRITE(ldtxt(8),*) '           bad value for cn_mpi_send = ', cn_mpi_send
-               nstop = nstop + 1
-            END SELECT
-
-#if ! defined key_agrif
-            CALL mpi_comm_dup( mpi_comm_world, mpi_comm_opa, code)
-            IF( code /= MPI_SUCCESS ) THEN
-               WRITE(*, cform_err)
-               WRITE(*, *) ' lib_mpp: Error in routine mpi_comm_dup'
-               CALL mpi_abort( mpi_comm_world, code, ierr )
-            ENDIF
-            !
-         ENDIF
-#endif
-#if defined key_agrif
-      ELSE
+      CALL mpi_initialized ( mpi_was_called, code )
+      IF( code /= MPI_SUCCESS ) THEN
+         WRITE(*, cform_err)
+         WRITE(*, *) 'lib_mpp: Error in routine mpi_initialized'
+         CALL mpi_abort( mpi_comm_world, code, ierr )
+      ENDIF
+
+      IF( mpi_was_called ) THEN
+         !
          SELECT CASE ( cn_mpi_send )
          CASE ( 'S' )                ! Standard mpi send (blocking)
@@ -259,4 +207,5 @@
          CASE ( 'B' )                ! Buffer mpi send (blocking)
             WRITE(ldtxt(7),*) '           Buffer blocking mpi send (bsend)'
+            CALL mpi_init_opa( ierr ) 
          CASE ( 'I' )                ! Immediate mpi send (non-blocking send)
             WRITE(ldtxt(7),*) '           Immediate non-blocking send (isend)'
@@ -267,8 +216,41 @@
             nstop = nstop + 1
          END SELECT
+      ELSE IF ( PRESENT(localComm) .and. .not. mpi_was_called ) THEN
+         WRITE(ldtxt(7),*) ' lib_mpp: You cannot provide a local communicator '
+         WRITE(ldtxt(8),*) '          without calling MPI_Init before ! '
+         nstop = nstop + 1
+      ELSE
+         SELECT CASE ( cn_mpi_send )
+         CASE ( 'S' )                ! Standard mpi send (blocking)
+            WRITE(ldtxt(7),*) '           Standard blocking mpi send (send)'
+            CALL mpi_init( ierr )
+         CASE ( 'B' )                ! Buffer mpi send (blocking)
+            WRITE(ldtxt(7),*) '           Buffer blocking mpi send (bsend)'
+            CALL mpi_init_opa( ierr )
+         CASE ( 'I' )                ! Immediate mpi send (non-blocking send)
+            WRITE(ldtxt(7),*) '           Immediate non-blocking send (isend)'
+            l_isend = .TRUE.
+            CALL mpi_init( ierr )
+         CASE DEFAULT
+            WRITE(ldtxt(7),cform_err)
+            WRITE(ldtxt(8),*) '           bad value for cn_mpi_send = ', cn_mpi_send
+            nstop = nstop + 1
+         END SELECT
+         !
       ENDIF
 
-      mpi_comm_opa = mpi_comm_world
-#endif
+      IF( PRESENT(localComm) ) THEN 
+         IF( Agrif_Root() ) THEN
+            mpi_comm_opa = localComm
+         ENDIF
+      ELSE
+         CALL mpi_comm_dup( mpi_comm_world, mpi_comm_opa, code)
+         IF( code /= MPI_SUCCESS ) THEN
+            WRITE(*, cform_err)
+            WRITE(*, *) ' lib_mpp: Error in routine mpi_comm_dup'
+            CALL mpi_abort( mpi_comm_world, code, ierr )
+         ENDIF
+      ENDIF 
+
       CALL mpi_comm_rank( mpi_comm_opa, mpprank, ierr )
       CALL mpi_comm_size( mpi_comm_opa, mppsize, ierr )
@@ -2067,4 +2049,5 @@
       ijpj   = 4
       ijpjm1 = 3
+      ztab(:,:,:) = 0.e0
       !
       DO jj = nlcj - ijpj +1, nlcj          ! put in znorthloc the last 4 jlines of pt3d
@@ -2132,4 +2115,5 @@
       ijpj   = 4
       ijpjm1 = 3
+      ztab(:,:) = 0.e0
       !
       DO jj = nlcj-ijpj+1, nlcj             ! put in znorthloc the last 4 jlines of pt2d
@@ -2197,4 +2181,5 @@
       !
       ijpj=4
+      ztab(:,:) = 0.e0
 
       ij=0
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/opa.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/opa.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/opa.F90	(revision 2007)
@@ -156,5 +156,7 @@
       CALL opa_closefile
 #if defined key_oasis3 || defined key_oasis4
-      CALL cpl_prism_finalize           ! end coupling and mpp communications with OASIS
+      IF( Agrif_Root() ) THEN
+         CALL cpl_prism_finalize           ! end coupling and mpp communications with OASIS
+     ENDIF 
 #else
       IF( lk_mpp )   CALL mppstop       ! end mpp communications
@@ -191,8 +193,12 @@
 #if defined key_iomput
 # if defined key_oasis3 || defined key_oasis4
-      CALL cpl_prism_init( ilocal_comm )      ! nemo local communicator given by oasis
-      CALL init_ioclient()                    ! io_server will get its communicators (if needed) from oasis (we don't see it)
+      IF( Agrif_Root() ) THEN
+         CALL cpl_prism_init( ilocal_comm )      ! nemo local communicator given by oasis
+         CALL init_ioclient()                    ! io_server will get its communicators (if needed) from oasis (we don't see it)
+      ENDIF
 # else
-      CALL init_ioclient( ilocal_comm )       ! nemo local communicator (used or not) given by the io_server
+      IF( Agrif_Root() ) THEN
+         CALL init_ioclient( ilocal_comm )       ! nemo local communicator (used or not) given by the io_server
+      ENDIF
 # endif
       narea = mynode( cltxt, ilocal_comm )    ! Nodes selection
@@ -200,5 +206,7 @@
 #else
 # if defined key_oasis3 || defined key_oasis4
-      CALL cpl_prism_init( ilocal_comm )      ! nemo local communicator given by oasis
+      IF( Agrif_Root() ) THEN
+         CALL cpl_prism_init( ilocal_comm )      ! nemo local communicator given by oasis
+      ENDIF
       narea = mynode( cltxt, ilocal_comm )    ! Nodes selection (control print return in cltxt)
 # else
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/par_POMME_R025.h90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/par_POMME_R025.h90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/par_POMME_R025.h90	(revision 2007)
@@ -22,9 +22,7 @@
       jp_cfg = 025  ,        &  !: resolution of the configuration (degrees)
       ! Original data size
-
       ! ORCA025 global grid size
       jpiglo_ORCA025 = 1442, &
       jpjglo_ORCA025 = 1021, &  ! not used currently
-
       ! POMME "global" domain localisation in the ORCA025 global grid
       jpi_iw    = 1059,      & 
@@ -32,9 +30,7 @@
       jpj_js    = 661,       &
       jpj_jn    = 700,       &
-
       jpidta  = ( jpi_ie - jpi_iw + 1 ), &   !: =30 first horizontal dimension > or = to jpi
       jpjdta  = ( jpj_jn - jpj_js + 1 ), &   !: =40 second                     > or = to jpj
       jpkdta  = 46 ,         &  !: number of levels           > or = to jpk
-
       ! total domain matrix size
       jpiglo  = jpidta,      &  !: first  dimension of global domain --> i
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 2007)
@@ -166,6 +166,9 @@
 #if defined key_agrif
       kstp = nit000 + Agrif_Nb_Step()
-!      IF ( Agrif_Root() .and. lwp) Write(*,*) '---'
-!      IF (lwp) Write(*,*) 'Grid Number',Agrif_Fixed(),' time step ',kstp
+!      IF( Agrif_Root() .and. lwp) Write(*,*) '---'
+!      IF(lwp) Write(*,*) 'Grid Number',Agrif_Fixed(),' time step ',kstp
+# if defined key_iomput
+      IF( Agrif_Nbstepint() == 0) CALL iom_swap
+# endif   
 #endif   
       indic = 1                                       ! reset to no error condition
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/trc_oce.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/trc_oce.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/trc_oce.F90	(revision 2007)
@@ -126,11 +126,11 @@
       zrgb(1,51) =  3.162   ;   zrgb(2,51) = 0.22703   ;   zrgb(3,51) = 0.16599   ;   zrgb(4,51) = 0.46601
       zrgb(1,52) =  3.548   ;   zrgb(2,52) = 0.24433   ;   zrgb(3,52) = 0.17334   ;   zrgb(4,52) = 0.47313
-      zrgb(1,53) =  3.981   ;   zrgb(2,53) = 0.26301   ;   zrgb(3,53) = 0.18126   ;   zrgb(4,54) = 0.48080
-      zrgb(1,54) =  4.467   ;   zrgb(2,54) = 0.28320   ;   zrgb(3,54) = 0.18981   ;   zrgb(4,55) = 0.48909
-      zrgb(1,55) =  5.012   ;   zrgb(2,55) = 0.30502   ;   zrgb(3,55) = 0.19903   ;   zrgb(4,56) = 0.49803
-      zrgb(1,56) =  5.623   ;   zrgb(2,56) = 0.32858   ;   zrgb(3,56) = 0.20898   ;   zrgb(4,57) = 0.50768
-      zrgb(1,57) =  6.310   ;   zrgb(2,57) = 0.35404   ;   zrgb(3,57) = 0.21971   ;   zrgb(4,58) = 0.51810
-      zrgb(1,58) =  7.079   ;   zrgb(2,58) = 0.38154   ;   zrgb(3,58) = 0.23129   ;   zrgb(4,59) = 0.52934
-      zrgb(1,59) =  7.943   ;   zrgb(2,59) = 0.41125   ;   zrgb(3,59) = 0.24378   ;   zrgb(4,50) = 0.54147
+      zrgb(1,53) =  3.981   ;   zrgb(2,53) = 0.26301   ;   zrgb(3,53) = 0.18126   ;   zrgb(4,53) = 0.48080
+      zrgb(1,54) =  4.467   ;   zrgb(2,54) = 0.28320   ;   zrgb(3,54) = 0.18981   ;   zrgb(4,54) = 0.48909
+      zrgb(1,55) =  5.012   ;   zrgb(2,55) = 0.30502   ;   zrgb(3,55) = 0.19903   ;   zrgb(4,55) = 0.49803
+      zrgb(1,56) =  5.623   ;   zrgb(2,56) = 0.32858   ;   zrgb(3,56) = 0.20898   ;   zrgb(4,56) = 0.50768
+      zrgb(1,57) =  6.310   ;   zrgb(2,57) = 0.35404   ;   zrgb(3,57) = 0.21971   ;   zrgb(4,57) = 0.51810
+      zrgb(1,58) =  7.079   ;   zrgb(2,58) = 0.38154   ;   zrgb(3,58) = 0.23129   ;   zrgb(4,58) = 0.52934
+      zrgb(1,59) =  7.943   ;   zrgb(2,59) = 0.41125   ;   zrgb(3,59) = 0.24378   ;   zrgb(4,59) = 0.54147
       zrgb(1,60) =  8.912   ;   zrgb(2,60) = 0.44336   ;   zrgb(3,60) = 0.25725   ;   zrgb(4,60) = 0.55457
       zrgb(1,61) = 10.000   ;   zrgb(2,61) = 0.47804   ;   zrgb(3,61) = 0.27178   ;   zrgb(4,61) = 0.56870
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/C14b/trclsm_c14b.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/C14b/trclsm_c14b.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/C14b/trclsm_c14b.F90	(revision 2007)
@@ -44,5 +44,5 @@
       INTEGER ::   numnatb
 
-#if defined key_trc_diaadd
+#if defined key_trc_diaadd && ! defined key_iomput
       ! definition of additional diagnostic as a structure
       INTEGER ::   jl, jn
@@ -58,5 +58,5 @@
       !!
       NAMELIST/namc14date/ ndate_beg_b, nyear_res_b
-#if defined key_trc_diaadd
+#if defined key_trc_diaadd && ! defined key_iomput
       NAMELIST/namc14dia/nwritedia, c14dia2d, c14dia3d     ! additional diagnostics
 #endif
@@ -81,5 +81,5 @@
       IF(lwp) WRITE(numout,*) '    initial year (aa)                  nyear_beg_b = ', nyear_beg_b
       !
-#if defined key_trc_diaadd
+#if defined key_trc_diaadd && ! defined key_iomput
 
       ! Namelist namc14dia
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/C14b/trcrst_c14b.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/C14b/trcrst_c14b.F90	(revision 2007)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/C14b/trcrst_c14b.F90	(revision 2007)
@@ -0,0 +1,91 @@
+MODULE trcrst_c14b
+   !!======================================================================
+   !!                       ***  MODULE trcrst_c14b  ***
+   !! TOP :   create, write, read the restart files of c14b tracer
+   !!======================================================================
+   !! History :   1.0  !  2010-01 (C. Ethe) Original
+   !!----------------------------------------------------------------------
+#if defined key_c14b
+   !!----------------------------------------------------------------------
+   !!   'key_c14b'                                               c14b tracers
+   !!----------------------------------------------------------------------
+   !!   trc_rst_read_c14b   : read  restart file
+   !!   trc_rst_wri_c14b    : write restart file
+   !!----------------------------------------------------------------------
+   USE oce_trc         ! Ocean variables
+   USE par_trc         ! TOP parameters
+   USE trc             ! TOP variables
+   USE trcsms_c14b          ! c14b sms trends
+   USE in_out_manager  ! I/O manager
+   USE iom
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC  trc_rst_read_c14b   ! called by trcini.F90 module
+   PUBLIC  trc_rst_wri_c14b   ! called by trcini.F90 module
+
+CONTAINS
+   
+   SUBROUTINE trc_rst_read_c14b( knum ) 
+      !!----------------------------------------------------------------------
+      !!                     ***  trc_rst_read_c14b  ***  
+      !!
+      !! ** Purpose : Read in restart file specific variables from c14b model
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in)  :: knum  ! unit of the restart file
+      INTEGER              :: jn    ! dummy loop indices
+      !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' trc_rst_read_c14b : Read specific variables from c14b model '
+      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
+      
+      DO jn = jp_c14b0, jp_c14b1
+         CALL iom_get( knum, jpdom_autoglo, 'qint_'//ctrcnm(jn), qint_c14(:,:,jn) ) 
+      END DO
+
+   END SUBROUTINE trc_rst_read_c14b
+
+   SUBROUTINE trc_rst_wri_c14b( kt, kitrst, knum )
+      !!----------------------------------------------------------------------
+      !!                     ***  trc_rst_read_c14b  ***
+      !!
+      !! ** Purpose : Read in restart file specific variables from c14b model
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in)  :: kt      ! time step
+      INTEGER, INTENT(in)  :: kitrst  ! time step of restart write
+      INTEGER, INTENT(in)  :: knum    ! unit of the restart file
+      INTEGER              :: jn      ! dummy loop indices
+      !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' trc_rst_wri_c14b : Write specific variables from c14b model '
+      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
+
+      DO jn = jp_c14b0, jp_c14b1
+         CALL iom_rstput( kt, kitrst, kum, 'qint_'//ctrcnm(jn), qint_c14(:,:,jn) )
+      END DO
+
+   END SUBROUTINE trc_rst_wri_c14b
+
+#else
+   !!----------------------------------------------------------------------
+   !!  Dummy module :                                     No passive tracer
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_rst_read_c14b( knum )
+     INTEGER, INTENT(in)  :: knum
+     WRITE(*,*) 'trc_rst_wri_c14b: You should not have seen this print! error?', knum
+   END SUBROUTINE trc_rst_read_c14b
+
+   SUBROUTINE trc_rst_wri_c14b( kt, kitrst, knum )
+     INTEGER, INTENT(in)  :: kt, kitrst, knum
+     WRITE(*,*) 'trc_rst_wri_c14b: You should not have seen this print! error?', kt, kitrst, knum
+   END SUBROUTINE trc_rst_wri_c14b
+#endif
+
+   !!======================================================================
+END MODULE trcrst_c14b
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/CFC/trcctl_cfc.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/CFC/trcctl_cfc.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/CFC/trcctl_cfc.F90	(revision 2007)
@@ -44,6 +44,5 @@
       IF( jp_cfc > 2) THEN 
           IF(lwp) THEN 
-              WRITE (numout,*) ' ===>>>> : w a r n i n g '
-              WRITE (numout,*) ' =======   ============= '
+              WRITE (numout,*) ' ===>>>> : w a r n i n g <<<<==='
               WRITE (numout,*)                             &
               &   ' STOP, change jp_cfc to 1 or 2 in par_CFC module '  
@@ -62,6 +61,5 @@
 
       IF(lwp) THEN
-         WRITE (numout,*) ' ===>>>> : w a r n i n g '
-         WRITE (numout,*) ' =======   ============= '
+         WRITE (numout,*) ' ===>>>> : w a r n i n g <<<<==='
          WRITE (numout,*) ' we force tracer names'
          DO jl = 1, jp_cfc
@@ -80,6 +78,5 @@
             ctrcun(jn) = 'mole/m3'
             IF(lwp) THEN
-               WRITE (numout,*) ' ===>>>> : w a r n i n g '
-               WRITE (numout,*) ' =======   ============= '
+               WRITE (numout,*) ' ===>>>> : w a r n i n g <<<<==='
                WRITE (numout,*) ' we force tracer unit'
                WRITE(numout,*) ' tracer  ',ctrcnm(jn), 'UNIT= ',ctrcun(jn)
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/CFC/trclsm_cfc.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/CFC/trclsm_cfc.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/CFC/trclsm_cfc.F90	(revision 2007)
@@ -43,5 +43,5 @@
       !!----------------------------------------------------------------------
       INTEGER ::   numnatc
-#if defined key_trc_diaadd
+#if defined key_trc_diaadd && ! defined key_iomput
       ! definition of additional diagnostic as a structure
       INTEGER :: jl, jn
@@ -56,5 +56,5 @@
       !!
       NAMELIST/namcfcdate/ ndate_beg, nyear_res
-#if defined key_trc_diaadd
+#if defined key_trc_diaadd && ! defined key_iomput
       NAMELIST/namcfcdia/nwritedia, cfcdia2d     ! additional diagnostics
 #endif
@@ -79,5 +79,5 @@
       IF(lwp) WRITE(numout,*) '    initial year (aa)                       nyear_beg = ', nyear_beg
       !
-#if defined key_trc_diaadd
+#if defined key_trc_diaadd && ! defined key_iomput
 
       ! Namelist namcfcdia
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/CFC/trcrst_cfc.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/CFC/trcrst_cfc.F90	(revision 2007)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/CFC/trcrst_cfc.F90	(revision 2007)
@@ -0,0 +1,91 @@
+MODULE trcrst_cfc
+   !!======================================================================
+   !!                       ***  MODULE trcrst_cfc  ***
+   !! TOP :   create, write, read the restart files of CFC tracer
+   !!======================================================================
+   !! History :   1.0  !  2010-01 (C. Ethe) Original
+   !!----------------------------------------------------------------------
+#if defined key_cfc
+   !!----------------------------------------------------------------------
+   !!   'key_cfc'                                               CFC tracers
+   !!----------------------------------------------------------------------
+   !!   trc_rst_read_cfc   : read  restart file
+   !!   trc_rst_wri_cfc    : write restart file
+   !!----------------------------------------------------------------------
+   USE oce_trc         ! Ocean variables
+   USE par_trc         ! TOP parameters
+   USE trc             ! TOP variables
+   USE trcsms_cfc          ! CFC sms trends
+   USE in_out_manager  ! I/O manager
+   USE iom
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC  trc_rst_read_cfc   ! called by trcini.F90 module
+   PUBLIC  trc_rst_wri_cfc   ! called by trcini.F90 module
+
+CONTAINS
+   
+   SUBROUTINE trc_rst_read_cfc( knum ) 
+      !!----------------------------------------------------------------------
+      !!                     ***  trc_rst_read_cfc  ***  
+      !!
+      !! ** Purpose : Read in restart file specific variables from CFC model
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in)  :: knum  ! unit of the restart file
+      INTEGER              :: jn    ! dummy loop indices
+      !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' trc_rst_read_cfc : Read specific variables from CFC model '
+      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
+      
+      DO jn = jp_cfc0, jp_cfc1
+         CALL iom_get( knum, jpdom_autoglo, 'qint_'//ctrcnm(jn), qint_cfc(:,:,jn) ) 
+      END DO
+
+   END SUBROUTINE trc_rst_read_cfc
+
+   SUBROUTINE trc_rst_wri_cfc( kt, kitrst, knum )
+      !!----------------------------------------------------------------------
+      !!                     ***  trc_rst_read_cfc  ***
+      !!
+      !! ** Purpose : Read in restart file specific variables from CFC model
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in)  :: kt      ! time step
+      INTEGER, INTENT(in)  :: kitrst  ! time step of restart write
+      INTEGER, INTENT(in)  :: knum    ! unit of the restart file
+      INTEGER              :: jn      ! dummy loop indices
+      !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' trc_rst_wri_cfc : Write specific variables from CFC model '
+      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
+
+      DO jn = jp_cfc0, jp_cfc1
+         CALL iom_rstput( kt, kitrst, knum, 'qint_'//ctrcnm(jn), qint_cfc(:,:,jn) )
+      END DO
+
+   END SUBROUTINE trc_rst_wri_cfc
+
+#else
+   !!----------------------------------------------------------------------
+   !!  Dummy module :                                     No passive tracer
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_rst_read_cfc( knum )
+     INTEGER, INTENT(in)  :: knum
+     WRITE(*,*) 'trc_rst_wri_cfc: You should not have seen this print! error?', knum
+   END SUBROUTINE trc_rst_read_cfc
+
+   SUBROUTINE trc_rst_wri_cfc( kt, kitrst, knum )
+     INTEGER, INTENT(in)  :: kt, kitrst, knum
+     WRITE(*,*) 'trc_rst_wri_cfc: You should not have seen this print! error?', kt, kitrst, knum
+   END SUBROUTINE trc_rst_wri_cfc
+#endif
+
+   !!======================================================================
+END MODULE trcrst_cfc
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/LOBSTER/trcbio.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/LOBSTER/trcbio.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/LOBSTER/trcbio.F90	(revision 2007)
@@ -482,4 +482,6 @@
       ENDIF
 
+      IF( l_trdtrc ) DEALLOCATE( ztrbio )
+
       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('bio')")
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/LOBSTER/trcexp.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/LOBSTER/trcexp.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/LOBSTER/trcexp.F90	(revision 2007)
@@ -164,4 +164,6 @@
       ENDIF
 
+      IF( l_trdtrc ) DEALLOCATE( ztrbio )
+
       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('exp')")
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/LOBSTER/trcini_lobster.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/LOBSTER/trcini_lobster.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/LOBSTER/trcini_lobster.F90	(revision 2007)
@@ -26,5 +26,4 @@
    PUBLIC   trc_ini_lobster   ! called by trcini.F90 module
 
-#  include "domzgr_substitute.h90"
 #  include "top_substitute.h90"
    !!----------------------------------------------------------------------
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/LOBSTER/trcopt.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/LOBSTER/trcopt.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/LOBSTER/trcopt.F90	(revision 2007)
@@ -28,5 +28,5 @@
 
    !!* Substitution
-#  include "domzgr_substitute.h90"
+#  include "top_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/LOBSTER/trcrst_lobster.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/LOBSTER/trcrst_lobster.F90	(revision 2007)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/LOBSTER/trcrst_lobster.F90	(revision 2007)
@@ -0,0 +1,88 @@
+MODULE trcrst_lobster
+   !!======================================================================
+   !!                       ***  MODULE trcrst_lobster  ***
+   !! TOP :   create, write, read the restart files of LOBSTER tracer
+   !!======================================================================
+   !! History :   1.0  !  2010-01 (C. Ethe) Original
+   !!----------------------------------------------------------------------
+#if defined key_lobster
+   !!----------------------------------------------------------------------
+   !!   'key_lobster'                                               lobster tracers
+   !!----------------------------------------------------------------------
+   !!   trc_rst_read_lobster   : read  restart file
+   !!   trc_rst_wri_lobster    : write restart file
+   !!----------------------------------------------------------------------
+   USE oce_trc         ! Ocean variables
+   USE par_trc         ! TOP parameters
+   USE trc             ! TOP variables
+   USE trcsms_lobster          ! lobster sms trends
+   USE sms_lobster          ! lobster sms trends
+   USE in_out_manager  ! I/O manager
+   USE iom
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC  trc_rst_read_lobster   ! called by trcini.F90 module
+   PUBLIC  trc_rst_wri_lobster   ! called by trcini.F90 module
+
+CONTAINS
+   
+   SUBROUTINE trc_rst_read_lobster( knum ) 
+      !!----------------------------------------------------------------------
+      !!                     ***  trc_rst_read_lobster  ***  
+      !!
+      !! ** Purpose : Read in restart file specific variables from lobster model
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in)  :: knum  ! unit of the restart file
+      !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' trc_rst_read_lobster : Read specific variables from lobster model '
+      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
+
+      CALL iom_get( knum, jpdom_autoglo, 'SEDB'//ctrcnm(jpdet), sedpocb(:,:) ) 
+      CALL iom_get( knum, jpdom_autoglo, 'SEDN'//ctrcnm(jpdet), sedpocn(:,:) ) 
+
+   END SUBROUTINE trc_rst_read_lobster
+
+   SUBROUTINE trc_rst_wri_lobster( kt, kitrst, knum )
+      !!----------------------------------------------------------------------
+      !!                     ***  trc_rst_read_lobster  ***
+      !!
+      !! ** Purpose : Read in restart file specific variables from lobster model
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in)  :: kt      ! time step
+      INTEGER, INTENT(in)  :: kitrst  ! time step of restart write
+      INTEGER, INTENT(in)  :: knum    ! unit of the restart file
+      !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' trc_rst_wri_lobster : Write specific variables from lobster model '
+      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
+
+      CALL iom_rstput( kt, kitrst, knum, 'SEDB'//ctrcnm(jpdet), sedpocb(:,:) )
+      CALL iom_rstput( kt, kitrst, knum, 'SEDN'//ctrcnm(jpdet), sedpocn(:,:) )
+
+   END SUBROUTINE trc_rst_wri_lobster
+
+#else
+   !!----------------------------------------------------------------------
+   !!  Dummy module :                                     No passive tracer
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_rst_read_lobster( knum )
+     INTEGER, INTENT(in)  :: knum
+     WRITE(*,*) 'trc_rst_wri_lobster: You should not have seen this print! error?',knum
+   END SUBROUTINE trc_rst_read_lobster
+
+   SUBROUTINE trc_rst_wri_lobster( kt, kitrst, knum )
+     INTEGER, INTENT(in)  :: kt, kitrst, knum
+     WRITE(*,*) 'trc_rst_wri_lobster: You should not have seen this print! error?', kt, kitrst, knum
+   END SUBROUTINE trc_rst_wri_lobster
+#endif
+
+   !!======================================================================
+END MODULE trcrst_lobster
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/LOBSTER/trcsed.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/LOBSTER/trcsed.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/LOBSTER/trcsed.F90	(revision 2007)
@@ -29,5 +29,5 @@
 
    !!* Substitution
-#  include "domzgr_substitute.h90"
+#  include "top_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
@@ -136,4 +136,6 @@
       ENDIF
 
+      IF( l_trdtrc ) DEALLOCATE( ztrbio )
+
       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('sed')")
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/MY_TRC/trcrst_my_trc.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/MY_TRC/trcrst_my_trc.F90	(revision 2007)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/MY_TRC/trcrst_my_trc.F90	(revision 2007)
@@ -0,0 +1,51 @@
+MODULE trcrst_my_trc
+   !!======================================================================
+   !!                       ***  MODULE trcrst_my_trc  ***
+   !! TOP :   create, write, read the restart files of MY_TRC tracer
+   !!======================================================================
+   !! History :   1.0  !  2010-01 (C. Ethe)  Original
+   !!----------------------------------------------------------------------
+#if defined key_my_trc
+   !!----------------------------------------------------------------------
+   !!   'key_my_trc'                                               CFC tracers
+   !!----------------------------------------------------------------------
+   !!   trc_rst_read_my_trc   : read  restart file
+   !!   trc_rst_wri_my_trc    : write restart file
+   !!----------------------------------------------------------------------
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC  trc_rst_read_my_trc   ! called by trcini.F90 module
+   PUBLIC  trc_rst_wri_my_trc   ! called by trcini.F90 module
+
+CONTAINS
+   
+   SUBROUTINE trc_rst_read_my_trc( knum ) 
+     INTEGER, INTENT(in)  :: knum
+     WRITE(*,*) 'trc_rst_read_my_trc: No specific variables to read on unit', knum
+   END SUBROUTINE trc_rst_read_my_trc
+
+   SUBROUTINE trc_rst_wri_my_trc( kt, kirst, knum )
+     INTEGER, INTENT(in)  :: kt, kirst, knum
+     WRITE(*,*) 'trc_rst_wri_my_trc: No specific variables to write on unit' ,knum, ' at time ', kt, kirst
+   END SUBROUTINE trc_rst_wri_my_trc
+
+#else
+   !!----------------------------------------------------------------------
+   !!  Dummy module :                                     No passive tracer
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_rst_read_my_trc( knum )
+     INTEGER, INTENT(in)  :: knum
+     WRITE(*,*) 'trc_rst_read_my_trc: You should not have seen this print! error?', knum
+   END SUBROUTINE trc_rst_read_my_trc
+
+   SUBROUTINE trc_rst_wri_my_trc( kt, kirst, knum )
+     INTEGER, INTENT(in)  :: kt, kirst, knum
+     WRITE(*,*) 'trc_rst_wri_my_trc: You should not have seen this print! error?', kt, kirst, knum
+   END SUBROUTINE trc_rst_wri_my_trc
+#endif
+
+   !!======================================================================
+END MODULE trcrst_my_trc
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zbio.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zbio.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zbio.F90	(revision 2007)
@@ -39,5 +39,5 @@
 
    !!* Substitution
-#  include "domzgr_substitute.h90"
+#  include "top_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zche.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zche.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zche.F90	(revision 2007)
@@ -147,5 +147,5 @@
 
    !!* Substitution
-#include "domzgr_substitute.h90"
+#include "top_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zflx.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zflx.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zflx.F90	(revision 2007)
@@ -52,5 +52,5 @@
 
    !!* Substitution
-#  include "domzgr_substitute.h90"
+#  include "top_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
@@ -204,12 +204,19 @@
           CALL mpp_sum( t_oce_co2_flx )   ! sum over the global domain
         ENDIF
+        ! Conversion in GtC/yr ; negative for outgoing from ocean
+        t_oce_co2_flx = (-1.) * t_oce_co2_flx  * 12. / 1.e15
+        !
         WRITE(numout,*) ' Atmospheric pCO2    :'
         WRITE(numout,*) '-------------------- : ',kt,'  ',t_atm_co2_flx
         WRITE(numout,*) '(ppm)'
-        WRITE(numout,*) 'Total Flux of Carbon :'
-        WRITE(numout,*) '-------------------- : ',t_oce_co2_flx * 12. / 1e15
-        WRITE(numout,*) '(GtC/an)'
+        WRITE(numout,*) 'Total Flux of Carbon out of the ocean :'
+        WRITE(numout,*) '-------------------- : ',t_oce_co2_flx
+        WRITE(numout,*) '(GtC/yr)'
         t_atm_co2_flx = 0.
         t_oce_co2_flx = 0.
+# if defined key_iomput
+        CALL iom_put( "tatpco2" , t_atm_co2_flx  )
+        CALL iom_put( "tco2flx" , t_oce_co2_flx  )
+#endif
       ENDIF
 #endif
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zlim.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zlim.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zlim.F90	(revision 2007)
@@ -41,5 +41,5 @@
 
    !!* Substitution
-#  include "domzgr_substitute.h90"
+#  include "top_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zlys.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zlys.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zlys.F90	(revision 2007)
@@ -67,5 +67,5 @@
 #if defined key_trc_dia3d && defined key_iomput
       REAL(wp) ::   zrfact2
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zcaldiss, zw3d
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zcaldiss
 #endif
       CHARACTER (len=25) :: charout
@@ -94,4 +94,7 @@
                   ! SET DUMMY VARIABLE FOR TOTAL BORATE
                   zbot  = borat(ji,jj,jk)
+
+                  ! SET DUMMY VARIABLE FOR TOTAL BORATE
+                  zbot  = borat(ji,jj,jk)
                   zfact = rhop (ji,jj,jk) / 1000. + rtrn
 
@@ -171,12 +174,8 @@
 #  else
       zrfact2 = 1.e3 * rfact2r
-      zw3d(:,:,:) = hi  (:,:,:)                    * tmask(:,:,:)
-      CALL iom_put( "PH", zw3d )
-      zw3d(:,:,:) = zco3(:,:,:)                    * tmask(:,:,:)
-      CALL iom_put( "CO3", zw3d )
-      zw3d(:,:,:) = aksp(:,:,:) / calcon           * tmask(:,:,:)
-      CALL iom_put( "CO3sat", zw3d )
-      zw3d(:,:,:) = zcaldiss(:,:,:) * zrfact2 * tmask(:,:,:)
-      CALL iom_put( "Dcal", zw3d )
+      CALL iom_put( "PH"    , hi      (:,:,:)           * tmask(:,:,:) )
+      CALL iom_put( "CO3"   , zco3    (:,:,:)           * tmask(:,:,:) )
+      CALL iom_put( "CO3sat", aksp    (:,:,:) / calcon  * tmask(:,:,:) )
+      CALL iom_put( "DCAL"  , zcaldiss(:,:,:) * zrfact2 * tmask(:,:,:) )
 #  endif
 # endif
@@ -232,5 +231,4 @@
    END SUBROUTINE p4z_lys
 #endif 
-
    !!======================================================================
 END MODULE  p4zlys
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zmeso.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zmeso.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zmeso.F90	(revision 2007)
@@ -45,5 +45,5 @@
 
    !!* Substitution
-#  include "domzgr_substitute.h90"
+#  include "top_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
@@ -76,5 +76,4 @@
 #if defined key_trc_diaadd && defined key_trc_dia3d && defined key_iomput
       REAL(wp) :: zrfact2
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d
 #endif
 
@@ -203,4 +202,10 @@
       END DO
       
+#if defined key_trc_dia3d
+      ! Total grazing ( grazing by microzoo is already computed in p4zmicro ) 
+      grazing(:,:,:) = grazing(:,:,:) + (  zgrazd  (:,:,:) + zgrazz  (:,:,:) + zgrazn(:,:,:) &
+                     &                   + zgrazpoc(:,:,:) + zgrazffe(:,:,:)  )
+#endif
+
 
       DO jk = 1,jpkm1
@@ -311,10 +316,12 @@
 #if defined key_trc_diaadd && defined key_trc_dia3d && defined key_iomput
       zrfact2 = 1.e3 * rfact2r
-      zw3d(:,:,:) = (     zgrazd(:,:,:) +   zgrazz(:,:,:) + zgrazn(:,:,:) &
-                    & + zgrazpoc(:,:,:) + zgrazffe(:,:,:)                 ) * zrfact2 * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "Graz2" , zw3d )
-
-      zw3d(:,:,:) = prodcal(:,:,:) * zrfact2 * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "Pcal"  , zw3d )
+      ! Total grazing of phyto by zoo
+      grazing(:,:,:) = grazing(:,:,:) * zrfact2 * tmask(:,:,:)
+      ! Calcite production
+      prodcal(:,:,:) = prodcal(:,:,:) * zrfact2 * tmask(:,:,:)
+      IF( jnt == nrdttrc ) then 
+         CALL iom_put( "GRAZ" , grazing  )  ! Total grazing of phyto by zooplankton
+         CALL iom_put( "PCAL" , prodcal  )  ! Calcite production
+      ENDIF
 #endif
 
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zmicro.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zmicro.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zmicro.F90	(revision 2007)
@@ -43,5 +43,5 @@
 
    !!* Substitution
-#  include "domzgr_substitute.h90"
+#  include "top_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
@@ -70,8 +70,4 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgrazmf, zgrazsf, zgrazpf
       CHARACTER (len=25) :: charout
-#if defined key_trc_diaadd && defined key_trc_dia3d && defined key_iomput
-      REAL(wp) :: zrfact2
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d
-#endif
 
       !!---------------------------------------------------------------------
@@ -88,4 +84,7 @@
       zgrazpf(:,:,:) = 0.
 
+#if defined key_trc_dia3d
+      grazing(:,:,:) = 0.  !: Initialisation of  grazing
+#endif
 
       zstep = rfact2 / rday      ! Time step duration for biology
@@ -156,4 +155,8 @@
       END DO
       
+#if defined key_trc_dia3d
+      ! Grazing by microzooplankton
+      grazing(:,:,:) = grazing(:,:,:) + zgrazp(:,:,:) + zgrazm(:,:,:) + zgrazsd(:,:,:) 
+#endif
 
       DO jk = 1,jpkm1
@@ -231,15 +234,9 @@
       END DO
       !
-#if defined key_trc_diaadd && defined key_trc_dia3d && defined key_iomput
-      zrfact2 = 1.e3 * rfact2r
-      zw3d(:,:,:) = ( zgrazp(:,:,:) + zgrazm(:,:,:) + zgrazsd(:,:,:) ) * zrfact2 * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "Graz" , zw3d )
-#endif
-
-       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('micro')")
          CALL prt_ctl_trc_info(charout)
          CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
-       ENDIF
+      ENDIF
 
    END SUBROUTINE p4z_micro
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zmort.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zmort.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zmort.F90	(revision 2007)
@@ -41,5 +41,5 @@
 
    !!* Substitution
-#  include "domzgr_substitute.h90"
+#  include "top_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zopt.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zopt.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zopt.F90	(revision 2007)
@@ -35,5 +35,5 @@
    
    !!* Substitution
-#  include "domzgr_substitute.h90"
+#  include "top_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
@@ -61,8 +61,4 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zekg, zekr, zekb
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze1 , ze2 , ze3, ze0
-#if defined key_trc_diaadd && defined key_iomput
-     REAL(wp), DIMENSION(jpi,jpj)      ::   zw2d
-     REAL(wp), DIMENSION(jpi,jpj,jpk)  ::   zw3d
-#endif
       !!---------------------------------------------------------------------
 
@@ -238,8 +234,8 @@
 # else
       ! write diagnostics 
-      zw2d(:,:  ) =  heup(:,:  ) * tmask(:,:,1)
-      zw3d(:,:,:) =  etot(:,:,:) * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "Heup", zw2d )               
-      IF( jnt == nrdttrc ) CALL iom_put( "PAR" , zw3d )
+      IF( jnt == nrdttrc ) then 
+         CALL iom_put( "Heup", heup(:,:  ) * tmask(:,:,1) )  ! euphotic layer deptht
+         CALL iom_put( "PAR" , etot(:,:,:) * tmask(:,:,:) )  ! Photosynthetically Available Radiation
+      ENDIF
 # endif
 #endif
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zprod.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zprod.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zprod.F90	(revision 2007)
@@ -53,5 +53,5 @@
 
    !!* Substitution
-#  include "domzgr_substitute.h90"
+#  include "top_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
@@ -81,7 +81,4 @@
 #if defined key_trc_diaadd && defined key_trc_dia3d
       REAL(wp) ::   zrfact2
-#if  defined key_iomput
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d
-#endif
 #endif
       REAL(wp), DIMENSION(jpi,jpj)     ::   zmixnano   , zmixdiat, zstrn
@@ -352,12 +349,11 @@
         WRITE(numout,*) 'Total PP :'
         WRITE(numout,*) '-------------------- : ',tpp * 12. / 1.E12
-        WRITE(numout,*) '(GtC/an)'
+        WRITE(numout,*) '(GtC/yr)'
         tpp = 0.
       ENDIF
 
-#if defined key_trc_diaadd && defined key_trc_dia3d
+#if defined key_trc_diaadd && defined key_trc_dia3d && ! defined key_iomput
+      !   Supplementary diagnostics
       zrfact2 = 1.e3 * rfact2r
-      !   Supplementary diagnostics
-#  if ! defined key_iomput
       trc3d(:,:,:,jp_pcs0_3d + 4)  = zprorca (:,:,:) * zrfact2 * tmask(:,:,:)
       trc3d(:,:,:,jp_pcs0_3d + 5)  = zprorcad(:,:,:) * zrfact2 * tmask(:,:,:)
@@ -366,24 +362,20 @@
       trc3d(:,:,:,jp_pcs0_3d + 8)  = zprorcad(:,:,:) * zrfact2 * tmask(:,:,:) * zysopt(:,:,:)
       trc3d(:,:,:,jp_pcs0_3d + 9)  = zprofed (:,:,:) * zrfact2 * tmask(:,:,:)
-#if ! defined key_kriest
+#  if ! defined key_kriest
       trc3d(:,:,:,jp_pcs0_3d + 10) = zprofen (:,:,:) * zrfact2 * tmask(:,:,:)
+#  endif
 #endif
 
-# else
-      zw3d(:,:,:) = zprorca (:,:,:) * zrfact2 * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "PPPHY" , zw3d )
-      zw3d(:,:,:) = zprorcad(:,:,:) * zrfact2 * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "PPPHY2", zw3d )
-      zw3d(:,:,:) = zpronew (:,:,:) * zrfact2 * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "PPNEWN" , zw3d )
-      zw3d(:,:,:) = zpronewd(:,:,:) * zrfact2 * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "PPNEWD", zw3d )
-      zw3d(:,:,:) = zprorcad(:,:,:) * zrfact2 * tmask(:,:,:) * zysopt(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "PBSi"  , zw3d )
-      zw3d(:,:,:) = zprofed (:,:,:) * zrfact2 * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "PFeD"  , zw3d )
-      zw3d(:,:,:) = zprofen (:,:,:) * zrfact2 * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "PFeN"  , zw3d )
-# endif
+#if defined key_trc_diaadd && defined key_trc_dia3d && defined key_iomput
+      zrfact2 = 1.e3 * rfact2r
+      IF ( jnt == nrdttrc ) then
+         CALL iom_put( "PPPHY" , zprorca (:,:,:) * zrfact2 * tmask(:,:,:) )  ! primary production by nanophyto
+         CALL iom_put( "PPPHY2", zprorcad(:,:,:) * zrfact2 * tmask(:,:,:) )  ! primary production by diatom
+         CALL iom_put( "PPNEWN", zpronew (:,:,:) * zrfact2 * tmask(:,:,:) )  ! new primary production by nanophyto
+         CALL iom_put( "PPNEWD", zpronewd(:,:,:) * zrfact2 * tmask(:,:,:) )  ! new primary production by diatom
+         CALL iom_put( "PBSi"  , zprorcad(:,:,:) * zrfact2 * tmask(:,:,:) * zysopt(:,:,:) ) ! biogenic silica production
+         CALL iom_put( "PFeD"  , zprofed (:,:,:) * zrfact2 * tmask(:,:,:) )  ! biogenic iron production by diatom
+         CALL iom_put( "PFeN"  , zprofen (:,:,:) * zrfact2 * tmask(:,:,:) )  ! biogenic iron production by nanophyto
+      ENDIF
 #endif
 
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zrem.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zrem.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zrem.F90	(revision 2007)
@@ -45,5 +45,5 @@
 
    !!* Substitution
-#  include "domzgr_substitute.h90"
+#  include "top_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zsed.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zsed.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zsed.F90	(revision 2007)
@@ -96,5 +96,4 @@
       REAL(wp) :: zrfact2
 # if defined key_iomput
-     REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d 
      REAL(wp), DIMENSION(jpi,jpj)    ::    zw2d 
 # endif
@@ -332,12 +331,11 @@
       trc2d(:,:,jp_pcs0_2d + 12) = znitrpot(:,:,1) * 1.e-7 * zrfact2 * fse3t(:,:,1) * tmask(:,:,1)
 # else
-      ! write diagnostics
-      zw2d(:,:)   =  ( zirondep(:,:,1) + ironsed(:,:,1) * rfact2 ) &
-      &            * zrfact2 * fse3t(:,:,1) * tmask(:,:,1)     
+      ! surface downward net flux of iron
+      zw2d(:,:)   =  ( zirondep(:,:,1) + ironsed(:,:,1) * rfact2 ) * zrfact2 * fse3t(:,:,1) * tmask(:,:,1) 
       IF( jnt == nrdttrc ) CALL iom_put( "Irondep", zw2d )
-      zw3d(:,:,:) = znitrpot(:,:,:) * 1.e-7 * zrfact2  * fse3t(:,:,:) * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "Nfix", zw3d  ) 
-# endif
-
+      ! nitrogen fixation at surface
+      zw2d(:,:)   =  znitrpot(:,:,1) * 1.e-7 * zrfact2  * fse3t(:,:,1) * tmask(:,:,1)
+      IF( jnt == nrdttrc ) CALL iom_put( "Nfix" , zw2d )
+# endif
 # endif
       !
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zsink.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zsink.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zsink.F90	(revision 2007)
@@ -69,5 +69,5 @@
 
    !!* Substitution
-#  include "domzgr_substitute.h90"
+#  include "top_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
@@ -99,8 +99,5 @@
 #if defined key_trc_diaadd
       REAL(wp) :: zrfact2
-      INTEGER  :: iksed1
-#if defined key_iomput
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d
-#endif
+      INTEGER  :: ik1
 #endif
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   znum3d
@@ -286,11 +283,11 @@
 #if defined key_trc_diaadd
       zrfact2 = 1.e3 * rfact2r
-      iksed1 = iksed + 1
+      ik1 = iksed + 1
 #  if ! defined key_iomput
-      trc2d(:,:  ,jp_pcs0_2d + 4)  = sinking (:,:,iksed1) * zrfact2 * tmask(:,:,1)
-      trc2d(:,:  ,jp_pcs0_2d + 5)  = sinking2(:,:,iksed1) * zrfact2 * tmask(:,:,1)
-      trc2d(:,:  ,jp_pcs0_2d + 6)  = sinkfer (:,:,iksed1) * zrfact2 * tmask(:,:,1)
-      trc2d(:,:  ,jp_pcs0_2d + 7)  = sinksil (:,:,iksed1) * zrfact2 * tmask(:,:,1)
-      trc2d(:,:  ,jp_pcs0_2d + 8)  = sinkcal (:,:,iksed1) * zrfact2 * tmask(:,:,1)
+      trc2d(:,:  ,jp_pcs0_2d + 4)  = sinking (:,:,ik1) * zrfact2 * tmask(:,:,1)
+      trc2d(:,:  ,jp_pcs0_2d + 5)  = sinking2(:,:,ik1) * zrfact2 * tmask(:,:,1)
+      trc2d(:,:  ,jp_pcs0_2d + 6)  = sinkfer (:,:,ik1) * zrfact2 * tmask(:,:,1)
+      trc2d(:,:  ,jp_pcs0_2d + 7)  = sinksil (:,:,ik1) * zrfact2 * tmask(:,:,1)
+      trc2d(:,:  ,jp_pcs0_2d + 8)  = sinkcal (:,:,ik1) * zrfact2 * tmask(:,:,1)
       trc3d(:,:,:,jp_pcs0_3d + 11) = sinking (:,:,:)      * zrfact2 * tmask(:,:,:)
       trc3d(:,:,:,jp_pcs0_3d + 12) = sinking2(:,:,:)      * zrfact2 * tmask(:,:,:)
@@ -301,28 +298,18 @@
       trc3d(:,:,:,jp_pcs0_3d + 17) = wsbio4  (:,:,:)                * tmask(:,:,:)
 #else
-      zw3d(:,:,:)  = sinking (:,:,:) * zrfact2 * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "PMO" , zw3d )
-      zw3d(:,:,:)  = sinking2(:,:,:) * zrfact2 * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "PMO2", zw3d )
-      zw3d(:,:,:)  = sinkfer (:,:,:) * zrfact2 * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "ExpFe1", zw3d )
-      zw3d(:,:,:)  = sinksil (:,:,:) * zrfact2 * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "ExpSi", zw3d )
-      zw3d(:,:,:)  = sinkcal (:,:,:) * zrfact2 * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "ExpCaCO3", zw3d )
-      zw3d(:,:,:)  = sinking (:,:,:) * zrfact2 * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "POCFlx", zw3d )
-      zw3d(:,:,:)  = sinking2(:,:,:) * zrfact2 * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "GOCFlx", zw3d )
-      zw3d(:,:,:)  = sinksil (:,:,:) * zrfact2 * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "SiFlx", zw3d )
-      zw3d(:,:,:)  = sinkcal (:,:,:) * zrfact2 * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "CaCO3Flx", zw3d )
-      zw3d(:,:,:)  = znum3d  (:,:,:)           * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "xnum", zw3d )
-      zw3d(:,:,:)  = wsbio3  (:,:,:)           * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "W1", zw3d )
-      zw3d(:,:,:)  = wsbio4  (:,:,:)           * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "W2", zw3d )
+      IF( jnt == nrdttrc ) then
+        CALL iom_put( "POCFlx"  , sinking (:,:,:)      * zrfact2 * tmask(:,:,:) )  ! POC export
+        CALL iom_put( "NumFlx"  , sinking2 (:,:,:)     * zrfact2 * tmask(:,:,:) )  ! Num export
+        CALL iom_put( "SiFlx"   , sinksil (:,:,:)      * zrfact2 * tmask(:,:,:) )  ! Silica export
+        CALL iom_put( "CaCO3Flx", sinkcal (:,:,:)      * zrfact2 * tmask(:,:,:) )  ! Calcite export
+        CALL iom_put( "xnum"    , znum3d  (:,:,:)                * tmask(:,:,:) )  ! Number of particles in aggregats
+        CALL iom_put( "W1"      , wsbio3  (:,:,:)                * tmask(:,:,:) )  ! sinking speed of POC
+        CALL iom_put( "W2"      , wsbio4  (:,:,:)                * tmask(:,:,:) )  ! sinking speed of aggregats
+        CALL iom_put( "PMO"     , sinking (:,:,ik1) * zrfact2 * tmask(:,:,1) )  ! POC export at 100m
+        CALL iom_put( "PMO2"    , sinking2(:,:,ik1) * zrfact2 * tmask(:,:,1) )  ! Num export at 100m
+        CALL iom_put( "ExpFe1"  , sinkfer (:,:,ik1) * zrfact2 * tmask(:,:,1) )  ! Export of iron at 100m
+        CALL iom_put( "ExpSi"   , sinksil (:,:,ik1) * zrfact2 * tmask(:,:,1) )  ! export of silica at 100m
+        CALL iom_put( "ExpCaCO3", sinkcal (:,:,ik1) * zrfact2 * tmask(:,:,1) )  ! export of calcite at 100m
+     ENDIF
 #  endif
 
@@ -489,8 +476,5 @@
 #if defined key_trc_dia3d
       REAL(wp) ::   zrfact2
-      INTEGER  ::   iksed1
-#endif
-#if defined key_iomput
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d
+      INTEGER  ::   ik1
 #endif
       CHARACTER (len=25) :: charout
@@ -613,26 +597,20 @@
 #if defined key_trc_diaadd
       zrfact2 = 1.e3 * rfact2r
-      iksed1 = iksed + 1
+      ik1  = iksed + 1
 #  if ! defined key_iomput
-      trc2d(:,:,jp_pcs0_2d + 4) = sinking (:,:,iksed1) * zrfact2 * tmask(:,:,1)
-      trc2d(:,:,jp_pcs0_2d + 5) = sinking2(:,:,iksed1) * zrfact2 * tmask(:,:,1)
-      trc2d(:,:,jp_pcs0_2d + 6) = sinkfer (:,:,iksed1) * zrfact2 * tmask(:,:,1)
-      trc2d(:,:,jp_pcs0_2d + 7) = sinkfer2(:,:,iksed1) * zrfact2 * tmask(:,:,1)
-      trc2d(:,:,jp_pcs0_2d + 8) = sinksil (:,:,iksed1) * zrfact2 * tmask(:,:,1)
-      trc2d(:,:,jp_pcs0_2d + 9) = sinkcal (:,:,iksed1) * zrfact2 * tmask(:,:,1)
+      trc2d(:,:,jp_pcs0_2d + 4) = sinking (:,:,ik1) * zrfact2 * tmask(:,:,1)
+      trc2d(:,:,jp_pcs0_2d + 5) = sinking2(:,:,ik1) * zrfact2 * tmask(:,:,1)
+      trc2d(:,:,jp_pcs0_2d + 6) = sinkfer (:,:,ik1) * zrfact2 * tmask(:,:,1)
+      trc2d(:,:,jp_pcs0_2d + 7) = sinkfer2(:,:,ik1) * zrfact2 * tmask(:,:,1)
+      trc2d(:,:,jp_pcs0_2d + 8) = sinksil (:,:,ik1) * zrfact2 * tmask(:,:,1)
+      trc2d(:,:,jp_pcs0_2d + 9) = sinkcal (:,:,ik1) * zrfact2 * tmask(:,:,1)
 #  else
-      zw3d(:,:,:)  = sinking (:,:,:) * zrfact2 * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "ExpPOC" , zw3d )
-      zw3d(:,:,:)  = sinking2(:,:,:) * zrfact2 * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "ExpGOC", zw3d )
-      zw3d(:,:,:)  = sinkfer (:,:,:) * zrfact2 * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "ExpFe1", zw3d )
-      zw3d(:,:,:)  = sinkfer2(:,:,:) * zrfact2 * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "ExpFe2", zw3d )
-      zw3d(:,:,:)  = sinksil (:,:,:) * zrfact2 * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "ExpSi", zw3d )
-      zw3d(:,:,:)  = sinkcal (:,:,:) * zrfact2 * tmask(:,:,:)
-      IF( jnt == nrdttrc ) CALL iom_put( "Expcal", zw3d )
-#  endif
+      IF( jnt == nrdttrc )  then
+         CALL iom_put( "EPC100"  , ( sinking(:,:,ik1) + sinking2(:,:,ik1) ) * zrfact2 * tmask(:,:,1) ) ! Export of carbon at 100m
+         CALL iom_put( "EPFE100" , ( sinkfer(:,:,ik1) + sinkfer2(:,:,ik1) ) * zrfact2 * tmask(:,:,1) ) ! Export of iron at 100m
+         CALL iom_put( "EPCAL100",   sinkcal(:,:,ik1)                       * zrfact2 * tmask(:,:,1) ) ! Export of calcite  at 100m
+         CALL iom_put( "EPSI100" ,   sinksil(:,:,ik1)                       * zrfact2 * tmask(:,:,1) ) ! Export of biogenic silica at 100m
+      ENDIF
+#endif
 #endif
       !
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90	(revision 2007)
@@ -38,4 +38,6 @@
    LOGICAL  ::   ln_pisdmp         !: relaxation or not of nutrients to a mean value
                                    !: when initialize from a restart file 
+   LOGICAL  ::   ln_pisclo         !: Restoring or not of nutrients to initial value
+                                   !: on close seas
 
    !!*  Biological fluxes for light
@@ -62,4 +64,5 @@
 #if defined key_trc_dia3d
    REAL(wp), DIMENSION(jpi,jpj,jpk) ::   prodcal    !: Calcite production
+   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   grazing    !: Total zooplankton grazing
 #endif
 
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90	(revision 2007)
@@ -38,5 +38,4 @@
       no3    =  31.04e-6 * 7.6
 
-#  include "domzgr_substitute.h90"
 #  include "top_substitute.h90"
    !!----------------------------------------------------------------------
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/trclsm_pisces.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/trclsm_pisces.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/trclsm_pisces.F90	(revision 2007)
@@ -67,5 +67,5 @@
       NAMELIST/nampisdia/ nwritedia, pisdia3d, pisdia2d     ! additional diagnostics
 #endif
-      NAMELIST/nampisdmp/ ln_pisdmp
+      NAMELIST/nampisdmp/ ln_pisdmp, ln_pisclo
 
       !!----------------------------------------------------------------------
@@ -188,5 +188,6 @@
          WRITE(numout,*)
          WRITE(numout,*) ' Namelist : nampisdmp'
-         WRITE(numout,*) '    Relaxation of tracer to glodap mean value    ln_pisdmp      =', ln_pisdmp
+         WRITE(numout,*) '    Relaxation of tracer to glodap mean value            ln_pisdmp      =', ln_pisdmp
+         WRITE(numout,*) '    Restoring of tracer to initial value  on closed seas  ln_pisclo      =', ln_pisclo
          WRITE(numout,*) ' '
       ENDIF
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/trcrst_pisces.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/trcrst_pisces.F90	(revision 2007)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/PISCES/trcrst_pisces.F90	(revision 2007)
@@ -0,0 +1,301 @@
+MODULE trcrst_pisces
+   !!======================================================================
+   !!                       ***  MODULE trcrst_pisces  ***
+   !! TOP :   create, write, read the restart files of PISCES tracer
+   !!======================================================================
+   !! History :   1.0  !  2010-01 (C. Ethe) Original
+   !!----------------------------------------------------------------------
+#if defined key_pisces
+   !!----------------------------------------------------------------------
+   !!   'key_pisces'                                               pisces tracers
+   !!----------------------------------------------------------------------
+   !!   trc_rst_read_pisces   : read  restart file
+   !!   trc_rst_wri_pisces    : write restart file
+   !!----------------------------------------------------------------------
+   USE oce_trc         ! Ocean variables
+   USE par_trc         ! TOP parameters
+   USE trc             ! TOP variables
+   USE trcsms_pisces          ! pisces sms trends
+   USE sms_pisces          ! pisces sms variables
+   USE in_out_manager  ! I/O manager
+   USE iom
+   USE trcdta
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC  trc_rst_read_pisces   ! called by trcini.F90 module
+   PUBLIC  trc_rst_wri_pisces   ! called by trcini.F90 module
+
+CONTAINS
+   
+   SUBROUTINE trc_rst_read_pisces( knum ) 
+      !!----------------------------------------------------------------------
+      !!                     ***  trc_rst_read_pisces  ***  
+      !!
+      !! ** Purpose : Read in restart file specific variables from pisces model
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in)  :: knum  ! unit of the restart file
+      INTEGER  ::  ji, jj, jk
+      REAL(wp) ::  zcaralk, zbicarb, zco3
+      REAL(wp) ::  ztmas, ztmas1
+      !!----------------------------------------------------------------------
+
+      !
+      IF( lk_dtatrc .AND. ln_pisclo ) CALL pis_dmp_clo  ! restoring of nutrients on close seas
+      IF( ln_pisdmp )                 CALL pis_dmp_ini  ! relaxation of some tracers
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' trc_rst_read_pisces : Read specific variables from pisces model '
+      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
+      ! 
+      IF( iom_varid( knum, 'PH', ldstop = .FALSE. ) > 0 ) THEN
+         CALL iom_get( knum, jpdom_autoglo, 'PH' , hi(:,:,:)  )
+      ELSE
+         ! Set PH from  total alkalinity, borat (???), akb3 (???) and ak23 (???)
+         ! --------------------------------------------------------
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztmas   = tmask(ji,jj,jk)
+                  ztmas1  = 1. - tmask(ji,jj,jk)
+                  zcaralk = trn(ji,jj,jk,jptal) - borat(ji,jj,jk) / (  1. + 1.E-8 / ( rtrn + akb3(ji,jj,jk) )  )
+                  zco3    = ( zcaralk - trn(ji,jj,jk,jpdic) ) * ztmas + 0.5e-3 * ztmas1
+                  zbicarb = ( 2. * trn(ji,jj,jk,jpdic) - zcaralk )
+                  hi(ji,jj,jk) = ( ak23(ji,jj,jk) * zbicarb / zco3 ) * ztmas + 1.e-9 * ztmas1
+               END DO
+            END DO
+         END DO
+      ENDIF
+      CALL iom_get( knum, jpdom_autoglo, 'Silicalim', xksi(:,:) ) 
+      IF( iom_varid( knum, 'Silicamax', ldstop = .FALSE. ) > 0 ) THEN
+         CALL iom_get( knum, jpdom_autoglo, 'Silicamax' , xksimax(:,:)  )
+      ELSE
+         xksimax(:,:) = xksi(:,:)
+      ENDIF
+
+   END SUBROUTINE trc_rst_read_pisces
+
+   SUBROUTINE trc_rst_wri_pisces( kt, kitrst, knum )
+      !!----------------------------------------------------------------------
+      !!                     ***  trc_rst_read_pisces  ***
+      !!
+      !! ** Purpose : Read in restart file specific variables from pisces model
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in)  :: kt      ! time step
+      INTEGER, INTENT(in)  :: kitrst  ! time step of restart write
+      INTEGER, INTENT(in)  :: knum    ! unit of the restart file
+      !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' trc_rst_wri_pisces : Write specific variables from pisces model '
+      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
+
+      CALL iom_rstput( kt, kitrst, knum, 'PH', hi(:,:,:) )
+      CALL iom_rstput( kt, kitrst, knum, 'Silicalim', xksi(:,:) ) 
+      CALL iom_rstput( kt, kitrst, knum, 'Silicamax', xksimax(:,:) )
+
+   END SUBROUTINE trc_rst_wri_pisces
+
+   SUBROUTINE pis_dmp_ini
+      !!----------------------------------------------------------------------
+      !!                    ***  pis_dmp_ini  ***
+      !!
+      !! ** purpose  : Relaxation of some tracers
+      !!----------------------------------------------------------------------
+      INTEGER  :: ji, jj, jk
+      REAL(wp) ::  &
+         alkmean = 2426. ,  & ! mean value of alkalinity ( Glodap ; for Goyet 2391. )
+         po4mean = 2.165 ,  & ! mean value of phosphates
+         no3mean = 30.90 ,  & ! mean value of nitrate
+         silmean = 91.51      ! mean value of silicate
+
+      REAL(wp) :: zarea, zvol, zalksum, zpo4sum, zno3sum, zsilsum
+
+
+      IF(lwp)  WRITE(numout,*)
+
+      IF( cp_cfg == "orca" .AND. .NOT. lk_trc_c1d ) THEN      ! ORCA condiguration (not 1D) !
+         !                                                    ! --------------------------- !
+         ! set total alkalinity, phosphate, nitrate & silicate
+
+         zalksum = 0.e0
+         zpo4sum = 0.e0
+         zno3sum = 0.e0
+         zsilsum = 0.e0
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zvol = cvol(ji,jj,jk)
+#  if defined key_off_degrad
+                  zvol = zvol * facvol(ji,jj,jk)
+#  endif
+                  zalksum = zalksum + trn(ji,jj,jk,jptal) * zvol
+                  zpo4sum = zpo4sum + trn(ji,jj,jk,jppo4) * zvol
+                  zno3sum = zno3sum + trn(ji,jj,jk,jpno3) * zvol
+                  zsilsum = zsilsum + trn(ji,jj,jk,jpsil) * zvol
+               END DO
+            END DO
+         END DO
+         IF( lk_mpp )   CALL mpp_sum( zalksum )     ! sum over the global domain
+         IF( lk_mpp )   CALL mpp_sum( zpo4sum )     ! sum over the global domain
+         IF( lk_mpp )   CALL mpp_sum( zno3sum )     ! sum over the global domain
+         IF( lk_mpp )   CALL mpp_sum( zsilsum )     ! sum over the global domain
+         zarea   = 1. / areatot * 1.e6
+         zalksum = zalksum * zarea
+         zpo4sum = zpo4sum * zarea / 122.
+         zno3sum = zno3sum * zarea / 7.6
+         zsilsum = zsilsum * zarea
+
+         IF(lwp) WRITE(numout,*) '       TALK mean : ', zalksum
+         trn(:,:,:,jptal) = trn(:,:,:,jptal) * alkmean / zalksum
+            
+         IF(lwp) WRITE(numout,*) '       PO4  mean : ', zpo4sum
+         trn(:,:,:,jppo4) = trn(:,:,:,jppo4) * po4mean / zpo4sum
+
+         IF(lwp) WRITE(numout,*) '       NO3  mean : ', zno3sum
+         trn(:,:,:,jpno3) = trn(:,:,:,jpno3) * no3mean / zno3sum
+
+         IF(lwp) WRITE(numout,*) '       SiO3 mean : ', zsilsum
+         trn(:,:,:,jpsil) = MIN( 400.e-6,trn(:,:,:,jpsil) * silmean / zsilsum )
+         !
+      ENDIF
+
+!#if defined key_kriest
+!     !! Initialize number of particles from a standart restart file
+!     !! The name of big organic particles jpgoc has been only change
+!     !! and replace by jpnum but the values here are concentration
+!     trn(:,:,:,jppoc) = trn(:,:,:,jppoc) + trn(:,:,:,jpnum) 
+!     trn(:,:,:,jpnum) = trn(:,:,:,jppoc) / ( 6. * xkr_massp )
+!#endif
+
+   END SUBROUTINE pis_dmp_ini
+
+   SUBROUTINE pis_dmp_clo   
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE pis_dmp_clo  ***
+      !!
+      !! ** Purpose :   Closed sea domain initialization
+      !!
+      !! ** Method  :   if a closed sea is located only in a model grid point
+      !!                we restore to initial data
+      !!
+      !! ** Action  :   ictsi1(), ictsj1() : south-west closed sea limits (i,j)
+      !!                ictsi2(), ictsj2() : north-east Closed sea limits (i,j)
+      !!----------------------------------------------------------------------
+      INTEGER, PARAMETER           ::   npicts   = 4       !: number of closed sea
+      INTEGER, DIMENSION(npicts)   ::   ictsi1, ictsj1     !: south-west closed sea limits (i,j)
+      INTEGER, DIMENSION(npicts)   ::   ictsi2, ictsj2     !: north-east closed sea limits (i,j)
+      INTEGER :: ji, jj, jk, jn, jc            ! dummy loop indices
+      !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*) 
+      IF(lwp) WRITE(numout,*)' pis_dmp_clo : closed seas '
+      IF(lwp) WRITE(numout,*)'~~~~~~~'
+
+      ! initial values
+      ictsi1(:) = 1  ;  ictsi2(:) = 1
+      ictsj1(:) = 1  ;  ictsj2(:) = 1
+
+      ! set the closed seas (in data domain indices)
+      ! -------------------
+
+      IF( cp_cfg == "orca" ) THEN
+         !
+         SELECT CASE ( jp_cfg )
+         !                                           ! =======================
+         CASE ( 2 )                                  !  ORCA_R2 configuration
+            !                                        ! =======================
+            !                                            ! Caspian Sea
+            ictsi1(1)   =  11  ;  ictsj1(1)   = 103
+            ictsi2(1)   =  17  ;  ictsj2(1)   = 112
+            !                                            ! Great North American Lakes
+            ictsi1(2)   =  97  ;  ictsj1(2)   = 107
+            ictsi2(2)   = 103  ;  ictsj2(2)   = 111
+            !                                            ! Black Sea 1 : west part of the Black Sea
+            ictsi1(3)   = 174  ; ictsj1(3)   = 107
+            ictsi2(3)   = 181  ; ictsj2(3)   = 112
+            !                                            ! Black Sea 2 : est part of the Black Sea
+            ictsi1(4)   =   2  ;  ictsj1(4)   = 107
+            ictsi2(4)   =   6  ;  ictsj2(4)   = 112
+            !                                        ! =======================
+         CASE ( 4 )                                  !  ORCA_R4 configuration
+            !                                        ! =======================
+            !                                            ! Caspian Sea
+            ictsi1(1)   =  4  ;  ictsj1(1)   = 53
+            ictsi2(1)   =  4  ;  ictsj2(1)   = 56
+            !                                            ! Great North American Lakes
+            ictsi1(2)   = 49  ;  ictsj1(2)   = 55
+            ictsi2(2)   = 51  ;  ictsj2(2)   = 56
+            !                                            ! Black Sea
+            ictsi1(3)   = 88  ;  ictsj1(3)   = 55
+            ictsi2(3)   = 91  ;  ictsj2(3)   = 56
+            !                                            ! Baltic Sea
+            ictsi1(4)   = 75  ;  ictsj1(4)   = 59
+            ictsi2(4)   = 76  ;  ictsj2(4)   = 61
+            !                                        ! =======================
+            !                                        ! =======================
+         CASE ( 025 )                                ! ORCA_R025 configuration
+            !                                        ! =======================
+                                                     ! Caspian + Aral sea
+            ictsi1(1)   = 1330 ; ictsj1(1)   = 645
+            ictsi2(1)   = 1400 ; ictsj2(1)   = 795
+            !                                        ! Azov Sea
+            ictsi1(2)   = 1284 ; ictsj1(2)   = 722
+            ictsi2(2)   = 1304 ; ictsj2(2)   = 747
+            !
+         END SELECT
+         !
+      ENDIF
+
+      ! convert the position in local domain indices
+      ! --------------------------------------------
+      DO jc = 1, npicts
+         ictsi1(jc)   = mi0( ictsi1(jc) )
+         ictsj1(jc)   = mj0( ictsj1(jc) )
+
+         ictsi2(jc)   = mi1( ictsi2(jc) )
+         ictsj2(jc)   = mj1( ictsj2(jc) )
+      END DO
+
+#if defined key_dtatrc
+      ! Restore close seas values to initial data
+      CALL trc_dta( nittrc000 ) 
+      DO jn = 1, jptra
+         IF( lutini(jn) ) THEN
+            DO jc = 1, npicts
+               DO jk = 1, jpkm1
+                  DO jj = ictsj1(jc), ictsj2(jc)
+                     DO ji = ictsi1(jc), ictsi2(jc)
+                        trn(ji,jj,jk,jn) = trdta(ji,jj,jk,jn) * tmask(ji,jj,jk) 
+                        trb(ji,jj,jk,jn) = trn(ji,jj,jk,jn)
+                     ENDDO
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDIF
+      ENDDO
+#endif
+   !
+   END SUBROUTINE pis_dmp_clo
+
+#else
+   !!----------------------------------------------------------------------
+   !!  Dummy module :                                     No passive tracer
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_rst_read_pisces( knum )
+      INTEGER, INTENT(in)  :: knum
+      WRITE(*,*) 'trc_rst_read_pisces: You should not have seen this print! error?', knum
+   END SUBROUTINE trc_rst_read_pisces
+
+   SUBROUTINE trc_rst_wri_pisces( kt, kitrst, knum )
+     INTEGER, INTENT(in)  :: kt, kitrst, knum
+     WRITE(*,*) 'trc_rst_wri_pisces: You should not have seen this print! error?', kt, kitrst, knum
+   END SUBROUTINE trc_rst_wri_pisces
+#endif
+
+   !!======================================================================
+END MODULE trcrst_pisces
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90	(revision 2007)
@@ -53,5 +53,5 @@
 
    !! * Substitutions
-#  include "domzgr_substitute.h90"
+#  include "top_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 1.0 , LOCEAN-IPSL (2005) 
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf_imp.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf_imp.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf_imp.F90	(revision 2007)
@@ -112,38 +112,40 @@
          rdttrc(:) =  rdttra(:) * FLOAT(ndttrc)      
       ENDIF
-     !                                                       ! ===========
+
+      ! Initialisation
+      zwd( 1 ,:,:) = 0.e0     ;     zwd(jpi,:,:) = 0.e0
+      zws( 1 ,:,:) = 0.e0     ;     zws(jpi,:,:) = 0.e0
+      zwi( 1 ,:,:) = 0.e0     ;     zwi(jpi,:,:) = 0.e0
+      !                                          
+      ! 0. Matrix construction 
+      ! ----------------------
+
+      ! Diagonal, inferior, superior
+      ! (including the bottom boundary condition via avs masked
+      DO jk = 1, jpkm1                    
+         DO jj = 2, jpjm1                                    
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zwi(ji,jj,jk) = - rdttrc(jk) * fstravs(ji,jj,jk  ) /( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
+               zws(ji,jj,jk) = - rdttrc(jk) * fstravs(ji,jj,jk+1) /( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
+               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zws(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+
+      ! Surface boudary conditions
+      DO jj = 2, jpjm1        
+         DO ji = fs_2, fs_jpim1
+            zwi(ji,jj,1) = 0.e0
+            zwd(ji,jj,1) = 1. - zws(ji,jj,1) 
+         END DO
+      END DO
+
+      !                                                       ! ===========
       DO jn = 1, jptra                                        ! tracer loop
          !                                                    ! ===========
          IF( l_trdtrc ) ztrtrd(:,:,:) = tra(:,:,:,jn)         ! ??? validation needed
 
-	 ! Initialisation     
-	 zwd( 1 ,:,:) = 0.e0     ;     zwd(jpi,:,:) = 0.e0
-	 zws( 1 ,:,:) = 0.e0     ;     zws(jpi,:,:) = 0.e0
-	 zwi( 1 ,:,:) = 0.e0     ;     zwi(jpi,:,:) = 0.e0
 	 zwt( 1 ,:,:) = 0.e0     ;     zwt(jpi,:,:) = 0.e0     
          zwt(  :,:,1) = 0.e0     ;     zwt(  :,:,jpk) = 0.e0
-         !                                          
-         ! 0. Matrix construction
-         ! ----------------------
-
-         ! Diagonal, inferior, superior
-         ! (including the bottom boundary condition via avs masked
-         DO jk = 1, jpkm1                                                     
-            DO jj = 2, jpjm1                                      
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zwi(ji,jj,jk) = - rdttrc(jk) * fstravs(ji,jj,jk  ) /( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
-                  zws(ji,jj,jk) = - rdttrc(jk) * fstravs(ji,jj,jk+1) /( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
-                  zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zws(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-
-         ! Surface boudary conditions
-         DO jj = 2, jpjm1        
-            DO ji = fs_2, fs_jpim1
-               zwi(ji,jj,1) = 0.e0
-               zwd(ji,jj,1) = 1. - zws(ji,jj,1)
-            END DO
-         END DO
          
          ! Second member construction
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf_iso.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf_iso.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf_iso.F90	(revision 2007)
@@ -182,6 +182,23 @@
 
 
-
-      DO jn = 1, jptra
+      ! 0.2 Update and save of avt (and avs if double diffusive mixing)
+      ! ---------------------------
+
+     DO jj = 2, jpjm1                                 !  Vertical slab
+        !                                             ! ===============
+         DO jk = 2, jpkm1
+            DO ji = 2, jpim1
+               zavi = fsahtw(ji,jj,jk)*( wslpi(ji,jj,jk)*wslpi(ji,jj,jk)   &
+                  &                     +wslpj(ji,jj,jk)*wslpj(ji,jj,jk) )
+               ! add isopycnal vertical coeff. to avs
+               fstravs(ji,jj,jk) = fstravs(ji,jj,jk) + zavi
+            END DO
+         END DO
+       !
+     END DO
+
+
+
+     DO jn = 1, jptra
 
          IF( l_trdtrc ) ztrtrd(:,:,:) = tra(:,:,:,jn)          ! save trends
@@ -262,19 +279,4 @@
             END DO
 
-
-            ! I.3  update and save of avt (and avs if double diffusive mixing)
-            ! ---------------------------
-
-            DO jk = 2, jpkm1
-               DO ji = 2, jpim1
-
-                  zavi = fsahtw(ji,jj,jk)*( wslpi(ji,jj,jk)*wslpi(ji,jj,jk)   &
-                     &                     +wslpj(ji,jj,jk)*wslpj(ji,jj,jk) )
-
-                  ! add isopycnal vertical coeff. to avs
-                  fstravs(ji,jj,jk) = fstravs(ji,jj,jk) + zavi
-
-               END DO
-            END DO
 
 #if defined key_trcldf_eiv
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf_iso_vopt.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf_iso_vopt.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf_iso_vopt.F90	(revision 2007)
@@ -154,5 +154,5 @@
                             zws   => va      ! workspace
       INTEGER, INTENT( in ) ::   kt          ! ocean time-step index
-      INTEGER ::   ji, jj, jk, jn            ! dummy loop indices
+      INTEGER  ::   ji, jj, jk, jn            ! dummy loop indices
       REAL(wp) ::   zavi, zrhs               ! temporary scalars
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
@@ -180,4 +180,89 @@
       ENDIF
 
+         
+      zwd  ( 1, :, : ) = 0.e0    ;     zwd  ( jpi, :,   : ) = 0.e0
+      zws  ( 1, :, : ) = 0.e0    ;     zws  ( jpi, :,   : ) = 0.e0
+      zwi  ( 1, :, : ) = 0.e0    ;     zwi  ( jpi, :,   : ) = 0.e0
+      zwt  ( 1, :, : ) = 0.e0    ;     zwt  ( jpi, :,   : ) = 0.e0
+      zwt  ( :, :, 1 ) = 0.e0    ;     zwt  (   :, :, jpk ) = 0.e0
+      zavsi( 1, :, : ) = 0.e0    ;     zavsi( jpi, :,   : ) = 0.e0 
+      zavsi( :, :, 1 ) = 0.e0    ;     zavsi(   :, :, jpk ) = 0.e0
+
+
+      ! II. Vertical trend associated with the vertical physics
+      !=======================================================
+      !     (including the vertical flux proportional to dk[t] associated
+      !      with the lateral mixing, through the avt update)
+      !     dk[ avt dk[ (t,s) ] ] diffusive trends
+
+      ! II.0 Matrix construction
+      ! ------------------------        
+      ! update and save of avt (and avs if double diffusive mixing)
+      DO jk = 2, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zavi = fsahtw(ji,jj,jk) * (                 &   ! vertical mixing coef. due to lateral mixing
+                  &                           wslpi(ji,jj,jk) * wslpi(ji,jj,jk)      &
+                  &                         + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
+               zavsi(ji,jj,jk) = fstravs(ji,jj,jk) + zavi        ! dd mixing: zavsi = total vertical mixing coef. on tracer
+            END DO
+         END DO
+      END DO
+
+      ! II.1 Vertical diffusion on tracer
+      ! ---------------------------------
+      ! Rebuild the Matrix as avt /= avs
+
+      ! Diagonal, inferior, superior  (including the bottom boundary condition via avs masked)
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zwi(ji,jj,jk) = - rdttrc(jk) * zavsi(ji,jj,jk  ) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
+               zws(ji,jj,jk) = - rdttrc(jk) * zavsi(ji,jj,jk+1) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
+               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zws(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+
+      ! Surface boudary conditions
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            zwi(ji,jj,1) = 0.e0
+            zwd(ji,jj,1) = 1. - zws(ji,jj,1)
+         END DO
+      END DO
+
+      !! Matrix inversion from the first level
+      !!----------------------------------------------------------------------
+      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
+      !
+      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
+      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
+      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
+      !        (        ...               )( ...  ) ( ...  )
+      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
+      !
+      !   m is decomposed in the product of an upper and lower triangular
+      !   matrix
+      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
+      !   The second member is in 2d array zwy
+      !   The solution is in 2d array zwx
+      !   The 3d arry zwt is a work space array
+      !   zwy is used and then used as a work space array : its value is modified!
+
+      ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1
+            zwt(ji,jj,1) = zwd(ji,jj,1)
+         END DO
+      END DO
+      DO jk = 2, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)/zwt(ji,jj,jk-1)
+            END DO
+         END DO
+      END DO
+
       IF( l_trdtrc ) ALLOCATE( ztrtrd(jpi,jpj,jpk) )
 
@@ -187,12 +272,4 @@
          
          IF( l_trdtrc ) ztrtrd(:,:,:) = tra(:,:,:,jn)          ! save trends
-         
-         zwd  ( 1, :, : ) = 0.e0    ;     zwd  ( jpi, :,   : ) = 0.e0
-         zws  ( 1, :, : ) = 0.e0    ;     zws  ( jpi, :,   : ) = 0.e0
-         zwi  ( 1, :, : ) = 0.e0    ;     zwi  ( jpi, :,   : ) = 0.e0
-         zwt  ( 1, :, : ) = 0.e0    ;     zwt  ( jpi, :,   : ) = 0.e0
-         zwt  ( :, :, 1 ) = 0.e0    ;     zwt  (   :, :, jpk ) = 0.e0
-         zavsi( 1, :, : ) = 0.e0    ;     zavsi( jpi, :,   : ) = 0.e0 
-         zavsi( :, :, 1 ) = 0.e0    ;     zavsi(   :, :, jpk ) = 0.e0
 
 #  if defined key_trc_diatrd
@@ -200,83 +277,4 @@
          ztrd(:,:,:) = tra(:,:,:,jn)
 #  endif
-
-         ! II. Vertical trend associated with the vertical physics
-         ! =======================================================
-         !     (including the vertical flux proportional to dk[t] associated
-         !      with the lateral mixing, through the avt update)
-         !     dk[ avt dk[ (t,s) ] ] diffusive trends
-
-
-         ! II.0 Matrix construction
-         ! ------------------------        
-         ! update and save of avt (and avs if double diffusive mixing)
-         DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zavi = fsahtw(ji,jj,jk) * (                 &   ! vertical mixing coef. due to lateral mixing
-                     &                           wslpi(ji,jj,jk) * wslpi(ji,jj,jk)      &
-                     &                         + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
-                  zavsi(ji,jj,jk) = fstravs(ji,jj,jk) + zavi        ! dd mixing: zavsi = total vertical mixing coef. on tracer
-
-               END DO
-            END DO
-         END DO
-
-
-         ! II.1 Vertical diffusion on tracer
-         ! ---------------------------------
-
-         ! Rebuild the Matrix as avt /= avs
-
-         ! Diagonal, inferior, superior  (including the bottom boundary condition via avs masked)
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zwi(ji,jj,jk) = - rdttrc(jk) * zavsi(ji,jj,jk  ) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
-                  zws(ji,jj,jk) = - rdttrc(jk) * zavsi(ji,jj,jk+1) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
-                  zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zws(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-
-         ! Surface boudary conditions
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zwi(ji,jj,1) = 0.e0
-               zwd(ji,jj,1) = 1. - zws(ji,jj,1)
-            END DO
-         END DO
-
-         !! Matrix inversion from the first level
-         !!----------------------------------------------------------------------
-         !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
-         !
-         !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
-         !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
-         !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
-         !        (        ...               )( ...  ) ( ...  )
-         !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
-         !
-         !   m is decomposed in the product of an upper and lower triangular
-         !   matrix
-         !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
-         !   The second member is in 2d array zwy
-         !   The solution is in 2d array zwx
-         !   The 3d arry zwt is a work space array
-         !   zwy is used and then used as a work space array : its value is modified!
-
-         ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-               zwt(ji,jj,1) = zwd(ji,jj,1)
-            END DO
-         END DO
-         DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                  zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1)
-               END DO
-            END DO
-         END DO
 
          ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/par_trc.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/par_trc.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/par_trc.F90	(revision 2007)
@@ -18,6 +18,6 @@
    USE par_lobster   ! LOBSTER model
    USE par_pisces    ! PISCES  model
+   USE par_c14b      ! C14 bomb tracer
    USE par_cfc       ! CFC 11 and 12 tracers
-   USE par_c14b      ! C14 bomb tracer 
    USE par_my_trc    ! user defined passive tracers
 
@@ -27,9 +27,9 @@
    ! Passive tracers : Total size
    ! ---------------               ! total number of passive tracers, of 2d and 3d output and trend arrays
-   INTEGER, PUBLIC, PARAMETER ::   jptra    =  jp_lobster    + jp_pisces     + jp_cfc     + jp_c14b     + jp_my_trc
-   INTEGER, PUBLIC, PARAMETER ::   jpdia2d  =  jp_lobster_2d + jp_pisces_2d  + jp_cfc_2d  + jp_c14b_2d  + jp_my_trc_2d
-   INTEGER, PUBLIC, PARAMETER ::   jpdia3d  =  jp_lobster_3d + jp_pisces_3d  + jp_cfc_3d  + jp_c14b_3d  + jp_my_trc_3d
+   INTEGER, PUBLIC, PARAMETER ::   jptra    =  jp_lobster    + jp_pisces     + jp_cfc     + jp_my_trc
+   INTEGER, PUBLIC, PARAMETER ::   jpdia2d  =  jp_lobster_2d + jp_pisces_2d  + jp_cfc_2d  + jp_my_trc_2d
+   INTEGER, PUBLIC, PARAMETER ::   jpdia3d  =  jp_lobster_3d + jp_pisces_3d  + jp_cfc_3d  + jp_my_trc_3d
    !                     ! total number of sms diagnostic arrays
-   INTEGER, PUBLIC, PARAMETER ::   jpdiabio = jp_lobster_trd + jp_pisces_trd + jp_cfc_trd + jp_c14b_trd + jp_my_trc_trd
+   INTEGER, PUBLIC, PARAMETER ::   jpdiabio = jp_lobster_trd + jp_pisces_trd + jp_cfc_trd + jp_my_trc_trd
    
    !  1D configuration ("key_c1d")
@@ -40,18 +40,34 @@
    LOGICAL, PUBLIC, PARAMETER ::   lk_trc_c1d   = .FALSE.  !: 1D pass. tracer configuration flag
 # endif
-
    ! Passive tracers : size for TRP trends diagnotics (used if 'key_trc_diatrd' defined)
-#if defined key_trcldf_eiv
-# if defined key_trcdmp
-   INTEGER, PARAMETER :: jpdiatrc = 11      !: trends: 3*(advection + diffusion + eiv ) + damping + sms
-# else
-   INTEGER, PARAMETER :: jpdiatrc = 10      !: trends: 3*(advection + diffusion + eiv )           + sms
-# endif
-#else
-# if defined key_trcdmp
-   INTEGER, PARAMETER :: jpdiatrc =  8      !: trends: 3*(advection + diffusion       ) + damping + sms
-# else
-   INTEGER, PARAMETER :: jpdiatrc =  7      !: trends: 3*(advection + diffusion       ) + damping + sms
-# endif
+# if defined key_trc_diatrd
+   ! Passive tracers : size for TRP trends diagnotics (used if 'key_trc_diatrd' defined)
+   INTEGER, PARAMETER ::   jptrc_xad     =  1   !: x- horizontal advection
+   INTEGER, PARAMETER ::   jptrc_yad     =  2   !: y- horizontal advection
+   INTEGER, PARAMETER ::   jptrc_zad     =  3   !: z- vertical   advection
+   INTEGER, PARAMETER ::   jptrc_xdf     =  4   !: lateral       diffusion
+   INTEGER, PARAMETER ::   jptrc_ydf     =  5   !: lateral       diffusion
+   INTEGER, PARAMETER ::   jptrc_zdf     =  6   !: vertical diffusion (Kz)
+   INTEGER, PARAMETER ::   jptrc_sbc     =  7   !: surface boundary condition
+#if ! defined key_trcldf_eiv && ! defined key_trcdmp
+   INTEGER, PARAMETER ::   jpdiatrc      =  7  !: trends: 3*(advection + diffusion       ) + sbc
+#endif
+#if defined key_trcldf_eiv && defined key_trcdmp
+   INTEGER, PARAMETER ::   jptrc_xei     =  8   !: x- horiz. EIV advection
+   INTEGER, PARAMETER ::   jptrc_yei     =  9   !: y- horiz. EIV advection
+   INTEGER, PARAMETER ::   jptrc_zei     = 10   !: z- vert.  EIV advection
+   INTEGER, PARAMETER ::   jptrc_dmp     = 11   !: damping
+   INTEGER, PARAMETER ::   jpdiatrc      = 11   !: trends: 3*(advection + diffusion + eiv ) + sbc + damping
+#endif
+#if defined key_trcldf_eiv && ! defined key_trcdmp
+   INTEGER, PARAMETER ::   jptrc_xei     =  8   !: x- horiz. EIV advection
+   INTEGER, PARAMETER ::   jptrc_yei     =  9   !: y- horiz. EIV advection
+   INTEGER, PARAMETER ::   jptrc_zei     = 10   !: z- vert.  EIV advection
+   INTEGER, PARAMETER ::   jpdiatrc      = 10   !: trends: 3*(advection + diffusion + eiv ) + sbc 
+#endif
+#if ! defined key_trcldf_eiv && defined key_trcdmp
+   INTEGER, PARAMETER ::   jptrc_dmp     =  8   !: damping
+   INTEGER, PARAMETER ::   jpdiatrc      =  8   !: trends: 3*(advection + diffusion       ) + sbc + damping
+#endif
 #endif
 
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/trcdia.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/trcdia.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/trcdia.F90	(revision 2007)
@@ -25,4 +25,5 @@
    USE trc
    USE trp_trc
+   USE par_trc
    USE trdmld_trc_oce, ONLY : luttrd
    USE dianam    ! build name of file (routine)
@@ -41,5 +42,5 @@
    INTEGER  ::   ndimt50   !: number of ocean points in index array
    INTEGER  ::   ndimt51   !: number of ocean points in index array
-   REAL(wp) ::   xjulian   !: ????   not DOCTOR !
+   REAL(wp) ::   zjulian   !: ????   not DOCTOR !
    INTEGER , DIMENSION (jpij*jpk) ::   ndext50   !: integer arrays for ocean 3D index
    INTEGER , DIMENSION (jpij)     ::   ndext51   !: integer arrays for ocean surface index
@@ -157,10 +158,10 @@
 
          ! Compute julian date from starting date of the run
-         CALL ymds2ju( nyear, nmonth, nday, rdt, xjulian )
-         xjulian = xjulian - adatrj   !   set calendar origin to the beginning of the experiment
+         CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
+         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
          IF(lwp)WRITE(numout,*)' '  
          IF(lwp)WRITE(numout,*)' Date 0 used :', nittrc000                         &
             &                 ,' YEAR ', nyear, ' MONTH ', nmonth, ' DAY ', nday   &
-            &                 ,'Julian day : ', xjulian  
+            &                 ,'Julian day : ', zjulian  
   
          IF(lwp) WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma,  &
@@ -171,5 +172,5 @@
          IF(lwp) THEN
             CALL dia_nam( clhstnam, nwritetrc,' ' )
-            CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
+            CALL ctlopn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', 1, numout, lwp, narea )
             WRITE(inum,*) clhstnam
             CLOSE(inum)
@@ -184,5 +185,5 @@
          CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,     &
             &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,         & 
-            &          nittrc000-ndttrc, xjulian, zdt, nhorit5, nit5 , domain_id=nidom)
+            &          nittrc000-ndttrc, zjulian, zdt, nhorit5, nit5 , domain_id=nidom)
 
          ! Vertical grid for tracer : gdept
@@ -258,5 +259,5 @@
       CHARACTER (len=80) ::   cltral
       CHARACTER (len=10) ::   csuff
-      INTEGER  ::   jn, jl
+      INTEGER  ::   jn, jl, ikn
       INTEGER  ::   iimi, iima, ijmi, ijma, ipk, it, itmod
       REAL(wp) ::   zsto, zout, zdt
@@ -313,5 +314,5 @@
                CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,       &
                   &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,   &
-                  &          nittrc000-ndttrc, xjulian, zdt, nhorit6(jn),  &
+                  &          nittrc000-ndttrc, zjulian, zdt, nhorit6(jn),  &
                   &          nit6(jn) , domain_id=nidom )
 
@@ -322,80 +323,68 @@
 
           ! Declare all the output fields as NETCDF variables
-
-          ! trends for tracer concentrations
           DO jn = 1, jptra
             IF( luttrd(jn) ) THEN
                 DO jl = 1, jpdiatrc
-                  IF( jl == 1 ) THEN
+                  IF( jl == jptrc_xad ) THEN
                       ! short and long title for x advection for tracer
                       WRITE (cltra,'("XAD_",16a)') ctrcnm(jn)
-                      WRITE (cltral,'("X advective trend for ",58a)')  &
-                         &      ctrcnl(jn)(1:58)
-                  END IF
-                  IF( jl == 2 ) THEN
+                      WRITE (cltral,'("X advective trend for ",58a)') ctrcnl(jn)(1:58)
+                  END IF
+                  IF( jl == jptrc_yad ) THEN
                       ! short and long title for y advection for tracer
                       WRITE (cltra,'("YAD_",16a)') ctrcnm(jn)
-                      WRITE (cltral,'("Y advective trend for ",58a)')  &
-                         &      ctrcnl(jn)(1:58)
-                  END IF
-                  IF( jl == 3 ) THEN
+                      WRITE (cltral,'("Y advective trend for ",58a)') ctrcnl(jn)(1:58)
+                  END IF
+                  IF( jl == jptrc_zad ) THEN
                       ! short and long title for Z advection for tracer
                       WRITE (cltra,'("ZAD_",16a)') ctrcnm(jn)
-                      WRITE (cltral,'("Z advective trend for ",58a)')  &
-                         &      ctrcnl(jn)(1:58)
-                  END IF
-                  IF( jl == 4 ) THEN
+                      WRITE (cltral,'("Z advective trend for ",58a)') ctrcnl(jn)(1:58)
+                  END IF
+                  IF( jl == jptrc_xdf ) THEN
                       ! short and long title for X diffusion for tracer
                       WRITE (cltra,'("XDF_",16a)') ctrcnm(jn)
-                      WRITE (cltral,'("X diffusion trend for ",58a)')  &
-                         &      ctrcnl(jn)(1:58)
-                  END IF
-                  IF( jl == 5 ) THEN
+                      WRITE (cltral,'("X diffusion trend for ",58a)') ctrcnl(jn)(1:58)
+                  END IF
+                  IF( jl == jptrc_ydf ) THEN
                       ! short and long title for Y diffusion for tracer
                       WRITE (cltra,'("YDF_",16a)') ctrcnm(jn)
-                      WRITE (cltral,'("Y diffusion trend for ",58a)')  &
-                         &      ctrcnl(jn)(1:58)
-                  END IF
-                  IF( jl == 6 ) THEN
+                      WRITE (cltral,'("Y diffusion trend for ",58a)') ctrcnl(jn)(1:58)
+                  END IF
+                  IF( jl == jptrc_zdf ) THEN
                       ! short and long title for Z diffusion for tracer
                       WRITE (cltra,'("ZDF_",16a)') ctrcnm(jn)
-                      WRITE (cltral,'("Z diffusion trend for ",58a)')  &
-                         &      ctrcnl(jn)(1:58)
+                      WRITE (cltral,'("Z diffusion trend for ",58a)') ctrcnl(jn)(1:58)
                   END IF
 # if defined key_trcldf_eiv
-                  IF( jl == 7 ) THEN
+                  IF( jl == jptrc_xei ) THEN
                       ! short and long title for x gent velocity for tracer
                       WRITE (cltra,'("XGV_",16a)') ctrcnm(jn)
-                      WRITE (cltral,'("X gent velocity trend for ",53a)')  &
-                         &      ctrcnl(jn)(1:53)
-                  END IF
-                  IF( jl == 8 ) THEN
+                      WRITE (cltral,'("X gent velocity trend for ",53a)') ctrcnl(jn)(1:53)
+                  END IF
+                  IF( jl == jptrc_yei ) THEN
                       ! short and long title for y gent velocity for tracer
                       WRITE (cltra,'("YGV_",16a)') ctrcnm(jn)
-                      WRITE (cltral,'("Y gent velocity trend for ",53a)')  &
-                         &      ctrcnl(jn)(1:53)
-                  END IF
-                  IF( jl == 9 ) THEN
+                      WRITE (cltral,'("Y gent velocity trend for ",53a)') ctrcnl(jn)(1:53)
+                  END IF
+                  IF( jl == jptrc_zei ) THEN
                       ! short and long title for Z gent velocity for tracer
                       WRITE (cltra,'("ZGV_",16a)') ctrcnm(jn)
-                      WRITE (cltral,'("Z gent velocity trend for ",53a)')  &
-                         &      ctrcnl(jn)(1:53)
+                      WRITE (cltral,'("Z gent velocity trend for ",53a)') ctrcnl(jn)(1:53)
                   END IF
 # endif
 # if defined key_trcdmp
-                  IF( jl == jpdiatrc - 1 ) THEN
+                  IF( jl == jptrc_dmp ) THEN
                       ! last trends for tracer damping : short and long title
                       WRITE (cltra,'("TDM_",16a)') ctrcnm(jn)
-                      WRITE (cltral,'("Tracer damping trend for ",55a)')  &
-                         &      ctrcnl(jn)(1:55)
-                  END IF
-# endif
-                  IF( jl == jpdiatrc ) THEN
+                      WRITE (cltral,'("Tracer damping trend for ",55a)') ctrcnl(jn)(1:55)
+                  END IF
+# endif
+                  IF( jl == jptrc_sbc ) THEN
                       ! last trends for tracer damping : short and long title
                       WRITE (cltra,'("SBC_",16a)') ctrcnm(jn)
-                      WRITE (cltral,'("Surface boundary flux ",58a)')  &
-                      &      ctrcnl(jn)(1:58)
-                  END IF
-
+                      WRITE (cltral,'("Surface boundary flux ",58a)') ctrcnl(jn)(1:55)
+                  END IF
+                      WRITE (cltral,'("Surface boundary flux ",58a)') ctrcnl(jn)(1:55)
+                  END IF
                   CALL FLUSH( numout )
                   cltrau = ctrcun(jn)      ! UNIT for tracer /trends
@@ -406,5 +395,4 @@
             END IF
          END DO
-
          ! CLOSE netcdf Files
           DO jn = 1, jptra
@@ -432,23 +420,24 @@
       DO jn = 1, jptra
          IF( luttrd(jn) ) THEN
+            ikn = ikeep(jn) 
             DO jl = 1, jpdiatrc
-               ! short titles 
-               IF( jl == 1)   WRITE (cltra,'("XAD_",16a)') ctrcnm(jn)      ! x advection for tracer
-               IF( jl == 2)   WRITE (cltra,'("YAD_",16a)') ctrcnm(jn)      ! z advection for tracer
-               IF( jl == 3)   WRITE (cltra,'("ZAD_",16a)') ctrcnm(jn)      ! z advection for tracer
-               IF( jl == 4)   WRITE (cltra,'("XDF_",16a)') ctrcnm(jn)      ! x diffusion for tracer
-               IF( jl == 5)   WRITE (cltra,'("YDF_",16a)') ctrcnm(jn)      ! y diffusion for tracer
-               IF( jl == 6)   WRITE (cltra,'("ZDF_",16a)') ctrcnm(jn)      ! z diffusion for tracer
+               ! short titles
+               IF( jl == jptrc_xad)   WRITE (cltra,'("XAD_",16a)') ctrcnm(jn)
+               IF( jl == jptrc_yad)   WRITE (cltra,'("YAD_",16a)') ctrcnm(jn)
+               IF( jl == jptrc_zad)   WRITE (cltra,'("ZAD_",16a)') ctrcnm(jn)
+               IF( jl == jptrc_xdf)   WRITE (cltra,'("XDF_",16a)') ctrcnm(jn)
+               IF( jl == jptrc_ydf)   WRITE (cltra,'("YDF_",16a)') ctrcnm(jn)
+               IF( jl == jptrc_zdf)   WRITE (cltra,'("ZDF_",16a)') ctrcnm(jn)
 # if defined key_trcldf_eiv
-               IF( jl == 7)   WRITE (cltra,'("XGV_",16a)') ctrcnm(jn)      ! x gent velocity for tracer
-               IF( jl == 8)   WRITE (cltra,'("YGV_",16a)') ctrcnm(jn)      ! y gent velocity for tracer
-               IF( jl == 9)   WRITE (cltra,'("ZGV_",16a)') ctrcnm(jn)      ! z gent velocity for tracer
+               IF( jl == jptrc_xei)   WRITE (cltra,'("XGV_",16a)') ctrcnm(jn)
+               IF( jl == jptrc_yei)   WRITE (cltra,'("YGV_",16a)') ctrcnm(jn)
+               IF( jl == jptrc_zei)   WRITE (cltra,'("ZGV_",16a)') ctrcnm(jn)
 # endif
 # if defined key_trcdmp
-               IF( jl == jpdiatrc - 1 )   WRITE (cltra,'("TDM_",16a)') ctrcnm(jn)      ! damping
-# endif
-               IF( jl == jpdiatrc )   WRITE (cltra,'("SBC_",a)') ctrcnm(jn)      ! surface boundary conditions
+               IF( jl == jptrc_dmp )  WRITE (cltra,'("TDM_",16a)') ctrcnm(jn)
+# endif
+               IF( jl == jptrc_sbc )  WRITE (cltra,'("SBC_",16a)') ctrcnm(jn)
                !
-               CALL histwrite(nit6(jn), cltra, it, trtrd(:,:,:,ikeep(jn),jl),ndimt50, ndext50)
+               CALL histwrite(nit6(jn), cltra, it, trtrd(:,:,:,ikn,jl),ndimt50, ndext50)
             END DO
          END IF
@@ -552,5 +541,5 @@
          CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,             &
             &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,         &
-            &          nittrc000-ndttrc, xjulian, zdt, nhoritd, nitd , domain_id=nidom )
+            &          nittrc000-ndttrc, zjulian, zdt, nhoritd, nitd , domain_id=nidom )
 
          ! Vertical grid for 2d and 3d arrays
@@ -700,5 +689,5 @@
          CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,      &
             &    iimi, iima-iimi+1, ijmi, ijma-ijmi+1,          &
-            &    nittrc000-ndttrc, xjulian, zdt, nhoritb, nitb , domain_id=nidom )
+            &    nittrc000-ndttrc, zjulian, zdt, nhoritb, nitb , domain_id=nidom )
          ! Vertical grid for biological trends
          CALL histvert(nitb, 'deptht', 'Vertical T levels', 'm', ipk, gdept_0, ndepitb)
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/trcdta.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/trcdta.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/trcdta.F90	(revision 2007)
@@ -25,4 +25,5 @@
    PUBLIC trc_dta   ! called in trcini.F90 and trcdmp.F90
 
+   LOGICAL , PUBLIC, PARAMETER ::   lk_dtatrc = .TRUE.   !: temperature data flag
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk,jptra) ::   trdta   !: tracer data at given time-step
 
@@ -62,5 +63,6 @@
       !!
       CHARACTER (len=39) ::   clname(jptra)
-      INTEGER, PARAMETER ::   jpmois  = 12        ! number of months
+      INTEGER, PARAMETER ::   &
+         jpmonth = 12    ! number of months
       INTEGER ::   ji, jj, jn, jl 
       INTEGER ::   imois, iman, i15, ik  ! temporary integers 
@@ -81,5 +83,5 @@
             ENDIF
             ! Initialization
-            iman = jpmois
+            iman = jpmonth
             i15  = nday / 16
             imois = nmonth + i15 -1
@@ -188,4 +190,5 @@
             ! Read init file only
             IF( kt == nittrc000  ) THEN
+               ntrc1(jn) = 1
                CALL iom_get ( numtr(jn), jpdom_data, ctrcnm(jn), trdta(:,:,:,jn), ntrc1(jn) )
                trdta(:,:,:,jn) = trdta(:,:,:,jn) * tmask(:,:,:)
@@ -204,4 +207,5 @@
    !!   Dummy module                              NO 3D passive tracer data
    !!----------------------------------------------------------------------
+   LOGICAL , PUBLIC, PARAMETER ::   lk_dtatrc = .FALSE.   !: temperature data flag
 CONTAINS
    SUBROUTINE trc_dta( kt )        ! Empty routine
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/trcini.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 2007)
@@ -122,8 +122,7 @@
          trb(:,:,:,:) = trn(:,:,:,:)
       ELSE
+         !
          CALL trc_rst_read      ! restart from a file
-#if defined key_off_tra
-         CALL day_init          ! calendar
-#endif
+         !
       ENDIF
 
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/trcrst.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 2007)
@@ -1,10 +1,19 @@
 MODULE trcrst
    !!======================================================================
-   !!                       ***  MODULE trcrst  ***
-   !! TOP :   create, write, read the restart files for passive tracers
+   !!                         ***  MODULE trcrst  ***
+   !! TOP :   Manage the passive tracer restart
    !!======================================================================
-   !! History :   1.0  !  2007-02 (C. Ethe) adaptation from the ocean
+   !! History :    -   !  1991-03  ()  original code
+   !!             1.0  !  2005-03 (O. Aumont, A. El Moussaoui) F90
+   !!              -   !  2005-10 (C. Ethe) print control
+   !!             2.0  !  2005-10 (C. Ethe, G. Madec) revised architecture
    !!----------------------------------------------------------------------
 #if defined key_top
+   !!----------------------------------------------------------------------
+   !!   'key_top'                                                TOP models
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   trc_rst :   Restart for passive tracer
+   !!----------------------------------------------------------------------
    !!----------------------------------------------------------------------
    !!   'key_top'                                                TOP models
@@ -16,22 +25,23 @@
    USE oce_trc
    USE trc
-   USE sms_lobster         ! LOBSTER variables
-   USE sms_pisces          ! PISCES variables
-   USE trcsms_cfc          ! CFC variables
-   USE trcsms_c14b         ! C14 variables
-   USE trcsms_my_trc       ! MY_TRC variables
-   USE trctrp_lec   
+   USE trctrp_lec
    USE lib_mpp
    USE iom
-   
+   USE trcrst_cfc      ! CFC      
+   USE trcrst_lobster  ! LOBSTER  restart
+   USE trcrst_pisces   ! PISCES   restart
+   USE trcrst_c14b     ! C14 bomb restart
+   USE trcrst_my_trc   ! MY_TRC   restart
+#if defined key_off_tra
+    USE daymod
+#endif
    IMPLICIT NONE
    PRIVATE
-   
+
    PUBLIC   trc_rst_opn       ! called by ???
    PUBLIC   trc_rst_read      ! called by ???
    PUBLIC   trc_rst_wri       ! called by ???
-   
+
    INTEGER, PUBLIC ::   numrtr, numrtw   !: logical unit for trc restart (read and write)
-
 
    !! * Substitutions
@@ -89,6 +99,5 @@
    END SUBROUTINE trc_rst_opn
 
-
-   SUBROUTINE trc_rst_read 
+   SUBROUTINE trc_rst_read
       !!----------------------------------------------------------------------
       !!                    ***  trc_rst_opn  ***
@@ -96,14 +105,9 @@
       !! ** purpose  :   read passive tracer fields in restart files
       !!----------------------------------------------------------------------
-      INTEGER  ::  jn  
-      INTEGER  ::  iarak0
+      INTEGER  ::  jn     
+      INTEGER  ::  iarak0 
       REAL(wp) ::  zarak0
       INTEGER  ::  jlibalt = jprstlib
       LOGICAL  ::  llok
-#if defined key_pisces 
-      INTEGER  ::  ji, jj, jk
-      REAL(wp) ::  zcaralk, zbicarb, zco3
-      REAL(wp) ::  ztmas, ztmas1
-#endif
 
       !!----------------------------------------------------------------------
@@ -115,10 +119,10 @@
       IF ( jprstlib == jprstdimg ) THEN
         ! eventually read netcdf file (monobloc)  for restarting on different number of processors
-        ! if {cn_trcrst_in}.nc exists, then set jlibalt to jpnf90
+        ! if {cn_trcrst_in}.nc exists, then set jlibalt to jpnf90 
         INQUIRE( FILE = TRIM(cn_trcrst_in)//'.nc', EXIST = llok )
-        IF ( llok ) THEN ; jlibalt = jpnf90  ; ELSE ; jlibalt = jprstlib ; ENDIF
-      ENDIF
-      
-      CALL iom_open( cn_trcrst_in, numrtr, kiolib = jlibalt )
+        IF ( llok ) THEN ; jlibalt = jpnf90  ; ELSE ; jlibalt = jprstlib ; ENDIF 
+      ENDIF
+
+      CALL iom_open( cn_trcrst_in, numrtr, kiolib = jlibalt ) 
 
       ! Time domain : restart
@@ -136,65 +140,24 @@
          & ' centered or euler '  )
       IF(lwp) WRITE(numout,*)
-
       IF(lwp) WRITE(numout,*) '    arakawa option      : ', NINT( zarak0 )
 
-
       ! READ prognostic variables and computes diagnostic variable
       DO jn = 1, jptra
-         CALL iom_get( numrtr, jpdom_autoglo, 'TRN'//ctrcnm(jn), trn(:,:,:,jn) ) 
-      END DO
-
-      DO jn = 1, jptra
-         CALL iom_get( numrtr, jpdom_autoglo, 'TRB'//ctrcnm(jn), trb(:,:,:,jn) ) 
-      END DO
-
-#if defined key_lobster
-      CALL iom_get( numrtr, jpdom_autoglo, 'SEDB'//ctrcnm(jpdet), sedpocb(:,:) ) 
-      CALL iom_get( numrtr, jpdom_autoglo, 'SEDN'//ctrcnm(jpdet), sedpocn(:,:) ) 
-#endif
-#if defined key_pisces
-      !
-      IF( ln_pisdmp ) CALL pis_dmp_ini  ! relaxation of some tracers
-      !
-      IF( iom_varid( numrtr, 'PH', ldstop = .FALSE. ) > 0 ) THEN
-         CALL iom_get( numrtr, jpdom_autoglo, 'PH' , hi(:,:,:)  )
-      ELSE
-         ! Set PH from  total alkalinity, borat (???), akb3 (???) and ak23 (???)
-         ! --------------------------------------------------------
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ztmas   = tmask(ji,jj,jk)
-                  ztmas1  = 1. - tmask(ji,jj,jk)
-                  zcaralk = trn(ji,jj,jk,jptal) - borat(ji,jj,jk) / (  1. + 1.E-8 / ( rtrn + akb3(ji,jj,jk) )  )
-                  zco3    = ( zcaralk - trn(ji,jj,jk,jpdic) ) * ztmas + 0.5e-3 * ztmas1
-                  zbicarb = ( 2. * trn(ji,jj,jk,jpdic) - zcaralk )
-                  hi(ji,jj,jk) = ( ak23(ji,jj,jk) * zbicarb / zco3 ) * ztmas + 1.e-9 * ztmas1
-               END DO
-            END DO
-         END DO
-      ENDIF
-      CALL iom_get( numrtr, jpdom_autoglo, 'Silicalim', xksi(:,:) ) 
-      IF( iom_varid( numrtr, 'Silicamax', ldstop = .FALSE. ) > 0 ) THEN
-         CALL iom_get( numrtr, jpdom_autoglo, 'Silicamax' , xksimax(:,:)  )
-      ELSE
-         xksimax(:,:) = xksi(:,:)
-      ENDIF
-#endif
-#if defined key_cfc
-      DO jn = jp_cfc0, jp_cfc1
-         CALL iom_get( numrtr, jpdom_autoglo, 'qint_'//ctrcnm(jn), qint_cfc(:,:,jn) ) 
-      END DO
-#endif
-#if defined key_c14b
-      CALL iom_get( numrtr, jpdom_autoglo, 'qint_'//ctrcnm(jn) , qint_c14(:,:) ) 
-#endif
-#if defined key_my_trc
-#endif
-      
+         CALL iom_get( numrtr, jpdom_autoglo, 'TRN'//ctrcnm(jn), trn(:,:,:,jn) )
+      END DO
+
+      DO jn = 1, jptra
+         CALL iom_get( numrtr, jpdom_autoglo, 'TRB'//ctrcnm(jn), trb(:,:,:,jn) )
+      END DO
+
+      IF( lk_lobster )   CALL trc_rst_read_lobster( numrtr )      ! LOBSTER bio-model
+      IF( lk_pisces  )   CALL trc_rst_read_pisces ( numrtr )      ! PISCES  bio-model
+      IF( lk_cfc     )   CALL trc_rst_read_cfc    ( numrtr )      ! CFC     tracers
+      IF( lk_c14b    )   CALL trc_rst_read_c14b   ( numrtr )      ! C14 bomb  tracer
+      IF( lk_my_trc  )   CALL trc_rst_read_my_trc ( numrtr )      ! MY_TRC  tracers
+
       CALL iom_close( numrtr )
       !
    END SUBROUTINE trc_rst_read
-
 
    SUBROUTINE trc_rst_wri( kt )
@@ -218,6 +181,6 @@
       CALL iom_rstput( kt, nitrst, numrtw, 'arak0', zarak0 )
 
-      ! prognostic variables
-      ! --------------------
+      ! prognostic variables 
+      ! -------------------- 
       DO jn = 1, jptra
          CALL iom_rstput( kt, nitrst, numrtw, 'TRN'//ctrcnm(jn), trn(:,:,:,jn) )
@@ -228,24 +191,10 @@
       END DO
 
-#if defined key_lobster
-         CALL iom_rstput( kt, nitrst, numrtw, 'SEDB'//ctrcnm(jpdet), sedpocb(:,:) )
-         CALL iom_rstput( kt, nitrst, numrtw, 'SEDN'//ctrcnm(jpdet), sedpocn(:,:) )
-#endif
-#if defined key_pisces 
-         CALL iom_rstput( kt, nitrst, numrtw, 'PH', hi(:,:,:) )
-         CALL iom_rstput( kt, nitrst, numrtw, 'Silicalim', xksi(:,:) )
-         CALL iom_rstput( kt, nitrst, numrtw, 'Silicamax', xksimax(:,:) )
-#endif
-#if defined key_cfc
-         DO jn = jp_cfc0, jp_cfc1
-            CALL iom_rstput( kt, nitrst, numrtw, 'qint_'//ctrcnm(jn), qint_cfc(:,:,jn) )
-         END DO
-#endif
-#if defined key_c14b
-         CALL iom_rstput( kt, nitrst, numrtw, 'qint_'//ctrcnm(jn), qint_c14(:,:) )
-#endif
-#if defined key_my_trc
-#endif
-      
+      IF( lk_lobster )   CALL trc_rst_wri_lobster( kt, nitrst, numrtw )      ! LOBSTER bio-model
+      IF( lk_pisces  )   CALL trc_rst_wri_pisces ( kt, nitrst, numrtw )      ! PISCES  bio-model
+      IF( lk_cfc     )   CALL trc_rst_wri_cfc    ( kt, nitrst, numrtw )      ! CFC     tracers
+      IF( lk_c14b    )   CALL trc_rst_wri_c14b   ( kt, nitrst, numrtw )      ! C14 bomb  tracer
+      IF( lk_my_trc  )   CALL trc_rst_wri_my_trc ( kt, nitrst, numrtw )      ! MY_TRC  tracers
+
       IF( kt == nitrst ) THEN
           CALL trc_rst_stat            ! statistics
@@ -256,5 +205,6 @@
       ENDIF
       !
-   END SUBROUTINE trc_rst_wri
+   END SUBROUTINE trc_rst_wri 
+
 
    SUBROUTINE trc_rst_cal( kt, cdrw )
@@ -329,4 +279,7 @@
            WRITE(numout,*)
          ENDIF
+         !
+         CALL day_init          ! compute calendar
+         !
 #endif
 
@@ -347,117 +300,4 @@
    END SUBROUTINE trc_rst_cal
 
-# if defined key_pisces 
-
-   SUBROUTINE pis_dmp_ini 
-      !!----------------------------------------------------------------------
-      !!                    ***  pis_dmp_ini  ***
-      !!
-      !! ** purpose  : Relaxation of some tracers
-      !!----------------------------------------------------------------------
-      INTEGER  :: ji, jj, jk  
-      REAL(wp) ::  &
-         alkmean = 2426. ,  & ! mean value of alkalinity ( Glodap ; for Goyet 2391. )
-         po4mean = 2.165 ,  & ! mean value of phosphates
-         no3mean = 30.90 ,  & ! mean value of nitrate
-         siomean = 91.51      ! mean value of silicate
-      
-      REAL(wp) ::   zvol, ztrasum
-
-
-      IF(lwp)  WRITE(numout,*)
-
-      IF( cp_cfg == "orca" .AND. .NOT. lk_trc_c1d ) THEN      ! ORCA condiguration (not 1D) !
-         !                                                    ! --------------------------- !
-         ! set total alkalinity, phosphate, NO3 & silicate
-
-         ! total alkalinity
-         ztrasum = 0.e0             
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zvol = cvol(ji,jj,jk)
-#  if defined key_off_degrad
-                  zvol = zvol * facvol(ji,jj,jk)
-#  endif
-                  ztrasum = ztrasum + trn(ji,jj,jk,jptal) * zvol
-               END DO
-            END DO
-         END DO
-         IF( lk_mpp )   CALL mpp_sum( ztrasum )     ! sum over the global domain  
-         
-         ztrasum = ztrasum / areatot * 1.e6
-         IF(lwp) WRITE(numout,*) '       TALK mean : ', ztrasum
-         trn(:,:,:,jptal) = trn(:,:,:,jptal) * alkmean / ztrasum
-            
-         ! phosphate
-         ztrasum = 0.e0
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zvol = cvol(ji,jj,jk)
-#  if defined key_off_degrad
-                  zvol = zvol * facvol(ji,jj,jk)
-#  endif
-                  ztrasum = ztrasum + trn(ji,jj,jk,jppo4) * zvol
-               END DO
-            END DO
-         END DO
-         IF( lk_mpp )   CALL mpp_sum( ztrasum )     ! sum over the global domain  
-         
-         ztrasum = ztrasum / areatot * 1.e6 / 122.
-         IF(lwp) WRITE(numout,*) '       PO4  mean : ', ztrasum
-         trn(:,:,:,jppo4) = trn(:,:,:,jppo4) * po4mean / ztrasum
-        
-         ! Nitrates          
-         ztrasum = 0.e0
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zvol = cvol(ji,jj,jk)
-#  if defined key_off_degrad
-                  zvol = zvol * facvol(ji,jj,jk)
-#  endif
-                  ztrasum = ztrasum + trn(ji,jj,jk,jpno3) * zvol
-               END DO
-            END DO
-         END DO
-         IF( lk_mpp )   CALL mpp_sum( ztrasum )     ! sum over the global domain  
-         
-         ztrasum = ztrasum / areatot * 1.e6 / 7.6
-         IF(lwp) WRITE(numout,*) '       NO3  mean : ', ztrasum
-         trn(:,:,:,jpno3) = trn(:,:,:,jpno3) * no3mean / ztrasum
-         
-         ! Silicate
-         ztrasum = 0.e0
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zvol = cvol(ji,jj,jk)
-#  if defined key_off_degrad
-                  zvol = zvol * facvol(ji,jj,jk)
-#  endif
-                  ztrasum = ztrasum + trn(ji,jj,jk,jpsil) * zvol
-               END DO
-            END DO
-         END DO
-         IF( lk_mpp )   CALL mpp_sum( ztrasum )     ! sum over the global domain  
-         ztrasum = ztrasum / areatot * 1.e6
-         IF(lwp) WRITE(numout,*) '       SiO3 mean : ', ztrasum
-         trn(:,:,:,jpsil) = MIN( 400.e-6,trn(:,:,:,jpsil) * siomean / ztrasum ) 
-         !
-      ENDIF
-
-!#if defined key_kriest
-!     !! Initialize number of particles from a standart restart file
-!     !! The name of big organic particles jpgoc has been only change
-!     !! and replace by jpnum but the values here are concentration
-!     trn(:,:,:,jppoc) = trn(:,:,:,jppoc) + trn(:,:,:,jpnum)
-!     trn(:,:,:,jpnum) = trn(:,:,:,jppoc) / ( 6. * xkr_massp )
-!#endif
-
-   END SUBROUTINE pis_dmp_ini
-
-#endif
-      !!----------------------------------------------------------------------
 
    SUBROUTINE trc_rst_stat
Index: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/trcwri.F90
===================================================================
--- branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/trcwri.F90	(revision 1987)
+++ branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/trcwri.F90	(revision 2007)
@@ -1,9 +1,9 @@
 MODULE trcwri
-   !!======================================================================
+   !!===================================================================================
    !!                       *** MODULE trcwri ***
-   !!    TOP :   Output of passive tracers 
-   !!======================================================================
-   !!             1.0  !  
-   !!                  !  2009-05 (C. Ethe )
+   !!    TOP :   Output of passive tracers
+   !!====================================================================================
+   !! History :   1.0  !  2009-05 (C. Ethe)  Original code
+   !!                  !  2010-03 (C. Ethe, R. Seferian ) Add the tracer transport trends
    !!----------------------------------------------------------------------
 #if defined key_top &&  defined key_iomput
@@ -11,7 +11,12 @@
    !!   'key_top' && 'key_iomput'                              TOP models
    !!----------------------------------------------------------------------
-   !! trc_wri     :  outputs of concentration fields
+   !! trc_wri_trc   :  outputs of concentration fields
+   !! trc_wri_trd   :  outputs of transport trends
    !!----------------------------------------------------------------------
+   USE dom_oce         ! ocean space and time domain variables
+   USE oce_trc
+   USE trp_trc
    USE trc
+   USE trdmld_trc_oce, ONLY : luttrd
    USE iom
 #if defined key_off_tra
@@ -35,7 +40,24 @@
 CONTAINS
 
-   SUBROUTINE trc_wri( kt )  
+   SUBROUTINE trc_wri( kt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE trc_wri  ***
+      !! 
+      !! ** Purpose :   output passive tracers fields and dynamical trends
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT( in ) :: kt
+      !!---------------------------------------------------------------------
+
+      !
+      CALL iom_setkt  ( kt + ndttrc - 1 )       ! set the passive tracer time step
+      CALL trc_wri_trc( kt              )       ! outputs for tracer concentration
+      CALL trc_wri_trd( kt              )       ! outputs for dynamical trends
+      CALL iom_setkt  ( kt              )       ! set the model time step
+      !
+   END SUBROUTINE trc_wri
+
+   SUBROUTINE trc_wri_trc( kt )  
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE trc_wri_trc  ***
       !!
       !! ** Purpose :   output passive tracers fields 
@@ -43,17 +65,11 @@
       INTEGER, INTENT( in ) :: kt       ! ocean time-step
       INTEGER               :: jn
-      CHARACTER (len=20)    :: cltra
+      CHARACTER (len=20)    :: cltra, cltras
 #if defined key_off_tra
       CHARACTER (len=40) :: clhstnam
       INTEGER ::   inum = 11            ! temporary logical unit
 #endif
-
       !!---------------------------------------------------------------------
  
-      ! Initialisation
-      ! --------------
-
-      CALL iom_setkt( kt + ndttrc - 1 ) ! set the passive tracer time step
-
 #if defined key_off_tra
       IF( kt == nittrc000 ) THEN
@@ -67,6 +83,4 @@
       ENDIF
 #endif
-
-
       ! write the tracer concentrations in the file
       ! ---------------------------------------
@@ -76,9 +90,53 @@
       END DO
       !
-      CALL iom_setkt( kt )       ! set the model time step
+   END SUBROUTINE trc_wri_trc
 
+# if defined key_trc_diatrd
+
+   SUBROUTINE trc_wri_trd( kt )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE trc_wri_trd  ***
+      !!
+      !! ** Purpose :   output of passive tracer : advection-diffusion trends
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt          ! ocean time-step
+      !!
+      CHARACTER (len=3) ::   cltra
+      INTEGER  ::   jn, jl, ikn
+      !!----------------------------------------------------------------------
+
+      DO jn = 1, jptra
+         IF( luttrd(jn) ) THEN
+            ikn = ikeep(jn)
+            DO jl = 1, jpdiatrc
+               IF( jl == jptrc_xad ) WRITE (cltra,"(3a)") 'XAD' ! x advection for tracer
+               IF( jl == jptrc_yad ) WRITE (cltra,"(3a)") 'YAD'  ! y advection for tracer
+               IF( jl == jptrc_zad ) WRITE (cltra,"(3a)") 'ZAD'  ! z advection for tracer
+               IF( jl == jptrc_xdf ) WRITE (cltra,"(3a)") 'XDF'  ! x diffusion for tracer
+               IF( jl == jptrc_ydf ) WRITE (cltra,"(3a)") 'YDF'  ! y diffusion for tracer
+               IF( jl == jptrc_zdf ) WRITE (cltra,"(3a)") 'ZDF'  ! z diffusion for tracer
+# if defined key_trcldf_eiv
+               IF( jl == jptrc_xei ) WRITE (cltra,"(3a)") 'XGV'  ! x gent velocity for tracer
+               IF( jl == jptrc_yei ) WRITE (cltra,"(3a)") 'YGV'  ! y gent velocity for tracer
+               IF( jl == jptrc_zei ) WRITE (cltra,"(3a)") 'ZGV'  ! z gent velocity for tracer
+# endif
+# if defined key_trcdmp
+               IF( jl == jptrc_dmp ) WRITE (cltra,"(3a)") 'DMP'  ! damping
+# endif
+               IF( jl == jptrc_sbc ) WRITE (cltra,"(3a)") 'SBC'  ! surface boundary conditions
+               ! write the trends
+               CALL iom_put( cltra, trtrd(:,:,:,ikn,jl) )
+            END DO
+         END IF
+      END DO
       !
-   END SUBROUTINE trc_wri
+   END SUBROUTINE trc_wri_trd
 
+# else
+   SUBROUTINE trc_wri_trd( kt )                      ! Dummy routine
+      INTEGER, INTENT ( in ) ::   kt
+   END SUBROUTINE trc_wri_trd
+#endif
 #else
    !!----------------------------------------------------------------------
@@ -90,5 +148,4 @@
    INTEGER, INTENT(in) :: kt
    END SUBROUTINE trc_wri
-
 #endif
 
