Index: branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 5260)
+++ branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 5989)
@@ -37,4 +37,5 @@
    USE limdyn          ! Ice dynamics
    USE limtrp          ! Ice transport
+   USE limhdf          ! Ice horizontal diffusion
    USE limthd          ! Ice thermodynamics
    USE limitd_me       ! Mechanics on ice thickness distribution
@@ -110,40 +111,84 @@
       REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky
       REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice          ! mean ice albedo (for coupled)
+      REAL(wp), POINTER, DIMENSION(:,:  )   ::   zutau_ice, zvtau_ice 
       !!----------------------------------------------------------------------
 
       IF( nn_timing == 1 )  CALL timing_start('sbc_ice_lim')
 
-      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !  Ice time-step only
-         !-----------------------!                                           
-         ! --- Bulk Formulae --- !                                           
-         !-----------------------!
-         u_oce(:,:) = ssu_m(:,:) * umask(:,:,1)      ! mean surface ocean current at ice velocity point
-         v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1)      ! (C-grid dynamics :  U- & V-points as the ocean)
+      !-----------------------!
+      ! --- Ice time step --- !
+      !-----------------------!
+      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
+
+         ! mean surface ocean current at ice velocity point (C-grid dynamics :  U- & V-points as the ocean)
+         u_oce(:,:) = ssu_m(:,:) * umask(:,:,1)
+         v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1)
          
          ! masked sea surface freezing temperature [Kelvin] (set to rt0 over land)
-         t_bo(:,:) = ( eos_fzp( sss_m ) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )  
-         !                                                                                      
-         ! Ice albedo
-         CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )
-         CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os )  ! cloud-sky and overcast-sky ice albedos
-
-         ! CORE and COUPLED bulk formulations
-         SELECT CASE( kblk )
-         CASE( jp_core , jp_cpl )
-
-            ! albedo depends on cloud fraction because of non-linear spectral effects
-            zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
-            ! In CLIO the cloud fraction is read in the climatology and the all-sky albedo 
-            ! (zalb_ice) is computed within the bulk routine
-            
-         END SELECT
-         
+         CALL eos_fzp( sss_m(:,:) , t_bo(:,:) )
+         t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
+          
          ! Mask sea ice surface temperature (set to rt0 over land)
          DO jl = 1, jpl
             t_su(:,:,jl) = t_su(:,:,jl) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
-         END DO
-     
-         ! Bulk formulae  - provides the following fields:
-         ! utau_ice, vtau_ice : surface ice stress                     (U- & V-points)   [N/m2]
+         END DO     
+         !
+         !------------------------------------------------!                                           
+         ! --- Dynamical coupling with the atmosphere --- !                                           
+         !------------------------------------------------!
+         ! It provides the following fields:
+         ! utau_ice, vtau_ice : surface ice stress (U- & V-points)   [N/m2]
+         !-----------------------------------------------------------------
+         SELECT CASE( kblk )
+         CASE( jp_clio    )   ;   CALL blk_ice_clio_tau                         ! CLIO bulk formulation            
+         CASE( jp_core    )   ;   CALL blk_ice_core_tau                         ! CORE bulk formulation
+         CASE( jp_purecpl )   ;   CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )   ! Coupled   formulation
+         END SELECT
+         
+         IF( ln_mixcpl) THEN   ! Case of a mixed Bulk/Coupled formulation
+            CALL wrk_alloc( jpi,jpj    , zutau_ice, zvtau_ice)
+            CALL sbc_cpl_ice_tau( zutau_ice , zvtau_ice )
+            utau_ice(:,:) = utau_ice(:,:) * xcplmask(:,:,0) + zutau_ice(:,:) * ( 1. - xcplmask(:,:,0) )
+            vtau_ice(:,:) = vtau_ice(:,:) * xcplmask(:,:,0) + zvtau_ice(:,:) * ( 1. - xcplmask(:,:,0) )
+            CALL wrk_dealloc( jpi,jpj  , zutau_ice, zvtau_ice)
+         ENDIF
+
+         !-------------------------------------------------------!
+         ! --- ice dynamics and transport (except in 1D case) ---!
+         !-------------------------------------------------------!
+         numit = numit + nn_fsbc                  ! Ice model time step
+         !                                                   
+         CALL sbc_lim_bef                         ! Store previous ice values
+         CALL sbc_lim_diag0                       ! set diag of mass, heat and salt fluxes to 0
+         CALL lim_rst_opn( kt )                   ! Open Ice restart file
+         !
+         IF( .NOT. lk_c1d ) THEN
+            !
+            CALL lim_dyn( kt )                    ! Ice dynamics    ( rheology/dynamics )   
+            !
+            CALL lim_trp( kt )                    ! Ice transport   ( Advection/diffusion )
+            !
+            IF( nn_monocat /= 2 ) CALL lim_itd_me ! Mechanical redistribution ! (ridging/rafting)
+            !
+#if defined key_bdy
+            CALL bdy_ice_lim( kt )                ! bdy ice thermo 
+            IF( ln_icectl )       CALL lim_prt( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )
+#endif
+            !
+            CALL lim_update1( kt )                ! Corrections
+            !
+         ENDIF
+         
+         ! previous lead fraction and ice volume for flux calculations
+         CALL sbc_lim_bef                        
+         CALL lim_var_glo2eqv                     ! ht_i and ht_s for ice albedo calculation
+         CALL lim_var_agg(1)                      ! at_i for coupling (via pfrld) 
+         pfrld(:,:)   = 1._wp - at_i(:,:)
+         phicif(:,:)  = vt_i(:,:)
+         
+         !------------------------------------------------------!                                           
+         ! --- Thermodynamical coupling with the atmosphere --- !                                           
+         !------------------------------------------------------!
+         ! It provides the following fields:
          ! qsr_ice , qns_ice  : solar & non solar heat flux over ice   (T-point)         [W/m2]
          ! qla_ice            : latent heat flux over ice              (T-point)         [W/m2]
@@ -151,85 +196,37 @@
          ! tprecip , sprecip  : total & solid precipitation            (T-point)         [Kg/m2/s]
          ! fr1_i0  , fr2_i0   : 1sr & 2nd fraction of qsr penetration in ice             [%]
-         !
+         !----------------------------------------------------------------------------------------
+         CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )
+         CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos
+
          SELECT CASE( kblk )
          CASE( jp_clio )                                       ! CLIO bulk formulation
-            CALL blk_ice_clio( t_su , zalb_cs    , zalb_os    , zalb_ice  ,               &
-               &                      utau_ice   , vtau_ice   , qns_ice   , qsr_ice   ,   &
-               &                      qla_ice    , dqns_ice   , dqla_ice  ,               &
-               &                      tprecip    , sprecip    ,                           &
-               &                      fr1_i0     , fr2_i0     , cp_ice_msh, jpl  )
-            !         
-            IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   &
-               &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx )
-
+            ! In CLIO the cloud fraction is read in the climatology and the all-sky albedo 
+            ! (zalb_ice) is computed within the bulk routine
+            CALL blk_ice_clio_flx( t_su, zalb_cs, zalb_os, zalb_ice )
+            IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su )
+            IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
          CASE( jp_core )                                       ! CORE bulk formulation
-            CALL blk_ice_core( t_su , u_ice     , v_ice     , zalb_ice   ,               &
-               &                      utau_ice  , vtau_ice  , qns_ice    , qsr_ice   ,   &
-               &                      qla_ice   , dqns_ice  , dqla_ice   ,               &
-               &                      tprecip   , sprecip   ,                            &
-               &                      fr1_i0    , fr2_i0    , cp_ice_msh, jpl  )
-               !
-            IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   &
-               &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx )
-            !
-         CASE ( jp_cpl )
-            
-            CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )
-
+            ! albedo depends on cloud fraction because of non-linear spectral effects
+            zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
+            CALL blk_ice_core_flx( t_su, zalb_ice )
+            IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su )
+            IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
+         CASE ( jp_purecpl )
+            ! albedo depends on cloud fraction because of non-linear spectral effects
+            zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
+                                 CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su )
+            ! clem: evap_ice is forced to 0 in coupled mode for now 
+            !       but it needs to be changed (along with modif in limthd_dh) once heat flux from evap will be avail. from atm. models
+            evap_ice  (:,:,:) = 0._wp   ;   devap_ice (:,:,:) = 0._wp
+            IF( nn_limflx == 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
          END SELECT
-         
-         !------------------------------!
-         ! --- LIM-3 main time-step --- !
-         !------------------------------!
-         numit = numit + nn_fsbc                     ! Ice model time step
-         !                                                   
-         CALL sbc_lim_bef                   ! Store previous ice values
-
-         CALL sbc_lim_diag0                 ! set diag of mass, heat and salt fluxes to 0
-         
-         CALL lim_rst_opn( kt )             ! Open Ice restart file
-         !
-         ! ----------------------------------------------
-         ! ice dynamics and transport (except in 1D case)
-         ! ----------------------------------------------
-         IF( .NOT. lk_c1d ) THEN
-            
-            CALL lim_dyn( kt )              ! Ice dynamics    ( rheology/dynamics )
-            
-            CALL lim_trp( kt )              ! Ice transport   ( Advection/diffusion )
-            
-            IF( nn_monocat /= 2 ) CALL lim_itd_me  ! Mechanical redistribution ! (ridging/rafting)
-
-#if defined key_bdy
-            CALL bdy_ice_lim( kt )         ! bdy ice thermo 
-            IF( ln_icectl )   CALL lim_prt( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )
-#endif
-            CALL lim_update1( kt )
-            
-         ENDIF
-         
-         CALL sbc_lim_bef                  ! Store previous ice values
- 
-         ! ----------------------------------------------
-         ! ice thermodynamics
-         ! ----------------------------------------------
-         CALL lim_var_agg(1)
-         
-         ! previous lead fraction and ice volume for flux calculations
-         pfrld(:,:)   = 1._wp - at_i(:,:)
-         phicif(:,:)  = vt_i(:,:)
-         
-         SELECT CASE( kblk )
-         CASE ( jp_cpl )
-            CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su    )
-            IF( nn_limflx == 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   &
-               &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx )
-            ! Latent heat flux is forced to 0 in coupled: it is included in qns (non-solar heat flux)
-            qla_ice  (:,:,:) = 0._wp
-            dqla_ice (:,:,:) = 0._wp
-         END SELECT
-         !
-         CALL lim_thd( kt )                         ! Ice thermodynamics 
-         
+         CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )
+
+         !----------------------------!
+         ! --- ice thermodynamics --- !
+         !----------------------------!
+         CALL lim_thd( kt )                         ! Ice thermodynamics      
+         !
          CALL lim_update2( kt )                     ! Corrections
          !
@@ -237,7 +234,7 @@
          !
          IF(ln_limdiaout) CALL lim_diahsb           ! Diagnostics and outputs 
-         
+         !
          CALL lim_wri( 1 )                          ! Ice outputs 
-         
+         !
          IF( kt == nit000 .AND. ln_rstart )   &
             &             CALL iom_close( numrir )  ! close input ice restart file
@@ -247,17 +244,14 @@
          IF( ln_icectl )  CALL lim_ctl( kt )        ! alerts in case of model crash
          !
-         CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )
-         !
       ENDIF   ! End sea-ice time step only
 
-      !--------------------------------!
-      ! --- at all ocean time step --- !
-      !--------------------------------!
-      ! Update surface ocean stresses (only in ice-dynamic case)
-      !    otherwise the atm.-ocean stresses are used everywhere
+      !-------------------------!
+      ! --- Ocean time step --- !
+      !-------------------------!
+      ! Update surface ocean stresses (only in ice-dynamic case) otherwise the atm.-ocean stresses are used everywhere
       IF( ln_limdyn )     CALL lim_sbc_tau( kt, ub(:,:,1), vb(:,:,1) )  ! using before instantaneous surf. currents
 !!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!!
       !
-      IF( nn_timing == 1 )  CALL timing_stop('sbc_ice_lim')
+      IF( nn_timing == 1 ) CALL timing_stop('sbc_ice_lim')
       !
    END SUBROUTINE sbc_ice_lim
@@ -300,4 +294,6 @@
       !
       CALL lim_itd_init                ! ice thickness distribution initialization
+      !
+      CALL lim_hdf_init                ! set ice horizontal diffusion computation parameters
       !
       CALL lim_thd_init                ! set ice thermodynics parameters
@@ -346,5 +342,5 @@
       !!-------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namicerun/ jpl, nlay_i, nlay_s, cn_icerst_in, cn_icerst_out,   &
+      NAMELIST/namicerun/ jpl, nlay_i, nlay_s, cn_icerst_in, cn_icerst_indir, cn_icerst_out, cn_icerst_outdir,  &
          &                ln_limdyn, rn_amax, ln_limdiahsb, ln_limdiaout, ln_icectl, iiceprt, jiceprt  
       !!-------------------------------------------------------------------
@@ -475,6 +471,5 @@
 
    
-   SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice,   &
-         &                          pdqn_ice, pqla_ice, pdql_ice, k_limflx )
+   SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_limflx )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE ice_lim_flx  ***
@@ -494,6 +489,6 @@
       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqsr_ice   ! net solar flux
       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdqn_ice   ! non solar flux sensitivity
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqla_ice   ! latent heat flux
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdql_ice   ! latent heat flux sensitivity
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pevap_ice  ! sublimation
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdevap_ice ! sublimation sensitivity
       !
       INTEGER  ::   jl      ! dummy loop index
@@ -504,7 +499,7 @@
       REAL(wp), POINTER, DIMENSION(:,:) :: z_qsr_m   ! Mean solar heat flux over all categories
       REAL(wp), POINTER, DIMENSION(:,:) :: z_qns_m   ! Mean non solar heat flux over all categories
-      REAL(wp), POINTER, DIMENSION(:,:) :: z_qla_m   ! Mean latent heat flux over all categories
+      REAL(wp), POINTER, DIMENSION(:,:) :: z_evap_m  ! Mean sublimation over all categories
       REAL(wp), POINTER, DIMENSION(:,:) :: z_dqn_m   ! Mean d(qns)/dT over all categories
-      REAL(wp), POINTER, DIMENSION(:,:) :: z_dql_m   ! Mean d(qla)/dT over all categories
+      REAL(wp), POINTER, DIMENSION(:,:) :: z_devap_m ! Mean d(evap)/dT over all categories
       !!----------------------------------------------------------------------
 
@@ -514,23 +509,23 @@
       SELECT CASE( k_limflx )                              !==  averaged on all ice categories  ==!
       CASE( 0 , 1 )
-         CALL wrk_alloc( jpi,jpj, z_qsr_m, z_qns_m, z_qla_m, z_dqn_m, z_dql_m)
-         !
-         z_qns_m(:,:) = fice_ice_ave ( pqns_ice (:,:,:) )
-         z_qsr_m(:,:) = fice_ice_ave ( pqsr_ice (:,:,:) )
-         z_dqn_m(:,:) = fice_ice_ave ( pdqn_ice (:,:,:) )
-         z_qla_m(:,:) = fice_ice_ave ( pqla_ice (:,:,:) )
-         z_dql_m(:,:) = fice_ice_ave ( pdql_ice (:,:,:) )
+         CALL wrk_alloc( jpi,jpj, z_qsr_m, z_qns_m, z_evap_m, z_dqn_m, z_devap_m)
+         !
+         z_qns_m  (:,:) = fice_ice_ave ( pqns_ice (:,:,:) )
+         z_qsr_m  (:,:) = fice_ice_ave ( pqsr_ice (:,:,:) )
+         z_dqn_m  (:,:) = fice_ice_ave ( pdqn_ice (:,:,:) )
+         z_evap_m (:,:) = fice_ice_ave ( pevap_ice (:,:,:) )
+         z_devap_m(:,:) = fice_ice_ave ( pdevap_ice (:,:,:) )
          DO jl = 1, jpl
-            pdqn_ice(:,:,jl) = z_dqn_m(:,:)
-            pdql_ice(:,:,jl) = z_dql_m(:,:)
+            pdqn_ice  (:,:,jl) = z_dqn_m(:,:)
+            pdevap_ice(:,:,jl) = z_devap_m(:,:)
          END DO
          !
          DO jl = 1, jpl
-            pqns_ice(:,:,jl) = z_qns_m(:,:)
-            pqsr_ice(:,:,jl) = z_qsr_m(:,:)
-            pqla_ice(:,:,jl) = z_qla_m(:,:)
+            pqns_ice (:,:,jl) = z_qns_m(:,:)
+            pqsr_ice (:,:,jl) = z_qsr_m(:,:)
+            pevap_ice(:,:,jl) = z_evap_m(:,:)
          END DO
          !
-         CALL wrk_dealloc( jpi,jpj, z_qsr_m, z_qns_m, z_qla_m, z_dqn_m, z_dql_m)
+         CALL wrk_dealloc( jpi,jpj, z_qsr_m, z_qns_m, z_evap_m, z_dqn_m, z_devap_m)
       END SELECT
 
@@ -542,7 +537,7 @@
          ztem_m(:,:) = fice_ice_ave ( ptn_ice  (:,:,:) ) 
          DO jl = 1, jpl
-            pqns_ice(:,:,jl) = pqns_ice(:,:,jl) + pdqn_ice(:,:,jl) * (ptn_ice(:,:,jl) - ztem_m(:,:))
-            pqla_ice(:,:,jl) = pqla_ice(:,:,jl) + pdql_ice(:,:,jl) * (ptn_ice(:,:,jl) - ztem_m(:,:))
-            pqsr_ice(:,:,jl) = pqsr_ice(:,:,jl) * ( 1._wp - palb_ice(:,:,jl) ) / ( 1._wp - zalb_m(:,:) ) 
+            pqns_ice (:,:,jl) = pqns_ice (:,:,jl) + pdqn_ice  (:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
+            pevap_ice(:,:,jl) = pevap_ice(:,:,jl) + pdevap_ice(:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
+            pqsr_ice (:,:,jl) = pqsr_ice (:,:,jl) * ( 1._wp - palb_ice(:,:,jl) ) / ( 1._wp - zalb_m(:,:) ) 
          END DO
          !
@@ -593,5 +588,4 @@
       wfx_spr(:,:) = 0._wp   ;   
       
-      hfx_in (:,:) = 0._wp   ;   hfx_out(:,:) = 0._wp
       hfx_thd(:,:) = 0._wp   ;   
       hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp
@@ -610,5 +604,6 @@
       
    END SUBROUTINE sbc_lim_diag0
-      
+
+     
    FUNCTION fice_cell_ave ( ptab )
       !!--------------------------------------------------------------------------
@@ -620,5 +615,4 @@
       
       fice_cell_ave (:,:) = 0.0_wp
-      
       DO jl = 1, jpl
          fice_cell_ave (:,:) = fice_cell_ave (:,:) + a_i (:,:,jl) * ptab (:,:,jl)
