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 3864)
+++ branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/CONFIG/ORCA2_LIM_CRS/EXP00/iodef.xml	(revision 3895)
@@ -349,5 +349,5 @@
              <field ref="sss"     name="sosaline" />
              <field ref="ssh"     name="sossheig" />
-             <field ref="eken"         name="energkin"  />
+             <field ref="eken"    name="energkin"  />
           </file>
 
@@ -362,51 +362,4 @@
           </file>
       
-<!-- 
-............................................................................................................
-           example of 3 types of 1d files
-............................................................................................................
--->
-	<!-- automatic definition of the name based on id="1d_grid_T" 
-	     => this name is used as the radical for all file names of this group
-	     => add a suffix to be sure that all files names of this group are different   -->
-	<group id="1d_grid_T" name="auto" description="ocean T grid variables" >
-	  <!-- mooring: automatic definition of the file name suffix based on id="0n180wT"  -->
-	  <file id="0n180wT" name_suffix="auto" >
-	    <group id="0n180wT" zoom_ref="0n180wT" >    <!-- group of variables sharing the same zoom -->
-	      <field ref="toce" name="votemper"  />
-	      <field ref="uoce" name="vozocrtx"  />     <!-- include a U-grid variable in the list => switch to T-grid -->
-              <field ref="eken"         name="energkin"  />
-
-	    </group>
-	  </file>
-	  <!-- Equatorial section: automatic definition of the file name suffix based on id="EqT" -->
-	  <file id="EqT" name_suffix="auto" >
-	    <group id="EqT" zoom_ref="EqT" >
-	      <field ref="toce" name="votemper"  />
-	    </group>
-	  </file>
-	  <!-- global file with different operations on data   -->
-	  <file id="global" >
-	    <field ref="sst"    name="sst_1d_ave"                       />     <!-- mean -->
-	    <field ref="sst"    name="sst_1d_max"  operation="t_max(X)" />     <!-- max -->
-            <field ref="M2x"     name="M2_x_elev"  />
-            <field ref="M2y"     name="M2_y_elev"  />
-	  </file>
-
-	</group>
-
-        <!-- variables available with key_float, instantaneous fields -->
-<!--
-        <file id="floats"  description="floats variables">
-            <field ref="traj_lon"   name="floats_longitude"   freq_op="86400" />
-            <field ref="traj_lat"   name="floats_latitude"    freq_op="86400" />
-            <field ref="traj_dep"   name="floats_depth"       freq_op="86400" />
-            <field ref="traj_temp"  name="floats_temperature" freq_op="86400" />
-            <field ref="traj_salt"  name="floats_salinity"    freq_op="86400" />
-            <field ref="traj_dens"  name="floats_density"     freq_op="86400" />
-            <field ref="traj_group" name="floats_group"       freq_op="86400" />
-        </file>
--->
-
       </group>
 
@@ -422,20 +375,4 @@
 	  <field ref="sss"          name="sosaline"  />
 	  <field ref="ssh"          name="sossheig"  />
-	  <field ref="empmr"        name="sowaflup"  />
-	  <field ref="qsr"          name="soshfldo"  />
-	  <field ref="empsmr"       name="sowaflcd"  />
-	  <field ref="qt"           name="sohefldo"  />
-	  <field ref="mldr10_1"     name="somxl010"  />
-	  <field ref="mldkz5"       name="somixhgt"  />
-	  <field ref="ice_cover"    name="soicecov"  />
-	  <field ref="wspd"         name="sowindsp"  />
-	  <field ref="qrp"          name="sohefldp"  />
-	  <field ref="erp"          name="sowafldp"  />
-	  <field ref="mlddzt"       name="sothedep"  />
-	  <field ref="20d"          name="so20chgt"  />
-	  <field ref="28d"          name="so28chgt"  />
-	  <field ref="hc300"        name="sohtc300"  />
-	  <field ref="ist_ipa"      name="soicetem"  />
-	  <field ref="icealb_cea"   name="soicealb"  />   
 	</file>
 
@@ -448,11 +385,8 @@
            <field ref="ssh_crs"     name="sossheig" />
            <field ref="hdiv_crs"    name="vohdiver" />
-           <field ref="eken_crs"    name="energkin" />
         </file>	
 
 	<file id="5d_grid_U" name="auto" description="ocean U grid variables" >
 	  <field ref="uoce"         name="vozocrtx"  />
-	  <field ref="uoce_eiv"     name="vozoeivu"  />
-	  <field ref="utau"         name="sozotaux"  />
 	</file>
 
@@ -465,6 +399,4 @@
 	<file id="5d_grid_V" name="auto" description="ocean V grid variables" >
 	  <field ref="voce"         name="vomecrty"  />
-	  <field ref="voce_eiv"     name="vomeeivv"  />
-	  <field ref="vtau"         name="sometauy"  />
 	</file>
 
@@ -477,12 +409,4 @@
 	<file id="5d_grid_W" name="auto" description="ocean W grid variables" >
 	  <field ref="woce"         name="vovecrtz" />
-	  <field ref="woce_eiv"     name="voveeivw" />
-	  <field ref="avt"          name="votkeavt" />
-	  <field ref="avt_evd"      name="votkeevd" />
-	  <field ref="avm"          name="votkeavm" />
-	  <field ref="avm_evd"      name="votkeevm" />
-	  <field ref="avs"          name="voddmavs" />
-	  <field ref="aht2d"        name="soleahtw" />
-	  <field ref="aht2d_eiv"    name="soleaeiw" />
 	</file>
 
@@ -491,28 +415,7 @@
         </file>	
 	
-	<file id="5d_icemod" name="auto" description="ice variables" >
-	  <field ref="ice_pres"                     />
-	  <field ref="snowthic_cea" name="isnowthi" />
-	  <field ref="icethic_cea"  name="iicethic" />
-	  <field ref="iceprod_cea"  name="iiceprod" />
-	  <field ref="ist_ipa"      name="iicetemp" />
-	  <field ref="ioceflxb"     name="ioceflxb" />
-	  <field ref="uice_ipa"     name="iicevelu" />
-	  <field ref="vice_ipa"     name="iicevelv" />
-	  <field ref="utau_ice"     name="iicestru" />
-	  <field ref="vtau_ice"     name="iicestrv" />
-	  <field ref="qsr_io_cea"   name="iicesflx" />
-	  <field ref="qns_io_cea"   name="iicenflx" />
-	  <field ref="snowpre"      name="isnowpre" />
-	</file>
-	
       </group>
       
       <group id="1m" output_freq="-1"     output_level="10" enabled=".FALSE.">                      <!-- real monthly files -->
-
-	<file id="1m_grid_T" name="auto" description="ocean T grid variables" >
-	  <field ref="sst"          name="sosstsst"  />
-	</file>   
-
       </group>
 
@@ -530,9 +433,4 @@
 
       <group id="1y" output_freq="-12"    output_level="10" enabled=".FALSE.">                      <!-- real yearly files -->
-
-	<file id="1y_grid_T" name="auto" description="ocean T grid variables" >
-	  <field ref="mldr10_1"     name="sobowlin"  operation="t_max(X)" />
-	</file>
-
       </group>
 
Index: branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/CONFIG/ORCA2_LIM_CRS/MY_SRC/lbcnfd.F90
===================================================================
--- branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/CONFIG/ORCA2_LIM_CRS/MY_SRC/lbcnfd.F90	(revision 3895)
+++ branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/CONFIG/ORCA2_LIM_CRS/MY_SRC/lbcnfd.F90	(revision 3895)
@@ -0,0 +1,365 @@
+MODULE lbcnfd
+   !!======================================================================
+   !!                       ***  MODULE  lbcnfd  ***
+   !! Ocean        : north fold  boundary conditions
+   !!======================================================================
+   !! History :  3.2  ! 2009-03  (R. Benshila)  Original code 
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   lbc_nfd       : generic interface for lbc_nfd_3d and lbc_nfd_2d routines
+   !!   lbc_nfd_3d    : lateral boundary condition: North fold treatment for a 3D arrays   (lbc_nfd)
+   !!   lbc_nfd_2d    : lateral boundary condition: North fold treatment for a 2D arrays   (lbc_nfd)
+   !!----------------------------------------------------------------------
+   USE dom_oce        ! ocean space and time domain 
+   USE in_out_manager ! I/O manager
+
+   IMPLICIT NONE
+   PRIVATE
+
+   INTERFACE lbc_nfd
+      MODULE PROCEDURE   lbc_nfd_3d, lbc_nfd_2d
+   END INTERFACE
+
+   PUBLIC   lbc_nfd   ! north fold conditions
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: lbcnfd.F90 3294 2012-01-28 16:44:18Z rblod $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE lbc_nfd_3d( pt3d, cd_type, psgn )
+      !!----------------------------------------------------------------------
+      !!                  ***  routine lbc_nfd_3d  ***
+      !!
+      !! ** Purpose :   3D lateral boundary condition : North fold treatment
+      !!              without processor exchanges. 
+      !!
+      !! ** Method  :   
+      !!
+      !! ** Action  :   pt3d with updated values along the north fold
+      !!----------------------------------------------------------------------
+      CHARACTER(len=1)          , INTENT(in   ) ::   cd_type   ! define the nature of ptab array grid-points
+      !                                                        !   = T , U , V , F , W points
+      REAL(wp)                  , INTENT(in   ) ::   psgn      ! control of the sign change
+      !                                                        !   = -1. , the sign is changed if north fold boundary
+      !                                                        !   =  1. , the sign is kept  if north fold boundary
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pt3d      ! 3D array on which the boundary condition is applied
+      !
+      INTEGER  ::   ji, jk
+      INTEGER  ::   ijt, iju, ijpj, ijpjm1
+      !!----------------------------------------------------------------------
+
+      SELECT CASE ( jpni )
+      CASE ( 1 )     ;   ijpj = nlcj      ! 1 proc only  along the i-direction
+      CASE DEFAULT   ;   ijpj = 4         ! several proc along the i-direction
+      END SELECT
+      ijpjm1 = ijpj-1
+
+      DO jk = 1, jpk
+         !
+         SELECT CASE ( npolj )
+         !
+         CASE ( 3 , 4 )                        ! *  North fold  T-point pivot
+            !                                        (the 2 North Poles are at (i,j) = (2,jpjglo-1) and (jpiglo/2+1,jpjglo-1)
+            SELECT CASE ( cd_type )
+            CASE ( 'T' , 'W' )                         ! T-, W-point
+               DO ji = 2, jpiglo
+                  ijt = jpiglo-ji+2
+                  pt3d(ji,ijpj,jk) = psgn * pt3d(ijt,ijpj-2,jk)
+               END DO
+               pt3d( 1 ,ijpj,jk) = psgn * pt3d( 3 ,ijpj-2,jk)
+               DO ji = jpiglo/2+1, jpiglo
+                  ijt = jpiglo-ji+2
+                  pt3d(ji,ijpjm1,jk) = psgn * pt3d(ijt,ijpjm1,jk)
+               END DO
+            CASE ( 'U' )                               ! U-point
+               DO ji = 1, jpiglo-1
+                  iju = jpiglo-ji+1
+                  pt3d(ji,ijpj,jk) = psgn * pt3d(iju,ijpj-2,jk)
+               END DO
+               pt3d(   1  ,ijpj,jk) = psgn * pt3d(    2   ,ijpj-2,jk)
+               pt3d(jpiglo,ijpj,jk) = psgn * pt3d(jpiglo-1,ijpj-2,jk) 
+               DO ji = jpiglo/2, jpiglo-1
+                  iju = jpiglo-ji+1
+                  pt3d(ji,ijpjm1,jk) = psgn * pt3d(iju,ijpjm1,jk)
+               END DO
+            CASE ( 'V' )                               ! V-point
+               DO ji = 2, jpiglo
+                  ijt = jpiglo-ji+2
+                  pt3d(ji,ijpjm1,jk) = psgn * pt3d(ijt,ijpj-2,jk)
+                  pt3d(ji,ijpj  ,jk) = psgn * pt3d(ijt,ijpj-3,jk)
+               END DO
+               pt3d( 1 ,ijpj,jk) = psgn * pt3d( 3 ,ijpj-3,jk) 
+            CASE ( 'F' )                               ! F-point
+               DO ji = 1, jpiglo-1
+                  iju = jpiglo-ji+1
+                  pt3d(ji,ijpj-1,jk) = psgn * pt3d(iju,ijpj-2,jk)
+                  pt3d(ji,ijpj  ,jk) = psgn * pt3d(iju,ijpj-3,jk)
+               END DO
+               pt3d(   1  ,ijpj,jk) = psgn * pt3d(    2   ,ijpj-3,jk)
+               pt3d(jpiglo,ijpj,jk) = psgn * pt3d(jpiglo-1,ijpj-3,jk) 
+            END SELECT
+            !
+         CASE ( 5 , 6 )                        ! *  North fold  F-point pivot
+            !
+            SELECT CASE ( cd_type )
+            CASE ( 'T' , 'W' )                         ! T-, W-point
+               DO ji = 1, jpiglo
+                  ijt = jpiglo-ji+1
+                  pt3d(ji,ijpj,jk) = psgn * pt3d(ijt,ijpj-1,jk)
+               END DO
+            CASE ( 'U' )                               ! U-point
+               DO ji = 1, jpiglo-1
+                  iju = jpiglo-ji
+                  pt3d(ji,ijpj,jk) = psgn * pt3d(iju,ijpj-1,jk)
+               END DO
+               pt3d(jpiglo,ijpj,jk) = psgn * pt3d(1,ijpj-1,jk)
+            CASE ( 'V' )                               ! V-point
+               DO ji = 1, jpiglo
+                  ijt = jpiglo-ji+1
+                  pt3d(ji,ijpj,jk) = psgn * pt3d(ijt,ijpj-2,jk)
+               END DO
+               DO ji = jpiglo/2+1, jpiglo
+                  ijt = jpiglo-ji+1
+                  pt3d(ji,ijpjm1,jk) = psgn * pt3d(ijt,ijpjm1,jk)
+               END DO
+            CASE ( 'F' )                               ! F-point
+               DO ji = 1, jpiglo-1
+                  iju = jpiglo-ji
+                  pt3d(ji,ijpj  ,jk) = psgn * pt3d(iju,ijpj-2,jk)
+               END DO
+               pt3d(jpiglo,ijpj,jk) = psgn * pt3d(1,ijpj-2,jk)
+               DO ji = jpiglo/2+1, jpiglo-1
+                  iju = jpiglo-ji
+                  pt3d(ji,ijpjm1,jk) = psgn * pt3d(iju,ijpjm1,jk)
+               END DO
+            END SELECT
+            !
+         CASE DEFAULT                           ! *  closed : the code probably never go through
+            !
+            SELECT CASE ( cd_type)
+            CASE ( 'T' , 'U' , 'V' , 'W' )             ! T-, U-, V-, W-points
+               pt3d(:, 1  ,jk) = 0.e0
+               pt3d(:,ijpj,jk) = 0.e0
+            CASE ( 'F' )                               ! F-point
+               pt3d(:,ijpj,jk) = 0.e0
+            END SELECT
+            !
+         END SELECT     !  npolj
+         !
+      END DO
+      !
+   END SUBROUTINE lbc_nfd_3d
+
+
+   SUBROUTINE lbc_nfd_2d( pt2d, cd_type, psgn, pr2dj )
+      !!----------------------------------------------------------------------
+      !!                  ***  routine lbc_nfd_2d  ***
+      !!
+      !! ** Purpose :   2D lateral boundary condition : North fold treatment
+      !!       without processor exchanges. 
+      !!
+      !! ** Method  :   
+      !!
+      !! ** Action  :   pt2d with updated values along the north fold
+      !!----------------------------------------------------------------------
+      CHARACTER(len=1)        , INTENT(in   ) ::   cd_type   ! define the nature of ptab array grid-points
+      !                                                      ! = T , U , V , F , W points
+      REAL(wp)                , INTENT(in   ) ::   psgn      ! control of the sign change
+      !                                                      !   = -1. , the sign is changed if north fold boundary
+      !                                                      !   =  1. , the sign is kept  if north fold boundary
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt2d      ! 2D array on which the boundary condition is applied
+      INTEGER , OPTIONAL      , INTENT(in   ) ::   pr2dj     ! number of additional halos
+      !
+      INTEGER  ::   ji, jl, ipr2dj
+      INTEGER  ::   ijt, iju, ijpj, ijpjm1
+      !!----------------------------------------------------------------------
+
+      SELECT CASE ( jpni )
+      CASE ( 1 )     ;   ijpj = nlcj      ! 1 proc only  along the i-direction
+      CASE DEFAULT   ;   ijpj = 4         ! several proc along the i-direction
+      END SELECT
+      !
+      IF( PRESENT(pr2dj) ) THEN           ! use of additional halos
+         ipr2dj = pr2dj
+         IF( jpni > 1 )   ijpj = ijpj + ipr2dj
+      ELSE
+         ipr2dj = 0 
+      ENDIF
+      !
+      ijpjm1 = ijpj-1
+
+      SELECT CASE ( npolj )
+      !
+      CASE ( 3, 4 )                       ! *  North fold  T-point pivot
+         !
+         SELECT CASE ( cd_type )
+         !
+         CASE ( 'T' , 'W' )                               ! T- , W-points
+            DO jl = 0, ipr2dj
+               DO ji = 2, jpiglo
+                  ijt=jpiglo-ji+2
+                  pt2d(ji,ijpj+jl) = psgn * pt2d(ijt,ijpj-2-jl)
+               END DO
+            END DO
+            pt2d( 1 ,ijpj)   = psgn * pt2d( 3 ,ijpj-2)
+            DO ji = jpiglo/2+1, jpiglo
+               ijt=jpiglo-ji+2
+               pt2d(ji,ijpj-1) = psgn * pt2d(ijt,ijpj-1)
+            END DO
+         CASE ( 'U' )                                     ! U-point
+            DO jl = 0, ipr2dj
+               DO ji = 2, jpiglo-1
+                  iju = jpiglo-ji+1
+                  pt2d(ji,ijpj+jl) = psgn * pt2d(iju,ijpj-2-jl)
+               END DO
+            END DO
+            pt2d(   1  ,ijpj) = psgn * pt2d(    2   ,ijpj-2)
+            pt2d(jpiglo,ijpj) = psgn * pt2d(jpiglo-1,ijpj-2)
+            pt2d(1,ijpj-1) = psgn * pt2d(jpiglo,ijpj-1)   
+            
+            DO ji = jpiglo/2, jpiglo-1
+               iju = jpiglo-ji+1
+               pt2d(ji,ijpjm1) = psgn * pt2d(iju,ijpjm1)
+            END DO
+         CASE ( 'V' )                                     ! V-point
+            DO jl = -1, ipr2dj
+               DO ji = 2, jpiglo
+                  ijt = jpiglo-ji+2
+                  pt2d(ji,ijpj+jl)   = psgn * pt2d(ijt,ijpj-3-jl)
+               END DO
+            END DO
+            pt2d( 1 ,ijpj)   = psgn * pt2d( 3 ,ijpj-3) 
+            
+         CASE ( 'F' )                                     ! F-point
+            DO jl = -1, ipr2dj
+               DO ji = 2, jpiglo-1
+                  iju = jpiglo-ji+1
+                  pt2d(ji,ijpj+jl) = psgn * pt2d(iju,ijpj-3-jl)
+               END DO
+            END DO
+            pt2d(   1  ,ijpj)   = psgn * pt2d(    2   ,ijpj-3)
+            pt2d(jpiglo,ijpj)   = psgn * pt2d(jpiglo-1,ijpj-3)
+            pt2d(jpiglo,ijpj-1) = psgn * pt2d(jpiglo-1,ijpj-2)      
+            pt2d(   1  ,ijpj-1) = psgn * pt2d(    2   ,ijpj-2)      
+         CASE ( 'I' )                                     ! ice U-V point (I-point)
+            DO jl = 0, ipr2dj
+               pt2d(2,ijpj+jl) = psgn * pt2d(3,ijpj-1+jl)
+               DO ji = 3, jpiglo
+                  iju = jpiglo - ji + 3
+                  pt2d(ji,ijpj+jl) = psgn * pt2d(iju,ijpj-1-jl)
+               END DO
+            END DO
+         CASE ( 'J' )                                     ! first ice U-V point
+            DO jl =0, ipr2dj
+               pt2d(2,ijpj+jl) = psgn * pt2d(3,ijpj-1+jl)
+               DO ji = 3, jpiglo
+                  iju = jpiglo - ji + 3
+                  pt2d(ji,ijpj+jl) = psgn * pt2d(iju,ijpj-1-jl)
+               END DO
+            END DO
+         CASE ( 'K' )                                     ! second ice U-V point
+            DO jl =0, ipr2dj
+               pt2d(2,ijpj+jl) = psgn * pt2d(3,ijpj-1+jl)
+               DO ji = 3, jpiglo
+                  iju = jpiglo - ji + 3
+                  pt2d(ji,ijpj+jl) = psgn * pt2d(iju,ijpj-1-jl)
+               END DO
+            END DO
+         END SELECT
+         !
+      CASE ( 5, 6 )                        ! *  North fold  F-point pivot
+         !
+         SELECT CASE ( cd_type )
+         CASE ( 'T' , 'W' )                               ! T-, W-point
+            DO jl = 0, ipr2dj
+               DO ji = 1, jpiglo
+                  ijt = jpiglo-ji+1
+                  pt2d(ji,ijpj+jl) = psgn * pt2d(ijt,ijpj-1-jl)
+               END DO
+            END DO
+         CASE ( 'U' )                                     ! U-point
+            DO jl = 0, ipr2dj
+               DO ji = 1, jpiglo-1
+                  iju = jpiglo-ji
+                  pt2d(ji,ijpj+jl) = psgn * pt2d(iju,ijpj-1-jl)
+               END DO
+            END DO
+            pt2d(jpiglo,ijpj) = psgn * pt2d(1,ijpj-1)
+         CASE ( 'V' )                                     ! V-point
+            DO jl = 0, ipr2dj
+               DO ji = 1, jpiglo
+                  ijt = jpiglo-ji+1
+                  pt2d(ji,ijpj+jl) = psgn * pt2d(ijt,ijpj-2-jl)
+               END DO
+            END DO
+            DO ji = jpiglo/2+1, jpiglo
+               ijt = jpiglo-ji+1
+               pt2d(ji,ijpjm1) = psgn * pt2d(ijt,ijpjm1)
+            END DO
+         CASE ( 'F' )                               ! F-point
+            DO jl = 0, ipr2dj
+               DO ji = 1, jpiglo-1
+                  iju = jpiglo-ji
+                  pt2d(ji,ijpj+jl) = psgn * pt2d(iju,ijpj-2-jl)
+               END DO
+            END DO
+            pt2d(jpiglo,ijpj) = psgn * pt2d(1,ijpj-2)
+            DO ji = jpiglo/2+1, jpiglo-1
+               iju = jpiglo-ji
+               pt2d(ji,ijpjm1) = psgn * pt2d(iju,ijpjm1)
+            END DO
+         CASE ( 'I' )                                  ! ice U-V point (I-point)
+            pt2d( 2 ,ijpj:ijpj+ipr2dj) = 0.e0
+            DO jl = 0, ipr2dj
+               DO ji = 2 , jpiglo-1
+                  ijt = jpiglo - ji + 2
+                  pt2d(ji,ijpj+jl)= 0.5 * ( pt2d(ji,ijpj-1-jl) + psgn * pt2d(ijt,ijpj-1-jl) )
+               END DO
+            END DO
+         CASE ( 'J' )                                  ! first ice U-V point
+            pt2d( 2 ,ijpj:ijpj+ipr2dj) = 0.e0
+            DO jl = 0, ipr2dj
+               DO ji = 2 , jpiglo-1
+                  ijt = jpiglo - ji + 2
+                  pt2d(ji,ijpj+jl)= pt2d(ji,ijpj-1-jl)
+               END DO
+            END DO
+         CASE ( 'K' )                                  ! second ice U-V point
+            pt2d( 2 ,ijpj:ijpj+ipr2dj) = 0.e0
+            DO jl = 0, ipr2dj
+               DO ji = 2 , jpiglo-1
+                  ijt = jpiglo - ji + 2
+                  pt2d(ji,ijpj+jl)= pt2d(ijt,ijpj-1-jl)
+               END DO
+            END DO
+         END SELECT
+         !
+      CASE DEFAULT                           ! *  closed : the code probably never go through
+         !
+         SELECT CASE ( cd_type)
+         CASE ( 'T' , 'U' , 'V' , 'W' )                 ! T-, U-, V-, W-points
+            pt2d(:, 1:1-ipr2dj     ) = 0.e0
+            pt2d(:,ijpj:ijpj+ipr2dj) = 0.e0
+         CASE ( 'F' )                                   ! F-point
+            pt2d(:,ijpj:ijpj+ipr2dj) = 0.e0
+         CASE ( 'I' )                                   ! ice U-V point
+            pt2d(:, 1:1-ipr2dj     ) = 0.e0
+            pt2d(:,ijpj:ijpj+ipr2dj) = 0.e0
+         CASE ( 'J' )                                   ! first ice U-V point
+            pt2d(:, 1:1-ipr2dj     ) = 0.e0
+            pt2d(:,ijpj:ijpj+ipr2dj) = 0.e0
+         CASE ( 'K' )                                   ! second ice U-V point
+            pt2d(:, 1:1-ipr2dj     ) = 0.e0
+            pt2d(:,ijpj:ijpj+ipr2dj) = 0.e0
+         END SELECT
+         !
+      END SELECT
+      !
+   END SUBROUTINE lbc_nfd_2d
+
+   !!======================================================================
+END MODULE lbcnfd
Index: branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/CONFIG/ORCA2_LIM_CRS/MY_SRC/lib_mpp.F90
===================================================================
--- branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/CONFIG/ORCA2_LIM_CRS/MY_SRC/lib_mpp.F90	(revision 3864)
+++ branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/CONFIG/ORCA2_LIM_CRS/MY_SRC/lib_mpp.F90	(revision 3895)
@@ -61,5 +61,5 @@
    PUBLIC   ctl_stop, ctl_warn, get_unit, ctl_opn
    PUBLIC   mynode, mppstop, mppsync, mpp_comm_free
-   PUBLIC   mpp_ini_north
+   PUBLIC   mpp_ini_north, mpp_lbc_north, mpp_lbc_north_e
    PUBLIC   mpp_min, mpp_max, mpp_sum, mpp_minloc, mpp_maxloc
    PUBLIC   mpp_lnk_3d, mpp_lnk_3d_gather, mpp_lnk_2d, mpp_lnk_2d_e
@@ -328,11 +328,11 @@
       REAL(wp) ::   zland
       INTEGER, DIMENSION(MPI_STATUS_SIZE) ::   ml_stat   ! for key_mpi_isend
-      REAL(wp), DIMENSION(:,:,:,:)  , ALLOCATABLE ::   zt3ns, zt3sn   ! 3d for north-south & south-north
-      REAL(wp), DIMENSION(:,:,:,:)  , ALLOCATABLE ::   zt3ew, zt3we   ! 3d for east-west & west-east
+      REAL(wp), DIMENSION(:,:,:,:)  , ALLOCATABLE ::   t3ns, t3sn   ! 3d for north-south & south-north
+      REAL(wp), DIMENSION(:,:,:,:)  , ALLOCATABLE ::   t3ew, t3we   ! 3d for east-west & west-east
 
       !!----------------------------------------------------------------------
       
-      ALLOCATE( zt3ns(jpi,jprecj,jpk,2)   , zt3sn(jpi,jprecj,jpk,2)   ,   &
-         &      zt3ew(jpj,jpreci,jpk,2)   , zt3we(jpj,jpreci,jpk,2))
+      ALLOCATE( t3ns(jpi,jprecj,jpk,2)   , t3sn(jpi,jprecj,jpk,2)   ,   &
+         &      t3ew(jpj,jpreci,jpk,2)   , t3we(jpj,jpreci,jpk,2))
          
       IF( PRESENT( pval ) ) THEN   ;   zland = pval      ! set land value
@@ -383,6 +383,6 @@
          iihom = nlci-nreci
          DO jl = 1, jpreci
-            zt3ew(:,jl,:,1) = ptab(jpreci+jl,:,:)
-            zt3we(:,jl,:,1) = ptab(iihom +jl,:,:)
+            t3ew(:,jl,:,1) = ptab(jpreci+jl,:,:)
+            t3we(:,jl,:,1) = ptab(iihom +jl,:,:)
          END DO
       END SELECT  
@@ -393,17 +393,17 @@
       SELECT CASE ( nbondi ) 
       CASE ( -1 )
-         CALL mppsend( 2, zt3we(1,1,1,1), imigr, noea, ml_req1 )
-         CALL mpprecv( 1, zt3ew(1,1,1,2), imigr, noea )
+         CALL mppsend( 2, t3we(1,1,1,1), imigr, noea, ml_req1 )
+         CALL mpprecv( 1, t3ew(1,1,1,2), imigr, noea )
          IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
       CASE ( 0 )
-         CALL mppsend( 1, zt3ew(1,1,1,1), imigr, nowe, ml_req1 )
-         CALL mppsend( 2, zt3we(1,1,1,1), imigr, noea, ml_req2 )
-         CALL mpprecv( 1, zt3ew(1,1,1,2), imigr, noea )
-         CALL mpprecv( 2, zt3we(1,1,1,2), imigr, nowe )
+         CALL mppsend( 1, t3ew(1,1,1,1), imigr, nowe, ml_req1 )
+         CALL mppsend( 2, t3we(1,1,1,1), imigr, noea, ml_req2 )
+         CALL mpprecv( 1, t3ew(1,1,1,2), imigr, noea )
+         CALL mpprecv( 2, t3we(1,1,1,2), imigr, nowe )
          IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
          IF(l_isend) CALL mpi_wait(ml_req2, ml_stat, ml_err)
       CASE ( 1 )
-         CALL mppsend( 1, zt3ew(1,1,1,1), imigr, nowe, ml_req1 )
-         CALL mpprecv( 2, zt3we(1,1,1,2), imigr, nowe )
+         CALL mppsend( 1, t3ew(1,1,1,1), imigr, nowe, ml_req1 )
+         CALL mpprecv( 2, t3we(1,1,1,2), imigr, nowe )
          IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
       END SELECT
@@ -415,14 +415,14 @@
       CASE ( -1 )
          DO jl = 1, jpreci
-            ptab(iihom+jl,:,:) = zt3ew(:,jl,:,2)
+            ptab(iihom+jl,:,:) = t3ew(:,jl,:,2)
          END DO
       CASE ( 0 ) 
          DO jl = 1, jpreci
-            ptab(jl      ,:,:) = zt3we(:,jl,:,2)
-            ptab(iihom+jl,:,:) = zt3ew(:,jl,:,2)
+            ptab(jl      ,:,:) = t3we(:,jl,:,2)
+            ptab(iihom+jl,:,:) = t3ew(:,jl,:,2)
          END DO
       CASE ( 1 )
          DO jl = 1, jpreci
-            ptab(jl      ,:,:) = zt3we(:,jl,:,2)
+            ptab(jl      ,:,:) = t3we(:,jl,:,2)
          END DO
       END SELECT
@@ -436,6 +436,6 @@
          ijhom = nlcj-nrecj
          DO jl = 1, jprecj
-            zt3sn(:,jl,:,1) = ptab(:,ijhom +jl,:)
-            zt3ns(:,jl,:,1) = ptab(:,jprecj+jl,:)
+            t3sn(:,jl,:,1) = ptab(:,ijhom +jl,:)
+            t3ns(:,jl,:,1) = ptab(:,jprecj+jl,:)
          END DO
       ENDIF
@@ -446,17 +446,17 @@
       SELECT CASE ( nbondj )     
       CASE ( -1 )
-         CALL mppsend( 4, zt3sn(1,1,1,1), imigr, nono, ml_req1 )
-         CALL mpprecv( 3, zt3ns(1,1,1,2), imigr, nono )
+         CALL mppsend( 4, t3sn(1,1,1,1), imigr, nono, ml_req1 )
+         CALL mpprecv( 3, t3ns(1,1,1,2), imigr, nono )
          IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
       CASE ( 0 )
-         CALL mppsend( 3, zt3ns(1,1,1,1), imigr, noso, ml_req1 )
-         CALL mppsend( 4, zt3sn(1,1,1,1), imigr, nono, ml_req2 )
-         CALL mpprecv( 3, zt3ns(1,1,1,2), imigr, nono )
-         CALL mpprecv( 4, zt3sn(1,1,1,2), imigr, noso )
+         CALL mppsend( 3, t3ns(1,1,1,1), imigr, noso, ml_req1 )
+         CALL mppsend( 4, t3sn(1,1,1,1), imigr, nono, ml_req2 )
+         CALL mpprecv( 3, t3ns(1,1,1,2), imigr, nono )
+         CALL mpprecv( 4, t3sn(1,1,1,2), imigr, noso )
          IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
          IF(l_isend) CALL mpi_wait(ml_req2, ml_stat, ml_err)
       CASE ( 1 ) 
-         CALL mppsend( 3, zt3ns(1,1,1,1), imigr, noso, ml_req1 )
-         CALL mpprecv( 4, zt3sn(1,1,1,2), imigr, noso )
+         CALL mppsend( 3, t3ns(1,1,1,1), imigr, noso, ml_req1 )
+         CALL mpprecv( 4, t3sn(1,1,1,2), imigr, noso )
          IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
       END SELECT
@@ -468,14 +468,14 @@
       CASE ( -1 )
          DO jl = 1, jprecj
-            ptab(:,ijhom+jl,:) = zt3ns(:,jl,:,2)
+            ptab(:,ijhom+jl,:) = t3ns(:,jl,:,2)
          END DO
       CASE ( 0 ) 
          DO jl = 1, jprecj
-            ptab(:,jl      ,:) = zt3sn(:,jl,:,2)
-            ptab(:,ijhom+jl,:) = zt3ns(:,jl,:,2)
+            ptab(:,jl      ,:) = t3sn(:,jl,:,2)
+            ptab(:,ijhom+jl,:) = t3ns(:,jl,:,2)
          END DO
       CASE ( 1 )
          DO jl = 1, jprecj
-            ptab(:,jl,:) = zt3sn(:,jl,:,2)
+            ptab(:,jl,:) = t3sn(:,jl,:,2)
          END DO
       END SELECT
@@ -494,5 +494,5 @@
       ENDIF
       !
-      DEALLOCATE( zt3ns   , zt3sn,  zt3ew   , zt3we )
+      DEALLOCATE( t3ns   , t3sn,  t3ew   , t3we )
    END SUBROUTINE mpp_lnk_3d
 
@@ -530,10 +530,10 @@
       REAL(wp) ::   zland
       INTEGER, DIMENSION(MPI_STATUS_SIZE) ::   ml_stat   ! for key_mpi_isend
-      REAL(wp), DIMENSION(:,:,:)    , ALLOCATABLE ::   zt2ns, zt2sn   ! 2d for north-south & south-north
-      REAL(wp), DIMENSION(:,:,:)    , ALLOCATABLE ::   zt2ew, zt2we   ! 2d for east-west & west-east
-
-      !!----------------------------------------------------------------------
-      ALLOCATE( zt2ns(jpi,jprecj    ,2)   , zt2sn(jpi,jprecj    ,2)   ,   &
-         &      zt2ew(jpj,jpreci    ,2)   , zt2we(jpj,jpreci    ,2))
+      REAL(wp), DIMENSION(:,:,:)    , ALLOCATABLE ::   t2ns, t2sn   ! 2d for north-south & south-north
+      REAL(wp), DIMENSION(:,:,:)    , ALLOCATABLE ::   t2ew, t2we   ! 2d for east-west & west-east
+
+      !!----------------------------------------------------------------------
+      ALLOCATE( t2ns(jpi,jprecj    ,2)   , t2sn(jpi,jprecj    ,2)   ,   &
+         &      t2ew(jpj,jpreci    ,2)   , t2we(jpj,jpreci    ,2))
       
       IF( PRESENT( pval ) ) THEN   ;   zland = pval      ! set land value
@@ -583,6 +583,6 @@
          iihom = nlci-nreci
          DO jl = 1, jpreci
-            zt2ew(:,jl,1) = pt2d(jpreci+jl,:)
-            zt2we(:,jl,1) = pt2d(iihom +jl,:)
+            t2ew(:,jl,1) = pt2d(jpreci+jl,:)
+            t2we(:,jl,1) = pt2d(iihom +jl,:)
          END DO
       END SELECT
@@ -593,17 +593,17 @@
       SELECT CASE ( nbondi )
       CASE ( -1 )
-         CALL mppsend( 2, zt2we(1,1,1), imigr, noea, ml_req1 )
-         CALL mpprecv( 1, zt2ew(1,1,2), imigr, noea )
+         CALL mppsend( 2, t2we(1,1,1), imigr, noea, ml_req1 )
+         CALL mpprecv( 1, t2ew(1,1,2), imigr, noea )
          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
       CASE ( 0 )
-         CALL mppsend( 1, zt2ew(1,1,1), imigr, nowe, ml_req1 )
-         CALL mppsend( 2, zt2we(1,1,1), imigr, noea, ml_req2 )
-         CALL mpprecv( 1, zt2ew(1,1,2), imigr, noea )
-         CALL mpprecv( 2, zt2we(1,1,2), imigr, nowe )
+         CALL mppsend( 1, t2ew(1,1,1), imigr, nowe, ml_req1 )
+         CALL mppsend( 2, t2we(1,1,1), imigr, noea, ml_req2 )
+         CALL mpprecv( 1, t2ew(1,1,2), imigr, noea )
+         CALL mpprecv( 2, t2we(1,1,2), imigr, nowe )
          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
          IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err)
       CASE ( 1 )
-         CALL mppsend( 1, zt2ew(1,1,1), imigr, nowe, ml_req1 )
-         CALL mpprecv( 2, zt2we(1,1,2), imigr, nowe )
+         CALL mppsend( 1, t2ew(1,1,1), imigr, nowe, ml_req1 )
+         CALL mpprecv( 2, t2we(1,1,2), imigr, nowe )
          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
       END SELECT
@@ -615,14 +615,14 @@
       CASE ( -1 )
          DO jl = 1, jpreci
-            pt2d(iihom+jl,:) = zt2ew(:,jl,2)
+            pt2d(iihom+jl,:) = t2ew(:,jl,2)
          END DO
       CASE ( 0 )
          DO jl = 1, jpreci
-            pt2d(jl      ,:) = zt2we(:,jl,2)
-            pt2d(iihom+jl,:) = zt2ew(:,jl,2)
+            pt2d(jl      ,:) = t2we(:,jl,2)
+            pt2d(iihom+jl,:) = t2ew(:,jl,2)
          END DO
       CASE ( 1 )
          DO jl = 1, jpreci
-            pt2d(jl      ,:) = zt2we(:,jl,2)
+            pt2d(jl      ,:) = t2we(:,jl,2)
          END DO
       END SELECT
@@ -636,6 +636,6 @@
          ijhom = nlcj-nrecj
          DO jl = 1, jprecj
-            zt2sn(:,jl,1) = pt2d(:,ijhom +jl)
-            zt2ns(:,jl,1) = pt2d(:,jprecj+jl)
+            t2sn(:,jl,1) = pt2d(:,ijhom +jl)
+            t2ns(:,jl,1) = pt2d(:,jprecj+jl)
          END DO
       ENDIF
@@ -646,17 +646,17 @@
       SELECT CASE ( nbondj )
       CASE ( -1 )
-         CALL mppsend( 4, zt2sn(1,1,1), imigr, nono, ml_req1 )
-         CALL mpprecv( 3, zt2ns(1,1,2), imigr, nono )
+         CALL mppsend( 4, t2sn(1,1,1), imigr, nono, ml_req1 )
+         CALL mpprecv( 3, t2ns(1,1,2), imigr, nono )
          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
       CASE ( 0 )
-         CALL mppsend( 3, zt2ns(1,1,1), imigr, noso, ml_req1 )
-         CALL mppsend( 4, zt2sn(1,1,1), imigr, nono, ml_req2 )
-         CALL mpprecv( 3, zt2ns(1,1,2), imigr, nono )
-         CALL mpprecv( 4, zt2sn(1,1,2), imigr, noso )
+         CALL mppsend( 3, t2ns(1,1,1), imigr, noso, ml_req1 )
+         CALL mppsend( 4, t2sn(1,1,1), imigr, nono, ml_req2 )
+         CALL mpprecv( 3, t2ns(1,1,2), imigr, nono )
+         CALL mpprecv( 4, t2sn(1,1,2), imigr, noso )
          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
          IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err)
       CASE ( 1 )
-         CALL mppsend( 3, zt2ns(1,1,1), imigr, noso, ml_req1 )
-         CALL mpprecv( 4, zt2sn(1,1,2), imigr, noso )
+         CALL mppsend( 3, t2ns(1,1,1), imigr, noso, ml_req1 )
+         CALL mpprecv( 4, t2sn(1,1,2), imigr, noso )
          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
       END SELECT
@@ -668,14 +668,14 @@
       CASE ( -1 )
          DO jl = 1, jprecj
-            pt2d(:,ijhom+jl) = zt2ns(:,jl,2)
+            pt2d(:,ijhom+jl) = t2ns(:,jl,2)
          END DO
       CASE ( 0 )
          DO jl = 1, jprecj
-            pt2d(:,jl      ) = zt2sn(:,jl,2)
-            pt2d(:,ijhom+jl) = zt2ns(:,jl,2)
+            pt2d(:,jl      ) = t2sn(:,jl,2)
+            pt2d(:,ijhom+jl) = t2ns(:,jl,2)
          END DO
       CASE ( 1 ) 
          DO jl = 1, jprecj
-            pt2d(:,jl      ) = zt2sn(:,jl,2)
+            pt2d(:,jl      ) = t2sn(:,jl,2)
          END DO
       END SELECT
@@ -694,5 +694,5 @@
       ENDIF
       !
-      DEALLOCATE( zt2ns   , zt2sn   , zt2ew   , zt2we)
+      DEALLOCATE( t2ns   , t2sn   , t2ew   , t2we)
    END SUBROUTINE mpp_lnk_2d
 
@@ -729,11 +729,10 @@
       INTEGER  ::   ml_req1, ml_req2, ml_err   ! for key_mpi_isend
       INTEGER, DIMENSION(MPI_STATUS_SIZE) ::   ml_stat   ! for key_mpi_isend
-
-      REAL(wp), DIMENSION(:,:,:,:,:), ALLOCATABLE ::   zt4ns, zt4sn   ! 2 x 3d for north-south & south-north
-      REAL(wp), DIMENSION(:,:,:,:,:), ALLOCATABLE ::   zt4ew, zt4we   ! 2 x 3d for east-west & west-east
-
-      !!----------------------------------------------------------------------
-      ALLOCATE( zt4ns(jpi,jprecj,jpk,2,2), zt4sn(jpi,jprecj,jpk,2,2) ,    &
-         &      zt4ew(jpj,jpreci,jpk,2,2), zt4we(jpj,jpreci,jpk,2,2) )
+      REAL(wp), DIMENSION(:,:,:,:,:), ALLOCATABLE ::   t4ns, t4sn   ! 2 x 3d for north-south & south-north
+      REAL(wp), DIMENSION(:,:,:,:,:), ALLOCATABLE ::   t4ew, t4we   ! 2 x 3d for east-west & west-east
+
+      !!----------------------------------------------------------------------
+      ALLOCATE(t4ns(jpi,jprecj,jpk,2,2) , t4sn(jpi,jprecj,jpk,2,2) ,    &
+         &     t4ew(jpj,jpreci,jpk,2,2) , t4we(jpj,jpreci,jpk,2,2))
 
 
@@ -770,8 +769,8 @@
          iihom = nlci-nreci
          DO jl = 1, jpreci
-            zt4ew(:,jl,:,1,1) = ptab1(jpreci+jl,:,:)
-            zt4we(:,jl,:,1,1) = ptab1(iihom +jl,:,:)
-            zt4ew(:,jl,:,2,1) = ptab2(jpreci+jl,:,:)
-            zt4we(:,jl,:,2,1) = ptab2(iihom +jl,:,:)
+            t4ew(:,jl,:,1,1) = ptab1(jpreci+jl,:,:)
+            t4we(:,jl,:,1,1) = ptab1(iihom +jl,:,:)
+            t4ew(:,jl,:,2,1) = ptab2(jpreci+jl,:,:)
+            t4we(:,jl,:,2,1) = ptab2(iihom +jl,:,:)
          END DO
       END SELECT
@@ -782,17 +781,17 @@
       SELECT CASE ( nbondi ) 
       CASE ( -1 )
-         CALL mppsend( 2, zt4we(1,1,1,1,1), imigr, noea, ml_req1 )
-         CALL mpprecv( 1, zt4ew(1,1,1,1,2), imigr, noea )
+         CALL mppsend( 2, t4we(1,1,1,1,1), imigr, noea, ml_req1 )
+         CALL mpprecv( 1, t4ew(1,1,1,1,2), imigr, noea )
          IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
       CASE ( 0 )
-         CALL mppsend( 1, zt4ew(1,1,1,1,1), imigr, nowe, ml_req1 )
-         CALL mppsend( 2, zt4we(1,1,1,1,1), imigr, noea, ml_req2 )
-         CALL mpprecv( 1, zt4ew(1,1,1,1,2), imigr, noea )
-         CALL mpprecv( 2, zt4we(1,1,1,1,2), imigr, nowe )
+         CALL mppsend( 1, t4ew(1,1,1,1,1), imigr, nowe, ml_req1 )
+         CALL mppsend( 2, t4we(1,1,1,1,1), imigr, noea, ml_req2 )
+         CALL mpprecv( 1, t4ew(1,1,1,1,2), imigr, noea )
+         CALL mpprecv( 2, t4we(1,1,1,1,2), imigr, nowe )
          IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
          IF(l_isend) CALL mpi_wait(ml_req2, ml_stat, ml_err)
       CASE ( 1 )
-         CALL mppsend( 1, zt4ew(1,1,1,1,1), imigr, nowe, ml_req1 )
-         CALL mpprecv( 2, zt4we(1,1,1,1,2), imigr, nowe )
+         CALL mppsend( 1, t4ew(1,1,1,1,1), imigr, nowe, ml_req1 )
+         CALL mpprecv( 2, t4we(1,1,1,1,2), imigr, nowe )
          IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
       END SELECT
@@ -804,18 +803,18 @@
       CASE ( -1 )
          DO jl = 1, jpreci
-            ptab1(iihom+jl,:,:) = zt4ew(:,jl,:,1,2)
-            ptab2(iihom+jl,:,:) = zt4ew(:,jl,:,2,2)
+            ptab1(iihom+jl,:,:) = t4ew(:,jl,:,1,2)
+            ptab2(iihom+jl,:,:) = t4ew(:,jl,:,2,2)
          END DO
       CASE ( 0 ) 
          DO jl = 1, jpreci
-            ptab1(jl      ,:,:) = zt4we(:,jl,:,1,2)
-            ptab1(iihom+jl,:,:) = zt4ew(:,jl,:,1,2)
-            ptab2(jl      ,:,:) = zt4we(:,jl,:,2,2)
-            ptab2(iihom+jl,:,:) = zt4ew(:,jl,:,2,2)
+            ptab1(jl      ,:,:) = t4we(:,jl,:,1,2)
+            ptab1(iihom+jl,:,:) = t4ew(:,jl,:,1,2)
+            ptab2(jl      ,:,:) = t4we(:,jl,:,2,2)
+            ptab2(iihom+jl,:,:) = t4ew(:,jl,:,2,2)
          END DO
       CASE ( 1 )
          DO jl = 1, jpreci
-            ptab1(jl      ,:,:) = zt4we(:,jl,:,1,2)
-            ptab2(jl      ,:,:) = zt4we(:,jl,:,2,2)
+            ptab1(jl      ,:,:) = t4we(:,jl,:,1,2)
+            ptab2(jl      ,:,:) = t4we(:,jl,:,2,2)
          END DO
       END SELECT
@@ -829,8 +828,8 @@
          ijhom = nlcj - nrecj
          DO jl = 1, jprecj
-            zt4sn(:,jl,:,1,1) = ptab1(:,ijhom +jl,:)
-            zt4ns(:,jl,:,1,1) = ptab1(:,jprecj+jl,:)
-            zt4sn(:,jl,:,2,1) = ptab2(:,ijhom +jl,:)
-            zt4ns(:,jl,:,2,1) = ptab2(:,jprecj+jl,:)
+            t4sn(:,jl,:,1,1) = ptab1(:,ijhom +jl,:)
+            t4ns(:,jl,:,1,1) = ptab1(:,jprecj+jl,:)
+            t4sn(:,jl,:,2,1) = ptab2(:,ijhom +jl,:)
+            t4ns(:,jl,:,2,1) = ptab2(:,jprecj+jl,:)
          END DO
       ENDIF
@@ -841,17 +840,17 @@
       SELECT CASE ( nbondj )     
       CASE ( -1 )
-         CALL mppsend( 4, zt4sn(1,1,1,1,1), imigr, nono, ml_req1 )
-         CALL mpprecv( 3, zt4ns(1,1,1,1,2), imigr, nono )
+         CALL mppsend( 4, t4sn(1,1,1,1,1), imigr, nono, ml_req1 )
+         CALL mpprecv( 3, t4ns(1,1,1,1,2), imigr, nono )
          IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
       CASE ( 0 )
-         CALL mppsend( 3, zt4ns(1,1,1,1,1), imigr, noso, ml_req1 )
-         CALL mppsend( 4, zt4sn(1,1,1,1,1), imigr, nono, ml_req2 )
-         CALL mpprecv( 3, zt4ns(1,1,1,1,2), imigr, nono )
-         CALL mpprecv( 4, zt4sn(1,1,1,1,2), imigr, noso )
+         CALL mppsend( 3, t4ns(1,1,1,1,1), imigr, noso, ml_req1 )
+         CALL mppsend( 4, t4sn(1,1,1,1,1), imigr, nono, ml_req2 )
+         CALL mpprecv( 3, t4ns(1,1,1,1,2), imigr, nono )
+         CALL mpprecv( 4, t4sn(1,1,1,1,2), imigr, noso )
          IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
          IF(l_isend) CALL mpi_wait(ml_req2, ml_stat, ml_err)
       CASE ( 1 ) 
-         CALL mppsend( 3, zt4ns(1,1,1,1,1), imigr, noso, ml_req1 )
-         CALL mpprecv( 4, zt4sn(1,1,1,1,2), imigr, noso )
+         CALL mppsend( 3, t4ns(1,1,1,1,1), imigr, noso, ml_req1 )
+         CALL mpprecv( 4, t4sn(1,1,1,1,2), imigr, noso )
          IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
       END SELECT
@@ -863,18 +862,18 @@
       CASE ( -1 )
          DO jl = 1, jprecj
-            ptab1(:,ijhom+jl,:) = zt4ns(:,jl,:,1,2)
-            ptab2(:,ijhom+jl,:) = zt4ns(:,jl,:,2,2)
+            ptab1(:,ijhom+jl,:) = t4ns(:,jl,:,1,2)
+            ptab2(:,ijhom+jl,:) = t4ns(:,jl,:,2,2)
          END DO
       CASE ( 0 ) 
          DO jl = 1, jprecj
-            ptab1(:,jl      ,:) = zt4sn(:,jl,:,1,2)
-            ptab1(:,ijhom+jl,:) = zt4ns(:,jl,:,1,2)
-            ptab2(:,jl      ,:) = zt4sn(:,jl,:,2,2)
-            ptab2(:,ijhom+jl,:) = zt4ns(:,jl,:,2,2)
+            ptab1(:,jl      ,:) = t4sn(:,jl,:,1,2)
+            ptab1(:,ijhom+jl,:) = t4ns(:,jl,:,1,2)
+            ptab2(:,jl      ,:) = t4sn(:,jl,:,2,2)
+            ptab2(:,ijhom+jl,:) = t4ns(:,jl,:,2,2)
          END DO
       CASE ( 1 )
          DO jl = 1, jprecj
-            ptab1(:,jl,:) = zt4sn(:,jl,:,1,2)
-            ptab2(:,jl,:) = zt4sn(:,jl,:,2,2)
+            ptab1(:,jl,:) = t4sn(:,jl,:,1,2)
+            ptab2(:,jl,:) = t4sn(:,jl,:,2,2)
          END DO
       END SELECT
@@ -896,6 +895,5 @@
       ENDIF
       !
-      DEALLOCATE(zt4ns , zt4sn , zt4ew , zt4we)
-      !
+      DEALLOCATE(t4ns , t4sn , t4ew , t4we)
    END SUBROUTINE mpp_lnk_3d_gather
 
@@ -932,12 +930,12 @@
       INTEGER  ::   ml_req1, ml_req2, ml_err   ! for key_mpi_isend
       INTEGER, DIMENSION(MPI_STATUS_SIZE) ::   ml_stat   ! for key_mpi_isend
-      REAL(wp), DIMENSION(:,:,:)    , ALLOCATABLE ::   ztr2ns, ztr2sn ! 2d for north-south & south-north + extra outer halo
-      REAL(wp), DIMENSION(:,:,:)    , ALLOCATABLE ::   ztr2ew, ztr2we ! 2d for east-west   & west-east   + extra outer halo
-
-      !!----------------------------------------------------------------------
-      ALLOCATE( ztr2ns(1-jpr2di:jpi+jpr2di,jprecj+jpr2dj,2) ,                                                     &
-         &      ztr2sn(1-jpr2di:jpi+jpr2di,jprecj+jpr2dj,2) ,                                                     &
-         &      ztr2ew(1-jpr2dj:jpj+jpr2dj,jpreci+jpr2di,2) ,                                                     &
-         &      ztr2we(1-jpr2dj:jpj+jpr2dj,jpreci+jpr2di,2))
+      REAL(wp), DIMENSION(:,:,:)    , ALLOCATABLE ::   tr2ns, tr2sn ! 2d for north-south & south-north + extra outer halo
+      REAL(wp), DIMENSION(:,:,:)    , ALLOCATABLE ::   tr2ew, tr2we ! 2d for east-west   & west-east   + extra outer halo
+
+      !!----------------------------------------------------------------------
+      ALLOCATE( tr2ns(1-jpr2di:jpi+jpr2di,jprecj+jpr2dj,2) ,                                                     &
+         &      tr2sn(1-jpr2di:jpi+jpr2di,jprecj+jpr2dj,2) ,                                                     &
+         &      tr2ew(1-jpr2dj:jpj+jpr2dj,jpreci+jpr2di,2) ,                                                     &
+         &      tr2we(1-jpr2dj:jpj+jpr2dj,jpreci+jpr2di,2))
          
       ipreci = jpreci + jpr2di      ! take into account outer extra 2D overlap area
@@ -984,6 +982,6 @@
          iihom = nlci-nreci-jpr2di
          DO jl = 1, ipreci
-            ztr2ew(:,jl,1) = pt2d(jpreci+jl,:)
-            ztr2we(:,jl,1) = pt2d(iihom +jl,:)
+            tr2ew(:,jl,1) = pt2d(jpreci+jl,:)
+            tr2we(:,jl,1) = pt2d(iihom +jl,:)
          END DO
       END SELECT
@@ -994,17 +992,17 @@
       SELECT CASE ( nbondi )
       CASE ( -1 )
-         CALL mppsend( 2, ztr2we(1-jpr2dj,1,1), imigr, noea, ml_req1 )
-         CALL mpprecv( 1, ztr2ew(1-jpr2dj,1,2), imigr, noea )
+         CALL mppsend( 2, tr2we(1-jpr2dj,1,1), imigr, noea, ml_req1 )
+         CALL mpprecv( 1, tr2ew(1-jpr2dj,1,2), imigr, noea )
          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
       CASE ( 0 )
-         CALL mppsend( 1, ztr2ew(1-jpr2dj,1,1), imigr, nowe, ml_req1 )
-         CALL mppsend( 2, ztr2we(1-jpr2dj,1,1), imigr, noea, ml_req2 )
-         CALL mpprecv( 1, ztr2ew(1-jpr2dj,1,2), imigr, noea )
-         CALL mpprecv( 2, ztr2we(1-jpr2dj,1,2), imigr, nowe )
+         CALL mppsend( 1, tr2ew(1-jpr2dj,1,1), imigr, nowe, ml_req1 )
+         CALL mppsend( 2, tr2we(1-jpr2dj,1,1), imigr, noea, ml_req2 )
+         CALL mpprecv( 1, tr2ew(1-jpr2dj,1,2), imigr, noea )
+         CALL mpprecv( 2, tr2we(1-jpr2dj,1,2), imigr, nowe )
          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
          IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err)
       CASE ( 1 )
-         CALL mppsend( 1, ztr2ew(1-jpr2dj,1,1), imigr, nowe, ml_req1 )
-         CALL mpprecv( 2, ztr2we(1-jpr2dj,1,2), imigr, nowe )
+         CALL mppsend( 1, tr2ew(1-jpr2dj,1,1), imigr, nowe, ml_req1 )
+         CALL mpprecv( 2, tr2we(1-jpr2dj,1,2), imigr, nowe )
          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
       END SELECT
@@ -1016,14 +1014,14 @@
       CASE ( -1 )
          DO jl = 1, ipreci
-            pt2d(iihom+jl,:) = ztr2ew(:,jl,2)
+            pt2d(iihom+jl,:) = tr2ew(:,jl,2)
          END DO
       CASE ( 0 )
          DO jl = 1, ipreci
-            pt2d(jl-jpr2di,:) = ztr2we(:,jl,2)
-            pt2d( iihom+jl,:) = ztr2ew(:,jl,2)
+            pt2d(jl-jpr2di,:) = tr2we(:,jl,2)
+            pt2d( iihom+jl,:) = tr2ew(:,jl,2)
          END DO
       CASE ( 1 )
          DO jl = 1, ipreci
-            pt2d(jl-jpr2di,:) = ztr2we(:,jl,2)
+            pt2d(jl-jpr2di,:) = tr2we(:,jl,2)
          END DO
       END SELECT
@@ -1037,6 +1035,6 @@
          ijhom = nlcj-nrecj-jpr2dj
          DO jl = 1, iprecj
-            ztr2sn(:,jl,1) = pt2d(:,ijhom +jl)
-            ztr2ns(:,jl,1) = pt2d(:,jprecj+jl)
+            tr2sn(:,jl,1) = pt2d(:,ijhom +jl)
+            tr2ns(:,jl,1) = pt2d(:,jprecj+jl)
          END DO
       ENDIF
@@ -1047,17 +1045,17 @@
       SELECT CASE ( nbondj )
       CASE ( -1 )
-         CALL mppsend( 4, ztr2sn(1-jpr2di,1,1), imigr, nono, ml_req1 )
-         CALL mpprecv( 3, ztr2ns(1-jpr2di,1,2), imigr, nono )
+         CALL mppsend( 4, tr2sn(1-jpr2di,1,1), imigr, nono, ml_req1 )
+         CALL mpprecv( 3, tr2ns(1-jpr2di,1,2), imigr, nono )
          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
       CASE ( 0 )
-         CALL mppsend( 3, ztr2ns(1-jpr2di,1,1), imigr, noso, ml_req1 )
-         CALL mppsend( 4, ztr2sn(1-jpr2di,1,1), imigr, nono, ml_req2 )
-         CALL mpprecv( 3, ztr2ns(1-jpr2di,1,2), imigr, nono )
-         CALL mpprecv( 4, ztr2sn(1-jpr2di,1,2), imigr, noso )
+         CALL mppsend( 3, tr2ns(1-jpr2di,1,1), imigr, noso, ml_req1 )
+         CALL mppsend( 4, tr2sn(1-jpr2di,1,1), imigr, nono, ml_req2 )
+         CALL mpprecv( 3, tr2ns(1-jpr2di,1,2), imigr, nono )
+         CALL mpprecv( 4, tr2sn(1-jpr2di,1,2), imigr, noso )
          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
          IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err)
       CASE ( 1 )
-         CALL mppsend( 3, ztr2ns(1-jpr2di,1,1), imigr, noso, ml_req1 )
-         CALL mpprecv( 4, ztr2sn(1-jpr2di,1,2), imigr, noso )
+         CALL mppsend( 3, tr2ns(1-jpr2di,1,1), imigr, noso, ml_req1 )
+         CALL mpprecv( 4, tr2sn(1-jpr2di,1,2), imigr, noso )
          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
       END SELECT
@@ -1069,17 +1067,17 @@
       CASE ( -1 )
          DO jl = 1, iprecj
-            pt2d(:,ijhom+jl) = ztr2ns(:,jl,2)
+            pt2d(:,ijhom+jl) = tr2ns(:,jl,2)
          END DO
       CASE ( 0 )
          DO jl = 1, iprecj
-            pt2d(:,jl-jpr2dj) = ztr2sn(:,jl,2)
-            pt2d(:,ijhom+jl ) = ztr2ns(:,jl,2)
+            pt2d(:,jl-jpr2dj) = tr2sn(:,jl,2)
+            pt2d(:,ijhom+jl ) = tr2ns(:,jl,2)
          END DO
       CASE ( 1 ) 
          DO jl = 1, iprecj
-            pt2d(:,jl-jpr2dj) = ztr2sn(:,jl,2)
+            pt2d(:,jl-jpr2dj) = tr2sn(:,jl,2)
          END DO
       END SELECT
-      DEALLOCATE( ztr2ns , ztr2sn , ztr2ew , ztr2we)
+      DEALLOCATE( tr2ns , tr2sn , tr2ew , tr2we)
    END SUBROUTINE mpp_lnk_2d_e
 
@@ -1751,10 +1749,10 @@
       INTEGER ::   ml_stat(MPI_STATUS_SIZE)    ! for key_mpi_isend
       REAL(wp), POINTER, DIMENSION(:,:) ::   ztab   ! temporary workspace
-      REAL(wp), DIMENSION(:,:,:)    , ALLOCATABLE ::   zt2ns, zt2sn   ! 2d for north-south & south-north
-      REAL(wp), DIMENSION(:,:,:)    , ALLOCATABLE ::   zt2ew, zt2we   ! 2d for east-west & west-east
-
-      !!----------------------------------------------------------------------
-      ALLOCATE( zt2ns(jpi,jprecj    ,2)   , zt2sn(jpi,jprecj    ,2)   ,   &
-         &      zt2ew(jpj,jpreci    ,2)   , zt2we(jpj,jpreci    ,2))
+      REAL(wp), DIMENSION(:,:,:)    , ALLOCATABLE ::   t2ns, t2sn   ! 2d for north-south & south-north
+      REAL(wp), DIMENSION(:,:,:)    , ALLOCATABLE ::   t2ew, t2we   ! 2d for east-west & west-east
+
+      !!----------------------------------------------------------------------
+      ALLOCATE( t2ns(jpi,jprecj    ,2)   , t2sn(jpi,jprecj    ,2)   ,   &
+         &      t2ew(jpj,jpreci    ,2)   , t2we(jpj,jpreci    ,2))
 
 
@@ -1811,6 +1809,6 @@
             iihom = nlci-nreci
             DO jl = 1, jpreci
-               zt2ew(:,jl,1) = ztab(jpreci+jl,:)
-               zt2we(:,jl,1) = ztab(iihom +jl,:)
+               t2ew(:,jl,1) = ztab(jpreci+jl,:)
+               t2we(:,jl,1) = ztab(iihom +jl,:)
             END DO
          ENDIF
@@ -1820,17 +1818,17 @@
          !
          IF( nbondi == -1 ) THEN
-            CALL mppsend( 2, zt2we(1,1,1), imigr, noea, ml_req1 )
-            CALL mpprecv( 1, zt2ew(1,1,2), imigr, noea )
+            CALL mppsend( 2, t2we(1,1,1), imigr, noea, ml_req1 )
+            CALL mpprecv( 1, t2ew(1,1,2), imigr, noea )
             IF(l_isend)   CALL mpi_wait( ml_req1, ml_stat, ml_err )
          ELSEIF( nbondi == 0 ) THEN
-            CALL mppsend( 1, zt2ew(1,1,1), imigr, nowe, ml_req1 )
-            CALL mppsend( 2, zt2we(1,1,1), imigr, noea, ml_req2 )
-            CALL mpprecv( 1, zt2ew(1,1,2), imigr, noea )
-            CALL mpprecv( 2, zt2we(1,1,2), imigr, nowe )
+            CALL mppsend( 1, t2ew(1,1,1), imigr, nowe, ml_req1 )
+            CALL mppsend( 2, t2we(1,1,1), imigr, noea, ml_req2 )
+            CALL mpprecv( 1, t2ew(1,1,2), imigr, noea )
+            CALL mpprecv( 2, t2we(1,1,2), imigr, nowe )
             IF(l_isend)   CALL mpi_wait( ml_req1, ml_stat, ml_err )
             IF(l_isend)   CALL mpi_wait( ml_req2, ml_stat, ml_err )
          ELSEIF( nbondi == 1 ) THEN
-            CALL mppsend( 1, zt2ew(1,1,1), imigr, nowe, ml_req1 )
-            CALL mpprecv( 2, zt2we(1,1,2), imigr, nowe )
+            CALL mppsend( 1, t2ew(1,1,1), imigr, nowe, ml_req1 )
+            CALL mpprecv( 2, t2we(1,1,2), imigr, nowe )
             IF(l_isend) CALL mpi_wait( ml_req1, ml_stat, ml_err )
          ENDIF
@@ -1841,10 +1839,10 @@
          IF( nbondi == 0 .OR. nbondi == 1 ) THEN
             DO jl = 1, jpreci
-               ztab(jl,:) = zt2we(:,jl,2)
+               ztab(jl,:) = t2we(:,jl,2)
             END DO
          ENDIF
          IF( nbondi == -1 .OR. nbondi == 0 ) THEN
             DO jl = 1, jpreci
-               ztab(iihom+jl,:) = zt2ew(:,jl,2)
+               ztab(iihom+jl,:) = t2ew(:,jl,2)
             END DO
          ENDIF
@@ -1857,6 +1855,6 @@
             ijhom = nlcj-nrecj
             DO jl = 1, jprecj
-               zt2sn(:,jl,1) = ztab(:,ijhom +jl)
-               zt2ns(:,jl,1) = ztab(:,jprecj+jl)
+               t2sn(:,jl,1) = ztab(:,ijhom +jl)
+               t2ns(:,jl,1) = ztab(:,jprecj+jl)
             END DO
          ENDIF
@@ -1866,17 +1864,17 @@
          !
          IF( nbondj == -1 ) THEN
-            CALL mppsend( 4, zt2sn(1,1,1), imigr, nono, ml_req1 )
-            CALL mpprecv( 3, zt2ns(1,1,2), imigr, nono )
+            CALL mppsend( 4, t2sn(1,1,1), imigr, nono, ml_req1 )
+            CALL mpprecv( 3, t2ns(1,1,2), imigr, nono )
             IF(l_isend) CALL mpi_wait( ml_req1, ml_stat, ml_err )
          ELSEIF( nbondj == 0 ) THEN
-            CALL mppsend( 3, zt2ns(1,1,1), imigr, noso, ml_req1 )
-            CALL mppsend( 4, zt2sn(1,1,1), imigr, nono, ml_req2 )
-            CALL mpprecv( 3, zt2ns(1,1,2), imigr, nono )
-            CALL mpprecv( 4, zt2sn(1,1,2), imigr, noso )
+            CALL mppsend( 3, t2ns(1,1,1), imigr, noso, ml_req1 )
+            CALL mppsend( 4, t2sn(1,1,1), imigr, nono, ml_req2 )
+            CALL mpprecv( 3, t2ns(1,1,2), imigr, nono )
+            CALL mpprecv( 4, t2sn(1,1,2), imigr, noso )
             IF( l_isend )   CALL mpi_wait( ml_req1, ml_stat, ml_err )
             IF( l_isend )   CALL mpi_wait( ml_req2, ml_stat, ml_err )
          ELSEIF( nbondj == 1 ) THEN
-            CALL mppsend( 3, zt2ns(1,1,1), imigr, noso, ml_req1 )
-            CALL mpprecv( 4, zt2sn(1,1,2), imigr, noso)
+            CALL mppsend( 3, t2ns(1,1,1), imigr, noso, ml_req1 )
+            CALL mpprecv( 4, t2sn(1,1,2), imigr, noso)
             IF( l_isend )   CALL mpi_wait( ml_req1, ml_stat, ml_err )
          ENDIF
@@ -1886,10 +1884,10 @@
          IF( nbondj == 0 .OR. nbondj == 1 ) THEN
             DO jl = 1, jprecj
-               ztab(:,jl) = zt2sn(:,jl,2)
+               ztab(:,jl) = t2sn(:,jl,2)
             END DO
          ENDIF
          IF( nbondj == 0 .OR. nbondj == -1 ) THEN
             DO jl = 1, jprecj
-               ztab(:,ijhom+jl) = zt2ns(:,jl,2)
+               ztab(:,ijhom+jl) = t2ns(:,jl,2)
             END DO
          ENDIF
@@ -1910,6 +1908,5 @@
       END DO
       !
-      DEALLOCATE( zt2ns   , zt2sn   , zt2ew   , zt2we)
-
+      DEALLOCATE( t2ns   , t2sn   , t2ew   , t2we)
       CALL wrk_dealloc( jpi,jpj, ztab )
       !
Index: branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/CONFIG/ORCA2_LIM_CRS/MY_SRC/mppini.F90
===================================================================
--- branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/CONFIG/ORCA2_LIM_CRS/MY_SRC/mppini.F90	(revision 3895)
+++ branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/CONFIG/ORCA2_LIM_CRS/MY_SRC/mppini.F90	(revision 3895)
@@ -0,0 +1,516 @@
+
+MODULE mppini
+   !!==============================================================================
+   !!                       ***  MODULE mppini   ***
+   !! Ocean initialization : distributed memory computing initialization
+   !!==============================================================================
+
+   !!----------------------------------------------------------------------
+   !!   mpp_init       : Lay out the global domain over processors
+   !!   mpp_init2      : Lay out the global domain over processors 
+   !!                    with land processor elimination
+   !!   mpp_init_ioispl: IOIPSL initialization in mpp
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE dom_oce         ! ocean space and time domain 
+   USE in_out_manager  ! I/O Manager
+   USE lib_mpp         ! distribued memory computing library
+   USE ioipsl
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC mpp_init       ! called by opa.F90
+   PUBLIC mpp_init2      ! called by opa.F90
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: mppini.F90 3294 2012-01-28 16:44:18Z rblod $ 
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+#if ! defined key_mpp_mpi
+   !!----------------------------------------------------------------------
+   !!   Default option :                            shared memory computing
+   !!----------------------------------------------------------------------
+
+   SUBROUTINE mpp_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE mpp_init  ***
+      !!
+      !! ** Purpose :   Lay out the global domain over processors.
+      !!
+      !! ** Method  :   Shared memory computing, set the local processor
+      !!      variables to the value of the global domain
+      !!
+      !! History :
+      !!   9.0  !  04-01  (G. Madec, J.M. Molines)  F90 : free form, north fold jpni >1
+      !!----------------------------------------------------------------------
+
+      ! No mpp computation
+      nimpp  = 1
+      njmpp  = 1
+      nlci   = jpi
+      nlcj   = jpj
+      nldi   = 1
+      nldj   = 1
+      nlei   = jpi
+      nlej   = jpj
+      nperio = jperio
+      nbondi = 2
+      nbondj = 2
+      nidom  = FLIO_DOM_NONE
+      npolj = jperio
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'mpp_init(2) : NO massively parallel processing'
+         WRITE(numout,*) '~~~~~~~~~~~: '
+         WRITE(numout,*) '         nperio = ', nperio
+         WRITE(numout,*) '         npolj  = ', npolj
+         WRITE(numout,*) '         nimpp  = ', nimpp
+         WRITE(numout,*) '         njmpp  = ', njmpp
+      ENDIF
+
+      IF(  jpni /= 1 .OR. jpnj /= 1 .OR. jpnij /= 1 ) &
+          CALL ctl_stop( 'equality  jpni = jpnj = jpnij = 1 is not satisfied',   &
+          &              'the domain is lay out for distributed memory computing! ' )
+
+   END SUBROUTINE mpp_init
+
+
+   SUBROUTINE mpp_init2 
+      CALL mpp_init                             ! same routine as mpp_init
+   END SUBROUTINE mpp_init2
+
+#else
+   !!----------------------------------------------------------------------
+   !!   'key_mpp_mpi'          OR         MPI massively parallel processing
+   !!----------------------------------------------------------------------
+
+   SUBROUTINE mpp_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE mpp_init  ***
+      !!                    
+      !! ** Purpose :   Lay out the global domain over processors.
+      !!
+      !! ** Method  :   Global domain is distributed in smaller local domains.
+      !!      Periodic condition is a function of the local domain position
+      !!      (global boundary or neighbouring domain) and of the global
+      !!      periodic
+      !!      Type :         jperio global periodic condition
+      !!                     nperio local  periodic condition
+      !!
+      !! ** Action  : - set domain parameters
+      !!                    nimpp     : longitudinal index 
+      !!                    njmpp     : latitudinal  index
+      !!                    nperio    : lateral condition type 
+      !!                    narea     : number for local area
+      !!                    nlci      : first dimension
+      !!                    nlcj      : second dimension
+      !!                    nbondi    : mark for "east-west local boundary"
+      !!                    nbondj    : mark for "north-south local boundary"
+      !!                    nproc     : number for local processor
+      !!                    noea      : number for local neighboring processor
+      !!                    nowe      : number for local neighboring processor
+      !!                    noso      : number for local neighboring processor
+      !!                    nono      : number for local neighboring processor
+      !!
+      !! History :
+      !!        !  94-11  (M. Guyon)  Original code
+      !!        !  95-04  (J. Escobar, M. Imbard)
+      !!        !  98-02  (M. Guyon)  FETI method
+      !!        !  98-05  (M. Imbard, J. Escobar, L. Colombet )  SHMEM and MPI versions
+      !!   8.5  !  02-08  (G. Madec)  F90 : free form
+      !!   3.4  !  11-11  (C. Harris) decomposition changes for running with CICE
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jn   ! dummy loop indices
+      INTEGER  ::   ii, ij, ifreq, il1, il2            ! local integers
+      INTEGER  ::   iresti, irestj, ijm1, imil, inum   !   -      -
+      REAL(wp) ::   zidom, zjdom                       ! local scalars
+      INTEGER, DIMENSION(jpni,jpnj) ::   iimppt, ijmppt, ilcit, ilcjt   ! local workspace
+      !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'mpp_init : Message Passing MPI'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~'
+
+
+      !  1. Dimension arrays for subdomains
+      ! -----------------------------------
+      !  Computation of local domain sizes ilcit() ilcjt()
+      !  These dimensions depend on global sizes jpni,jpnj and jpiglo,jpjglo
+      !  The subdomains are squares leeser than or equal to the global
+      !  dimensions divided by the number of processors minus the overlap
+      !  array (cf. par_oce.F90).
+      
+      nreci  = 2 * jpreci
+      nrecj  = 2 * jprecj
+      iresti = MOD( jpiglo - nreci , jpni )
+      irestj = MOD( jpjglo - nrecj , jpnj )
+
+      IF(  iresti == 0 )   iresti = jpni
+
+#if defined key_nemocice_decomp
+      ! In order to match CICE the size of domains in NEMO has to be changed
+      ! The last line of blocks (west) will have fewer points
+
+      DO jj = 1, jpnj
+         DO ji=1, jpni-1
+            ilcit(ji,jj) = jpi
+         END DO
+         ilcit(jpni,jj) = jpiglo - (jpni - 1) * (jpi - nreci)
+      END DO
+
+#else
+
+      DO jj = 1, jpnj
+         DO ji = 1, iresti
+            ilcit(ji,jj) = jpi
+         END DO
+         DO ji = iresti+1, jpni
+            ilcit(ji,jj) = jpi -1
+         END DO
+      END DO
+      
+#endif
+      IF( irestj == 0 )   irestj = jpnj
+
+#if defined key_nemocice_decomp
+      ! Same change to domains in North-South direction as in East-West. 
+      DO ji=1,jpni
+         DO jj=1,jpnj-1
+            ilcjt(ji,jj) = jpj
+         END DO
+         ilcjt(ji,jpnj) = jpjglo - (jpnj - 1) * (jpj - nrecj)
+      END DO
+
+#else
+
+      DO ji = 1, jpni
+         DO jj = 1, irestj
+            ilcjt(ji,jj) = jpj
+         END DO
+         DO jj = irestj+1, jpnj
+            ilcjt(ji,jj) = jpj -1
+         END DO
+      END DO
+      
+#endif
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) '           defines mpp subdomains'
+         WRITE(numout,*) '           ----------------------'
+         WRITE(numout,*) '           iresti=',iresti,' irestj=',irestj
+         WRITE(numout,*) '           jpni  =',jpni  ,' jpnj  =',jpnj
+         ifreq = 4
+         il1   = 1
+         DO jn = 1, (jpni-1)/ifreq+1
+            il2 = MIN( jpni, il1+ifreq-1 )
+            WRITE(numout,*)
+            WRITE(numout,9200) ('***',ji = il1,il2-1)
+            DO jj = jpnj, 1, -1
+               WRITE(numout,9203) ('   ',ji = il1,il2-1)
+               WRITE(numout,9202) jj, ( ilcit(ji,jj),ilcjt(ji,jj),ji = il1,il2 )
+               WRITE(numout,9203) ('   ',ji = il1,il2-1)
+               WRITE(numout,9200) ('***',ji = il1,il2-1)
+            END DO
+            WRITE(numout,9201) (ji,ji = il1,il2)
+            il1 = il1+ifreq
+         END DO
+ 9200    FORMAT('     ***',20('*************',a3))
+ 9203    FORMAT('     *     ',20('         *   ',a3))
+ 9201    FORMAT('        ',20('   ',i3,'          '))
+ 9202    FORMAT(' ',i3,' *  ',20(i3,'  x',i3,'   *   '))
+      ENDIF
+
+      zidom = nreci
+      DO ji = 1, jpni
+         zidom = zidom + ilcit(ji,1) - nreci
+      END DO
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*)' sum ilcit(i,1) = ', zidom, ' jpiglo = ', jpiglo
+      
+      zjdom = nrecj
+      DO jj = 1, jpnj
+         zjdom = zjdom + ilcjt(1,jj) - nrecj
+      END DO
+      IF(lwp) WRITE(numout,*)' sum ilcit(1,j) = ', zjdom, ' jpjglo = ', jpjglo
+      IF(lwp) WRITE(numout,*)
+      
+
+      !  2. Index arrays for subdomains
+      ! -------------------------------
+      
+      iimppt(:,:) = 1
+      ijmppt(:,:) = 1
+      
+      IF( jpni > 1 ) THEN
+         DO jj = 1, jpnj
+            DO ji = 2, jpni
+               iimppt(ji,jj) = iimppt(ji-1,jj) + ilcit(ji-1,jj) - nreci
+            END DO
+         END DO
+      ENDIF
+
+      IF( jpnj > 1 ) THEN
+         DO jj = 2, jpnj
+            DO ji = 1, jpni
+               ijmppt(ji,jj) = ijmppt(ji,jj-1)+ilcjt(ji,jj-1)-nrecj
+            END DO
+         END DO
+      ENDIF
+      
+      ! 3. Subdomain description
+      ! ------------------------
+
+      DO jn = 1, jpnij
+         ii = 1 + MOD( jn-1, jpni )
+         ij = 1 + (jn-1) / jpni
+         nimppt(jn) = iimppt(ii,ij)
+         njmppt(jn) = ijmppt(ii,ij)
+         nlcit (jn) = ilcit (ii,ij)     
+         nlci       = nlcit (jn)     
+         nlcjt (jn) = ilcjt (ii,ij)     
+         nlcj       = nlcjt (jn)
+         nbondj = -1                                   ! general case
+         IF( jn   >  jpni          )   nbondj = 0      ! first row of processor
+         IF( jn   >  (jpnj-1)*jpni )   nbondj = 1      ! last  row of processor
+         IF( jpnj == 1             )   nbondj = 2      ! one processor only in j-direction
+         ibonjt(jn) = nbondj
+         
+         nbondi = 0                                    ! 
+         IF( MOD( jn, jpni ) == 1 )   nbondi = -1      !
+         IF( MOD( jn, jpni ) == 0 )   nbondi =  1      !
+         IF( jpni            == 1 )   nbondi =  2      ! one processor only in i-direction
+         ibonit(jn) = nbondi
+         
+         nldi =  1   + jpreci
+         nlei = nlci - jpreci
+         IF( nbondi == -1 .OR. nbondi == 2 )   nldi = 1
+         IF( nbondi ==  1 .OR. nbondi == 2 )   nlei = nlci
+         nldj =  1   + jprecj
+         nlej = nlcj - jprecj
+         IF( nbondj == -1 .OR. nbondj == 2 )   nldj = 1
+         IF( nbondj ==  1 .OR. nbondj == 2 )   nlej = nlcj
+         nldit(jn) = nldi
+         nleit(jn) = nlei
+         nldjt(jn) = nldj
+         nlejt(jn) = nlej
+   !      if(lwp) then
+            write(numout,*) ' proc : ', jn, ' nlcjt/njmppt : ', nlcjt(jn),njmppt(jn)
+   !      endif 
+      END DO
+      
+
+      ! 4. From global to local
+      ! -----------------------
+
+      nperio = 0
+      IF( jperio == 2 .AND. nbondj == -1 )   nperio = 2
+
+
+      ! 5. Subdomain neighbours
+      ! ----------------------
+
+      nproc = narea - 1
+      noso  = nproc - jpni
+      nowe  = nproc - 1
+      noea  = nproc + 1
+      nono  = nproc + jpni
+      ! great neighbours
+      npnw = nono - 1
+      npne = nono + 1
+      npsw = noso - 1
+      npse = noso + 1
+      nbsw = 1
+      nbnw = 1
+      IF( MOD( nproc, jpni ) == 0 ) THEN
+         nbsw = 0
+         nbnw = 0
+      ENDIF
+      nbse = 1
+      nbne = 1
+      IF( MOD( nproc, jpni ) == jpni-1 ) THEN
+         nbse = 0
+         nbne = 0
+      ENDIF
+      IF(nproc < jpni) THEN
+         nbsw = 0
+         nbse = 0
+      ENDIF
+      IF( nproc >= (jpnj-1)*jpni ) THEN
+         nbnw = 0
+         nbne = 0
+      ENDIF
+      nlcj = nlcjt(narea)  
+      nlci = nlcit(narea)  
+      nldi = nldit(narea)
+      nlei = nleit(narea)
+      nldj = nldjt(narea)
+      nlej = nlejt(narea)
+      nbondi = ibonit(narea)
+      nbondj = ibonjt(narea)
+      nimpp  = nimppt(narea)  
+      njmpp  = njmppt(narea)  
+
+     ! Save processor layout in layout.dat file 
+       IF (lwp) THEN
+        CALL ctl_opn( inum, 'layout.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE., narea )
+        WRITE(inum,'(a)') '   jpnij     jpi     jpj     jpk  jpiglo  jpjglo'
+        WRITE(inum,'(6i8)') jpnij,jpi,jpj,jpk,jpiglo,jpjglo
+        WRITE(inum,'(a)') 'NAREA nlci nlcj nldi nldj nlei nlej nimpp njmpp'
+
+        DO  jn = 1, jpnij
+         WRITE(inum,'(9i5)') jn, nlcit(jn), nlcjt(jn), &
+                                      nldit(jn), nldjt(jn), &
+                                      nleit(jn), nlejt(jn), &
+                                      nimppt(jn), njmppt(jn)
+        END DO
+        CLOSE(inum)   
+      END IF
+
+
+      ! w a r n i n g  narea (zone) /= nproc (processors)!
+
+      IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
+         IF( jpni == 1 )THEN
+            nbondi = 2
+            nperio = 1
+         ELSE
+            nbondi = 0
+         ENDIF
+         IF( MOD( narea, jpni ) == 0 ) THEN
+            noea = nproc-(jpni-1)
+            npne = npne-jpni
+            npse = npse-jpni
+         ENDIF
+         IF( MOD( narea, jpni ) == 1 ) THEN
+            nowe = nproc+(jpni-1)
+            npnw = npnw+jpni
+            npsw = npsw+jpni
+         ENDIF
+         nbsw = 1
+         nbnw = 1
+         nbse = 1
+         nbne = 1
+         IF( nproc < jpni ) THEN
+            nbsw = 0
+            nbse = 0
+         ENDIF
+         IF( nproc >= (jpnj-1)*jpni ) THEN
+            nbnw = 0
+            nbne = 0
+         ENDIF
+      ENDIF
+      npolj = 0
+      IF( jperio == 3 .OR. jperio == 4 ) THEN
+         ijm1 = jpni*(jpnj-1)
+         imil = ijm1+(jpni+1)/2
+         IF( narea > ijm1 ) npolj = 3
+         IF( MOD(jpni,2) == 1 .AND. narea == imil ) npolj = 4
+         IF( npolj == 3 ) nono = jpni*jpnj-narea+ijm1
+      ENDIF
+      IF( jperio == 5 .OR. jperio == 6 ) THEN
+          ijm1 = jpni*(jpnj-1)
+          imil = ijm1+(jpni+1)/2
+          IF( narea > ijm1) npolj = 5
+          IF( MOD(jpni,2) == 1 .AND. narea == imil ) npolj = 6
+          IF( npolj == 5 ) nono = jpni*jpnj-narea+ijm1
+      ENDIF
+
+      ! Periodicity : no corner if nbondi = 2 and nperio != 1
+
+      IF(lwp) THEN
+         WRITE(numout,*) ' nproc  = ', nproc
+         WRITE(numout,*) ' nowe   = ', nowe  , ' noea   =  ', noea
+         WRITE(numout,*) ' nono   = ', nono  , ' noso   =  ', noso
+         WRITE(numout,*) ' nbondi = ', nbondi
+         WRITE(numout,*) ' nbondj = ', nbondj
+         WRITE(numout,*) ' npolj  = ', npolj
+         WRITE(numout,*) ' nperio = ', nperio
+         WRITE(numout,*) ' nlci   = ', nlci
+         WRITE(numout,*) ' nlcj   = ', nlcj
+         WRITE(numout,*) ' nimpp  = ', nimpp
+         WRITE(numout,*) ' njmpp  = ', njmpp
+         WRITE(numout,*) ' nbse   = ', nbse  , ' npse   = ', npse
+         WRITE(numout,*) ' nbsw   = ', nbsw  , ' npsw   = ', npsw
+         WRITE(numout,*) ' nbne   = ', nbne  , ' npne   = ', npne
+         WRITE(numout,*) ' nbnw   = ', nbnw  , ' npnw   = ', npnw
+      ENDIF
+
+      IF( nperio == 1 .AND. jpni /= 1 ) CALL ctl_stop( ' mpp_init: error on cyclicity' )
+
+      ! Prepare mpp north fold
+
+      IF (jperio >= 3 .AND. jperio <= 6 .AND. jpni > 1 ) THEN
+         CALL mpp_ini_north
+      END IF
+
+      ! Prepare NetCDF output file (if necessary)
+      CALL mpp_init_ioipsl
+
+   END SUBROUTINE mpp_init
+
+#  include "mppini_2.h90"
+
+# if defined key_dimgout
+   !!----------------------------------------------------------------------
+   !!   'key_dimgout'                  NO use of NetCDF files
+   !!----------------------------------------------------------------------
+   SUBROUTINE mpp_init_ioipsl       ! Dummy routine
+   END SUBROUTINE mpp_init_ioipsl  
+# else
+   SUBROUTINE mpp_init_ioipsl
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE mpp_init_ioipsl  ***
+      !!
+      !! ** Purpose :   
+      !!
+      !! ** Method  :   
+      !!
+      !! History :
+      !!   9.0  !  04-03  (G. Madec )  MPP-IOIPSL 
+      !!   " "  !  08-12  (A. Coward)  addition in case of jpni*jpnj < jpnij
+      !!----------------------------------------------------------------------
+      INTEGER, DIMENSION(2) ::   iglo, iloc, iabsf, iabsl, ihals, ihale, idid
+      !!----------------------------------------------------------------------
+
+      ! The domain is split only horizontally along i- or/and j- direction
+      ! So we need at the most only 1D arrays with 2 elements.
+      ! Set idompar values equivalent to the jpdom_local_noextra definition
+      ! used in IOM. This works even if jpnij .ne. jpni*jpnj.
+      iglo(1) = jpiglo
+      iglo(2) = jpjglo
+      iloc(1) = nlci
+      iloc(2) = nlcj
+      iabsf(1) = nimppt(narea)
+      iabsf(2) = njmppt(narea)
+      iabsl(:) = iabsf(:) + iloc(:) - 1
+      ihals(1) = nldi - 1
+      ihals(2) = nldj - 1
+      ihale(1) = nlci - nlei
+      ihale(2) = nlcj - nlej
+      idid(1) = 1
+      idid(2) = 2
+
+      IF(lwp) THEN
+          WRITE(numout,*)
+          WRITE(numout,*) 'mpp_init_ioipsl :   iloc  = ', iloc (1), iloc (2)
+          WRITE(numout,*) '~~~~~~~~~~~~~~~     iabsf = ', iabsf(1), iabsf(2)
+          WRITE(numout,*) '                    ihals = ', ihals(1), ihals(2)
+          WRITE(numout,*) '                    ihale = ', ihale(1), ihale(2)
+      ENDIF
+      !
+      CALL flio_dom_set ( jpnij, nproc, idid, iglo, iloc, iabsf, iabsl, ihals, ihale, 'BOX', nidom)
+      !
+   END SUBROUTINE mpp_init_ioipsl  
+
+# endif
+#endif
+
+   !!======================================================================
+END MODULE mppini
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 3864)
+++ branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crs.F90	(revision 3895)
@@ -28,6 +28,6 @@
                    jpj_full                     !: 2nd dimension of local parent grid domain
 
-      INTEGER  ::  nistart, njstart
-      INTEGER  ::  niend  , njend
+      INTEGER  ::  nistr , njstr
+      INTEGER  ::  niend , njend
 
       INTEGER  ::  jpi_crsm1, jpj_crsm1         !: loop indices      
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 3864)
+++ branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsdiawri.F90	(revision 3895)
@@ -4,18 +4,5 @@
    !! Ocean diagnostics :  write ocean output files
    !!=====================================================================
-   !! History :  OPA  ! 1991-03  (M.-A. Foujols)  Original code
-   !!            4.0  ! 1991-11  (G. Madec)
-   !!                 ! 1992-06  (M. Imbard)  correction restart file
-   !!                 ! 1992-07  (M. Imbard)  split into diawri and rstwri
-   !!                 ! 1993-03  (M. Imbard)  suppress writibm
-   !!                 ! 1998-01  (C. Levy)  NETCDF format using ioipsl INTERFACE
-   !!                 ! 1999-02  (E. Guilyardi)  name of netCDF files + variables
-   !!            8.2  ! 2000-06  (M. Imbard)  Original code (diabort.F)
-   !!   NEMO     1.0  ! 2002-06  (A.Bozec, E. Durand)  Original code (diainit.F)
-   !!             -   ! 2002-09  (G. Madec)  F90: Free form and module
-   !!             -   ! 2002-12  (G. Madec)  merge of diabort and diainit, F90
-   !!                 ! 2005-11  (V. Garnier) Surface pressure gradient organization
-   !!            3.2  ! 2008-11  (B. Lemaire) creation from old diawri
-   !!                 ! 2012-07  (J. Simeon, G. Madec, C. Ethe, C. Calone) Modified for coarsened output
+   !!   2012-07  (J. Simeon, C. Calone, G. Madec, C. Ethe)
    !!----------------------------------------------------------------------
 
@@ -131,6 +118,5 @@
       !  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=zt, p_cfield3d=zt_crs )
+      CALL crs_dom_ope( zt, 'VOL', 'T', tmask, zt_crs, p_e12=e1e2t, p_e3=zfse3t )
       tsn_crs(:,:,:,jp_tem) = zt_crs(:,:,:)
       CALL crs_iom_put( "toce_crs", pv_r3d=tsn_crs(:,:,:,jp_tem) )    ! temp
@@ -138,7 +124,6 @@
 
       !  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=zt, p_cfield3d=zt_crs )
+      zs(:,:,:) = tsn(:,:,:,jp_sal)  ;      zs_crs(:,:,:) = 0._wp
+      CALL crs_dom_ope( zs, 'VOL', 'T', tmask, zs_crs, p_e12=e1e2t, p_e3=zfse3t )
       tsn_crs(:,:,:,jp_sal) = zt_crs(:,:,:)
       CALL crs_iom_put( "soce_crs" , pv_r3d=tsn_crs(:,:,:,jp_sal) )    ! sal
@@ -146,6 +131,5 @@
 
       !  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=un_crs )
+      CALL crs_dom_ope( un, 'SUM', 'U', umask, un_crs, p_e12=e2u, p_e3=zfse3u, p_surf_crs=e2e3u_msk )
       CALL crs_iom_put( "uoce_crs"  , pv_r3d=un_crs  )   ! i-current 
       !
@@ -159,15 +143,12 @@
          END DO
       END DO
-      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_dom_ope( zt, 'SUM', 'U', umask, zt_crs, p_e12=e2u, p_e3=zfse3u, p_surf_crs=e2e3u_msk )
       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_dom_ope( zs, 'SUM', 'U', umask, zs_crs, p_e12=e2u, p_e3=zfse3u, p_surf_crs=e2e3u_msk )
       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=vn_crs )
+      CALL crs_dom_ope( vn, 'SUM', 'V', vmask, vn_crs, p_e12=e1v, p_e3=zfse3v, p_surf_crs=e1e3v_msk )
       CALL crs_iom_put( "voce_crs"  , pv_r3d=vn_crs  )   ! v-current 
       !                                                                                 
@@ -181,21 +162,14 @@
          END DO
       END DO
-      CALL crsfun( p_e1_e2=e1v, cd_type='V', psgn=-1.0, p_pmask=vmask, &
-         &         p_fse3=zfse3v, p_pfield=zt, p_cfield3d=zt_crs )
+      CALL crs_dom_ope( zt, 'SUM', 'V', vmask, zt_crs, p_e12=e1v, p_e3=zfse3v, p_surf_crs=e1e3v_msk )
       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=zs, p_cfield3d=zs_crs )
+      CALL crs_dom_ope( zs, 'SUM', 'V', vmask, zs_crs, p_e12=e1v, p_e3=zfse3v, p_surf_crs=e1e3v_msk )
       CALL crs_iom_put( "voces_crs" , pv_r3d=zs_crs )   ! vS
       
       !  Kinetic energy
-      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=rke, p_cfield3d=rke_crs )
-      rke_crs(:,:,:) = rke_crs(:,:,:) * tmask_crs(:,:,:)
+      CALL crs_dom_ope( rke, 'VOL', 'T', tmask, rke_crs, p_e12=e1e2t, p_e3=zfse3t )
       CALL crs_iom_put( "eken_crs", pv_r3d=rke_crs )
    
-      !
-
-
       !  Horizontal divergence ( following OPA_SRC/DYN/divcur.F90 ) 
       DO jk = 1, jpkm1
@@ -213,22 +187,18 @@
          ENDDO
       ENDDO
-      CALL crs_lbc_lnk( hdivn_crs,'T', 1.0     )
+      CALL crs_lbc_lnk( hdivn_crs, 'T', 1.0  )
       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_dom_ope( sshn, 'VOL', 'T', tmask, sshn_crs, p_e12=e1e2t, p_e3=zfse3t )
       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 
-      
-
+      CALL crs_dom_ope( wn, 'VOL', 'W', tmask, wn_crs, p_e12=e1e2t, p_e3=zfse3w )
+      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 )
Index: branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsdom.F90
===================================================================
--- branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsdom.F90	(revision 3864)
+++ branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsdom.F90	(revision 3895)
@@ -28,5 +28,5 @@
    !!               Climate Dynamics, 14:101-116.
    !!  History: 
-   !!       Original.   May 2012.  (J. Simeon, G. Madec, C. Ethe, C. Calone)
+   !!       Original.   May 2012.  (J. Simeon, C. Calone, G. Madec, C. Ethe)
    !!===================================================================
 
@@ -43,9 +43,10 @@
    PRIVATE 
 
-   PUBLIC crsfun, crsfun_wgt
-   PUBLIC crs_dom_e3_max, crs_dom_sfc, crs_dom_msk, crs_dom_def, crs_dom_hgr, crs_dom_coordinates, crs_dom_bat
-
-   INTERFACE crsfun
-      MODULE PROCEDURE crsfun_UV, crsfun_TW
+   PUBLIC crs_dom_ope
+   PUBLIC crs_dom_e3_max, crs_dom_sfc, crs_dom_msk, crs_dom_hgr, crs_dom_coordinates
+   PUBLIC crs_dom_facvol, crs_dom_def, crs_dom_bat
+
+   INTERFACE crs_dom_ope
+      MODULE PROCEDURE crs_dom_ope_3d, crs_dom_ope_2d
    END INTERFACE
 
@@ -58,7 +59,7 @@
    SUBROUTINE crs_dom_msk
       
-      INTEGER  ::  ji, jj, jk, ijpk                   ! dummy loop indices
-      INTEGER  :: ijie,ijis,ijje,ijjs
-      REAL(wp) ::   zmask
+      INTEGER  ::  ji, jj, jk                   ! dummy loop indices
+      INTEGER  ::  ijie,ijis,ijje,ijjs
+      REAL(wp) ::  zmask
       
       ! Initialize
@@ -69,23 +70,23 @@
       fmask_crs(:,:,:) = 0.0
       
-      DO jk = 1, jpk
+      DO jk = 1, jpkm1
          DO ji = 2, nlei_crs  
             ijie = mie_crs(ji)
             ijis = mis_crs(ji)
-            DO jj = njstart, njend 
+            DO jj = nldj_crs, nlej_crs 
                ijje = mje_crs(jj) 
                ijjs = mjs_crs(jj) 
                                  
-               zmask = 0
+               zmask = 0.0
                zmask = SUM( tmask(ijis:ijie,ijjs:ijje,jk) ) 
-               IF ( zmask > 0 ) tmask_crs(ji,jj,jk) = 1
+               IF ( zmask > 0.0 ) tmask_crs(ji,jj,jk) = 1.0
                
-               zmask = 0
+               zmask = 0.0
                zmask = SUM( vmask(ijis:ijie,ijje     ,jk) )  
-               IF ( zmask > 0 ) vmask_crs(ji,jj,jk) = 1
+               IF ( zmask > 0.0 ) vmask_crs(ji,jj,jk) = 1.0
                
-               zmask = 0
+               zmask = 0.0
                zmask = SUM( umask(ijie     ,ijjs:ijje,jk) )  
-               IF ( zmask > 0 ) umask_crs(ji,jj,jk) = 1
+               IF ( zmask > 0.0 ) umask_crs(ji,jj,jk) = 1.0
                
                fmask_crs(ji,jj,jk) = fmask(ijie,ijje,jk)  
@@ -95,12 +96,12 @@
       ENDDO
       !
-      CALL crs_lbc_lnk( tmask_crs(:,:,:), 'T', psgn = 1.0 )
-      CALL crs_lbc_lnk( vmask_crs(:,:,:), 'V', psgn = 1.0 )
-      CALL crs_lbc_lnk( umask_crs(:,:,:), 'U', psgn = 1.0 )
-      CALL crs_lbc_lnk( fmask_crs(:,:,:), 'F', psgn = 1.0 )
+      CALL crs_lbc_lnk( tmask_crs, 'T', 1.0 )
+      CALL crs_lbc_lnk( vmask_crs, 'V', 1.0 )
+      CALL crs_lbc_lnk( umask_crs, 'U', 1.0 )
+      CALL crs_lbc_lnk( fmask_crs, 'F', 1.0 )
       !
    END SUBROUTINE crs_dom_msk
 
-   SUBROUTINE crs_dom_coordinates( p_pgphi, p_pglam, cd_type, p_cgphi, p_cglam )
+   SUBROUTINE crs_dom_coordinates( p_gphi, p_glam, cd_type, p_gphi_crs, p_glam_crs )
       !!----------------------------------------------------------------
       !!               *** SUBROUTINE crs_coordinates ***
@@ -116,70 +117,90 @@
       !!               the central f-corner.
       !!
-      !! ** Input   :  p_pgphi = parent grid gphi[t|u|v|f]
-      !!               p_pglam = parent grid glam[t|u|v|f]
+      !! ** Input   :  p_gphi = parent grid gphi[t|u|v|f]
+      !!               p_glam = parent grid glam[t|u|v|f]
       !!               cd_type  = grid type (T,U,V,F)
-      !! ** Output  :  p_cgphi = coarse grid gphi[t|u|v|f]
-      !!               p_cglam = coarse grid glam[t|u|v|f]
+      !! ** Output  :  p_gphi_crs = coarse grid gphi[t|u|v|f]
+      !!               p_glam_crs = coarse grid glam[t|u|v|f]
       !!              
       !! History. 1 Jun.
       !!----------------------------------------------------------------
       !! Arguments
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: p_pgphi  ! Parent grid latitude
-      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), 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
+      REAL(wp), DIMENSION(jpi,jpj)        , INTENT(in)  :: p_gphi  ! Parent grid latitude
+      REAL(wp), DIMENSION(jpi,jpj)        , INTENT(in)  :: p_glam  ! Parent grid longitude
+      CHARACTER(len=1),                     INTENT(in)  :: cd_type   ! grid type (T,U,V,F) 
+      REAL(wp), DIMENSION(jpi_crs,jpj_crs), INTENT(out) :: p_gphi_crs  ! Coarse grid latitude
+      REAL(wp), DIMENSION(jpi_crs,jpj_crs), INTENT(out) :: p_glam_crs  ! Coarse grid longitude
 
       !! Local variables
-      INTEGER ::  ji, jj, jk                   ! dummy loop indices
-      INTEGER :: ijie,ijis,ijje,ijjs,ijpk
+      INTEGER :: ji, jj, jk                   ! dummy loop indices
+      INTEGER :: ijis, ijjs
 
   
-      !! Initialize output fields
-      p_cgphi(:,:) = 0.e0
-      p_cglam(:,:) = 0.e0
-
-       DO ji = 2, nlei_crs
-
-         IF ( cd_type == 'T' .OR. cd_type == 'V' )  ijis = mis_crs(ji) + mxbinctr 
-         IF ( cd_type == 'U' .OR. cd_type == 'F' )  ijis = mie_crs(ji)
-
-         DO jj = njstart, njend
-     
-            IF ( cd_type == 'T' .OR. cd_type == 'U' ) ijjs = mjs_crs(jj) + mybinctr                  
-            IF ( cd_type == 'V' .OR. cd_type == 'F' ) ijjs = mje_crs(jj)
-
-            p_cgphi(ji,jj) = p_pgphi(ijis,ijjs)
-            p_cglam(ji,jj) = p_pglam(ijis,ijjs)
-             
-         ENDDO
-
-      ENDDO
-
-
-       ! Retroactively add back the boundary halo cells.
-
-         CALL crs_lbc_lnk( p_cgphi(:,:),cd_type,1.0 )
-         CALL crs_lbc_lnk( p_cglam(:,:),cd_type,1.0 )
+      SELECT CASE ( cd_type )
+         CASE ( 'T' )
+            DO jj =  nldj_crs, nlej_crs
+               ijjs = mjs_crs(jj) + mybinctr
+               DO ji = 2, nlei_crs
+                  ijis = mis_crs(ji) + mxbinctr 
+                  p_gphi_crs(ji,jj) = p_gphi(ijis,ijjs)
+                  p_glam_crs(ji,jj) = p_glam(ijis,ijjs)
+               ENDDO
+            ENDDO
+         CASE ( 'U' )
+            DO jj =  nldj_crs, nlej_crs
+               ijjs = mjs_crs(jj) + mybinctr                  
+               DO ji = 2, nlei_crs
+                  ijis = mis_crs(ji)
+                  p_gphi_crs(ji,jj) = p_gphi(ijis,ijjs)
+                  p_glam_crs(ji,jj) = p_glam(ijis,ijjs)
+               ENDDO
+            ENDDO
+         CASE ( 'V' )
+            DO jj =  nldj_crs, nlej_crs
+               ijjs = mjs_crs(jj)
+               DO ji = 2, nlei_crs
+                  ijis = mis_crs(ji) + mxbinctr 
+                  p_gphi_crs(ji,jj) = p_gphi(ijis,ijjs)
+                  p_glam_crs(ji,jj) = p_glam(ijis,ijjs)
+               ENDDO
+            ENDDO
+         CASE ( 'F' )
+            DO jj =  nldj_crs, nlej_crs
+               ijjs = mjs_crs(jj)
+               DO ji = 2, nlei_crs
+                  ijis = mis_crs(ji)
+                  p_gphi_crs(ji,jj) = p_gphi(ijis,ijjs)
+                  p_glam_crs(ji,jj) = p_glam(ijis,ijjs)
+               ENDDO
+            ENDDO
+      END SELECT
+
+      ! Retroactively add back the boundary halo cells.
+      CALL crs_lbc_lnk( p_gphi_crs, cd_type, 1.0 )
+      CALL crs_lbc_lnk( p_glam_crs, cd_type, 1.0 )
          
       ! Fill up jrow=1 which is zeroed out or not handled by lbc_lnk and lbc_nfd
-
-    DO ji = 2, nlei_crs
-
-         IF ( cd_type == 'T' .OR. cd_type == 'V' )  ijis = mis_crs(ji) + mxbinctr 
-         IF ( cd_type == 'U' .OR. cd_type == 'F' )  ijis = mie_crs(ji)
-
-         p_cgphi(ji,1) = p_pgphi(ijis,1)
-         p_cglam(ji,1) = p_pglam(ijis,1)
-
-      ENDDO
-
-      ! Fill i=1, i=jpi at j=1
-      p_cgphi(1,1) = p_cgphi(jpi_crsm1,1)
-      p_cglam(1,1) = p_cglam(jpi_crsm1,1)
-      p_cgphi(jpi_crs,1) = p_cgphi(2,1)
-      p_cglam(jpi_crs,1) = p_cglam(2,1)
-      ! Fill upper-right corner i=1, j=jpj_crs
- 
+      SELECT CASE ( cd_type )
+         CASE ( 'T', 'V' )
+            DO ji = 2, nlei_crs
+               ijis = mis_crs(ji) + mxbinctr 
+               p_gphi_crs(ji,1) = p_gphi(ijis,1)
+               p_glam_crs(ji,1) = p_glam(ijis,1)
+            ENDDO
+         CASE ( 'U', 'F' )
+            DO ji = 2, nlei_crs
+               ijis = mis_crs(ji) 
+               p_gphi_crs(ji,1) = p_gphi(ijis,1)
+               p_glam_crs(ji,1) = p_glam(ijis,1)
+            ENDDO
+      END SELECT
+
+      !                                             ! Fill i=1, i=jpi at j=1
+      p_gphi_crs(1      ,1) = p_gphi(jpi_crsm1,1)
+      p_glam_crs(1      ,1) = p_glam(jpi_crsm1,1)
+      !                                             ! Fill upper-right corner i=1, j=jpj_crs
+      p_gphi_crs(jpi_crs,1) = p_gphi(2        ,1)
+      p_glam_crs(jpi_crs,1) = p_glam(2        ,1)
+      !                                            
    END SUBROUTINE crs_dom_coordinates
 
@@ -211,6 +232,6 @@
 
       !! Local variables
-      INTEGER                                 :: ji, jj, jk     ! dummy loop indices
-      INTEGER                                 :: ijie,ijis,ijje,ijjs,ijrs
+      INTEGER :: ji, jj, jk     ! dummy loop indices
+      INTEGER :: ijie,ijis,ijje,ijjs,ijrs
   
       !!----------------------------------------------------------------  
@@ -221,5 +242,5 @@
             ijie = mie_crs(ji)
             ijis = mis_crs(ji)
-            DO jj = njstart, njend
+            DO jj = nldj_crs, nlej_crs
                ijje = mje_crs(jj) 
                ijjs = mjs_crs(jj)                   
@@ -262,6 +283,5 @@
 
 
-   SUBROUTINE crsfun_wgt( cd_type, cd_op, p_pmask, p_e1, p_e2, p_fse3, &
-      &                   p_cfield2d_1, p_cfield2d_2, p_cfield3d_1, p_cfield3d_2 )
+   SUBROUTINE crs_dom_facvol( p_mask, cd_type, p_e1, p_e2, p_e3, p_fld1_crs, p_fld2_crs )
       !!----------------------------------------------------------------
       !!               *** SUBROUTINE crsfun_wgt ***
@@ -304,247 +324,125 @@
       !! 
       !!  Arguments
-      CHARACTER(len=1),                 INTENT(in) :: cd_type  ! grid type U,V 
-      CHARACTER(len=3),                 INTENT(in) :: cd_op    ! operation sum or average
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: p_pmask  ! Parent grid U,V mask
-      REAL(wp), DIMENSION(jpi,jpj),     INTENT(in) :: p_e1     ! Parent grid U,V scale factors (e1)
-      REAL(wp), DIMENSION(jpi,jpj),     INTENT(in) :: p_e2     ! Parent grid U,V scale factors (e2)
-
-      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)
+      CHARACTER(len=1),                         INTENT(in)  :: cd_type  ! grid type U,V 
+      REAL(wp), DIMENSION(jpi,jpj,jpk)        , INTENT(in)  :: p_mask  ! Parent grid U,V mask
+      REAL(wp), DIMENSION(jpi,jpj)            , INTENT(in)  :: p_e1     ! Parent grid U,V scale factors (e1)
+      REAL(wp), DIMENSION(jpi,jpj)            , INTENT(in)  :: p_e2     ! Parent grid U,V scale factors (e2)
+      REAL(wp), DIMENSION(jpi,jpj,jpk)        , INTENT(in)  :: p_e3     ! Parent grid vertical level thickness (fse3u, fse3v)
+
+      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(out) :: p_fld1_crs ! Coarse grid box 3D quantity 
+      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(out) :: p_fld2_crs ! Coarse grid box 3D quantity 
 
       !! Local variables
-      INTEGER                                 ::  ji, jj, jk, jii, jjj     ! dummy loop indices
-      INTEGER                                 :: ijie,ijis,ijje,ijjs,ijpk
-      REAL(wp)                                :: zdAm                      ! masked face area
-      REAL(wp), DIMENSION(:,:),   POINTER     :: ze1, ze2
-      REAL(wp), DIMENSION(:,:,:), POINTER     :: ze3       
-      REAL(wp), DIMENSION(:,:),   POINTER     :: zcfield2d_1, zcfield2d_2
-      REAL(wp), DIMENSION(:,:,:), POINTER     :: zcfield3d_1, zcfield3d_2
-  
+      REAL(wp)                                :: zdAm
+      INTEGER                                 :: ji, jj, jk       ! dummy loop indices
+      INTEGER                                 :: ii, ij, ijie,ijje
+
+      REAL(wp), DIMENSION(:,:,:), POINTER     :: zvol      
       !!----------------------------------------------------------------  
-      ! Initialize      
-
-      ! Arrays, scalars initialization 
-      CALL wrk_alloc(jpi    , jpj         , ze1, ze2 )
-      CALL wrk_alloc(jpi    , jpj    , jpk, ze3 )
-      CALL wrk_alloc(jpi_crs, jpj_crs,      zcfield2d_1,  zcfield2d_2 )
-      CALL wrk_alloc(jpi_crs, jpj_crs, jpk, zcfield3d_1,  zcfield3d_2 )
-
-      ze1(:,:) = 1.0
-      ze2(:,:) = 1.0
-      ze3(:,:,:) = 1.0
-      zcfield2d_1(:,:)  = 0.0
-      zcfield2d_2(:,:)  = 0.0
-      zcfield3d_1(:,:,:) = 0.0
-      zcfield3d_2(:,:,:) = 0.0
-    
-      ijpk = jpk
-
-     ! Control of arguments
-      ze1(:,:) = p_e1(:,:)
-      ze2(:,:) = p_e2(:,:)
-
-      IF ( PRESENT(p_cfield2d_1) ) p_cfield2d_1(:,:) = 0.0
-      IF ( PRESENT(p_cfield2d_2) ) p_cfield2d_2(:,:) = 0.0
-      IF ( PRESENT(p_cfield3d_1) ) p_cfield3d_1(:,:,:) = 0.0
-      IF ( PRESENT(p_cfield3d_2) ) p_cfield3d_2(:,:,:) = 0.0
-
-      IF ( PRESENT(p_fse3) ) ze3(:,:,:) = p_fse3(:,:,:)
-
-
-          DO jk = 1, ijpk    
-
-             zcfield2d_1(:,:) = 0.0 ; zcfield2d_2(:,:) = 0.0        
-            DO ji = 2, nlei_crs
-                ijie = mie_crs(ji)
-                ijis = mis_crs(ji)
-
-
-                DO jj = njstart, njend
-                   ijje = mje_crs(jj) 
-                   ijjs = mjs_crs(jj)                   
-                   
-                   IF ( cd_op == 'POS' ) THEN      !cc
-                   
-                      IF ( nn_factx == 3 .AND. nn_facty == 3) THEN
-                      
-                         SELECT CASE ( cd_type )
-                         
-                            CASE ( 'T' )
-
-                               SELECT CASE ( mje_crs(jj)-mjs_crs(jj) )
-                                     
-                                  CASE( 0, 1 )    ! Si à la frontière sud on a pas assez de maille de la grille mère
-                                      
-                                     zcfield2d_1(ji,jj) = ze1(ijie-1,ijje  ) * nn_factx
-                                     zcfield2d_2(ji,jj) = ze2(ijie-1,ijje  ) * nn_facty
-                                       
-                                  CASE DEFAULT
    
-                                     zcfield2d_1(ji,jj) = ze1(ijie-1,ijje-1) * nn_factx
-                                     zcfield2d_2(ji,jj) = ze2(ijie-1,ijje-1) * nn_facty
- 
-                                  END SELECT
-                             
-                            CASE ( 'U' )
-                            
-                               SELECT CASE ( mje_crs(jj)-mjs_crs(jj) )
-                               
-                                  CASE( 0, 1 )    ! Si à la frontière sud on a pas assez de maille de la grille mère
-                                  
-                                     zcfield2d_1(ji,jj) = ze1(ijie  ,ijje  ) * nn_factx                            
-                                     zcfield2d_2(ji,jj) = ze2(ijie  ,ijje  ) * nn_facty
-                                     
-                                  CASE DEFAULT
-                                     
-                                     zcfield2d_1(ji,jj) = ze1(ijie  ,ijje-1) * nn_factx                            
-                                     zcfield2d_2(ji,jj) = ze2(ijie  ,ijje-1) * nn_facty
-                                     
-                               END SELECT
-                                
-                            CASE ( 'V' )
-
-                               zcfield2d_1(ji,jj) = ze1(ijie-1,ijje  ) * nn_factx
-                               zcfield2d_2(ji,jj) = ze2(ijie-1,ijje  ) * nn_facty
-                                                              
-                            CASE ( 'F' )
-
-                               zcfield2d_1(ji,jj) = ze1(ijie  ,ijje  ) * nn_factx                         
-                               zcfield2d_2(ji,jj) = ze2(ijie  ,ijje  ) * nn_facty
-                               
-                         END SELECT
-                      ENDIF
-                   ENDIF                    
-                   
-
-                   IF ( cd_op == 'WGT' ) THEN
-
-                      zdAm = 0.0
-
-                      IF ( cd_type == 'V' ) THEN
-                         ! 
-                         DO jii = ijis, ijie
-                            zdAm = zdAm + ( ze1(jii,ijje) * ze3(jii,ijje,jk) * p_pmask(jii,ijje,jk) ) 
-                         ENDDO
-                         IF ( zdAm /= 0 ) zcfield3d_1(ji,jj,jk) = zdAm 
-    
-                      ELSEIF ( cd_type == 'U') THEN
-                         DO jjj = ijjs, ijje
-                            zdAm = zdAm + ( ze2(ijie,jjj) * ze3(ijie,jjj,jk) * p_pmask(ijie,jjj,jk) )
-                         ENDDO
-                         IF ( zdAm /= 0 ) zcfield3d_1(ji,jj,jk) = zdAm 
-
-                      ELSEIF ( cd_type == 'W' ) THEN
-                         DO jii = ijis, ijie
-                            DO jjj = ijjs, ijje
-                               zdAm = zdAm + ( ze1(jii,jjj) * ze2(jii,jjj) * p_pmask(jii,jjj,jk) )
-                            ENDDO
-                         ENDDO
-                         IF ( zdAm /= 0 ) zcfield3d_1(ji,jj,jk) = zdAm
-
-                      ELSEIF ( cd_type == 'T' ) THEN
-                         DO jii = ijis, ijie
-                            DO jjj = ijjs, ijje
-                               zdAm = zdAm + ( ze1(jii,jjj) * ze2(jii,jjj) * ze3(jii,jjj,jk) * p_pmask(jii,jjj,jk) )
-                            ENDDO
-                         ENDDO
-                         IF ( zdAm /= 0 ) zcfield3d_1(ji,jj,jk) = zdAm
-
-                      ELSE
-  
-                         ! jes. Add a stop?
-
-                      ENDIF
-
-                   ENDIF
-
-                   IF ( cd_op == 'VOL' ) THEN
-
-                      zdAm = 0.0
-
-                      IF ( cd_type == 'W' .OR. cd_type == 'T' ) THEN
-
-                         DO jii = ijis, ijie
-                            DO jjj = ijjs, ijje
-                               zcfield3d_1(ji,jj,jk) = zcfield3d_1(ji,jj,jk) + ( ze1(jii,jjj) * ze2(jii,jjj) * ze3(jii,jjj,jk) )
-                               zdAm = zdAm + ( ze1(jii,jjj) * ze2(jii,jjj) * ze3(jii,jjj,jk) * p_pmask(jii,jjj,jk) )
-                            ENDDO
-                         ENDDO
-                         IF ( zcfield3d_1(ji,jj,jk) /= 0 ) zcfield3d_2(ji,jj,jk) = zdAm / zcfield3d_1(ji,jj,jk)
-
-                      ELSE
-                          ! jes. add a stop?
-                      ENDIF
-
-                   ENDIF
-                                    
-                ENDDO
-             ENDDO
-
-          ENDDO
-
-! Retroactively add back the boundary halo cells.
-
-
-
-         ! Take care of the 2D arrays
-         IF ( cd_op == 'SUM' .OR.  cd_op == 'POS') THEN
-            IF ( PRESENT(p_cfield2d_1) ) THEN
-               p_cfield2d_1(:,:) = zcfield2d_1(:,:)
-               CALL crs_lbc_lnk( p_cfield2d_1(:,:),cd_type,1.0, pval=1.0)
-            ENDIF
-            
-            IF ( PRESENT(p_cfield2d_2) ) THEN
-               p_cfield2d_2(:,:) = zcfield2d_2(:,:)
-               CALL crs_lbc_lnk( p_cfield2d_2(:,:),cd_type,1.0, pval=1.0 )
-               
-               IF ( cd_op == 'SUM') THEN 
-                  DO jii = 1 , jpiglo_crs
-                  p_cfield2d_2(jii,1) = p_cfield2d_2(jii,1) * 3
-                  ENDDO
-               ENDIF
-            ENDIF
-
-         ELSE
-
-            CALL crs_lbc_lnk( zcfield2d_1(:,:),cd_type,1.0 ) 
-            IF ( PRESENT(p_cfield2d_1) ) p_cfield2d_1(:,:) = zcfield2d_1(:,:)
-            CALL crs_lbc_lnk( zcfield2d_2(:,:),cd_type,1.0 ) 
-            IF ( PRESENT(p_cfield2d_2) ) p_cfield2d_2(:,:) = zcfield2d_2(:,:)
-
-         ENDIF
-
-         ! Take care now of 3d arrays
-         IF ( cd_op == 'SUM' .OR. cd_op == 'VOL' .OR. cd_op == 'POS'  ) THEN
-            CALL crs_lbc_lnk( zcfield3d_1(:,:,:),cd_type,1.0 ) 
-            IF ( PRESENT(p_cfield3d_1) ) p_cfield3d_1(:,:,:) = zcfield3d_1(:,:,:)
-            CALL crs_lbc_lnk( zcfield3d_2(:,:,:),cd_type,1.0 ) 
-            IF ( PRESENT(p_cfield3d_2) ) p_cfield3d_2(:,:,:) = zcfield3d_2(:,:,:)
-         ELSE
-            p_cfield3d_1(:,:,:) = zcfield3d_1(:,:,:)
-            CALL crs_lbc_lnk( p_cfield3d_1(:,:,:),cd_type,1.0 )
-
-            ! Fill upper-right corner i=1, j=jpj_crs
-            IF ( nperio == 4 ) THEN
-               p_cfield3d_1(1,jpj_crs,:) = p_cfield3d_1(jpi_crsm1,jpj_crs-2,:)
-            ELSEIF ( nperio == 6 ) THEN
-               p_cfield3d_1(1,jpj_crs,:) = p_cfield3d_1(jpi_crs,jpj_crsm1,:)
-            ENDIF
-
-         ENDIF
-
-      CALL wrk_dealloc(jpi    , jpj         , ze1, ze2 )
-      CALL wrk_dealloc(jpi    , jpj    , jpk, ze3 )
-      CALL wrk_dealloc(jpi_crs, jpj_crs,      zcfield2d_1,  zcfield2d_2 )
-      CALL wrk_dealloc(jpi_crs, jpj_crs, jpk, zcfield3d_1,  zcfield3d_2 )
-
-   END SUBROUTINE crsfun_wgt
-
-
-   SUBROUTINE crsfun_UV( p_e1_e2, cd_type, psgn, p_pmask, p_fse3, p_pfield, p_surf_crs, p_cfield3d )
+      CALL wrk_alloc( jpi, jpj, jpk, zvol )
+
+      DO jk = 1, jpk
+         zvol(:,:,jk) =  p_e1(:,:) * p_e2(:,:) * p_e3(:,:,jk) 
+      ENDDO
+
+      DO jk = 1, jpk           
+         DO ji = nistr, niend, nn_factx
+            DO jj   = njstr, njend, nn_facty
+               ii   = ( ji - mis_crs(2) ) * rfactx_r + 2                 ! cordinate in parent grid
+               ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+               ijje = mje_crs(ij)
+               ijie = mie_crs(ii)
+               !
+               p_fld1_crs(ii,ij,jk) =  zvol(ji,jj  ,jk) + zvol(ji+1,jj  ,jk) + zvol(ji+2,jj  ,jk)  &
+                   &                 + zvol(ji,jj+1,jk) + zvol(ji+1,jj+1,jk) + zvol(ji+2,jj+1,jk)  &
+                   &                 + zvol(ji,jj+2,jk) + zvol(ji+1,jj+2,jk) + zvol(ji+2,jj+2,jk)  
+            ENDDO
+         ENDDO
+      ENDDO
+
+      IF( cd_type == 'T' ) THEN
+         DO jk = 1, jpk           
+            DO ji = nistr, niend, nn_factx
+               DO jj   = njstr, njend, nn_facty
+                  ii   = ( ji - mis_crs(2) ) * rfactx_r + 2                 ! cordinate in parent grid
+                  ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                  ijje = mje_crs(ij)
+                  ijie = mie_crs(ii)
+                  !
+                  zdAm =  zvol(ji  ,jj  ,jk) * p_mask(ji  ,jj  ,jk)  &
+                    &   + zvol(ji+1,jj  ,jk) * p_mask(ji+1,jj  ,jk)  &
+                    &   + zvol(ji+2,jj  ,jk) * p_mask(ji+2,jj  ,jk)  &
+                    &   + zvol(ji  ,jj+1,jk) * p_mask(ji  ,jj+1,jk)  &
+                    &   + zvol(ji+1,jj+1,jk) * p_mask(ji+1,jj+1,jk)  &
+                    &   + zvol(ji+2,jj+1,jk) * p_mask(ji+2,jj+1,jk)  &
+                    &   + zvol(ji  ,jj+2,jk) * p_mask(ji  ,jj+2,jk)  &
+                    &   + zvol(ji+1,jj+2,jk) * p_mask(ji+1,jj+2,jk)  &
+                    &   + zvol(ji+2,jj+2,jk) * p_mask(ji+2,jj+2,jk)
+                    ! 
+                   IF( p_fld1_crs(ii,ij,jk) /= 0._wp ) p_fld2_crs(ii,ij,jk) = zdAm / p_fld1_crs(ii,ij,jk) 
+                   ! 
+               ENDDO
+            ENDDO
+         ENDDO
+      ENDIF
+      !
+      IF( cd_type == 'W' ) THEN
+         DO jk = 2, jpk           
+            DO ji = nistr, niend, nn_factx
+               DO jj   = njstr, njend, nn_facty
+                  ii   = ( ji - mis_crs(2) ) * rfactx_r + 2                 ! cordinate in parent grid
+                  ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                  ijje = mje_crs(ij)
+                  ijie = mie_crs(ii)
+                  !
+                  zdAm =  zvol(ji  ,jj  ,jk) * p_mask(ji  ,jj  ,jk-1)  &
+                    &   + zvol(ji+1,jj  ,jk) * p_mask(ji+1,jj  ,jk-1)  &
+                    &   + zvol(ji+2,jj  ,jk) * p_mask(ji+2,jj  ,jk-1)  &
+                    &   + zvol(ji  ,jj+1,jk) * p_mask(ji  ,jj+1,jk-1)  &
+                    &   + zvol(ji+1,jj+1,jk) * p_mask(ji+1,jj+1,jk-1)  &
+                    &   + zvol(ji+2,jj+1,jk) * p_mask(ji+2,jj+1,jk-1)  &
+                    &   + zvol(ji  ,jj+2,jk) * p_mask(ji  ,jj+2,jk-1)  &
+                    &   + zvol(ji+1,jj+2,jk) * p_mask(ji+1,jj+2,jk-1)  &
+                    &   + zvol(ji+2,jj+2,jk) * p_mask(ji+2,jj+2,jk-1)
+                  !
+                  IF( p_fld1_crs(ii,ij,jk) /= 0._wp ) p_fld2_crs(ii,ij,jk) = zdAm / p_fld1_crs(ii,ij,jk) 
+                  !
+               ENDDO
+            ENDDO
+         ENDDO
+         DO ji = nistr, niend, nn_factx
+            DO jj   = njstr, njend, nn_facty
+               ii   = ( ji - mis_crs(2) ) * rfactx_r + 2                 ! cordinate in parent grid
+               ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+               ijje = mje_crs(ij)
+               ijie = mie_crs(ii)
+               !
+               zdAm =  zvol(ji  ,jj  ,1) * p_mask(ji  ,jj  ,1)  &
+                 &   + zvol(ji+1,jj  ,1) * p_mask(ji+1,jj  ,1)  &
+                 &   + zvol(ji+2,jj  ,1) * p_mask(ji+2,jj  ,1)  &
+                 &   + zvol(ji  ,jj+1,1) * p_mask(ji  ,jj+1,1)  &
+                 &   + zvol(ji+1,jj+1,1) * p_mask(ji+1,jj+1,1)  &
+                 &   + zvol(ji+2,jj+1,1) * p_mask(ji+2,jj+1,1)  &
+                 &   + zvol(ji  ,jj+2,1) * p_mask(ji  ,jj+2,1)  &
+                 &   + zvol(ji+1,jj+2,1) * p_mask(ji+1,jj+2,1)  &
+                 &   + zvol(ji+2,jj+2,1) * p_mask(ji+2,jj+2,1)
+               ! 
+               IF( p_fld1_crs(ii,ij,1) /= 0._wp ) p_fld2_crs(ii,ij,1) = zdAm / p_fld2_crs(ii,ij,1) 
+               ! 
+            ENDDO
+         ENDDO
+      ENDIF
+      
+      !                                             !  Retroactively add back the boundary halo cells.
+      CALL crs_lbc_lnk( p_fld1_crs, cd_type, 1.0 ) 
+      CALL crs_lbc_lnk( p_fld2_crs, cd_type, 1.0 ) 
+      !
+      CALL wrk_dealloc( jpi, jpj, jpk, zvol )
+      !
+   END SUBROUTINE crs_dom_facvol
+
+
+   SUBROUTINE crs_dom_ope_3d( p_fld, cd_op, cd_type, p_mask, p_fld_crs, p_e12, p_e3, p_surf_crs, p_mask_crs )
       !!----------------------------------------------------------------
       !!               *** SUBROUTINE crsfun_UV ***
@@ -571,334 +469,826 @@
       !! 
       !!  Arguments
-      REAL(wp), DIMENSION(jpi,jpj),     INTENT(in)     :: p_e1_e2    ! Parent grid U,V scale factors (e1 or e2)
-      CHARACTER(len=1),                 INTENT(in)     :: cd_type    ! grid type U,V 
-      REAL(wp),                         INTENT(in)     :: psgn       ! sign change option across north fold
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)     :: p_pmask    ! Parent grid U,V mask
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)     :: p_fse3     ! Parent grid vertical level thickness (fse3u, fse3v)
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)     :: p_pfield   ! U or V on parent grid
+      REAL(wp), DIMENSION(jpi,jpj,jpk),         INTENT(in)           :: p_fld   ! T, U, V or W on parent grid
+      CHARACTER(len=3),                         INTENT(in)           :: cd_op    ! Operation SUM, MAX or MIN
+      CHARACTER(len=1),                         INTENT(in)           :: cd_type    ! grid type U,V 
+      REAL(wp), DIMENSION(jpi,jpj,jpk),         INTENT(in)           :: p_mask    ! Parent grid T,U,V mask
+      REAL(wp), DIMENSION(jpi,jpj),             INTENT(in), OPTIONAL :: p_e12    ! Parent grid T,U,V scale factors (e1 or e2)
+      REAL(wp), DIMENSION(jpi,jpj,jpk),         INTENT(in), OPTIONAL :: p_e3     ! Parent grid vertical level thickness (fse3u, fse3v)
       REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(in), OPTIONAL :: p_surf_crs ! Coarse grid area-weighting denominator    
-
-      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(out)  :: p_cfield3d ! Coarse grid box 3D quantity 
+      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(in), OPTIONAL :: p_mask_crs    ! Coarse grid T,U,V mask
+
+      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(out)          :: p_fld_crs ! Coarse grid box 3D quantity 
 
       !! Local variables
-      INTEGER  :: ji, jj, jk , jii, jjj                  ! dummy loop indices
-      INTEGER  :: ijie, ijis, ijje, ijjs
-      REAL(wp), DIMENSION(:,:,:), POINTER :: zsurfcrs   
+      INTEGER  :: ji, jj, jk                 ! dummy loop indices
+      INTEGER  :: ijie, ijje, ii, ij
+      REAL(wp) :: zflcrs, zsfcrs   
+      REAL(wp) :: zeps = 1.e20    
+      REAL(wp), DIMENSION(:,:,:), POINTER :: zsurf   
 
       !!----------------------------------------------------------------  
-
-      CALL wrk_alloc( jpi_crs, jpj_crs, jpk, zsurfcrs )
-      zsurfcrs(:,:,:) = 1.0
-      IF ( PRESENT(p_surf_crs) ) THEN
-         WHERE ( p_surf_crs /= 0 ) zsurfcrs(:,:,:) = 1.0/p_surf_crs(:,:,:)
-      ENDIF
-
-      DO jk = 1, jpk    
-
-        DO ji = 2, nlei_crs
-            ijie = mie_crs(ji)
-            ijis = mis_crs(ji)
-
-            DO jj = njstart, njend
-               ijje = mje_crs(jj) 
-               ijjs = mjs_crs(jj)                   
-
-               IF ( cd_type == 'V' ) THEN
-
-                  DO jii = ijis, ijie
-                     p_cfield3d(ji,jj,jk) = p_cfield3d(ji,jj,jk) &
-                       & + ( p_pfield(jii,ijje,jk) * p_e1_e2(jii,ijje) * p_fse3(jii,ijje,jk) * p_pmask(jii,ijje,jk) ) 
+   
+
+      SELECT CASE ( cd_op )
+      
+         CASE ( 'VOL' )
+      
+            CALL wrk_alloc( jpi, jpj, jpk, zsurf )
+            DO jk = 1, jpk
+              zsurf(:,:,jk) =  p_e12(:,:) * p_e3(:,:,jk)
+            ENDDO
+         
+            SELECT CASE ( cd_type )
+            
+               CASE( 'T' )
+         
+                  DO jk = 1, jpk
+                    
+                     DO ji = nistr, niend, nn_factx
+                        DO jj = njstr, njend, nn_facty
+                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                           ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                           ijje = mje_crs(ij)
+                           ijie = mie_crs(ii)                  
+ 
+                           zflcrs =  p_fld(ji  ,jj  ,jk) * zsurf(ji  ,jj  ,jk) * p_mask(ji  ,jj  ,jk) &
+                             &     + p_fld(ji+1,jj  ,jk) * zsurf(ji+1,jj  ,jk) * p_mask(ji+1,jj  ,jk) &
+                             &     + p_fld(ji+2,jj  ,jk) * zsurf(ji+2,jj  ,jk) * p_mask(ji+2,jj  ,jk) &
+                             &     + p_fld(ji  ,jj+1,jk) * zsurf(ji  ,jj+1,jk) * p_mask(ji  ,jj+1,jk) &
+                             &     + p_fld(ji+1,jj+1,jk) * zsurf(ji+1,jj+1,jk) * p_mask(ji+1,jj+1,jk) &
+                             &     + p_fld(ji+2,jj+1,jk) * zsurf(ji+2,jj+1,jk) * p_mask(ji+2,jj+1,jk) &
+                             &     + p_fld(ji  ,jj+2,jk) * zsurf(ji  ,jj+2,jk) * p_mask(ji  ,jj+2,jk) &
+                             &     + p_fld(ji+1,jj+2,jk) * zsurf(ji+1,jj+2,jk) * p_mask(ji+1,jj+2,jk) &
+                             &     + p_fld(ji+2,jj+2,jk) * zsurf(ji+2,jj+2,jk) * p_mask(ji+2,jj+2,jk) 
+ 
+                           zsfcrs =  zsurf(ji  ,jj  ,jk) * p_mask(ji  ,jj  ,jk) &
+                             &     + zsurf(ji+1,jj  ,jk) * p_mask(ji+1,jj  ,jk) &
+                             &     + zsurf(ji+2,jj  ,jk) * p_mask(ji+2,jj  ,jk) &
+                             &     + zsurf(ji  ,jj+1,jk) * p_mask(ji  ,jj+1,jk) &
+                             &     + zsurf(ji+1,jj+1,jk) * p_mask(ji+1,jj+1,jk) &
+                             &     + zsurf(ji+2,jj+1,jk) * p_mask(ji+2,jj+1,jk) &
+                             &     + zsurf(ji  ,jj+2,jk) * p_mask(ji  ,jj+2,jk) &
+                             &     + zsurf(ji+1,jj+2,jk) * p_mask(ji+1,jj+2,jk) &
+                             &     + zsurf(ji+2,jj+2,jk) * p_mask(ji+2,jj+2,jk) 
+                           !
+                           p_fld_crs(ii,ij,jk) = zflcrs
+                           IF( zsfcrs /= 0.0 )  p_fld_crs(ii,ij,jk) = zflcrs / zsfcrs
+
+                        ENDDO      
+                     ENDDO
+                  ENDDO   
+            
+               CASE( 'W' )
+         
+                  DO jk = 2, jpk
+                    
+                     DO ji = nistr, niend, nn_factx
+                        DO jj = njstr, njend, nn_facty
+                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                           ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                           ijje = mje_crs(ij)
+                           ijie = mie_crs(ii)                  
+ 
+                           zflcrs =  p_fld(ji  ,jj  ,jk) * zsurf(ji  ,jj  ,jk) * p_mask(ji  ,jj  ,jk-1) &
+                             &     + p_fld(ji+1,jj  ,jk) * zsurf(ji+1,jj  ,jk) * p_mask(ji+1,jj  ,jk-1) &
+                             &     + p_fld(ji+2,jj  ,jk) * zsurf(ji+2,jj  ,jk) * p_mask(ji+2,jj  ,jk-1) &
+                             &     + p_fld(ji  ,jj+1,jk) * zsurf(ji  ,jj+1,jk) * p_mask(ji  ,jj+1,jk-1) &
+                             &     + p_fld(ji+1,jj+1,jk) * zsurf(ji+1,jj+1,jk) * p_mask(ji+1,jj+1,jk-1) &
+                             &     + p_fld(ji+2,jj+1,jk) * zsurf(ji+2,jj+1,jk) * p_mask(ji+2,jj+1,jk-1) &
+                             &     + p_fld(ji  ,jj+2,jk) * zsurf(ji  ,jj+2,jk) * p_mask(ji  ,jj+2,jk-1) &
+                             &     + p_fld(ji+1,jj+2,jk) * zsurf(ji+1,jj+2,jk) * p_mask(ji+1,jj+2,jk-1) &
+                             &     + p_fld(ji+2,jj+2,jk) * zsurf(ji+2,jj+2,jk) * p_mask(ji+2,jj+2,jk-1) 
+ 
+                           zsfcrs =  zsurf(ji  ,jj  ,jk) * p_mask(ji  ,jj  ,jk-1) &
+                             &     + zsurf(ji+1,jj  ,jk) * p_mask(ji+1,jj  ,jk-1) &
+                             &     + zsurf(ji+2,jj  ,jk) * p_mask(ji+2,jj  ,jk-1) &
+                             &     + zsurf(ji  ,jj+1,jk) * p_mask(ji  ,jj+1,jk-1) &
+                             &     + zsurf(ji+1,jj+1,jk) * p_mask(ji+1,jj+1,jk-1) &
+                             &     + zsurf(ji+2,jj+1,jk) * p_mask(ji+2,jj+1,jk-1) &
+                             &     + zsurf(ji  ,jj+2,jk) * p_mask(ji  ,jj+2,jk-1) &
+                             &     + zsurf(ji+1,jj+2,jk) * p_mask(ji+1,jj+2,jk-1) &
+                             &     + zsurf(ji+2,jj+2,jk) * p_mask(ji+2,jj+2,jk-1) 
+                           !
+                           p_fld_crs(ii,ij,jk) = zflcrs
+                           IF( zsfcrs /= 0.0 )  p_fld_crs(ii,ij,jk) = zflcrs / zsfcrs
+
+                        ENDDO      
+                     ENDDO
+                  ENDDO   
+
+                  DO ji = nistr, niend, nn_factx
+                     DO jj = njstr, njend, nn_facty
+                        ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                        ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                        ijje = mje_crs(ij)
+                        ijie = mie_crs(ii)                  
+
+                          zflcrs =   p_fld(ji  ,jj  ,1) * zsurf(ji  ,jj  ,1) * p_mask(ji  ,jj  ,1) &
+                             &     + p_fld(ji+1,jj  ,1) * zsurf(ji+1,jj  ,1) * p_mask(ji+1,jj  ,1) &
+                             &     + p_fld(ji+2,jj  ,1) * zsurf(ji+2,jj  ,1) * p_mask(ji+2,jj  ,1) &
+                             &     + p_fld(ji  ,jj+1,1) * zsurf(ji  ,jj+1,1) * p_mask(ji  ,jj+1,1) &
+                             &     + p_fld(ji+1,jj+1,1) * zsurf(ji+1,jj+1,1) * p_mask(ji+1,jj+1,1) &
+                             &     + p_fld(ji+2,jj+1,1) * zsurf(ji+2,jj+1,1) * p_mask(ji+2,jj+1,1) &
+                             &     + p_fld(ji  ,jj+2,1) * zsurf(ji  ,jj+2,1) * p_mask(ji  ,jj+2,1) &
+                             &     + p_fld(ji+1,jj+2,1) * zsurf(ji+1,jj+2,1) * p_mask(ji+1,jj+2,1) &
+                             &     + p_fld(ji+2,jj+2,1) * zsurf(ji+2,jj+2,1) * p_mask(ji+2,jj+2,1) 
+ 
+                           zsfcrs =  zsurf(ji  ,jj  ,1) * p_mask(ji  ,jj  ,1) &
+                             &     + zsurf(ji+1,jj  ,1) * p_mask(ji+1,jj  ,1) &
+                             &     + zsurf(ji+2,jj  ,1) * p_mask(ji+2,jj  ,1) &
+                             &     + zsurf(ji  ,jj+1,1) * p_mask(ji  ,jj+1,1) &
+                             &     + zsurf(ji+1,jj+1,1) * p_mask(ji+1,jj+1,1) &
+                             &     + zsurf(ji+2,jj+1,1) * p_mask(ji+2,jj+1,1) &
+                             &     + zsurf(ji  ,jj+2,1) * p_mask(ji  ,jj+2,1) &
+                             &     + zsurf(ji+1,jj+2,1) * p_mask(ji+1,jj+2,1) &
+                             &     + zsurf(ji+2,jj+2,1) * p_mask(ji+2,jj+2,1) 
+                         
+                        p_fld_crs(ii,ij,1) = zflcrs
+                        IF( zsfcrs /= 0.0 )  p_fld_crs(ii,ij,1) = zflcrs / zsfcrs
+
+                     ENDDO      
                   ENDDO
-                  p_cfield3d(ji,jj,jk) = p_cfield3d(ji,jj,jk) * zsurfcrs(ji,jj,jk) 
-
-               ELSEIF ( cd_type == 'U') THEN
-
-                  DO jjj = ijjs, ijje
-                     p_cfield3d(ji,jj,jk) = p_cfield3d(ji,jj,jk) &
-                       & + ( p_pfield(ijie,jjj,jk) * p_e1_e2(ijie,jjj) * p_fse3(ijie,jjj,jk) * p_pmask(ijie,jjj,jk) )
+
+              END SELECT
+
+              CALL wrk_dealloc( jpi, jpj, jpk, zsurf )
+
+         CASE ( 'SUM' )
+         
+            CALL wrk_alloc( jpi, jpj, jpk, zsurf )
+            DO jk = 1, jpk
+              zsurf(:,:,jk) =  p_e12(:,:) * p_e3(:,:,jk)
+            ENDDO
+         
+            SELECT CASE ( cd_type )
+            
+               CASE( 'T' )
+         
+                  DO jk = 1, jpk
+                     DO ji = nistr, niend, nn_factx
+                        DO jj = njstr, njend, nn_facty
+                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                           ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                           ijje = mje_crs(ij)
+                           ijie = mie_crs(ii)                  
+ 
+                           zflcrs  =  p_fld(ji  ,jj  ,jk) * zsurf(ji  ,jj  ,jk) * p_mask(ji  ,jj  ,jk) &
+                             &      + p_fld(ji+1,jj  ,jk) * zsurf(ji+1,jj  ,jk) * p_mask(ji+1,jj  ,jk) &
+                             &      + p_fld(ji+2,jj  ,jk) * zsurf(ji+2,jj  ,jk) * p_mask(ji+2,jj  ,jk) &
+                             &      + p_fld(ji  ,jj+1,jk) * zsurf(ji  ,jj+1,jk) * p_mask(ji  ,jj+1,jk) &
+                             &      + p_fld(ji+1,jj+1,jk) * zsurf(ji+1,jj+1,jk) * p_mask(ji+1,jj+1,jk) &
+                             &      + p_fld(ji+2,jj+1,jk) * zsurf(ji+2,jj+1,jk) * p_mask(ji+2,jj+1,jk) &
+                             &      + p_fld(ji  ,jj+2,jk) * zsurf(ji  ,jj+2,jk) * p_mask(ji  ,jj+2,jk) &
+                             &      + p_fld(ji+1,jj+2,jk) * zsurf(ji+1,jj+2,jk) * p_mask(ji+1,jj+2,jk) &
+                             &      + p_fld(ji+2,jj+2,jk) * zsurf(ji+2,jj+2,jk) * p_mask(ji+2,jj+2,jk) 
+                           !
+                           p_fld_crs(ii,ij,jk) = zflcrs
+                           !
+                        ENDDO      
+                     ENDDO
+                  ENDDO   
+            
+               CASE( 'W' )
+         
+                  DO jk = 2, jpk
+                     DO ji = nistr, niend, nn_factx
+                        DO jj = njstr, njend, nn_facty
+                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                           ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                           ijje = mje_crs(ij)
+                           ijie = mie_crs(ii)
+                           !                  
+                           zflcrs  =  p_fld(ji  ,jj  ,jk) * zsurf(ji  ,jj  ,jk) * p_mask(ji  ,jj  ,jk-1) &
+                             &      + p_fld(ji+1,jj  ,jk) * zsurf(ji+1,jj  ,jk) * p_mask(ji+1,jj  ,jk-1) &
+                             &      + p_fld(ji+2,jj  ,jk) * zsurf(ji+2,jj  ,jk) * p_mask(ji+2,jj  ,jk-1) &
+                             &      + p_fld(ji  ,jj+1,jk) * zsurf(ji  ,jj+1,jk) * p_mask(ji  ,jj+1,jk-1) &
+                             &      + p_fld(ji+1,jj+1,jk) * zsurf(ji+1,jj+1,jk) * p_mask(ji+1,jj+1,jk-1) &
+                             &      + p_fld(ji+2,jj+1,jk) * zsurf(ji+2,jj+1,jk) * p_mask(ji+2,jj+1,jk-1) &
+                             &      + p_fld(ji  ,jj+2,jk) * zsurf(ji  ,jj+2,jk) * p_mask(ji  ,jj+2,jk-1) &
+                             &      + p_fld(ji+1,jj+2,jk) * zsurf(ji+1,jj+2,jk) * p_mask(ji+1,jj+2,jk-1) &
+                             &      + p_fld(ji+2,jj+2,jk) * zsurf(ji+2,jj+2,jk) * p_mask(ji+2,jj+2,jk-1) 
+                           !
+                           p_fld_crs(ii,ij,jk) = zflcrs
+                           !
+                        ENDDO      
+                     ENDDO
+                  ENDDO   
+
+                  DO ji = nistr, niend, nn_factx
+                     DO jj = njstr, njend, nn_facty
+                        ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                        ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                        ijje = mje_crs(ij)
+                        ijie = mie_crs(ii)                  
+                        !
+                        zflcrs  =   p_fld(ji  ,jj  ,1) * zsurf(ji  ,jj  ,1) * p_mask(ji  ,jj  ,1) &
+                           &      + p_fld(ji+1,jj  ,1) * zsurf(ji+1,jj  ,1) * p_mask(ji+1,jj  ,1) &
+                           &      + p_fld(ji+2,jj  ,1) * zsurf(ji+2,jj  ,1) * p_mask(ji+2,jj  ,1) &
+                           &      + p_fld(ji  ,jj+1,1) * zsurf(ji  ,jj+1,1) * p_mask(ji  ,jj+1,1) &
+                           &      + p_fld(ji+1,jj+1,1) * zsurf(ji+1,jj+1,1) * p_mask(ji+1,jj+1,1) &
+                           &      + p_fld(ji+2,jj+1,1) * zsurf(ji+2,jj+1,1) * p_mask(ji+2,jj+1,1) &
+                           &      + p_fld(ji  ,jj+2,1) * zsurf(ji  ,jj+2,1) * p_mask(ji  ,jj+2,1) &
+                           &      + p_fld(ji+1,jj+2,1) * zsurf(ji+1,jj+2,1) * p_mask(ji+1,jj+2,1) &
+                           &      + p_fld(ji+2,jj+2,1) * zsurf(ji+2,jj+2,1) * p_mask(ji+2,jj+2,1) 
+                        !
+                        p_fld_crs(ii,ij,1) = zflcrs
+                        !
+                     ENDDO      
                   ENDDO
-                  p_cfield3d(ji,jj,jk) = p_cfield3d(ji,jj,jk) * zsurfcrs(ji,jj,jk)
-
-               ENDIF
-                 
-            ENDDO
-         ENDDO
-      ENDDO
-
-! Retroactively add back the boundary halo cells.
-
-      CALL crs_lbc_lnk( p_cfield3d(:,:,:),cd_type,psgn )
-      
-      CALL wrk_dealloc( jpi_crs, jpj_crs, jpk, zsurfcrs )
-
-   END SUBROUTINE crsfun_UV
-
-   SUBROUTINE crsfun_TW( p_e1e2t, cd_type, cd_op, p_cmask, p_ptmask, psgn, p_pfield2d, p_pfield3d_1, p_pfield3d_2, &
-      &                  p_cfield2d, p_cfield3d)
+           
+               CASE( 'V' )
+         
+                  DO jk = 1, jpk
+                     DO ji = nistr, niend, nn_factx
+                        DO jj = njstr, njend, nn_facty
+                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                           ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                           ijje = mje_crs(ij)
+                           ijie = mie_crs(ii)
+                           !                  
+                           zflcrs  =  p_fld(ji  ,ijje,jk) * zsurf(ji  ,ijje,jk) * p_mask(ji  ,ijje,jk) &
+                             &      + p_fld(ji+1,ijje,jk) * zsurf(ji+1,ijje,jk) * p_mask(ji+1,ijje,jk) &
+                             &      + p_fld(ji+2,ijje,jk) * zsurf(ji+2,ijje,jk) * p_mask(ji+2,ijje,jk) 
+                           !
+                           p_fld_crs(ii,ij,jk) = zflcrs
+                           !
+                         ENDDO      
+                     ENDDO
+                  ENDDO   
+
+            
+               CASE( 'U' )
+         
+                  DO jk = 1, jpk
+                     DO ji = nistr, niend, nn_factx
+                        DO jj = njstr, njend, nn_facty
+                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                           ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                           ijje = mje_crs(ij)
+                           ijie = mie_crs(ii)
+                           !                  
+                           zflcrs =  p_fld(ijie,jj  ,jk) * zsurf(ijie,jj  ,jk) * p_mask(ijie,jj  ,jk) &
+                             &     + p_fld(ijie,jj+1,jk) * zsurf(ijie,jj+1,jk) * p_mask(ijie,jj+1,jk) &
+                             &     + p_fld(ijie,jj+2,jk) * zsurf(ijie,jj+2,jk) * p_mask(ijie,jj+2,jk)
+                           !
+                           p_fld_crs(ii,ij,jk) = zflcrs
+                           !
+                        ENDDO      
+                     ENDDO
+                  ENDDO   
+
+              END SELECT
+
+              IF( PRESENT( p_surf_crs ) ) THEN
+                 WHERE ( p_surf_crs /= 0.0 ) p_fld_crs(:,:,:) = p_fld_crs(:,:,:) / p_surf_crs(:,:,:)
+              ENDIF
+
+              CALL wrk_dealloc( jpi, jpj, jpk, zsurf )
+
+         CASE ( 'MAX' )
+         
+            SELECT CASE ( cd_type )
+            
+               CASE( 'T' )
+         
+                  DO jk = 1, jpk
+                     DO ji = nistr, niend, nn_factx
+                        DO jj = njstr, njend, nn_facty
+                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                           ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                           ijje = mje_crs(ij)
+                           ijie = mie_crs(ii)                  
+ 
+                           zflcrs =  MAX( p_fld(ji  ,jj  ,jk) * p_mask(ji  ,jj  ,jk),  &
+                             &            p_fld(ji+1,jj  ,jk) * p_mask(ji+1,jj  ,jk),  &
+                             &            p_fld(ji+2,jj  ,jk) * p_mask(ji+2,jj  ,jk),  &
+                             &            p_fld(ji  ,jj+1,jk) * p_mask(ji  ,jj+1,jk),  &
+                             &            p_fld(ji+1,jj+1,jk) * p_mask(ji+1,jj+1,jk),  &
+                             &            p_fld(ji+2,jj+1,jk) * p_mask(ji+2,jj+1,jk),  &
+                             &            p_fld(ji  ,jj+2,jk) * p_mask(ji  ,jj+2,jk),  &
+                             &            p_fld(ji+1,jj+2,jk) * p_mask(ji+1,jj+2,jk),  &
+                             &            p_fld(ji+2,jj+2,jk) * p_mask(ji+2,jj+2,jk)   )
+                           !
+                           p_fld_crs(ii,ij,jk) = zflcrs
+                           !
+                        ENDDO      
+                     ENDDO
+                  ENDDO   
+            
+               CASE( 'W' )
+         
+                  DO jk = 2, jpk
+                     DO ji = nistr, niend, nn_factx
+                        DO jj = njstr, njend, nn_facty
+                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                           ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                           ijje = mje_crs(ij)
+                           ijie = mie_crs(ii)
+                           !                  
+                           zflcrs =  MAX( p_fld(ji  ,jj  ,jk) * p_mask(ji  ,jj  ,jk-1),  &
+                             &            p_fld(ji+1,jj  ,jk) * p_mask(ji+1,jj  ,jk-1),  &
+                             &            p_fld(ji+2,jj  ,jk) * p_mask(ji+2,jj  ,jk-1),  &
+                             &            p_fld(ji  ,jj+1,jk) * p_mask(ji  ,jj+1,jk-1),  &
+                             &            p_fld(ji+1,jj+1,jk) * p_mask(ji+1,jj+1,jk-1),  &
+                             &            p_fld(ji+2,jj+1,jk) * p_mask(ji+2,jj+1,jk-1),  &
+                             &            p_fld(ji  ,jj+2,jk) * p_mask(ji  ,jj+2,jk-1),  &
+                             &            p_fld(ji+1,jj+2,jk) * p_mask(ji+1,jj+2,jk-1),  &
+                             &            p_fld(ji+2,jj+2,jk) * p_mask(ji+2,jj+2,jk-1) )
+                           !
+                           p_fld_crs(ii,ij,jk) = zflcrs
+                           !
+                        ENDDO      
+                     ENDDO
+                  ENDDO   
+
+                  DO ji = nistr, niend, nn_factx
+                     DO jj = njstr, njend, nn_facty
+                        ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                        ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                        ijje = mje_crs(ij)
+                        ijie = mie_crs(ii)                  
+                        !
+                        zflcrs = MAX( p_fld(ji  ,jj  ,1) * p_mask(ji  ,jj  ,1),  &
+                           &          p_fld(ji+1,jj  ,1) * p_mask(ji+1,jj  ,1),  &
+                           &          p_fld(ji+2,jj  ,1) * p_mask(ji+2,jj  ,1),  &
+                           &          p_fld(ji  ,jj+1,1) * p_mask(ji  ,jj+1,1),  &
+                           &          p_fld(ji+1,jj+1,1) * p_mask(ji+1,jj+1,1),  &
+                           &          p_fld(ji+2,jj+1,1) * p_mask(ji+2,jj+1,1),  &
+                           &          p_fld(ji  ,jj+2,1) * p_mask(ji  ,jj+2,1),  &
+                           &          p_fld(ji+1,jj+2,1) * p_mask(ji+1,jj+2,1),  &
+                           &          p_fld(ji+2,jj+2,1) * p_mask(ji+2,jj+2,1)  ) 
+                         !
+                         p_fld_crs(ii,ij,1) = zflcrs
+                         !
+                     ENDDO      
+                  ENDDO
+           
+               CASE( 'V' )
+         
+                  DO jk = 1, jpk
+                     DO ji = nistr, niend, nn_factx
+                        DO jj = njstr, njend, nn_facty
+                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                           ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                           ijje = mje_crs(ij)
+                           ijie = mie_crs(ii)
+                           !                  
+                           zflcrs = MAX( p_fld(ji  ,ijje,jk) * p_mask(ji  ,ijje,jk),  &
+                             &           p_fld(ji+1,ijje,jk) * p_mask(ji+1,ijje,jk),  &
+                             &           p_fld(ji+2,ijje,jk) * p_mask(ji+2,ijje,jk) )
+                           !
+                           p_fld_crs(ii,ij,jk) = zflcrs
+                           !
+                        ENDDO      
+                     ENDDO
+                  ENDDO   
+
+            
+               CASE( 'U' )
+         
+                  DO jk = 1, jpk
+                     DO ji = nistr, niend, nn_factx
+                        DO jj = njstr, njend, nn_facty
+                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                           ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                           ijje = mje_crs(ij)
+                           ijie = mie_crs(ii)
+                           !                  
+                           Zflcrs = MAX( p_fld(ijie,jj  ,jk) * p_mask(ijie,jj  ,jk),  &
+                             &           p_fld(ijie,jj+1,jk) * p_mask(ijie,jj+1,jk),  &
+                             &           p_fld(ijie,jj+2,jk) * p_mask(ijie,jj+2,jk) )
+                           !
+                           p_fld_crs(ii,ij,jk) = zflcrs
+                           !
+                        ENDDO      
+                     ENDDO
+                  ENDDO   
+
+              END SELECT
+
+         CASE ( 'MIN' )
+            !   Search the min of masked grid cells
+            SELECT CASE ( cd_type )
+            
+               CASE( 'T' )
+         
+                  DO jk = 1, jpk
+                     DO ji = nistr, niend, nn_factx
+                        DO jj = njstr, njend, nn_facty
+                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                           ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                           ijje = mje_crs(ij)
+                           ijie = mie_crs(ii)                  
+                           
+                           zflcrs =  MIN( p_fld(ji  ,jj  ,jk) * ( 1. + ( 1. - p_mask(ji  ,jj  ,jk) ) * zeps ),  &
+                             &            p_fld(ji+1,jj  ,jk) * ( 1. + ( 1. - p_mask(ji+1,jj  ,jk) ) * zeps ),  &
+                             &            p_fld(ji+2,jj  ,jk) * ( 1. + ( 1. - p_mask(ji+2,jj  ,jk) ) * zeps ),  &
+                             &            p_fld(ji  ,jj+1,jk) * ( 1. + ( 1. - p_mask(ji  ,jj  ,jk) ) * zeps ),  &
+                             &            p_fld(ji+1,jj+1,jk) * ( 1. + ( 1. - p_mask(ji+1,jj+1,jk) ) * zeps ),  &
+                             &            p_fld(ji+2,jj+1,jk) * ( 1. + ( 1. - p_mask(ji+2,jj+2,jk) ) * zeps ),  &
+                             &            p_fld(ji  ,jj+2,jk) * ( 1. + ( 1. - p_mask(ji  ,jj  ,jk) ) * zeps ),  &
+                             &            p_fld(ji+1,jj+2,jk) * ( 1. + ( 1. - p_mask(ji+1,jj+1,jk) ) * zeps ),  &
+                             &            p_fld(ji+2,jj+2,jk) * ( 1. + ( 1. - p_mask(ji+2,jj+2,jk) ) * zeps )   )
+                           !
+                           p_fld_crs(ii,ij,jk) = zflcrs * p_mask_crs(ii,ij,jk)
+                           !
+                        ENDDO      
+                     ENDDO
+                  ENDDO   
+           
+               CASE( 'W' )
+         
+                  DO jk = 2, jpk
+                     DO ji = nistr, niend, nn_factx
+                        DO jj = njstr, njend, nn_facty
+                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                           ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                           ijje = mje_crs(ij)
+                           ijie = mie_crs(ii)                  
+                           
+                           zflcrs =  MIN( p_fld(ji  ,jj  ,jk) * ( 1. + ( 1. - p_mask(ji  ,jj  ,jk-1) ) * zeps ),  &
+                             &            p_fld(ji+1,jj  ,jk) * ( 1. + ( 1. - p_mask(ji+1,jj  ,jk-1) ) * zeps ),  &
+                             &            p_fld(ji+2,jj  ,jk) * ( 1. + ( 1. - p_mask(ji+2,jj  ,jk-1) ) * zeps ),  &
+                             &            p_fld(ji  ,jj+1,jk) * ( 1. + ( 1. - p_mask(ji  ,jj  ,jk-1) ) * zeps ),  &
+                             &            p_fld(ji+1,jj+1,jk) * ( 1. + ( 1. - p_mask(ji+1,jj+1,jk-1) ) * zeps ),  &
+                             &            p_fld(ji+2,jj+1,jk) * ( 1. + ( 1. - p_mask(ji+2,jj+2,jk-1) ) * zeps ),  &
+                             &            p_fld(ji  ,jj+2,jk) * ( 1. + ( 1. - p_mask(ji  ,jj  ,jk-1) ) * zeps ),  &
+                             &            p_fld(ji+1,jj+2,jk) * ( 1. + ( 1. - p_mask(ji+1,jj+1,jk-1) ) * zeps ),  &
+                             &            p_fld(ji+2,jj+2,jk) * ( 1. + ( 1. - p_mask(ji+2,jj+2,jk-1) ) * zeps )   )
+                           !
+                           p_fld_crs(ii,ij,jk) = zflcrs * p_mask_crs(ii,ij,jk)
+                           !
+                        ENDDO      
+                     ENDDO
+                  ENDDO  
+ 
+                  DO ji = nistr, niend, nn_factx
+                     DO jj = njstr, njend, nn_facty
+                        ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                        ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                        ijje = mje_crs(ij)
+                        ijie = mie_crs(ii)                  
+                        
+                        zflcrs =  MIN( p_fld(ji  ,jj  ,1) * ( 1. + ( 1. - p_mask(ji  ,jj  ,1) ) * zeps ),  &
+                          &            p_fld(ji+1,jj  ,1) * ( 1. + ( 1. - p_mask(ji+1,jj  ,1) ) * zeps ),  &
+                          &            p_fld(ji+2,jj  ,1) * ( 1. + ( 1. - p_mask(ji+2,jj  ,1) ) * zeps ),  &
+                          &            p_fld(ji  ,jj+1,1) * ( 1. + ( 1. - p_mask(ji  ,jj  ,1) ) * zeps ),  &
+                          &            p_fld(ji+1,jj+1,1) * ( 1. + ( 1. - p_mask(ji+1,jj+1,1) ) * zeps ),  &
+                          &            p_fld(ji+2,jj+1,1) * ( 1. + ( 1. - p_mask(ji+2,jj+2,1) ) * zeps ),  &
+                          &            p_fld(ji  ,jj+2,1) * ( 1. + ( 1. - p_mask(ji  ,jj  ,1) ) * zeps ),  &
+                          &            p_fld(ji+1,jj+2,1) * ( 1. + ( 1. - p_mask(ji+1,jj+1,1) ) * zeps ),  &
+                          &            p_fld(ji+2,jj+2,1) * ( 1. + ( 1. - p_mask(ji+2,jj+2,1) ) * zeps )   )
+                        !
+                        p_fld_crs(ii,ij,1) = zflcrs * p_mask_crs(ii,ij,1)
+                        !
+                     ENDDO      
+                  ENDDO
+
+               CASE( 'V' )
+         
+                  DO jk = 1, jpk
+                     DO ji = nistr, niend, nn_factx
+                        DO jj = njstr, njend, nn_facty
+                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                           ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                           ijje = mje_crs(ij)
+                           ijie = mie_crs(ii)                  
+                           
+                           zflcrs =  MIN( p_fld(ji  ,ijje,jk) * ( 1. + ( 1. - p_mask(ji  ,ijje,jk) ) * zeps ),  &
+                             &            p_fld(ji+1,ijje,jk) * ( 1. + ( 1. - p_mask(ji+1,ijje,jk) ) * zeps ),  &
+                             &            p_fld(ji+2,ijje,jk) * ( 1. + ( 1. - p_mask(ji+2,ijje,jk) ) * zeps )   )
+                           !
+                           p_fld_crs(ii,ij,jk) = zflcrs * p_mask_crs(ii,ij,jk)
+                           !
+                        ENDDO      
+                     ENDDO
+                  ENDDO  
+
+
+               CASE( 'U' )
+         
+                  DO jk = 1, jpk
+                     DO ji = nistr, niend, nn_factx
+                        DO jj = njstr, njend, nn_facty
+                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                           ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                           ijje = mje_crs(ij)
+                           ijie = mie_crs(ii)                  
+                           
+                           zflcrs =  MIN( p_fld(ijie,jj  ,jk) * ( 1. + ( 1. - p_mask(ijie,jj  ,jk) ) * zeps ),  &
+                             &            p_fld(ijie,jj+1,jk) * ( 1. + ( 1. - p_mask(ijie,jj+1,jk) ) * zeps ),  &
+                             &            p_fld(ijie,jj+2,jk) * ( 1. + ( 1. - p_mask(ijie,jj+2,jk) ) * zeps )   )
+                           !
+                           p_fld_crs(ii,ij,jk) = zflcrs * p_mask_crs(ii,ij,jk)
+                           !
+                        ENDDO      
+                     ENDDO
+                  ENDDO  
+            END SELECT
+            !
+         END SELECT
+         !
+         CALL crs_lbc_lnk( p_fld_crs, cd_type, 1.0 )
+         !
+    END SUBROUTINE crs_dom_ope_3d
+
+    SUBROUTINE crs_dom_ope_2d( p_fld, cd_op, cd_type, p_mask, p_fld_crs, p_e12, p_e3, p_surf_crs, p_mask_crs )
       !!----------------------------------------------------------------
-      !!               *** SUBROUTINE crsfun_TW ***
-      !! ** Purpose :  Five applications.
-      !!               1) Maximum surface quantity 
-      !!                  - Vertical scale factors (fse3t or fse3w) 
-      !!                    max thickness of the parent grid for coarse grid scale factors.
-      !!                  - or diffusion test
-      !!               2) Area-weighted mean quantity: w, or other 3D or 2D quantity
-      !!               3) Volume-weighted mean quantity: tracer
-      !!               4) Minimum surface quantity (diffusion test)
-      !!               5) Area- or Volume- weighted sum.
-      !! ** Method  :  1) - cd_op = 'MAX'. Determines the max vertical thickness of grid boxes
-      !!                    including partial steps for at the bottom
-      !!                    for the coarsened grid, where within the subset of 
-      !!                    the parent grid cells the maximum thickness is taken.
-      !!                    Valid arguments: p_e1e2t, cd_type, cd_op, p_cmask, p_pfield3d_1
-      !!                    Where, normally p_pfield3d_1 is e3t.
-      !!                  - cd_op = 'MAX'. May also be used for say, determining the maximum value of Kz, 
-      !!                    thus p_pfield3d_1 is set to the 3D field, Kz.
-      !!                    Valid arguments: p_e1e2t, cd_type, cd_op, p_cmask, p_pfield3d_1
-      !!               2) - cd_op = 'ARE'. Calculate the area-weighted average (surface e1t*e2t)  
-      !!                    of vertical velocity, w.
-      !!                    Valid arguments: p_e1e2t, cd_type, cd_op, p_cmask, p_pfield3d_1
-      !!                  - cd_op = 'ARE'. Calculate area-weighted average of a 2D quantity (e.g. emp)
-      !!                    Valid arguments: p_e1e2t, cd_type, cd_op, p_cmask, p_pfield2d
-      !!               3) - cd_op = 'VOL'. Calculate the ocean volume (e1e2t*[fse3t|fse3w]) 
-      !!                    Valid arguments: p_e1e2t, cd_type, cd_op, p_cmask, p_pfield3d_1
-      !!                  - cd_op = 'VOL'. Calculate volume-weighted average (volume e1t*e2t*fse3t) of a quantity.
-      !!                    Valid arguments: p_e1e2t, cd_type, cd_op, p_cmask, p_pfield3d_1, p_pfield3d_2
-      !!               4) - cd_op = 'MIN'. Calculate the minimum value on surface e1t*e2t for 3D variables
-      !!                  Valid arguments: p_e1e2t, cd_type, cd_op, p_cmask, p_pfield3d_1
-      !!               5) - cd_op = 'SUM'. Calculate a dimesionally-weighted sum.  This could be area-weighted
-      !!                  or volume-weighted sum. 
-      !! ** Inputs  : p_e1e2t      = parent grid top face surface area, e1t*e2t 
-      !!              cd_type      = grid type T, W (U, V, F) 
-      !!              cd_op        = MAX, ARE, VOL, MIN, SUM
-      !!              p_cmask      =  coarse grid mask
-      !!              p_ptmask     =  parent grid tmask     
-      !!              psgn         = (Optional) sign for lbc_lnk  
-      !!              p_pfield2d   = (Optional) 2D field on parent grid
-      !!              p_pfield3d_1 = (Optional) parent grid fse3t or fse3w
-      !!              p_pfield3d_2 = (Optional) 3D field on parent grid
-      !! ** Outputs : p_cfield2d   = (Optional) 2D field on coarse grid
-      !!              p_cfield3d   = (Optional) 3D field on coarse grid
+      !!               *** SUBROUTINE crsfun_UV ***
+      !! ** Purpose :  Average, area-weighted, of U or V on the east and north faces 
       !!
-      !! 
-      !! History.  30 May.  Editing.  To decide later: Keep all functionality or separate out the mean function.
-      !!            7 Jun   TODO. Need to fix up the parent grid mask to optional like crsfun_UV...
+      !! ** Method  :  The U and V velocities (3D) are determined as the area-weighted averages
+      !!               on the east and north faces, respectively,
+      !!               of the parent grid subset comprising the coarse grid box. 
+      !!               In the case of the V and F grid, the last jrow minus 1 is spurious.
+      !! ** Inputs  : p_e1_e2     = parent grid e1 or e2 (t,u,v,f)
+      !!              cd_type     = grid type (T,U,V,F) for scale factors; for velocities (U or V)
+      !!              psgn        = sign change over north fold (See lbclnk.F90)
+      !!              p_pmask     = parent grid mask (T,U,V,F) for scale factors; 
+      !!                                       for velocities (U or V)
+      !!              p_fse3      = parent grid vertical level thickness (fse3u or fse3v)
+      !!              p_pfield    = U or V on the parent grid
+      !!              p_surf_crs  = (Optional) Coarse grid weight for averaging
+      !! ** Outputs : p_cfield3d = 3D field on coarse grid
+      !!
+      !! History.  29 May.  completed draft.
+      !!            4 Jun.  Revision for WGT
+      !!            5 Jun.  Streamline for area-weighted average only ; separate scale factor and weights.
       !!----------------------------------------------------------------
       !! 
       !!  Arguments
-      REAL(wp), DIMENSION(jpi,jpj),               INTENT(in) :: p_e1e2t      ! Parent grid T surface (e1*e2)
-      CHARACTER(len=1),                           INTENT(in) :: cd_type      ! grid type T, W ( U, V, F)
-      CHARACTER(len=3),                           INTENT(in) :: cd_op        ! operation max, min, area-average, volume-average
-      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk),   INTENT(in) :: p_cmask      ! Coarse grid T mask
-      REAL(wp), DIMENSION(jpi,jpj,jpk),           INTENT(in) :: p_ptmask     ! Parent grid T mask
-      REAL(wp),                         OPTIONAL, INTENT(in) :: psgn         ! sign for lbc_lnk
-      REAL(wp), DIMENSION(jpi,jpj),     OPTIONAL, INTENT(in) :: p_pfield2d   ! 2D quantity on parent grid, e.g. ssh
-      REAL(wp), DIMENSION(jpi,jpj,jpk), OPTIONAL, INTENT(in) :: p_pfield3d_1 ! Normally parent grid vertical level thickness
-      REAL(wp), DIMENSION(jpi,jpj,jpk), OPTIONAL, INTENT(in) :: p_pfield3d_2 ! 3D tracer or W on parent grid
-
-      REAL(wp), DIMENSION(jpi_crs,jpj_crs),     OPTIONAL, INTENT(out):: p_cfield2d ! Coarse grid box east or north face quantity
-      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), OPTIONAL, INTENT(out):: p_cfield3d ! Coarse grid box east or north face quantity 
+      REAL(wp), DIMENSION(jpi,jpj),             INTENT(in)           :: p_fld   ! T, U, V or W on parent grid
+      CHARACTER(len=3),                         INTENT(in)           :: cd_op    ! Operation SUM, MAX or MIN
+      CHARACTER(len=1),                         INTENT(in)           :: cd_type    ! grid type U,V 
+      REAL(wp), DIMENSION(jpi,jpj,jpk),         INTENT(in)           :: p_mask    ! Parent grid T,U,V mask
+      REAL(wp), DIMENSION(jpi,jpj),             INTENT(in), OPTIONAL :: p_e12    ! Parent grid T,U,V scale factors (e1 or e2)
+      REAL(wp), DIMENSION(jpi,jpj,jpk),         INTENT(in), OPTIONAL :: p_e3     ! Parent grid vertical level thickness (fse3u, fse3v)
+      REAL(wp), DIMENSION(jpi_crs,jpj_crs)    , INTENT(in), OPTIONAL :: p_surf_crs ! Coarse grid area-weighting denominator    
+      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(in), OPTIONAL :: p_mask_crs    ! Coarse grid T,U,V mask
+
+      REAL(wp), DIMENSION(jpi_crs,jpj_crs)    , INTENT(out)          :: p_fld_crs ! Coarse grid box 3D quantity 
 
       !! Local variables
-      INTEGER ::  ji, jj, jk                   ! dummy loop indices
-      INTEGER :: ijie,ijis,ijje,ijjs,ijpk,jii,jjj
-      INTEGER, DIMENSION(3) :: idims
-      REAL(wp), POINTER, DIMENSION(:,:)   :: ze1e2, zpfield2d, zcfield2d
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3, zpfield3d, zcfield3d, zcmask, zpmask  
-      REAL(wp)                            :: zdAm, zsgn
+      INTEGER  :: ji, jj, jk                 ! dummy loop indices
+      INTEGER  :: ijie, ijje, ii, ij
+      REAL(wp) :: zflcrs, zsfcrs   
+      REAL(wp) :: zeps = 1.e20    
+      REAL(wp), DIMENSION(:,:), POINTER :: zsurf   
+
       !!----------------------------------------------------------------  
-      ! Initialize
-
-      CALL wrk_alloc(jpi    , jpj         , ze1e2, zpfield2d )
-      CALL wrk_alloc(jpi    , jpj    , jpk,  ze3 , zpfield3d, zpmask )
-      CALL wrk_alloc(jpi_crs, jpj_crs,      zcfield2d )
-      CALL wrk_alloc(jpi_crs, jpj_crs, jpk, zcfield3d,  zcmask )
-
-
-      ! Arrays, scalars initialization 
-      zpfield2d(:,:)   = 0.0
-      zpfield3d(:,:,:) = 0.0
-      zcfield2d(:,:)   = 0.0
-      zcfield3d(:,:,:) = 0.0
-      zpmask(:,:,:)    = 1.0
-      idims(:)         = 1
-
-      zcmask(:,:,:)    = p_cmask(:,:,:)     
-      zpmask(:,:,:)    = p_ptmask(:,:,:)
-
-      ijpk             = jpk
-
-
-      ! Control of optional arguments
-      !
-      IF ( PRESENT(psgn) ) THEN
-         zsgn = psgn
-      ELSE
-         zsgn = 1.0 
-      ENDIF
-      !
-      IF ( TRIM(cd_op) == 'MAX' ) THEN
-         ! Find the maximum thickness in each parent grid subset
-         IF ( PRESENT(p_pfield3d_1) )  THEN
-            zpfield3d(:,:,:) = p_pfield3d_1(:,:,:)
-            ze3(:,:,:) = 0.0
-            ze1e2(:,:) = 0.0
-         ELSE 
-            WRITE(numout,*) 'crsfun_TW. MAX only 3D arrays supported' 
-         ENDIF 
-      ELSEIF ( TRIM(cd_op) == 'VOL' ) THEN
-         IF ( PRESENT(p_pfield3d_1) ) THEN
-            ze3(:,:,:)       = p_pfield3d_1(:,:,:)
-            IF ( PRESENT(p_pfield3d_2) ) THEN
-!              ! Prep to calculate a volume-averaged mean
-               zpfield3d(:,:,:) = p_pfield3d_2(:,:,:)
-               ze1e2(:,:)       = p_e1e2t(:,:)     
-            ELSE
-               WRITE(numout,*) 'crsfun_TW. WARNING. Supply both e3t and the field for volume-averaged field'
-            ENDIF
-         ELSE
-            WRITE(numout,*) 'crsfun_TW. VOL only 3D arrays supported'  
-         ENDIF
-      ELSEIF ( TRIM(cd_op) == 'ARE' ) THEN
-         ze1e2(:,:) = p_e1e2t(:,:)
-         ze3(:,:,:) = 1.0
-         IF ( PRESENT(p_pfield3d_2) ) THEN
-            ! Prep to do the area-weighted average of (3D) W
-            zpfield3d(:,:,:) = p_pfield3d_2(:,:,:)
-         ELSEIF ( PRESENT(p_pfield2d) ) THEN
-            ! Prep to do the area-weighted average of some 2D quantity 
-            zpfield2d(:,:) = p_pfield2d(:,:) 
-            ijpk=1
-         ENDIF
-      ELSEIF ( TRIM(cd_op) == 'MIN' ) THEN
-         IF ( PRESENT(p_pfield3d_1) ) THEN
-            ! Prep to do get the minimum diffusion on the top face
-            zpfield3d(:,:,:) = p_pfield3d_1(:,:,:)
-            ze3(:,:,:) = 0.0
-            ze1e2(:,:) = 0.0
-         ELSE
-            WRITE(numout,*) 'crsfun_TW. MIN Only 3D arrays supported' 
-         ENDIF
-      ELSEIF ( TRIM(cd_op) == 'SUM' ) THEN
-         ze1e2(:,:) = p_e1e2t(:,:)
-         zpmask(:,:,:) = p_ptmask(:,:,:)
-         ze3(:,:,:) = 1.0
-        IF ( PRESENT(p_pfield3d_1) ) THEN
-            IF ( PRESENT(p_pfield3d_2) ) THEN
-!              ! Prep to calculate a volume-weighted sum
-               zpfield3d(:,:,:) = p_pfield3d_2(:,:,:)
-               ze3(:,:,:) = p_pfield3d_1(:,:,:)
-            ELSE
-               ! Prep to do the area-weighted sum of (3D) W
-               zpfield3d(:,:,:) = p_pfield3d_1(:,:,:)
-            ENDIF
-         ELSEIF ( PRESENT(p_pfield2d) ) THEN
-            ! Prep to do the area-weighted sum of some 2D quantity 
-            zpfield2d(:,:) = p_pfield2d(:,:) 
-            ijpk=1
-         ENDIF 
-      ELSE
-         WRITE(numout,*) 'crsfun_TW. valid cd_op are MAX, MIN, ARE, VOL, SUM' 
-      ENDIF
-
-      ! Determine output 
-      DO jk = 1, ijpk    
-
-         IF ( ijpk == jpk ) zpfield2d(:,:) = zpfield3d(:,:,jk) 
-         zcfield2d(:,:) = 0.0 
-
-           DO ji = 2, nlei_crs
-               ijie = mie_crs(ji)
-               ijis = mis_crs(ji)
-
-             DO jj = njstart, njend
-                  ijje = mje_crs(jj) 
-                  ijjs = mjs_crs(jj)                   
-
-                  ! First determine weighted sums 
-                  IF ( TRIM(cd_op) == 'SUM' .OR. TRIM(cd_op) == 'ARE' .OR. TRIM(cd_op) == 'VOL' ) THEN
-                     ! Area- or volume- weighted sum
-                     ! Accumulate to sum in parent grid subset 
-                    DO jii = ijis, ijie
-                        DO jjj = ijjs, ijje
-                           zcfield2d(ji,jj) = zcfield2d(ji,jj) &
-                              &              + ( zpfield2d(jii,jjj)      &
-                              &                  *   ze1e2(jii,jjj)      &
-                              &                  *     ze3(jii,jjj,jk)   &
-                              &                  *  zpmask(jii,jjj,jk) ) 
-
-                        ENDDO
+   
+
+      SELECT CASE ( cd_op )
+      
+         CASE ( 'VOL' )
+      
+           CALL wrk_alloc( jpi, jpj, zsurf )
+           zsurf(:,:) =  p_e12(:,:) * p_e3(:,:,1) * p_mask(:,:,1)
+                    
+           DO ji = nistr, niend, nn_factx
+              DO jj = njstr, njend, nn_facty
+                 ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                 ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                 ijje = mje_crs(ij)
+                 ijie = mie_crs(ii)                  
+ 
+                 zflcrs =  p_fld(ji  ,jj  ) * zsurf(ji  ,jj  )  &
+                   &     + p_fld(ji+1,jj  ) * zsurf(ji+1,jj  )  &
+                   &     + p_fld(ji+2,jj  ) * zsurf(ji+2,jj  )  &
+                   &     + p_fld(ji  ,jj+1) * zsurf(ji  ,jj+1)  &
+                   &     + p_fld(ji+1,jj+1) * zsurf(ji+1,jj+1)  &
+                   &     + p_fld(ji+2,jj+1) * zsurf(ji+2,jj+1)  &
+                   &     + p_fld(ji  ,jj+2) * zsurf(ji  ,jj+2)  &
+                   &     + p_fld(ji+1,jj+2) * zsurf(ji+1,jj+2)  &
+                   &     + p_fld(ji+2,jj+2) * zsurf(ji+2,jj+2)  
+
+                 zsfcrs =  zsurf(ji,jj  ) + zsurf(ji+1,jj  ) + zsurf(ji+2,jj  )  &
+                   &     + zsurf(ji,jj+1) + zsurf(ji+1,jj+1) + zsurf(ji+2,jj+1)  &
+                   &     + zsurf(ji,jj+2) + zsurf(ji+1,jj+2) + zsurf(ji+2,jj+2)  
+                 !
+                 p_fld_crs(ii,ij) = zflcrs
+                 IF( zsfcrs /= 0.0 )  p_fld_crs(ii,ij) = zflcrs / zsfcrs
+
+              ENDDO      
+           ENDDO
+            
+           CALL wrk_dealloc( jpi, jpj, zsurf )
+
+         CASE ( 'SUM' )
+         
+            CALL wrk_alloc( jpi, jpj, zsurf )
+            zsurf(:,:) =  p_e12(:,:) * p_e3(:,:,1) * p_mask(:,:,1)
+         
+            SELECT CASE ( cd_type )
+            
+               CASE( 'T', 'W' )
+         
+                  DO ji = nistr, niend, nn_factx
+                     DO jj = njstr, njend, nn_facty
+                        ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                        ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                        ijje = mje_crs(ij)
+                        ijie = mie_crs(ii)                  
+
+                        zflcrs  =  p_fld(ji  ,jj  ) * zsurf(ji  ,jj  )  &
+                          &      + p_fld(ji+1,jj  ) * zsurf(ji+1,jj  )  &
+                          &      + p_fld(ji+2,jj  ) * zsurf(ji+2,jj  )  &
+                          &      + p_fld(ji  ,jj+1) * zsurf(ji  ,jj+1)  &
+                          &      + p_fld(ji+1,jj+1) * zsurf(ji+1,jj+1)  &
+                          &      + p_fld(ji+2,jj+1) * zsurf(ji+2,jj+1)  &
+                          &      + p_fld(ji  ,jj+2) * zsurf(ji  ,jj+2)  &
+                          &      + p_fld(ji+1,jj+2) * zsurf(ji+1,jj+2)  &
+                          &      + p_fld(ji+2,jj+2) * zsurf(ji+2,jj+2)  
+                        !
+                        p_fld_crs(ii,ij) = zflcrs
+                        !
+                     ENDDO      
+                  ENDDO
+            
+               CASE( 'V' )
+         
+                  DO jk = 1, jpk
+                     DO ji = nistr, niend, nn_factx
+                        DO jj = njstr, njend, nn_facty
+                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                           ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                           ijje = mje_crs(ij)
+                           ijie = mie_crs(ii)
+                           !                  
+                           zflcrs  =  p_fld(ji  ,ijje) * zsurf(ji  ,ijje)  &
+                             &      + p_fld(ji+1,ijje) * zsurf(ji+1,ijje)  &
+                             &      + p_fld(ji+2,ijje) * zsurf(ji+2,ijje)  
+                           !
+                           p_fld_crs(ii,ij) = zflcrs
+                           !
+                         ENDDO      
                      ENDDO
-
-                  ENDIF
-
-                  ! Calculate weighted average if desired
-                  IF ( TRIM(cd_op) == 'ARE' .OR. TRIM(cd_op) == 'VOL' ) THEN
-
-                     ! Area- or volume- weighted mean
-                     ! Sum first parent grid subset  
-                     zdAm = 0.0
-                     DO jii = ijis, ijie
-                        DO jjj = ijjs, ijje
-                           zdAm = zdAm + (   ze1e2(jii,jjj)               &
-                              &           *    ze3(jii,jjj,jk)            &
-                              &           * zpmask(jii,jjj,jk) )
-                        ENDDO
+                  ENDDO   
+
+            
+               CASE( 'U' )
+         
+                  DO jk = 1, jpk
+                     DO ji = nistr, niend, nn_factx
+                        DO jj = njstr, njend, nn_facty
+                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                           ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                           ijje = mje_crs(ij)
+                           ijie = mie_crs(ii)
+                           !                  
+                           zflcrs =  p_fld(ijie,jj  ) * zsurf(ijie,jj  )  &
+                             &     + p_fld(ijie,jj+1) * zsurf(ijie,jj+1)  &
+                             &     + p_fld(ijie,jj+2) * zsurf(ijie,jj+2) 
+                           !
+                           p_fld_crs(ii,ij) = zflcrs
+                           !
+                        ENDDO      
                      ENDDO
-                     
-                     IF ( zdAm /= 0 )  zcfield2d(ji,jj) = zcfield2d(ji,jj) / zdAm
-
-                  ENDIF
-
-
-                  IF ( TRIM(cd_op) == 'MAX' ) THEN
-                     ! Find max in parent grid subset 
-                     DO jii = ijis, ijie
-                        DO jjj = ijjs, ijje
-                           zcfield2d(ji,jj) = MAX( zcfield2d(ji,jj), zpfield3d(jii,jjj,jk)*zpmask(jii,jjj,jk) ) 
-                        ENDDO
-                     ENDDO
-                  ENDIF
-
-                  IF ( TRIM(cd_op) == 'MIN' ) THEN
-                     ! Find min in parent grid subset  
-                     DO jii = ijis, ijie
-                        DO jjj = ijjs, ijje
-                           IF ( zpmask(jii,jjj,jk) == 1 ) THEN
-                              zcfield2d(ji,jj) = MIN( zcfield2d(ji,jj), zpfield3d(jii,jjj,jk) ) 
-                           ENDIF
-                        ENDDO
-                     ENDDO
-                  ENDIF
-
-               ENDDO
-            ENDDO
-
-            IF ( ijpk == 1 ) THEN
-               IF ( PRESENT(p_cfield2d) ) p_cfield2d(:,:) = zcfield2d(:,:) * zcmask(:,:,jk)
-            ELSE
-               IF ( PRESENT(p_cfield3d) ) p_cfield3d(:,:,jk) = zcfield2d(:,:) * zcmask(:,:,jk)
-            ENDIF    
-
-      ENDDO
-
-
-! Retroactively add back the boundary halo cells.
-
-         IF ( ijpk == 1 ) THEN
-            IF ( PRESENT(p_cfield2d) ) CALL crs_lbc_lnk( p_cfield2d(:,:),cd_type,zsgn )
-         ELSE
-            IF ( PRESENT(p_cfield3d) ) THEN
-               CALL crs_lbc_lnk( p_cfield3d(:,:,:),cd_type,zsgn )
-            ENDIF
-         ENDIF
-     
-      CALL wrk_dealloc(jpi    , jpj         , ze1e2, zpfield2d )
-      CALL wrk_dealloc(jpi    , jpj    , jpk,  ze3 , zpfield3d, zpmask )
-      CALL wrk_dealloc(jpi_crs, jpj_crs,      zcfield2d )
-      CALL wrk_dealloc(jpi_crs, jpj_crs, jpk, zcfield3d,  zcmask )
+                  ENDDO   
+
+              END SELECT
+
+              IF( PRESENT( p_surf_crs ) ) THEN
+                 WHERE ( p_surf_crs /= 0.0 ) p_fld_crs(:,:) = p_fld_crs(:,:) / p_surf_crs(:,:)
+              ENDIF
+
+              CALL wrk_dealloc( jpi, jpj, zsurf )
+
+         CASE ( 'MAX' )
+         
+            SELECT CASE ( cd_type )
+            
+               CASE( 'T', 'W' )
+         
+                  DO ji = nistr, niend, nn_factx
+                     DO jj = njstr, njend, nn_facty
+                        ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                        ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                        ijje = mje_crs(ij)
+                        ijie = mie_crs(ii)                  
+
+                        zflcrs =  MAX( p_fld(ji  ,jj  ) * p_mask(ji  ,jj  ,1),  &
+                          &            p_fld(ji+1,jj  ) * p_mask(ji+1,jj  ,1),  &
+                          &            p_fld(ji+2,jj  ) * p_mask(ji+2,jj  ,1),  &
+                          &            p_fld(ji  ,jj+1) * p_mask(ji  ,jj+1,1),  &
+                          &            p_fld(ji+1,jj+1) * p_mask(ji+1,jj+1,1),  &
+                          &            p_fld(ji+2,jj+1) * p_mask(ji+2,jj+1,1),  &
+                          &            p_fld(ji  ,jj+2) * p_mask(ji  ,jj+2,1),  &
+                          &            p_fld(ji+1,jj+2) * p_mask(ji+1,jj+2,1),  &
+                          &            p_fld(ji+2,jj+2) * p_mask(ji+2,jj+2,1)   )
+                        !
+                        p_fld_crs(ii,ij) = zflcrs
+                        !
+                     ENDDO      
+                  ENDDO
+            
+               CASE( 'V' )
+         
+                  DO ji = nistr, niend, nn_factx
+                     DO jj = njstr, njend, nn_facty
+                        ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                        ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                        ijje = mje_crs(ij)
+                        ijie = mie_crs(ii)
+                        !                  
+                        zflcrs = MAX( p_fld(ji  ,ijje) * p_mask(ji  ,ijje,1),  &
+                          &           p_fld(ji+1,ijje) * p_mask(ji+1,ijje,1),  &
+                          &           p_fld(ji+2,ijje) * p_mask(ji+2,ijje,1) )
+                        !
+                        p_fld_crs(ii,ij) = zflcrs
+                        !
+                     ENDDO      
+                  ENDDO
+            
+               CASE( 'U' )
+         
+                  DO ji = nistr, niend, nn_factx
+                     DO jj = njstr, njend, nn_facty
+                        ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                        ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                        ijje = mje_crs(ij)
+                        ijie = mie_crs(ii)
+                        !                  
+                        zflcrs = MAX( p_fld(ijie,jj  ) * p_mask(ijie,jj  ,1),  &
+                          &           p_fld(ijie,jj+1) * p_mask(ijie,jj+1,1),  &
+                          &           p_fld(ijie,jj+2) * p_mask(ijie,jj+2,1) )
+                        !
+                        p_fld_crs(ii,ij) = zflcrs
+                        !
+                     ENDDO      
+                  ENDDO
+
+              END SELECT
+
+         CASE ( 'MIN' )
+            !   Search the min of masked grid cells
+            SELECT CASE ( cd_type )
+            
+               CASE( 'T', 'W' )
+         
+                  DO ji = nistr, niend, nn_factx
+                     DO jj = njstr, njend, nn_facty
+                        ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                        ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                        ijje = mje_crs(ij)
+                        ijie = mie_crs(ii)                  
+                        
+                        zflcrs =  MIN( p_fld(ji  ,jj  ) * ( 1. + ( 1. - p_mask(ji  ,jj  ,1) ) * zeps ),  &
+                          &            p_fld(ji+1,jj  ) * ( 1. + ( 1. - p_mask(ji+1,jj  ,1) ) * zeps ),  &
+                          &            p_fld(ji+2,jj  ) * ( 1. + ( 1. - p_mask(ji+2,jj  ,1) ) * zeps ),  &
+                          &            p_fld(ji  ,jj+1) * ( 1. + ( 1. - p_mask(ji  ,jj  ,1) ) * zeps ),  &
+                          &            p_fld(ji+1,jj+1) * ( 1. + ( 1. - p_mask(ji+1,jj+1,1) ) * zeps ),  &
+                          &            p_fld(ji+2,jj+1) * ( 1. + ( 1. - p_mask(ji+2,jj+2,1) ) * zeps ),  &
+                          &            p_fld(ji  ,jj+2) * ( 1. + ( 1. - p_mask(ji  ,jj  ,1) ) * zeps ),  &
+                          &            p_fld(ji+1,jj+2) * ( 1. + ( 1. - p_mask(ji+1,jj+1,1) ) * zeps ),  &
+                          &            p_fld(ji+2,jj+2) * ( 1. + ( 1. - p_mask(ji+2,jj+2,1) ) * zeps )   )
+                        !
+                        p_fld_crs(ii,ij) = zflcrs * p_mask_crs(ii,ij,1)
+                        !
+                     ENDDO      
+                  ENDDO
  
-
-   END SUBROUTINE crsfun_TW
+               CASE( 'V' )
+         
+                  DO ji = nistr, niend, nn_factx
+                     DO jj = njstr, njend, nn_facty
+                        ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                        ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                        ijje = mje_crs(ij)
+                        ijie = mie_crs(ii)                  
+                        
+                        zflcrs =  MIN( p_fld(ji  ,ijje) * ( 1. + ( 1. - p_mask(ji  ,ijje,1) ) * zeps ),  &
+                          &            p_fld(ji+1,ijje) * ( 1. + ( 1. - p_mask(ji+1,ijje,1) ) * zeps ),  &
+                          &            p_fld(ji+2,ijje) * ( 1. + ( 1. - p_mask(ji+2,ijje,1) ) * zeps )   )
+                        !
+                        p_fld_crs(ii,ij) = zflcrs * p_mask_crs(ii,ij,1)
+                        !
+                     ENDDO      
+                  ENDDO
+
+               CASE( 'U' )
+         
+                  DO ji = nistr, niend, nn_factx
+                     DO jj = njstr, njend, nn_facty
+                        ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                        ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                        ijje = mje_crs(ij)
+                        ijie = mie_crs(ii)                  
+                        
+                        zflcrs =  MIN( p_fld(ijie,jj  ) * ( 1. + ( 1. - p_mask(ijie,jj  ,1) ) * zeps ),  &
+                          &            p_fld(ijie,jj+1) * ( 1. + ( 1. - p_mask(ijie,jj+1,1) ) * zeps ),  &
+                          &            p_fld(ijie,jj+2) * ( 1. + ( 1. - p_mask(ijie,jj+2,1) ) * zeps )   )
+                        !
+                        p_fld_crs(ii,ij) = zflcrs * p_mask_crs(ii,ij,1)
+                        !
+                     ENDDO      
+                  ENDDO
+            END SELECT
+            !
+         END SELECT
+         !
+         CALL crs_lbc_lnk( p_fld_crs, cd_type, 1.0 )
+         !
+   END SUBROUTINE crs_dom_ope_2d
 
 
@@ -914,5 +1304,6 @@
       !! Local variables
       INTEGER ::  ji, jj, jk                   ! dummy loop indices
-      INTEGER :: ijie,ijis,ijje,ijjs,jii,jjj
+      INTEGER ::  ijie, ijje, ii, ij
+      REAL(wp) :: ze3crs  
       !!----------------------------------------------------------------  
       ! Initialize
@@ -923,19 +1314,23 @@
          
             DO jk = 1 , jpk
-              
-               DO ji = 2, nlei_crs        ! ji = 1 et jpi_crs definit par cyclique est-ouest et pivot T
-                   ijie = mie_crs(ji)
-                   ijis = mis_crs(ji)
-
-                   DO jj =  njstart, njend  ! jj = jpj_crs definit par pivot T 
-                      ijje = mje_crs(jj) 
-                      ijjs = mjs_crs(jj)  
-                      
-                      DO jii = ijis, ijie
-                         DO jjj = ijjs, ijje
-                            p_e3_crs(ji,jj,jk) = max( p_e3_crs(ji,jj,jk), p_e3(jii,jjj,jk) * p_mask(jii,jjj,jk)  )
-                         ENDDO
-                      ENDDO
-                   ENDDO
+               DO ji = nistr, niend, nn_factx
+                  DO jj = njstr, njend, nn_facty
+                     ii   = ( ji - mis_crs(2) ) * rfactx_r + 2                 ! cordinate in parent grid
+                     ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                     ijje = mje_crs(ij)
+                     ijie = mie_crs(ii)
+                     !  
+                     ze3crs = MAX( p_e3(ji  ,jj  ,jk) * p_mask(ji  ,jj  ,jk),  &
+                        &          p_e3(ji+1,jj  ,jk) * p_mask(ji+1,jj  ,jk),  &
+                        &          p_e3(ji+2,jj  ,jk) * p_mask(ji+2,jj  ,jk),  &
+                        &          p_e3(ji  ,jj+1,jk) * p_mask(ji  ,jj+1,jk),  &
+                        &          p_e3(ji+1,jj+1,jk) * p_mask(ji+1,jj+1,jk),  &
+                        &          p_e3(ji+2,jj+1,jk) * p_mask(ji+2,jj+1,jk),  &
+                        &          p_e3(ji  ,jj+2,jk) * p_mask(ji  ,jj+2,jk),  &
+                        &          p_e3(ji+1,jj+2,jk) * p_mask(ji+1,jj+2,jk),  &
+                        &          p_e3(ji+2,jj+2,jk) * p_mask(ji+2,jj+2,jk) )
+                    
+                     p_e3_crs(ii,ij,jk) = ze3crs
+                  ENDDO
                ENDDO
             ENDDO
@@ -944,159 +1339,200 @@
          
             DO jk = 2 , jpk
-            
-                                          ! ji = 1 et jpi_crs definit par cyclique est-ouest et pivot T
-               DO ji = 2, nlei_crs
-                   ijie = mie_crs(ji)
-                   ijis = mis_crs(ji)
-
-                   DO jj = njstart, njend ! jj = jpj_crs definit par pivot T 
-                      ijje = mje_crs(jj) 
-                      ijjs = mjs_crs(jj)  
-                   
-                      DO jii = ijis, ijie
-                         DO jjj = ijjs, ijje
-                            p_e3_crs(ji,jj,jk) = max( p_e3_crs(ji,jj,jk), p_e3(jii,jjj,jk) * p_mask(jii,jjj,jk-1)  )
-                         ENDDO
-                      ENDDO
-                   ENDDO
-               ENDDO
-            ENDDO
-             
-            jk = 1                                           ! cas particulier car zpmask(jii,jjj,0) n'existe pas      
-          
-            DO ji = 2, nlei_crs
-               ijie = mie_crs(ji)
-               ijis = mis_crs(ji)
-
-               DO jj = njstart, njend            
-                  ijje = mje_crs(jj) 
-                  ijjs = mjs_crs(jj)  
-                   
-                  DO jii = ijis, ijie
-                     DO jjj = ijjs, ijje
-                        p_e3_crs(ji,jj,jk) = max( p_e3_crs(ji,jj,jk), p_e3(jii,jjj,jk) * p_mask(jii,jjj,jk)  )
-                     ENDDO
+               DO ji = nistr, niend, nn_factx
+                  DO jj = njstr, njend, nn_facty
+                     ii   = ( ji - mis_crs(2) ) * rfactx_r + 2                 ! cordinate in parent grid
+                     ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                     ijje = mje_crs(ij)
+                     ijie = mie_crs(ii)
+                     !  
+                     ze3crs = MAX( p_e3(ji  ,jj  ,jk) * p_mask(ji  ,jj  ,jk-1),  &
+                        &          p_e3(ji+1,jj  ,jk) * p_mask(ji+1,jj  ,jk-1),  &
+                        &          p_e3(ji+2,jj  ,jk) * p_mask(ji+2,jj  ,jk-1),  &
+                        &          p_e3(ji  ,jj+1,jk) * p_mask(ji  ,jj+1,jk-1),  &
+                        &          p_e3(ji+1,jj+1,jk) * p_mask(ji+1,jj+1,jk-1),  &
+                        &          p_e3(ji+2,jj+1,jk) * p_mask(ji+2,jj+1,jk-1),  &
+                        &          p_e3(ji  ,jj+2,jk) * p_mask(ji  ,jj+2,jk-1),  &
+                        &          p_e3(ji+1,jj+2,jk) * p_mask(ji+1,jj+2,jk-1),  &
+                        &          p_e3(ji+2,jj+2,jk) * p_mask(ji+2,jj+2,jk-1) )
+                    
+                     p_e3_crs(ii,ij,jk) = ze3crs
                   ENDDO
                ENDDO
-            ENDDO  
+            ENDDO
+             
+            DO ji = nistr, niend, nn_factx
+               DO jj = njstr, njend, nn_facty
+                  ii   = ( ji - mis_crs(2) ) * rfactx_r + 2                 ! cordinate in parent grid
+                  ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                  ijje = mje_crs(ij)
+                  ijie = mie_crs(ii)
+                  !  
+                  ze3crs = MAX( p_e3(ji  ,jj  ,1) * p_mask(ji  ,jj  ,1),  &
+                     &          p_e3(ji+1,jj  ,1) * p_mask(ji+1,jj  ,1),  &
+                     &          p_e3(ji+2,jj  ,1) * p_mask(ji+2,jj  ,1),  &
+                     &          p_e3(ji  ,jj+1,1) * p_mask(ji  ,jj+1,1),  &
+                     &          p_e3(ji+1,jj+1,1) * p_mask(ji+1,jj+1,1),  &
+                     &          p_e3(ji+2,jj+1,1) * p_mask(ji+2,jj+1,1),  &
+                     &          p_e3(ji  ,jj+2,1) * p_mask(ji  ,jj+2,1),  &
+                     &          p_e3(ji+1,jj+2,1) * p_mask(ji+1,jj+2,1),  &
+                     &          p_e3(ji+2,jj+2,1) * p_mask(ji+2,jj+2,1) )
+                 
+                  p_e3_crs(ii,ij,1) = ze3crs
+               ENDDO
+            ENDDO
             
          END SELECT 
-                  
-         CALL crs_lbc_lnk( p_e3_crs(:,:,:),cd_type, 1.0, pval=1.0 )  
-     
-         WRITE(numout,*) 'crs_e3_max : end of subroutine '
- 
-
+         !              
+         CALL crs_lbc_lnk( p_e3_crs, cd_type, 1.0, pval=1.0 )  
+         !              
    END SUBROUTINE crs_dom_e3_max
 
-
-   SUBROUTINE crs_dom_sfc( p_e1, p_e2, p_e3, cd_type, p_mask, p_surf_crs, p_surf_crs_msk )
+   SUBROUTINE crs_dom_sfc( p_mask, cd_type, p_surf_crs, p_surf_crs_msk,  p_e1, p_e2, p_e3 )
 
       !!  Arguments
       CHARACTER(len=1),                 INTENT(in) :: cd_type      ! grid type T, W ( U, V, F)
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: p_mask       ! Parent grid T mask
-      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in) :: p_e1, p_e2         ! 3D tracer T or W on parent grid
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: p_e3         ! 3D tracer T or W on parent grid
-      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(inout):: p_surf_crs ! Coarse grid box east or north face quantity 
-      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(inout):: p_surf_crs_msk ! Coarse grid box east or north face quantity 
+      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in), OPTIONAL :: p_e1, p_e2         ! 3D tracer T or W on parent grid
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL :: p_e3         ! 3D tracer T or W on parent grid
+      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(out):: p_surf_crs ! Coarse grid box east or north face quantity 
+      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(out):: p_surf_crs_msk ! Coarse grid box east or north face quantity 
 
       !! Local variables
-      INTEGER  ::  ji, jj, jk                   ! dummy loop indices
-      INTEGER  :: ijie,ijis,ijje,ijjs,ijpk,jii,jjj
-      REAL(wp) :: zsfc
+      INTEGER  :: ji, jj, jk                   ! dummy loop indices
+      INTEGER  :: ijie, ijje, ii, ij
+      REAL(wp), DIMENSION(:,:)  , POINTER :: zsurf   
+      REAL(wp), DIMENSION(:,:,:), POINTER :: zsurf3d   
+      REAL(wp) :: zsfcrs, zsfcrs_msk
       !!----------------------------------------------------------------  
       ! Initialize
 
-      p_surf_crs    (:,:,:) = 0._wp
-      p_surf_crs_msk(:,:,:) = 0._wp
+
       !
       SELECT CASE ( cd_type )
       
          CASE ('W')
-
-            DO jk = 2 , jpk
-                DO ji = 2, nlei_crs
-                   ijie = mie_crs(ji)
-                   ijis = mis_crs(ji)
-                   DO jj = njstart, njend
-                      ijje = mje_crs(jj) 
-                      ijjs = mjs_crs(jj)
-                      DO jii = ijis, ijie
-                         DO jjj = ijjs, ijje
-                            zsfc = p_e1(jii,jjj) * p_e2(jii,jjj)
-                            p_surf_crs(ji,jj,jk)     = p_surf_crs(ji,jj,jk) + zsfc
-                            p_surf_crs_msk(ji,jj,jk) = p_surf_crs_msk(ji,jj,jk) + zsfc * p_mask(jii,jjj,jk-1)               
-                         ENDDO
-                      ENDDO
-                      IF( njstart == 1 ) p_surf_crs(ji,jj,jk) = p_surf_crs(ji,jj,jk) * nn_facty
-                   ENDDO                  
+      
+           CALL wrk_alloc( jpi, jpj, zsurf )
+           zsurf(:,:) =  p_e1(:,:) * p_e2(:,:)
+        
+            DO ji = nistr, niend, nn_factx
+               DO jj   = njstr, njend, nn_facty
+                  ii   = ( ji - mis_crs(2) ) * rfactx_r + 2                 ! cordinate in parent grid
+                  ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                  ijje = mje_crs(ij)
+                  ijie = mie_crs(ii)
+                  !
+                  zsfcrs     =  zsurf(ji,jj  ) + zsurf(ji+1,jj  ) + zsurf(ji+2,jj  )  &
+                    &         + zsurf(ji,jj+1) + zsurf(ji+1,jj+1) + zsurf(ji+2,jj+1)  &
+                    &         + zsurf(ji,jj+2) + zsurf(ji+1,jj+2) + zsurf(ji+2,jj+2)  
+                  !
+                  zsfcrs_msk =  zsurf(ji  ,jj  ) * p_mask(ji  ,jj  ,1)  &
+                    &         + zsurf(ji+1,jj  ) * p_mask(ji+1,jj  ,1)  &
+                    &         + zsurf(ji+2,jj  ) * p_mask(ji+2,jj  ,1)  &
+                    &         + zsurf(ji  ,jj+1) * p_mask(ji  ,jj+1,1)  &
+                    &         + zsurf(ji+1,jj+1) * p_mask(ji+1,jj+1,1)  &
+                    &         + zsurf(ji+2,jj+1) * p_mask(ji+2,jj+1,1)  &
+                    &         + zsurf(ji  ,jj+2) * p_mask(ji  ,jj+2,1)  &
+                    &         + zsurf(ji+1,jj+2) * p_mask(ji+1,jj+2,1)  &
+                    &         + zsurf(ji+2,jj+2) * p_mask(ji+2,jj+2,1)
+                  !
+                  p_surf_crs    (ii,ij,1) = zsfcrs
+                  p_surf_crs_msk(ii,ij,1) = zsfcrs_msk
+                  !
                ENDDO
             ENDDO
+            DO jk = 2, jpk
+               !
+               p_surf_crs(:,:,jk) = p_surf_crs(:,:,1)
+               !
+               DO ji = nistr, niend, nn_factx
+                  DO jj   = njstr, njend, nn_facty
+                     ii   = ( ji - mis_crs(2) ) * rfactx_r + 2                 ! cordinate in parent grid
+                     ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                     ijje = mje_crs(ij)
+                     ijie = mie_crs(ii)
+                     !
+                     zsfcrs_msk =  zsurf(ji  ,jj  ) * p_mask(ji  ,jj  ,jk-1)  &
+                       &         + zsurf(ji+1,jj  ) * p_mask(ji+1,jj  ,jk-1)  &
+                       &         + zsurf(ji+2,jj  ) * p_mask(ji+2,jj  ,jk-1)  &
+                       &         + zsurf(ji  ,jj+1) * p_mask(ji  ,jj+1,jk-1)  &
+                       &         + zsurf(ji+1,jj+1) * p_mask(ji+1,jj+1,jk-1)  &
+                       &         + zsurf(ji+2,jj+1) * p_mask(ji+2,jj+1,jk-1)  &
+                       &         + zsurf(ji  ,jj+2) * p_mask(ji  ,jj+2,jk-1)  &
+                       &         + zsurf(ji+1,jj+2) * p_mask(ji+1,jj+2,jk-1)  &
+                       &         + zsurf(ji+2,jj+2) * p_mask(ji+2,jj+2,jk-1)
+                       !
+                       p_surf_crs_msk(ii,ij,jk) = zsfcrs_msk
+                       !
+                   ENDDO
+                ENDDO
+            ENDDO
                         
-            ! surface   ; ji = 1 et jpi_crs definit par cyclique est-ouest et pivot T
-            DO ji = 2, nlei_crs
-               ijie = mie_crs(ji)
-               ijis = mis_crs(ji)
-               DO jj = njstart, njend                           
-                  ijje = mje_crs(jj) 
-                  ijjs = mjs_crs(jj)
-                  DO jii = ijis, ijie 
-                     DO jjj = ijjs, ijje
-                        zsfc = p_e1(jii,jjj) * p_e2(jii,jjj)
-                        p_surf_crs(ji,jj,1)     = p_surf_crs(ji,jj,1) + zsfc
-                        p_surf_crs_msk(ji,jj,1) = p_surf_crs_msk(ji,jj,1) * zsfc * p_mask(jii,jjj,1)               
-                     ENDDO  
-                  ENDDO
-                  IF( njstart == 1 ) p_surf_crs(ji,jj,1) = p_surf_crs(ji,jj,1) * nn_facty               
-               ENDDO
-            ENDDO
+            CALL wrk_dealloc( jpi, jpj, zsurf )
+          
+         CASE( 'V' )
+           
+           CALL wrk_alloc( jpi, jpj, jpk, zsurf3d )
+           DO jk = 1, jpk
+              zsurf3d(:,:,jk) =  p_e1(:,:) * p_e3(:,:,jk)
+           ENDDO
+         
+           DO jk = 1, jpk
+              DO ji = nistr, niend, nn_factx
+                 DO jj = njstr, njend, nn_facty
+                    ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                    ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                    ijje = mje_crs(ij)
+                    ijie = mie_crs(ii)
+                    !                  
+                    zsfcrs      =  zsurf3d(ji,ijje,jk) + zsurf3d(ji+1,ijje,jk) + zsurf3d(ji+2,ijje,jk)
+                    !
+                    zsfcrs_msk  =  zsurf3d(ji  ,ijje,jk) * p_mask(ji  ,ijje,jk) &
+                      &          + zsurf3d(ji+1,ijje,jk) * p_mask(ji+1,ijje,jk) &
+                      &          + zsurf3d(ji+2,ijje,jk) * p_mask(ji+2,ijje,jk) 
+                    !
+                    p_surf_crs    (ii,ij,jk) = zsfcrs
+                    p_surf_crs_msk(ii,ij,jk) = zsfcrs_msk
+                    !
+                  ENDDO      
+              ENDDO
+           ENDDO   
+
+           CALL wrk_dealloc( jpi, jpj, jpk, zsurf3d )
             
-       CASE ('U')
-       
-          DO jk = 1 , jpk
-                                             ! ji = 1 et jpi_crs definit par cyclique est-ouest et pivot T
-             DO ji = 2, nlei_crs
-                ijie = mie_crs(ji)
-                ijis = mis_crs(ji)
-                DO jj = njstart, njend 
-                   ijje = mje_crs(jj)  
-                   ijjs = mjs_crs(jj)
-                   DO jii = ijis, ijie
-                      DO jjj = ijjs, ijje   
-                         zsfc =  p_e3(jii,jjj,jk) * p_e2(jii,jjj)  
-                         p_surf_crs(ji,jj,jk)     = p_surf_crs(ji,jj,jk) + zsfc
-                         p_surf_crs_msk(ji,jj,jk) = p_surf_crs_msk(ji,jj,jk) + zsfc * p_mask(jii,jjj,jk)               
-                      ENDDO
-                   ENDDO
-                  IF( njstart == 1 ) p_surf_crs(ji,jj,jk) = p_surf_crs(ji,jj,jk) * nn_facty
-               ENDDO
-            ENDDO
-         ENDDO
-       
-      CASE ('V')
-          
-         DO jk = 1 , ijpk
-                                                                   ! ji = 1 et jpi_crs definit par cyclique est-ouest et pivot T
-            DO ji = 2, nlei_crs
-               ijie = mie_crs(ji)
-               ijis = mis_crs(ji)
-               DO jj = njstart, njend                  ! jj = jpj_crs definit par pivot T 
-                  ijje = mje_crs(jj)   
-                  ijjs = mjs_crs(jj)  
-                  DO jii = ijis, ijie
-                     DO jjj = ijjs, ijje
-                        zsfc = p_e3(jii,jjj,jk) * p_e1(jii,jjj) 
-                        p_surf_crs(ji,jj,jk)     = p_surf_crs(ji,jj,jk) + zsfc
-                        p_surf_crs_msk(ji,jj,jk) = p_surf_crs_msk(ji,jj,jk) + zsfc * p_mask(jii,jjj,jk)   
-                     ENDDO
-                  ENDDO
-               ENDDO
-            ENDDO
-         ENDDO
+         CASE( 'U' )
+         
+           CALL wrk_alloc( jpi, jpj, jpk, zsurf3d )
+           DO jk = 1, jpk
+              zsurf3d(:,:,jk) =  p_e2(:,:) * p_e3(:,:,jk)
+           ENDDO
+         
+           DO jk = 1, jpk
+              DO ji = nistr, niend, nn_factx
+                 DO jj = njstr, njend, nn_facty
+                    ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
+                    ij   = ( jj - mjs_crs(2) ) * rfacty_r + 2
+                    ijje = mje_crs(ij)
+                    ijie = mie_crs(ii)
+                    !                  
+                    zsfcrs      =  zsurf3d(ijie,jj,jk) + zsurf3d(ijie,jj+1,jk) + zsurf3d(ijie,jj+2,jk)
+                    !
+                    zsfcrs_msk  =  zsurf3d(ijie  ,jj,jk) * p_mask(ijie,jj  ,jk) &
+                      &          + zsurf3d(ijie,jj+1,jk) * p_mask(ijie,jj+1,jk) &
+                      &          + zsurf3d(ijie,jj+2,jk) * p_mask(ijie,jj+2,jk) 
+                    !
+                    p_surf_crs    (ii,ij,jk) = zsfcrs
+                    p_surf_crs_msk(ii,ij,jk) = zsfcrs_msk
+                    !
+                  ENDDO      
+              ENDDO
+           ENDDO  
+ 
+           CALL wrk_dealloc( jpi, jpj, jpk, zsurf3d )
+
       END SELECT
- 
-      
-      CALL crs_lbc_lnk( p_surf_crs    (:,:,:), cd_type, 1.0, pval=1.0 )
-      CALL crs_lbc_lnk( p_surf_crs_msk(:,:,:), cd_type, 1.0, pval=1.0 )
+      
+      CALL crs_lbc_lnk( p_surf_crs    , cd_type, 1.0, pval=1.0 )
+      CALL crs_lbc_lnk( p_surf_crs_msk, cd_type, 1.0, pval=1.0 )
 
 
@@ -1155,8 +1591,6 @@
          nlej_crs   = jpj_crs
          
-        !!! Calculs suivant une découpage en j
-     
+        ! Calculs suivant une découpage en j
         DO jn = 1, jpnij, jpni
-         
            IF( jn < (jpnij-jpni + 1)) THEN
               nlejt_crs(jn) = AINT( REAL( ( jpjglo - (njmppt(jn     ) - 1) ) / nn_facty, wp ) ) &
@@ -1167,5 +1601,4 @@
            
            SELECT CASE( ibonjt(jn) )
-           
               CASE ( -1 )
                 IF( MOD( jpjglo - njmppt(jn), nn_facty) > 0 )    nlejt_crs(jn) = nlejt_crs(jn) + 1
@@ -1198,11 +1631,8 @@
            ENDDO
         ENDDO 
-                 WRITE(numout,*) ' njmppt_crs', njmppt_crs
-         nlej_crs  = nlejt_crs(nproc + 1) 
-         nlcj_crs  = nlcjt_crs(nproc + 1)
-         nldj_crs  = nldjt_crs(nproc + 1)
-         njmpp_crs = njmppt_crs(nproc + 1)
-
-
+        nlej_crs  = nlejt_crs(nproc + 1) 
+        nlcj_crs  = nlcjt_crs(nproc + 1)
+        nldj_crs  = nldjt_crs(nproc + 1)
+        njmpp_crs = njmppt_crs(nproc + 1)
 
          !!!! Calcul suivant un decoupage en i
@@ -1218,5 +1648,4 @@
            
                  CASE ( -1 )
-                 
                    nleit_crs(jn) = nleit_crs(jn) + jpreci           
                    nlcit_crs(jn) = nleit_crs(jn) + jpreci
@@ -1224,5 +1653,4 @@
               
                  CASE ( 0 )
-              
                    nleit_crs(jn) = nleit_crs(jn) + jpreci
                    nlcit_crs(jn) = nleit_crs(jn) + jpreci
@@ -1230,5 +1658,4 @@
                 
                  CASE ( 1, 2 )
-                 
                    IF( MOD( jpiglo - nimppt(jn), nn_factx) > 0 )    nleit_crs(jn) = nleit_crs(jn) + 1
                    nleit_crs(jn) = nleit_crs(jn) + jpreci
@@ -1254,7 +1681,5 @@
          nimpp_crs = nimppt_crs(nproc + 1)
 
-
-         
-         !!! rajouter la condition stop
+         ! rajouter la condition stop
          IF( jpizoom /= 1 .OR. jpjzoom /= 1)    STOP 
          DO ji = 1, jpi_crs
@@ -1265,6 +1690,4 @@
          ENDDO
 
-         
-          
       ENDIF
       
@@ -1434,59 +1857,42 @@
       END SELECT
 
-
-        ! Pad the boundaries, do not know if it is necessary
-         mis2_crs(1) = 1            ; mis2_crs(jpiglo_crs) = mie2_crs(jpiglo_crs - 1) + 1   
-         mie2_crs(1) = nn_factx     ; mie2_crs(jpiglo_crs) = jpiglo                         
-         mje2_crs(1) = mjs2_crs(2)-1; mje2_crs(jpjglo_crs) = jpjglo
-         mjs2_crs(1) = 1            ; mjs2_crs(jpjglo_crs) = mje2_crs(jpjglo_crs) - nn_facty + 1 
- 
-
-         
+     ! Pad the boundaries, do not know if it is necessary
+      mis2_crs(1) = 1             ; mis2_crs(jpiglo_crs) = mie2_crs(jpiglo_crs - 1) + 1   
+      mie2_crs(1) = nn_factx      ; mie2_crs(jpiglo_crs) = jpiglo                         
+      mje2_crs(1) = mjs2_crs(2)-1 ; mje2_crs(jpjglo_crs) = jpjglo
+      mjs2_crs(1) = 1             ; mjs2_crs(jpjglo_crs) = mje2_crs(jpjglo_crs) - nn_facty + 1 
  
       IF( .NOT. lk_mpp ) THEN     
-         njstart = 1     ;    njend =  jpj_crsm1
-         mis_crs(:) = mis2_crs(:) 
-         mie_crs(:) = mie2_crs(:)
-         mjs_crs(:) = mjs2_crs(:) 
-         mje_crs(:) = mje2_crs(:)
+        mis_crs(:) = mis2_crs(:) 
+        mie_crs(:) = mie2_crs(:)
+        mjs_crs(:) = mjs2_crs(:) 
+        mje_crs(:) = mje2_crs(:) 
       ELSE
-         !
-         IF( nldj==1 )  THEN  ;   njstart = 1
-         ELSE                 ;   njstart = 2
-         ENDIF
-         !
-         IF(nlcj == nlej) THEN
-            njend = nlej_crs - 1
-         ELSE
-            njend = nlej_crs
-         ENDIF
-         
-         DO jj = 1, nlej_crs
-            mjs_crs(jj) = mjs2_crs(mjg_crs(jj)) - njmpp + 1
-            mje_crs(jj) = mje2_crs(mjg_crs(jj)) - njmpp + 1
-         ENDDO
-         DO ji = 1, nlei_crs
-            mis_crs(ji) = mis2_crs(mig_crs(ji)) - nimpp + 1
-            mie_crs(ji) = mie2_crs(mig_crs(ji)) - nimpp + 1
-         ENDDO
-
-         !
+        DO jj = 1, nlej_crs
+           mjs_crs(jj) = mjs2_crs(mjg_crs(jj)) - njmpp + 1
+           mje_crs(jj) = mje2_crs(mjg_crs(jj)) - njmpp + 1
+        ENDDO
+        DO ji = 1, nlei_crs
+           mis_crs(ji) = mis2_crs(mig_crs(ji)) - nimpp + 1
+           mie_crs(ji) = mie2_crs(mig_crs(ji)) - nimpp + 1
+        ENDDO
       ENDIF
-         WRITE(numout,*) 'crs_init. coarse grid bounds on parent grid'
-         WRITE(numout,*) 'mjg_crs='  , mjg_crs
-         WRITE(numout,*) 'mig_crs='  , mig_crs
-         WRITE(numout,*) 'nimpp_crs=', nimpp_crs
-         WRITE(numout,*) 'njmpp_crs=', njmpp_crs
-         WRITE(numout,*) 'njend='    , njend
-         WRITE(numout,*) 'mis_crs='  , mis_crs
-         WRITE(numout,*) 'mie_crs='  , mie_crs
-         WRITE(numout,*) 'mjs_crs='  , mjs_crs
-         WRITE(numout,*) 'mje_crs='  , mje_crs
-   
-   
+      njstr = mjs_crs(2)  ;   njend = mjs_crs(nlcj_crs - 1)
+      nistr = mis_crs(2)  ;   niend = mis_crs(nlci_crs - 1)
+      WRITE(numout,*) 'crs_init. coarse grid bounds on parent grid'
+      WRITE(numout,*) 'mjg_crs='  , mjg_crs
+      WRITE(numout,*) 'mig_crs='  , mig_crs
+      WRITE(numout,*) 'nimpp_crs=', nimpp_crs
+      WRITE(numout,*) 'njmpp_crs=', njmpp_crs
+      WRITE(numout,*) 'njend='    , njend
+      WRITE(numout,*) 'mis_crs='  , mis_crs
+      WRITE(numout,*) 'mie_crs='  , mie_crs
+      WRITE(numout,*) 'mjs_crs='  , mjs_crs
+      WRITE(numout,*) 'mje_crs='  , mje_crs
+      !
    END SUBROUTINE crs_dom_def
    
    SUBROUTINE crs_dom_bat
-    !!----------------------------------------------------------------
+      !!----------------------------------------------------------------
       !!               *** SUBROUTINE crs_dom_bat ***
       !! ** Purpose :  coarsenig bathy
@@ -1494,10 +1900,9 @@
       !! 
       !!  local variables
-
+      INTEGER  :: ji,jj,jk      ! dummy indices
+      REAL(wp), DIMENSION(:,:)  , POINTER :: zmbk
+      !!----------------------------------------------------------------
    
-   INTEGER  :: ji,jj,jk      ! dummy indices
-   REAL(wp), DIMENSION(:,:)  , POINTER :: zmbk
-   
-   CALL wrk_alloc( jpi_crs, jpj_crs, zmbk )
+      CALL wrk_alloc( jpi_crs, jpj_crs, zmbk )
    
       mbathy_crs(:,:) = jpkm1
@@ -1538,8 +1943,5 @@
       END DO
 
-      WRITE(numout,*) 'crsini. Set mbku, mkbv'
-
       ! convert into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
-  
       zmbk(:,:) = 1.e0;    
       zmbk(:,:) = REAL( mbku_crs(:,:), wp )   ;   CALL crs_lbc_lnk(zmbk,'U',1.0) ; mbku_crs  (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
@@ -1547,6 +1949,5 @@
       !
       CALL wrk_dealloc( jpi_crs, jpj_crs, zmbk )
-      WRITE(numout,*) 'crs_init   = finished section 3d.1 jpi=', jpi, 'jpj=',jpj, ' jpk=', jpk
-
+      !
    END SUBROUTINE crs_dom_bat
 
Index: branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsdomwri.F90
===================================================================
--- branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsdomwri.F90	(revision 3864)
+++ branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsdomwri.F90	(revision 3895)
@@ -21,4 +21,5 @@
    USE lib_mpp         ! MPP library
 !   USE wrk_nemo        ! Memory allocation
+   USE iom_def
    USE iom
    USE crs         ! coarse grid domain
@@ -295,12 +296,12 @@
       SELECT CASE ( MOD(nn_msh_crs, 3) )
       CASE ( 1 )                
-         CALL iom_close( inum0 )
+         CALL crs_iom_close( inum0 )
       CASE ( 2 )
-         CALL iom_close( inum1 )
-         CALL iom_close( inum2 )
+         CALL crs_iom_close( inum1 )
+         CALL crs_iom_close( inum2 )
       CASE ( 0 )
-         CALL iom_close( inum2 )
-         CALL iom_close( inum3 )
-         CALL iom_close( inum4 )
+         CALL crs_iom_close( inum2 )
+         CALL crs_iom_close( inum3 )
+         CALL crs_iom_close( inum4 )
       END SELECT
       !
Index: branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsini.F90
===================================================================
--- branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsini.F90	(revision 3864)
+++ branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsini.F90	(revision 3895)
@@ -15,5 +15,4 @@
    USE in_out_manager
    USE par_kind, ONLY: wp
-   USE crs
    USE crsdom
    USE crsdomwri
@@ -208,16 +207,15 @@
 
      !    3.d.3   Vertical depth (meters)
-
-     CALL crsfun( p_e1e2t=e1e2t, cd_type='T', cd_op='MAX', p_cmask=tmask_crs, p_ptmask=tmask, &
-          &       p_pfield3d_1=gdept, p_cfield3d=gdept_crs ) 
-     CALL crsfun( p_e1e2t=e1e2t, cd_type='W', cd_op='MAX', p_cmask=tmask_crs, p_ptmask=tmask, &
-          &       p_pfield3d_1=gdepw, p_cfield3d=gdepw_crs )
+     CALL crs_dom_ope( gdept, 'MAX', 'T', tmask, gdept_crs, p_e3=zfse3t ) 
+     CALL crs_dom_ope( gdepw, 'MAX', 'W', tmask, gdepw_crs, p_e3=zfse3w )
+
 
      !    3.d.4   Surfaces 
-      
-     CALL crs_dom_sfc( e1t, e2t, zfse3w, 'W', tmask, e1e2w, e1e2w_msk )
-     CALL crs_dom_sfc( e1u, e2u, zfse3u, 'U', umask, e2e3u, e2e3u_msk )
-     CALL crs_dom_sfc( e1v, e2v, zfse3v, 'V', vmask, e1e3v, e1e3v_msk )
+     CALL crs_dom_sfc( tmask, 'T', e1e2w, e1e2w_msk, p_e1=e1t, p_e2=e2t )
+     CALL crs_dom_sfc( umask, 'U', e2e3u, e2e3u_msk, p_e2=e2u, p_e3=zfse3u )
+     CALL crs_dom_sfc( vmask, 'V', e1e3v, e1e3v_msk, p_e1=e1v, p_e3=zfse3v )
    
+     facsurfu(:,:,:) = umask_crs(:,:,:) * e2e3u_msk(:,:,:) / e2e3u(:,:,:)
+     facsurfv(:,:,:) = vmask_crs(:,:,:) * e1e3v_msk(:,:,:) / e1e3v(:,:,:)
 
      !---------------------------------------------------------
@@ -227,17 +225,13 @@
 
 !!     ! jes. May not need ocean_volume_crs_t, ocean_volume_crs_w as calculated already in trc_init as cvol
-      CALL crsfun_wgt( cd_type='T', cd_op='VOL', p_pmask=tmask, p_e1=e1t, p_e2=e2t, p_fse3=zfse3t, &
-         &             p_cfield3d_1=ocean_volume_crs_t, p_cfield3d_2=facvol_t )
-      
-      r1_bt_crs(:,:,:) = 0._wp 
-      bt_crs(:,:,:) = ocean_volume_crs_t(:,:,:)* facvol_t(:,:,:)
-      WHERE( bt_crs /= 0._wp ) r1_bt_crs(:,:,:) = 1._wp/bt_crs(:,:,:)
-
-      CALL crsfun_wgt( cd_type='W', cd_op='VOL', p_pmask=tmask, p_e1=e1t, p_e2=e2t, p_fse3=zfse3w, &
-         &             p_cfield3d_1=ocean_volume_crs_w, p_cfield3d_2=facvol_w )
-
-     ! 4.c. T volume weights
-      CALL crsfun_wgt( cd_type='T', cd_op='WGT', p_pmask=tmask, p_e1=e1t, p_e2=e2t, p_fse3=zfse3t, p_cfield3d_1=crs_volt_wgt )  
-
+     CALL crs_dom_facvol( tmask, 'T', e1t, e2t, zfse3t, ocean_volume_crs_t, facvol_t )
+     !
+     bt_crs(:,:,:) = ocean_volume_crs_t(:,:,:) * facvol_t(:,:,:)
+     !
+     r1_bt_crs(:,:,:) = 0._wp 
+     WHERE( bt_crs /= 0._wp ) r1_bt_crs(:,:,:) = 1._wp / bt_crs(:,:,:)
+
+     CALL crs_dom_facvol( tmask, 'W', e1t, e2t, zfse3w, ocean_volume_crs_w, facvol_w )
+     !
      !---------------------------------------------------------
      ! 5.  Write out coarse meshmask  (see OPA_SRC/DOM/domwri.F90 for ideas later)
Index: branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsiom.F90
===================================================================
--- branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsiom.F90	(revision 3864)
+++ branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsiom.F90	(revision 3895)
@@ -31,8 +31,13 @@
    !! jes. 28 Jun 2012. TODO. make sure of variable declarations to be placed here or crs_dom.F90
    !!--------------------------------------------------------------------
+   USE timing
+   USE crs
    USE dom_oce         ! ocean space and time domain
-   USE crs
+   USE iom_def         ! iom variables definitions
+   USE netcdf          ! NetCDF library
    USE in_out_manager  ! I/O manager
+   USE lib_mpp         ! MPP library
    USE iom             ! I/O library
+   USE par_kind, ONLY: wp   
 
 
@@ -40,7 +45,22 @@
    PRIVATE
 
-   PUBLIC crs_iom_open, crs_iom_rstput, crs_iom_put
+   PUBLIC crs_iom_open, crs_iom_close, crs_iom_rstput, crs_iom_put
+
+!   PUBLIC crs_iom_varid, crs_iom_get, crs_iom_gettime
+
+
+   INTEGER, PARAMETER ::   jpdomcrs_data          = 1   !: ( 1  :jpiglo_crs, 1  :jpjglo_crs)
+   INTEGER, PARAMETER ::   jpdomcrs_global        = 2   !: ( 1  :jpiglo_crs, 1  :jpjglo_crs)
+   INTEGER, PARAMETER ::   jpdomcrs_local         = 3   !: One of the 3 following cases
+   INTEGER, PARAMETER ::   jpdomcrs_local_full    = 4   !: ( 1  :jpi_crs   , 1  :jpj_crs   )
+   INTEGER, PARAMETER ::   jpdomcrs_local_noextra = 5   !: ( 1  :nlci_crs  , 1  :nlcj_crs  )
+   INTEGER, PARAMETER ::   jpdomcrs_local_noovlap = 6   !: (nldi_crs:nlei_crs  ,nldj_crs:nlej_crs  )
+   INTEGER, PARAMETER ::   jpdomcrs_unknown       = 7   !: No dimension checking
+   INTEGER, PARAMETER ::   jpdomcrs_autoglo       = 8   !:
+   INTEGER, PARAMETER ::   jpdomcrs_autodta       = 9   !:
+   INTEGER            ::   ipdomcrs_local_noovlap_crs, ipdomcrs_local_full_crs, idomcrs_local_noextra_crs
 
    INTEGER                 ::   idomcrs     ! Type of domain to be written (default = jpdom_local_noovlap)
+   INTEGER, DIMENSION(2,5) ::   idompar_crs ! domain parameters:    
    LOGICAL                 ::   llnoov      ! local definition to read overlap
 
@@ -65,6 +85,6 @@
 
       llnoov = (jpni * jpnj ) == jpnij .AND. .NOT. lk_agrif 
-      IF( llnoov ) THEN   ;   idomcrs = jpdom_local_noovlap   ! default definition
-      ELSE                ;   idomcrs = jpdom_local_full      ! default definition
+      IF( llnoov ) THEN   ;   idomcrs = jpdomcrs_local_noovlap   ! default definition
+      ELSE                ;   idomcrs = jpdomcrs_local_full      ! default definition
       ENDIF
       IF ( PRESENT(kdom) ) idomcrs = kdom
@@ -72,7 +92,30 @@
       CALL iom_open( cdname, kiomid, ldwrt, idomcrs, kiolib )
 
+      WRITE(numout,*) 'crs_iom_open. after iom_open call kiomid=', kiomid
+
       CALL dom_grid_glo   ! Return to parent grid domain
 
    END SUBROUTINE crs_iom_open
+
+
+   SUBROUTINE crs_iom_close( kiomid )
+      !!--------------------------------------------------------------------
+      !!                       ***  MODULE crs_iom_open  ***
+      !!
+      !! ** Purpose : open an input file with NF90 on coarsened grid
+      !!---------------------------------------------------------------------
+      !! Arguments
+      INTEGER                , INTENT(inout)           ::   kiomid      ! nf90 identifier of the opened file
+      !! Local variable
+      CHARACTER(LEN=100)  ::   clinfo   ! info character
+      !---------------------------------------------------------------------
+      !
+      WRITE(numout,*) 'crs_iom_close. kiomid=', kiomid
+
+      CALL iom_close( kiomid )
+      WRITE(numout,*) 'crs_iom_close. after iom_open call kiomid=', kiomid
+
+      !    
+   END SUBROUTINE crs_iom_close
 
 
@@ -137,5 +180,4 @@
       ELSEIF( PRESENT(pv_r3d) ) THEN   ;  CALL iom_put( cdvar, pv_r3d )
       ENDIF
-
       CALL dom_grid_glo   ! Return to parent grid domain
 
