Changeset 7487 for branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/NEMO/LIM_SRC_3

Timestamp:

2016-12-12T11:53:26+01:00 (8 years ago)

Author:

timgraham

Message:

Merged head of v3.6 Stable into NEMOGCM/NEMO

Location:

branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/NEMO/LIM_SRC_3

Files:

• ice.F90 (modified) (9 diffs)
• limcons.F90 (modified) (1 diff)
• limdiahsb.F90 (modified) (7 diffs)
• limhdf.F90 (modified) (10 diffs)
• limistate.F90 (modified) (6 diffs)
• limitd_me.F90 (modified) (3 diffs)
• limrhg.F90 (modified) (13 diffs)
• limsbc.F90 (modified) (1 diff)
• limthd_dh.F90 (modified) (3 diffs)
• limthd_sal.F90 (modified) (2 diffs)
• limtrp.F90 (modified) (8 diffs)
• limvar.F90 (modified) (7 diffs)
• limwri.F90 (modified) (9 diffs)

branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

@@ -234 +234 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   u_oce, v_oce   !: surface ocean velocity used in ice dynamics
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ahiu , ahiv    !: hor. diffusivity coeff. at U- and V-points [m2/s]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   pahu , pahv    !: ice hor. eddy diffusivity coef. at U- and V-points
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ust2s, hicol   !: friction velocity, ice collection thickness accreted in leads
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   strp1, strp2   !: strength at previous time steps
@@ -244 +243 @@
    !
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sist        !: Average Sea-Ice Surface Temperature [Kelvin]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   icethi      !: total ice thickness (for all categories) (diag only)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_bo        !: Sea-Ice bottom temperature [Kelvin]     
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   frld        !: Leads fraction = 1 - ice fraction
@@ -303 +301 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_res     !: residual heat flux due to correction of ice thickness [W.m-2]
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ftr_ice   !: transmitted solar radiation under ice
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ftr_ice   !: transmitted solar radiation under ice   
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   pahu3D , pahv3D
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rn_amax_2d  !: maximum ice concentration 2d array
 
@@ -320 +319 @@
    !                                                                  !  this is an extensive variable that has to be transported
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   o_i     !: Sea-Ice Age (days)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ov_i    !: Sea-Ice Age times volume per area (days.m)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   oa_i    !: Sea-Ice Age times ice area (days)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   bv_i    !: brine volume
 
    !! Variables summed over all categories, or associated to all the ice in a single grid cell
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   u_ice, v_ice   !: components of the ice velocity (m/s)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tio_u, tio_v   !: components of the ice-ocean stress (N/m2)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vt_i , vt_s    !: ice and snow total volume per unit area (m)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   at_i           !: ice total fractional area (ice concentration)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ato_i          !: =1-at_i ; total open water fractional area
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   et_i , et_s    !: ice and snow total heat content
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ot_i           !: mean age over all categories
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tm_i           !: mean ice temperature over all categories
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bv_i           !: brine volume averaged over all categories
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   smt_i          !: mean sea ice salinity averaged over all categories [PSU]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tm_i         !: mean ice temperature over all categories
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bvm_i        !: brine volume averaged over all categories
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   smt_i        !: mean sea ice salinity averaged over all categories [PSU]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tm_su        !: mean surface temperature over all categories
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htm_i        !: mean ice  thickness over all categories
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htm_s        !: mean snow thickness over all categories
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   om_i         !: mean ice age over all categories
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   t_s        !: Snow temperatures [K]
@@ -429 +430 @@
       ALLOCATE( u_oce    (jpi,jpj) , v_oce    (jpi,jpj) ,                           &
          &      ahiu     (jpi,jpj) , ahiv     (jpi,jpj) ,                           &
-         &      pahu     (jpi,jpj) , pahv     (jpi,jpj) ,                           &
          &      ust2s    (jpi,jpj) , hicol    (jpi,jpj) ,                           &
          &      strp1    (jpi,jpj) , strp2    (jpi,jpj) , strength  (jpi,jpj) ,     &
@@ -436 +436 @@
 
       ii = ii + 1
-      ALLOCATE( sist   (jpi,jpj) , icethi (jpi,jpj) , t_bo   (jpi,jpj) ,                        &
+      ALLOCATE( sist   (jpi,jpj) , t_bo   (jpi,jpj) ,                        &
          &      frld   (jpi,jpj) , pfrld  (jpi,jpj) , phicif (jpi,jpj) ,                        &
          &      wfx_snw(jpi,jpj) , wfx_ice(jpi,jpj) , wfx_sub(jpi,jpj) ,                        &
@@ -442 +442 @@
          &      wfx_res(jpi,jpj) , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) ,     &
          &      afx_tot(jpi,jpj) , afx_thd(jpi,jpj),  afx_dyn(jpi,jpj) ,                        &
-         &      fhtur  (jpi,jpj) , ftr_ice(jpi,jpj,jpl), qlead  (jpi,jpj) ,                     &
-         &      rn_amax_2d(jpi,jpj),                                                            &
-         &      sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) , sfx_sub(jpi,jpj) ,                       &
+         &      fhtur  (jpi,jpj) , ftr_ice(jpi,jpj,jpl), pahu3D(jpi,jpj,jpl+1), pahv3D(jpi,jpj,jpl+1),            &
+         &      qlead  (jpi,jpj) , rn_amax_2d(jpi,jpj),                                         &
+         &      sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) , sfx_sub(jpi,jpj),      &
          &      sfx_bog(jpi,jpj) , sfx_bom(jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) ,    &
          &      hfx_res(jpi,jpj) , hfx_snw(jpi,jpj) , hfx_sub(jpi,jpj) , hfx_err(jpi,jpj) ,     & 
@@ -457 +457 @@
          &      v_s  (jpi,jpj,jpl) , ht_s (jpi,jpj,jpl) , t_su (jpi,jpj,jpl) ,     &
          &      sm_i (jpi,jpj,jpl) , smv_i(jpi,jpj,jpl) , o_i  (jpi,jpj,jpl) ,     &
-         &      ov_i (jpi,jpj,jpl) , oa_i (jpi,jpj,jpl)                      , STAT=ierr(ii) )
-      ii = ii + 1
-      ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) , tio_u(jpi,jpj) , tio_v(jpi,jpj) ,     &
+         &      oa_i (jpi,jpj,jpl) , bv_i (jpi,jpj,jpl) , STAT=ierr(ii) )
+      ii = ii + 1
+      ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) ,      &
          &      vt_i (jpi,jpj) , vt_s (jpi,jpj) , at_i (jpi,jpj) , ato_i(jpi,jpj) ,     &
-         &      et_i (jpi,jpj) , et_s (jpi,jpj) , ot_i (jpi,jpj) , tm_i (jpi,jpj) ,     &
-         &      bv_i (jpi,jpj) , smt_i(jpi,jpj)                                   , STAT=ierr(ii) )
+         &      et_i (jpi,jpj) , et_s (jpi,jpj) , tm_i (jpi,jpj) , bvm_i(jpi,jpj) ,     &
+         &      smt_i(jpi,jpj) , tm_su(jpi,jpj) , htm_i(jpi,jpj) , htm_s(jpi,jpj) ,     &
+         &      om_i (jpi,jpj)                              , STAT=ierr(ii) )
       ii = ii + 1
       ALLOCATE( t_s(jpi,jpj,nlay_s,jpl) , e_s(jpi,jpj,nlay_s,jpl) , STAT=ierr(ii) )
@@ -514 +515 @@
    !!======================================================================
 END MODULE ice
+

branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90

@@ -288 +288 @@
 #if ! defined key_bdy
       ! heat flux
-      zhfx  = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es - SUM( qevap_ice * a_i_b, dim=3 ) )  &
-         &              * e12t * tmask(:,:,1) * zconv ) 
+      zhfx  = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es  &
+      !  &              - SUM( qevap_ice * a_i_b, dim=3 )                           & !!clem: I think this line must be commented (but need check)
+         &              ) * e12t * tmask(:,:,1) * zconv ) 
       ! salt flux
       zsfx  = glob_sum( ( sfx + diag_smvi ) * e12t * tmask(:,:,1) * zconv ) * rday

branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90

@@ -31 +31 @@
 
    PUBLIC   lim_diahsb        ! routine called by ice_step.F90
-
-   real(wp) ::   frc_sal, frc_vol   ! global forcing trends
-   real(wp) ::   bg_grme            ! global ice growth+melt trends
-
+   PUBLIC   lim_diahsb_init   ! routine called in sbcice_lim.F90
+
+   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   vol_loc_ini, sal_loc_ini, tem_loc_ini ! initial volume, salt and heat contents
+   REAL(wp)                              ::   frc_sal, frc_voltop, frc_volbot, frc_temtop, frc_tembot  ! global forcing trends
+   
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -46 +47 @@
 CONTAINS
 
-   SUBROUTINE lim_diahsb
+   SUBROUTINE lim_diahsb( kt )
       !!---------------------------------------------------------------------------
       !!                  ***  ROUTINE lim_diahsb  ***
@@ -53 +54 @@
       !! 
       !!---------------------------------------------------------------------------
+      INTEGER, INTENT(in) :: kt    ! number of iteration
       !!
-      real(wp)   ::   zbg_ivo, zbg_svo, zbg_are, zbg_sal ,zbg_tem ,zbg_ihc ,zbg_shc
-      real(wp)   ::   zbg_sfx, zbg_sfx_bri, zbg_sfx_bog, zbg_sfx_bom, zbg_sfx_sum, zbg_sfx_sni,   &
-      &               zbg_sfx_opw, zbg_sfx_res, zbg_sfx_dyn, zbg_sfx_sub 
-      real(wp)   ::   zbg_vfx, zbg_vfx_bog, zbg_vfx_opw, zbg_vfx_sni, zbg_vfx_dyn
-      real(wp)   ::   zbg_vfx_bom, zbg_vfx_sum, zbg_vfx_res, zbg_vfx_spr, zbg_vfx_snw, zbg_vfx_sub  
-      real(wp)   ::   zbg_hfx_dhc, zbg_hfx_spr
-      real(wp)   ::   zbg_hfx_res, zbg_hfx_sub, zbg_hfx_dyn, zbg_hfx_thd, zbg_hfx_snw, zbg_hfx_out, zbg_hfx_in   
-      real(wp)   ::   zbg_hfx_sum, zbg_hfx_bom, zbg_hfx_bog, zbg_hfx_dif, zbg_hfx_opw
-      real(wp)   ::   z_frc_vol, z_frc_sal, z_bg_grme 
-      real(wp)   ::   z1_area                     !    -     -
-      REAL(wp)   ::   ztmp
+      real(wp)   ::   zbg_ivol, zbg_svol, zbg_area, zbg_isal, zbg_item ,zbg_stem
+      REAL(wp)   ::   z_frc_voltop, z_frc_volbot, z_frc_sal, z_frc_temtop, z_frc_tembot  
+      REAL(wp)   ::   zdiff_vol, zdiff_sal, zdiff_tem  
       !!---------------------------------------------------------------------------
       IF( nn_timing == 1 )   CALL timing_start('lim_diahsb')
 
-      IF( numit == nstart ) CALL lim_diahsb_init 
-
-      ! 1/area
-      z1_area = 1._wp / MAX( glob_sum( e12t(:,:) * tmask(:,:,1) ), epsi06 )
-
-      rswitch = MAX( 0._wp , SIGN( 1._wp , glob_sum( e12t(:,:) * tmask(:,:,1) ) - epsi06 ) )
-      ! ----------------------- !
-      ! 1 -  Content variations !
-      ! ----------------------- !
-      zbg_ivo = glob_sum( vt_i(:,:) * e12t(:,:) * tmask(:,:,1) ) ! volume ice 
-      zbg_svo = glob_sum( vt_s(:,:) * e12t(:,:) * tmask(:,:,1) ) ! volume snow
-      zbg_are = glob_sum( at_i(:,:) * e12t(:,:) * tmask(:,:,1) ) ! area
-      zbg_sal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) )       ! mean salt content
-      zbg_tem = glob_sum( ( tm_i(:,:) - rt0 ) * vt_i(:,:) * e12t(:,:) * tmask(:,:,1) )  ! mean temp content
-
-      !zbg_ihc = glob_sum( et_i(:,:) * e12t(:,:) * tmask(:,:,1) ) / MAX( zbg_ivo,epsi06 ) ! ice heat content
-      !zbg_shc = glob_sum( et_s(:,:) * e12t(:,:) * tmask(:,:,1) ) / MAX( zbg_svo,epsi06 ) ! snow heat content
-
-      ! Volume
-      ztmp = rswitch * z1_area * r1_rau0 * rday
-      zbg_vfx     = ztmp * glob_sum(     emp(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_vfx_bog = ztmp * glob_sum( wfx_bog(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_vfx_opw = ztmp * glob_sum( wfx_opw(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_vfx_sni = ztmp * glob_sum( wfx_sni(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_vfx_dyn = ztmp * glob_sum( wfx_dyn(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_vfx_bom = ztmp * glob_sum( wfx_bom(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_vfx_sum = ztmp * glob_sum( wfx_sum(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_vfx_res = ztmp * glob_sum( wfx_res(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_vfx_spr = ztmp * glob_sum( wfx_spr(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_vfx_snw = ztmp * glob_sum( wfx_snw(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_vfx_sub = ztmp * glob_sum( wfx_sub(:,:) * e12t(:,:) * tmask(:,:,1) )
-
-      ! Salt
-      zbg_sfx     = ztmp * glob_sum(     sfx(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_sfx_bri = ztmp * glob_sum( sfx_bri(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_sfx_res = ztmp * glob_sum( sfx_res(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_sfx_dyn = ztmp * glob_sum( sfx_dyn(:,:) * e12t(:,:) * tmask(:,:,1) )
-
-      zbg_sfx_bog = ztmp * glob_sum( sfx_bog(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_sfx_opw = ztmp * glob_sum( sfx_opw(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_sfx_sni = ztmp * glob_sum( sfx_sni(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_sfx_bom = ztmp * glob_sum( sfx_bom(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_sfx_sum = ztmp * glob_sum( sfx_sum(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_sfx_sub = ztmp * glob_sum( sfx_sub(:,:) * e12t(:,:) * tmask(:,:,1) )
-
-      ! Heat budget
-      zbg_ihc      = glob_sum( et_i(:,:) * e12t(:,:) * 1.e-20 ) ! ice heat content  [1.e20 J]
-      zbg_shc      = glob_sum( et_s(:,:) * e12t(:,:) * 1.e-20 ) ! snow heat content [1.e20 J]
-      zbg_hfx_dhc  = glob_sum( diag_heat(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-      zbg_hfx_spr  = glob_sum( hfx_spr(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-
-      zbg_hfx_thd  = glob_sum( hfx_thd(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-      zbg_hfx_dyn  = glob_sum( hfx_dyn(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-      zbg_hfx_res  = glob_sum( hfx_res(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-      zbg_hfx_sub  = glob_sum( hfx_sub(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-      zbg_hfx_snw  = glob_sum( hfx_snw(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-      zbg_hfx_sum  = glob_sum( hfx_sum(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-      zbg_hfx_bom  = glob_sum( hfx_bom(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-      zbg_hfx_bog  = glob_sum( hfx_bog(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-      zbg_hfx_dif  = glob_sum( hfx_dif(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-      zbg_hfx_opw  = glob_sum( hfx_opw(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-      zbg_hfx_out  = glob_sum( hfx_out(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-      zbg_hfx_in   = glob_sum(  hfx_in(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-    
-      ! --------------------------------------------- !
-      ! 2 - Trends due to forcing and ice growth/melt !
-      ! --------------------------------------------- !
-      z_frc_vol = r1_rau0 * glob_sum( - emp(:,:) * e12t(:,:) * tmask(:,:,1) ) ! volume fluxes
-      z_frc_sal = r1_rau0 * glob_sum(   sfx(:,:) * e12t(:,:) * tmask(:,:,1) ) ! salt fluxes
-      z_bg_grme = glob_sum( - ( wfx_bog(:,:) + wfx_opw(:,:) + wfx_sni(:,:) + wfx_dyn(:,:) + &
-                          &     wfx_bom(:,:) + wfx_sum(:,:) + wfx_res(:,:) + wfx_snw(:,:) + &
-                          &     wfx_sub(:,:) ) * e12t(:,:) * tmask(:,:,1) ) ! volume fluxes
-      !
-      frc_vol  = frc_vol  + z_frc_vol  * rdt_ice
-      frc_sal  = frc_sal  + z_frc_sal  * rdt_ice
-      bg_grme  = bg_grme  + z_bg_grme  * rdt_ice
+      ! ----------------------- !
+      ! 1 -  Contents !
+      ! ----------------------- !
+      zbg_ivol = glob_sum( vt_i(:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-9 )                  ! ice volume (km3)
+      zbg_svol = glob_sum( vt_s(:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-9 )                  ! snow volume (km3)
+      zbg_area = glob_sum( at_i(:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-6 )                  ! area (km2)
+      zbg_isal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) * 1.e-9 ) ! salt content (pss*km3)
+      zbg_item = glob_sum( et_i * e12t(:,:) * tmask(:,:,1) * 1.e-20 )                      ! heat content (1.e20 J)
+      zbg_stem = glob_sum( et_s * e12t(:,:) * tmask(:,:,1) * 1.e-20 )                      ! heat content (1.e20 J)
       
-      ! difference
-      !frc_vol = zbg_ivo - frc_vol
-      !frc_sal = zbg_sal - frc_sal
-      
-      ! ----------------------- !
-      ! 3 - Diagnostics writing !
-      ! ----------------------- !
-      rswitch = MAX( 0._wp , SIGN( 1._wp , zbg_ivo - epsi06 ) )
-      !
-      IF( iom_use('ibgvoltot') )   &
-      CALL iom_put( 'ibgvoltot' , zbg_ivo * rhoic * r1_rau0 * 1.e-9        )   ! ice volume (km3 equivalent liquid)         
-      IF( iom_use('sbgvoltot') )   &
-      CALL iom_put( 'sbgvoltot' , zbg_svo * rhosn * r1_rau0 * 1.e-9        )   ! snw volume (km3 equivalent liquid)       
-      IF( iom_use('ibgarea') )   &
-      CALL iom_put( 'ibgarea'   , zbg_are * 1.e-6                          )   ! ice area   (km2)
-      IF( iom_use('ibgsaline') )   &
-      CALL iom_put( 'ibgsaline' , rswitch * zbg_sal / MAX( zbg_ivo, epsi06 ) )   ! ice saline (psu)
-      IF( iom_use('ibgtemper') )   &
-      CALL iom_put( 'ibgtemper' , rswitch * 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)
-      IF( iom_use('ibgsaltco') )   &
-      CALL iom_put( 'ibgsaltco' , zbg_sal * rhoic * r1_rau0 * 1.e-9        )   ! ice salt content (psu*km3 equivalent liquid)        
-
-      CALL iom_put( 'ibgvfx'    , zbg_vfx                                  )   ! volume flux emp (m/day liquid)
-      CALL iom_put( 'ibgvfxbog' , zbg_vfx_bog                              )   ! volume flux bottom growth     -(m/day equivalent liquid)
-      CALL iom_put( 'ibgvfxopw' , zbg_vfx_opw                              )   ! volume flux open water growth -
-      CALL iom_put( 'ibgvfxsni' , zbg_vfx_sni                              )   ! volume flux snow ice growth   -
-      CALL iom_put( 'ibgvfxdyn' , zbg_vfx_dyn                              )   ! volume flux dynamic growth    -
-      CALL iom_put( 'ibgvfxbom' , zbg_vfx_bom                              )   ! volume flux bottom melt       -
-      CALL iom_put( 'ibgvfxsum' , zbg_vfx_sum                              )   ! volume flux surface melt      -
-      CALL iom_put( 'ibgvfxres' , zbg_vfx_res                              )   ! volume flux resultant         -
-      CALL iom_put( 'ibgvfxspr' , zbg_vfx_spr                              )   ! volume flux from snow precip         -
-      CALL iom_put( 'ibgvfxsnw' , zbg_vfx_snw                              )   ! volume flux from snow melt         -
-      CALL iom_put( 'ibgvfxsub' , zbg_vfx_sub                              )   ! volume flux from sublimation         -
-          
-      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( 'ibgsfxdyn' , zbg_sfx_dyn                              )   ! salt flux dynamic -    
-      CALL iom_put( 'ibgsfxres' , zbg_sfx_res                              )   ! salt flux result  -    
-      CALL iom_put( 'ibgsfxbog' , zbg_sfx_bog                              )   ! salt flux bottom growth   
-      CALL iom_put( 'ibgsfxopw' , zbg_sfx_opw                              )   ! salt flux open water growth -
-      CALL iom_put( 'ibgsfxsni' , zbg_sfx_sni                              )   ! salt flux snow ice growth   -
-      CALL iom_put( 'ibgsfxbom' , zbg_sfx_bom                              )   ! salt flux bottom melt       -
-      CALL iom_put( 'ibgsfxsum' , zbg_sfx_sum                              )   ! salt flux surface melt      -
-      CALL iom_put( 'ibgsfxsub' , zbg_sfx_sub                              )   ! salt flux sublimation      -
-
-      CALL iom_put( 'ibghfxdhc' , zbg_hfx_dhc                              )   ! Heat content variation in snow and ice [W]
-      CALL iom_put( 'ibghfxspr' , zbg_hfx_spr                              )   ! Heat content of snow precip [W]
-
-      CALL iom_put( 'ibghfxres' , zbg_hfx_res                              )   ! 
-      CALL iom_put( 'ibghfxsub' , zbg_hfx_sub                              )   ! 
-      CALL iom_put( 'ibghfxdyn' , zbg_hfx_dyn                              )   ! 
-      CALL iom_put( 'ibghfxthd' , zbg_hfx_thd                              )   ! 
-      CALL iom_put( 'ibghfxsnw' , zbg_hfx_snw                              )   ! 
-      CALL iom_put( 'ibghfxsum' , zbg_hfx_sum                              )   ! 
-      CALL iom_put( 'ibghfxbom' , zbg_hfx_bom                              )   ! 
-      CALL iom_put( 'ibghfxbog' , zbg_hfx_bog                              )   ! 
-      CALL iom_put( 'ibghfxdif' , zbg_hfx_dif                              )   ! 
-      CALL iom_put( 'ibghfxopw' , zbg_hfx_opw                              )   ! 
-      CALL iom_put( 'ibghfxout' , zbg_hfx_out                              )   ! 
-      CALL iom_put( 'ibghfxin'  , zbg_hfx_in                               )   ! 
-
-      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)   
-      IF( iom_use('ibgvolgrm') )   &
-      CALL iom_put( 'ibgvolgrm' , bg_grme * r1_rau0 * 1.e-9                )   ! vol growth + melt (km3 equivalent liquid)         
-
+      ! ---------------------------!
+      ! 2 - Trends due to forcing  !
+      ! ---------------------------!
+      z_frc_volbot = r1_rau0 * glob_sum( - ( wfx_ice(:,:) + wfx_snw(:,:) + wfx_err_sub(:,:) ) * e12t(:,:) * tmask(:,:,1) * 1.e-9 )  ! freshwater flux ice/snow-ocean 
+      z_frc_voltop = r1_rau0 * glob_sum( - ( wfx_sub(:,:) + wfx_spr(:,:) ) * e12t(:,:) * tmask(:,:,1) * 1.e-9 )                     ! freshwater flux ice/snow-atm
+      z_frc_sal    = r1_rau0 * glob_sum( - sfx(:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-9 )                                            ! salt fluxes ice/snow-ocean
+      z_frc_tembot =           glob_sum( hfx_out(:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-20 )                                         ! heat on top of ocean (and below ice)
+      z_frc_temtop =           glob_sum( hfx_in (:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-20 )                                         ! heat on top of ice-coean
+      !
+      frc_voltop  = frc_voltop  + z_frc_voltop  * rdt_ice ! km3
+      frc_volbot  = frc_volbot  + z_frc_volbot  * rdt_ice ! km3
+      frc_sal     = frc_sal     + z_frc_sal     * rdt_ice ! km3*pss
+      frc_temtop  = frc_temtop  + z_frc_temtop  * rdt_ice ! 1.e20 J
+      frc_tembot  = frc_tembot  + z_frc_tembot  * rdt_ice ! 1.e20 J
+            
+      ! ----------------------- !
+      ! 3 -  Content variations !
+      ! ----------------------- !
+      zdiff_vol = r1_rau0 * glob_sum( ( rhoic * vt_i(:,:) + rhosn * vt_s(:,:) - vol_loc_ini(:,:)  &  ! freshwater trend (km3) 
+         &                            ) * e12t(:,:) * tmask(:,:,1) * 1.e-9 ) 
+      zdiff_sal = r1_rau0 * glob_sum( ( rhoic * SUM( smv_i(:,:,:), dim=3 ) - sal_loc_ini(:,:)     &  ! salt content trend (km3*pss)
+         &                            ) * e12t(:,:) * tmask(:,:,1) * 1.e-9 )
+      zdiff_tem =           glob_sum( ( et_i(:,:) + et_s(:,:) - tem_loc_ini(:,:)                  &  ! heat content trend (1.e20 J)
+      !  &                            + SUM( qevap_ice * a_i_b, dim=3 ) &     !! clem: I think this line should be commented (but needs a check)
+         &                            ) * e12t(:,:) * tmask(:,:,1) * 1.e-20 )
+
+      ! ----------------------- !
+      ! 4 -  Drifts             !
+      ! ----------------------- !
+      zdiff_vol = zdiff_vol - ( frc_voltop + frc_volbot )
+      zdiff_sal = zdiff_sal - frc_sal
+      zdiff_tem = zdiff_tem - ( frc_tembot - frc_temtop )
+
+      ! ----------------------- !
+      ! 5 - Diagnostics writing !
+      ! ----------------------- !
+      !
+      IF( iom_use('ibgvolume') )  CALL iom_put( 'ibgvolume' , zdiff_vol        )   ! ice/snow volume  drift            (km3 equivalent ocean water)         
+      IF( iom_use('ibgsaltco') )  CALL iom_put( 'ibgsaltco' , zdiff_sal        )   ! ice salt content drift            (psu*km3 equivalent ocean water)
+      IF( iom_use('ibgheatco') )  CALL iom_put( 'ibgheatco' , zdiff_tem        )   ! ice/snow heat content drift       (1.e20 J)
+      IF( iom_use('ibgheatfx') )  CALL iom_put( 'ibgheatfx' , zdiff_tem /      &   ! ice/snow heat flux drift          (W/m2)
+         &                                                    glob_sum( e12t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) )
+
+      IF( iom_use('ibgfrcvoltop') )  CALL iom_put( 'ibgfrcvoltop' , frc_voltop )   ! vol  forcing ice/snw-atm          (km3 equivalent ocean water) 
+      IF( iom_use('ibgfrcvolbot') )  CALL iom_put( 'ibgfrcvolbot' , frc_volbot )   ! vol  forcing ice/snw-ocean        (km3 equivalent ocean water) 
+      IF( iom_use('ibgfrcsal') )     CALL iom_put( 'ibgfrcsal'    , frc_sal    )   ! sal - forcing                     (psu*km3 equivalent ocean water)   
+      IF( iom_use('ibgfrctemtop') )  CALL iom_put( 'ibgfrctemtop' , frc_temtop )   ! heat on top of ice/snw/ocean      (1.e20 J)   
+      IF( iom_use('ibgfrctembot') )  CALL iom_put( 'ibgfrctembot' , frc_tembot )   ! heat on top of ocean(below ice)   (1.e20 J)   
+      IF( iom_use('ibgfrchfxtop') )  CALL iom_put( 'ibgfrchfxtop' , frc_temtop / & ! heat on top of ice/snw/ocean      (W/m2) 
+         &                                                    glob_sum( e12t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) )
+      IF( iom_use('ibgfrchfxbot') )  CALL iom_put( 'ibgfrchfxbot' , frc_tembot / & ! heat on top of ocean(below ice)   (W/m2) 
+         &                                                    glob_sum( e12t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) )
+
+      IF( iom_use('ibgvol_tot' ) )  CALL iom_put( 'ibgvol_tot'  , zbg_ivol     )   ! ice volume                        (km3)
+      IF( iom_use('sbgvol_tot' ) )  CALL iom_put( 'sbgvol_tot'  , zbg_svol     )   ! snow volume                       (km3)
+      IF( iom_use('ibgarea_tot') )  CALL iom_put( 'ibgarea_tot' , zbg_area     )   ! ice area                          (km2)
+      IF( iom_use('ibgsalt_tot') )  CALL iom_put( 'ibgsalt_tot' , zbg_isal     )   ! ice salinity content              (pss*km3)
+      IF( iom_use('ibgheat_tot') )  CALL iom_put( 'ibgheat_tot' , zbg_item     )   ! ice heat content                  (1.e20 J)
+      IF( iom_use('sbgheat_tot') )  CALL iom_put( 'sbgheat_tot' , zbg_stem     )   ! snow heat content                 (1.e20 J)
       !
       IF( lrst_ice )   CALL lim_diahsb_rst( numit, 'WRITE' )
       !
       IF( nn_timing == 1 )   CALL timing_stop('lim_diahsb')
-!
+      !
    END SUBROUTINE lim_diahsb
 
@@ -233 +151 @@
       !!             - Compute coefficients for conversion
       !!---------------------------------------------------------------------------
-      INTEGER            ::   jk       ! dummy loop indice
       INTEGER            ::   ierror   ! local integer
       !!
@@ -247 +164 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
       ENDIF
-      !
+      !      
+      ALLOCATE( vol_loc_ini(jpi,jpj), sal_loc_ini(jpi,jpj), tem_loc_ini(jpi,jpj), STAT=ierror )
+      IF( ierror > 0 )  THEN
+         CALL ctl_stop( 'lim_diahsb: unable to allocate vol_loc_ini' )
+         RETURN
+      ENDIF
+
       CALL lim_diahsb_rst( nstart, 'READ' )  !* read or initialize all required files
       !
@@ -263 +186 @@
      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
      !
-     INTEGER ::   id1, id2, id3   ! local integers
      !!----------------------------------------------------------------------
      !
      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
         IF( ln_rstart ) THEN                   !* Read the restart file
-           !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_grme', bg_grme )
+           IF(lwp) WRITE(numout,*) ' lim_diahsb_rst read at it= ', kt,' date= ', ndastp
+           IF(lwp) WRITE(numout,*) '~~~~~~~'
+           CALL iom_get( numrir, 'frc_voltop' , frc_voltop  )
+           CALL iom_get( numrir, 'frc_volbot' , frc_volbot  )
+           CALL iom_get( numrir, 'frc_temtop' , frc_temtop  )
+           CALL iom_get( numrir, 'frc_tembot' , frc_tembot  )
+           CALL iom_get( numrir, 'frc_sal'    , frc_sal     )
+           CALL iom_get( numrir, jpdom_autoglo, 'vol_loc_ini', vol_loc_ini )
+           CALL iom_get( numrir, jpdom_autoglo, 'tem_loc_ini', tem_loc_ini )
+           CALL iom_get( numrir, jpdom_autoglo, 'sal_loc_ini', sal_loc_ini )
         ELSE
            IF(lwp) WRITE(numout,*) '~~~~~~~'
            IF(lwp) WRITE(numout,*) ' lim_diahsb at initial state '
            IF(lwp) WRITE(numout,*) '~~~~~~~'
-           frc_vol  = 0._wp                                          
-           frc_sal  = 0._wp                                                 
-           bg_grme  = 0._wp                                       
+           ! set trends to 0
+           frc_voltop  = 0._wp                                          
+           frc_volbot  = 0._wp                                          
+           frc_temtop  = 0._wp                                                 
+           frc_tembot  = 0._wp                                                 
+           frc_sal     = 0._wp                                                 
+           ! record initial ice volume, salt and temp
+           vol_loc_ini(:,:) = rhoic * vt_i(:,:) + rhosn * vt_s(:,:)  ! ice/snow volume (kg/m2)
+           tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:)                  ! ice/snow heat content (J)
+           sal_loc_ini(:,:) = rhoic * SUM( smv_i(:,:,:), dim=3 )     ! ice salt content (pss*kg/m2)
+           
        ENDIF
 
@@ -288 +222 @@
         !                                   ! -------------------
         IF(lwp) WRITE(numout,*) '~~~~~~~'
-        IF(lwp) WRITE(numout,*) ' lim_diahsb_rst at it= ', kt,' date= ', ndastp
+        IF(lwp) WRITE(numout,*) ' lim_diahsb_rst write 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_grme'   , bg_grme     )
+        CALL iom_rstput( kt, nitrst, numriw, 'frc_voltop' , frc_voltop  )
+        CALL iom_rstput( kt, nitrst, numriw, 'frc_volbot' , frc_volbot  )
+        CALL iom_rstput( kt, nitrst, numriw, 'frc_temtop' , frc_temtop  )
+        CALL iom_rstput( kt, nitrst, numriw, 'frc_tembot' , frc_tembot  )
+        CALL iom_rstput( kt, nitrst, numriw, 'frc_sal'    , frc_sal     )
+        CALL iom_rstput( kt, nitrst, numriw, 'vol_loc_ini', vol_loc_ini )
+        CALL iom_rstput( kt, nitrst, numriw, 'tem_loc_ini', tem_loc_ini )
+        CALL iom_rstput( kt, nitrst, numriw, 'sal_loc_ini', sal_loc_ini )
         !
      ENDIF

branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90

@@ -7 +7 @@
    !!             -   !  2001-05 (G. Madec, R. Hordoir) opa norm
    !!            1.0  !  2002-08 (C. Ethe)  F90, free form
+   !!            3.0  !  2015-08 (O. Tintó and M. Castrillo)  added lim_hdf (multiple)
    !!----------------------------------------------------------------------
 #if defined key_lim3
@@ -27 +28 @@
    PRIVATE
 
-   PUBLIC   lim_hdf         ! called by lim_trp
+   PUBLIC   lim_hdf ! called by lim_trp
    PUBLIC   lim_hdf_init    ! called by sbc_lim_init
 
@@ -43 +44 @@
 CONTAINS
 
-   SUBROUTINE lim_hdf( ptab )
+   SUBROUTINE lim_hdf( ptab , ihdf_vars , jpl , nlay_i )
       !!-------------------------------------------------------------------
       !!                  ***  ROUTINE lim_hdf  ***
@@ -54 +55 @@
       !! ** Action  :    update ptab with the diffusive contribution
       !!-------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) ::   ptab    ! Field on which the diffusion is applied
-      !
-      INTEGER                           ::  ji, jj                    ! dummy loop indices
+      INTEGER                           :: jpl, nlay_i, isize, ihdf_vars
+      REAL(wp),  DIMENSION(:,:,:), INTENT( inout ),TARGET ::   ptab    ! Field on which the diffusion is applied
+      !
+      INTEGER                           ::  ji, jj, jk, jl , jm               ! dummy loop indices
       INTEGER                           ::  iter, ierr           ! local integers
-      REAL(wp)                          ::  zrlxint, zconv     ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:) ::  zrlx, zflu, zflv, zdiv0, zdiv, ztab0
+      REAL(wp)                          ::  zrlxint     ! local scalars
+      REAL(wp), POINTER , DIMENSION ( : )        :: zconv     ! local scalars
+      REAL(wp), POINTER , DIMENSION(:,:,:) ::  zrlx,zdiv0, ztab0
+      REAL(wp), POINTER , DIMENSION(:,:) ::  zflu, zflv, zdiv
       CHARACTER(lc)                     ::  charout                   ! local character
       REAL(wp), PARAMETER               ::  zrelax = 0.5_wp           ! relaxation constant for iterative procedure
@@ -65 +69 @@
       INTEGER , PARAMETER               ::  its    = 100              ! Maximum number of iteration
       !!-------------------------------------------------------------------
+      TYPE(arrayptr)   , ALLOCATABLE, DIMENSION(:) ::   pt2d_array, zrlx_array
+      CHARACTER(len=1) , ALLOCATABLE, DIMENSION(:) ::   type_array ! define the nature of ptab array grid-points
+      !                                                            ! = T , U , V , F , W and I points
+      REAL(wp)        , ALLOCATABLE, DIMENSION(:)  ::   psgn_array    ! =-1 the sign change across the north fold boundary
+
+     !!--------------------------------------------------------------------- 
+
+      !                       !==  Initialisation  ==!
+      ! +1 open water diffusion
+      isize = jpl*(ihdf_vars+nlay_i)+1
+      ALLOCATE( zconv (isize) )
+      ALLOCATE( pt2d_array(isize) , zrlx_array(isize) )
+      ALLOCATE( type_array(isize) )
+      ALLOCATE( psgn_array(isize) )
       
-      CALL wrk_alloc( jpi, jpj, zrlx, zflu, zflv, zdiv0, zdiv, ztab0 )
-
-      !                       !==  Initialisation  ==!
+      CALL wrk_alloc( jpi, jpj, isize, zrlx, zdiv0, ztab0 )
+      CALL wrk_alloc( jpi, jpj, zflu, zflv, zdiv )
+
+      DO jk= 1 , isize
+         pt2d_array(jk)%pt2d=>ptab(:,:,jk)
+         zrlx_array(jk)%pt2d=>zrlx(:,:,jk)
+         type_array(jk)='T'
+         psgn_array(jk)=1.
+      END DO
+
       !
       IF( linit ) THEN              ! Metric coefficient (compute at the first call and saved in efact)
@@ -74 +99 @@
          IF( lk_mpp    )   CALL mpp_sum( ierr )
          IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'lim_hdf : unable to allocate arrays' )
-         DO jj = 2, jpjm1  
+         DO jj = 2, jpjm1
             DO ji = fs_2 , fs_jpim1   ! vector opt.
                efact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj) + e1v(ji,jj) + e1v(ji,jj-1) ) * r1_e12t(ji,jj)
@@ -83 +108 @@
       !                             ! Time integration parameters
       !
-      ztab0(:, : ) = ptab(:,:)      ! Arrays initialization
-      zdiv0(:, 1 ) = 0._wp
-      zdiv0(:,jpj) = 0._wp
-      zflu (jpi,:) = 0._wp   
-      zflv (jpi,:) = 0._wp
-      zdiv0(1,  :) = 0._wp
-      zdiv0(jpi,:) = 0._wp
+      zflu (jpi,: ) = 0._wp
+      zflv (jpi,: ) = 0._wp
+
+      DO jk=1 , isize
+         ztab0(:, : , jk ) = ptab(:,:,jk)      ! Arrays initialization
+         zdiv0(:, 1 , jk ) = 0._wp
+         zdiv0(:,jpj, jk ) = 0._wp
+         zdiv0(1,  :, jk ) = 0._wp
+         zdiv0(jpi,:, jk ) = 0._wp
+      END DO
 
       zconv = 1._wp           !==  horizontal diffusion using a Crant-Nicholson scheme  ==!
       iter  = 0
       !
-      DO WHILE( zconv > ( 2._wp * 1.e-04 ) .AND. iter <= its )   ! Sub-time step loop
+      DO WHILE( MAXVAL(zconv(:)) > ( 2._wp * 1.e-04 ) .AND. iter <= its )   ! Sub-time step loop
          !
          iter = iter + 1                                 ! incrementation of the sub-time step number
          !
+         DO jk = 1 , isize
+            jl = (jk-1) /( ihdf_vars+nlay_i)+1
+            IF (zconv(jk) > ( 2._wp * 1.e-04 )) THEN
+               DO jj = 1, jpjm1                                ! diffusive fluxes in U- and V- direction
+                  DO ji = 1 , fs_jpim1   ! vector opt.
+                     zflu(ji,jj) = pahu3D(ji,jj,jl) * e2u(ji,jj) * r1_e1u(ji,jj) * ( ptab(ji+1,jj,jk) - ptab(ji,jj,jk) )
+                     zflv(ji,jj) = pahv3D(ji,jj,jl) * e1v(ji,jj) * r1_e2v(ji,jj) * ( ptab(ji,jj+1,jk) - ptab(ji,jj,jk) )
+                  END DO
+               END DO
+               !
+               DO jj= 2, jpjm1                                 ! diffusive trend : divergence of the fluxes
+                  DO ji = fs_2 , fs_jpim1   ! vector opt. 
+                     zdiv(ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj) + zflv(ji,jj) - zflv(ji,jj-1) ) * r1_e12t(ji,jj)
+                  END DO
+               END DO
+               !
+               IF( iter == 1 )   zdiv0(:,:,jk) = zdiv(:,:)        ! save the 1st evaluation of the diffusive trend in zdiv0
+               !
+               DO jj = 2, jpjm1                                ! iterative evaluation
+                  DO ji = fs_2 , fs_jpim1   ! vector opt.
+                     zrlxint = (   ztab0(ji,jj,jk)    &
+                        &       +  rdt_ice * (           zalfa   * ( zdiv(ji,jj) + efact(ji,jj) * ptab(ji,jj,jk) )   &
+                        &                      + ( 1.0 - zalfa ) *   zdiv0(ji,jj,jk) )                               &
+                        &      ) / ( 1.0 + zalfa * rdt_ice * efact(ji,jj) )
+                     zrlx(ji,jj,jk) = ptab(ji,jj,jk) + zrelax * ( zrlxint - ptab(ji,jj,jk) )
+                  END DO
+               END DO
+            END IF
+
+         END DO
+
+         CALL lbc_lnk_multi( zrlx_array, type_array , psgn_array , isize ) ! Multiple interchange of all the variables
+         !
+         
+         IF ( MOD( iter-1 , nn_convfrq ) == 0 )  THEN   !Convergence test every nn_convfrq iterations (perf. optimization ) 
+            DO jk=1,isize
+               zconv(jk) = 0._wp                                   ! convergence test
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1
+                     zconv(jk) = MAX( zconv(jk), ABS( zrlx(ji,jj,jk) - ptab(ji,jj,jk) )  )
+                  END DO
+               END DO
+            END DO
+            IF( lk_mpp ) CALL mpp_max_multiple( zconv , isize )            ! max over the global domain for all the variables
+         ENDIF
+         !
+         DO jk=1,isize
+            ptab(:,:,jk) = zrlx(:,:,jk)
+         END DO
+         !
+      END DO                                       ! end of sub-time step loop
+
+     ! -----------------------
+      !!! final step (clem) !!!
+      DO jk = 1, isize
+         jl = (jk-1) /( ihdf_vars+nlay_i)+1
          DO jj = 1, jpjm1                                ! diffusive fluxes in U- and V- direction
             DO ji = 1 , fs_jpim1   ! vector opt.
-               zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) * r1_e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) )
-               zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) * r1_e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) )
+               zflu(ji,jj) = pahu3D(ji,jj,jl) * e2u(ji,jj) * r1_e1u(ji,jj) * ( ptab(ji+1,jj,jk) - ptab(ji,jj,jk) )
+               zflv(ji,jj) = pahv3D(ji,jj,jl) * e1v(ji,jj) * r1_e2v(ji,jj) * ( ptab(ji,jj+1,jk) - ptab(ji,jj,jk) )
             END DO
          END DO
@@ -108 +192 @@
             DO ji = fs_2 , fs_jpim1   ! vector opt. 
                zdiv(ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj) + zflv(ji,jj) - zflv(ji,jj-1) ) * r1_e12t(ji,jj)
-            END DO
-         END DO
-         !
-         IF( iter == 1 )   zdiv0(:,:) = zdiv(:,:)        ! save the 1st evaluation of the diffusive trend in zdiv0
-         !
-         DO jj = 2, jpjm1                                ! iterative evaluation
-            DO ji = fs_2 , fs_jpim1   ! vector opt.
-               zrlxint = (   ztab0(ji,jj)    &
-                  &       +  rdt_ice * (           zalfa   * ( zdiv(ji,jj) + efact(ji,jj) * ptab(ji,jj) )   &
-                  &                      + ( 1.0 - zalfa ) *   zdiv0(ji,jj) )                               & 
-                  &      ) / ( 1.0 + zalfa * rdt_ice * efact(ji,jj) )
-               zrlx(ji,jj) = ptab(ji,jj) + zrelax * ( zrlxint - ptab(ji,jj) )
-            END DO
-         END DO
-         CALL lbc_lnk( zrlx, 'T', 1. )                   ! lateral boundary condition
-         !
-         IF ( MOD( iter, nn_convfrq ) == 0 )  THEN    ! convergence test every nn_convfrq iterations (perf. optimization)
-            zconv = 0._wp
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                  zconv = MAX( zconv, ABS( zrlx(ji,jj) - ptab(ji,jj) )  )
-               END DO
-            END DO
-            IF( lk_mpp )   CALL mpp_max( zconv )      ! max over the global domain
-         ENDIF
-         !
-         ptab(:,:) = zrlx(:,:)
-         !
-      END DO                                       ! end of sub-time step loop
-
-      ! -----------------------
-      !!! final step (clem) !!!
-      DO jj = 1, jpjm1                                ! diffusive fluxes in U- and V- direction
-         DO ji = 1 , fs_jpim1   ! vector opt.
-            zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) * r1_e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) )
-            zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) * r1_e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) )
+               ptab(ji,jj,jk) = ztab0(ji,jj,jk) + 0.5 * ( zdiv(ji,jj) + zdiv0(ji,jj,jk) )
+            END DO
          END DO
       END DO
-      !
-      DO jj= 2, jpjm1                                 ! diffusive trend : divergence of the fluxes
-         DO ji = fs_2 , fs_jpim1   ! vector opt. 
-            zdiv(ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj) + zflv(ji,jj) - zflv(ji,jj-1) ) * r1_e12t(ji,jj)
-            ptab(ji,jj) = ztab0(ji,jj) + 0.5 * ( zdiv(ji,jj) + zdiv0(ji,jj) )
-         END DO
-      END DO
-      CALL lbc_lnk( ptab, 'T', 1. )                   ! lateral boundary condition
+
+      CALL lbc_lnk_multi( pt2d_array, type_array , psgn_array , isize ) ! Multiple interchange of all the variables
+
       !!! final step (clem) !!!
       ! -----------------------
 
       IF(ln_ctl)   THEN
-         zrlx(:,:) = ptab(:,:) - ztab0(:,:)
-         WRITE(charout,FMT="(' lim_hdf  : zconv =',D23.16, ' iter =',I4,2X)") zconv, iter
-         CALL prt_ctl( tab2d_1=zrlx, clinfo1=charout )
-      ENDIF
-      !
-      CALL wrk_dealloc( jpi, jpj, zrlx, zflu, zflv, zdiv0, zdiv, ztab0 )
+         DO jk = 1 , isize
+            zrlx(:,:,jk) = ptab(:,:,jk) - ztab0(:,:,jk)
+            WRITE(charout,FMT="(' lim_hdf  : zconv =',D23.16, ' iter =',I4,2X)") zconv, iter
+            CALL prt_ctl( tab2d_1=zrlx(:,:,jk), clinfo1=charout )
+         END DO
+      ENDIF
+      !
+      CALL wrk_dealloc( jpi, jpj, isize, zrlx, zdiv0, ztab0 )
+      CALL wrk_dealloc( jpi, jpj, zflu, zflv, zdiv )
+
+      DEALLOCATE( zconv )
+      DEALLOCATE( pt2d_array , zrlx_array )
+      DEALLOCATE( type_array )
+      DEALLOCATE( psgn_array )
       !
    END SUBROUTINE lim_hdf
+
 
    
@@ -179 +233 @@
       !!-------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namicehdf/ nn_convfrq
+      NAMELIST/namicehdf/ nn_convfrq 
       !!-------------------------------------------------------------------
       !
@@ -212 +266 @@
    !!======================================================================
 END MODULE limhdf
+

branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90

@@ -24 +24 @@
    USE par_oce          ! ocean parameters
    USE dom_ice          ! sea-ice domain
+   USE limvar           ! lim_var_salprof
    USE in_out_manager   ! I/O manager
    USE lib_mpp          ! MPP library
@@ -246 +247 @@
                ztest_1 = 1
             ELSE 
-              ! this write is useful
-              IF(lwp)  WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(i_hemis) 
                ztest_1 = 0
             ENDIF
@@ -258 +257 @@
                ztest_2 = 1
             ELSE
-              ! this write is useful
-              IF(lwp)  WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, &
-                            ' zvt_i_ini = ', zvt_i_ini(i_hemis)
                ztest_2 = 0
             ENDIF
@@ -268 +264 @@
                ztest_3 = 1
             ELSE
-               ! this write is useful
-               IF(lwp) WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zh_i_ini(i_fill,i_hemis) = ', &
-               zh_i_ini(i_fill,i_hemis), ' hi_max(jpl-1) = ', hi_max(i_fill-1)
                ztest_3 = 0
             ENDIF
@@ -278 +271 @@
             DO jl = 1, jpl
                IF ( za_i_ini(jl,i_hemis) .LT. 0._wp ) THEN 
-                  ! this write is useful
-                  IF(lwp) WRITE(numout,*) ' * TEST 4 POSITIVITY NOT OK FOR CAT ', jl, ' WITH A = ', za_i_ini(jl,i_hemis)
                   ztest_4 = 0
                ENDIF
@@ -337 +328 @@
          END DO
       END DO
+
+      ! for constant salinity in time
+      IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
+         CALL lim_var_salprof
+         smv_i = sm_i * v_i
+      ENDIF
 
       ! Snow temperature and heat content

branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90

@@ -651 +651 @@
             wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ij) * rhoic * r1_rdtice   ! increase in ice volume due to seawater frozen in voids
             
+             ! virtual salt flux to keep salinity constant
+            IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
+               srdg2(ij)      = srdg2(ij) - vsw(ij) * ( sss_m(ji,jj) - sm_i(ji,jj,jl1) )           ! ridge salinity = sm_i
+               sfx_bri(ji,jj) = sfx_bri(ji,jj) + sss_m(ji,jj)    * vsw(ij) * rhoic * r1_rdtice  &  ! put back sss_m into the ocean
+                  &                            - sm_i(ji,jj,jl1) * vsw(ij) * rhoic * r1_rdtice     ! and get  sm_i  from the ocean 
+            ENDIF
+
             !------------------------------------------            
             ! 3.7 Put the snow somewhere in the ocean
@@ -664 +671 @@
             hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnowrdg )         & 
                &                                - esrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice        ! heat sink for ocean (<0, W.m-2)
-
+               
             !-----------------------------------------------------------------
             ! 3.8 Compute quantities used to apportion ice among categories
@@ -859 +866 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0)))
+               strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bvm_i(ji,jj),0.0)))
             END DO
          END DO

branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

@@ -10 +10 @@
    !!            3.4  !  2011-01  (A. Porter)  dynamical allocation 
    !!            3.5  !  2012-08  (R. Benshila)  AGRIF 
+   !!            3.6  !  2016-06  (C. Rousset) Rewriting (conserves energy)
    !!----------------------------------------------------------------------
 #if defined key_lim3 || (  defined key_lim2 && ! defined key_lim2_vp )
@@ -95 +96 @@
       !!                 coriolis terms of the momentum equation
       !!              3) Solve the momentum equation (iterative procedure)
-      !!              4) Prevent high velocities if the ice is thin
-      !!              5) Recompute invariants of the strain rate tensor
+      !!              4) Recompute invariants of the strain rate tensor
       !!                 which are inputs of the ITD, store stress
       !!                 for the next time step
-      !!              6) Control prints of residual (convergence)
+      !!              5) Control prints of residual (convergence)
       !!                 and charge ellipse.
       !!                 The user should make sure that the parameters
@@ -106 +106 @@
       !!                 e.g. in the Canadian Archipelago
       !!
+      !! ** Notes   : Boundary condition for ice is chosen no-slip 
+      !!              but can be adjusted with param rn_shlat
+      !!
       !! References : Hunke and Dukowicz, JPO97
       !!              Bouillon et al., Ocean Modelling 2009
@@ -115 +118 @@
       INTEGER ::   jter     ! local integers
       CHARACTER (len=50) ::   charout
-      REAL(wp) ::   zt11, zt12, zt21, zt22, ztagnx, ztagny, delta                         !
-      REAL(wp) ::   za, zstms          ! local scalars
-      REAL(wp) ::   zc1, zc2, zc3      ! ice mass
-
-      REAL(wp) ::   dtevp , z1_dtevp              ! time step for subcycling
-      REAL(wp) ::   dtotel, z1_dtotel, ecc2, ecci ! square of yield ellipse eccenticity
-      REAL(wp) ::   z0, zr, zcca, zccb            ! temporary scalars
-      REAL(wp) ::   zu_ice2, zv_ice1              !
-      REAL(wp) ::   zddc, zdtc                    ! delta on corners and on centre
-      REAL(wp) ::   zdst                          ! shear at the center of the grid point
-      REAL(wp) ::   zdsshx, zdsshy                ! term for the gradient of ocean surface
-      REAL(wp) ::   sigma1, sigma2                ! internal ice stress
-
-      REAL(wp) ::   zresm         ! Maximal error on ice velocity
-      REAL(wp) ::   zintb, zintn  ! dummy argument
-
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zpresh           ! temporary array for ice strength
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zpreshc          ! Ice strength on grid cell corners (zpreshc)
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zfrld1, zfrld2   ! lead fraction on U/V points
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zmass1, zmass2   ! ice/snow mass on U/V points
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zcorl1, zcorl2   ! coriolis parameter on U/V points
-      REAL(wp), POINTER, DIMENSION(:,:) ::   za1ct , za2ct    ! temporary arrays
-      REAL(wp), POINTER, DIMENSION(:,:) ::   v_oce1           ! ocean u/v component on U points                           
-      REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce2           ! ocean u/v component on V points
-      REAL(wp), POINTER, DIMENSION(:,:) ::   u_ice2, v_ice1   ! ice u/v component on V/U point
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zf1   , zf2      ! arrays for internal stresses
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zmask            ! mask ocean grid points
+
+      REAL(wp) ::   zdtevp, z1_dtevp                                         ! time step for subcycling
+      REAL(wp) ::   ecc2, z1_ecc2                                            ! square of yield ellipse eccenticity
+      REAL(wp) ::   zbeta, zalph1, z1_alph1, zalph2, z1_alph2                ! alpha and beta from Bouillon 2009 and 2013
+      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                            ! ice/snow mass
+      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2            ! temporary scalars
+      REAL(wp) ::   zTauO, zTauE, zCor                                       ! temporary scalars
+
+      REAL(wp) ::   zsig1, zsig2                                             ! internal ice stress
+      REAL(wp) ::   zresm                                                    ! Maximal error on ice velocity
+      REAL(wp) ::   zintb, zintn                                             ! dummy argument
       
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zdt              ! tension at centre of grid cells
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zds              ! Shear on northeast corner of grid cells
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zs1   , zs2      ! Diagonal stress tensor components zs1 and zs2 
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zs12             ! Non-diagonal stress tensor component zs12
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_ice, zv_ice, zresr   ! Local error on velocity
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zpice            ! array used for the calculation of ice surface slope:
-                                                              !   ocean surface (ssh_m) if ice is not embedded
-                                                              !   ice top surface if ice is embedded   
-
-      REAL(wp), PARAMETER               ::   zepsi = 1.0e-20_wp ! tolerance parameter
-      REAL(wp), PARAMETER               ::   zvmin = 1.0e-03_wp ! ice volume below which ice velocity equals ocean velocity
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zpresh                          ! temporary array for ice strength
+      REAL(wp), POINTER, DIMENSION(:,:) ::   z1_e1t0, z1_e2t0                ! scale factors
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zp_delt                         ! P/delta at T points
+      !
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zaU   , zaV                     ! ice fraction on U/V points
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zmU_t, zmV_t                    ! ice/snow mass/dt on U/V points
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zmf                             ! coriolis parameter at T points
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zspgU , zspgV                   ! surface pressure gradient at U/V points
+      REAL(wp), POINTER, DIMENSION(:,:) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                           
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zfU   , zfV                     ! internal stresses
+      
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zds                             ! shear
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zs1, zs2, zs12                  ! stress tensor components
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_ice, zv_ice, zresr           ! check convergence
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zpice                           ! array used for the calculation of ice surface slope:
+                                                                             !   ocean surface (ssh_m) if ice is not embedded
+                                                                             !   ice top surface if ice is embedded   
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zswitchU, zswitchV              ! dummy arrays
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zmaskU, zmaskV                  ! mask for ice presence
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zfmask, zwf                     ! mask at F points for the ice
+
+      REAL(wp), PARAMETER               ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
+      REAL(wp), PARAMETER               ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity
+      REAL(wp), PARAMETER               ::   zshlat = 2._wp                  ! boundary condition for sea-ice velocity (2=no slip ; 0=free slip)
       !!-------------------------------------------------------------------
 
-      CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
-      CALL wrk_alloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask               )
-      CALL wrk_alloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  )
-      CALL wrk_alloc( jpi,jpj, zs1   , zs2   , zs12   , zresr , zpice                 )
+      CALL wrk_alloc( jpi,jpj, zpresh, z1_e1t0, z1_e2t0, zp_delt )
+      CALL wrk_alloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia )
+      CALL wrk_alloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV )
+      CALL wrk_alloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
+      CALL wrk_alloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf )
 
 #if  defined key_lim2 && ! defined key_lim2_vp
@@ -176 +178 @@
       !
       !------------------------------------------------------------------------------!
-      ! 1) Ice strength (zpresh)                                !
-      !------------------------------------------------------------------------------!
-      !
-      ! Put every vector to 0
-      delta_i(:,:) = 0._wp   ;
-      zpresh (:,:) = 0._wp   ;  
-      zpreshc(:,:) = 0._wp
-      u_ice2 (:,:) = 0._wp   ;   v_ice1(:,:) = 0._wp
-      divu_i (:,:) = 0._wp   ;   zdt   (:,:) = 0._wp   ;   zds(:,:) = 0._wp
-      shear_i(:,:) = 0._wp
-
+      ! 0) mask at F points for the ice (on the whole domain, not only k_j1,k_jpj) 
+      !------------------------------------------------------------------------------!
+      ! ocean/land mask
+      DO jj = 1, jpjm1
+         DO ji = 1, jpim1      ! NO vector opt.
+            zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
+         END DO
+      END DO
+      CALL lbc_lnk( zfmask, 'F', 1._wp )
+
+      ! Lateral boundary conditions on velocity (modify zfmask)
+      zwf(:,:) = zfmask(:,:)
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            IF( zfmask(ji,jj) == 0._wp ) THEN
+               zfmask(ji,jj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) )
+            ENDIF
+         END DO
+      END DO
+      DO jj = 2, jpjm1
+         IF( zfmask(1,jj) == 0._wp ) THEN
+            zfmask(1  ,jj) = zshlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
+         ENDIF
+         IF( zfmask(jpi,jj) == 0._wp ) THEN
+            zfmask(jpi,jj) = zshlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
+         ENDIF
+      END DO
+      DO ji = 2, jpim1
+         IF( zfmask(ji,1) == 0._wp ) THEN
+            zfmask(ji,1  ) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
+         ENDIF
+         IF( zfmask(ji,jpj) == 0._wp ) THEN
+            zfmask(ji,jpj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
+         ENDIF
+      END DO
+      CALL lbc_lnk( zfmask, 'F', 1._wp )
+
+      !------------------------------------------------------------------------------!
+      ! 1) define some variables and initialize arrays
+      !------------------------------------------------------------------------------!
+      ! ecc2: square of yield ellipse eccenticrity
+      ecc2    = rn_ecc * rn_ecc
+      z1_ecc2 = 1._wp / ecc2
+
+      ! Time step for subcycling
+      zdtevp   = rdt_ice / REAL( nn_nevp )
+      z1_dtevp = 1._wp / zdtevp
+
+      ! alpha parameters (Bouillon 2009)
 #if defined key_lim3
-      CALL lim_itd_me_icestrength( nn_icestr )      ! LIM-3: Ice strength on T-points
-#endif
-
-      DO jj = k_j1 , k_jpj       ! Ice mass and temp variables
-         DO ji = 1 , jpi
+      zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
+#else
+      zalph1 = ( 2._wp * telast ) * z1_dtevp
+#endif
+      zalph2 = zalph1 * z1_ecc2
+
+      z1_alph1 = 1._wp / ( zalph1 + 1._wp )
+      z1_alph2 = 1._wp / ( zalph2 + 1._wp )
+
+      ! Initialise stress tensor 
+      zs1 (:,:) = stress1_i (:,:) 
+      zs2 (:,:) = stress2_i (:,:)
+      zs12(:,:) = stress12_i(:,:)
+
+      ! Ice strength
 #if defined key_lim3
-            zpresh(ji,jj) = tmask(ji,jj,1) *  strength(ji,jj)
-#endif
-#if defined key_lim2
-            zpresh(ji,jj) = tmask(ji,jj,1) *  pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) )
-#endif
-            ! zmask = 1 where there is ice or on land
-            zmask(ji,jj)    = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - zepsi ) ) ) * tmask(ji,jj,1)
+      CALL lim_itd_me_icestrength( nn_icestr )
+      zpresh(:,:) = tmask(:,:,1) *  strength(:,:)
+#else
+      zpresh(:,:) = tmask(:,:,1) *  pstar * vt_i(:,:) * EXP( -c_rhg * (1. - at_i(:,:) ) )
+#endif
+
+      ! scale factors
+      DO jj = k_j1+1, k_jpj-1
+         DO ji = fs_2, fs_jpim1
+            z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj  ) + e1t(ji,jj  ) )
+            z1_e2t0(ji,jj) = 1._wp / ( e2t(ji  ,jj+1) + e2t(ji,jj  ) )
          END DO
       END DO
-
-      ! Ice strength on grid cell corners (zpreshc)
-      ! needed for calculation of shear stress 
-      DO jj = k_j1+1, k_jpj-1
-         DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1)
-            zstms          =  tmask(ji+1,jj+1,1) * wght(ji+1,jj+1,2,2) + tmask(ji,jj+1,1) * wght(ji+1,jj+1,1,2) +   &
-               &              tmask(ji+1,jj,1)   * wght(ji+1,jj+1,2,1) + tmask(ji,jj,1)   * wght(ji+1,jj+1,1,1)
-            zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) +   &
-               &               zpresh(ji+1,jj)   * wght(ji+1,jj+1,2,1) + zpresh(ji,jj)   * wght(ji+1,jj+1,1,1)     &
-               &             ) / MAX( zstms, zepsi )
-         END DO
-      END DO
-      CALL lbc_lnk( zpreshc(:,:), 'F', 1. )
+            
       !
       !------------------------------------------------------------------------------!
       ! 2) Wind / ocean stress, mass terms, coriolis terms
       !------------------------------------------------------------------------------!
-      !
-      !  Wind stress, coriolis and mass terms on the sides of the squares        
-      !  zfrld1: lead fraction on U-points                                      
-      !  zfrld2: lead fraction on V-points                                     
-      !  zmass1: ice/snow mass on U-points                                    
-      !  zmass2: ice/snow mass on V-points                                   
-      !  zcorl1: Coriolis parameter on U-points                             
-      !  zcorl2: Coriolis parameter on V-points                            
-      !  (ztagnx,ztagny): wind stress on U/V points                       
-      !  v_oce1: ocean v component on u points                          
-      !  u_oce2: ocean u component on v points                         
 
       IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==!
@@ -242 +273 @@
          zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
          !
-         zpice(:,:) = ssh_m(:,:) + (  zintn * snwice_mass(:,:) +  zintb * snwice_mass_b(:,:)  ) * r1_rau0
+         zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0
          !
       ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==!
@@ -251 +282 @@
          DO ji = fs_2, fs_jpim1
 
-            zc1 = tmask(ji  ,jj  ,1) * ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) )
-            zc2 = tmask(ji+1,jj  ,1) * ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) )
-            zc3 = tmask(ji  ,jj+1,1) * ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) )
-
-            zt11 = tmask(ji  ,jj,1) * e1t(ji  ,jj)
-            zt12 = tmask(ji+1,jj,1) * e1t(ji+1,jj)
-            zt21 = tmask(ji,jj  ,1) * e2t(ji,jj  )
-            zt22 = tmask(ji,jj+1,1) * e2t(ji,jj+1)
-
-            ! Leads area.
-            zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + zepsi )
-            zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + zepsi )
-
-            ! Mass, coriolis coeff. and currents
-            zmass1(ji,jj) = ( zt12 * zc1 + zt11 * zc2 ) / ( zt11 + zt12 + zepsi )
-            zmass2(ji,jj) = ( zt22 * zc1 + zt21 * zc3 ) / ( zt21 + zt22 + zepsi )
-            zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj) * fcor(ji,jj) + e1t(ji,jj) * fcor(ji+1,jj) )   &
-               &                          / ( e1t(ji,jj) + e1t(ji+1,jj) + zepsi )
-            zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1) * fcor(ji,jj) + e2t(ji,jj) * fcor(ji,jj+1) )   &
-               &                          / ( e2t(ji,jj+1) + e2t(ji,jj) + zepsi )
-            !
-            ! Ocean has no slip boundary condition
-            v_oce1(ji,jj)  = 0.5 * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji,jj)      &
-               &                   + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji+1,jj) )  &
-               &                   / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)  
-
-            u_oce2(ji,jj)  = 0.5 * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj)      &
-               &                   + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj+1) )  &
-               &                   / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
-
-            ! Wind stress at U,V-point
-            ztagnx = ( 1. - zfrld1(ji,jj) ) * utau_ice(ji,jj)
-            ztagny = ( 1. - zfrld2(ji,jj) ) * vtau_ice(ji,jj)
-
-            ! Computation of the velocity field taking into account the ice internal interaction.
-            ! Terms that are independent of the velocity field.
-
-            ! SB On utilise maintenant le gradient de la pente de l'ocean
-            ! include it later
-
-            zdsshx =  ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
-            zdsshy =  ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
-
-            za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx
-            za2ct(ji,jj) = ztagny - zmass2(ji,jj) * grav * zdsshy
+            ! ice fraction at U-V points
+            zaU(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e12t(ji,jj) + at_i(ji+1,jj) * e12t(ji+1,jj) ) * r1_e12u(ji,jj) * umask(ji,jj,1)
+            zaV(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e12t(ji,jj) + at_i(ji,jj+1) * e12t(ji,jj+1) ) * r1_e12v(ji,jj) * vmask(ji,jj,1)
+
+            ! Ice/snow mass at U-V points
+            zm1 = ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) )
+            zm2 = ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) )
+            zm3 = ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) )
+            zmassU = 0.5_wp * ( zm1 * e12t(ji,jj) + zm2 * e12t(ji+1,jj) ) * r1_e12u(ji,jj) * umask(ji,jj,1)
+            zmassV = 0.5_wp * ( zm1 * e12t(ji,jj) + zm3 * e12t(ji,jj+1) ) * r1_e12v(ji,jj) * vmask(ji,jj,1)
+
+            ! Ocean currents at U-V points
+            v_oceU(ji,jj)   = 0.5_wp * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji+1,jj)    &
+               &                       + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
+            
+            u_oceV(ji,jj)   = 0.5_wp * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj+1)    &
+               &                       + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
+
+            ! Coriolis at T points (m*f)
+            zmf(ji,jj)      = zm1 * fcor(ji,jj)
+
+            ! m/dt
+            zmU_t(ji,jj)    = zmassU * z1_dtevp
+            zmV_t(ji,jj)    = zmassV * z1_dtevp
+
+            ! Drag ice-atm.
+            zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
+            zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
+
+            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
+            zspgU(ji,jj)    = - zmassU * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
+            zspgV(ji,jj)    = - zmassV * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
+
+            ! masks
+            zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
+            zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
+
+            ! switches
+            zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin
+            zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin
 
          END DO
       END DO
-
+      CALL lbc_lnk( zmf, 'T', 1. )
       !
       !------------------------------------------------------------------------------!
@@ -305 +331 @@
       !------------------------------------------------------------------------------!
       !
-      ! Time step for subcycling
-      dtevp  = rdt_ice / nn_nevp
-#if defined key_lim3
-      dtotel = dtevp / ( 2._wp * rn_relast * rdt_ice )
-#else
-      dtotel = dtevp / ( 2._wp * telast )
-#endif
-      z1_dtotel = 1._wp / ( 1._wp + dtotel )
-      z1_dtevp  = 1._wp / dtevp
-      !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter)
-      ecc2 = rn_ecc * rn_ecc
-      ecci = 1. / ecc2
-
-      !-Initialise stress tensor 
-      zs1 (:,:) = stress1_i (:,:) 
-      zs2 (:,:) = stress2_i (:,:)
-      zs12(:,:) = stress12_i(:,:)
-
       !                                               !----------------------!
       DO jter = 1 , nn_nevp                           !    loop over jter    !
          !                                            !----------------------!        
-         DO jj = k_j1, k_jpj-1
-            zu_ice(:,jj) = u_ice(:,jj)    ! velocity at previous time step
-            zv_ice(:,jj) = v_ice(:,jj)
-         END DO
-
-         DO jj = k_j1+1, k_jpj-1
-            DO ji = fs_2, fs_jpim1   !RB bug no vect opt due to zmask
-
-               !  
-               !- Divergence, tension and shear (Section a. Appendix B of Hunke & Dukowicz, 2002)
-               !- divu_i(:,:), zdt(:,:): divergence and tension at centre of grid cells
-               !- zds(:,:): shear on northeast corner of grid cells
-               !
-               !- IMPORTANT REMINDER: Dear Gurvan, note that, the way these terms are coded, 
-               !                      there are many repeated calculations. 
-               !                      Speed could be improved by regrouping terms. For
-               !                      the moment, however, the stress is on clarity of coding to avoid
-               !                      bugs (Martin, for Miguel).
-               !
-               !- ALSO: arrays zdt, zds and delta could 
-               !  be removed in the future to minimise memory demand.
-               !
-               !- MORE NOTES: Note that we are calculating deformation rates and stresses on the corners of
-               !              grid cells, exactly as in the B grid case. For simplicity, the indexation on
-               !              the corners is the same as in the B grid.
-               !
-               !
-               divu_i(ji,jj) = (  e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
-                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
-                  &            ) * r1_e12t(ji,jj)
-
-               zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
-                  &         - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
-                  &         ) * r1_e12t(ji,jj)
-
-               !
+         IF(ln_ctl) THEN   ! Convergence test
+            DO jj = k_j1, k_jpj-1
+               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
+               zv_ice(:,jj) = v_ice(:,jj)
+            END DO
+         ENDIF
+
+         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
+         DO jj = k_j1, k_jpj-1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points
+            DO ji = 1, jpim1
+
+               ! shear at F points
                zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
                   &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
-                  &         ) * r1_e12f(ji,jj) * ( 2._wp - fmask(ji,jj,1) )   &
-                  &         * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1)
-
-
-               v_ice1(ji,jj)  = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
-                  &                      + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) )   &
-                  &                      / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 
-
-               u_ice2(ji,jj)  = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
-                  &                      + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) )   &
-                  &                      / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
-            END DO
-         END DO
-
-         CALL lbc_lnk_multi( v_ice1, 'U', -1., u_ice2, 'V', -1. )      ! lateral boundary cond.
-         
+                  &         ) * r1_e12f(ji,jj) * zfmask(ji,jj)
+
+            END DO
+         END DO
+         CALL lbc_lnk( zds, 'F', 1. )
+
          DO jj = k_j1+1, k_jpj-1
-            DO ji = fs_2, fs_jpim1
-
-               !- Calculate Delta at centre of grid cells
-               zdst          = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )   &
-                  &            + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1)   &
-                  &            ) * r1_e12t(ji,jj)
-
-               delta          = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 )  
-               delta_i(ji,jj) = delta + rn_creepl
-
-               !- Calculate Delta on corners
-               zddc  = (  ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  &
-                  &     + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  &
-                  &    ) * r1_e12f(ji,jj)
-
-               zdtc  = (- ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  &
-                  &     + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  &
-                  &    ) * r1_e12f(ji,jj)
-
-               zddc = SQRT( zddc**2 + ( zdtc**2 + zds(ji,jj)**2 ) * usecc2 ) + rn_creepl
-
-               !-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 * ( divu_i(ji,jj) - delta ) / delta_i(ji,jj)   * zpresh(ji,jj)   &
-                  &          ) * z1_dtotel
-               zs2(ji,jj)  = ( zs2 (ji,jj) + dtotel *         ecci * zdt(ji,jj) / delta_i(ji,jj)   * zpresh(ji,jj)   &
-                  &          ) * z1_dtotel
-               !-Calculate stress tensor component zs12 at corners
-               zs12(ji,jj) = ( zs12(ji,jj) + dtotel *         ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj)  &
-                  &          ) * z1_dtotel 
-
-            END DO
-         END DO
-
-         CALL lbc_lnk_multi( zs1 , 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
+            DO ji = 2, jpim1 ! no vector loop
+
+               ! shear**2 at T points (doc eq. A16)
+               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e12f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e12f(ji-1,jj  )  &
+                  &   + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1)  &
+                  &   ) * 0.25_wp * r1_e12t(ji,jj)
+              
+               ! divergence at T points
+               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
+                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
+                  &    ) * r1_e12t(ji,jj)
+               zdiv2 = zdiv * zdiv
+               
+               ! tension at T points
+               zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
+                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
+                  &   ) * r1_e12t(ji,jj)
+               zdt2 = zdt * zdt
+               
+               ! delta at T points
+               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * usecc2 )  
+
+               ! P/delta at T points
+               zp_delt(ji,jj) = zpresh(ji,jj) / ( zdelta + rn_creepl )
+               
+               ! stress at T points
+               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1
+               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2
+             
+            END DO
+         END DO
+         CALL lbc_lnk( zp_delt, 'T', 1. )
+
+         DO jj = k_j1, k_jpj-1
+            DO ji = 1, jpim1
+
+               ! P/delta at F points
+               zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) )
+               
+               ! stress at F points
+               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2
+
+            END DO
+         END DO
+         CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
  
-         ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002)
+         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
          DO jj = k_j1+1, k_jpj-1
-            DO ji = fs_2, fs_jpim1
-               !- contribution of zs1, zs2 and zs12 to zf1
-               zf1(ji,jj) = 0.5 * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)  &
-                  &             + ( zs2(ji+1,jj) * e2t(ji+1,jj)**2 - zs2(ji,jj) * e2t(ji,jj)**2 ) * r1_e2u(ji,jj)          &
-                  &             + 2.0 * ( zs12(ji,jj) * e1f(ji,jj)**2 - zs12(ji,jj-1) * e1f(ji,jj-1)**2 ) * r1_e1u(ji,jj)  &
-                  &                ) * r1_e12u(ji,jj)
-               ! contribution of zs1, zs2 and zs12 to zf2
-               zf2(ji,jj) = 0.5 * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)  &
-                  &             - ( zs2(ji,jj+1) * e1t(ji,jj+1)**2 - zs2(ji,jj) * e1t(ji,jj)**2 ) * r1_e1v(ji,jj)          &
-                  &             + 2.0 * ( zs12(ji,jj) * e2f(ji,jj)**2 - zs12(ji-1,jj) * e2f(ji-1,jj)**2 ) * r1_e2v(ji,jj)  &
-                  &               )  * r1_e12v(ji,jj)
+            DO ji = fs_2, fs_jpim1               
+
+               ! U points
+               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
+                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
+                  &                    ) * r1_e2u(ji,jj)                                                                      &
+                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
+                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
+                  &                  ) * r1_e12u(ji,jj)
+
+               ! V points
+               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
+                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
+                  &                    ) * r1_e1v(ji,jj)                                                                      &
+                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
+                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
+                  &                  ) * r1_e12v(ji,jj)
+
+               ! u_ice at V point
+               u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
+                  &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
+               
+               ! v_ice at U point
+               v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
+                  &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
+
             END DO
          END DO
          !
-         ! Computation of ice velocity
-         !
-         ! Both the Coriolis term and the ice-ocean drag are solved semi-implicitly.
-         !
-         IF (MOD(jter,2).eq.0) THEN 
-
+         ! --- Computation of ice velocity --- !
+         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp
+         !  Bouillon et al. 2009 (eq 34-35) => stable
+         IF( MOD(jter,2) .EQ. 0 ) THEN ! even iterations
+            
             DO jj = k_j1+1, k_jpj-1
                DO ji = fs_2, fs_jpim1
-                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1)
-                  z0           = zmass1(ji,jj) * z1_dtevp
-
-                  ! SB modif because ocean has no slip boundary condition
-                  zv_ice1      = 0.5 * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji  ,jj)     &
-                     &                 + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) )   &
-                     &                 / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)
-                  za           = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 +  &
-                     &                         ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) )
-                  zr           = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj)
-                  zcca         = z0 + za
-                  zccb         = zcorl1(ji,jj)
-                  u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 
+
+                  ! tau_io/(v_oce - v_ice)
+                  zTauO = zaV(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
+                     &                             + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
+
+                  ! Coriolis at V-points (energy conserving formulation)
+                  zCor  = - 0.25_wp * r1_e2v(ji,jj) *  &
+                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
+                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
+
+                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
+                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
+                  
+                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  v_ice(ji,jj) = ( ( zmV_t(ji,jj) * v_ice(ji,jj) + zTauE + zTauO * v_ice(ji,jj)  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &             ) / MAX( zepsi, zmV_t(ji,jj) + zTauO ) * zswitchV(ji,jj)      &  ! m/dt + tau_io(only ice part)
+                     &             + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
+                     &           ) * zmaskV(ji,jj)
                END DO
             END DO
-
-            CALL lbc_lnk( u_ice(:,:), 'U', -1. )
+            CALL lbc_lnk( v_ice, 'V', -1. )
+            
+#if defined key_agrif && defined key_lim2
+            CALL agrif_rhg_lim2( jter, nn_nevp, 'V' )
+#endif
+#if defined key_bdy
+            CALL bdy_ice_lim_dyn( 'V' )
+#endif         
+
+            DO jj = k_j1+1, k_jpj-1
+               DO ji = fs_2, fs_jpim1
+                               
+                  ! tau_io/(u_oce - u_ice)
+                  zTauO = zaU(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
+                     &                             + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
+
+                  ! Coriolis at U-points (energy conserving formulation)
+                  zCor  =   0.25_wp * r1_e1u(ji,jj) *  &
+                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
+                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
+                  
+                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
+                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
+
+                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  u_ice(ji,jj) = ( ( zmU_t(ji,jj) * u_ice(ji,jj) + zTauE + zTauO * u_ice(ji,jj)  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &             ) / MAX( zepsi, zmU_t(ji,jj) + zTauO ) * zswitchU(ji,jj)      &  ! m/dt + tau_io(only ice part)
+                     &             + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
+                     &           ) * zmaskU(ji,jj)
+               END DO
+            END DO
+            CALL lbc_lnk( u_ice, 'U', -1. )
+            
 #if defined key_agrif && defined key_lim2
             CALL agrif_rhg_lim2( jter, nn_nevp, 'U' )
 #endif
 #if defined key_bdy
-         CALL bdy_ice_lim_dyn( 'U' )
+            CALL bdy_ice_lim_dyn( 'U' )
 #endif         
+
+         ELSE ! odd iterations
 
             DO jj = k_j1+1, k_jpj-1
                DO ji = fs_2, fs_jpim1
-
-                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1)
-                  z0           = zmass2(ji,jj) * z1_dtevp
-                  ! SB modif because ocean has no slip boundary condition
-                  zu_ice2      = 0.5 * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj)       &
-                     &                 + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) )   &
-                     &                 / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
-                  za           = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 +  & 
-                     &                         ( v_ice(ji,jj) - v_oce(ji,jj))**2 ) * ( 1.0 - zfrld2(ji,jj) )
-                  zr           = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj)
-                  zcca         = z0 + za
-                  zccb         = zcorl2(ji,jj)
-                  v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch
+                               
+                  ! tau_io/(u_oce - u_ice)
+                  zTauO = zaU(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
+                     &                             + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
+
+                  ! Coriolis at U-points (energy conserving formulation)
+                  zCor  =   0.25_wp * r1_e1u(ji,jj) *  &
+                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
+                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
+                  
+                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
+                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
+
+                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  u_ice(ji,jj) = ( ( zmU_t(ji,jj) * u_ice(ji,jj) + zTauE + zTauO * u_ice(ji,jj)  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &             ) / MAX( zepsi, zmU_t(ji,jj) + zTauO ) * zswitchU(ji,jj)      &  ! m/dt + tau_io(only ice part)
+                     &             + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
+                     &           ) * zmaskU(ji,jj)
                END DO
             END DO
-
-            CALL lbc_lnk( v_ice(:,:), 'V', -1. )
+            CALL lbc_lnk( u_ice, 'U', -1. )
+            
+#if defined key_agrif && defined key_lim2
+            CALL agrif_rhg_lim2( jter, nn_nevp, 'U' )
+#endif
+#if defined key_bdy
+            CALL bdy_ice_lim_dyn( 'U' )
+#endif         
+
+           DO jj = k_j1+1, k_jpj-1
+               DO ji = fs_2, fs_jpim1
+
+                  ! tau_io/(v_oce - v_ice)
+                  zTauO = zaV(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
+                     &                             + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
+
+                  ! Coriolis at V-points (energy conserving formulation)
+                  zCor  = - 0.25_wp * r1_e2v(ji,jj) *  &
+                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
+                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
+
+                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
+                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
+                  
+                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  v_ice(ji,jj) = ( ( zmV_t(ji,jj) * v_ice(ji,jj) + zTauE + zTauO * v_ice(ji,jj)  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &             ) / MAX( zepsi, zmV_t(ji,jj) + zTauO ) * zswitchV(ji,jj)      &  ! m/dt + tau_io(only ice part)
+                     &             + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
+                     &           ) * zmaskV(ji,jj)
+               END DO
+            END DO
+            CALL lbc_lnk( v_ice, 'V', -1. )
+            
 #if defined key_agrif && defined key_lim2
             CALL agrif_rhg_lim2( jter, nn_nevp, 'V' )
 #endif
 #if defined key_bdy
-         CALL bdy_ice_lim_dyn( 'V' )
+            CALL bdy_ice_lim_dyn( 'V' )
 #endif         
 
-         ELSE 
+         ENDIF
+         
+         IF(ln_ctl) THEN   ! Convergence test
             DO jj = k_j1+1, k_jpj-1
-               DO ji = fs_2, fs_jpim1
-                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1)
-                  z0           = zmass2(ji,jj) * z1_dtevp
-                  ! SB modif because ocean has no slip boundary condition
-                  zu_ice2      = 0.5 * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj)       &
-                     &                  +( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) )   &
-                     &                 / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)   
-
-                  za           = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 +  &
-                     &                         ( v_ice(ji,jj) - v_oce(ji,jj) )**2 ) * ( 1.0 - zfrld2(ji,jj) )
-                  zr           = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj)
-                  zcca         = z0 + za
-                  zccb         = zcorl2(ji,jj)
-                  v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch
-               END DO
-            END DO
-
-            CALL lbc_lnk( v_ice(:,:), 'V', -1. )
-#if defined key_agrif && defined key_lim2
-            CALL agrif_rhg_lim2( jter, nn_nevp, 'V' )
-#endif
-#if defined key_bdy
-         CALL bdy_ice_lim_dyn( 'V' )
-#endif         
-
-            DO jj = k_j1+1, k_jpj-1
-               DO ji = fs_2, fs_jpim1
-                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1)
-                  z0           = zmass1(ji,jj) * z1_dtevp
-                  zv_ice1      = 0.5 * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji,jj)       &
-                     &                 + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) )   &
-                     &                 / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)
-
-                  za           = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 +  &
-                     &                         ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) )
-                  zr           = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj)
-                  zcca         = z0 + za
-                  zccb         = zcorl1(ji,jj)
-                  u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 
-               END DO
-            END DO
-
-            CALL lbc_lnk( u_ice(:,:), 'U', -1. )
-#if defined key_agrif && defined key_lim2
-            CALL agrif_rhg_lim2( jter, nn_nevp, 'U' )
-#endif
-#if defined key_bdy
-         CALL bdy_ice_lim_dyn( 'U' )
-#endif         
-
-         ENDIF
-         
-         IF(ln_ctl) THEN
-            !---  Convergence test.
-            DO jj = k_j1+1 , k_jpj-1
                zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
             END DO
@@ -552 +572 @@
             IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
          ENDIF
-
+         !
          !                                                ! ==================== !
       END DO                                              !  end loop over jter  !
@@ -558 +578 @@
       !
       !------------------------------------------------------------------------------!
-      ! 4) Prevent ice velocities when the ice is thin
-      !------------------------------------------------------------------------------!
-      ! If the ice volume is below zvmin then ice velocity should equal the
-      ! ocean velocity. This prevents high velocity when ice is thin
-      DO jj = k_j1+1, k_jpj-1
-         DO ji = fs_2, fs_jpim1
-            IF ( vt_i(ji,jj) <= zvmin ) THEN
-               u_ice(ji,jj) = u_oce(ji,jj)
-               v_ice(ji,jj) = v_oce(ji,jj)
-            ENDIF
+      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution) 
+      !------------------------------------------------------------------------------!
+      DO jj = k_j1, k_jpj-1 
+         DO ji = 1, jpim1
+
+            ! shear at F points
+            zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
+               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
+               &         ) * r1_e12f(ji,jj) * zfmask(ji,jj)
+
+         END DO
+      END DO           
+      CALL lbc_lnk( zds, 'F', 1. )
+      
+      DO jj = k_j1+1, k_jpj-1 
+         DO ji = 2, jpim1 ! no vector loop
+            
+            ! tension**2 at T points
+            zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
+               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
+               &   ) * r1_e12t(ji,jj)
+            zdt2 = zdt * zdt
+            
+            ! shear**2 at T points (doc eq. A16)
+            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e12f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e12f(ji-1,jj  )  &
+               &   + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1)  &
+               &   ) * 0.25_wp * r1_e12t(ji,jj)
+            
+            ! shear at T points
+            shear_i(ji,jj) = SQRT( zdt2 + zds2 )
+
+            ! divergence at T points
+            divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
+               &            + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
+               &            ) * r1_e12t(ji,jj)
+            
+            ! delta at T points
+            zdelta         = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * usecc2 )  
+            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
+            delta_i(ji,jj) = zdelta + rn_creepl * rswitch
+
          END DO
       END DO
-
-      CALL lbc_lnk_multi( u_ice(:,:), 'U', -1., v_ice(:,:), 'V', -1. )
-
-#if defined key_agrif && defined key_lim2
-      CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'U' )
-      CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'V' )
-#endif
-#if defined key_bdy
-      CALL bdy_ice_lim_dyn( 'U' )
-      CALL bdy_ice_lim_dyn( 'V' )
-#endif         
-
-      DO jj = k_j1+1, k_jpj-1 
-         DO ji = fs_2, fs_jpim1
-            IF ( vt_i(ji,jj) <= zvmin ) THEN
-               v_ice1(ji,jj)  = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji,  jj-1) ) * e1t(ji+1,jj)     &
-                  &                      + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) )   &
-                  &                      / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)
-
-               u_ice2(ji,jj)  = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
-                  &                      + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) )   &
-                  &                      / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
-            ENDIF 
-         END DO
-      END DO
-
-      CALL lbc_lnk_multi( u_ice2(:,:), 'V', -1., v_ice1(:,:), 'U', -1. )
-
-      ! Recompute delta, shear and div, inputs for mechanical redistribution 
-      DO jj = k_j1+1, k_jpj-1
-         DO ji = fs_2, jpim1   !RB bug no vect opt due to zmask
-            !- divu_i(:,:), zdt(:,:): divergence and tension at centre 
-            !- zds(:,:): shear on northeast corner of grid cells
-            IF ( vt_i(ji,jj) <= zvmin ) THEN
-
-               divu_i(ji,jj) = (  e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj  ) * u_ice(ji-1,jj  )   &
-                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji  ,jj-1) * v_ice(ji  ,jj-1)   &
-                  &            ) * r1_e12t(ji,jj)
-
-               zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)  &
-                  &          -( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)  &
-                  &         ) * r1_e12t(ji,jj)
-               !
-               ! SB modif because ocean has no slip boundary condition 
-               zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  &
-                  &          +( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  &
-                  &         ) * r1_e12f(ji,jj) * ( 2.0 - fmask(ji,jj,1) )                                     &
-                  &         * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1)
-
-               zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )    &
-                  &   + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1) ) * r1_e12t(ji,jj)
-
-               delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 )  
-               delta_i(ji,jj) = delta + rn_creepl
-            
-            ENDIF
-         END DO
-      END DO
-      !
-      !------------------------------------------------------------------------------!
-      ! 5) Store stress tensor and its invariants
-      !------------------------------------------------------------------------------!
-      ! * Invariants of the stress tensor are required for limitd_me
-      !   (accelerates convergence and improves stability)
-      DO jj = k_j1+1, k_jpj-1
-         DO ji = fs_2, fs_jpim1
-            zdst           = (  e2u(ji,jj) * v_ice1(ji,jj) - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)  &   
-               &              + e1v(ji,jj) * u_ice2(ji,jj) - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1) ) * r1_e12t(ji,jj) 
-            shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst )
-         END DO
-      END DO
-
-      ! Lateral boundary condition
-      CALL lbc_lnk_multi( divu_i (:,:), 'T', 1., delta_i(:,:), 'T', 1.,  shear_i(:,:), 'T', 1. )
-
-      ! * Store the stress tensor for the next time step
+      CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. )
+      
+      ! --- Store the stress tensor for the next time step --- !
       stress1_i (:,:) = zs1 (:,:)
       stress2_i (:,:) = zs2 (:,:)
@@ -652 +630 @@
       !
       !------------------------------------------------------------------------------!
-      ! 6) Control prints of residual and charge ellipse
+      ! 5) Control prints of residual and charge ellipse
       !------------------------------------------------------------------------------!
       !
@@ -675 +653 @@
                DO ji = 2, jpim1
                   IF (zpresh(ji,jj) > 1.0) THEN
-                     sigma1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 
-                     sigma2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) )
+                     zsig1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 
+                     zsig2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) )
                      WRITE(charout,FMT="('lim_rhg  :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)")
                      CALL prt_ctl_info(charout)
@@ -687 +665 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
-      CALL wrk_dealloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask               )
-      CALL wrk_dealloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  )
-      CALL wrk_dealloc( jpi,jpj, zs1   , zs2   , zs12   , zresr , zpice                 )
+      CALL wrk_dealloc( jpi,jpj, zpresh, z1_e1t0, z1_e2t0, zp_delt )
+      CALL wrk_dealloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia )
+      CALL wrk_dealloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV )
+      CALL wrk_dealloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
+      CALL wrk_dealloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf )
 
    END SUBROUTINE lim_rhg

branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90

@@ -110 +110 @@
       !!---------------------------------------------------------------------
 
-      ! make calls for heat fluxes before it is modified
-      ! pfrld is the lead fraction at the previous time step (actually between TRP and THD)
-      IF( iom_use('qsr_oce') )   CALL iom_put( "qsr_oce" , qsr_oce(:,:) * pfrld(:,:) )                                   !     solar flux at ocean surface
-      IF( iom_use('qns_oce') )   CALL iom_put( "qns_oce" , qns_oce(:,:) * pfrld(:,:) + qemp_oce(:,:) )                   ! non-solar flux at ocean surface
-      IF( iom_use('qsr_ice') )   CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux at ice surface
-      IF( iom_use('qns_ice') )   CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) ! non-solar flux at ice surface
-      IF( iom_use('qtr_ice') )   CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux transmitted thru ice
-      IF( iom_use('qt_oce' ) )   CALL iom_put( "qt_oce"  , ( qsr_oce(:,:) + qns_oce(:,:) ) * pfrld(:,:) + qemp_oce(:,:) )  
-      IF( iom_use('qt_ice' ) )   CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) )   &
-         &                                                      * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) )
-      IF( iom_use('qemp_oce') )  CALL iom_put( "qemp_oce" , qemp_oce(:,:) )  
-      IF( iom_use('qemp_ice') )  CALL iom_put( "qemp_ice" , qemp_ice(:,:) )  
-      IF( iom_use('emp_oce' ) )  CALL iom_put( "emp_oce"  , emp_oce(:,:) )   ! emp over ocean (taking into account the snow blown away from the ice)
-      IF( iom_use('emp_ice' ) )  CALL iom_put( "emp_ice"  , emp_ice(:,:) )   ! emp over ice   (taking into account the snow blown away from the ice)
-
-      ! albedo output
+      ! make call for albedo output before it is modified
       CALL wrk_alloc( jpi,jpj, zalb )    
 

branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

@@ -116 +116 @@
 
       ! Discriminate between varying salinity (nn_icesal=2) and prescribed cases (other values)
-      SELECT CASE( nn_icesal )                       ! varying salinity or not
-         CASE( 1, 3, 4 ) ;   zswitch_sal = 0       ! prescribed salinity profile
-         CASE( 2 )       ;   zswitch_sal = 1       ! varying salinity profile
+      SELECT CASE( nn_icesal )                  ! varying salinity or not
+         CASE( 1, 3 ) ;   zswitch_sal = 0       ! prescribed salinity profile
+         CASE( 2 )    ;   zswitch_sal = 1       ! varying salinity profile
       END SELECT
 
@@ -651 +651 @@
 
          ! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0)
-         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
          zfmdt          = ( rhosn - rhoic ) * MAX( dh_snowice(ji), 0._wp )    ! <0
          zsstK          = sst_m(ii,ij) + rt0                                
@@ -662 +661 @@
          ! Contribution to salt flux
          sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_1d(ji) * zfmdt * r1_rdtice 
+
+         ! virtual salt flux to keep salinity constant
+         IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
+            sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_m(ii,ij) * a_i_1d(ji) * zfmdt                  * r1_rdtice  & ! put back sss_m into the ocean
+               &                            - sm_i_1d(ji)  * a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice    ! and get  sm_i  from the ocean 
+         ENDIF
           
          ! Contribution to mass flux

branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90

@@ -62 +62 @@
       END DO
  
-      !------------------------------------------------------------------------------|
-      ! 1) Constant salinity, constant in time                                       |
-      !------------------------------------------------------------------------------|
-!!gm comment: if nn_icesal = 1 s_i_new, s_i_1d and sm_i_1d can be set to rn_icesal one for all in the initialisation phase !!
-!!gm           ===>>>   simplification of almost all test on nn_icesal value
-      IF(  nn_icesal == 1  ) THEN
-            s_i_1d (kideb:kiut,1:nlay_i) =  rn_icesal
-            sm_i_1d(kideb:kiut)          =  rn_icesal 
-            s_i_new(kideb:kiut)          =  rn_icesal
-      ENDIF
+      !--------------------------------------------------------------------|
+      ! 1) salinity constant in time                                       |
+      !--------------------------------------------------------------------|
+      ! do nothing
 
-      !------------------------------------------------------------------------------|
-      !  Module 2 : Constant salinity varying in time                                |
-      !------------------------------------------------------------------------------|
+      !----------------------------------------------------------------------|
+      !  2) salinity varying in time                                         |
+      !----------------------------------------------------------------------|
       IF(  nn_icesal == 2  ) THEN
 
@@ -113 +107 @@
 
       !------------------------------------------------------------------------------|
-      !  Module 3 : Profile of salinity, constant in time                            |
+      !  3) vertical profile of salinity, constant in time                           |
       !------------------------------------------------------------------------------|
       IF(  nn_icesal == 3  )   CALL lim_var_salprof1d( kideb, kiut )

branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90

@@ -63 +63 @@
       INTEGER, INTENT(in) ::   kt           ! number of iteration
       !
-      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices
+      INTEGER  ::   ji, jj, jk, jm , jl, jt      ! dummy loop indices
       INTEGER  ::   initad                  ! number of sub-timestep for the advection
       REAL(wp) ::   zcfl , zusnit           !   -      -
@@ -75 +75 @@
       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zhimax                   ! old ice thickness
       REAL(wp), POINTER, DIMENSION(:,:)      ::   zatold, zeiold, zesold   ! old concentration, enthalpies
+      REAL(wp), POINTER, DIMENSION(:,:,:)             ::   zhdfptab
       REAL(wp) ::    zdv, zvi, zvs, zsmv, zes, zei
       REAL(wp) ::    zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
+      !!---------------------------------------------------------------------
+      INTEGER                                ::  ihdf_vars  = 6  !!Number of variables in which we apply horizontal diffusion
+                                                                   !!  inside limtrp for each ice category , not counting the 
+                                                                   !!  variables corresponding to ice_layers 
       !!---------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('limtrp')
@@ -85 +90 @@
       CALL wrk_alloc( jpi,jpj,nlay_i,jpl, z0ei )
       CALL wrk_alloc( jpi,jpj,jpl,        zhimax, zviold, zvsold, zsmvold )
+      CALL wrk_alloc( jpi,jpj,jpl*(ihdf_vars + nlay_i)+1,zhdfptab)
 
       IF( numit == nstart .AND. lwp ) THEN
@@ -170 +176 @@
             z0oi (:,:,jl)   = oa_i (:,:,jl) * e12t(:,:)    ! Age content
             z0es (:,:,jl)   = e_s  (:,:,1,jl) * e12t(:,:)  ! Snow heat content
-            DO jk = 1, nlay_i
+           DO jk = 1, nlay_i
                z0ei  (:,:,jk,jl) = e_i  (:,:,jk,jl) * e12t(:,:) ! Ice  heat content
             END DO
@@ -284 +290 @@
          ! Diffusion of Ice fields                  
          !------------------------------------------------------------------------------!
-
+         !------------------------------------
+         !  Diffusion of other ice variables
+         !------------------------------------
+         jm=1
+         DO jl = 1, jpl
+         !                             ! Masked eddy diffusivity coefficient at ocean U- and V-points
+         !   DO jj = 1, jpjm1                 ! NB: has not to be defined on jpj line and jpi row
+         !      DO ji = 1 , fs_jpim1   ! vector opt.
+         !         pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji  ,jj,jl) ) ) )   &
+         !            &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj,jl) ) ) ) * ahiu(ji,jj)
+         !         pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji,jj  ,jl) ) ) )   &
+         !            &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- a_i(ji,jj+1,jl) ) ) ) * ahiv(ji,jj)
+         !      END DO
+         !   END DO
+            DO jj = 1, jpjm1                 ! NB: has not to be defined on jpj line and jpi row
+               DO ji = 1 , fs_jpim1   ! vector opt.
+                  pahu3D(ji,jj,jl) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji  ,jj,  jl ) ) ) )   &
+                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj,  jl ) ) ) ) * ahiu(ji,jj)
+                  pahv3D(ji,jj,jl) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji,  jj,  jl ) ) ) )   &
+                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- a_i(ji,  jj+1,jl ) ) ) ) * ahiv(ji,jj)
+               END DO
+            END DO
+
+            zhdfptab(:,:,jm)= a_i  (:,:,  jl); jm = jm + 1    
+            zhdfptab(:,:,jm)= v_i  (:,:,  jl); jm = jm + 1
+            zhdfptab(:,:,jm)= v_s  (:,:,  jl); jm = jm + 1 
+            zhdfptab(:,:,jm)= smv_i(:,:,  jl); jm = jm + 1
+            zhdfptab(:,:,jm)= oa_i (:,:,  jl); jm = jm + 1
+            zhdfptab(:,:,jm)= e_s  (:,:,1,jl); jm = jm + 1
+         ! Sample of adding more variables to apply lim_hdf using lim_hdf optimization---
+         !   zhdfptab(:,:,jm) = variable_1 (:,:,1,jl); jm = jm + 1  
+         !   zhdfptab(:,:,jm) = variable_2 (:,:,1,jl); jm = jm + 1 
+         !
+         ! and in this example the parameter ihdf_vars musb be changed to 8 (necessary for allocation)
+         !----------------------------------------------------------------------------------------
+            DO jk = 1, nlay_i
+              zhdfptab(:,:,jm)=e_i(:,:,jk,jl); jm= jm+1
+            END DO
+         END DO
          !
          !--------------------------------
@@ -290 +334 @@
          !--------------------------------
          !                             ! Masked eddy diffusivity coefficient at ocean U- and V-points
+         !DO jj = 1, jpjm1                    ! NB: has not to be defined on jpj line and jpi row
+         !   DO ji = 1 , fs_jpim1   ! vector opt.
+         !      pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji  ,jj) ) ) )   &
+         !         &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj)
+         !      pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj  ) ) ) )   &
+         !         &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj)
+         !   END DO
+         !END DO
+         
          DO jj = 1, jpjm1                    ! NB: has not to be defined on jpj line and jpi row
             DO ji = 1 , fs_jpim1   ! vector opt.
-               pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji  ,jj) ) ) )   &
-                  &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj)
-               pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj  ) ) ) )   &
-                  &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj)
+               pahu3D(ji,jj,jpl+1) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji  ,jj) ) ) )   &
+                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj)
+               pahv3D(ji,jj,jpl+1) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj  ) ) ) )   &
+                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj)
             END DO
          END DO
          !
-         CALL lim_hdf( ato_i (:,:) )
-
-         !------------------------------------
-         !  Diffusion of other ice variables
-         !------------------------------------
-         DO jl = 1, jpl
-         !                             ! Masked eddy diffusivity coefficient at ocean U- and V-points
-            DO jj = 1, jpjm1                 ! NB: has not to be defined on jpj line and jpi row
-               DO ji = 1 , fs_jpim1   ! vector opt.
-                  pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji  ,jj,jl) ) ) )   &
-                     &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj,jl) ) ) ) * ahiu(ji,jj)
-                  pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji,jj  ,jl) ) ) )   &
-                     &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- a_i(ji,jj+1,jl) ) ) ) * ahiv(ji,jj)
-               END DO
-            END DO
-
-            CALL lim_hdf( v_i  (:,:,  jl) )
-            CALL lim_hdf( v_s  (:,:,  jl) )
-            CALL lim_hdf( smv_i(:,:,  jl) )
-            CALL lim_hdf( oa_i (:,:,  jl) )
-            CALL lim_hdf( a_i  (:,:,  jl) )
-            CALL lim_hdf( e_s  (:,:,1,jl) )
+         zhdfptab(:,:,jm)= ato_i  (:,:);
+         CALL lim_hdf( zhdfptab, ihdf_vars, jpl, nlay_i) 
+
+         jm=1
+         DO jl = 1, jpl
+            a_i  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1      
+            v_i  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1 
+            v_s  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1 
+            smv_i(:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1 
+            oa_i (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1 
+            e_s  (:,:,1,jl) = zhdfptab(:,:,jm); jm = jm + 1 
+         ! Sample of adding more variables to apply lim_hdf---------
+         !   variable_1  (:,:,1,jl) = zhdfptab(:,:, jm  ) ; jm + 1 
+         !   variable_2  (:,:,1,jl) = zhdfptab(:,:, jm  ) ; jm + 1
+         !-----------------------------------------------------------
             DO jk = 1, nlay_i
-               CALL lim_hdf( e_i(:,:,jk,jl) )
-            END DO
-         END DO
+               e_i(:,:,jk,jl) = zhdfptab(:,:,jm);jm= jm + 1 
+            END DO
+         END DO
+
+         ato_i  (:,:) = zhdfptab(:,:,jm)
 
          !------------------------------------------------------------------------------!
@@ -464 +512 @@
       CALL wrk_dealloc( jpi,jpj,nlay_i,jpl, z0ei )
       CALL wrk_dealloc( jpi,jpj,jpl,        zviold, zvsold, zhimax, zsmvold )
+      CALL wrk_dealloc( jpi,jpj,jpl*(ihdf_vars+nlay_i)+1,zhdfptab)
       !
       IF( nn_timing == 1 )  CALL timing_stop('limtrp')
@@ -479 +528 @@
    !!======================================================================
 END MODULE limtrp
+

branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90

@@ -54 +54 @@
    PUBLIC   lim_var_eqv2glo      
    PUBLIC   lim_var_salprof      
-   PUBLIC   lim_var_icetm        
    PUBLIC   lim_var_bv           
    PUBLIC   lim_var_salprof1d    
@@ -89 +88 @@
       ! Compute variables
       !--------------------
-      vt_i (:,:) = 0._wp
-      vt_s (:,:) = 0._wp
-      at_i (:,:) = 0._wp
-      ato_i(:,:) = 1._wp
-      !
-      DO jl = 1, jpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               !
-               vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume
-               vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume
-               at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration
-               !
-               rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 
-               icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch  ! ice thickness
-            END DO
-         END DO
-      END DO
-
+      ! integrated values
+      vt_i (:,:) = SUM( v_i, dim=3 )
+      vt_s (:,:) = SUM( v_s, dim=3 )
+      at_i (:,:) = SUM( a_i, dim=3 )
+      et_s(:,:)  = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 )
+      et_i(:,:)  = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 )
+      !
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -115 +102 @@
 
       IF( kn > 1 ) THEN
-         et_s (:,:) = 0._wp
-         ot_i (:,:) = 0._wp
+         !
+         ! mean ice/snow thickness
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               rswitch      = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 
+               htm_i(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch
+               htm_s(ji,jj) = vt_s(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch
+            ENDDO
+         ENDDO
+
+         ! mean temperature (K), salinity and age
          smt_i(:,:) = 0._wp
-         et_i (:,:) = 0._wp
-         !
+         tm_i(:,:)  = 0._wp
+         tm_su(:,:) = 0._wp
+         om_i (:,:) = 0._wp
          DO jl = 1, jpl
+            
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  et_s(ji,jj)  = et_s(ji,jj)  + e_s(ji,jj,1,jl)                                           ! snow heat content
-                  rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi20 ) ) 
-                  smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi20 ) * rswitch   ! ice salinity
-                  rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi20 ) ) 
-                  ot_i(ji,jj)  = ot_i(ji,jj)  + oa_i(ji,jj,jl)  / MAX( at_i(ji,jj) , epsi20 ) * rswitch   ! ice age
-               END DO
-            END DO
-         END DO
-         !
-         DO jl = 1, jpl
+                  rswitch      = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )
+                  tm_su(ji,jj) = tm_su(ji,jj) + rswitch * ( t_su(ji,jj,jl) - rt0 ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi10 )
+                  om_i (ji,jj) = om_i (ji,jj) + rswitch *   oa_i(ji,jj,jl)                         / MAX( at_i(ji,jj) , epsi10 )
+               END DO
+            END DO
+            
             DO jk = 1, nlay_i
-               et_i(:,:) = et_i(:,:) + e_i(:,:,jk,jl)       ! ice heat content
-            END DO
-         END DO
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )
+                     tm_i(ji,jj)  = tm_i(ji,jj)  + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl)  &
+                        &            / MAX( vt_i(ji,jj) , epsi10 )
+                     smt_i(ji,jj) = smt_i(ji,jj) + r1_nlay_i * rswitch * s_i(ji,jj,jk,jl) * v_i(ji,jj,jl)  &
+                        &            / MAX( vt_i(ji,jj) , epsi10 )
+                  END DO
+               END DO
+            END DO
+         END DO
+         tm_i  = tm_i  + rt0
+         tm_su = tm_su + rt0
          !
       ENDIF
@@ -246 +250 @@
       ! Mean temperature
       !-------------------
-      vt_i (:,:) = 0._wp
-      DO jl = 1, jpl
-         vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl)
-      END DO
+      ! integrated values
+      vt_i (:,:) = SUM( v_i, dim=3 )
+      vt_s (:,:) = SUM( v_s, dim=3 )
+      at_i (:,:) = SUM( a_i, dim=3 )
 
       tm_i(:,:) = 0._wp
@@ -314 +318 @@
       ! Vertically constant, constant in time
       !---------------------------------------
-      IF(  nn_icesal == 1  )   s_i(:,:,:,:) = rn_icesal
+      IF(  nn_icesal == 1  )  THEN
+         s_i (:,:,:,:) = rn_icesal
+         sm_i(:,:,:)   = rn_icesal
+      ENDIF
 
       !-----------------------------------
@@ -394 +401 @@
    END SUBROUTINE lim_var_salprof
 
-
-   SUBROUTINE lim_var_icetm
-      !!------------------------------------------------------------------
-      !!                ***  ROUTINE lim_var_icetm ***
-      !!
-      !! ** Purpose :   computes mean sea ice temperature
+   SUBROUTINE lim_var_bv
+      !!------------------------------------------------------------------
+      !!                ***  ROUTINE lim_var_bv ***
+      !!
+      !! ** Purpose :   computes mean brine volume (%) in sea ice
+      !!
+      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
+      !!
+      !! References : Vancoppenolle et al., JGR, 2007
       !!------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
       !!------------------------------------------------------------------
-
-      ! Mean sea ice temperature
-      vt_i (:,:) = 0._wp
-      DO jl = 1, jpl
-         vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl)
-      END DO
-
-      tm_i(:,:) = 0._wp
+      !
+      bvm_i(:,:)   = 0._wp
+      bv_i (:,:,:) = 0._wp
       DO jl = 1, jpl
          DO jk = 1, nlay_i
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )
-                  tm_i(ji,jj) = tm_i(ji,jj) + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl)  &
-                     &            / MAX( vt_i(ji,jj) , epsi10 )
-               END DO
-            END DO
-         END DO
-      END DO
-      tm_i = tm_i + rt0
-
-   END SUBROUTINE lim_var_icetm
-
-
-   SUBROUTINE lim_var_bv
-      !!------------------------------------------------------------------
-      !!                ***  ROUTINE lim_var_bv ***
-      !!
-      !! ** Purpose :   computes mean brine volume (%) in sea ice
-      !!
-      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
-      !!
-      !! References : Vancoppenolle et al., JGR, 2007
-      !!------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
-      REAL(wp) ::   zbvi             ! local scalars
-      !!------------------------------------------------------------------
-      !
-      vt_i (:,:) = 0._wp
-      DO jl = 1, jpl
-         vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl)
-      END DO
-
-      bv_i(:,:) = 0._wp
-      DO jl = 1, jpl
-         DO jk = 1, nlay_i
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt0) + epsi10 ) )  )
-                  zbvi  = - rswitch * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 )   &
-                     &                   * v_i(ji,jj,jl) * r1_nlay_i
-                  rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi20 ) )  )
-                  bv_i(ji,jj) = bv_i(ji,jj) + rswitch * zbvi  / MAX( vt_i(ji,jj) , epsi20 )
-               END DO
+                  rswitch        = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt0) + epsi10 ) )  )
+                  bv_i(ji,jj,jl) = bv_i(ji,jj,jl) - rswitch * tmut * s_i(ji,jj,jk,jl) * r1_nlay_i  &
+                     &                            / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 )
+               END DO
+            END DO
+         END DO
+         
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               rswitch      = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )
+               bvm_i(ji,jj) = bvm_i(ji,jj) + rswitch * bv_i(ji,jj,jl) * v_i(ji,jj,jl) / MAX( vt_i(ji,jj), epsi10 )
             END DO
          END DO
@@ -712 +686 @@
             zht_i(ji,1:jpl) = 0._wp
             za_i (ji,1:jpl) = 0._wp
-            
+            itest(:)        = 0      
+      
             ! *** case very thin ice: fill only category 1
             IF ( i_fill == 1 ) THEN

branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90

@@ -17 +17 @@
    USE sbc_oce         ! Surface boundary condition: ocean fields
    USE sbc_ice         ! Surface boundary condition: ice fields
-   USE dom_ice
    USE ice
    USE limvar
@@ -40 +39 @@
    !!----------------------------------------------------------------------
 CONTAINS
-
-#if defined key_dimgout
-# include "limwri_dimg.h90"
-#else
 
    SUBROUTINE lim_wri( kindic )
@@ -59 +54 @@
       INTEGER  ::  ji, jj, jk, jl  ! dummy loop indices
       REAL(wp) ::  z1_365
-      REAL(wp) ::  ztmp
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zoi, zei, zt_i, zt_s
-      REAL(wp), POINTER, DIMENSION(:,:)   ::  z2d, z2da, z2db, zswi    ! 2D workspace
+      REAL(wp) ::  z2da, z2db, ztmp
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zswi2
+      REAL(wp), POINTER, DIMENSION(:,:)   ::  z2d, zswi    ! 2D workspace
       !!-------------------------------------------------------------------
 
       IF( nn_timing == 1 )  CALL timing_start('limwri')
 
-      CALL wrk_alloc( jpi, jpj, jpl, zoi, zei, zt_i, zt_s )
-      CALL wrk_alloc( jpi, jpj     , z2d, z2da, z2db, zswi )
+      CALL wrk_alloc( jpi, jpj, jpl, zswi2 )
+      CALL wrk_alloc( jpi, jpj     , z2d, zswi )
 
       !-----------------------------
@@ -74 +69 @@
       z1_365 = 1._wp / 365._wp
 
-      CALL lim_var_icetm      ! mean sea ice temperature
-
-      CALL lim_var_bv         ! brine volume
-
-      DO jj = 1, jpj          ! presence indicator of ice
+      ! brine volume
+      CALL lim_var_bv 
+
+      ! tresholds for outputs
+      DO jj = 1, jpj
          DO ji = 1, jpi
             zswi(ji,jj)  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) )
          END DO
       END DO
-      !
-      !
-      !                                             
-      IF ( iom_use( "icethic_cea" ) ) THEN                       ! mean ice thickness
-         DO jj = 1, jpj 
+      DO jl = 1, jpl
+         DO jj = 1, jpj
             DO ji = 1, jpi
-               z2d(ji,jj)  = vt_i(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zswi(ji,jj)
+               zswi2(ji,jj,jl)  = MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - epsi06 ) )
             END DO
          END DO
-         CALL iom_put( "icethic_cea"  , z2d              )
-      ENDIF
-
-      IF ( iom_use( "snowthic_cea" ) ) THEN                      ! snow thickness = mean snow thickness over the cell 
-         DO jj = 1, jpj                                            
-            DO ji = 1, jpi
-               z2d(ji,jj)  = vt_s(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zswi(ji,jj)
-            END DO
-         END DO
-         CALL iom_put( "snowthic_cea" , z2d              )       
-      ENDIF
+      END DO
       !
+      ! fluxes 
+      ! pfrld is the lead fraction at the previous time step (actually between TRP and THD)
+      IF( iom_use('qsr_oce') )   CALL iom_put( "qsr_oce" , qsr_oce(:,:) * pfrld(:,:) )                                   !     solar flux at ocean surface
+      IF( iom_use('qns_oce') )   CALL iom_put( "qns_oce" , qns_oce(:,:) * pfrld(:,:) + qemp_oce(:,:) )                   ! non-solar flux at ocean surface
+      IF( iom_use('qsr_ice') )   CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux at ice surface
+      IF( iom_use('qns_ice') )   CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) ! non-solar flux at ice surface
+      IF( iom_use('qtr_ice') )   CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux transmitted thru ice
+      IF( iom_use('qt_oce' ) )   CALL iom_put( "qt_oce"  , ( qsr_oce(:,:) + qns_oce(:,:) ) * pfrld(:,:) + qemp_oce(:,:) )  
+      IF( iom_use('qt_ice' ) )   CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) )   &
+         &                                                      * a_i_b(:,:,:),dim=3 ) + qemp_ice(:,:) )
+      IF( iom_use('qemp_oce') )  CALL iom_put( "qemp_oce" , qemp_oce(:,:) )  
+      IF( iom_use('qemp_ice') )  CALL iom_put( "qemp_ice" , qemp_ice(:,:) )  
+      IF( iom_use('emp_oce' ) )  CALL iom_put( "emp_oce"  , emp_oce(:,:) )   !emp over ocean (taking into account the snow blown away from the ice)
+      IF( iom_use('emp_ice' ) )  CALL iom_put( "emp_ice"  , emp_ice(:,:) )   !emp over ice   (taking into account the snow blown away from the ice)
+
+      ! velocity
       IF ( iom_use( "uice_ipa" ) .OR. iom_use( "vice_ipa" ) .OR. iom_use( "icevel" ) ) THEN 
          DO jj = 2 , jpjm1
             DO ji = 2 , jpim1
-               z2da(ji,jj)  = (  u_ice(ji,jj) * umask(ji,jj,1) + u_ice(ji-1,jj) * umask(ji-1,jj,1) ) * 0.5_wp
-               z2db(ji,jj)  = (  v_ice(ji,jj) * vmask(ji,jj,1) + v_ice(ji,jj-1) * vmask(ji,jj-1,1) ) * 0.5_wp
+               z2da  = ( u_ice(ji,jj) * umask(ji,jj,1) + u_ice(ji-1,jj) * umask(ji-1,jj,1) ) * 0.5_wp
+               z2db  = ( v_ice(ji,jj) * vmask(ji,jj,1) + v_ice(ji,jj-1) * vmask(ji,jj-1,1) ) * 0.5_wp
+               z2d(ji,jj) = SQRT( z2da * z2da + z2db * z2db )
            END DO
          END DO
-         CALL lbc_lnk( z2da, 'T', -1. )
-         CALL lbc_lnk( z2db, 'T', -1. )
-         CALL iom_put( "uice_ipa"     , z2da             )       ! ice velocity u component
-         CALL iom_put( "vice_ipa"     , z2db             )       ! ice velocity v component
-         DO jj = 1, jpj                                 
-            DO ji = 1, jpi
-               z2d(ji,jj)  = SQRT( z2da(ji,jj) * z2da(ji,jj) + z2db(ji,jj) * z2db(ji,jj) ) 
-            END DO
-         END DO
-         CALL iom_put( "icevel"       , z2d              )       ! ice velocity module
+         CALL lbc_lnk( z2d, 'T', 1. )
+         CALL iom_put( "uice_ipa"     , u_ice      )       ! ice velocity u component
+         CALL iom_put( "vice_ipa"     , v_ice      )       ! ice velocity v component
+         CALL iom_put( "icevel"       , z2d        )       ! ice velocity module
       ENDIF
       !
-      IF ( iom_use( "miceage" ) ) THEN 
-         z2d(:,:) = 0.e0
-         DO jl = 1, jpl
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  rswitch    = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.1 ) )
-                  z2d(ji,jj) = z2d(ji,jj) + rswitch * oa_i(ji,jj,jl) / MAX( at_i(ji,jj), 0.1 )
-               END DO
-            END DO
-         END DO
-         CALL iom_put( "miceage"     , z2d * z1_365      )        ! mean ice age
-      ENDIF
-
-      IF ( iom_use( "micet" ) ) THEN 
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               z2d(ji,jj) = ( tm_i(ji,jj) - rt0 ) * zswi(ji,jj)
-            END DO
-         END DO
-         CALL iom_put( "micet"       , z2d               )        ! mean ice temperature
-      ENDIF
+      IF ( iom_use( "miceage" ) )       CALL iom_put( "miceage"     , om_i * zswi * z1_365   )  ! mean ice age
+      IF ( iom_use( "icethic_cea" ) )   CALL iom_put( "icethic_cea" , htm_i * zswi           )  ! ice thickness mean
+      IF ( iom_use( "snowthic_cea" ) )  CALL iom_put( "snowthic_cea", htm_s * zswi           )  ! snow thickness mean
+      IF ( iom_use( "micet" ) )         CALL iom_put( "micet"       , ( tm_i  - rt0 ) * zswi )  ! ice mean    temperature
+      IF ( iom_use( "icest" ) )         CALL iom_put( "icest"       , ( tm_su - rt0 ) * zswi )  ! ice surface temperature
+      IF ( iom_use( "icecolf" ) )       CALL iom_put( "icecolf"     , hicol                  )  ! frazil ice collection thickness
       !
-      IF ( iom_use( "icest" ) ) THEN 
-         z2d(:,:) = 0.e0
-         DO jl = 1, jpl
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  z2d(ji,jj) = z2d(ji,jj) + zswi(ji,jj) * ( t_su(ji,jj,jl) - rt0 ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi06 )
-               END DO
-            END DO
-         END DO
-         CALL iom_put( "icest"       , z2d              )        ! ice surface temperature
-      ENDIF
-
-      IF ( iom_use( "icecolf" ) )   CALL iom_put( "icecolf", hicol )  ! frazil ice collection thickness
-
       CALL iom_put( "isst"        , sst_m               )        ! sea surface temperature
       CALL iom_put( "isss"        , sss_m               )        ! sea surface salinity
-      CALL iom_put( "iceconc"     , at_i                )        ! ice concentration
-      CALL iom_put( "icevolu"     , vt_i                )        ! ice volume = mean ice thickness over the cell
-      CALL iom_put( "icehc"       , et_i                )        ! ice total heat content
-      CALL iom_put( "isnowhc"     , et_s                )        ! snow total heat content
-      CALL iom_put( "ibrinv"      , bv_i * 100._wp      )        ! brine volume
+      CALL iom_put( "iceconc"     , at_i  * zswi        )        ! ice concentration
+      CALL iom_put( "icevolu"     , vt_i  * zswi        )        ! ice volume = mean ice thickness over the cell
+      CALL iom_put( "icehc"       , et_i  * zswi        )        ! ice total heat content
+      CALL iom_put( "isnowhc"     , et_s  * zswi        )        ! snow total heat content
+      CALL iom_put( "ibrinv"      , bvm_i * zswi * 100. )        ! brine volume
       CALL iom_put( "utau_ice"    , utau_ice            )        ! wind stress over ice along i-axis at I-point
       CALL iom_put( "vtau_ice"    , vtau_ice            )        ! wind stress over ice along j-axis at I-point
       CALL iom_put( "snowpre"     , sprecip * 86400.    )        ! snow precipitation 
-      CALL iom_put( "micesalt"    , smt_i               )        ! mean ice salinity
-
-      CALL iom_put( "icestr"      , strength * 0.001    )        ! ice strength
-      CALL iom_put( "idive"       , divu_i * 1.0e8      )        ! divergence
-      CALL iom_put( "ishear"      , shear_i * 1.0e8     )        ! shear
-      CALL iom_put( "snowvol"     , vt_s                )        ! snow volume
+      CALL iom_put( "micesalt"    , smt_i   * zswi      )        ! mean ice salinity
+
+      CALL iom_put( "icestr"      , strength * zswi )    ! ice strength
+      CALL iom_put( "idive"       , divu_i * 1.0e8      )    ! divergence
+      CALL iom_put( "ishear"      , shear_i * 1.0e8     )    ! shear
+      CALL iom_put( "snowvol"     , vt_s   * zswi       )        ! snow volume
       
       CALL iom_put( "icetrp"      , diag_trp_vi * rday  )        ! ice volume transport
@@ -183 +147 @@
 
       CALL iom_put( "sfxbog"      , sfx_bog * rday      )        ! salt flux from bottom growth
-      CALL iom_put( "sfxbom"      , sfx_bom * rday      )        ! salt flux from bottom melt
-      CALL iom_put( "sfxsum"      , sfx_sum * rday      )        ! salt flux from surface melt
+      CALL iom_put( "sfxbom"      , sfx_bom * rday      )        ! salt flux from bottom melting
+      CALL iom_put( "sfxsum"      , sfx_sum * rday      )        ! salt flux from surface melting
       CALL iom_put( "sfxsni"      , sfx_sni * rday      )        ! salt flux from snow ice formation
       CALL iom_put( "sfxopw"      , sfx_opw * rday      )        ! salt flux from open water formation
       CALL iom_put( "sfxdyn"      , sfx_dyn * rday      )        ! salt flux from ridging rafting
-      CALL iom_put( "sfxres"      , sfx_res * rday      )        ! salt flux from residual
+      CALL iom_put( "sfxres"      , sfx_res * rday      )        ! salt flux from limupdate (resultant)
       CALL iom_put( "sfxbri"      , sfx_bri * rday      )        ! salt flux from brines
       CALL iom_put( "sfxsub"      , sfx_sub * rday      )        ! salt flux from sublimation
@@ -202 +166 @@
       CALL iom_put( "vfxbom"     , wfx_bom * ztmp       )        ! bottom melt 
       CALL iom_put( "vfxice"     , wfx_ice * ztmp       )        ! total ice growth/melt 
+
+      IF ( iom_use( "vfxthin" ) ) THEN   ! ice production for open water + thin ice (<20cm) => comparable to observations  
+         WHERE( htm_i(:,:) < 0.2 .AND. htm_i(:,:) > 0. ) ; z2d = wfx_bog
+         ELSEWHERE                                       ; z2d = 0._wp
+         END WHERE
+         CALL iom_put( "vfxthin", ( wfx_opw + z2d ) * ztmp )
+      ENDIF
+
+      ztmp = rday / rhosn
+      CALL iom_put( "vfxspr"     , wfx_spr * ztmp       )        ! precip (snow)
       CALL iom_put( "vfxsnw"     , wfx_snw * ztmp       )        ! total snw growth/melt 
-      CALL iom_put( "vfxsub"     , wfx_sub * ztmp       )        ! sublimation (snow) 
-      CALL iom_put( "vfxspr"     , wfx_spr * ztmp       )        ! precip (snow)
-      
+      CALL iom_put( "vfxsub"     , wfx_sub * ztmp       )        ! sublimation (snow/ice) 
+      CALL iom_put( "vfxsub_err" , wfx_err_sub * ztmp   )        ! "excess" of sublimation sent to ocean      
+ 
       CALL iom_put( "afxtot"     , afx_tot * rday       )        ! concentration tendency (total)
       CALL iom_put( "afxdyn"     , afx_dyn * rday       )        ! concentration tendency (dynamics)
@@ -225 +199 @@
       CALL iom_put ('hfxdif'     , hfx_dif(:,:)         )   !  
       CALL iom_put ('hfxopw'     , hfx_opw(:,:)         )   !  
-      CALL iom_put ('hfxtur'     , fhtur(:,:) * SUM(a_i_b(:,:,:), dim=3) ) ! turbulent heat flux at ice base 
+      CALL iom_put ('hfxtur'     , fhtur(:,:) * SUM( a_i_b(:,:,:), dim=3 ) ) ! turbulent heat flux at ice base 
       CALL iom_put ('hfxdhc'     , diag_heat(:,:)       )   ! Heat content variation in snow and ice 
       CALL iom_put ('hfxspr'     , hfx_spr(:,:)         )   ! Heat content of snow precip 
 
-
-      IF ( iom_use( "vfxthin" ) ) THEN   ! ice production for open water + thin ice (<20cm) => comparable to observations  
-         DO jj = 1, jpj 
-            DO ji = 1, jpi
-               z2d(ji,jj)  = vt_i(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zswi(ji,jj) ! mean ice thickness
-            END DO
-         END DO
-         WHERE( z2d(:,:) < 0.2 .AND. z2d(:,:) > 0. ) ; z2da = wfx_bog
-         ELSEWHERE                                   ; z2da = 0._wp
-         END WHERE
-         CALL iom_put( "vfxthin", ( wfx_opw + z2da ) * ztmp )
-      ENDIF
       
       !--------------------------------
       ! Output values for each category
       !--------------------------------
-      CALL iom_put( "iceconc_cat"      , a_i         )        ! area for categories
-      CALL iom_put( "icethic_cat"      , ht_i        )        ! thickness for categories
-      CALL iom_put( "snowthic_cat"     , ht_s        )        ! snow depth for categories
-      CALL iom_put( "salinity_cat"     , sm_i        )        ! salinity for categories
-
+      IF ( iom_use( "iceconc_cat"  ) )  CALL iom_put( "iceconc_cat"      , a_i   * zswi2   )        ! area for categories
+      IF ( iom_use( "icethic_cat"  ) )  CALL iom_put( "icethic_cat"      , ht_i  * zswi2   )        ! thickness for categories
+      IF ( iom_use( "snowthic_cat" ) )  CALL iom_put( "snowthic_cat"     , ht_s  * zswi2   )        ! snow depth for categories
+      IF ( iom_use( "salinity_cat" ) )  CALL iom_put( "salinity_cat"     , sm_i  * zswi2   )        ! salinity for categories
       ! ice temperature
-      IF ( iom_use( "icetemp_cat" ) ) THEN 
-         zt_i(:,:,:) = SUM( t_i(:,:,:,:), dim=3 ) * r1_nlay_i
-         CALL iom_put( "icetemp_cat"   , zt_i - rt0  )
-      ENDIF
-      
+      IF ( iom_use( "icetemp_cat"  ) )  CALL iom_put( "icetemp_cat", ( SUM( t_i(:,:,:,:), dim=3 ) * r1_nlay_i - rt0 ) * zswi2 )
       ! snow temperature
-      IF ( iom_use( "snwtemp_cat" ) ) THEN 
-         zt_s(:,:,:) = SUM( t_s(:,:,:,:), dim=3 ) * r1_nlay_s
-         CALL iom_put( "snwtemp_cat"   , zt_s - rt0  )
-      ENDIF
-
-      ! Compute ice age
-      IF ( iom_use( "iceage_cat" ) ) THEN 
-         DO jl = 1, jpl 
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.1 ) )
-                  rswitch = rswitch * MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - 0.1 ) )
-                  zoi(ji,jj,jl) = oa_i(ji,jj,jl)  / MAX( a_i(ji,jj,jl) , 0.1 ) * rswitch
-               END DO
-            END DO
-         END DO
-         CALL iom_put( "iceage_cat"   , zoi * z1_365 )        ! ice age for categories
-      ENDIF
-
-      ! Compute brine volume
-      IF ( iom_use( "brinevol_cat" ) ) THEN 
-         zei(:,:,:) = 0._wp
-         DO jl = 1, jpl 
-            DO jk = 1, nlay_i
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     rswitch = MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - epsi06 ) )
-                     zei(ji,jj,jl) = zei(ji,jj,jl) + 100.0 *  &
-                        ( - tmut * s_i(ji,jj,jk,jl) / MIN( ( t_i(ji,jj,jk,jl) - rt0 ), - epsi06 ) ) * &
-                        rswitch * r1_nlay_i
-                  END DO
-               END DO
-            END DO
-         END DO
-         CALL iom_put( "brinevol_cat"     , zei      )        ! brine volume for categories
-      ENDIF
+      IF ( iom_use( "snwtemp_cat"  ) )  CALL iom_put( "snwtemp_cat", ( SUM( t_s(:,:,:,:), dim=3 ) * r1_nlay_s - rt0 ) * zswi2 )
+      ! ice age
+      IF ( iom_use( "iceage_cat"   ) )  CALL iom_put( "iceage_cat" , o_i * zswi2 * z1_365 )
+      ! brine volume
+      IF ( iom_use( "brinevol_cat" ) )  CALL iom_put( "brinevol_cat", bv_i * 100. * zswi2 )
 
       !     !  Create an output files (output.lim.abort.nc) if S < 0 or u > 20 m/s
@@ -298 +224 @@
       !     not yet implemented
       
-      CALL wrk_dealloc( jpi, jpj, jpl, zoi, zei, zt_i, zt_s )
-      CALL wrk_dealloc( jpi, jpj     , z2d, zswi, z2da, z2db )
+      CALL wrk_dealloc( jpi, jpj, jpl, zswi2 )
+      CALL wrk_dealloc( jpi, jpj     , z2d, zswi )
 
       IF( nn_timing == 1 )  CALL timing_stop('limwri')
       
    END SUBROUTINE lim_wri
-#endif
 
  
@@ -319 +244 @@
       !!   4.0  !  2013-06  (C. Rousset)
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt               ! ocean time-step index)
-      INTEGER, INTENT( in ) ::   kid , kh_i       
+      INTEGER, INTENT( in )   ::   kt               ! ocean time-step index)
+      INTEGER, INTENT( in )   ::   kid , kh_i
+      INTEGER                 ::   nz_i, jl
+      REAL(wp), DIMENSION(jpl) :: jcat
       !!----------------------------------------------------------------------
-
-      CALL histdef( kid, "iicethic", "Ice thickness"           , "m"      ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iiceconc", "Ice concentration"       , "%"      ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicetemp", "Ice temperature"         , "C"      ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicevelu", "i-Ice speed (I-point)"   , "m/s"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicevelv", "j-Ice speed (I-point)"   , "m/s"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicestru", "i-Wind stress over ice (I-pt)", "Pa",   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicestrv", "j-Wind stress over ice (I-pt)", "Pa",   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicesflx", "Solar flux over ocean"     , "w/m2" ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicenflx", "Non-solar flux over ocean" , "w/m2" ,   &
+      DO jl = 1, jpl
+         jcat(jl) = REAL(jl)
+      ENDDO
+      
+      CALL histvert( kid, "ncatice", "Ice Categories","", jpl, jcat, nz_i, "up")
+
+      CALL histdef( kid, "sithic", "Ice thickness"           , "m"      ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "siconc", "Ice concentration"       , "%"      ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "sitemp", "Ice temperature"         , "C"      ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "sivelu", "i-Ice speed "            , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "sivelv", "j-Ice speed "            , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+      CALL histdef( kid, "sistru", "i-Wind stress over ice " , "Pa"     ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "sistrv", "j-Wind stress over ice " , "Pa"     ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+      CALL histdef( kid, "sisflx", "Solar flux over ocean"     , "w/m2" ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+      CALL histdef( kid, "sinflx", "Non-solar flux over ocean" , "w/m2" ,   &
       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
       CALL histdef( kid, "isnowpre", "Snow precipitation"      , "kg/m2/s",   &
       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicesali", "Ice salinity"            , "PSU"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicevolu", "Ice volume"              , "m"      ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicedive", "Ice divergence"          , "10-8s-1",   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicebopr", "Ice bottom production"   , "m/s"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicedypr", "Ice dynamic production"  , "m/s"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicelapr", "Ice open water prod"     , "m/s"    ,   &
+      CALL histdef( kid, "sisali", "Ice salinity"            , "PSU"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+      CALL histdef( kid, "sivolu", "Ice volume"              , "m"      ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+      CALL histdef( kid, "sidive", "Ice divergence"          , "10-8s-1",   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+
+      CALL histdef( kid, "vfxbog", "Ice bottom production"   , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "vfxdyn", "Ice dynamic production"  , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "vfxopw", "Ice open water prod"     , "m/s"    ,   &
       &       jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicesipr", "Snow ice production "    , "m/s"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicerepr", "Ice prod from limupdate" , "m/s"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicebome", "Ice bottom melt"         , "m/s"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicesume", "Ice surface melt"        , "m/s"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iisfxdyn", "Salt flux from dynmics"  , ""       ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iisfxres", "Salt flux from limupdate", ""       ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "vfxsni", "Snow ice production "    , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "vfxres", "Ice prod from limupdate" , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "vfxbom", "Ice bottom melt"         , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "vfxsum", "Ice surface melt"        , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+
+      CALL histdef( kid, "sithicat", "Ice thickness"         , "m"      ,   &
+      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "siconcat", "Ice concentration"     , "%"      ,   &
+      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "sisalcat", "Ice salinity"           , ""      ,   &
+      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "sitemcat", "Ice temperature"       , "C"      ,   &
+      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "snthicat", "Snw thickness"         , "m"      ,   &
+      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "sntemcat", "Snw temperature"       , "C"      ,   &
+      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )
 
       CALL histend( kid, snc4set )   ! end of the file definition
 
-      CALL histwrite( kid, "iicethic", kt, icethi        , jpi*jpj, (/1/) )    
-      CALL histwrite( kid, "iiceconc", kt, at_i          , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicetemp", kt, tm_i - rt0    , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicevelu", kt, u_ice          , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicevelv", kt, v_ice          , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicestru", kt, utau_ice       , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicestrv", kt, vtau_ice       , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicesflx", kt, qsr , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicenflx", kt, qns , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sithic", kt, htm_i         , jpi*jpj, (/1/) )    
+      CALL histwrite( kid, "siconc", kt, at_i          , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sitemp", kt, tm_i - rt0    , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sivelu", kt, u_ice          , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sivelv", kt, v_ice          , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sistru", kt, utau_ice       , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sistrv", kt, vtau_ice       , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sisflx", kt, qsr , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sinflx", kt, qns , jpi*jpj, (/1/) )
       CALL histwrite( kid, "isnowpre", kt, sprecip        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicesali", kt, smt_i          , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicevolu", kt, vt_i           , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicedive", kt, divu_i*1.0e8   , jpi*jpj, (/1/) )
-
-      CALL histwrite( kid, "iicebopr", kt, wfx_bog        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicedypr", kt, wfx_dyn        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicelapr", kt, wfx_opw        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicesipr", kt, wfx_sni        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicerepr", kt, wfx_res        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicebome", kt, wfx_bom        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicesume", kt, wfx_sum        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iisfxdyn", kt, sfx_dyn        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iisfxres", kt, sfx_res        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sisali", kt, smt_i          , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sivolu", kt, vt_i           , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sidive", kt, divu_i*1.0e8   , jpi*jpj, (/1/) )
+
+      CALL histwrite( kid, "vfxbog", kt, wfx_bog        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "vfxdyn", kt, wfx_dyn        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "vfxopw", kt, wfx_opw        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "vfxsni", kt, wfx_sni        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "vfxres", kt, wfx_res        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "vfxbom", kt, wfx_bom        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "vfxsum", kt, wfx_sum        , jpi*jpj, (/1/) )
+
+      CALL histwrite( kid, "sithicat", kt, ht_i        , jpi*jpj*jpl, (/1/) )    
+      CALL histwrite( kid, "siconcat", kt, a_i         , jpi*jpj*jpl, (/1/) )    
+      CALL histwrite( kid, "sisalcat", kt, sm_i        , jpi*jpj*jpl, (/1/) )    
+      CALL histwrite( kid, "sitemcat", kt, tm_i - rt0  , jpi*jpj*jpl, (/1/) )    
+      CALL histwrite( kid, "snthicat", kt, ht_s        , jpi*jpj*jpl, (/1/) )    
+      CALL histwrite( kid, "sntemcat", kt, tm_su - rt0 , jpi*jpj*jpl, (/1/) )    
 
       ! Close the file