Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/dom_ice.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/dom_ice.F90	(revision 5058)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/dom_ice.F90	(revision 5059)
@@ -21,5 +21,4 @@
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fcor       !: coriolis coefficient
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   covrai     !: sine of geographic latitude
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   area       !: surface of grid cell 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tms, tmi   !: temperature mask, mask for stress
@@ -43,6 +42,5 @@
       !!-------------------------------------------------------------------
       !
-      ALLOCATE( fcor(jpi,jpj)                   ,      &
-         &      covrai(jpi,jpj) , area(jpi,jpj) ,      &
+      ALLOCATE( fcor(jpi,jpj)   , area(jpi,jpj) ,      &
          &      tms   (jpi,jpj) , tmi (jpi,jpj) ,      &
          &      tmu   (jpi,jpj) , tmv (jpi,jpj) ,      &
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 5058)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 5059)
@@ -380,14 +380,4 @@
    ! thd refers to changes induced by thermodynamics
    ! trp   ''         ''     ''       advection (transport of ice)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   d_a_i_thd  , d_a_i_trp                 !: icefractions                  
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   d_v_s_thd  , d_v_s_trp                 !: snow volume
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   d_v_i_thd  , d_v_i_trp                 !: ice  volume
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   d_smv_i_thd, d_smv_i_trp               !:     
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   d_oa_i_thd , d_oa_i_trp                !:
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   d_e_s_thd  , d_e_s_trp     !:
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   d_e_i_thd  , d_e_i_trp     !:
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   d_u_ice_dyn, d_v_ice_dyn   !: ice velocity 
-
    LOGICAL , PUBLIC                                        ::   ln_limdiahsb  !: flag for ice diag (T) or not (F)
    LOGICAL , PUBLIC                                        ::   ln_limdiaout  !: flag for ice diag (T) or not (F)
@@ -414,5 +404,5 @@
       INTEGER :: ice_alloc
       !
-      INTEGER :: ierr(19), ii
+      INTEGER :: ierr(17), ii
       !!-----------------------------------------------------------------
 
@@ -485,14 +475,4 @@
          &      oa_i_b (jpi,jpj,jpl)                                                        ,     &
          &      u_ice_b(jpi,jpj)     , v_ice_b(jpi,jpj)                                   , STAT=ierr(ii) )
-
-      ! * Increment of global variables
-      ii = ii + 1
-      ALLOCATE( d_a_i_thd(jpi,jpj,jpl) , d_a_i_trp (jpi,jpj,jpl) , d_v_s_thd  (jpi,jpj,jpl) , d_v_s_trp  (jpi,jpj,jpl) ,   &
-         &      d_v_i_thd(jpi,jpj,jpl) , d_v_i_trp (jpi,jpj,jpl) , d_smv_i_thd(jpi,jpj,jpl) , d_smv_i_trp(jpi,jpj,jpl) ,   &     
-         &      d_oa_i_thd(jpi,jpj,jpl) , d_oa_i_trp (jpi,jpj,jpl) ,   &
-         &     STAT=ierr(ii) )
-      ii = ii + 1
-      ALLOCATE( d_e_s_thd(jpi,jpj,nlay_s,jpl) , d_e_i_thd(jpi,jpj,nlay_i,jpl) , d_u_ice_dyn(jpi,jpj) ,     &
-         &      d_e_s_trp(jpi,jpj,nlay_s,jpl) , d_e_i_trp(jpi,jpj,nlay_i,jpl) , d_v_ice_dyn(jpi,jpj) , STAT=ierr(ii) )
       
       ! * Ice thickness distribution variables
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90	(revision 5058)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90	(revision 5059)
@@ -6,5 +6,5 @@
    !! History :   -   ! Original code from William H. Lipscomb, LANL
    !!            3.0  ! 2004-06  (M. Vancoppenolle)   Energy Conservation 
-   !!            4.0  ! 2011-02  (G. Madec)  add mpp considerations
+   !!            3.5  ! 2011-02  (G. Madec)  add mpp considerations
    !!             -   ! 2014-05  (C. Rousset) add lim_cons_hsm
    !!----------------------------------------------------------------------
@@ -207,9 +207,9 @@
             IF ( ABS( zsmv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (',cd_routine,') = ',(zsmv * rday)
             IF ( ABS( zei    ) >  1.    ) WRITE(numout,*) 'violation enthalpy [1e9 J]    (',cd_routine,') = ',(zei)
-            IF ( zvmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',(zvmin)
+            IF ( zvmin <  1.e-10        ) WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',(zvmin)
             IF( cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' .AND. zamax > amax+1.e-10 ) THEN
                                           WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax
             ENDIF
-            IF ( zamin <  0.            ) WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin
+            IF ( zamin <  1.e-10        ) WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin
          ENDIF
 
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90	(revision 5058)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90	(revision 5059)
@@ -77,6 +77,4 @@
                   !WRITE(numout,*) ' a_i *** a_i_b   ', a_i      (ji,jj,jl), a_i_b  (ji,jj,jl)
                   !WRITE(numout,*) ' v_i *** v_i_b   ', v_i      (ji,jj,jl), v_i_b  (ji,jj,jl)
-                  !WRITE(numout,*) ' d_a_i_thd/trp   ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl)
-                  !WRITE(numout,*) ' d_v_i_thd/trp   ', d_v_i_thd(ji,jj,jl), d_v_i_trp(ji,jj,jl)
                   inb_alp(ialert_id) = inb_alp(ialert_id) + 1
                ENDIF
@@ -187,7 +185,5 @@
                   !WRITE(numout,*) ' Category no: ', jl
                   !WRITE(numout,*) ' a_i        : ', a_i      (ji,jj,jl) , ' a_i_b      : ', a_i_b  (ji,jj,jl)   
-                  !WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl) , ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl) 
                   !WRITE(numout,*) ' v_i        : ', v_i      (ji,jj,jl) , ' v_i_b      : ', v_i_b  (ji,jj,jl)   
-                  !WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl) , ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl) 
                   !WRITE(numout,*) ' '
                !END DO
@@ -363,5 +359,4 @@
                WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ji,jj)
                WRITE(numout,*) ' strength      : ', strength(ji,jj)
-               WRITE(numout,*) ' d_u_ice_dyn   : ', d_u_ice_dyn(ji,jj), ' d_v_ice_dyn   : ', d_v_ice_dyn(ji,jj)
                WRITE(numout,*) ' u_ice_b       : ', u_ice_b(ji,jj)    , ' v_ice_b       : ', v_ice_b(ji,jj)  
                WRITE(numout,*)
@@ -375,19 +370,11 @@
                   WRITE(numout,*) ' sm_i       : ', sm_i(ji,jj,jl)             , ' o_i        : ', o_i(ji,jj,jl)
                   WRITE(numout,*) ' a_i        : ', a_i(ji,jj,jl)              , ' a_i_b      : ', a_i_b(ji,jj,jl)   
-                  WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl)        , ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl) 
                   WRITE(numout,*) ' v_i        : ', v_i(ji,jj,jl)              , ' v_i_b      : ', v_i_b(ji,jj,jl)   
-                  WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl)        , ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl) 
                   WRITE(numout,*) ' v_s        : ', v_s(ji,jj,jl)              , ' v_s_b      : ', v_s_b(ji,jj,jl)  
-                  WRITE(numout,*) ' d_v_s_trp  : ', d_v_s_trp(ji,jj,jl)        , ' d_v_s_thd  : ', d_v_s_thd(ji,jj,jl)
                   WRITE(numout,*) ' e_i1       : ', e_i(ji,jj,1,jl)/1.0e9      , ' ei1        : ', e_i_b(ji,jj,1,jl)/1.0e9 
-                  WRITE(numout,*) ' de_i1_trp  : ', d_e_i_trp(ji,jj,1,jl)/1.0e9, ' de_i1_thd  : ', d_e_i_thd(ji,jj,1,jl)/1.0e9
                   WRITE(numout,*) ' e_i2       : ', e_i(ji,jj,2,jl)/1.0e9      , ' ei2_b      : ', e_i_b(ji,jj,2,jl)/1.0e9  
-                  WRITE(numout,*) ' de_i2_trp  : ', d_e_i_trp(ji,jj,2,jl)/1.0e9, ' de_i2_thd  : ', d_e_i_thd(ji,jj,2,jl)/1.0e9
                   WRITE(numout,*) ' e_snow     : ', e_s(ji,jj,1,jl)            , ' e_snow_b   : ', e_s_b(ji,jj,1,jl) 
-                  WRITE(numout,*) ' d_e_s_trp  : ', d_e_s_trp(ji,jj,1,jl)      , ' d_e_s_thd  : ', d_e_s_thd(ji,jj,1,jl)
                   WRITE(numout,*) ' smv_i      : ', smv_i(ji,jj,jl)            , ' smv_i_b    : ', smv_i_b(ji,jj,jl)   
-                  WRITE(numout,*) ' d_smv_i_trp: ', d_smv_i_trp(ji,jj,jl)      , ' d_smv_i_thd: ', d_smv_i_thd(ji,jj,jl) 
                   WRITE(numout,*) ' oa_i       : ', oa_i(ji,jj,jl)             , ' oa_i_b     : ', oa_i_b(ji,jj,jl)
-                  WRITE(numout,*) ' d_oa_i_trp : ', d_oa_i_trp(ji,jj,jl)       , ' d_oa_i_thd : ', d_oa_i_thd(ji,jj,jl)
                END DO !jl
                
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 5058)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 5059)
@@ -6,5 +6,5 @@
    !! history :  1.0  ! 2002-08  (C. Ethe, G. Madec)  original VP code 
    !!            3.0  ! 2007-03  (MA Morales Maqueda, S. Bouillon, M. Vancoppenolle)  LIM3: EVP-Cgrid
-   !!            4.0  ! 2011-02  (G. Madec) dynamical allocation
+   !!            3.5  ! 2011-02  (G. Madec) dynamical allocation
    !!----------------------------------------------------------------------
 #if defined key_lim3
@@ -244,4 +244,6 @@
       INTEGER  ::   ios                 ! Local integer output status for namelist read
       NAMELIST/namicedyn/ cw, pstar, c_rhg, creepl, ecc, ahi0, nevp, relast
+      INTEGER  ::   ji, jj
+      REAL(wp) ::   za00, zd_max
       !!-------------------------------------------------------------------
 
@@ -264,5 +266,5 @@
          WRITE(numout,*) '   creep limit                                      creepl = ', creepl
          WRITE(numout,*) '   eccentricity of the elliptical yield curve       ecc    = ', ecc
-         WRITE(numout,*) '   horizontal diffusivity coeff. for sea-ice        ahi0   = ', ahi0
+         WRITE(numout,*) '   horizontal diffusivity coeff. (orca2 grid)       ahi0   = ', ahi0
          WRITE(numout,*) '   number of iterations for subcycling              nevp   = ', nevp
          WRITE(numout,*) '   ratio of elastic timescale over ice time step    relast = ', relast
@@ -271,9 +273,22 @@
       usecc2 = 1._wp / ( ecc * ecc )
       rhoco  = rau0  * cw
-
+      !
       !  Diffusion coefficients.
-      ahiu(:,:) = ahi0 * umask(:,:,1)
-      ahiv(:,:) = ahi0 * vmask(:,:,1)
-      !
+      zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
+      IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
+           
+      za00 = ahi0 / ( 1.e05_wp )              ! 1.e05 = 100km = max grid space at 60° latitude in orca2
+                                              !                 (i.e. 60° = min latitude for ice cover)  
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1)
+            ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1)
+         END DO
+      END DO
+      !
+      IF(lwp) WRITE(numout,*) ''
+      IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to e1'
+      IF(lwp) WRITE(numout,*) '   maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max
+ 
    END SUBROUTINE lim_dyn_init
 
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 5058)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 5059)
@@ -61,10 +61,10 @@
       !! ** action :
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! number of iteration
+      INTEGER, INTENT(in) ::   kt           ! number of iteration
       !
-      INTEGER  ::   ji, jj, jk, jl, jt   ! dummy loop indices
+      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices
       INTEGER  ::   initad                  ! number of sub-timestep for the advection
       REAL(wp) ::   zcfl , zusnit           !   -      -
-      CHARACTER(len=80) :: cltmp
+      CHARACTER(len=80) ::   cltmp
       !
       REAL(wp), POINTER, DIMENSION(:,:)      ::   zsm, zs0at
@@ -75,7 +75,6 @@
       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zhimax                   ! old ice thickness
       REAL(wp), POINTER, DIMENSION(:,:)      ::   zatold, zeiold, zesold   ! old concentration, enthalpies
-      REAL(wp) :: zdv, zvi, zvs, zsmv, zes, zei
-      !
-      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
+      REAL(wp) ::    zdv, zvi, zvs, zsmv, zes, zei
+      REAL(wp) ::    zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
       !!---------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('limtrp')
@@ -122,8 +121,4 @@
                DO ji = 2, jpim1
                   zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) )
-                  !zhimax(ji,jj,jl) = ( ht_i(ji  ,jj  ,jl) * tmask(ji,  jj  ,1) + ht_i(ji-1,jj-1,jl) * tmask(ji-1,jj-1,1) + ht_i(ji+1,jj+1,jl) * tmask(ji+1,jj+1,1) &
-                  !     &             + ht_i(ji-1,jj  ,jl) * tmask(ji-1,jj  ,1) + ht_i(ji  ,jj-1,jl) * tmask(ji  ,jj-1,1) &
-                  !     &             + ht_i(ji+1,jj  ,jl) * tmask(ji+1,jj  ,1) + ht_i(ji  ,jj+1,jl) * tmask(ji  ,jj+1,1) &
-                  !     &             + ht_i(ji-1,jj+1,jl) * tmask(ji-1,jj+1,1) + ht_i(ji+1,jj-1,jl) * tmask(ji+1,jj-1,1) )
                END DO
             END DO
@@ -139,16 +134,14 @@
          zcfl  = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice / e2v(:,:) ) )
          IF(lk_mpp )   CALL mpp_max( zcfl )
-!!gm more readability:
-!         IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp
-!         ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp
-!         ENDIF
-!!gm end
-         initad = 1 + NINT( MAX( 0._wp, SIGN( 1._wp, zcfl-0.5 ) ) )
-         zusnit = 1.0 / REAL( initad ) 
+
+         IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp
+         ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp
+         ENDIF
+
          IF( zcfl > 0.5_wp .AND. lwp )   ncfl = ncfl + 1
          IF( numit == nlast .AND. lwp ) THEN
             IF( ncfl > 0 ) THEN   
-             WRITE(cltmp,'(i6.1)') ncfl
-             CALL ctl_stop('STOP',TRIM(cltmp) )
+               WRITE(cltmp,'(i6.1)') ncfl
+               CALL ctl_stop('STOP',TRIM(cltmp) )
                CALL ctl_warn( 'lim_trp: ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ')
             ELSE
@@ -160,5 +153,5 @@
          ! transported fields                                        
          !-------------------------
-         zs0ow(:,:,1) = ato_i(:,:) * area(:,:)               ! Open water area 
+         zs0ow(:,:,1) = ato_i(:,:) * area(:,:)              ! Open water area 
          DO jl = 1, jpl
             zs0sn (:,:,jl)   = v_s  (:,:,jl) * area(:,:)    ! Snow volume
@@ -277,5 +270,5 @@
 
          !------------------------------------------------------------------------------!
-         ! 4) Diffusion of Ice fields                  
+         ! Diffusion of Ice fields                  
          !------------------------------------------------------------------------------!
 
@@ -322,11 +315,7 @@
 
          !------------------------------------------------------------------------------!
-         ! 5) Update and limit ice properties after transport                           
+         ! limit ice properties after transport                           
          !------------------------------------------------------------------------------!
-
 !!gm & cr   :  MAX should not be active if adv scheme is positive !
-         !--------------------------------------------------
-         ! 5.1) Recover mean values over the grid squares.
-         !--------------------------------------------------
          DO jl = 1, jpl
             DO jj = 1, jpj
@@ -340,6 +329,5 @@
                END DO
             END DO
-         END DO
-         DO jl = 1, jpl
+
             DO jk = 1, nlay_i
                DO jj = 1, jpj
@@ -458,37 +446,7 @@
 
       ! -------------------------------------------------
-      IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx,-1, ' - ice dyn & trp - ' )   ! control print
+      ! control prints
       ! -------------------------------------------------
-      IF(ln_ctl) THEN   ! Control print
-         CALL prt_ctl_info(' ')
-         CALL prt_ctl_info(' - Cell values : ')
-         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
-         CALL prt_ctl(tab2d_1=area , clinfo1=' lim_trp  : cell area :')
-         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_trp  : at_i      :')
-         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_trp  : vt_i      :')
-         CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_trp  : vt_s      :')
-         DO jl = 1, jpl
-            CALL prt_ctl_info(' ')
-            CALL prt_ctl_info(' - Category : ', ivar1=jl)
-            CALL prt_ctl_info('   ~~~~~~~~~~')
-            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_trp  : a_i      : ')
-            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_trp  : ht_i     : ')
-            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_trp  : ht_s     : ')
-            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_trp  : v_i      : ')
-            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_trp  : v_s      : ')
-            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_trp  : e_s      : ')
-            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_trp  : t_su     : ')
-            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_trp  : t_snow   : ')
-            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_trp  : sm_i     : ')
-            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_trp  : smv_i    : ')
-            DO jk = 1, nlay_i
-               CALL prt_ctl_info(' ')
-               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
-               CALL prt_ctl_info('   ~~~~~~~')
-               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_trp  : t_i      : ')
-               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_trp  : e_i      : ')
-            END DO
-         END DO
-      ENDIF
+      IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx,-1, ' - ice dyn & trp - ' )
       !
       CALL wrk_dealloc( jpi,jpj,           zsm, zs0at, zatold, zeiold, zesold )
@@ -499,4 +457,5 @@
       !
       IF( nn_timing == 1 )  CALL timing_stop('limtrp')
+
    END SUBROUTINE lim_trp
 
@@ -509,5 +468,4 @@
    END SUBROUTINE lim_trp
 #endif
-
    !!======================================================================
 END MODULE limtrp
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 5058)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 5059)
@@ -5,5 +5,5 @@
    !!======================================================================
    !! History :  3.0  !  2006-04  (M. Vancoppenolle) Original code
-   !!            3.6  !  2014-06  (C. Rousset)       Complete rewriting/cleaning
+   !!            3.5  !  2014-06  (C. Rousset)       Complete rewriting/cleaning
    !!----------------------------------------------------------------------
 #if defined key_lim3
@@ -13,9 +13,4 @@
    !!    lim_update1   : computes update of sea-ice global variables from trend terms
    !!----------------------------------------------------------------------
-   USE limrhg          ! ice rheology
-
-   USE dom_oce
-   USE oce             ! dynamics and tracers variables
-   USE in_out_manager
    USE sbc_oce         ! Surface boundary condition: ocean fields
    USE sbc_ice         ! Surface boundary condition: ice fields
@@ -23,29 +18,19 @@
    USE phycst          ! physical constants
    USE ice
-   USE limdyn
-   USE limtrp
-   USE limthd
-   USE limsbc
-   USE limdiahsb
-   USE limwri
-   USE limrst
    USE thd_ice         ! LIM thermodynamic sea-ice variables
    USE limitd_th
-   USE limitd_me
    USE limvar
-   USE prtctl           ! Print control
-   USE lbclnk           ! lateral boundary condition - MPP exchanges
-   USE wrk_nemo         ! work arrays
-   USE lib_fortran     ! glob_sum
-   USE in_out_manager   ! I/O manager
-   USE iom              ! I/O manager
-   USE lib_mpp          ! MPP library
+   USE prtctl          ! Print control
+   USE wrk_nemo        ! work arrays
    USE timing          ! Timing
-   USE limcons        ! conservation tests
+   USE limcons         ! conservation tests
+   USE lib_mpp         ! MPP library
+   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE in_out_manager  ! I/O manager
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   lim_update1   ! routine called by ice_step
+   PUBLIC   lim_update1
 
    !! * Substitutions
@@ -66,10 +51,8 @@
       !!                
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) :: kt    ! number of iteration
+      INTEGER, INTENT(in) ::   kt    ! number of iteration
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
-      INTEGER  ::   i_ice_switch
       REAL(wp) ::   zsal
-      !
-      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
+      REAL(wp) ::   zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
       !!-------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('limupdate1')
@@ -141,6 +124,6 @@
                   smv_i(ji,jj,jl) = sm_i(ji,jj,jl) * v_i(ji,jj,jl)
                   ! salinity stays in bounds
-                  i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) )
-                  smv_i(ji,jj,jl) = i_ice_switch * MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) )
+                  rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) )
+                  smv_i(ji,jj,jl) = rswitch * MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) )
                   ! associated salt flux
                   sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
@@ -150,4 +133,7 @@
       ENDIF
 
+      ! conservation test
+      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+
       ! -------------------------------------------------
       ! Diagnostics
@@ -157,18 +143,16 @@
       END DO
 
-      d_u_ice_dyn(:,:)     = u_ice(:,:)     - u_ice_b(:,:)
-      d_v_ice_dyn(:,:)     = v_ice(:,:)     - v_ice_b(:,:)
-      d_a_i_trp  (:,:,:)   = a_i  (:,:,:)   - a_i_b  (:,:,:)
-      d_v_s_trp  (:,:,:)   = v_s  (:,:,:)   - v_s_b  (:,:,:)  
-      d_v_i_trp  (:,:,:)   = v_i  (:,:,:)   - v_i_b  (:,:,:)   
-      d_e_s_trp  (:,:,:,:) = e_s  (:,:,:,:) - e_s_b  (:,:,:,:)  
-      d_e_i_trp  (:,:,1:nlay_i,:) = e_i  (:,:,1:nlay_i,:) - e_i_b(:,:,1:nlay_i,:)
-      d_oa_i_trp (:,:,:)   = oa_i (:,:,:)   - oa_i_b (:,:,:)
-      d_smv_i_trp(:,:,:)   = 0._wp
-      IF(  num_sal == 2  ) d_smv_i_trp(:,:,:) = smv_i(:,:,:) - smv_i_b(:,:,:)
-
-      ! conservation test
-      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-
+      ! heat content variation (W.m-2)
+      DO jj = 1, jpj
+         DO ji = 1, jpi            
+            diag_heat_dhc(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  & 
+               &                     SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    &
+               &                   ) * unit_fac * r1_rdtice / area(ji,jj)   
+         END DO
+      END DO
+
+      ! -------------------------------------------------
+      ! control prints
+      ! -------------------------------------------------
       IF(ln_ctl) THEN   ! Control print
          CALL prt_ctl_info(' ')
@@ -181,5 +165,4 @@
          CALL prt_ctl(tab2d_1=strength   , clinfo1=' lim_update1  : strength    :')
          CALL prt_ctl(tab2d_1=u_ice      , clinfo1=' lim_update1  : u_ice       :', tab2d_2=v_ice      , clinfo2=' v_ice       :')
-         CALL prt_ctl(tab2d_1=d_u_ice_dyn, clinfo1=' lim_update1  : d_u_ice_dyn :', tab2d_2=d_v_ice_dyn, clinfo2=' d_v_ice_dyn :')
          CALL prt_ctl(tab2d_1=u_ice_b    , clinfo1=' lim_update1  : u_ice_b     :', tab2d_2=v_ice_b    , clinfo2=' v_ice_b     :')
 
@@ -196,26 +179,18 @@
             CALL prt_ctl(tab2d_1=a_i        (:,:,jl)        , clinfo1= ' lim_update1  : a_i         : ')
             CALL prt_ctl(tab2d_1=a_i_b      (:,:,jl)        , clinfo1= ' lim_update1  : a_i_b       : ')
-            CALL prt_ctl(tab2d_1=d_a_i_trp  (:,:,jl)        , clinfo1= ' lim_update1  : d_a_i_trp   : ')
             CALL prt_ctl(tab2d_1=v_i        (:,:,jl)        , clinfo1= ' lim_update1  : v_i         : ')
             CALL prt_ctl(tab2d_1=v_i_b      (:,:,jl)        , clinfo1= ' lim_update1  : v_i_b       : ')
-            CALL prt_ctl(tab2d_1=d_v_i_trp  (:,:,jl)        , clinfo1= ' lim_update1  : d_v_i_trp   : ')
             CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' lim_update1  : v_s         : ')
             CALL prt_ctl(tab2d_1=v_s_b      (:,:,jl)        , clinfo1= ' lim_update1  : v_s_b       : ')
-            CALL prt_ctl(tab2d_1=d_v_s_trp  (:,:,jl)        , clinfo1= ' lim_update1  : d_v_s_trp   : ')
             CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : e_i1        : ')
             CALL prt_ctl(tab2d_1=e_i_b      (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : e_i1_b      : ')
-            CALL prt_ctl(tab2d_1=d_e_i_trp  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : de_i1_trp   : ')
             CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : e_i2        : ')
             CALL prt_ctl(tab2d_1=e_i_b      (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : e_i2_b      : ')
-            CALL prt_ctl(tab2d_1=d_e_i_trp  (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : de_i2_trp   : ')
             CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' lim_update1  : e_snow      : ')
             CALL prt_ctl(tab2d_1=e_s_b      (:,:,1,jl)      , clinfo1= ' lim_update1  : e_snow_b    : ')
-            CALL prt_ctl(tab2d_1=d_e_s_trp  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : d_e_s_trp   : ')
             CALL prt_ctl(tab2d_1=smv_i      (:,:,jl)        , clinfo1= ' lim_update1  : smv_i       : ')
             CALL prt_ctl(tab2d_1=smv_i_b    (:,:,jl)        , clinfo1= ' lim_update1  : smv_i_b     : ')
-            CALL prt_ctl(tab2d_1=d_smv_i_trp(:,:,jl)        , clinfo1= ' lim_update1  : d_smv_i_trp : ')
             CALL prt_ctl(tab2d_1=oa_i       (:,:,jl)        , clinfo1= ' lim_update1  : oa_i        : ')
             CALL prt_ctl(tab2d_1=oa_i_b     (:,:,jl)        , clinfo1= ' lim_update1  : oa_i_b      : ')
-            CALL prt_ctl(tab2d_1=d_oa_i_trp (:,:,jl)        , clinfo1= ' lim_update1  : d_oa_i_trp  : ')
 
             DO jk = 1, nlay_i
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 5058)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 5059)
@@ -5,5 +5,5 @@
    !!======================================================================
    !! History :  3.0  !  2006-04  (M. Vancoppenolle) Original code
-   !!            3.6  !  2014-06  (C. Rousset)       Complete rewriting/cleaning
+   !!            3.5  !  2014-06  (C. Rousset)       Complete rewriting/cleaning
    !!----------------------------------------------------------------------
 #if defined key_lim3
@@ -13,9 +13,4 @@
    !!    lim_update2   : computes update of sea-ice global variables from trend terms
    !!----------------------------------------------------------------------
-   USE limrhg          ! ice rheology
-
-   USE dom_oce
-   USE oce             ! dynamics and tracers variables
-   USE in_out_manager
    USE sbc_oce         ! Surface boundary condition: ocean fields
    USE sbc_ice         ! Surface boundary condition: ice fields
@@ -23,22 +18,16 @@
    USE phycst          ! physical constants
    USE ice
-   USE limdyn
-   USE limtrp
-   USE limthd
-   USE limsbc
-   USE limdiahsb
-   USE limwri
-   USE limrst
    USE thd_ice         ! LIM thermodynamic sea-ice variables
    USE limitd_th
-   USE limitd_me
    USE limvar
-   USE prtctl           ! Print control
-   USE lbclnk           ! lateral boundary condition - MPP exchanges
-   USE wrk_nemo         ! work arrays
-   USE lib_fortran     ! glob_sum
+   USE prtctl          ! Print control
+   USE lbclnk          ! lateral boundary condition - MPP exchanges
+   USE wrk_nemo        ! work arrays
    USE timing          ! Timing
-   USE limcons        ! conservation tests
+   USE limcons         ! conservation tests
    USE limctl
+   USE lib_mpp         ! MPP library
+   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE in_out_manager
 
    IMPLICIT NONE
@@ -64,10 +53,8 @@
       !!
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) :: kt    ! number of iteration
-      INTEGER  ::   ji, jj, jk, jl    ! dummy loop indices
-      INTEGER  ::   i_ice_switch
+      INTEGER, INTENT(in) ::   kt    ! number of iteration
+      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
       REAL(wp) ::   zh, zsal
-      !
-      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
+      REAL(wp) ::   zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
       !!-------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('limupdate2')
@@ -142,5 +129,5 @@
 
       !---------------------
-      ! 2.11) Ice salinity
+      ! Ice salinity
       !---------------------
       IF (  num_sal == 2  ) THEN 
@@ -151,6 +138,6 @@
                   smv_i(ji,jj,jl) = sm_i(ji,jj,jl) * v_i(ji,jj,jl)
                   ! salinity stays in bounds
-                  i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) )
-                  smv_i(ji,jj,jl) = i_ice_switch * MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) ) !+ s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl)
+                  rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) )
+                  smv_i(ji,jj,jl) = rswitch * MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) ) !+ s_i_min * ( 1._wp - rswitch ) * v_i(ji,jj,jl)
                   ! associated salt flux
                   sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
@@ -161,5 +148,5 @@
 
       !------------------------------------------------------------------------------
-      ! 2) Corrections to avoid wrong values                                        |
+      ! Corrections to avoid wrong values                                        |
       !------------------------------------------------------------------------------
       ! Ice drift
@@ -186,4 +173,7 @@
       CALL lim_var_agg(2)             ! aggregate ice thickness categories
 
+      ! conservation test
+      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+
       ! -------------------------------------------------
       ! Diagnostics
@@ -193,26 +183,18 @@
       END DO
       afx_tot = afx_thd + afx_dyn
-
-      d_a_i_thd(:,:,:)   = a_i(:,:,:)   - a_i_b(:,:,:) 
-      d_v_s_thd(:,:,:)   = v_s(:,:,:)   - v_s_b(:,:,:)
-      d_v_i_thd(:,:,:)   = v_i(:,:,:)   - v_i_b(:,:,:)  
-      d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - e_s_b(:,:,:,:) 
-      d_e_i_thd(:,:,1:nlay_i,:) = e_i(:,:,1:nlay_i,:) - e_i_b(:,:,1:nlay_i,:)
-      !?? d_oa_i_thd(:,:,:)  = oa_i (:,:,:) - oa_i_b (:,:,:)
-      d_smv_i_thd(:,:,:) = 0._wp
-      IF( num_sal == 2 )   d_smv_i_thd(:,:,:) = smv_i(:,:,:) - smv_i_b(:,:,:)
 
       ! heat content variation (W.m-2)
       DO jj = 1, jpj
          DO ji = 1, jpi            
-            diag_heat_dhc(ji,jj) = ( SUM( d_e_i_trp(ji,jj,1:nlay_i,:) + d_e_i_thd(ji,jj,1:nlay_i,:) ) +  & 
-               &                     SUM( d_e_s_trp(ji,jj,1:nlay_s,:) + d_e_s_thd(ji,jj,1:nlay_s,:) )    &
+            diag_heat_dhc(ji,jj) = diag_heat_dhc(ji,jj)+   &
+               &                   ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  & 
+               &                     SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    &
                &                   ) * unit_fac * r1_rdtice / area(ji,jj)   
          END DO
       END DO
 
-      ! conservation test
-      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-
+      ! -------------------------------------------------
+      ! control prints
+      ! -------------------------------------------------
       IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 2, ' - Final state - ' )   ! control print
 
@@ -241,26 +223,18 @@
             CALL prt_ctl(tab2d_1=a_i        (:,:,jl)        , clinfo1= ' lim_update2  : a_i         : ')
             CALL prt_ctl(tab2d_1=a_i_b      (:,:,jl)        , clinfo1= ' lim_update2  : a_i_b       : ')
-            CALL prt_ctl(tab2d_1=d_a_i_thd  (:,:,jl)        , clinfo1= ' lim_update2  : d_a_i_thd   : ')
             CALL prt_ctl(tab2d_1=v_i        (:,:,jl)        , clinfo1= ' lim_update2  : v_i         : ')
             CALL prt_ctl(tab2d_1=v_i_b      (:,:,jl)        , clinfo1= ' lim_update2  : v_i_b       : ')
-            CALL prt_ctl(tab2d_1=d_v_i_thd  (:,:,jl)        , clinfo1= ' lim_update2  : d_v_i_thd   : ')
             CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' lim_update2  : v_s         : ')
             CALL prt_ctl(tab2d_1=v_s_b      (:,:,jl)        , clinfo1= ' lim_update2  : v_s_b       : ')
-            CALL prt_ctl(tab2d_1=d_v_s_thd  (:,:,jl)        , clinfo1= ' lim_update2  : d_v_s_thd   : ')
             CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : e_i1        : ')
             CALL prt_ctl(tab2d_1=e_i_b      (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : e_i1_b      : ')
-            CALL prt_ctl(tab2d_1=d_e_i_thd  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : de_i1_thd   : ')
             CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2  : e_i2        : ')
             CALL prt_ctl(tab2d_1=e_i_b      (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2  : e_i2_b      : ')
-            CALL prt_ctl(tab2d_1=d_e_i_thd  (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2  : de_i2_thd   : ')
             CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' lim_update2  : e_snow      : ')
             CALL prt_ctl(tab2d_1=e_s_b      (:,:,1,jl)      , clinfo1= ' lim_update2  : e_snow_b    : ')
-            CALL prt_ctl(tab2d_1=d_e_s_thd  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : d_e_s_thd   : ')
             CALL prt_ctl(tab2d_1=smv_i      (:,:,jl)        , clinfo1= ' lim_update2  : smv_i       : ')
             CALL prt_ctl(tab2d_1=smv_i_b    (:,:,jl)        , clinfo1= ' lim_update2  : smv_i_b     : ')
-            CALL prt_ctl(tab2d_1=d_smv_i_thd(:,:,jl)        , clinfo1= ' lim_update2  : d_smv_i_thd : ')
             CALL prt_ctl(tab2d_1=oa_i       (:,:,jl)        , clinfo1= ' lim_update2  : oa_i        : ')
             CALL prt_ctl(tab2d_1=oa_i_b     (:,:,jl)        , clinfo1= ' lim_update2  : oa_i_b      : ')
-            CALL prt_ctl(tab2d_1=d_oa_i_thd (:,:,jl)        , clinfo1= ' lim_update2  : d_oa_i_thd  : ')
 
             DO jk = 1, nlay_i
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 5058)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 5059)
@@ -183,19 +183,10 @@
          !------------------------------!
          numit = numit + nn_fsbc                     ! Ice model time step
-         !
-         !                                           ! Store previous ice values
-         a_i_b  (:,:,:)   = a_i  (:,:,:)     ! ice area
-         e_i_b  (:,:,:,:) = e_i  (:,:,:,:)   ! ice thermal energy
-         v_i_b  (:,:,:)   = v_i  (:,:,:)     ! ice volume
-         v_s_b  (:,:,:)   = v_s  (:,:,:)     ! snow volume 
-         e_s_b  (:,:,:,:) = e_s  (:,:,:,:)   ! snow thermal energy
-         smv_i_b(:,:,:)   = smv_i(:,:,:)     ! salt content
-         oa_i_b (:,:,:)   = oa_i (:,:,:)     ! areal age content
-         u_ice_b(:,:)     = u_ice(:,:)
-         v_ice_b(:,:)     = v_ice(:,:)
-
-                          CALL sbc_lim_diag0         ! set diag of mass, heat and salt fluxes to 0
-
-                          CALL lim_rst_opn( kt )     ! Open Ice restart file
+         !                                                   
+         CALL sbc_lim_update                ! Store previous ice values
+
+         CALL sbc_lim_diag0                 ! set diag of mass, heat and salt fluxes to 0
+         
+         CALL lim_rst_opn( kt )             ! Open Ice restart file
          !
          ! ----------------------------------------------
@@ -203,84 +194,71 @@
          ! ----------------------------------------------
          IF( .NOT. lk_c1d ) THEN
-
-                          CALL lim_dyn( kt )              ! Ice dynamics    ( rheology/dynamics )
-
-                          CALL lim_trp( kt )              ! Ice transport   ( Advection/diffusion )
-
-         IF( nn_monocat /= 2 )   &
-            &             CALL lim_itd_me                 ! Mechanical redistribution ! (ridging/rafting)
+            
+            CALL lim_dyn( kt )              ! Ice dynamics    ( rheology/dynamics )
+            
+            CALL lim_trp( kt )              ! Ice transport   ( Advection/diffusion )
+            
+            IF( nn_monocat /= 2 ) CALL lim_itd_me  ! Mechanical redistribution ! (ridging/rafting)
 
 #if defined key_bdy
-                          ! bdy ice thermo 
-                          CALL lim_var_glo2eqv            ! equivalent variables
-                          CALL bdy_ice_lim( kt )
-                          CALL lim_var_zapsmall
-                          CALL lim_var_agg(1)
-         IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 1, ' - ice thermo bdy - ' )   ! control print
+            CALL lim_var_glo2eqv
+            CALL bdy_ice_lim( kt )         ! bdy ice thermo 
+            CALL lim_var_zapsmall
+            CALL lim_var_agg(1)
+            IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 1, ' - ice thermo bdy - ' )
 #endif
-
-                          CALL lim_update1( kt )
-
+            CALL lim_update1( kt )
+            
          ENDIF
-
-         !- Change old values for new values
-         u_ice_b(:,:)     = u_ice(:,:)
-         v_ice_b(:,:)     = v_ice(:,:)
-         a_i_b  (:,:,:)   = a_i  (:,:,:)
-         v_s_b  (:,:,:)   = v_s  (:,:,:)
-         v_i_b  (:,:,:)   = v_i  (:,:,:)
-         e_s_b  (:,:,:,:) = e_s  (:,:,:,:)
-         e_i_b  (:,:,:,:) = e_i  (:,:,:,:)
-         oa_i_b (:,:,:)   = oa_i (:,:,:)
-         smv_i_b(:,:,:)   = smv_i(:,:,:)
+         
+         CALL sbc_lim_update                ! Store previous ice values
  
          ! ----------------------------------------------
          ! ice thermodynamics
          ! ----------------------------------------------
-                          CALL lim_var_glo2eqv            ! equivalent variables
-                          CALL lim_var_agg(1)             ! aggregate ice categories
-
-                          ! previous lead fraction and ice volume for flux calculations
-                          pfrld(:,:)   = 1._wp - at_i(:,:)
-                          phicif(:,:)  = vt_i(:,:)
-
-                          SELECT CASE( kblk )
-                             CASE ( jp_cpl )
-                             CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su    )
-                             IF( nn_limflx == 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   &
-                          &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx )
-                             ! Latent heat flux is forced to 0 in coupled: it is included in qns (non-solar heat flux)
-                             qla_ice  (:,:,:) = 0._wp
-                             dqla_ice (:,:,:) = 0._wp
-                          END SELECT
-                          !
-                          CALL lim_thd( kt )              ! Ice thermodynamics 
-
-                          CALL lim_update2( kt )          ! Global variables update
-         !
-                          CALL lim_sbc_flx( kt )          ! Update surface ocean mass, heat and salt fluxes
-         !
-         IF(ln_limdiaout) CALL lim_diahsb                 ! Diagnostics and outputs 
-
-                          CALL lim_wri( 1 )               ! Ice outputs 
-
+         CALL lim_var_glo2eqv
+         CALL lim_var_agg(1)
+         
+         ! previous lead fraction and ice volume for flux calculations
+         pfrld(:,:)   = 1._wp - at_i(:,:)
+         phicif(:,:)  = vt_i(:,:)
+         
+         SELECT CASE( kblk )
+         CASE ( jp_cpl )
+            CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su    )
+            IF( nn_limflx == 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   &
+               &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx )
+            ! Latent heat flux is forced to 0 in coupled: it is included in qns (non-solar heat flux)
+            qla_ice  (:,:,:) = 0._wp
+            dqla_ice (:,:,:) = 0._wp
+         END SELECT
+         !
+         CALL lim_thd( kt )                         ! Ice thermodynamics 
+         
+         CALL lim_update2( kt )                     ! Global variables update
+         !
+         CALL lim_sbc_flx( kt )                     ! Update surface ocean mass, heat and salt fluxes
+         !
+         IF(ln_limdiaout) CALL lim_diahsb           ! Diagnostics and outputs 
+         
+         CALL lim_wri( 1 )                          ! Ice outputs 
+         
          IF( kt == nit000 .AND. ln_rstart )   &
-            &             CALL iom_close( numrir )        ! close input ice restart file
-         !
-         IF( lrst_ice )   CALL lim_rst_write( kt )        ! Ice restart file 
-                          CALL lim_var_glo2eqv            ! ???
-         !
-         IF( ln_nicep )   CALL lim_ctl( kt )              ! alerts in case of model crash
+            &             CALL iom_close( numrir )  ! close input ice restart file
+         !
+         IF( lrst_ice )   CALL lim_rst_write( kt )  ! Ice restart file 
+         CALL lim_var_glo2eqv                       ! ???
+         !
+         IF( ln_nicep )   CALL lim_ctl( kt )        ! alerts in case of model crash
          !
          CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )
          !
-      ENDIF                                    ! End sea-ice time step only
-
-      !                                        !--------------------------!
-      !                                        !  at all ocean time step  !
-      !                                        !--------------------------!
-      !                                               
-      !                                              ! Update surface ocean stresses (only in ice-dynamic case)
-      !                                                   ! otherwise the atm.-ocean stresses are used everywhere
+      ENDIF   ! End sea-ice time step only
+
+      !--------------------------------!
+      ! --- at all ocean time step --- !
+      !--------------------------------!
+      ! Update surface ocean stresses (only in ice-dynamic case)
+      !    otherwise the atm.-ocean stresses are used everywhere
       IF( ln_limdyn )     CALL lim_sbc_tau( kt, ub(:,:,1), vb(:,:,1) )  ! using before instantaneous surf. currents
 !!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!!
@@ -504,5 +482,5 @@
          &                          pdqn_ice, pqla_ice, pdql_ice, k_limflx )
       !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE sbc_ice_lim  ***
+      !!                  ***  ROUTINE ice_lim_flx  ***
       !!                   
       !! ** Purpose :   update the ice surface boundary condition by averaging and / or
@@ -580,7 +558,25 @@
    END SUBROUTINE ice_lim_flx
 
+   SUBROUTINE sbc_lim_update
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_lim_update  ***
+      !!
+      !! ** purpose :  store ice variables at "before" time step 
+      !!----------------------------------------------------------------------
+      a_i_b  (:,:,:)   = a_i  (:,:,:)     ! ice area
+      e_i_b  (:,:,:,:) = e_i  (:,:,:,:)   ! ice thermal energy
+      v_i_b  (:,:,:)   = v_i  (:,:,:)     ! ice volume
+      v_s_b  (:,:,:)   = v_s  (:,:,:)     ! snow volume 
+      e_s_b  (:,:,:,:) = e_s  (:,:,:,:)   ! snow thermal energy
+      smv_i_b(:,:,:)   = smv_i(:,:,:)     ! salt content
+      oa_i_b (:,:,:)   = oa_i (:,:,:)     ! areal age content
+      u_ice_b(:,:)     = u_ice(:,:)
+      v_ice_b(:,:)     = v_ice(:,:)
+      
+   END SUBROUTINE sbc_lim_update
+
    SUBROUTINE sbc_lim_diag0
       !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE sbc_lim_flx0  ***
+      !!                  ***  ROUTINE sbc_lim_diag0  ***
       !!
       !! ** purpose :  set ice-ocean and ice-atm. fluxes to zeros at the beggining
@@ -612,4 +608,6 @@
       afx_tot(:,:) = 0._wp   ;
       afx_dyn(:,:) = 0._wp   ;   afx_thd(:,:) = 0._wp
+
+      diag_heat_dhc(:,:) = 0._wp ;
       
    END SUBROUTINE sbc_lim_diag0
