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 3777)
+++ /branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/CONFIG/ORCA2_LIM_CRS/EXP00/iodef.xml	(revision 3778)
@@ -291,4 +291,5 @@
        	<field id="soce_crs"   description="salinity"                unit="psu"  axis_ref="deptht" />
 	<field id="ssh_crs"    description="sea surface height"      unit="m"                      />
+        <field id="ssh2_crs"    description="sea surface height"      unit="m"                      />
 	<field id="sst_crs"    description="sea surface temperature" unit="degC"                   />
 	<field id="sss_crs"    description="sea surface salinity"    unit="psu"                    />
@@ -311,4 +312,5 @@
       <group id="gcrs_W"  axis_ref="depthw" grid_ref="grid_W_crs">
 	<field id="woce_crs"   description="ocean vertical velocity"   unit="m/s" />
+        <field id="woce2_crs"   description="ocean vertical velocity"   unit="m/s" />
 	<field id="wocet_crs"  description="ocean vertical velocity times temperature" unit="degC.m/s" />
 	<field id="woces_crs"  description="ocean vertical velocity times salinity"    unit="psu.m/s" />
@@ -347,4 +349,5 @@
              <field ref="soce_crs"    name="vosaline" />
              <field ref="ssh_crs"     name="sossheig" />
+             <field ref="ssh2_crs"    name="sossheig2" />
              <field ref="hdiv_crs"    name="vohdiver" />
       	     <field ref="sst_crs"     name="sosstsst" />
@@ -431,4 +434,5 @@
            <field ref="soce_crs"    name="vosaline" />
            <field ref="ssh_crs"     name="sossheig" />
+           <field ref="ssh2_crs"     name="sossheig2" />
 	   <field ref="sst_crs"     name="sosstsst" />
 	   <field ref="sss_crs"     name="sosaline" />
@@ -474,4 +478,5 @@
         <file id="5d_gcrs_W" name="auto" description="ocean U grid coarsened variables" >
            <field ref="woce_crs"    name="vovecrtz"  />
+           <field ref="woce2_crs"   name="vovecrtz2"  />
            <field ref="wocet_crs"   name="vovewoct" />
            <field ref="woces_crs"   name="vovewocs" />
Index: /branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/CONFIG/ORCA2_LIM_CRS/MY_SRC/iom.F90
===================================================================
--- /branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/CONFIG/ORCA2_LIM_CRS/MY_SRC/iom.F90	(revision 3777)
+++ /branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/CONFIG/ORCA2_LIM_CRS/MY_SRC/iom.F90	(revision 3778)
@@ -114,10 +114,12 @@
       CALL set_grid( "grid_V", glamv, gphiv )
       CALL set_grid( "grid_W", glamt, gphit )
-
+      !cc
+      IF(ln_crs) THEN
       ! horizontal coarse grid definition
-      CALL setgrid_crs( "grid_T_crs", glamt_crs, gphit_crs )
-      CALL setgrid_crs( "grid_U_crs", glamu_crs, gphiu_crs )
-      CALL setgrid_crs( "grid_V_crs", glamv_crs, gphiv_crs )
-      CALL setgrid_crs( "grid_W_crs", glamt_crs, gphit_crs )
+         CALL setgrid_crs( "grid_T_crs", glamt_crs, gphit_crs )
+         CALL setgrid_crs( "grid_U_crs", glamu_crs, gphiu_crs )
+         CALL setgrid_crs( "grid_V_crs", glamv_crs, gphiv_crs )
+         CALL setgrid_crs( "grid_W_crs", glamt_crs, gphit_crs )
+      ENDIF
 
       ! vertical grid definition
Index: /branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/CONFIG/ORCA2_LIM_CRS/MY_SRC/lbclnk.F90
===================================================================
--- /branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/CONFIG/ORCA2_LIM_CRS/MY_SRC/lbclnk.F90	(revision 3777)
+++ /branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/CONFIG/ORCA2_LIM_CRS/MY_SRC/lbclnk.F90	(revision 3778)
@@ -161,6 +161,5 @@
             END SELECT
             !                                          ! North fold
-            pt3d( 1 ,jpj,:) = zland
-            pt3d(jpi,jpj,:) = zland
+          
             CALL lbc_nfd( pt3d(:,:,:), cd_type, psgn )
             !
@@ -249,8 +248,5 @@
                pt2d(:, 1 ) = zland
             END SELECT
-            !                                          ! North fold
-            pt2d( 1 ,1  ) = zland 
-            pt2d( 1 ,jpj) = zland 
-            pt2d(jpi,jpj) = zland
+            
             CALL lbc_nfd( pt2d(:,:), cd_type, psgn )
             !
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 3777)
+++ /branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crs.F90	(revision 3778)
@@ -43,5 +43,5 @@
    PRIVATE 
 
-   PUBLIC crsfun
+   PUBLIC crsfun, crs_e3_max, crs_surf
 
    INTERFACE crsfun
@@ -282,11 +282,11 @@
       p_cglam(jpi_crs,1) = p_cglam(2,1)
       ! Fill upper-right corner i=1, j=jpj_crs
-      IF ( nperio == 4 ) THEN
-         p_cgphi(1,jpj_crs) = p_cgphi(jpi_crsm1,jpj_crs-2)
-         p_cglam(1,jpj_crs) = p_cglam(jpi_crsm1,jpj_crs-2)
-      ELSEIF ( nperio == 6 ) THEN
-         p_cgphi(1,jpj_crs) = p_cgphi(jpi_crs,jpj_crsm1)
-         p_cglam(1,jpj_crs) = p_cglam(jpi_crs,jpj_crsm1)
-      ENDIF
+     !cc IF ( nperio == 4 ) THEN
+     !cc    p_cgphi(1,jpj_crs) = p_cgphi(jpi_crsm1,jpj_crs-2)
+     !cc    p_cglam(1,jpj_crs) = p_cglam(jpi_crsm1,jpj_crs-2)
+     !cc ELSEIF ( nperio == 6 ) THEN
+     !cc    p_cgphi(1,jpj_crs) = p_cgphi(jpi_crs,jpj_crsm1)
+     !cc    p_cglam(1,jpj_crs) = p_cglam(jpi_crs,jpj_crsm1)
+     !cc ENDIF
 
       WRITE(numout,*) 'crsfun_coordinates. done'
@@ -398,9 +398,11 @@
 
              zcfield2d_1(:,:) = 0.0 ; zcfield2d_2(:,:) = 0.0        
-             DO ji = 2, jpi_crsm1
+             ! DO ji = 2, jpi_crsm1
+             DO ji = 1, jpi_crs
                 ijie = mie_crs(ji)
                 ijis = mis_crs(ji)
 
-                DO jj = 1, jpj_crsm1
+             !   DO jj = 1, jpj_crsm1
+                DO jj = 1, jpj_crs
                    ijje = mje_crs(jj) 
                    ijjs = mjs_crs(jj)                   
@@ -430,15 +432,72 @@
 
                       ELSE
-
-                         ! Calculate e1 scale factor or if present ze3, unmasked surface area
-                         DO jii = ijis, ijie
-                            zcfield2d_1(ji,jj) = zcfield2d_1(ji,jj) +  ( ze1(jii,ijje) * ze3(jii,ijje,jk) ) 
-                         ENDDO
-
-                         ! Calculate e2 scale factor
-                         DO jjj = ijjs, ijje
-                            zcfield2d_2(ji,jj) = zcfield2d_2(ji,jj) +  ( ze2(ijie,jjj) * ze3(ijie,jjj,jk) )
-                         ENDDO
- 
+                         SELECT CASE ( cd_type )
+                         
+                         CASE ( 'T' )
+                            IF ( nn_factx == 3 ) THEN
+                               ! Calculate e1 scale factor or if present ze3, unmasked surface area
+                               IF (jj == 1) THEN
+                                  DO jii = ijis, ijie
+                                  zcfield2d_1(ji,jj) = zcfield2d_1(ji,jj) +  ( ze1(jii,ijje) * ze3(jii,ijje-1,jk) ) 
+                                  ENDDO
+                               ELSE
+                                  DO jii = ijis, ijie
+                                     zcfield2d_1(ji,jj) = zcfield2d_1(ji,jj) +  ( ze1(jii,ijje-1) * ze3(jii,ijje-1,jk) ) 
+                                  ENDDO
+                               ENDIF
+                            
+                               ! Calculate e2 scale factor
+                               DO jjj = ijjs, ijje
+                                zcfield2d_2(ji,jj) = zcfield2d_2(ji,jj) +  ( ze2(ijie-1,jjj) * ze3(ijie-1,jjj,jk) )
+                               ENDDO
+                            ENDIF
+                            
+                         CASE ( 'U' )
+                            IF ( nn_factx == 3 ) THEN
+                               ! Calculate e1 scale factor or if present ze3, unmasked surface area
+                               IF (jj == 1) THEN
+                                  DO jii = ijis, ijie
+                                  zcfield2d_1(ji,jj) = zcfield2d_1(ji,jj) +  ( ze1(jii+1,ijje) * ze3(jii+1,ijje-1,jk) ) 
+                                  ENDDO
+                               ELSE
+                                  DO jii = ijis, ijie
+                                     zcfield2d_1(ji,jj) = zcfield2d_1(ji,jj) +  ( ze1(jii+1,ijje-1) * ze3(jii+1,ijje-1,jk) ) 
+                                  ENDDO
+                               ENDIF                            
+                               ! Calculate e2 scale factor
+                               DO jjj = ijjs, ijje
+                                zcfield2d_2(ji,jj) = zcfield2d_2(ji,jj) +  ( ze2(ijie,jjj) * ze3(ijie,jjj,jk) )
+                               ENDDO
+                            ENDIF
+                            
+                         CASE ( 'V' )
+                            IF ( nn_factx == 3 ) THEN
+                               ! Calculate e1 scale factor or if present ze3, unmasked surface area
+                               DO jii = ijis, ijie
+                                  zcfield2d_1(ji,jj) = zcfield2d_1(ji,jj) +  ( ze1(jii,ijje) * ze3(jii,ijje,jk) ) 
+                               ENDDO
+                            
+                               ! Calculate e2 scale factor
+                               DO jjj = ijjs, ijje
+                                zcfield2d_2(ji,jj) = zcfield2d_2(ji,jj) +  ( ze2(ijie-1,jjj+1) * ze3(ijie-1,jjj+1,jk) )
+                               ENDDO
+                            ENDIF
+                            
+                            CASE ( 'F' )
+                            IF ( nn_factx == 3 ) THEN
+                               ! Calculate e1 scale factor or if present ze3, unmasked surface area
+                               DO jii = ijis, ijie
+                                  zcfield2d_1(ji,jj) = zcfield2d_1(ji,jj) +  ( ze1(jii+1,ijje) * ze3(jii+1,ijje,jk) ) 
+                               ENDDO
+                            
+                               ! Calculate e2 scale factor
+                               DO jjj = ijjs, ijje
+                                zcfield2d_2(ji,jj) = zcfield2d_2(ji,jj) +  ( ze2(ijie,jjj+1) * ze3(ijie,jjj+1,jk) )
+                               ENDDO
+                            ENDIF
+                         END SELECT
+                         
+                         
+                         
                          IF ( PRESENT(p_cfield3d_1) ) THEN
 
@@ -468,4 +527,80 @@
 
                    ENDIF
+                   
+                   IF ( cd_op == 'POS' ) THEN      !cc
+                   
+                      IF ( nn_factx == 3 .AND. nn_facty == 3) THEN
+                      
+                         SELECT CASE ( cd_type )
+                         
+                            CASE ( 'T' )
+                            
+                               IF ((jj == 1) .AND. (ji == 1)) THEN
+                                  ! Calculate e1 scale factor or if present ze3, unmasked surface area
+                                  zcfield2d_1(ji,jj) = ( ze1(ijie,ijje  ) * ze3(ijie,ijje,jk) ) * nn_factx
+                            
+                                  ! Calculate e2 scale factor
+                                  zcfield2d_2(ji,jj) = ( ze2(ijie,ijje  ) * ze3(ijie,ijje,jk) ) * nn_facty
+                               ELSEIF (jj == 1) THEN
+                                  ! Calculate e1 scale factor or if present ze3, unmasked surface area
+                                  zcfield2d_1(ji,jj) = ( ze1(ijie-1,ijje  ) * ze3(ijie-1,ijje,jk) ) * nn_factx
+                            
+                                  ! Calculate e2 scale factor
+                                  zcfield2d_2(ji,jj) = ( ze2(ijie-1,ijje  ) * ze3(ijie-1,ijje,jk) ) * nn_facty
+                               ELSEIF (ji == 1) THEN 
+                                  ! Calculate e1 scale factor or if present ze3, unmasked surface area
+                                  zcfield2d_1(ji,jj) = ( ze1(ijie,ijje-1) * ze3(ijie,ijje-1,jk) ) * nn_factx
+                            
+                                  ! Calculate e2 scale factor
+                                  zcfield2d_2(ji,jj) = ( ze2(ijie,ijje-1) * ze3(ijie,ijje-1,jk) ) * nn_facty
+                               ELSE
+                                 ! Calculate e1 scale factor or if present ze3, unmasked surface area
+                                  zcfield2d_1(ji,jj) = ( ze1(ijie-1,ijje-1) * ze3(ijie-1,ijje-1,jk) ) * nn_factx
+                            
+                                  ! Calculate e2 scale factor
+                                  zcfield2d_2(ji,jj) = ( ze2(ijie-1,ijje-1) * ze3(ijie-1,ijje-1,jk) ) * nn_facty
+                               ENDIF
+                               
+                            CASE ( 'U' )
+                               IF (jj == 1) THEN
+                                  ! Calculate e1 scale factor or if present ze3, unmasked surface area
+                                  zcfield2d_1(ji,jj) = ( ze1(ijie  ,ijje  ) * ze3(ijie  ,ijje,jk) ) * nn_factx
+                            
+                                  ! Calculate e2 scale factor
+                                  zcfield2d_2(ji,jj) = ( ze2(ijie  ,ijje  ) * ze3(ijie  ,ijje,jk) ) * nn_facty
+                               ELSE
+                                  ! Calculate e1 scale factor or if present ze3, unmasked surface area
+                                  zcfield2d_1(ji,jj) = ( ze1(ijie  ,ijje-1) * ze3(ijie  ,ijje-1,jk) ) * nn_factx
+                            
+                                  ! Calculate e2 scale factor
+                                  zcfield2d_2(ji,jj) = ( ze2(ijie  ,ijje-1) * ze3(ijie  ,ijje-1,jk) ) * nn_facty
+                              ENDIF
+                                
+                            CASE ( 'V' )
+                               IF (ji == 1) THEN 
+                               ! Calculate e1 scale factor or if present ze3, unmasked surface area
+                               zcfield2d_1(ji,jj) = ( ze1(ijie,ijje  ) * ze3(ijie,ijje  ,jk) ) * nn_factx
+                            
+                               ! Calculate e2 scale factor
+                               zcfield2d_2(ji,jj) = ( ze2(ijie,ijje  ) * ze3(ijie,ijje  ,jk) ) * nn_facty
+                               ELSE
+                               ! Calculate e1 scale factor or if present ze3, unmasked surface area
+                               zcfield2d_1(ji,jj) = ( ze1(ijie-1,ijje  ) * ze3(ijie-1,ijje  ,jk) ) * nn_factx
+                            
+                               ! Calculate e2 scale factor
+                               zcfield2d_2(ji,jj) = ( ze2(ijie-1,ijje  ) * ze3(ijie-1,ijje  ,jk) ) * nn_facty
+                               ENDIF
+                               
+                            CASE ( 'F' )
+                               ! Calculate e1 scale factor or if present ze3, unmasked surface area
+                               zcfield2d_1(ji,jj) = ( ze1(ijie  ,ijje  ) * ze3(ijie  ,ijje  ,jk) ) * nn_factx
+                            
+                               ! Calculate e2 scale factor
+                               zcfield2d_2(ji,jj) = ( ze2(ijie  ,ijje  ) * ze3(ijie  ,ijje  ,jk) ) * nn_facty
+                               
+                         END SELECT
+                      ENDIF
+                   ENDIF                    !cc
+                   
 
                    IF ( cd_op == 'WGT' ) THEN
@@ -540,5 +675,5 @@
 
          ! Take care of the 2D arrays
-         IF ( cd_op == 'SUM' ) THEN
+         IF ( cd_op == 'SUM' .OR.  cd_op == 'POS') THEN
             IF ( PRESENT(p_cfield2d_1) ) THEN
                p_cfield2d_1(:,:) = zcfield2d_1(:,:)
@@ -546,15 +681,17 @@
 
                ! Fill up jrow=1 which is zeroed out or not handled by lbc_lnk and lbc_nfd
-               p_cfield2d_1(:,1) = zcfield2d_1(:,1)
+               p_cfield2d_1(:,1) = zcfield2d_1(:,1)  !cc 
                ! Fill i=1, i=jpi at j=1
                p_cfield2d_1(1,1) = p_cfield2d_1(jpi_crsm1,1)
                p_cfield2d_1(jpi_crs,1) = p_cfield2d_1(2,1)
-
-               ! Fill upper-right corner i=1, j=jpj_crs
-               IF ( nperio == 4 ) THEN
-                  p_cfield2d_1(1,jpj_crs) = p_cfield2d_1(jpi_crsm1,jpj_crs-2)
-               ELSEIF ( nperio == 6 ) THEN
-                  p_cfield2d_1(1,jpj_crs) = p_cfield2d_1(jpi_crs,jpj_crsm1)
-               ENDIF
+               
+             !cc  p_cfield2d_1(1,jpj_crs-1) = p_cfield2d_1(3,jpj_crs-1)
+
+              ! Fill upper-right corner i=1, j=jpj_crs
+              !cc IF ( nperio == 4 ) THEN on écrase les valeurs limites calculées dans crs_lbc_lnk
+              !cc   p_cfield2d_1(1,jpj_crs) = p_cfield2d_1(jpi_crsm1,jpj_crs-2)
+              !cc ELSEIF ( nperio == 6 ) THEN
+              !cc    p_cfield2d_1(1,jpj_crs) = p_cfield2d_1(jpi_crs,jpj_crsm1)
+              !cc ENDIF
 
             ENDIF
@@ -562,5 +699,5 @@
                p_cfield2d_2(:,:) = zcfield2d_2(:,:)
                CALL crs_lbc_lnk( cd_type,1.0,pt2d=p_cfield2d_2 )
-
+               
                ! Fill up jrow=1 which is zeroed out or not handled by lbc_lnk and lbc_nfd
                p_cfield2d_2(:,1) = zcfield2d_2(:,1)
@@ -569,11 +706,15 @@
                p_cfield2d_2(1,1) = p_cfield2d_2(jpi_crsm1,1)
                p_cfield2d_2(jpi_crs,1) = p_cfield2d_2(2,1)
-
+               IF ( cd_op == 'SUM') THEN 
+                  DO jii = 1 , jpiglo_crs
+                  p_cfield2d_2(jii,1) = p_cfield2d_2(jii,1) * 3
+                  ENDDO
+               ENDIF
                ! Fill upper-right corner i=1, j=jpj_crs
-               IF ( nperio == 4 ) THEN
-                  p_cfield2d_2(1,jpj_crs) = p_cfield2d_2(jpi_crsm1,jpj_crs-2)
-               ELSEIF ( nperio == 6 ) THEN
-                  p_cfield2d_2(1,jpj_crs) = p_cfield2d_2(jpi_crs,jpj_crsm1)
-               ENDIF
+           !cc     IF ( nperio == 4 ) THEN
+           !cc       p_cfield2d_2(1,jpj_crs) = p_cfield2d_2(jpi_crsm1,jpj_crs-2)
+           !cc    ELSEIF ( nperio == 6 ) THEN
+           !cc       p_cfield2d_2(1,jpj_crs) = p_cfield2d_2(jpi_crs,jpj_crsm1)
+           !cc ENDIF
             ENDIF
 
@@ -588,5 +729,5 @@
 
          ! Take care now of 3d arrays
-         IF ( cd_op == 'SUM' .OR. cd_op == 'VOL' ) THEN
+         IF ( cd_op == 'SUM' .OR. cd_op == 'VOL' .OR. cd_op == 'POS'  ) THEN
             CALL crs_lbc_lnk( cd_type,1.0,pt3d1=zcfield3d_1 ) 
             IF ( PRESENT(p_cfield3d_1) ) p_cfield3d_1(:,:,:) = zcfield3d_1(:,:,:)
@@ -607,5 +748,5 @@
 
       ELSE
-         IF ( cd_op == 'SUM' ) THEN      
+         IF ( cd_op == 'SUM' .OR. cd_op == 'POS'  ) THEN      
             IF ( PRESENT(p_cfield2d_1) ) THEN
                p_cfield2d_1(:,:) = zcfield2d_1(:,:)
@@ -918,5 +1059,6 @@
                ijis = mis_crs(ji)
 
-               DO jj = 2, jpj_crsm1
+            ! DO jj = 2, jpj_crsm1
+              DO jj = 1, jpj_crsm1
                   ijje = mje_crs(jj) 
                   ijjs = mjs_crs(jj)                   
@@ -1027,4 +1169,377 @@
    END SUBROUTINE crsfun_TW
 
+SUBROUTINE crs_e3_max( p_e3, cd_type, p_mask, p_e3_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
+      !!
+      !! 
+      !! 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...
+      !!----------------------------------------------------------------
+      !! 
+      !!  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,jpk), OPTIONAL, INTENT(in) :: p_e3         ! 3D tracer T or W on parent grid
+      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), OPTIONAL, INTENT(out):: p_e3_crs ! 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), ALLOCATABLE, DIMENSION(:,:,:) :: ze3, ze3_crs, zpmask  
+      !!----------------------------------------------------------------  
+      ! Initialize
+
+      ALLOCATE( ze3(jpi,jpj,jpk), zpmask(jpi,jpj,jpk) )
+      ALLOCATE( ze3_crs(jpi_crs,jpj_crs,jpk) )
+
+      ! Arrays, scalars initialization 
+      ze3(:,:,:)       = p_e3(:,:,:)
+      ze3_crs(:,:,:)   = 0.0
+      zpmask(:,:,:)    = p_mask(:,:,:)
+      ijpk             = jpk
+
+      SELECT CASE ( cd_type )
+      
+         CASE ('T')
+         
+            DO jk = 1 , ijpk
+            
+               DO ji = 1, jpi_crs                       ! ji = 1 et jpi_crs definit par cyclique est-ouest et pivot T
+                   ijie = mie_crs(ji)
+                   ijis = mis_crs(ji)
+
+                   DO jj = 1, jpj_crs                   ! jj = jpj_crs definit par pivot T 
+                      ijje = mje_crs(jj) 
+                      ijjs = mjs_crs(jj)  
+                   
+                      DO jii = ijis, ijie
+                         DO jjj = ijjs, ijje
+                            ze3_crs(ji,jj,jk) = max( ze3_crs(ji,jj,jk), ze3(jii,jjj,jk) * zpmask(jii,jjj,jk)  )
+                         ENDDO
+                      ENDDO
+                   ENDDO
+               ENDDO
+            ENDDO
+            
+         CASE ('W')
+         
+            DO jk = 2 , ijpk
+            
+               DO ji = 1, jpi_crs                       ! ji = 1 et jpi_crs definit par cyclique est-ouest et pivot T
+                   ijie = mie_crs(ji)
+                   ijis = mis_crs(ji)
+
+                   DO jj = 1, jpj_crs                   ! jj = jpj_crs definit par pivot T 
+                      ijje = mje_crs(jj) 
+                      ijjs = mjs_crs(jj)  
+                   
+                      DO jii = ijis, ijie
+                         DO jjj = ijjs, ijje
+                            ze3_crs(ji,jj,jk) = max( ze3_crs(ji,jj,jk), ze3(jii,jjj,jk) * zpmask(jii,jjj,jk-1)  )
+                         ENDDO
+                      ENDDO
+                   ENDDO
+               ENDDO
+            ENDDO
+             
+            jk = 1                                           ! cas particulier car zpmask(jii,jjj,0) n'existe pas      
+          
+            DO ji = 1, jpi_crs                       
+               ijie = mie_crs(ji)
+               ijis = mis_crs(ji)
+
+               DO jj = 1, jpj_crs                 
+                  ijje = mje_crs(jj) 
+                  ijjs = mjs_crs(jj)  
+                   
+                  DO jii = ijis, ijie
+                     DO jjj = ijjs, ijje
+                        ze3_crs(ji,jj,jk) = max( ze3_crs(ji,jj,jk), ze3(jii,jjj,jk) * zpmask(jii,jjj,jk)  )
+                     ENDDO
+                  ENDDO
+               ENDDO
+            ENDDO  
+            
+         END SELECT 
+         
+         p_e3_crs(:,:,:) = ze3_crs(:,:,:)
+         
+         CALL crs_lbc_lnk( cd_type, 1.0, pt3d1=p_e3_crs )
+
+         ! lbcnlk met la ligne jpj = 1 a 0 donc il faut la remettre en ne pas oubliant le cyclique est-ouest
+                   
+         p_e3_crs(   :   ,1,:) = ze3_crs(    :    ,1,:)   
+         p_e3_crs(   1   ,1,:) = ze3_crs(jpi_crsm1,1,:)
+         p_e3_crs(jpi_crs,1,:) = ze3_crs(    2    ,1,:)       
+ 
+
+      DEALLOCATE( ze3 , zpmask )
+      DEALLOCATE( ze3_crs )
+ 
+
+   END SUBROUTINE crs_e3_max
+
+
+SUBROUTINE crs_surf(p_e1, p_e2 ,p_e3, cd_type, p_mask, surf_crs, surf_msk_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
+      !!
+      !! 
+      !! 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...
+      !!----------------------------------------------------------------
+      !! 
+      !!  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,jpk), OPTIONAL, INTENT(in) :: p_e1, p_e2, p_e3         ! 3D tracer T or W on parent grid
+      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), OPTIONAL, INTENT(out):: surf_crs, surf_msk_crs ! 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), ALLOCATABLE, DIMENSION(:,:,:) :: ze1, ze2, ze3, zsurf_crs, zsurf_msk_crs, zpmask  
+      !!----------------------------------------------------------------  
+      ! Initialize
+
+      ALLOCATE( ze1(jpi,jpj,jpk), ze2(jpi,jpj,jpk), ze3(jpi,jpj,jpk), zpmask(jpi,jpj,jpk) )
+      ALLOCATE( zsurf_crs(jpi_crs,jpj_crs,jpk), zsurf_msk_crs(jpi_crs,jpj_crs,jpk) )
+
+      ! Arrays, scalars initialization 
+      ze1(:,:,:)           = p_e1(:,:,:)
+      ze2(:,:,:)           = p_e2(:,:,:)
+      ze3(:,:,:)           = p_e3(:,:,:)
+      zsurf_crs(:,:,:)     = 0.0
+      zsurf_msk_crs(:,:,:) = 0.0
+      zpmask(:,:,:)        = p_mask(:,:,:)
+      ijpk                 = jpk
+
+      SELECT CASE ( cd_type )
+      
+         CASE ('W')
+
+            DO jk = 2 , ijpk
+            
+               DO ji = 1, jpi_crs                       ! ji = 1 et jpi_crs definit par cyclique est-ouest et pivot T
+                   ijie = mie_crs(ji)
+                   ijis = mis_crs(ji)
+                   jj   = 1
+                   ijje = mje_crs(jj) 
+                   ijjs = mjs_crs(jj)
+                   
+                   DO jii = ijis, ijie
+                      DO jjj = ijjs, ijje
+                         zsurf_crs(ji,jj,jk) = zsurf_crs(ji,jj,jk) + ze1(ji,jj,jk) * ze2(jii,jjj,jk)
+                         zsurf_msk_crs(ji,jj,jk) = zsurf_msk_crs(ji,jj,jk) + ( ze1(ji,jj,jk) * ze2(jii,jjj,jk) ) * zpmask(jii,jjj,jk-1)               
+                      ENDDO
+                   ENDDO
+                   
+                   zsurf_crs(ji,jj,jk) = zsurf_crs(ji,jj,jk) * 3
+                   
+                   DO jj = 2, jpj_crs                   ! jj = jpj_crs definit par pivot T 
+                      ijje = mje_crs(jj) 
+                      ijjs = mjs_crs(jj)  
+                   
+                      DO jii = ijis, ijie
+                         DO jjj = ijjs, ijje
+                            zsurf_crs(ji,jj,jk) = zsurf_crs(ji,jj,jk) + ze1(ji,jj,jk) * ze2(jii,jjj,jk)
+                            zsurf_msk_crs(ji,jj,jk) = zsurf_msk_crs(ji,jj,jk) + ( ze1(ji,jj,jk) * ze2(jii,jjj,jk) ) * zpmask(jii,jjj,jk-1)   
+                         ENDDO
+                      ENDDO
+                   ENDDO
+               ENDDO
+            ENDDO
+            
+            jk = 1                                      !cas particulier ou on est en surface
+            
+            DO ji = 1, jpi_crs                       ! ji = 1 et jpi_crs definit par cyclique est-ouest et pivot T
+               ijie = mie_crs(ji)
+               ijis = mis_crs(ji)
+               jj   = 1                           
+               ijje = mje_crs(jj) 
+               ijjs = mjs_crs(jj)
+                   
+               DO jii = ijis, ijie 
+                  DO jjj = ijjs, ijje
+                     zsurf_crs(ji,jj,jk) = zsurf_crs(ji,jj,jk) + ze1(ji,jj,jk) * ze2(jii,jjj,jk)
+                     zsurf_msk_crs(ji,jj,jk) = zsurf_msk_crs(ji,jj,jk) + ( ze1(ji,jj,jk) * ze2(jii,jjj,jk) ) * zpmask(jii,jjj,jk)               
+                  ENDDO  
+               ENDDO
+               zsurf_crs(ji,jj,jk) = zsurf_crs(ji,jj,jk) * 3                
+               DO jj = 2, jpj_crs                   ! jj = jpj_crs definit par pivot T 
+                  ijje = mje_crs(jj) 
+                  ijjs = mjs_crs(jj)              
+                  DO jii = ijis, ijie
+                     DO jjj = ijjs, ijje   
+                        zsurf_crs(ji,jj,jk) = zsurf_crs(ji,jj,jk) + ze1(ji,jj,jk) * ze2(jii,jjj,jk)
+                        zsurf_msk_crs(ji,jj,jk) = zsurf_msk_crs(ji,jj,jk) + ( ze1(ji,jj,jk) * ze2(jii,jjj,jk) ) * zpmask(jii,jjj,jk)   
+                     ENDDO   
+                  ENDDO
+               ENDDO
+            ENDDO
+            
+       CASE ('U')
+       
+          DO jk = 1 , ijpk
+            
+             DO ji = 1, jpi_crs                       ! ji = 1 et jpi_crs definit par cyclique est-ouest et pivot T
+                ijie = mie_crs(ji)
+                ijis = mis_crs(ji)
+                jj   = 1
+                ijje = mje_crs(jj)  
+                ijjs = mjs_crs(jj)
+                   
+                DO jii = ijis, ijie
+                   DO jjj = ijjs, ijje      
+                      zsurf_crs(ji,jj,jk) = zsurf_crs(ji,jj,jk) + ze3(ji,jj,jk) * ze2(jii,jjj,jk)
+                      zsurf_msk_crs(ji,jj,jk) = zsurf_msk_crs(ji,jj,jk) + ( ze3(ji,jj,jk) * ze2(jii,jjj,jk) ) * zpmask(jii,jjj,jk)               
+                   ENDDO
+                ENDDO
+                   
+                zsurf_crs(ji,jj,jk) = zsurf_crs(ji,jj,jk) * 3
+                   
+                DO jj = 2, jpj_crs                   ! jj = jpj_crs definit par pivot T 
+                   ijje = mje_crs(jj) 
+                   ijjs = mjs_crs(jj)  
+               
+                   DO jii = ijis, ijie
+                      DO jjj = ijjs, ijje
+                         zsurf_crs(ji,jj,jk) = zsurf_crs(ji,jj,jk) + ze3(ji,jj,jk) * ze2(jii,jjj,jk)
+                         zsurf_msk_crs(ji,jj,jk) = zsurf_msk_crs(ji,jj,jk) + ( ze3(ji,jj,jk) * ze2(jii,jjj,jk) ) * zpmask(jii,jjj,jk)
+                      ENDDO         
+                   ENDDO
+                ENDDO
+            ENDDO
+         ENDDO
+       
+      CASE ('V')
+          
+         DO jk = 1 , ijpk
+            
+            DO ji = 1, jpi_crs                       ! ji = 1 et jpi_crs definit par cyclique est-ouest et pivot T
+               ijie = mie_crs(ji)
+               ijis = mis_crs(ji)
+                   
+               DO jj = 1, jpj_crs                   ! jj = jpj_crs definit par pivot T 
+                  ijje = mje_crs(jj)   
+                  ijjs = mjs_crs(jj)  
+                   
+                  DO jii = ijis, ijie
+                     DO jjj = ijjs, ijje
+                        zsurf_crs(ji,jj,jk) = zsurf_crs(ji,jj,jk) + ze3(ji,jj,jk) * ze1(jii,jjj,jk) 
+                        zsurf_msk_crs(ji,jj,jk) = zsurf_msk_crs(ji,jj,jk) + ( ze3(ji,jj,jk) * ze1(jii,jjj,jk) ) * zpmask(jii,jjj,jk)   
+                     ENDDO
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDDO
+      END SELECT
+ 
+      surf_crs(:,:,:) = zsurf_crs(:,:,:)
+      CALL crs_lbc_lnk( cd_type, 1.0, pt3d1=surf_crs )
+         ! lbcnlk met la ligne jpj = 1 a 0 donc il faut la remettre en ne pas oubliant le cyclique est-ouest 
+         ! a faire un case pour cyclique est-ouest ou non.      
+         surf_crs(   :   ,1,:) = zsurf_crs(    :    ,1,:)   
+         surf_crs(   1   ,1,:) = zsurf_crs(jpi_crsm1,1,:)
+         surf_crs(jpi_crs,1,:) = zsurf_crs(    2    ,1,:)   
+         
+         surf_msk_crs(:,:,:) = zsurf_msk_crs(:,:,:)
+      CALL crs_lbc_lnk( cd_type, 1.0, pt3d1=surf_msk_crs )
+         ! lbcnlk met la ligne jpj = 1 a 0 donc il faut la remettre en ne pas oubliant le cyclique est-ouest        
+         surf_msk_crs(   :   ,1,:) = zsurf_msk_crs(    :    ,1,:)   
+         surf_msk_crs(   1   ,1,:) = zsurf_msk_crs(jpi_crsm1,1,:)
+         surf_msk_crs(jpi_crs,1,:) = zsurf_msk_crs(    2    ,1,:) 
+
+      DEALLOCATE( ze3 , ze2, ze1, zpmask )
+      DEALLOCATE( zsurf_msk_crs, zsurf_crs )
+ 
+
+   END SUBROUTINE crs_surf
+
 
 END MODULE crs
Index: /branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crs_dom.F90
===================================================================
--- /branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crs_dom.F90	(revision 3777)
+++ /branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crs_dom.F90	(revision 3778)
@@ -13,5 +13,6 @@
    PUBLIC
 
-   PUBLIC crs_dom_alloc  ! Called from crsini.F90
+   
+      PUBLIC crs_dom_alloc  ! Called from crsini.F90
    PUBLIC dom_grid_glo   
    PUBLIC dom_grid_crs   
@@ -30,4 +31,5 @@
       INTEGER  ::  npolj_full, npolj_crs        !: north fold mark
       INTEGER  ::  jpiglo_full, jpjglo_full     !: jpiglo / jpjglo
+      INTEGER  ::  npiglo, npjglo      !: jpjglo
       INTEGER  ::  nlci_full, nlcj_full         !: i-, j-dimension of local or sub domain on parent grid
       INTEGER  ::  nldi_full, nldj_full         !: starting indices of internal sub-domain on parent grid
@@ -59,5 +61,7 @@
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: fse3t_crs, fse3u_crs, fse3v_crs, fse3f_crs, fse3w_crs
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: fse3t_n_crs, fse3t_b_crs, fse3t_a_crs
-
+      
+      ! Surface
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: e1e2w_msk, e2e3u_msk, e1e3v_msk, e1e2w, e2e3u, e1e3v
                                                                   ! vertical scale factors 
       ! Coordinates
@@ -70,5 +74,5 @@
       ! Weights
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: facsurfv, facsurfu, facvol_t, facvol_w
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ocean_volume_crs_t, ocean_volume_crs_w
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ocean_volume_crs_t, ocean_volume_crs_w, bt_crs, r1_bt_crs
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: crs_surfu_wgt, crs_surfv_wgt, crs_surfw_wgt, crs_volt_wgt
 
@@ -192,10 +196,14 @@
          &      e3f_crs(jpi_crs,jpj_crs,jpk)    , fse3f_crs(jpi_crs,jpj_crs,jpk) , & 
          &      fse3t_b_crs(jpi_crs,jpj_crs,jpk), fse3t_n_crs(jpi_crs,jpj_crs,jpk),&
-         &      fse3t_a_crs(jpi_crs,jpj_crs,jpk), STAT=ierr(6))
+         &      fse3t_a_crs(jpi_crs,jpj_crs,jpk), e1e2w_msk(jpi_crs,jpj_crs,jpk) , &
+         &      e2e3u_msk(jpi_crs,jpj_crs,jpk)  , e1e3v_msk(jpi_crs,jpj_crs,jpk) , &
+         &      e1e2w(jpi_crs,jpj_crs,jpk)      , e2e3u(jpi_crs,jpj_crs,jpk)     , &
+         &      e1e3v(jpi_crs,jpj_crs,jpk)      , STAT=ierr(6))
 
 
       ALLOCATE( facsurfv(jpi_crs,jpj_crs,jpk) , facsurfu(jpi_crs,jpj_crs,jpk) , & 
          &      facvol_t(jpi_crs,jpj_crs,jpk) , facvol_w(jpi_crs,jpj_crs,jpk) , &
-         &      ocean_volume_crs_t(jpi_crs,jpj_crs,jpk) , ocean_volume_crs_w(jpi_crs,jpj_crs,jpk) , STAT=ierr(7))
+         &      ocean_volume_crs_t(jpi_crs,jpj_crs,jpk) , ocean_volume_crs_w(jpi_crs,jpj_crs,jpk), &
+         &      bt_crs(jpi_crs,jpj_crs,jpk)   , r1_bt_crs(jpi_crs,jpj_crs,jpk) , STAT=ierr(7))
 
 
Index: /branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crs_iom.F90
===================================================================
--- /branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crs_iom.F90	(revision 3777)
+++ /branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crs_iom.F90	(revision 3778)
@@ -179,5 +179,4 @@
       ELSEIF( PRESENT(pv_r3d) ) THEN   ;  CALL iom_put( cdvar, pv_r3d )
       ENDIF
-
       CALL dom_grid_glo   ! Return to parent grid domain
 
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 3777)
+++ /branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsdiawri.F90	(revision 3778)
@@ -103,9 +103,12 @@
       REAL(wp), POINTER, DIMENSION(:,:)   :: ze1e2u ! 2D workspace
       REAL(wp), POINTER, DIMENSION(:,:)   :: ze1e2v ! 2D workspace
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: z3dcrs ! 3D workspace for coarse grid
-      REAL(wp), POINTER, DIMENSION(:,:)   :: z2dcrs ! 2D workspace for coarse grid
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: z3dcrs, zw ! 3D workspace for coarse grid !cc
+      REAL(wp), POINTER, DIMENSION(:,:)   :: z2dcrs, ssh_crs2 ! 2D workspace for coarse grid
       INTEGER        :: ialloc, iiki, iikn
       INTEGER, SAVE  :: itsct
       REAL(wp)       :: zdtj, zrtsct
+!!cc1
+      REAL(wp)       :: z2dcrsu, z2dcrsv
+!!cc1
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), SAVE :: ztsnm
 
@@ -125,8 +128,8 @@
       REAL(wp)                                      :: zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0, zraur
       REAL(wp)                                      :: zij, zip1j, zijp1
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)         :: z2dcrsu, z2dcrsv, z2dcrsf
+!!cc1      REAL(wp), ALLOCATABLE, DIMENSION(:,:)         :: z2dcrsu, z2dcrsv, z2dcrsf, zhdivbt
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)         :: z2dcrsf, zhdivbt
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)         :: zsshub, zsshua, zsshvb, zsshva   ! temp work arrays for instantaneous fields
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)       :: zhdiv, zvolt_wgt, zrhd, zrhop, zavt
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)         :: zhdivbt
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)         :: zee_t, zee_f, zee_u, zee_v
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)       :: zmut, zmuf
@@ -150,6 +153,7 @@
       CALL wrk_alloc( jpi , jpj , jpk , zfse3u, zfse3v )
       CALL wrk_alloc( jpi , jpj , jpk , zfse3t, zfse3w )
-      CALL wrk_alloc( jpi_crs , jpj_crs , jpk , z3dcrs )
-      CALL wrk_alloc( jpi_crs , jpj_crs , z2dcrs )
+      CALL wrk_alloc( jpi_crs , jpj_crs , jpk , z3dcrs, zw ) !cc
+      CALL wrk_alloc( jpi_crs , jpj_crs , z2dcrs, ssh_crs2 )
+      
 
       IF ( .NOT. ALLOCATED(ztsnm) ) THEN
@@ -187,6 +191,8 @@
       ENDIF
 
-      IF ( .NOT. ALLOCATED(z2dcrsu) ) THEN
-         ALLOCATE( z2dcrsu(jpi_crs,jpj_crs) , z2dcrsv(jpi_crs,jpj_crs) , &
+      IF ( .NOT. ALLOCATED(z2dcrsf) ) THEN
+!!cc1      IF ( .NOT. ALLOCATED(z2dcrsu) ) THEN
+!!cc1         ALLOCATE( z2dcrsu(jpi_crs,jpj_crs) , z2dcrsv(jpi_crs,jpj_crs) , &
+         ALLOCATE(   &
             &      z2dcrsf(jpi_crs,jpj_crs) , zhdivbt(jpi_crs,jpj_crs) , &
             &      zsshub(jpi_crs,jpj_crs)  , zsshvb(jpi_crs,jpj_crs)  , &
@@ -213,5 +219,5 @@
       ENDIF
 
-
+      zw(:,:,:)=0.0
       ! generic work arrays
       z2d(:,:)      = 0.0
@@ -320,5 +326,5 @@
       z3d(:,:,:) = tsn(:,:,:,jp_tem)
       CALL crsfun( p_e1e2t=e1e2t, cd_type='T', cd_op='VOL', p_cmask=tmask_crs, p_ptmask=tmask, &
-         &         p_pfield3d_1=zfse3t, p_pfield3d_2=z3d, p_cfield3d=z3dcrs )
+         &         p_pfield3d_1=zfse3t, p_pfield3d_2=z3d, p_cfield3d=z3dcrs ) 
       ztsn(:,:,:,1) = z3dcrs(:,:,:)
       ztsnm(:,:,:,1) = ztsnm(:,:,:,1) + ztsn(:,:,:,1) 
@@ -392,4 +398,12 @@
          &         p_fse3=zfse3v, p_pfield=z3d1, p_cfield3d=z3dcrs )
       zvsm(:,:,:) = zvsm(:,:,:) + z3dcrs(:,:,:)
+      
+      
+      ! Vitesse vertical !cc
+      z3dcrs(:,:,:) = 0.0
+      CALL crsfun( p_e1e2t=e1e2t, cd_type='T', cd_op='ARE', p_cmask=tmask_crs, p_ptmask=tmask, &
+        &           p_pfield3d_1=wn, p_cfield3d=z3dcrs)
+      zw(:,:,:) = z3dcrs(:,:,:)
+
 
       ! 2.5. Surface boundary conditions
@@ -495,5 +509,5 @@
       zhdivbt(:,:) = 0.0
       DO jk = 1, jpkm1
-         z2dcrsu(:,:) = 0.0; z2dcrsv(:,:) = 0.0   
+!!cc1         z2dcrsu(:,:) = 0.0; z2dcrsv(:,:) = 0.0   
          DO ji = 2, jpi_crsm1
             DO jj = 2, jpj_crsm1
@@ -501,10 +515,13 @@
                ! Horizontal divergence ( following OPA_SRC/DYN/divcur.F90 ) 
                ! with partial steps and/or variable layer thicknesses for W
-               z2dcrsu(ji,jj) =  ( zum(ji,  jj,jk) * crs_surfu_wgt(ji,  jj,jk)  ) &
+!!cc1               z2dcrsu(ji,jj) =  ( zum(ji,  jj,jk) * crs_surfu_wgt(ji,  jj,jk)  ) &
+               z2dcrsu =  ( zum(ji,  jj,jk) * crs_surfu_wgt(ji,  jj,jk)  ) &
                   &            - ( zum(ji-1,jj,jk) * crs_surfu_wgt(ji-1,jj,jk)  )
-               z2dcrsv(ji,jj) =  ( zvm(ji,jj,  jk) * crs_surfv_wgt(ji,jj  ,jk)  ) &
+!!cc1               z2dcrsv(ji,jj) =  ( zvm(ji,jj,  jk) * crs_surfv_wgt(ji,jj  ,jk)  ) &
+               z2dcrsv =  ( zvm(ji,jj,  jk) * crs_surfv_wgt(ji,jj  ,jk)  ) &
                   &            - ( zvm(ji,jj-1,jk) * crs_surfv_wgt(ji,jj-1,jk)  )
 
-               zhdiv(ji,jj,jk) = ( z2dcrsu(ji,jj) + z2dcrsv(ji,jj) ) * zvolt_wgt(ji,jj,jk) 
+!!cc1               zhdiv(ji,jj,jk) = ( z2dcrsu(ji,jj) + z2dcrsv(ji,jj) ) * zvolt_wgt(ji,jj,jk) 
+               zhdiv(ji,jj,jk) = ( z2dcrsu + z2dcrsv ) * zvolt_wgt(ji,jj,jk) 
  
             ENDDO
@@ -523,7 +540,61 @@
       zhdivbtm(:,:) = zhdivbtm(:,:) + zhdivbt
 
-
-      !
+!      !!! Calcul de l'energie cinétique   !cc  !! DECLARE LES VARIABLES
+!      
+!      zun2(:,:,:) = un(:,:,:) * un(:,:,:)    ! U^2
+!      zvn2(:,:,:) = vn(:,:,:) * vn(:,:,:)    ! V^2
+!      
+!!      ! moyenne sur i de U^2
+!      
+!      DO ji = 1, jpiglo-1
+!         zun2(ji,:,:) = 0.5 * (zun2(ji,:,:) + zun2(ji+1,:,:))
+!      END
+!      ji = jpiglo
+!      zun2(ji,:,:) = 0.5 * zun2(ji,:,:) 
+!      uun2(:,:,:) = zun2(:,:,:)
+!      
+!      
+!         CALL crs_lbc_lnk( cd_type='T', 1.0, pt3d1=uun2)
+!         ! lbcnlk met la ligne jpj = 1 a 0 donc il faut la remettre en ne pas oubliant le cyclique est-ouest 
+!         ! a faire un case pour cyclique est-ouest ou non.      
+!         uun2(   :   ,1,:) = zun2(    :    ,1,:)   
+!         
+!         IF ((jperio==4) .OR. (jperio==6)) THEN
+!            uun2(   1   ,1,:) = zun2(jpi_crsm1,1,:)
+!            uun2(jpi_crs,1,:) = zun2(    2    ,1,:)
+!         ENDIF  
+!
+!
+ !     DO jj = 1, jpjglo-1
+ !        zvn2(:,jj,:) = 0.5 * (zvn2(:,jj,:) + zvn2(:,jj+1,:))
+ !     END
+ !     jj = jpjglo 
+ !     zvn2(:,jj,:) = 0.5 * zvn2(:,jj,:)
+ !     vvn2(:,:,:) = zvn2(:,:,:)
+ !     
+ !     CALL crs_lbc_lnk( cd_type='T', 1.0, pt3d1=vvn2)
+ !        ! lbcnlk met la ligne jpj = 1 a 0 donc il faut la remettre en ne pas oubliant le cyclique est-ouest 
+ !        ! a faire un case pour cyclique est-ouest ou non.      
+ !        vvn2(   :   ,1,:) = zvn2(    :    ,1,:)   
+ !        
+ !        IF ((jperio==4) .OR. (jperio==6)) THEN
+ !           vvn2(   1   ,1,:) = zvn2(jpi_crsm1,1,:)
+ !           vvn2(jpi_crs,1,:) = zvn2(    2    ,1,:)
+ !        ENDIF  
+
+
+
+      
       !    2.6.3. Sea-surface Height  ( See DOM/istate.F90: ssh init; OPA_SRC/DYN/sshwzv.F90: ssh update )
+      !cc
+      z2dcrs(:,:) = 0.0
+      CALL crsfun( p_e1e2t=e1e2t, cd_type='T', cd_op='ARE', p_cmask=tmask_crs, p_ptmask=tmask, p_pfield2d=sshn, &
+      &                  p_cfield2d=z2dcrs )
+      ssh_crs2(:,:) = z2dcrs(:,:)
+  !    WRITE(numout,*) 'test', sshn(:,:)
+  !    WRITE(numout,*) 'test', ssh_crs2(:,:)
+      
+     !cc    
+      
       !       set some temp variables
       z2dcrs(:,:) = 0.0; z2dcrsf(:,:) = 0.0
@@ -696,4 +767,5 @@
          WRITE(numout,*) 'crsdiawri.', clmode
          zrtsct = 1.0/REAL(itsct, wp)
+
          !
          !    3.1.2 Weights for spatial averages
@@ -710,4 +782,5 @@
          CALL crs_iom_put( "toce_crs"   , pv_r3d=tsn_crs(:,:,:,1)  )    ! temperature
          CALL crs_iom_put( "sst_crs"    , pv_r2d=z2dcrs            )    ! sst
+     
          !
          !    3.1.4 Salinity
@@ -725,5 +798,5 @@
          us_crs(:,:,:) = zusm(:,:,:) * zsuru(:,:,:) ! area-weighted- , time- mean      
          CALL crs_iom_put( "uoces_crs"     , pv_r3d=us_crs             )   ! uS
-
+              
          !    3.1.6 V-velocity
          vn_crs(:,:,:) = zvm(:,:,:) * zrtsct ! area-weighted- , time- mean 
@@ -737,8 +810,11 @@
          wn_crs(:,:,:) = zwm(:,:,:) * zrtsct ! area-weighted- , time- mean 
          CALL crs_iom_put( "woce_crs"   , pv_r3d=wn_crs                   )    ! W-velocity
+         CALL crs_iom_put( "woce2_crs"   , pv_r3d=zw                      )    ! cc
 
          !    3.1.8 Horizontal divergence
          hdivn_crs(:,:,:) = zhdivnm(:,:,:) * zrtsct
-         CALL crs_iom_put( "hdivn_crs"   , pv_r3d=hdivn_crs            )    ! hdiv
+      !   CALL crs_iom_put( "hdivn_crs"   , pv_r3d=hdivn_crs            )   ! hdiv
+         CALL crs_iom_put( "hdiv_crs"   , pv_r3d=hdivn_crs            )  
+
          hdivbt_crs(:,:) = zhdivbtm(:,:) * zrtsct
 
@@ -758,5 +834,5 @@
          ENDIF
          CALL crs_iom_put( "ssh_crs"   , pv_r2d=sshn_crs                  )   ! ssh output 
-
+         CALL crs_iom_put( "ssh2_crs"  , pv_r2d=ssh_crs2                  )   !cc
 
          !    3.1.10 Potential density
@@ -866,8 +942,9 @@
       CALL wrk_dealloc( jpi , jpj , jpk , zfse3u, zfse3v )
       CALL wrk_dealloc( jpi , jpj , jpk , zfse3t, zfse3w )
-      CALL wrk_dealloc( jpi_crs , jpj_crs , jpk , z3dcrs )
-      CALL wrk_dealloc( jpi_crs , jpj_crs , z2dcrs )
-
-      DEALLOCATE( z2dcrsu, z2dcrsv, z2dcrsf, zhdivbt )
+      CALL wrk_dealloc( jpi_crs , jpj_crs , jpk , z3dcrs, zw ) !cc
+      CALL wrk_dealloc( jpi_crs , jpj_crs , z2dcrs, ssh_crs2 )
+
+!!cc1      DEALLOCATE( z2dcrsu, z2dcrsv, z2dcrsf, zhdivbt ) ! probleme de malloc au 65 eme pas de temps
+      DEALLOCATE( z2dcrsf, zhdivbt ) ! probleme de malloc au 65 eme pas de temps
       DEALLOCATE( zsshub, zsshua, zsshvb, zsshva )
       DEALLOCATE( zee_t, zee_f, zee_u, zee_v )
@@ -876,4 +953,5 @@
       DEALLOCATE( zhdiv, zvolt_wgt) 
       !
+                        
       IF( nn_timing == 1 )   CALL timing_stop('crs_dia_wri')
       !
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 3777)
+++ /branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsdomwri.F90	(revision 3778)
@@ -137,23 +137,23 @@
       !========================================================
       !                                                         ! horizontal mesh (inum3)
-      CALL crs_iom_rstput( 0, 0, inum3, 'glamt_crs', pv_r2d=glamt_crs, ktype = jp_r4 )     !    ! latitude
-      CALL crs_iom_rstput( 0, 0, inum3, 'glamu_crs', pv_r2d=glamu_crs, ktype = jp_r4 )
-      CALL crs_iom_rstput( 0, 0, inum3, 'glamv_crs', pv_r2d=glamv_crs, ktype = jp_r4 )
-      CALL crs_iom_rstput( 0, 0, inum3, 'glamf_crs', pv_r2d=glamf_crs, ktype = jp_r4 )
-      
-      CALL crs_iom_rstput( 0, 0, inum3, 'gphit_crs', pv_r2d=gphit_crs, ktype = jp_r4 )     !    ! longitude
-      CALL crs_iom_rstput( 0, 0, inum3, 'gphiu_crs', pv_r2d=gphiu_crs, ktype = jp_r4 )
-      CALL crs_iom_rstput( 0, 0, inum3, 'gphiv_crs', pv_r2d=gphiv_crs, ktype = jp_r4 )
-      CALL crs_iom_rstput( 0, 0, inum3, 'gphif_crs', pv_r2d=gphif_crs, ktype = jp_r4 )
-      
-      CALL crs_iom_rstput( 0, 0, inum3, 'e1t_crs', pv_r2d=e1t_crs, ktype = jp_r8 )         !    ! e1 scale factors
-      CALL crs_iom_rstput( 0, 0, inum3, 'e1u_crs', pv_r2d=e1u_crs, ktype = jp_r8 )
-      CALL crs_iom_rstput( 0, 0, inum3, 'e1v_crs', pv_r2d=e1v_crs, ktype = jp_r8 )
-      CALL crs_iom_rstput( 0, 0, inum3, 'e1f_crs', pv_r2d=e1f_crs, ktype = jp_r8 )
-      
-      CALL crs_iom_rstput( 0, 0, inum3, 'e2t_crs', pv_r2d=e2t_crs, ktype = jp_r8 )         !    ! e2 scale factors
-      CALL crs_iom_rstput( 0, 0, inum3, 'e2u_crs', pv_r2d=e2u_crs, ktype = jp_r8 )
-      CALL crs_iom_rstput( 0, 0, inum3, 'e2v_crs', pv_r2d=e2v_crs, ktype = jp_r8 )
-      CALL crs_iom_rstput( 0, 0, inum3, 'e2f_crs', pv_r2d=e2f_crs, ktype = jp_r8 )
+      CALL crs_iom_rstput( 0, 0, inum3, 'glamt', pv_r2d=glamt_crs, ktype = jp_r4 )     !    ! latitude
+      CALL crs_iom_rstput( 0, 0, inum3, 'glamu', pv_r2d=glamu_crs, ktype = jp_r4 )
+      CALL crs_iom_rstput( 0, 0, inum3, 'glamv', pv_r2d=glamv_crs, ktype = jp_r4 )
+      CALL crs_iom_rstput( 0, 0, inum3, 'glamf', pv_r2d=glamf_crs, ktype = jp_r4 )
+      
+      CALL crs_iom_rstput( 0, 0, inum3, 'gphit', pv_r2d=gphit_crs, ktype = jp_r4 )     !    ! longitude
+      CALL crs_iom_rstput( 0, 0, inum3, 'gphiu', pv_r2d=gphiu_crs, ktype = jp_r4 )
+      CALL crs_iom_rstput( 0, 0, inum3, 'gphiv', pv_r2d=gphiv_crs, ktype = jp_r4 )
+      CALL crs_iom_rstput( 0, 0, inum3, 'gphif', pv_r2d=gphif_crs, ktype = jp_r4 )
+      
+      CALL crs_iom_rstput( 0, 0, inum3, 'e1t', pv_r2d=e1t_crs, ktype = jp_r8 )         !    ! e1 scale factors
+      CALL crs_iom_rstput( 0, 0, inum3, 'e1u', pv_r2d=e1u_crs, ktype = jp_r8 )
+      CALL crs_iom_rstput( 0, 0, inum3, 'e1v', pv_r2d=e1v_crs, ktype = jp_r8 )
+      CALL crs_iom_rstput( 0, 0, inum3, 'e1f', pv_r2d=e1f_crs, ktype = jp_r8 )
+      
+      CALL crs_iom_rstput( 0, 0, inum3, 'e2t', pv_r2d=e2t_crs, ktype = jp_r8 )         !    ! e2 scale factors
+      CALL crs_iom_rstput( 0, 0, inum3, 'e2u', pv_r2d=e2u_crs, ktype = jp_r8 )
+      CALL crs_iom_rstput( 0, 0, inum3, 'e2v', pv_r2d=e2v_crs, ktype = jp_r8 )
+      CALL crs_iom_rstput( 0, 0, inum3, 'e2f', pv_r2d=e2f_crs, ktype = jp_r8 )
       
       CALL crs_iom_rstput( 0, 0, inum3, 'ff_crs', pv_r2d=ff_crs, ktype = jp_r8 )           !    ! coriolis factor
@@ -167,9 +167,10 @@
       IF( ln_zps ) THEN                       ! z-coordinate - partial steps
 
+            
          IF ( nn_msh_crs <= 6 ) THEN
-            CALL crs_iom_rstput( 0, 0, inum4, 'e3t_crs', pv_r3d=fse3t_crs )      
-            CALL crs_iom_rstput( 0, 0, inum4, 'e3w_crs', pv_r3d=fse3w_crs )      
-            CALL crs_iom_rstput( 0, 0, inum4, 'e3u_crs', pv_r3d=fse3u_crs )      
-            CALL crs_iom_rstput( 0, 0, inum4, 'e3v_crs', pv_r3d=fse3v_crs )      
+            CALL crs_iom_rstput( 0, 0, inum4, 'e3t', pv_r3d=e3t_crs )      
+            CALL crs_iom_rstput( 0, 0, inum4, 'e3w', pv_r3d=e3w_crs )      
+            CALL crs_iom_rstput( 0, 0, inum4, 'e3u', pv_r3d=e3u_crs )      
+            CALL crs_iom_rstput( 0, 0, inum4, 'e3v', pv_r3d=e3v_crs )      
          ELSE
             DO jj = 1,jpj_crs   
@@ -188,5 +189,5 @@
 
          IF ( nn_msh_crs <= 3 ) THEN
-            CALL crs_iom_rstput( 0, 0, inum4, 'gdept_crs', pv_r3d=gdept_crs, ktype = jp_r4 )     
+            CALL crs_iom_rstput( 0, 0, inum4, 'gdept_crs', pv_r3d=gdept_crs, ktype = jp_r4 ) 
             DO jk = 1,jpk   
                DO jj = 1, jpj_crsm1   
@@ -219,8 +220,16 @@
 
          CALL crs_iom_rstput(  0, 0, inum4, 'ocean_volume_crs_t', pv_r3d=ocean_volume_crs_t ) 
-         CALL crs_iom_rstput(  0, 0, inum4, 'facvol_t', pv_r3d=facvol_t ) 
-         CALL crs_iom_rstput(  0, 0, inum4, 'facvol_w', pv_r3d=facvol_w ) 
-         CALL crs_iom_rstput(  0, 0, inum4, 'facsurfu', pv_r3d=facsurfu ) 
-         CALL crs_iom_rstput(  0, 0, inum4, 'facsurfv', pv_r3d=facsurfv ) 
+         CALL crs_iom_rstput(  0, 0, inum4, 'facvol_t' , pv_r3d=facvol_t ) 
+         CALL crs_iom_rstput(  0, 0, inum4, 'facvol_w' , pv_r3d=facvol_w ) 
+         CALL crs_iom_rstput(  0, 0, inum4, 'facsurfu' , pv_r3d=facsurfu ) 
+         CALL crs_iom_rstput(  0, 0, inum4, 'facsurfv' , pv_r3d=facsurfv ) 
+         CALL crs_iom_rstput(  0, 0, inum4, 'e1e2w_msk', pv_r3d=e1e2w_msk) 
+         CALL crs_iom_rstput(  0, 0, inum4, 'e2e3u_msk', pv_r3d=e2e3u_msk) 
+         CALL crs_iom_rstput(  0, 0, inum4, 'e1e3v_msk', pv_r3d=e1e3v_msk)
+         CALL crs_iom_rstput(  0, 0, inum4, 'e1e2w'    , pv_r3d=e1e2w    ) 
+         CALL crs_iom_rstput(  0, 0, inum4, 'e2e3u'    , pv_r3d=e2e3u    ) 
+         CALL crs_iom_rstput(  0, 0, inum4, 'e1e3v'    , pv_r3d=e1e3v    )
+         CALL crs_iom_rstput(  0, 0, inum4, 'bt_crs'   , pv_r3d=bt_crs   )
+         CALL crs_iom_rstput(  0, 0, inum4, 'r1_bt_crs', pv_r3d=r1_bt_crs)
 
          CALL crs_iom_rstput(  0, 0, inum4, 'crs_surfu_wgt', pv_r3d=crs_surfu_wgt) 
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 3777)
+++ /branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsini.F90	(revision 3778)
@@ -193,9 +193,10 @@
 
             DO ji = 2, jpiglo_crsm1
-               ijie = (ji*nn_factx)-nn_factx+1
+               !cc ijie = (ji*nn_factx)-nn_factx+1
+               ijie = (ji*nn_factx)-nn_factx   !cc
                ijis = ijie-nn_factx+1
 
                IF ( ji == jpiglo_crsm1 ) THEN
-                  IF ( ((jpiglo-1)-ijie) <= nn_factx )   ijie = jpiglo-1
+                  IF ( ((jpiglo-1)-ijie) <= nn_factx ) ijie = jpiglo-2  ! ijie = jpiglo-1 !cc
                ENDIF
 
@@ -207,6 +208,7 @@
                   ENDIF
 
-                  DO jj = 2, jpjglo_crsm1
-                     ijje = ijjgloT-nn_facty*(jj-2)
+                 DO jj = 2, jpjglo_crsm1
+                ! cc ijje = ijjgloT-nn_facty*(jj-2)
+                     ijje = ijjgloT-nn_facty*(jj-2) - 1
                      ijjs = ijje-nn_facty+1                   
                   
@@ -271,8 +273,8 @@
 
         ! Pad the boundaries, do not know if it is necessary
-         mis_crs(1) = 1; mis_crs(jpiglo_crs) = jpiglo
-         mie_crs(1) = 1; mie_crs(jpiglo_crs) = jpiglo
-         mjs_crs(1) = 1; mjs_crs(jpjglo_crs) = jpjglo
-         mje_crs(1) = 1; mje_crs(jpjglo_crs) = jpjglo 
+         mis_crs(1) = 1           ; mis_crs(jpiglo_crs) = mie_crs(jpiglo_crs - 1) + 1    !cc
+         mie_crs(1) = nn_factx    ; mie_crs(jpiglo_crs) = jpiglo                         !cc
+         mjs_crs(1) = 1           ; mjs_crs(jpjglo_crs) = mje_crs(jpjglo_crs - 1) + 1
+         mje_crs(1) = mjs_crs(2)-1; mje_crs(jpjglo_crs) = jpjglo 
 
 !         WRITE(numout,*) 'crs_init. coarse grid bounds on parent grid'
@@ -313,17 +315,17 @@
         WRITE(numout,*) 'crsini. count 1'
 
-!        CALL crsfun( gphiu, glamu, 'U', gphiu_crs, glamu_crs )
-!        WRITE(numout,*) 'crsini. gphiu_crs(15,15)', gphiu_crs(15,15)
-!        WRITE(numout,*) 'crsini. glamu_crs(15,15)', glamu_crs(15,15)
+        CALL crsfun( gphiu, glamu, 'U', gphiu_crs, glamu_crs )       !cc
+        WRITE(numout,*) 'crsini. gphiu_crs(15,15)', gphiu_crs(15,15) !cc
+        WRITE(numout,*) 'crsini. glamu_crs(15,15)', glamu_crs(15,15) !cc
         WRITE(numout,*) 'crsini. count 2'
  
-!        CALL crsfun( p_pgphi=gphiv, p_pglam=glamv, cd_type='V', p_cgphi=gphiv_crs, p_cglam=glamv_crs )
-!        WRITE(numout,*) 'crsini. gphiv_crs(15,15)', gphiv_crs(15,15)
-!        WRITE(numout,*) 'crsini. glamv_crs(15,15)', glamv_crs(15,15)
+        CALL crsfun( p_pgphi=gphiv, p_pglam=glamv, cd_type='V', p_cgphi=gphiv_crs, p_cglam=glamv_crs ) !cc
+        WRITE(numout,*) 'crsini. gphiv_crs(15,15)', gphiv_crs(15,15) !cc
+        WRITE(numout,*) 'crsini. glamv_crs(15,15)', glamv_crs(15,15) !cc
 
         WRITE(numout,*) 'crsini. count 3'
-!        CALL crsfun( p_pgphi=gphif, p_pglam=glamf, cd_type='F', p_cgphi=gphif_crs, p_cglam=glamf_crs )
-!        WRITE(numout,*) 'crsini. gphif_crs(15,15)', gphif_crs(15,15)
-!        WRITE(numout,*) 'crsini. glamf_crs(15,15)', glamf_crs(15,15)
+        CALL crsfun( p_pgphi=gphif, p_pglam=glamf, cd_type='F', p_cgphi=gphif_crs, p_cglam=glamf_crs ) !cc
+        WRITE(numout,*) 'crsini. gphif_crs(15,15)', gphif_crs(15,15) !cc
+        WRITE(numout,*) 'crsini. glamf_crs(15,15)', glamf_crs(15,15) !cc
 
         WRITE(numout,*) 'crsini. count 4'
@@ -350,8 +352,12 @@
 
      !      3.c.1 Horizontal scale factors
-     CALL crsfun( cd_type='T', cd_op='SUM', p_pmask=tmask, p_e1=e1t, p_e2=e2t, p_cfield2d_1=e1t_crs, p_cfield2d_2=e2t_crs )
-     CALL crsfun( cd_type='U', cd_op='SUM', p_pmask=umask, p_e1=e1u, p_e2=e2u, p_cfield2d_1=e1u_crs, p_cfield2d_2=e2u_crs )
-     CALL crsfun( cd_type='V', cd_op='SUM', p_pmask=vmask, p_e1=e1v, p_e2=e2v, p_cfield2d_1=e1v_crs, p_cfield2d_2=e2v_crs )
-     CALL crsfun( cd_type='F', cd_op='SUM', p_pmask=fmask, p_e1=e1f, p_e2=e2f, p_cfield2d_1=e1f_crs, p_cfield2d_2=e2f_crs )
+  !   CALL crsfun( cd_type='T', cd_op='SUM', p_pmask=tmask, p_e1=e1t, p_e2=e2t, p_cfield2d_1=e1t_crs, p_cfield2d_2=e2t_crs )
+  !   CALL crsfun( cd_type='U', cd_op='SUM', p_pmask=umask, p_e1=e1u, p_e2=e2u, p_cfield2d_1=e1u_crs, p_cfield2d_2=e2u_crs )
+  !   CALL crsfun( cd_type='V', cd_op='SUM', p_pmask=vmask, p_e1=e1v, p_e2=e2v, p_cfield2d_1=e1v_crs, p_cfield2d_2=e2v_crs )
+  !   CALL crsfun( cd_type='F', cd_op='SUM', p_pmask=fmask, p_e1=e1f, p_e2=e2f, p_cfield2d_1=e1f_crs, p_cfield2d_2=e2f_crs )
+     CALL crsfun( cd_type='T', cd_op='POS', p_pmask=tmask, p_e1=e1t, p_e2=e2t, p_cfield2d_1=e1t_crs, p_cfield2d_2=e2t_crs )
+     CALL crsfun( cd_type='U', cd_op='POS', p_pmask=umask, p_e1=e1u, p_e2=e2u, p_cfield2d_1=e1u_crs, p_cfield2d_2=e2u_crs )
+     CALL crsfun( cd_type='V', cd_op='POS', p_pmask=vmask, p_e1=e1v, p_e2=e2v, p_cfield2d_1=e1v_crs, p_cfield2d_2=e2v_crs )
+     CALL crsfun( cd_type='F', cd_op='POS', p_pmask=fmask, p_e1=e1f, p_e2=e2f, p_cfield2d_1=e1f_crs, p_cfield2d_2=e2f_crs )
 
      e1e2t_crs(:,:) = e1t_crs(:,:) * e2t_crs(:,:)
@@ -440,11 +446,15 @@
      zfse3f(:,:,:) = fse3f(:,:,:)
      zfse3w(:,:,:) = fse3w(:,:,:)
-
-     CALL crsfun( p_e1e2t=e1e2t, cd_type='T', cd_op='MAX', p_cmask=tmask_crs, p_ptmask=tmask, p_pfield3d_1=zfse3t, p_cfield3d=e3t_crs )
-     CALL crsfun( p_e1e2t=e1e2t, cd_type='W', cd_op='MAX', p_cmask=tmask_crs, p_ptmask=tmask, p_pfield3d_1=zfse3w, p_cfield3d=e3w_crs )
-     CALL crsfun( p_e1e2t=e1e2t, cd_type='U', cd_op='MIN', p_cmask=umask_crs, p_ptmask=umask, p_pfield3d_1=zfse3u, p_cfield3d=e3u_crs )
-     CALL crsfun( p_e1e2t=e1e2t, cd_type='V', cd_op='MIN', p_cmask=vmask_crs, p_ptmask=vmask, p_pfield3d_1=zfse3v, p_cfield3d=e3v_crs )
-     CALL crsfun( p_e1e2t=e1e2t, cd_type='F', cd_op='MIN', p_cmask=fmask_crs, p_ptmask=fmask, p_pfield3d_1=zfse3f, p_cfield3d=e3f_crs )
-
+     
+      
+
+     !CALL crsfun( p_e1e2t=e1e2t, cd_type='T', cd_op='MAX', p_cmask=tmask_crs, p_ptmask=tmask, p_pfield3d_1=zfse3t, p_cfield3d=e3t_crs )
+     !CALL crsfun( p_e1e2t=e1e2t, cd_type='W', cd_op='MAX', p_cmask=tmask_crs, p_ptmask=tmask, p_pfield3d_1=zfse3w, p_cfield3d=e3w_crs )
+     !CALL crsfun( p_e1e2t=e1e2t, cd_type='U', cd_op='MIN', p_cmask=umask_crs, p_ptmask=umask, p_pfield3d_1=zfse3u, p_cfield3d=e3u_crs )
+     !CALL crsfun( p_e1e2t=e1e2t, cd_type='V', cd_op='MIN', p_cmask=vmask_crs, p_ptmask=vmask, p_pfield3d_1=zfse3v, p_cfield3d=e3v_crs )
+     !CALL crsfun( p_e1e2t=e1e2t, cd_type='F', cd_op='MIN', p_cmask=fmask_crs, p_ptmask=fmask, p_pfield3d_1=zfse3f, p_cfield3d=e3f_crs )
+     CALL crs_e3_max( p_e3=zfse3t, cd_type='T', p_mask=tmask, p_e3_crs=e3t_crs)
+     CALL crs_e3_max( p_e3=zfse3w, cd_type='W', p_mask=tmask, p_e3_crs=e3w_crs)
+    
      ! Reset 0 to e3t_0 or e3w_0
      DO jk = 1, jpk
@@ -465,5 +475,9 @@
      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 )
 
-
+     !    3.d.4   Surfaces 
+     
+     CALL crs_surf(p_e1=e1t, p_e2=e2t ,p_e3=zfse3w, cd_type='W', p_mask=tmask, surf_crs=e1e2w, surf_msk_crs=e1e2w_msk)
+     CALL crs_surf(p_e1=e1u, p_e2=e2u ,p_e3=zfse3u, cd_type='U', p_mask=umask, surf_crs=e2e3u, surf_msk_crs=e2e3u_msk)
+     CALL crs_surf(p_e1=e1v, p_e2=e2v ,p_e3=zfse3v, cd_type='V', p_mask=vmask, surf_crs=e1e3v, surf_msk_crs=e1e3v_msk)
 
 
@@ -476,4 +490,8 @@
       CALL crsfun( 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( cd_type='W', cd_op='VOL', p_pmask=tmask, p_e1=e1t, p_e2=e2t, p_fse3=zfse3w, &
@@ -495,5 +513,4 @@
      ! 5.  Write out coarse meshmask  (see OPA_SRC/DOM/domwri.F90 for ideas later)
      !---------------------------------------------------------
-     
      IF ( nn_msh_crs > 0 ) CALL crs_dom_wri
 
