Changeset 8813

Timestamp:

2017-11-24T17:56:51+01:00 (6 years ago)

Author:

gm

Message:

#1911 (ENHANCE-09): PART I.3 - phasing with updated branch dev_r8183_ICEMODEL revision 8787

Location:

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM

Files:

• CONFIG/SHARED/namelist_ice_lim3_ref (modified) (4 diffs)
• NEMO/LIM_SRC_3/ice.F90 (modified) (2 diffs)
• NEMO/LIM_SRC_3/ice1d.F90 (modified) (7 diffs)
• NEMO/LIM_SRC_3/icedyn_rhg.F90 (modified) (4 diffs)
• NEMO/LIM_SRC_3/icedyn_rhg_evp.F90 (modified) (23 diffs)
• NEMO/LIM_SRC_3/iceforcing.F90 (modified) (4 diffs)
• NEMO/LIM_SRC_3/iceistate.F90 (modified) (28 diffs)
• NEMO/LIM_SRC_3/iceitd.F90 (modified) (24 diffs)
• NEMO/LIM_SRC_3/icestp.F90 (modified) (2 diffs)
• NEMO/LIM_SRC_3/icethd.F90 (modified) (4 diffs)
• NEMO/LIM_SRC_3/icethd_dh.F90 (modified) (1 diff)
• NEMO/LIM_SRC_3/icethd_zdf.F90 (modified) (23 diffs)
• NEMO/LIM_SRC_3/icevar.F90 (modified) (20 diffs)
• NEMO/OPA_SRC/BDY/bdydta.F90 (modified) (45 diffs)
• NEMO/OPA_SRC/SBC/sbc_ice.F90 (modified) (6 diffs)
• NEMO/OPA_SRC/SBC/sbcblk.F90 (modified) (18 diffs)
• NEMO/OPA_SRC/SBC/sbcblk_algo_coare.F90 (modified) (1 diff)
• NEMO/OPA_SRC/SBC/sbccpl.F90 (modified) (17 diffs)
• NEMO/OPA_SRC/SBC/sbcice_cice.F90 (modified) (1 diff)
• NEMO/OPA_SRC/SBC/sbcisf.F90 (modified) (4 diffs)

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref

@@ -41 +41 @@
 &namitd         !   Ice discretization
 !------------------------------------------------------------------------------
-   rn_himean        =   2.0           !  expected domain-average ice thickness (m)
+   ln_cat_hfn       = .true.          !  ice categories are defined by a function following rn_himean**(-0.05)
+      rn_himean     =   2.0           !  expected domain-average ice thickness (m)
+   ln_cat_usr       = .false.         !  ice categories are defined by rn_catbnd below (m)
+      rn_catbnd     =   0.,0.45,1.1,2.1,3.7,6.0  
    rn_himin         =   0.1           !  minimum ice thickness (m) used in remapping
 /
@@ -90 +93 @@
 !------------------------------------------------------------------------------
    ln_rhg_EVP       = .true.          !  EVP rheology
-      rn_creepl     =   1.0e-12       !     creep limit [1/s]
+      ln_aEVP       = .false.         !     adaptive rheology (Kimmritz et al. 2016 & 2017)
+      rn_creepl     =   2.0e-9        !     creep limit [1/s]
       rn_ecc        =   2.0           !     eccentricity of the elliptical yield curve          
       nn_nevp       = 120             !     number of EVP subcycles                             
@@ -118 +122 @@
                                       !     = 2  Redistribute a single flux over categories
                                       !          ==> coupled mode only
+   nice_jules       = 1               !  Jules coupling (0=OFF, 1=EMULATED, 2=ACTIVE)
 /
 !------------------------------------------------------------------------------
@@ -136 +141 @@
                                       !     Obs: 0.1-0.5 (Lecomte et al, JAMES 2013)
    rn_kappa_i       =   1.0           !  radiation attenuation coefficient in sea ice [1/m]
-   ln_dqns_i        = .true.          !  change the surface non-solar flux with surface temperature (T) or not (F)
 /
 !------------------------------------------------------------------------------

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

@@ -178 +178 @@
    !
    !                                     !!** ice-rheology namelist (namrhg) **
+   LOGICAL , PUBLIC ::   ln_aEVP          !: using adaptive EVP (T or F) 
    REAL(wp), PUBLIC ::   rn_creepl        !: creep limit : has to be under 1.0e-9
    REAL(wp), PUBLIC ::   rn_ecc           !: eccentricity of the elliptical yield curve
@@ -194 +195 @@
    !                                      !                                   using T-ice and albedo sensitivity
    !                                      !   = 2  Redistribute a single flux over categories
+
+   INTEGER , PUBLIC            ::   nice_jules            !: Choice of jules coupling 
+   !                                                      !  Associated indices
+   INTEGER , PUBLIC, PARAMETER ::   np_jules_OFF    = 0   !: no Jules coupling (ice thermodynamics forced via qsr and qns)
+   INTEGER , PUBLIC, PARAMETER ::   np_jules_EMULE  = 1   !: emulated Jules coupling via icethd_zdf.F90 (BL99) (1st round compute qcn and qsr_tr, 2nd round use it)
+   INTEGER , PUBLIC, PARAMETER ::   np_jules_ACTIVE = 2   !: active Jules coupling                      (SM0L) (compute qcn and qsr_tr via sbcblk.F90 or sbccpl.F90)
 
    !                                     !!** ice-salinity namelist (namthd_sal) **

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/ice1d.F90

@@ -33 +33 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ftr_ice_1d   
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qsr_ice_1d  
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fr1_i0_1d   
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fr2_i0_1d   
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qns_ice_1d  
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t_bo_1d     
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rn_amax_1d
+   
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qml_ice_1d     !: heat available for snow / ice surface melting [W/m2] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qcn_ice_1d     !: heat available for snow / ice surface sublimation¬ [W/m2] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qsr_ice_tr_1d  !: solar flux transmitted below the ice surface¬ [W/m2] 
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_sum_1d
@@ -96 +98 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   evap_ice_1d   !: <==> the 2D  evap_ice
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qprec_ice_1d  !: <==> the 2D  qprec_ice
-   !                                                     ! to reintegrate longwave flux inside the ice thermodynamics
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   i0            !: fraction of radiation transmitted to the ice
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t_su_1d       !: <==> the 2D  t_su
@@ -178 +178 @@
       ALLOCATE( nptidx   (jpij) ,   &
          &      qlead_1d (jpij) , ftr_ice_1d(jpij) , qsr_ice_1d (jpij) ,   &
-         &      fr1_i0_1d(jpij) , fr2_i0_1d (jpij) , qns_ice_1d(jpij)  ,   &
+         &      qns_ice_1d(jpij)  ,   &
+         &      qml_ice_1d(jpij), qcn_ice_1d(jpij) , qsr_ice_tr_1d(jpij) , &
          &      t_bo_1d   (jpij) ,                                         &
          &      hfx_sum_1d(jpij) , hfx_bom_1d(jpij) ,hfx_bog_1d(jpij) ,    & 
@@ -194 +195 @@
          &      wfx_snw_sub_1d(jpij), wfx_snw_dyn_1d(jpij), wfx_ice_sub_1d(jpij), wfx_err_sub_1d(jpij) ,   &
          &      wfx_lam_1d(jpij)  , wfx_dyn_1d(jpij), wfx_pnd_1d(jpij),dqns_ice_1d(jpij) , evap_ice_1d (jpij),   &
-         &      qprec_ice_1d(jpij), i0         (jpij) ,                                         &  
+         &      qprec_ice_1d(jpij),                                                             &  
          &      sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) , sfx_sum_1d (jpij),  &
          &      sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , sfx_sub_1d (jpij),  &
@@ -200 +201 @@
       !
       ii = ii + 1
-      ALLOCATE( t_su_1d  (jpij) , t_si_1d   (jpij) , a_i_1d  (jpij) , a_ib_1d(jpij) ,                &
-         &      h_i_1d  (jpij) , h_ib_1d  (jpij) , h_s_1d (jpij) , fc_su  (jpij) , fc_bo_i(jpij) ,   &    
-         &      dh_s_tot (jpij) , dh_i_surf (jpij) , dh_i_sub(jpij) ,                                &    
+      ALLOCATE( t_su_1d  (jpij) , t_si_1d   (jpij) , a_i_1d  (jpij) , a_ib_1d(jpij) ,                  &
+         &      h_i_1d  (jpij) , h_ib_1d  (jpij) , h_s_1d (jpij) , fc_su  (jpij) , fc_bo_i(jpij) ,  &    
+         &      dh_s_tot (jpij) , dh_i_surf (jpij) , dh_i_sub(jpij) ,                                  &    
          &      dh_i_bott(jpij) , dh_s_mlt(jpij), dh_snowice(jpij) , s_i_1d (jpij) , s_i_new(jpij) , &
          &      a_ip_1d  (jpij) , v_ip_1d   (jpij) , v_i_1d  (jpij) , v_s_1d (jpij) , &
@@ -210 +211 @@
       ii = ii + 1
       ALLOCATE( t_s_1d  (jpij,nlay_s)     , t_i_1d (jpij,nlay_i)     , sz_i_1d(jpij,nlay_i) ,  &            
-         &      e_i_1d  (jpij,nlay_i)     , e_s_1d (jpij,nlay_s)     ,                         &
+         &      e_i_1d  (jpij,nlay_i)     , e_s_1d (jpij,nlay_s)     ,                        &
          &      eh_i_old(jpij,0:nlay_i+1) , h_i_old(jpij,0:nlay_i+1) , STAT=ierr(ii) )
       !
@@ -220 +221 @@
       !
       ii = ii + 1
-      ALLOCATE( a_i_2d (jpij,jpl) , a_ib_2d(jpij,jpl) , h_i_2d (jpij,jpl) , h_ib_2d(jpij,jpl) ,   &
-         &      v_i_2d (jpij,jpl) , v_s_2d (jpij,jpl) , oa_i_2d(jpij,jpl) , sv_i_2d(jpij,jpl) ,   &
-         &      a_ip_2d(jpij,jpl) , v_ip_2d(jpij,jpl) , t_su_2d(jpij,jpl) ,                   STAT=ierr(ii) )
+      ALLOCATE( a_i_2d(jpij,jpl) , a_ib_2d(jpij,jpl) , h_i_2d(jpij,jpl) , h_ib_2d(jpij,jpl) , &
+         &      v_i_2d(jpij,jpl) ,v_s_2d(jpij,jpl) ,oa_i_2d(jpij,jpl) ,sv_i_2d(jpij,jpl) ,  &
+         &      a_ip_2d(jpij,jpl) ,v_ip_2d(jpij,jpl) ,t_su_2d(jpij,jpl) ,  &
+         &      STAT=ierr(ii) )
 
       ice1D_alloc = MAXVAL( ierr(:) )

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rhg.F90

@@ -78 +78 @@
       CASE( np_rhgEVP )                ! Elasto-Viscous-Plastic !
          !                             !------------------------!
-         CALL ice_dyn_rhg_evp( kt, stress1_i, stress2_i, stress12_i, u_ice, v_ice, shear_i, divu_i, delta_i )
+         CALL ice_dyn_rhg_evp( kt, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i )
          !         
       END SELECT
@@ -107 +107 @@
       INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
       !!
-      NAMELIST/namdyn_rhg/  ln_rhg_EVP, rn_creepl, rn_ecc , nn_nevp, rn_relast
+      NAMELIST/namdyn_rhg/  ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast
       !!-------------------------------------------------------------------
       !
@@ -125 +125 @@
          WRITE(numout,*) '   Namelist namdyn_rhg:'
          WRITE(numout,*) '      rheology EVP (icedyn_rhg_evp)                        ln_rhg_EVP = ', ln_rhg_EVP
+         WRITE(numout,*) '         use adaptive EVP (aEVP)                           ln_aEVP    = ', ln_aEVP
          WRITE(numout,*) '         creep limit                                       rn_creepl  = ', rn_creepl
          WRITE(numout,*) '         eccentricity of the elliptical yield curve        rn_ecc     = ', rn_ecc
@@ -137 +138 @@
       IF( ioptio /= 1 )   CALL ctl_stop( 'ice_dyn_rhg_init: choose one and only one ice rheology' )
       !
-      IF( ln_rhg_EVP )   CALL rhg_evp_rst( 'READ' )  !* read or initialize all required files
+      IF( ln_rhg_EVP  )   CALL rhg_evp_rst( 'READ' )  !* read or initialize all required files
       !
    END SUBROUTINE ice_dyn_rhg_init

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rhg_evp.F90

@@ -16 +16 @@
    !!   'key_lim3'                                       ESIM sea-ice model
    !!----------------------------------------------------------------------
-   !!   ice_dyn_rhg_evp       : computes ice velocities from EVP rheology
+   !!   ice_dyn_rhg_evp : computes ice velocities from EVP rheology
+   !!   rhg_evp_rst     : read/write EVP fields in ice restart
    !!----------------------------------------------------------------------
    USE phycst         ! Physical constant
@@ -24 +25 @@
    USE ice            ! sea-ice: ice variables
    USE icedyn_rdgrft  ! sea-ice: ice strength
+   USE bdy_oce , ONLY : ln_bdy 
+   USE bdyice 
+#if defined key_agrif
+   USE agrif_lim3_interp
+#endif
    !
    USE in_out_manager ! I/O manager
@@ -31 +37 @@
    USE lbclnk         ! lateral boundary conditions (or mpp links)
    USE prtctl         ! Print control
-   !
-#if defined key_agrif
-   USE agrif_lim3_interp
-#endif
-   USE bdy_oce   , ONLY: ln_bdy 
-   USE bdyice 
 
    IMPLICIT NONE
@@ -53 +53 @@
 CONTAINS
 
-   SUBROUTINE ice_dyn_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, u_ice, v_ice, pshear_i, pdivu_i, pdelta_i )
+   SUBROUTINE ice_dyn_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i )
       !!-------------------------------------------------------------------
       !!                 ***  SUBROUTINE ice_dyn_rhg_evp  ***
-      !!                          EVP-C-grid
+      !!                             EVP-C-grid
       !!
       !! ** purpose : determines sea ice drift from wind stress, ice-ocean
@@ -95 +95 @@
       !!              5) Diagnostics including charge ellipse
       !!
-      !! ** Notes   : There is the possibility to use mEVP from Bouillon 2013
-      !!              (by uncommenting some lines in part 3 and changing alpha and beta parameters)
-      !!              but this solution appears very unstable (see Kimmritz et al 2016)
+      !! ** Notes   : There is the possibility to use aEVP from the nice work of Kimmritz et al. (2016 & 2017)
+      !!              by setting up ln_aEVP=T (i.e. changing alpha and beta parameters).
+      !!              This is an upgraded version of mEVP from Bouillon et al. 2013
+      !!              (i.e. more stable and better convergence)
       !!
       !! References : Hunke and Dukowicz, JPO97
       !!              Bouillon et al., Ocean Modelling 2009
       !!              Bouillon et al., Ocean Modelling 2013
+      !!              Kimmritz et al., Ocean Modelling 2016 & 2017
       !!-------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt      ! time step
-      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i
-      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   u_ice, v_ice 
-      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i, pdivu_i, pdelta_i 
+      INTEGER                 , INTENT(in   ) ::   kt                                    ! time step
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i   !
+      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      !
       !!
       INTEGER ::   ji, jj       ! dummy loop indices
       INTEGER ::   jter         ! local integers
-
+      !
       REAL(wp) ::   zrhoco                                                   ! rau0 * rn_cio
       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) ::   zalph1, z1_alph1, zalph2, z1_alph2                       ! alpha coef from Bouillon 2009 or Kimmritz 2017
       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, zTauB, zTauE, zvel                                ! temporary scalars
-
+      !
       REAL(wp) ::   zresm                                                    ! Maximal error on ice velocity
       REAL(wp) ::   zintb, zintn                                             ! dummy argument
       REAL(wp) ::   zfac_x, zfac_y
       REAL(wp) ::   zshear, zdum1, zdum2
-      
+      !
       REAL(wp), DIMENSION(jpi,jpj) ::   z1_e1t0, z1_e2t0                ! scale factors
       REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points
-      !
+      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017
+      !
+      REAL(wp), DIMENSION(jpi,jpj) ::   zdt_m                           ! (dt / ice-snow_mass) on T points
       REAL(wp), DIMENSION(jpi,jpj) ::   zaU   , zaV                     ! ice fraction on U/V points
-      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! ice/snow mass/dt on U/V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! (ice-snow_mass / dt) on U/V points
       REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points
       REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points
@@ -134 +137 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                           
       REAL(wp), DIMENSION(jpi,jpj) ::   zfU   , zfV                     ! internal stresses
-      
+      !
       REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear
       REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components
       REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice, zresr           ! check convergence
       REAL(wp), DIMENSION(jpi,jpj) ::   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   
+      !                                                                      !   ocean surface (ssh_m) if ice is not embedded
+      !                                                                       !   ice top surface if ice is embedded   
       REAL(wp), DIMENSION(jpi,jpj) ::   zCorx, zCory                    ! Coriolis stress array
       REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! Ocean-to-ice stress array
-
+      !
       REAL(wp), DIMENSION(jpi,jpj) ::   zswitchU, zswitchV              ! dummy arrays
       REAL(wp), DIMENSION(jpi,jpj) ::   zmaskU, zmaskV                  ! mask for ice presence
@@ -179 +182 @@
 #endif
       !
+!!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization....
       !------------------------------------------------------------------------------!
       ! 0) mask at F points for the ice
@@ -231 +235 @@
 
       ! alpha parameters (Bouillon 2009)
-      zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
-      zalph2 = zalph1 * z1_ecc2
-
-      ! alpha and beta parameters (Bouillon 2013)
-      !!zalph1 = 40.
-      !!zalph2 = 40.
-      !!zbeta  = 3000.
-      !!zbeta = REAL( nn_nevp )   ! close to classical EVP of Hunke (2001)
-
-      z1_alph1 = 1._wp / ( zalph1 + 1._wp )
-      z1_alph2 = 1._wp / ( zalph2 + 1._wp )
-
+      IF( .NOT. ln_aEVP ) THEN
+         zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
+         zalph2 = zalph1 * z1_ecc2
+
+         z1_alph1 = 1._wp / ( zalph1 + 1._wp )
+         z1_alph2 = 1._wp / ( zalph2 + 1._wp )
+      ENDIF
+         
       ! Initialise stress tensor 
       zs1 (:,:) = pstress1_i (:,:) 
@@ -266 +266 @@
       IF( ln_ice_embd ) THEN             !== embedded sea ice: compute representative ice top surface ==!
          !                                            
-         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
-         !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
+         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
+         !                                               = (1/nn_fsbc)^2 * {SUM[n]        , n=0,nn_fsbc-1}
          zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp     
          !
-         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
+         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   *    {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
          !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
          zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
@@ -304 +304 @@
             zmf(ji,jj)      = zm1 * ff_t(ji,jj)
 
+            ! dt/m at T points (for alpha and beta coefficients)
+            zdt_m(ji,jj)    = zdtevp / MAX( zm1, zmmin )
+            
             ! 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)
@@ -326 +329 @@
          END DO
       END DO
-      CALL lbc_lnk( zmf, 'T', 1. )
+      CALL lbc_lnk_multi( zmf, 'T', 1., zdt_m, 'T', 1. )
       !
       !------------------------------------------------------------------------------!
@@ -380 +383 @@
                ! P/delta at T points
                zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
+
+               ! alpha & beta for aEVP
+               !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m
+               !   alpha = beta = sqrt(4*gamma)
+               IF( ln_aEVP ) THEN
+                  zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
+                  z1_alph1 = 1._wp / ( zalph1 + 1._wp )
+                  zalph2   = zalph1
+                  z1_alph2 = z1_alph1
+               ENDIF
                
                ! stress at T points
@@ -392 +405 @@
             DO ji = 1, jpim1
 
+               ! alpha & beta for aEVP
+               IF( ln_aEVP ) THEN
+                  zalph2   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
+                  z1_alph2 = 1._wp / ( zalph2 + 1._wp )
+                  zbeta(ji,jj) = zalph2
+               ENDIF
+               
                ! 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) )
@@ -411 +431 @@
                   &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
                   &                  ) * r1_e1e2u(ji,jj)
-                  !
-                  !                !--- V points
+               !
+               !                !--- 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)    &
@@ -419 +439 @@
                   &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
                   &                  ) * r1_e1e2v(ji,jj)
-                  !
-                  !                !--- u_ice at V point
+               !
+               !                !--- 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_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 --- !
-         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp
+         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta vary as in Kimmritz 2016 & 2017
          !  Bouillon et al. 2009 (eq 34-35) => stable
          IF( MOD(jter,2) == 0 ) THEN ! even iterations
@@ -438 +457 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
-                     !                 !--- tau_io/(v_oce - v_ice)
+                  !                 !--- tau_io/(v_oce - v_ice)
                   zTauO = zaV(ji,jj) * zrhoco * 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) ) )
-                     !                 !--- Ocean-to-Ice stress
+                  !                 !--- Ocean-to-Ice stress
                   ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
-                     !
-                     !                 !--- tau_bottom/v_ice
+                  !
+                  !                 !--- tau_bottom/v_ice
                   zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
                   zTauB = - tau_icebfr(ji,jj) / zvel
-                     !
-                     !                 !--- Coriolis at V-points (energy conserving formulation)
+                  !
+                  !                 !--- Coriolis at V-points (energy conserving formulation)
                   zCory(ji,jj)  = - 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
+                  !
+                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
                   zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
-                     !
-                     !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
+                  !
+                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
-                     !
-                     !                 !--- ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  !
+                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
+                  v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
+                     &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                     &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
+                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin
+                     &           ) * zmaskV(ji,jj)
+                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
                   v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             &  ! previous velocity
                      &                                     + zTauE + zTauO * v_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
@@ -466 +492 @@
                      &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
                      &            ) * zmaskV(ji,jj)
-                     !
-                  ! Bouillon 2013
-                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
-                  !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
-                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
-                  !
+                  ENDIF
                END DO
             END DO
@@ -483 +504 @@
             !
             DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                               
-                  ! tau_io/(u_oce - u_ice)
+               DO ji = fs_2, fs_jpim1          
+                  !                 !--- tau_io/(u_oce - u_ice)
                   zTauO = zaU(ji,jj) * zrhoco * 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) ) )
-
-                  ! Ocean-to-Ice stress
+                  !                 !--- Ocean-to-Ice stress
                   ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
-
-                  ! tau_bottom/u_ice
+                  !
+                  !                 !--- tau_bottom/u_ice
                   zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
                   zTauB = - tau_icebfr(ji,jj) / zvel
-
-                  ! Coriolis at U-points (energy conserving formulation)
+                  !
+                  !                 !--- Coriolis at U-points (energy conserving formulation)
                   zCorx(ji,jj)  =   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
+                  !
+                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
                   zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
-
-                  ! landfast switch => 0 = static friction ; 1 = sliding friction
+                  !
+                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
-
-                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  !
+                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
+                  u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
+                     &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0
+                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin 
+                     &            ) * zmaskU(ji,jj)
+                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
                   u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                             &  ! previous velocity
                      &                                     + zTauE + zTauO * u_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
@@ -514 +540 @@
                      &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
                      &            ) * zmaskU(ji,jj)
-                  ! Bouillon 2013
-                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
-                  !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
-                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
+                  ENDIF
                END DO
             END DO
@@ -532 +555 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
-
-                  ! tau_io/(u_oce - u_ice)
+                  !                 !--- tau_io/(u_oce - u_ice)
                   zTauO = zaU(ji,jj) * zrhoco * 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) ) )
-
-                  ! Ocean-to-Ice stress
+                  !                 !--- Ocean-to-Ice stress
                   ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
-
-                  ! tau_bottom/u_ice
+                  !
+                  !                 !--- tau_bottom/u_ice
                   zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
                   zTauB = - tau_icebfr(ji,jj) / zvel
-
-                  ! Coriolis at U-points (energy conserving formulation)
+                  !
+                  !                 !--- Coriolis at U-points (energy conserving formulation)
                   zCorx(ji,jj)  =   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
+                  !
+                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
                   zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
-
-                  ! landfast switch => 0 = static friction ; 1 = sliding friction
+                  !
+                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
-
-                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  !
+                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
+                  u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
+                     &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0
+                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin 
+                     &            ) * zmaskU(ji,jj)
+                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
                   u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                             &  ! previous velocity
                      &                                     + zTauE + zTauO * u_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                      &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
                      &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
-                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
+                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
                      &            ) * zmaskU(ji,jj)
-                  ! Bouillon 2013
-                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
-                  !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
-                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
+                  ENDIF
                END DO
             END DO
@@ -578 +603 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
-         
-                  ! tau_io/(v_oce - v_ice)
+                  !                 !--- tau_io/(v_oce - v_ice)
                   zTauO = zaV(ji,jj) * zrhoco * 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) ) )
-
-                  ! Ocean-to-Ice stress
+                  !                 !--- Ocean-to-Ice stress
                   ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
-
-                  ! tau_bottom/v_ice
+                  !
+                  !                 !--- tau_bottom/v_ice
                   zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
                   ztauB = - tau_icebfr(ji,jj) / zvel
-                  
-                  ! Coriolis at V-points (energy conserving formulation)
+                  !
+                  !                 !--- Coriolis at v-points (energy conserving formulation)
                   zCory(ji,jj)  = - 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
+                  !
+                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
                   zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
-
-                  ! landfast switch => 0 = static friction (tau_icebfr > zTauE); 1 = sliding friction
+                  !
+                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
-                  
-                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  !
+                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
+                  v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
+                     &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                     &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
+                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin
+                     &           ) * zmaskV(ji,jj)
+                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
                   v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             &  ! previous velocity
                      &                                     + zTauE + zTauO * v_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
@@ -608 +638 @@
                      &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
                      &            ) * zmaskV(ji,jj)
-                  ! Bouillon 2013
-                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
-                  !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
-                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
+                  ENDIF
                END DO
             END DO
@@ -820 +847 @@
    END SUBROUTINE ice_dyn_rhg_evp
 
+
    SUBROUTINE rhg_evp_rst( cdrw, kt )
       !!---------------------------------------------------------------------

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90

@@ -99 +99 @@
       !!
       !! ** Action  : It provides the following fields used in sea ice model:
-      !!                fr1_i0  , fr2_i0                         = 1sr & 2nd fraction of qsr penetration in ice  [%]
       !!                emp_oce , emp_ice                        = E-P over ocean and sea ice                    [Kg/m2/s]
       !!                sprecip                                  = solid precipitation                           [Kg/m2/s]
@@ -142 +141 @@
          !
       CASE( jp_blk )                !--- bulk formulation
-                                CALL blk_ice_flx( t_su, alb_ice )    ! 
-         IF( ln_mixcpl      )   CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su )
+                                CALL blk_ice_flx( t_su, h_s, h_i, alb_ice )    ! 
+         IF( ln_mixcpl      )   CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )
          IF( nn_iceflx /= 2 )   CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_iceflx )
          !
       CASE ( jp_purecpl )           !--- coupled formulation
-                                CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su )
+                                CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )
          IF( nn_iceflx == 2 )   CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_iceflx )
          !
@@ -268 +267 @@
       INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
       !!
-      NAMELIST/namforcing/ rn_cio, rn_blow_s, nn_iceflx
+      NAMELIST/namforcing/ rn_cio, rn_blow_s, nn_iceflx, nice_jules
       !!-------------------------------------------------------------------
       !
@@ -285 +284 @@
          WRITE(numout,*) '~~~~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namforcing:'
-         WRITE(numout,*) '      drag coefficient for oceanic stress              rn_cio    = ', rn_cio
-         WRITE(numout,*) '      coefficient for ice-lead partition of snowfall   rn_blow_s = ', rn_blow_s
-         WRITE(numout,*) '      Multicategory heat flux formulation              nn_iceflx = ', nn_iceflx
+         WRITE(numout,*) '      drag coefficient for oceanic stress              rn_cio     = ', rn_cio
+         WRITE(numout,*) '      coefficient for ice-lead partition of snowfall   rn_blow_s  = ', rn_blow_s
+         WRITE(numout,*) '      Multicategory heat flux formulation              nn_iceflx  = ', nn_iceflx
+         WRITE(numout,*) '      Jules coupling (0=no, 1=emulated, 2=active)      nice_jules = ', nice_jules
       ENDIF
       !

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90

@@ -48 +48 @@
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si  ! structure of input fields (file informations, fields read)
    !
-   ! ** namelist (namini) **
+   !                             !! ** namelist (namini) **
    LOGICAL  ::   ln_iceini        ! initialization or not
    LOGICAL  ::   ln_iceini_file   ! Ice initialization state from 2D netcdf file
@@ -91 +91 @@
       !!              where there is no ice (clem: I do not know why, is it mandatory?) 
       !!--------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk, jl             ! dummy loop indices
+      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices
+      INTEGER  ::   i_hemis, i_fill, jl0   ! local integers
       REAL(wp) ::   ztmelts, zdh
-      INTEGER  ::   i_hemis, i_fill, jl0
       REAL(wp) ::   zarg, zV, zconv, zdv, zfac
       INTEGER , DIMENSION(4)           ::   itest
@@ -100 +100 @@
       REAL(wp), DIMENSION(jpi,jpj)     ::   zht_i_ini, zat_i_ini, zvt_i_ini            !data from namelist or nc file
       REAL(wp), DIMENSION(jpi,jpj)     ::   zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zh_i_ini, za_i_ini                         !data by cattegories to fill
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zh_i_ini , za_i_ini                        !data by cattegories to fill
       !--------------------------------------------------------------------
 
@@ -116 +116 @@
          tn_ice(:,:,jl) = rt0 * tmask(:,:,1)
       END DO
-
-      ! init basal temperature (considered at freezing point)
+      !
+      ! init basal temperature (considered at freezing point)   [Kelvin]
       CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
       t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 
@@ -142 +142 @@
             zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:)
             !
-         !                             !---------------!
+            !                          !---------------!
          ELSE                          ! Read namelist !
             !                          !---------------!
-
-           ! no ice if sst <= t-freez + ttest
-            WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp 
-            ELSEWHERE                                                                  ; zswitch(:,:) = tmask(:,:,1)
+            ! no ice if sst <= t-freez + ttest
+            WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst )   ;   zswitch(:,:) = 0._wp 
+            ELSEWHERE                                                                    ;   zswitch(:,:) = tmask(:,:,1)
             END WHERE
-
+            !
             ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
             WHERE( ff_t(:,:) >= 0._wp )
@@ -242 +241 @@
                         !
                      ENDIF
-
+                     !
                      ! Compatibility tests
                      zconv = ABS( zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:jpl) ) )           ! Test 1: area conservation
                      IF ( zconv < epsi06 ) itest(1) = 1
-                     
+                     !
                      zconv = ABS(       zat_i_ini(ji,jj)       * zht_i_ini(ji,jj)   &         ! Test 2: volume conservation
                         &        - SUM( za_i_ini (ji,jj,1:jpl) * zh_i_ini (ji,jj,1:jpl) ) )
                      IF ( zconv < epsi06 ) itest(2) = 1
-                     
+                     !
                      IF ( zh_i_ini(ji,jj,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1           ! Test 3: thickness of the last category is in-bounds ?
-                     
+                     !
                      itest(4) = 1
                      DO jl = 1, i_fill
@@ -260 +259 @@
                   END DO                                                        ! end iteration on categories
                   !                                                             !----------------------------
-                  !
                   IF( lwp .AND. SUM(itest) /= 4 ) THEN 
                      WRITE(numout,*)
@@ -270 +268 @@
                      WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj)
                   ENDIF
-               
+                  !
                ENDIF
                !
@@ -279 +277 @@
          ! 4) Fill in sea ice arrays
          !---------------------------------------------------------------------
-
+         !
          ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature
          DO jl = 1, jpl ! loop over categories
@@ -289 +287 @@
                   o_i(ji,jj,jl)  = 0._wp                                                     ! age (0 day)
                   t_su(ji,jj,jl) = zswitch(ji,jj) * zts_u_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp
-
+                  !
                   IF( zht_i_ini(ji,jj) > 0._wp )THEN
                     h_s(ji,jj,jl)= h_i(ji,jj,jl) * ( zht_s_ini(ji,jj) / zht_i_ini(ji,jj) )  ! snow depth
@@ -295 +293 @@
                     h_s(ji,jj,jl)= 0._wp
                   ENDIF
-
+                  !
                   ! This case below should not be used if (h_s/h_i) is ok in namelist
                   ! In case snow load is in excess that would lead to transformation from snow to ice
@@ -303 +301 @@
                   h_i(ji,jj,jl) = MIN( hi_max(jl), h_i(ji,jj,jl) + zdh )
                   h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoic * r1_rhosn )
-
+                  !
                   ! ice volume, salt content, age content
                   v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume
@@ -312 +310 @@
             END DO
          END DO
-
-         ! for constant salinity in time
-         IF( nn_icesal /= 2 )  THEN
+         !
+         IF( nn_icesal /= 2 )  THEN         ! for constant salinity in time
             CALL ice_var_salprof
             sv_i = s_i * v_i
          ENDIF
-            
+         !  
          ! Snow temperature and heat content
          DO jk = 1, nlay_s
@@ -327 +324 @@
                      ! Snow energy of melting
                      e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )
-
+                     !
                      ! Mutliply by volume, and divide by number of layers to get heat content in J/m2
                      e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s
@@ -334 +331 @@
             END DO
          END DO
-
+         !
          ! Ice salinity, temperature and heat content
          DO jk = 1, nlay_i
@@ -343 +340 @@
                      sz_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rn_simin
                      ztmelts          = - tmut * sz_i(ji,jj,jk,jl) + rt0 !Melting temperature in K
-
+                     !
                      ! heat content per unit volume
-                     e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
-                        +   lfus    * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) &
-                        -   rcp     * ( ztmelts - rt0 ) )
-
+                     e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) )         &
+                        &             + lfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0) , -epsi20 )  )   &
+                        &             - rcp  * ( ztmelts - rt0 ) )
+                     !
                      ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
                      e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i
@@ -355 +352 @@
             END DO
          END DO
-
+         !
          tn_ice (:,:,:) = t_su (:,:,:)
-
+         !
          ! Melt pond volume and fraction
          IF ( ln_pnd_CST .OR. ln_pnd_H12 ) THEN   ;   zfac = 1._wp
@@ -368 +365 @@
          a_ip(:,:,:) = a_ip_frac(:,:,:) * a_i (:,:,:) 
          v_ip(:,:,:) = h_ip     (:,:,:) * a_ip(:,:,:)
-
+         !
       ELSE ! if ln_iceini=false
          a_i  (:,:,:) = 0._wp
@@ -379 +376 @@
          s_i  (:,:,:) = 0._wp
          o_i  (:,:,:) = 0._wp
-
+         !
          e_i(:,:,:,:) = 0._wp
          e_s(:,:,:,:) = 0._wp
-
+         !
          DO jl = 1, jpl
             DO jk = 1, nlay_i
@@ -391 +388 @@
             END DO
          END DO
-
+         !
          a_ip(:,:,:)      = 0._wp
          v_ip(:,:,:)      = 0._wp
          a_ip_frac(:,:,:) = 0._wp
          h_ip     (:,:,:) = 0._wp
-
+         !
       ENDIF ! ln_iceini
-      
+      !
       at_i (:,:) = 0.0_wp
       DO jl = 1, jpl
@@ -405 +402 @@
       !
       ! --- set ice velocities --- !
-      u_ice (:,:)     = 0._wp
-      v_ice (:,:)     = 0._wp
+      u_ice (:,:) = 0._wp
+      v_ice (:,:) = 0._wp
       !
       !----------------------------------------------
@@ -415 +412 @@
       !
       IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area
-
+         !
          sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
          sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
-
+         !
          IF( .NOT.ln_linssh ) THEN
-            
+            !
             WHERE( ht_0(:,:) > 0 )   ;   z2d(:,:) = 1._wp + sshn(:,:)*tmask(:,:,1) / ht_0(:,:)
             ELSEWHERE                ;   z2d(:,:) = 1._wp   ;   END WHERE
-         
+            !
             DO jk = 1,jpkm1                     ! adjust initial vertical scale factors                
                e3t_n(:,:,jk) = e3t_0(:,:,jk) * z2d(:,:)
@@ -429 +426 @@
                e3t_a(:,:,jk) = e3t_n(:,:,jk)
             END DO
-            
+            !
             ! Reconstruction of all vertical scale factors at now and before time-steps
             ! =========================================================================
@@ -478 +475 @@
 !!      CALL dia_wri_state( 'output.init', nit000 )
 !!!      
-
+      !
    END SUBROUTINE ice_istate
+
 
    SUBROUTINE ice_istate_init
@@ -485 +483 @@
       !!                   ***  ROUTINE ice_istate_init  ***
       !!        
-      !! ** Purpose : Definition of initial state of the ice 
-      !!
-      !! ** Method : Read the namini namelist and check the parameter 
-      !!       values called at the first timestep (nit000)
-      !!
-      !! ** input : 
-      !!        Namelist namini
-      !!
-      !! history :
-      !!  8.5  ! 03-08 (C. Ethe) original code 
-      !!  8.5  ! 07-11 (M. Vancoppenolle) rewritten initialization
+      !! ** Purpose :   Definition of initial state of the ice 
+      !!
+      !! ** Method  :   Read the namini namelist and check the parameter 
+      !!              values called at the first timestep (nit000)
+      !!
+      !! ** input   :  Namelist namini
+      !!
       !!-----------------------------------------------------------------------------
-      !
-      INTEGER :: ios,ierr,inum_ice                 ! Local integer output status for namelist read
-      INTEGER :: ji,jj
-      INTEGER :: ifpr, ierror
+      INTEGER ::   ji, jj
+      INTEGER ::   ios, ierr, inum_ice   ! Local integer output status for namelist read
+      INTEGER ::   ifpr, ierror
       !
       CHARACTER(len=256) ::  cn_dir          ! Root directory for location of ice files
@@ -515 +508 @@
       READ  ( numnam_ice_ref, namini, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namini in reference namelist', lwp )
-
+      !
       REWIND( numnam_ice_cfg )              ! Namelist namini in configuration namelist : Ice initial state
       READ  ( numnam_ice_cfg, namini, IOSTAT = ios, ERR = 902 )
 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namini in configuration namelist', lwp )
       IF(lwm) WRITE ( numoni, namini )
-
+      !
       slf_i(jp_hti) = sn_hti  ;  slf_i(jp_hts) = sn_hts
       slf_i(jp_ati) = sn_ati  ;  slf_i(jp_tsu) = sn_tsu
@@ -545 +538 @@
          WRITE(numout,*) '      initial  ice/snw temp  in the south                    rn_tmi_ini_s    = ', rn_tmi_ini_s
       ENDIF
-
+      !
       IF( ln_iceini_file ) THEN                      ! Ice initialization using input file
          !
@@ -553 +546 @@
             CALL ctl_stop( 'Ice_ini in iceistate: unable to allocate si structure' )   ;   RETURN
          ENDIF
-
+         !
          DO ifpr = 1, jpfldi
             ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
             ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
          END DO
-
+         !
          ! fill si with slf_i and control print
          CALL fld_fill( si, slf_i, cn_dir, 'ice_istate', 'ice istate ini', 'numnam_ice' )
-
+         !
          CALL fld_read( nit000, 1, si )                ! input fields provided at the current time-step
-
+         !
       ENDIF
-
+      !
    END SUBROUTINE ice_istate_init
 

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90

@@ -13 +13 @@
    !!   'key_lim3'                                       ESIM sea-ice model
    !!----------------------------------------------------------------------
-   !!   ice_itd_init   : read ice thicknesses mean and min from namelist
-   !!   ice_itd_rem    : redistribute ice thicknesses after thermo growth and melt
-   !!   ice_itd_reb    : rebin ice thicknesses into bounded categories
+   !!   ice_itd_rem   : redistribute ice thicknesses after thermo growth and melt
+   !!   itd_glinear   : build g(h) satisfying area and volume constraints
+   !!   itd_shiftice  : shift ice across category boundaries, conserving everything
+   !!   ice_itd_reb   : rebin ice thicknesses into bounded categories
+   !!   ice_itd_init  : read ice thicknesses mean and min from namelist
    !!----------------------------------------------------------------------
    USE dom_oce        ! ocean domain
@@ -36 +38 @@
    PUBLIC   ice_itd_reb   ! called in icecor
 
-   ! ** namelist (namitd) **
-   REAL(wp) ::   rn_himean   ! mean thickness of the domain
-
+   INTEGER            ::   nice_catbnd     ! choice of the type of ice category function
+   !                                       ! associated indices:
+   INTEGER, PARAMETER ::   np_cathfn = 1   ! categories defined by a function
+   INTEGER, PARAMETER ::   np_catusr = 2   ! categories defined by the user
+   !
+   !                                           !! ** namelist (namitd) **
+   LOGICAL                    ::   ln_cat_hfn   ! ice categories are defined by function like rn_himean**(-0.05)
+   REAL(wp)                   ::   rn_himean    ! mean thickness of the domain
+   LOGICAL                    ::   ln_cat_usr   ! ice categories are defined by rn_catbnd
+   REAL(wp), DIMENSION(0:100) ::   rn_catbnd    ! ice categories bounds
+   !
    !!----------------------------------------------------------------------
    !! NEMO/ICE 4.0 , NEMO Consortium (2017)
@@ -65 +75 @@
       REAL(wp) ::   zx3        
       REAL(wp) ::   zslope          ! used to compute local thermodynamic "speeds"
-
+      !
       INTEGER , DIMENSION(jpij)       ::   iptidx          ! compute remapping or not
       INTEGER , DIMENSION(jpij,jpl-1) ::   jdonor          ! donor category index
@@ -97 +107 @@
       !-----------------------------------------------------------------------------------------------
       IF( npti > 0 ) THEN
-
+         !
          zdhice(:,:) = 0._wp
          zhbnew(:,:) = 0._wp
-
+         !
          CALL tab_3d_2d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i   )
          CALL tab_3d_2d( npti, nptidx(1:npti), h_ib_2d(1:npti,1:jpl), h_i_b )
          CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i    )
          CALL tab_3d_2d( npti, nptidx(1:npti), a_ib_2d (1:npti,1:jpl), a_i_b  )
-
+         !
          DO jl = 1, jpl
             ! Compute thickness change in each ice category
@@ -112 +122 @@
             END DO
          END DO
-         
+         !
          ! --- New boundaries for category 1:jpl-1 --- !
          DO jl = 1, jpl - 1
@@ -135 +145 @@
                ! 1) hn(t+1)+espi < Hn* < hn+1(t+1)-epsi               
                !    Note: hn(t+1) must not be too close to either HR or HL otherwise a division by nearly 0 is possible 
-               !          in ice_itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
+               !          in itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
                IF( a_i_2d(ji,jl  ) > epsi10 .AND. h_i_2d(ji,jl  ) > ( zhbnew(ji,jl) - epsi10 ) )   nptidx(ji) = 0
                IF( a_i_2d(ji,jl+1) > epsi10 .AND. h_i_2d(ji,jl+1) < ( zhbnew(ji,jl) + epsi10 ) )   nptidx(ji) = 0
-               
+               !
                ! 2) Hn-1 < Hn* < Hn+1  
                IF( zhbnew(ji,jl) < hi_max(jl-1) )   nptidx(ji) = 0
                IF( zhbnew(ji,jl) > hi_max(jl+1) )   nptidx(ji) = 0
-               
+               !
             END DO
          END DO
@@ -153 +163 @@
                zhbnew(ji,jpl) = hi_max(jpl)  
             ENDIF
-            
+            !
             ! --- 1 additional condition for remapping (1st category) --- !
             ! H0+epsi < h1(t) < H1-epsi 
             !    h1(t) must not be too close to either HR or HL otherwise a division by nearly 0 is possible 
-            !    in ice_itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
+            !    in itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
             IF( h_ib_2d(ji,1) < ( hi_max(0) + epsi10 ) )   nptidx(ji) = 0
             IF( h_ib_2d(ji,1) > ( hi_max(1) - epsi10 ) )   nptidx(ji) = 0
@@ -165 +175 @@
          !  3) Identify cells where remapping
          !-----------------------------------------------------------------------------------------------
-         ipti = 0 ; iptidx(:) = 0
+         ipti = 0   ;   iptidx(:) = 0
          DO ji = 1, npti
             IF( nptidx(ji) /= 0 ) THEN
@@ -197 +207 @@
                !  
                ! --- g(h) for category 1 --- !
-               CALL ice_itd_glinear( zhb0(1:npti)  , zhb1(1:npti)  , h_ib_1d(1:npti)  , a_i_1d(1:npti)  ,  &   ! in
-                  &                  g0  (1:npti,1), g1  (1:npti,1), hL      (1:npti,1), hR    (1:npti,1)   )   ! out
+               CALL itd_glinear( zhb0(1:npti)  , zhb1(1:npti)  , h_ib_1d(1:npti)  , a_i_1d(1:npti)  ,  &   ! in
+                  &              g0  (1:npti,1), g1  (1:npti,1), hL      (1:npti,1), hR    (1:npti,1)   )   ! out
                   !
                ! Area lost due to melting of thin ice
@@ -240 +250 @@
             !
             ! --- g(h) for each thickness category --- !  
-            CALL ice_itd_glinear( zhbnew(1:npti,jl-1), zhbnew(1:npti,jl), h_i_1d(1:npti)   , a_i_1d(1:npti)   ,  &   ! in
-               &                  g0    (1:npti,jl  ), g1    (1:npti,jl), hL     (1:npti,jl), hR    (1:npti,jl)   )   ! out
+            CALL itd_glinear( zhbnew(1:npti,jl-1), zhbnew(1:npti,jl), h_i_1d(1:npti)   , a_i_1d(1:npti)   ,  &   ! in
+               &              g0    (1:npti,jl  ), g1    (1:npti,jl), hL     (1:npti,jl), hR    (1:npti,jl)   )   ! out
             !
          END DO
@@ -281 +291 @@
          ! 6) Shift ice between categories
          !----------------------------------------------------------------------------------------------
-         CALL ice_itd_shiftice ( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )
+         CALL itd_shiftice ( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )
          
          !----------------------------------------------------------------------------------------------
@@ -289 +299 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,1)    )
          CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip(:,:,1)   )
-         
+         !
          DO ji = 1, npti
             IF ( a_i_1d(ji) > epsi10 .AND. h_i_1d(ji) < rn_himin ) THEN
@@ -309 +319 @@
 
 
-   SUBROUTINE ice_itd_glinear( HbL, Hbr, phice, paice, pg0, pg1, phL, phR )
-      !!------------------------------------------------------------------
-      !!                ***  ROUTINE ice_itd_glinear ***
+   SUBROUTINE itd_glinear( HbL, Hbr, phice, paice, pg0, pg1, phL, phR )
+      !!------------------------------------------------------------------
+      !!                ***  ROUTINE itd_glinear ***
       !!
       !! ** Purpose :   build g(h) satisfying area and volume constraints (Eq. 6 and 9)
@@ -367 +377 @@
       END DO
       !
-   END SUBROUTINE ice_itd_glinear
-
-
-   SUBROUTINE ice_itd_shiftice( kdonor, pdaice, pdvice )
-      !!------------------------------------------------------------------
-      !!                ***  ROUTINE ice_itd_shiftice ***
+   END SUBROUTINE itd_glinear
+
+
+   SUBROUTINE itd_shiftice( kdonor, pdaice, pdvice )
+      !!------------------------------------------------------------------
+      !!                ***  ROUTINE itd_shiftice ***
       !!
       !! ** Purpose :   shift ice across category boundaries, conserving everything
@@ -486 +496 @@
             END DO
          END DO
-
+         !
          DO jk = 1, nlay_i         !--- Ice heat content
             DO ji = 1, npti
@@ -529 +539 @@
       CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su )
       !
-   END SUBROUTINE ice_itd_shiftice
+   END SUBROUTINE itd_shiftice
    
 
@@ -550 +560 @@
       !
       IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution' 
-
+      !
       jdonor(:,:) = 0
       zdaice(:,:) = 0._wp
@@ -587 +597 @@
          !
          IF( npti > 0 ) THEN
-            CALL ice_itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl=>jl+1
+            CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl=>jl+1
             ! Reset shift parameters
             jdonor(1:npti,jl) = 0
@@ -618 +628 @@
          !
          IF( npti > 0 ) THEN
-            CALL ice_itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl+1=>jl
+            CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl+1=>jl
             ! Reset shift parameters
             jdonor(1:npti,jl) = 0
@@ -629 +639 @@
    END SUBROUTINE ice_itd_reb
 
+
    SUBROUTINE ice_itd_init
       !!------------------------------------------------------------------
@@ -637 +648 @@
       !! ** input   :   Namelist namitd
       !!-------------------------------------------------------------------
-      INTEGER  ::   jl    ! dummy loop index
-      INTEGER  ::   ios   ! Local integer output status for namelist read
+      INTEGER  ::   jl            ! dummy loop index
+      INTEGER  ::   ios, ioptio   ! Local integer output status for namelist read
       REAL(wp) ::   zhmax, znum, zden, zalpha   !   -      -
-      !!
-      NAMELIST/namitd/ rn_himean, rn_himin
+      !
+      NAMELIST/namitd/ ln_cat_hfn, rn_himean, ln_cat_usr, rn_catbnd, rn_himin
       !!------------------------------------------------------------------
       !
@@ -647 +658 @@
       READ  ( numnam_ice_ref, namitd, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namitd in reference namelist', lwp )
-
+      !
       REWIND( numnam_ice_cfg )      ! Namelist namitd in configuration namelist : Parameters for ice
       READ  ( numnam_ice_cfg, namitd, IOSTAT = ios, ERR = 902 )
@@ -658 +669 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namitd: '
-         WRITE(numout,*) '      mean ice thickness in the domain               rn_himean = ', rn_himean
-         WRITE(numout,*) '      minimum ice thickness                          rn_himin  = ', rn_himin 
+         WRITE(numout,*) '      Ice categories are defined by a function of rn_himean**(-0.05)    ln_cat_hfn = ', ln_cat_hfn
+         WRITE(numout,*) '         mean ice thickness in the domain                               rn_himean  = ', rn_himean
+         WRITE(numout,*) '      Ice categories are defined by rn_catbnd                           ln_cat_usr = ', ln_cat_usr
+         WRITE(numout,*) '         ice categories boundaries (m)                                  rn_catbnd  = ', rn_catbnd 
+         WRITE(numout,*) '      minimum ice thickness                                             rn_himin   = ', rn_himin 
       ENDIF
       !
@@ -665 +679 @@
       !  Thickness categories boundaries  !
       !-----------------------------------!
-      !
-      zalpha = 0.05_wp              ! max of each category (from h^(-alpha) function)
-      zhmax  = 3._wp * rn_himean
-      DO jl = 1, jpl
-         znum = jpl * ( zhmax+1 )**zalpha
-         zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp )
-         hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
-      END DO
+      !                             !== set the choice of ice categories ==!
+      ioptio = 0 
+      IF( ln_cat_hfn ) THEN   ;   ioptio = ioptio + 1   ;   nice_catbnd = np_cathfn    ;   ENDIF
+      IF( ln_cat_usr ) THEN   ;   ioptio = ioptio + 1   ;   nice_catbnd = np_catusr    ;   ENDIF
+      IF( ioptio /= 1 )   CALL ctl_stop( 'ice_itd_init: choose one and only one ice categories boundaries' )
+      !
+      SELECT CASE( nice_catbnd )
+      !                                !------------------------!
+      CASE( np_cathfn )                ! h^(-alpha) function
+         !                             !------------------------!
+         zalpha = 0.05_wp
+         zhmax  = 3._wp * rn_himean
+         DO jl = 1, jpl
+            znum = jpl * ( zhmax+1 )**zalpha
+            zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp )
+            hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
+         END DO
+         !                             !------------------------!
+      CASE( np_catusr )                ! user defined
+         !                             !------------------------!
+         DO jl = 0, jpl
+            hi_max(jl) = rn_catbnd(jl)
+         END DO
+         !
+      END SELECT
       !
       DO jl = 1, jpl                ! mean thickness by category

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90

@@ -153 +153 @@
          !------------------------------------------------------!
          ! It provides the following fields used in sea ice model:
-         !    fr1_i0  , fr2_i0     = 1sr & 2nd fraction of qsr penetration in ice  [%]
          !    emp_oce , emp_ice    = E-P over ocean and sea ice                    [Kg/m2/s]
          !    sprecip              = solid precipitation                           [Kg/m2/s]
@@ -320 +319 @@
          IF(lwp) WRITE(numout,*) '   nn_monocat forced to 0 as jpl>1, i.e. multi-category case is chosen'
       ENDIF
-      IF ( jpl == 1 .AND. nn_monocat == 0 ) THEN
-         CALL ctl_stop( 'STOP', 'par_init : if jpl=1 then nn_monocat should be between 1 and 4' )
-      ENDIF
+!     IF ( jpl == 1 .AND. nn_monocat == 0 ) THEN
+!        CALL ctl_stop( 'STOP', 'par_init : if jpl=1 then nn_monocat should be between 1 and 4' )
+!     ENDIF
       !
       IF( ln_bdy .AND. ln_icediachk )   CALL ctl_warn('par_init: online conservation check does not work with BDY')

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90

@@ -26 +26 @@
    USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, qns_tot, qsr_tot, sprecip, ln_cpl
    USE sbc_ice , ONLY : qsr_oce, qns_oce, qemp_oce, qsr_ice, qns_ice, dqns_ice, evap_ice, qprec_ice, qevap_ice, &
-      &                 fr1_i0, fr2_i0
+      &                 qml_ice, qcn_ice, qsr_ice_tr
    USE ice1D          ! sea-ice: thermodynamics variables
    USE icethd_zdf     ! sea-ice: vertical heat diffusion
@@ -383 +383 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d(1:npti), qprec_ice        )
          CALL tab_2d_1d( npti, nptidx(1:npti), qsr_ice_1d  (1:npti), qsr_ice (:,:,kl) )
-         CALL tab_2d_1d( npti, nptidx(1:npti), fr1_i0_1d   (1:npti), fr1_i0           )
-         CALL tab_2d_1d( npti, nptidx(1:npti), fr2_i0_1d   (1:npti), fr2_i0           )
          CALL tab_2d_1d( npti, nptidx(1:npti), qns_ice_1d  (1:npti), qns_ice (:,:,kl) )
          CALL tab_2d_1d( npti, nptidx(1:npti), ftr_ice_1d  (1:npti), ftr_ice (:,:,kl) )
@@ -393 +391 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), fhtur_1d    (1:npti), fhtur            )
          CALL tab_2d_1d( npti, nptidx(1:npti), fhld_1d     (1:npti), fhld             )
+         
+         CALL tab_2d_1d( npti, nptidx(1:npti), qml_ice_1d   (1:npti), qml_ice      (:,:,kl)  )
+         CALL tab_2d_1d( npti, nptidx(1:npti), qcn_ice_1d   (1:npti), qcn_ice      (:,:,kl) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), qsr_ice_tr_1d(1:npti), qsr_ice_tr   (:,:,kl) )
          !
          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni   )
@@ -498 +500 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr        )
          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam        )
-         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd        )
          !
          CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog        )

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90

@@ -130 +130 @@
       !------------------------------------------------------------------------------!
       !
-      DO ji = 1, npti
-         zdum       = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 
+      SELECT CASE ( nice_jules )
+      !
+      CASE ( np_jules_ACTIVE )
+         !
+         DO ji = 1, npti
+            zq_su(ji)        =   MAX( 0._wp, qml_ice_1d(ji) * rdt_ice )
+         END DO
+         !
+      CASE ( np_jules_EMULE )
+         !
+         DO ji = 1, npti
+            zdum             =   ( qns_ice_1d(ji) + qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) - fc_su(ji) ) 
+            qml_ice_1d(ji)   =   zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
+            zq_su(ji)        =   MAX( 0._wp, qml_ice_1d(ji) * rdt_ice )
+         END DO
+         !
+      CASE ( np_jules_OFF ) 
+         !
+         DO ji = 1, npti
+            zdum             =   ( qns_ice_1d(ji) + qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) - fc_su(ji) ) 
+            zdum             =   zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
+            zq_su(ji)        =   MAX( 0._wp, zdum             * rdt_ice )
+         END DO
+         !
+      END SELECT
+      !
+      DO ji = 1, npti
          zf_tt(ji)  = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 
-
-         zq_su (ji) = MAX( 0._wp, zdum      * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
          zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice )
       END DO

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90

@@ -15 +15 @@
    !!   'key_lim3'                                       ESIM sea-ice model
    !!----------------------------------------------------------------------
+   !!  ice_thd_zdf      : select the appropriate routine for vertical heat diffusion calculation
+   !!  ice_thd_zdf_BL99 : 
+   !!  ice_thd_enmelt   :
+   !!  ice_thd_zdf_init :
+   !!----------------------------------------------------------------------
    USE dom_oce        ! ocean space and time domain
    USE phycst         ! physical constants (ocean directory) 
@@ -34 +39 @@
    LOGICAL  ::   ln_cndi_U64      ! thermal conductivity: Untersteiner (1964)
    LOGICAL  ::   ln_cndi_P07      ! thermal conductivity: Pringle et al (2007)
-   REAL(wp) ::   rn_cnd_s         ! thermal conductivity of the snow [W/m/K]
    REAL(wp) ::   rn_kappa_i       ! coef. for the extinction of radiation Grenfell et al. (2006) [1/m]
-   LOGICAL  ::   ln_dqns_i        ! change non-solar surface flux with changing surface temperature (T) or not (F)
-
+   REAL(wp), PUBLIC ::   rn_cnd_s ! thermal conductivity of the snow [W/m/K]
+
+   INTEGER  ::   nice_zdf         ! Choice of the type of vertical heat diffusion formulation
+   !                                           ! associated indices:
+   INTEGER, PARAMETER ::   np_BL99    = 1      ! Bitz and Lipscomb (1999)
+
+   INTEGER, PARAMETER ::   np_zdf_jules_OFF    = 0   !   compute all temperatures from qsr and qns
+   INTEGER, PARAMETER ::   np_zdf_jules_SND    = 1   !   compute conductive heat flux and surface temperature from qsr and qns
+   INTEGER, PARAMETER ::   np_zdf_jules_RCV    = 2   !   compute snow and inner ice temperatures from qcnd
+   
    !!----------------------------------------------------------------------
    !! NEMO/ICE 4.0 , NEMO Consortium (2017)
@@ -48 +60 @@
       !!-------------------------------------------------------------------
       !!                ***  ROUTINE ice_thd_zdf  ***
-      !! ** Purpose :
-      !!           This routine determines the time evolution of snow and sea-ice 
-      !!           temperature profiles.
-      !! ** Method  :
-      !!           This is done by solving the heat equation diffusion with
-      !!           a Neumann boundary condition at the surface and a Dirichlet one
-      !!           at the bottom. Solar radiation is partially absorbed into the ice.
-      !!           The specific heat and thermal conductivities depend on ice salinity
-      !!           and temperature to take into account brine pocket melting. The 
-      !!           numerical
-      !!           scheme is an iterative Crank-Nicolson on a non-uniform multilayer grid 
-      !!           in the ice and snow system.
+      !!
+      !! ** Purpose :   select the appropriate routine for the computation
+      !!              of vertical diffusion
+      !!-------------------------------------------------------------------
+   
+      SELECT CASE ( nice_zdf )      ! Choose the vertical heat diffusion solver
+      !
+      !                             !-------------      
+      CASE( np_BL99 )               ! BL99 solver
+         !                          !-------------
+         SELECT CASE( nice_jules )
+         CASE( np_jules_OFF    )   ;   CALL ice_thd_zdf_BL99 ( np_zdf_jules_OFF )   ! No Jules coupler     
+         CASE( np_jules_EMULE  )   ;   CALL ice_thd_zdf_BL99 ( np_zdf_jules_SND )   ! Jules coupler is emulated (send/recieve)
+                                       CALL ice_thd_zdf_BL99 ( np_zdf_jules_RCV )
+         CASE( np_jules_ACTIVE )   ;   CALL ice_thd_zdf_BL99 ( np_zdf_jules_RCV )   ! Jules coupler is active (receive only)
+         END SELECT
+      END SELECT
+      
+   END SUBROUTINE ice_thd_zdf
+   
+   
+   SUBROUTINE ice_thd_zdf_BL99( k_jules )
+      !!-------------------------------------------------------------------
+      !!                ***  ROUTINE ice_thd_zdf_BL99  ***
+      !!
+      !! ** Purpose :   computes the time evolution of snow and sea-ice temperature
+      !!              profiles, using the original Bitz and Lipscomb (1999) algorithm
+      !!
+      !! ** Method  :   solves the heat equation diffusion with a Neumann boundary
+      !!              condition at the surface and a Dirichlet one at the bottom. 
+      !!              Solar radiation is partially absorbed into the ice.
+      !!              The specific heat and thermal conductivities depend on ice 
+      !!              salinity and temperature to take into account brine pocket   
+      !!              melting. The numerical scheme is an iterative Crank-Nicolson
+      !!              on a non-uniform multilayer grid in the ice and snow system.
       !!
       !!           The successive steps of this routine are
@@ -83 +118 @@
       !!           total ice/snow thickness : h_i_1d, h_s_1d
       !!-------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   k_jules     ! Jules coupling (0=OFF, 1=RECEIVE, 2=SEND)
+      !
       INTEGER ::   ji, jk         ! spatial loop index
       INTEGER ::   jm             ! current reference number of equation
@@ -88 +125 @@
       INTEGER ::   iconv          ! number of iterations in iterative procedure
       INTEGER ::   iconv_max = 50 ! max number of iterations in iterative procedure
-      
+      !
       INTEGER, DIMENSION(jpij) ::   jm_min    ! reference number of top equation
       INTEGER, DIMENSION(jpij) ::   jm_max    ! reference number of bottom equation
-      
+      !
       REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system
       REAL(wp) ::   zg1       =  2._wp        !
@@ -101 +138 @@
       REAL(wp) ::   zdti_bnd  =  1.e-4_wp     ! maximal authorized error on temperature 
       REAL(wp) ::   ztmelt_i                  ! ice melting temperature
-      REAL(wp) ::   z1_hsu
       REAL(wp) ::   zdti_max                  ! current maximal error on temperature 
       REAL(wp) ::   zcpi                      ! Ice specific heat
       REAL(wp) ::   zhfx_err, zdq             ! diag errors on heat
       REAL(wp) ::   zfac                      ! dummy factor
-      
+      !
       REAL(wp), DIMENSION(jpij)     ::   isnow       ! switch for presence (1) or absence (0) of snow
       REAL(wp), DIMENSION(jpij)     ::   ztsub       ! surface temperature at previous iteration
       REAL(wp), DIMENSION(jpij)     ::   zh_i, z1_h_i ! ice layer thickness
       REAL(wp), DIMENSION(jpij)     ::   zh_s, z1_h_s ! snow layer thickness
-      REAL(wp), DIMENSION(jpij)     ::   zfsw        ! solar radiation absorbed at the surface
       REAL(wp), DIMENSION(jpij)     ::   zqns_ice_b  ! solar radiation absorbed at the surface
       REAL(wp), DIMENSION(jpij)     ::   zfnet       ! surface flux function
       REAL(wp), DIMENSION(jpij)     ::   zdqns_ice_b ! derivative of the surface flux function
-      REAL(wp), DIMENSION(jpij)     ::   zftrice     ! solar radiation transmitted through the ice
-      
+      !
+      REAL(wp), DIMENSION(jpij       )     ::   ztsuold     ! Old surface temperature in the ice
       REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztiold      ! Old temperature in the ice
       REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsold      ! Old temperature in the snow
@@ -136 +171 @@
       REAL(wp), DIMENSION(jpij)            ::   zq_ini      ! diag errors on heat
       REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat
-      
+      !
+      REAL(wp) ::   zfr1, zfr2, zfrqsr_tr_i   ! dummy factor
+      !
       ! Mono-category
       REAL(wp) :: zepsilon      ! determines thres. above which computation of G(h) is done
@@ -162 +199 @@
       ELSEWHERE                         ;   z1_h_i(1:npti) = 0._wp
       END WHERE
-
+      !
       WHERE( zh_s(1:npti) >= epsi10 )   ;   z1_h_s(1:npti) = 1._wp / zh_s(1:npti)
       ELSEWHERE                         ;   z1_h_s(1:npti) = 0._wp
       END WHERE
       !
-      ! temperatures
-      ztsub  (1:npti)   = t_su_1d(1:npti)  ! temperature at the previous iteration
+      ! Store initial temperatures and non solar heat fluxes
+      IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND  ) THEN ! OFF or SND mode
+         !
+         ztsub  (1:npti)       =   t_su_1d(1:npti)                          ! surface temperature at iteration n-1
+         ztsuold(1:npti)       =   t_su_1d(1:npti)                          ! surface temperature initial value
+         zdqns_ice_b(1:npti)   =   dqns_ice_1d(1:npti)                      ! derivative of incoming nonsolar flux 
+         zqns_ice_b (1:npti)   =   qns_ice_1d(1:npti)                       ! store previous qns_ice_1d value
+         !
+         t_su_1d(1:npti)       =   MIN( t_su_1d(1:npti), rt0 - ztsu_err )   ! required to leave the choice between melting or not
+         !
+      ENDIF   
+      !
       ztsold (1:npti,:) = t_s_1d(1:npti,:)   ! Old snow temperature
       ztiold (1:npti,:) = t_i_1d(1:npti,:)   ! Old ice temperature
-      t_su_1d(1:npti)   = MIN( t_su_1d(1:npti), rt0 - ztsu_err )   ! necessary
-      !
+
       !-------------
       ! 2) Radiation
       !-------------
-      z1_hsu = 1._wp / 0.1_wp ! threshold for the computation of i0
-      DO ji = 1, npti
-         ! --- Computation of i0 --- !
-         ! i0 describes the fraction of solar radiation which does not contribute
-         ! to the surface energy budget but rather penetrates inside the ice.
-         ! We assume that no radiation is transmitted through the snow
-         ! If there is no no snow
-         ! zfsw    = (1-i0).qsr_ice   is absorbed at the surface 
-         ! zftrice = io.qsr_ice       is below the surface 
-         ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice 
-         ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover
-         zfac = MAX( 0._wp , 1._wp - ( h_i_1d(ji) * z1_hsu ) )     
-         i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zfac * fr2_i0_1d(ji) )
-
-         ! --- Solar radiation absorbed / transmitted at the surface --- !
-         !     Derivative of the non solar flux
-         zfsw   (ji)     =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface
-         zftrice(ji)     =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer
-         zdqns_ice_b(ji) = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux 
-         zqns_ice_b (ji) =  qns_ice_1d(ji)                    ! store previous qns_ice_1d value
-      END DO
-
       ! --- Transmission/absorption of solar radiation in the ice --- !
-      zradtr_s(1:npti,0) = zftrice(1:npti)
+!     zfr1 = ( 0.18 * ( 1.0 - 0.81 ) + 0.35 * 0.81 )            !   standard value
+!     zfr2 = ( 0.82 * ( 1.0 - 0.81 ) + 0.65 * 0.81 )            !   zfr2 such that zfr1 + zfr2 to equal 1
+!
+!     DO ji = 1, npti
+!
+!        zfac = MAX( 0._wp , 1._wp - ( h_i_1d(ji) * 10._wp ) )     
+!
+!              zfrqsr_tr_i = zfr1 + zfac * zfr2                         !   below 10 cm, linearly increase zfrqsr_tr_i until 1 at zero thickness
+!              IF ( h_s_1d(ji) >= 0.0_wp )   zfrqsr_tr_i = 0._wp     !   snow fully opaque
+!
+!              qsr_ice_tr_1d(ji) = zfrqsr_tr_i * qsr_ice_1d(ji)   !   transmitted solar radiation 
+!
+!              zfsw(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji)
+!              zftrice(ji) = qsr_ice_tr_1d(ji)
+!              i0(ji) = zfrqsr_tr_i
+!
+!     END DO
+
+      zradtr_s(1:npti,0) = qsr_ice_tr_1d(1:npti)
       DO jk = 1, nlay_s
          DO ji = 1, npti
@@ -208 +250 @@
          END DO
       END DO
-         
-      zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + zftrice(1:npti) * ( 1._wp - isnow(1:npti) )
+      !
+      zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + qsr_ice_tr_1d(1:npti) * ( 1._wp - isnow(1:npti) )
       DO jk = 1, nlay_i 
          DO ji = 1, npti
@@ -218 +260 @@
          END DO
       END DO
-
+      !
       ftr_ice_1d(1:npti) = zradtr_i(1:npti,nlay_i)   ! record radiation transmitted below the ice
       !
@@ -273 +315 @@
          !
          SELECT CASE ( nn_monocat )
-
+         !
          CASE ( 1 , 3 )
-
+            !
             zepsilon = 0.1_wp
             DO ji = 1, npti
@@ -284 +326 @@
                ENDIF
             END DO
-
+            !
          END SELECT
          !
@@ -330 +372 @@
             END DO
          END DO
-         !
-         !----------------------------
-         ! 6) surface flux computation
-         !----------------------------
-         IF ( ln_dqns_i ) THEN 
-            DO ji = 1, npti
-               ! update of the non solar flux according to the update in T_su
+         
+         !
+         !----------------------------------------!
+         !                                        !
+         !       JULES IF (OFF or SND MODE)       !
+         !                                        !
+         !----------------------------------------!
+         !
+         
+         IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND  ) THEN ! OFF or SND mode
+         
+            !
+            ! In OFF mode the original BL99 temperature computation is used
+            ! (with qsr_ice, qns_ice and dqns_ice as inputs)
+            !
+            ! In SND mode, the computation is required to compute the conduction fluxes
+            !
+         
+            !----------------------------
+            ! 6) surface flux computation
+            !----------------------------
+
+            DO ji = 1, npti
+            ! update of the non solar flux according to the update in T_su
                qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) )
             END DO
-         ENDIF
-
-         DO ji = 1, npti
-            zfnet(ji) = zfsw(ji) + qns_ice_1d(ji) ! incoming = net absorbed solar radiation + non solar total flux (LWup, LWdw, SH, LH)
-         END DO
-         !
-         !----------------------------
-         ! 7) tridiagonal system terms
-         !----------------------------
-         !!layer denotes the number of the layer in the snow or in the ice
-         !!jm denotes the reference number of the equation in the tridiagonal
-         !!system, terms of tridiagonal system are indexed as following :
-         !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
-
-         !!ice interior terms (top equation has the same form as the others)
-         ztrid   (1:npti,:,:) = 0._wp
-         zindterm(1:npti,:)   = 0._wp
-         zindtbis(1:npti,:)   = 0._wp
-         zdiagbis(1:npti,:)   = 0._wp
-
-         DO jm = nlay_s + 2, nlay_s + nlay_i 
+
+            DO ji = 1, npti
+               zfnet(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) + qns_ice_1d(ji) ! net heat flux = net solar - transmitted solar + non solar
+            END DO
+            !
+            !----------------------------
+            ! 7) tridiagonal system terms
+            !----------------------------
+            !!layer denotes the number of the layer in the snow or in the ice
+            !!jm denotes the reference number of the equation in the tridiagonal
+            !!system, terms of tridiagonal system are indexed as following :
+            !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
+
+            !!ice interior terms (top equation has the same form as the others)
+            ztrid   (1:npti,:,:) = 0._wp
+            zindterm(1:npti,:)   = 0._wp
+            zindtbis(1:npti,:)   = 0._wp
+            zdiagbis(1:npti,:)   = 0._wp
+
+            DO jm = nlay_s + 2, nlay_s + nlay_i 
             DO ji = 1, npti
                jk = jm - nlay_s - 1
@@ -367 +425 @@
                zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
             END DO
-         ENDDO
-
-         jm =  nlay_s + nlay_i + 1
-         DO ji = 1, npti
+            END DO
+
+            jm =  nlay_s + nlay_i + 1
+            DO ji = 1, npti
             !!ice bottom term
             ztrid(ji,jm,1)  = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)   
@@ -377 +435 @@
             zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
                &            ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) ) 
-         ENDDO
-
-
-         DO ji = 1, npti
+            END DO
+
+
+            DO ji = 1, npti
             !                               !---------------------!
-            IF ( h_s_1d(ji) > 0.0 ) THEN   !  snow-covered cells !
+            IF ( h_s_1d(ji) > 0.0 ) THEN    !  snow-covered cells !
                !                            !---------------------!
                ! snow interior terms (bottom equation has the same form as the others)
@@ -428 +486 @@
                      &           ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
                ENDIF
-            !                               !---------------------!
+               !                            !---------------------!
             ELSE                            ! cells without snow  !
                !                            !---------------------!
@@ -453 +511 @@
                      ztrid(ji,jm_min(ji),1)    = 0.0
                      ztrid(ji,jm_min(ji),2)    = zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0
-                     ztrid(ji,jm_min(ji),3)    = zkappa_i(ji,0) * 2.0
+                     ztrid(ji,jm_min(ji),3)    =  zkappa_i(ji,0) * 2.0
                      ztrid(ji,jm_min(ji)+1,1)  = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1)
                      ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
@@ -488 +546 @@
             zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2)
             !
-         END DO
-         !
-         !------------------------------
-         ! 8) tridiagonal system solving
-         !------------------------------
-         ! Solve the tridiagonal system with Gauss elimination method.
-         ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
-         jm_maxt = 0
-         jm_mint = nlay_i+5
-         DO ji = 1, npti
-            jm_mint = MIN(jm_min(ji),jm_mint)
-            jm_maxt = MAX(jm_max(ji),jm_maxt)
-         END DO
-
-         DO jk = jm_mint+1, jm_maxt
+            END DO
+            !
+            !------------------------------
+            ! 8) tridiagonal system solving
+            !------------------------------
+            ! Solve the tridiagonal system with Gauss elimination method.
+            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
+            jm_maxt = 0
+            jm_mint = nlay_i+5
+            DO ji = 1, npti
+               jm_mint = MIN(jm_min(ji),jm_mint)
+               jm_maxt = MAX(jm_max(ji),jm_maxt)
+            END DO
+
+            DO jk = jm_mint+1, jm_maxt
             DO ji = 1, npti
                jm = min(max(jm_min(ji)+1,jk),jm_max(ji))
@@ -508 +566 @@
                zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1)
             END DO
-         END DO
-
-         DO ji = 1, npti
+            END DO
+
+            DO ji = 1, npti
             ! ice temperatures
             t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
-         END DO
-
-         DO jm = nlay_i + nlay_s, nlay_s + 2, -1
+            END DO
+
+            DO jm = nlay_i + nlay_s, nlay_s + 2, -1
             DO ji = 1, npti
                jk    =  jm - nlay_s - 1
                t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
             END DO
-         END DO
-
-         DO ji = 1, npti
+            END DO
+
+            DO ji = 1, npti
             ! snow temperatures      
             IF( h_s_1d(ji) > 0._wp ) THEN
                t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  &
-                  &                / zdiagbis(ji,nlay_s+1)
+               &                / zdiagbis(ji,nlay_s+1)
             ENDIF
             ! surface temperature
@@ -534 +592 @@
                   &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))
             ENDIF
-         END DO
-         !
-         !--------------------------------------------------------------
-         ! 9) Has the scheme converged ?, end of the iterative procedure
-         !--------------------------------------------------------------
-         ! check that nowhere it has started to melt
-         ! zdti_max is a measure of error, it has to be under zdti_bnd
-         zdti_max = 0._wp
-         DO ji = 1, npti
+            END DO
+            !
+            !--------------------------------------------------------------
+            ! 9) Has the scheme converged ?, end of the iterative procedure
+            !--------------------------------------------------------------
+            ! check that nowhere it has started to melt
+            ! zdti_max is a measure of error, it has to be under zdti_bnd
+            zdti_max = 0._wp
+            DO ji = 1, npti
             t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp )
             zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )     
-         END DO
-
-         DO jk  =  1, nlay_s
+            END DO
+
+            DO jk  =  1, nlay_s
             DO ji = 1, npti
                t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )
                zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )
             END DO
-         END DO
-
-         DO jk  =  1, nlay_i
+            END DO
+
+            DO jk  =  1, nlay_i
             DO ji = 1, npti
                ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0 
@@ -560 +618 @@
                zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
             END DO
-         END DO
-
-         ! Compute spatial maximum over all errors
-         ! note that this could be optimized substantially by iterating only the non-converging points
-         IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice )
-
+            END DO
+
+            ! Compute spatial maximum over all errors
+            ! note that this could be optimized substantially by iterating only the non-converging points
+            IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice )
+         !
+         !----------------------------------------!
+         !                                        !
+         !       JULES IF (RCV MODE)              !
+         !                                        !
+         !----------------------------------------!
+         !
+
+         ELSE IF ( k_jules == np_zdf_jules_RCV  ) THEN ! RCV mode
+         
+            !
+            ! In RCV mode, we use a modified BL99 solver 
+            ! with conduction flux (qcn_ice) as forcing term
+            !
+            !----------------------------
+            ! 7) tridiagonal system terms
+            !----------------------------
+            !!layer denotes the number of the layer in the snow or in the ice
+            !!jm denotes the reference number of the equation in the tridiagonal
+            !!system, terms of tridiagonal system are indexed as following :
+            !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
+
+            !!ice interior terms (top equation has the same form as the others)
+            ztrid   (1:npti,:,:) = 0._wp
+            zindterm(1:npti,:)   = 0._wp
+            zindtbis(1:npti,:)   = 0._wp
+            zdiagbis(1:npti,:)   = 0._wp
+
+            DO jm = nlay_s + 2, nlay_s + nlay_i 
+            DO ji = 1, npti
+               jk = jm - nlay_s - 1
+               ztrid(ji,jm,1)  = - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
+               ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
+               ztrid(ji,jm,3)  = - zeta_i(ji,jk) * zkappa_i(ji,jk)
+               zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
+            END DO
+            ENDDO
+
+            jm =  nlay_s + nlay_i + 1
+            DO ji = 1, npti
+            !!ice bottom term
+            ztrid(ji,jm,1)  = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)   
+            ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) )
+            ztrid(ji,jm,3)  = 0.0
+            zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
+               &            ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) ) 
+            ENDDO
+
+
+            DO ji = 1, npti
+            !                              !---------------------!
+            IF ( h_s_1d(ji) > 0.0 ) THEN   !  snow-covered cells !
+               !                           !---------------------!
+               ! snow interior terms (bottom equation has the same form as the others)
+               DO jm = 3, nlay_s + 1
+               jk = jm - 1
+               ztrid(ji,jm,1)  = - zeta_s(ji,jk) * zkappa_s(ji,jk-1)
+               ztrid(ji,jm,2)  = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
+               ztrid(ji,jm,3)  = - zeta_s(ji,jk)*zkappa_s(ji,jk)
+               zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
+               END DO
+
+               ! case of only one layer in the ice (ice equation is altered)
+               IF ( nlay_i == 1 ) THEN
+               ztrid(ji,nlay_s+2,3)  = 0.0
+               zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 
+               ENDIF
+
+               jm_min(ji) = 2
+               jm_max(ji) = nlay_i + nlay_s + 1
+
+               ! first layer of snow equation
+               ztrid(ji,2,1)  = 0.0
+               ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * zkappa_s(ji,1)
+               ztrid(ji,2,3)  = - zeta_s(ji,1)*zkappa_s(ji,1) 
+               zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) *  &
+               &           ( zradab_s(ji,1) + qcn_ice_1d(ji) ) 
+
+            !                               !---------------------!
+            ELSE                            ! cells without snow  !
+            !                               !---------------------!
+
+               jm_min(ji)    =  nlay_s + 2
+               jm_max(ji)    =  nlay_i + nlay_s + 1
+
+               ! first layer of ice equation
+               ztrid(ji,jm_min(ji),1)  = 0.0
+               ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1)  
+               ztrid(ji,jm_min(ji),3)  = - zeta_i(ji,1) * zkappa_i(ji,1)
+               zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) *  &
+               &                    ( zradab_i(ji,1) + qcn_ice_1d(ji) )
+
+               ! case of only one layer in the ice (surface & ice equations are altered)
+               IF ( nlay_i == 1 ) THEN
+               ztrid(ji,jm_min(ji),1)  = 0.0
+               ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1)
+               ztrid(ji,jm_min(ji),3)  = 0.0
+               zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji)    &
+               &                    + qcn_ice_1d(ji) )
+
+               ENDIF
+
+            ENDIF
+            !
+            zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
+            zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2)
+            !
+            END DO
+            !
+            !------------------------------
+            ! 8) tridiagonal system solving
+            !------------------------------
+            ! Solve the tridiagonal system with Gauss elimination method.
+            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
+            jm_maxt = 0
+            jm_mint = nlay_i+5
+            DO ji = 1, npti
+            jm_mint = MIN(jm_min(ji),jm_mint)
+            jm_maxt = MAX(jm_max(ji),jm_maxt)
+            END DO
+
+            DO jk = jm_mint+1, jm_maxt
+            DO ji = 1, npti
+               jm = min(max(jm_min(ji)+1,jk),jm_max(ji))
+               zdiagbis(ji,jm) = ztrid(ji,jm,2)  - ztrid(ji,jm,1) * ztrid(ji,jm-1,3)  / zdiagbis(ji,jm-1)
+               zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1)
+            END DO
+            END DO
+
+            DO ji = 1, npti
+            ! ice temperatures
+            t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
+            END DO
+
+            DO jm = nlay_i + nlay_s, nlay_s + 2, -1
+            DO ji = 1, npti
+               jk    =  jm - nlay_s - 1
+               t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
+            END DO
+            END DO
+
+            DO ji = 1, npti
+            ! snow temperatures      
+            IF( h_s_1d(ji) > 0._wp ) THEN
+               t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  &
+               &                / zdiagbis(ji,nlay_s+1)
+            ENDIF
+            END DO
+            !
+            !--------------------------------------------------------------
+            ! 9) Has the scheme converged ?, end of the iterative procedure
+            !--------------------------------------------------------------
+            ! check that nowhere it has started to melt
+            ! zdti_max is a measure of error, it has to be under zdti_bnd
+            zdti_max = 0._wp
+
+            DO jk  =  1, nlay_s
+            DO ji = 1, npti
+               t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )
+               zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )
+            END DO
+            END DO
+
+            DO jk  =  1, nlay_i
+            DO ji = 1, npti
+               ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0 
+               t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp )
+               zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
+            END DO
+            END DO
+
+            ! Compute spatial maximum over all errors
+            ! note that this could be optimized substantially by iterating only the non-converging points
+            IF( lk_mpp )   CALL mpp_max( zdti_max, kcom=ncomm_ice )
+
+         ENDIF ! k_jules
+         
       END DO  ! End of the do while iterative procedure
-
+      
       IF( ln_icectl .AND. lwp ) THEN
          WRITE(numout,*) ' zdti_max : ', zdti_max
          WRITE(numout,*) ' iconv    : ', iconv
       ENDIF
+      
       !
       !-----------------------------
       ! 10) Fluxes at the interfaces
       !-----------------------------
+      !
+      ! --- update conduction fluxes
+      !
       DO ji = 1, npti
-         !                                ! surface ice conduction flux
+      !                                ! surface ice conduction flux
          fc_su(ji)   =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   &
             &           - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji))
-         !                                ! bottom ice conduction flux
+      !                                ! bottom ice conduction flux
          fc_bo_i(ji) =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) )
       END DO
-
-      ! --- computes sea ice energy of melting compulsory for icethd_dh --- !
-      CALL ice_thd_enmelt
-
-      ! --- diagnose the change in non-solar flux due to surface temperature change --- !
-      IF ( ln_dqns_i ) THEN
+      
+      !
+      ! --- Diagnose the heat loss due to changing non-solar / conduction flux --- !
+      !
+      IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND  ) THEN  ! OFF or SND mode
          DO ji = 1, npti
-            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji)  - zqns_ice_b(ji) ) * a_i_1d(ji) 
-         END DO
-      END IF
-
-      ! --- diag conservation imbalance on heat diffusion - PART 2 --- !
-      !     hfx_dif = Heat flux used to warm/cool ice in W.m-2
-      !     zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
-      DO ji = 1, npti
-         zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  &
-            &                   SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s )
-
-         IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
-            zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 
-         ELSE                          ! case T_su = 0degC
-            zhfx_err = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji)
-         ENDIF
-         hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji)
-
-         ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation)
-         hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err
-
-      END DO   
-
-      ! --- Diagnostics SIMIP --- !
-      DO ji = 1, npti
-         !--- Conduction fluxes (positive downwards)
-         diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji)
-         diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su(ji)   * a_i_1d(ji) / at_i_1d(ji)
-
-         !--- Snow-ice interfacial temperature (diagnostic SIMIP)
-         zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji)
-         IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN
-            t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   &
-               &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac
-         ELSE
-            t_si_1d(ji) = t_su_1d(ji)
-         ENDIF
-      END DO
-      !
-   END SUBROUTINE ice_thd_zdf
+            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji) 
+         END DO
+      ELSE        ! RCV mode
+         DO ji = 1, npti
+            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( fc_su(ji)      - qcn_ice_1d(ji) ) * a_i_1d(ji) 
+         END DO
+      ENDIF
+      
+      !
+      ! --- Diagnose the heat loss due to non-fully converged temperature solution (should not be above 10-4 W-m2) --- !
+      !
+
+      IF ( ( k_jules == np_zdf_jules_OFF ) .OR. ( k_jules == np_zdf_jules_RCV ) ) THEN ! OFF
+      
+         CALL ice_thd_enmelt       
+
+         !     zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
+         DO ji = 1, npti
+            zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  &
+               &                   SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s )
+               
+            IF ( ( k_jules == np_zdf_jules_OFF ) ) THEN
+               
+                  IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
+                     zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji)    - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 
+                  ELSE                          ! case T_su = 0degC
+                     zhfx_err = ( fc_su(ji)      + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji)
+                  ENDIF
+            
+            ELSE ! RCV CASE
+            
+               zhfx_err = ( fc_su(ji) + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji)
+            
+            ENDIF
+            !
+            ! total heat sink to be sent to the ocean
+            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err
+            !
+            ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2   
+            hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji)
+            !
+         END DO
+         !
+         ! --- SIMIP diagnostics
+         !
+         DO ji = 1, npti
+            !--- Conduction fluxes (positive downwards)
+            diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji)
+            diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su(ji)   * a_i_1d(ji) / at_i_1d(ji)
+   
+            !--- Snow-ice interfacial temperature (diagnostic SIMIP)
+            zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji)
+            IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN
+               t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   &
+                  &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac
+            ELSE
+               t_si_1d(ji) = t_su_1d(ji)
+            ENDIF
+         END DO
+         !
+      ENDIF
+      !
+      !---------------------------------------------------------------------------------------
+      ! 11) Jules coupling: reset inner snow and ice temperatures, update conduction fluxes
+      !---------------------------------------------------------------------------------------
+      !
+      IF ( k_jules == np_zdf_jules_SND ) THEN   ! --- Jules coupling in "SND" mode
+         !
+         ! Restore temperatures to their initial values
+         t_s_1d(1:npti,:)        = ztsold (1:npti,:)
+         t_i_1d(1:npti,:)        = ztiold (1:npti,:)
+         qcn_ice_1d(1:npti)      = fc_su(1:npti)
+         !  
+      ENDIF
+      !
+   END SUBROUTINE ice_thd_zdf_BL99
 
 
@@ -675 +948 @@
       !! ** input   :   Namelist namthd_zdf
       !!-------------------------------------------------------------------
-      INTEGER  ::   ios   ! Local integer output status for namelist read
+      INTEGER  ::   ios, ioptio   ! Local integer
       !!
-      NAMELIST/namthd_zdf/ ln_zdf_BL99, ln_cndi_U64, ln_cndi_P07, rn_cnd_s, rn_kappa_i, ln_dqns_i
+      NAMELIST/namthd_zdf/ ln_zdf_BL99, ln_cndi_U64, ln_cndi_P07, rn_cnd_s, rn_kappa_i
       !!-------------------------------------------------------------------
       !
@@ -694 +967 @@
          WRITE(numout,*) '~~~~~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namthd_zdf:'
-         WRITE(numout,*) '      Diffusion follows a Bitz and Lipscomb (1999)            ln_zdf_BL99  = ', ln_zdf_BL99
+         WRITE(numout,*) '      Bitz and Lipscomb (1999) formulation                    ln_zdf_BL99  = ', ln_zdf_BL99
          WRITE(numout,*) '      thermal conductivity in the ice (Untersteiner 1964)     ln_cndi_U64  = ', ln_cndi_U64
          WRITE(numout,*) '      thermal conductivity in the ice (Pringle et al 2007)    ln_cndi_P07  = ', ln_cndi_P07
          WRITE(numout,*) '      thermal conductivity in the snow                        rn_cnd_s     = ', rn_cnd_s
          WRITE(numout,*) '      extinction radiation parameter in sea ice               rn_kappa_i   = ', rn_kappa_i
-         WRITE(numout,*) '      change the surface non-solar flux with Tsu or not       ln_dqns_i    = ', ln_dqns_i
      ENDIF
+
       !
       IF ( ( ln_cndi_U64 .AND. ln_cndi_P07 ) .OR. ( .NOT.ln_cndi_U64 .AND. .NOT.ln_cndi_P07 ) ) THEN
          CALL ctl_stop( 'ice_thd_zdf_init: choose one and only one formulation for thermal conduction (ln_cndi_U64 or ln_cndi_P07)' )
       ENDIF
+
+      !                             !== set the choice of ice vertical thermodynamic formulation ==!
+      ioptio = 0 
+      !      !--- BL99 thermo dynamics                               (linear liquidus + constant thermal properties)
+      IF( ln_zdf_BL99 ) THEN   ;   ioptio = ioptio + 1   ;   nice_zdf = np_BL99   ;   ENDIF
+      IF( ioptio /= 1 )    CALL ctl_stop( 'ice_thd_init: one and only one ice thermo option has to be defined ' )
       !
    END SUBROUTINE ice_thd_zdf_init

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90

@@ -44 +44 @@
    !!   ice_var_salprof1d : salinity profile in the ice 1D
    !!   ice_var_zapsmall  : remove very small area and volume
-   !!   ice_var_itd       : convert 1-cat to multiple cat
+   !!   ice_var_itd       : convert 1-cat to jpl-cat
+   !!   ice_var_itd2      : convert N-cat to jpl-cat
    !!   ice_var_bv        : brine volume
    !!----------------------------------------------------------------------
@@ -67 +68 @@
    PUBLIC   ice_var_zapsmall
    PUBLIC   ice_var_itd
+   PUBLIC   ice_var_itd2
    PUBLIC   ice_var_bv           
 
@@ -96 +98 @@
       et_s(:,:)  = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 )
       et_i(:,:)  = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 )
-
+      !
       at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds
       vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 )
-
-      ato_i(:,:) = 1._wp - at_i(:,:)    ! open water fraction  
+      !
+      ato_i(:,:) = 1._wp - at_i(:,:)         ! open water fraction  
 
       IF( kn > 1 ) THEN
@@ -238 +240 @@
          END WHERE
       END DO
-
+      !
       ! integrated values 
       vt_i (:,:) = SUM( v_i, dim=3 )
@@ -322 +324 @@
             END DO
          END DO
-
+         !
          ! Computation of the profile
          DO jl = 1, jpl
@@ -362 +364 @@
    END SUBROUTINE ice_var_salprof
 
+
    SUBROUTINE ice_var_salprof1d
       !!-------------------------------------------------------------------
@@ -457 +460 @@
          ELSEWHERE                                                                           ;   zswitch(:,:) = 1._wp
          END WHERE
-         
+         !
          DO jk = 1, nlay_i
             DO jj = 1 , jpj
@@ -468 +471 @@
             END DO
          END DO
-
+         !
          DO jj = 1 , jpj
             DO ji = 1 , jpi
@@ -484 +487 @@
                t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
                e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zswitch(ji,jj)
-
+               !
                !-----------------------------------------------------------------
                ! zap ice and snow volume, add water and salt to ocean
@@ -495 +498 @@
                oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj)
                sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj)
-
+               !
                h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj)
                h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj)
-
+               !
                a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
                v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
-
-            END DO
-         END DO
-        
+               !
+            END DO
+         END DO
+         !
       END DO 
 
@@ -548 +551 @@
       !!-------------------------------------------------------------------
       INTEGER  :: ji, jk, jl             ! dummy loop indices
-      INTEGER  :: ijpij, i_fill, jl0  
+      INTEGER  :: idim, i_fill, jl0  
       REAL(wp) :: zarg, zV, zconv, zdh, zdv
       REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
@@ -561 +564 @@
       ! then we check whether the distribution fullfills
       ! volume and area conservation, positivity and ice categories bounds
-      ijpij = SIZE( zhti , 1 )
-      zh_i(1:ijpij,1:jpl) = 0._wp
-      zh_s(1:ijpij,1:jpl) = 0._wp
-      za_i (1:ijpij,1:jpl) = 0._wp
-
-      DO ji = 1, ijpij
-         
+      idim = SIZE( zhti , 1 )
+      zh_i(1:idim,1:jpl) = 0._wp
+      zh_s(1:idim,1:jpl) = 0._wp
+      za_i(1:idim,1:jpl) = 0._wp
+      !
+      DO ji = 1, idim
+         !
          IF( zhti(ji) > 0._wp ) THEN
-
+            !
             ! find which category (jl0) the input ice thickness falls into
             jl0 = jpl
@@ -578 +581 @@
                ENDIF
             END DO
-
+            !
             itest(:) = 0
             i_fill   = jpl + 1                                            !------------------------------------
@@ -586 +589 @@
                !
                zh_i(ji,1:jpl) = 0._wp
-               za_i (ji,1:jpl) = 0._wp
-               itest(:)        = 0      
-               
+               za_i(ji,1:jpl) = 0._wp
+               itest(:)       = 0      
+               !
                IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1
                   zh_i(ji,1) = zhti(ji)
@@ -597 +600 @@
                      zh_i(ji,jl) = hi_mean(jl)
                   END DO
-                  
+                  !
                   ! concentration
                   za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl))
@@ -606 +609 @@
                      ENDIF
                   END DO
-                  
+                  !
                   ! last category
                   za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) )
                   zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) )
                   zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) 
-                  
+                  !
                   ! clem: correction if concentration of upper cat is greater than lower cat
                   !       (it should be a gaussian around jl0 but sometimes it is not)
@@ -622 +625 @@
                            za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
                         END IF
-                     ENDDO
+                     END DO
                   ENDIF
-               
+                  !
                ENDIF
-            
+               !
                ! Compatibility tests
                zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) ) 
-               IF ( zconv < epsi06 ) itest(1) = 1                                        ! Test 1: area conservation
-            
+               IF ( zconv < epsi06 )   itest(1) = 1                                        ! Test 1: area conservation
+               !
                zconv = ABS( zhti(ji)*zati(ji) - SUM( za_i(ji,1:jpl)*zh_i(ji,1:jpl) ) )
-               IF ( zconv < epsi06 ) itest(2) = 1                                        ! Test 2: volume conservation
-               
-               IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1                  ! Test 3: thickness of the last category is in-bounds ?
-               
+               IF ( zconv < epsi06 )   itest(2) = 1                                        ! Test 2: volume conservation
+               !
+               IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) )   itest(3) = 1                  ! Test 3: thickness of the last category is in-bounds ?
+               !
                itest(4) = 1
                DO jl = 1, i_fill
@@ -642 +645 @@
                !                                         !----------------------------
             END DO                                       ! end iteration on categories
-               !                                         !----------------------------
+            !                                            !----------------------------
          ENDIF
       END DO
@@ -648 +651 @@
       ! Add Snow in each category where za_i is not 0
       DO jl = 1, jpl
-         DO ji = 1, ijpij
+         DO ji = 1, idim
             IF( za_i(ji,jl) > 0._wp ) THEN
                zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) )
@@ -661 +664 @@
       END DO
       !
-    END SUBROUTINE ice_var_itd
-
-
-    SUBROUTINE ice_var_bv
+   END SUBROUTINE ice_var_itd
+
+
+   SUBROUTINE ice_var_itd2( zhti, zhts, zati, zh_i, zh_s, za_i )
+      !!-------------------------------------------------------------------
+      !!                ***  ROUTINE ice_var_itd2   ***
+      !!
+      !! ** Purpose :  converting N-cat ice to jpl ice categories
+      !!
+      !!                  ice thickness distribution follows a gaussian law
+      !!               around the concentration of the most likely ice thickness
+      !!                           (similar as iceistate.F90)
+      !!
+      !! ** Method:   Iterative procedure
+      !!                
+      !!               1) Fill ice cat that correspond to input thicknesses
+      !!                  Find the lowest(jlmin) and highest(jlmax) cat that are filled
+      !!
+      !!               2) Expand the filling to the cat jlmin-1 and jlmax+1
+      !!                   by removing 25% ice area from jlmin and jlmax (resp.) 
+      !!              
+      !!               3) Expand the filling to the empty cat between jlmin and jlmax 
+      !!                   by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)
+      !!                      b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
+      !!
+      !! ** Arguments : zhti: N-cat ice thickness
+      !!                zhts: N-cat snow depth
+      !!                zati: N-cat ice concentration
+      !!
+      !! ** Output    : jpl-cat 
+      !!
+      !!  (Example of application: BDY forcings when inputs have N-cat /= jpl)  
+      !!-------------------------------------------------------------------
+      INTEGER  ::   ji, jl, jl1, jl2             ! dummy loop indices
+      INTEGER  ::   idim, icat  
+      INTEGER, PARAMETER ::   ztrans = 0.25_wp
+      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables
+      INTEGER , DIMENSION(:,:), ALLOCATABLE   ::   jlfil, jlfil2
+      INTEGER , DIMENSION(:)  , ALLOCATABLE   ::   jlmax, jlmin
+      !!-------------------------------------------------------------------
+      !
+      idim = SIZE( zhti, 1 )
+      icat = SIZE( zhti, 2 )
+      !
+      ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) )       ! allocate arrays
+      ALLOCATE( jlmin(idim), jlmax(idim) )
+
+      ! --- initialize output fields to 0 --- !
+      zh_i(1:idim,1:jpl) = 0._wp
+      zh_s(1:idim,1:jpl) = 0._wp
+      za_i(1:idim,1:jpl) = 0._wp
+      !
+      ! --- fill the categories --- !
+      !     find where cat-input = cat-output and fill cat-output fields  
+      jlmax(:) = 0
+      jlmin(:) = 999
+      jlfil(:,:) = 0
+      DO jl1 = 1, jpl
+         DO jl2 = 1, icat
+            DO ji = 1, idim
+               IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN
+                  ! fill the right category
+                  zh_i(ji,jl1) = zhti(ji,jl2)
+                  zh_s(ji,jl1) = zhts(ji,jl2)
+                  za_i(ji,jl1) = zati(ji,jl2)
+                  ! record categories that are filled
+                  jlmax(ji) = MAX( jlmax(ji), jl1 )
+                  jlmin(ji) = MIN( jlmin(ji), jl1 )
+                  jlfil(ji,jl1) = jl1
+               ENDIF
+            END DO
+         END DO
+      END DO
+      !
+      ! --- fill the gaps between categories --- !  
+      !     transfer from categories filled at the previous step to the empty ones in between
+      DO ji = 1, idim
+         jl1 = jlmin(ji)
+         jl2 = jlmax(ji)
+         IF( jl1 > 1 ) THEN
+            ! fill the lower cat (jl1-1)
+            za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)
+            zh_i(ji,jl1-1) = hi_mean(jl1-1)
+            ! remove from cat jl1
+            za_i(ji,jl1  ) = ( 1._wp - ztrans ) * za_i(ji,jl1)
+         ENDIF
+         IF( jl2 < jpl ) THEN
+            ! fill the upper cat (jl2+1)
+            za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)
+            zh_i(ji,jl2+1) = hi_mean(jl2+1)
+            ! remove from cat jl2
+            za_i(ji,jl2  ) = ( 1._wp - ztrans ) * za_i(ji,jl2)
+         ENDIF
+      END DO
+      !
+      jlfil2(:,:) = jlfil(:,:) 
+      ! fill categories from low to high
+      DO jl = 2, jpl-1
+         DO ji = 1, idim
+            IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
+               ! fill high
+               za_i(ji,jl) = ztrans * za_i(ji,jl-1)
+               zh_i(ji,jl) = hi_mean(jl)
+               jlfil(ji,jl) = jl
+               ! remove low
+               za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)
+            ENDIF
+         END DO
+      END DO
+      !
+      ! fill categories from high to low
+      DO jl = jpl-1, 2, -1
+         DO ji = 1, idim
+            IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
+               ! fill low
+               za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)
+               zh_i(ji,jl) = hi_mean(jl) 
+               jlfil2(ji,jl) = jl
+               ! remove high
+               za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)
+            ENDIF
+         END DO
+      END DO
+      !
+      DEALLOCATE( jlfil, jlfil2 )      ! deallocate arrays
+      DEALLOCATE( jlmin, jlmax )
+      !
+   END SUBROUTINE ice_var_itd2
+
+
+   SUBROUTINE ice_var_bv
       !!-------------------------------------------------------------------
       !!                ***  ROUTINE ice_var_bv ***

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90

@@ -51 +51 @@
 
 #if defined key_lim3
-   LOGICAL :: ll_bdylim3                  ! determine whether ice input is 1cat (F) or Xcat (T) type
-   INTEGER :: jfld_hti, jfld_hts, jfld_ai ! indices of ice thickness, snow thickness and concentration in bf structure
+   INTEGER ::   nice_cat                      ! number of categories in the input file
+   INTEGER ::   jfld_hti, jfld_hts, jfld_ai   ! indices of ice thickness, snow thickness and concentration in bf structure
 #endif
 
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
    !! $Id$ 
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
@@ -80 +80 @@
       !                                               ! etc.
       !
-      INTEGER ::  ib_bdy, jfld, jstart, jend, ib, ii, ij, ik, igrd, jl  ! local indices
+      INTEGER ::  jbdy, jfld, jstart, jend, ib, jl  ! dummy loop indices
+      INTEGER ::  ii, ij, ik, igrd                  ! local integers
       INTEGER,          DIMENSION(jpbgrd) ::   ilen1 
       INTEGER, POINTER, DIMENSION(:)      ::   nblen, nblenrim  ! short cuts
@@ -95 +96 @@
          !-----------------------------
          
-         DO ib_bdy = 1, nb_bdy
+         DO jbdy = 1, nb_bdy
             !
-            nblen => idx_bdy(ib_bdy)%nblen
-            nblenrim => idx_bdy(ib_bdy)%nblenrim
-            dta => dta_bdy(ib_bdy)
-
-            IF( nn_dyn2d_dta(ib_bdy) == 0 ) THEN 
+            nblen    => idx_bdy(jbdy)%nblen
+            nblenrim => idx_bdy(jbdy)%nblenrim
+            dta      => dta_bdy(jbdy)
+            !
+            IF( nn_dyn2d_dta(jbdy) == 0 ) THEN 
                ilen1(:) = nblen(:)
                IF( dta%ll_ssh ) THEN 
                   igrd = 1
                   DO ib = 1, ilen1(igrd)
-                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                     dta_bdy(ib_bdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1)         
+                     ii = idx_bdy(jbdy)%nbi(ib,igrd)
+                     ij = idx_bdy(jbdy)%nbj(ib,igrd)
+                     dta_bdy(jbdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1)         
                   END DO 
-               END IF
+               ENDIF
                IF( dta%ll_u2d ) THEN 
                   igrd = 2
                   DO ib = 1, ilen1(igrd)
-                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                     dta_bdy(ib_bdy)%u2d(ib) = un_b(ii,ij) * umask(ii,ij,1)         
+                     ii = idx_bdy(jbdy)%nbi(ib,igrd)
+                     ij = idx_bdy(jbdy)%nbj(ib,igrd)
+                     dta_bdy(jbdy)%u2d(ib) = un_b(ii,ij) * umask(ii,ij,1)         
                   END DO 
-               END IF
+               ENDIF
                IF( dta%ll_v2d ) THEN 
                   igrd = 3
                   DO ib = 1, ilen1(igrd)
-                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                     dta_bdy(ib_bdy)%v2d(ib) = vn_b(ii,ij) * vmask(ii,ij,1)         
+                     ii = idx_bdy(jbdy)%nbi(ib,igrd)
+                     ij = idx_bdy(jbdy)%nbj(ib,igrd)
+                     dta_bdy(jbdy)%v2d(ib) = vn_b(ii,ij) * vmask(ii,ij,1)         
                   END DO 
-               END IF
-            ENDIF
-
-            IF( nn_dyn3d_dta(ib_bdy) == 0 ) THEN 
+               ENDIF
+            ENDIF
+            !
+            IF( nn_dyn3d_dta(jbdy) == 0 ) THEN 
                ilen1(:) = nblen(:)
                IF( dta%ll_u3d ) THEN 
@@ -135 +136 @@
                   DO ib = 1, ilen1(igrd)
                      DO ik = 1, jpkm1
-                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                        dta_bdy(ib_bdy)%u3d(ib,ik) =  ( un(ii,ij,ik) - un_b(ii,ij) ) * umask(ii,ij,ik)         
+                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
+                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
+                        dta_bdy(jbdy)%u3d(ib,ik) =  ( un(ii,ij,ik) - un_b(ii,ij) ) * umask(ii,ij,ik)         
                      END DO
                   END DO 
-               END IF
+               ENDIF
                IF( dta%ll_v3d ) THEN 
                   igrd = 3 
                   DO ib = 1, ilen1(igrd)
                      DO ik = 1, jpkm1
-                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                        dta_bdy(ib_bdy)%v3d(ib,ik) =  ( vn(ii,ij,ik) - vn_b(ii,ij) ) * vmask(ii,ij,ik)         
+                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
+                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
+                        dta_bdy(jbdy)%v3d(ib,ik) =  ( vn(ii,ij,ik) - vn_b(ii,ij) ) * vmask(ii,ij,ik)         
                         END DO
                   END DO 
-               END IF
-            ENDIF
-
-            IF( nn_tra_dta(ib_bdy) == 0 ) THEN 
+               ENDIF
+            ENDIF
+
+            IF( nn_tra_dta(jbdy) == 0 ) THEN 
                ilen1(:) = nblen(:)
                IF( dta%ll_tem ) THEN
@@ -159 +160 @@
                   DO ib = 1, ilen1(igrd)
                      DO ik = 1, jpkm1
-                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                        dta_bdy(ib_bdy)%tem(ib,ik) = tsn(ii,ij,ik,jp_tem) * tmask(ii,ij,ik)         
+                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
+                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
+                        dta_bdy(jbdy)%tem(ib,ik) = tsn(ii,ij,ik,jp_tem) * tmask(ii,ij,ik)         
                      END DO
                   END DO 
-               END IF
+               ENDIF
                IF( dta%ll_sal ) THEN
                   igrd = 1 
                   DO ib = 1, ilen1(igrd)
                      DO ik = 1, jpkm1
-                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                        dta_bdy(ib_bdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_sal) * tmask(ii,ij,ik)         
+                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
+                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
+                        dta_bdy(jbdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_sal) * tmask(ii,ij,ik)         
                      END DO
                   END DO 
-               END IF
-            ENDIF
-
-#if defined key_lim3
-            IF( nn_ice_lim_dta(ib_bdy) == 0 ) THEN 
+               ENDIF
+            ENDIF
+
+#if defined key_lim3
+            IF( nn_ice_lim_dta(jbdy) == 0 ) THEN    ! set ice to initial values
                ilen1(:) = nblen(:)
                IF( dta%ll_a_i ) THEN
@@ -184 +185 @@
                   DO jl = 1, jpl
                      DO ib = 1, ilen1(igrd)
-                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                        dta_bdy(ib_bdy)%a_i (ib,jl) =  a_i(ii,ij,jl) * tmask(ii,ij,1) 
+                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
+                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
+                        dta_bdy(jbdy)%a_i (ib,jl) =  a_i(ii,ij,jl) * tmask(ii,ij,1) 
                      END DO
                   END DO
@@ -194 +195 @@
                   DO jl = 1, jpl
                      DO ib = 1, ilen1(igrd)
-                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                        dta_bdy(ib_bdy)%h_i (ib,jl) =  h_i(ii,ij,jl) * tmask(ii,ij,1) 
+                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
+                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
+                        dta_bdy(jbdy)%h_i (ib,jl) =  h_i(ii,ij,jl) * tmask(ii,ij,1) 
                      END DO
                   END DO
@@ -204 +205 @@
                   DO jl = 1, jpl
                      DO ib = 1, ilen1(igrd)
-                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                        dta_bdy(ib_bdy)%h_s (ib,jl) =  h_s(ii,ij,jl) * tmask(ii,ij,1) 
+                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
+                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
+                        dta_bdy(jbdy)%h_s (ib,jl) =  h_s(ii,ij,jl) * tmask(ii,ij,1) 
                      END DO
                   END DO
@@ -212 +213 @@
             ENDIF
 #endif
-         END DO ! ib_bdy
+         END DO ! jbdy
          !
       ENDIF ! kt == nit000
@@ -220 +221 @@
      
       jstart = 1
-      DO ib_bdy = 1, nb_bdy   
-         dta => dta_bdy(ib_bdy)
-         IF( nn_dta(ib_bdy) == 1 ) THEN ! skip this bit if no external data required
+      DO jbdy = 1, nb_bdy   
+         dta => dta_bdy(jbdy)
+         IF( nn_dta(jbdy) == 1 ) THEN ! skip this bit if no external data required
       
             IF( PRESENT(jit) ) THEN
                ! Update barotropic boundary conditions only
                ! jit is optional argument for fld_read and bdytide_update
-               IF( cn_dyn2d(ib_bdy) /= 'none' ) THEN
-                  IF( nn_dyn2d_dta(ib_bdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
+               IF( cn_dyn2d(jbdy) /= 'none' ) THEN
+                  IF( nn_dyn2d_dta(jbdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
                      IF( dta%ll_ssh ) dta%ssh(:) = 0._wp
                      IF( dta%ll_u2d ) dta%u2d(:) = 0._wp
                      IF( dta%ll_u3d ) dta%v2d(:) = 0._wp
                   ENDIF
-                  IF (cn_tra(ib_bdy) /= 'runoff') THEN
-                     IF( nn_dyn2d_dta(ib_bdy) == 1 .OR. nn_dyn2d_dta(ib_bdy) == 3 ) THEN
+                  IF (cn_tra(jbdy) /= 'runoff') THEN
+                     IF( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 ) THEN
 
                         jend = jstart + dta%nread(2) - 1
-                        IF( ln_full_vel_array(ib_bdy) ) THEN
+                        IF( ln_full_vel_array(jbdy) ) THEN
                            CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  &
-                                     & kit=jit, kt_offset=time_offset , jpk_bdy=nb_jpk_bdy, fvl=ln_full_vel_array(ib_bdy)  )
+                                     & kit=jit, kt_offset=time_offset , jpk_bdy=nb_jpk_bdy, fvl=ln_full_vel_array(jbdy)  )
                         ELSE
                            CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  &
@@ -246 +247 @@
 
                         ! If full velocities in boundary data then extract barotropic velocities from 3D fields
-                        IF( ln_full_vel_array(ib_bdy) .AND.                                             &
-                          &    ( nn_dyn2d_dta(ib_bdy) == 1 .OR. nn_dyn2d_dta(ib_bdy) == 3 .OR.  &
-                          &      nn_dyn3d_dta(ib_bdy) == 1 ) )THEN
+                        IF( ln_full_vel_array(jbdy) .AND.                                             &
+                          &    ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 .OR.  &
+                          &      nn_dyn3d_dta(jbdy) == 1 ) )THEN
 
                            igrd = 2                      ! zonal velocity
                            dta%u2d(:) = 0._wp
-                           DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
-                              ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                              ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                           DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
+                              ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+                              ij   = idx_bdy(jbdy)%nbj(ib,igrd)
                               DO ik = 1, jpkm1
                                  dta%u2d(ib) = dta%u2d(ib) &
@@ -263 +264 @@
                            igrd = 3                      ! meridional velocity
                            dta%v2d(:) = 0._wp
-                           DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
-                              ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                              ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                           DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
+                              ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+                              ij   = idx_bdy(jbdy)%nbj(ib,igrd)
                               DO ik = 1, jpkm1
                                  dta%v2d(ib) = dta%v2d(ib) &
@@ -274 +275 @@
                         ENDIF                    
                      ENDIF
-                     IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing
-                        CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta, td=tides(ib_bdy),   & 
+                     IF( nn_dyn2d_dta(jbdy) .ge. 2 ) THEN ! update tidal harmonic forcing
+                        CALL bdytide_update( kt=kt, idx=idx_bdy(jbdy), dta=dta, td=tides(jbdy),   & 
                           &                 jit=jit, time_offset=time_offset )
                      ENDIF
@@ -281 +282 @@
                ENDIF
             ELSE
-               IF (cn_tra(ib_bdy) == 'runoff') then      ! runoff condition
-                  jend = nb_bdy_fld(ib_bdy)
+               IF (cn_tra(jbdy) == 'runoff') then      ! runoff condition
+                  jend = nb_bdy_fld(jbdy)
                   CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend),  &
                                & map=nbmap_ptr(jstart:jend), kt_offset=time_offset )
                   !
                   igrd = 2                      ! zonal velocity
-                  DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
-                     ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                     ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                  DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
+                     ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+                     ij   = idx_bdy(jbdy)%nbj(ib,igrd)
                      dta%u2d(ib) = dta%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) )
                   END DO
                   !
                   igrd = 3                      ! meridional velocity
-                  DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
-                     ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                     ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                  DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
+                     ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+                     ij   = idx_bdy(jbdy)%nbj(ib,igrd)
                      dta%v2d(ib) = dta%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) )
                   END DO
                ELSE
-                  IF( nn_dyn2d_dta(ib_bdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
+                  IF( nn_dyn2d_dta(jbdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
                      IF( dta%ll_ssh ) dta%ssh(:) = 0._wp
                      IF( dta%ll_u2d ) dta%u2d(:) = 0._wp
@@ -308 +309 @@
                      jend = jstart + dta%nread(1) - 1
                      CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), &
-                                  & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy, fvl=ln_full_vel_array(ib_bdy) )
+                                  & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy, fvl=ln_full_vel_array(jbdy) )
                   ENDIF
                   ! If full velocities in boundary data then split into barotropic and baroclinic data
-                  IF( ln_full_vel_array(ib_bdy) .and.                                             &
-                    & ( nn_dyn2d_dta(ib_bdy) == 1 .OR. nn_dyn2d_dta(ib_bdy) == 3 .OR. &
-                    &   nn_dyn3d_dta(ib_bdy) == 1 ) ) THEN
+                  IF( ln_full_vel_array(jbdy) .and.                                             &
+                    & ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 .OR. &
+                    &   nn_dyn3d_dta(jbdy) == 1 ) ) THEN
                      igrd = 2                      ! zonal velocity
                      dta%u2d(:) = 0._wp
-                     DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
-                        ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                        ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                     DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
+                        ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+                        ij   = idx_bdy(jbdy)%nbj(ib,igrd)
                         DO ik = 1, jpkm1
                            dta%u2d(ib) = dta%u2d(ib) &
@@ -330 +331 @@
                      igrd = 3                      ! meridional velocity
                      dta%v2d(:) = 0._wp
-                     DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
-                        ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                        ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                     DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
+                        ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+                        ij   = idx_bdy(jbdy)%nbj(ib,igrd)
                         DO ik = 1, jpkm1
                            dta%v2d(ib) = dta%v2d(ib) &
@@ -346 +347 @@
                ENDIF
 #if defined key_lim3
-               IF( .NOT. ll_bdylim3 .AND. cn_ice_lim(ib_bdy) /= 'none' .AND. nn_ice_lim_dta(ib_bdy) == 1 ) THEN ! bdy ice input (case input is 1cat)
-                CALL ice_var_itd ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), &
-                                  & dta_bdy(ib_bdy)%h_i,     dta_bdy(ib_bdy)%h_s,     dta_bdy(ib_bdy)%a_i     )
+               IF( cn_ice_lim(jbdy) /= 'none' .AND. nn_ice_lim_dta(jbdy) == 1 ) THEN
+                  IF( nice_cat == 1 ) THEN ! case input cat = 1
+                     CALL ice_var_itd ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), &
+                        &               dta_bdy(jbdy)%h_i     , dta_bdy(jbdy)%h_s     , dta_bdy(jbdy)%a_i    )
+                  ELSEIF( nice_cat /= 1 .AND. nice_cat /= jpl ) THEN ! case input cat /=1 and /=jpl
+                     CALL ice_var_itd2( bf(jfld_hti)%fnow(:,1,:), bf(jfld_hts)%fnow(:,1,:), bf(jfld_ai)%fnow(:,1,:), &
+                        &               dta_bdy(jbdy)%h_i     , dta_bdy(jbdy)%h_s     , dta_bdy(jbdy)%a_i    )
+                  ENDIF
                ENDIF
 #endif
             ENDIF
             jstart = jstart + dta%nread(1)
-         END IF ! nn_dta(ib_bdy) = 1
-      END DO  ! ib_bdy
+         ENDIF    ! nn_dta(jbdy) = 1
+      END DO  ! jbdy
 
       IF ( ln_tide ) THEN
          IF (ln_dynspg_ts) THEN      ! Fill temporary arrays with slow-varying bdy data                           
-            DO ib_bdy = 1, nb_bdy    ! Tidal component added in ts loop
-               IF ( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN
-                  nblen => idx_bdy(ib_bdy)%nblen
-                  nblenrim => idx_bdy(ib_bdy)%nblenrim
-                  IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN; ilen1(:)=nblen(:) ; ELSE ; ilen1(:)=nblenrim(:) ; ENDIF 
-                  IF ( dta_bdy(ib_bdy)%ll_ssh ) dta_bdy_s(ib_bdy)%ssh(1:ilen1(1)) = dta_bdy(ib_bdy)%ssh(1:ilen1(1))
-                  IF ( dta_bdy(ib_bdy)%ll_u2d ) dta_bdy_s(ib_bdy)%u2d(1:ilen1(2)) = dta_bdy(ib_bdy)%u2d(1:ilen1(2))
-                  IF ( dta_bdy(ib_bdy)%ll_v2d ) dta_bdy_s(ib_bdy)%v2d(1:ilen1(3)) = dta_bdy(ib_bdy)%v2d(1:ilen1(3))
+            DO jbdy = 1, nb_bdy    ! Tidal component added in ts loop
+               IF ( nn_dyn2d_dta(jbdy) .ge. 2 ) THEN
+                  nblen => idx_bdy(jbdy)%nblen
+                  nblenrim => idx_bdy(jbdy)%nblenrim
+                  IF( cn_dyn2d(jbdy) == 'frs' ) THEN; ilen1(:)=nblen(:) ; ELSE ; ilen1(:)=nblenrim(:) ; ENDIF 
+                  IF ( dta_bdy(jbdy)%ll_ssh ) dta_bdy_s(jbdy)%ssh(1:ilen1(1)) = dta_bdy(jbdy)%ssh(1:ilen1(1))
+                  IF ( dta_bdy(jbdy)%ll_u2d ) dta_bdy_s(jbdy)%u2d(1:ilen1(2)) = dta_bdy(jbdy)%u2d(1:ilen1(2))
+                  IF ( dta_bdy(jbdy)%ll_v2d ) dta_bdy_s(jbdy)%v2d(1:ilen1(3)) = dta_bdy(jbdy)%v2d(1:ilen1(3))
                ENDIF
             END DO
@@ -375 +381 @@
 
       IF ( ln_apr_obc ) THEN
-         DO ib_bdy = 1, nb_bdy
-            IF (cn_tra(ib_bdy) /= 'runoff')THEN
+         DO jbdy = 1, nb_bdy
+            IF (cn_tra(jbdy) /= 'runoff')THEN
                igrd = 1                      ! meridional velocity
-               DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
-                  ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                  ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                  dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + ssh_ib(ii,ij)
+               DO ib = 1, idx_bdy(jbdy)%nblenrim(igrd)
+                  ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+                  ij   = idx_bdy(jbdy)%nbj(ib,igrd)
+                  dta_bdy(jbdy)%ssh(ib) = dta_bdy(jbdy)%ssh(ib) + ssh_ib(ii,ij)
                END DO
             ENDIF
@@ -402 +408 @@
       !!                
       !!----------------------------------------------------------------------
-      INTEGER ::   ib_bdy, jfld, jstart, jend, ierror, ios     ! Local integers
+      INTEGER ::   jbdy, jfld, jstart, jend, ierror, ios     ! Local integers
       !
       CHARACTER(len=100)                     ::   cn_dir        ! Root directory for location of data files
@@ -408 +414 @@
       CHARACTER(len = 256)::   clname                           ! temporary file name
       LOGICAL                                ::   ln_full_vel   ! =T => full velocities in 3D boundary data
-                                                                ! =F => baroclinic velocities in 3D boundary data
+      !                                                         ! =F => baroclinic velocities in 3D boundary data
       INTEGER                                ::   ilen_global   ! Max length required for global bdy dta arrays
       INTEGER, ALLOCATABLE, DIMENSION(:)     ::   ilen1, ilen3  ! size of 1st and 3rd dimensions of local arrays
@@ -416 +422 @@
       TYPE(OBC_DATA), POINTER                ::   dta           ! short cut
 #if defined key_lim3
-      INTEGER               ::   zndims   ! number of dimensions in an array (i.e. 3 = wo ice cat; 4 = w ice cat)
+      INTEGER               ::   kndims   ! number of dimensions in an array (i.e. 3 = wo ice cat; 4 = w ice cat)
+      INTEGER, DIMENSION(4) ::   kdimsz   ! size   of dimensions
       INTEGER               ::   inum,id1 ! local integer
 #endif
@@ -440 +447 @@
 
       ! Set nn_dta
-      DO ib_bdy = 1, nb_bdy
-         nn_dta(ib_bdy) = MAX(  nn_dyn2d_dta(ib_bdy)       &
-                               ,nn_dyn3d_dta(ib_bdy)       &
-                               ,nn_tra_dta(ib_bdy)         &
-#if defined key_lim3
-                              ,nn_ice_lim_dta(ib_bdy)    &
+      DO jbdy = 1, nb_bdy
+         nn_dta(jbdy) = MAX(   nn_dyn2d_dta  (jbdy)    &
+            &                , nn_dyn3d_dta  (jbdy)    &
+            &                , nn_tra_dta    (jbdy)    &
+#if defined key_lim3
+            &                , nn_ice_lim_dta(jbdy)    &
 #endif
                               )
-         IF(nn_dta(ib_bdy) > 1) nn_dta(ib_bdy) = 1
+         IF(nn_dta(jbdy) > 1)   nn_dta(jbdy) = 1
       END DO
 
@@ -455 +462 @@
       ALLOCATE( nb_bdy_fld(nb_bdy) )
       nb_bdy_fld(:) = 0
-      DO ib_bdy = 1, nb_bdy         
-         IF( cn_dyn2d(ib_bdy) /= 'none' .and. ( nn_dyn2d_dta(ib_bdy) == 1 .or. nn_dyn2d_dta(ib_bdy) == 3 ) ) THEN
-            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 3
-         ENDIF
-         IF( cn_dyn3d(ib_bdy) /= 'none' .and. nn_dyn3d_dta(ib_bdy) == 1 ) THEN
-            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 2
-         ENDIF
-         IF( cn_tra(ib_bdy) /= 'none' .and. nn_tra_dta(ib_bdy) == 1  ) THEN
-            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 2
-         ENDIF
-#if defined key_lim3
-         IF( cn_ice_lim(ib_bdy) /= 'none' .and. nn_ice_lim_dta(ib_bdy) == 1  ) THEN
-            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 3
+      DO jbdy = 1, nb_bdy         
+         IF( cn_dyn2d(jbdy) /= 'none' .AND. ( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) ) THEN
+            nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 3
+         ENDIF
+         IF( cn_dyn3d(jbdy) /= 'none' .AND. nn_dyn3d_dta(jbdy) == 1 ) THEN
+            nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 2
+         ENDIF
+         IF( cn_tra(jbdy) /= 'none' .AND. nn_tra_dta(jbdy) == 1  ) THEN
+            nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 2
+         ENDIF
+#if defined key_lim3
+         IF( cn_ice_lim(jbdy) /= 'none' .AND. nn_ice_lim_dta(jbdy) == 1  ) THEN
+            nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 3
          ENDIF
 #endif               
-         IF(lwp) WRITE(numout,*) 'Maximum number of files to open =',nb_bdy_fld(ib_bdy)
+         IF(lwp) WRITE(numout,*) 'Maximum number of files to open =', nb_bdy_fld(jbdy)
       END DO            
 
@@ -496 +503 @@
       REWIND(numnam_cfg)
       jfld = 0 
-      DO ib_bdy = 1, nb_bdy         
-         IF( nn_dta(ib_bdy) == 1 ) THEN
+      DO jbdy = 1, nb_bdy         
+         IF( nn_dta(jbdy) == 1 ) THEN
             READ  ( numnam_ref, nambdy_dta, IOSTAT = ios, ERR = 901)
 901         IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_dta in reference namelist', lwp )
@@ -505 +512 @@
             IF(lwm) WRITE( numond, nambdy_dta )
 
-            cn_dir_array(ib_bdy) = cn_dir
-            ln_full_vel_array(ib_bdy) = ln_full_vel
-
-            nblen => idx_bdy(ib_bdy)%nblen
-            nblenrim => idx_bdy(ib_bdy)%nblenrim
-            dta => dta_bdy(ib_bdy)
+            cn_dir_array(jbdy) = cn_dir
+            ln_full_vel_array(jbdy) = ln_full_vel
+
+            nblen => idx_bdy(jbdy)%nblen
+            nblenrim => idx_bdy(jbdy)%nblenrim
+            dta => dta_bdy(jbdy)
             dta%nread(2) = 0
 
             ! Only read in necessary fields for this set.
             ! Important that barotropic variables come first.
-            IF( nn_dyn2d_dta(ib_bdy) == 1 .or. nn_dyn2d_dta(ib_bdy) == 3 ) THEN 
+            IF( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) THEN 
 
                IF( dta%ll_ssh ) THEN 
@@ -521 +528 @@
                   jfld = jfld + 1
                   blf_i(jfld) = bn_ssh
-                  ibdy(jfld) = ib_bdy
+                  ibdy(jfld) = jbdy
                   igrid(jfld) = 1
                   ilen1(jfld) = nblen(igrid(jfld))
@@ -528 +535 @@
                ENDIF
 
-               IF( dta%ll_u2d .and. .not. ln_full_vel_array(ib_bdy) ) THEN
+               IF( dta%ll_u2d .and. .not. ln_full_vel_array(jbdy) ) THEN
                   if(lwp) write(numout,*) '++++++ reading in u2d field'
                   jfld = jfld + 1
                   blf_i(jfld) = bn_u2d
-                  ibdy(jfld) = ib_bdy
+                  ibdy(jfld) = jbdy
                   igrid(jfld) = 2
                   ilen1(jfld) = nblen(igrid(jfld))
@@ -539 +546 @@
                ENDIF
 
-               IF( dta%ll_v2d .and. .not. ln_full_vel_array(ib_bdy) ) THEN
+               IF( dta%ll_v2d .and. .not. ln_full_vel_array(jbdy) ) THEN
                   if(lwp) write(numout,*) '++++++ reading in v2d field'
                   jfld = jfld + 1
                   blf_i(jfld) = bn_v2d
-                  ibdy(jfld) = ib_bdy
+                  ibdy(jfld) = jbdy
                   igrid(jfld) = 3
                   ilen1(jfld) = nblen(igrid(jfld))
@@ -554 +561 @@
             ! read 3D velocities if baroclinic velocities require OR if
             ! barotropic velocities required and ln_full_vel set to .true.
-            IF( nn_dyn3d_dta(ib_bdy) == 1 .OR. &
-           &  ( ln_full_vel_array(ib_bdy) .AND. ( nn_dyn2d_dta(ib_bdy) == 1 .or. nn_dyn2d_dta(ib_bdy) == 3 ) ) ) THEN
-
-               IF( dta%ll_u3d .OR. ( ln_full_vel_array(ib_bdy) .and. dta%ll_u2d ) ) THEN
+            IF( nn_dyn3d_dta(jbdy) == 1 .OR. &
+           &  ( ln_full_vel_array(jbdy) .AND. ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 ) ) ) THEN
+
+               IF( dta%ll_u3d .OR. ( ln_full_vel_array(jbdy) .and. dta%ll_u2d ) ) THEN
                   if(lwp) write(numout,*) '++++++ reading in u3d field'
                   jfld = jfld + 1
                   blf_i(jfld) = bn_u3d
-                  ibdy(jfld) = ib_bdy
+                  ibdy(jfld) = jbdy
                   igrid(jfld) = 2
                   ilen1(jfld) = nblen(igrid(jfld))
                   ilen3(jfld) = jpk
-                  IF( ln_full_vel_array(ib_bdy) .and. dta%ll_u2d ) dta%nread(2) = dta%nread(2) + 1
-               ENDIF
-
-               IF( dta%ll_v3d .OR. ( ln_full_vel_array(ib_bdy) .and. dta%ll_v2d ) ) THEN
+                  IF( ln_full_vel_array(jbdy) .and. dta%ll_u2d ) dta%nread(2) = dta%nread(2) + 1
+               ENDIF
+
+               IF( dta%ll_v3d .OR. ( ln_full_vel_array(jbdy) .and. dta%ll_v2d ) ) THEN
                   if(lwp) write(numout,*) '++++++ reading in v3d field'
                   jfld = jfld + 1
                   blf_i(jfld) = bn_v3d
-                  ibdy(jfld) = ib_bdy
+                  ibdy(jfld) = jbdy
                   igrid(jfld) = 3
                   ilen1(jfld) = nblen(igrid(jfld))
                   ilen3(jfld) = jpk
-                  IF( ln_full_vel_array(ib_bdy) .and. dta%ll_v2d ) dta%nread(2) = dta%nread(2) + 1
+                  IF( ln_full_vel_array(jbdy) .and. dta%ll_v2d ) dta%nread(2) = dta%nread(2) + 1
                ENDIF
 
@@ -582 +589 @@
 
             ! temperature and salinity
-            IF( nn_tra_dta(ib_bdy) == 1 ) THEN
+            IF( nn_tra_dta(jbdy) == 1 ) THEN
 
                IF( dta%ll_tem ) THEN
@@ -588 +595 @@
                   jfld = jfld + 1
                   blf_i(jfld) = bn_tem
-                  ibdy(jfld) = ib_bdy
+                  ibdy(jfld) = jbdy
                   igrid(jfld) = 1
                   ilen1(jfld) = nblen(igrid(jfld))
@@ -598 +605 @@
                   jfld = jfld + 1
                   blf_i(jfld) = bn_sal
-                  ibdy(jfld) = ib_bdy
+                  ibdy(jfld) = jbdy
                   igrid(jfld) = 1
                   ilen1(jfld) = nblen(igrid(jfld))
@@ -608 +615 @@
 #if defined key_lim3
             ! sea ice
-            IF( nn_ice_lim_dta(ib_bdy) == 1 ) THEN
+            IF( nn_ice_lim_dta(jbdy) == 1 ) THEN
                ! Test for types of ice input (1cat or Xcat) 
                ! Build file name to find dimensions 
@@ -622 +629 @@
                !
                CALL iom_open  ( clname, inum )
-               id1 = iom_varid( inum, bn_a_i%clvar, kndims=zndims, ldstop = .FALSE. )
+               id1 = iom_varid( inum, bn_a_i%clvar, kdimsz=kdimsz, kndims=kndims, ldstop = .FALSE. )
                CALL iom_close ( inum )
 
-                IF ( zndims == 4 ) THEN
-                 ll_bdylim3 = .TRUE.   ! Xcat input
+                IF ( kndims == 4 ) THEN
+                 nice_cat = kdimsz(4)   ! Xcat input
                ELSE
-                 ll_bdylim3 = .FALSE.  ! 1cat input      
+                 nice_cat = 1           ! 1cat input      
                ENDIF
                ! End test
@@ -635 +642 @@
                   jfld = jfld + 1
                   blf_i(jfld) = bn_a_i
-                  ibdy(jfld) = ib_bdy
+                  ibdy(jfld)  = jbdy
                   igrid(jfld) = 1
                   ilen1(jfld) = nblen(igrid(jfld))
-                  IF ( ll_bdylim3 ) THEN ; ilen3(jfld)=jpl ; ELSE ; ilen3(jfld)=1 ; ENDIF
+                  ilen3(jfld) = nice_cat
                ENDIF
 
@@ -644 +651 @@
                   jfld = jfld + 1
                   blf_i(jfld) = bn_h_i
-                  ibdy(jfld) = ib_bdy
+                  ibdy(jfld)  = jbdy
                   igrid(jfld) = 1
                   ilen1(jfld) = nblen(igrid(jfld))
-                  IF ( ll_bdylim3 ) THEN ; ilen3(jfld)=jpl ; ELSE ; ilen3(jfld)=1 ; ENDIF
+                  ilen3(jfld) = nice_cat
                ENDIF
 
                IF( dta%ll_h_s ) THEN
                   jfld = jfld + 1
-                   blf_i(jfld) = bn_h_s
-                  ibdy(jfld) = ib_bdy
+                  blf_i(jfld) = bn_h_s
+                  ibdy(jfld)  = jbdy
                   igrid(jfld) = 1
                   ilen1(jfld) = nblen(igrid(jfld))
-                  IF ( ll_bdylim3 ) THEN ; ilen3(jfld)=jpl ; ELSE ; ilen3(jfld)=1 ; ENDIF
+                  ilen3(jfld) = nice_cat
                ENDIF
 
@@ -663 +670 @@
             ! Recalculate field counts
             !-------------------------
-            IF( ib_bdy == 1 ) THEN 
+            IF( jbdy == 1 ) THEN 
                nb_bdy_fld_sum = 0
-               nb_bdy_fld(ib_bdy) = jfld
+               nb_bdy_fld(jbdy) = jfld
                nb_bdy_fld_sum     = jfld              
             ELSE
-               nb_bdy_fld(ib_bdy) = jfld - nb_bdy_fld_sum
-               nb_bdy_fld_sum = nb_bdy_fld_sum + nb_bdy_fld(ib_bdy)
-            ENDIF
-
-            dta%nread(1) = nb_bdy_fld(ib_bdy)
+               nb_bdy_fld(jbdy) = jfld - nb_bdy_fld_sum
+               nb_bdy_fld_sum = nb_bdy_fld_sum + nb_bdy_fld(jbdy)
+            ENDIF
+
+            dta%nread(1) = nb_bdy_fld(jbdy)
 
          ENDIF ! nn_dta == 1
-      ENDDO ! ib_bdy
+      ENDDO ! jbdy
 
       DO jfld = 1, nb_bdy_fld_sum
@@ -687 +694 @@
       !-------------------------------------
       jstart = 1
-      DO ib_bdy = 1, nb_bdy
-         jend = jstart - 1 + nb_bdy_fld(ib_bdy) 
-         CALL fld_fill( bf(jstart:jend), blf_i(jstart:jend), cn_dir_array(ib_bdy), 'bdy_dta',   &
+      DO jbdy = 1, nb_bdy
+         jend = jstart - 1 + nb_bdy_fld(jbdy) 
+         CALL fld_fill( bf(jstart:jend), blf_i(jstart:jend), cn_dir_array(jbdy), 'bdy_dta',   &
          &              'open boundary conditions', 'nambdy_dta' )
          jstart = jend + 1
@@ -700 +707 @@
 
       jfld = 0
-      DO ib_bdy=1, nb_bdy
-
-         nblen => idx_bdy(ib_bdy)%nblen
-         dta => dta_bdy(ib_bdy)
+      DO jbdy=1, nb_bdy
+
+         nblen => idx_bdy(jbdy)%nblen
+         dta => dta_bdy(jbdy)
 
          if(lwp) then
@@ -715 +722 @@
          endif
 
-         IF ( nn_dyn2d_dta(ib_bdy) == 0 .or. nn_dyn2d_dta(ib_bdy) == 2 ) THEN
+         IF ( nn_dyn2d_dta(jbdy) == 0 .or. nn_dyn2d_dta(jbdy) == 2 ) THEN
             if(lwp) write(numout,*) '++++++ dta%ssh/u2d/u3d allocated space'
             IF( dta%ll_ssh ) ALLOCATE( dta%ssh(nblen(1)) )
@@ -721 +728 @@
             IF( dta%ll_v2d ) ALLOCATE( dta%v2d(nblen(3)) )
          ENDIF
-         IF ( nn_dyn2d_dta(ib_bdy) == 1 .or. nn_dyn2d_dta(ib_bdy) == 3 ) THEN
+         IF ( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) THEN
             IF( dta%ll_ssh ) THEN
                if(lwp) write(numout,*) '++++++ dta%ssh pointing to fnow'
@@ -728 +735 @@
             ENDIF
             IF ( dta%ll_u2d ) THEN
-               IF ( ln_full_vel_array(ib_bdy) ) THEN
+               IF ( ln_full_vel_array(jbdy) ) THEN
                   if(lwp) write(numout,*) '++++++ dta%u2d allocated space'
                   ALLOCATE( dta%u2d(nblen(2)) )
@@ -738 +745 @@
             ENDIF
             IF ( dta%ll_v2d ) THEN
-               IF ( ln_full_vel_array(ib_bdy) ) THEN
+               IF ( ln_full_vel_array(jbdy) ) THEN
                   if(lwp) write(numout,*) '++++++ dta%v2d allocated space'
                   ALLOCATE( dta%v2d(nblen(3)) )
@@ -749 +756 @@
          ENDIF
 
-         IF ( nn_dyn3d_dta(ib_bdy) == 0 ) THEN
+         IF ( nn_dyn3d_dta(jbdy) == 0 ) THEN
             if(lwp) write(numout,*) '++++++ dta%u3d/v3d allocated space'
-            IF( dta%ll_u3d ) ALLOCATE( dta_bdy(ib_bdy)%u3d(nblen(2),jpk) )
-            IF( dta%ll_v3d ) ALLOCATE( dta_bdy(ib_bdy)%v3d(nblen(3),jpk) )
-         ENDIF
-         IF ( nn_dyn3d_dta(ib_bdy) == 1 .or. &
-           &  ( ln_full_vel_array(ib_bdy) .and. ( nn_dyn2d_dta(ib_bdy) == 1 .or. nn_dyn2d_dta(ib_bdy) == 3 ) ) ) THEN
-            IF ( dta%ll_u3d .or. ( ln_full_vel_array(ib_bdy) .and. dta%ll_u2d ) ) THEN
+            IF( dta%ll_u3d ) ALLOCATE( dta_bdy(jbdy)%u3d(nblen(2),jpk) )
+            IF( dta%ll_v3d ) ALLOCATE( dta_bdy(jbdy)%v3d(nblen(3),jpk) )
+         ENDIF
+         IF ( nn_dyn3d_dta(jbdy) == 1 .or. &
+           &  ( ln_full_vel_array(jbdy) .and. ( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) ) ) THEN
+            IF ( dta%ll_u3d .or. ( ln_full_vel_array(jbdy) .and. dta%ll_u2d ) ) THEN
                if(lwp) write(numout,*) '++++++ dta%u3d pointing to fnow'
                jfld = jfld + 1
-               dta_bdy(ib_bdy)%u3d => bf(jfld)%fnow(:,1,:)
-            ENDIF
-            IF ( dta%ll_v3d .or. ( ln_full_vel_array(ib_bdy) .and. dta%ll_v2d ) ) THEN
+               dta_bdy(jbdy)%u3d => bf(jfld)%fnow(:,1,:)
+            ENDIF
+            IF ( dta%ll_v3d .or. ( ln_full_vel_array(jbdy) .and. dta%ll_v2d ) ) THEN
                if(lwp) write(numout,*) '++++++ dta%v3d pointing to fnow'
                jfld = jfld + 1
-               dta_bdy(ib_bdy)%v3d => bf(jfld)%fnow(:,1,:)
-            ENDIF
-         ENDIF
-
-         IF( nn_tra_dta(ib_bdy) == 0 ) THEN
+               dta_bdy(jbdy)%v3d => bf(jfld)%fnow(:,1,:)
+            ENDIF
+         ENDIF
+
+         IF( nn_tra_dta(jbdy) == 0 ) THEN
             if(lwp) write(numout,*) '++++++ dta%tem/sal allocated space'
-            IF( dta%ll_tem ) ALLOCATE( dta_bdy(ib_bdy)%tem(nblen(1),jpk) )
-            IF( dta%ll_sal ) ALLOCATE( dta_bdy(ib_bdy)%sal(nblen(1),jpk) )
+            IF( dta%ll_tem ) ALLOCATE( dta_bdy(jbdy)%tem(nblen(1),jpk) )
+            IF( dta%ll_sal ) ALLOCATE( dta_bdy(jbdy)%sal(nblen(1),jpk) )
          ELSE
             IF( dta%ll_tem ) THEN
                if(lwp) write(numout,*) '++++++ dta%tem pointing to fnow'
                jfld = jfld + 1
-               dta_bdy(ib_bdy)%tem => bf(jfld)%fnow(:,1,:)
+               dta_bdy(jbdy)%tem => bf(jfld)%fnow(:,1,:)
             ENDIF
             IF( dta%ll_sal ) THEN 
                if(lwp) write(numout,*) '++++++ dta%sal pointing to fnow'
                jfld = jfld + 1
-               dta_bdy(ib_bdy)%sal => bf(jfld)%fnow(:,1,:)
-            ENDIF
-         ENDIF
-
-#if defined key_lim3
-         IF (cn_ice_lim(ib_bdy) /= 'none') THEN
-            IF( nn_ice_lim_dta(ib_bdy) == 0 ) THEN
-               ALLOCATE( dta_bdy(ib_bdy)%a_i(nblen(1),jpl) )
-               ALLOCATE( dta_bdy(ib_bdy)%h_i(nblen(1),jpl) )
-               ALLOCATE( dta_bdy(ib_bdy)%h_s(nblen(1),jpl) )
+               dta_bdy(jbdy)%sal => bf(jfld)%fnow(:,1,:)
+            ENDIF
+         ENDIF
+
+#if defined key_lim3
+         IF (cn_ice_lim(jbdy) /= 'none') THEN
+            IF( nn_ice_lim_dta(jbdy) == 0 ) THEN
+               ALLOCATE( dta_bdy(jbdy)%a_i(nblen(1),jpl) )
+               ALLOCATE( dta_bdy(jbdy)%h_i(nblen(1),jpl) )
+               ALLOCATE( dta_bdy(jbdy)%h_s(nblen(1),jpl) )
             ELSE
-               IF ( ll_bdylim3 ) THEN ! case input is Xcat
-                  jfld = jfld + 1
-                  dta_bdy(ib_bdy)%a_i => bf(jfld)%fnow(:,1,:)
-                  jfld = jfld + 1
-                  dta_bdy(ib_bdy)%h_i => bf(jfld)%fnow(:,1,:)
-                  jfld = jfld + 1
-                  dta_bdy(ib_bdy)%h_s => bf(jfld)%fnow(:,1,:)
-               ELSE ! case input is 1cat
+               IF ( nice_cat == jpl ) THEN ! case input cat = jpl
+                  jfld = jfld + 1
+                  dta_bdy(jbdy)%a_i => bf(jfld)%fnow(:,1,:)
+                  jfld = jfld + 1
+                  dta_bdy(jbdy)%h_i => bf(jfld)%fnow(:,1,:)
+                  jfld = jfld + 1
+                  dta_bdy(jbdy)%h_s => bf(jfld)%fnow(:,1,:)
+               ELSE                        ! case input cat = 1 OR (/=1 and /=jpl)
                   jfld_ai  = jfld + 1
                   jfld_hti = jfld + 2
                   jfld_hts = jfld + 3
                   jfld     = jfld + 3
-                  ALLOCATE( dta_bdy(ib_bdy)%a_i(nblen(1),jpl) )
-                  ALLOCATE( dta_bdy(ib_bdy)%h_i(nblen(1),jpl) )
-                  ALLOCATE( dta_bdy(ib_bdy)%h_s(nblen(1),jpl) )
-                  dta_bdy(ib_bdy)%a_i(:,:) = 0._wp
-                  dta_bdy(ib_bdy)%h_i(:,:) = 0._wp
-                  dta_bdy(ib_bdy)%h_s(:,:) = 0._wp
+                  ALLOCATE( dta_bdy(jbdy)%a_i(nblen(1),jpl) )
+                  ALLOCATE( dta_bdy(jbdy)%h_i(nblen(1),jpl) )
+                  ALLOCATE( dta_bdy(jbdy)%h_s(nblen(1),jpl) )
+                  dta_bdy(jbdy)%a_i(:,:) = 0._wp
+                  dta_bdy(jbdy)%h_i(:,:) = 0._wp
+                  dta_bdy(jbdy)%h_s(:,:) = 0._wp
                ENDIF
 
@@ -816 +823 @@
 #endif
          !
-      END DO ! ib_bdy 
+      END DO ! jbdy 
       !
       IF( nn_timing == 1 ) CALL timing_stop('bdy_dta_init')

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90

@@ -48 +48 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   alb_ice        !: ice albedo                                       [-]
 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qml_ice        !: heat available for snow / ice surface melting     [W/m2] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qcn_ice        !: heat conduction flux in the layer below surface   [W/m2] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qsr_ice_tr     !: solar flux transmitted below the ice surface      [W/m2]
+
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   utau_ice       !: atmos-ice u-stress. VP: I-pt ; EVP: U,V-pts   [N/m2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   vtau_ice       !: atmos-ice v-stress. VP: I-pt ; EVP: U,V-pts   [N/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr1_i0         !: Solar surface transmission parameter, thick ice  [-]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr2_i0         !: Solar surface transmission parameter, thin ice   [-]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   emp_ice        !: sublimation - precip over sea ice          [kg/m2/s]
 
@@ -123 +125 @@
       ALLOCATE( qns_ice (jpi,jpj,jpl) , qsr_ice (jpi,jpj,jpl) ,     &
          &      qla_ice (jpi,jpj,jpl) , dqla_ice(jpi,jpj,jpl) ,     &
-         &      dqns_ice(jpi,jpj,jpl) , tn_ice  (jpi,jpj,jpl) , alb_ice (jpi,jpj,jpl) ,   &
+         &      dqns_ice(jpi,jpj,jpl) , tn_ice  (jpi,jpj,jpl)  , alb_ice (jpi,jpj,jpl) ,   &
+         &      qml_ice(jpi,jpj,jpl)  , qcn_ice(jpi,jpj,jpl)   , qsr_ice_tr(jpi,jpj,jpl),  &
          &      utau_ice(jpi,jpj)     , vtau_ice(jpi,jpj)     , wndm_ice(jpi,jpj)     ,   &
-         &      fr1_i0  (jpi,jpj)     , fr2_i0  (jpi,jpj)     ,     &
          &      evap_ice(jpi,jpj,jpl) , devap_ice(jpi,jpj,jpl) , qprec_ice(jpi,jpj) ,   &
          &      qemp_ice(jpi,jpj)     , qevap_ice(jpi,jpj,jpl) , qemp_oce (jpi,jpj) ,   &
@@ -139 +141 @@
                 a_i(jpi,jpj,ncat)     , topmelt(jpi,jpj,ncat) , botmelt(jpi,jpj,ncat) , &
                 STAT= ierr(2) )
-      IF( ln_cpl )   ALLOCATE( u_ice(jpi,jpj)        , fr1_i0(jpi,jpj)       , tn_ice (jpi,jpj,1)    , &
-         &                     v_ice(jpi,jpj)        , fr2_i0(jpi,jpj)       , alb_ice(jpi,jpj,1)    , &
+      IF( ln_cpl )   ALLOCATE( u_ice(jpi,jpj)        , tn_ice (jpi,jpj,1)    , &
+         &                     v_ice(jpi,jpj)        , alb_ice(jpi,jpj,1)    , &
          &                     emp_ice(jpi,jpj)      , qns_ice(jpi,jpj,1)    , dqns_ice(jpi,jpj,1)   , &
          &                     STAT= ierr(3) )      
@@ -161 +163 @@
    PRIVATE
 
-   PUBLIC sbc_ice_alloc
+   PUBLIC   sbc_ice_alloc   !
 
    LOGICAL         , PUBLIC, PARAMETER ::   lk_lim3    = .FALSE.  !: no LIM-3 ice model
@@ -168 +170 @@
    REAL(wp)        , PUBLIC, PARAMETER ::   cldf_ice = 0.81       !: cloud fraction over sea ice, summer CLIO value   [-]
    INTEGER         , PUBLIC, PARAMETER ::   jpl = 1 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice,fr1_i0,fr2_i0          ! jpi, jpj
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice                        ! jpi, jpj
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tn_ice, alb_ice, qns_ice, dqns_ice  ! (jpi,jpj,jpl)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_i
@@ -176 +178 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   topmelt, botmelt
    !
-   !! arrays relating to embedding ice in the ocean. These arrays need to be declared 
-   !! even if no ice model is required. In the no ice model or traditional levitating 
-   !! ice cases they contain only zeros
+   !! arrays related to embedding ice in the ocean. 
+   !! These arrays need to be declared even if no ice model is required. 
+   !! In the no ice model or traditional levitating ice cases they contain only zeros
    !! ---------------------
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   snwice_mass        !: mass of snow and ice at current  ice time step   [Kg/m2]

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90

@@ -15 +15 @@
    !!            3.7  !  2014-06  (L. Brodeau)  simplification and optimization of CORE bulk
    !!            4.0  !  2016-06  (L. Brodeau)  sbcblk_core becomes sbcblk and is not restricted to the CORE algorithm anymore
-   !!                                          ==> based on AeroBulk (http://aerobulk.sourceforge.net/)
+   !!                 !                        ==> based on AeroBulk (http://aerobulk.sourceforge.net/)
    !!            4.0  !  2016-10  (G. Madec)  introduce a sbc_blk_init routine
+   !!            4.0  !  2016-10  (M. Vancoppenolle)  Introduce Jules emulator (M. Vancoppenolle) 
    !!----------------------------------------------------------------------
 
@@ -23 +24 @@
    !!   sbc_blk       : bulk formulation as ocean surface boundary condition
    !!   blk_oce       : computes momentum, heat and freshwater fluxes over ocean
-   !!   blk_ice       : computes momentum, heat and freshwater fluxes over sea ice
    !!   rho_air       : density of (moist) air (depends on T_air, q_air and SLP
    !!   cp_air        : specific heat of (moist) air (depends spec. hum. q_air)
    !!   q_sat         : saturation humidity as a function of SLP and temperature
    !!   L_vap         : latent heat of vaporization of water as a function of temperature
+   !!             sea-ice case only : 
+   !!   blk_ice_tau   : provide the air-ice stress
+   !!   blk_ice_flx   : provide the heat and mass fluxes at air-ice interface
+   !!   blk_ice_qcn   : provide ice surface temperature and snow/ice conduction flux (emulating JULES coupler)
+   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
+   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers
@@ -42 +48 @@
    USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, tm_su
    USE icethd_dh      ! for CALL ice_thd_snwblow
+   USE icethd_zdf,ONLY:   rn_cnd_s  ! for blk_ice_qcn
 #endif
    USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009) 
@@ -63 +70 @@
    PUBLIC   blk_ice_tau   ! routine called in icestp module
    PUBLIC   blk_ice_flx   ! routine called in icestp module
-#endif
+   PUBLIC   blk_ice_qcn   ! routine called in icestp module
+#endif 
 
 !!Lolo: should ultimately be moved in the module with all physical constants ?
@@ -111 +119 @@
    LOGICAL  ::   ln_Cd_L15 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2015)
    !
-!!cr   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_oce   ! air-ocean drag (clem)
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_atm                    ! transfer coefficient for momentum      (tau)
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ch_atm                    ! transfer coefficient for sensible heat (Q_sens)
@@ -139 +146 @@
       !!             ***  ROUTINE sbc_blk_alloc ***
       !!-------------------------------------------------------------------
-!!cr      ALLOCATE( Cd_oce(jpi,jpj) , STAT=sbc_blk_alloc )
       ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), &
          &      cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc )
@@ -147 +153 @@
    END FUNCTION sbc_blk_alloc
 
+
    SUBROUTINE sbc_blk_init
       !!---------------------------------------------------------------------
@@ -154 +161 @@
       !!
       !! ** Method  : 
-      !!
-      !!      C A U T I O N : never mask the surface stress fields
-      !!                      the stress is assumed to be in the (i,j) mesh referential
-      !!
-      !! ** Action  :   
       !!
       !!----------------------------------------------------------------------
@@ -285 +287 @@
       !!
       !! ** Purpose :   provide at each time step the surface ocean fluxes
-      !!      (momentum, heat, freshwater and runoff)
+      !!              (momentum, heat, freshwater and runoff)
       !!
       !! ** Method  : (1) READ each fluxes in NetCDF files:
@@ -566 +568 @@
    END SUBROUTINE blk_oce
 
+
+
+   FUNCTION rho_air( ptak, pqa, pslp )
+      !!-------------------------------------------------------------------------------
+      !!                           ***  FUNCTION rho_air  ***
+      !!
+      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
+      !!-------------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa]
+      REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air   ! density of moist air   [kg/m^3]
+      !!-------------------------------------------------------------------------------
+      !
+      rho_air = pslp / (  R_dry*ptak * ( 1._wp + rctv0*pqa )  )
+      !
+   END FUNCTION rho_air
+
+
+   FUNCTION cp_air( pqa )
+      !!-------------------------------------------------------------------------------
+      !!                           ***  FUNCTION cp_air  ***
+      !!
+      !! ** Purpose : Compute specific heat (Cp) of moist air
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
+      !!-------------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg]
+      REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air   ! specific heat of moist air   [J/K/kg]
+      !!-------------------------------------------------------------------------------
+      !
+      Cp_air = Cp_dry + Cp_vap * pqa
+      !
+   END FUNCTION cp_air
+
+
+   FUNCTION q_sat( ptak, pslp )
+      !!----------------------------------------------------------------------------------
+      !!                           ***  FUNCTION q_sat  ***
+      !!
+      !! ** Purpose : Specific humidity at saturation in [kg/kg]
+      !!              Based on accurate estimate of "e_sat"
+      !!              aka saturation water vapor (Goff, 1957)
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
+      !!----------------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak    ! air temperature                       [K]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp    ! sea level atmospheric pressure       [Pa]
+      REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg]
+      !
+      INTEGER  ::   ji, jj         ! dummy loop indices
+      REAL(wp) ::   ze_sat, ztmp   ! local scalar
+      !!----------------------------------------------------------------------------------
+      !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            !
+            ztmp = rt0 / ptak(ji,jj)
+            !
+            ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957)
+            ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0)        &
+               &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) )  &
+               &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614  )
+               !
+            q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat )   ! 0.01 because SLP is in [Pa]
+            !
+         END DO
+      END DO
+      !
+   END FUNCTION q_sat
+
+
+   FUNCTION gamma_moist( ptak, pqa )
+      !!----------------------------------------------------------------------------------
+      !!                           ***  FUNCTION gamma_moist  ***
+      !!
+      !! ** Purpose : Compute the moist adiabatic lapse-rate.
+      !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
+      !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
+      !!----------------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg]
+      REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist   ! moist adiabatic lapse-rate
+      !
+      INTEGER  ::   ji, jj         ! dummy loop indices
+      REAL(wp) :: zrv, ziRT        ! local scalar
+      !!----------------------------------------------------------------------------------
+      !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zrv = pqa(ji,jj) / (1. - pqa(ji,jj))
+            ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT
+            gamma_moist(ji,jj) = grav * ( 1. + cevap*zrv*ziRT ) / ( Cp_dry + cevap*cevap*zrv*reps0*ziRT/ptak(ji,jj) )
+         END DO
+      END DO
+      !
+   END FUNCTION gamma_moist
+
+
+   FUNCTION L_vap( psst )
+      !!---------------------------------------------------------------------------------
+      !!                           ***  FUNCTION L_vap  ***
+      !!
+      !! ** Purpose : Compute the latent heat of vaporization of water from temperature
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
+      !!----------------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K]
+      !!----------------------------------------------------------------------------------
+      !
+      L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6
+      !
+   END FUNCTION L_vap
+
 #if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3'                                       ESIM sea-ice model
+   !!----------------------------------------------------------------------
+   !!   blk_ice_tau : provide the air-ice stress
+   !!   blk_ice_flx : provide the heat and mass fluxes at air-ice interface
+   !!   blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating JULES coupler)
+   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
+   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 
+   !!----------------------------------------------------------------------
 
    SUBROUTINE blk_ice_tau
@@ -688 +818 @@
 
 
-   SUBROUTINE blk_ice_flx( ptsu, palb )
+   SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE blk_ice_flx  ***
       !!
-      !! ** Purpose :   provide the surface boundary condition over sea-ice
+      !! ** Purpose :   provide the heat and mass fluxes at air-ice interface
       !!
       !! ** Method  :   compute heat and freshwater exchanged
@@ -701 +831 @@
       !!---------------------------------------------------------------------
       REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu   ! sea ice surface temperature
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness
       REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb   ! ice albedo (all skies)
       !!
@@ -707 +839 @@
       REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      -
       REAL(wp) ::   zztmp, z1_lsub           !   -      -
+      REAL(wp) ::   zfrqsr_tr_i              ! sea ice shortwave fraction transmitted below through the ice
+      REAL(wp) ::   zfr1, zfr2, zfac         ! local variables
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qlw         ! long wave heat flux over ice
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qsb         ! sensible  heat flux over ice
@@ -717 +851 @@
       IF( ln_timing )   CALL timing_start('blk_ice_flx')
       !
-      !
-      ! local scalars
-      zcoef_dqlw   = 4.0 * 0.95 * Stef
-      zcoef_dqla   = -Ls * 11637800. * (-5897.8)
+      zcoef_dqlw = 4.0 * 0.95 * Stef             ! local scalars
+      zcoef_dqla = -Ls * 11637800. * (-5897.8)
       !
       zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
@@ -815 +947 @@
       END DO
 
-      !--------------------------------------------------------------------
-      ! FRACTIONs of net shortwave radiation which is not absorbed in the
-      ! thin surface layer and penetrates inside the ice cover
-      ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 )
-      !
-      fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
-      fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
+      ! --- absorbed and transmitted shortwave radiation (W/m2) --- !
+      !
+      ! former coding was
+      !     fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
+      !     fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
+
+      ! --- surface transmission parameter (i0, Grenfell Maykut 77) --- !
+      zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            !   standard value
+      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            !   zfr2 such that zfr1 + zfr2 to equal 1
+
+      qsr_ice_tr(:,:,:) = 0._wp
+
+      DO jl = 1, jpl
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               !
+               zfac = MAX( 0._wp , 1._wp - phi(ji,jj,jl) * 10._wp )     !   linear weighting factor: =0 for phi=0, =1 for phi = 0.1
+               zfrqsr_tr_i = zfr1 + zfac * zfr2                         !   below 10 cm, linearly increase zfrqsr_tr_i until 1 at zero thickness
+               !
+               IF ( phs(ji,jj,jl) <= 0._wp ) THEN   ;    zfrqsr_tr_i  = zfr1 + zfac * zfr2
+               ELSE                                 ;    zfrqsr_tr_i  = 0._wp                                  !   snow fully opaque
+               ENDIF
+               !
+               qsr_ice_tr(ji,jj,jl) = zfrqsr_tr_i * qsr_ice(ji,jj,jl)   !   transmitted solar radiation 
+               !
+            END DO
+         END DO
+      END DO
       !
       IF(ln_ctl) THEN
@@ -836 +989 @@
    END SUBROUTINE blk_ice_flx
    
-#endif
-
-   FUNCTION rho_air( ptak, pqa, pslp )
-      !!-------------------------------------------------------------------------------
-      !!                           ***  FUNCTION rho_air  ***
-      !!
-      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
-      !!-------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa]
-      REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air   ! density of moist air   [kg/m^3]
-      !!-------------------------------------------------------------------------------
-      !
-      rho_air = pslp / (  R_dry*ptak * ( 1._wp + rctv0*pqa )  )
-      !
-   END FUNCTION rho_air
-
-
-   FUNCTION cp_air( pqa )
-      !!-------------------------------------------------------------------------------
-      !!                           ***  FUNCTION cp_air  ***
-      !!
-      !! ** Purpose : Compute specific heat (Cp) of moist air
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
-      !!-------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg]
-      REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air   ! specific heat of moist air   [J/K/kg]
-      !!-------------------------------------------------------------------------------
-      !
-      Cp_air = Cp_dry + Cp_vap * pqa
-      !
-   END FUNCTION cp_air
-
-
-   FUNCTION q_sat( ptak, pslp )
-      !!----------------------------------------------------------------------------------
-      !!                           ***  FUNCTION q_sat  ***
-      !!
-      !! ** Purpose : Specific humidity at saturation in [kg/kg]
-      !!              Based on accurate estimate of "e_sat"
-      !!              aka saturation water vapor (Goff, 1957)
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
-      !!----------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak    ! air temperature                       [K]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp    ! sea level atmospheric pressure       [Pa]
-      REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg]
-      !
-      INTEGER  ::   ji, jj         ! dummy loop indices
-      REAL(wp) ::   ze_sat, ztmp   ! local scalar
-      !!----------------------------------------------------------------------------------
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            !
-            ztmp = rt0 / ptak(ji,jj)
-            !
-            ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957)
-            ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0)        &
-               &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) )  &
-               &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614  )
+
+   SUBROUTINE blk_ice_qcn( k_monocat, ptsu, ptb, phs, phi )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE blk_ice_qcn  ***
+      !!
+      !! ** Purpose :   Compute surface temperature and snow/ice conduction flux
+      !!                to force sea ice / snow thermodynamics
+      !!                in the case JULES coupler is emulated
+      !!                
+      !! ** Method  :   compute surface energy balance assuming neglecting heat storage
+      !!                following the 0-layer Semtner (1976) approach
+      !!
+      !! ** Outputs : - ptsu    : sea-ice / snow surface temperature (K)
+      !!              - qcn_ice : surface inner conduction flux (W/m2)
+      !!
+      !!---------------------------------------------------------------------
+      INTEGER                   , INTENT(in   ) ::   k_monocat   ! single-category option
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptsu        ! sea ice / snow surface temperature
+      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   ptb         ! sea ice base temperature
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phs         ! snow thickness
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phi         ! sea ice thickness
+      !
+      INTEGER , PARAMETER ::   nit = 10                  ! number of iterations
+      REAL(wp), PARAMETER ::   zepsilon = 0.1_wp         ! characteristic thickness for enhanced conduction
+      !
+      INTEGER  ::   ji, jj, jl           ! dummy loop indices
+      INTEGER  ::   iter                 ! local integer
+      REAL(wp) ::   zfac, zfac2, zfac3   ! local scalars
+      REAL(wp) ::   zkeff_h, ztsu        !
+      REAL(wp) ::   zqc, zqnet           !
+      REAL(wp) ::   zhe, zqa0            !
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zgfac   ! enhanced conduction factor
+      !!---------------------------------------------------------------------
+
+      IF( nn_timing == 1 )  CALL timing_start('blk_ice_qcn')
+      
+      ! -------------------------------------!
+      !      I   Enhanced conduction factor  !
+      ! -------------------------------------!
+      ! Emulates the enhancement of conduction by unresolved thin ice (k_monocat = 1/3)
+      ! Fichefet and Morales Maqueda, JGR 1997
+      !
+      zgfac(:,:,:) = 1._wp
+      
+      SELECT CASE ( k_monocat ) 
+      !
+      CASE ( 1 , 3 )
+         !
+         zfac    =   1._wp /  ( rn_cnd_s + rcdic )
+         zfac2   =   EXP(1._wp) * 0.5_wp * zepsilon
+         zfac3   =   2._wp / zepsilon
+         !   
+         DO jl = 1, jpl                
+            DO jj = 1 , jpj
+               DO ji = 1, jpi
+                  !                                ! Effective thickness
+                  zhe               =   ( rn_cnd_s * phi(ji,jj,jl) + rcdic * phs(ji,jj,jl) ) * zfac
+                  !                                ! Enhanced conduction factor
+                  IF( zhe >=  zfac2 )  &
+                     zgfac(ji,jj,jl)   =   MIN( 2._wp, ( 0.5_wp + 0.5 * LOG( zhe * zfac3 ) ) )
+               END DO
+            END DO
+         END DO
+         !      
+      END SELECT
+      
+      ! -------------------------------------------------------------!
+      !      II   Surface temperature and conduction flux            !
+      ! -------------------------------------------------------------!
+      !
+      zfac   =   rcdic * rn_cnd_s 
+      !                             ! ========================== !
+      DO jl = 1, jpl                !  Loop over ice categories  !
+         !                          ! ========================== !
+         DO jj = 1 , jpj
+            DO ji = 1, jpi
+               !                    ! Effective conductivity of the snow-ice system divided by thickness
+               zkeff_h =   zfac * zgfac(ji,jj,jl) / ( rcdic * phs(ji,jj,jl) + rn_cnd_s * phi(ji,jj,jl) )
+               !                    ! Store initial surface temperature
+               ztsu    =   ptsu(ji,jj,jl)
+               !                    ! Net initial atmospheric heat flux
+               zqa0    =   qsr_ice(ji,jj,jl) - qsr_ice_tr(ji,jj,jl) + qns_ice(ji,jj,jl)
                !
-            q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat )   ! 0.01 because SLP is in [Pa]
-            !
-         END DO
-      END DO
-      !
-   END FUNCTION q_sat
-
-
-   FUNCTION gamma_moist( ptak, pqa )
-      !!----------------------------------------------------------------------------------
-      !!                           ***  FUNCTION gamma_moist  ***
-      !!
-      !! ** Purpose : Compute the moist adiabatic lapse-rate.
-      !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
-      !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
-      !!----------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg]
-      REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist   ! moist adiabatic lapse-rate
-      !
-      INTEGER  ::   ji, jj         ! dummy loop indices
-      REAL(wp) :: zrv, ziRT        ! local scalar
-      !!----------------------------------------------------------------------------------
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zrv = pqa(ji,jj) / (1. - pqa(ji,jj))
-            ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT
-            gamma_moist(ji,jj) = grav * ( 1. + cevap*zrv*ziRT ) / ( Cp_dry + cevap*cevap*zrv*reps0*ziRT/ptak(ji,jj) )
-         END DO
-      END DO
-      !
-   END FUNCTION gamma_moist
-
-
-   FUNCTION L_vap( psst )
-      !!---------------------------------------------------------------------------------
-      !!                           ***  FUNCTION L_vap  ***
-      !!
-      !! ** Purpose : Compute the latent heat of vaporization of water from temperature
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
-      !!----------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K]
-      !!----------------------------------------------------------------------------------
-      !
-      L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6
-      !
-   END FUNCTION L_vap
-
-#if defined key_lim3
+               DO iter = 1, nit     ! --- Iteration loop
+                   !                      ! Conduction heat flux through snow-ice system (>0 downwards)
+                   zqc   =   zkeff_h * ( ztsu - ptb(ji,jj) )
+                   !                      ! Surface energy budget
+                   zqnet =   zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc
+                   !                      ! Temperature update
+                   ztsu  =   ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h )
+               END DO
+               !
+               ptsu   (ji,jj,jl) = MIN( rt0, ztsu )
+               qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) )
+               !
+            END DO
+         END DO
+         !
+      END DO 
+      !      
+      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_qcn')
+      !
+   END SUBROUTINE blk_ice_qcn
+   
 
    SUBROUTINE Cdn10_Lupkes2012( Cd )
@@ -961 +1098 @@
       !!                      ***  ROUTINE  Cdn10_Lupkes2012  ***
       !!
-      !! ** Purpose :    Recompute the ice-atm drag at 10m height to make
-      !!                 it dependent on edges at leads, melt ponds and flows.
+      !! ** Purpose :    Recompute the neutral air-ice drag referenced at 10m 
+      !!                 to make it dependent on edges at leads, melt ponds and flows.
       !!                 After some approximations, this can be resumed to a dependency
       !!                 on ice concentration.
@@ -1012 +1149 @@
       !!                      ***  ROUTINE  Cdn10_Lupkes2015  ***
       !!
-      !! ** pUrpose :    1lternative turbulent transfert coefficients formulation
+      !! ** pUrpose :    Alternative turbulent transfert coefficients formulation
       !!                 between sea-ice and atmosphere with distinct momentum 
       !!                 and heat coefficients depending on sea-ice concentration 

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_coare.F90

@@ -49 +49 @@
    PUBLIC ::   TURB_COARE   ! called by sbcblk.F90
 
-   !                                           !! COARE own values for given constants:
-   REAL(wp), PARAMETER ::   zi0     = 600.      ! scale height of the atmospheric boundary layer...1
-   REAL(wp), PARAMETER ::   Beta0   =   1.25    ! gustiness parameter
-   REAL(wp), PARAMETER ::   rctv0   =   0.608   ! constant to obtain virtual temperature...
+   !                                              !! COARE own values for given constants:
+   REAL(wp), PARAMETER ::   zi0     = 600._wp      ! scale height of the atmospheric boundary layer...
+   REAL(wp), PARAMETER ::   Beta0   =   1.250_wp   ! gustiness parameter
+   REAL(wp), PARAMETER ::   rctv0   =   0.608_wp   ! constant to obtain virtual temperature...
 
    !!----------------------------------------------------------------------

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

@@ -971 +971 @@
       !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case)
       !!----------------------------------------------------------------------
-      USE zdf_oce,  ONLY : ln_zdfswm
-
-      IMPLICIT NONE
-
-      INTEGER, INTENT(in)           ::   kt          ! ocean model time step index
-      INTEGER, INTENT(in)           ::   k_fsbc      ! frequency of sbc (-> ice model) computation 
-      INTEGER, INTENT(in)           ::   k_ice       ! ice management in the sbc (=0/1/2/3)
+      USE zdf_oce,  ONLY :   ln_zdfswm
+      !
+      INTEGER, INTENT(in) ::   kt          ! ocean model time step index
+      INTEGER, INTENT(in) ::   k_fsbc      ! frequency of sbc (-> ice model) computation 
+      INTEGER, INTENT(in) ::   k_ice       ! ice management in the sbc (=0/1/2/3)
       !!
       LOGICAL  ::   llnewtx, llnewtau      ! update wind stress components and module??
@@ -1170 +1168 @@
          ! Calculate the 3D Stokes drift both in coupled and not fully uncoupled mode
          IF( srcv(jpr_sdrftx)%laction .OR. srcv(jpr_sdrfty)%laction .OR. srcv(jpr_wper)%laction &
-                                                                    .OR. srcv(jpr_hsig)%laction ) THEN
+            &                                                       .OR. srcv(jpr_hsig)%laction ) THEN
             CALL sbc_stokes()
          ENDIF
@@ -1525 +1523 @@
    
 
-   SUBROUTINE sbc_cpl_ice_flx( picefr, palbi, psst, pist )
+   SUBROUTINE sbc_cpl_ice_flx( picefr, palbi, psst, pist, phs, phi )
       !!----------------------------------------------------------------------
       !!             ***  ROUTINE sbc_cpl_ice_flx  ***
@@ -1576 +1574 @@
       !!----------------------------------------------------------------------
       REAL(wp), INTENT(in), DIMENSION(:,:)             ::   picefr     ! ice fraction                [0 to 1]
-      ! optional arguments, used only in 'mixed oce-ice' case
+      !                                                !!           ! optional arguments, used only in 'mixed oce-ice' case
       REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   palbi      ! all skies ice albedo 
       REAL(wp), INTENT(in), DIMENSION(:,:  ), OPTIONAL ::   psst       ! sea surface temperature     [Celsius]
       REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   pist       ! ice surface temperature     [Kelvin]
-      !
-      INTEGER ::   jl         ! dummy loop index
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice
+      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   phs        ! snow depth                  [m]
+      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   phi        ! ice thickness               [m]
+      !
+      INTEGER  ::   ji, jj, jl   ! dummy loop index
+      REAL(wp) ::   ztri         ! local scalar
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip  , zevap_oce, zevap_ice, zdevap_ice
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice    !!gm , zfrqsr_tr_i
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_flx')
       !
-      CALL wrk_alloc( jpi,jpj,     zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw )
-      CALL wrk_alloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice )
-      CALL wrk_alloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice )
-      CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice )
-
       IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
-      ziceld(:,:) = 1. - picefr(:,:)
-      zcptn(:,:) = rcp * sst_m(:,:)
+      ziceld(:,:) = 1._wp - picefr(:,:)
+      zcptn (:,:) = rcp * sst_m(:,:)
       !
       !                                                      ! ========================= !
@@ -1622 +1618 @@
 #if defined key_lim3
       ! zsnw = snow fraction over ice after wind blowing (=picefr if no blowing)
-      zsnw(:,:) = 0._wp  ;  CALL ice_thd_snwblow( ziceld, zsnw )
+      zsnw(:,:) = 0._wp   ;   CALL ice_thd_snwblow( ziceld, zsnw )
       
       ! --- evaporation minus precipitation corrected (because of wind blowing on snow) --- !
@@ -1659 +1655 @@
          sprecip(:,:) = sprecip(:,:) * xcplmask(:,:,0) + zsprecip(:,:) * zmsk(:,:)
          tprecip(:,:) = tprecip(:,:) * xcplmask(:,:,0) + ztprecip(:,:) * zmsk(:,:)
-         DO jl=1,jpl
+         DO jl = 1, jpl
             evap_ice (:,:,jl) = evap_ice (:,:,jl) * xcplmask(:,:,0) + zevap_ice (:,:) * zmsk(:,:)
             devap_ice(:,:,jl) = devap_ice(:,:,jl) * xcplmask(:,:,0) + zdevap_ice(:,:) * zmsk(:,:)
-         ENDDO
+         END DO
       ELSE
-         emp_tot(:,:) =         zemp_tot(:,:)
-         emp_ice(:,:) =         zemp_ice(:,:)
-         emp_oce(:,:) =         zemp_oce(:,:)     
-         sprecip(:,:) =         zsprecip(:,:)
-         tprecip(:,:) =         ztprecip(:,:)
-         DO jl=1,jpl
+         emp_tot(:,:) = zemp_tot(:,:)
+         emp_ice(:,:) = zemp_ice(:,:)
+         emp_oce(:,:) = zemp_oce(:,:)     
+         sprecip(:,:) = zsprecip(:,:)
+         tprecip(:,:) = ztprecip(:,:)
+         DO jl = 1, jpl
             evap_ice (:,:,jl) = zevap_ice (:,:)
             devap_ice(:,:,jl) = zdevap_ice(:,:)
-         ENDDO
+         END DO
       ENDIF
 
@@ -1691 +1687 @@
         fwfisf(:,:) = - frcv(jpr_isf)%z3(:,:,1)
       ENDIF
-
+      !
       IF( ln_mixcpl ) THEN
          emp_tot(:,:) = emp_tot(:,:) * xcplmask(:,:,0) + zemp_tot(:,:) * zmsk(:,:)
@@ -1703 +1699 @@
          tprecip(:,:) =                                  ztprecip(:,:)
       ENDIF
-
+      !
 #endif
+
       ! outputs
 !!      IF( srcv(jpr_rnf)%laction )   CALL iom_put( 'runoffs' , rnf(:,:) * tmask(:,:,1)                                 )  ! runoff
@@ -1730 +1727 @@
             zqns_ice(:,:,1:jpl) = frcv(jpr_qnsice)%z3(:,:,1:jpl)
          ELSE
-            DO jl=1,jpl
+            DO jl = 1, jpl
                zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1) ! Set all category values equal
-            ENDDO
+            END DO
          ENDIF
       CASE( 'oce and ice' )      ! the total flux is computed from ocean and ice fluxes
@@ -1743 +1740 @@
          ELSE
             qns_tot(:,:) = qns_tot(:,:) + picefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
-            DO jl=1,jpl
+            DO jl = 1, jpl
                zqns_tot(:,:   ) = zqns_tot(:,:) + picefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
                zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
-            ENDDO
+            END DO
          ENDIF
       CASE( 'mixed oce-ice' )    ! the ice flux is cumputed from the total flux, the SST and ice informations
@@ -1766 +1763 @@
       !         note: ziceld cannot be = 0 since we limit the ice concentration to amax
       zqns_oce = 0._wp
-      WHERE( ziceld /= 0._wp )  zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / ziceld(:,:)
+      WHERE( ziceld /= 0._wp )   zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / ziceld(:,:)
 
       ! Heat content per unit mass of snow (J/kg)
@@ -1859 +1856 @@
          ELSE
             ! Set all category values equal for the moment
-            DO jl=1,jpl
+            DO jl = 1, jpl
                zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
-            ENDDO
+            END DO
          ENDIF
          zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
@@ -1868 +1865 @@
          zqsr_tot(:,:  ) =  ziceld(:,:) * frcv(jpr_qsroce)%z3(:,:,1)
          IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
-            DO jl=1,jpl
+            DO jl = 1, jpl
                zqsr_tot(:,:   ) = zqsr_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qsrice)%z3(:,:,jl)   
                zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,jl)
-            ENDDO
+            END DO
          ELSE
             qsr_tot(:,:   ) = qsr_tot(:,:) + picefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
-            DO jl=1,jpl
+            DO jl = 1, jpl
                zqsr_tot(:,:   ) = zqsr_tot(:,:) + picefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
                zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
-            ENDDO
+            END DO
          ENDIF
       CASE( 'mixed oce-ice' )
@@ -1890 +1887 @@
       IF( ln_dm2dc .AND. ln_cpl ) THEN   ! modify qsr to include the diurnal cycle
          zqsr_tot(:,:  ) = sbc_dcy( zqsr_tot(:,:  ) )
-         DO jl=1,jpl
+         DO jl = 1, jpl
             zqsr_ice(:,:,jl) = sbc_dcy( zqsr_ice(:,:,jl) )
-         ENDDO
+         END DO
       ENDIF
 
@@ -1908 +1905 @@
          qsr_tot(:,:) = qsr(:,:) * ziceld(:,:) + SUM( qsr_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk
          qsr_tot(:,:) = qsr_tot(:,:) * xcplmask(:,:,0) +  zqsr_tot(:,:)* zmsk(:,:)
-         DO jl=1,jpl
+         DO jl = 1, jpl
             qsr_ice(:,:,jl) = qsr_ice(:,:,jl) * xcplmask(:,:,0) +  zqsr_ice(:,:,jl)* zmsk(:,:)
-         ENDDO
+         END DO
       ELSE
          qsr_tot(:,:  ) = zqsr_tot(:,:  )
@@ -1946 +1943 @@
       END SELECT
 
-      ! Surface transimission parameter io (Maykut Untersteiner , 1971 ; Ebert and Curry, 1993 )
-      ! Used for LIM3
-      ! Coupled case: since cloud cover is not received from atmosphere 
-      !               ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
-      fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
-      fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
-
-      CALL wrk_dealloc( jpi,jpj,     zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw )
-      CALL wrk_dealloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice )
-      CALL wrk_dealloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice )
-      CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice )
+      !                                                      ! ========================= !
+      !                                                      !      Transmitted Qsr      !   [W/m2]
+      !                                                      ! ========================= !
+      SELECT CASE( nice_jules )
+      CASE( np_jules_OFF    )       !==  No Jules coupler  ==!
+         !
+!!gm         ! former coding was    
+!!gm         ! Coupled case: since cloud cover is not received from atmosphere 
+!!gm         !               ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
+!!gm         !     fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
+!!gm         !     fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
+!!gm         
+!!gm         ! to retrieve that coding, we needed to access h_i & h_s from here
+!!gm         ! we could even retrieve cloud fraction from the coupler
+!!gm         !
+!!gm         zfrqsr_tr_i(:,:,:) = 0._wp   !   surface transmission parameter
+!!gm         !
+!!gm         DO jl = 1, jpl
+!!gm            DO jj = 1 , jpj
+!!gm               DO ji = 1, jpi
+!!gm                  !              !--- surface transmission parameter (Grenfell Maykut 77) --- !
+!!gm                  zfrqsr_tr_i(ji,jj,jl) = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice 
+!!gm                  !
+!!gm                  !              ! --- influence of snow and thin ice --- !
+!!gm                  IF ( phs(ji,jj,jl) >= 0.0_wp )   zfrqsr_tr_i(ji,jj,jl) = 0._wp   !   snow fully opaque
+!!gm                  IF ( phi(ji,jj,jl) <= 0.1_wp )   zfrqsr_tr_i(ji,jj,jl) = 1._wp   !   thin ice transmits all solar radiation
+!!gm               END DO
+!!gm            END DO
+!!gm         END DO
+!!gm         !
+!!gm         qsr_ice_tr(:,:,:) =   zfrqsr_tr_i(:,:,:) * qsr_ice(:,:,:)               !   transmitted solar radiation 
+!!gm         !
+!!gm better coding of the above calculation:
+         !
+         !                    ! ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
+         ztri = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice    ! surface transmission parameter (Grenfell Maykut 77)
+         !
+         qsr_ice_tr(:,:,:) = ztri * qsr_ice(:,:,:)
+         WHERE( phs(:,:,:) >= 0.0_wp )   qsr_ice_tr(:,:,:) = 0._wp            ! snow fully opaque
+         WHERE( phi(:,:,:) <= 0.1_wp )   qsr_ice_tr(:,:,:) = qsr_ice(:,:,:)   ! thin ice transmits all solar radiation
+!!gm end
+         !     
+      CASE( np_jules_ACTIVE )       !==  Jules coupler is active  ==!
+         !
+         !                    ! ===> here we must receive the qsr_ice_tr array from the coupler
+         !                           for now just assume zero (fully opaque ice)
+         qsr_ice_tr(:,:,:) = 0._wp
+         !
+      END SELECT
       !
       IF( nn_timing == 1 )   CALL timing_stop('sbc_cpl_ice_flx')
@@ -2081 +2116 @@
       IF( ssnd(jps_albmix)%laction ) THEN                         ! mixed ice-ocean
          ztmp1(:,:) = alb_oce_mix(:,:) * zfr_l(:,:)
-         DO jl=1,jpl
+         DO jl = 1, jpl
             ztmp1(:,:) = ztmp1(:,:) + alb_ice(:,:,jl) * a_i(:,:,jl)
          END DO

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90

@@ -710 +710 @@
       !                                         !   Prepare Coupling fields   !
       !                                         ! =========================== !
-
-! x and y comp of ice velocity
-
+      !
+      ! x and y comp of ice velocity
+      !
       CALL cice2nemo(uvel,u_ice,'F', -1. )
       CALL cice2nemo(vvel,v_ice,'F', -1. )
-
-! Ice concentration (CO_1) = a_i calculated at end of cice_sbc_out  
-
-! Snow and ice thicknesses (CO_2 and CO_3)
-
-      DO jl = 1,ncat
+      !
+      ! Ice concentration (CO_1) = a_i calculated at end of cice_sbc_out  
+      !
+      ! Snow and ice thicknesses (CO_2 and CO_3)
+      !
+      DO jl = 1, ncat
          CALL cice2nemo( vsnon(:,:,jl,:), h_s(:,:,jl),'T', 1. )
          CALL cice2nemo( vicen(:,:,jl,:), h_i(:,:,jl),'T', 1. )
-      ENDDO
+      END DO
       !
       IF( nn_timing == 1 )  CALL timing_stop('cice_sbc_hadgam')

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90

@@ -21 +21 @@
    !
    USE in_out_manager ! I/O manager
-   USE iom            ! I/O manager library
+   USE iom            ! I/O library
    USE fldread        ! read input field at current time step
    USE lbclnk         !
@@ -41 +41 @@
    REAL(wp), PUBLIC ::   rn_gammas0                  !: salinity    exchange coeficient []
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   risf_tsc_b, risf_tsc  !: before and now T & S isf contents [K.m/s & PSU.m/s]  
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf                  !: net heat flux from ice shelf      [W/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  rzisf_tbl              !:depth of calving front (shallowest point) nn_isf ==2/3
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  rhisf_tbl, rhisf_tbl_0 !:thickness of tbl  [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  r1_hisf_tbl            !:1/thickness of tbl
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  ralpha                 !:proportion of bottom cell influenced by tbl 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  risfLeff               !:effective length (Leff) BG03 nn_isf==2
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  ttbl, stbl, utbl, vtbl !:top boundary layer variable at T point
-   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  misfkt, misfkb         !:Level of ice shelf base
-
-   LOGICAL, PUBLIC ::   l_isfcpl = .false.       ! isf recieved from oasis
-
-   REAL(wp), PUBLIC, SAVE ::   rcpi     = 2000.0_wp     ! specific heat of ice shelf             [J/kg/K]
-   REAL(wp), PUBLIC, SAVE ::   rkappa   = 1.54e-6_wp    ! heat diffusivity through the ice-shelf [m2/s]
-   REAL(wp), PUBLIC, SAVE ::   rhoisf   = 920.0_wp      ! volumic mass of ice shelf              [kg/m3]
-   REAL(wp), PUBLIC, SAVE ::   tsurf    = -20.0_wp      ! air temperature on top of ice shelf    [C]
-   REAL(wp), PUBLIC, SAVE ::   rlfusisf = 0.334e6_wp    ! latent heat of fusion of ice shelf     [J/kg]
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   misfkt   , misfkb        !: Level of ice shelf base
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rzisf_tbl                !: depth of calving front (shallowest point) nn_isf ==2/3
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rhisf_tbl, rhisf_tbl_0   !: thickness of tbl  [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   r1_hisf_tbl              !: 1/thickness of tbl
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ralpha                   !: proportion of bottom cell influenced by tbl 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   risfLeff                 !: effective length (Leff) BG03 nn_isf==2
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ttbl, stbl, utbl, vtbl   !: top boundary layer variable at T point
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf                     !: net heat flux from ice shelf      [W/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   risf_tsc_b, risf_tsc     !: before and now T & S isf contents [K.m/s & PSU.m/s]  
+
+   LOGICAL, PUBLIC ::   l_isfcpl = .false.       !: isf recieved from oasis
+
+   REAL(wp), PUBLIC, SAVE ::   rcpi     = 2000.0_wp     !: specific heat of ice shelf             [J/kg/K]
+   REAL(wp), PUBLIC, SAVE ::   rkappa   = 1.54e-6_wp    !: heat diffusivity through the ice-shelf [m2/s]
+   REAL(wp), PUBLIC, SAVE ::   rhoisf   = 920.0_wp      !: volumic mass of ice shelf              [kg/m3]
+   REAL(wp), PUBLIC, SAVE ::   tsurf    = -20.0_wp      !: air temperature on top of ice shelf    [C]
+   REAL(wp), PUBLIC, SAVE ::   rlfusisf = 0.334e6_wp    !: latent heat of fusion of ice shelf     [J/kg]
 
 !: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3)
@@ -70 +70 @@
    
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.7 , LOCEAN-IPSL (2015)
+   !! NEMO/OPA 4.0 , LOCEAN-IPSL (2017)
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -91 +91 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time step
       !
-      INTEGER :: ji, jj, jk   ! loop index
-      INTEGER :: ikt, ikb     ! local integers
+      INTEGER ::   ji, jj, jk   ! loop index
+      INTEGER ::   ikt, ikb     ! local integers
       REAL(wp), DIMENSION(jpi,jpj) ::   zt_frz, zdep   ! freezing temperature (zt_frz) at depth (zdep) 
       REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   zqhcisf2d