Index: /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90
===================================================================
--- /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90	(revision 3962)
+++ /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90	(revision 3963)
@@ -127,10 +127,10 @@
       CALL iom_put( 'ibgfrcemp',frc_vol * 1.e-9 )           !vol - forcing : km3 (equivalent liquid)         
       CALL iom_put( 'ibgfrcemps',frc_sal * 1.e-9 )          !sal - forcing : psu*km3   
-      CALL iom_put( 'ibgemps',zbg_emps *86400)            
-      CALL iom_put( 'ibgemp',zbg_emp *86400)              
-      CALL iom_put( 'ibgfsbri',zbg_fsbri *86400)         
-      CALL iom_put( 'ibgfseqv',zbg_fseqv *86400)         
-      CALL iom_put( 'ibgfsaltres',zbg_fsalt_res *86400)         
-      CALL iom_put( 'ibgfsaltrpo',zbg_fsalt_rpo *86400)         
+      CALL iom_put( 'ibgemps',zbg_emps *86400.)            
+      CALL iom_put( 'ibgemp',zbg_emp *86400.)              
+      CALL iom_put( 'ibgfsbri',zbg_fsbri *86400.)         
+      CALL iom_put( 'ibgfseqv',zbg_fseqv *86400.)         
+      CALL iom_put( 'ibgfsaltres',zbg_fsalt_res *86400.)         
+      CALL iom_put( 'ibgfsaltrpo',zbg_fsalt_rpo *86400.)         
       CALL iom_put( 'ibggrpme',bg_grpme * rhoic / rau0 * 1.e-9 ) ! km3 (equivalent liquid)         
       !
Index: /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
===================================================================
--- /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 3962)
+++ /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 3963)
@@ -66,4 +66,5 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   zu_io, zv_io   ! ice-ocean velocity
       REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset)
+      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)
      !!---------------------------------------------------------------------
 
@@ -234,9 +235,14 @@
          zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic )
 
+         zchk_vmin = glob_min(v_i)
+         zchk_amax = glob_max(SUM(a_i,dim=3))
+         zchk_amin = glob_min(a_i)
+
          IF(lwp) THEN
-            IF (    ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limdyn) = ',(zchk_v_i * 86400.)
-            IF (    ABS( zchk_smv   ) >  1.e-4  ) WRITE(numout,*) 'violation saline [psu*m3/day] (limdyn) = ',(zchk_smv * 86400.)
-            IF ( MINVAL( v_i(:,:,:) ) <  0.    ) WRITE(numout,*) 'violation v_i<0  [mm]         (limdyn) = ',(MINVAL(v_i) * 1.e-3)
-            IF ( MAXVAL( SUM(a_i(:,:,:),dim=3) ) > amax+1.e-10 ) WRITE(numout,*) 'violation a_i>amax    (limdyn) = ',MAXVAL(SUM(a_i,dim=3))
+            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limdyn) = ',(zchk_v_i * 86400.)
+            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limdyn) = ',(zchk_smv * 86400.)
+            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limdyn) = ',(zchk_vmin * 1.e-3)
+            !IF ( zchk_amax >  amax+1.e-10   ) WRITE(numout,*) 'violation a_i>amax            (limdyn) = ',zchk_amax
+            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limdyn) = ',zchk_amin
          ENDIF
          !CALL ctl_opn( numhsb , cl_name , 'UNKNOWN' , 'FORMATTED' , 'SEQUENTIAL' , 1 , numout , lwp , 1 )
Index: /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
===================================================================
--- /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 3962)
+++ /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 3963)
@@ -142,4 +142,5 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   vt_i_init, vt_i_final  !  ice volume summed over categories
       REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset)
+      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)
       ! mass and salt flux (clem)
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zviold, zvsold, zsmvold   ! old ice volume...
@@ -473,10 +474,14 @@
          zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic )
 
+         zchk_vmin = glob_min(v_i)
+         zchk_amax = glob_max(SUM(a_i,dim=3))
+         zchk_amin = glob_min(a_i)
+       
          IF(lwp) THEN
-            IF (    ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limitd_me) = ',(zchk_v_i * 86400.)
-            IF (    ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_me) = ',(zchk_smv * 86400.)
-            IF ( MINVAL( v_i(:,:,:) ) <  0.    ) WRITE(numout,*) 'violation v_i<0  [mm]         (limitd_me) = ',(MINVAL(v_i) * 1.e-3)
-            IF ( MAXVAL( SUM(a_i(:,:,:),dim=3) ) > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax            (limitd_me) = ',MAXVAL(SUM(a_i,dim=3))
-            IF ( MINVAL( a_i(:,:,:) ) <  0.    ) WRITE(numout,*) 'violation a_i<0               (limitd_me) = ',MINVAL(a_i)
+            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limitd_me) = ',(zchk_v_i * 86400.)
+            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_me) = ',(zchk_smv * 86400.)
+            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limitd_me) = ',(zchk_vmin * 1.e-3)
+            IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limitd_me) = ',zchk_amax
+            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limitd_me) = ',zchk_amin
          ENDIF
       ENDIF
Index: /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
===================================================================
--- /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 3962)
+++ /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 3963)
@@ -68,4 +68,5 @@
       INTEGER ::   jl, ja, jm, jbnd1, jbnd2   ! ice types    dummy loop index         
       REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset)
+      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)
       !!------------------------------------------------------------------
       ! -------------------------------
@@ -165,9 +166,14 @@
          zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic )
 
+         zchk_vmin = glob_min(v_i)
+         zchk_amax = glob_max(SUM(a_i,dim=3))
+         zchk_amin = glob_min(a_i)
+
          IF(lwp) THEN
-            IF (    ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limitd_th) = ',(zchk_v_i * 86400.)
-            IF (    ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_th) = ',(zchk_smv * 86400.)
-            IF ( MINVAL( v_i(:,:,:) ) <  0.    ) WRITE(numout,*) 'violation v_i<0  [mm]         (limitd_th) = ',(MINVAL(v_i) * 1.e-3)
-            IF ( MAXVAL( SUM(a_i(:,:,:),dim=3) ) >  amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax    (limitd_th) = ',MAXVAL(SUM(a_i,dim=3))
+            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limitd_th) = ',(zchk_v_i * 86400.)
+            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_th) = ',(zchk_smv * 86400.)
+            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limitd_th) = ',(zchk_vmin * 1.e-3)
+            IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limitd_th) = ',zchk_amax
+            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limitd_th) = ',zchk_amin
          ENDIF
        ENDIF
@@ -206,5 +212,5 @@
       INTEGER  ::   zji, zjj, nd   ! local integer
       REAL(wp) ::   zx1, zwk1, zdh0, zetamin, zdamax   ! local scalars
-      REAL(wp) ::   zx2, zwk2, zda0, zetamax, zhimin   !   -      -
+      REAL(wp) ::   zx2, zwk2, zda0, zetamax           !   -      -
       REAL(wp) ::   zx3,             zareamin, zindb   !   -      -
       CHARACTER (len = 15) :: fieldid
@@ -242,5 +248,4 @@
       CALL wrk_alloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final )
 
-      zhimin   = 0.1      !minimum ice thickness tolerated by the model
       zareamin = epsi10   !minimum area in thickness categories tolerated by the conceptors of the model
 
@@ -527,10 +532,8 @@
          zji = nind_i(ji)
          zjj = nind_j(ji)
-         IF ( ( zhimin .GT. 0.0 ) .AND. & 
-            ( ( a_i(zji,zjj,1) .GT. epsi10 ) .AND. ( ht_i(zji,zjj,1) .LT. zhimin ) ) &
-            ) THEN
-            a_i(zji,zjj,1)  = a_i(zji,zjj,1) * ht_i(zji,zjj,1) / zhimin 
-            ht_i(zji,zjj,1) = zhimin
-            v_i(zji,zjj,1)  = a_i(zji,zjj,1)*ht_i(zji,zjj,1) !clem@useless
+         IF ( ( a_i(zji,zjj,1) .GT. epsi10 ) .AND. ( ht_i(zji,zjj,1) .LT. hiclim ) ) THEN
+            a_i(zji,zjj,1)  = a_i(zji,zjj,1) * ht_i(zji,zjj,1) / hiclim 
+            ht_i(zji,zjj,1) = hiclim
+            v_i(zji,zjj,1)  = a_i(zji,zjj,1) * ht_i(zji,zjj,1) !clem@useless
          ENDIF
       END DO !ji
Index: /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 3962)
+++ /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 3963)
@@ -94,4 +94,5 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   zqlbsbq   ! link with lead energy budget qldif
       REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset)
+      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)
       !!-------------------------------------------------------------------
 
@@ -476,9 +477,14 @@
          zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic )
 
+         zchk_vmin = glob_min(v_i)
+         zchk_amax = glob_max(SUM(a_i,dim=3))
+         zchk_amin = glob_min(a_i)
+
          IF(lwp) THEN
-            IF (    ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limthd) = ',(zchk_v_i * 86400.)
-            IF (    ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limthd) = ',(zchk_smv * 86400.)
-            IF ( MINVAL( v_i(:,:,:) ) <  0.    ) WRITE(numout,*) 'violation v_i<0  [mm]         (limthd) = ',(MINVAL(v_i) * 1.e-3)
-            IF ( MAXVAL( SUM(a_i(:,:,:),dim=3) ) >  amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax    (limthd) = ',MAXVAL(SUM(a_i,dim=3))
+            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limthd) = ',(zchk_v_i * 86400.)
+            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limthd) = ',(zchk_smv * 86400.)
+            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limthd) = ',(zchk_vmin * 1.e-3)
+            IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limthd) = ',zchk_amax
+            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limthd) = ',zchk_amin
          ENDIF
       ENDIF
Index: /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 3962)
+++ /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 3963)
@@ -121,5 +121,5 @@
 
       ! mass and salt flux (clem)
-      REAL(wp) :: zdvres
+      REAL(wp) :: zdvres, zdvsur, zdvbot
       REAL(wp), POINTER, DIMENSION(:) ::   zviold, zvsold   ! old ice volume...
 
@@ -255,4 +255,6 @@
          zhn            =  1.0 - MAX( zzero , SIGN( zone , - zhsnew ) )
          ht_s_b(ji)     =  MAX( zzero , zhsnew )
+         ! we recompute dh_s_tot (clem) 
+         dh_s_tot (ji)  =  ht_s_b(ji) - zhsold(ji)
          ! Volume and mass variations of snow
          ! dvsbq_1d  (ji) =  a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_mel(ji) )
@@ -283,9 +285,8 @@
             zdq_i    (ji)    = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) / rdt_ice
             !
-            ! IOVINO
-            !zfsalt_melt(ji) = zfsalt_melt(ji) - sm_i_b(ji) * a_i_b(ji)    &
-            !   &                              * MIN( zdeltah(ji,jk) , 0.e0 ) * rhoic / rdt_ice
+            ! clem
             fseqv_1d(ji) = fseqv_1d(ji) - sm_i_b(ji) * a_i_b(ji)    &
                &                              * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic / rdt_ice
+!            fseqv_1d(ji) = fseqv_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic / rdt_ice
          END DO
       END DO
@@ -449,6 +450,6 @@
                ! Salinity update
                ! entrapment during bottom growth
-               dsm_i_se_1d(ji) = ( s_i_new(ji) * dh_i_bott(ji) + sm_i_b(ji) * ht_i_b(ji) )    &
-                  &            / MAX( ht_i_b(ji) + dh_i_bott(ji) ,epsi13 ) - sm_i_b(ji)
+               !clem:movedown dsm_i_se_1d(ji) = ( s_i_new(ji) * dh_i_bott(ji) + sm_i_b(ji) * ht_i_b(ji) )    &
+               !clem:movedown   &            / MAX( ht_i_b(ji) + dh_i_bott(ji) ,epsi13 ) - sm_i_b(ji)
                fseqv_1d(ji) = fseqv_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic / rdt_ice
             ENDIF ! heat budget
@@ -489,9 +490,8 @@
                   zdq_i(ji)       = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) / rdt_ice
                ENDIF
-               ! IOVINO contribution to salt flux
-               !zfsalt_melt(ji) = zfsalt_melt(ji) - sm_i_b(ji) * a_i_b(ji)   &
-               !     &                              * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic / rdt_ice
+               ! clem
                fseqv_1d(ji) = fseqv_1d(ji) - sm_i_b(ji) * a_i_b(ji)    &
                     &                              * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic / rdt_ice
+!               fseqv_1d(ji) = fseqv_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic / rdt_ice
             ENDIF
          END DO ! ji
@@ -549,8 +549,4 @@
          !                     ! excessive energy is sent to lateral ablation
          fsup     (ji) =  rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi13 ) * zdvres / rdt_ice
-         !                     !since ice volume is only used for outputs, we keep it global for all categories
-         dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji)
-         !                     !new ice thickness
-         zhgnew   (ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji)
       END DO
 
@@ -564,9 +560,19 @@
          ! Adapt the remaining energy if too much ice melts
          !--------------------------------------------------
-         zdvres     = MAX( 0._wp, - zhgnew(ji) )
-         zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice
-         dh_i_bott (ji) = dh_i_bott(ji) + zdvres ! clem@bug
-         zhgnew    (ji) =         zihgnew   * zhgnew(ji)      ! ice thickness is put to 0
-         ! remaining heat
+         zdvres     = MAX( 0._wp, - ht_i_b(ji) - dh_i_surf(ji) - dh_i_bott(ji) )
+         zdvsur     = MIN( 0._wp, dh_i_surf(ji) + zdvres ) - dh_i_surf(ji) ! fill the surface first
+         zdvbot     = MAX( 0._wp, zdvres - zdvsur ) ! then the bottom
+         dh_i_surf (ji) = dh_i_surf(ji) + zdvsur ! clem
+         dh_i_bott (ji) = dh_i_bott(ji) + zdvbot ! clem
+
+         ! new ice thickness (clem)
+         zhgnew(ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji)
+         zihgnew    = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice
+         zhgnew(ji) = zihgnew * zhgnew(ji)      ! ice thickness is put to 0
+ 
+         !                     !since ice volume is only used for outputs, we keep it global for all categories
+         dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji)
+
+        ! remaining heat
          zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) +  zqfont_bo(ji) ) 
 
@@ -582,4 +588,6 @@
          ht_s_b(ji)     =  MAX( zzero , zhnfi )
          zqt_s(ji)      =  zqt_s(ji) * ht_s_b(ji)
+         ! we recompute dh_s_tot (clem)
+         dh_s_tot (ji)  =  ht_s_b(ji) - zhsold(ji)
 
          ! Mass variations of ice and snow
@@ -598,6 +606,16 @@
          !---------------------------------
          focea(ji)  = - zfdt_final(ji) / rdt_ice         ! focea is in W.m-2 * dt
-         ! salt flux
-         fseqv_1d(ji)  = fseqv_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic / rdt_ice
+
+         ! residual salt flux (clem)
+         !--------------------------
+         ! surface
+         fseqv_1d(ji)  = fseqv_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvsur * rhoic / rdt_ice
+         ! bottom
+         IF     ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0._wp ) THEN ! melting
+            fseqv_1d(ji)  = fseqv_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvbot * rhoic / rdt_ice
+         ELSE ! growth
+            fseqv_1d(ji)  = fseqv_1d(ji) - s_i_new(ji) * a_i_b(ji) * zdvbot * rhoic / rdt_ice
+         ENDIF
+         !
          !                     ! diagnostic ( bottom ice growth )
          zji = MOD( npb(ji) - 1, jpi ) + 1
@@ -700,6 +718,4 @@
          !ENDIF
          ! entrapment during snow ice formation
-         i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - ht_i_b(ji) + epsi13 ) )
-         isnowic      = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - dh_snowice(ji) ) ) * i_ice_switch
 
          !clem IF(  num_sal == 2  .OR.  num_sal == 4  )   &
@@ -707,7 +723,15 @@
          !clem   &               + sm_i_b(ji) * ht_i_b(ji) / MAX( ht_i_b(ji) + dh_snowice(ji), epsi13)   &
          !clem   &               - sm_i_b(ji) ) * isnowic
-         IF (  num_sal == 2  .OR.  num_sal == 4  )   &
-            & dsm_i_si_1d(ji) = ( ( zsm_snowice * dh_snowice(ji) + sm_i_b(ji) * ht_i_b(ji) ) &
-            &                    / MAX( ht_i_b(ji) + dh_snowice(ji), epsi13 ) - sm_i_b(ji) ) * isnowic     
+
+         ! clem: new salinity difference stored (to be used in limthd_ent.F90)
+         IF (  num_sal == 2  .OR.  num_sal == 4  ) THEN
+            i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - zhgnew(ji) + epsi13 ) )
+            ! salinity dif due to snow-ice formation
+            dsm_i_si_1d(ji) = ( zsm_snowice - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi13 ) * i_ice_switch     
+            ! salinity dif due to bottom growth 
+            IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) < 0._wp ) THEN
+               dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( zhgnew(ji), epsi13 ) * i_ice_switch
+            ENDIF
+         ENDIF
 
          !  Actualize new snow and ice thickness.
Index: /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
===================================================================
--- /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 3962)
+++ /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 3963)
@@ -81,5 +81,5 @@
       INTEGER ::   layer, nbpac     ! local integers 
       INTEGER ::   zji, zjj, iter   !   -       -
-      REAL(wp)  ::   ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zde  ! local scalars
+      REAL(wp)  ::   ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zinda, zde  ! local scalars
       REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new        !   -      -
       REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      -
@@ -454,5 +454,5 @@
             zji                  = MOD( npac(ji) - 1, jpi ) + 1
             zjj                  = ( npac(ji) - 1 ) / jpi + 1
-            diag_lat_gr(zji,zjj) = zv_newice(ji) / rdt_ice
+            diag_lat_gr(zji,zjj) = diag_lat_gr(zji,zjj) + zv_newice(ji) / rdt_ice ! clem
           END DO !ji
 
@@ -547,5 +547,6 @@
             DO ji = 1, nbpac
                zindb = MAX( 0._wp, SIGN( 1._wp , zdv_res(ji) ) )
-               zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi06 )
+               zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_lev(ji) - epsi06 ) )  ! clem
+               zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zinda * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi06 )
             END DO
          END DO
@@ -623,13 +624,13 @@
          ! Update salinity
          !-----------------
-         IF(  num_sal == 2  .OR.  num_sal == 4  ) THEN
+         !clem IF(  num_sal == 2  .OR.  num_sal == 4  ) THEN
             DO jl = 1, jpl
                DO ji = 1, nbpac
                   zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) )  ! 0 if no ice and 1 if yes
                   zdv   = zv_i_ac(ji,jl) - zv_old(ji,jl)
-                  zsmv_i_ac(ji,jl) = ( zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) ) * zindb
+                  zsmv_i_ac(ji,jl) = zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) * zindb ! clem modif
                END DO
             END DO   
-         ENDIF
+         !clem ENDIF
 
          !--------------------------------
@@ -640,5 +641,5 @@
                zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) )  ! 0 if no ice and 1 if yes
                zdv   = zv_i_ac(ji,jl) - zv_old(ji,jl)
-               rdmicif_1d(ji) = rdmicif_1d(ji) + zdv * rhoic !* zindb
+               rdmicif_1d(ji) = rdmicif_1d(ji) + zdv * rhoic * zindb
                fseqv_1d(ji)   =   fseqv_1d(ji) - zdv * rhoic * zs_newice(ji) / rdt_ice * zindb
            END DO
@@ -652,5 +653,5 @@
             CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_ac (1:nbpac,jl), jpi, jpj )
             CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_ac(1:nbpac,jl), jpi, jpj )
-            IF (  num_sal == 2  .OR.  num_sal == 4  )   &
+            !clem IF (  num_sal == 2  .OR.  num_sal == 4  )   &
                CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_ac(1:nbpac,jl) , jpi, jpj )
             DO jk = 1, nlay_i
Index: /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90
===================================================================
--- /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90	(revision 3962)
+++ /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90	(revision 3963)
@@ -127,5 +127,17 @@
             i_ice_switch = 1._wp - MAX ( 0._wp, SIGN( 1._wp , - ht_i_b(ji) ) )
             sm_i_b(ji)   = i_ice_switch * sm_i_b(ji) + s_i_min * ( 1._wp - i_ice_switch )
-         END DO ! ji
+
+            !----------------------------
+            ! Heat flux - brine drainage
+            !----------------------------
+            fhbri_1d(ji) = 0._wp
+
+            !----------------------------
+            ! Salt flux - brine drainage
+            !----------------------------
+            fsbri_1d(ji) = fsbri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji) * ( sm_i_b(ji) - zsiold(ji) ) / rdt_ice
+            IF( num_sal == 4 ) fsbri_1d(ji) = 0._wp
+
+        END DO ! ji
 
          ! Salinity profile
@@ -136,29 +148,29 @@
          !----------------------------
 
-         DO ji = kideb, kiut
+!clem:move         DO ji = kideb, kiut
 !!gm useless
             ! iflush  : 1 if summer 
-            iflush  =  MAX( 0._wp , SIGN ( 1._wp , t_su_b(ji) - rtt ) ) 
+!clem            iflush  =  MAX( 0._wp , SIGN ( 1._wp , t_su_b(ji) - rtt ) ) 
             ! igravdr : 1 if t_su lt t_bo
-            igravdr =  MAX( 0._wp , SIGN ( 1._wp , t_bo_b(ji) - t_su_b(ji) ) ) 
+!clem            igravdr =  MAX( 0._wp , SIGN ( 1._wp , t_bo_b(ji) - t_su_b(ji) ) ) 
             ! iaccrbo : 1 if bottom accretion
-            iaccrbo =  MAX( 0._wp , SIGN ( 1._wp , dh_i_bott(ji) ) )
+!clem            iaccrbo =  MAX( 0._wp , SIGN ( 1._wp , dh_i_bott(ji) ) )
 !!gm end useless
             !
-            fhbri_1d(ji) = 0._wp
-         END DO ! ji
+!clem:move            fhbri_1d(ji) = 0._wp
+!clem:move         END DO ! ji
 
          !----------------------------
          ! Salt flux - brine drainage
          !----------------------------
-          DO ji = kideb, kiut
-            i_ice_switch = 1._wp - MAX ( 0._wp, SIGN( 1._wp , - ht_i_b(ji) ) )
-            fsbri_1d(ji) = fsbri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji) * ( sm_i_b(ji) - zsiold(ji) ) / rdt_ice
+!clem:move          DO ji = kideb, kiut
+!clem:move            i_ice_switch = 1._wp - MAX ( 0._wp, SIGN( 1._wp , - ht_i_b(ji) ) )
+!clem:move            fsbri_1d(ji) = fsbri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji) * ( sm_i_b(ji) - zsiold(ji) ) / rdt_ice
             !i_ice_switch = 1._wp - MAX ( 0._wp, SIGN( 1._wp , - zhiold(ji) ) )
             !fsbri_1d(ji) = fsbri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * zhiold(ji) * ( sm_i_b(ji) - zsiold(ji) ) / rdt_ice
             !clem fsbri_1d(ji) = fsbri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji)         &
             !clem     &         * ( MAX(dsm_i_gd_1d(ji) + dsm_i_fl_1d(ji), sm_i_b(ji) - zsiold(ji) ) ) / rdt_ice
-            IF( num_sal == 4 ) fsbri_1d(ji) = 0._wp
-          END DO ! ji
+!clem:move            IF( num_sal == 4 ) fsbri_1d(ji) = 0._wp
+!clem:move          END DO ! ji
 
          ! Only necessary for conservation check since salinity is modified
Index: /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 3962)
+++ /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 3963)
@@ -81,4 +81,5 @@
       REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::   zs0e
       REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset)
+      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)
       ! mass and salt flux (clem)
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zviold   ! old ice volume...
@@ -170,5 +171,5 @@
          zusnit = 1.0 / REAL( initad ) 
          IF( zcfl > 0.5 .AND. lwp )   &
-            WRITE(numout,*) 'lim_trp_2 : CFL violation at day ', nday, ', cfl = ', zcfl,   &
+            WRITE(numout,*) 'lim_trp   : CFL violation at day ', nday, ', cfl = ', zcfl,   &
                &                        ': the ice time stepping is split in two'
 
@@ -430,5 +431,5 @@
                DO ji = 1, jpi
 
-                  IF ( v_i(ji,jj,jl) > 0. ) THEN
+                  IF ( v_i(ji,jj,jl) > 0._wp ) THEN
                      zvi = v_i(ji,jj,jl)
                      zvs = v_s(ji,jj,jl)
@@ -437,5 +438,5 @@
                      
                      IF ( ( zdv > 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) .AND. SUM( zaiold(ji,jj,1:jpl) ) < 0.80 ) .OR. &
-                          ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN                                          
+                        & ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN                                          
                         ht_i(ji,jj,jl) = MIN( zhimax(ji,jj,jl), hi_max(jl) )
                         zindh   =  MAX( rzero, SIGN( rone, ht_i(ji,jj,jl) - epsi10 ) )
@@ -447,9 +448,9 @@
                      ENDIF
 
-                     !zindh   =  MAX( rzero, SIGN( rone, ht_i(ji,jj,jl) - epsi10 ) )
-                     v_i(ji,jj,jl) = a_i(ji,jj,jl) * ht_i(ji,jj,jl)
-                     v_s(ji,jj,jl) = a_i(ji,jj,jl) * ht_s(ji,jj,jl)
-
-                     ! Update mass fluxes (clem)
+                     ! small correction due to *zindh for a_i
+                     v_i(ji,jj,jl) = zindh * v_i(ji,jj,jl)
+                     v_s(ji,jj,jl) = zindh * v_s(ji,jj,jl)
+
+                     ! Update mass fluxes
                      rdmicif(ji,jj) = rdmicif(ji,jj) + ( v_i(ji,jj,jl) - zvi ) * rhoic
                      rdmsnif(ji,jj) = rdmsnif(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn
@@ -587,9 +588,13 @@
          zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic )
 
+         zchk_vmin = glob_min(v_i)
+         zchk_amax = glob_max(SUM(a_i,dim=3))
+         zchk_amin = glob_min(a_i)
+
          IF(lwp) THEN
-            IF (    ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limtrp) = ',(zchk_v_i * 86400.)
-            IF (    ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limtrp) = ',(zchk_smv * 86400.)
-            IF ( MINVAL( v_i(:,:,:) ) <  0.    ) WRITE(numout,*) 'violation v_i<0  [mm]         (limtrp) = ',(MINVAL(v_i) * 1.e-3)
-            IF ( MINVAL( a_i(:,:,:) ) <  0.    ) WRITE(numout,*) 'violation a_i<0               (limtrp) = ',MINVAL(a_i)
+            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limtrp) = ',(zchk_v_i * 86400.)
+            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limtrp) = ',(zchk_smv * 86400.)
+            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limtrp) = ',(zchk_vmin * 1.e-3)
+            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limtrp) = ',zchk_amin
          ENDIF
       ENDIF
Index: /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
===================================================================
--- /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 3962)
+++ /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 3963)
@@ -86,6 +86,6 @@
       INTEGER ::   i_ice_switch
       INTEGER ::   ind_im, layer      ! indices for internal melt
-      REAL(wp) ::   zweight, zesum, z_da_i
-      REAL(wp) ::   zindb, zindsn, zindic
+      REAL(wp) ::   zweight, zesum, z_da_i, zhimax
+      REAL(wp) ::   zinda, zindb, zindsn, zindic
       REAL(wp) ::   zindg, zh, zdvres, zviold2
       REAL(wp) ::   zbigvalue, zvsold2, z_da_ex
@@ -94,4 +94,5 @@
       REAL(wp), POINTER, DIMENSION(:) ::   zthick0, zqm0      ! thickness of the layers and heat contents for
       REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset)
+      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)
       ! mass and salt flux (clem)
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zviold, zvsold, zsmvold   ! old ice volume...
@@ -182,4 +183,70 @@
       ! 2. Review of all pathological cases
       !--------------------------------------
+
+      !-------------------------------------------
+      ! 2.1) Advection of ice in an ice-free cell
+      !-------------------------------------------
+      ! should be removed since it is treated after dynamics now
+
+!      !IF( ln_nicep ) THEN  
+!         WRITE(numout,*) ' limupdate1 - before h correction '
+!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl)
+!         WRITE(numout,*) ' at_i : ', at_i(12,44)
+!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl)
+!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl)
+!      !ENDIF
+!
+      zhimax = 1._wp
+      ! first category
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            !--- the thickness of such an ice is often out of bounds
+            !--- thus we recompute a new area while conserving ice volume
+            zat_i_old = SUM( old_a_i(ji,jj,:) )
+            zindb     = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_trp(ji,jj,1)) - epsi10 ) ) 
+            IF ( ( ABS(d_v_i_trp(ji,jj,1))/MAX(ABS(d_a_i_trp(ji,jj,1)),epsi10)*zindb .GT. zhimax) &
+                 .AND.( ( v_i(ji,jj,1)/MAX(a_i(ji,jj,1),epsi10)*zindb) .GT. zhimax ) &
+                 .AND.( zat_i_old .LT. 1.e-6 ) )  THEN ! new line
+               ht_i(ji,jj,1) = hi_max(1) / 2.0
+               a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1)
+            ENDIF
+         END DO
+      END DO
+
+
+!      !IF( ln_nicep ) THEN  
+!         at_i(:,:) = 0._wp
+!         DO jl = 1, jpl
+!            at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
+!         END DO
+!         WRITE(numout,*) ' limupdate1 - after h correction 1 '
+!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl)
+!         WRITE(numout,*) ' at_i : ', at_i(12,44)
+!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl)
+!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl)
+!      !ENDIF
+
+      zhimax = 10._wp
+      ! other categories
+      DO jl = 2, jpl
+         jm = ice_types(jl)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zindb =  MAX( rzero, SIGN( rone, ABS(d_a_i_trp(ji,jj,jl)) - epsi10 ) ) 
+               ! this correction is very tricky... sometimes, advection gets wrong i don't know why
+               ! it makes problems when the advected volume and concentration do not seem to be 
+               ! related with each other
+               ! the new thickness is sometimes very big!
+               ! and sometimes d_a_i_trp and d_v_i_trp have different sign
+               ! which of course is plausible
+               ! but fuck! it fucks everything up :)
+               IF ( (ABS(d_v_i_trp(ji,jj,jl))/MAX(ABS(d_a_i_trp(ji,jj,jl)),epsi10)*zindb .GT. zhimax) &
+                    .AND.(v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl),epsi10)*zindb) .GT. zhimax ) THEN
+                  ht_i(ji,jj,jl) = ( hi_max_typ(jl-ice_cat_bounds(jm,1),jm) + hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) / 2.0
+                  a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl)
+               ENDIF
+            END DO ! ji
+         END DO !jj
+      END DO !jl
      
       at_i(:,:) = 0._wp
@@ -187,4 +254,12 @@
          at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
       END DO
+
+!      !IF( ln_nicep ) THEN  
+!         WRITE(numout,*) ' limupdate1 - after h correction 2 '
+!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl)
+!         WRITE(numout,*) ' at_i : ', at_i(12,44)
+!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl)
+!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl)
+!      !ENDIF
 
       !----------------------------------------------------
@@ -247,6 +322,9 @@
                !----------------------------------------------------------------------------
                zindg          = tms(ji,jj) *  MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) ) )
-               v_i(ji,jj,jl)  = v_i(ji,jj,jl) + zindg * rhosn * v_s(ji,jj,jl) / rau0
-               v_s(ji,jj,jl)  = ( 1._wp - zindg ) * v_s(ji,jj,jl) + zindg * v_i(ji,jj,jl) * ( rau0 - rhoic ) / rhosn
+               zdvres         = zindg * rhosn * v_s(ji,jj,jl) / rau0
+               v_i(ji,jj,jl)  = v_i(ji,jj,jl)  + zdvres
+
+               zdvres         = zindsn*zindb * ( - zindg * v_s(ji,jj,jl) + zindg * v_i(ji,jj,jl) * ( rau0 - rhoic ) / rhosn )
+               v_s(ji,jj,jl)  = v_s(ji,jj,jl)  + zdvres
 
                !--- 2.7 Correction to ice concentrations 
@@ -286,19 +364,25 @@
 
       !--- 2.13 ice concentration should not exceed amax 
+      !         (it should not be the case) 
       !-----------------------------------------------------
       DO jj = 1, jpj
          DO ji = 1, jpi
             z_da_ex =  MAX( at_i(ji,jj) - amax , 0.0 )        
+            zindb   =  MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi06 ) ) 
             DO jl  = 1, jpl
-               zindb   =  MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) ) 
                z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi06 ) * zindb
-               a_i(ji,jj,jl) = a_i(ji,jj,jl) - z_da_i
-           END DO
-         END DO 
-      END DO 
+               a_i(ji,jj,jl) = MAX( 0._wp, a_i(ji,jj,jl) - z_da_i )
+               !
+               zinda   =  MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) ) 
+               ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi06 ) * zinda
+               !v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! makes ice shrinken but should not be used
+            END DO
+         END DO
+      END DO
       at_i(:,:) = a_i(:,:,1)
       DO jl = 2, jpl
          at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
       END DO
+
 
       ! Final thickness distribution rebinning
@@ -311,4 +395,5 @@
          ENDIF
       END DO
+
 
       !---------------------
@@ -356,4 +441,5 @@
       !- check conservation (C Rousset)
       IF (ln_limdiahsb) THEN
+
          zchk_fs  = glob_sum( ( fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b
          zchk_fw  = glob_sum( rdmicif(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b
@@ -362,10 +448,14 @@
          zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic )
 
+         zchk_vmin = glob_min(v_i)
+         zchk_amax = glob_max(SUM(a_i,dim=3))
+         zchk_amin = glob_min(a_i)
+
          IF(lwp) THEN
-            IF (    ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limupdate1) = ',(zchk_v_i * 86400.)
-            IF (    ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limupdate1) = ',(zchk_smv * 86400.)
-            IF ( MINVAL( v_i(:,:,:) ) <  0.    ) WRITE(numout,*) 'violation v_i<0  [mm]         (limupdate1) = ',(MINVAL(v_i) * 1.e-3)
-            IF ( MAXVAL( SUM(a_i(:,:,:),dim=3) ) >  amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax    (limupdate1) = ',MAXVAL(SUM(a_i,dim=3))
-            IF ( MINVAL( a_i(:,:,:) ) <  0.    ) WRITE(numout,*) 'violation a_i<0               (limupdate1) = ',MINVAL(a_i)
+            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limupdate1) = ',(zchk_v_i * 86400.)
+            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limupdate1) = ',(zchk_smv * 86400.)
+            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limupdate1) = ',(zchk_vmin * 1.e-3)
+            IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limupdate1) = ',zchk_amax
+            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limupdate1) = ',zchk_amin
          ENDIF
       ENDIF
Index: /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
===================================================================
--- /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 3962)
+++ /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 3963)
@@ -83,5 +83,5 @@
       INTEGER ::   ind_im, layer      ! indices for internal melt
       REAL(wp) ::   zweight, zesum, zhimax, z_da_i
-      REAL(wp) ::   zindb, zindsn, zindic
+      REAL(wp) ::   zinda, zindb, zindsn, zindic
       REAL(wp) ::   zindg, zh, zdvres, zviold2
       REAL(wp) ::   zbigvalue, zvsold2, z_da_ex
@@ -91,5 +91,5 @@
       REAL(wp), POINTER, DIMENSION(:) ::   zthick0, zqm0      ! thickness of the layers and heat contents for
       REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset)
-
+      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)
       ! mass and salt flux (clem)
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zviold, zvsold, zsmvold   ! old ice volume...
@@ -177,9 +177,84 @@
       END DO
 
-
       !--------------------------------------
       ! 2. Review of all pathological cases
       !--------------------------------------
-
+      !-------------------------------------------
+      ! 2.1) Advection of ice in an ice-free cell
+      !-------------------------------------------
+      ! should be removed since it is treated after dynamics now
+
+!      !IF( ln_nicep ) THEN  
+!         WRITE(numout,*) ' limupdate2 - before h correction '
+!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl)
+!         WRITE(numout,*) ' at_i : ', at_i(12,44)
+!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl)
+!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl)
+!      !ENDIF
+
+      zhimax = 1._wp
+      ! first category
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            !--- the thickness of such an ice is often out of bounds
+            !--- thus we recompute a new area while conserving ice volume
+            zat_i_old = SUM( old_a_i(ji,jj,:) )
+            zindb     = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_thd(ji,jj,1)) - epsi10 ) ) 
+            IF ( ( ABS(d_v_i_thd(ji,jj,1))/MAX(ABS(d_a_i_thd(ji,jj,1)),epsi10)*zindb .GT. zhimax) &
+                 .AND.( ( v_i(ji,jj,1)/MAX(a_i(ji,jj,1),epsi10)*zindb) .GT. zhimax ) &
+                 .AND.( zat_i_old .LT. 1.e-6 ) )  THEN ! new line
+               ht_i(ji,jj,1) = hi_max(1) / 2.0
+               a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1)
+            ENDIF
+         END DO
+      END DO
+
+!      !IF( ln_nicep ) THEN  
+!         at_i(:,:) = 0._wp
+!         DO jl = 1, jpl
+!            at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
+!         END DO
+!         WRITE(numout,*) ' limupdate2 - after h correction 1 '
+!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl)
+!         WRITE(numout,*) ' at_i : ', at_i(12,44)
+!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl)
+!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl)
+!      !ENDIF
+! 
+      zhimax = 10._wp
+      ! other categories
+      DO jl = 2, jpl
+         jm = ice_types(jl)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zindb =  MAX( rzero, SIGN( rone, ABS(d_a_i_thd(ji,jj,jl)) - epsi10 ) ) 
+              ! this correction is very tricky... sometimes, advection gets wrong i don't know why
+               ! it makes problems when the advected volume and concentration do not seem to be 
+               ! related with each other
+               ! the new thickness is sometimes very big!
+               ! and sometimes d_a_i_trp and d_v_i_trp have different sign
+               ! which of course is plausible
+               ! but fuck! it fucks everything up :)
+               IF ( (ABS(d_v_i_thd(ji,jj,jl))/MAX(ABS(d_a_i_thd(ji,jj,jl)),epsi10)*zindb .GT. zhimax) &
+                    .AND.(v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl),epsi10)*zindb) .GT. zhimax ) THEN
+                  ht_i(ji,jj,jl) = ( hi_max_typ(jl-ice_cat_bounds(jm,1),jm) + hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) / 2.0
+                  a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl)
+               ENDIF
+            END DO ! ji
+         END DO !jj
+      END DO !jl
+     
+      at_i(:,:) = 0._wp
+      DO jl = 1, jpl
+         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
+      END DO
+
+!      !IF( ln_nicep ) THEN  
+!         WRITE(numout,*) ' limupdate2 - after h correction 2 '
+!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl)
+!         WRITE(numout,*) ' at_i : ', at_i(12,44)
+!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl)
+!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl)
+!      !ENDIF
 
       !----------------------------------------------------
@@ -363,5 +438,5 @@
                v_i(ji,jj,jl)  = v_i(ji,jj,jl)  + zdvres
 
-               zdvres         = - zindg * v_s(ji,jj,jl) + zindg * v_i(ji,jj,jl) * ( rau0 - rhoic ) / rhosn
+               zdvres         = zindsn*zindb * ( - zindg * v_s(ji,jj,jl) + zindg * v_i(ji,jj,jl) * ( rau0 - rhoic ) / rhosn )
                v_s(ji,jj,jl)  = v_s(ji,jj,jl)  + zdvres
 
@@ -421,12 +496,17 @@
       
       !--- 2.13 ice concentration should not exceed amax 
+      !         (it should not be the case) 
       !-----------------------------------------------------
       DO jj = 1, jpj
          DO ji = 1, jpi
             z_da_ex =  MAX( at_i(ji,jj) - amax , 0.0 )        
+            zindb   =  MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi06 ) ) 
             DO jl  = 1, jpl
-               zindb   =  MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) ) 
                z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi06 ) * zindb
-               a_i(ji,jj,jl) = a_i(ji,jj,jl) - z_da_i
+               a_i(ji,jj,jl) = MAX( 0._wp, a_i(ji,jj,jl) - z_da_i )
+               !
+               zinda   =  MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) ) 
+               ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi06 ) * zinda
+               !v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! makes ice shrinken but should not be used
             END DO
          END DO
@@ -517,10 +597,14 @@
          zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic )
 
+         zchk_vmin = glob_min(v_i)
+         zchk_amax = glob_max(SUM(a_i,dim=3))
+         zchk_amin = glob_min(a_i)
+
          IF(lwp) THEN
-            IF (    ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limupdate2) = ',(zchk_v_i * 86400.)
-            IF (    ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limupdate2) = ',(zchk_smv * 86400.)
-            IF ( MINVAL( v_i(:,:,:) ) <  0.    ) WRITE(numout,*) 'violation v_i<0  [mm]         (limupdate2) = ',(MINVAL(v_i) * 1.e-3)
-            IF ( MAXVAL( SUM(a_i(:,:,:),dim=3) ) > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax    (limupdate2) = ',MAXVAL(SUM(a_i,dim=3))
-            IF ( MINVAL( a_i(:,:,:) ) <  0.    ) WRITE(numout,*) 'violation a_i<0               (limupdate2) = ',MINVAL(a_i)
+            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limupdate2) = ',(zchk_v_i * 86400.)
+            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limupdate2) = ',(zchk_smv * 86400.)
+            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limupdate2) = ',(zchk_vmin * 1.e-3)
+            IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limupdate2) = ',zchk_amax
+            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limupdate2) = ',zchk_amin
          ENDIF
       ENDIF
Index: /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
===================================================================
--- /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 3962)
+++ /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 3963)
@@ -180,4 +180,5 @@
                ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb
                o_i(ji,jj,jl)  = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb
+               a_i(ji,jj,jl) = a_i (ji,jj,jl) * zindb ! clem correction
             END DO
          END DO
Index: /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
===================================================================
--- /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 3962)
+++ /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 3963)
@@ -122,33 +122,34 @@
          ! Normal file
          !-------------
-
-         zsto     = rdt_ice
-         IF( ln_mskland )   THEN   ;   clop = "ave(only(x))"   ! put 1.e+20 on land (very expensive!!)
-         ELSE                      ;   clop = "ave(x)"         ! no use of the mask value (require less cpu time)
-         ENDIF
-         zout     = nwrite * rdt_ice / nn_fsbc
          niter    = ( nit000 - 1 ) / nn_fsbc
-         zdept(1) = 0.
-
          CALL ymds2ju ( nyear, nmonth, nday, rdt, zjulian )
          zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
-         CALL dia_nam ( clhstnam, nwrite, 'icemod_old' )
-         CALL histbeg ( clhstnam, jpi, glamt, jpj, gphit, 1, jpi, 1, jpj, niter, zjulian, rdt_ice,   &
-            &           nhorid, nice, domain_id=nidom, snc4chunks=snc4set )
-         CALL histvert( nice, "deptht", "Vertical T levels", "m", 1, zdept, ndepid, "down")
-         CALL wheneq  ( jpij , tmask(:,:,1), 1, 1., ndex51, ndim)
-
-         DO jf = 1 , noumef
-            IF(lwp) WRITE(numout,*) 'jf', jf
-            IF ( nc(jf) == 1 ) THEN
-               CALL histdef( nice, nam(jf), titn(jf), uni(jf), jpi, jpj &
-                  , nhorid, 1, 1, 1, -99, 32, clop, zsto, zout )
-               IF(lwp) WRITE(numout,*) 'nice, nam(jf), titn(jf), uni(jf), nhorid, clop, zsto, zout'
-               IF(lwp) WRITE(numout,*)  nice, nam(jf), titn(jf), uni(jf), nhorid, clop, zsto, zout 
-            ENDIF
-         END DO
-
-         CALL histend(nice, snc4set)
-
+!clem
+!         zsto     = rdt_ice
+!         IF( ln_mskland )   THEN   ;   clop = "ave(only(x))"   ! put 1.e+20 on land (very expensive!!)
+!         ELSE                      ;   clop = "ave(x)"         ! no use of the mask value (require less cpu time)
+!         ENDIF
+!         zout     = nwrite * rdt_ice / nn_fsbc
+!         zdept(1) = 0.
+!
+!         CALL dia_nam ( clhstnam, nwrite, 'icemod_old' )
+!         CALL histbeg ( clhstnam, jpi, glamt, jpj, gphit, 1, jpi, 1, jpj, niter, zjulian, rdt_ice,   &
+!            &           nhorid, nice, domain_id=nidom, snc4chunks=snc4set )
+!         CALL histvert( nice, "deptht", "Vertical T levels", "m", 1, zdept, ndepid, "down")
+!         CALL wheneq  ( jpij , tmask(:,:,1), 1, 1., ndex51, ndim)
+!
+!         DO jf = 1 , noumef
+!            IF(lwp) WRITE(numout,*) 'jf', jf
+!            IF ( nc(jf) == 1 ) THEN
+!               CALL histdef( nice, nam(jf), titn(jf), uni(jf), jpi, jpj &
+!                  , nhorid, 1, 1, 1, -99, 32, clop, zsto, zout )
+!               IF(lwp) WRITE(numout,*) 'nice, nam(jf), titn(jf), uni(jf), nhorid, clop, zsto, zout'
+!               IF(lwp) WRITE(numout,*)  nice, nam(jf), titn(jf), uni(jf), nhorid, clop, zsto, zout 
+!            ENDIF
+!         END DO
+!
+!         CALL histend(nice, snc4set)
+!clem
+         !
          !-----------------
          ! ITD file output
@@ -280,27 +281,30 @@
       ! ecriture d'un fichier netcdf
       !
+ 
       niter = niter + 1
-      DO jf = 1 , noumef
-         !
-         zfield(:,:) = zcmo(:,:,jf) * cmulti(jf) + cadd(jf)
-         !
-         IF( jf == 7  .OR. jf == 8  .OR. jf == 15 .OR. jf == 16 ) THEN   ;   CALL lbc_lnk( zfield, 'T', -1. )
-         ELSE                                                            ;   CALL lbc_lnk( zfield, 'T',  1. )
-         ENDIF
-         !
-         IF( ln_nicep ) THEN 
-            WRITE(numout,*)
-            WRITE(numout,*) 'nc(jf), nice, nam(jf), niter, ndim'
-            WRITE(numout,*) nc(jf), nice, nam(jf), niter, ndim
-         ENDIF
-         IF( nc(jf) == 1 ) CALL histwrite( nice, nam(jf), niter, zfield, ndim, ndex51 )
-         !
-      END DO
-
-      IF( ( nn_fsbc * niter ) >= nitend .OR. kindic < 0 ) THEN
-         IF( lwp) WRITE(numout,*) ' Closing the icemod file '
-         CALL histclo( nice )
-      ENDIF
-
+!clem
+!      DO jf = 1 , noumef
+!         !
+!         zfield(:,:) = zcmo(:,:,jf) * cmulti(jf) + cadd(jf)
+!         !
+!         IF( jf == 7  .OR. jf == 8  .OR. jf == 15 .OR. jf == 16 ) THEN   ;   CALL lbc_lnk( zfield, 'T', -1. )
+!         ELSE                                                            ;   CALL lbc_lnk( zfield, 'T',  1. )
+!         ENDIF
+!         !
+!         IF( ln_nicep ) THEN 
+!            WRITE(numout,*)
+!            WRITE(numout,*) 'nc(jf), nice, nam(jf), niter, ndim'
+!            WRITE(numout,*) nc(jf), nice, nam(jf), niter, ndim
+!         ENDIF
+!         IF( nc(jf) == 1 ) CALL histwrite( nice, nam(jf), niter, zfield, ndim, ndex51 )
+!         !
+!      END DO
+!
+!      IF( ( nn_fsbc * niter ) >= nitend .OR. kindic < 0 ) THEN
+!         IF( lwp) WRITE(numout,*) ' Closing the icemod file '
+!         CALL histclo( nice )
+!      ENDIF
+!clem
+      !
        CALL iom_put ('iceconc', zcmo(:,:,1) )          ! field1: ice concentration
        CALL iom_put ('icethic_cea', zcmo(:,:,2) )      ! field2: ice thickness (i.e. icethi(:,:))
@@ -328,6 +332,6 @@
        CALL iom_put ('isnowhc', zcmo(:,:,26) )         ! field 26: snow total heat content
        CALL iom_put ('icest', zcmo(:,:,27) )           ! field 27: ice surface temperature
-       CALL iom_put ('fsbri', zcmo(:,:,28) * 86400 )           ! field 28: brine salt flux
-       CALL iom_put ('fseqv', zcmo(:,:,29) * 86400 )           ! field 29: equivalent FW salt flux
+       CALL iom_put ('fsbri', zcmo(:,:,28) * 86400. )           ! field 28: brine salt flux
+       CALL iom_put ('fseqv', zcmo(:,:,29) * 86400. )           ! field 29: equivalent FW salt flux
        CALL iom_put ('ibrinv', zcmo(:,:,30) *100 )     ! field 30: brine volume
        CALL iom_put ('icecolf', zcmo(:,:,31) )         ! field 31: frazil ice collection thickness
@@ -341,6 +345,6 @@
        CALL iom_put ('icevolu', zcmo(:,:,39) ) ! field 39: ice volume
        CALL iom_put ('snowvol', zcmo(:,:,40) ) ! field 40: snow volume
-       CALL iom_put ('fsrpo', zcmo(:,:,41) * 86400 )           ! field 41: salt flux from ridging rafting
-       CALL iom_put ('fsres', zcmo(:,:,42) * 86400 )           ! field 42: salt flux from limupdate (resultant)
+       CALL iom_put ('fsrpo', zcmo(:,:,41) * 86400. )           ! field 41: salt flux from ridging rafting
+       CALL iom_put ('fsres', zcmo(:,:,42) * 86400. )           ! field 42: salt flux from limupdate (resultant)
        CALL iom_put ('icetrp', zcmo(:,:,43) )    ! field 43: ice volume transport
 
Index: /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri_dimg.h90
===================================================================
--- /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri_dimg.h90	(revision 3962)
+++ /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri_dimg.h90	(revision 3963)
@@ -14,4 +14,5 @@
    !!  modif : 03/06/98
    !!-------------------------------------------------------------------
+   USE lib_fortran     ! to use key_nosignedzero
    USE  diawri, ONLY : dia_wri_dimg
    REAL(wp),DIMENSION(1) ::   zdept
Index: /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 3962)
+++ /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 3963)
@@ -41,4 +41,5 @@
    USE sbc_ice         ! Surface boundary condition: ice fields
 #endif
+   USE lib_fortran     ! to use key_nosignedzero
 
    IMPLICIT NONE
@@ -530,7 +531,7 @@
             DO ji = fs_2, fs_jpim1   ! vect. opt.
                p_taui(ji,jj) = zcoef_wnorm2 * ( z_wnds_t(ji+1,jj  ) + z_wnds_t(ji,jj) )                          &
-                  &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - pui(ji,jj) )
+                  &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * pui(ji,jj) )
                p_tauj(ji,jj) = zcoef_wnorm2 * ( z_wnds_t(ji,jj+1  ) + z_wnds_t(ji,jj) )                          &
-                  &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - pvi(ji,jj) )
+                  &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * pvi(ji,jj) )
             END DO
          END DO
@@ -560,7 +561,7 @@
                ! iovino
                IF( ff(ji,jj) .GT. 0._wp ) THEN
-                  z_qlw(ji,jj,jl) = ( 0.95 * sf(jp_qlw)%fnow(ji,jj,1) - Stef * pst(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) ! MV test 2 garde juste cette ligne ci
+                  z_qlw(ji,jj,jl) = ( 0.95 * sf(jp_qlw)%fnow(ji,jj,1) - Stef * pst(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
                ELSE
-                  z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * pst(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) ! MV test 1 garde juste cette ligne ci
+                  z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * pst(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
                ENDIF
                ! lw sensitivity
@@ -610,6 +611,6 @@
 !CDIR COLLAPSE
       p_spr(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac      ! solid precipitation [kg/m2/s]
-      CALL iom_put( 'snowpre', p_spr * 86400 )                  ! Snow precipitation 
-      CALL iom_put( 'precip', p_tpr * 86400 )                   ! Total precipitation 
+      CALL iom_put( 'snowpre', p_spr * 86400. )                  ! Snow precipitation 
+      CALL iom_put( 'precip', p_tpr * 86400. )                   ! Total precipitation 
       !
       IF(ln_ctl) THEN
Index: /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 3962)
+++ /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 3963)
@@ -97,9 +97,9 @@
 
       !                          ! overwrite namelist parameter using CPP key information
-      !IF( Agrif_Root() ) THEN                ! AGRIF zoom
-      !  IF( lk_lim2 )   nn_ice      = 2
-      !  IF( lk_lim3 )   nn_ice      = 3
-      !  IF( lk_cice )   nn_ice      = 4
-      !ENDIF
+      IF( Agrif_Root() ) THEN                ! AGRIF zoom
+        IF( lk_lim2 )   nn_ice      = 2
+        IF( lk_lim3 )   nn_ice      = 3
+        IF( lk_cice )   nn_ice      = 4
+      ENDIF
       IF( cp_cfg == 'gyre' ) THEN            ! GYRE configuration
           ln_ana      = .TRUE.   
Index: /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/lib_fortran.F90
===================================================================
--- /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/lib_fortran.F90	(revision 3962)
+++ /branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/lib_fortran.F90	(revision 3963)
@@ -5,4 +5,6 @@
    !!======================================================================
    !! History :  3.2  !  2010-05  (M. Dunphy, R. Benshila)  Original code
+   !!            3.4  !  2013-06  (C. Rousset)  add glob_min, glob_max 
+   !!                                           + 3d dim. of input is fexible (jpk, jpl...) 
    !!----------------------------------------------------------------------
 
@@ -23,4 +25,5 @@
 
    PUBLIC glob_sum
+   PUBLIC glob_min, glob_max
 #if defined key_nosignedzero
    PUBLIC SIGN
@@ -29,4 +32,10 @@
    INTERFACE glob_sum
       MODULE PROCEDURE glob_sum_2d, glob_sum_3d,glob_sum_2d_a, glob_sum_3d_a 
+   END INTERFACE
+   INTERFACE glob_min
+      MODULE PROCEDURE glob_min_2d, glob_min_3d,glob_min_2d_a, glob_min_3d_a 
+   END INTERFACE
+   INTERFACE glob_max
+      MODULE PROCEDURE glob_max_2d, glob_max_3d,glob_max_2d_a, glob_max_3d_a 
    END INTERFACE
 
@@ -47,4 +56,6 @@
 
 #if ! defined key_mpp_rep
+
+   ! --- SUM ---
    FUNCTION glob_sum_2d( ptab ) 
       !!-----------------------------------------------------------------------
@@ -61,6 +72,5 @@
       !
    END FUNCTION glob_sum_2d
-   
-   
+  
    FUNCTION glob_sum_3d( ptab ) 
       !!-----------------------------------------------------------------------
@@ -73,8 +83,11 @@
       !!
       INTEGER :: jk
-      !!-----------------------------------------------------------------------
+      INTEGER :: zjpk ! local variable: size of the 3d dimension of ptab
+      !!-----------------------------------------------------------------------
+      !
+      zjpk = SIZE(ptab,3)
       !
       glob_sum_3d = 0.e0
-      DO jk = 1, jpk
+      DO jk = 1, zjpk
          glob_sum_3d = glob_sum_3d + SUM( ptab(:,:,jk)*tmask_i(:,:) )
       END DO
@@ -111,8 +124,11 @@
       !!
       INTEGER :: jk
-      !!-----------------------------------------------------------------------
+      INTEGER :: zjpk ! local variable: size of the 3d dimension of ptab
+      !!-----------------------------------------------------------------------
+      !
+      zjpk = SIZE(ptab1,3)
       !
       glob_sum_3d_a(:) = 0.e0
-      DO jk = 1, jpk
+      DO jk = 1, zjpk
          glob_sum_3d_a(1) = glob_sum_3d_a(1) + SUM( ptab1(:,:,jk)*tmask_i(:,:) )
          glob_sum_3d_a(2) = glob_sum_3d_a(2) + SUM( ptab2(:,:,jk)*tmask_i(:,:) )
@@ -121,4 +137,168 @@
       !
    END FUNCTION glob_sum_3d_a
+ 
+
+   ! --- MIN ---
+   FUNCTION glob_min_2d( ptab ) 
+      !!-----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_min_2D  ***
+      !!
+      !! ** Purpose : perform a masked min on the inner global domain of a 2D array
+      !!-----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab          ! input 2D array
+      REAL(wp)                             ::   glob_min_2d   ! global masked min
+      !!-----------------------------------------------------------------------
+      !
+      glob_min_2d = MINVAL( ptab(:,:)*tmask_i(:,:) )
+      IF( lk_mpp )   CALL mpp_min( glob_min_2d )
+      !
+   END FUNCTION glob_min_2d
+ 
+   FUNCTION glob_min_3d( ptab ) 
+      !!-----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_min_3D  ***
+      !!
+      !! ** Purpose : perform a masked min on the inner global domain of a 3D array
+      !!-----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab          ! input 3D array
+      REAL(wp)                               ::   glob_min_3d   ! global masked min
+      !!
+      INTEGER :: jk
+      INTEGER :: zjpk ! local variable: size of the 3d dimension of ptab
+      !!-----------------------------------------------------------------------
+      !
+      zjpk = SIZE(ptab,3)
+      !
+      glob_min_3d = 0.e0
+      DO jk = 1, zjpk
+         glob_min_3d = glob_min_3d + MINVAL( ptab(:,:,jk)*tmask_i(:,:) )
+      END DO
+      IF( lk_mpp )   CALL mpp_min( glob_min_3d )
+      !
+   END FUNCTION glob_min_3d
+
+
+   FUNCTION glob_min_2d_a( ptab1, ptab2 ) 
+      !!-----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_min_2D _a ***
+      !!
+      !! ** Purpose : perform a masked min on the inner global domain of two 2D array
+      !!-----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab1, ptab2    ! input 2D array
+      REAL(wp)            , DIMENSION(2)   ::   glob_min_2d_a   ! global masked min
+      !!-----------------------------------------------------------------------
+      !             
+      glob_min_2d_a(1) = MINVAL( ptab1(:,:)*tmask_i(:,:) )
+      glob_min_2d_a(2) = MINVAL( ptab2(:,:)*tmask_i(:,:) )
+      IF( lk_mpp )   CALL mpp_min( glob_min_2d_a, 2 )
+      !
+   END FUNCTION glob_min_2d_a
+ 
+ 
+   FUNCTION glob_min_3d_a( ptab1, ptab2 ) 
+      !!-----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_min_3D_a ***
+      !!
+      !! ** Purpose : perform a masked min on the inner global domain of two 3D array
+      !!-----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab1, ptab2    ! input 3D array
+      REAL(wp)            , DIMENSION(2)     ::   glob_min_3d_a   ! global masked min
+      !!
+      INTEGER :: jk
+      INTEGER :: zjpk ! local variable: size of the 3d dimension of ptab
+      !!-----------------------------------------------------------------------
+      !
+      zjpk = SIZE(ptab1,3)
+      !
+      glob_min_3d_a(:) = 0.e0
+      DO jk = 1, zjpk
+         glob_min_3d_a(1) = glob_min_3d_a(1) + MINVAL( ptab1(:,:,jk)*tmask_i(:,:) )
+         glob_min_3d_a(2) = glob_min_3d_a(2) + MINVAL( ptab2(:,:,jk)*tmask_i(:,:) )
+      END DO
+      IF( lk_mpp )   CALL mpp_min( glob_min_3d_a, 2 )
+      !
+   END FUNCTION glob_min_3d_a
+
+   ! --- MAX ---
+   FUNCTION glob_max_2d( ptab ) 
+      !!-----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_max_2D  ***
+      !!
+      !! ** Purpose : perform a masked max on the inner global domain of a 2D array
+      !!-----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab          ! input 2D array
+      REAL(wp)                             ::   glob_max_2d   ! global masked max
+      !!-----------------------------------------------------------------------
+      !
+      glob_max_2d = MAXVAL( ptab(:,:)*tmask_i(:,:) )
+      IF( lk_mpp )   CALL mpp_max( glob_max_2d )
+      !
+   END FUNCTION glob_max_2d
+ 
+   FUNCTION glob_max_3d( ptab ) 
+      !!-----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_max_3D  ***
+      !!
+      !! ** Purpose : perform a masked max on the inner global domain of a 3D array
+      !!-----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab          ! input 3D array
+      REAL(wp)                               ::   glob_max_3d   ! global masked max
+      !!
+      INTEGER :: jk
+      INTEGER :: zjpk ! local variable: size of the 3d dimension of ptab
+      !!-----------------------------------------------------------------------
+      !
+      zjpk = SIZE(ptab,3)
+      !
+      glob_max_3d = 0.e0
+      DO jk = 1, zjpk
+         glob_max_3d = glob_max_3d + MAXVAL( ptab(:,:,jk)*tmask_i(:,:) )
+      END DO
+      IF( lk_mpp )   CALL mpp_max( glob_max_3d )
+      !
+   END FUNCTION glob_max_3d
+
+
+   FUNCTION glob_max_2d_a( ptab1, ptab2 ) 
+      !!-----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_max_2D _a ***
+      !!
+      !! ** Purpose : perform a masked max on the inner global domain of two 2D array
+      !!-----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab1, ptab2    ! input 2D array
+      REAL(wp)            , DIMENSION(2)   ::   glob_max_2d_a   ! global masked max
+      !!-----------------------------------------------------------------------
+      !             
+      glob_max_2d_a(1) = MAXVAL( ptab1(:,:)*tmask_i(:,:) )
+      glob_max_2d_a(2) = MAXVAL( ptab2(:,:)*tmask_i(:,:) )
+      IF( lk_mpp )   CALL mpp_max( glob_max_2d_a, 2 )
+      !
+   END FUNCTION glob_max_2d_a
+ 
+ 
+   FUNCTION glob_max_3d_a( ptab1, ptab2 ) 
+      !!-----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_max_3D_a ***
+      !!
+      !! ** Purpose : perform a masked max on the inner global domain of two 3D array
+      !!-----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab1, ptab2    ! input 3D array
+      REAL(wp)            , DIMENSION(2)     ::   glob_max_3d_a   ! global masked max
+      !!
+      INTEGER :: jk
+      INTEGER :: zjpk ! local variable: size of the 3d dimension of ptab
+      !!-----------------------------------------------------------------------
+      !
+      zjpk = SIZE(ptab1,3)
+      !
+      glob_max_3d_a(:) = 0.e0
+      DO jk = 1, zjpk
+         glob_max_3d_a(1) = glob_max_3d_a(1) + MAXVAL( ptab1(:,:,jk)*tmask_i(:,:) )
+         glob_max_3d_a(2) = glob_max_3d_a(2) + MAXVAL( ptab2(:,:,jk)*tmask_i(:,:) )
+      END DO
+      IF( lk_mpp )   CALL mpp_max( glob_max_3d_a, 2 )
+      !
+   END FUNCTION glob_max_3d_a
+
 
 #else  
@@ -127,4 +307,5 @@
    !!----------------------------------------------------------------------
    
+   ! --- SUM ---
    FUNCTION glob_sum_2d( ptab ) 
       !!----------------------------------------------------------------------
@@ -133,16 +314,19 @@
       !! ** Purpose : perform a sum in calling DDPDD routine
       !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in), DIMENSION(jpi,jpj) ::   ptab
-      REAL(wp)                                 ::   glob_sum_2d   ! global masked sum
+      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab
+      REAL(wp)                             ::   glob_sum_2d   ! global masked sum
       !!
       COMPLEX(wp)::   ctmp
       REAL(wp)   ::   ztmp
       INTEGER    ::   ji, jj   ! dummy loop indices
-      !!-----------------------------------------------------------------------
-      !
-      ztmp = 0.e0
-      ctmp = CMPLX( 0.e0, 0.e0, wp )
-      DO jj = 1, jpj
-         DO ji =1, jpi
+      INTEGER    ::   zjpi, zjpj ! local variables: size of ptab
+      !!-----------------------------------------------------------------------
+      zjpi = SIZE(ptab,1)
+      zjpj = SIZE(ptab,2)
+      !
+      ztmp = 0.e0
+      ctmp = CMPLX( 0.e0, 0.e0, wp )
+      DO jj = 1, zjpj
+         DO ji =1, zjpi
          ztmp =  ptab(ji,jj) * tmask_i(ji,jj)
          CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
@@ -161,17 +345,22 @@
       !! ** Purpose : perform a sum on a 3D array in calling DDPDD routine
       !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   ptab
-      REAL(wp)                                     ::   glob_sum_3d   ! global masked sum
+      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab
+      REAL(wp)                               ::   glob_sum_3d   ! global masked sum
       !!
       COMPLEX(wp)::   ctmp
       REAL(wp)   ::   ztmp
       INTEGER    ::   ji, jj, jk   ! dummy loop indices
-      !!-----------------------------------------------------------------------
-      !
-      ztmp = 0.e0
-      ctmp = CMPLX( 0.e0, 0.e0, wp )
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji =1, jpi
+      INTEGER    ::   zjpi, zjpj, zjpk ! local variables: size of ptab
+      !!-----------------------------------------------------------------------
+      !
+      zjpi = SIZE(ptab,1)
+      zjpj = SIZE(ptab,2)
+      zjpk = SIZE(ptab,3)
+      !
+      ztmp = 0.e0
+      ctmp = CMPLX( 0.e0, 0.e0, wp )
+      DO jk = 1, zjpk
+         DO jj = 1, zjpj
+            DO ji =1, zjpi
             ztmp =  ptab(ji,jj,jk) * tmask_i(ji,jj)
             CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
@@ -191,16 +380,20 @@
       !! ** Purpose : perform a sum on two 2D arrays in calling DDPDD routine
       !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in), DIMENSION(jpi,jpj) ::   ptab1, ptab2
-      REAL(wp)                                 ::   glob_sum_2d_a   ! global masked sum
+      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab1, ptab2
+      REAL(wp)                             ::   glob_sum_2d_a   ! global masked sum
       !!
       COMPLEX(wp)::   ctmp
       REAL(wp)   ::   ztmp
       INTEGER    ::   ji, jj   ! dummy loop indices
-      !!-----------------------------------------------------------------------
-      !
-      ztmp = 0.e0
-      ctmp = CMPLX( 0.e0, 0.e0, wp )
-      DO jj = 1, jpj
-         DO ji =1, jpi
+      INTEGER    ::   zjpi, zjpj ! local variables: size of ptab
+      !!-----------------------------------------------------------------------
+      !
+      zjpi = SIZE(ptab1,1)
+      zjpj = SIZE(ptab1,2)
+      !
+      ztmp = 0.e0
+      ctmp = CMPLX( 0.e0, 0.e0, wp )
+      DO jj = 1, zjpj
+         DO ji =1, zjpi
          ztmp =  ptab1(ji,jj) * tmask_i(ji,jj)
          CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
@@ -221,17 +414,22 @@
       !! ** Purpose : perform a sum on two 3D array in calling DDPDD routine
       !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   ptab1, ptab2
-      REAL(wp)                                     ::   glob_sum_3d_a   ! global masked sum
+      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab1, ptab2
+      REAL(wp)                               ::   glob_sum_3d_a   ! global masked sum
       !!
       COMPLEX(wp)::   ctmp
       REAL(wp)   ::   ztmp
       INTEGER    ::   ji, jj, jk   ! dummy loop indices
-      !!-----------------------------------------------------------------------
-      !
-      ztmp = 0.e0
-      ctmp = CMPLX( 0.e0, 0.e0, wp )
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji =1, jpi
+      INTEGER    ::   zjpi, zjpj, zjpk ! local variables: size of ptab
+      !!-----------------------------------------------------------------------
+      !
+      zjpi = SIZE(ptab1,1)
+      zjpj = SIZE(ptab1,2)
+      zjpk = SIZE(ptab1,3)
+      !
+      ztmp = 0.e0
+      ctmp = CMPLX( 0.e0, 0.e0, wp )
+      DO jk = 1, zjpk
+         DO jj = 1, zjpj
+            DO ji =1, zjpi
             ztmp =  ptab1(ji,jj,jk) * tmask_i(ji,jj)
             CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
@@ -247,4 +445,280 @@
 
 
+   ! --- MIN ---
+   FUNCTION glob_min_2d( ptab ) 
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_min_2d ***
+      !!
+      !! ** Purpose : perform a min in calling DDPDD routine
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab
+      REAL(wp)                             ::   glob_min_2d   ! global masked min
+      !!
+      COMPLEX(wp)::   ctmp
+      REAL(wp)   ::   ztmp
+      INTEGER    ::   ji, jj   ! dummy loop indices
+      INTEGER    ::   zjpi, zjpj ! local variables: size of ptab
+      !!-----------------------------------------------------------------------
+      zjpi = SIZE(ptab,1)
+      zjpj = SIZE(ptab,2)
+      !
+      ztmp = 0.e0
+      ctmp = CMPLX( 0.e0, 0.e0, wp )
+      DO jj = 1, zjpj
+         DO ji =1, zjpi
+         ztmp =  ptab(ji,jj) * tmask_i(ji,jj)
+         CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
+         END DO
+      END DO
+      IF( lk_mpp )   CALL mpp_min( ctmp )   ! min over the global domain
+      glob_min_2d = REAL(ctmp,wp)
+      !
+   END FUNCTION glob_min_2d   
+
+
+   FUNCTION glob_min_3d( ptab ) 
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_min_3d ***
+      !!
+      !! ** Purpose : perform a min on a 3D array in calling DDPDD routine
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab
+      REAL(wp)                               ::   glob_min_3d   ! global masked min
+      !!
+      COMPLEX(wp)::   ctmp
+      REAL(wp)   ::   ztmp
+      INTEGER    ::   ji, jj, jk   ! dummy loop indices
+      INTEGER    ::   zjpi, zjpj, zjpk ! local variables: size of ptab
+      !!-----------------------------------------------------------------------
+      !
+      zjpi = SIZE(ptab,1)
+      zjpj = SIZE(ptab,2)
+      zjpk = SIZE(ptab,3)
+      !
+      ztmp = 0.e0
+      ctmp = CMPLX( 0.e0, 0.e0, wp )
+      DO jk = 1, zjpk
+         DO jj = 1, zjpj
+            DO ji =1, zjpi
+            ztmp =  ptab(ji,jj,jk) * tmask_i(ji,jj)
+            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
+            END DO
+         END DO    
+      END DO
+      IF( lk_mpp )   CALL mpp_min( ctmp )   ! min over the global domain
+      glob_min_3d = REAL(ctmp,wp)
+      !
+   END FUNCTION glob_min_3d   
+
+
+   FUNCTION glob_min_2d_a( ptab1, ptab2 ) 
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_min_2d_a ***
+      !!
+      !! ** Purpose : perform a min on two 2D arrays in calling DDPDD routine
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab1, ptab2
+      REAL(wp)                             ::   glob_min_2d_a   ! global masked min
+      !!
+      COMPLEX(wp)::   ctmp
+      REAL(wp)   ::   ztmp
+      INTEGER    ::   ji, jj   ! dummy loop indices
+      INTEGER    ::   zjpi, zjpj ! local variables: size of ptab
+      !!-----------------------------------------------------------------------
+      !
+      zjpi = SIZE(ptab1,1)
+      zjpj = SIZE(ptab1,2)
+      !
+      ztmp = 0.e0
+      ctmp = CMPLX( 0.e0, 0.e0, wp )
+      DO jj = 1, zjpj
+         DO ji =1, zjpi
+         ztmp =  ptab1(ji,jj) * tmask_i(ji,jj)
+         CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
+         ztmp =  ptab2(ji,jj) * tmask_i(ji,jj)
+         CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
+         END DO
+      END DO
+      IF( lk_mpp )   CALL mpp_min( ctmp )   ! min over the global domain
+      glob_min_2d_a = REAL(ctmp,wp)
+      !
+   END FUNCTION glob_min_2d_a   
+
+
+   FUNCTION glob_min_3d_a( ptab1, ptab2 ) 
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_min_3d_a ***
+      !!
+      !! ** Purpose : perform a min on two 3D array in calling DDPDD routine
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab1, ptab2
+      REAL(wp)                               ::   glob_min_3d_a   ! global masked min
+      !!
+      COMPLEX(wp)::   ctmp
+      REAL(wp)   ::   ztmp
+      INTEGER    ::   ji, jj, jk   ! dummy loop indices
+      INTEGER    ::   zjpi, zjpj, zjpk ! local variables: size of ptab
+      !!-----------------------------------------------------------------------
+      !
+      zjpi = SIZE(ptab1,1)
+      zjpj = SIZE(ptab1,2)
+      zjpk = SIZE(ptab1,3)
+      !
+      ztmp = 0.e0
+      ctmp = CMPLX( 0.e0, 0.e0, wp )
+      DO jk = 1, zjpk
+         DO jj = 1, zjpj
+            DO ji =1, zjpi
+            ztmp =  ptab1(ji,jj,jk) * tmask_i(ji,jj)
+            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
+            ztmp =  ptab2(ji,jj,jk) * tmask_i(ji,jj)
+            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
+            END DO
+         END DO    
+      END DO
+      IF( lk_mpp )   CALL mpp_min( ctmp )   ! min over the global domain
+      glob_min_3d_a = REAL(ctmp,wp)
+      !
+   END FUNCTION glob_min_3d_a   
+
+ 
+   ! --- MAX ---
+   FUNCTION glob_max_2d( ptab ) 
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_max_2d ***
+      !!
+      !! ** Purpose : perform a max in calling DDPDD routine
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab
+      REAL(wp)                             ::   glob_max_2d   ! global masked max
+      !!
+      COMPLEX(wp)::   ctmp
+      REAL(wp)   ::   ztmp
+      INTEGER    ::   ji, jj   ! dummy loop indices
+      INTEGER    ::   zjpi, zjpj ! local variables: size of ptab
+      !!-----------------------------------------------------------------------
+      zjpi = SIZE(ptab,1)
+      zjpj = SIZE(ptab,2)
+      !
+      ztmp = 0.e0
+      ctmp = CMPLX( 0.e0, 0.e0, wp )
+      DO jj = 1, zjpj
+         DO ji =1, zjpi
+         ztmp =  ptab(ji,jj) * tmask_i(ji,jj)
+         CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
+         END DO
+      END DO
+      IF( lk_mpp )   CALL mpp_max( ctmp )   ! max over the global domain
+      glob_max_2d = REAL(ctmp,wp)
+      !
+   END FUNCTION glob_max_2d   
+
+
+   FUNCTION glob_max_3d( ptab ) 
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_max_3d ***
+      !!
+      !! ** Purpose : perform a max on a 3D array in calling DDPDD routine
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab
+      REAL(wp)                               ::   glob_max_3d   ! global masked max
+      !!
+      COMPLEX(wp)::   ctmp
+      REAL(wp)   ::   ztmp
+      INTEGER    ::   ji, jj, jk   ! dummy loop indices
+      INTEGER    ::   zjpi, zjpj, zjpk ! local variables: size of ptab
+      !!-----------------------------------------------------------------------
+      !
+      zjpi = SIZE(ptab,1)
+      zjpj = SIZE(ptab,2)
+      zjpk = SIZE(ptab,3)
+      !
+      ztmp = 0.e0
+      ctmp = CMPLX( 0.e0, 0.e0, wp )
+      DO jk = 1, zjpk
+         DO jj = 1, zjpj
+            DO ji =1, zjpi
+            ztmp =  ptab(ji,jj,jk) * tmask_i(ji,jj)
+            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
+            END DO
+         END DO    
+      END DO
+      IF( lk_mpp )   CALL mpp_max( ctmp )   ! max over the global domain
+      glob_max_3d = REAL(ctmp,wp)
+      !
+   END FUNCTION glob_max_3d   
+
+
+   FUNCTION glob_max_2d_a( ptab1, ptab2 ) 
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_max_2d_a ***
+      !!
+      !! ** Purpose : perform a max on two 2D arrays in calling DDPDD routine
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab1, ptab2
+      REAL(wp)                             ::   glob_max_2d_a   ! global masked max
+      !!
+      COMPLEX(wp)::   ctmp
+      REAL(wp)   ::   ztmp
+      INTEGER    ::   ji, jj   ! dummy loop indices
+      INTEGER    ::   zjpi, zjpj ! local variables: size of ptab
+      !!-----------------------------------------------------------------------
+      !
+      zjpi = SIZE(ptab1,1)
+      zjpj = SIZE(ptab1,2)
+      !
+      ztmp = 0.e0
+      ctmp = CMPLX( 0.e0, 0.e0, wp )
+      DO jj = 1, zjpj
+         DO ji =1, zjpi
+         ztmp =  ptab1(ji,jj) * tmask_i(ji,jj)
+         CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
+         ztmp =  ptab2(ji,jj) * tmask_i(ji,jj)
+         CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
+         END DO
+      END DO
+      IF( lk_mpp )   CALL mpp_max( ctmp )   ! max over the global domain
+      glob_max_2d_a = REAL(ctmp,wp)
+      !
+   END FUNCTION glob_max_2d_a   
+
+
+   FUNCTION glob_max_3d_a( ptab1, ptab2 ) 
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_max_3d_a ***
+      !!
+      !! ** Purpose : perform a max on two 3D array in calling DDPDD routine
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab1, ptab2
+      REAL(wp)                               ::   glob_max_3d_a   ! global masked max
+      !!
+      COMPLEX(wp)::   ctmp
+      REAL(wp)   ::   ztmp
+      INTEGER    ::   ji, jj, jk   ! dummy loop indices
+      INTEGER    ::   zjpi, zjpj, zjpk ! local variables: size of ptab
+      !!-----------------------------------------------------------------------
+      !
+      zjpi = SIZE(ptab1,1)
+      zjpj = SIZE(ptab1,2)
+      zjpk = SIZE(ptab1,3)
+      !
+      ztmp = 0.e0
+      ctmp = CMPLX( 0.e0, 0.e0, wp )
+      DO jk = 1, zjpk
+         DO jj = 1, zjpj
+            DO ji =1, zjpi
+            ztmp =  ptab1(ji,jj,jk) * tmask_i(ji,jj)
+            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
+            ztmp =  ptab2(ji,jj,jk) * tmask_i(ji,jj)
+            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
+            END DO
+         END DO    
+      END DO
+      IF( lk_mpp )   CALL mpp_max( ctmp )   ! max over the global domain
+      glob_max_3d_a = REAL(ctmp,wp)
+      !
+   END FUNCTION glob_max_3d_a   
+
+
    SUBROUTINE DDPDD( ydda, yddb )
       !!----------------------------------------------------------------------
