Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis3.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis3.F90	(revision 5891)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis3.F90	(revision 5901)
@@ -4,14 +4,15 @@
    !! Coupled O/A : coupled ocean-atmosphere case using OASIS3-MCT
    !!=====================================================================
-   !! History :   
-   !!   9.0  !  04-06  (R. Redler, NEC Laboratories Europe, Germany) Original code
-   !!   " "  !  04-11  (R. Redler, NEC Laboratories Europe; N. Keenlyside, W. Park, IFM-GEOMAR, Germany) revision
-   !!   " "  !  04-11  (V. Gayler, MPI M&D) Grid writing
-   !!   " "  !  05-08  (R. Redler, W. Park) frld initialization, paral(2) revision
-   !!   " "  !  05-09  (R. Redler) extended to allow for communication over root only
-   !!   " "  !  06-01  (W. Park) modification of physical part
-   !!   " "  !  06-02  (R. Redler, W. Park) buffer array fix for root exchange
-   !!   3.4  !  11-11  (C. Harris) Changes to allow mutiple category fields
-   !!----------------------------------------------------------------------
+   !! History :  1.0  !  2004-06  (R. Redler, NEC Laboratories Europe, Germany) Original code
+   !!             -   !  2004-11  (R. Redler, NEC Laboratories Europe; N. Keenlyside, W. Park, IFM-GEOMAR, Germany) revision
+   !!             -   !  2004-11  (V. Gayler, MPI M&D) Grid writing
+   !!            2.0  !  2005-08  (R. Redler, W. Park) frld initialization, paral(2) revision
+   !!             -   !  2005-09  (R. Redler) extended to allow for communication over root only
+   !!             -   !  2006-01  (W. Park) modification of physical part
+   !!             -   !  2006-02  (R. Redler, W. Park) buffer array fix for root exchange
+   !!            3.4  !  2011-11  (C. Harris) Changes to allow mutiple category fields
+   !!            3.6  !  2014-11  (S. Masson) OASIS3-MCT
+   !!----------------------------------------------------------------------
+   
    !!----------------------------------------------------------------------
    !!   'key_oasis3'                    coupled Ocean/Atmosphere via OASIS3-MCT
@@ -20,6 +21,6 @@
    !!   cpl_init     : initialization of coupled mode communication
    !!   cpl_define   : definition of grid and fields
-   !!   cpl_snd     : snd out fields in coupled mode
-   !!   cpl_rcv     : receive fields in coupled mode
+   !!   cpl_snd      : snd out fields in coupled mode
+   !!   cpl_rcv      : receive fields in coupled mode
    !!   cpl_finalize : finalize the coupled mode communication
    !!----------------------------------------------------------------------
@@ -99,6 +100,6 @@
       !! ** Method  :   OASIS3 MPI communication 
       !!--------------------------------------------------------------------
-      CHARACTER(len = *), INTENT(in) ::   cd_modname   ! model name as set in namcouple file
-      INTEGER          , INTENT(out) ::   kl_comm      ! local communicator of the model
+      CHARACTER(len = *), INTENT(in   ) ::   cd_modname   ! model name as set in namcouple file
+      INTEGER           , INTENT(  out) ::   kl_comm      ! local communicator of the model
       !!--------------------------------------------------------------------
 
@@ -163,5 +164,4 @@
          CALL oasis_abort ( ncomp_id, 'cpl_define', 'nsnd is larger than nmaxfld, increase nmaxfld')   ;   RETURN
       ENDIF
-
       !
       ! ... Define the shape for the area that excludes the halo
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90	(revision 5891)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90	(revision 5901)
@@ -195,5 +195,4 @@
 
       DO jj = 2, jpjm1
-!CDIR NOVERRCHK
          DO ji = fs_2, jpi   ! vector opt.
 
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 5891)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 5901)
@@ -80,5 +80,4 @@
    INTEGER , PUBLIC, PARAMETER ::   jp_mfs     = 6        !: MFS  bulk                     formulation
    INTEGER , PUBLIC, PARAMETER ::   jp_none    = 7        !: for OPA when doing coupling via SAS module
-   INTEGER , PUBLIC, PARAMETER ::   jp_esopa   = -1       !: esopa test, ALL formulations
    
    !!----------------------------------------------------------------------
@@ -200,7 +199,5 @@
       !!---------------------------------------------------------------------
       zcoef = 0.5 / ( zrhoa * zcdrag ) 
-!CDIR NOVERRCHK
       DO jj = 2, jpjm1
-!CDIR NOVERRCHK
          DO ji = fs_2, fs_jpim1   ! vect. opt.
             ztx = utau(ji-1,jj  ) + utau(ji,jj) 
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcana.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcana.F90	(revision 5891)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcana.F90	(revision 5901)
@@ -279,7 +279,5 @@
       ! module of wind stress and wind speed at T-point
       zcoef = 1. / ( zrhoa * zcdrag ) 
-!CDIR NOVERRCHK
       DO jj = 2, jpjm1
-!CDIR NOVERRCHK
          DO ji = fs_2, fs_jpim1   ! vect. opt.
             ztx = utau(ji-1,jj  ) + utau(ji,jj) 
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcapr.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcapr.F90	(revision 5891)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcapr.F90	(revision 5901)
@@ -62,5 +62,5 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in)::   kt   ! ocean time step
-      !!
+      !
       INTEGER            ::   ierror  ! local integer 
       INTEGER            ::   ios     ! Local integer output status for namelist read
@@ -71,5 +71,4 @@
       NAMELIST/namsbc_apr/ cn_dir, sn_apr, ln_ref_apr, rn_pref, ln_apr_obc
       !!----------------------------------------------------------------------
-      !
       !
       !                                         ! -------------------- !
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90	(revision 5891)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90	(revision 5901)
@@ -243,7 +243,5 @@
       !   momentum fluxes  (utau, vtau )   !
       !------------------------------------!
-!CDIR COLLAPSE
       utau(:,:) = sf(jp_utau)%fnow(:,:,1)
-!CDIR COLLAPSE
       vtau(:,:) = sf(jp_vtau)%fnow(:,:,1)
 
@@ -251,7 +249,5 @@
       !   wind stress module (taum )       !
       !------------------------------------!
-!CDIR NOVERRCHK
       DO jj = 2, jpjm1
-!CDIR NOVERRCHK
          DO ji = fs_2, fs_jpim1   ! vector opt.
             ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
@@ -268,5 +264,4 @@
       !   store the wind speed  (wndm )    !
       !------------------------------------!
-!CDIR COLLAPSE
       wndm(:,:) = sf(jp_wndm)%fnow(:,:,1)
       wndm(:,:) = wndm(:,:) * tmask(:,:,1)
@@ -281,8 +276,5 @@
       !   Other ocean fluxes   !
       !------------------------!
-!CDIR NOVERRCHK
-!CDIR COLLAPSE
       DO jj = 1, jpj
-!CDIR NOVERRCHK
          DO ji = 1, jpi
             !
@@ -375,9 +367,7 @@
       zcprec = rcp /  rday     ! convert prec ( mm/day ==> m/s)  ==> W/m2
 
-!CDIR COLLAPSE
       emp(:,:) = zqla(:,:) / cevap                                        &   ! freshwater flux
          &     - sf(jp_prec)%fnow(:,:,1) / rday * tmask(:,:,1)
       !
-!CDIR COLLAPSE
       qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                        &   ! Downward Non Solar flux
          &     - zqla(:,:)             * pst(:,:) * zcevap                &   ! remove evap.   heat content at SST in Celcius
@@ -415,4 +405,5 @@
 
 # if defined key_lim2 || defined key_lim3
+
    SUBROUTINE blk_ice_clio_tau
       !!---------------------------------------------------------------------------
@@ -429,15 +420,15 @@
       !
       IF( nn_timing == 1 )  CALL timing_start('blk_ice_clio_tau')
-
+      !
       SELECT CASE( cp_ice_msh )
-
+      !
       CASE( 'C' )                          ! C-grid ice dynamics
-
+         !
          zcoef  = cai / cao                         ! Change from air-sea stress to air-ice stress
          utau_ice(:,:) = zcoef * utau(:,:)
          vtau_ice(:,:) = zcoef * vtau(:,:)
-
+         !
       CASE( 'I' )                          ! I-grid ice dynamics:  I-point (i.e. F-point lower-left corner)
-
+         !
          zcoef  = 0.5_wp * cai / cao                ! Change from air-sea stress to air-ice stress
          DO jj = 2, jpj         ! stress from ocean U- and V-points to ice U,V point
@@ -447,19 +438,21 @@
             END DO
          END DO
-
+         !
          CALL lbc_lnk( utau_ice(:,:), 'I', -1. )   ;   CALL lbc_lnk( vtau_ice(:,:), 'I', -1. )   ! I-point
-
+         !
       END SELECT
-
+      !
       IF(ln_ctl) THEN
          CALL prt_ctl(tab2d_1=utau_ice , clinfo1=' blk_ice_clio: utau_ice : ', tab2d_2=vtau_ice , clinfo2=' vtau_ice : ')
       ENDIF
-
+      !
       IF( nn_timing == 1 )  CALL timing_stop('blk_ice_clio_tau')
-
+      !
    END SUBROUTINE blk_ice_clio_tau
+   
 #endif
 
 # if defined key_lim2 || defined key_lim3
+
    SUBROUTINE blk_ice_clio_flx(  ptsu , palb_cs, palb_os, palb )
       !!---------------------------------------------------------------------------
@@ -520,8 +513,5 @@
       !--------------------------------------------------------------------------------
 
-!CDIR NOVERRCHK
-!CDIR COLLAPSE
       DO jj = 1, jpj
-!CDIR NOVERRCHK
          DO ji = 1, jpi
             ztatm (ji,jj) = sf(jp_tair)%fnow(ji,jj,1)                ! air temperature in Kelvins 
@@ -573,10 +563,7 @@
 
       !                                     ! ========================== !
-      DO jl = 1, jpl                       !  Loop over ice categories  !
+      DO jl = 1, jpl                        !  Loop over ice categories  !
          !                                  ! ========================== !
-!CDIR NOVERRCHK
-!CDIR COLLAPSE
          DO jj = 1 , jpj
-!CDIR NOVERRCHK
             DO ji = 1, jpi
                !-------------------------------------------!
@@ -636,7 +623,5 @@
       ! ----------------------------------------------------------------------------- !
       !
-!CDIR COLLAPSE
       qns_ice(:,:,:) = z_qlw (:,:,:) - z_qsb (:,:,:) - qla_ice (:,:,:)      ! Downward Non Solar flux
-!CDIR COLLAPSE
       tprecip(:,:)   = sf(jp_prec)%fnow(:,:,1) / rday                     ! total precipitation [kg/m2/s]
       !
@@ -644,5 +629,4 @@
       !    Correct the OCEAN non solar flux with the existence of solid precipitation !
       ! ---------------=====--------------------------------------------------------- !
-!CDIR COLLAPSE
       qns(:,:) = qns(:,:)                                                           &   ! update the non-solar heat flux with:
          &     - sprecip(:,:) * lfus                                                  &   ! remove melting solid precip
@@ -782,8 +766,5 @@
       ! Saturated water vapour and vapour pressure
       ! ------------------------------------------
-!CDIR NOVERRCHK
-!CDIR COLLAPSE
       DO jj = 1, jpj
-!CDIR NOVERRCHK
          DO ji = 1, jpi
             ztamr = sf(jp_tair)%fnow(ji,jj,1) - rtt
@@ -814,7 +795,5 @@
       zdaycor  = 1.0 + 0.0013 * SIN( zdist ) + 0.0342 * COS( zdist )
 
-!CDIR NOVERRCHK
       DO jj = 1, jpj
-!CDIR NOVERRCHK
          DO ji = 1, jpi
             !  product of sine (cosine) of latitude and sine (cosine) of solar declination
@@ -837,11 +816,7 @@
 
       ! compute and sum ocean qsr over the daylight (i.e. between sunrise and sunset)
-!CDIR NOVERRCHK   
       DO jt = 1, jp24
          zcoef = FLOAT( jt ) - 0.5
-!CDIR NOVERRCHK     
-!CDIR COLLAPSE
          DO jj = 1, jpj
-!CDIR NOVERRCHK
             DO ji = 1, jpi
                zlha = COS(  zlsrise(ji,jj) - zcoef * zdlha(ji,jj)  )                  ! local hour angle
@@ -862,5 +837,4 @@
       ! Taking into account the ellipsity of the earth orbit, the clouds AND masked if sea-ice cover > 0%
       zcoef1 = srgamma * zdaycor / ( 2. * rpi )
-!CDIR COLLAPSE
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -920,8 +894,5 @@
       ! Saturated water vapour and vapour pressure
       ! ------------------------------------------
-!CDIR NOVERRCHK
-!CDIR COLLAPSE
       DO jj = 1, jpj
-!CDIR NOVERRCHK
          DO ji = 1, jpi           
             ztamr = sf(jp_tair)%fnow(ji,jj,1) - rtt           
@@ -952,7 +923,5 @@
       zdaycor  = 1.0 + 0.0013 * SIN( zdist ) + 0.0342 * COS( zdist )
 
-!CDIR NOVERRCHK
       DO jj = 1, jpj
-!CDIR NOVERRCHK
          DO ji = 1, jpi
             !  product of sine (cosine) of latitude and sine (cosine) of solar declination
@@ -979,11 +948,7 @@
       DO jl = 1, ijpl      !  loop over ice categories  !
          !                 !----------------------------! 
-!CDIR NOVERRCHK   
          DO jt = 1, jp24   
             zcoef = FLOAT( jt ) - 0.5
-!CDIR NOVERRCHK     
-!CDIR COLLAPSE
             DO jj = 1, jpj
-!CDIR NOVERRCHK
                DO ji = 1, jpi
                   zlha = COS(  zlsrise(ji,jj) - zcoef * zdlha(ji,jj)  )                  ! local hour angle
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90	(revision 5891)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90	(revision 5901)
@@ -251,7 +251,5 @@
          ! for basin budget and cooerence
          !--------------------------------------------------
-!CDIR COLLAPSE
-           emp (:,:) = evap(:,:) - sf(jp_prec)%fnow(:,:,1) * tmask(:,:,1)
-!CDIR COLLAPSE
+         emp(:,:) = evap(:,:) - sf(jp_prec)%fnow(:,:,1) * tmask(:,:,1)
 
          CALL iom_put( "qlw_oce",   qbw  )                 ! output downward longwave heat over the ocean
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 5891)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 5901)
@@ -931,7 +931,5 @@
          ! => need to be done only when otx1 was changed
          IF( llnewtx ) THEN
-!CDIR NOVERRCHK
             DO jj = 2, jpjm1
-!CDIR NOVERRCHK
                DO ji = fs_2, fs_jpim1   ! vect. opt.
                   zzx = frcv(jpr_otx1)%z3(ji-1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1)
@@ -961,7 +959,5 @@
          IF( llnewtau ) THEN 
             zcoef = 1. / ( zrhoa * zcdrag ) 
-!CDIR NOVERRCHK
             DO jj = 1, jpj
-!CDIR NOVERRCHK
                DO ji = 1, jpi 
                   frcv(jpr_w10m)%z3(ji,jj,1) = SQRT( frcv(jpr_taum)%z3(ji,jj,1) * zcoef )
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 5891)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 5901)
@@ -131,5 +131,4 @@
          ELSE                  ;   qsr(:,:) =          sf(jp_qsr)%fnow(:,:,1)
          ENDIF
-!CDIR COLLAPSE
          DO jj = 1, jpj                                           ! set the ocean fluxes from read fields
             DO ji = 1, jpi
@@ -145,7 +144,5 @@
          !                                                        ! module of wind stress and wind speed at T-point
          zcoef = 1. / ( zrhoa * zcdrag )
-!CDIR NOVERRCHK
          DO jj = 2, jpjm1
-!CDIR NOVERRCHK
             DO ji = fs_2, fs_jpim1   ! vect. opt.
                ztx = utau(ji-1,jj  ) + utau(ji,jj) 
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 5891)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 5901)
@@ -108,5 +108,5 @@
          !
          IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
-            z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) -  snwice_fmass(:,:) ) ) / area   ! sum over the global domain
+            z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) - snwice_fmass(:,:) ) ) / area   ! sum over the global domain
             zcoef = z_fwf * rcp
             emp(:,:) = emp(:,:) - z_fwf              * tmask(:,:,1)
@@ -162,5 +162,5 @@
             zsurf_pos = glob_sum( e1e2t(:,:)*ztmsk_pos(:,:) )
             !                                                  ! fwf global mean (excluding ocean to ice/snow exchanges) 
-            z_fwf     = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) - snwice_fmass(:,:) ) ) / area
+            z_fwf     = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) - snwice_fmass(:,:) ) ) / area
             !            
             IF( z_fwf < 0._wp ) THEN         ! spread out over >0 erp area to increase evaporation
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 5891)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 5901)
@@ -67,5 +67,4 @@
    PRIVATE
 
-   !! * Routine accessibility
    PUBLIC cice_sbc_init   ! routine called by sbc_init
    PUBLIC cice_sbc_final  ! routine called by sbc_final
@@ -95,6 +94,9 @@
    !! * Substitutions
 #  include "domzgr_substitute.h90"
-
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.7 , NEMO-consortium (2015) 
    !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
@@ -154,9 +156,11 @@
    END SUBROUTINE sbc_ice_cice
 
-   SUBROUTINE cice_sbc_init (ksbc)
+
+   SUBROUTINE cice_sbc_init( ksbc )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE cice_sbc_init  ***
       !! ** Purpose: Initialise ice related fields for NEMO and coupling
       !!
+      !!---------------------------------------------------------------------
       INTEGER, INTENT( in  ) ::   ksbc                ! surface forcing type
       REAL(wp), DIMENSION(:,:), POINTER :: ztmp1, ztmp2
@@ -289,5 +293,5 @@
 
    
-   SUBROUTINE cice_sbc_in (kt, ksbc)
+   SUBROUTINE cice_sbc_in( kt, ksbc )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE cice_sbc_in  ***
@@ -296,5 +300,5 @@
       INTEGER, INTENT(in   ) ::   kt   ! ocean time step
       INTEGER, INTENT(in   ) ::   ksbc ! surface forcing type
-
+      !
       INTEGER  ::   ji, jj, jl                   ! dummy loop indices      
       REAL(wp), DIMENSION(:,:), POINTER :: ztmp, zpice
@@ -490,21 +494,19 @@
 ! x comp and y comp of sea surface slope (on F points)
 ! T point to F point
-      DO jj=1,jpjm1
-         DO ji=1,jpim1
-            ztmp(ji,jj)=0.5 * (  (zpice(ji+1,jj  )-zpice(ji,jj  ))/e1u(ji,jj  )   &
-                               + (zpice(ji+1,jj+1)-zpice(ji,jj+1))/e1u(ji,jj+1) ) & 
-                            *  fmask(ji,jj,1)
-         ENDDO
-      ENDDO
-      CALL nemo2cice(ztmp,ss_tltx,'F', -1. )
+      DO jj = 1, jpjm1
+         DO ji = 1, jpim1
+            ztmp(ji,jj)=0.5 * (  (zpice(ji+1,jj  )-zpice(ji,jj  )) * r1_e1u(ji,jj  )    &
+               &               + (zpice(ji+1,jj+1)-zpice(ji,jj+1)) * r1_e1u(ji,jj+1)  ) * fmask(ji,jj,1)
+         END DO
+      END DO
+      CALL nemo2cice( ztmp,ss_tltx,'F', -1. )
 
 ! T point to F point
-      DO jj=1,jpjm1
-         DO ji=1,jpim1
-            ztmp(ji,jj)=0.5 * (  (zpice(ji  ,jj+1)-zpice(ji  ,jj))/e2v(ji  ,jj)   &
-                               + (zpice(ji+1,jj+1)-zpice(ji+1,jj))/e2v(ji+1,jj) ) &
-                            *  fmask(ji,jj,1)
-         ENDDO
-      ENDDO
+      DO jj = 1, jpjm1
+         DO ji = 1, jpim1
+            ztmp(ji,jj)=0.5 * (  (zpice(ji  ,jj+1)-zpice(ji  ,jj)) * r1_e2v(ji  ,jj)    &
+               &               + (zpice(ji+1,jj+1)-zpice(ji+1,jj)) * r1_e2v(ji+1,jj)  ) *  fmask(ji,jj,1)
+         END DO
+      END DO
       CALL nemo2cice(ztmp,ss_tlty,'F', -1. )
 
@@ -517,5 +519,5 @@
 
 
-   SUBROUTINE cice_sbc_out (kt,ksbc)
+   SUBROUTINE cice_sbc_out( kt, ksbc )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE cice_sbc_out  ***
@@ -575,5 +577,5 @@
 ! Update taum with modulus of ice-ocean stress 
 ! strocnxT and strocnyT are not weighted by ice fraction in CICE so must be done here 
-taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1**2. + ztmp2**2.) 
+taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2) 
 
 ! Freshwater fluxes 
@@ -888,5 +890,4 @@
 #endif
       !!---------------------------------------------------------------------
-
       CHARACTER(len=1), INTENT( in ) ::   &
           cd_type       ! nature of pn grid-point
@@ -908,4 +909,5 @@
 
       INTEGER  ::   ji, jj, jn                      ! dummy loop indices
+      !!---------------------------------------------------------------------
 
 !     A. Ensure all haloes are filled in NEMO field (pn)
@@ -1096,5 +1098,4 @@
    !!   Default option           Dummy module         NO CICE sea-ice model
    !!----------------------------------------------------------------------
-   !! $Id$
 CONTAINS
 
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 5891)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 5901)
@@ -18,19 +18,18 @@
    USE eosbn2          ! equation of state
    USE sbc_oce         ! surface boundary condition: ocean fields
+   USE zdfbfr          !
+   !
+   USE in_out_manager  ! I/O manager
+   USE iom             ! I/O manager library
+   USE fldread         ! read input field at current time step
    USE lbclnk          !
-   USE iom             ! I/O manager library
-   USE in_out_manager  ! I/O manager
    USE wrk_nemo        ! Memory allocation
    USE timing          ! Timing
    USE lib_fortran     ! glob_sum
-   USE zdfbfr
-   USE fldread         ! read input field at current time step
-
-
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   sbc_isf, sbc_isf_div, sbc_isf_alloc  ! routine called in sbcmod and divcur
+   PUBLIC   sbc_isf, sbc_isf_div, sbc_isf_alloc  ! routine called in sbcmod and divhor
 
    ! public in order to be able to output then 
@@ -53,12 +52,5 @@
    REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  risfLeff               !:effective length (Leff) BG03 nn_isf==2
    REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  ttbl, stbl, utbl, vtbl !:top boundary layer variable at T point
-#if defined key_agrif
-   ! AGRIF can not handle these arrays as integers. The reason is a mystery but problems avoided by declaring them as reals
-   REAL(wp),    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  misfkt, misfkb         !:Level of ice shelf base
-                                                                                          !: (first wet level and last level include in the tbl)
-#else
    INTEGER,    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  misfkt, misfkb         !:Level of ice shelf base
-#endif
-
 
    REAL(wp), PUBLIC, SAVE ::   rcpi   = 2000.0_wp     ! phycst ?
@@ -79,26 +71,29 @@
 #  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)
+   !! NEMO/OPA 3.7 , LOCEAN-IPSL (2015)
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
  
-  SUBROUTINE sbc_isf(kt)
-    INTEGER, INTENT(in)          ::   kt         ! ocean time step
-    INTEGER                      ::   ji, jj, jk, ijkmin, inum, ierror
-    INTEGER                      ::   ikt, ikb   ! top and bottom level of the isf boundary layer
-    REAL(wp)                     ::   rmin
-    REAL(wp)                     ::   zhk
-    CHARACTER(len=256)           ::   cfisf, cvarzisf, cvarhisf   ! name for isf file
-    CHARACTER(LEN=256)           :: cnameis                     ! name of iceshelf file
-    CHARACTER (LEN=32)           :: cvarLeff                    ! variable name for efficient Length scale
-    INTEGER           ::   ios           ! Local integer output status for namelist read
-      !
+   SUBROUTINE sbc_isf(kt)
       !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE sbc_isf  ***
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in)          ::   kt         ! ocean time step
+      !
+      INTEGER                      ::   ji, jj, jk, ijkmin, inum, ierror
+      INTEGER                      ::   ikt, ikb   ! top and bottom level of the isf boundary layer
+      REAL(wp)                     ::   rmin
+      REAL(wp)                     ::   zhk
+      REAL(wp)                     ::   zt_frz, zpress
+      CHARACTER(len=256)           ::   cfisf , cvarzisf, cvarhisf   ! name for isf file
+      CHARACTER(LEN=256)           :: cnameis                     ! name of iceshelf file
+      CHARACTER (LEN=32)           :: cvarLeff                    ! variable name for efficient Length scale
+      INTEGER           ::   ios           ! Local integer output status for namelist read
+      !!
       NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, ln_divisf, ln_conserve, rn_gammat0, rn_gammas0, nn_gammablk, &
-                         & sn_fwfisf, sn_qisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf
-      !
+         &                 sn_fwfisf, sn_qisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf
+      !!---------------------------------------------------------------------
       !
       !                                         ! ====================== !
@@ -113,5 +108,4 @@
 902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp )
          IF(lwm) WRITE ( numond, namsbc_isf )
-
 
          IF ( lwp ) WRITE(numout,*)
@@ -194,7 +188,6 @@
          END IF
          
+         ! compute bottom level of isf tbl and thickness of tbl below the ice shelf
          rhisf_tbl_0(:,:) = rhisf_tbl(:,:)
-
-         ! compute bottom level of isf tbl and thickness of tbl below the ice shelf
          DO jj = 1,jpj
             DO ji = 1,jpi
@@ -217,5 +210,5 @@
             END DO
          END DO
-         
+         !
       END IF
 
@@ -270,10 +263,19 @@
          END IF
          ! compute tsc due to isf
-         ! WARNING water add at temp = 0C, correction term is added in trasbc, maybe better here but need a 3D variable).
-         risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp !
+         ! WARNING water add at temp = 0C, correction term is added, maybe better here but need a 3D variable).
+!         zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04
+         zt_frz = -1.9 !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress )
+         risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - rdivisf * fwfisf(:,:) * zt_frz * r1_rau0 !
          
          ! salt effect already take into account in vertical advection
          risf_tsc(:,:,jp_sal) = (1.0_wp-rdivisf) * fwfisf(:,:) * stbl(:,:) * r1_rau0
-          
+
+         ! output
+         IF( iom_use('qisf'  ) )   CALL iom_put('qisf'  , qisf)
+         IF( iom_use('fwfisf') )   CALL iom_put('fwfisf', fwfisf * stbl(:,:) / soce )
+
+         ! if apply only on the trend and not as a volume flux (rdivisf = 0), fwfisf have to be set to 0 now
+         fwfisf(:,:) = rdivisf * fwfisf(:,:)         
+ 
          ! lbclnk
          CALL lbc_lnk(risf_tsc(:,:,jp_tem),'T',1.)
@@ -295,10 +297,8 @@
          ENDIF
          ! 
-         ! output
-         CALL iom_put('qisf'  , qisf)
-         IF( iom_use('fwfisf') )   CALL iom_put('fwfisf', fwfisf * stbl(:,:) / soce )
       END IF
-  
+      !  
   END SUBROUTINE sbc_isf
+
 
   INTEGER FUNCTION sbc_isf_alloc()
@@ -321,26 +321,26 @@
   END FUNCTION
 
-  SUBROUTINE sbc_isf_bg03(kt)
-   !!==========================================================================
-   !!                 *** SUBROUTINE sbcisf_bg03  ***
-   !! add net heat and fresh water flux from ice shelf melting
-   !! into the adjacent ocean using the parameterisation by
-   !! Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean
-   !!     interaction for climate models", Ocean Modelling 5(2003) 157-170.
-   !!  (hereafter BG)
-   !!==========================================================================
-   !!----------------------------------------------------------------------
-   !!   sbc_isf_bg03      : routine called from sbcmod
-   !!----------------------------------------------------------------------
-   !!
-   !! ** Purpose   :   Add heat and fresh water fluxes due to ice shelf melting
-   !! ** Reference :   Beckmann et Goosse, 2003, Ocean Modelling
-   !!
-   !! History :
-   !!      !  06-02  (C. Wang) Original code
-   !!----------------------------------------------------------------------
-
-    INTEGER, INTENT ( in ) :: kt
-
+
+   SUBROUTINE sbc_isf_bg03(kt)
+      !!==========================================================================
+      !!                 *** SUBROUTINE sbcisf_bg03  ***
+      !! add net heat and fresh water flux from ice shelf melting
+      !! into the adjacent ocean using the parameterisation by
+      !! Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean
+      !!     interaction for climate models", Ocean Modelling 5(2003) 157-170.
+      !!  (hereafter BG)
+      !!==========================================================================
+      !!----------------------------------------------------------------------
+      !!   sbc_isf_bg03      : routine called from sbcmod
+      !!----------------------------------------------------------------------
+      !!
+      !! ** Purpose   :   Add heat and fresh water fluxes due to ice shelf melting
+      !! ** Reference :   Beckmann et Goosse, 2003, Ocean Modelling
+      !!
+      !! History :
+      !!      !  06-02  (C. Wang) Original code
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT ( in ) :: kt
+      !
     INTEGER :: ji, jj, jk, jish  !temporary integer
     INTEGER :: ijkmin
@@ -386,9 +386,11 @@
              qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp
           END IF
-       ENDDO
-    ENDDO
+       END DO
+    END DO
     !
     IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_bg03')
+      !
   END SUBROUTINE sbc_isf_bg03
+
 
    SUBROUTINE sbc_isf_cav( kt )
@@ -439,5 +441,4 @@
       !
       !
-!CDIR COLLAPSE
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -472,9 +473,6 @@
 
                      nit = nit + 1
-                     IF (nit .GE. 100) THEN
-                        !WRITE(numout,*) "sbcisf : too many iteration ... ", zhtflx, zhtflx_b,zgammat, rn_gammat0, rn_tfri2, nn_gammablk, ji,jj
-                        !WRITE(numout,*) "sbcisf : too many iteration ... ", (zhtflx - zhtflx_b)/zhtflx
-                        CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
-                     END IF
+                     IF (nit .GE. 100) CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
+
 ! save gammat and compute zhtflx_b
                      zgammat2d(ji,jj)=zgammat
@@ -496,5 +494,4 @@
 
 ! More complicated 3 equation thermodynamics as in MITgcm
-!CDIR COLLAPSE
          DO jj = 2, jpj
             DO ji = 2, jpi
@@ -565,6 +562,7 @@
       !
       IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_cav')
-
+      !
    END SUBROUTINE sbc_isf_cav
+
 
    SUBROUTINE sbc_isf_gammats(gt, gs, zqhisf, zqwisf, ji, jj, lit )
@@ -693,6 +691,7 @@
                END IF
       END IF
-
+      !
    END SUBROUTINE
+
 
    SUBROUTINE sbc_isf_tbl( varin, varout, cptin )
@@ -756,5 +755,5 @@
       IF (cptin == 'T') CALL lbc_lnk(varout,'T',1.)
       IF (cptin == 'U' .OR. cptin == 'V') CALL lbc_lnk(varout,'T',-1.)
-
+      !
    END SUBROUTINE sbc_isf_tbl
       
@@ -794,5 +793,5 @@
                ! test on tmask useless ?????
                DO jk = ikt, mbkt(ji,jj)
-!                  IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
+                  IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
                END DO
                rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
@@ -823,5 +822,6 @@
       !
    END SUBROUTINE sbc_isf_div
-                        
+
+
    FUNCTION tinsitu( ptem, psal, ppress ) RESULT( pti )
       !!----------------------------------------------------------------------
@@ -874,5 +874,6 @@
       !
    END FUNCTION tinsitu
-   !
+
+
    FUNCTION fsatg( pfps, pfpt, pfphp )
       !!----------------------------------------------------------------------
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 5891)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 5901)
@@ -28,5 +28,4 @@
    USE sbcdcy           ! surface boundary condition: diurnal cycle
    USE sbcssm           ! surface boundary condition: sea-surface mean variables
-   USE sbcapr           ! surface boundary condition: atmospheric pressure
    USE sbcana           ! surface boundary condition: analytical formulation
    USE sbcflx           ! surface boundary condition: flux formulation
@@ -133,5 +132,5 @@
          WRITE(numout,*) '              forced-coupled mixed formulation           ln_mixcpl   = ', ln_mixcpl
          WRITE(numout,*) '              OASIS coupling (with atm or sas)           lk_oasis    = ', lk_oasis
-         WRITE(numout,*) '              components of your executable              nn_components = ', nn_components
+         WRITE(numout,*) '              components of your executable            nn_components = ', nn_components
          WRITE(numout,*) '              Multicategory heat flux formulation (LIM3) nn_limflx   = ', nn_limflx
          WRITE(numout,*) '           Misc. options of sbc : '
@@ -176,11 +175,11 @@
 
       !                              ! allocate sbc arrays
-      IF( sbc_oce_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_oce arrays' )
+      IF( sbc_oce_alloc() /= 0 )   CALL ctl_stop( 'sbc_init : unable to allocate sbc_oce arrays' )
 
       !                          ! Checks:
       IF( nn_isf .EQ. 0 ) THEN                      ! variable initialisation if no ice shelf 
-         IF( sbc_isf_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_isf arrays' )
-         fwfisf  (:,:) = 0.0_wp
-         fwfisf_b(:,:) = 0.0_wp
+         IF( sbc_isf_alloc() /= 0 )   CALL ctl_stop( 'sbc_init : unable to allocate sbc_isf arrays' )
+         fwfisf  (:,:)   = 0.0_wp ; fwfisf_b  (:,:)   = 0.0_wp
+         risf_tsc(:,:,:) = 0.0_wp ; risf_tsc_b(:,:,:) = 0.0_wp
          rdivisf       = 0.0_wp
       END IF
@@ -224,7 +223,7 @@
          ENDIF
       ELSE
-      IF ( ln_cdgw .OR. ln_sdw  )                                         & 
-         &   CALL ctl_stop('Not Activated Wave Module (ln_wave=F) but     &
-         & asked coupling with drag coefficient (ln_cdgw =T) or Stokes drift (ln_sdw=T) ')
+      IF ( ln_cdgw .OR. ln_sdw  )                                                           & 
+         &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    &
+         &                  'with drag coefficient (ln_cdgw =T) or Stokes drift (ln_sdw=T) ')
       ENDIF 
       !                          ! Choice of the Surface Boudary Condition (set nsbc)
@@ -241,16 +240,9 @@
       IF( nn_components == jp_iam_opa )   &
          &                  THEN   ;   nsbc = jp_none    ; icpt = icpt + 1   ;   ENDIF       ! opa coupling via SAS module
-      IF( lk_esopa        )            nsbc = jp_esopa                                       ! esopa test, ALL formulations
-      !
-      IF( icpt /= 1 .AND. .NOT.lk_esopa ) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) '           E R R O R in setting the sbc, one and only one namelist/CPP key option '
-         WRITE(numout,*) '                     must be choosen. You choose ', icpt, ' option(s)'
-         WRITE(numout,*) '                     We stop'
-         nstop = nstop + 1
-      ENDIF
+      !
+      IF( icpt /= 1 )    CALL ctl_stop( 'sbc_init: choose ONE and only ONE sbc option' )
+      !
       IF(lwp) THEN
          WRITE(numout,*)
-         IF( nsbc == jp_esopa   )   WRITE(numout,*) '              ESOPA test All surface boundary conditions'
          IF( nsbc == jp_gyre    )   WRITE(numout,*) '              GYRE analytical formulation'
          IF( nsbc == jp_ana     )   WRITE(numout,*) '              analytical formulation'
@@ -267,10 +259,9 @@
       !
       IF( lk_oasis )   CALL sbc_cpl_init (nn_ice)   ! OASIS initialisation. must be done before: (1) first time step
-      !                                                     !                                            (2) the use of nn_fsbc
+      !                                             !                                            (2) the use of nn_fsbc
 
 !     nn_fsbc initialization if OPA-SAS coupling via OASIS
 !     sas model time step has to be declared in OASIS (mandatory) -> nn_fsbc has to be modified accordingly
       IF ( nn_components /= jp_iam_nemo ) THEN
-
          IF ( nn_components == jp_iam_opa ) nn_fsbc = cpl_freq('O_SFLX') / NINT(rdt)
          IF ( nn_components == jp_iam_sas ) nn_fsbc = cpl_freq('I_SFLX') / NINT(rdt)
@@ -344,9 +335,4 @@
       !                                            !        forcing field computation         !
       !                                            ! ---------------------------------------- !
-      !
-      IF ( .NOT. lk_bdy ) then
-         IF( ln_apr_dyn ) CALL sbc_apr( kt )                ! atmospheric pressure provided at kt+0.5*nn_fsbc
-      ENDIF
-                                                         ! (caution called before sbc_ssm)
       !
       IF( nn_components /= jp_iam_sas )   CALL sbc_ssm( kt )   ! ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m)
@@ -373,11 +359,4 @@
          IF( nn_components == jp_iam_opa ) &
                              CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! OPA-SAS coupling: OPA receiving fields from SAS
-      CASE( jp_esopa )                                
-                             CALL sbc_ana     ( kt )                    ! ESOPA, test ALL the formulations
-                             CALL sbc_gyre    ( kt )                    !
-                             CALL sbc_flx     ( kt )                    !
-                             CALL sbc_blk_clio( kt )                    !
-                             CALL sbc_blk_core( kt )                    !
-                             CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   !
       END SELECT
 
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 5891)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 5901)
@@ -31,8 +31,9 @@
    PRIVATE
 
-   PUBLIC   sbc_rnf       ! routine call in sbcmod module
-   PUBLIC   sbc_rnf_div   ! routine called in divcurl module
-   PUBLIC   sbc_rnf_alloc ! routine call in sbcmod module
-   PUBLIC   sbc_rnf_init  ! (PUBLIC for TAM)
+   PUBLIC   sbc_rnf       ! routine called in sbcmod module
+   PUBLIC   sbc_rnf_div   ! routine called in divhor module
+   PUBLIC   sbc_rnf_alloc ! routine called in sbcmod module
+   PUBLIC   sbc_rnf_init  ! routine called in sbcmod module
+   
    !                                                     !!* namsbc_rnf namelist *
    CHARACTER(len=100)         ::   cn_dir          !: Root directory for location of rnf files
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90	(revision 5891)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90	(revision 5901)
@@ -107,5 +107,4 @@
             IF( nn_sssr == 1 ) THEN                                   !* Salinity damping term (salt flux only (sfx))
                zsrp = rn_deds / rday                                  ! from [mm/day] to [kg/m2/s]
-!CDIR COLLAPSE
                DO jj = 1, jpj
                   DO ji = 1, jpi
@@ -121,5 +120,4 @@
                zsrp = rn_deds / rday                                  ! from [mm/day] to [kg/m2/s]
                zerp_bnd = rn_sssr_bnd / rday                          !       -              -    
-!CDIR COLLAPSE
                DO jj = 1, jpj
                   DO ji = 1, jpi                            
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 5891)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 5901)
@@ -4,19 +4,21 @@
    !! Wave module 
    !!======================================================================
-   !! History :  3.3.1  !   2011-09  (Adani M)  Original code: Drag Coefficient 
-   !!         :  3.4    !   2012-10  (Adani M)                 Stokes Drift 
+   !! History :  3.3  !   2011-09  (Adani M)  Original code: Drag Coefficient 
+   !!         :  3.4  !   2012-10  (Adani M)                 Stokes Drift 
    !!----------------------------------------------------------------------
-   USE iom             ! I/O manager library
-   USE in_out_manager  ! I/O manager
-   USE lib_mpp         ! distribued memory computing library
+
+   !!----------------------------------------------------------------------
+   !!   sbc_wave      : read drag coefficient from wave model in netcdf files 
+   !!----------------------------------------------------------------------
+   USE oce            ! 
+   USE sbc_oce	       ! Surface boundary condition: ocean fields
+   USE bdy_oce        !
+   USE domvvl         !
+   !
+   USE iom            ! I/O manager library
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! distribued memory computing library
    USE fldread	       ! read input fields
-   USE oce
-   USE sbc_oce	       ! Surface boundary condition: ocean fields
-   USE domvvl
-
-   
-   !!----------------------------------------------------------------------
-   !!   sbc_wave       : read drag coefficient from wave model in netcdf files 
-   !!----------------------------------------------------------------------
+   USE wrk_nemo       !
 
    IMPLICIT NONE
@@ -25,18 +27,21 @@
    PUBLIC   sbc_wave    ! routine called in sbc_blk_core or sbc_blk_mfs
    
-   INTEGER , PARAMETER ::   jpfld  = 3           ! maximum number of files to read for srokes drift
-   INTEGER , PARAMETER ::   jp_usd = 1           ! index of stokes drift  (i-component) (m/s)    at T-point
-   INTEGER , PARAMETER ::   jp_vsd = 2           ! index of stokes drift  (j-component) (m/s)    at T-point
-   INTEGER , PARAMETER ::   jp_wn  = 3           ! index of wave number                 (1/m)    at T-point
+   INTEGER , PARAMETER ::   jpfld  = 3   ! maximum number of files to read for srokes drift
+   INTEGER , PARAMETER ::   jp_usd = 1   ! index of stokes drift  (i-component) (m/s)    at T-point
+   INTEGER , PARAMETER ::   jp_vsd = 2   ! index of stokes drift  (j-component) (m/s)    at T-point
+   INTEGER , PARAMETER ::   jp_wn  = 3   ! index of wave number                 (1/m)    at T-point
+
    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_cd	  ! structure of input fields (file informations, fields read) Drag Coefficient
    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_sd	  ! structure of input fields (file informations, fields read) Stokes Drift
-   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: cdn_wave 
-   REAL(wp),ALLOCATABLE,DIMENSION (:,:)              :: usd2d,vsd2d,uwavenum,vwavenum 
-   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:,:)     :: usd3d,vsd3d,wsd3d 
+
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION (:,:)   :: cdn_wave 
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION (:,:,:) :: usd3d, vsd3d, wsd3d 
+   REAL(wp),         ALLOCATABLE, DIMENSION (:,:)   :: usd2d, vsd2d, uwavenum, vwavenum 
 
    !! * Substitutions
 #  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
+   !! NEMO/OPA 3.7 , NEMO Consortium (2014) 
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -56,27 +61,17 @@
       !!              - Compute 3d stokes drift using monochromatic
       !! ** action  :   
-      !!               
       !!---------------------------------------------------------------------
-      USE oce,  ONLY : un,vn,hdivn,rotn
-      USE divcur
-      USE wrk_nemo
-#if defined key_bdy
-      USE bdy_oce, ONLY : bdytmask
-#endif
-      INTEGER, INTENT( in  ) ::  kt       ! ocean time step
-      INTEGER                ::  ierror   ! return error code
-      INTEGER                ::  ifpr, jj,ji,jk 
+      INTEGER, INTENT( in  ) ::   kt       ! ocean time step
+      !
+      INTEGER                ::   ierror   ! return error code
+      INTEGER                ::   ifpr, jj,ji,jk 
       INTEGER                ::   ios     ! Local integer output status for namelist read
-      REAL(wp),DIMENSION(:,:,:),POINTER             ::  udummy,vdummy,hdivdummy,rotdummy
-      REAL                                          ::  z2dt,z1_2dt
       TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i     ! array of namelist informations on the fields to read
       CHARACTER(len=100)     ::  cn_dir                          ! Root directory for location of drag coefficient files
       TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd, sn_wn   ! informations about the fields to be read
-      !!---------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:), POINTER ::   zusd_t, zvsd_t, ze3hdiv   ! 3D workspace
+      !!
       NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_wn
       !!---------------------------------------------------------------------
-
-      !!----------------------------------------------------------------------
-      !
       !
       !                                         ! -------------------- !
@@ -92,5 +87,4 @@
          IF(lwm) WRITE ( numond, namsbc_wave )
          !
-
          IF ( ln_cdgw ) THEN
             ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg
@@ -102,5 +96,5 @@
             ALLOCATE( cdn_wave(jpi,jpj) )
             cdn_wave(:,:) = 0.0
-        ENDIF
+         ENDIF
          IF ( ln_sdw ) THEN
             slf_i(jp_usd) = sn_usd ; slf_i(jp_vsd) = sn_vsd; slf_i(jp_wn) = sn_wn
@@ -113,80 +107,72 @@
             END DO
             CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' )
-            ALLOCATE( usd2d(jpi,jpj),vsd2d(jpi,jpj),uwavenum(jpi,jpj),vwavenum(jpi,jpj) )
+            ALLOCATE( usd2d(jpi,jpj) , vsd2d(jpi,jpj) , uwavenum(jpi,jpj) , vwavenum(jpi,jpj) )
             ALLOCATE( usd3d(jpi,jpj,jpk),vsd3d(jpi,jpj,jpk),wsd3d(jpi,jpj,jpk) )
-            usd2d(:,:) = 0.0 ;  vsd2d(:,:) = 0.0 ; uwavenum(:,:) = 0.0 ; vwavenum(:,:) = 0.0
-            usd3d(:,:,:) = 0.0 ;vsd3d(:,:,:) = 0.0 ; wsd3d(:,:,:) = 0.0
+            usd3d(:,:,:) = 0._wp   ;   usd2d(:,:) = 0._wp   ;    uwavenum(:,:) = 0._wp
+            vsd3d(:,:,:) = 0._wp   ;   vsd2d(:,:) = 0._wp   ;    vwavenum(:,:) = 0._wp
+            wsd3d(:,:,:) = 0._wp
          ENDIF
       ENDIF
+      !
+      IF ( ln_cdgw ) THEN              !==  Neutral drag coefficient  ==!
+         CALL fld_read( kt, nn_fsbc, sf_cd )      ! read from external forcing
+         cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1)
+      ENDIF
+      !
+      IF ( ln_sdw )  THEN              !==  Computation of the 3d Stokes Drift  ==!
+         !
+         CALL fld_read( kt, nn_fsbc, sf_sd )    !* read drag coefficient from external forcing
          !
          !
-      IF ( ln_cdgw ) THEN
-         CALL fld_read( kt, nn_fsbc, sf_cd )      !* read drag coefficient from external forcing
-         cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1)
-      ENDIF
-      IF ( ln_sdw )  THEN
-          CALL fld_read( kt, nn_fsbc, sf_sd )      !* read drag coefficient from external forcing
-
-         ! Interpolate wavenumber, stokes drift into the grid_V and grid_V
-         !-------------------------------------------------
-
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               uwavenum(ji,jj)=0.5 * ( 2. - umask(ji,jj,1) ) * ( sf_sd(3)%fnow(ji,jj,1) * tmask(ji,jj,1) &
-               &                                + sf_sd(3)%fnow(ji+1,jj,1) * tmask(ji+1,jj,1) )
-
-               vwavenum(ji,jj)=0.5 * ( 2. - vmask(ji,jj,1) ) * ( sf_sd(3)%fnow(ji,jj,1) * tmask(ji,jj,1) &
-               &                                + sf_sd(3)%fnow(ji,jj+1,1) * tmask(ji,jj+1,1) )
-
-               usd2d(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( sf_sd(1)%fnow(ji,jj,1) * tmask(ji,jj,1) &
-               &                                + sf_sd(1)%fnow(ji+1,jj,1) * tmask(ji+1,jj,1) )
-
-               vsd2d(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( sf_sd(2)%fnow(ji,jj,1) * tmask(ji,jj,1) &
-               &                                + sf_sd(2)%fnow(ji,jj+1,1) * tmask(ji,jj+1,1) )
+         CALL wrk_alloc( jpi,jpj,jpk,   zusd_t, zvsd_t, ze3hdiv )
+         !                                      !* distribute it on the vertical
+         DO jk = 1, jpkm1
+            zusd_t(:,:,jk) = sf_sd(jp_usd)%fnow(:,:,1) * EXP( -2._wp * sf_sd(jp_wn)%fnow(:,:,1) * fsdept_n(:,:,jk) )
+            zvsd_t(:,:,jk) = sf_sd(jp_vsd)%fnow(:,:,1) * EXP( -2._wp * sf_sd(jp_wn)%fnow(:,:,1) * fsdept_n(:,:,jk) )
+         END DO
+         !                                      !* interpolate the stokes drift from t-point to u- and v-points
+         DO jk = 1, jpkm1
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1
+                   usd3d(ji,jj,jk) = 0.5_wp * ( zusd_t(ji  ,jj,jk) + zusd_t(ji+1,jj,jk) ) * umask(ji,jj,jk)
+                   vsd3d(ji,jj,jk) = 0.5_wp * ( zvsd_t(ji  ,jj,jk) + zvsd_t(ji,jj+1,jk) ) * vmask(ji,jj,jk)
+               END DO
             END DO
          END DO
-
-          !Computation of the 3d Stokes Drift
-          DO jk = 1, jpk
-             DO jj = 1, jpj-1
-                DO ji = 1, jpi-1
-                   usd3d(ji,jj,jk) = usd2d(ji,jj)*exp(2.0*uwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj  ,jk))))
-                   vsd3d(ji,jj,jk) = vsd2d(ji,jj)*exp(2.0*vwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji  ,jj+1,jk))))
-                END DO
-             END DO
-             usd3d(jpi,:,jk) = usd2d(jpi,:)*exp( 2.0*uwavenum(jpi,:)*(-gdept_0(jpi,:,jk)) )
-             vsd3d(:,jpj,jk) = vsd2d(:,jpj)*exp( 2.0*vwavenum(:,jpj)*(-gdept_0(:,jpj,jk)) )
-          END DO
-
-          CALL wrk_alloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy)
-          
-          udummy(:,:,:)=un(:,:,:)
-          vdummy(:,:,:)=vn(:,:,:)
-          hdivdummy(:,:,:)=hdivn(:,:,:)
-          rotdummy(:,:,:)=rotn(:,:,:)
-          un(:,:,:)=usd3d(:,:,:)
-          vn(:,:,:)=vsd3d(:,:,:)
-          CALL div_cur(kt)
-      !                                           !------------------------------!
-      !                                           !     Now Vertical Velocity    !
-      !                                           !------------------------------!
-          z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
-
-          z1_2dt = 1.e0 / z2dt
-          DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
-             ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
-             wsd3d(:,:,jk) = wsd3d(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
-                &                      - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
-                &                         * tmask(:,:,jk) * z1_2dt
+         CALL lbc_lnk( usd3d(:,:,:), 'U', -1. )
+         CALL lbc_lnk( vsd3d(:,:,:), 'V', -1. )
+         !
+         DO jk = 1, jpkm1                       !* e3t * Horizontal divergence  ==!
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ze3hdiv(ji,jj,jk) = (  e2u(ji  ,jj) * fse3u_n(ji  ,jj,jk) * usd3d(ji  ,jj,jk)     &
+                     &                 - e2u(ji-1,jj) * fse3u_n(ji-1,jj,jk) * usd3d(ji-1,jj,jk)     &
+                     &                 + e1v(ji,jj  ) * fse3v_n(ji,jj  ,jk) * vsd3d(ji,jj  ,jk)     &
+                     &                 - e1v(ji,jj-1) * fse3v_n(ji,jj-1,jk) * vsd3d(ji,jj-1,jk)   ) * r1_e1e2t(ji,jj)
+               END DO  
+            END DO  
+            IF( .NOT. AGRIF_Root() ) THEN
+               IF( nbondi ==  1 .OR. nbondi == 2 )   ze3hdiv(nlci-1,   :  ,jk) = 0._wp      ! east
+               IF( nbondi == -1 .OR. nbondi == 2 )   ze3hdiv(  2   ,   :  ,jk) = 0._wp      ! west
+               IF( nbondj ==  1 .OR. nbondj == 2 )   ze3hdiv(  :   ,nlcj-1,jk) = 0._wp      ! north
+               IF( nbondj == -1 .OR. nbondj == 2 )   ze3hdiv(  :   ,  2   ,jk) = 0._wp      ! south
+            ENDIF
+         END DO
+         CALL lbc_lnk( ze3hdiv, 'T', 1. ) 
+         !
+         DO jk = jpkm1, 1, -1                   !* integrate from the bottom the e3t * hor. divergence
+            wsd3d(:,:,jk) = wsd3d(:,:,jk+1) - ze3hdiv(:,:,jk)
+         END DO
 #if defined key_bdy
-             wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:)
+         IF( lk_bdy ) THEN
+            DO jk = 1, jpkm1
+               wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:)
+            END DO
+         ENDIF
 #endif
-          END DO
-          hdivn(:,:,:)=hdivdummy(:,:,:)
-          rotn(:,:,:)=rotdummy(:,:,:)
-          vn(:,:,:)=vdummy(:,:,:)
-          un(:,:,:)=udummy(:,:,:)
-          CALL wrk_dealloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy)
+         CALL wrk_dealloc( jpi,jpj,jpk,   zusd_t, zvsd_t, ze3hdiv )
+         ! 
       ENDIF
+      !
    END SUBROUTINE sbc_wave
       
