Index: branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limsbc_2.F90
===================================================================
--- branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 1858)
+++ branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 1859)
@@ -4,8 +4,9 @@
    !!           computation of the flux at the sea ice/ocean interface
    !!======================================================================
-   !! History : LIM  ! 2000-01 (H. Goosse) Original code
-   !!           2.0  ! 2002-07 (C. Ethe, G. Madec) re-writing F90
-   !!            -   ! 2006-07 (G. Madec) surface module
-   !!           2.1  ! 2010-05  (Y. Aksenov, M. Vancoppenolle, G. Madec) add heat content exchanges
+   !! History :  1.0  !  2000-01  (H. Goosse) Original code
+   !!            2.0  !  2002-07  (C. Ethe, G. Madec) re-writing F90
+   !!             -   !  2006-07  (G. Madec) surface module
+   !!             -   !  2008-07  (C. Talandier,G.  Madec) 2D fields for soce and sice
+   !!            2.1  !  2010-05  (Y. Aksenov G. Madec) salt flux + heat associated with emp
    !!----------------------------------------------------------------------
 #if defined key_lim2
@@ -17,15 +18,15 @@
    USE par_oce          ! ocean parameters
    USE dom_oce          ! ocean domain
-   USE sbc_ice          ! surface boundary condition
-   USE sbc_oce          ! surface boundary condition
+   USE sbc_ice          ! ice   surface boundary condition
+   USE sbc_oce          ! ocean surface boundary condition
    USE phycst           ! physical constants
-   USE ice_2            ! LIM sea-ice variables
+   USE albedo           ! albedo parameters
+   USE ice_2            ! LIM-2 sea-ice variables
 
    USE lbclnk           ! ocean lateral boundary condition
    USE in_out_manager   ! I/O manager
+   USE iom              ! 
+   USE prtctl           ! Print control
    USE diaar5, ONLY :   lk_diaar5
-   USE iom              ! 
-   USE albedo           ! albedo parameters
-   USE prtctl           ! Print control
    USE cpl_oasis3, ONLY : lk_cpl
 
@@ -33,16 +34,15 @@
    PRIVATE
 
-   PUBLIC lim_sbc_2     ! called by sbc_ice_lim_2
+   PUBLIC   lim_sbc_2   ! called by sbc_ice_lim_2
 
    REAL(wp)  ::   epsi16 = 1.e-16  ! constant values
    REAL(wp)  ::   rzero  = 0.e0    
    REAL(wp)  ::   rone   = 1.e0
-   REAL(wp), DIMENSION(jpi,jpj)  ::   soce_r
-   REAL(wp), DIMENSION(jpi,jpj)  ::   sice_r
+   REAL(wp), DIMENSION(jpi,jpj)  ::   soce_r, sice_r   ! ocean and ice 2D constant salinity fields (used if lk_vvl=F)
 
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006) 
+   !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -62,8 +62,8 @@
       !!              - Update 
       !!     
-      !! ** Outputs : - qsr     : sea heat flux:     solar 
-      !!              - qns     : sea heat flux: non solar
-      !!              - emp     : freshwater budget: volume flux 
-      !!              - emps    : freshwater budget: concentration/dillution 
+      !! ** Outputs : - qsr     : solar heat flux
+      !!              - qns     : non-solar heat flux including heat content of mass flux
+      !!              - emp     : mass flux
+      !!              - emps    : salt flux due to Freezing/Melting 
       !!              - utau    : sea surface i-stress (ocean referential)
       !!              - vtau    : sea surface j-stress (ocean referential)
@@ -75,60 +75,62 @@
       !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108.
       !!---------------------------------------------------------------------
-      INTEGER ::   kt    ! number of iteration
+      INTEGER, INTENT(in) ::   kt    ! number of iteration
       !!
-      INTEGER  ::   ji, jj           ! dummy loop indices
-      INTEGER  ::   ifvt, i1mfr, idfr               ! some switches
-      INTEGER  ::   iflt, ial, iadv, ifral, ifrdv
-      INTEGER  ::   ii0, ii1, ij0, ij1  ! temporary integers
-      REAL(wp) ::   zrdtir           ! 1. / rdt_ice
-      REAL(wp) ::   zqsr  , zqns     ! solar & non solar heat flux
-      REAL(wp) ::   zinda            ! switch for testing the values of ice concentration
-      REAL(wp) ::   zfons            ! salt exchanges at the ice/ocean interface
-      REAL(wp) ::   zemp             ! freshwater exchanges at the ice/ocean interface
-      REAL(wp) ::   zfrldu, zfrldv   ! lead fraction at U- & V-points
-      REAL(wp) ::   zutau , zvtau    ! lead fraction at U- & V-points
-      REAL(wp) ::   zu_io , zv_io    ! 2 components of the ice-ocean velocity
-! interface 2D --> 3D
-      REAL(wp), DIMENSION(jpi,jpj,1) ::   zalb     ! albedo of ice under overcast sky
-      REAL(wp), DIMENSION(jpi,jpj,1) ::   zalbp    ! albedo of ice under clear sky
+      INTEGER  ::   ji, jj                     ! dummy loop indices
+      INTEGER  ::   ifvt, idfr , iadv, i1mfr   ! local integers
+      INTEGER  ::   iflt, ifrdv, ial , ifral   !   -      -
+      INTEGER  ::   ii0, ii1, ij0, ij1         !   -      -
+      REAL(wp) ::   zqsr, zqns, zqhc, zemp     ! local scalars
+      REAL(wp) ::   zinda, zswitch, zcd        !   -      -
+      REAL(wp) ::   zfrldu, zutau, zu_io       !   -      -
+      REAL(wp) ::   zfrldv, zvtau, zv_io       !   -      -
+      REAL(wp) ::   zemp_snw, zfmm, zfsalt     !   -      -
       REAL(wp) ::   zsang, zmod, zztmp, zfm
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztio_u, ztio_v   ! component of ocean stress below sea-ice at I-point
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztiomi           ! module    of ocean stress below sea-ice at I-point
-      REAL(wp), DIMENSION(jpi,jpj) ::   zqnsoce          ! save qns before its modification by ice model
-
+      REAL(wp), DIMENSION(jpi,jpj,1) ::   zalb, zalbp    ! 3D workspace
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztio_u, ztio_v   ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztiomi, zqnsoce  !  -      -
       !!---------------------------------------------------------------------
-     
-      zrdtir = 1. / rdt_ice
-      
+           
       IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'lim_sbc_2 : LIM 2.0 sea-ice - surface boundary condition'
          IF(lwp) WRITE(numout,*) '~~~~~~~~~   '
-
-         soce_r(:,:) = soce
-         sice_r(:,:) = sice
-         !
-         IF( cp_cfg == "orca"  .AND. jp_cfg == 2 ) THEN
-            !                                        ! =======================
-            !                                        !  ORCA_R2 configuration
-            !                                        ! =======================
-            ii0 = 145   ;   ii1 = 180        ! Baltic Sea
+         !                              ! 2D fields for constant ice and ocean salinities
+         soce_r(:,:) = soce             !    in order to use different value in the Baltic sea
+         sice_r(:,:) = sice             !    which is much less salty than polar regions
+         !
+         IF( cp_cfg == "orca" ) THEN    ! ORCA configuration
+            IF( jp_cfg == 2       ) THEN     !  ORCA_R2 configuration
+            ii0 = 145   ;   ii1 = 180              ! Baltic Sea
             ij0 = 113   ;   ij1 = 130   ;   soce_r(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0
-                                            sice_r(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 2.e0
+                                            sice_r(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
+!!gm add here the R1 R05 and R025  cases
+!!          ELSEIF( jp_cfg == 1   ) THEN           !  ORCA_R1   configuration
+!!          ELSEIF( jp_cfg == 05  ) THEN           !  ORCA_R05  configuration
+!!          ELSEIF( jp_cfg == 025 ) THEN           !  ORCA_R025 configuration
+!!
+!!gm or better introduce the baltic change as a function of lat/lon of the baltic sea
+!!end gm
+            ENDIF
          ENDIF
          !
       ENDIF
 
-      !------------------------------------------!
-      !      heat flux at the ocean surface      !
-      !------------------------------------------!
-
-!!gm
-!!gm CAUTION   
-!!gm re-verifies the non solar expression, especially over open ocen
-!!gm
-      zqnsoce(:,:) = qns(:,:)
+      zqnsoce(:,:) = qns(:,:)      ! save non-solar flux prior to its modification by ice-ocean fluxes (diag.)
+      !
+      zswitch = 1       ! standard levitating sea-ice : salt exchange only
+      !
+!!gm ice embedment
+!      SELECT CASE( nn_ice_embd )       ! levitating/embedded sea-ice
+!      CASE( 0    )   ;   zswitch = 1       ! standard levitating sea-ice : salt exchange only
+!      CASE( 1, 2 )   ;   zswitch = 0       ! other levitating sea-ice or embedded sea-ice : salt and volume fluxes
+!      END SELECT                           !    
+!!gm end embedment
+
       DO jj = 1, jpj
          DO ji = 1, jpi
+            !                          !------------------------------------------!
+            !                          !      heat flux at the ocean surface      !
+            !                          !------------------------------------------!
             zinda   = 1.0   - MAX( rzero , SIGN( rone, - ( 1.0 - pfrld(ji,jj) )   ) )
             ifvt    = zinda * MAX( rzero , SIGN( rone,  - phicif(ji,jj)           ) )
@@ -138,42 +140,42 @@
             ial     = ifvt   * i1mfr + ( 1 - ifvt ) * idfr
             iadv    = ( 1  - i1mfr ) * zinda
-            ifral   = ( 1  - i1mfr * ( 1 - ial ) )   
-            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv 
-
-!!$            zinda   = 1.0 - AINT( pfrld(ji,jj) )                   !   = 0. if pure ocean else 1. (at previous time)
-!!$
-!!$            i1mfr   = 1.0 - AINT(  frld(ji,jj) )                   !   = 0. if pure ocean else 1. (at current  time)
-!!$
-!!$            IF( phicif(ji,jj) <= 0. ) THEN   ;   ifvt = zinda      !   = 1. if (snow and no ice at previous time) else 0. ???
-!!$            ELSE                             ;   ifvt = 0.
-!!$            ENDIF
-!!$
-!!$            IF( frld(ji,jj) >= pfrld(ji,jj) ) THEN   ;   idfr = 0.  !   = 0. if lead fraction increases from previous to current
-!!$            ELSE                                     ;   idfr = 1.   
-!!$            ENDIF
-!!$
-!!$            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )    !   = 1. if ice (not only snow) at previous and pure ocean at current
-!!$
+            ifral   = ( 1  - i1mfr * ( 1 - ial ) )
+            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv
+
+!!gm  attempt to understand and comment the tricky flags used here....
+!
+!gm      zinda   = 1.0 - AINT( pfrld(ji,jj) )    ! = 0. free-ice ocean else 1. (after ice adv, but before ice thermo)
+!gm      i1mfr   = 1.0 - AINT(  frld(ji,jj) )    ! = 0. free-ice ocean else 1. (after ice thermo)
+!
+!gm      IF( phicif(ji,jj) <= 0. ) THEN   ;   ifvt = zinda   ! = 1. if (snow and no ice at previous time) else 0. ???
+!gm      ELSE                             ;   ifvt = 0.      ! correspond to a overmelting of snow in surface ablation
+!gm      ENDIF                                               ! 
+!
+!gm      IF( frld(ji,jj) >= pfrld(ji,jj) ) THEN   ;   idfr = 0.  !   = 0. if lead fraction increases due to ice thermo
+!gm      ELSE                                     ;   idfr = 1.   
+!gm      ENDIF
+!
+!!$      iflt    = zinda  * (1 - i1mfr) * (1 - ifvt ) ! = 1. if ice (not only snow) at previous and pure ocean at current
+!
 !!$            ial     = ifvt   * i1mfr    +    ( 1 - ifvt ) * idfr
 !!$!                 snow no ice   ice         ice or nothing  lead fraction increases
 !!$!                 at previous   now           at previous
 !!$!                -> ice aera increases  ???         -> ice aera decreases ???
-!!$
+!
 !!$            iadv    = ( 1  - i1mfr ) * zinda  
 !!$!                     pure ocean      ice at
 !!$!                     at current      previous
 !!$!                        -> = 1. if ice disapear between previous and current
-!!$
+!
 !!$            ifral   = ( 1  - i1mfr * ( 1 - ial ) )  
 !!$!                            ice at     ???
 !!$!                            current         
 !!$!                         -> ???
-!!$ 
+!
 !!$            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv 
 !!$!                                                    ice disapear                           
-!!$
-!!$
-
-            !   computation the solar flux at ocean surface
+!
+            !
+            ! - computation the solar flux at ocean surface
 #if defined key_coupled 
             zqsr = qsr_tot(ji,jj) + ( fstric(ji,jj) - qsr_ice(ji,jj,1) ) * ( 1.0 - pfrld(ji,jj) )
@@ -181,60 +183,56 @@
             zqsr = pfrld(ji,jj) * qsr(ji,jj)  + ( 1.  - pfrld(ji,jj) ) * fstric(ji,jj)
 #endif            
-            !  computation the non solar heat flux at ocean surface
+            !
+            ! - computation the non solar heat flux at ocean surface
             zqns    =  - ( 1. - thcm(ji,jj) ) * zqsr   &   ! part of the solar energy used in leads
-               &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                            &
-               &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * zrdtir    &
-               &       + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) ) * zrdtir
-
-            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     ! ???
-            
-            qsr  (ji,jj) = zqsr                                          ! solar heat flux 
-            qns  (ji,jj) = zqns - fdtcn(ji,jj)                           ! non solar heat flux
+               &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                                &
+               &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdt_ice    &
+               &       + ifrdv   * (       qfvbq(ji,jj) +             qdtcn(ji,jj) ) * r1_rdt_ice
+
+            ! - store residual heat flux (put in the ocean at the next time-step)
+            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)   ! ???
+            !
+            ! - heat content of mass exchanged between ocean and sea-ice
+            zqhc = ( rdq_snw(ji,jj) + rdq_ice(ji,jj) ) * r1_rdt_ice    ! heat flux due to sown & ice heat content exchanges
+            !            
+            qsr(ji,jj) = zqsr                                          ! solar heat flux 
+            qns(ji,jj) = zqns - fdtcn(ji,jj) + zqhc                    ! non solar heat flux
+  
+            !                          !------------------------------------------!
+            !                          !      mass flux at the ocean surface      !
+            !                          !------------------------------------------!
+            !
+            ! mass flux at the ocean-atmosphere interface (open ocean fraction = leads area)
+#if defined key_coupled
+            !                                                       ! coupled mode: 
+            zemp = + emp_tot(ji,jj)                              &       ! net mass flux over the grid cell (ice+ocean area)
+               &   - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )              ! minus the mass flux intercepted by sea-ice
+#else
+            !                                                       ! forced  mode: 
+            zemp = + emp(ji,jj)     *         frld(ji,jj)      &         ! mass flux over open ocean fraction 
+               &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )    &         ! liquid precip. over ice reaches directly the ocean
+               &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )    &         ! snow is intercepted by sea-ice (previous frld)
+#endif            
+            !
+            ! mass flux at the ocean/ice interface (sea ice fraction)
+            zemp_snw = rdm_snw(ji,jj) * r1_rdt_ice                  ! snow melting = pure water that enters the ocean
+            zfmm     = rdm_ice(ji,jj) * r1_rdt_ice                  ! Freezing minus Melting (F-M)
+
+            ! salt flux at the ice/ocean interface (sea ice fraction) [PSU*kg/m2/s]
+            zfsalt = - sice_r(ji,jj) * zfmm                         ! F-M salt exchange
+            zcd    =   soce_r(ji,jj) * zfmm                         ! concentration/dilution term due to F-M
+            !
+            ! salt flux only       : add concentration dilution term in salt flux  and no  F-M term in volume flux
+            ! salt and mass fluxes : non concentartion dilution term in salt flux  and add F-M term in volume flux
+            emps(ji,jj) = zfsalt +                  zswitch  * zcd   ! salt flux (+ C/D if no ice/ocean mass exchange)
+            emp (ji,jj) = zemp   + zemp_snw + ( 1.- zswitch) * zfmm  ! mass flux (- F/M mass flux if no ice/ocean mass exchange)
+            !
          END DO
       END DO
+
 
       CALL iom_put( 'hflx_ice_cea', - fdtcn(:,:) )      
       CALL iom_put( 'qns_io_cea', qns(:,:) - zqnsoce(:,:) * pfrld(:,:) )      
       CALL iom_put( 'qsr_io_cea', fstric(:,:) * (1. - pfrld(:,:)) )
-
-      !------------------------------------------!
-      !      mass flux at the ocean surface      !
-      !------------------------------------------!
-
-!!gm
-!!gm CAUTION   
-!!gm re-verifies the emp & emps expression, especially the absence of 1-frld on zfm
-!!gm
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            
-#if defined key_coupled
-          zemp = emp_tot(ji,jj) - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )    &   ! 
-             &   + rdm_snw(ji,jj) * zrdtir                                      !  freshwaterflux due to snow melting 
-#else
-!!$            !  computing freshwater exchanges at the ice/ocean interface
-!!$            zpme = - evap(ji,jj)    *   frld(ji,jj)           &   !  evaporation over oceanic fraction
-!!$               &   + tprecip(ji,jj)                           &   !  total precipitation
-!!$               &   - sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )   &   !  remov. snow precip over ice
-!!$               &   - rdm_snw(ji,jj) / rdt_ice                     !  freshwaterflux due to snow melting 
-            !  computing freshwater exchanges at the ice/ocean interface
-            zemp = + emp(ji,jj)     *         frld(ji,jj)      &   !  e-p budget over open ocean fraction 
-               &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )    &   !  liquid precipitation reaches directly the ocean
-               &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !  taking into account change in ice cover within the time step
-               &   + rdm_snw(ji,jj) * zrdtir                       !  freshwaterflux due to snow melting 
-               !                                                   !  ice-covered fraction:
-#endif            
-
-            !  computing salt exchanges at the ice/ocean interface
-            zfons =  ( soce_r(ji,jj) - sice_r(ji,jj) ) * ( rdm_ice(ji,jj) * zrdtir ) 
-            
-            !  converting the salt flux from ice to a freshwater flux from ocean
-            zfm  = zfons / ( sss_m(ji,jj) + epsi16 )
-            
-            emps(ji,jj) = zemp + zfm      ! surface ocean concentration/dilution effect (use on SSS evolution)
-            emp (ji,jj) = zemp            ! surface ocean volume flux (use on sea-surface height evolution)
-
-         END DO
-      END DO
 
       IF( lk_diaar5 ) THEN
@@ -250,5 +248,5 @@
       IF ( ln_limdyn ) THEN                        ! Update the stress over ice-over area (only in ice-dynamic case)
          !                                         ! otherwise the atmosphere-ocean stress is used everywhere
-
+         !
          ! ... ice stress over ocean with a ice-ocean rotation angle (at I-point)
 !CDIR NOVERRCHK
@@ -290,22 +288,22 @@
             END DO
          END DO
-
-         ! boundary condition on the stress (utau,vtau,taum)
-         CALL lbc_lnk( utau, 'U', -1. )
-         CALL lbc_lnk( vtau, 'V', -1. )
+         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )         ! lateral boundary conditions
          CALL lbc_lnk( taum, 'T',  1. )
-
+         !
       ENDIF
 
       !-----------------------------------------------!
-      !   Coupling variables                          !
+      !   Storing the transmitted variables           !
       !-----------------------------------------------!
-
-      IF ( lk_cpl ) THEN           
-         ! Ice surface temperature 
+!!gm where this is done ?????   ==>>> limthd_2   not logic ???
+!!gm      fr_i(:,:) = 1.0 - frld(:,:)       ! sea-ice fraction
+!!gm
+
+      IF ( lk_cpl ) THEN      ! coupled mode :
          tn_ice(:,:,1) = sist(:,:)          ! sea-ice surface temperature       
-         ! Computation of snow/ice and ocean albedo
+         !                                  ! snow/ice and ocean albedo
          CALL albedo_ice( tn_ice, reshape( hicif, (/jpi,jpj,1/) ), reshape( hsnif, (/jpi,jpj,1/) ), zalbp, zalb )
          alb_ice(:,:,1) =  0.5 * ( zalbp(:,:,1) + zalb (:,:,1) )   ! Ice albedo (mean clear and overcast skys)
+         !
          CALL iom_put( "icealb_cea", alb_ice(:,:,1) * fr_i(:,:) )  ! ice albedo
       ENDIF
@@ -318,6 +316,6 @@
          CALL prt_ctl(tab2d_1=fr_i  , clinfo1=' lim_sbc: fr_i   : ', tab2d_2=tn_ice(:,:,1), clinfo2=' tn_ice  : ')
       ENDIF 
-   
-    END SUBROUTINE lim_sbc_2
+      !
+   END SUBROUTINE lim_sbc_2
 
 #else
