Index: branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/CONFIG/ORCA2_LIM_CRS/EXP00/iodef.xml
===================================================================
--- branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/CONFIG/ORCA2_LIM_CRS/EXP00/iodef.xml	(revision 3779)
+++ branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/CONFIG/ORCA2_LIM_CRS/EXP00/iodef.xml	(revision 3790)
@@ -291,5 +291,4 @@
        	<field id="soce_crs"   description="salinity"                unit="psu"  axis_ref="deptht" />
 	<field id="ssh_crs"    description="sea surface height"      unit="m"                      />
-        <field id="ssh2_crs"    description="sea surface height"      unit="m"                      />
 	<field id="sst_crs"    description="sea surface temperature" unit="degC"                   />
 	<field id="sss_crs"    description="sea surface salinity"    unit="psu"                    />
@@ -312,7 +311,4 @@
       <group id="gcrs_W"  axis_ref="depthw" grid_ref="grid_W_crs">
 	<field id="woce_crs"   description="ocean vertical velocity"   unit="m/s" />
-        <field id="woce2_crs"   description="ocean vertical velocity"   unit="m/s" />
-	<field id="wocet_crs"  description="ocean vertical velocity times temperature" unit="degC.m/s" />
-	<field id="woces_crs"  description="ocean vertical velocity times salinity"    unit="psu.m/s" />
       </group>
 
@@ -348,9 +344,8 @@
              <field ref="toce_crs"    name="votemper" />
              <field ref="soce_crs"    name="vosaline" />
-             <field ref="ssh_crs"     name="sossheig" />
-             <field ref="ssh2_crs"    name="sossheig2" />
-             <field ref="hdiv_crs"    name="vohdiver" />
       	     <field ref="sst_crs"     name="sosstsst" />
 	     <field ref="sss_crs"     name="sosaline" />
+             <field ref="ssh_crs"     name="sossheig" />
+             <field ref="hdiv_crs"    name="vohdiver" />
           </file>
 <!-- 
@@ -433,8 +428,7 @@
            <field ref="toce_crs"    name="votemper" />
            <field ref="soce_crs"    name="vosaline" />
-           <field ref="ssh_crs"     name="sossheig" />
-           <field ref="ssh2_crs"     name="sossheig2" />
 	   <field ref="sst_crs"     name="sosstsst" />
 	   <field ref="sss_crs"     name="sosaline" />
+           <field ref="ssh_crs"     name="sossheig" />
            <field ref="hdiv_crs"    name="vohdiver" />
         </file>	
@@ -478,7 +472,4 @@
         <file id="5d_gcrs_W" name="auto" description="ocean U grid coarsened variables" >
            <field ref="woce_crs"    name="vovecrtz"  />
-           <field ref="woce2_crs"   name="vovecrtz2"  />
-           <field ref="wocet_crs"   name="vovewoct" />
-           <field ref="woces_crs"   name="vovewocs" />
         </file>	
 	
Index: branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/CONFIG/ORCA2_LIM_CRS/MY_SRC/oce_trc.F90
===================================================================
--- branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/CONFIG/ORCA2_LIM_CRS/MY_SRC/oce_trc.F90	(revision 3779)
+++ branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/CONFIG/ORCA2_LIM_CRS/MY_SRC/oce_trc.F90	(revision 3790)
@@ -12,5 +12,4 @@
    !!----------------------------------------------------------------------
 
-   IF ( ln_crs ) THEN
 
    !* Domain size *
@@ -544,5 +543,4 @@
 # endif
 
-   ENDIF
 
 #else
Index: branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crs.F90
===================================================================
--- branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crs.F90	(revision 3779)
+++ branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crs.F90	(revision 3790)
@@ -54,5 +54,5 @@
 CONTAINS
 
-   SUBROUTINE crsfun_mask( p_ptmask, cd_type, p_pmask2d, p_cmask, p_cmask2d )
+   SUBROUTINE crsfun_mask( p_ptmask, cd_type, p_cmask )
       !!----------------------------------------------------------------
       !!               *** SUBROUTINE crsfun_mask ***
@@ -77,9 +77,7 @@
       !!----------------------------------------------------------------
       ! Arguments
-      REAL, DIMENSION(jpi,jpj,jpk)                 , INTENT(in)  :: p_ptmask  ! Parent grid tmask
-      CHARACTER(len=1)                             , INTENT(in)  :: cd_type   ! grid type (T,U,V,F)
-      REAL, DIMENSION(:,:),   OPTIONAL, ALLOCATABLE, INTENT(in)  :: p_pmask2d ! 2D mask, such as rnfmsk       
-      REAL, DIMENSION(:,:,:), OPTIONAL, ALLOCATABLE, INTENT(out) :: p_cmask   ! Coarse grid [t,u,v,f]mask
-      REAL, DIMENSION(:,:),   OPTIONAL, ALLOCATABLE, INTENT(out) :: p_cmask2d ! Coarse grid [t,u,v,f]mask
+      REAL, DIMENSION(jpi,jpj,jpk)        , INTENT(in)  :: p_ptmask  ! Parent grid tmask
+      CHARACTER(len=1)                    , INTENT(in)  :: cd_type   ! grid type (T,U,V,F)
+      REAL, DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(out) :: p_cmask   ! Coarse grid [t,u,v,f]mask
 
 
@@ -89,18 +87,5 @@
 
       ! Initialize
-      IF ( PRESENT( p_cmask ) ) THEN
-         ALLOCATE( p_cmask(jpi_crs,jpj_crs,jpk) )
-         p_cmask(:,:,:) = 0.0
-         ijpk = jpk
-      ENDIF
-
-      IF ( PRESENT( p_cmask2d ) ) THEN
-         ALLOCATE( p_cmask2d(jpi_crs,jpj_crs) )
-         p_cmask2d(:,:) = 0.0
-         ijpk = 1
-      ENDIF
- 
-
-      DO jk = 1, ijpk
+      DO jk = 1, jpk
          DO ji = 2, jpi_crsm1
             ijie = mie_crs(ji)
@@ -111,6 +96,4 @@
                ijjs = mjs_crs(jj)                   
 
-               IF ( PRESENT( p_cmask ) ) THEN
-               
                   IF ( cd_type == 'T' ) THEN  
                      p_cmask(ji,jj,jk) = SUM( p_ptmask(ijis:ijie,ijjs:ijje,jk) ) 
@@ -128,24 +111,4 @@
                   ENDIF
 
-               ENDIF
-
-               IF ( PRESENT( p_cmask2d ) ) THEN
-                  IF ( cd_type == 'T' ) THEN  
-                     p_cmask2d(ji,jj) = SUM( p_pmask2d(ijis:ijie,ijjs:ijje) ) 
-                     IF ( p_cmask2d(ji,jj) > 0 )  p_cmask2d(ji,jj) = 1
-                  ELSEIF ( cd_type == 'V' ) THEN               
-                     p_cmask2d(ji,jj) = SUM( p_pmask2d(ijis:ijie,ijje) * p_pmask2d(ijis:ijie,ijje+1) )
-                     IF ( p_cmask2d(ji,jj) > 0 ) p_cmask2d(ji,jj) = 1
-                  ELSEIF ( cd_type == 'U' ) THEN
-                     p_cmask2d(ji,jj) = SUM( p_pmask2d(ijie,ijjs:ijje) * p_pmask2d(ijie+1,ijjs:ijje) ) 
-                     IF ( p_cmask2d(ji,jj) > 0 ) p_cmask2d(ji,jj) = 1
-                  ELSE   ! fmask
-                     p_cmask2d(ji,jj) =  p_pmask2d(ijie,ijje) + p_pmask2d(ijie+1,ijje)        &
-                         &            + p_pmask2d(ijie,ijje+1) + p_pmask2d(ijie+1,ijje+1) 
-                     IF ( p_cmask2d(ji,jj) > 0 ) p_cmask2d(ji,jj) = 1
-                  ENDIF
-
-               ENDIF
-
             ENDDO
          ENDDO
@@ -155,26 +118,12 @@
 ! Retroactively add back the halo cells j=1, j=jpj_crs and i=1, i=jpi_crs.
 
-      IF ( PRESENT( p_cmask ) ) THEN
-         IF( nperio /= 0 ) THEN
-            CALL crs_lbc_lnk( cd_type,1.0,pt3d1=p_cmask )
-         ELSE
-            p_cmask( 1  ,:   ,:) = p_cmask(     2,:   ,:)          ! all points
-            p_cmask(jpi_crs,:,:) = p_cmask(jpi_crsm1,:,:)
-            p_cmask(   :,1   ,:) = p_cmask(     :,2   ,:)
-            p_cmask(:,jpj_crs,:) = p_cmask(:,jpj_crsm1,:)
-         ENDIF
+      IF( nperio /= 0 ) THEN
+          CALL crs_lbc_lnk( cd_type,1.0,pt3d1=p_cmask )
+      ELSE
+         p_cmask( 1  ,:   ,:) = p_cmask(     2,:   ,:)          ! all points
+         p_cmask(jpi_crs,:,:) = p_cmask(jpi_crsm1,:,:)
+         p_cmask(   :,1   ,:) = p_cmask(     :,2   ,:)
+         p_cmask(:,jpj_crs,:) = p_cmask(:,jpj_crsm1,:)
       ENDIF
-
-      IF ( PRESENT( p_cmask2d ) ) THEN
-         IF( nperio /= 0 ) THEN
-            CALL crs_lbc_lnk( cd_type,1.0,pt2d=p_cmask2d )
-         ELSE
-            p_cmask2d( 1  ,:   ) = p_cmask2d(     2,:   )          ! all points
-            p_cmask2d(jpi_crs,:) = p_cmask2d(jpi_crsm1,:)
-            p_cmask2d(   :,1   ) = p_cmask2d(     :,2   )
-            p_cmask2d(:,jpj_crs) = p_cmask2d(:,jpj_crsm1)
-         ENDIF
-      ENDIF
-
 
 
@@ -207,6 +156,6 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: p_pglam  ! Parent grid longitude
       CHARACTER(len=1),             INTENT(in)  :: cd_type   ! grid type (T,U,V,F) 
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:), INTENT(out) :: p_cgphi  ! Coarse grid latitude
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:), INTENT(out) :: p_cglam  ! Coarse grid longitude
+      REAL(wp), DIMENSION(jpi_crs,jpj_crs), INTENT(out) :: p_cgphi  ! Coarse grid latitude
+      REAL(wp), DIMENSION(jpi_crs,jpj_crs), INTENT(out) :: p_cglam  ! Coarse grid longitude
 
       !! Local variables
@@ -217,7 +166,4 @@
 
       !! Initialize output fields
-      IF ( .NOT. ALLOCATED( p_cgphi ) ) THEN
-        ALLOCATE( p_cgphi(jpi_crs,jpj_crs), p_cglam(jpi_crs,jpj_crs) )
-      ENDIF
       p_cgphi(:,:) = 0.e0
       p_cglam(:,:) = 0.e0
@@ -343,8 +289,8 @@
       REAL(wp), DIMENSION(jpi,jpj),     INTENT(in) :: p_e2     ! Parent grid U,V scale factors (e2)
 
-      REAL(wp), DIMENSION(:,:),   INTENT(out), ALLOCATABLE, OPTIONAL :: p_cfield2d_1 ! Coarse grid box 2D quantity
-      REAL(wp), DIMENSION(:,:),   INTENT(out), ALLOCATABLE, OPTIONAL :: p_cfield2d_2 ! Coarse grid box 2D quantity
-      REAL(wp), DIMENSION(:,:,:), INTENT(out), ALLOCATABLE, OPTIONAL :: p_cfield3d_1 ! Coarse grid box 3D quantity 
-      REAL(wp), DIMENSION(:,:,:), INTENT(out), ALLOCATABLE, OPTIONAL :: p_cfield3d_2 ! Coarse grid box 3D quantity 
+      REAL(wp), DIMENSION(:,:),   INTENT(out), OPTIONAL :: p_cfield2d_1 ! Coarse grid box 2D quantity
+      REAL(wp), DIMENSION(:,:),   INTENT(out), OPTIONAL :: p_cfield2d_2 ! Coarse grid box 2D quantity
+      REAL(wp), DIMENSION(:,:,:), INTENT(out), OPTIONAL :: p_cfield3d_1 ! Coarse grid box 3D quantity 
+      REAL(wp), DIMENSION(:,:,:), INTENT(out), OPTIONAL :: p_cfield3d_2 ! Coarse grid box 3D quantity 
 
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL :: p_fse3     ! Parent grid vertical level thickness (fse3u, fse3v)
@@ -381,9 +327,4 @@
       ze1(:,:) = p_e1(:,:)
       ze2(:,:) = p_e2(:,:)
-
-      IF ( PRESENT(p_cfield2d_1) ) ALLOCATE( p_cfield2d_1(jpi_crs,jpj_crs) )
-      IF ( PRESENT(p_cfield2d_2) ) ALLOCATE( p_cfield2d_2(jpi_crs,jpj_crs) ) 
-      IF ( PRESENT(p_cfield3d_1) ) ALLOCATE( p_cfield3d_1(jpi_crs,jpj_crs,jpk) ) 
-      IF ( PRESENT(p_cfield3d_2) ) ALLOCATE( p_cfield3d_2(jpi_crs,jpj_crs,jpk) ) 
 
       IF ( PRESENT(p_cfield2d_1) ) p_cfield2d_1(:,:) = 0.0
@@ -1012,7 +953,7 @@
          ze1e2(:,:) = p_e1e2t(:,:)
          ze3(:,:,:) = 1.0
-         IF ( PRESENT(p_pfield3d_1) ) THEN
+         IF ( PRESENT(p_pfield3d_2) ) THEN
             ! Prep to do the area-weighted average of (3D) W
-            zpfield3d(:,:,:) = p_pfield3d_1(:,:,:)
+            zpfield3d(:,:,:) = p_pfield3d_2(:,:,:)
          ELSEIF ( PRESENT(p_pfield2d) ) THEN
             ! Prep to do the area-weighted average of some 2D quantity 
@@ -1125,5 +1066,4 @@
             ENDDO
 
-
             IF ( ijpk == 1 ) THEN
                IF ( PRESENT(p_cfield2d) ) p_cfield2d(:,:) = zcfield2d(:,:) * zcmask(:,:,jk)
@@ -1171,5 +1111,5 @@
    END SUBROUTINE crsfun_TW
 
-SUBROUTINE crs_e3_max( p_e3, cd_type, p_mask, p_e3_crs)
+   SUBROUTINE crs_e3_max( p_e3, cd_type, p_mask, p_e3_crs)
       !!----------------------------------------------------------------
       !!               *** SUBROUTINE crsfun_TW ***
@@ -1326,5 +1266,5 @@
 
 
-SUBROUTINE crs_surf(p_e1, p_e2 ,p_e3, cd_type, p_mask, surf_crs, surf_msk_crs)
+   SUBROUTINE crs_surf(p_e1, p_e2 ,p_e3, cd_type, p_mask, surf_crs, surf_msk_crs)
       !!----------------------------------------------------------------
       !!               *** SUBROUTINE crsfun_TW ***
@@ -1386,5 +1326,6 @@
       INTEGER ::  ji, jj, jk                   ! dummy loop indices
       INTEGER :: ijie,ijis,ijje,ijjs,ijpk,jii,jjj
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze1, ze2, ze3, zsurf_crs, zsurf_msk_crs, zpmask  
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze1, ze2, ze3  
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zsurf_crs, zsurf_msk_crs, zpmask  
       !!----------------------------------------------------------------  
       ! Initialize
Index: branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crs_dom.F90
===================================================================
--- branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crs_dom.F90	(revision 3779)
+++ branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crs_dom.F90	(revision 3790)
@@ -108,39 +108,11 @@
       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE      :: tsn_crs
       REAL(wp), DIMENSION(:,:,:)  , ALLOCATABLE      :: un_crs, vn_crs, wn_crs
-      REAL(wp), DIMENSION(:,:,:)  , ALLOCATABLE      :: ut_crs, vt_crs, wt_crs, us_crs, vs_crs, ws_crs 
-      REAL(wp), DIMENSION(:,:,:)  , ALLOCATABLE      :: rhd_crs, rhop_crs, hdivn_crs    
-      REAL(wp), DIMENSION(:,:)    , ALLOCATABLE      :: sshb_crs, sshn_crs, ssha_crs    
-      REAL(wp), DIMENSION(:,:)    , ALLOCATABLE      :: sshun_crs, sshvn_crs, sshfn_crs 
-      REAL(wp), DIMENSION(:,:)    , ALLOCATABLE      :: sshub_crs, sshvb_crs, sshua_crs, sshva_crs 
-      REAL(wp), DIMENSION(:,:)    , ALLOCATABLE      :: hu_crs, hv_crs 
-      REAL(wp), DIMENSION(:,:)    , ALLOCATABLE      :: hdivbt_crs
-      REAL(wp), DIMENSION(:,:)    , ALLOCATABLE      :: ssh_b_crs, ssh_a_crs, ssh_un_crs, ssh_vn_crs        ! instantaneous fields
-
+      REAL(wp), DIMENSION(:,:,:)  , ALLOCATABLE      :: hdivn_crs    
+      REAL(wp), DIMENSION(:,:)    , ALLOCATABLE      :: sshn_crs    
 
       ! 
       ! Surface fluxes to pass to TOP
-      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE      ::  utau_crs, vtau_crs, wndm_crs, qsr_crs
-      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE      ::  del_emp_crs, sum_emp_crs
-      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE      ::  emp_crs, emp_b_crs, emps_crs 
-      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE      ::  rnf_crs, fr_i_crs, h_rnf_crs
-    
-      !
-      ! Lateral diffusivity (tracers) to pass to TOP
-      REAL(wp)                                     ::  rldf_crs, rn_aht_0_crs, aht0_crs, ahtb0_crs 
-
-#if defined key_traldf_c3d
-      REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) ::   ahtt_crs, ahtu_crs, ahtv_crs, ahtw_crs   !: 3D coefficients at T-,U-,V-,W-points
-      REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) ::   aeiu_crs, aeiv_crs, aeiw_crs   
-#elif defined key_traldf_c2d
-      REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   ahtt_crs, ahtu_crs, ahtv_crs, ahtw_crs   !: 2D coefficients at T-,U-,V-,W-points
-      REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   aeiu_crs, aeiv_crs, aeiw_crs   
-#elif defined key_traldf_c1d
-      REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:)     ::   ahtt_crs, ahtu_crs, ahtv_crs, ahtw_crs   !: 1D coefficients at T-,U-,V-,W-points
-      REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:)     ::   aeiu_crs, aeiv_crs, aeiw_crs   
-#else
-      REAL(wp), PUBLIC                                ::   ahtt_crs, ahtu_crs, ahtv_crs, ahtw_crs   !: scalar coefficients at T-,U-,V-,W-points
-#endif
-
- 
+      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE        ::  wndm_crs, qsr_crs, emp_crs, emps_crs
+
       ! Vertical diffusion
       REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)  ::  avt_crs           !: vert. diffusivity coef. [m2/s] at w-point for temp  
@@ -153,5 +125,4 @@
 
       ! Direction of lateral diffusion
-
 
 
@@ -164,5 +135,5 @@
       !!-------------------------------------------------------------------
       !! Local variables
-      INTEGER, DIMENSION(14) :: ierr
+      INTEGER, DIMENSION(15) :: ierr
 
       ierr(:) = 0
@@ -219,21 +190,10 @@
 
 
-      ALLOCATE( un_crs(jpi_crs,jpj_crs,jpk) ,   vn_crs(jpi_crs,jpj_crs,jpk) , &
-         &      wn_crs(jpi_crs,jpj_crs,jpk) ,                                 &
-         &      ut_crs(jpi_crs,jpj_crs,jpk) ,   vt_crs(jpi_crs,jpj_crs,jpk) , &
-         &      us_crs(jpi_crs,jpj_crs,jpk) ,   vs_crs(jpi_crs,jpj_crs,jpk) , &
-         &      wt_crs(jpi_crs,jpj_crs,jpk) ,   ws_crs(jpi_crs,jpj_crs,jpk) , &
-         &      rhd_crs(jpi_crs,jpj_crs,jpk) , rhop_crs(jpi_crs,jpj_crs,jpk) , &
-         &      hdivn_crs(jpi_crs,jpj_crs,jpk), &
-         &      STAT=ierr(11))
-
-      ALLOCATE( sshb_crs(jpi_crs,jpj_crs)   , sshn_crs(jpi_crs,jpj_crs)   , &
-         &      sshun_crs(jpi_crs,jpj_crs)  , sshvn_crs(jpi_crs,jpj_crs)  , &
-         &      sshfn_crs(jpi_crs,jpj_crs)  , emp_crs(jpi_crs,jpj_crs)    , &
-         &      del_emp_crs(jpi_crs,jpj_crs), sum_emp_crs(jpi_crs,jpj_crs), &
-         &      emp_b_crs(jpi_crs,jpj_crs)  , emps_crs(jpi_crs,jpj_crs)   , &
-         &      ssh_b_crs(jpi_crs,jpj_crs)  , ssh_a_crs(jpi_crs,jpj_crs)  , &
-         &      ssh_un_crs(jpi_crs,jpj_crs) , ssh_vn_crs(jpi_crs,jpj_crs)  , &
-         &      STAT=ierr(12)  )
+      ALLOCATE( un_crs(jpi_crs,jpj_crs,jpk) , vn_crs(jpi_crs,jpj_crs,jpk) , &
+         &      wn_crs(jpi_crs,jpj_crs,jpk) , hdivn_crs(jpi_crs,jpj_crs,jpk), STAT=ierr(11))
+
+      ALLOCATE( sshn_crs(jpi_crs,jpj_crs),  emp_crs(jpi_crs,jpj_crs)    , &
+         &      qsr_crs(jpi_crs,jpj_crs) ,  wndm_crs(jpi_crs,jpj_crs)    , &
+         &      emps_crs(jpi_crs,jpj_crs),        STAT=ierr(12)  )
  
       ALLOCATE( tsn_crs(jpi_crs,jpj_crs,jpk,jpts), avt_crs(jpi_crs,jpj_crs,jpk),    &
@@ -245,4 +205,5 @@
       ALLOCATE( nmln_crs(jpi_crs,jpj_crs) , hmld_crs(jpi_crs,jpj_crs) , &
          &      hmlp_crs(jpi_crs,jpj_crs) , hmlpt_crs(jpi_crs,jpj_crs) , STAT=ierr(14) )
+
 
       crs_dom_alloc = MAXVAL(ierr)
Index: branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsdiawri.F90
===================================================================
--- branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsdiawri.F90	(revision 3779)
+++ branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsdiawri.F90	(revision 3790)
@@ -92,136 +92,22 @@
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
       !!
-      INTEGER                      ::   ji, jj, jk              ! dummy loop indices
-      !!
-      REAL(wp), POINTER, DIMENSION(:,:)   :: z2d    ! 2D workspace
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d    ! 3D workspace
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d1   ! 3D workspace
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zfse3t ! 3D workspace for e3
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zfse3u ! 3D workspace for e3
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zfse3v ! 3D workspace for e3
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zfse3w ! 3D workspace for e3
-      REAL(wp), POINTER, DIMENSION(:,:)   :: ze1e2u ! 2D workspace
-      REAL(wp), POINTER, DIMENSION(:,:)   :: ze1e2v ! 2D workspace
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: z3dcrs, zw ! 3D workspace for coarse grid !cc
-      REAL(wp), POINTER, DIMENSION(:,:)   :: z2dcrs, ssh_crs2 ! 2D workspace for coarse grid
-      INTEGER        :: ialloc, iiki, iikn
-      INTEGER, SAVE  :: itsct
-      REAL(wp)       :: zdtj, zrtsct
-!!cc1
+      INTEGER               ::   ji, jj, jk              ! dummy loop indices
+      !!
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zfse3t, zfse3u, zfse3v, zfse3w ! 3D workspace for e3
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zt, zs 
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zt_crs, zs_crs !
       REAL(wp)       :: z2dcrsu, z2dcrsv
-!!cc1
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), SAVE :: ztsnm
-
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: zum, zvm, zwm  ! 3d work arrays for
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: zutm, zusm, zvtm, zvsm, zwtm, zwsm              ! accumulating sums for averaging
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: zvolw, zsuru, zsurv, zsurw, zhdivnm
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: zrn_absm, zrhopm, zrhdm, zrotbm
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: zavtm, zavsm, zuslpm,zvslpm, zwslpim, zwslpjm
-
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ), SAVE :: zsshbm, zsshm, zssham, zsshfm                    ! 2d work arrays 
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ), SAVE :: zsshubm, zsshum, zsshuam, zsurtaux, zutaum       ! for accumulating sums for averaging
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ), SAVE :: zsshvbm, zsshvm, zsshvam, zsurtauy, zvtaum       
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ), SAVE :: zhdivbtm, zempm, zemp_bm, zempsm
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ), SAVE :: zwndmm, zqsrm, zrnfm, zfr_im, zh_rnfm
-
-      !!
-      REAL(wp)                                      :: zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0, zraur
-      REAL(wp)                                      :: zij, zip1j, zijp1
-!!cc1      REAL(wp), ALLOCATABLE, DIMENSION(:,:)         :: z2dcrsu, z2dcrsv, z2dcrsf, zhdivbt
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)         :: z2dcrsf, zhdivbt
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)         :: zsshub, zsshua, zsshvb, zsshva   ! temp work arrays for instantaneous fields
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)       :: zhdiv, zvolt_wgt, zrhd, zrhop, zavt
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)         :: zee_t, zee_f, zee_u, zee_v
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)       :: zmut, zmuf
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)         :: zcrs_surtaux_wgt , zcrs_surtauy_wgt
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:)     :: ztsn      
-      REAL(wp) ::   zrho_c = 0.01_wp    ! density criterion for mixed layer depth
-      REAL(wp) ::   zavt_c = 5.e-4_wp   ! Kz criterion for the turbocline depth
-
-      INTEGER , ALLOCATABLE, DIMENSION(:,:  )       :: imld
-
-      CHARACTER(LEN=80) :: clmode
-      !!----------------------------------------------------------------------
+      !!
+       !!----------------------------------------------------------------------
       ! 
 
       IF( nn_timing == 1 )   CALL timing_start('crs_dia_wri')
-      ! 
-      ! 1. Initialize arrays
-      ! --------------------------------------------------------
-      CALL wrk_alloc( jpi , jpj , z2d, ze1e2u, ze1e2v )
-      CALL wrk_alloc( jpi , jpj , jpk , z3d, z3d1 )
-      CALL wrk_alloc( jpi , jpj , jpk , zfse3u, zfse3v )
-      CALL wrk_alloc( jpi , jpj , jpk , zfse3t, zfse3w )
-      CALL wrk_alloc( jpi_crs , jpj_crs , jpk , z3dcrs, zw ) !cc
-      CALL wrk_alloc( jpi_crs , jpj_crs , z2dcrs, ssh_crs2 )
-      
-
-      IF ( .NOT. ALLOCATED(ztsnm) ) THEN
-         ALLOCATE( ztsnm(jpi_crs,jpj_crs,jpk,jpts), STAT=ialloc )
-         IF ( ialloc /= 0 ) CALL ctl_warn('crsdiawri : failed to allocate cumulative 4d array')
-      ENDIF
-
-      IF ( .NOT. ALLOCATED(zum) ) THEN
-         ALLOCATE( zum(jpi_crs,jpj_crs,jpk) ,   zvm(jpi_crs,jpj_crs,jpk) ,  &
-            &    zavtm(jpi_crs,jpj_crs,jpk) , zvolw(jpi_crs,jpj_crs,jpk) ,  &
-            &    zsuru(jpi_crs,jpj_crs,jpk) , zsurv(jpi_crs,jpj_crs,jpk) ,  &
-            &    zsurw(jpi_crs,jpj_crs,jpk) ,   zwm(jpi_crs,jpj_crs,jpk) ,  &
-            &     zutm(jpi_crs,jpj_crs,jpk) ,  zusm(jpi_crs,jpj_crs,jpk) ,  &
-            &     zvtm(jpi_crs,jpj_crs,jpk) ,  zvsm(jpi_crs,jpj_crs,jpk) ,  &
-            &     zwtm(jpi_crs,jpj_crs,jpk) ,  zwsm(jpi_crs,jpj_crs,jpk) ,  &
-            &     zrhopm(jpi_crs,jpj_crs,jpk), zrhdm(jpi_crs,jpj_crs,jpk) ,  &
-            &   zhdivnm(jpi_crs,jpj_crs,jpk) ,  STAT=ialloc )
-         IF( ialloc /= 0 )   CALL ctl_warn('crsdiawri : failed to allocate cumulative 3d arrays')
-      ENDIF
-
-      IF ( .NOT. ALLOCATED(zsshm) ) THEN
-         ALLOCATE( zsshm(jpi_crs,jpj_crs) ,   zsshfm(jpi_crs,jpj_crs) , &
-            &     zsshbm(jpi_crs,jpj_crs) ,   zssham(jpi_crs,jpj_crs) , &
-            &     zsshum(jpi_crs,jpj_crs) ,   zsshvm(jpi_crs,jpj_crs) , &
-            &   zhdivbtm(jpi_crs,jpj_crs) ,                             &
-            &      zempm(jpi_crs,jpj_crs) ,  zemp_bm(jpi_crs,jpj_crs) , &
-            &     zempsm(jpi_crs,jpj_crs) ,                             &
-            &     zvtaum(jpi_crs,jpj_crs) ,   zutaum(jpi_crs,jpj_crs) , &
-            &   zsurtaux(jpi_crs,jpj_crs) , zsurtauy(jpi_crs,jpj_crs) , &
-            &     zwndmm(jpi_crs,jpj_crs) ,    zqsrm(jpi_crs,jpj_crs) , &
-            &      zrnfm(jpi_crs,jpj_crs) ,   zfr_im(jpi_crs,jpj_crs) , &
-            &    zh_rnfm(jpi_crs,jpj_crs) ,                             &
-            &     STAT=ialloc )
-         IF( ialloc /= 0 )   CALL ctl_warn('crsdiawri : failed to allocate cumulative 2d arrays')
-      ENDIF
-
-      IF ( .NOT. ALLOCATED(z2dcrsf) ) THEN
-!!cc1      IF ( .NOT. ALLOCATED(z2dcrsu) ) THEN
-!!cc1         ALLOCATE( z2dcrsu(jpi_crs,jpj_crs) , z2dcrsv(jpi_crs,jpj_crs) , &
-         ALLOCATE(   &
-            &      z2dcrsf(jpi_crs,jpj_crs) , zhdivbt(jpi_crs,jpj_crs) , &
-            &      zsshub(jpi_crs,jpj_crs)  , zsshvb(jpi_crs,jpj_crs)  , &
-            &      zsshua(jpi_crs,jpj_crs)  , zsshva(jpi_crs,jpj_crs)  , &
-            &      zee_t(jpi_crs,jpj_crs)   , zee_u(jpi_crs,jpj_crs)    , &
-            &      zee_v(jpi_crs,jpj_crs)   , zee_f(jpi_crs,jpj_crs)    , &
-            &      zmut(jpi_crs,jpj_crs,jpk), zmuf(jpi_crs,jpj_crs,jpk)    , &
-            &      zcrs_surtaux_wgt(jpi_crs,jpj_crs) , zcrs_surtauy_wgt(jpi_crs,jpj_crs) , &
-            &        imld(jpi_crs,jpj_crs)  ,                            &
-            &      STAT=ialloc )
-         IF( ialloc /= 0 )   CALL ctl_warn('crsdiawri : failed to allocate temp 2d arrays')
-      ENDIF
- 
-      IF ( .NOT. ALLOCATED(zhdiv) ) THEN
-         ALLOCATE( zhdiv(jpi_crs,jpj_crs,jpk) , zvolt_wgt(jpi_crs,jpj_crs,jpk) , &
-            &       zrhd(jpi_crs,jpj_crs,jpk) ,     zrhop(jpi_crs,jpj_crs,jpk) , &
-            &       zavt(jpi_crs,jpj_crs,jpk) , STAT=ialloc )
-         IF( ialloc /= 0 )   CALL ctl_warn('crsdiawri : failed to allocate temp 3d arrays')
-      ENDIF
-
-      IF ( .NOT. ALLOCATED(ztsn) ) THEN
-         ALLOCATE( ztsn(jpi_crs,jpj_crs,jpk,jpts),    STAT=ialloc )
-         IF( ialloc /= 0 )   CALL ctl_warn('crsdiawri : failed to allocate temp 4d arrays')
-      ENDIF
-
-      zw(:,:,:)=0.0
-      ! generic work arrays
-      z2d(:,:)      = 0.0
-      z3d(:,:,:)    = 0.0
-      z3d1(:,:,:)   = 0.0
+      
+      !  Initialize arrays
+      CALL wrk_alloc( jpi, jpj, jpk, zfse3t, zfse3w )
+      CALL wrk_alloc( jpi, jpj, jpk, zfse3u, zfse3v )
+      CALL wrk_alloc( jpi, jpj, jpk, zt, zs       )
+      !
+      CALL wrk_alloc( jpi_crs, jpj_crs, jpk, zt_crs, zs_crs )
 
       ! Depth work arrrays
@@ -231,727 +117,118 @@
       zfse3w(:,:,:) = fse3w(:,:,:)
 
-      ! Horizontal divergence work array
-      zhdiv(:,:,:) = 0._wp
-      zhdivbt(:,:) = 0._wp
-
-      zvolt_wgt(:,:,:) = 1._wp
-
-      ! Sea-surface height work arrays
-      zee_t(:,:) = 0._wp    ;   zee_f(:,:) = 0._wp
-      zee_u(:,:) = 0._wp    ;   zee_v(:,:) = 0._wp
-      zmut(:,:,:) = 0._wp   ;   zmuf(:,:,:) = 0._wp
-
-      ! work arrays
-      ze1e2u(:,:) = e1u(:,:) * e2u(:,:)
-      ze1e2v(:,:) = e1v(:,:) * e2v(:,:)
- 
       IF( kt == nit000  ) THEN
-         ! reset work arrays for average computation
-         itsct = 0
-         !
-         zum(:,:,:)   = 0._wp
-         zvm(:,:,:)   = 0._wp
-         ztsnm(:,:,:,:) = 0._wp
-         zwm(:,:,:)   = 0._wp
-         zutm(:,:,:)  = 0._wp
-         zusm(:,:,:)  = 0._wp
-         zvtm(:,:,:)  = 0._wp
-         zvsm(:,:,:)  = 0._wp
-         zwtm(:,:,:)  = 0._wp
-         zwsm(:,:,:)  = 0._wp
-         zavtm(:,:,:)  = 0._wp
-
-         zvolw(:,:,:) = 0._wp   ! t-centered grid box volume, masked
-         zsuru(:,:,:) = 0._wp   ! east face area, masked
-         zsurv(:,:,:) = 0._wp   ! north face area, masked
-         zsurw(:,:,:) = 0._wp   ! top face area, masked on gridT
-         zsurtaux(:,:) = 0._wp  ! top face area, masked on gridU
-         zsurtauy(:,:) = 0._wp  ! top face area, masked on gridV
-
-         zsshm(:,:)   = 0._wp
-         zsshum(:,:)  = 0._wp
-         zsshvm(:,:)  = 0._wp
-         zsshfm(:,:)  = 0._wp
-         zsshbm(:,:) = 0._wp
-         zsshubm(:,:) = 0._wp
-         zsshvbm(:,:) = 0._wp
-         zssham(:,:) = 0._wp
-         zsshuam(:,:) = 0._wp
-         zsshvam(:,:) = 0._wp
-         zwndmm(:,:) = 0._wp
-         zutaum(:,:) = 0._wp
-         zvtaum(:,:) = 0._wp
-         zqsrm(:,:) = 0._wp
-         zrnfm(:,:) = 0._wp
-         zh_rnfm(:,:) = 0._wp
-         zfr_im(:,:) = 0._wp
-
-         tsn_crs(:,:,:,:) = 0._wp               ! passive tracer array, now 
-         un_crs(:,:,:) = 0._wp                  ! u-velocity
-         vn_crs(:,:,:) = 0._wp                  ! v-velocity
-         wn_crs(:,:,:) = 0._wp                  ! bottom boundary condition: w(:,:,jpk)=0 (set once for all) 
-         sshn_crs(:,:) = 0._wp                  ! sea-surface height now
-         ut_crs(:,:,:) = 0._wp                  ! U*T
-         us_crs(:,:,:) = 0._wp                  ! U*S
-         vt_crs(:,:,:) = 0._wp                  ! V*T
-         vs_crs(:,:,:) = 0._wp                  ! V*S
-         avt_crs(:,:,:) = 0._wp                  ! avt
-
-         zcrs_surtaux_wgt(:,:) = ze1e2u(:,:) * umask(:,:,1)
-         zcrs_surtauy_wgt(:,:) = ze1e2v(:,:) * vmask(:,:,1)
-
+         tsn_crs(:,:,:,:) = 0._wp    ! temp/sal  array, now 
+         un_crs (:,:,:  ) = 0._wp    ! u-velocity
+         vn_crs (:,:,:  ) = 0._wp    ! v-velocity
+         wn_crs (:,:,:  ) = 0._wp    ! w
+         avt_crs(:,:,:  ) = 0._wp    ! avt
+         sshn_crs(:,:   ) = 0._wp    ! ssh
       ENDIF
 
-      ! 2. Coarsen fields and accumulate sums at each time step
+
+      ! 2. Coarsen fields at each time step
       ! --------------------------------------------------------
-! jes, TODO. 31 aug 2012
-! need to pass the 
-! 2D. utau, vtau, wndm, qsr, rnf, fi_i(:,:), h_rnf , nmln, hmld, hmlp, hmlpt
-! 3D. rn_abs(:,:,:), rhop, rhd=(rho-rau0)/rau0,  rotb, avt, avs, slopes, 
-! DIM?. aht0, ahtu, ahtv, ahtw, ahtt, aeiu, aeiv, aeiw
-
-
-      itsct = itsct+1     ! counter for the number of time steps 
-
-      ! 2.0  Weights
-      zvolw(:,:,:) = zvolw(:,:,:) + crs_volt_wgt(:,:,:)   ! masked ocean volume
-      zsurv(:,:,:) = zsurv(:,:,:) + crs_surfv_wgt(:,:,:)  ! masked north face surface area
-      zsuru(:,:,:) = zsuru(:,:,:) + crs_surfu_wgt(:,:,:)  ! masked east face surface area
-      zsurw(:,:,:) = zsurw(:,:,:) + crs_surfw_wgt(:,:,:)  ! masked top face surface area of gridT
-      zsurtaux(:,:) = zsurtaux(:,:) + zcrs_surtaux_wgt(:,:) ! masked top face surface area of gridU
-      zsurtauy(:,:) = zsurtauy(:,:) + zcrs_surtauy_wgt(:,:) ! masked top face surface area of gridV
-
-      ! 2.1 Temperature
-      !    2.1.1. now-temperature
-      z3d(:,:,:) = tsn(:,:,:,jp_tem)
+
+      !  Temperature
+      zt(:,:,:) = tsn(:,:,:,jp_tem)  ;      zt_crs(:,:,:) = 0._wp
       CALL crsfun( p_e1e2t=e1e2t, cd_type='T', cd_op='VOL', p_cmask=tmask_crs, p_ptmask=tmask, &
-         &         p_pfield3d_1=zfse3t, p_pfield3d_2=z3d, p_cfield3d=z3dcrs ) 
-      ztsn(:,:,:,1) = z3dcrs(:,:,:)
-      ztsnm(:,:,:,1) = ztsnm(:,:,:,1) + ztsn(:,:,:,1) 
-
-      ! 2.2 Salinity
-      !    2.2.1. now-salinity
-      z3d(:,:,:) = tsn(:,:,:,jp_sal)   ; z3dcrs(:,:,:) = 0.0 
+         &         p_pfield3d_1=zfse3t, p_pfield3d_2=zt, p_cfield3d=zt_crs )
+      tsn_crs(:,:,:,jp_tem) = zt_crs(:,:,:)
+      CALL crs_iom_put( "toce_crs", pv_r3d=tsn_crs(:,:,:,jp_tem) )    ! temp
+      CALL crs_iom_put( "sst_crs" , pv_r2d=tsn_crs(:,:,1,jp_tem) )    ! sst
+
+      !  Salinity
+      zt(:,:,:) = tsn(:,:,:,jp_sal)  ;      zt_crs(:,:,:) = 0._wp
       CALL crsfun( p_e1e2t=e1e2t, cd_type='T', cd_op='VOL', p_cmask=tmask_crs, p_ptmask=tmask, &
-         &         p_pfield3d_1=zfse3t, p_pfield3d_2=z3d, p_cfield3d=z3dcrs )
-      ztsn(:,:,:,2) = z3dcrs(:,:,:)
-      ztsnm(:,:,:,2) = ztsnm(:,:,:,2) + ztsn(:,:,:,2) 
-
-      !    2.5.9 Vertical eddy diffusion coefficient (avt) on gridW
-      z3d(:,:,:) = avt(:,:,:)   ;  z3dcrs(:,:,:) = 0.0 
-      CALL crsfun( p_e1e2t=e1e2t, cd_type='T', cd_op='SUM', p_cmask=tmask_crs, p_ptmask=tmask, &
-         &         p_pfield3d_1=zfse3t, p_pfield3d_2=z3d, p_cfield3d=z3dcrs )
-      zavtm(:,:,:) = zavtm(:,:,:) + z3dcrs(:,:,:)       
-
-      ! 2.3. U-Velocity
-      !    2.3.1. U
-      z3dcrs(:,:,:) = 0.0
+         &         p_pfield3d_1=zfse3t, p_pfield3d_2=zt, p_cfield3d=zt_crs )
+      tsn_crs(:,:,:,jp_sal) = zt_crs(:,:,:)
+      CALL crs_iom_put( "soce_crs" , pv_r3d=tsn_crs(:,:,:,jp_sal) )    ! sal
+      CALL crs_iom_put( "sss_crs"  , pv_r2d=tsn_crs(:,:,1,jp_sal) )    ! sss
+
+      !  U-velocity
       CALL crsfun( p_e1_e2=e2u, cd_type='U', psgn=-1.0, p_pmask=umask, &
-         &         p_fse3=zfse3u, p_pfield=un, p_surf_crs=crs_surfu_wgt, p_cfield3d=z3dcrs )
-      zum(:,:,:) = z3dcrs(:,:,:)
-      !    2.3.2.  UT, US
-      z3d(:,:,:) = 0.e0 ; z3d1(:,:,:) = 0.0
+         &         p_fse3=zfse3u, p_pfield=un, p_surf_crs=crs_surfu_wgt, p_cfield3d=un_crs )
+      CALL crs_iom_put( "uoce_crs"  , pv_r3d=un_crs  )   ! i-current 
+      !
+      zt(:,:,:) = 0._wp     ;    zs(:,:,:) = 0._wp  ;   zt_crs(:,:,:) = 0._wp   ;    zs_crs(:,:,:) = 0._wp
       DO jk = 1, jpkm1
          DO jj = 2, jpjm1
             DO ji = 2, jpim1   
-               z3d(ji,jj,jk)  = un(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) ) 
-               z3d1(ji,jj,jk) = un(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) ) 
+               zt(ji,jj,jk)  = un(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) ) 
+               zs(ji,jj,jk)  = un(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) ) 
             END DO
          END DO
       END DO
-      !   UT
-      z3dcrs(:,:,:) = 0.0
-      CALL crsfun(  p_e1_e2=e2u, cd_type='U', psgn=-1.0, p_pmask=umask, &
-         &         p_fse3=zfse3u, p_pfield=z3d, p_cfield3d=z3dcrs )
-      zutm(:,:,:) = zutm(:,:,:) + z3dcrs(:,:,:)
-      !   US
-      z3dcrs(:,:,:) = 0.0
-      CALL crsfun(  p_e1_e2=e2u, cd_type='U', psgn=-1.0, p_pmask=umask, &
-         &         p_fse3=zfse3u, p_pfield=z3d1, p_cfield3d=z3dcrs )
-      zusm(:,:,:) = zusm(:,:,:) + z3dcrs(:,:,:)
-
-
-      ! 2.4. V-Velocity
-      !    2.4.1. V
-      z3dcrs(:,:,:) = 0.0
+      CALL crsfun( p_e1_e2=e2u, cd_type='U', psgn=-1.0, p_pmask=umask, &
+         &         p_fse3=zfse3u, p_pfield=zt,  p_cfield3d=zt_crs )
+      CALL crs_iom_put( "uocet_crs" , pv_r3d=zt_crs )   ! uT
+      CALL crsfun( p_e1_e2=e2u, cd_type='U', psgn=-1.0, p_pmask=umask, &
+         &         p_fse3=zfse3u, p_pfield=zs,  p_cfield3d=zs_crs )
+      CALL crs_iom_put( "uoces_crs" , pv_r3d=zs_crs )   ! uS
+
+
+      !  V-velocity
       CALL crsfun( p_e1_e2=e1v, cd_type='V', psgn=-1.0, p_pmask=vmask, &
-         &         p_fse3=zfse3v, p_pfield=vn, p_surf_crs=crs_surfv_wgt, p_cfield3d=z3dcrs )
-      zvm(:,:,:) = z3dcrs(:,:,:)
-      !    2.4.2  VT, VS
-      z3d(:,:,:) = 0.e0 ; z3d1(:,:,:) = 0.0
+         &         p_fse3=zfse3v, p_pfield=vn, p_surf_crs=crs_surfv_wgt, p_cfield3d=vn_crs )
+      CALL crs_iom_put( "voce_crs"  , pv_r3d=vn_crs  )   ! v-current 
+      !                                                                                 
+      zt(:,:,:) = 0._wp     ;    zs(:,:,:) = 0._wp  ;   zt_crs(:,:,:) = 0._wp   ;    zs_crs(:,:,:) = 0._wp
       DO jk = 1, jpkm1
          DO jj = 2, jpjm1
             DO ji = 2, jpim1   
-               z3d(ji,jj,jk)  = vn(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) ) 
-               z3d1(ji,jj,jk) = vn(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) ) 
+               zt(ji,jj,jk)  = vn(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) ) 
+               zs(ji,jj,jk)  = vn(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) ) 
             END DO
          END DO
       END DO
-      !   VT
-      z3dcrs(:,:,:) = 0.0
       CALL crsfun( p_e1_e2=e1v, cd_type='V', psgn=-1.0, p_pmask=vmask, &
-         &         p_fse3=zfse3v, p_pfield=z3d, p_cfield3d=z3dcrs )
-      zvtm(:,:,:) = zvtm(:,:,:) + z3dcrs(:,:,:)
-      !   VS
-      z3dcrs(:,:,:) = 0.0
+         &         p_fse3=zfse3v, p_pfield=zt, p_cfield3d=zt_crs )
+      CALL crs_iom_put( "vocet_crs" , pv_r3d=zt_crs )   ! vT
+
       CALL crsfun( p_e1_e2=e1v, cd_type='V', psgn=-1.0, p_pmask=vmask, &
-         &         p_fse3=zfse3v, p_pfield=z3d1, p_cfield3d=z3dcrs )
-      zvsm(:,:,:) = zvsm(:,:,:) + z3dcrs(:,:,:)
-      
-      
-      ! Vitesse vertical !cc
-      z3dcrs(:,:,:) = 0.0
-      CALL crsfun( p_e1e2t=e1e2t, cd_type='T', cd_op='ARE', p_cmask=tmask_crs, p_ptmask=tmask, &
-        &           p_pfield3d_1=wn, p_cfield3d=z3dcrs)
-      zw(:,:,:) = z3dcrs(:,:,:)
-
-
-      ! 2.5. Surface boundary conditions
-      !    2.5.1 Evaporation minus Precipitation (emp)
-      z2d(:,:) = emp_b(:,:)   ; z2dcrs(:,:) = 0._wp
-      CALL crsfun( p_e1e2t=e1e2t, cd_type='T', cd_op='ARE', p_cmask=tmask_crs, p_ptmask=tmask, &
-         &         p_pfield2d=z2d, p_cfield2d=z2dcrs )   
-      zemp_bm(:,:) = zemp_bm(:,:) + z2dcrs     ! summation of emp_b_crs
-      !
-      z2d(:,:) = emp(:,:)   ; z2dcrs(:,:) = 0._wp
-      CALL crsfun( p_e1e2t=e1e2t, cd_type='T', cd_op='ARE', p_cmask=tmask_crs, p_ptmask=tmask, &
-         &         p_pfield2d=z2d, p_cfield2d=z2dcrs )
-      zempm(:,:) = zempm(:,:) + z2dcrs     ! summation of emp_crs
-      !
-      z2d(:,:) = emps(:,:)   ; z2dcrs(:,:) = 0._wp
-      CALL crsfun( p_e1e2t=e1e2t, cd_type='T', cd_op='ARE', p_cmask=tmask_crs, p_ptmask=tmask, &
-         &         p_pfield2d=z2d, p_cfield2d=z2dcrs )
-      zempsm(:,:) = zempsm(:,:) + z2dcrs     ! summation of emps_crs
-      !         terms needed for ssh calculation
-      z2d(:,:) = emp_b(:,:)-emp(:,:)
-      CALL crsfun( p_e1e2t=e1e2t, cd_type='T', cd_op='ARE', p_cmask=tmask_crs, p_ptmask=tmask, &
-         &         p_pfield2d=z2d, p_cfield2d=del_emp_crs )
-      z2d(:,:) = emp_b(:,:)+emp(:,:)
-      CALL crsfun( p_e1e2t=e1e2t, cd_type='T', cd_op='ARE', p_cmask=tmask_crs, p_ptmask=tmask, &
-         &         p_pfield2d=z2d, p_cfield2d=sum_emp_crs )
-      !
-      !    2.5.2 Wind Stress (utau, vtau) 
-      !          Taux
-      z2d(:,:) = utau(:,:)   ; z2dcrs(:,:) = 0._wp
-      CALL crsfun( p_e1e2t=ze1e2u, cd_type='U', cd_op='SUM', p_cmask=umask_crs, p_ptmask=umask, &
-         &         p_pfield2d=z2d, p_cfield2d=z2dcrs )
-      zutaum(:,:) = zutaum(:,:) + z2dcrs(:,:)      
-      !          Tauy
-      z2d(:,:) = vtau(:,:)   ; z2dcrs(:,:) = 0._wp
-      CALL crsfun( p_e1e2t=ze1e2v, cd_type='V', cd_op='SUM', p_cmask=vmask_crs, p_ptmask=vmask, &
-         &         p_pfield2d=z2d, p_cfield2d=z2dcrs )
-      zvtaum(:,:) = zvtaum(:,:) + z2dcrs(:,:)
- 
-      !    2.5.3 Scalar Wind (wndm)
-      z2d(:,:) = wndm(:,:)   ; z2dcrs(:,:) = 0._wp
-      CALL crsfun( p_e1e2t=e1e2t, cd_type='T', cd_op='SUM', p_cmask=tmask_crs, p_ptmask=tmask, &
-         &         p_pfield2d=z2d, p_cfield2d=z2dcrs )
-      zwndmm(:,:) = zwndmm(:,:) + z2dcrs(:,:)
-
-      !    2.5.4 Surface shortwave radiative flux (qsr)
-      z2d(:,:) = qsr(:,:)   ; z2dcrs(:,:) = 0._wp
-      CALL crsfun( p_e1e2t=e1e2t, cd_type='T', cd_op='SUM', p_cmask=tmask_crs, p_ptmask=tmask, &
-         &         p_pfield2d=z2d, p_cfield2d=z2dcrs )
-      zqsrm(:,:) = zqsrm(:,:) + z2dcrs(:,:)
-
-      !    2.5.5 Runoff freshwater flux
-      z2d(:,:) = rnf(:,:)   ; z2dcrs(:,:) = 0._wp
-      CALL crsfun( p_e1e2t=e1e2t, cd_type='T', cd_op='SUM', p_cmask=tmask_crs, p_ptmask=tmask, &
-         &         p_pfield2d=z2d, p_cfield2d=z2dcrs )
-      zrnfm(:,:) = zrnfm(:,:) + z2dcrs(:,:)
-
-      !    2.5.6 Depth of runoff in meters h_rnf(:,:))
-      z2d(:,:) = h_rnf(:,:) ; z2dcrs(:,:) = 0._wp
-      CALL crsfun(p_e1e2t=e1e2t, cd_type='T', cd_op='SUM', p_cmask=tmask_crs, p_ptmask=tmask, &
-         &         p_pfield2d=z2d, p_cfield2d=z2dcrs )
-      zh_rnfm(:,:) = zh_rnfm(:,:) + z2dcrs(:,:)
-
-      !    2.5.7 Sea-ice fraction (fr_i)
-      z2d(:,:) = fr_i(:,:)   ; z2dcrs(:,:) = 0._wp
-      CALL crsfun( p_e1e2t=e1e2t, cd_type='T', cd_op='SUM', p_cmask=tmask_crs, p_ptmask=tmask, &
-         &         p_pfield2d=z2d, p_cfield2d=z2dcrs )
-      zfr_im(:,:) = zfr_im(:,:) + z2dcrs(:,:)
-
-
-
-      !    2.5.8 Potential density (rhop) 
-! jes need to declare rho cumulative arrays
-  !    CALL crs_eos( ztsn, zrhd, zrhop )
-
-      CALL dom_grid_crs   ! Save the parent grid information  & Switch to coarse grid domain
-
-      CALL eos( ztsn, zrhd, zrhop )
-
-      CALL dom_grid_glo   ! Return to parent grid domain
-
-      zrhdm(:,:,:) = zrhdm(:,:,:) + zrhd
-      zrhopm(:,:,:) = zrhopm(:,:,:) + zrhop
-
-      ! 2.6. Diagnostic variables: horizontal divergence and vertical velocity
-      !    2.6.1 ocean depth at u- v- points
-      ! hu_crs(:,:) = ssh_un_crs(:,:)                ! now ocean depth (at u- and v-points)
-      ! hv_crs(:,:) = ssh_vn_crs(:,:)
-      !                                            ! now masked inverse of the ocean depth (at u- and v-points)
-      ! jes. not hur_crs, huv_crs are needed. remove?
-      !hur_crs(:,:) = umask_crs(:,:,1) / ( hu_crs(:,:) + 1._wp - umask_crs(:,:,1) )
-      !hvr_crs(:,:) = vmask_crs(:,:,1) / ( hv_crs(:,:) + 1._wp - vmask_crs(:,:,1) )
-      !
-      !    2.6.2. Horizontal divergence
-      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
-      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
-      z1_2dt = 0.e0 / (z2dt)
-      !
-      WHERE ( crs_volt_wgt /= 0 ) zvolt_wgt(:,:,:) = 1.0/crs_volt_wgt(:,:,:)
-      !
-      fse3t_crs(:,:,:) = e3t_crs(:,:,:)
-      IF ( lk_vvl ) fse3t_crs(:,:,:) = fse3t_n_crs(:,:,:)
-
-      zhdivbt(:,:) = 0.0
+         &         p_fse3=zfse3v, p_pfield=zs, p_cfield3d=zs_crs )
+      CALL crs_iom_put( "voces_crs" , pv_r3d=zs_crs )   ! vS
+      !
+
+
+      !  Horizontal divergence ( following OPA_SRC/DYN/divcur.F90 ) 
       DO jk = 1, jpkm1
-!!cc1         z2dcrsu(:,:) = 0.0; z2dcrsv(:,:) = 0.0   
          DO ji = 2, jpi_crsm1
             DO jj = 2, jpj_crsm1
-
-               ! Horizontal divergence ( following OPA_SRC/DYN/divcur.F90 ) 
-               ! with partial steps and/or variable layer thicknesses for W
-!!cc1               z2dcrsu(ji,jj) =  ( zum(ji,  jj,jk) * crs_surfu_wgt(ji,  jj,jk)  ) &
-               z2dcrsu =  ( zum(ji,  jj,jk) * crs_surfu_wgt(ji,  jj,jk)  ) &
-                  &            - ( zum(ji-1,jj,jk) * crs_surfu_wgt(ji-1,jj,jk)  )
-!!cc1               z2dcrsv(ji,jj) =  ( zvm(ji,jj,  jk) * crs_surfv_wgt(ji,jj  ,jk)  ) &
-               z2dcrsv =  ( zvm(ji,jj,  jk) * crs_surfv_wgt(ji,jj  ,jk)  ) &
-                  &            - ( zvm(ji,jj-1,jk) * crs_surfv_wgt(ji,jj-1,jk)  )
-
-!!cc1               zhdiv(ji,jj,jk) = ( z2dcrsu(ji,jj) + z2dcrsv(ji,jj) ) * zvolt_wgt(ji,jj,jk) 
-               zhdiv(ji,jj,jk) = ( z2dcrsu + z2dcrsv ) * zvolt_wgt(ji,jj,jk) 
- 
+               IF( tmask_crs(ji,jj,jk ) > 0 ) THEN
+                   z2dcrsu =  ( un_crs(ji  ,jj  ,jk) * crs_surfu_wgt(ji  ,jj  ,jk) ) &
+                      &     - ( un_crs(ji-1,jj  ,jk) * crs_surfu_wgt(ji-1,jj  ,jk) )
+                   z2dcrsv =  ( vn_crs(ji  ,jj  ,jk) * crs_surfv_wgt(ji  ,jj  ,jk) ) &
+                      &     - ( vn_crs(ji  ,jj-1,jk) * crs_surfv_wgt(ji  ,jj-1,jk) )
+                   !
+                   hdivn_crs(ji,jj,jk) = ( z2dcrsu + z2dcrsv ) / crs_volt_wgt(ji,jj,jk) 
+               ENDIF
             ENDDO
          ENDDO
-
-         ! Horizontal barotropic divergence for ssh ( following OPA_SRC/DYN/sshwzv.F90 )
-
-         zhdivbt(:,:) =  zhdivbt(:,:) + fse3t_crs(:,:,jk) * zhdiv(:,:,jk) 
-
       ENDDO
-
-      CALL crs_lbc_lnk( 'T', 1.0, pt3d1=zhdiv )
-      CALL crs_lbc_lnk( 'T', 1.0, pt3d1=zhdivbt )
-      ! accumulate 
-      zhdivnm(:,:,:) = zhdivnm(:,:,:) + zhdiv
-      zhdivbtm(:,:) = zhdivbtm(:,:) + zhdivbt
-
-!      !!! Calcul de l'energie cinétique   !cc  !! DECLARE LES VARIABLES
-!      
-!      zun2(:,:,:) = un(:,:,:) * un(:,:,:)    ! U^2
-!      zvn2(:,:,:) = vn(:,:,:) * vn(:,:,:)    ! V^2
-!      
-!!      ! moyenne sur i de U^2
-!      
-!      DO ji = 1, jpiglo-1
-!         zun2(ji,:,:) = 0.5 * (zun2(ji,:,:) + zun2(ji+1,:,:))
-!      END
-!      ji = jpiglo
-!      zun2(ji,:,:) = 0.5 * zun2(ji,:,:) 
-!      uun2(:,:,:) = zun2(:,:,:)
-!      
-!      
-!         CALL crs_lbc_lnk( cd_type='T', 1.0, pt3d1=uun2)
-!         ! lbcnlk met la ligne jpj = 1 a 0 donc il faut la remettre en ne pas oubliant le cyclique est-ouest 
-!         ! a faire un case pour cyclique est-ouest ou non.      
-!         uun2(   :   ,1,:) = zun2(    :    ,1,:)   
-!         
-!         IF ((jperio==4) .OR. (jperio==6)) THEN
-!            uun2(   1   ,1,:) = zun2(jpi_crsm1,1,:)
-!            uun2(jpi_crs,1,:) = zun2(    2    ,1,:)
-!         ENDIF  
-!
-!
- !     DO jj = 1, jpjglo-1
- !        zvn2(:,jj,:) = 0.5 * (zvn2(:,jj,:) + zvn2(:,jj+1,:))
- !     END
- !     jj = jpjglo 
- !     zvn2(:,jj,:) = 0.5 * zvn2(:,jj,:)
- !     vvn2(:,:,:) = zvn2(:,:,:)
- !     
- !     CALL crs_lbc_lnk( cd_type='T', 1.0, pt3d1=vvn2)
- !        ! lbcnlk met la ligne jpj = 1 a 0 donc il faut la remettre en ne pas oubliant le cyclique est-ouest 
- !        ! a faire un case pour cyclique est-ouest ou non.      
- !        vvn2(   :   ,1,:) = zvn2(    :    ,1,:)   
- !        
- !        IF ((jperio==4) .OR. (jperio==6)) THEN
- !           vvn2(   1   ,1,:) = zvn2(jpi_crsm1,1,:)
- !           vvn2(jpi_crs,1,:) = zvn2(    2    ,1,:)
- !        ENDIF  
-
-
-
-      
-      !    2.6.3. Sea-surface Height  ( See DOM/istate.F90: ssh init; OPA_SRC/DYN/sshwzv.F90: ssh update )
-      !cc
-      z2dcrs(:,:) = 0.0
-      CALL crsfun( p_e1e2t=e1e2t, cd_type='T', cd_op='ARE', p_cmask=tmask_crs, p_ptmask=tmask, p_pfield2d=sshn, &
-      &                  p_cfield2d=z2dcrs )
-      ssh_crs2(:,:) = z2dcrs(:,:)
-  !    WRITE(numout,*) 'test', sshn(:,:)
-  !    WRITE(numout,*) 'test', ssh_crs2(:,:)
-      
-     !cc    
-      
-      !       set some temp variables
-      z2dcrs(:,:) = 0.0; z2dcrsf(:,:) = 0.0
-      zcoefu = 0.0; zcoefv = 0.0; zcoeff = 0.0; 
-      zraur = 1._wp / rau0
-      z1_rau0 = 0.5 / rau0
-      !       sshn, sshb to accumulate
-      DO jj = 2, jpj_crsm1                                      ! leap-frog on ssh_e
-         DO ji = 2, jpi_crsm1
-            z2dcrs(ji,jj) = ( ssh_b_crs(ji,jj) - z2dt * ( zraur * ( del_emp_crs(ji,jj) ) + zhdivbt(ji,jj) ) ) * tmask_crs(ji,jj,1) 
-         END DO
-      END DO
-      CALL crs_lbc_lnk('T',1.0,pt2d=z2dcrs)   
-      zsshm(:,:)  = zsshm(:,:)  + z2dcrs(:,:)
-      zsshbm(:,:) = zsshbm(:,:) + ssh_b_crs(:,:)
-
-      !
-      !       sshu, sshv, sshf, for vvl case. Calculate then accumulate the fields for time-mean
-      !
-            DO jj = 1, jpj_crsm1
-               DO ji = 1, jpi_crsm1                    
-                  zcoefu = 0.5  * umask_crs(ji,jj,1) / ( e1u_crs(ji,jj) * e2u_crs(ji,jj) )
-                  zcoefv = 0.5  * vmask_crs(ji,jj,1) / ( e1v_crs(ji,jj) * e2v_crs(ji,jj) )
-                  zij   =  e1e2t_crs(ji  ,jj  ) * z2dcrs(ji  ,jj  ) 
-                  zip1j =  e1e2t_crs(ji+1,jj  ) * z2dcrs(ji+1,jj  ) 
-                  zijp1 =  e1e2t_crs(ji  ,jj+1) * z2dcrs(ji  ,jj+1) 
-              
-                  ssh_un_crs(ji,jj) = zcoefu * ( zij + zip1j )
-                  ssh_vn_crs(ji,jj) = zcoefv * ( zij + zijp1 )
-         IF( lk_vvl ) THEN
-                  zsshub(ji,jj) = zcoefu * ( e1e2t_crs(ji,  jj) * ssh_b_crs(ji  ,jj)   &
-                     &                        +  e1e2t_crs(ji+1,jj) * ssh_b_crs(ji+1,jj) )
-                  zsshvb(ji,jj) = zcoefv * ( e1e2t_crs(ji,  jj) * ssh_b_crs(ji,  jj)   &
-                     &                        +  e1e2t_crs(ji,jj+1) * ssh_b_crs(ji,jj+1) )
-         ENDIF
-               END DO
-            END DO
-
-            DO jj = 1, jpj_crsm1
-               DO ji = 1, jpi_crsm1     
-                  z2dcrsf(ji,jj) = 0.5  * umask_crs(ji,jj,1) * umask_crs(ji,jj+1,1)                   &
-                       &               / ( e1f_crs(ji,jj  ) * e2f_crs(ji,jj  ) )                     &
-                       &               * ( e1u_crs(ji,jj  ) * e2u_crs(ji,jj  ) * ssh_un_crs(ji,jj  )     &
-                       &                 + e1u_crs(ji,jj+1) * e2u_crs(ji,jj+1) * ssh_un_crs(ji,jj+1) )
-               END DO
-            END DO
-      !     accumulate fields in time
-            CALL crs_lbc_lnk( 'U', 1., pt2d=ssh_un_crs )   ! sshun_crs
-            CALL crs_lbc_lnk( 'V', 1., pt2d=ssh_vn_crs )   ! sshvn_crs
-            CALL crs_lbc_lnk('F', 1., pt2d=z2dcrsf )
-            zsshum(:,:)  = zsshum(:,:)  + ssh_un_crs(:,:)
-            zsshvm(:,:)  = zsshvm(:,:)  + ssh_vn_crs(:,:)
-            zsshfm(:,:)  = zsshfm(:,:)  + z2dcrsf(:,:)
-         IF( lk_vvl ) THEN
-            CALL crs_lbc_lnk( 'U', 1., pt2d=zsshub )
-            CALL crs_lbc_lnk( 'V', 1., pt2d=zsshvb )
-            zsshubm(:,:) = zsshubm(:,:) + zsshub(:,:) 
-            zsshvbm(:,:) = zsshvbm(:,:) + zsshvb(:,:) 
-         ENDIF
-      !       update sshb to sshn
-      ssh_b_crs(:,:) = z2dcrs(:,:) 
-      !       update ssha from sshn
-      ssh_a_crs(:,:) = (  z2dcrs(:,:) - z2dt * ( z1_rau0 * ( sum_emp_crs(:,:) ) + zhdivbt(:,:) )  ) * tmask(:,:,1)
-      CALL crs_lbc_lnk('T', 1.0, pt2d=ssh_a_crs)
-      zssham(:,:) = zssham(:,:) + ssh_a_crs(:,:)
-      !       accumulate fields in time and update sshua from sshun
-      IF ( lk_vvl ) THEN
-         DO jj = 1, jpj_crsm1
-            DO ji = 1, jpi_crsm1     
-               zsshua(ji,jj) = 0.5  * umask_crs(ji,jj,1) / ( e1u_crs(ji  ,jj) * e2u_crs(ji  ,jj) )                   &
-                  &                                  * ( e1t_crs(ji  ,jj) * e2t_crs(ji  ,jj) * ssh_a_crs(ji  ,jj)     &
-                  &                                    + e1t_crs(ji+1,jj) * e2t_crs(ji+1,jj) * ssh_a_crs(ji+1,jj) )
-               zsshva(ji,jj) = 0.5  * vmask_crs(ji,jj,1) / ( e1v_crs(ji,jj  ) * e2v_crs(ji,jj  ) )                   &
-                  &                                  * ( e1t_crs(ji,jj  ) * e2t_crs(ji,jj  ) * ssh_a_crs(ji,jj  )     &
-                  &                                    + e1t_crs(ji,jj+1) * e2t_crs(ji,jj+1) * ssh_a_crs(ji,jj+1) )
-            END DO
-         END DO
-
-         CALL crs_lbc_lnk( 'U', 1., pt2d=zsshua )
-         CALL crs_lbc_lnk( 'V', 1., pt2d=zsshva )
-
-         zsshuam(:,:) = zsshuam(:,:) + zsshua(:,:) 
-         zsshvam(:,:) = zsshvam(:,:) + zsshva(:,:)     
-      ENDIF
-      !
-      !    2.5.2.  W-velocity ( See OPA_SRC/DYN/sshwzv.F90 )
-      IF ( lk_vvl ) THEN
-         !                                 !==  mu computation  ==!
-         zee_t(:,:) = e3t_crs(:,:,1)       ! Lower bound : thickness of the first model level
-         zee_f(:,:) = e3f_crs(:,:,1)      
-         zee_u(:,:) = e3u_crs(:,:,1)      
-         zee_v(:,:) = e3v_crs(:,:,1)      
-         !
-         DO jk = 2, jpkm1                  ! Sum of the masked vertical scale factors
-            zee_t(:,:) = zee_t(:,:) + e3t_crs(:,:,jk) * tmask_crs(:,:,jk)
-            zee_u(:,:) = zee_u(:,:) + e3u_crs(:,:,jk) * umask_crs(:,:,jk)
-            zee_v(:,:) = zee_v(:,:) + e3v_crs(:,:,jk) * vmask_crs(:,:,jk)
-            DO jj = 1, jpj_crsm1                      ! f-point : fmask=shlat at coasts, use the product of umask
-               zee_f(:,jj) = zee_f(:,jj) + e3f_crs(:,jj,jk) *  umask_crs(:,jj,jk) * umask_crs(:,jj+1,jk)
-            END DO
-         END DO
-         !                                 ! Compute and mask the inverse of the local depth at T, U, V and F points
-         zee_t(:,:) = 1._wp / zee_t(:,:) * tmask_crs(:,:,1)
-         zee_u(:,:) = 1._wp / zee_u(:,:) * umask_crs(:,:,1)
-         zee_v(:,:) = 1._wp / zee_v(:,:) * vmask_crs(:,:,1)
-         DO jj = 1, jpj_crsm1                               ! f-point case fmask cannot be used
-            zee_f(:,jj) = 1._wp / zee_f(:,jj) * umask_crs(:,jj,1) * umask_crs(:,jj+1,1)
-         END DO
-         CALL crs_lbc_lnk( 'F', 1., pt2d=zee_f )                 ! lateral boundary condition on ee_f
-         !
-         DO jk = 1, jpk                    ! mu coefficients
-            zmut(:,:,jk) = zee_t(:,:) * tmask_crs(:,:,jk)     ! T-point at T levels
-         END DO
-         DO jk = 1, jpk                                 ! F-point : fmask=shlat at coasts, use the product of umask
-            DO jj = 1, jpj_crsm1
-               zmuf(:,jj,jk) = zee_f(:,jj) * umask_crs(:,jj,jk) * umask_crs(:,jj+1,jk)   ! at T levels
-            END DO
-            zmuf(:,jpj,jk) = 0._wp
-         END DO
-         CALL crs_lbc_lnk( 'F', 1., pt3d1=zmuf )                   ! lateral boundary condition
-
-         DO jk = 1, jpk                    ! mu coefficients
-            fse3t_crs(:,:,jk)   = e3t_crs(:,:,jk)
-            fse3t_n_crs(:,:,jk) = e3t_crs(:,:,jk) + z2dcrs(:,:)*zmut(:,:,jk)
-            fse3t_b_crs(:,:,jk) = fse3t_b_crs(:,:,jk)
-            fse3t_a_crs(:,:,jk) = e3t_crs(:,:,jk) + ssh_a_crs(:,:)*zmut(:,:,jk)
-         ENDDO
-! jes. 30 aug 2012.  Need to add here fse3u, fse3v, fse3f - n,b,a
-      ELSE
-         fse3t_crs = e3t_crs
-!         fse3w_crs = e3w_crs
-!         fse3u_crs = e3u_crs
-!         fse3v_crs = e3v_crs
-      ENDIF
-      !
-      z3dcrs(:,:,:) = 0.0      
-      IF ( lk_vvl ) THEN
-         DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
-            z3dcrs(:,:,jk) = z3dcrs(:,:,jk+1) -   fse3t_n_crs(:,:,jk) * zhdiv(:,:,jk)        &
-               &                      - ( fse3t_a_crs(:,:,jk) - fse3t_b_crs(:,:,jk) )    &
-               &                         * tmask_crs(:,:,jk) * z1_2dt
-         END DO
-      ELSE
-         DO jk = jpkm1, 1, -1                             ! calculate from the bottom-up
-            z3dcrs(:,:,jk) = z3dcrs(:,:,jk+1) - fse3t_crs(:,:,jk) * zhdiv(:,:,jk)        
-         ENDDO
-      ENDIF
-      CALL crs_lbc_lnk( 'T', 1.0, pt3d1=z3dcrs )
-      zwm(:,:,:) = zwm(:,:,:) + z3dcrs(:,:,:)    
-
-      ! 2.5.3 Slopes  for slopes look at OPA_SRC/DYN/dynldf_iso.F90 and OPA_SRC/LDF/ldfslp.F90
-
-
-      ! 3. Output of some dynamics and tracer fields and selected fields
-      !   * Within the nn_fcrs frequence for iom 
-      ! --------------------------------------------------------
-      !
-      CALL iom_setkt( kt + nn_fcrs - 1 )    ! in crs, iom_put is called every nn_fcrs time step(s)
- 
-      IF( MOD( kt-1, nn_fcrs ) == 0 ) THEN
-
-         ! It is time to dump some fields
-         ! compute average by scaling the sum by spatial dimension then time
-         ! OR by time only
-
-         ! 3.1  Ocean fields
-         !    3.1.1 Weight for time average
-         zdtj=rdt/86400.   ! time step in days
-         WRITE(clmode,'(f5.1,a)' ) zdtj,' days average'
-         WRITE(numout,*) 'crsdiawri.', clmode
-         zrtsct = 1.0/REAL(itsct, wp)
-
-         !
-         !    3.1.2 Weights for spatial averages
-         WHERE ( zvolw /= 0 ) zvolw(:,:,:) = 1.0/zvolw(:,:,:)
-         WHERE ( zsurv /= 0 ) zsurv(:,:,:) = 1.0/zsurv(:,:,:)
-         WHERE ( zsuru /= 0 ) zsuru(:,:,:) = 1.0/zsuru(:,:,:)
-         WHERE ( zsurw /= 0 ) zsurw(:,:,:) = 1.0/zsurw(:,:,:)
-         WHERE ( zsurtaux /= 0 ) zsurtaux(:,:) = 1.0 / zsurtaux(:,:)
-         WHERE ( zsurtauy /= 0 ) zsurtauy(:,:) = 1.0 / zsurtauy(:,:)
-
-         !    3.1.3 Temperature
-         tsn_crs(:,:,:,1) = ztsnm(:,:,:,1) * zrtsct ! volume-weighted- , time- mean  
-         z2dcrs(:,:) = tsn_crs(:,:,1,1)
-         CALL crs_iom_put( "toce_crs"   , pv_r3d=tsn_crs(:,:,:,1)  )    ! temperature
-         CALL crs_iom_put( "sst_crs"    , pv_r2d=z2dcrs            )    ! sst
-     
-         !
-         !    3.1.4 Salinity
-         tsn_crs(:,:,:,2) = ztsnm(:,:,:,2) * zrtsct ! volume-weighted- , time- mean 
-         z2dcrs(:,:) = tsn_crs(:,:,1,2)
-         CALL crs_iom_put( "soce_crs"   , pv_r3d=tsn_crs(:,:,:,2)   )   !  salinity
-         CALL crs_iom_put( "sss_crs"    , pv_r2d=z2dcrs             )   !  sss
-         !
-
-         !    3.1.5 U-velocity
-         un_crs(:,:,:) = zum(:,:,:) * zrtsct ! area-weighted- , time- mean 
-         CALL crs_iom_put( "uoce_crs"   , pv_r3d=un_crs             )   ! i-current 
-         ut_crs(:,:,:) = zutm(:,:,:) * zsuru(:,:,:) ! area-weighted- , time- mean      
-         CALL crs_iom_put( "uocet_crs"     , pv_r3d=ut_crs             )   ! uT
-         us_crs(:,:,:) = zusm(:,:,:) * zsuru(:,:,:) ! area-weighted- , time- mean      
-         CALL crs_iom_put( "uoces_crs"     , pv_r3d=us_crs             )   ! uS
-              
-         !    3.1.6 V-velocity
-         vn_crs(:,:,:) = zvm(:,:,:) * zrtsct ! area-weighted- , time- mean 
-         CALL crs_iom_put( "voce_crs"   , pv_r3d=vn_crs                   )    ! j-current
-         vt_crs(:,:,:) = zvtm(:,:,:) * zsurv(:,:,:) ! area-weighted- , time- mean      
-         CALL crs_iom_put( "vocet_crs"     , pv_r3d=vt_crs             )   ! vT
-         vs_crs(:,:,:) = zvsm(:,:,:) * zsurv(:,:,:) ! area-weighted- , time- mean      
-         CALL crs_iom_put( "voces_crs"     , pv_r3d=vs_crs             )   ! vS
-
-         !    3.1.7 W-velocity
-         wn_crs(:,:,:) = zwm(:,:,:) * zrtsct ! area-weighted- , time- mean 
-         CALL crs_iom_put( "woce_crs"   , pv_r3d=wn_crs                   )    ! W-velocity
-         CALL crs_iom_put( "woce2_crs"   , pv_r3d=zw                      )    ! cc
-
-         !    3.1.8 Horizontal divergence
-         hdivn_crs(:,:,:) = zhdivnm(:,:,:) * zrtsct
-      !   CALL crs_iom_put( "hdivn_crs"   , pv_r3d=hdivn_crs            )   ! hdiv
-         CALL crs_iom_put( "hdiv_crs"   , pv_r3d=hdivn_crs            )  
-
-         hdivbt_crs(:,:) = zhdivbtm(:,:) * zrtsct
-
-         !    3.1.9 Sea-surface Height
-         sshn_crs(:,:) = zsshm(:,:)  * zrtsct
-         sshb_crs(:,:) = zsshbm(:,:) * zrtsct
-         ssha_crs(:,:) = zssham(:,:) * zrtsct
-
-         IF( lk_vvl ) THEN                    ! now Sea SSH at u-, v-, f-points (vvl case only)
-            sshun_crs(:,:) = zsshum(:,:)  * zrtsct                              ! sshun 
-            sshub_crs(:,:) = zsshubm(:,:) * zrtsct                              ! sshub 
-            sshua_crs(:,:) = zsshuam(:,:) * zrtsct                              ! sshua 
-            sshvn_crs(:,:) = zsshvm(:,:)  * zrtsct                              ! sshvn 
-            sshvb_crs(:,:) = zsshvbm(:,:) * zrtsct                              ! sshvb 
-            sshva_crs(:,:) = zsshvam(:,:) * zrtsct                              ! sshva 
-            sshfn_crs(:,:) = zsshfm(:,:) * zrtsct                               ! sshfn 
-         ENDIF
-         CALL crs_iom_put( "ssh_crs"   , pv_r2d=sshn_crs                  )   ! ssh output 
-         CALL crs_iom_put( "ssh2_crs"  , pv_r2d=ssh_crs2                  )   !cc
-
-         !    3.1.10 Potential density
-         rhop_crs(:,:,:) = zrhopm(:,:,:) * zrtsct
-         rhd_crs(:,:,:)  =  zrhdm(:,:,:) * zrtsct 
-
-         !    3.1.10 Vertical eddy diffusivity
-         avt_crs(:,:,:) = zavtm(:,:,:) * zsurw(:,:,:) 
-
-         !    3.1.11  Mixed-layer diagnostics
-         !            Recalculated following OPA_SRC/ZDF/zdfmxl.F90
-         !   
-         ! w-level of the mixing and mixed layers
-         nmln_crs(:,:) = mbkt_crs(:,:) + 1        ! Initialization to the number of w ocean point
-         imld(:,:) = mbkt_crs(:,:) + 1
-         DO jk = jpkm1, nlb10, -1         ! from the bottom to nlb10
-            DO jj = 1, jpj_crs
-               DO ji = 1, jpi_crs
-                  IF( rhop_crs(ji,jj,jk) > rhop_crs(ji,jj,nla10) + zrho_c )   nmln_crs(ji,jj) = jk      ! Mixed layer
-                  IF(  avt_crs(ji,jj,jk) < zavt_c                     )           imld(ji,jj) = jk      ! Turbocline
-               END DO
-            END DO
-         END DO
-         ! depth of the mixing and mixed layers
-         DO jj = 1, jpj_crs
-            DO ji = 1, jpi_crs
-               iiki = imld(ji,jj)
-               iikn = nmln_crs(ji,jj)
-               ! jes. TODO.  update here gdep for vvl case
-               hmld_crs (ji,jj) = gdepw_crs(ji,jj,iiki  ) * tmask_crs(ji,jj,1)    ! Turbocline depth
-               hmlp_crs (ji,jj) = gdepw_crs(ji,jj,iikn  ) * tmask_crs(ji,jj,1)    ! Mixed layer depth
-               hmlpt_crs(ji,jj) = gdept_crs(ji,jj,iikn-1)                         ! depth of the last T-point inside the mixed layer
-            END DO
-         END DO
-
-
-         ! 3.2. Surface fluxes
-         !    3.2.1. Evaporation minus Precipitation
-         emp_crs(:,:)   = zempm(:,:)   * zrtsct                     
-         emp_b_crs(:,:) = zemp_bm(:,:) * zrtsct                     
-         emps_crs(:,:)  = zempsm(:,:)  * zrtsct                     
-         !    3.2.2. Wind stresses, Scalar Wind
-         utau_crs(:,:) = zutaum(:,:) * zsurtaux(:,:)
-         vtau_crs(:,:) = zvtaum(:,:) * zsurtauy(:,:)
-         wndm_crs(:,:) = zwndmm(:,:) * zsurw(:,:,1)
-         !    3.2.3. Shortwave radiative flux
-         qsr_crs(:,:) = zqsrm(:,:) * zsurw(:,:,1)
-         !    3.2.4. Runoff freshwater flux
-         rnf_crs(:,:) = zrnfm(:,:) * zsurw(:,:,1)
-         !    3.2.4. Depth of freshwater flux input
-         h_rnf_crs(:,:) = zh_rnfm(:,:) * zrtsct
-         !    3.2.5. Ice fraction
-         fr_i_crs(:,:) = zfr_im(:,:) * zsurw(:,:,1)
-
-
-         ! Reset work arrays and counter to 0 after writing
-         !
-         itsct = 0
-         !
-         zum(:,:,:)   = 0._wp
-         zvm(:,:,:)   = 0._wp
-         ztsnm(:,:,:,:) = 0._wp
-         zwm(:,:,:)   = 0._wp
-         zutm(:,:,:)  = 0._wp
-         zusm(:,:,:)  = 0._wp
-         zvtm(:,:,:)  = 0._wp
-         zvsm(:,:,:)  = 0._wp
-         zwtm(:,:,:)  = 0._wp
-         zwsm(:,:,:)  = 0._wp
-
-         zvolw(:,:,:) = 0._wp
-         zsuru(:,:,:) = 0._wp
-         zsurv(:,:,:) = 0._wp
-         zsurw(:,:,:) = 0._wp   
-         zsurtaux(:,:) = 0._wp
-         zsurtauy(:,:) = 0._wp
-
-         zsshm(:,:) = 0._wp
-         zsshum(:,:) = 0._wp
-         zsshvm(:,:) = 0._wp
-         zsshfm(:,:) = 0._wp
-         zsshbm(:,:) = 0._wp
-         zsshubm(:,:) = 0._wp
-         zsshvbm(:,:) = 0._wp
-         zssham(:,:) = 0._wp
-         zsshuam(:,:) = 0._wp
-         zsshvam(:,:) = 0._wp
-         zwndmm(:,:) = 0._wp
-         zutaum(:,:) = 0._wp
-         zvtaum(:,:) = 0._wp
-         zqsrm(:,:) = 0._wp
-         zrnfm(:,:) = 0._wp
-         zh_rnfm(:,:) = 0._wp
-         zfr_im(:,:) = 0._wp
-
-      ENDIF
-
-      CALL iom_setkt( kt )                  ! iom_put outside of crs is called at every time step
-
-
-
-      ! 5. Clean-up
-      ! --------------------------------------------------------
-      !
-      CALL wrk_dealloc( jpi , jpj , z2d, ze1e2u, ze1e2v )
-      CALL wrk_dealloc( jpi , jpj , jpk , z3d, z3d1 )
-      CALL wrk_dealloc( jpi , jpj , jpk , zfse3u, zfse3v )
-      CALL wrk_dealloc( jpi , jpj , jpk , zfse3t, zfse3w )
-      CALL wrk_dealloc( jpi_crs , jpj_crs , jpk , z3dcrs, zw ) !cc
-      CALL wrk_dealloc( jpi_crs , jpj_crs , z2dcrs, ssh_crs2 )
-
-!!cc1      DEALLOCATE( z2dcrsu, z2dcrsv, z2dcrsf, zhdivbt ) ! probleme de malloc au 65 eme pas de temps
-      DEALLOCATE( z2dcrsf, zhdivbt ) ! probleme de malloc au 65 eme pas de temps
-      DEALLOCATE( zsshub, zsshua, zsshvb, zsshva )
-      DEALLOCATE( zee_t, zee_f, zee_u, zee_v )
-      DEALLOCATE( zmut, zmuf )
-      DEALLOCATE( zcrs_surtaux_wgt, zcrs_surtauy_wgt )
-      DEALLOCATE( zhdiv, zvolt_wgt) 
-      !
-                        
+      CALL crs_lbc_lnk( 'T', 1.0,    pt3d1=hdivn_crs )
+      CALL crs_iom_put( "hdiv_crs", pv_r3d=hdivn_crs )  
+
+      
+      !  Sea-surface Height 
+      CALL crsfun( p_e1e2t=e1e2t, cd_type='T', cd_op='ARE', p_cmask=tmask_crs, p_ptmask=tmask, &
+           &       p_pfield2d=sshn, p_cfield2d=sshn_crs )
+      CALL crs_iom_put( "ssh_crs" , pv_r2d=sshn_crs  )   ! ssh output 
+
+      !  W-velocity
+      CALL crsfun( p_e1e2t=e1e2t, cd_type='T', cd_op='ARE', p_cmask=tmask_crs, p_ptmask=tmask, &
+         &         p_pfield3d_2=wn, p_cfield3d=wn_crs )
+      CALL crs_iom_put( "woce_crs"  , pv_r3d=wn_crs  )   ! i-current 
+      
+
+
+      !  Clean-up
+      
+      CALL wrk_dealloc( jpi, jpj, jpk, zfse3t, zfse3w )
+      CALL wrk_dealloc( jpi, jpj, jpk, zfse3u, zfse3v )
+      CALL wrk_dealloc( jpi, jpj, jpk, zt, zs       )
+      !
+      CALL wrk_dealloc( jpi_crs, jpj_crs, jpk, zt_crs, zs_crs )
+      !
       IF( nn_timing == 1 )   CALL timing_stop('crs_dia_wri')
       !
