Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 3609)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 3625)
@@ -12,4 +12,5 @@
    !!             -   ! 2010-10  (J. Chanut, C. Bricaud, G. Madec)  add the surface pressure forcing
    !!            3.4  ! 2011-11  (C. Harris) CICE added as an option
+   !!            3.5  ! 2012-11  (A. Coward, G. Madec) Rethink of heat, mass and salt surface fluxes
    !!----------------------------------------------------------------------
 
@@ -84,6 +85,7 @@
       INTEGER ::   icpt   ! local integer
       !!
-      NAMELIST/namsbc/ nn_fsbc   , ln_ana , ln_flx  , ln_blk_clio, ln_blk_core, ln_cpl,   &
-         &             ln_blk_mfs, ln_apr_dyn, nn_ice , ln_dm2dc, ln_rnf, ln_ssr     , nn_fwb, ln_cdgw
+      NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx,  ln_blk_clio, ln_blk_core, ln_cpl,   &
+         &             ln_blk_mfs, ln_apr_dyn, nn_ice,  nn_ice_embd, ln_dm2dc   , ln_rnf,   &
+         &             ln_ssr    , nn_fwb    , ln_cdgw
       !!----------------------------------------------------------------------
 
@@ -121,4 +123,5 @@
          WRITE(numout,*) '              Patm gradient added in ocean & ice Eqs.    ln_apr_dyn  = ', ln_apr_dyn
          WRITE(numout,*) '              ice management in the sbc (=0/1/2/3)       nn_ice      = ', nn_ice 
+         WRITE(numout,*) '              ice-ocean embedded/levitating (=0/1/2)     nn_ice_embd = ', nn_ice_embd
          WRITE(numout,*) '              daily mean to diurnal cycle qsr            ln_dm2dc    = ', ln_dm2dc 
          WRITE(numout,*) '              runoff / runoff mouths                     ln_rnf      = ', ln_rnf
@@ -136,9 +139,12 @@
          IF( sbc_rnf_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_rnf arrays' )
          nkrnf         = 0
-         rnf     (:,:) = 0.e0
-         rnfmsk  (:,:) = 0.e0
-         rnfmsk_z(:)   = 0.e0
+         rnf     (:,:) = 0.0_wp
+         rnfmsk  (:,:) = 0.0_wp
+         rnfmsk_z(:)   = 0.0_wp
       ENDIF
       IF( nn_ice == 0  )   fr_i(:,:) = 0.e0        ! no ice in the domain, ice fraction is always zero
+
+      sfx(:,:) = 0.0_wp                            ! the salt flux due to freezing/melting will be computed (i.e. will be non-zero) 
+                                                   ! only if sea-ice is present
 
       !                                            ! restartability   
@@ -157,6 +163,8 @@
       IF( nn_ice == 4 .AND. .NOT.( ln_blk_core .OR. lk_cpl ) )   &
          &   CALL ctl_stop( 'CICE sea-ice model requires ln_blk_core or lk_cpl' )
-      IF( nn_ice == 4 .AND. ( .NOT. ( cp_cfg == 'orca' ) .OR. lk_agrif ) )   &
-         &   CALL ctl_stop( 'CICE sea-ice model currently only available in a global ORCA configuration without AGRIF' )
+      IF( nn_ice == 4 .AND. lk_agrif )   &
+         &   CALL ctl_stop( 'CICE sea-ice model not currently available with AGRIF' )
+      IF( ( nn_ice == 3 .OR. nn_ice == 4 ) .AND. nn_ice_embd == 0 )   &
+         &   CALL ctl_stop( 'LIM3 and CICE sea-ice models require nn_ice_embd = 2 or 3' )
       
       IF( ln_dm2dc )   nday_qsr = -1   ! initialisation flag
@@ -226,6 +234,6 @@
       !! ** Action  : - set the ocean surface boundary condition at before and now 
       !!                time step, i.e.  
-      !!                utau_b, vtau_b, qns_b, qsr_b, emp_n, emps_b, qrp_b, erp_b
-      !!                utau  , vtau  , qns  , qsr  , emp  , emps  , qrp  , erp
+      !!                utau_b, vtau_b, qns_b, qsr_b, emp_n, sfx_b, qrp_b, erp_b
+      !!                utau  , vtau  , qns  , qsr  , emp  , sfx  , qrp  , erp
       !!              - updte the ice fraction : fr_i
       !!----------------------------------------------------------------------
@@ -243,6 +251,6 @@
          ! The 3D heat content due to qsr forcing is treated in traqsr
          ! qsr_b (:,:) = qsr (:,:)
-         emp_b (:,:) = emp (:,:)
-         emps_b(:,:) = emps(:,:)
+         emp_b(:,:) = emp(:,:)
+         sfx_b(:,:) = sfx(:,:)
       ENDIF
       !                                            ! ---------------------------------------- !
@@ -262,5 +270,5 @@
                                                             
       SELECT CASE( nsbc )                                ! Compute ocean surface boundary condition
-      !                                                  ! (i.e. utau,vtau, qns, qsr, emp, emps)
+      !                                                  ! (i.e. utau,vtau, qns, qsr, emp, sfx)
       CASE(  0 )   ;   CALL sbc_gyre    ( kt )                    ! analytical formulation : GYRE configuration
       CASE(  1 )   ;   CALL sbc_ana     ( kt )                    ! analytical formulation : uniform sbc
@@ -314,7 +322,12 @@
             CALL iom_get( numror, jpdom_autoglo, 'qns_b' , qns_b  )   ! before non solar heat flux (T-point)
             ! The 3D heat content due to qsr forcing is treated in traqsr
-            ! CALL iom_get( numror, jpdom_autoglo, 'qsr_b' , qsr_b  )   ! before     solar heat flux (T-point)
-            CALL iom_get( numror, jpdom_autoglo, 'emp_b' , emp_b  )   ! before     freshwater flux (T-point)
-            CALL iom_get( numror, jpdom_autoglo, 'emps_b', emps_b )   ! before C/D freshwater flux (T-point)
+            ! CALL iom_get( numror, jpdom_autoglo, 'qsr_b' , qsr_b  ) ! before     solar heat flux (T-point)
+            CALL iom_get( numror, jpdom_autoglo, 'emp_b', emp_b  )    ! before     freshwater flux (T-point)
+            ! To ensure restart capability with 3.3x/3.4 restart files    !! to be removed in v3.6
+            IF( iom_varid( numror, 'sfx_b', ldstop = .FALSE. ) > 0 ) THEN
+               CALL iom_get( numror, jpdom_autoglo, 'sfx_b', sfx_b )  ! before salt flux (T-point)
+            ELSE
+               sfx_b (:,:) = sfx(:,:)
+            ENDIF
          ELSE                                                   !* no restart: set from nit000 values
             IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields set to nit000'
@@ -322,7 +335,6 @@
             vtau_b(:,:) = vtau(:,:)
             qns_b (:,:) = qns (:,:)
-            ! qsr_b (:,:) = qsr (:,:)
-            emp_b (:,:) = emp (:,:)
-            emps_b(:,:) = emps(:,:)
+            emp_b (:,:) = emp(:,:)
+            sfx_b (:,:) = sfx(:,:)
          ENDIF
       ENDIF
@@ -340,5 +352,5 @@
          ! CALL iom_rstput( kt, nitrst, numrow, 'qsr_b'  , qsr  )
          CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp  )
-         CALL iom_rstput( kt, nitrst, numrow, 'emps_b' , emps )
+         CALL iom_rstput( kt, nitrst, numrow, 'sfx_b' , sfx )
       ENDIF
 
@@ -348,5 +360,7 @@
       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
          CALL iom_put( "empmr" , emp  - rnf )                   ! upward water flux
-         CALL iom_put( "empsmr", emps - rnf )                   ! c/d water flux
+         CALL iom_put( "saltflx", sfx  )                        ! downward salt flux  
+                                                                ! (includes virtual salt flux beneath ice 
+                                                                ! in linear free surface case)
          CALL iom_put( "qt"    , qns  + qsr )                   ! total heat flux 
          CALL iom_put( "qns"   , qns        )                   ! solar heat flux
@@ -365,5 +379,5 @@
          CALL prt_ctl(tab2d_1=fr_i             , clinfo1=' fr_i     - : ', mask1=tmask, ovlap=1 )
          CALL prt_ctl(tab2d_1=(emp-rnf)        , clinfo1=' emp-rnf  - : ', mask1=tmask, ovlap=1 )
-         CALL prt_ctl(tab2d_1=(emps-rnf)       , clinfo1=' emps-rnf - : ', mask1=tmask, ovlap=1 )
+         CALL prt_ctl(tab2d_1=(sfx-rnf)        , clinfo1=' sfx-rnf  - : ', mask1=tmask, ovlap=1 )
          CALL prt_ctl(tab2d_1=qns              , clinfo1=' qns      - : ', mask1=tmask, ovlap=1 )
          CALL prt_ctl(tab2d_1=qsr              , clinfo1=' qsr      - : ', mask1=tmask, ovlap=1 )
