Index: NEMO/trunk/src/ICE/ice.F90
===================================================================
--- NEMO/trunk/src/ICE/ice.F90	(revision 10534)
+++ NEMO/trunk/src/ICE/ice.F90	(revision 10535)
@@ -150,5 +150,5 @@
    LOGICAL , PUBLIC ::   ln_adv_UMx       !: Ultimate-Macho advection scheme
    !
-   !                                     !!** ice-surface forcing namelist (namforcing) **
+   !                                     !!** ice-surface boundary conditions namelist (namsbc) **
                                           ! -- icethd_dh -- !
    REAL(wp), PUBLIC ::   rn_blow_s        !: coef. for partitioning of snowfall between leads and sea ice
Index: NEMO/trunk/src/ICE/icealb.F90
===================================================================
--- NEMO/trunk/src/ICE/icealb.F90	(revision 10534)
+++ NEMO/trunk/src/ICE/icealb.F90	(revision 10535)
@@ -27,5 +27,5 @@
 
    PUBLIC   ice_alb_init   ! called in icestp
-   PUBLIC   ice_alb        ! called in iceforcing.F90 and iceupdate.F90
+   PUBLIC   ice_alb        ! called in icesbc.F90 and iceupdate.F90
 
    REAL(wp), PUBLIC, PARAMETER ::   rn_alb_oce = 0.066   !: ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001)
Index: NEMO/trunk/src/ICE/iceforcing.F90
===================================================================
--- NEMO/trunk/src/ICE/iceforcing.F90	(revision 10534)
+++ 	(revision )
@@ -1,317 +1,0 @@
-MODULE iceforcing
-   !!======================================================================
-   !!                       ***  MODULE  iceforcing  ***
-   !! Sea-Ice :   air-ice forcing fields
-   !!=====================================================================
-   !! History :  4.0  !  2017-08  (C. Rousset)       Original code
-   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
-   !!----------------------------------------------------------------------
-#if defined key_si3
-   !!----------------------------------------------------------------------
-   !!   'key_si3' :                                     SI3 sea-ice model
-   !!----------------------------------------------------------------------
-   USE oce            ! ocean dynamics and tracers
-   USE dom_oce        ! ocean space and time domain
-   USE ice            ! sea-ice: variables
-   USE sbc_oce        ! Surface boundary condition: ocean fields
-   USE sbc_ice        ! Surface boundary condition: ice   fields
-   USE usrdef_sbc     ! Surface boundary condition: user defined
-   USE sbcblk         ! Surface boundary condition: bulk
-   USE sbccpl         ! Surface boundary condition: coupled interface
-   USE icealb         ! sea-ice: albedo
-   !
-   USE in_out_manager ! I/O manager
-   USE iom            ! I/O manager library
-   USE lib_mpp        ! MPP library
-   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
-   USE lbclnk         ! lateral boundary conditions (or mpp links)
-   USE timing         ! Timing
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC ice_forcing_tau   ! called by icestp.F90
-   PUBLIC ice_forcing_flx   ! called by icestp.F90
-   PUBLIC ice_forcing_init  ! called by icestp.F90
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
-   !! $Id$
-   !! Software governed by the CeCILL license (see ./LICENSE)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE ice_forcing_tau( kt, ksbc, utau_ice, vtau_ice )
-      !!-------------------------------------------------------------------
-      !!                  ***  ROUTINE ice_forcing_tau  ***
-      !!
-      !! ** Purpose : provide surface boundary condition for sea ice (momentum)
-      !!
-      !! ** Action  : It provides the following fields:
-      !!              utau_ice, vtau_ice : surface ice stress (U- & V-points) [N/m2]
-      !!-------------------------------------------------------------------
-      INTEGER                     , INTENT(in   ) ::   kt                   ! ocean time step
-      INTEGER                     , INTENT(in   ) ::   ksbc                 ! type of sbc flux
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   utau_ice, vtau_ice   ! air-ice stress   [N/m2]
-      !!
-      INTEGER  ::   ji, jj                 ! dummy loop index
-      REAL(wp), DIMENSION(jpi,jpj) ::   zutau_ice, zvtau_ice 
-      !!-------------------------------------------------------------------
-      !
-      IF( ln_timing )   CALL timing_start('ice_forcing')
-      !
-      IF( kt == nit000 .AND. lwp ) THEN
-         WRITE(numout,*)
-         WRITE(numout,*)'ice_forcing_tau: Surface boundary condition for sea ice (momentum)'
-         WRITE(numout,*)'~~~~~~~~~~~~~~~'
-      ENDIF
-      !
-      SELECT CASE( ksbc )
-         CASE( jp_usr     )   ;    CALL usrdef_sbc_ice_tau( kt )                 ! user defined formulation
-         CASE( jp_blk     )   ;    CALL blk_ice_tau                              ! 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 sbc_cpl_ice_tau( zutau_ice , zvtau_ice )
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               utau_ice(ji,jj) = utau_ice(ji,jj) * xcplmask(ji,jj,0) + zutau_ice(ji,jj) * ( 1. - xcplmask(ji,jj,0) )
-               vtau_ice(ji,jj) = vtau_ice(ji,jj) * xcplmask(ji,jj,0) + zvtau_ice(ji,jj) * ( 1. - xcplmask(ji,jj,0) )
-            END DO
-         END DO
-         CALL lbc_lnk_multi( 'iceforcing', utau_ice, 'U', -1., vtau_ice, 'V', -1. )
-      ENDIF
-      !
-      IF( ln_timing )   CALL timing_stop('ice_forcing')
-      !
-   END SUBROUTINE ice_forcing_tau
-
-   
-   SUBROUTINE ice_forcing_flx( kt, ksbc )
-      !!-------------------------------------------------------------------
-      !!                  ***  ROUTINE ice_forcing_flx  ***
-      !!
-      !! ** Purpose : provide surface boundary condition for sea ice (flux)
-      !!
-      !! ** Action  : It provides the following fields used in sea ice model:
-      !!                emp_oce , emp_ice                        = E-P over ocean and sea ice                    [Kg/m2/s]
-      !!                sprecip                                  = solid precipitation                           [Kg/m2/s]
-      !!                evap_ice                                 = sublimation                                   [Kg/m2/s]
-      !!                qsr_tot , qns_tot                        = solar & non solar heat flux (total)           [W/m2]
-      !!                qsr_ice , qns_ice                        = solar & non solar heat flux over ice          [W/m2]
-      !!                dqns_ice                                 = non solar  heat sensistivity                  [W/m2]
-      !!                qemp_oce, qemp_ice, qprec_ice, qevap_ice = sensible heat (associated with evap & precip) [W/m2]
-      !!            + some fields that are not used outside this module: 
-      !!                qla_ice                                  = latent heat flux over ice                     [W/m2]
-      !!                dqla_ice                                 = latent heat sensistivity                      [W/m2]
-      !!                tprecip                                  = total  precipitation                          [Kg/m2/s]
-      !!                alb_ice                                  = albedo above sea ice
-      !!-------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt     ! ocean time step
-      INTEGER, INTENT(in) ::   ksbc   ! flux formulation (user defined, bulk or Pure Coupled)
-      !
-      INTEGER  ::   ji, jj, jl                                ! dummy loop index
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zalb              ! 2D workspace
-      !!--------------------------------------------------------------------
-      !
-      IF( ln_timing )   CALL timing_start('ice_forcing_flx')
-
-      IF( kt == nit000 .AND. lwp ) THEN
-         WRITE(numout,*)
-         WRITE(numout,*)'ice_forcing_flx: Surface boundary condition for sea ice (flux)'
-         WRITE(numout,*)'~~~~~~~~~~~~~~~'
-      ENDIF
-
-      ! --- cloud-sky and overcast-sky ice albedos --- !
-      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os )
-
-      ! albedo depends on cloud fraction because of non-linear spectral effects
-!!gm cldf_ice is a real, DOCTOR naming rule: start with cd means CHARACTER passed in argument !
-      alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
-      !
-      !
-      SELECT CASE( ksbc )   !== fluxes over sea ice ==!
-      !
-      CASE( jp_usr )              !--- user defined formulation
-                                  CALL usrdef_sbc_ice_flx( kt, h_s, h_i )
-      CASE( jp_blk )              !--- bulk formulation
-                                  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_flxdist /= -1 )   CALL ice_flx_dist   ( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist )
-         !                        !    compute conduction flux and surface temperature (as in Jules surface module)
-         IF( ln_cndflx .AND. .NOT.ln_cndemulate ) &
-            &                     CALL blk_ice_qcn    ( ln_virtual_itd, t_su, t_bo, h_s, h_i )
-      CASE ( jp_purecpl )         !--- coupled formulation
-                                  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_flxdist /= -1 )   CALL ice_flx_dist   ( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist )
-      END SELECT
-
-      !--- output ice albedo and surface albedo ---!
-      IF( iom_use('icealb') ) THEN
-         WHERE( at_i_b <= epsi06 )   ;   zalb(:,:) = rn_alb_oce
-         ELSEWHERE                   ;   zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) / at_i_b
-         END WHERE
-         CALL iom_put( "icealb" , zalb(:,:) )
-      ENDIF
-      IF( iom_use('albedo') ) THEN
-         zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) + rn_alb_oce * ( 1._wp - at_i_b )
-         CALL iom_put( "albedo" , zalb(:,:) )
-      ENDIF
-      !
-      IF( ln_timing )   CALL timing_stop('ice_forcing_flx')
-      !
-   END SUBROUTINE ice_forcing_flx
-
-
-   SUBROUTINE ice_flx_dist( ptn_ice, palb_ice, pqns_ice, pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_flxdist )
-      !!-------------------------------------------------------------------
-      !!                  ***  ROUTINE ice_flx_dist  ***
-      !!
-      !! ** Purpose :   update the ice surface boundary condition by averaging 
-      !!              and/or redistributing fluxes on ice categories
-      !!
-      !! ** Method  :   average then redistribute
-      !!
-      !! ** Action  :   depends on k_flxdist
-      !!                = -1  Do nothing (needs N(cat) fluxes)
-      !!                =  0  Average N(cat) fluxes then apply the average over the N(cat) ice 
-      !!                =  1  Average N(cat) fluxes then redistribute over the N(cat) ice
-      !!                                                 using T-ice and albedo sensitivity
-      !!                =  2  Redistribute a single flux over categories
-      !!-------------------------------------------------------------------
-      INTEGER                   , INTENT(in   ) ::   k_flxdist  ! redistributor
-      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptn_ice    ! ice surface temperature
-      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   palb_ice   ! ice albedo
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqns_ice   ! non solar flux
-      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) ::   pevap_ice  ! sublimation
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdevap_ice ! sublimation sensitivity
-      !
-      INTEGER  ::   jl      ! dummy loop index
-      !
-      REAL(wp), DIMENSION(jpi,jpj) ::   z1_at_i   ! inverse of concentration
-      !
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_qsr_m   ! Mean solar heat flux over all categories
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_qns_m   ! Mean non solar heat flux over all categories
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_evap_m  ! Mean sublimation over all categories
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_dqn_m   ! Mean d(qns)/dT over all categories
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_devap_m ! Mean d(evap)/dT over all categories
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zalb_m    ! Mean albedo over all categories
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   ztem_m    ! Mean temperature over all categories
-      !!----------------------------------------------------------------------
-      !
-      WHERE ( at_i (:,:) > 0._wp )   ; z1_at_i(:,:) = 1._wp / at_i (:,:)
-      ELSEWHERE                      ; z1_at_i(:,:) = 0._wp
-      END WHERE
-      
-      SELECT CASE( k_flxdist )       !==  averaged on all ice categories  ==!
-      !
-      CASE( 0 , 1 )
-         !
-         ALLOCATE( z_qns_m(jpi,jpj), z_qsr_m(jpi,jpj), z_dqn_m(jpi,jpj), z_evap_m(jpi,jpj), z_devap_m(jpi,jpj) )  
-         !
-         z_qns_m  (:,:) = SUM( a_i(:,:,:) * pqns_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:)
-         z_qsr_m  (:,:) = SUM( a_i(:,:,:) * pqsr_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:)
-         z_dqn_m  (:,:) = SUM( a_i(:,:,:) * pdqn_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:)
-         z_evap_m (:,:) = SUM( a_i(:,:,:) * pevap_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
-         z_devap_m(:,:) = SUM( a_i(:,:,:) * pdevap_ice(:,:,:) , dim=3 ) * z1_at_i(:,:)
-         DO jl = 1, jpl
-            pqns_ice  (:,:,jl) = z_qns_m (:,:)
-            pqsr_ice  (:,:,jl) = z_qsr_m (:,:)
-            pdqn_ice  (:,:,jl) = z_dqn_m  (:,:)
-            pevap_ice (:,:,jl) = z_evap_m(:,:)
-            pdevap_ice(:,:,jl) = z_devap_m(:,:)
-         END DO
-         !
-         DEALLOCATE( z_qns_m, z_qsr_m, z_dqn_m, z_evap_m, z_devap_m )  
-         !
-      END SELECT
-      !
-      SELECT CASE( k_flxdist )       !==  redistribution on all ice categories  ==!
-      !
-      CASE( 1 , 2 )
-         !
-         ALLOCATE( zalb_m(jpi,jpj), ztem_m(jpi,jpj) )  
-         !
-         zalb_m(:,:) = SUM( a_i(:,:,:) * palb_ice(:,:,:) , dim=3 ) * z1_at_i(:,:)
-         ztem_m(:,:) = SUM( a_i(:,:,:) * ptn_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
-         DO jl = 1, jpl
-            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
-         !
-         DEALLOCATE( zalb_m, ztem_m )  
-         !
-      END SELECT
-      !
-   END SUBROUTINE ice_flx_dist
-
-
-   SUBROUTINE ice_forcing_init
-      !!-------------------------------------------------------------------
-      !!                  ***  ROUTINE ice_forcing_init  ***
-      !!
-      !! ** Purpose :   Physical constants and parameters linked to the ice dynamics
-      !!      
-      !! ** Method  :   Read the namforcing namelist and check the ice-dynamic
-      !!              parameter values called at the first timestep (nit000)
-      !!
-      !! ** input   :   Namelist namforcing
-      !!-------------------------------------------------------------------
-      INTEGER ::   ios, ioptio   ! Local integer
-      !!
-      NAMELIST/namforcing/ rn_cio, rn_blow_s, nn_flxdist, ln_cndflx, ln_cndemulate
-      !!-------------------------------------------------------------------
-      !
-      REWIND( numnam_ice_ref )         ! Namelist namforcing in reference namelist : Ice dynamics
-      READ  ( numnam_ice_ref, namforcing, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namforcing in reference namelist', lwp )
-      REWIND( numnam_ice_cfg )         ! Namelist namforcing in configuration namelist : Ice dynamics
-      READ  ( numnam_ice_cfg, namforcing, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namforcing in configuration namelist', lwp )
-      IF(lwm) WRITE( numoni, namforcing )
-      !
-      IF(lwp) THEN                     ! control print
-         WRITE(numout,*)
-         WRITE(numout,*) 'ice_forcing_init: ice parameters for ice dynamics '
-         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_flxdist    = ', nn_flxdist
-         WRITE(numout,*) '      Use conduction flux as surface condition         ln_cndflx     = ', ln_cndflx
-         WRITE(numout,*) '         emulate conduction flux                       ln_cndemulate = ', ln_cndemulate
-      ENDIF
-      !
-      IF(lwp) WRITE(numout,*)
-      SELECT CASE( nn_flxdist )         ! SI3 Multi-category heat flux formulation
-      CASE( -1  )
-         IF(lwp) WRITE(numout,*) '   SI3: use per-category fluxes (nn_flxdist = -1) '
-      CASE(  0  )
-         IF(lwp) WRITE(numout,*) '   SI3: use average per-category fluxes (nn_flxdist = 0) '
-      CASE(  1  )
-         IF(lwp) WRITE(numout,*) '   SI3: use average then redistribute per-category fluxes (nn_flxdist = 1) '
-         IF( ln_cpl )         CALL ctl_stop( 'ice_thd_init: the chosen nn_flxdist for SI3 in coupled mode must be /=1' )
-      CASE(  2  )
-         IF(lwp) WRITE(numout,*) '   SI3: Redistribute a single flux over categories (nn_flxdist = 2) '
-         IF( .NOT. ln_cpl )   CALL ctl_stop( 'ice_thd_init: the chosen nn_flxdist for SI3 in forced mode must be /=2' )
-      CASE DEFAULT
-         CALL ctl_stop( 'ice_thd_init: SI3 option, nn_flxdist, should be between -1 and 2' )
-      END SELECT
-      !
-   END SUBROUTINE ice_forcing_init
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option :         Empty module         NO SI3 sea-ice model
-   !!----------------------------------------------------------------------
-#endif
-
-   !!======================================================================
-END MODULE iceforcing
Index: NEMO/trunk/src/ICE/icesbc.F90
===================================================================
--- NEMO/trunk/src/ICE/icesbc.F90	(revision 10535)
+++ NEMO/trunk/src/ICE/icesbc.F90	(revision 10535)
@@ -0,0 +1,317 @@
+MODULE icesbc
+   !!======================================================================
+   !!                       ***  MODULE  icesbc  ***
+   !! Sea-Ice :   air-ice sbc fields
+   !!=====================================================================
+   !! History :  4.0  !  2017-08  (C. Rousset)       Original code
+   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
+   !!----------------------------------------------------------------------
+#if defined key_si3
+   !!----------------------------------------------------------------------
+   !!   'key_si3' :                                     SI3 sea-ice model
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers
+   USE dom_oce        ! ocean space and time domain
+   USE ice            ! sea-ice: variables
+   USE sbc_oce        ! Surface boundary condition: ocean fields
+   USE sbc_ice        ! Surface boundary condition: ice   fields
+   USE usrdef_sbc     ! Surface boundary condition: user defined
+   USE sbcblk         ! Surface boundary condition: bulk
+   USE sbccpl         ! Surface boundary condition: coupled interface
+   USE icealb         ! sea-ice: albedo
+   !
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O manager library
+   USE lib_mpp        ! MPP library
+   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
+   USE lbclnk         ! lateral boundary conditions (or mpp links)
+   USE timing         ! Timing
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC ice_sbc_tau   ! called by icestp.F90
+   PUBLIC ice_sbc_flx   ! called by icestp.F90
+   PUBLIC ice_sbc_init  ! called by icestp.F90
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE ice_sbc_tau( kt, ksbc, utau_ice, vtau_ice )
+      !!-------------------------------------------------------------------
+      !!                  ***  ROUTINE ice_sbc_tau  ***
+      !!
+      !! ** Purpose : provide surface boundary condition for sea ice (momentum)
+      !!
+      !! ** Action  : It provides the following fields:
+      !!              utau_ice, vtau_ice : surface ice stress (U- & V-points) [N/m2]
+      !!-------------------------------------------------------------------
+      INTEGER                     , INTENT(in   ) ::   kt                   ! ocean time step
+      INTEGER                     , INTENT(in   ) ::   ksbc                 ! type of sbc flux
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   utau_ice, vtau_ice   ! air-ice stress   [N/m2]
+      !!
+      INTEGER  ::   ji, jj                 ! dummy loop index
+      REAL(wp), DIMENSION(jpi,jpj) ::   zutau_ice, zvtau_ice 
+      !!-------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('ice_sbc')
+      !
+      IF( kt == nit000 .AND. lwp ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*)'ice_sbc_tau: Surface boundary condition for sea ice (momentum)'
+         WRITE(numout,*)'~~~~~~~~~~~~~~~'
+      ENDIF
+      !
+      SELECT CASE( ksbc )
+         CASE( jp_usr     )   ;    CALL usrdef_sbc_ice_tau( kt )                 ! user defined formulation
+         CASE( jp_blk     )   ;    CALL blk_ice_tau                              ! 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 sbc_cpl_ice_tau( zutau_ice , zvtau_ice )
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               utau_ice(ji,jj) = utau_ice(ji,jj) * xcplmask(ji,jj,0) + zutau_ice(ji,jj) * ( 1. - xcplmask(ji,jj,0) )
+               vtau_ice(ji,jj) = vtau_ice(ji,jj) * xcplmask(ji,jj,0) + zvtau_ice(ji,jj) * ( 1. - xcplmask(ji,jj,0) )
+            END DO
+         END DO
+         CALL lbc_lnk_multi( 'icesbc', utau_ice, 'U', -1., vtau_ice, 'V', -1. )
+      ENDIF
+      !
+      IF( ln_timing )   CALL timing_stop('ice_sbc')
+      !
+   END SUBROUTINE ice_sbc_tau
+
+   
+   SUBROUTINE ice_sbc_flx( kt, ksbc )
+      !!-------------------------------------------------------------------
+      !!                  ***  ROUTINE ice_sbc_flx  ***
+      !!
+      !! ** Purpose : provide surface boundary condition for sea ice (flux)
+      !!
+      !! ** Action  : It provides the following fields used in sea ice model:
+      !!                emp_oce , emp_ice                        = E-P over ocean and sea ice                    [Kg/m2/s]
+      !!                sprecip                                  = solid precipitation                           [Kg/m2/s]
+      !!                evap_ice                                 = sublimation                                   [Kg/m2/s]
+      !!                qsr_tot , qns_tot                        = solar & non solar heat flux (total)           [W/m2]
+      !!                qsr_ice , qns_ice                        = solar & non solar heat flux over ice          [W/m2]
+      !!                dqns_ice                                 = non solar  heat sensistivity                  [W/m2]
+      !!                qemp_oce, qemp_ice, qprec_ice, qevap_ice = sensible heat (associated with evap & precip) [W/m2]
+      !!            + some fields that are not used outside this module: 
+      !!                qla_ice                                  = latent heat flux over ice                     [W/m2]
+      !!                dqla_ice                                 = latent heat sensistivity                      [W/m2]
+      !!                tprecip                                  = total  precipitation                          [Kg/m2/s]
+      !!                alb_ice                                  = albedo above sea ice
+      !!-------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt     ! ocean time step
+      INTEGER, INTENT(in) ::   ksbc   ! flux formulation (user defined, bulk or Pure Coupled)
+      !
+      INTEGER  ::   ji, jj, jl                                ! dummy loop index
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zalb              ! 2D workspace
+      !!--------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('ice_sbc_flx')
+
+      IF( kt == nit000 .AND. lwp ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*)'ice_sbc_flx: Surface boundary condition for sea ice (flux)'
+         WRITE(numout,*)'~~~~~~~~~~~~~~~'
+      ENDIF
+
+      ! --- cloud-sky and overcast-sky ice albedos --- !
+      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os )
+
+      ! albedo depends on cloud fraction because of non-linear spectral effects
+!!gm cldf_ice is a real, DOCTOR naming rule: start with cd means CHARACTER passed in argument !
+      alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
+      !
+      !
+      SELECT CASE( ksbc )   !== fluxes over sea ice ==!
+      !
+      CASE( jp_usr )              !--- user defined formulation
+                                  CALL usrdef_sbc_ice_flx( kt, h_s, h_i )
+      CASE( jp_blk )              !--- bulk formulation
+                                  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_flxdist /= -1 )   CALL ice_flx_dist   ( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist )
+         !                        !    compute conduction flux and surface temperature (as in Jules surface module)
+         IF( ln_cndflx .AND. .NOT.ln_cndemulate ) &
+            &                     CALL blk_ice_qcn    ( ln_virtual_itd, t_su, t_bo, h_s, h_i )
+      CASE ( jp_purecpl )         !--- coupled formulation
+                                  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_flxdist /= -1 )   CALL ice_flx_dist   ( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist )
+      END SELECT
+
+      !--- output ice albedo and surface albedo ---!
+      IF( iom_use('icealb') ) THEN
+         WHERE( at_i_b <= epsi06 )   ;   zalb(:,:) = rn_alb_oce
+         ELSEWHERE                   ;   zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) / at_i_b
+         END WHERE
+         CALL iom_put( "icealb" , zalb(:,:) )
+      ENDIF
+      IF( iom_use('albedo') ) THEN
+         zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) + rn_alb_oce * ( 1._wp - at_i_b )
+         CALL iom_put( "albedo" , zalb(:,:) )
+      ENDIF
+      !
+      IF( ln_timing )   CALL timing_stop('ice_sbc_flx')
+      !
+   END SUBROUTINE ice_sbc_flx
+
+
+   SUBROUTINE ice_flx_dist( ptn_ice, palb_ice, pqns_ice, pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_flxdist )
+      !!-------------------------------------------------------------------
+      !!                  ***  ROUTINE ice_flx_dist  ***
+      !!
+      !! ** Purpose :   update the ice surface boundary condition by averaging 
+      !!              and/or redistributing fluxes on ice categories
+      !!
+      !! ** Method  :   average then redistribute
+      !!
+      !! ** Action  :   depends on k_flxdist
+      !!                = -1  Do nothing (needs N(cat) fluxes)
+      !!                =  0  Average N(cat) fluxes then apply the average over the N(cat) ice 
+      !!                =  1  Average N(cat) fluxes then redistribute over the N(cat) ice
+      !!                                                 using T-ice and albedo sensitivity
+      !!                =  2  Redistribute a single flux over categories
+      !!-------------------------------------------------------------------
+      INTEGER                   , INTENT(in   ) ::   k_flxdist  ! redistributor
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptn_ice    ! ice surface temperature
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   palb_ice   ! ice albedo
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqns_ice   ! non solar flux
+      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) ::   pevap_ice  ! sublimation
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdevap_ice ! sublimation sensitivity
+      !
+      INTEGER  ::   jl      ! dummy loop index
+      !
+      REAL(wp), DIMENSION(jpi,jpj) ::   z1_at_i   ! inverse of concentration
+      !
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_qsr_m   ! Mean solar heat flux over all categories
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_qns_m   ! Mean non solar heat flux over all categories
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_evap_m  ! Mean sublimation over all categories
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_dqn_m   ! Mean d(qns)/dT over all categories
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_devap_m ! Mean d(evap)/dT over all categories
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zalb_m    ! Mean albedo over all categories
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   ztem_m    ! Mean temperature over all categories
+      !!----------------------------------------------------------------------
+      !
+      WHERE ( at_i (:,:) > 0._wp )   ; z1_at_i(:,:) = 1._wp / at_i (:,:)
+      ELSEWHERE                      ; z1_at_i(:,:) = 0._wp
+      END WHERE
+      
+      SELECT CASE( k_flxdist )       !==  averaged on all ice categories  ==!
+      !
+      CASE( 0 , 1 )
+         !
+         ALLOCATE( z_qns_m(jpi,jpj), z_qsr_m(jpi,jpj), z_dqn_m(jpi,jpj), z_evap_m(jpi,jpj), z_devap_m(jpi,jpj) )  
+         !
+         z_qns_m  (:,:) = SUM( a_i(:,:,:) * pqns_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:)
+         z_qsr_m  (:,:) = SUM( a_i(:,:,:) * pqsr_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:)
+         z_dqn_m  (:,:) = SUM( a_i(:,:,:) * pdqn_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:)
+         z_evap_m (:,:) = SUM( a_i(:,:,:) * pevap_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
+         z_devap_m(:,:) = SUM( a_i(:,:,:) * pdevap_ice(:,:,:) , dim=3 ) * z1_at_i(:,:)
+         DO jl = 1, jpl
+            pqns_ice  (:,:,jl) = z_qns_m (:,:)
+            pqsr_ice  (:,:,jl) = z_qsr_m (:,:)
+            pdqn_ice  (:,:,jl) = z_dqn_m  (:,:)
+            pevap_ice (:,:,jl) = z_evap_m(:,:)
+            pdevap_ice(:,:,jl) = z_devap_m(:,:)
+         END DO
+         !
+         DEALLOCATE( z_qns_m, z_qsr_m, z_dqn_m, z_evap_m, z_devap_m )  
+         !
+      END SELECT
+      !
+      SELECT CASE( k_flxdist )       !==  redistribution on all ice categories  ==!
+      !
+      CASE( 1 , 2 )
+         !
+         ALLOCATE( zalb_m(jpi,jpj), ztem_m(jpi,jpj) )  
+         !
+         zalb_m(:,:) = SUM( a_i(:,:,:) * palb_ice(:,:,:) , dim=3 ) * z1_at_i(:,:)
+         ztem_m(:,:) = SUM( a_i(:,:,:) * ptn_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
+         DO jl = 1, jpl
+            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
+         !
+         DEALLOCATE( zalb_m, ztem_m )  
+         !
+      END SELECT
+      !
+   END SUBROUTINE ice_flx_dist
+
+
+   SUBROUTINE ice_sbc_init
+      !!-------------------------------------------------------------------
+      !!                  ***  ROUTINE ice_sbc_init  ***
+      !!
+      !! ** Purpose :   Physical constants and parameters linked to the ice dynamics
+      !!      
+      !! ** Method  :   Read the namsbc namelist and check the ice-dynamic
+      !!              parameter values called at the first timestep (nit000)
+      !!
+      !! ** input   :   Namelist namsbc
+      !!-------------------------------------------------------------------
+      INTEGER ::   ios, ioptio   ! Local integer
+      !!
+      NAMELIST/namsbc/ rn_cio, rn_blow_s, nn_flxdist, ln_cndflx, ln_cndemulate
+      !!-------------------------------------------------------------------
+      !
+      REWIND( numnam_ice_ref )         ! Namelist namsbc in reference namelist : Ice dynamics
+      READ  ( numnam_ice_ref, namsbc, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc in reference namelist', lwp )
+      REWIND( numnam_ice_cfg )         ! Namelist namsbc in configuration namelist : Ice dynamics
+      READ  ( numnam_ice_cfg, namsbc, IOSTAT = ios, ERR = 902 )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc in configuration namelist', lwp )
+      IF(lwm) WRITE( numoni, namsbc )
+      !
+      IF(lwp) THEN                     ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'ice_sbc_init: ice parameters for ice dynamics '
+         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namsbc:'
+         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_flxdist    = ', nn_flxdist
+         WRITE(numout,*) '      Use conduction flux as surface condition         ln_cndflx     = ', ln_cndflx
+         WRITE(numout,*) '         emulate conduction flux                       ln_cndemulate = ', ln_cndemulate
+      ENDIF
+      !
+      IF(lwp) WRITE(numout,*)
+      SELECT CASE( nn_flxdist )         ! SI3 Multi-category heat flux formulation
+      CASE( -1  )
+         IF(lwp) WRITE(numout,*) '   SI3: use per-category fluxes (nn_flxdist = -1) '
+      CASE(  0  )
+         IF(lwp) WRITE(numout,*) '   SI3: use average per-category fluxes (nn_flxdist = 0) '
+      CASE(  1  )
+         IF(lwp) WRITE(numout,*) '   SI3: use average then redistribute per-category fluxes (nn_flxdist = 1) '
+         IF( ln_cpl )         CALL ctl_stop( 'ice_thd_init: the chosen nn_flxdist for SI3 in coupled mode must be /=1' )
+      CASE(  2  )
+         IF(lwp) WRITE(numout,*) '   SI3: Redistribute a single flux over categories (nn_flxdist = 2) '
+         IF( .NOT. ln_cpl )   CALL ctl_stop( 'ice_thd_init: the chosen nn_flxdist for SI3 in forced mode must be /=2' )
+      CASE DEFAULT
+         CALL ctl_stop( 'ice_thd_init: SI3 option, nn_flxdist, should be between -1 and 2' )
+      END SELECT
+      !
+   END SUBROUTINE ice_sbc_init
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option :         Empty module         NO SI3 sea-ice model
+   !!----------------------------------------------------------------------
+#endif
+
+   !!======================================================================
+END MODULE icesbc
Index: NEMO/trunk/src/ICE/icestp.F90
===================================================================
--- NEMO/trunk/src/ICE/icestp.F90	(revision 10534)
+++ NEMO/trunk/src/ICE/icestp.F90	(revision 10535)
@@ -52,5 +52,5 @@
    USE sbc_ice        ! Surface boundary condition: ice   fields
    !
-   USE iceforcing     ! sea-ice: Surface boundary condition       !!gm why not icesbc module name
+   USE icesbc         ! sea-ice: Surface boundary conditions
    USE icedyn         ! sea-ice: dynamics
    USE icethd         ! sea-ice: thermodynamics
@@ -152,5 +152,5 @@
          !    utau_ice, vtau_ice = surface ice stress [N/m2]
          !------------------------------------------------!
-                                        CALL ice_forcing_tau( kt, ksbc, utau_ice, vtau_ice )          
+                                        CALL ice_sbc_tau( kt, ksbc, utau_ice, vtau_ice )          
          !-------------------------------------!
          ! --- ice dynamics and advection  --- !
@@ -183,5 +183,5 @@
          !    qprec_ice, qevap_ice
          !------------------------------------------------------!
-                                        CALL ice_forcing_flx( kt, ksbc )
+                                        CALL ice_sbc_flx( kt, ksbc )
          !----------------------------!
          ! --- ice thermodynamics --- !
@@ -241,5 +241,5 @@
       !                                ! Allocate the ice arrays (sbc_ice already allocated in sbc_init)
       ierr =        ice_alloc        ()      ! ice variables
-      ierr = ierr + sbc_ice_alloc    ()      ! surface forcing 
+      ierr = ierr + sbc_ice_alloc    ()      ! surface boundary conditions 
       ierr = ierr + ice1D_alloc      ()      ! thermodynamics
       !
@@ -261,5 +261,5 @@
       CALL ice_var_agg(1)
       !
-      CALL ice_forcing_init            ! set ice-ocean and ice-atm. coupling parameters
+      CALL ice_sbc_init                ! set ice-ocean and ice-atm. coupling parameters
       !
       CALL ice_dyn_init                ! set ice dynamics parameters
