Index: branches/devmercator2010_1/NEMO/LIM_SRC_2/limsbc_2.F90
===================================================================
--- branches/devmercator2010_1/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 2092)
+++ branches/devmercator2010_1/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 2135)
@@ -4,12 +4,12 @@
    !!           computation of the flux at the sea ice/ocean interface
    !!======================================================================
-   !! History : 00-01 (H. Goosse) Original code
-   !!           02-07 (C. Ethe, G. Madec) re-writing F90
-   !!           06-07 (G. Madec) surface module
+   !! History :  LIM  ! 2000-01 (H. Goosse) Original code
+   !!            1.0  ! 2002-07 (C. Ethe, G. Madec) re-writing F90
+   !!            3.0  ! 2006-07 (G. Madec) surface module
+   !!            3.3  ! 2009-05 (G.Garric, C. Bricaud) addition of the lim2_evp case
    !!----------------------------------------------------------------------
 #if defined key_lim2
    !!----------------------------------------------------------------------
    !!   'key_lim2'                                    LIM 2.0 sea-ice model
-   !!----------------------------------------------------------------------
    !!----------------------------------------------------------------------
    !!   lim_sbc_2  : flux at the ice / ocean interface
@@ -17,13 +17,13 @@
    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          ! surface boundary condition: ice
+   USE sbc_oce          ! surface boundary condition: ocean
    USE phycst           ! physical constants
-   USE ice_2            ! LIM sea-ice variables
-
-   USE lbclnk           ! ocean lateral boundary condition
+   USE ice_2            ! LIM2: ice variables
+
+   USE lbclnk           ! ocean lateral boundary condition - MPP exchanges
    USE in_out_manager   ! I/O manager
    USE diaar5, ONLY :   lk_diaar5
-   USE iom              ! 
+   USE iom              ! IOM library
    USE albedo           ! albedo parameters
    USE prtctl           ! Print control
@@ -33,20 +33,20 @@
    PRIVATE
 
-   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
+   PUBLIC   lim_sbc_2   ! called by sbc_ice_lim_2
+
+   REAL(wp)  ::   r1_rdtice                    ! constant values
+   REAL(wp)  ::   epsi16 = 1.e-16              !     -      -
+   REAL(wp)  ::   rzero  = 0.e0                !     -      -
+   REAL(wp)  ::   rone   = 1.e0                !     -      -
+   !
+   REAL(wp), DIMENSION(jpi,jpj) ::   soce_r, sice_r   ! constant SSS and ice salinity used in levitating sea-ice case
 
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006) 
+   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010)
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -78,26 +78,16 @@
       !!
       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
-      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
-
+      INTEGER  ::   ii0, ii1, ij0, ij1         ! local integers
+      INTEGER  ::   ifvt, i1mfr, idfr, iflt    !   -       -
+      INTEGER  ::   ial, iadv, ifral, ifrdv    !   -       -
+      REAL(wp) ::   zqsr, zqns, zsang, zmod, zfm   ! local scalars
+      REAL(wp) ::   zinda, zfons, zemp, zztmp      !   -      -
+      REAL(wp) ::   zfrldu, zutau, zu_io           !   -      -
+      REAL(wp) ::   zfrldv, zvtau, zv_io           !   -      -
+      REAL(wp), DIMENSION(jpi,jpj)   ::   ztio_u, ztio_v    ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj)   ::   ztiomi, zqnsoce   !  -     -
+      REAL(wp), DIMENSION(jpi,jpj,1) ::   zalb, zalbp   ! 2D/3D workspace
       !!---------------------------------------------------------------------
      
-      zrdtir = 1. / rdt_ice
       
       IF( kt == nit000 ) THEN
@@ -105,12 +95,11 @@
          IF(lwp) WRITE(numout,*) 'lim_sbc_2 : LIM 2.0 sea-ice - surface boundary condition'
          IF(lwp) WRITE(numout,*) '~~~~~~~~~   '
-
+         !
+         r1_rdtice = 1. / rdt_ice
+         !
          soce_r(:,:) = soce
          sice_r(:,:) = sice
          !
-         IF( cp_cfg == "orca"  .AND. jp_cfg == 2 ) THEN
-            !                                        ! =======================
-            !                                        !  ORCA_R2 configuration
-            !                                        ! =======================
+         IF( cp_cfg == "orca"  .AND. 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
@@ -175,5 +164,5 @@
 !!$
 
-            !   computation the solar flux at ocean surface
+            ! 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,12 +170,12 @@
             zqsr = pfrld(ji,jj) * qsr(ji,jj)  + ( 1.  - pfrld(ji,jj) ) * fstric(ji,jj)
 #endif            
-            !  computation the non solar heat flux at ocean surface
+            ! 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
-
+               &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                              &
+               &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice   &
+               &       + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) )                   * r1_rdtice
+            !
             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
@@ -194,52 +183,39 @@
       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(:,:)) )
+      CALL iom_put( 'hflx_ice_cea', - fdtcn(:,:)                         )      
+      CALL iom_put( 'qns_io_cea'  , qns(:,:) - zqnsoce(:,:) * pfrld(:,:) )      
+      CALL iom_put( 'qsr_io_cea'  , fstric(:,:) * ( 1.e0 - 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) )    &   ! 
-             &   + rdmsnif(ji,jj) * zrdtir                                      !  freshwaterflux due to snow melting 
+            ! freshwater exchanges at the ice-atmosphere / ocean interface (coupled mode)
+            zemp = emp_tot(ji,jj) - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )   &   ! atmosphere-ocean freshwater flux
+               &                  + rdmsnif(ji,jj) * r1_rdtice                   ! freshwater flux 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
-!!$               &   - rdmsnif(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
-               &   + rdmsnif(ji,jj) * zrdtir                       !  freshwaterflux due to snow melting 
-               !                                                   !  ice-covered fraction:
+            ! freshwater exchanges at the ice-atmosphere / ocean interface (forced mode)
+            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) )   &   ! (account for change in ice cover within the timestep
+               &   + rdmsnif(ji,jj) * r1_rdtice                   ! freshwaterflux due to snow melting 
 #endif            
-
-            !  computing salt exchanges at the ice/ocean interface
-            zfons =  ( soce_r(ji,jj) - sice_r(ji,jj) ) * ( rdmicif(ji,jj) * zrdtir ) 
-            
-            !  converting the salt flux from ice to a freshwater flux from ocean
+            ! salt exchanges at the ice/ocean interface
+            zfons =  ( soce_r(ji,jj) - sice_r(ji,jj) ) * ( rdmicif(ji,jj) * r1_rdtice ) 
+            !
+            ! convert the salt flux from ice into 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
-         CALL iom_put( 'isnwmlt_cea'  ,                 rdmsnif(:,:) * zrdtir )
-         CALL iom_put( 'fsal_virt_cea',   soce_r(:,:) * rdmicif(:,:) * zrdtir )
-         CALL iom_put( 'fsal_real_cea', - sice_r(:,:) * rdmicif(:,:) * zrdtir )
+         CALL iom_put( 'isnwmlt_cea'  ,                 rdmsnif(:,:) * r1_rdtice )
+         CALL iom_put( 'fsal_virt_cea',   soce_r(:,:) * rdmicif(:,:) * r1_rdtice )
+         CALL iom_put( 'fsal_real_cea', - sice_r(:,:) * rdmicif(:,:) * r1_rdtice )
       ENDIF
 
@@ -275,6 +251,11 @@
             DO ji = 2, jpim1   ! NO vector opt.
                ! ... components of ice-ocean stress at U and V-points  (from I-point values)
-               zutau  = 0.5 * ( ztio_u(ji+1,jj) + ztio_u(ji+1,jj+1) )
+#if defined key_lim2_vp
+               zutau  = 0.5 * ( ztio_u(ji+1,jj) + ztio_u(ji+1,jj+1) )      ! VP rheology
                zvtau  = 0.5 * ( ztio_v(ji,jj+1) + ztio_v(ji+1,jj+1) )
+#else
+               zutau  = ztio_u(ji,jj)                                      ! EVP rheology
+               zvtau  = ztio_v(ji,jj)
+#endif
                ! ... open-ocean (lead) fraction at U- & V-points (from T-point values)
                zfrldu = 0.5 * ( frld(ji,jj) + frld(ji+1,jj  ) )
@@ -290,20 +271,15 @@
             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 condition 
          CALL lbc_lnk( taum, 'T',  1. )
 
       ENDIF
 
+      IF( lk_cpl ) THEN           
       !-----------------------------------------------!
       !   Coupling variables                          !
       !-----------------------------------------------!
-
-      IF ( lk_cpl ) THEN           
-         ! Ice surface temperature 
-         tn_ice(:,:,1) = sist(:,:)          ! sea-ice surface temperature       
-         ! Computation of snow/ice and ocean albedo
+         tn_ice(:,:,1) = sist(:,:)           ! sea-ice surface temperature       
+         !                                   ! 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)
@@ -318,6 +294,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
