Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90	(revision 4344)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90	(revision 4346)
@@ -34,6 +34,6 @@
    !!PUBLIC   lim_diahsb_rst   ! routine called by ice_init.F90
 
-   REAL(wp), SAVE   ::   frc_sal, frc_vol   ! global forcing trends
-   REAL(wp), SAVE   ::   bg_grpme           ! global ice growth+melt trends
+   REAL(dp) ::   frc_sal, frc_vol   ! global forcing trends
+   REAL(dp) ::   bg_grme            ! global ice growth+melt trends
    REAL(wp) ::   epsi06 = 1.e-6_wp  ! small number
    REAL(wp) ::   epsi03 = 1.e-3_wp  ! small number
@@ -57,14 +57,12 @@
       !! ** Purpose: Compute the ice global heat content, salt content and volume conservation
       !!	
-      !! ** Method : - Compute the deviation of heat content, salt content and volume
-      !!	            at the current time step from their values at nit000
-      !!	            - Compute the contribution of forcing and remove it from these deviations
+      !!---------------------------------------------------------------------------
       !!
-      !!---------------------------------------------------------------------------
-      !!
-      REAL(wp)   ::   zbg_ivo, zbg_svo, zbg_are, zbg_sal ,zbg_tem ,zbg_ihc ,zbg_shc
-      REAL(wp)   ::   zbg_sfx, zbg_emp, zbg_sfx_bri, zbg_sfx_thd, zbg_sfx_res, zbg_sfx_mec 
-      REAL(wp)   ::   z_frc_vol, z_frc_sal, z_bg_grpme 
-      REAL(wp)   ::   zdeltat, z1_area                     !    -     -
+      REAL(dp)   ::   zbg_ivo, zbg_svo, zbg_are, zbg_sal ,zbg_tem ,zbg_ihc ,zbg_shc
+      REAL(dp)   ::   zbg_sfx, zbg_sfx_bri, zbg_sfx_thd, zbg_sfx_res, zbg_sfx_mec 
+      REAL(dp)   ::   zbg_emp, zbg_emp_bog, zbg_emp_lag, zbg_emp_sig, zbg_emp_dyg, zbg_emp_bom, zbg_emp_sum, zbg_emp_res 
+      REAL(dp)   ::   z_frc_vol, z_frc_sal, z_bg_grme 
+      REAL(dp)   ::   z1_area                     !    -     -
+      REAL(dp)   ::   zinda, zindb
       !!---------------------------------------------------------------------------
       IF( nn_timing == 1 )   CALL timing_start('lim_diahsb')
@@ -73,5 +71,7 @@
 
       ! 1/area
-      z1_area = 1.0 / MAX( glob_sum( area(:,:) * tms(:,:) ),epsi06 )
+      z1_area = 1.d0 / MAX( glob_sum( area(:,:) * tms(:,:) ), epsi06 )
+
+      zinda = MAX( 0.d0 , SIGN( 1.d0 , glob_sum( area(:,:) * tms(:,:) ) - epsi06 ) )
       ! ----------------------- !
       ! 1 -  Content variations !
@@ -80,31 +80,39 @@
       zbg_svo = glob_sum( vt_s(:,:) * area(:,:) * tms(:,:) ) ! volume snow
       zbg_are = glob_sum( at_i(:,:) * area(:,:) * tms(:,:) ) ! area
-      zbg_sal = glob_sum( SUM(smv_i(:,:,:),dim=3) * area(:,:) * tms(:,:) ) ! mean salt content
-      zbg_tem = glob_sum( (tm_i(:,:)-rtt) * vt_i(:,:) * area(:,:) * tms(:,:) )  ! mean temp content
+      zbg_sal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )       ! mean salt content
+      zbg_tem = glob_sum( ( tm_i(:,:) - rtt ) * vt_i(:,:) * area(:,:) * tms(:,:) )  ! mean temp content
 
       !zbg_ihc = glob_sum( et_i(:,:) * area(:,:) * tms(:,:) ) / MAX( zbg_ivo,epsi06 ) ! ice heat content
       !zbg_shc = glob_sum( et_s(:,:) * area(:,:) * tms(:,:) ) / MAX( zbg_svo,epsi06 ) ! snow heat content
 
-      zbg_ihc = glob_sum( et_i(:,:) ) ! ice heat content  [10^9 J]
-      zbg_shc = glob_sum( et_s(:,:) ) ! snow heat content [10^9 J]
-
-      zbg_sfx      = glob_sum(      sfx(:,:) * area(:,:) * tms(:,:) ) * z1_area
-      zbg_emp       = glob_sum(       emp(:,:) * area(:,:) * tms(:,:) ) * z1_area
-      zbg_sfx_bri     = glob_sum(     sfx_bri(:,:) * area(:,:) * tms(:,:) ) * z1_area
-      zbg_sfx_thd     = glob_sum(     sfx_thd(:,:) * area(:,:) * tms(:,:) ) * z1_area
-      zbg_sfx_res = glob_sum( sfx_res(:,:) * area(:,:) * tms(:,:) ) * z1_area
-      zbg_sfx_mec = glob_sum( sfx_mec(:,:) * area(:,:) * tms(:,:) ) * z1_area
+      zbg_ihc = glob_sum( et_i(:,:) * 1.e-11 ) ! ice heat content  [10^9*1.e-11 J]
+      zbg_shc = glob_sum( et_s(:,:) * 1.e-11 ) ! snow heat content [10^9*1.e-11 J]
+
+      zbg_emp     = zinda * glob_sum(         emp(:,:) * area(:,:) * tms(:,:) ) * z1_area         * r1_rau0 * rday
+      zbg_emp_bog = zinda * glob_sum( diag_bot_gr(:,:) * area(:,:) * tms(:,:) ) * z1_area * rhoic * r1_rau0 * rday
+      zbg_emp_lag = zinda * glob_sum( diag_lat_gr(:,:) * area(:,:) * tms(:,:) ) * z1_area * rhoic * r1_rau0 * rday
+      zbg_emp_sig = zinda * glob_sum( diag_sni_gr(:,:) * area(:,:) * tms(:,:) ) * z1_area * rhoic * r1_rau0 * rday
+      zbg_emp_dyg = zinda * glob_sum( diag_dyn_gr(:,:) * area(:,:) * tms(:,:) ) * z1_area * rhoic * r1_rau0 * rday
+      zbg_emp_bom = zinda * glob_sum( diag_bot_me(:,:) * area(:,:) * tms(:,:) ) * z1_area * rhoic * r1_rau0 * rday
+      zbg_emp_sum = zinda * glob_sum( diag_sur_me(:,:) * area(:,:) * tms(:,:) ) * z1_area * rhoic * r1_rau0 * rday
+      zbg_emp_res = zinda * glob_sum( diag_res_pr(:,:) * area(:,:) * tms(:,:) ) * z1_area * rhoic * r1_rau0 * rday
+
+      zbg_sfx     = zinda * glob_sum(     sfx(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
+      zbg_sfx_bri = zinda * glob_sum( sfx_bri(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
+      zbg_sfx_thd = zinda * glob_sum( sfx_thd(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
+      zbg_sfx_res = zinda * glob_sum( sfx_res(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
+      zbg_sfx_mec = zinda * glob_sum( sfx_mec(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
       
       ! --------------------------------------------- !
       ! 2 - Trends due to forcing and ice growth/melt !
       ! --------------------------------------------- !
-      z_frc_vol = ( 1.e0 / rau0 ) * glob_sum( - emp(:,:) * area(:,:) * tms(:,:) ) ! volume fluxes
-      z_frc_sal = ( 1.e0 / rau0 ) * glob_sum(  sfx(:,:) * area(:,:) * tms(:,:) ) ! salt fluxes
-      z_bg_grpme = glob_sum( ( diag_bot_gr(:,:) + diag_lat_gr(:,:) + diag_sni_gr(:,:) + diag_dyn_gr(:,:) + &
+      z_frc_vol = r1_rau0 * glob_sum( - emp(:,:) * area(:,:) * tms(:,:) ) ! volume fluxes
+      z_frc_sal = r1_rau0 * glob_sum(   sfx(:,:) * area(:,:) * tms(:,:) ) ! salt fluxes
+      z_bg_grme = glob_sum( ( diag_bot_gr(:,:) + diag_lat_gr(:,:) + diag_sni_gr(:,:) + diag_dyn_gr(:,:) + &
                           &    diag_bot_me(:,:) + diag_sur_me(:,:) + diag_res_pr(:,:) ) * area(:,:) * tms(:,:) ) ! volume fluxes
       !
       frc_vol  = frc_vol  + z_frc_vol  * rdt_ice
       frc_sal  = frc_sal  + z_frc_sal  * rdt_ice
-      bg_grpme = bg_grpme + z_bg_grpme * rdt_ice
+      bg_grme  = bg_grme  + z_bg_grme  * rdt_ice
       
       ! difference
@@ -115,23 +123,33 @@
       ! 3 - Diagnostics writing !
       ! ----------------------- !
-      !zdeltat  = 1.e0 / ( ( kt - nit000 + 1 ) * rdt_ice )
-      !
-      CALL iom_put( 'ibgvoltot',zbg_ivo * rhoic / rau0 * 1.e-9 ) ! m3 (equivalent liquid)         
-      CALL iom_put( 'sbgvoltot',zbg_svo * rhosn / rau0 * 1.e-9 )  ! m3 (equivalent liquid)       
-      CALL iom_put( 'ibgarea',zbg_are * 1.e-6 )         
-      CALL iom_put( 'ibgsaline',zbg_sal / MAX( zbg_ivo,epsi03 ) )         
-      CALL iom_put( 'ibgtemper',zbg_tem / MAX( zbg_ivo,epsi03 ) )         
-      CALL iom_put( 'ibgheatco',zbg_ihc )           ! 10^9 J        
-      CALL iom_put( 'ibgsaltco',zbg_sal * rhoic / rau0 * 1.e-9 ) ! psu*km3        
-      CALL iom_put( 'sbgheatco',zbg_shc )           ! 10^9 J
-      CALL iom_put( 'ibgfrcemp',frc_vol * 1.e-9 )           !vol - forcing : km3 (equivalent liquid)         
-      CALL iom_put( 'ibgfrcsfx',frc_sal * 1.e-9 )          !sal - forcing : psu*km3   
-      CALL iom_put( 'ibgsfx',zbg_sfx *rday)            
-      CALL iom_put( 'ibgemp',zbg_emp *rday)              
-      CALL iom_put( 'ibgsfxbri',zbg_sfx_bri *rday)         
-      CALL iom_put( 'ibgsfxthd',zbg_sfx_thd *rday)         
-      CALL iom_put( 'ibgsfxres',zbg_sfx_res *rday)         
-      CALL iom_put( 'ibgsfxmec',zbg_sfx_mec *rday)         
-      CALL iom_put( 'ibggrpme',bg_grpme * rhoic / rau0 * 1.e-9 ) ! km3 (equivalent liquid)         
+      zindb = MAX( 0.d0 , SIGN( 1.d0 , zbg_ivo - epsi06 ) )
+      !
+      CALL iom_put( 'ibgvoltot' , zbg_ivo * rhoic * r1_rau0 * 1.e-9        )   ! ice volume (km3 equivalent liquid)         
+      CALL iom_put( 'sbgvoltot' , zbg_svo * rhosn * r1_rau0 * 1.e-9        )   ! snw volume (km3 equivalent liquid)       
+      CALL iom_put( 'ibgarea'   , zbg_are * 1.e-6                          )   ! ice area   (km2)
+      CALL iom_put( 'ibgsaline' , zindb * zbg_sal / MAX( zbg_ivo, epsi06 ) )   ! ice saline (psu)
+      CALL iom_put( 'ibgtemper' , zindb * zbg_tem / MAX( zbg_ivo, epsi06 ) )   ! ice temper (C)
+      CALL iom_put( 'ibgheatco' , zbg_ihc                                  )   ! ice heat content (1.e20 J)        
+      CALL iom_put( 'sbgheatco' , zbg_shc                                  )   ! snw heat content (1.e20 J)
+      CALL iom_put( 'ibgsaltco' , zbg_sal * rhoic * r1_rau0 * 1.e-9        )   ! ice salt content (psu*km3 equivalent liquid)        
+
+      CALL iom_put( 'ibgemp'    , zbg_emp                                  )   ! volume flux emp (m/day liquid)
+      CALL iom_put( 'ibgempbog' , zbg_emp_bog                              )   ! volume flux bottom growth     -(m/day equivalent liquid)
+      CALL iom_put( 'ibgemplag' , zbg_emp_lag                              )   ! volume flux open water growth -
+      CALL iom_put( 'ibgempsig' , zbg_emp_sig                              )   ! volume flux snow ice growth   -
+      CALL iom_put( 'ibgempdyg' , zbg_emp_dyg                              )   ! volume flux dynamic growth    -
+      CALL iom_put( 'ibgempbom' , zbg_emp_bom                              )   ! volume flux bottom melt       -
+      CALL iom_put( 'ibgempsum' , zbg_emp_sum                              )   ! volume flux surface melt      -
+      CALL iom_put( 'ibgempres' , zbg_emp_res                              )   ! volume flux resultant         -
+          
+      CALL iom_put( 'ibgsfx'    , zbg_sfx                                  )   ! salt flux         -(psu*m/day equivalent liquid)       
+      CALL iom_put( 'ibgsfxbri' , zbg_sfx_bri                              )   ! salt flux brines  -      
+      CALL iom_put( 'ibgsfxthd' , zbg_sfx_thd                              )   ! salt flux thermo  -    
+      CALL iom_put( 'ibgsfxmec' , zbg_sfx_mec                              )   ! salt flux dynamic -    
+      CALL iom_put( 'ibgsfxres' , zbg_sfx_res                              )   ! salt flux result  -    
+
+      CALL iom_put( 'ibgfrcvol' , frc_vol * 1.e-9                          )   ! vol - forcing     (km3 equivalent liquid) 
+      CALL iom_put( 'ibgfrcsfx' , frc_sal * 1.e-9                          )   ! sal - forcing     (psu*km3 equivalent liquid)   
+      CALL iom_put( 'ibggrme'   , bg_grme * rhoic * r1_rau0 * 1.e-9        )   ! vol growth + melt (km3 equivalent liquid)         
       !
       IF( lrst_ice )   CALL lim_diahsb_rst( numit, 'WRITE' )
@@ -172,7 +190,7 @@
       ! 2 - initial conservation variables !
       ! ---------------------------------- !
-      !frc_vol = 0._wp                                           ! volume       trend due to forcing
-      !frc_sal = 0._wp                                           ! salt content   -    -   -    -         
-      !bg_grpme = 0._wp                                          ! ice growth + melt volume trend
+      !frc_vol = 0.d0                                           ! volume       trend due to forcing
+      !frc_sal = 0.d0                                           ! salt content   -    -   -    -         
+      !bg_grme = 0.d0                                           ! ice growth + melt volume trend
       !
       CALL lim_diahsb_rst( nstart, 'READ' )  !* read or initialize all required files
@@ -196,21 +214,29 @@
      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
         IF( ln_rstart ) THEN                   !* Read the restart file
-           id1 = iom_varid( numrir, 'frc_vol'  , ldstop = .TRUE. )
+           !id1 = iom_varid( numrir, 'frc_vol'  , ldstop = .TRUE. )
            !
+           IF(lwp) WRITE(numout,*) '~~~~~~~'
+           IF(lwp) WRITE(numout,*) ' lim_diahsb_rst at it= ', kt,' date= ', ndastp
+           IF(lwp) WRITE(numout,*) '~~~~~~~'
            CALL iom_get( numrir, 'frc_vol', frc_vol )
            CALL iom_get( numrir, 'frc_sal', frc_sal )
-           CALL iom_get( numrir, 'bg_grpme', bg_grpme )
-       ELSE
-       	   frc_vol = 0._wp                                           
-       	   frc_sal = 0._wp                                                  
-      	   bg_grpme = 0._wp                                        
+           CALL iom_get( numrir, 'bg_grme', bg_grme )
+        ELSE
+           IF(lwp) WRITE(numout,*) '~~~~~~~'
+           IF(lwp) WRITE(numout,*) ' lim_diahsb at initial state '
+           IF(lwp) WRITE(numout,*) '~~~~~~~'
+           frc_vol  = 0.d0                                           
+           frc_sal  = 0.d0                                                  
+           bg_grme  = 0.d0                                        
        ENDIF	
 
      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
         !                                   ! -------------------
-        IF(lwp) WRITE(numout,*) '---- limdia-rst ----'
+        IF(lwp) WRITE(numout,*) '~~~~~~~'
+        IF(lwp) WRITE(numout,*) ' lim_diahsb_rst at it= ', kt,' date= ', ndastp
+        IF(lwp) WRITE(numout,*) '~~~~~~~'
         CALL iom_rstput( kt, nitrst, numriw, 'frc_vol'   , frc_vol     )
         CALL iom_rstput( kt, nitrst, numriw, 'frc_sal'   , frc_sal     )
-        CALL iom_rstput( kt, nitrst, numriw, 'bg_grpme'  , bg_grpme     )
+        CALL iom_rstput( kt, nitrst, numriw, 'bg_grme'   , bg_grme     )
         !
      ENDIF
Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 4344)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 4346)
@@ -1484,5 +1484,5 @@
             END DO
          END DO
-         zmask_glo = glob_sum(zmask)
+         !zmask_glo = glob_sum(zmask)
          !IF( ln_nicep .AND. lwp ) WRITE(numout,*) zmask_glo, ' cells of ice zapped in the ocean '
 
Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 4344)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 4346)
@@ -427,17 +427,17 @@
                !-Calculate stress tensor components zs1 and zs2 
                !-at centre of grid cells (see section 3.5 of CICE user's guide).
-               zs1(ji,jj) = ( zs1(ji,jj) - dtotel*( ( 1._wp - alphaevp) * zs1(ji,jj) +   &
-                  &          ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) * zpresh(ji,jj) ) )  &       
-                  &          / ( 1._wp + alphaevp * dtotel )
-
-               zs2(ji,jj) = ( zs2(ji,jj) - dtotel * ( ( 1._wp - alphaevp ) * ecc2 * zs2(ji,jj) -   &
-                             zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )   &
-                  &          / ( 1._wp + alphaevp * ecc2 * dtotel )
+               !zs1(ji,jj) = ( zs1(ji,jj) - dtotel*( ( 1._wp - alphaevp) * zs1(ji,jj) +   &
+               !   &          ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) * zpresh(ji,jj) ) )  &       
+               !   &          / ( 1._wp + alphaevp * dtotel )
+
+               !zs2(ji,jj) = ( zs2(ji,jj) - dtotel * ( ( 1._wp - alphaevp ) * ecc2 * zs2(ji,jj) -   &
+               !              zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )   &
+               !   &          / ( 1._wp + alphaevp * ecc2 * dtotel )
 
                ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp)
-               !zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( zdd(ji,jj) / deltat(ji,jj) - delta / deltat(ji,jj) )  &
-               !   &         * zpresh(ji,jj) ) ) / ( 1._wp + dtotel )
-               !zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )  &
-               !   &         / ( 1._wp + dtotel )
+               zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( zdd(ji,jj) / deltat(ji,jj) - delta / deltat(ji,jj) )  &
+                  &         * zpresh(ji,jj) ) ) / ( 1._wp + dtotel )
+               zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )  &
+                  &         / ( 1._wp + dtotel )
 
             END DO
@@ -473,12 +473,12 @@
 
                !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide).
-               zs12(ji,jj) = ( zs12(ji,jj) - dtotel * ( (1.0-alphaevp) * ecc2 * zs12(ji,jj) - zds(ji,jj) /  &
-                  &          ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  &
-                  &          / ( 1._wp + alphaevp * ecc2 * dtotel ) 
+               !zs12(ji,jj) = ( zs12(ji,jj) - dtotel * ( (1.0-alphaevp) * ecc2 * zs12(ji,jj) - zds(ji,jj) /  &
+               !   &          ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  &
+               !   &          / ( 1._wp + alphaevp * ecc2 * dtotel ) 
 
                ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp)
-               !zs12(ji,jj) = ( zs12(ji,jj) + dtotel *  &
-               !   &          ( ecci * zds(ji,jj) / ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  &
-               !   &          / ( 1.0 + dtotel ) 
+               zs12(ji,jj) = ( zs12(ji,jj) + dtotel *  &
+                  &          ( ecci * zds(ji,jj) / ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  &
+                  &          / ( 1.0 + dtotel ) 
 
             END DO ! ji
