Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90	(revision 6121)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90	(revision 6140)
@@ -4,9 +4,7 @@
    !! Ocean forcing:  read input field for surface boundary condition
    !!=====================================================================
-   !! History :  2.0  !  06-2006  (S. Masson, G. Madec) Original code
-   !!                 !  05-2008  (S. Alderson) Modified for Interpolation in memory
-   !!                 !                         from input grid to model grid
-   !!                 !  10-2013  (D. Delrosso, P. Oddo) implement suppression of 
-   !!                 !                         land point prior to interpolation
+   !! History :  2.0  !  06-2006  (S. Masson, G. Madec)  Original code
+   !!                 !  05-2008  (S. Alderson)  Modified for Interpolation in memory from input grid to model grid
+   !!                 !  10-2013  (D. Delrosso, P. Oddo)  suppression of land point prior to interpolation
    !!----------------------------------------------------------------------
 
@@ -15,15 +13,16 @@
    !!                   surface boundary condition
    !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers
-   USE dom_oce         ! ocean space and time domain
-   USE phycst          ! ???
-   USE in_out_manager  ! I/O manager
-   USE iom             ! I/O manager library
-   USE geo2ocean       ! for vector rotation on to model grid
-   USE lib_mpp         ! MPP library
-   USE wrk_nemo        ! work arrays
-   USE lbclnk          ! ocean lateral boundary conditions (C1D case)
-   USE ioipsl, ONLY :   ymds2ju, ju2ymds   ! for calendar
-   USE sbc_oce
+   USE oce            ! ocean dynamics and tracers
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constant
+   USE sbc_oce        ! surface boundary conditions : fields
+   USE geo2ocean      ! for vector rotation on to model grid
+   !
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O manager library
+   USE ioipsl  , ONLY : ymds2ju, ju2ymds   ! for calendar
+   USE lib_mpp        ! MPP library
+   USE wrk_nemo       ! work arrays
+   USE lbclnk         ! ocean lateral boundary conditions (C1D case)
    
    IMPLICIT NONE
@@ -60,6 +59,6 @@
       INTEGER , DIMENSION(2)          ::   nrec_b       ! before record (1: index, 2: second since Jan. 1st 00h of nit000 year)
       INTEGER , DIMENSION(2)          ::   nrec_a       ! after  record (1: index, 2: second since Jan. 1st 00h of nit000 year)
-      REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:  ) ::   fnow       ! input fields interpolated to now time step
-      REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:,:) ::   fdta       ! 2 consecutive record of input fields
+      REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:  ) ::   fnow   ! input fields interpolated to now time step
+      REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:,:) ::   fdta   ! 2 consecutive record of input fields
       CHARACTER(len = 256)            ::   wgtname      ! current name of the NetCDF weight file acting as a key
       !                                                 ! into the WGTLIST structure
@@ -133,22 +132,21 @@
       INTEGER  , INTENT(in   ), OPTIONAL     ::   kit       ! subcycle timestep for timesplitting option
       INTEGER  , INTENT(in   ), OPTIONAL     ::   kt_offset ! provide fields at time other than "now"
-                                                            !   kt_offset = -1 => fields at "before" time level
-                                                            !   kt_offset = +1 => fields at "after"  time level
-                                                            !   etc.
-      !!
-      INTEGER  ::   itmp       ! temporary variable
-      INTEGER  ::   imf        ! size of the structure sd
-      INTEGER  ::   jf         ! dummy indices
-      INTEGER  ::   isecend    ! number of second since Jan. 1st 00h of nit000 year at nitend
-      INTEGER  ::   isecsbc    ! number of seconds between Jan. 1st 00h of nit000 year and the middle of sbc time step
-      INTEGER  ::   it_offset  ! local time offset variable
-      LOGICAL  ::   llnxtyr    ! open next year  file?
-      LOGICAL  ::   llnxtmth   ! open next month file?
-      LOGICAL  ::   llstop     ! stop is the file does not exist
+      !                                                     !   kt_offset = -1 => fields at "before" time level
+      !                                                     !   kt_offset = +1 => fields at "after"  time level
+      !                                                     !   etc.
+      INTEGER  ::   itmp         ! local variable
+      INTEGER  ::   imf          ! size of the structure sd
+      INTEGER  ::   jf           ! dummy indices
+      INTEGER  ::   isecend      ! number of second since Jan. 1st 00h of nit000 year at nitend
+      INTEGER  ::   isecsbc      ! number of seconds between Jan. 1st 00h of nit000 year and the middle of sbc time step
+      INTEGER  ::   it_offset    ! local time offset variable
+      LOGICAL  ::   llnxtyr      ! open next year  file?
+      LOGICAL  ::   llnxtmth     ! open next month file?
+      LOGICAL  ::   llstop       ! stop is the file does not exist
       LOGICAL  ::   ll_firstcall ! true if this is the first call to fld_read for this set of fields
-      REAL(wp) ::   ztinta     ! ratio applied to after  records when doing time interpolation
-      REAL(wp) ::   ztintb     ! ratio applied to before records when doing time interpolation
-      CHARACTER(LEN=1000) ::   clfmt   ! write format
-      TYPE(MAP_POINTER) ::   imap   ! global-to-local mapping indices
+      REAL(wp) ::   ztinta       ! ratio applied to after  records when doing time interpolation
+      REAL(wp) ::   ztintb       ! ratio applied to before records when doing time interpolation
+      CHARACTER(LEN=1000) ::   clfmt  ! write format
+      TYPE(MAP_POINTER)   ::   imap   ! global-to-local mapping indices
       !!---------------------------------------------------------------------
       ll_firstcall = kt == nit000
@@ -166,5 +164,5 @@
          isecsbc = nsec_year + nsec1jan000 + (kit+it_offset)*NINT( rdt/REAL(nn_baro,wp) )
       ELSE                      ! middle of sbc time step
-         isecsbc = nsec_year + nsec1jan000 + NINT(0.5 * REAL(kn_fsbc - 1,wp) * rdttra(1)) + it_offset * NINT(rdttra(1))
+         isecsbc = nsec_year + nsec1jan000 + NINT(0.5 * REAL(kn_fsbc - 1,wp) * rdt) + it_offset * NINT(rdt)
       ENDIF
       imf = SIZE( sd )
@@ -193,5 +191,5 @@
                CALL fld_rec( kn_fsbc, sd(jf), kt_offset = it_offset, kit = kit )    ! update after record informations
 
-               ! if kn_fsbc*rdttra is larger than nfreqh (which is kind of odd),
+               ! if kn_fsbc*rdt is larger than nfreqh (which is kind of odd),
                ! it is possible that the before value is no more the good one... we have to re-read it
                ! if before is not the last record of the file currently opened and after is the first record to be read
@@ -214,5 +212,5 @@
                IF( sd(jf)%ln_tint ) THEN
                   
-                  ! if kn_fsbc*rdttra is larger than nfreqh (which is kind of odd),
+                  ! if kn_fsbc*rdt is larger than nfreqh (which is kind of odd),
                   ! it is possible that the before value is no more the good one... we have to re-read it
                   ! if before record is not just just before the after record...
@@ -245,5 +243,5 @@
                         ! year/month/week/day file to be not present. If the run continue further than the current
                         ! year/month/week/day, next year/month/week/day file must exist
-                        isecend = nsec_year + nsec1jan000 + (nitend - kt) * NINT(rdttra(1))   ! second at the end of the run 
+                        isecend = nsec_year + nsec1jan000 + (nitend - kt) * NINT(rdt)   ! second at the end of the run
                         llstop = isecend > sd(jf)%nrec_a(2)                                   ! read more than 1 record of next year
                         ! we suppose that the date of next file is next day (should be ok even for weekly files...)
@@ -299,7 +297,5 @@
          END DO                                    ! --- end loop over field --- !
          !
-         !                                         ! ====================================== !
-      ENDIF                                        ! update field at each kn_fsbc time-step !
-      !                                            ! ====================================== !
+      ENDIF
       !
    END SUBROUTINE fld_read
@@ -333,12 +329,12 @@
       llprevday  = .FALSE.
       isec_week  = 0
-            
+      !
       ! define record informations
       CALL fld_rec( kn_fsbc, sdjf, ldbefore = .TRUE. )  ! return before values in sdjf%nrec_a (as we will swap it later)
-
+      !
       ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar 
-
+      !
       IF( sdjf%ln_tint ) THEN ! we need to read the previous record and we will put it in the current record structure
-
+         !
          IF( sdjf%nrec_a(1) == 0  ) THEN   ! we redefine record sdjf%nrec_a(1) with the last record of previous year file
             IF    ( sdjf%nfreqh == -12 ) THEN   ! yearly mean
@@ -391,5 +387,5 @@
          !
          CALL fld_clopn( sdjf, iyear, imonth, iday, .NOT. llprev )
-
+         !
          ! if previous year/month/day file does not exist, we switch to the current year/month/day
          IF( llprev .AND. sdjf%num <= 0 ) THEN
@@ -401,5 +397,5 @@
             CALL fld_clopn( sdjf )
          ENDIF
-         
+         !
          IF( llprev ) THEN   ! check if the record sdjf%nrec_a(1) exists in the file
             idvar = iom_varid( sdjf%num, sdjf%clvar )                                        ! id of the variable sdjf%clvar
@@ -408,11 +404,10 @@
             sdjf%nrec_a(1) = MIN( sdjf%nrec_a(1), inrec )   ! make sure we select an existing record
          ENDIF
-
-         ! read before data in after arrays(as we will swap it later)
-         CALL fld_get( sdjf, map )
-
+         !
+         CALL fld_get( sdjf, map )         ! read before data in after arrays(as we will swap it later)
+         !
          clfmt = "('fld_init : time-interpolation for ', a, ' read previous record = ', i6, ' at time = ', f7.2, ' days')"
          IF(lwp) WRITE(numout, clfmt) TRIM(sdjf%clvar), sdjf%nrec_a(1), REAL(sdjf%nrec_a(2),wp)/rday
-
+         !
       ENDIF
       !
@@ -435,8 +430,8 @@
       LOGICAL  , INTENT(in   ), OPTIONAL ::   ldbefore  ! sent back before record values (default = .FALSE.)
       INTEGER  , INTENT(in   ), OPTIONAL ::   kit       ! index of barotropic subcycle
-                                                        ! used only if sdjf%ln_tint = .TRUE.
+      !                                                 ! used only if sdjf%ln_tint = .TRUE.
       INTEGER  , INTENT(in   ), OPTIONAL ::   kt_offset ! Offset of required time level compared to "now"
-                                                        !   time level in units of time steps.
-      !!
+      !                                                 !   time level in units of time steps.
+      !
       LOGICAL  ::   llbefore    ! local definition of ldbefore
       INTEGER  ::   iendrec     ! end of this record (in seconds)
@@ -459,5 +454,5 @@
       IF( PRESENT(kt_offset) )   it_offset = kt_offset
       IF( PRESENT(kit) ) THEN   ;   it_offset = ( kit + it_offset ) * NINT( rdt/REAL(nn_baro,wp) )
-      ELSE                      ;   it_offset =         it_offset   * NINT(       rdttra(1)      )
+      ELSE                      ;   it_offset =         it_offset   * NINT(       rdt            )
       ENDIF
       !
@@ -536,6 +531,6 @@
          ELSE                                           ;   ztmp = REAL(nsec_year ,wp)  ! since 00h on Jan 1 of the current year
          ENDIF
-         ztmp = ztmp + 0.5 * REAL(kn_fsbc - 1, wp) * rdttra(1) + REAL( it_offset, wp )  ! centrered in the middle of sbc time step
-         ztmp = ztmp + 0.01 * rdttra(1)                                                 ! avoid truncation error 
+         ztmp = ztmp + 0.5 * REAL(kn_fsbc - 1, wp) * rdt + REAL( it_offset, wp )        ! centrered in the middle of sbc time step
+         ztmp = ztmp + 0.01 * rdt                                                       ! avoid truncation error 
          IF( sdjf%ln_tint ) THEN                ! time interpolation, shift by 1/2 record
             !
@@ -592,13 +587,13 @@
       !! ** Purpose :   read the data
       !!----------------------------------------------------------------------
-      TYPE(FLD), INTENT(inout) ::   sdjf   ! input field related variables
-      TYPE(MAP_POINTER),INTENT(in) ::   map   ! global-to-local mapping indices
-      !!
-      INTEGER                  ::   ipk    ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk )
-      INTEGER                  ::   iw     ! index into wgts array
-      INTEGER                  ::   ipdom  ! index of the domain
-      INTEGER                  ::   idvar  ! variable ID
-      INTEGER                  ::   idmspc ! number of spatial dimensions
-      LOGICAL                  ::   lmoor  ! C1D case: point data
+      TYPE(FLD)        , INTENT(inout) ::   sdjf   ! input field related variables
+      TYPE(MAP_POINTER), INTENT(in   ) ::   map    ! global-to-local mapping indices
+      !
+      INTEGER ::   ipk      ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk )
+      INTEGER ::   iw       ! index into wgts array
+      INTEGER ::   ipdom    ! index of the domain
+      INTEGER ::   idvar    ! variable ID
+      INTEGER ::   idmspc   ! number of spatial dimensions
+      LOGICAL ::   lmoor    ! C1D case: point data
       !!---------------------------------------------------------------------
       !
@@ -611,19 +606,19 @@
       ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN
          CALL wgt_list( sdjf, iw )
-         IF( sdjf%ln_tint ) THEN   ;   CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk  , sdjf%fdta(:,:,:,2), & 
-              & sdjf%nrec_a(1), sdjf%lsmname )
-         ELSE                      ;   CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk  , sdjf%fnow(:,:,:  ), &
-              & sdjf%nrec_a(1), sdjf%lsmname )
+         IF( sdjf%ln_tint ) THEN   ;   CALL fld_interp( sdjf%num, sdjf%clvar, iw, ipk, sdjf%fdta(:,:,:,2),          & 
+            &                                                                          sdjf%nrec_a(1), sdjf%lsmname )
+         ELSE                      ;   CALL fld_interp( sdjf%num, sdjf%clvar, iw, ipk, sdjf%fnow(:,:,:  ),          &
+            &                                                                          sdjf%nrec_a(1), sdjf%lsmname )
          ENDIF
       ELSE
-         IF( SIZE(sdjf%fnow, 1) == jpi ) THEN  ;  ipdom = jpdom_data
-         ELSE                                  ;  ipdom = jpdom_unknown
+         IF( SIZE(sdjf%fnow, 1) == jpi ) THEN   ;   ipdom = jpdom_data
+         ELSE                                   ;   ipdom = jpdom_unknown
          ENDIF
          ! C1D case: If product of spatial dimensions == ipk, then x,y are of
          ! size 1 (point/mooring data): this must be read onto the central grid point
          idvar  = iom_varid( sdjf%num, sdjf%clvar )
-         idmspc = iom_file( sdjf%num )%ndims( idvar )
+         idmspc = iom_file ( sdjf%num )%ndims( idvar )
          IF( iom_file( sdjf%num )%luld( idvar ) )   idmspc = idmspc - 1
-         lmoor  = (idmspc == 0 .OR. PRODUCT( iom_file( sdjf%num )%dimsz( 1:MAX(idmspc,1) ,idvar ) ) == ipk)
+         lmoor  = (  idmspc == 0 .OR. PRODUCT( iom_file( sdjf%num )%dimsz( 1:MAX(idmspc,1) ,idvar ) ) == ipk  )
          !
          SELECT CASE( ipk )
@@ -660,6 +655,7 @@
       !
       sdjf%rotn(2) = .false.   ! vector not yet rotated
-
+      !
    END SUBROUTINE fld_get
+
 
    SUBROUTINE fld_map( num, clvar, dta, nrec, map )
@@ -688,9 +684,9 @@
       REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read  ! work space for global data
       !!---------------------------------------------------------------------
-            
+      !
       ipi = SIZE( dta, 1 )
       ipj = 1
       ipk = SIZE( dta, 3 )
-
+      !
       idvar   = iom_varid( num, clvar )
       ilendta = iom_file(num)%dimsz(1,idvar)
@@ -698,14 +694,14 @@
 #if defined key_bdy
       ipj = iom_file(num)%dimsz(2,idvar)
-      IF ( map%ll_unstruc) THEN ! unstructured open boundary data file
+      IF( map%ll_unstruc) THEN   ! unstructured open boundary data file
          dta_read => dta_global
-      ELSE                      ! structured open boundary data file
+      ELSE                       ! structured open boundary data file
          dta_read => dta_global2
       ENDIF
 #endif
 
-      IF(lwp) WRITE(numout,*) 'Dim size for ',TRIM(clvar),' is ', ilendta
+      IF(lwp) WRITE(numout,*) 'Dim size for ',        TRIM(clvar),' is ', ilendta
       IF(lwp) WRITE(numout,*) 'Number of levels for ',TRIM(clvar),' is ', ipk
-
+      !
       SELECT CASE( ipk )
       CASE(1)        ;   CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1    ), nrec )
@@ -713,5 +709,5 @@
       END SELECT
       !
-      IF ( map%ll_unstruc ) THEN ! unstructured open boundary data file
+      IF( map%ll_unstruc ) THEN  ! unstructured open boundary data file
          DO ib = 1, ipi
             DO ik = 1, ipk
@@ -728,5 +724,5 @@
          END DO
       ENDIF
-
+      !
    END SUBROUTINE fld_map
 
@@ -738,17 +734,17 @@
       !! ** Purpose :   Vector fields may need to be rotated onto the local grid direction
       !!----------------------------------------------------------------------
-      INTEGER  , INTENT(in   )               ::   kt        ! ocean time step
-      TYPE(FLD), INTENT(inout), DIMENSION(:) ::   sd        ! input field related variables
-      !!
-      INTEGER                           ::   ju,jv,jk,jn  ! loop indices
-      INTEGER                           ::   imf          ! size of the structure sd
-      INTEGER                           ::   ill          ! character length
-      INTEGER                           ::   iv           ! indice of V component
+      INTEGER                , INTENT(in   ) ::   kt   ! ocean time step
+      TYPE(FLD), DIMENSION(:), INTENT(inout) ::   sd   ! input field related variables
+      !
+      INTEGER ::   ju, jv, jk, jn  ! loop indices
+      INTEGER ::   imf             ! size of the structure sd
+      INTEGER ::   ill             ! character length
+      INTEGER ::   iv              ! indice of V component
+      CHARACTER (LEN=100)               ::   clcomp       ! dummy weight name
       REAL(wp), POINTER, DIMENSION(:,:) ::   utmp, vtmp   ! temporary arrays for vector rotation
-      CHARACTER (LEN=100)               ::   clcomp       ! dummy weight name
-      !!---------------------------------------------------------------------
-
-      CALL wrk_alloc( jpi,jpj, utmp, vtmp )
-
+      !!---------------------------------------------------------------------
+      !
+      CALL wrk_alloc( jpi,jpj,   utmp, vtmp )
+      !
       !! (sga: following code should be modified so that pairs arent searched for each time
       !
@@ -786,5 +782,5 @@
        END DO
       !
-      CALL wrk_dealloc( jpi,jpj, utmp, vtmp )
+      CALL wrk_dealloc( jpi,jpj,   utmp, vtmp )
       !
    END SUBROUTINE fld_rot
@@ -802,5 +798,5 @@
       INTEGER, OPTIONAL, INTENT(in   ) ::   kday     ! day value
       LOGICAL, OPTIONAL, INTENT(in   ) ::   ldstop   ! stop if open to read a non-existing file (default = .TRUE.)
-      !!
+      !
       LOGICAL :: llprevyr              ! are we reading previous year  file?
       LOGICAL :: llprevmth             ! are we reading previous month file?
@@ -853,9 +849,9 @@
       !
       IF( TRIM(clname) /= TRIM(sdjf%clname) .OR. sdjf%num == 0 ) THEN   ! new file to be open 
-
+         !
          sdjf%clname = TRIM(clname)
          IF( sdjf%num /= 0 )   CALL iom_close( sdjf%num )   ! close file if already open
          CALL iom_open( sdjf%clname, sdjf%num, ldstop = ldstop, ldiof =  LEN(TRIM(sdjf%wgtname)) > 0 )
-
+         !
          ! find the last record to be read -> update sdjf%nreclast
          indexyr = iyear - nyear + 1
@@ -866,5 +862,5 @@
          CASE ( 2 )   ;   imonth_len = 31   ! next     year -> imonth = 1
          END SELECT
-         
+         !
          ! last record to be read in the current file
          IF    ( sdjf%nfreqh == -12 ) THEN                 ;   sdjf%nreclast = 1    !  yearly mean
@@ -880,5 +876,5 @@
             ENDIF
          ENDIF
-         
+         !
       ENDIF
       !
@@ -901,5 +897,5 @@
       INTEGER  ::   jf       ! dummy indices
       !!---------------------------------------------------------------------
-
+      !
       DO jf = 1, SIZE(sdf)
          sdf(jf)%clrootname = TRIM( cdir )//TRIM( sdf_n(jf)%clname )
@@ -923,5 +919,5 @@
          sdf(jf)%nreclast = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn
       END DO
-
+      !
       IF(lwp) THEN      ! control print
          WRITE(numout,*)
@@ -943,5 +939,5 @@
          END DO
       ENDIF
-      
+      !
    END SUBROUTINE fld_fill
 
@@ -958,5 +954,5 @@
       TYPE( FLD ), INTENT(in   ) ::   sd        ! field with name of weights file
       INTEGER    , INTENT(inout) ::   kwgt      ! index of weights
-      !!
+      !
       INTEGER ::   kw, nestid   ! local integer
       LOGICAL ::   found        ! local logical
@@ -966,5 +962,5 @@
       !! weights filename is either present or we hit the end of the list
       found = .FALSE.
-
+      !
       !! because agrif nest part of filenames are now added in iom_open
       !! to distinguish between weights files on the different grids, need to track
@@ -1028,18 +1024,18 @@
       TYPE( FLD ), INTENT(in) ::   sd   ! field with name of weights file
       !!
-      INTEGER                           ::   jn            ! dummy loop indices
-      INTEGER                           ::   inum          ! temporary logical unit
-      INTEGER                           ::   id            ! temporary variable id
-      INTEGER                           ::   ipk           ! temporary vertical dimension
-      CHARACTER (len=5)                 ::   aname
+      INTEGER ::   jn         ! dummy loop indices
+      INTEGER ::   inum       ! local logical unit
+      INTEGER ::   id         ! local variable id
+      INTEGER ::   ipk        ! local vertical dimension
+      INTEGER ::   zwrap      ! local integer
+      LOGICAL ::   cyclical   ! 
+      CHARACTER (len=5) ::   aname   !
       INTEGER , DIMENSION(:), ALLOCATABLE ::   ddims
       INTEGER , POINTER, DIMENSION(:,:) ::   data_src
       REAL(wp), POINTER, DIMENSION(:,:) ::   data_tmp
-      LOGICAL                           ::   cyclical
-      INTEGER                           ::   zwrap      ! local integer
-      !!----------------------------------------------------------------------
-      !
-      CALL wrk_alloc( jpi,jpj, data_src )   ! integer
-      CALL wrk_alloc( jpi,jpj, data_tmp )
+      !!----------------------------------------------------------------------
+      !
+      CALL wrk_alloc( jpi,jpj,   data_src )   ! integer
+      CALL wrk_alloc( jpi,jpj,   data_tmp )
       !
       IF( nxt_wgt > tot_wgts ) THEN
@@ -1151,7 +1147,7 @@
          ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3 ,ipk) )
          IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col(1,ref_wgts(nxt_wgt)%jpjwgt+3,ipk) )
-
+         !
          nxt_wgt = nxt_wgt + 1
-
+         !
       ELSE 
          CALL ctl_stop( '    fld_weight : unable to read the file ' )
@@ -1166,6 +1162,6 @@
 
 
-   SUBROUTINE apply_seaoverland(clmaskfile,zfieldo,jpi1_lsm,jpi2_lsm,jpj1_lsm, &
-                          &      jpj2_lsm,itmpi,itmpj,itmpz,rec1_lsm,recn_lsm)
+   SUBROUTINE apply_seaoverland( clmaskfile, zfieldo, jpi1_lsm, jpi2_lsm, jpj1_lsm,   &
+                          &      jpj2_lsm, itmpi, itmpj, itmpz, rec1_lsm, recn_lsm )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE apply_seaoverland  ***
@@ -1176,67 +1172,61 @@
       !!      D. Delrosso INGV          
       !!---------------------------------------------------------------------- 
-      INTEGER                                   :: inum,jni,jnj,jnz,jc                  ! temporary indices
-      INTEGER,                   INTENT(in)     :: itmpi,itmpj,itmpz                    ! lengths
-      INTEGER,                   INTENT(in)     :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm  ! temporary indices
-      INTEGER, DIMENSION(3),     INTENT(in)     :: rec1_lsm,recn_lsm                    ! temporary arrays for start and length
-      REAL(wp),DIMENSION (:,:,:),INTENT(inout)  :: zfieldo                              ! input/output array for seaoverland application
-      REAL(wp),DIMENSION (:,:,:),ALLOCATABLE    :: zslmec1                              ! temporary array for land point detection
-      REAL(wp),DIMENSION (:,:),  ALLOCATABLE    :: zfieldn                              ! array of forcing field with undeff for land points
-      REAL(wp),DIMENSION (:,:),  ALLOCATABLE    :: zfield                               ! array of forcing field
-      CHARACTER (len=100),       INTENT(in)     :: clmaskfile                           ! land/sea mask file name
-      !!---------------------------------------------------------------------
-      ALLOCATE ( zslmec1(itmpi,itmpj,itmpz) )
-      ALLOCATE ( zfieldn(itmpi,itmpj) )
-      ALLOCATE ( zfield(itmpi,itmpj) )
-
+      INTEGER,                   INTENT(in   ) :: itmpi,itmpj,itmpz                    ! lengths
+      INTEGER,                   INTENT(in   ) :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm  ! temporary indices
+      INTEGER, DIMENSION(3),     INTENT(in   ) :: rec1_lsm,recn_lsm                    ! temporary arrays for start and length
+      REAL(wp),DIMENSION (:,:,:),INTENT(inout) :: zfieldo                              ! input/output array for seaoverland application
+      CHARACTER (len=100),       INTENT(in   ) :: clmaskfile                           ! land/sea mask file name
+      !
+      INTEGER :: inum,jni,jnj,jnz,jc   ! local indices
+      REAL(wp),DIMENSION (:,:,:),ALLOCATABLE :: zslmec1             ! local array for land point detection
+      REAL(wp),DIMENSION (:,:),  ALLOCATABLE :: zfieldn   ! array of forcing field with undeff for land points
+      REAL(wp),DIMENSION (:,:),  ALLOCATABLE :: zfield    ! array of forcing field
+      !!---------------------------------------------------------------------
+      !
+      ALLOCATE ( zslmec1(itmpi,itmpj,itmpz), zfieldn(itmpi,itmpj), zfield(itmpi,itmpj) )
+      !
       ! Retrieve the land sea mask data
       CALL iom_open( clmaskfile, inum )
       SELECT CASE( SIZE(zfieldo(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) )
       CASE(1)
-           CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), 1, rec1_lsm, recn_lsm)
+         CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), 1, rec1_lsm, recn_lsm)
       CASE DEFAULT
-           CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), 1, rec1_lsm, recn_lsm)
+         CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), 1, rec1_lsm, recn_lsm)
       END SELECT
       CALL iom_close( inum )
-
-      DO jnz=1,rec1_lsm(3)                            !! Loop over k dimension
-
-         DO jni=1,itmpi                               !! copy the original field into a tmp array
-            DO jnj=1,itmpj                            !! substituting undeff over land points
-            zfieldn(jni,jnj) = zfieldo(jni,jnj,jnz)
-               IF ( zslmec1(jni,jnj,jnz) == 1. ) THEN
-                  zfieldn(jni,jnj) = undeff_lsm
-               ENDIF
+      !
+      DO jnz=1,rec1_lsm(3)             !! Loop over k dimension
+         !
+         DO jni = 1, itmpi                               !! copy the original field into a tmp array
+            DO jnj = 1, itmpj                            !! substituting undeff over land points
+               zfieldn(jni,jnj) = zfieldo(jni,jnj,jnz)
+               IF( zslmec1(jni,jnj,jnz) == 1. )   zfieldn(jni,jnj) = undeff_lsm
             END DO
          END DO
-  
-      CALL seaoverland(zfieldn,itmpi,itmpj,zfield)
-      DO jc=1,nn_lsm
-         CALL seaoverland(zfield,itmpi,itmpj,zfield)
-      END DO
-
-      !   Check for Undeff and substitute original values
-      IF(ANY(zfield==undeff_lsm)) THEN
-         DO jni=1,itmpi
-            DO jnj=1,itmpj
-               IF (zfield(jni,jnj)==undeff_lsm) THEN
-                  zfield(jni,jnj) = zfieldo(jni,jnj,jnz)
-               ENDIF
-            ENDDO
-         ENDDO
-      ENDIF
-
-      zfieldo(:,:,jnz)=zfield(:,:)
-
-      END DO                          !! End Loop over k dimension
-
-      DEALLOCATE ( zslmec1 )
-      DEALLOCATE ( zfieldn )
-      DEALLOCATE ( zfield )
-
+         !
+         CALL seaoverland( zfieldn, itmpi, itmpj, zfield )
+         DO jc = 1, nn_lsm
+            CALL seaoverland( zfield, itmpi, itmpj, zfield )
+         END DO
+         !
+         !   Check for Undeff and substitute original values
+         IF( ANY(zfield==undeff_lsm) ) THEN
+            DO jni = 1, itmpi
+               DO jnj = 1, itmpj
+                  IF( zfield(jni,jnj)==undeff_lsm )   zfield(jni,jnj) = zfieldo(jni,jnj,jnz)
+               END DO
+            END DO
+         ENDIF
+         !
+         zfieldo(:,:,jnz) = zfield(:,:)
+         !
+      END DO                           !! End Loop over k dimension
+      !
+      DEALLOCATE ( zslmec1, zfieldn, zfield )
+      !
    END SUBROUTINE apply_seaoverland 
 
 
-   SUBROUTINE seaoverland(zfieldn,ileni,ilenj,zfield)
+   SUBROUTINE seaoverland( zfieldn, ileni, ilenj, zfield )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE seaoverland  ***
@@ -1245,29 +1235,31 @@
       !!      D. Delrosso INGV
       !!---------------------------------------------------------------------- 
-      INTEGER,INTENT(in)                       :: ileni,ilenj              ! lengths 
-      REAL,DIMENSION (ileni,ilenj),INTENT(in)  :: zfieldn                  ! array of forcing field with undeff for land points
-      REAL,DIMENSION (ileni,ilenj),INTENT(out) :: zfield                   ! array of forcing field
-      REAL,DIMENSION (ileni,ilenj)             :: zmat1,zmat2,zmat3,zmat4  ! temporary arrays for seaoverland application
-      REAL,DIMENSION (ileni,ilenj)             :: zmat5,zmat6,zmat7,zmat8  ! temporary arrays for seaoverland application
-      REAL,DIMENSION (ileni,ilenj)             :: zlsm2d                   ! temporary arrays for seaoverland application
-      REAL,DIMENSION (ileni,ilenj,8)           :: zlsm3d                   ! temporary arrays for seaoverland application
-      LOGICAL,DIMENSION (ileni,ilenj,8)        :: ll_msknan3d              ! logical mask for undeff detection
-      LOGICAL,DIMENSION (ileni,ilenj)          :: ll_msknan2d              ! logical mask for undeff detection
+      INTEGER                      , INTENT(in   ) :: ileni,ilenj   ! lengths 
+      REAL, DIMENSION (ileni,ilenj), INTENT(in   ) :: zfieldn       ! array of forcing field with undeff for land points
+      REAL, DIMENSION (ileni,ilenj), INTENT(  out) :: zfield        ! array of forcing field
+      !
+      REAL   , DIMENSION (ileni,ilenj)   :: zmat1, zmat2, zmat3, zmat4  ! local arrays 
+      REAL   , DIMENSION (ileni,ilenj)   :: zmat5, zmat6, zmat7, zmat8  !   -     - 
+      REAL   , DIMENSION (ileni,ilenj)   :: zlsm2d                      !   -     - 
+      REAL   , DIMENSION (ileni,ilenj,8) :: zlsm3d                      !   -     -
+      LOGICAL, DIMENSION (ileni,ilenj,8) :: ll_msknan3d                 ! logical mask for undeff detection
+      LOGICAL, DIMENSION (ileni,ilenj)   :: ll_msknan2d                 ! logical mask for undeff detection
       !!---------------------------------------------------------------------- 
-      zmat8 = eoshift(zfieldn   ,  SHIFT=-1, BOUNDARY = (/zfieldn(:,1)/)    ,DIM=2)
-      zmat1 = eoshift(zmat8     ,  SHIFT=-1, BOUNDARY = (/zmat8(1,:)/)      ,DIM=1)
-      zmat2 = eoshift(zfieldn   ,  SHIFT=-1, BOUNDARY = (/zfieldn(1,:)/)    ,DIM=1)
-      zmat4 = eoshift(zfieldn   ,  SHIFT= 1, BOUNDARY = (/zfieldn(:,ilenj)/),DIM=2)
-      zmat3 = eoshift(zmat4     ,  SHIFT=-1, BOUNDARY = (/zmat4(1,:)/)      ,DIM=1)
-      zmat5 = eoshift(zmat4     ,  SHIFT= 1, BOUNDARY = (/zmat4(ileni,:)/)  ,DIM=1)
-      zmat6 = eoshift(zfieldn   ,  SHIFT= 1, BOUNDARY = (/zfieldn(ileni,:)/),DIM=1)
-      zmat7 = eoshift(zmat8     ,  SHIFT= 1, BOUNDARY = (/zmat8(ileni,:)/)  ,DIM=1)
-
+      zmat8 = eoshift( zfieldn , SHIFT=-1 , BOUNDARY = (/zfieldn(:,1)/)     , DIM=2 )
+      zmat1 = eoshift( zmat8   , SHIFT=-1 , BOUNDARY = (/zmat8(1,:)/)       , DIM=1 )
+      zmat2 = eoshift( zfieldn , SHIFT=-1 , BOUNDARY = (/zfieldn(1,:)/)     , DIM=1 )
+      zmat4 = eoshift( zfieldn , SHIFT= 1 , BOUNDARY = (/zfieldn(:,ilenj)/) , DIM=2 )
+      zmat3 = eoshift( zmat4   , SHIFT=-1 , BOUNDARY = (/zmat4(1,:)/)       , DIM=1 )
+      zmat5 = eoshift( zmat4   , SHIFT= 1 , BOUNDARY = (/zmat4(ileni,:)/)   , DIM=1 )
+      zmat6 = eoshift( zfieldn , SHIFT= 1 , BOUNDARY = (/zfieldn(ileni,:)/) , DIM=1 )
+      zmat7 = eoshift( zmat8   , SHIFT= 1 , BOUNDARY = (/zmat8(ileni,:)/)   , DIM=1 )
+      !
       zlsm3d  = RESHAPE( (/ zmat1, zmat2, zmat3, zmat4, zmat5, zmat6, zmat7, zmat8 /), (/ ileni, ilenj, 8 /))
-      ll_msknan3d = .not.(zlsm3d==undeff_lsm)
-      ll_msknan2d = .not.(zfieldn==undeff_lsm)  ! FALSE where is Undeff (land)
-      zlsm2d = (SUM ( zlsm3d, 3 , ll_msknan3d ) )/(MAX(1,(COUNT( ll_msknan3d , 3 ))   ))
-      WHERE ((COUNT( ll_msknan3d , 3 )) == 0.0_wp)  zlsm2d = undeff_lsm
-      zfield = MERGE (zfieldn,zlsm2d,ll_msknan2d)
+      ll_msknan3d = .NOT.( zlsm3d  == undeff_lsm )
+      ll_msknan2d = .NOT.( zfieldn == undeff_lsm )  ! FALSE where is Undeff (land)
+      zlsm2d = SUM( zlsm3d, 3 , ll_msknan3d ) / MAX( 1 , COUNT( ll_msknan3d , 3 ) )
+      WHERE( COUNT( ll_msknan3d , 3 ) == 0._wp )   zlsm2d = undeff_lsm
+      zfield = MERGE( zfieldn, zlsm2d, ll_msknan2d )
+      !
    END SUBROUTINE seaoverland
 
@@ -1288,19 +1280,18 @@
       INTEGER                   , INTENT(in   ) ::   nrec    ! record number to read (ie time slice)
       CHARACTER(LEN=*)          , INTENT(in   ) ::   lsmfile ! land sea mask file name
-      !! 
-      REAL(wp),DIMENSION(:,:,:),ALLOCATABLE     ::   ztmp_fly_dta,zfieldo                  ! temporary array of values on input grid
-      INTEGER, DIMENSION(3)                     ::   rec1,recn                             ! temporary arrays for start and length
-      INTEGER, DIMENSION(3)                     ::   rec1_lsm,recn_lsm                     ! temporary arrays for start and length in case of seaoverland
-      INTEGER                                   ::   ii_lsm1,ii_lsm2,ij_lsm1,ij_lsm2       ! temporary indices
-      INTEGER                                   ::   jk, jn, jm, jir, jjr                  ! loop counters
-      INTEGER                                   ::   ni, nj                                ! lengths
-      INTEGER                                   ::   jpimin,jpiwid                         ! temporary indices
-      INTEGER                                   ::   jpimin_lsm,jpiwid_lsm                 ! temporary indices
-      INTEGER                                   ::   jpjmin,jpjwid                         ! temporary indices
-      INTEGER                                   ::   jpjmin_lsm,jpjwid_lsm                 ! temporary indices
-      INTEGER                                   ::   jpi1,jpi2,jpj1,jpj2                   ! temporary indices
-      INTEGER                                   ::   jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm   ! temporary indices
-      INTEGER                                   ::   itmpi,itmpj,itmpz                     ! lengths
-      
+      !
+      INTEGER, DIMENSION(3) ::   rec1, recn           ! temporary arrays for start and length
+      INTEGER, DIMENSION(3) ::   rec1_lsm, recn_lsm   ! temporary arrays for start and length in case of seaoverland
+      INTEGER ::   ii_lsm1,ii_lsm2,ij_lsm1,ij_lsm2    ! temporary indices
+      INTEGER ::   jk, jn, jm, jir, jjr               ! loop counters
+      INTEGER ::   ni, nj                             ! lengths
+      INTEGER ::   jpimin,jpiwid                      ! temporary indices
+      INTEGER ::   jpimin_lsm,jpiwid_lsm              ! temporary indices
+      INTEGER ::   jpjmin,jpjwid                      ! temporary indices
+      INTEGER ::   jpjmin_lsm,jpjwid_lsm              ! temporary indices
+      INTEGER ::   jpi1,jpi2,jpj1,jpj2                ! temporary indices
+      INTEGER ::   jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm   ! temporary indices
+      INTEGER ::   itmpi,itmpj,itmpz                     ! lengths
+      REAL(wp),DIMENSION(:,:,:), ALLOCATABLE ::   ztmp_fly_dta                 ! local array of values on input grid     
       !!----------------------------------------------------------------------
       !
@@ -1355,6 +1346,6 @@
 
 
-         itmpi=SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),1)
-         itmpj=SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),2)
+         itmpi=jpi2_lsm-jpi1_lsm+1
+         itmpj=jpj2_lsm-jpj1_lsm+1
          itmpz=kk
          ALLOCATE(ztmp_fly_dta(itmpi,itmpj,itmpz))
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90	(revision 6121)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90	(revision 6140)
@@ -5,39 +5,41 @@
    !!======================================================================
    !! History :  OPA  !  07-1996  (O. Marti)  Original code
-   !!   NEMO     1.0  !  02-2008  (G. Madec)  F90: Free form
-   !!            3.0  !  
+   !!   NEMO     1.0  !  06-2006  (G. Madec )  Free form, F90 + opt.
+   !!                 !  04-2007  (S. Masson)  angle: Add T, F points and bugfix in cos lateral boundary
+   !!            3.0  !  07-2008  (G. Madec)  geo2oce suppress lon/lat agruments
+   !!            3.7  !  11-2015  (G. Madec)  remove the unused repere and repcmo routines
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
-   !!   repcmo      : 
-   !!   angle       :
-   !!   geo2oce     :
-   !!   repere      :   old routine suppress it ???
+   !!   rot_rep       : Rotate the Repere: geographic grid <==> stretched coordinates grid
+   !!   angle         :
+   !!   geo2oce       :
+   !!   oce2geo       :
    !!----------------------------------------------------------------------
-   USE dom_oce         ! mesh and scale factors
-   USE phycst          ! physical constants
-   USE in_out_manager  ! I/O manager
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE lib_mpp         ! MPP library
+   USE dom_oce        ! mesh and scale factors
+   USE phycst         ! physical constants
+   !
+   USE in_out_manager ! I/O manager
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp        ! MPP library
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   rot_rep, repcmo, repere, geo2oce, oce2geo   ! only rot_rep should be used
-                                             ! repcmo and repere are keep only for compatibility.
-                                             ! they are only a useless overlay of rot_rep
-
-   PUBLIC   obs_rot
-
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   &
-      gsint, gcost,   &  ! cos/sin between model grid lines and NP direction at T point
-      gsinu, gcosu,   &  ! cos/sin between model grid lines and NP direction at U point
-      gsinv, gcosv,   &  ! cos/sin between model grid lines and NP direction at V point
-      gsinf, gcosf       ! cos/sin between model grid lines and NP direction at F point
+   PUBLIC   rot_rep   ! called in sbccpl, fldread, and cyclone
+   PUBLIC   geo2oce   ! called in sbccpl
+   PUBLIC   oce2geo   ! called in sbccpl
+   PUBLIC   obs_rot   ! called in obs_rot_vel and obs_write
+
+   !                                         ! cos/sin between model grid lines and NP direction
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gsint, gcost   ! at T point
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gsinu, gcosu   ! at U point
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gsinv, gcosv   ! at V point
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gsinf, gcosf   ! at F point
 
    LOGICAL ,              SAVE, DIMENSION(4)     ::   linit = .FALSE.
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gsinlon, gcoslon, gsinlat, gcoslat
 
-   LOGICAL ::   lmust_init = .TRUE.        !: used to initialize the cos/sin variables (se above)
+   LOGICAL ::   lmust_init = .TRUE.        !: used to initialize the cos/sin variables (see above)
 
    !! * Substitutions
@@ -50,31 +52,4 @@
 CONTAINS
 
-   SUBROUTINE repcmo ( pxu1, pyu1, pxv1, pyv1,   &
-                       px2 , py2 )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE repcmo  ***
-      !!
-      !! ** Purpose :   Change vector componantes from a geographic grid to a
-      !!      stretched coordinates grid.
-      !!
-      !! ** Method  :   Initialization of arrays at the first call.
-      !!
-      !! ** Action  : - px2 : first  componante (defined at u point)
-      !!              - py2 : second componante (defined at v point)
-      !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   pxu1, pyu1   ! geographic vector componantes at u-point
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   pxv1, pyv1   ! geographic vector componantes at v-point
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   px2          ! i-componante (defined at u-point)
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   py2          ! j-componante (defined at v-point)
-      !!----------------------------------------------------------------------
-      
-      ! Change from geographic to stretched coordinate
-      ! ----------------------------------------------
-      CALL rot_rep( pxu1, pyu1, 'U', 'en->i',px2 )
-      CALL rot_rep( pxv1, pyv1, 'V', 'en->j',py2 )
-      
-   END SUBROUTINE repcmo
-
-
    SUBROUTINE rot_rep ( pxin, pyin, cd_type, cdtodo, prot )
       !!----------------------------------------------------------------------
@@ -83,33 +58,27 @@
       !! ** Purpose :   Rotate the Repere: Change vector componantes between
       !!                geographic grid <--> stretched coordinates grid.
-      !!
-      !! History :
-      !!   9.2  !  07-04  (S. Masson)  
-      !!                  (O. Marti ) Original code (repere and repcmo)
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT( IN ) ::   pxin, pyin   ! vector componantes
-      CHARACTER(len=1),             INTENT( IN ) ::   cd_type      ! define the nature of pt2d array grid-points
-      CHARACTER(len=5),             INTENT( IN ) ::   cdtodo       ! specify the work to do:
-      !!                                                           ! 'en->i' east-north componantes to model i componante
-      !!                                                           ! 'en->j' east-north componantes to model j componante
-      !!                                                           ! 'ij->e' model i-j componantes to east componante
-      !!                                                           ! 'ij->n' model i-j componantes to east componante
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) ::   prot      
-      !!----------------------------------------------------------------------
-
-      ! Initialization of gsin* and gcos* at first call
-      ! -----------------------------------------------
-
-      IF( lmust_init ) THEN
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pxin, pyin   ! vector componantes
+      CHARACTER(len=1),             INTENT(in   ) ::   cd_type      ! define the nature of pt2d array grid-points
+      CHARACTER(len=5),             INTENT(in   ) ::   cdtodo       ! type of transpormation:
+      !                                                             ! 'en->i' = east-north to i-component
+      !                                                             ! 'en->j' = east-north to j-component
+      !                                                             ! 'ij->e' = (i,j) components to east
+      !                                                             ! 'ij->n' = (i,j) components to north
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   prot      
+      !!----------------------------------------------------------------------
+      !
+      IF( lmust_init ) THEN      ! at 1st call only: set  gsin. & gcos.
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) ' rot_rep : geographic <--> stretched'
-         IF(lwp) WRITE(numout,*) ' ~~~~~    coordinate transformation'
+         IF(lwp) WRITE(numout,*) ' rot_rep: coordinate transformation : geographic <==> model (i,j)-components'
+         IF(lwp) WRITE(numout,*) ' ~~~~~~~~    '
          !
          CALL angle       ! initialization of the transformation
          lmust_init = .FALSE.
       ENDIF
-      
-      SELECT CASE (cdtodo)
-      CASE ('en->i')      ! 'en->i' est-north componantes to model i componante
+      !
+      SELECT CASE( cdtodo )      ! type of rotation
+      !
+      CASE( 'en->i' )                  ! east-north to i-component
          SELECT CASE (cd_type)
          CASE ('T')   ;   prot(:,:) = pxin(:,:) * gcost(:,:) + pyin(:,:) * gsint(:,:)
@@ -119,5 +88,5 @@
          CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
          END SELECT
-      CASE ('en->j')      ! 'en->j' est-north componantes to model j componante
+      CASE ('en->j')                   ! east-north to j-component
          SELECT CASE (cd_type)
          CASE ('T')   ;   prot(:,:) = pyin(:,:) * gcost(:,:) - pxin(:,:) * gsint(:,:)
@@ -127,5 +96,5 @@
          CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
          END SELECT
-      CASE ('ij->e')      ! 'ij->e' model i-j componantes to est componante
+      CASE ('ij->e')                   ! (i,j)-components to east
          SELECT CASE (cd_type)
          CASE ('T')   ;   prot(:,:) = pxin(:,:) * gcost(:,:) - pyin(:,:) * gsint(:,:)
@@ -135,5 +104,5 @@
          CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
          END SELECT
-      CASE ('ij->n')      ! 'ij->n' model i-j componantes to est componante
+      CASE ('ij->n')                   ! (i,j)-components to north 
          SELECT CASE (cd_type)
          CASE ('T')   ;   prot(:,:) = pyin(:,:) * gcost(:,:) + pxin(:,:) * gsint(:,:)
@@ -144,6 +113,7 @@
          END SELECT
       CASE DEFAULT   ;   CALL ctl_stop( 'rot_rep: Syntax Error in the definition of cdtodo' )
+      !
       END SELECT
-      
+      !
    END SUBROUTINE rot_rep
 
@@ -155,31 +125,24 @@
       !! ** Purpose :   Compute angles between model grid lines and the North direction
       !!
-      !! ** Method  :
-      !!
-      !! ** Action  :   Compute (gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf) arrays:
-      !!      sinus and cosinus of the angle between the north-south axe and the 
-      !!      j-direction at t, u, v and f-points
-      !!
-      !! History :
-      !!   7.0  !  96-07  (O. Marti )  Original code
-      !!   8.0  !  98-06  (G. Madec )
-      !!   8.5  !  98-06  (G. Madec )  Free form, F90 + opt.
-      !!   9.2  !  07-04  (S. Masson)  Add T, F points and bugfix in cos lateral boundary
-      !!----------------------------------------------------------------------
-      INTEGER ::   ji, jj   ! dummy loop indices
-      INTEGER ::   ierr     ! local integer
-      REAL(wp) ::   &
-         zlam, zphi,            &  ! temporary scalars
-         zlan, zphh,            &  !    "         "
-         zxnpt, zynpt, znnpt,   &  ! x,y components and norm of the vector: T point to North Pole
-         zxnpu, zynpu, znnpu,   &  ! x,y components and norm of the vector: U point to North Pole
-         zxnpv, zynpv, znnpv,   &  ! x,y components and norm of the vector: V point to North Pole
-         zxnpf, zynpf, znnpf,   &  ! x,y components and norm of the vector: F point to North Pole
-         zxvvt, zyvvt, znvvt,   &  ! x,y components and norm of the vector: between V points below and above a T point
-         zxffu, zyffu, znffu,   &  ! x,y components and norm of the vector: between F points below and above a U point
-         zxffv, zyffv, znffv,   &  ! x,y components and norm of the vector: between F points left  and right a V point
-         zxuuf, zyuuf, znuuf       ! x,y components and norm of the vector: between U points below and above a F point
-      !!----------------------------------------------------------------------
-
+      !! ** Method  :   sinus and cosinus of the angle between the north-south axe 
+      !!              and the j-direction at t, u, v and f-points
+      !!                dot and cross products are used to obtain cos and sin, resp.
+      !!
+      !! ** Action  : - gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji, jj   ! dummy loop indices
+      INTEGER  ::   ierr     ! local integer
+      REAL(wp) ::   zlam, zphi            ! local scalars
+      REAL(wp) ::   zlan, zphh            !   -      -
+      REAL(wp) ::   zxnpt, zynpt, znnpt   ! x,y components and norm of the vector: T point to North Pole
+      REAL(wp) ::   zxnpu, zynpu, znnpu   ! x,y components and norm of the vector: U point to North Pole
+      REAL(wp) ::   zxnpv, zynpv, znnpv   ! x,y components and norm of the vector: V point to North Pole
+      REAL(wp) ::   zxnpf, zynpf, znnpf   ! x,y components and norm of the vector: F point to North Pole
+      REAL(wp) ::   zxvvt, zyvvt, znvvt   ! x,y components and norm of the vector: between V points below and above a T point
+      REAL(wp) ::   zxffu, zyffu, znffu   ! x,y components and norm of the vector: between F points below and above a U point
+      REAL(wp) ::   zxffv, zyffv, znffv   ! x,y components and norm of the vector: between F points left  and right a V point
+      REAL(wp) ::   zxuuf, zyuuf, znuuf   ! x,y components and norm of the vector: between U points below and above a F point
+      !!----------------------------------------------------------------------
+      !
       ALLOCATE( gsint(jpi,jpj), gcost(jpi,jpj),   & 
          &      gsinu(jpi,jpj), gcosu(jpi,jpj),   & 
@@ -187,44 +150,39 @@
          &      gsinf(jpi,jpj), gcosf(jpi,jpj), STAT=ierr )
       IF(lk_mpp)   CALL mpp_sum( ierr )
-      IF( ierr /= 0 )   CALL ctl_stop('angle: unable to allocate arrays' )
-
+      IF( ierr /= 0 )   CALL ctl_stop( 'angle: unable to allocate arrays' )
+      !
       ! ============================= !
       ! Compute the cosinus and sinus !
       ! ============================= !
       ! (computation done on the north stereographic polar plane)
-
+      !
       DO jj = 2, jpjm1
          DO ji = fs_2, jpi   ! vector opt.
-
-            ! north pole direction & modulous (at t-point)
-            zlam = glamt(ji,jj)
+            !                  
+            zlam = glamt(ji,jj)     ! north pole direction & modulous (at t-point)
             zphi = gphit(ji,jj)
             zxnpt = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
             zynpt = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
             znnpt = zxnpt*zxnpt + zynpt*zynpt
-
-            ! north pole direction & modulous (at u-point)
-            zlam = glamu(ji,jj)
+            !
+            zlam = glamu(ji,jj)     ! north pole direction & modulous (at u-point)
             zphi = gphiu(ji,jj)
             zxnpu = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
             zynpu = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
             znnpu = zxnpu*zxnpu + zynpu*zynpu
-
-            ! north pole direction & modulous (at v-point)
-            zlam = glamv(ji,jj)
+            !
+            zlam = glamv(ji,jj)     ! north pole direction & modulous (at v-point)
             zphi = gphiv(ji,jj)
             zxnpv = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
             zynpv = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
             znnpv = zxnpv*zxnpv + zynpv*zynpv
-
-            ! north pole direction & modulous (at f-point)
-            zlam = glamf(ji,jj)
+            !
+            zlam = glamf(ji,jj)     ! north pole direction & modulous (at f-point)
             zphi = gphif(ji,jj)
             zxnpf = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
             zynpf = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
             znnpf = zxnpf*zxnpf + zynpf*zynpf
-
-            ! j-direction: v-point segment direction (around t-point)
-            zlam = glamv(ji,jj  )
+            !
+            zlam = glamv(ji,jj  )   ! j-direction: v-point segment direction (around t-point)
             zphi = gphiv(ji,jj  )
             zlan = glamv(ji,jj-1)
@@ -236,7 +194,6 @@
             znvvt = SQRT( znnpt * ( zxvvt*zxvvt + zyvvt*zyvvt )  )
             znvvt = MAX( znvvt, 1.e-14 )
-
-            ! j-direction: f-point segment direction (around u-point)
-            zlam = glamf(ji,jj  )
+            !
+            zlam = glamf(ji,jj  )   ! j-direction: f-point segment direction (around u-point)
             zphi = gphif(ji,jj  )
             zlan = glamf(ji,jj-1)
@@ -248,7 +205,6 @@
             znffu = SQRT( znnpu * ( zxffu*zxffu + zyffu*zyffu )  )
             znffu = MAX( znffu, 1.e-14 )
-
-            ! i-direction: f-point segment direction (around v-point)
-            zlam = glamf(ji  ,jj)
+            !
+            zlam = glamf(ji  ,jj)   ! i-direction: f-point segment direction (around v-point)
             zphi = gphif(ji  ,jj)
             zlan = glamf(ji-1,jj)
@@ -260,7 +216,6 @@
             znffv = SQRT( znnpv * ( zxffv*zxffv + zyffv*zyffv )  )
             znffv = MAX( znffv, 1.e-14 )
-
-            ! j-direction: u-point segment direction (around f-point)
-            zlam = glamu(ji,jj+1)
+            !
+            zlam = glamu(ji,jj+1)   ! j-direction: u-point segment direction (around f-point)
             zphi = gphiu(ji,jj+1)
             zlan = glamu(ji,jj  )
@@ -272,19 +227,18 @@
             znuuf = SQRT( znnpf * ( zxuuf*zxuuf + zyuuf*zyuuf )  )
             znuuf = MAX( znuuf, 1.e-14 )
-
-            ! cosinus and sinus using scalar and vectorial products
+            !
+            !                       ! cosinus and sinus using dot and cross products
             gsint(ji,jj) = ( zxnpt*zyvvt - zynpt*zxvvt ) / znvvt
             gcost(ji,jj) = ( zxnpt*zxvvt + zynpt*zyvvt ) / znvvt
-
+            !
             gsinu(ji,jj) = ( zxnpu*zyffu - zynpu*zxffu ) / znffu
             gcosu(ji,jj) = ( zxnpu*zxffu + zynpu*zyffu ) / znffu
-
+            !
             gsinf(ji,jj) = ( zxnpf*zyuuf - zynpf*zxuuf ) / znuuf
             gcosf(ji,jj) = ( zxnpf*zxuuf + zynpf*zyuuf ) / znuuf
-
-            ! (caution, rotation of 90 degres)
+            !
             gsinv(ji,jj) = ( zxnpv*zxffv + zynpv*zyffv ) / znffv
-            gcosv(ji,jj) =-( zxnpv*zyffv - zynpv*zxffv ) / znffv
-
+            gcosv(ji,jj) =-( zxnpv*zyffv - zynpv*zxffv ) / znffv     ! (caution, rotation of 90 degres)
+            !
          END DO
       END DO
@@ -318,16 +272,14 @@
       ! Lateral boundary conditions !
       ! =========================== !
-
-      ! lateral boundary cond.: T-, U-, V-, F-pts, sgn
+      !           ! lateral boundary cond.: T-, U-, V-, F-pts, sgn
       CALL lbc_lnk( gcost, 'T', -1. )   ;   CALL lbc_lnk( gsint, 'T', -1. )
       CALL lbc_lnk( gcosu, 'U', -1. )   ;   CALL lbc_lnk( gsinu, 'U', -1. )
       CALL lbc_lnk( gcosv, 'V', -1. )   ;   CALL lbc_lnk( gsinv, 'V', -1. )
       CALL lbc_lnk( gcosf, 'F', -1. )   ;   CALL lbc_lnk( gsinf, 'F', -1. )
-
+      !
    END SUBROUTINE angle
 
 
-   SUBROUTINE geo2oce ( pxx, pyy, pzz, cgrid,     &
-                        pte, ptn )
+   SUBROUTINE geo2oce ( pxx, pyy, pzz, cgrid, pte, ptn )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE geo2oce  ***
@@ -335,19 +287,11 @@
       !! ** Purpose :
       !!
-      !! ** Method  :   Change wind stress from geocentric to east/north
-      !!
-      !! History :
-      !!        !         (O. Marti)  Original code
-      !!        !  91-03  (G. Madec)
-      !!        !  92-07  (M. Imbard)
-      !!        !  99-11  (M. Imbard) NetCDF format with IOIPSL
-      !!        !  00-08  (D. Ludicone) Reduced section at Bab el Mandeb
-      !!   8.5  !  02-06  (G. Madec)  F90: Free form
-      !!   3.0  !  07-08  (G. Madec)  geo2oce suppress lon/lat agruments
+      !! ** Method  :   Change a vector from geocentric to east/north 
+      !!
       !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::  pxx, pyy, pzz
       CHARACTER(len=1)            , INTENT(in   ) ::  cgrid
       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::  pte, ptn
-      !!
+      !
       REAL(wp), PARAMETER :: rpi = 3.141592653e0
       REAL(wp), PARAMETER :: rad = rpi / 180.e0
@@ -355,5 +299,5 @@
       INTEGER ::   ierr   ! local integer
       !!----------------------------------------------------------------------
-
+      !
       IF( .NOT. ALLOCATED( gsinlon ) ) THEN
          ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) ,   &
@@ -361,4 +305,81 @@
          IF( lk_mpp    )   CALL mpp_sum( ierr )
          IF( ierr /= 0 )   CALL ctl_stop('geo2oce: unable to allocate arrays' )
+      ENDIF
+      !
+      SELECT CASE( cgrid)
+      CASE ( 'T' )   
+         ig = 1
+         IF( .NOT. linit(ig) ) THEN 
+            gsinlon(:,:,ig) = SIN( rad * glamt(:,:) )
+            gcoslon(:,:,ig) = COS( rad * glamt(:,:) )
+            gsinlat(:,:,ig) = SIN( rad * gphit(:,:) )
+            gcoslat(:,:,ig) = COS( rad * gphit(:,:) )
+            linit(ig) = .TRUE.
+         ENDIF
+      CASE ( 'U' )   
+         ig = 2
+         IF( .NOT. linit(ig) ) THEN 
+            gsinlon(:,:,ig) = SIN( rad * glamu(:,:) )
+            gcoslon(:,:,ig) = COS( rad * glamu(:,:) )
+            gsinlat(:,:,ig) = SIN( rad * gphiu(:,:) )
+            gcoslat(:,:,ig) = COS( rad * gphiu(:,:) )
+            linit(ig) = .TRUE.
+         ENDIF
+      CASE ( 'V' )   
+         ig = 3
+         IF( .NOT. linit(ig) ) THEN 
+            gsinlon(:,:,ig) = SIN( rad * glamv(:,:) )
+            gcoslon(:,:,ig) = COS( rad * glamv(:,:) )
+            gsinlat(:,:,ig) = SIN( rad * gphiv(:,:) )
+            gcoslat(:,:,ig) = COS( rad * gphiv(:,:) )
+            linit(ig) = .TRUE.
+         ENDIF
+      CASE ( 'F' )   
+         ig = 4
+         IF( .NOT. linit(ig) ) THEN 
+            gsinlon(:,:,ig) = SIN( rad * glamf(:,:) )
+            gcoslon(:,:,ig) = COS( rad * glamf(:,:) )
+            gsinlat(:,:,ig) = SIN( rad * gphif(:,:) )
+            gcoslat(:,:,ig) = COS( rad * gphif(:,:) )
+            linit(ig) = .TRUE.
+         ENDIF
+      CASE default   
+         WRITE(ctmp1,*) 'geo2oce : bad grid argument : ', cgrid
+         CALL ctl_stop( ctmp1 )
+      END SELECT
+      !
+      pte = - gsinlon(:,:,ig) * pxx + gcoslon(:,:,ig) * pyy
+      ptn = - gcoslon(:,:,ig) * gsinlat(:,:,ig) * pxx    &
+         &  - gsinlon(:,:,ig) * gsinlat(:,:,ig) * pyy    &
+         &  + gcoslat(:,:,ig) * pzz
+      !
+   END SUBROUTINE geo2oce
+
+
+   SUBROUTINE oce2geo ( pte, ptn, cgrid, pxx , pyy , pzz )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE oce2geo  ***
+      !!      
+      !! ** Purpose :
+      !!
+      !! ** Method  :   Change vector from east/north to geocentric
+      !!
+      !! History :     ! (A. Caubel)  oce2geo - Original code
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT( IN    ) ::  pte, ptn
+      CHARACTER(len=1)            , INTENT( IN    ) ::  cgrid
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(   OUT ) ::  pxx , pyy , pzz
+      !!
+      REAL(wp), PARAMETER :: rpi = 3.141592653E0
+      REAL(wp), PARAMETER :: rad = rpi / 180.e0
+      INTEGER ::   ig     !
+      INTEGER ::   ierr   ! local integer
+      !!----------------------------------------------------------------------
+
+      IF( .NOT. ALLOCATED( gsinlon ) ) THEN
+         ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) ,   &
+            &      gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr )
+         IF( lk_mpp    )   CALL mpp_sum( ierr )
+         IF( ierr /= 0 )   CALL ctl_stop('oce2geo: unable to allocate arrays' )
       ENDIF
 
@@ -404,139 +425,13 @@
             CALL ctl_stop( ctmp1 )
       END SELECT
-      
-      pte = - gsinlon(:,:,ig) * pxx + gcoslon(:,:,ig) * pyy
-      ptn = - gcoslon(:,:,ig) * gsinlat(:,:,ig) * pxx    &
-            - gsinlon(:,:,ig) * gsinlat(:,:,ig) * pyy    &
-            + gcoslat(:,:,ig) * pzz
-!!$   ptv =   gcoslon(:,:,ig) * gcoslat(:,:,ig) * pxx    &
-!!$         + gsinlon(:,:,ig) * gcoslat(:,:,ig) * pyy    &
-!!$         + gsinlat(:,:,ig) * pzz
-      !
-   END SUBROUTINE geo2oce
-
-   SUBROUTINE oce2geo ( pte, ptn, cgrid,     &
-                        pxx , pyy , pzz )
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE oce2geo  ***
-      !!      
-      !! ** Purpose :
-      !!
-      !! ** Method  :   Change vector from east/north to geocentric
-      !!
-      !! History :
-      !!        !         (A. Caubel)  oce2geo - Original code
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT( IN    ) ::  pte, ptn
-      CHARACTER(len=1)            , INTENT( IN    ) ::  cgrid
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(   OUT ) ::  pxx , pyy , pzz
-      !!
-      REAL(wp), PARAMETER :: rpi = 3.141592653E0
-      REAL(wp), PARAMETER :: rad = rpi / 180.e0
-      INTEGER ::   ig     !
-      INTEGER ::   ierr   ! local integer
-      !!----------------------------------------------------------------------
-
-      IF( .NOT. ALLOCATED( gsinlon ) ) THEN
-         ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) ,   &
-            &      gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr )
-         IF( lk_mpp    )   CALL mpp_sum( ierr )
-         IF( ierr /= 0 )   CALL ctl_stop('oce2geo: unable to allocate arrays' )
-      ENDIF
-
-      SELECT CASE( cgrid)
-         CASE ( 'T' )   
-            ig = 1
-            IF( .NOT. linit(ig) ) THEN 
-               gsinlon(:,:,ig) = SIN( rad * glamt(:,:) )
-               gcoslon(:,:,ig) = COS( rad * glamt(:,:) )
-               gsinlat(:,:,ig) = SIN( rad * gphit(:,:) )
-               gcoslat(:,:,ig) = COS( rad * gphit(:,:) )
-               linit(ig) = .TRUE.
-            ENDIF
-         CASE ( 'U' )   
-            ig = 2
-            IF( .NOT. linit(ig) ) THEN 
-               gsinlon(:,:,ig) = SIN( rad * glamu(:,:) )
-               gcoslon(:,:,ig) = COS( rad * glamu(:,:) )
-               gsinlat(:,:,ig) = SIN( rad * gphiu(:,:) )
-               gcoslat(:,:,ig) = COS( rad * gphiu(:,:) )
-               linit(ig) = .TRUE.
-            ENDIF
-         CASE ( 'V' )   
-            ig = 3
-            IF( .NOT. linit(ig) ) THEN 
-               gsinlon(:,:,ig) = SIN( rad * glamv(:,:) )
-               gcoslon(:,:,ig) = COS( rad * glamv(:,:) )
-               gsinlat(:,:,ig) = SIN( rad * gphiv(:,:) )
-               gcoslat(:,:,ig) = COS( rad * gphiv(:,:) )
-               linit(ig) = .TRUE.
-            ENDIF
-         CASE ( 'F' )   
-            ig = 4
-            IF( .NOT. linit(ig) ) THEN 
-               gsinlon(:,:,ig) = SIN( rad * glamf(:,:) )
-               gcoslon(:,:,ig) = COS( rad * glamf(:,:) )
-               gsinlat(:,:,ig) = SIN( rad * gphif(:,:) )
-               gcoslat(:,:,ig) = COS( rad * gphif(:,:) )
-               linit(ig) = .TRUE.
-            ENDIF
-         CASE default   
-            WRITE(ctmp1,*) 'geo2oce : bad grid argument : ', cgrid
-            CALL ctl_stop( ctmp1 )
-      END SELECT
-
-       pxx = - gsinlon(:,:,ig) * pte - gcoslon(:,:,ig) * gsinlat(:,:,ig) * ptn 
-       pyy =   gcoslon(:,:,ig) * pte - gsinlon(:,:,ig) * gsinlat(:,:,ig) * ptn
-       pzz =   gcoslat(:,:,ig) * ptn
-
-      
+      !
+      pxx = - gsinlon(:,:,ig) * pte - gcoslon(:,:,ig) * gsinlat(:,:,ig) * ptn 
+      pyy =   gcoslon(:,:,ig) * pte - gsinlon(:,:,ig) * gsinlat(:,:,ig) * ptn
+      pzz =   gcoslat(:,:,ig) * ptn
+      !
    END SUBROUTINE oce2geo
 
 
-   SUBROUTINE repere ( px1, py1, px2, py2, kchoix, cd_type )
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE repere  ***
-      !!        
-      !! ** Purpose :   Change vector componantes between a geopgraphic grid 
-      !!      and a stretched coordinates grid.
-      !!
-      !! ** Method  :   
-      !!
-      !! ** Action  :
-      !!
-      !! History :
-      !!        !  89-03  (O. Marti)  original code
-      !!        !  92-02  (M. Imbard)
-      !!        !  93-03  (M. Guyon)  symetrical conditions
-      !!        !  98-05  (B. Blanke)
-      !!   8.5  !  02-08  (G. Madec)  F90: Free form
-      !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   px1, py1   ! two horizontal components to be rotated
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   px2, py2   ! the two horizontal components in the model repere
-      INTEGER , INTENT(in   )                     ::   kchoix     ! type of transformation
-      !                                                           ! = 1 change from geographic to model grid.
-      !                                                           ! =-1 change from model to geographic grid
-      CHARACTER(len=1), INTENT(in   ), OPTIONAL   ::   cd_type    ! define the nature of pt2d array grid-points
-      !
-      CHARACTER(len=1) ::   cl_type      ! define the nature of pt2d array grid-points (T point by default)
-      !!----------------------------------------------------------------------
-
-      cl_type = 'T'
-      IF( PRESENT(cd_type) )   cl_type = cd_type
-         !
-      SELECT CASE (kchoix)
-      CASE ( 1)      ! change from geographic to model grid.
-         CALL rot_rep( px1, py1, cl_type, 'en->i', px2 )
-         CALL rot_rep( px1, py1, cl_type, 'en->j', py2 )
-      CASE (-1)      ! change from model to geographic grid
-         CALL rot_rep( px1, py1, cl_type, 'ij->e', px2 )
-         CALL rot_rep( px1, py1, cl_type, 'ij->n', py2 )
-      CASE DEFAULT   ;   CALL ctl_stop( 'repere: Syntax Error in the definition of kchoix (1 OR -1' )
-      END SELECT
-      
-   END SUBROUTINE repere
-
-
-   SUBROUTINE obs_rot ( psinu, pcosu, psinv, pcosv )
+   SUBROUTINE obs_rot( psinu, pcosu, psinv, pcosv )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE obs_rot  ***
@@ -546,28 +441,24 @@
       !!                current at observation points
       !!
-      !! History :
-      !!   9.2  !  09-02  (K. Mogensen)
+      !! History :  9.2  !  09-02  (K. Mogensen)
       !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj), INTENT( OUT )::   psinu, pcosu, psinv, pcosv   ! copy of data
       !!----------------------------------------------------------------------
-
+      !
       ! Initialization of gsin* and gcos* at first call
       ! -----------------------------------------------
-
       IF( lmust_init ) THEN
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) ' obs_rot : geographic <--> stretched'
          IF(lwp) WRITE(numout,*) ' ~~~~~~~   coordinate transformation'
-
          CALL angle       ! initialization of the transformation
          lmust_init = .FALSE.
-
       ENDIF
-
+      !
       psinu(:,:) = gsinu(:,:)
       pcosu(:,:) = gcosu(:,:)
       psinv(:,:) = gsinv(:,:)
       pcosv(:,:) = gcosv(:,:)
-
+      !
    END SUBROUTINE obs_rot
 
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 6121)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 6140)
@@ -44,8 +44,8 @@
    LOGICAL , PUBLIC ::   ln_dm2dc       !: Daily mean to Diurnal Cycle short wave (qsr)
    LOGICAL , PUBLIC ::   ln_rnf         !: runoffs / runoff mouths
+   LOGICAL , PUBLIC ::   ln_isf         !: ice shelf melting
    LOGICAL , PUBLIC ::   ln_ssr         !: Sea Surface restoring on SST and/or SSS      
    LOGICAL , PUBLIC ::   ln_apr_dyn     !: Atmospheric pressure forcing used on dynamics (ocean & ice)
    INTEGER , PUBLIC ::   nn_ice         !: flag for ice in the surface boundary condition (=0/1/2/3)
-   INTEGER , PUBLIC ::   nn_isf         !: flag for isf in the surface boundary condition (=0/1/2/3/4) 
    INTEGER , PUBLIC ::   nn_ice_embd    !: flag for levitating/embedding sea-ice in the ocean
    !                                             !: =0 levitating ice (no mass exchange, concentration/dilution effect)
@@ -172,7 +172,5 @@
          &      ssv_m  (jpi,jpj) , sss_m(jpi,jpj) , ssh_m(jpi,jpj) , STAT=ierr(4) )
          !
-#if defined key_vvl
       ALLOCATE( e3t_m(jpi,jpj) , STAT=ierr(5) )
-#endif
          !
       sbc_oce_alloc = MAXVAL( ierr )
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcana.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcana.F90	(revision 6121)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcana.F90	(revision 6140)
@@ -37,5 +37,4 @@
    
    !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcapr.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcapr.F90	(revision 6121)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcapr.F90	(revision 6140)
@@ -13,4 +13,5 @@
    USE sbc_oce         ! surface boundary condition
    USE phycst          ! physical constants
+   !
    USE fldread         ! read input fields
    USE in_out_manager  ! I/O manager
@@ -38,6 +39,4 @@
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_apr   ! structure of input fields (file informations, fields read)
 
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
@@ -112,5 +111,5 @@
          ENDIF
 !jc: stop below should rather be a warning 
-         IF( ( ln_apr_obc ) .AND. .NOT. ln_apr_dyn   )   &
+         IF( ln_apr_obc .AND. .NOT.ln_apr_dyn   )   &
             CALL ctl_stop( 'sbc_apr: use inverse barometer ssh at open boundary ONLY requires ln_apr_dyn=T' )
       ENDIF
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 6121)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 6140)
@@ -19,36 +19,37 @@
 
    !!----------------------------------------------------------------------
-   !!   sbc_blk_core    : bulk formulation as ocean surface boundary condition (forced mode, CORE bulk formulea)
-   !!   blk_oce_core    : computes momentum, heat and freshwater fluxes over ocean
-   !!   blk_ice_core    : computes momentum, heat and freshwater fluxes over ice
-   !!   turb_core_2z    : Computes turbulent transfert coefficients
-   !!   cd_neutral_10m  : Estimate of the neutral drag coefficient at 10m
-   !!   psi_m           : universal profile stability function for momentum
-   !!   psi_h           : universal profile stability function for temperature and humidity
+   !!   sbc_blk_core  : bulk formulation as ocean surface boundary condition (forced mode, CORE bulk formulea)
+   !!   blk_oce_core  : computes momentum, heat and freshwater fluxes over ocean
+   !!   blk_ice_core  : computes momentum, heat and freshwater fluxes over ice
+   !!   turb_core_2z  : Computes turbulent transfert coefficients
+   !!   cd_neutral_10m: Estimate of the neutral drag coefficient at 10m
+   !!   psi_m         : universal profile stability function for momentum
+   !!   psi_h         : universal profile stability function for temperature and humidity
    !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers
-   USE dom_oce         ! ocean space and time domain
-   USE phycst          ! physical constants
-   USE fldread         ! read input fields
-   USE sbc_oce         ! Surface boundary condition: ocean fields
-   USE cyclone         ! Cyclone 10m wind form trac of cyclone centres
-   USE sbcdcy          ! surface boundary condition: diurnal cycle
-   USE iom             ! I/O manager library
-   USE in_out_manager  ! I/O manager
-   USE lib_mpp         ! distribued memory computing library
-   USE wrk_nemo        ! work arrays
-   USE timing          ! Timing
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE prtctl          ! Print control
-   USE sbcwave, ONLY   :  cdn_wave ! wave module
-   USE sbc_ice         ! Surface boundary condition: ice fields
-   USE lib_fortran     ! to use key_nosignedzero
+   USE oce            ! ocean dynamics and tracers
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   USE fldread        ! read input fields
+   USE sbc_oce        ! Surface boundary condition: ocean fields
+   USE cyclone        ! Cyclone 10m wind form trac of cyclone centres
+   USE sbcdcy         ! surface boundary condition: diurnal cycle
+   USE sbcwave , ONLY :   cdn_wave ! wave module
+   USE sbc_ice        ! Surface boundary condition: ice fields
+   USE lib_fortran    ! to use key_nosignedzero
 #if defined key_lim3
-   USE ice, ONLY       : u_ice, v_ice, jpl, pfrld, a_i_b
-   USE limthd_dh       ! for CALL lim_thd_snwblow
+   USE ice     , ONLY :   u_ice, v_ice, jpl, pfrld, a_i_b
+   USE limthd_dh      ! for CALL lim_thd_snwblow
 #elif defined key_lim2
-   USE ice_2, ONLY     : u_ice, v_ice
-   USE par_ice_2
+   USE ice_2   , ONLY :   u_ice, v_ice
+   USE par_ice_2      ! LIM-2 parameters
 #endif
+   !
+   USE iom            ! I/O manager library
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! distribued memory computing library
+   USE wrk_nemo       ! work arrays
+   USE timing         ! Timing
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE prtctl         ! Print control
 
    IMPLICIT NONE
@@ -84,5 +85,5 @@
    REAL(wp), PARAMETER ::   albo =    0.066       ! ocean albedo assumed to be constant
 
-   !                                  !!* Namelist namsbc_core : CORE bulk parameters
+   !                        !!* Namelist namsbc_core : CORE bulk parameters
    LOGICAL  ::   ln_taudif   ! logical flag to use the "mean of stress module - module of mean stress" data
    REAL(wp) ::   rn_pfac     ! multiplication factor for precipitation
@@ -93,5 +94,4 @@
 
    !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -149,7 +149,7 @@
       TYPE(FLD_N) ::   sn_tdif                                 !   "                                 "
       NAMELIST/namsbc_core/ cn_dir , ln_taudif, rn_pfac, rn_efac, rn_vfac,  &
-         &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr ,           &
-         &                  sn_qlw , sn_tair, sn_prec  , sn_snow,           &
-         &                  sn_tdif, rn_zqt,  rn_zu
+         &                  sn_wndi, sn_wndj  , sn_humi, sn_qsr ,           &
+         &                  sn_qlw , sn_tair  , sn_prec, sn_snow,           &
+         &                  sn_tdif, rn_zqt   ,  rn_zu
       !!---------------------------------------------------------------------
       !
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90	(revision 6121)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90	(revision 6140)
@@ -8,21 +8,22 @@
 
    !!----------------------------------------------------------------------
-   !!   sbc_blk_mfs  : bulk formulation as ocean surface boundary condition
+   !!   sbc_blk_mfs   : bulk formulation as ocean surface boundary condition
    !!                   (forced mode, mfs bulk formulae)
-   !!   blk_oce_mfs  : ocean: computes momentum, heat and freshwater fluxes
+   !!   blk_oce_mfs   : ocean: computes momentum, heat and freshwater fluxes
    !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers
-   USE dom_oce         ! ocean space and time domain
-   USE phycst          ! physical constants
-   USE fldread         ! read input fields
-   USE sbc_oce         ! Surface boundary condition: ocean fields
-   USE iom             ! I/O manager library
-   USE in_out_manager  ! I/O manager
-   USE lib_mpp         ! distribued memory computing library
-   USE wrk_nemo        ! work arrays
-   USE timing          ! Timing
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE prtctl          ! Print control
-   USE sbcwave,ONLY : cdn_wave !wave module
+   USE oce            ! ocean dynamics and tracers
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   USE fldread        ! read input fields
+   USE sbc_oce        ! Surface boundary condition: ocean fields
+   USE sbcwave  ,ONLY :   cdn_wave !wave module
+   !
+   USE iom            ! I/O manager library
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! distribued memory computing library
+   USE wrk_nemo       ! work arrays
+   USE timing         ! Timing
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE prtctl         ! Print control
 
    IMPLICIT NONE
@@ -42,5 +43,4 @@
          
    !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -49,7 +49,5 @@
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
-
 
    SUBROUTINE sbc_blk_mfs( kt )
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 6121)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 6140)
@@ -18,29 +18,23 @@
    !!   sbc_cpl_snd     : send     fields to the atmosphere
    !!----------------------------------------------------------------------
-   USE dom_oce         ! ocean space and time domain
-   USE sbc_oce         ! Surface boundary condition: ocean fields
-   USE sbc_ice         ! Surface boundary condition: ice fields
-   USE sbcapr
-   USE sbcdcy          ! surface boundary condition: diurnal cycle
-   USE phycst          ! physical constants
+   USE dom_oce        ! ocean space and time domain
+   USE sbc_oce        ! Surface boundary condition: ocean fields
+   USE sbc_ice        ! Surface boundary condition: ice fields
+   USE sbcapr         ! Stochastic param. : ???
+   USE sbcdcy         ! surface boundary condition: diurnal cycle
+   USE phycst         ! physical constants
 #if defined key_lim3
-   USE ice             ! ice variables
+   USE ice            ! ice variables
 #endif
 #if defined key_lim2
-   USE par_ice_2       ! ice parameters
-   USE ice_2           ! ice variables
+   USE par_ice_2      ! ice parameters
+   USE ice_2          ! ice variables
 #endif
-   USE cpl_oasis3      ! OASIS3 coupling
-   USE geo2ocean       ! 
+   USE cpl_oasis3     ! OASIS3 coupling
+   USE geo2ocean      ! 
    USE oce   , ONLY : tsn, un, vn, sshn, ub, vb, sshb, fraqsr_1lev
-   USE albedo          !
-   USE in_out_manager  ! I/O manager
-   USE iom             ! NetCDF library
-   USE lib_mpp         ! distribued memory computing library
-   USE wrk_nemo        ! work arrays
-   USE timing          ! Timing
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE eosbn2
-   USE sbcrnf   , ONLY : l_rnfcpl
+   USE albedo         ! 
+   USE eosbn2         ! 
+   USE sbcrnf  , ONLY : l_rnfcpl
 #if defined key_cpl_carbon_cycle
    USE p4zflx, ONLY : oce_co2
@@ -50,127 +44,132 @@
 #endif
 #if defined key_lim3
-   USE limthd_dh       ! for CALL lim_thd_snwblow
+   USE limthd_dh      ! for CALL lim_thd_snwblow
 #endif
+   !
+   USE in_out_manager ! I/O manager
+   USE iom            ! NetCDF library
+   USE lib_mpp        ! distribued memory computing library
+   USE wrk_nemo       ! work arrays
+   USE timing         ! Timing
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   sbc_cpl_init       ! routine called by sbcmod.F90
-   PUBLIC   sbc_cpl_rcv        ! routine called by sbc_ice_lim(_2).F90
-   PUBLIC   sbc_cpl_snd        ! routine called by step.F90
-   PUBLIC   sbc_cpl_ice_tau    ! routine called by sbc_ice_lim(_2).F90
-   PUBLIC   sbc_cpl_ice_flx    ! routine called by sbc_ice_lim(_2).F90
-   PUBLIC   sbc_cpl_alloc      ! routine called in sbcice_cice.F90
-
-   INTEGER, PARAMETER ::   jpr_otx1   =  1            ! 3 atmosphere-ocean stress components on grid 1
-   INTEGER, PARAMETER ::   jpr_oty1   =  2            ! 
-   INTEGER, PARAMETER ::   jpr_otz1   =  3            ! 
-   INTEGER, PARAMETER ::   jpr_otx2   =  4            ! 3 atmosphere-ocean stress components on grid 2
-   INTEGER, PARAMETER ::   jpr_oty2   =  5            ! 
-   INTEGER, PARAMETER ::   jpr_otz2   =  6            ! 
-   INTEGER, PARAMETER ::   jpr_itx1   =  7            ! 3 atmosphere-ice   stress components on grid 1
-   INTEGER, PARAMETER ::   jpr_ity1   =  8            ! 
-   INTEGER, PARAMETER ::   jpr_itz1   =  9            ! 
-   INTEGER, PARAMETER ::   jpr_itx2   = 10            ! 3 atmosphere-ice   stress components on grid 2
-   INTEGER, PARAMETER ::   jpr_ity2   = 11            ! 
-   INTEGER, PARAMETER ::   jpr_itz2   = 12            ! 
-   INTEGER, PARAMETER ::   jpr_qsroce = 13            ! Qsr above the ocean
-   INTEGER, PARAMETER ::   jpr_qsrice = 14            ! Qsr above the ice
+   PUBLIC   sbc_cpl_init      ! routine called by sbcmod.F90
+   PUBLIC   sbc_cpl_rcv       ! routine called by sbc_ice_lim(_2).F90
+   PUBLIC   sbc_cpl_snd       ! routine called by step.F90
+   PUBLIC   sbc_cpl_ice_tau   ! routine called by sbc_ice_lim(_2).F90
+   PUBLIC   sbc_cpl_ice_flx   ! routine called by sbc_ice_lim(_2).F90
+   PUBLIC   sbc_cpl_alloc     ! routine called in sbcice_cice.F90
+
+   INTEGER, PARAMETER ::   jpr_otx1   =  1   ! 3 atmosphere-ocean stress components on grid 1
+   INTEGER, PARAMETER ::   jpr_oty1   =  2   ! 
+   INTEGER, PARAMETER ::   jpr_otz1   =  3   ! 
+   INTEGER, PARAMETER ::   jpr_otx2   =  4   ! 3 atmosphere-ocean stress components on grid 2
+   INTEGER, PARAMETER ::   jpr_oty2   =  5   ! 
+   INTEGER, PARAMETER ::   jpr_otz2   =  6   ! 
+   INTEGER, PARAMETER ::   jpr_itx1   =  7   ! 3 atmosphere-ice   stress components on grid 1
+   INTEGER, PARAMETER ::   jpr_ity1   =  8   ! 
+   INTEGER, PARAMETER ::   jpr_itz1   =  9   ! 
+   INTEGER, PARAMETER ::   jpr_itx2   = 10   ! 3 atmosphere-ice   stress components on grid 2
+   INTEGER, PARAMETER ::   jpr_ity2   = 11   ! 
+   INTEGER, PARAMETER ::   jpr_itz2   = 12   ! 
+   INTEGER, PARAMETER ::   jpr_qsroce = 13   ! Qsr above the ocean
+   INTEGER, PARAMETER ::   jpr_qsrice = 14   ! Qsr above the ice
    INTEGER, PARAMETER ::   jpr_qsrmix = 15 
-   INTEGER, PARAMETER ::   jpr_qnsoce = 16            ! Qns above the ocean
-   INTEGER, PARAMETER ::   jpr_qnsice = 17            ! Qns above the ice
+   INTEGER, PARAMETER ::   jpr_qnsoce = 16   ! Qns above the ocean
+   INTEGER, PARAMETER ::   jpr_qnsice = 17   ! Qns above the ice
    INTEGER, PARAMETER ::   jpr_qnsmix = 18
-   INTEGER, PARAMETER ::   jpr_rain   = 19            ! total liquid precipitation (rain)
-   INTEGER, PARAMETER ::   jpr_snow   = 20            ! solid precipitation over the ocean (snow)
-   INTEGER, PARAMETER ::   jpr_tevp   = 21            ! total evaporation
-   INTEGER, PARAMETER ::   jpr_ievp   = 22            ! solid evaporation (sublimation)
-   INTEGER, PARAMETER ::   jpr_sbpr   = 23            ! sublimation - liquid precipitation - solid precipitation
-   INTEGER, PARAMETER ::   jpr_semp   = 24            ! solid freshwater budget (sublimation - snow)
-   INTEGER, PARAMETER ::   jpr_oemp   = 25            ! ocean freshwater budget (evap - precip)
-   INTEGER, PARAMETER ::   jpr_w10m   = 26            ! 10m wind
-   INTEGER, PARAMETER ::   jpr_dqnsdt = 27            ! d(Q non solar)/d(temperature)
-   INTEGER, PARAMETER ::   jpr_rnf    = 28            ! runoffs
-   INTEGER, PARAMETER ::   jpr_cal    = 29            ! calving
-   INTEGER, PARAMETER ::   jpr_taum   = 30            ! wind stress module
+   INTEGER, PARAMETER ::   jpr_rain   = 19   ! total liquid precipitation (rain)
+   INTEGER, PARAMETER ::   jpr_snow   = 20   ! solid precipitation over the ocean (snow)
+   INTEGER, PARAMETER ::   jpr_tevp   = 21   ! total evaporation
+   INTEGER, PARAMETER ::   jpr_ievp   = 22   ! solid evaporation (sublimation)
+   INTEGER, PARAMETER ::   jpr_sbpr   = 23   ! sublimation - liquid precipitation - solid precipitation
+   INTEGER, PARAMETER ::   jpr_semp   = 24   ! solid freshwater budget (sublimation - snow)
+   INTEGER, PARAMETER ::   jpr_oemp   = 25   ! ocean freshwater budget (evap - precip)
+   INTEGER, PARAMETER ::   jpr_w10m   = 26   ! 10m wind
+   INTEGER, PARAMETER ::   jpr_dqnsdt = 27   ! d(Q non solar)/d(temperature)
+   INTEGER, PARAMETER ::   jpr_rnf    = 28   ! runoffs
+   INTEGER, PARAMETER ::   jpr_cal    = 29   ! calving
+   INTEGER, PARAMETER ::   jpr_taum   = 30   ! wind stress module
    INTEGER, PARAMETER ::   jpr_co2    = 31
-   INTEGER, PARAMETER ::   jpr_topm   = 32            ! topmeltn
-   INTEGER, PARAMETER ::   jpr_botm   = 33            ! botmeltn
-   INTEGER, PARAMETER ::   jpr_sflx   = 34            ! salt flux
-   INTEGER, PARAMETER ::   jpr_toce   = 35            ! ocean temperature
-   INTEGER, PARAMETER ::   jpr_soce   = 36            ! ocean salinity
-   INTEGER, PARAMETER ::   jpr_ocx1   = 37            ! ocean current on grid 1
-   INTEGER, PARAMETER ::   jpr_ocy1   = 38            !
-   INTEGER, PARAMETER ::   jpr_ssh    = 39            ! sea surface height
-   INTEGER, PARAMETER ::   jpr_fice   = 40            ! ice fraction          
-   INTEGER, PARAMETER ::   jpr_e3t1st = 41            ! first T level thickness 
-   INTEGER, PARAMETER ::   jpr_fraqsr = 42            ! fraction of solar net radiation absorbed in the first ocean level
-   INTEGER, PARAMETER ::   jprcv      = 42            ! total number of fields received
-
-   INTEGER, PARAMETER ::   jps_fice   =  1            ! ice fraction sent to the atmosphere
-   INTEGER, PARAMETER ::   jps_toce   =  2            ! ocean temperature
-   INTEGER, PARAMETER ::   jps_tice   =  3            ! ice   temperature
-   INTEGER, PARAMETER ::   jps_tmix   =  4            ! mixed temperature (ocean+ice)
-   INTEGER, PARAMETER ::   jps_albice =  5            ! ice   albedo
-   INTEGER, PARAMETER ::   jps_albmix =  6            ! mixed albedo
-   INTEGER, PARAMETER ::   jps_hice   =  7            ! ice  thickness
-   INTEGER, PARAMETER ::   jps_hsnw   =  8            ! snow thickness
-   INTEGER, PARAMETER ::   jps_ocx1   =  9            ! ocean current on grid 1
-   INTEGER, PARAMETER ::   jps_ocy1   = 10            !
-   INTEGER, PARAMETER ::   jps_ocz1   = 11            !
-   INTEGER, PARAMETER ::   jps_ivx1   = 12            ! ice   current on grid 1
-   INTEGER, PARAMETER ::   jps_ivy1   = 13            !
-   INTEGER, PARAMETER ::   jps_ivz1   = 14            !
+   INTEGER, PARAMETER ::   jpr_topm   = 32   ! topmeltn
+   INTEGER, PARAMETER ::   jpr_botm   = 33   ! botmeltn
+   INTEGER, PARAMETER ::   jpr_sflx   = 34   ! salt flux
+   INTEGER, PARAMETER ::   jpr_toce   = 35   ! ocean temperature
+   INTEGER, PARAMETER ::   jpr_soce   = 36   ! ocean salinity
+   INTEGER, PARAMETER ::   jpr_ocx1   = 37   ! ocean current on grid 1
+   INTEGER, PARAMETER ::   jpr_ocy1   = 38   !
+   INTEGER, PARAMETER ::   jpr_ssh    = 39   ! sea surface height
+   INTEGER, PARAMETER ::   jpr_fice   = 40   ! ice fraction          
+   INTEGER, PARAMETER ::   jpr_e3t1st = 41   ! first T level thickness 
+   INTEGER, PARAMETER ::   jpr_fraqsr = 42   ! fraction of solar net radiation absorbed in the first ocean level
+   INTEGER, PARAMETER ::   jprcv      = 42   ! total number of fields received
+
+   INTEGER, PARAMETER ::   jps_fice   =  1   ! ice fraction sent to the atmosphere
+   INTEGER, PARAMETER ::   jps_toce   =  2   ! ocean temperature
+   INTEGER, PARAMETER ::   jps_tice   =  3   ! ice   temperature
+   INTEGER, PARAMETER ::   jps_tmix   =  4   ! mixed temperature (ocean+ice)
+   INTEGER, PARAMETER ::   jps_albice =  5   ! ice   albedo
+   INTEGER, PARAMETER ::   jps_albmix =  6   ! mixed albedo
+   INTEGER, PARAMETER ::   jps_hice   =  7   ! ice  thickness
+   INTEGER, PARAMETER ::   jps_hsnw   =  8   ! snow thickness
+   INTEGER, PARAMETER ::   jps_ocx1   =  9   ! ocean current on grid 1
+   INTEGER, PARAMETER ::   jps_ocy1   = 10   !
+   INTEGER, PARAMETER ::   jps_ocz1   = 11   !
+   INTEGER, PARAMETER ::   jps_ivx1   = 12   ! ice   current on grid 1
+   INTEGER, PARAMETER ::   jps_ivy1   = 13   !
+   INTEGER, PARAMETER ::   jps_ivz1   = 14   !
    INTEGER, PARAMETER ::   jps_co2    = 15
-   INTEGER, PARAMETER ::   jps_soce   = 16            ! ocean salinity
-   INTEGER, PARAMETER ::   jps_ssh    = 17            ! sea surface height
-   INTEGER, PARAMETER ::   jps_qsroce = 18            ! Qsr above the ocean
-   INTEGER, PARAMETER ::   jps_qnsoce = 19            ! Qns above the ocean
-   INTEGER, PARAMETER ::   jps_oemp   = 20            ! ocean freshwater budget (evap - precip)
-   INTEGER, PARAMETER ::   jps_sflx   = 21            ! salt flux
-   INTEGER, PARAMETER ::   jps_otx1   = 22            ! 2 atmosphere-ocean stress components on grid 1
-   INTEGER, PARAMETER ::   jps_oty1   = 23            ! 
-   INTEGER, PARAMETER ::   jps_rnf    = 24            ! runoffs
-   INTEGER, PARAMETER ::   jps_taum   = 25            ! wind stress module
-   INTEGER, PARAMETER ::   jps_fice2  = 26            ! ice fraction sent to OPA (by SAS when doing SAS-OPA coupling)
-   INTEGER, PARAMETER ::   jps_e3t1st = 27            ! first level depth (vvl)
-   INTEGER, PARAMETER ::   jps_fraqsr = 28            ! fraction of solar net radiation absorbed in the first ocean level
-   INTEGER, PARAMETER ::   jpsnd      = 28            ! total number of fields sended
-
-   !                                                         !!** namelist namsbc_cpl **
-   TYPE ::   FLD_C
-      CHARACTER(len = 32) ::   cldes                  ! desciption of the coupling strategy
-      CHARACTER(len = 32) ::   clcat                  ! multiple ice categories strategy
-      CHARACTER(len = 32) ::   clvref                 ! reference of vector ('spherical' or 'cartesian')
-      CHARACTER(len = 32) ::   clvor                  ! orientation of vector fields ('eastward-northward' or 'local grid')
-      CHARACTER(len = 32) ::   clvgrd                 ! grids on which is located the vector fields
+   INTEGER, PARAMETER ::   jps_soce   = 16   ! ocean salinity
+   INTEGER, PARAMETER ::   jps_ssh    = 17   ! sea surface height
+   INTEGER, PARAMETER ::   jps_qsroce = 18   ! Qsr above the ocean
+   INTEGER, PARAMETER ::   jps_qnsoce = 19   ! Qns above the ocean
+   INTEGER, PARAMETER ::   jps_oemp   = 20   ! ocean freshwater budget (evap - precip)
+   INTEGER, PARAMETER ::   jps_sflx   = 21   ! salt flux
+   INTEGER, PARAMETER ::   jps_otx1   = 22   ! 2 atmosphere-ocean stress components on grid 1
+   INTEGER, PARAMETER ::   jps_oty1   = 23   ! 
+   INTEGER, PARAMETER ::   jps_rnf    = 24   ! runoffs
+   INTEGER, PARAMETER ::   jps_taum   = 25   ! wind stress module
+   INTEGER, PARAMETER ::   jps_fice2  = 26   ! ice fraction sent to OPA (by SAS when doing SAS-OPA coupling)
+   INTEGER, PARAMETER ::   jps_e3t1st = 27   ! first level depth (vvl)
+   INTEGER, PARAMETER ::   jps_fraqsr = 28   ! fraction of solar net radiation absorbed in the first ocean level
+   INTEGER, PARAMETER ::   jpsnd      = 28   ! total number of fields sended
+
+   !                                  !!** namelist namsbc_cpl **
+   TYPE ::   FLD_C                     !   
+      CHARACTER(len = 32) ::   cldes      ! desciption of the coupling strategy
+      CHARACTER(len = 32) ::   clcat      ! multiple ice categories strategy
+      CHARACTER(len = 32) ::   clvref     ! reference of vector ('spherical' or 'cartesian')
+      CHARACTER(len = 32) ::   clvor      ! orientation of vector fields ('eastward-northward' or 'local grid')
+      CHARACTER(len = 32) ::   clvgrd     ! grids on which is located the vector fields
    END TYPE FLD_C
-   ! Send to the atmosphere                           !
+   !                                   ! Send to the atmosphere  
    TYPE(FLD_C) ::   sn_snd_temp, sn_snd_alb, sn_snd_thick, sn_snd_crt, sn_snd_co2                        
-   ! Received from the atmosphere                     !
+   !                                   ! Received from the atmosphere
    TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr, sn_rcv_qns, sn_rcv_emp, sn_rcv_rnf
    TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2                        
-   ! Other namelist parameters                        !
-   INTEGER     ::   nn_cplmodel            ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
-   LOGICAL     ::   ln_usecplmask          !  use a coupling mask file to merge data received from several models
-                                           !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
+   !                                   ! Other namelist parameters
+   INTEGER     ::   nn_cplmodel           ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
+   LOGICAL     ::   ln_usecplmask         !  use a coupling mask file to merge data received from several models
+                                         !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
    TYPE ::   DYNARR     
       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z3   
    END TYPE DYNARR
 
-   TYPE( DYNARR ), SAVE, DIMENSION(jprcv) ::   frcv                      ! all fields recieved from the atmosphere
-
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   albedo_oce_mix     ! ocean albedo sent to atmosphere (mix clear/overcast sky)
-
-   INTEGER , ALLOCATABLE, SAVE, DIMENSION(    :) ::   nrcvinfo           ! OASIS info argument
+   TYPE( DYNARR ), SAVE, DIMENSION(jprcv) ::   frcv                     ! all fields recieved from the atmosphere
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   albedo_oce_mix    ! ocean albedo sent to atmosphere (mix clear/overcast sky)
+
+   INTEGER , ALLOCATABLE, SAVE, DIMENSION(    :) ::   nrcvinfo          ! OASIS info argument
 
    !! Substitution
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
   
@@ -209,9 +208,8 @@
       !!              * initialise the OASIS coupler
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   k_ice       ! ice management in the sbc (=0/1/2/3)
-      !!
-      INTEGER ::   jn   ! dummy loop index
-      INTEGER ::   ios  ! Local integer output status for namelist read
-      INTEGER ::   inum 
+      INTEGER, INTENT(in) ::   k_ice   ! ice management in the sbc (=0/1/2/3)
+      !
+      INTEGER ::   jn          ! dummy loop index
+      INTEGER ::   ios, inum   ! Local integer
       REAL(wp), POINTER, DIMENSION(:,:) ::   zacs, zaos
       !!
@@ -222,21 +220,21 @@
       !!---------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_init')
-      !
-      CALL wrk_alloc( jpi,jpj, zacs, zaos )
+      IF( nn_timing == 1 )   CALL timing_start('sbc_cpl_init')
+      !
+      CALL wrk_alloc( jpi,jpj,   zacs, zaos )
 
       ! ================================ !
       !      Namelist informations       !
       ! ================================ !
-
+      !
       REWIND( numnam_ref )              ! Namelist namsbc_cpl in reference namelist : Variables for OASIS coupling
       READ  ( numnam_ref, namsbc_cpl, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cpl in reference namelist', lwp )
-
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_cpl in reference namelist', lwp )
+      !
       REWIND( numnam_cfg )              ! Namelist namsbc_cpl in configuration namelist : Variables for OASIS coupling
       READ  ( numnam_cfg, namsbc_cpl, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cpl in configuration namelist', lwp )
+902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_cpl in configuration namelist', lwp )
       IF(lwm) WRITE ( numond, namsbc_cpl )
-
+      !
       IF(lwp) THEN                        ! control print
          WRITE(numout,*)
@@ -374,5 +372,5 @@
          srcv(jpr_ity1)%clgrid = 'V'                  ! i.e. it is always at U- & V-points for i- & j-comp. resp.
       ENDIF
-       
+      !
       !                                                      ! ------------------------- !
       !                                                      !    freshwater budget      !   E-P
@@ -396,5 +394,5 @@
       CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_emp%cldes' )
       END SELECT
-
+      !
       !                                                      ! ------------------------- !
       !                                                      !     Runoffs & Calving     !   
@@ -410,5 +408,5 @@
       !
       srcv(jpr_cal   )%clname = 'OCalving'   ;   IF( TRIM( sn_rcv_cal%cldes ) == 'coupled' )   srcv(jpr_cal)%laction = .TRUE.
-
+      !
       !                                                      ! ------------------------- !
       !                                                      !    non solar radiation    !   Qns
@@ -535,5 +533,5 @@
          IF( .NOT. ln_cpl ) srcv(:)%nsgn    = 1.        ! force default definition in case of opa <-> sas coupling
          srcv( (/jpr_toce, jpr_soce, jpr_ssh, jpr_fraqsr, jpr_ocx1, jpr_ocy1/) )%laction = .TRUE.
-         srcv( jpr_e3t1st )%laction = lk_vvl
+         srcv( jpr_e3t1st )%laction = .NOT.ln_linssh
          srcv(jpr_ocx1)%clgrid = 'U'        ! oce components given at U-point
          srcv(jpr_ocy1)%clgrid = 'V'        !           and           V-point
@@ -701,5 +699,5 @@
          ssnd(:)%laction = .FALSE.   ! force default definition in case of opa <-> sas coupling
          ssnd( (/jps_toce, jps_soce, jps_ssh, jps_fraqsr, jps_ocx1, jps_ocy1/) )%laction = .TRUE.
-         ssnd( jps_e3t1st )%laction = lk_vvl
+         ssnd( jps_e3t1st )%laction = .NOT.ln_linssh
          ! vector definition: not used but cleaner...
          ssnd(jps_ocx1)%clgrid  = 'U'        ! oce components given at U-point
@@ -785,7 +783,7 @@
       ncpl_qsr_freq = 86400 / ncpl_qsr_freq
 
-      CALL wrk_dealloc( jpi,jpj, zacs, zaos )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_init')
+      CALL wrk_dealloc( jpi,jpj,   zacs, zaos )
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('sbc_cpl_init')
       !
    END SUBROUTINE sbc_cpl_init
@@ -837,12 +835,12 @@
       !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case)
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in)           ::   kt          ! ocean model time step index
-      INTEGER, INTENT(in)           ::   k_fsbc      ! frequency of sbc (-> ice model) computation 
-      INTEGER, INTENT(in)           ::   k_ice       ! ice management in the sbc (=0/1/2/3)
+      INTEGER, INTENT(in) ::   kt       ! ocean model time step index
+      INTEGER, INTENT(in) ::   k_fsbc   ! frequency of sbc (-> ice model) computation 
+      INTEGER, INTENT(in) ::   k_ice    ! ice management in the sbc (=0/1/2/3)
 
       !!
       LOGICAL  ::   llnewtx, llnewtau      ! update wind stress components and module??
       INTEGER  ::   ji, jj, jn             ! dummy loop indices
-      INTEGER  ::   isec                   ! number of seconds since nit000 (assuming rdttra did not change since nit000)
+      INTEGER  ::   isec                   ! number of seconds since nit000 (assuming rdt did not change since nit000)
       REAL(wp) ::   zcumulneg, zcumulpos   ! temporary scalars     
       REAL(wp) ::   zcoef                  ! temporary scalar
@@ -853,7 +851,7 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_rcv')
-      !
-      CALL wrk_alloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr )
+      IF( nn_timing == 1 )   CALL timing_start('sbc_cpl_rcv')
+      !
+      CALL wrk_alloc( jpi,jpj,   ztx, zty, zmsk, zemp, zqns, zqsr )
       !
       IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
@@ -862,5 +860,5 @@
       !                                                      ! Receive all the atmos. fields (including ice information)
       !                                                      ! ======================================================= !
-      isec = ( kt - nit000 ) * NINT( rdttra(1) )                ! date of exchanges
+      isec = ( kt - nit000 ) * NINT( rdt )                      ! date of exchanges
       DO jn = 1, jprcv                                          ! received fields sent by the atmosphere
          IF( srcv(jn)%laction )   CALL cpl_rcv( jn, isec, frcv(jn)%z3, xcplmask(:,:,1:nn_cplmodel), nrcvinfo(jn) )
@@ -1104,10 +1102,9 @@
          IF( srcv(jpr_fice )%laction )   fr_i(:,:) = frcv(jpr_fice )%z3(:,:,1)
          !
-
-      ENDIF
-      !
-      CALL wrk_dealloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_rcv')
+      ENDIF
+      !
+      CALL wrk_dealloc( jpi,jpj,   ztx, zty, zmsk, zemp, zqns, zqsr )
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('sbc_cpl_rcv')
       !
    END SUBROUTINE sbc_cpl_rcv
@@ -1150,12 +1147,12 @@
       REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_tauj   ! at I-point (B-grid) or U & V-point (C-grid)
       !!
-      INTEGER ::   ji, jj                          ! dummy loop indices
-      INTEGER ::   itx                             ! index of taux over ice
+      INTEGER ::   ji, jj   ! dummy loop indices
+      INTEGER ::   itx      ! index of taux over ice
       REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty 
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_tau')
-      !
-      CALL wrk_alloc( jpi,jpj, ztx, zty )
+      IF( nn_timing == 1 )   CALL timing_start('sbc_cpl_ice_tau')
+      !
+      CALL wrk_alloc( jpi,jpj,   ztx, zty )
 
       IF( srcv(jpr_itx1)%laction ) THEN   ;   itx =  jpr_itx1   
@@ -1165,5 +1162,4 @@
       ! do something only if we just received the stress from atmosphere
       IF(  nrcvinfo(itx) == OASIS_Rcv ) THEN
-
          !                                                      ! ======================= !
          IF( srcv(jpr_itx1)%laction ) THEN                      !   ice stress received   !
@@ -1318,7 +1314,7 @@
       ENDIF
       !   
-      CALL wrk_dealloc( jpi,jpj, ztx, zty )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_ice_tau')
+      CALL wrk_dealloc( jpi,jpj,   ztx, zty )
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('sbc_cpl_ice_tau')
       !
    END SUBROUTINE sbc_cpl_ice_tau
@@ -1365,11 +1361,11 @@
       !!                   sprecip             solid precipitation over the ocean  
       !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in   ), DIMENSION(:,:)   ::   p_frld     ! lead fraction                [0 to 1]
+      REAL(wp), INTENT(in   ), DIMENSION(:,:)             ::   p_frld  ! lead fraction            [0 to 1]
       ! optional arguments, used only in 'mixed oce-ice' case
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   palbi      ! all skies ice albedo 
-      REAL(wp), INTENT(in   ), DIMENSION(:,:  ), OPTIONAL ::   psst       ! sea surface temperature     [Celsius]
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   pist       ! ice surface temperature     [Kelvin]
-      !
-      INTEGER ::   jl         ! dummy loop index
+      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   palbi   ! all skies ice albedo 
+      REAL(wp), INTENT(in   ), DIMENSION(:,:  ), OPTIONAL ::   psst    ! sea surface temperature  [Celsius]
+      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   pist    ! ice surface temperature  [Kelvin]
+      !
+      INTEGER ::   jl   ! dummy loop index
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcptn, ztmp, zicefr, zmsk
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot
@@ -1378,8 +1374,8 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_flx')
-      !
-      CALL wrk_alloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot )
-      CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice )
+      IF( nn_timing == 1 )   CALL timing_start('sbc_cpl_ice_flx')
+      !
+      CALL wrk_alloc( jpi,jpj,       zcptn, ztmp, zicefr, zmsk, zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot )
+      CALL wrk_alloc( jpi,jpj,jpl,   zqns_ice, zqsr_ice, zdqns_ice )
 
       IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
@@ -1554,6 +1550,6 @@
       CALL wrk_dealloc( jpi,jpj, zevap, zsnw, zqns_oce, zqprec_ice, zqemp_oce ) 
 #else
-
-      ! clem: this formulation is certainly wrong... but better than it was...
+      !
+      ! clem: this formulation is certainly wrong... but better than it was before...
       zqns_tot(:,:) = zqns_tot(:,:)                       &            ! zqns_tot update over free ocean with:
          &          - ztmp(:,:)                           &            ! remove the latent heat flux of solid precip. melting
@@ -1571,7 +1567,6 @@
          qns_ice(:,:,:) = zqns_ice(:,:,:)
       ENDIF
-
+      !
 #endif
-
       !                                                      ! ========================= !
       SELECT CASE( TRIM( sn_rcv_qsr%cldes ) )                !      solar heat fluxes    !   (qsr)
@@ -1682,8 +1677,8 @@
       fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
 
-      CALL wrk_dealloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot )
-      CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_ice_flx')
+      CALL wrk_dealloc( jpi,jpj,       zcptn, ztmp, zicefr, zmsk, zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot )
+      CALL wrk_dealloc( jpi,jpj,jpl,   zqns_ice, zqsr_ice, zdqns_ice )
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('sbc_cpl_ice_flx')
       !
    END SUBROUTINE sbc_cpl_ice_flx
@@ -1708,10 +1703,10 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_snd')
-      !
-      CALL wrk_alloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
-      CALL wrk_alloc( jpi,jpj,jpl, ztmp3, ztmp4 )
-
-      isec = ( kt - nit000 ) * NINT(rdttra(1))        ! date of exchanges
+      IF( nn_timing == 1 )   CALL timing_start('sbc_cpl_snd')
+      !
+      CALL wrk_alloc( jpi,jpj,       zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
+      CALL wrk_alloc( jpi,jpj,jpl,   ztmp3, ztmp4 )
+
+      isec = ( kt - nit000 ) * NINT( rdt )        ! date of exchanges
 
       zfr_l(:,:) = 1.- fr_i(:,:)
@@ -2031,5 +2026,5 @@
       !                                                        ! first T level thickness 
       IF( ssnd(jps_e3t1st )%laction )  THEN
-         CALL cpl_snd( jps_e3t1st, isec, RESHAPE ( fse3t_n(:,:,1)   , (/jpi,jpj,1/) ), info )
+         CALL cpl_snd( jps_e3t1st, isec, RESHAPE ( e3t_n(:,:,1)   , (/jpi,jpj,1/) ), info )
       ENDIF
       !                                                        ! Qsr fraction
@@ -2049,8 +2044,8 @@
       IF( ssnd(jps_taum  )%laction )  CALL cpl_snd( jps_taum  , isec, RESHAPE ( taum, (/jpi,jpj,1/) ), info )
 
-      CALL wrk_dealloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
-      CALL wrk_dealloc( jpi,jpj,jpl, ztmp3, ztmp4 )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_snd')
+      CALL wrk_dealloc( jpi,jpj,       zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
+      CALL wrk_dealloc( jpi,jpj,jpl,   ztmp3, ztmp4 )
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('sbc_cpl_snd')
       !
    END SUBROUTINE sbc_cpl_snd
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcdcy.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcdcy.F90	(revision 6121)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcdcy.F90	(revision 6140)
@@ -90,6 +90,6 @@
 
       ! When are we during the day (from 0 to 1)
-      zlo = ( REAL(nsec_day, wp) - 0.5_wp * rdttra(1) ) / rday
-      zup = zlo + ( REAL(nn_fsbc, wp)     * rdttra(1) ) / rday
+      zlo = ( REAL(nsec_day, wp) - 0.5_wp * rdt ) / rday
+      zup = zlo + ( REAL(nn_fsbc, wp)     * rdt ) / rday
       !                                          
       IF( nday_qsr == -1 ) THEN       ! first time step only  
@@ -189,5 +189,5 @@
          END DO  
          !
-         ztmp = rday / ( rdttra(1) * REAL(nn_fsbc, wp) )
+         ztmp = rday / ( rdt * REAL(nn_fsbc, wp) )
          rscal(:,:) = rscal(:,:) * ztmp
          !
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 6121)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 6140)
@@ -17,4 +17,5 @@
    USE sbcdcy          ! surface boundary condition: diurnal cycle on qsr
    USE phycst          ! physical constants
+   !
    USE fldread         ! read input fields
    USE iom             ! IOM library
@@ -37,5 +38,4 @@
 
    !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -165,7 +165,5 @@
                WRITE(numout,*) 
                WRITE(numout,*) ' day: ', ndastp , TRIM(sf(jf)%clvar), ' * ', zfact
-               CALL prihre( sf(jf)%fnow, jpi, jpj, 1, jpi, 20, 1, jpj, 10, zfact, numout )
             END DO
-            CALL FLUSH(numout)
          ENDIF
          !
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 6121)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 6140)
@@ -12,20 +12,20 @@
 
    !!----------------------------------------------------------------------
-   !!   sbc_fwb      : freshwater budget for global ocean configurations
-   !!                  in free surface and forced mode
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers
-   USE dom_oce         ! ocean space and time domain
-   USE sbc_oce         ! surface ocean boundary condition
-   USE phycst          ! physical constants
-   USE sbcrnf          ! ocean runoffs
-   USE sbcisf          ! ice shelf melting contribution
-   USE sbcssr          ! SS damping terms
-   USE in_out_manager  ! I/O manager
-   USE lib_mpp         ! distribued memory computing library
-   USE wrk_nemo        ! work arrays
-   USE timing          ! Timing
-   USE lbclnk          ! ocean lateral boundary conditions
-   USE lib_fortran
+   !!   sbc_fwb       : freshwater budget for global ocean configurations (free surface & forced mode)
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers
+   USE dom_oce        ! ocean space and time domain
+   USE sbc_oce        ! surface ocean boundary condition
+   USE phycst         ! physical constants
+   USE sbcrnf         ! ocean runoffs
+   USE sbcisf         ! ice shelf melting contribution
+   USE sbcssr         ! Sea-Surface damping terms
+   !
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! distribued memory computing library
+   USE wrk_nemo       ! work arrays
+   USE timing         ! Timing
+   USE lbclnk         ! ocean lateral boundary conditions
+   USE lib_fortran    ! 
 
    IMPLICIT NONE
@@ -40,5 +40,4 @@
 
    !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -129,5 +128,5 @@
          ENDIF   
          !                                         ! Update fwfold if new year start
-         ikty = 365 * 86400 / rdttra(1)    !!bug  use of 365 days leap year or 360d year !!!!!!!
+         ikty = 365 * 86400 / rdt                  !!bug  use of 365 days leap year or 360d year !!!!!!!
          IF( MOD( kt, ikty ) == 0 ) THEN
             a_fwb_b = a_fwb                           ! mean sea level taking into account the ice+snow
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 6121)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 6140)
@@ -9,6 +9,4 @@
    !!----------------------------------------------------------------------
    !!   sbc_ice_cice  : sea-ice model time-stepping and update ocean sbc over ice-covered area
-   !!   
-   !!   
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
@@ -92,6 +90,4 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:), PRIVATE ::   png     ! local array used in sbc_cice_ice
 
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.7 , NEMO-consortium (2015) 
@@ -244,46 +240,48 @@
          snwice_mass_b(:,:) = 0.0_wp         ! no mass exchanges
       ENDIF
-      IF( .NOT. ln_rstart ) THEN
+      IF( .NOT.ln_rstart ) THEN
          IF( nn_ice_embd == 2 ) THEN            ! full embedment (case 2) deplete the initial ssh below sea-ice area
             sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
             sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
-#if defined key_vvl            
-           ! key_vvl necessary? clem: yes for compilation purpose
-            DO jk = 1,jpkm1                     ! adjust initial vertical scale factors
-               fse3t_n(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshn(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
-               fse3t_b(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshb(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
-            ENDDO
-            fse3t_a(:,:,:) = fse3t_b(:,:,:)
-            ! Reconstruction of all vertical scale factors at now and before time
-            ! steps
-            ! =============================================================================
-            ! Horizontal scale factor interpolations
-            ! --------------------------------------
-            CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
-            CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
-            CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
-            CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
-            CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
-            ! Vertical scale factor interpolations
-            ! ------------------------------------
-            CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
-            CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
-            CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
-            CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
-            CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
-            ! t- and w- points depth
-            ! ----------------------
-            fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
-            fsdepw_n(:,:,1) = 0.0_wp
-            fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
-            DO jk = 2, jpk
-               fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk)
-               fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1)
-               fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:)
-            END DO
-#endif
+
+!!gm This should be put elsewhere....   (same remark for limsbc)
+!!gm especially here it is assumed zstar coordinate, but it can be ztilde....
+            IF( .NOT.ln_linssh ) THEN
+               !
+               DO jk = 1,jpkm1                     ! adjust initial vertical scale factors
+                  e3t_n(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshn(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
+                  e3t_b(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshb(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
+               ENDDO
+               e3t_a(:,:,:) = e3t_b(:,:,:)
+               ! Reconstruction of all vertical scale factors at now and before time-steps
+               ! =============================================================================
+               ! Horizontal scale factor interpolations
+               ! --------------------------------------
+               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
+               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
+               CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
+               CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
+               CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
+               ! Vertical scale factor interpolations
+               ! ------------------------------------
+               CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )
+               CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
+               CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
+               CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
+               CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
+               ! t- and w- points depth
+               ! ----------------------
+               gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
+               gdepw_n(:,:,1) = 0.0_wp
+               gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
+               DO jk = 2, jpk
+                  gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk)
+                  gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1)
+                  gde3w_n(:,:,jk) = gdept_n(:,:,jk  ) - sshn   (:,:)
+               END DO
+            ENDIF
          ENDIF
       ENDIF
- 
+      !
       CALL wrk_dealloc( jpi,jpj, ztmp1, ztmp2 )
       !
@@ -448,5 +446,5 @@
 ! Freezing/melting potential
 ! Calculated over NEMO leapfrog timestep (hence 2*dt)
-      nfrzmlt(:,:)=rau0*rcp*fse3t_m(:,:)*(Tocnfrz-sst_m(:,:))/(2.0*dt)
+      nfrzmlt(:,:) = rau0 * rcp * e3t_m(:,:) * ( Tocnfrz-sst_m(:,:) ) / ( 2.0*dt )
 
       ztmp(:,:) = nfrzmlt(:,:)
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_if.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_if.F90	(revision 6121)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_if.F90	(revision 6140)
@@ -34,6 +34,4 @@
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_ice   ! structure of input ice-cover (file informations, fields read)
    
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 6121)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 6140)
@@ -25,5 +25,5 @@
    USE thd_ice         ! LIM-3: thermodynamical variables
    USE dom_ice         ! LIM-3: ice domain
-
+   !
    USE sbc_oce         ! Surface boundary condition: ocean fields
    USE sbc_ice         ! Surface boundary condition: ice   fields
@@ -32,5 +32,5 @@
    USE sbccpl          ! Surface boundary condition: coupled interface
    USE albedo          ! ocean & ice albedo
-
+   !
    USE phycst          ! Define parameters for the routines
    USE eosbn2          ! equation of state
@@ -47,19 +47,18 @@
    USE limupdate2      ! update of global variables
    USE limvar          ! Ice variables switch
-
+   USE limctl          ! 
    USE limmsh          ! LIM mesh
    USE limistate       ! LIM initial state
    USE limthd_sal      ! LIM ice thermodynamics: salinity
-
+   !
    USE c1d             ! 1D vertical configuration
+   USE in_out_manager  ! I/O manager
+   USE iom             ! I/O manager library
+   USE prtctl          ! Print control
+   USE lib_fortran     ! 
    USE lbclnk          ! lateral boundary condition - MPP link
    USE lib_mpp         ! MPP library
    USE wrk_nemo        ! work arrays
    USE timing          ! Timing
-   USE iom             ! I/O manager library
-   USE in_out_manager  ! I/O manager
-   USE prtctl          ! Print control
-   USE lib_fortran     ! 
-   USE limctl
 
 #if defined key_bdy 
@@ -74,5 +73,4 @@
    
    !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -82,6 +80,4 @@
    !!----------------------------------------------------------------------
 CONTAINS
-
-   !!======================================================================
 
    SUBROUTINE sbc_ice_lim( kt, kblk )
@@ -270,9 +266,9 @@
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   via Louvain la Neuve Ice Model (LIM-3) time stepping'
       !
-                                       ! Open the reference and configuration namelist files and namelist output file 
+      !                                ! Open the reference and configuration namelist files and namelist output file 
       CALL ctl_opn( numnam_ice_ref, 'namelist_ice_ref',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp ) 
       CALL ctl_opn( numnam_ice_cfg, 'namelist_ice_cfg',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
       IF(lwm) CALL ctl_opn( numoni, 'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 )
-
+      !
       CALL ice_run                     ! set some ice run parameters
       !
@@ -348,11 +344,10 @@
       REWIND( numnam_ice_ref )              ! Namelist namicerun in reference namelist : Parameters for ice
       READ  ( numnam_ice_ref, namicerun, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in reference namelist', lwp )
-
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namicerun in reference namelist', lwp )
+      !
       REWIND( numnam_ice_cfg )              ! Namelist namicerun in configuration namelist : Parameters for ice
       READ  ( numnam_ice_cfg, namicerun, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in configuration namelist', lwp )
-      IF(lwm) WRITE ( numoni, namicerun )
-      !
+902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namicerun in configuration namelist', lwp )
+      IF(lwm) WRITE( numoni, namicerun )
       !
       IF(lwp) THEN                        ! control print
@@ -373,5 +368,5 @@
       !
       ! sea-ice timestep and inverse
-      rdt_ice   = nn_fsbc * rdttra(1)  
+      rdt_ice   = nn_fsbc * rdt  
       r1_rdtice = 1._wp / rdt_ice 
 
@@ -405,11 +400,10 @@
       REWIND( numnam_ice_ref )              ! Namelist namiceitd in reference namelist : Parameters for ice
       READ  ( numnam_ice_ref, namiceitd, IOSTAT = ios, ERR = 903)
-903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in reference namelist', lwp )
-
+903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namiceitd in reference namelist', lwp )
+      !
       REWIND( numnam_ice_cfg )              ! Namelist namiceitd in configuration namelist : Parameters for ice
       READ  ( numnam_ice_cfg, namiceitd, IOSTAT = ios, ERR = 904 )
-904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in configuration namelist', lwp )
-      IF(lwm) WRITE ( numoni, namiceitd )
-      !
+904   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namiceitd in configuration namelist', lwp )
+      IF(lwm) WRITE( numoni, namiceitd )
       !
       IF(lwp) THEN                        ! control print
@@ -417,8 +411,8 @@
          WRITE(numout,*) 'ice_itd : ice cat distribution'
          WRITE(numout,*) ' ~~~~~~'
-         WRITE(numout,*) '   shape of ice categories distribution                          nn_catbnd = ', nn_catbnd
-         WRITE(numout,*) '   mean ice thickness in the domain (only active if nn_catbnd=2) rn_himean = ', rn_himean
+         WRITE(numout,*) '   shape of ice categories distribution                     nn_catbnd = ', nn_catbnd
+         WRITE(numout,*) '   mean ice thickness in the domain (used if nn_catbnd=2)   rn_himean = ', rn_himean
       ENDIF
-
+      !
       !----------------------------------
       !- Thickness categories boundaries 
@@ -427,42 +421,35 @@
       IF(lwp) WRITE(numout,*) 'lim_itd_init : Initialization of ice cat distribution '
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
-
+      !
       hi_max(:) = 0._wp
-
-      SELECT CASE ( nn_catbnd  )       
-                                   !----------------------
-         CASE (1)                  ! tanh function (CICE)
-                                   !----------------------
+      !
+      SELECT CASE ( nn_catbnd  )    ! type of ice categories distribution
+      !
+      CASE (1)                            !==  tanh function (CICE)  ==!
          zc1 =  3._wp / REAL( jpl, wp )
          zc2 = 10._wp * zc1
          zc3 =  3._wp
-
          DO jl = 1, jpl
             zx1 = REAL( jl-1, wp ) / REAL( jpl, wp )
             hi_max(jl) = hi_max(jl-1) + zc1 + zc2 * (1._wp + TANH( zc3 * (zx1 - 1._wp ) ) )
          END DO
-
-                                   !----------------------
-         CASE (2)                  ! h^(-alpha) function
-                                   !----------------------
-         zalpha = 0.05             ! exponent of the transform function
-
-         zhmax  = 3.*rn_himean
-
+         !
+      CASE (2)                            !==  h^(-alpha) function  ==!
+         zalpha = 0.05_wp
+         zhmax  = 3._wp * rn_himean
          DO jl = 1, jpl 
             znum = jpl * ( zhmax+1 )**zalpha
-            zden = ( jpl - jl ) * ( zhmax+1 )**zalpha + jl
+            zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp )
             hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
          END DO
-
+         !
       END SELECT
-
-      DO jl = 1, jpl
+      !
+      DO jl = 1, jpl                ! mean thickness by category
          hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
       END DO
-
-      ! Set hi_max(jpl) to a big value to ensure that all ice is thinner than hi_max(jpl)
-      hi_max(jpl) = 99._wp
-
+      !
+      hi_max(jpl) = 99._wp          ! set to a big value to ensure that all ice is thinner than hi_max(jpl)
+      !
       IF(lwp) WRITE(numout,*) ' Thickness category boundaries '
       IF(lwp) WRITE(numout,*) ' hi_max ', hi_max(0:jpl)
@@ -471,5 +458,6 @@
 
    
-   SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_limflx )
+   SUBROUTINE ice_lim_flx( ptn_ice , palb_ice, pqns_ice ,    &
+      &                    pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_limflx )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE ice_lim_flx  ***
@@ -483,5 +471,5 @@
       !!---------------------------------------------------------------------
       INTEGER                   , INTENT(in   ) ::   k_limflx   ! =-1 do nothing; =0 average ; 
-                                                                ! =1 average and redistribute ; =2 redistribute
+      !                                                         ! =1 average and redistribute ; =2 redistribute
       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptn_ice    ! ice surface temperature 
       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   palb_ice   ! ice albedo
@@ -503,7 +491,6 @@
       REAL(wp), POINTER, DIMENSION(:,:) :: z_devap_m ! Mean d(evap)/dT over all categories
       !!----------------------------------------------------------------------
-
+      !
       IF( nn_timing == 1 )  CALL timing_start('ice_lim_flx')
-      !
       !
       SELECT CASE( k_limflx )                              !==  averaged on all ice categories  ==!
@@ -529,5 +516,5 @@
          CALL wrk_dealloc( jpi,jpj, z_qsr_m, z_qns_m, z_evap_m, z_dqn_m, z_devap_m)
       END SELECT
-
+      !
       SELECT CASE( k_limflx )                              !==  redistribution on all ice categories  ==!
       CASE( 1 , 2 )
@@ -548,4 +535,5 @@
       !
    END SUBROUTINE ice_lim_flx
+
 
    SUBROUTINE sbc_lim_bef
@@ -564,6 +552,7 @@
       u_ice_b(:,:)     = u_ice(:,:)
       v_ice_b(:,:)     = v_ice(:,:)
-      
+      !      
    END SUBROUTINE sbc_lim_bef
+
 
    SUBROUTINE sbc_lim_diag0
@@ -580,5 +569,5 @@
       sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp
       sfx_res(:,:) = 0._wp
-      
+      !
       wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp
       wfx_sni(:,:) = 0._wp   ;   wfx_opw(:,:) = 0._wp
@@ -587,5 +576,5 @@
       wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp
       wfx_spr(:,:) = 0._wp   ;   
-      
+      !
       hfx_thd(:,:) = 0._wp   ;   
       hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp
@@ -596,11 +585,11 @@
       hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp
       hfx_err_dif(:,:) = 0._wp   ;
-
+      !
       afx_tot(:,:) = 0._wp   ;
       afx_dyn(:,:) = 0._wp   ;   afx_thd(:,:) = 0._wp
-
+      !
       diag_heat(:,:) = 0._wp ;   diag_smvi(:,:) = 0._wp ;
       diag_vice(:,:) = 0._wp ;   diag_vsnw(:,:) = 0._wp ;
-      
+      !
    END SUBROUTINE sbc_lim_diag0
 
@@ -634,5 +623,4 @@
    END FUNCTION fice_ice_ave
 
-
 #else
    !!----------------------------------------------------------------------
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90	(revision 6121)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90	(revision 6140)
@@ -64,5 +64,4 @@
 
    !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 6121)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 6140)
@@ -35,39 +35,37 @@
    ! public in order to be able to output then 
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   risf_tsc_b, risf_tsc   
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf              !: net heat flux from ice shelf
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   risf_tsc_b, risf_tsc  !: before and now T & S isf contents [K.m/s & PSU.m/s]  
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf                  !: net heat flux from ice shelf      [W/m2]
    REAL(wp), PUBLIC ::   rn_hisf_tbl                 !: thickness of top boundary layer [m]
-   LOGICAL , PUBLIC ::   ln_divisf                   !: flag to correct divergence 
-   INTEGER , PUBLIC ::   nn_isfblk                   !: 
-   INTEGER , PUBLIC ::   nn_gammablk                 !:
-   LOGICAL , PUBLIC ::   ln_conserve                 !:
-   REAL(wp), PUBLIC ::   rn_gammat0                  !: temperature exchange coeficient
-   REAL(wp), PUBLIC ::   rn_gammas0                  !: salinity    exchange coeficient 
-   REAL(wp), PUBLIC ::   rdivisf                     !: flag to test if fwf apply on divergence
+   INTEGER , PUBLIC ::   nn_isf                      !: flag to choose between explicit/param/specified  
+   INTEGER , PUBLIC ::   nn_isfblk                   !: flag to choose the bulk formulation to compute the ice shelf melting
+   INTEGER , PUBLIC ::   nn_gammablk                 !: flag to choose how the exchange coefficient is computed
+   REAL(wp), PUBLIC ::   rn_gammat0                  !: temperature exchange coeficient []
+   REAL(wp), PUBLIC ::   rn_gammas0                  !: salinity    exchange coeficient []
 
    REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  rzisf_tbl              !:depth of calving front (shallowest point) nn_isf ==2/3
-   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  rhisf_tbl, rhisf_tbl_0 !:thickness of tbl
+   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  rhisf_tbl, rhisf_tbl_0 !:thickness of tbl  [m]
    REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  r1_hisf_tbl            !:1/thickness of tbl
    REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  ralpha                 !:proportion of bottom cell influenced by tbl 
    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
-   INTEGER,    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  misfkt, misfkb         !:Level of ice shelf base
-
-   REAL(wp), PUBLIC, SAVE ::   rcpi   = 2000.0_wp     ! phycst ?
-   REAL(wp), PUBLIC, SAVE ::   kappa  = 1.54e-6_wp    ! phycst ?
-   REAL(wp), PUBLIC, SAVE ::   rhoisf = 920.0_wp      ! phycst ?
-   REAL(wp), PUBLIC, SAVE ::   tsurf  = -20.0_wp      ! phycst ?
-   REAL(wp), PUBLIC, SAVE ::   lfusisf= 0.334e6_wp    ! phycst ?
+   INTEGER,    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)      ::  misfkt, misfkb         !:Level of ice shelf base
+
+   REAL(wp), PUBLIC, SAVE ::   rcpi     = 2000.0_wp     ! specific heat of ice shelf             [J/kg/K]
+   REAL(wp), PUBLIC, SAVE ::   rkappa   = 1.54e-6_wp    ! heat diffusivity through the ice-shelf [m2/s]
+   REAL(wp), PUBLIC, SAVE ::   rhoisf   = 920.0_wp      ! volumic mass of ice shelf              [kg/m3]
+   REAL(wp), PUBLIC, SAVE ::   tsurf    = -20.0_wp      ! air temperature on top of ice shelf    [C]
+   REAL(wp), PUBLIC, SAVE ::   rlfusisf = 0.334e6_wp    ! latent heat of fusion of ice shelf     [J/kg]
 
 !: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3)
-   CHARACTER(len=100), PUBLIC ::   cn_dirisf  = './'    !: Root directory for location of ssr files
-   TYPE(FLD_N)       , PUBLIC ::   sn_qisf, sn_fwfisf     !: information about the runoff file to be read
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_qisf, sf_fwfisf
-   TYPE(FLD_N)       , PUBLIC ::   sn_rnfisf              !: information about the runoff file to be read
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_rnfisf           
-   TYPE(FLD_N)       , PUBLIC ::   sn_depmax_isf, sn_depmin_isf, sn_Leff_isf     !: information about the runoff file to be read
+   CHARACTER(len=100), PUBLIC           :: cn_dirisf  = './' !: Root directory for location of ssr files
+   TYPE(FLD_N)       , PUBLIC           :: sn_fwfisf         !: information about the isf melting file to be read
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_fwfisf
+   TYPE(FLD_N)       , PUBLIC           :: sn_rnfisf         !: information about the isf melting param.   file to be read
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnfisf           
+   TYPE(FLD_N)       , PUBLIC           :: sn_depmax_isf     !: information about the grounding line depth file to be read
+   TYPE(FLD_N)       , PUBLIC           :: sn_depmin_isf     !: information about the calving   line depth file to be read
+   TYPE(FLD_N)       , PUBLIC           :: sn_Leff_isf       !: information about the effective length     file to be read
    
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.7 , LOCEAN-IPSL (2015)
@@ -77,22 +75,21 @@
 CONTAINS
  
-   SUBROUTINE sbc_isf(kt)
+  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
+      !!                  ***  ROUTINE sbc_isf  ***
+      !!
+      !! ** Purpose : Compute Salt and Heat fluxes related to ice_shelf 
+      !!              melting and freezing 
+      !!
+      !! ** Method  :  4 parameterizations are available according to nn_isf 
+      !!               nn_isf = 1 : Realistic ice_shelf formulation
+      !!                        2 : Beckmann & Goose parameterization
+      !!                        3 : Specified runoff in deptht (Mathiot & al. )
+      !!                        4 : specified fwf and heat flux forcing beneath the ice shelf
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) :: kt                   ! ocean time step
+      !
+      INTEGER               :: ji, jj               ! loop index
+      REAL(wp), DIMENSION (:,:), POINTER :: zt_frz, zdep ! freezing temperature (zt_frz) at depth (zdep) 
       !!---------------------------------------------------------------------
       !
@@ -100,189 +97,72 @@
       IF( kt == nit000 ) THEN                   !  First call kt=nit000  !
          !                                      ! ====================== !
-         REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs 
-         READ  ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901)
-901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp )
-
-         REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs
-         READ  ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 )
-902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp )
-         IF(lwm) WRITE ( numond, namsbc_isf )
-
-         IF ( lwp ) WRITE(numout,*)
-         IF ( lwp ) WRITE(numout,*) 'sbc_isf: heat flux of the ice shelf'
-         IF ( lwp ) WRITE(numout,*) '~~~~~~~~~'
-         IF ( lwp ) WRITE(numout,*) 'sbcisf :' 
-         IF ( lwp ) WRITE(numout,*) '~~~~~~~~'
-         IF ( lwp ) WRITE(numout,*) '        nn_isf      = ', nn_isf
-         IF ( lwp ) WRITE(numout,*) '        nn_isfblk   = ', nn_isfblk
-         IF ( lwp ) WRITE(numout,*) '        rn_hisf_tbl = ', rn_hisf_tbl
-         IF ( lwp ) WRITE(numout,*) '        ln_divisf   = ', ln_divisf 
-         IF ( lwp ) WRITE(numout,*) '        nn_gammablk = ', nn_gammablk 
-         IF ( lwp ) WRITE(numout,*) '        rn_tfri2    = ', rn_tfri2 
-         IF (ln_divisf) THEN       ! keep it in the namelist ??? used true anyway as for runoff ? (PM)
-            rdivisf = 1._wp
-         ELSE
-            rdivisf = 0._wp
-         END IF
-         !
-         ! Allocate public variable
-         IF ( sbc_isf_alloc()  /= 0 )         CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' )
-         !
-         ! initialisation
-         qisf(:,:)        = 0._wp  ; fwfisf(:,:) = 0._wp
-         risf_tsc(:,:,:)  = 0._wp
-         !
-         ! define isf tbl tickness, top and bottom indice
-         IF      (nn_isf == 1) THEN
-            rhisf_tbl(:,:) = rn_hisf_tbl
-            misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
-         ELSE IF ((nn_isf == 3) .OR. (nn_isf == 2)) THEN
-            ALLOCATE( sf_rnfisf(1), STAT=ierror )
-            ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) )
-            CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
-
-            !: read effective lenght (BG03)
-            IF (nn_isf == 2) THEN
-               ! Read Data and save some integral values
-               CALL iom_open( sn_Leff_isf%clname, inum )
-               cvarLeff  = 'soLeff'               !: variable name for Efficient Length scale
-               CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1)
-               CALL iom_close(inum)
-               !
-               risfLeff = risfLeff*1000           !: convertion in m
-            END IF
-
-           ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth)
-            CALL iom_open( sn_depmax_isf%clname, inum )
-            cvarhisf = TRIM(sn_depmax_isf%clvar)
-            CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base
-            CALL iom_close(inum)
-            !
-            CALL iom_open( sn_depmin_isf%clname, inum )
-            cvarzisf = TRIM(sn_depmin_isf%clvar)
-            CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base
-            CALL iom_close(inum)
-            !
-            rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:)        !: tickness isf boundary layer
-
-           !! compute first level of the top boundary layer
-           DO ji = 1, jpi
-              DO jj = 1, jpj
-                  jk = 2
-                  DO WHILE ( jk .LE. mbkt(ji,jj) .AND. fsdepw(ji,jj,jk) < rzisf_tbl(ji,jj) ) ;  jk = jk + 1 ;  END DO
-                  misfkt(ji,jj) = jk-1
-               END DO
-            END DO
-
-         ELSE IF ( nn_isf == 4 ) THEN
-            ! as in nn_isf == 1
-            rhisf_tbl(:,:) = rn_hisf_tbl
-            misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
-            
-            ! load variable used in fldread (use for temporal interpolation of isf fwf forcing)
-            ALLOCATE( sf_fwfisf(1), sf_qisf(1), STAT=ierror )
-            ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) )
-            ALLOCATE( sf_qisf(1)%fnow(jpi,jpj,1), sf_qisf(1)%fdta(jpi,jpj,1,2) )
-            CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
-            !CALL fld_fill( sf_qisf  , (/ sn_qisf   /), cn_dirisf, 'sbc_isf_init', 'read heat flux isf data'       , 'namsbc_isf' )
-         END IF
-         
-         ! compute bottom level of isf tbl and thickness of tbl below the ice shelf
-         rhisf_tbl_0(:,:) = rhisf_tbl(:,:)
-         DO jj = 1,jpj
-            DO ji = 1,jpi
-               ikt = misfkt(ji,jj)
-               ikb = misfkt(ji,jj)
-               ! thickness of boundary layer at least the top level thickness
-               rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t_n(ji,jj,ikt))
-
-               ! determine the deepest level influenced by the boundary layer
-               ! test on tmask useless ?????
-               DO jk = ikt, mbkt(ji,jj)
-                  IF ( (SUM(fse3t_n(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_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
-               misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl
-               r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
-
-               zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1
-               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
-            END DO
-         END DO
+         CALL sbc_isf_init
+      !                                         ! ---------------------------------------- !
+      ELSE                                      !          Swap of forcing fields          !
+         !                                      ! ---------------------------------------- !
+         fwfisf_b  (:,:  ) = fwfisf  (:,:  )    ! Swap the ocean forcing fields except at nit000
+         risf_tsc_b(:,:,:) = risf_tsc(:,:,:)    ! where before fields are set at the end of the routine
          !
       END IF
 
-      !                                            ! ---------------------------------------- !
-      IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          !
-         !                                         ! ---------------------------------------- !
-         fwfisf_b  (:,:  ) = fwfisf  (:,:  )               ! Swap the ocean forcing fields except at nit000
-         risf_tsc_b(:,:,:) = risf_tsc(:,:,:)               ! where before fields are set at the end of the routine
-         !
-      ENDIF
-
       IF( MOD( kt-1, nn_fsbc) == 0 ) THEN
-
-
-         ! compute salf and heat flux
-         IF (nn_isf == 1) THEN
-            ! realistic ice shelf formulation
+         ! allocation
+         CALL wrk_alloc( jpi,jpj, zt_frz, zdep  )
+
+         ! compute salt and heat flux
+         SELECT CASE ( nn_isf )
+         CASE ( 1 )    ! realistic ice shelf formulation
             ! compute T/S/U/V for the top boundary layer
             CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T')
             CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T')
-            CALL sbc_isf_tbl(un(:,:,:),utbl(:,:),'U')
-            CALL sbc_isf_tbl(vn(:,:,:),vtbl(:,:),'V')
+            CALL sbc_isf_tbl(un(:,:,:)        ,utbl(:,:),'U')
+            CALL sbc_isf_tbl(vn(:,:,:)        ,vtbl(:,:),'V')
             ! iom print
             CALL iom_put('ttbl',ttbl(:,:))
             CALL iom_put('stbl',stbl(:,:))
-            CALL iom_put('utbl',utbl(:,:))
-            CALL iom_put('vtbl',vtbl(:,:))
+            CALL iom_put('utbl',utbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:))
+            CALL iom_put('vtbl',vtbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:))
             ! compute fwf and heat flux
             CALL sbc_isf_cav (kt)
 
-         ELSE IF (nn_isf == 2) THEN
-            ! Beckmann and Goosse parametrisation 
+         CASE ( 2 )    ! Beckmann and Goosse parametrisation 
             stbl(:,:)   = soce
             CALL sbc_isf_bg03(kt)
 
-         ELSE IF (nn_isf == 3) THEN
-            ! specified runoff in depth (Mathiot et al., XXXX in preparation)
+         CASE ( 3 )    ! specified runoff in depth (Mathiot et al., XXXX in preparation)
             CALL fld_read ( kt, nn_fsbc, sf_rnfisf   )
-            fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting) 
-            qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux
+            fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fwf  flux from the isf (fwfisf <0 mean melting) 
+            qisf(:,:)   = fwfisf(:,:) * rlfusisf             ! heat flux
             stbl(:,:)   = soce
 
-         ELSE IF (nn_isf == 4) THEN
-            ! specified fwf and heat flux forcing beneath the ice shelf
+         CASE ( 4 )    ! specified fwf and heat flux forcing beneath the ice shelf
             CALL fld_read ( kt, nn_fsbc, sf_fwfisf   )
-            !CALL fld_read ( kt, nn_fsbc, sf_qisf   )
-            fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1)            ! fwf
-            qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux
-            !qisf(:,:)   = sf_qisf(1)%fnow(:,:,1)              ! heat flux
+            fwfisf(:,:) = - sf_fwfisf(1)%fnow(:,:,1)           ! fwf  flux from the isf (fwfisf <0 mean melting)
+            qisf(:,:)   = fwfisf(:,:) * rlfusisf               ! heat flux
             stbl(:,:)   = soce
 
-         END IF
+         END SELECT
+
          ! compute tsc due to isf
-         ! 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 !
+         ! isf melting implemented as a volume flux and we assume that melt water is at 0 PSU.
+         ! WARNING water add at temp = 0C, need to add a correction term (fwfisf * tfreez / rau0).
+         ! compute freezing point beneath ice shelf (or top cell if nn_isf = 3)
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               zdep(ji,jj)=gdepw_n(ji,jj,misfkt(ji,jj))
+            END DO
+         END DO
+         CALL eos_fzp( stbl(:,:), zt_frz(:,:), zdep(:,:) )
          
-         ! 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(:,:)         
- 
+         risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - fwfisf(:,:) * zt_frz(:,:) * r1_rau0 !
+         risf_tsc(:,:,jp_sal) = 0.0_wp
+
          ! lbclnk
          CALL lbc_lnk(risf_tsc(:,:,jp_tem),'T',1.)
          CALL lbc_lnk(risf_tsc(:,:,jp_sal),'T',1.)
-         CALL lbc_lnk(fwfisf(:,:)   ,'T',1.)
-         CALL lbc_lnk(qisf(:,:)     ,'T',1.)
-
-         IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
+         CALL lbc_lnk(fwfisf(:,:)         ,'T',1.)
+         CALL lbc_lnk(qisf(:,:)           ,'T',1.)
+
+         IF( kt == nit000 ) THEN                         !   set the forcing field at nit000 - 1    !
             IF( ln_rstart .AND.    &                     ! Restart: read in restart file
                  & iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN
@@ -295,6 +175,12 @@
                risf_tsc_b(:,:,:)= risf_tsc(:,:,:)
             END IF
-         ENDIF
+         END IF
          ! 
+         ! output
+         CALL iom_put('qisf'  , qisf)
+         CALL iom_put('fwfisf', fwfisf)
+
+         ! deallocation
+         CALL wrk_dealloc( jpi,jpj, zt_frz, zdep  )
       END IF
       !  
@@ -315,84 +201,205 @@
                &    STAT= sbc_isf_alloc )
          !
-         IF( lk_mpp                  )   CALL mpp_sum ( sbc_isf_alloc )
+         IF( lk_mpp             )   CALL mpp_sum ( sbc_isf_alloc )
          IF( sbc_isf_alloc /= 0 )   CALL ctl_warn('sbc_isf_alloc: failed to allocate arrays.')
          !
-      ENDIF
+      END IF
   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
-      !!
+  SUBROUTINE sbc_isf_init
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_isf_init  ***
+      !!
+      !! ** Purpose : Initialisation of variables for iceshelf fluxes formulation
+      !!
+      !! ** Method  :  4 parameterizations are available according to nn_isf 
+      !!               nn_isf = 1 : Realistic ice_shelf formulation
+      !!                        2 : Beckmann & Goose parameterization
+      !!                        3 : Specified runoff in deptht (Mathiot & al. )
+      !!                        4 : specified fwf and heat flux forcing beneath the ice shelf
+      !!----------------------------------------------------------------------
+      INTEGER               :: ji, jj, jk           ! loop index
+      INTEGER               :: ik                   ! current level index
+      INTEGER               :: ikt, ikb             ! top and bottom level of the isf boundary layer
+      INTEGER               :: inum, ierror
+      INTEGER               :: ios                  ! Local integer output status for namelist read
+      REAL(wp)              :: zhk
+      CHARACTER(len=256)    :: cvarzisf, cvarhisf   ! name for isf file
+      CHARACTER(LEN=32 )    :: cvarLeff             ! variable name for efficient Length scale
+      !!----------------------------------------------------------------------
+      NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, rn_gammat0, rn_gammas0, nn_gammablk, nn_isf, &
+                         & sn_fwfisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf
+      !!----------------------------------------------------------------------
+
+      REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs 
+      READ  ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs
+      READ  ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namsbc_isf )
+
+      IF ( lwp ) WRITE(numout,*)
+      IF ( lwp ) WRITE(numout,*) 'sbc_isf: heat flux of the ice shelf'
+      IF ( lwp ) WRITE(numout,*) '~~~~~~~~~'
+      IF ( lwp ) WRITE(numout,*) 'sbcisf :' 
+      IF ( lwp ) WRITE(numout,*) '~~~~~~~~'
+      IF ( lwp ) WRITE(numout,*) '        nn_isf      = ', nn_isf
+      IF ( lwp ) WRITE(numout,*) '        nn_isfblk   = ', nn_isfblk
+      IF ( lwp ) WRITE(numout,*) '        rn_hisf_tbl = ', rn_hisf_tbl
+      IF ( lwp ) WRITE(numout,*) '        nn_gammablk = ', nn_gammablk 
+      IF ( lwp ) WRITE(numout,*) '        rn_gammat0  = ', rn_gammat0  
+      IF ( lwp ) WRITE(numout,*) '        rn_gammas0  = ', rn_gammas0  
+      IF ( lwp ) WRITE(numout,*) '        rn_tfri2    = ', rn_tfri2 
+      !
+      ! Allocate public variable
+      IF ( sbc_isf_alloc()  /= 0 )         CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' )
+      !
+      ! initialisation
+      qisf(:,:)        = 0._wp  ; fwfisf  (:,:) = 0._wp
+      risf_tsc(:,:,:)  = 0._wp  ; fwfisf_b(:,:) = 0._wp
+      !
+      ! define isf tbl tickness, top and bottom indice
+      SELECT CASE ( nn_isf )
+      CASE ( 1 ) 
+         rhisf_tbl(:,:) = rn_hisf_tbl
+         misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
+
+      CASE ( 2 , 3 )
+         ALLOCATE( sf_rnfisf(1), STAT=ierror )
+         ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) )
+         CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
+
+         !  read effective lenght (BG03)
+         IF (nn_isf == 2) THEN
+            CALL iom_open( sn_Leff_isf%clname, inum )
+            cvarLeff = TRIM(sn_Leff_isf%clvar)
+            CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1)
+            CALL iom_close(inum)
+            !
+            risfLeff = risfLeff*1000.0_wp           !: convertion in m
+         END IF
+
+         ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth)
+         CALL iom_open( sn_depmax_isf%clname, inum )
+         cvarhisf = TRIM(sn_depmax_isf%clvar)
+         CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base
+         CALL iom_close(inum)
+         !
+         CALL iom_open( sn_depmin_isf%clname, inum )
+         cvarzisf = TRIM(sn_depmin_isf%clvar)
+         CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base
+         CALL iom_close(inum)
+         !
+         rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:)        !: tickness isf boundary layer
+
+         !! compute first level of the top boundary layer
+         DO ji = 1, jpi
+            DO jj = 1, jpj
+                ik = 2
+                DO WHILE ( ik <= mbkt(ji,jj) .AND. gdepw_n(ji,jj,ik) < rzisf_tbl(ji,jj) ) ;  ik = ik + 1 ;  END DO
+                misfkt(ji,jj) = ik-1
+            END DO
+         END DO
+
+      CASE ( 4 ) 
+         ! as in nn_isf == 1
+         rhisf_tbl(:,:) = rn_hisf_tbl
+         misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
+         
+         ! load variable used in fldread (use for temporal interpolation of isf fwf forcing)
+         ALLOCATE( sf_fwfisf(1), STAT=ierror )
+         ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) )
+         CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
+
+      END SELECT
+         
+      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
+            ikt = misfkt(ji,jj)
+            ikb = misfkt(ji,jj)
+            ! thickness of boundary layer at least the top level thickness
+            rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt))
+
+            ! determine the deepest level influenced by the boundary layer
+            DO jk = ikt+1, mbkt(ji,jj)
+               IF ( (SUM(e3t_n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
+            END DO
+            rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.
+            misfkb(ji,jj) = ikb                                                   ! last wet level of the tbl
+            r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
+
+            zhk           = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1
+            ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
+         END DO
+      END DO
+
+  END SUBROUTINE sbc_isf_init
+
+  SUBROUTINE sbc_isf_bg03(kt)
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_isf_bg03  ***
+      !!
+      !! ** Purpose : add net heat and fresh water flux from ice shelf melting
+      !!          into the adjacent ocean
+      !!
+      !! ** Method  :   See reference
+      !!
+      !! ** Reference : Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean
+      !!         interaction for climate models", Ocean Modelling 5(2003) 157-170.
+      !!         (hereafter BG)
       !! History :
-      !!      !  06-02  (C. Wang) Original code
+      !!         06-02  (C. Wang) Original code
       !!----------------------------------------------------------------------
       INTEGER, INTENT ( in ) :: kt
       !
-    INTEGER :: ji, jj, jk, jish  !temporary integer
-    INTEGER :: ijkmin
-    INTEGER :: ii, ij, ik 
-    INTEGER :: inum
-
-    REAL(wp) :: zt_sum      ! sum of the temperature between 200m and 600m
-    REAL(wp) :: zt_ave      ! averaged temperature between 200m and 600m
-    REAL(wp) :: zt_frz      ! freezing point temperature at depth z
-    REAL(wp) :: zpress      ! pressure to compute the freezing point in depth
-    
-    !!----------------------------------------------------------------------
-    IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03')
-     !
-
-    ! This test is false only in the very first time step of a run (JMM ???- Initialy build to skip 1rst year of run )
-    DO ji = 1, jpi
-       DO jj = 1, jpj
-          ik = misfkt(ji,jj)
-          !! Initialize arrays to 0 (each step)
-          zt_sum = 0.e0_wp
-          IF ( ik .GT. 1 ) THEN
-    ! 3. -----------the average temperature between 200m and 600m ---------------------
-             DO jk = misfkt(ji,jj),misfkb(ji,jj)
-             ! freezing point temperature  at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!)
-             ! after verif with UNESCO, wrong sign in BG eq. 2
-             ! Calculate freezing temperature
-                zpress = grav*rau0*fsdept(ji,jj,ik)*1.e-04 
-                CALL eos_fzp(tsb(ji,jj,ik,jp_sal), zt_frz, zpress) 
-                zt_sum = zt_sum + (tsn(ji,jj,ik,jp_tem)-zt_frz) * fse3t(ji,jj,ik) * tmask(ji,jj,ik)  ! sum temp
-             ENDDO
-             zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value
-    
-    ! 4. ------------Net heat flux and fresh water flux due to the ice shelf
-          ! For those corresponding to zonal boundary    
-             qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave  &
-                         & / (e1t(ji,jj) * e2t(ji,jj)) * tmask(ji,jj,ik) 
+      INTEGER  :: ji, jj, jk ! dummy loop index
+      INTEGER  :: ik         ! current level
+      REAL(wp) :: zt_sum     ! sum of the temperature between 200m and 600m
+      REAL(wp) :: zt_ave     ! averaged temperature between 200m and 600m
+      REAL(wp) :: zt_frz     ! freezing point temperature at depth z
+      REAL(wp) :: zpress     ! pressure to compute the freezing point in depth
+      !!----------------------------------------------------------------------
+
+      IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03')
+      !
+      DO ji = 1, jpi
+         DO jj = 1, jpj
+            ik = misfkt(ji,jj)
+            !! Initialize arrays to 0 (each step)
+            zt_sum = 0.e0_wp
+            IF ( ik > 1 ) THEN
+               ! 1. -----------the average temperature between 200m and 600m ---------------------
+               DO jk = misfkt(ji,jj),misfkb(ji,jj)
+                  ! freezing point temperature  at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!)
+                  ! after verif with UNESCO, wrong sign in BG eq. 2
+                  ! Calculate freezing temperature
+                  CALL eos_fzp(stbl(ji,jj), zt_frz, zpress) 
+                  zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * e3t_n(ji,jj,jk) * tmask(ji,jj,jk)  ! sum temp
+               END DO
+               zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value
+               ! 2. ------------Net heat flux and fresh water flux due to the ice shelf
+               ! For those corresponding to zonal boundary    
+               qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave  &
+                           & * r1_e1e2t(ji,jj) * tmask(ji,jj,jk)
              
-             fwfisf(ji,jj) = qisf(ji,jj) / lfusisf          !fresh water flux kg/(m2s)                  
-             fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) )
-             !add to salinity trend
-          ELSE
-             qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp
-          END IF
-       END DO
-    END DO
-    !
-    IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_bg03')
+               fwfisf(ji,jj) = qisf(ji,jj) / rlfusisf          !fresh water flux kg/(m2s)                  
+               fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) )
+               !add to salinity trend
+            ELSE
+               qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp
+            END IF
+         END DO
+      END DO
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_bg03')
       !
   END SUBROUTINE sbc_isf_bg03
 
-
-   SUBROUTINE sbc_isf_cav( kt )
+  SUBROUTINE sbc_isf_cav( kt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE sbc_isf_cav  ***
@@ -409,155 +416,119 @@
       INTEGER, INTENT(in)          ::   kt         ! ocean time step
       !
-      LOGICAL :: ln_isomip = .true.
-      REAL(wp), DIMENSION(:,:), POINTER       ::   zfrz,zpress,zti
-      REAL(wp), DIMENSION(:,:), POINTER       ::   zgammat2d, zgammas2d 
-      !REAL(wp), DIMENSION(:,:), POINTER ::   zqisf, zfwfisf
+      INTEGER  ::   ji, jj     ! dummy loop indices
+      INTEGER  ::   nit
       REAL(wp) ::   zlamb1, zlamb2, zlamb3
       REAL(wp) ::   zeps1,zeps2,zeps3,zeps4,zeps6,zeps7
       REAL(wp) ::   zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac
-      REAL(wp) ::   zfwflx, zhtflx, zhtflx_b
-      REAL(wp) ::   zgammat, zgammas
-      REAL(wp) ::   zeps   =  -1.e-20_wp        !==   Local constant initialization   ==!
-      INTEGER  ::   ji, jj     ! dummy loop indices
-      INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integers
-      INTEGER  ::   ierror     ! return error code
-      LOGICAL  ::   lit=.TRUE.
-      INTEGER  ::   nit
+      REAL(wp) ::   zeps = 1.e-20_wp        
+      REAL(wp) ::   zerr
+      REAL(wp), DIMENSION(:,:), POINTER ::   zfrz
+      REAL(wp), DIMENSION(:,:), POINTER ::   zgammat, zgammas 
+      REAL(wp), DIMENSION(:,:), POINTER ::   zfwflx, zhtflx, zhtflx_b
+      LOGICAL  ::   lit
       !!---------------------------------------------------------------------
-      !
-      ! coeficient for linearisation of tfreez
-      zlamb1=-0.0575
-      zlamb2=0.0901
-      zlamb3=-7.61e-04
+      ! coeficient for linearisation of potential tfreez
+      ! Crude approximation for pressure (but commonly used)
+      zlamb1 =-0.0573_wp
+      zlamb2 = 0.0832_wp
+      zlamb3 =-7.53e-08_wp * grav * rau0
       IF( nn_timing == 1 )  CALL timing_start('sbc_isf_cav')
       !
-      CALL wrk_alloc( jpi,jpj, zfrz,zpress,zti, zgammat2d, zgammas2d )
-
-      zcfac=0.0_wp 
-      IF (ln_conserve)  zcfac=1.0_wp
-      zpress(:,:)=0.0_wp
-      zgammat2d(:,:)=0.0_wp
-      zgammas2d(:,:)=0.0_wp
-      !
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ! Crude approximation for pressure (but commonly used)
-            ! 1e-04 to convert from Pa to dBar
-            zpress(ji,jj)=grav*rau0*fsdepw(ji,jj,mikt(ji,jj))*1.e-04
-            !
-         END DO
+      CALL wrk_alloc( jpi,jpj, zfrz  , zgammat, zgammas  )
+      CALL wrk_alloc( jpi,jpj, zfwflx, zhtflx , zhtflx_b )
+
+      ! initialisation
+      zgammat(:,:) = rn_gammat0 ; zgammas (:,:) = rn_gammas0
+      zhtflx (:,:) = 0.0_wp     ; zhtflx_b(:,:) = 0.0_wp    
+      zfwflx (:,:) = 0.0_wp
+
+      ! compute ice shelf melting
+      nit = 1 ; lit = .TRUE.
+      DO WHILE ( lit )    ! maybe just a constant number of iteration as in blk_core is fine
+         SELECT CASE ( nn_isfblk )
+         CASE ( 1 )   !  ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006)
+            ! Calculate freezing temperature
+            CALL eos_fzp( stbl(:,:), zfrz(:,:), risfdep(:,:) )
+
+            ! compute gammat every where (2d)
+            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx)
+            
+            ! compute upward heat flux zhtflx and upward water flux zwflx
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zhtflx(ji,jj) =   zgammat(ji,jj)*rcp*rau0*(ttbl(ji,jj)-zfrz(ji,jj))
+                  zfwflx(ji,jj) = - zhtflx(ji,jj)/rlfusisf
+               END DO
+            END DO
+
+            ! Compute heat flux and upward fresh water flux
+            qisf  (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
+            fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
+
+         CASE ( 2 )  ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015)
+            ! compute gammat every where (2d)
+            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx)
+
+            ! compute upward heat flux zhtflx and upward water flux zwflx
+            ! Resolution of a 2d equation from equation 21, 22 and 23 to find Sb (Asay-Davis et al., 2015)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ! compute coeficient to solve the 2nd order equation
+                  zeps1 = rcp*rau0*zgammat(ji,jj)
+                  zeps2 = rlfusisf*rau0*zgammas(ji,jj)
+                  zeps3 = rhoisf*rcpi*rkappa/MAX(risfdep(ji,jj),zeps)
+                  zeps4 = zlamb2+zlamb3*risfdep(ji,jj)
+                  zeps6 = zeps4-ttbl(ji,jj)
+                  zeps7 = zeps4-tsurf
+                  zaqe  = zlamb1 * (zeps1 + zeps3)
+                  zaqer = 0.5_wp/MIN(zaqe,-zeps)
+                  zbqe  = zeps1*zeps6+zeps3*zeps7-zeps2
+                  zcqe  = zeps2*stbl(ji,jj)
+                  zdis  = zbqe*zbqe-4.0_wp*zaqe*zcqe               
+
+                  ! Presumably zdis can never be negative because gammas is very small compared to gammat
+                  ! compute s freeze
+                  zsfrz=(-zbqe-SQRT(zdis))*zaqer
+                  IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer
+
+                  ! compute t freeze (eq. 22)
+                  zfrz(ji,jj)=zeps4+zlamb1*zsfrz
+  
+                  ! zfwflx is upward water flux
+                  ! zhtflx is upward heat flux (out of ocean)
+                  ! compute the upward water and heat flux (eq. 28 and eq. 29)
+                  zfwflx(ji,jj) = rau0 * zgammas(ji,jj) * (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps)
+                  zhtflx(ji,jj) = zgammat(ji,jj) * rau0 * rcp * (ttbl(ji,jj) - zfrz(ji,jj) ) 
+               END DO
+            END DO
+
+            ! compute heat and water flux
+            qisf  (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
+            fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
+
+         END SELECT
+
+         ! define if we need to iterate (nn_gammablk 0/1 do not need iteration)
+         IF ( nn_gammablk <  2 ) THEN ; lit = .FALSE.
+         ELSE                           
+            ! check total number of iteration
+            IF (nit >= 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
+            ELSE                 ; nit = nit + 1
+            END IF
+
+            ! compute error between 2 iterations
+            ! if needed save gammat and compute zhtflx_b for next iteration
+            zerr = MAXVAL(ABS(zhtflx-zhtflx_b))
+            IF ( zerr <= 0.01_wp ) THEN ; lit = .FALSE.
+            ELSE                        ; zhtflx_b(:,:) = zhtflx(:,:)
+            END IF
+         END IF
       END DO
-
-! Calculate in-situ temperature (ref to surface)
-      zti(:,:)=tinsitu( ttbl, stbl, zpress )
-! Calculate freezing temperature
-      CALL eos_fzp( sss_m(:,:), zfrz(:,:), zpress )
-
-      
-      zhtflx=0._wp ; zfwflx=0._wp
-      IF (nn_isfblk == 1) THEN
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF (mikt(ji,jj) > 1 ) THEN
-                  nit = 1; lit = .TRUE.; zgammat=rn_gammat0; zgammas=rn_gammas0; zhtflx_b=0._wp
-                  DO WHILE ( lit )
-! compute gamma
-                     CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, ji, jj, lit)
-! zhtflx is upward heat flux (out of ocean)
-                     zhtflx = zgammat*rcp*rau0*(zti(ji,jj)-zfrz(ji,jj))
-! zwflx is upward water flux
-                     zfwflx = - zhtflx/lfusisf
-! test convergence and compute gammat
-                     IF ( (zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE.
-
-                     nit = nit + 1
-                     IF (nit .GE. 100) CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
-
-! save gammat and compute zhtflx_b
-                     zgammat2d(ji,jj)=zgammat
-                     zhtflx_b = zhtflx
-                  END DO
-
-                  qisf(ji,jj) = - zhtflx
-! For genuine ISOMIP protocol this should probably be something like
-                  fwfisf(ji,jj) = zfwflx  * ( soce / MAX(stbl(ji,jj),zeps))
-               ELSE
-                  fwfisf(ji,jj) = 0._wp
-                  qisf(ji,jj)   = 0._wp
-               END IF
-            !
-            END DO
-         END DO
-
-      ELSE IF (nn_isfblk == 2 ) THEN
-
-! More complicated 3 equation thermodynamics as in MITgcm
-         DO jj = 2, jpj
-            DO ji = 2, jpi
-               IF (mikt(ji,jj) > 1 ) THEN
-                  nit=1; lit=.TRUE.; zgammat=rn_gammat0; zgammas=rn_gammas0; zhtflx_b=0._wp; zhtflx=0._wp
-                  DO WHILE ( lit )
-                     CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, ji, jj, lit)
-
-                     zeps1=rcp*rau0*zgammat
-                     zeps2=lfusisf*rau0*zgammas
-                     zeps3=rhoisf*rcpi*kappa/risfdep(ji,jj)
-                     zeps4=zlamb2+zlamb3*risfdep(ji,jj)
-                     zeps6=zeps4-zti(ji,jj)
-                     zeps7=zeps4-tsurf
-                     zaqe=zlamb1 * (zeps1 + zeps3)
-                     zaqer=0.5/zaqe
-                     zbqe=zeps1*zeps6+zeps3*zeps7-zeps2
-                     zcqe=zeps2*stbl(ji,jj)
-                     zdis=zbqe*zbqe-4.0*zaqe*zcqe               
-! Presumably zdis can never be negative because gammas is very small compared to gammat
-                     zsfrz=(-zbqe-SQRT(zdis))*zaqer
-                     IF (zsfrz .lt. 0.0) zsfrz=(-zbqe+SQRT(zdis))*zaqer
-                     zfrz(ji,jj)=zeps4+zlamb1*zsfrz
-  
-! zfwflx is upward water flux
-                     zfwflx= rau0 * zgammas * ( (zsfrz-stbl(ji,jj)) / zsfrz )
-! zhtflx is upward heat flux (out of ocean)
-! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value
-                     zhtflx = ( zgammat*rau0 - zcfac*zfwflx ) * rcp * (zti(ji,jj) - zfrz(ji,jj) ) 
-! zwflx is upward water flux
-! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz
-                     zfwflx = ( zgammas*rau0 - zcfac*zfwflx ) * (zsfrz - stbl(ji,jj)) / stbl(ji,jj)
-! test convergence and compute gammat
-                     IF (( zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE.
-
-                     nit = nit + 1
-                     IF (nit .GE. 51) THEN
-                        WRITE(numout,*) "sbcisf : too many iteration ... ", &
-                            &  zhtflx, zhtflx_b, zgammat, zgammas, nn_gammablk, ji, jj, mikt(ji,jj), narea
-                        CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
-                     END IF
-! save gammat and compute zhtflx_b
-                     zgammat2d(ji,jj)=zgammat
-                     zgammas2d(ji,jj)=zgammas
-                     zhtflx_b = zhtflx
-
-                  END DO
-! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value
-                  qisf(ji,jj) = - zhtflx 
-! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz
-                  fwfisf(ji,jj) = zfwflx 
-               ELSE
-                  fwfisf(ji,jj) = 0._wp
-                  qisf(ji,jj)   = 0._wp
-               ENDIF
-               !
-            END DO
-         END DO
-      ENDIF
-      ! lbclnk
-      CALL lbc_lnk(zgammas2d(:,:),'T',1.)
-      CALL lbc_lnk(zgammat2d(:,:),'T',1.)
-      ! output
-      CALL iom_put('isfgammat', zgammat2d)
-      CALL iom_put('isfgammas', zgammas2d)
-      !
-      CALL wrk_dealloc( jpi,jpj, zfrz,zpress,zti, zgammat2d, zgammas2d )
+      !
+      CALL iom_put('isfgammat', zgammat)
+      CALL iom_put('isfgammas', zgammas)
+      ! 
+      CALL wrk_dealloc( jpi,jpj, zfrz  , zgammat, zgammas  )
+      CALL wrk_dealloc( jpi,jpj, zfwflx, zhtflx , zhtflx_b )
       !
       IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_cav')
@@ -565,6 +536,5 @@
    END SUBROUTINE sbc_isf_cav
 
-
-   SUBROUTINE sbc_isf_gammats(gt, gs, zqhisf, zqwisf, ji, jj, lit )
+   SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf )
       !!----------------------------------------------------------------------
       !! ** Purpose    : compute the coefficient echange for heat flux
@@ -575,10 +545,10 @@
       !!                Jenkins et al., 2010, JPO, p2298-2312
       !!---------------------------------------------------------------------
-      REAL(wp), INTENT(inout) :: gt, gs, zqhisf, zqwisf
-      INTEGER , INTENT(in)    :: ji,jj
-      LOGICAL , INTENT(inout) :: lit
-
-      INTEGER  :: ikt                 ! loop index
-      REAL(wp) :: zut, zvt, zustar           ! U, V at T point and friction velocity
+      REAL(wp), DIMENSION(:,:), INTENT(out) :: pgt, pgs
+      REAL(wp), DIMENSION(:,:), INTENT(in ) :: pqhisf, pqwisf
+      !
+      INTEGER  :: ikt                        
+      INTEGER  :: ji, jj                     ! loop index
+      REAL(wp), DIMENSION(:,:), POINTER :: zustar           ! U, V at T point and friction velocity
       REAL(wp) :: zdku, zdkv                 ! U, V shear 
       REAL(wp) :: zPr, zSc, zRc              ! Prandtl, Scmidth and Richardson number 
@@ -590,169 +560,211 @@
       REAL(wp) :: zcoef                      ! temporary coef
       REAL(wp) :: zdep
-      REAL(wp), PARAMETER :: zxsiN = 0.052   ! dimensionless constant
-      REAL(wp), PARAMETER :: epsln = 1.0e-20 ! a small positive number
-      REAL(wp), PARAMETER :: znu   = 1.95e-6 ! kinamatic viscosity of sea water (m2.s-1)
-      REAL(wp) ::   rcs      = 1.0e-3_wp        ! conversion: mm/s ==> m/s
+      REAL(wp) :: zeps = 1.0e-20_wp    
+      REAL(wp), PARAMETER :: zxsiN = 0.052_wp   ! dimensionless constant
+      REAL(wp), PARAMETER :: znu   = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1)
       REAL(wp), DIMENSION(2) :: zts, zab
       !!---------------------------------------------------------------------
-      !
-      IF( nn_gammablk == 0 ) THEN
-      !! gamma is constant (specified in namelist)
-         gt = rn_gammat0
-         gs = rn_gammas0
-         lit = .FALSE.
-      ELSE IF ( nn_gammablk == 1 ) THEN
-      !! gamma is assume to be proportional to u* 
-      !! WARNING in case of Losh 2008 tbl parametrization, 
-      !! you have to used the mean value of u in the boundary layer) 
-      !! not yet coded
-      !! Jenkins et al., 2010, JPO, p2298-2312
-         ikt = mikt(ji,jj)
-      !! Compute U and V at T points
-   !      zut = 0.5 * ( utbl(ji-1,jj  ) + utbl(ji,jj) )
-   !      zvt = 0.5 * ( vtbl(ji  ,jj-1) + vtbl(ji,jj) )
-          zut = utbl(ji,jj)
-          zvt = vtbl(ji,jj)
-
-      !! compute ustar
-         zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) )
-      !! Compute mean value over the TBL
-
-      !! Compute gammats
-         gt = zustar * rn_gammat0
-         gs = zustar * rn_gammas0
-         lit = .FALSE.
-      ELSE IF ( nn_gammablk == 2 ) THEN
-      !! gamma depends of stability of boundary layer
-      !! WARNING in case of Losh 2008 tbl parametrization, 
-      !! you have to used the mean value of u in the boundary layer) 
-      !! not yet coded
-      !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
-      !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO)
+      CALL wrk_alloc( jpi,jpj, zustar )
+      !
+      SELECT CASE ( nn_gammablk )
+      CASE ( 0 ) ! gamma is constant (specified in namelist)
+         !! ISOMIP formulation (Hunter et al, 2006)
+         pgt(:,:) = rn_gammat0
+         pgs(:,:) = rn_gammas0
+
+      CASE ( 1 ) ! gamma is assume to be proportional to u*
+         !! Jenkins et al., 2010, JPO, p2298-2312
+         !! Adopted by Asay-Davis et al. (2015)
+
+         !! compute ustar (eq. 24)
+         zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) )
+
+         !! Compute gammats
+         pgt(:,:) = zustar(:,:) * rn_gammat0
+         pgs(:,:) = zustar(:,:) * rn_gammas0
+      
+      CASE ( 2 ) ! gamma depends of stability of boundary layer
+         !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
+         !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO)
+         !! compute ustar
+         zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) )
+
+         !! compute Pr and Sc number (can be improved)
+         zPr =   13.8_wp
+         zSc = 2432.0_wp
+
+         !! compute gamma mole
+         zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp
+         zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp
+
+         !! compute gamma
+         DO ji=2,jpi
+            DO jj=2,jpj
                ikt = mikt(ji,jj)
 
-      !! Compute U and V at T points
-               zut = 0.5 * ( utbl(ji-1,jj  ) + utbl(ji,jj) )
-               zvt = 0.5 * ( vtbl(ji  ,jj-1) + vtbl(ji,jj) )
-
-      !! compute ustar
-               zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) )
-               IF (zustar == 0._wp) THEN           ! only for kt = 1 I think
-                 gt = rn_gammat0
-                 gs = rn_gammas0
+               IF (zustar(ji,jj) == 0._wp) THEN           ! only for kt = 1 I think
+                  pgt = rn_gammat0
+                  pgs = rn_gammas0
                ELSE
-      !! compute Rc number (as done in zdfric.F90)
-               zcoef = 0.5 / fse3w(ji,jj,ikt)
-               !                                            ! shear of horizontal velocity
-               zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )   &
-                  &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  )
-               zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )   &
-                  &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  )
-               !                                            ! richardson number (minimum value set to zero)
-               zRc = rn2(ji,jj,ikt+1) / ( zdku*zdku + zdkv*zdkv + 1.e-20 )
-
-      !! compute Pr and Sc number (can be improved)
-               zPr =   13.8
-               zSc = 2432.0
-
-      !! compute gamma mole
-               zgmolet = 12.5 * zPr ** (2.0/3.0) - 6.0
-               zgmoles = 12.5 * zSc ** (2.0/3.0) -6.0
-
-      !! compute bouyancy 
-               zts(jp_tem) = ttbl(ji,jj)
-               zts(jp_sal) = stbl(ji,jj)
-               zdep        = fsdepw(ji,jj,ikt)
-               !
-               CALL eos_rab( zts, zdep, zab )
+                  !! compute Rc number (as done in zdfric.F90)
+                  zcoef = 0.5_wp / e3w_n(ji,jj,ikt)
+                  !                                            ! shear of horizontal velocity
+                  zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )  &
+                     &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  )
+                  zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )  &
+                     &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  )
+                  !                                            ! richardson number (minimum value set to zero)
+                  zRc = rn2(ji,jj,ikt+1) / MAX( zdku*zdku + zdkv*zdkv, zeps )
+
+                  !! compute bouyancy 
+                  zts(jp_tem) = ttbl(ji,jj)
+                  zts(jp_sal) = stbl(ji,jj)
+                  zdep        = gdepw_n(ji,jj,ikt)
                   !
-      !! compute length scale 
-               zbuofdep = grav * ( zab(jp_tem) * zqhisf - zab(jp_sal) * zqwisf )  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-      !! compute Monin Obukov Length
-               ! Maximum boundary layer depth
-               zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) -0.001
-               ! Compute Monin obukhov length scale at the surface and Ekman depth:
-               zmob   = zustar ** 3 / (vkarmn * (zbuofdep + epsln))
-               zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
-
-      !! compute eta* (stability parameter)
-               zetastar = 1 / ( SQRT(1 + MAX(zxsiN * zustar / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0)))
-
-      !! compute the sublayer thickness
-               zhnu = 5 * znu / zustar
-      !! compute gamma turb
-               zgturb = 1/vkarmn * LOG(zustar * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) &
-               &      + 1 / ( 2 * zxsiN * zetastar ) - 1/vkarmn
-
-      !! compute gammats
-               gt = zustar / (zgturb + zgmolet)
-               gs = zustar / (zgturb + zgmoles)
+                  CALL eos_rab( zts, zdep, zab )
+                  !
+                  !! compute length scale 
+                  zbuofdep = grav * ( zab(jp_tem) * pqhisf(ji,jj) - zab(jp_sal) * pqwisf(ji,jj) )  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+                  !! compute Monin Obukov Length
+                  ! Maximum boundary layer depth
+                  zhmax = gdept_n(ji,jj,mbkt(ji,jj)) - gdepw_n(ji,jj,mikt(ji,jj)) -0.001_wp
+                  ! Compute Monin obukhov length scale at the surface and Ekman depth:
+                  zmob   = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps))
+                  zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
+
+                  !! compute eta* (stability parameter)
+                  zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0_wp)))
+
+                  !! compute the sublayer thickness
+                  zhnu = 5 * znu / zustar(ji,jj)
+
+                  !! compute gamma turb
+                  zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) &
+                  &      + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn
+
+                  !! compute gammats
+                  pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet)
+                  pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles)
                END IF
-      END IF
-      !
-   END SUBROUTINE
-
-
-   SUBROUTINE sbc_isf_tbl( varin, varout, cptin )
+            END DO
+         END DO
+         CALL lbc_lnk(pgt(:,:),'T',1.)
+         CALL lbc_lnk(pgs(:,:),'T',1.)
+      END SELECT
+      CALL wrk_dealloc( jpi,jpj, zustar )
+      !
+   END SUBROUTINE sbc_isf_gammats
+
+   SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE sbc_isf_tbl  ***
       !!
-      !! ** Purpose : compute mean T/S/U/V in the boundary layer 
-      !!
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:), INTENT(in) :: varin
-      REAL(wp), DIMENSION(:,:)  , INTENT(out):: varout
+      !! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:), INTENT( in  ) :: pvarin
+      REAL(wp), DIMENSION(:,:)  , INTENT( out ) :: pvarout
+      CHARACTER(len=1),           INTENT( in  ) :: cd_ptin ! point of variable in/out
+      !
+      REAL(wp) :: ze3, zhk
+      REAL(wp), DIMENSION(:,:), POINTER :: zhisf_tbl ! thickness of the tbl
+
+      INTEGER :: ji, jj, jk                  ! loop index
+      INTEGER :: ikt, ikb                    ! top and bottom index of the tbl
+      !!----------------------------------------------------------------------
+      ! allocation
+      CALL wrk_alloc( jpi,jpj, zhisf_tbl)
       
-      CHARACTER(len=1), INTENT(in) :: cptin ! point of variable in/out
-
-      REAL(wp) :: ze3, zhk
-      REAL(wp), DIMENSION(:,:), POINTER :: zikt
-
-      INTEGER :: ji,jj,jk
-      INTEGER :: ikt,ikb
-      INTEGER, DIMENSION(:,:), POINTER :: mkt, mkb
-
-      CALL wrk_alloc( jpi,jpj, mkt, mkb  )
-      CALL wrk_alloc( jpi,jpj, zikt )
-
-      ! get first and last level of tbl
-      mkt(:,:) = misfkt(:,:)
-      mkb(:,:) = misfkb(:,:)
-
-      varout(:,:)=0._wp
-      DO jj = 2,jpj
-         DO ji = 2,jpi
-            IF (ssmask(ji,jj) == 1) THEN
-               ikt = mkt(ji,jj)
-               ikb = mkb(ji,jj)
+      ! initialisation
+      pvarout(:,:)=0._wp
+   
+      SELECT CASE ( cd_ptin )
+      CASE ( 'U' ) ! compute U in the top boundary layer at T- point 
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               ikt = miku(ji,jj) ; ikb = miku(ji,jj)
+               ! thickness of boundary layer at least the top level thickness
+               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3u_n(ji,jj,ikt))
+
+               ! determine the deepest level influenced by the boundary layer
+               DO jk = ikt+1, mbku(ji,jj)
+                  IF ( (SUM(e3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk
+               END DO
+               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3u_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
 
                ! level fully include in the ice shelf boundary layer
                DO jk = ikt, ikb - 1
-                  ze3 = fse3t_n(ji,jj,jk)
-                  IF (cptin == 'T' ) varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3
-                  IF (cptin == 'U' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji-1,jj,jk)) &
-                     &                                                       * r1_hisf_tbl(ji,jj) * ze3
-                  IF (cptin == 'V' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji,jj-1,jk)) &
-                     &                                                       * r1_hisf_tbl(ji,jj) * ze3
+                  ze3 = e3u_n(ji,jj,jk)
+                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
                END DO
 
                ! level partially include in ice shelf boundary layer 
-               zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)
-               IF (cptin == 'T') &
-                   &  varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk)
-               IF (cptin == 'U') &
-                   &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji-1,jj,ikb)) * (1._wp - zhk)
-               IF (cptin == 'V') &
-                   &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji,jj-1,ikb)) * (1._wp - zhk)
-            END IF
-         END DO
-      END DO
-
-      CALL wrk_dealloc( jpi,jpj, mkt, mkb )      
-      CALL wrk_dealloc( jpi,jpj, zikt ) 
-
-      IF (cptin == 'T') CALL lbc_lnk(varout,'T',1.)
-      IF (cptin == 'U' .OR. cptin == 'V') CALL lbc_lnk(varout,'T',-1.)
+               zhk = SUM( e3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
+               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
+            END DO
+         END DO
+         DO jj = 2,jpj
+            DO ji = 2,jpi
+               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji-1,jj))
+            END DO
+         END DO
+         CALL lbc_lnk(pvarout,'T',-1.)
+      
+      CASE ( 'V' ) ! compute V in the top boundary layer at T- point 
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               ikt = mikv(ji,jj) ; ikb = mikv(ji,jj)
+               ! thickness of boundary layer at least the top level thickness
+               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3v_n(ji,jj,ikt))
+
+               ! determine the deepest level influenced by the boundary layer
+               DO jk = ikt+1, mbkv(ji,jj)
+                  IF ( (SUM(e3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk
+               END DO
+               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3v_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
+
+               ! level fully include in the ice shelf boundary layer
+               DO jk = ikt, ikb - 1
+                  ze3 = e3v_n(ji,jj,jk)
+                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
+               END DO
+
+               ! level partially include in ice shelf boundary layer 
+               zhk = SUM( e3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
+               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
+            END DO
+         END DO
+         DO jj = 2,jpj
+            DO ji = 2,jpi
+               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji,jj-1))
+            END DO
+         END DO
+         CALL lbc_lnk(pvarout,'T',-1.)
+
+      CASE ( 'T' ) ! compute T in the top boundary layer at T- point 
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               ikt = misfkt(ji,jj)
+               ikb = misfkb(ji,jj)
+
+               ! level fully include in the ice shelf boundary layer
+               DO jk = ikt, ikb - 1
+                  ze3 = e3t_n(ji,jj,jk)
+                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3
+               END DO
+
+               ! level partially include in ice shelf boundary layer 
+               zhk = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)
+               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
+            END DO
+         END DO
+      END SELECT
+
+      ! mask mean tbl value
+      pvarout(:,:) = pvarout(:,:) * ssmask(:,:)
+
+      ! deallocation
+      CALL wrk_dealloc( jpi,jpj, zhisf_tbl )      
       !
    END SUBROUTINE sbc_isf_tbl
@@ -771,16 +783,15 @@
       !! ** Action  :   phdivn   decreased by the runoff inflow
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
-      !!
+      REAL(wp), DIMENSION(:,:,:), INTENT( inout ) ::   phdivn   ! horizontal divergence
+      ! 
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       INTEGER  ::   ikt, ikb 
-      INTEGER  ::   nk_isf
-      REAL(wp)     ::   zhk, z1_hisf_tbl, zhisf_tbl
-      REAL(wp)     ::   zfact     ! local scalar
+      REAL(wp) ::   zhk
+      REAL(wp) ::   zfact     ! local scalar
       !!----------------------------------------------------------------------
       !
       zfact   = 0.5_wp
       !
-      IF (lk_vvl) THEN     ! need to re compute level distribution of isf fresh water
+      IF(.NOT.ln_linssh ) THEN     ! need to re compute level distribution of isf fresh water
          DO jj = 1,jpj
             DO ji = 1,jpi
@@ -788,21 +799,21 @@
                ikb = misfkt(ji,jj)
                ! thickness of boundary layer at least the top level thickness
-               rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t(ji,jj,ikt))
+               rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt))
 
                ! determine the deepest level influenced by the boundary layer
-               ! 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
-               END DO
-               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
+                  IF ( (SUM(e3t_n(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(e3t_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
                misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl
                r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
 
-               zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1
-               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
-            END DO
-         END DO
-      END IF ! vvl case
-      !
+               zhk           = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1
+               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
+            END DO
+         END DO
+      END IF 
+      !
+      !==   ice shelf melting distributed over several levels   ==!
       DO jj = 1,jpj
          DO ji = 1,jpi
@@ -812,91 +823,13 @@
                DO jk = ikt, ikb - 1
                   phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) &
-                    &               * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact
+                    &              * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact
                END DO
                ! level partially include in ice shelf boundary layer 
                phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) &
-                  &             + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 
-            !==   ice shelf melting mass distributed over several levels   ==!
+                    &            + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 
          END DO
       END DO
       !
    END SUBROUTINE sbc_isf_div
-
-
-   FUNCTION tinsitu( ptem, psal, ppress ) RESULT( pti )
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE eos_init  ***
-      !!
-      !! ** Purpose :   Compute the in-situ temperature [Celcius]
-      !!
-      !! ** Method  :   
-      !!
-      !! Reference  :   Bryden,h.,1973,deep-sea res.,20,401-408
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptem   ! potential temperature [Celcius]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity             [psu]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ppress ! pressure             [dBar]
-      REAL(wp), DIMENSION(:,:), POINTER           ::   pti    ! in-situ temperature [Celcius]
-!      REAL(wp) :: fsatg
-!      REAL(wp) :: pfps, pfpt, pfphp 
-      REAL(wp) :: zt, zs, zp, zh, zq, zxk
-      INTEGER  :: ji, jj            ! dummy loop indices
-      !
-      CALL wrk_alloc( jpi,jpj, pti  )
-      ! 
-      DO jj=1,jpj
-         DO ji=1,jpi
-            zh = ppress(ji,jj)
-! Theta1
-            zt = ptem(ji,jj)
-            zs = psal(ji,jj)
-            zp = 0.0
-            zxk= zh * fsatg( zs, zt, zp )
-            zt = zt + 0.5 * zxk
-            zq = zxk
-! Theta2
-            zp = zp + 0.5 * zh
-            zxk= zh*fsatg( zs, zt, zp )
-            zt = zt + 0.29289322 * ( zxk - zq )
-            zq = 0.58578644 * zxk + 0.121320344 * zq
-! Theta3
-            zxk= zh * fsatg( zs, zt, zp )
-            zt = zt + 1.707106781 * ( zxk - zq )
-            zq = 3.414213562 * zxk - 4.121320344 * zq
-! Theta4
-            zp = zp + 0.5 * zh
-            zxk= zh * fsatg( zs, zt, zp )
-            pti(ji,jj) = zt + ( zxk - 2.0 * zq ) / 6.0
-         END DO
-      END DO
-      !
-      CALL wrk_dealloc( jpi,jpj, pti  )
-      !
-   END FUNCTION tinsitu
-
-
-   FUNCTION fsatg( pfps, pfpt, pfphp )
-      !!----------------------------------------------------------------------
-      !!                 ***  FUNCTION fsatg  ***
-      !!
-      !! ** Purpose    :   Compute the Adiabatic laspse rate [Celcius].[decibar]^-1
-      !!
-      !! ** Reference  :   Bryden,h.,1973,deep-sea res.,20,401-408
-      !! 
-      !! ** units      :   pressure        pfphp    decibars
-      !!                   temperature     pfpt     deg celsius (ipts-68)
-      !!                   salinity        pfps     (ipss-78)
-      !!                   adiabatic       fsatg    deg. c/decibar
-      !!----------------------------------------------------------------------
-      REAL(wp) :: pfps, pfpt, pfphp 
-      REAL(wp) :: fsatg
-      !
-      fsatg = (((-2.1687e-16*pfpt+1.8676e-14)*pfpt-4.6206e-13)*pfphp         &
-        &    +((2.7759e-12*pfpt-1.1351e-10)*(pfps-35.)+((-5.4481e-14*pfpt    &
-        &    +8.733e-12)*pfpt-6.7795e-10)*pfpt+1.8741e-8))*pfphp             &
-        &    +(-4.2393e-8*pfpt+1.8932e-6)*(pfps-35.)                         &
-        &    +((6.6228e-10*pfpt-6.836e-8)*pfpt+8.5258e-6)*pfpt+3.5803e-5
-      !
-    END FUNCTION fsatg
-    !!======================================================================
+   !!======================================================================
 END MODULE sbcisf
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 6121)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 6140)
@@ -17,40 +17,44 @@
 
    !!----------------------------------------------------------------------
-   !!   sbc_init       : read namsbc namelist
-   !!   sbc            : surface ocean momentum, heat and freshwater boundary conditions
+   !!   sbc_init      : read namsbc namelist
+   !!   sbc           : surface ocean momentum, heat and freshwater boundary conditions
    !!----------------------------------------------------------------------
-   USE oce              ! ocean dynamics and tracers
-   USE dom_oce          ! ocean space and time domain
-   USE phycst           ! physical constants
-   USE sbc_oce          ! Surface boundary condition: ocean fields
-   USE trc_oce          ! shared ocean-passive tracers variables
-   USE sbc_ice          ! Surface boundary condition: ice fields
-   USE sbcdcy           ! surface boundary condition: diurnal cycle
-   USE sbcssm           ! surface boundary condition: sea-surface mean variables
-   USE sbcana           ! surface boundary condition: analytical formulation
-   USE sbcflx           ! surface boundary condition: flux formulation
-   USE sbcblk_clio      ! surface boundary condition: bulk formulation : CLIO
-   USE sbcblk_core      ! surface boundary condition: bulk formulation : CORE
-   USE sbcblk_mfs       ! surface boundary condition: bulk formulation : MFS
-   USE sbcice_if        ! surface boundary condition: ice-if sea-ice model
-   USE sbcice_lim       ! surface boundary condition: LIM 3.0 sea-ice model
-   USE sbcice_lim_2     ! surface boundary condition: LIM 2.0 sea-ice model
-   USE sbcice_cice      ! surface boundary condition: CICE    sea-ice model
-   USE sbccpl           ! surface boundary condition: coupled florulation
-   USE cpl_oasis3       ! OASIS routines for coupling
-   USE sbcssr           ! surface boundary condition: sea surface restoring
-   USE sbcrnf           ! surface boundary condition: runoffs
-   USE sbcisf           ! surface boundary condition: ice shelf
-   USE sbcfwb           ! surface boundary condition: freshwater budget
-   USE closea           ! closed sea
-   USE icbstp           ! Icebergs!
-
-   USE prtctl           ! Print control                    (prt_ctl routine)
-   USE iom              ! IOM library
-   USE in_out_manager   ! I/O manager
-   USE lib_mpp          ! MPP library
-   USE timing           ! Timing
-   USE sbcwave          ! Wave module
-   USE bdy_par          ! Require lk_bdy
+   USE oce            ! ocean dynamics and tracers
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   USE sbc_oce        ! Surface boundary condition: ocean fields
+   USE trc_oce        ! shared ocean-passive tracers variables
+   USE sbc_ice        ! Surface boundary condition: ice fields
+   USE sbcdcy         ! surface boundary condition: diurnal cycle
+   USE sbcssm         ! surface boundary condition: sea-surface mean variables
+   USE sbcana         ! surface boundary condition: analytical formulation
+   USE sbcflx         ! surface boundary condition: flux formulation
+   USE sbcblk_clio    ! surface boundary condition: bulk formulation : CLIO
+   USE sbcblk_core    ! surface boundary condition: bulk formulation : CORE
+   USE sbcblk_mfs     ! surface boundary condition: bulk formulation : MFS
+   USE sbcice_if      ! surface boundary condition: ice-if sea-ice model
+   USE sbcice_lim     ! surface boundary condition: LIM 3.0 sea-ice model
+   USE sbcice_lim_2   ! surface boundary condition: LIM 2.0 sea-ice model
+   USE sbcice_cice    ! surface boundary condition: CICE    sea-ice model
+   USE sbccpl         ! surface boundary condition: coupled florulation
+   USE cpl_oasis3     ! OASIS routines for coupling
+   USE sbcssr         ! surface boundary condition: sea surface restoring
+   USE sbcrnf         ! surface boundary condition: runoffs
+   USE sbcisf         ! surface boundary condition: ice shelf
+   USE sbcfwb         ! surface boundary condition: freshwater budget
+   USE closea         ! closed sea
+   USE icbstp         ! Icebergs
+   USE traqsr         ! active tracers: light penetration
+   USE sbcwave        ! Wave module
+   USE bdy_par        ! Require lk_bdy
+   !
+   USE prtctl         ! Print control                    (prt_ctl routine)
+   USE iom            ! IOM library
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   USE timing         ! Timing
+
+   USE diurnal_bulk, ONLY: &
+      & ln_diurnal_only 
 
    IMPLICIT NONE
@@ -62,6 +66,4 @@
    INTEGER ::   nsbc   ! type of surface boundary condition (deduced from namsbc informations)
       
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO-consortium (2011) 
@@ -85,13 +87,16 @@
       INTEGER ::   icpt   ! local integer
       !!
-      NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx, ln_blk_clio, ln_blk_core, ln_mixcpl,   &
-         &             ln_blk_mfs, ln_apr_dyn, nn_ice, nn_ice_embd, ln_dm2dc   , ln_rnf   ,   &
-         &             ln_ssr    , nn_isf    , nn_fwb, ln_cdgw    , ln_wave    , ln_sdw   ,   &
-         &             nn_lsm    , nn_limflx , nn_components, ln_cpl
+      NAMELIST/namsbc/ nn_fsbc  , ln_ana   , ln_flx, ln_blk_clio, ln_blk_core, ln_blk_mfs,   &
+         &             ln_cpl   , ln_mixcpl, nn_components      , nn_limflx  ,               &
+         &             ln_traqsr, ln_dm2dc ,                                                 &  
+         &             nn_ice   , nn_ice_embd,                                               &
+         &             ln_rnf   , ln_ssr   , ln_isf   , nn_fwb    , ln_apr_dyn,              &
+         &             ln_wave  ,                                                            &
+         &             nn_lsm   
       INTEGER  ::   ios
       INTEGER  ::   ierr, ierr0, ierr1, ierr2, ierr3, jpm
       LOGICAL  ::   ll_purecpl
       !!----------------------------------------------------------------------
-
+      !
       IF(lwp) THEN
          WRITE(numout,*)
@@ -99,14 +104,14 @@
          WRITE(numout,*) '~~~~~~~~ '
       ENDIF
-
+      !
       REWIND( numnam_ref )              ! Namelist namsbc in reference namelist : Surface boundary
       READ  ( numnam_ref, namsbc, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc in reference namelist', lwp )
-
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc in reference namelist', lwp )
+      !
       REWIND( numnam_cfg )              ! Namelist namsbc in configuration namelist : Parameters of the run
       READ  ( numnam_cfg, namsbc, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc in configuration namelist', lwp )
-      IF(lwm) WRITE ( numond, namsbc )
-
+902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc in configuration namelist', lwp )
+      IF(lwm) WRITE( numond, namsbc )
+      !
       !                          ! overwrite namelist parameter using CPP key information
       IF( Agrif_Root() ) THEN                ! AGRIF zoom
@@ -119,52 +124,54 @@
           nn_ice      =   0
       ENDIF
-
+      !
       IF(lwp) THEN               ! Control print
          WRITE(numout,*) '        Namelist namsbc (partly overwritten with CPP key setting)'
          WRITE(numout,*) '           frequency update of sbc (and ice)             nn_fsbc     = ', nn_fsbc
-         WRITE(numout,*) '           Type of sbc : '
-         WRITE(numout,*) '              analytical formulation                     ln_ana      = ', ln_ana
-         WRITE(numout,*) '              flux       formulation                     ln_flx      = ', ln_flx
-         WRITE(numout,*) '              CLIO bulk  formulation                     ln_blk_clio = ', ln_blk_clio
-         WRITE(numout,*) '              CORE bulk  formulation                     ln_blk_core = ', ln_blk_core
-         WRITE(numout,*) '              MFS  bulk  formulation                     ln_blk_mfs  = ', ln_blk_mfs
-         WRITE(numout,*) '              ocean-atmosphere coupled formulation       ln_cpl      = ', ln_cpl
-         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,*) '              Multicategory heat flux formulation (LIM3) nn_limflx   = ', nn_limflx
+         WRITE(numout,*) '           Type of air-sea fluxes : '
+         WRITE(numout,*) '              analytical formulation                     ln_ana        = ', ln_ana
+         WRITE(numout,*) '              flux       formulation                     ln_flx        = ', ln_flx
+         WRITE(numout,*) '              CLIO bulk  formulation                     ln_blk_clio   = ', ln_blk_clio
+         WRITE(numout,*) '              CORE bulk  formulation                     ln_blk_core   = ', ln_blk_core
+         WRITE(numout,*) '              MFS  bulk  formulation                     ln_blk_mfs    = ', ln_blk_mfs
+         WRITE(numout,*) '           Type of coupling (Ocean/Ice/Atmosphere) : '
+         WRITE(numout,*) '              ocean-atmosphere coupled formulation       ln_cpl        = ', ln_cpl
+         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,*) '              Multicategory heat flux formulation (LIM3) nn_limflx     = ', nn_limflx
+         WRITE(numout,*) '           Sea-ice : '
+         WRITE(numout,*) '              ice management in the sbc (=0/1/2/3)       nn_ice        = ', nn_ice 
+         WRITE(numout,*) '              ice-ocean embedded/levitating (=0/1/2)     nn_ice_embd   = ', nn_ice_embd
          WRITE(numout,*) '           Misc. options of sbc : '
-         WRITE(numout,*) '              Patm gradient added in ocean & ice Eqs.    ln_apr_dyn  = ', ln_apr_dyn
-         WRITE(numout,*) '              ice management in the sbc (=0/1/2/3)       nn_ice      = ', nn_ice 
-         WRITE(numout,*) '              ice-ocean embedded/levitating (=0/1/2)     nn_ice_embd = ', nn_ice_embd
-         WRITE(numout,*) '              daily mean to diurnal cycle qsr            ln_dm2dc    = ', ln_dm2dc 
-         WRITE(numout,*) '              runoff / runoff mouths                     ln_rnf      = ', ln_rnf
-         WRITE(numout,*) '              iceshelf formulation                       nn_isf      = ', nn_isf
-         WRITE(numout,*) '              Sea Surface Restoring on SST and/or SSS    ln_ssr      = ', ln_ssr
-         WRITE(numout,*) '              FreshWater Budget control  (=0/1/2)        nn_fwb      = ', nn_fwb
-         WRITE(numout,*) '              closed sea (=0/1) (set in namdom)          nn_closea   = ', nn_closea
-         WRITE(numout,*) '              n. of iterations if land-sea-mask applied  nn_lsm      = ', nn_lsm
-      ENDIF
-
-      ! LIM3 Multi-category heat flux formulation
-      SELECT CASE ( nn_limflx)
-      CASE ( -1 )
-         IF(lwp) WRITE(numout,*) '              Use of per-category fluxes (nn_limflx = -1) '
-      CASE ( 0  )
-         IF(lwp) WRITE(numout,*) '              Average per-category fluxes (nn_limflx = 0) ' 
-      CASE ( 1  )
-         IF(lwp) WRITE(numout,*) '              Average then redistribute per-category fluxes (nn_limflx = 1) '
-      CASE ( 2  )
-         IF(lwp) WRITE(numout,*) '              Redistribute a single flux over categories (nn_limflx = 2) '
-      END SELECT
-      !
-      IF ( nn_components /= jp_iam_nemo .AND. .NOT. lk_oasis )   &
-         &      CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but key_oasis3 disabled' )
-      IF ( nn_components == jp_iam_opa .AND. ln_cpl )   &
-         &      CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but ln_cpl = T in OPA' )
-      IF ( nn_components == jp_iam_opa .AND. ln_mixcpl )   &
-         &      CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but ln_mixcpl = T in OPA' )
-      IF ( ln_cpl .AND. .NOT. lk_oasis )    &
-         &      CALL ctl_stop( 'STOP', 'sbc_init : OASIS-coupled atmosphere model, but key_oasis3 disabled' )
+         WRITE(numout,*) '              Light penetration in temperature Eq.       ln_traqsr     = ', ln_traqsr
+         WRITE(numout,*) '                 daily mean to diurnal cycle qsr            ln_dm2dc   = ', ln_dm2dc 
+         WRITE(numout,*) '              Sea Surface Restoring on SST and/or SSS    ln_ssr        = ', ln_ssr
+         WRITE(numout,*) '              FreshWater Budget control  (=0/1/2)        nn_fwb        = ', nn_fwb
+         WRITE(numout,*) '              Patm gradient added in ocean & ice Eqs.    ln_apr_dyn    = ', ln_apr_dyn
+         WRITE(numout,*) '              runoff / runoff mouths                     ln_rnf        = ', ln_rnf
+         WRITE(numout,*) '              iceshelf formulation                       ln_isf        = ', ln_isf
+         WRITE(numout,*) '              closed sea (=0/1) (set in namdom)          nn_closea     = ', nn_closea
+         WRITE(numout,*) '              nb of iterations if land-sea-mask applied  nn_lsm        = ', nn_lsm
+         WRITE(numout,*) '              surface wave                               ln_wave       = ', ln_wave  
+      ENDIF
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         SELECT CASE ( nn_limflx )        ! LIM3 Multi-category heat flux formulation
+         CASE ( -1 )   ;   WRITE(numout,*) '   LIM3: use per-category fluxes (nn_limflx = -1) '
+         CASE ( 0  )   ;   WRITE(numout,*) '   LIM3: use average per-category fluxes (nn_limflx = 0) ' 
+         CASE ( 1  )   ;   WRITE(numout,*) '   LIM3: use average then redistribute per-category fluxes (nn_limflx = 1) '
+         CASE ( 2  )   ;   WRITE(numout,*) '   LIM3: Redistribute a single flux over categories (nn_limflx = 2) '
+         END SELECT
+      ENDIF
+      !
+      IF( nn_components /= jp_iam_nemo .AND. .NOT. lk_oasis )   &
+         &      CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but key_oasis3 disabled' )
+      IF( nn_components == jp_iam_opa .AND. ln_cpl )   &
+         &      CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but ln_cpl = T in OPA' )
+      IF( nn_components == jp_iam_opa .AND. ln_mixcpl )   &
+         &      CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but ln_mixcpl = T in OPA' )
+      IF( ln_cpl .AND. .NOT. lk_oasis )    &
+         &      CALL ctl_stop( 'sbc_init : OASIS-coupled atmosphere model, but key_oasis3 disabled' )
       IF( ln_mixcpl .AND. .NOT. lk_oasis )    &
          &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl) requires the cpp key key_oasis3' )
@@ -178,18 +185,17 @@
 
       !                          ! Checks:
-      IF( nn_isf .EQ. 0 ) THEN                      ! variable initialisation if no ice shelf 
-         IF( sbc_isf_alloc() /= 0 )   CALL ctl_stop( 'sbc_init : unable to allocate sbc_isf arrays' )
+      IF( .NOT. ln_isf ) 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
          risf_tsc(:,:,:) = 0.0_wp ; risf_tsc_b(:,:,:) = 0.0_wp
-         rdivisf       = 0.0_wp
       END IF
-      IF( nn_ice == 0 .AND. nn_components /= jp_iam_opa )   fr_i(:,:) = 0.e0 ! no ice in the domain, ice fraction is always zero
-
-      sfx(:,:) = 0.0_wp                            ! the salt flux due to freezing/melting will be computed (i.e. will be non-zero) 
+      IF( nn_ice == 0 .AND. nn_components /= jp_iam_opa )   fr_i(:,:) = 0._wp    ! no ice in the domain, ice fraction is always zero
+
+      sfx(:,:) = 0._wp                             ! the salt flux due to freezing/melting will be computed (i.e. will be non-zero) 
                                                    ! only if sea-ice is present
  
-      fmmflx(:,:) = 0.0_wp                        ! freezing-melting array initialisation
+      fmmflx(:,:) = 0._wp                          ! freezing-melting array initialisation
       
-      taum(:,:) = 0.0_wp                           ! Initialise taum for use in gls in case of reduced restart
+      taum(:,:) = 0._wp                            ! Initialise taum for use in gls in case of reduced restart
 
       !                                            ! restartability   
@@ -214,17 +220,4 @@
          &   CALL ctl_stop( 'diurnal cycle into qsr field from daily values requires a flux or core-bulk formulation' )
       
-      IF ( ln_wave ) THEN
-      !Activated wave module but neither drag nor stokes drift activated
-         IF ( .NOT.(ln_cdgw .OR. ln_sdw) )   THEN
-            CALL ctl_warn( 'Ask for wave coupling but nor drag coefficient (ln_cdgw=F) neither stokes drift activated (ln_sdw=F)' )
-      !drag coefficient read from wave model definable only with mfs bulk formulae and core 
-         ELSEIF (ln_cdgw .AND. .NOT.(ln_blk_mfs .OR. ln_blk_core) )       THEN       
-             CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core')
-         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) ')
-      ENDIF 
       !                          ! Choice of the Surface Boudary Condition (set nsbc)
       ll_purecpl = ln_cpl .AND. .NOT. ln_mixcpl
@@ -245,25 +238,26 @@
       IF(lwp) THEN
          WRITE(numout,*)
-         IF( nsbc == jp_gyre    )   WRITE(numout,*) '              GYRE analytical formulation'
-         IF( nsbc == jp_ana     )   WRITE(numout,*) '              analytical formulation'
-         IF( nsbc == jp_flx     )   WRITE(numout,*) '              flux formulation'
-         IF( nsbc == jp_clio    )   WRITE(numout,*) '              CLIO bulk formulation'
-         IF( nsbc == jp_core    )   WRITE(numout,*) '              CORE bulk formulation'
-         IF( nsbc == jp_purecpl )   WRITE(numout,*) '              pure coupled formulation'
-         IF( nsbc == jp_mfs     )   WRITE(numout,*) '              MFS Bulk formulation'
-         IF( nsbc == jp_none    )   WRITE(numout,*) '              OPA coupled to SAS via oasis'
-         IF( ln_mixcpl          )   WRITE(numout,*) '              + forced-coupled mixed formulation'
+         SELECT CASE( nsbc )
+         CASE( jp_gyre    )   ;   WRITE(numout,*) '   GYRE analytical formulation'
+         CASE( jp_ana     )   ;   WRITE(numout,*) '   analytical formulation'
+         CASE( jp_flx     )   ;   WRITE(numout,*) '   flux formulation'
+         CASE( jp_clio    )   ;   WRITE(numout,*) '   CLIO bulk formulation'
+         CASE( jp_core    )   ;   WRITE(numout,*) '   CORE bulk formulation'
+         CASE( jp_purecpl )   ;   WRITE(numout,*) '   pure coupled formulation'
+         CASE( jp_mfs     )   ;   WRITE(numout,*) '   MFS Bulk formulation'
+         CASE( jp_none    )   ;   WRITE(numout,*) '   OPA coupled to SAS via oasis'
+            IF( ln_mixcpl )       WRITE(numout,*) '       + forced-coupled mixed formulation'
+         END SELECT
          IF( nn_components/= jp_iam_nemo )  &
-            &                       WRITE(numout,*) '              + OASIS coupled SAS'
+            &                     WRITE(numout,*) '       + OASIS coupled SAS'
       ENDIF
       !
       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
-
-!     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)
+      !     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)
          !
          IF(lwp)THEN
@@ -273,5 +267,5 @@
          ENDIF
       ENDIF
-
+      !
       IF( MOD( nitend - nit000 + 1, nn_fsbc) /= 0 .OR.   &
           MOD( nstock             , nn_fsbc) /= 0 ) THEN 
@@ -286,15 +280,15 @@
       IF( ln_dm2dc .AND. ( ( NINT(rday) / ( nn_fsbc * NINT(rdt) ) )  < 8 ) )   &
          &   CALL ctl_warn( 'diurnal cycle for qsr: the sampling of the diurnal cycle is too small...' )
-
-                               CALL sbc_ssm_init               ! Sea-surface mean fields initialisation
-      !
-      IF( ln_ssr           )   CALL sbc_ssr_init               ! Sea-Surface Restoring initialisation
-      !
-                               CALL sbc_rnf_init               ! Runof initialisation
-      !
-      IF( nn_ice == 3      )   CALL sbc_lim_init               ! LIM3 initialisation
-
-      IF( nn_ice == 4      )   CALL cice_sbc_init( nsbc )      ! CICE initialisation
-      
+      !
+                          CALL sbc_ssm_init               ! Sea-surface mean fields initialisation
+      !
+      IF( ln_ssr      )   CALL sbc_ssr_init               ! Sea-Surface Restoring initialisation
+      !
+                          CALL sbc_rnf_init               ! Runof initialisation
+      !
+      IF( nn_ice == 3 )   CALL sbc_lim_init               ! LIM3 initialisation
+      !
+      IF( nn_ice == 4 )   CALL cice_sbc_init( nsbc )      ! CICE initialisation
+      !
    END SUBROUTINE sbc_init
 
@@ -327,8 +321,6 @@
          vtau_b(:,:) = vtau(:,:)                         ! (except at nit000 where before fields
          qns_b (:,:) = qns (:,:)                         !  are set at the end of the routine)
-         ! The 3D heat content due to qsr forcing is treated in traqsr
-         ! qsr_b (:,:) = qsr (:,:)
-         emp_b(:,:) = emp(:,:)
-         sfx_b(:,:) = sfx(:,:)
+         emp_b (:,:) = emp (:,:)
+         sfx_b (:,:) = sfx (:,:)
       ENDIF
       !                                            ! ---------------------------------------- !
@@ -336,8 +328,9 @@
       !                                            ! ---------------------------------------- !
       !
-      IF( nn_components /= jp_iam_sas )   CALL sbc_ssm( kt )   ! ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m)
+      IF( nn_components /= jp_iam_sas )   CALL sbc_ssm ( kt )  ! ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m)
       !                                                        ! averaged over nf_sbc time-step
-
-      IF (ln_wave) CALL sbc_wave( kt )
+      IF( ln_wave                     )   CALL sbc_wave( kt )  ! surface waves
+      
+      
                                                    !==  sbc formulation  ==!
                                                             
@@ -357,13 +350,13 @@
       CASE( jp_mfs   )   ;   CALL sbc_blk_mfs ( kt )                    ! bulk formulation : MFS for the ocean
       CASE( jp_none  ) 
-         IF( nn_components == jp_iam_opa ) &
-                             CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! OPA-SAS coupling: OPA receiving fields from SAS
+         IF( nn_components == jp_iam_opa )   &
+            &                CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! OPA-SAS coupling: OPA receiving fields from SAS
       END SELECT
 
       IF( ln_mixcpl )        CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! forced-coupled mixed formulation after forcing
 
-
+      !
       !                                            !==  Misc. Options  ==!
-      
+      !
       SELECT CASE( nn_ice )                                       ! Update heat and freshwater fluxes over sea-ice areas
       CASE(  1 )   ;         CALL sbc_ice_if   ( kt )                ! Ice-cover climatology ("Ice-if" model)
@@ -375,5 +368,5 @@
       IF( ln_icebergs    )   CALL icb_stp( kt )                   ! compute icebergs
 
-      IF( nn_isf   /= 0  )   CALL sbc_isf( kt )                    ! compute iceshelves
+      IF( ln_isf         )   CALL sbc_isf( kt )                   ! compute iceshelves
 
       IF( ln_rnf         )   CALL sbc_rnf( kt )                   ! add runoffs to fresh water fluxes
@@ -383,6 +376,9 @@
       IF( nn_fwb    /= 0 )   CALL sbc_fwb( kt, nn_fwb, nn_fsbc )  ! control the freshwater budget
 
-      IF( nn_closea == 1 )   CALL sbc_clo( kt )                   ! treatment of closed sea in the model domain 
-      !                                                           ! (update freshwater fluxes)
+      ! treatment of closed sea in the model domain 
+      ! (update freshwater fluxes)
+      ! Should not be ran if ln_diurnal_only
+      IF( .NOT.(ln_diurnal_only) .AND. (nn_closea == 1) )   CALL sbc_clo( kt )   
+
 !RBbug do not understand why see ticket 667
 !clem: it looks like it is necessary for the north fold (in certain circumstances). Don't know why.
@@ -430,5 +426,4 @@
          CALL iom_rstput( kt, nitrst, numrow, 'sfx_b'  , sfx  )
       ENDIF
-
       !                                                ! ---------------------------------------- !
       !                                                !        Outputs and control print         !
@@ -452,7 +447,7 @@
       !
       IF(ln_ctl) THEN         ! print mean trends (used for debugging)
-         CALL prt_ctl(tab2d_1=fr_i              , clinfo1=' fr_i     - : ', mask1=tmask, ovlap=1 )
-         CALL prt_ctl(tab2d_1=(emp-rnf + fwfisf), clinfo1=' emp-rnf  - : ', mask1=tmask, ovlap=1 )
-         CALL prt_ctl(tab2d_1=(sfx-rnf + fwfisf), clinfo1=' sfx-rnf  - : ', mask1=tmask, ovlap=1 )
+         CALL prt_ctl(tab2d_1=fr_i              , clinfo1=' fr_i    - : ', mask1=tmask, ovlap=1 )
+         CALL prt_ctl(tab2d_1=(emp-rnf + fwfisf), clinfo1=' emp-rnf - : ', mask1=tmask, ovlap=1 )
+         CALL prt_ctl(tab2d_1=(sfx-rnf + fwfisf), clinfo1=' sfx-rnf - : ', mask1=tmask, ovlap=1 )
          CALL prt_ctl(tab2d_1=qns              , clinfo1=' qns      - : ', mask1=tmask, ovlap=1 )
          CALL prt_ctl(tab2d_1=qsr              , clinfo1=' qsr      - : ', mask1=tmask, ovlap=1 )
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 6121)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 6140)
@@ -12,50 +12,51 @@
 
    !!----------------------------------------------------------------------
-   !!   sbc_rnf      : monthly runoffs read in a NetCDF file
-   !!   sbc_rnf_init : runoffs initialisation
-   !!   rnf_mouth    : set river mouth mask
+   !!   sbc_rnf       : monthly runoffs read in a NetCDF file
+   !!   sbc_rnf_init  : runoffs initialisation
+   !!   rnf_mouth     : set river mouth mask
    !!----------------------------------------------------------------------
-   USE dom_oce         ! ocean space and time domain
-   USE phycst          ! physical constants
-   USE sbc_oce         ! surface boundary condition variables
-   USE sbcisf          ! PM we could remove it I think
-   USE closea          ! closed seas
-   USE fldread         ! read input field at current time step
-   USE in_out_manager  ! I/O manager
-   USE iom             ! I/O module
-   USE lib_mpp         ! MPP library
-   USE eosbn2
-   USE wrk_nemo        ! Memory allocation
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   USE sbc_oce        ! surface boundary condition variables
+   USE sbcisf         ! PM we could remove it I think
+   USE closea         ! closed seas
+   USE eosbn2         ! Equation Of State
+   !
+   USE in_out_manager ! I/O manager
+   USE fldread        ! read input field at current time step
+   USE iom            ! I/O module
+   USE lib_mpp        ! MPP library
+   USE wrk_nemo       ! Memory allocation
 
    IMPLICIT NONE
    PRIVATE
 
-   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
+   PUBLIC   sbc_rnf       ! called in sbcmod module
+   PUBLIC   sbc_rnf_div   ! called in divhor module
+   PUBLIC   sbc_rnf_alloc ! called in sbcmod module
+   PUBLIC   sbc_rnf_init  ! called in sbcmod module
    
-   !                                                     !!* namsbc_rnf namelist *
-   CHARACTER(len=100)         ::   cn_dir          !: Root directory for location of rnf files
+   !                                                !!* namsbc_rnf namelist *
+   CHARACTER(len=100)         ::   cn_dir            !: Root directory for location of rnf files
    LOGICAL                    ::   ln_rnf_depth      !: depth       river runoffs attribute specified in a file
-   LOGICAL                    ::   ln_rnf_depth_ini  !: depth       river runoffs  computed at the initialisation
-   REAL(wp)                   ::   rn_rnf_max        !: maximum value of the runoff climatologie ( ln_rnf_depth_ini = .true )
-   REAL(wp)                   ::   rn_dep_max        !: depth over which runoffs is spread ( ln_rnf_depth_ini = .true )
-   INTEGER                    ::   nn_rnf_depth_file !: create (=1) a runoff depth file or not (=0)
-   LOGICAL                    ::   ln_rnf_tem      !: temperature river runoffs attribute specified in a file
-   LOGICAL           , PUBLIC ::   ln_rnf_sal      !: salinity    river runoffs attribute specified in a file
-   TYPE(FLD_N)       , PUBLIC ::   sn_rnf          !: information about the runoff file to be read
-   TYPE(FLD_N)                ::   sn_cnf          !: information about the runoff mouth file to be read
-   TYPE(FLD_N)                ::   sn_s_rnf        !: information about the salinities of runoff file to be read
-   TYPE(FLD_N)                ::   sn_t_rnf        !: information about the temperatures of runoff file to be read
-   TYPE(FLD_N)                ::   sn_dep_rnf      !: information about the depth which river inflow affects
-   LOGICAL           , PUBLIC ::   ln_rnf_mouth    !: specific treatment in mouths vicinity
-   REAL(wp)                   ::   rn_hrnf         !: runoffs, depth over which enhanced vertical mixing is used
-   REAL(wp)          , PUBLIC ::   rn_avt_rnf      !: runoffs, value of the additional vertical mixing coef. [m2/s]
-   REAL(wp)                   ::   rn_rfact        !: multiplicative factor for runoff
-
-   LOGICAL           , PUBLIC ::   l_rnfcpl = .false.       ! runoffs recieved from oasis
-
-   INTEGER , PUBLIC  ::   nkrnf = 0         !: nb of levels over which Kz is increased at river mouths
+   LOGICAL                    ::      ln_rnf_depth_ini  !: depth       river runoffs  computed at the initialisation
+   REAL(wp)                   ::      rn_rnf_max        !: maximum value of the runoff climatologie (ln_rnf_depth_ini =T)
+   REAL(wp)                   ::      rn_dep_max        !: depth over which runoffs is spread       (ln_rnf_depth_ini =T)
+   INTEGER                    ::      nn_rnf_depth_file !: create (=1) a runoff depth file or not (=0)
+   LOGICAL                    ::   ln_rnf_tem        !: temperature river runoffs attribute specified in a file
+   LOGICAL           , PUBLIC ::   ln_rnf_sal        !: salinity    river runoffs attribute specified in a file
+   TYPE(FLD_N)       , PUBLIC ::   sn_rnf            !: information about the runoff file to be read
+   TYPE(FLD_N)                ::   sn_cnf            !: information about the runoff mouth file to be read
+   TYPE(FLD_N)                ::   sn_s_rnf          !: information about the salinities of runoff file to be read
+   TYPE(FLD_N)                ::   sn_t_rnf          !: information about the temperatures of runoff file to be read
+   TYPE(FLD_N)                ::   sn_dep_rnf        !: information about the depth which river inflow affects
+   LOGICAL           , PUBLIC ::   ln_rnf_mouth      !: specific treatment in mouths vicinity
+   REAL(wp)                   ::   rn_hrnf           !: runoffs, depth over which enhanced vertical mixing is used
+   REAL(wp)          , PUBLIC ::   rn_avt_rnf        !: runoffs, value of the additional vertical mixing coef. [m2/s]
+   REAL(wp)          , PUBLIC ::   rn_rfact          !: multiplicative factor for runoff
+
+   LOGICAL , PUBLIC ::   l_rnfcpl = .false.   !: runoffs recieved from oasis
+   INTEGER , PUBLIC ::   nkrnf = 0            !: nb of levels over which Kz is increased at river mouths
+   
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rnfmsk              !: river mouth mask (hori.)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)     ::   rnfmsk_z            !: river mouth mask (vert.)
@@ -68,6 +69,4 @@
    TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_t_rnf     ! structure: river runoff temperature (file information, fields read)  
  
-   !! * Substitutions  
-#  include "domzgr_substitute.h90"  
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -127,12 +126,4 @@
       IF(   ln_rnf_sal   )   CALL fld_read ( kt, nn_fsbc, sf_s_rnf )    ! idem for runoffs salinity    if required
       !
-      ! Runoff reduction only associated to the ORCA2_LIM configuration
-      ! when reading the NetCDF file runoff_1m_nomask.nc
-      IF( cp_cfg == 'orca' .AND. jp_cfg == 2 .AND. .NOT. l_rnfcpl )   THEN
-         WHERE( 40._wp < gphit(:,:) .AND. gphit(:,:) < 65._wp )
-            sf_rnf(1)%fnow(:,:,1) = 0.85 * sf_rnf(1)%fnow(:,:,1)
-         END WHERE
-      ENDIF
-      !
       IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
          !
@@ -142,10 +133,10 @@
          IF( ln_rnf_tem ) THEN                                       ! use runoffs temperature data
             rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0
+            CALL eos_fzp( sss_m(:,:), ztfrz(:,:) )
             WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999._wp )             ! if missing data value use SST as runoffs temperature
                rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0
             END WHERE
             WHERE( sf_t_rnf(1)%fnow(:,:,1) == -222._wp )             ! where fwf comes from melting of ice shelves or iceberg
-               ztfrz(:,:) = -1.9 !tfreez( sss_m(:,:) ) !PM to be discuss (trouble if sensitivity study)
-               rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * lfusisf * r1_rau0_rcp
+               rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * rlfusisf * r1_rau0_rcp
             END WHERE
          ELSE                                                        ! use SST as runoffs temperature
@@ -211,10 +202,18 @@
       !
       IF( ln_rnf_depth .OR. ln_rnf_depth_ini ) THEN      !==   runoff distributed over several levels   ==!
-         IF( lk_vvl ) THEN             ! variable volume case
+         IF( ln_linssh ) THEN    !* constant volume case : just apply the runoff input flow
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  DO jk = 1, nk_rnf(ji,jj)
+                     phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * r1_rau0 / h_rnf(ji,jj)
+                  END DO
+               END DO
+            END DO
+         ELSE                    !* variable volume case
             DO jj = 1, jpj                   ! update the depth over which runoffs are distributed
                DO ji = 1, jpi
                   h_rnf(ji,jj) = 0._wp
                   DO jk = 1, nk_rnf(ji,jj)                           ! recalculates h_rnf to be the depth in metres
-                     h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk)   ! to the bottom of the relevant grid box
+                     h_rnf(ji,jj) = h_rnf(ji,jj) + e3t_n(ji,jj,jk)   ! to the bottom of the relevant grid box
                   END DO
                   !                          ! apply the runoff input flow
@@ -224,18 +223,8 @@
                END DO
             END DO
-         ELSE                          ! constant volume case : just apply the runoff input flow
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  DO jk = 1, nk_rnf(ji,jj)
-                     phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * r1_rau0 / h_rnf(ji,jj)
-                  END DO
-               END DO
-            END DO
          ENDIF
       ELSE                       !==   runoff put only at the surface   ==!
-         IF( lk_vvl ) THEN              ! variable volume case
-            h_rnf(:,:) = fse3t(:,:,1)   ! recalculate h_rnf to be depth of top box
-         ENDIF
-         phdivn(:,:,1) = phdivn(:,:,1) - ( rnf(:,:) + rnf_b(:,:) ) * zfact * r1_rau0 / fse3t(:,:,1)
+         h_rnf (:,:)   = e3t_n (:,:,1)        ! update h_rnf to be depth of top box
+         phdivn(:,:,1) = phdivn(:,:,1) - ( rnf(:,:) + rnf_b(:,:) ) * zfact * r1_rau0 / e3t_n(:,:,1)
       ENDIF
       !
@@ -377,5 +366,5 @@
                h_rnf(ji,jj) = 0._wp
                DO jk = 1, nk_rnf(ji,jj)
-                  h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk)
+                  h_rnf(ji,jj) = h_rnf(ji,jj) + e3t_n(ji,jj,jk)
                END DO
             END DO
@@ -435,5 +424,5 @@
                h_rnf(ji,jj) = 0._wp
                DO jk = 1, nk_rnf(ji,jj)
-                  h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk)
+                  h_rnf(ji,jj) = h_rnf(ji,jj) + e3t_n(ji,jj,jk)
                END DO
             END DO
@@ -448,5 +437,5 @@
       ELSE                                       ! runoffs applied at the surface
          nk_rnf(:,:) = 1
-         h_rnf (:,:) = fse3t(:,:,1)
+         h_rnf (:,:) = e3t_n(:,:,1)
       ENDIF
       !
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90	(revision 6121)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90	(revision 6140)
@@ -6,31 +6,29 @@
    !! History :  9.0  ! 2006-07  (G. Madec)  Original code
    !!            3.3  ! 2010-10  (C. Bricaud, G. Madec)  add the Patm forcing for sea-ice
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   sbc_ssm        : calculate sea surface mean currents, temperature,  
-   !!                    and salinity over nn_fsbc time-step
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers
-   USE dom_oce         ! ocean space and time domain
-   USE sbc_oce         ! surface boundary condition: ocean fields
-   USE sbcapr          ! surface boundary condition: atmospheric pressure
-   USE eosbn2          ! equation of state and related derivatives
+   !!            3.7  ! 2015-11  (G. Madec)  non linear free surface by default: e3t_m always computed
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   sbc_ssm       : calculate sea surface mean currents, temperature,  
+   !!                   and salinity over nn_fsbc time-step
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers
+   USE dom_oce        ! ocean space and time domain
+   USE sbc_oce        ! surface boundary condition: ocean fields
+   USE sbcapr         ! surface boundary condition: atmospheric pressure
+   USE eosbn2         ! equation of state and related derivatives
    !
-   USE in_out_manager  ! I/O manager
-   USE prtctl          ! Print control
-   USE iom             ! IOM library
+   USE in_out_manager ! I/O manager
+   USE prtctl         ! Print control
+   USE iom            ! IOM library
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   sbc_ssm         ! routine called by step.F90
-   PUBLIC   sbc_ssm_init    ! routine called by sbcmod.F90
-
-   LOGICAL, SAVE  ::   l_ssm_mean = .FALSE.       ! keep track of whether means have been read
-                                                  ! from restart file
+   PUBLIC   sbc_ssm        ! routine called by step.F90
+   PUBLIC   sbc_ssm_init   ! routine called by sbcmod.F90
+
+   LOGICAL, SAVE ::   l_ssm_mean = .FALSE.   ! keep track of whether means have been read from restart file
    
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -59,5 +57,5 @@
       REAL(wp), DIMENSION(jpi,jpj,jpts) :: zts
       !!---------------------------------------------------------------------
-
+      !
       !                                        !* surface T-, U-, V- ocean level variables (T, S, depth, velocity)
       DO jj = 1, jpj
@@ -81,5 +79,5 @@
          ENDIF
          !
-         IF( lk_vvl )   e3t_m(:,:) = fse3t_n(:,:,1)
+         e3t_m(:,:) = e3t_n(:,:,1)
          !
          frq_m(:,:) = fraqsr_1lev(:,:)
@@ -103,5 +101,5 @@
             ENDIF
             !
-            IF( lk_vvl )   e3t_m(:,:) = zcoef * fse3t_n(:,:,1)
+            e3t_m(:,:) = zcoef * e3t_n(:,:,1)
             !
             frq_m(:,:) = zcoef * fraqsr_1lev(:,:)
@@ -109,11 +107,11 @@
          ELSEIF( MOD( kt - 2 , nn_fsbc ) == 0 ) THEN      !   Initialisation: New mean computation   !
             !                                             ! ---------------------------------------- !
-            ssu_m(:,:) = 0.e0      ! reset to zero ocean mean sbc fields
-            ssv_m(:,:) = 0.e0
-            sst_m(:,:) = 0.e0
-            sss_m(:,:) = 0.e0
-            ssh_m(:,:) = 0.e0
-            IF( lk_vvl )   e3t_m(:,:) = 0.e0
-            frq_m(:,:) = 0.e0
+            ssu_m(:,:) = 0._wp     ! reset to zero ocean mean sbc fields
+            ssv_m(:,:) = 0._wp
+            sst_m(:,:) = 0._wp
+            sss_m(:,:) = 0._wp
+            ssh_m(:,:) = 0._wp
+            e3t_m(:,:) = 0._wp
+            frq_m(:,:) = 0._wp
          ENDIF
          !                                                ! ---------------------------------------- !
@@ -131,7 +129,7 @@
          ENDIF
          !
-         IF( lk_vvl )   e3t_m(:,:) = fse3t_m(:,:) + fse3t_n(:,:,1)
-         !
-         frq_m(:,:) =   frq_m(:,:) + fraqsr_1lev(:,:)
+         e3t_m(:,:) = e3t_m(:,:) + e3t_n(:,:,1)
+         !
+         frq_m(:,:) = frq_m(:,:) + fraqsr_1lev(:,:)
 
          !                                                ! ---------------------------------------- !
@@ -139,11 +137,11 @@
             !                                             ! ---------------------------------------- !
             zcoef = 1. / REAL( nn_fsbc, wp )
-            sst_m(:,:) = sst_m(:,:) * zcoef           ! mean SST             [Celcius]
-            sss_m(:,:) = sss_m(:,:) * zcoef           ! mean SSS             [psu]
-            ssu_m(:,:) = ssu_m(:,:) * zcoef           ! mean suface current  [m/s]
-            ssv_m(:,:) = ssv_m(:,:) * zcoef           !
-            ssh_m(:,:) = ssh_m(:,:) * zcoef           ! mean SSH             [m]
-            IF( lk_vvl )   e3t_m(:,:) = fse3t_m(:,:) * zcoef   ! mean vertical scale factor [m]
-            frq_m(:,:) = frq_m(:,:) * zcoef   ! mean fraction of solar net radiation absorbed in the 1st T level [-]
+            sst_m(:,:) = sst_m(:,:) * zcoef     ! mean SST             [Celcius]
+            sss_m(:,:) = sss_m(:,:) * zcoef     ! mean SSS             [psu]
+            ssu_m(:,:) = ssu_m(:,:) * zcoef     ! mean suface current  [m/s]
+            ssv_m(:,:) = ssv_m(:,:) * zcoef     !
+            ssh_m(:,:) = ssh_m(:,:) * zcoef     ! mean SSH             [m]
+            e3t_m(:,:) = e3t_m(:,:) * zcoef     ! mean vertical scale factor [m]
+            frq_m(:,:) = frq_m(:,:) * zcoef     ! mean fraction of solar net radiation absorbed in the 1st T level [-]
             !
          ENDIF
@@ -162,5 +160,5 @@
             CALL iom_rstput( kt, nitrst, numrow, 'sss_m'  , sss_m  )
             CALL iom_rstput( kt, nitrst, numrow, 'ssh_m'  , ssh_m  )
-            IF( lk_vvl )   CALL iom_rstput( kt, nitrst, numrow, 'e3t_m'  , e3t_m  )
+            CALL iom_rstput( kt, nitrst, numrow, 'e3t_m'  , e3t_m  )
             CALL iom_rstput( kt, nitrst, numrow, 'frq_m'  , frq_m  )
             !
@@ -175,9 +173,10 @@
          CALL iom_put( 'sss_m', sss_m )
          CALL iom_put( 'ssh_m', ssh_m )
-         IF( lk_vvl )   CALL iom_put( 'e3t_m', e3t_m )
+         CALL iom_put( 'e3t_m', e3t_m )
          CALL iom_put( 'frq_m', frq_m )
       ENDIF
       !
    END SUBROUTINE sbc_ssm
+
 
    SUBROUTINE sbc_ssm_init
@@ -189,7 +188,7 @@
       !! ** Action  : - read parameters
       !!----------------------------------------------------------------------
-      REAL(wp) ::   zcoef, zf_sbc       ! local scalar
+      REAL(wp) ::   zcoef, zf_sbc   ! local scalar
       !!----------------------------------------------------------------------
-
+      !
       IF( nn_fsbc == 1 ) THEN
          !
@@ -206,11 +205,11 @@
          IF( ln_rstart .AND. iom_varid( numror, 'nn_fsbc', ldstop = .FALSE. ) > 0 ) THEN
             l_ssm_mean = .TRUE.
-            CALL iom_get( numror               , 'nn_fsbc', zf_sbc )   ! sbc frequency of previous run
-            CALL iom_get( numror, jpdom_autoglo, 'ssu_m'  , ssu_m  )   ! sea surface mean velocity    (T-point)
-            CALL iom_get( numror, jpdom_autoglo, 'ssv_m'  , ssv_m  )   !   "         "    velocity    (V-point)
-            CALL iom_get( numror, jpdom_autoglo, 'sst_m'  , sst_m  )   !   "         "    temperature (T-point)
-            CALL iom_get( numror, jpdom_autoglo, 'sss_m'  , sss_m  )   !   "         "    salinity    (T-point)
-            CALL iom_get( numror, jpdom_autoglo, 'ssh_m'  , ssh_m  )   !   "         "    height      (T-point)
-            IF( lk_vvl ) CALL iom_get( numror, jpdom_autoglo, 'e3t_m', e3t_m )
+            CALL iom_get( numror               , 'nn_fsbc', zf_sbc )    ! sbc frequency of previous run
+            CALL iom_get( numror, jpdom_autoglo, 'ssu_m'  , ssu_m  )    ! sea surface mean velocity    (U-point)
+            CALL iom_get( numror, jpdom_autoglo, 'ssv_m'  , ssv_m  )    !   "         "    velocity    (V-point)
+            CALL iom_get( numror, jpdom_autoglo, 'sst_m'  , sst_m  )    !   "         "    temperature (T-point)
+            CALL iom_get( numror, jpdom_autoglo, 'sss_m'  , sss_m  )    !   "         "    salinity    (T-point)
+            CALL iom_get( numror, jpdom_autoglo, 'ssh_m'  , ssh_m  )    !   "         "    height      (T-point)
+            CALL iom_get( numror, jpdom_autoglo, 'e3t_m'  , e3t_m  )    ! 1st level thickness          (T-point)
             ! fraction of solar net radiation absorbed in 1st T level
             IF( iom_varid( numror, 'frq_m', ldstop = .FALSE. ) > 0 ) THEN
@@ -229,5 +228,5 @@
                sss_m(:,:) = zcoef * sss_m(:,:)
                ssh_m(:,:) = zcoef * ssh_m(:,:)
-               IF( lk_vvl )   e3t_m(:,:) = zcoef * fse3t_m(:,:)
+               e3t_m(:,:) = zcoef * e3t_m(:,:)
                frq_m(:,:) = zcoef * frq_m(:,:)
             ELSE
@@ -245,7 +244,7 @@
          ELSE                    ;   sst_m(:,:) = tsn(:,:,1,jp_tem)
          ENDIF
-         sss_m(:,:) = tsn(:,:,1,jp_sal)
-         ssh_m(:,:) = sshn(:,:)
-         IF( lk_vvl )   e3t_m(:,:) = fse3t_n(:,:,1)
+         sss_m(:,:) = tsn  (:,:,1,jp_sal)
+         ssh_m(:,:) = sshn (:,:)
+         e3t_m(:,:) = e3t_n(:,:,1)
          frq_m(:,:) = 1._wp
          !
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90	(revision 6121)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90	(revision 6140)
@@ -19,6 +19,6 @@
    !
    USE fldread        ! read input fields
+   USE in_out_manager ! I/O manager
    USE iom            ! I/O manager
-   USE in_out_manager ! I/O manager
    USE lib_mpp        ! distribued memory computing library
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
@@ -47,6 +47,4 @@
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_sss   ! structure of input SSS (file informations, fields read)
 
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO Consortium (2011)
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbctide.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbctide.F90	(revision 6121)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbctide.F90	(revision 6140)
@@ -6,14 +6,14 @@
    !! History :  9.0  !  2007  (O. Le Galloudec)  Original code
    !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers variables
-   USE dom_oce         ! ocean space and time domain
-   USE phycst
-   USE daymod
-   USE tideini
+   USE oce            ! ocean dynamics and tracers variables
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constant
+   USE daymod         ! calandar
+   USE tideini        ! 
    !
-   USE iom
-   USE in_out_manager  ! I/O units
-   USE ioipsl          ! NetCDF IPSL library
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE in_out_manager ! I/O units
+   USE iom            ! xIOs server
+   USE ioipsl         ! NetCDF IPSL library
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
 
    IMPLICIT NONE
@@ -46,7 +46,8 @@
       INTEGER, INTENT( in ) ::   kt     ! ocean time-step
       INTEGER               ::   jk     ! dummy loop index
+      INTEGER               ::   nsec_day_orig     ! Temporary variable
       !!----------------------------------------------------------------------
-
-      IF( nsec_day == NINT(0.5_wp * rdttra(1)) ) THEN      ! start a new day
+      
+      IF( nsec_day == NINT(0.5_wp * rdt) .OR. kt == nit000 ) THEN      ! start a new day
          !
          IF( kt == nit000 ) THEN
@@ -59,7 +60,16 @@
          pot_astro(:,:) = 0._wp
          !
+         ! If the run does not start from midnight then need to initialise tides
+         ! at the start of the current day (only occurs when kt==nit000)
+         ! Temporarily set nsec_day to beginning of day.
+         nsec_day_orig = nsec_day
+         IF ( nsec_day /= NINT(0.5_wp * rdt) ) THEN 
+            kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt
+            nsec_day = NINT(0.5_wp * rdt)
+         ELSE
+            kt_tide = kt 
+         ENDIF
          CALL tide_harmo( omega_tide, v0tide, utide, ftide, ntide, nb_harmo )
          !
-         kt_tide = kt
          !
          IF(lwp) THEN
@@ -74,4 +84,6 @@
          IF( ln_tide_pot )   CALL tide_init_potential
          !
+         ! Reset nsec_day
+         nsec_day = nsec_day_orig 
       ENDIF
       !
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 6121)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 6140)
@@ -40,5 +40,4 @@
 
    !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -72,5 +71,5 @@
       REAL(wp), DIMENSION(:,:,:), POINTER ::   zusd_t, zvsd_t, ze3hdiv   ! 3D workspace
       !!
-      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_wn
+      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_wn, ln_cdgw , ln_sdw
       !!---------------------------------------------------------------------
       !
@@ -81,5 +80,5 @@
          READ  ( numnam_ref, namsbc_wave, IOSTAT = ios, ERR = 901)
 901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in reference namelist', lwp )
-
+         !
          REWIND( numnam_cfg )              ! Namelist namsbc_wave in configuration namelist : File for drag coeff. from wave model
          READ  ( numnam_cfg, namsbc_wave, IOSTAT = ios, ERR = 902 )
@@ -87,5 +86,16 @@
          IF(lwm) WRITE ( numond, namsbc_wave )
          !
-         IF ( ln_cdgw ) THEN
+         IF(lwp) THEN               ! Control print
+            WRITE(numout,*) '        Namelist namsbc_wave : surface wave setting' 
+            WRITE(numout,*) '           wave drag coefficient                      ln_cdgw  = ', ln_cdgw  
+            WRITE(numout,*) '           wave stokes drift                          ln_sdw   = ', ln_sdw
+         ENDIF
+         !
+         IF( .NOT.( ln_cdgw .OR. ln_sdw ) )    &
+            &  CALL ctl_warn( 'ln_sbcwave=T but nor drag coefficient (ln_cdgw=F) neither stokes drift activated (ln_sdw=F)' )
+         IF( ln_cdgw .AND. .NOT.(ln_blk_mfs .OR. ln_blk_core) )   &       
+            &  CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core')
+         !
+         IF( ln_cdgw ) THEN
             ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg
             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' )
@@ -97,5 +107,5 @@
             cdn_wave(:,:) = 0.0
          ENDIF
-         IF ( ln_sdw ) THEN
+         IF( ln_sdw ) THEN
             slf_i(jp_usd) = sn_usd ; slf_i(jp_vsd) = sn_vsd; slf_i(jp_wn) = sn_wn
             ALLOCATE( sf_sd(3), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg
@@ -115,22 +125,20 @@
       ENDIF
       !
-      IF ( ln_cdgw ) THEN              !==  Neutral drag coefficient  ==!
+      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  ==!
+      IF( ln_sdw )  THEN               !==  Computation of the 3d Stokes Drift  ==!
+         !
+         CALL wrk_alloc( jpi,jpj,jpk,   zusd_t, zvsd_t, ze3hdiv )
          !
          CALL fld_read( kt, nn_fsbc, sf_sd )    !* read drag coefficient from external forcing
          !
-         !
-         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) )
+         DO jk = 1, jpkm1                       !* distribute it on the vertical
+            zusd_t(:,:,jk) = sf_sd(jp_usd)%fnow(:,:,1) * EXP( -2._wp * sf_sd(jp_wn)%fnow(:,:,1) * gdept_n(:,:,jk) )
+            zvsd_t(:,:,jk) = sf_sd(jp_vsd)%fnow(:,:,1) * EXP( -2._wp * sf_sd(jp_wn)%fnow(:,:,1) * gdept_n(:,:,jk) )
          END DO
-         !                                      !* interpolate the stokes drift from t-point to u- and v-points
-         DO jk = 1, jpkm1
+         DO jk = 1, jpkm1                       !* interpolate the stokes drift from t-point to u- and v-points
             DO jj = 1, jpjm1
                DO ji = 1, jpim1
@@ -146,8 +154,8 @@
             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)
+                  ze3hdiv(ji,jj,jk) = (  e2u(ji  ,jj) * e3u_n(ji  ,jj,jk) * usd3d(ji  ,jj,jk)     &
+                     &                 - e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) * usd3d(ji-1,jj,jk)     &
+                     &                 + e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * vsd3d(ji,jj  ,jk)     &
+                     &                 - e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vsd3d(ji,jj-1,jk)   ) * r1_e1e2t(ji,jj)
                END DO  
             END DO  
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/tideini.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/tideini.F90	(revision 6121)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/tideini.F90	(revision 6140)
@@ -6,14 +6,14 @@
    !! History :  1.0  !  2007  (O. Le Galloudec)  Original code
    !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers variables
-   USE dom_oce         ! ocean space and time domain
-   USE phycst
-   USE daymod
-   USE tide_mod
+   USE oce            ! ocean dynamics and tracers variables
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constant
+   USE daymod         ! calandar
+   USE tide_mod       ! 
    !
-   USE iom
-   USE in_out_manager  ! I/O units
-   USE ioipsl          ! NetCDF IPSL library
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE in_out_manager ! I/O units
+   USE iom            ! xIOs server
+   USE ioipsl         ! NetCDF IPSL library
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
 
    IMPLICIT NONE
@@ -27,7 +27,7 @@
    LOGICAL , PUBLIC ::   ln_tide_pot     !:
    LOGICAL , PUBLIC ::   ln_tide_ramp    !:
-   INTEGER , PUBLIC ::   nb_harmo                 !:
-   INTEGER , PUBLIC ::   kt_tide                  !:
-   REAL(wp), PUBLIC ::   rdttideramp              !:
+   INTEGER , PUBLIC ::   nb_harmo        !:
+   INTEGER , PUBLIC ::   kt_tide         !:
+   REAL(wp), PUBLIC ::   rdttideramp     !:
    
    INTEGER , PUBLIC, ALLOCATABLE, DIMENSION(:) ::   ntide   !:
@@ -40,75 +40,68 @@
 CONTAINS
    
-  SUBROUTINE tide_init ( kt )
-    !!----------------------------------------------------------------------
-    !!                 ***  ROUTINE tide_init  ***
-    !!----------------------------------------------------------------------      
-    !! * Local declarations
-    INTEGER  :: ji, jk
-    INTEGER, INTENT( in ) ::   kt     ! ocean time-step
-    CHARACTER(LEN=4), DIMENSION(jpmax_harmo) :: clname
-    INTEGER  ::   ios                 ! Local integer output status for namelist read
-    !
-    NAMELIST/nam_tide/ln_tide_pot, ln_tide_ramp, rdttideramp, clname
-    !!----------------------------------------------------------------------
-
-    IF ( kt == nit000 ) THEN
-       !
-       IF(lwp) THEN
-          WRITE(numout,*)
-          WRITE(numout,*) 'tide_init : Initialization of the tidal components'
-          WRITE(numout,*) '~~~~~~~~~ '
-       ENDIF
-       !
-       CALL tide_init_Wave
-       !
-       ! Read Namelist nam_tide
-       REWIND( numnam_ref )              ! Namelist nam_tide in reference namelist : Tides
-       READ  ( numnam_ref, nam_tide, IOSTAT = ios, ERR = 901)
-901    IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_tide in reference namelist', lwp )
-
-       REWIND( numnam_cfg )              ! Namelist nam_tide in configuration namelist : Tides
-       READ  ( numnam_cfg, nam_tide, IOSTAT = ios, ERR = 902 )
-902    IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_tide in configuration namelist', lwp )
-       IF(lwm) WRITE ( numond, nam_tide )
-       !
-       nb_harmo=0
-       DO jk = 1, jpmax_harmo
-          DO ji = 1,jpmax_harmo
-             IF( TRIM(clname(jk)) == Wave(ji)%cname_tide )   nb_harmo = nb_harmo + 1
-          END DO
-       END DO
-       !       
-       ! Ensure that tidal components have been set in namelist_cfg
-       IF( nb_harmo .EQ. 0 ) CALL ctl_stop( 'tide_init : No tidal components set in nam_tide' )
-       !
-       IF(lwp) THEN
-          WRITE(numout,*) '   Namelist nam_tide'
-          WRITE(numout,*) '      Apply astronomical potential : ln_tide_pot  =', ln_tide_pot
-          WRITE(numout,*) '                                     nb_harmo     = ', nb_harmo
-          WRITE(numout,*) '                                     ln_tide_ramp = ', ln_tide_ramp 
-          WRITE(numout,*) '                                     rdttideramp  = ', rdttideramp
-       ENDIF
-       IF( ln_tide_ramp.AND.((nitend-nit000+1)*rdt/rday < rdttideramp) )   &
-          &   CALL ctl_stop('rdttideramp must be lower than run duration')
-       IF( ln_tide_ramp.AND.(rdttideramp<0.) ) &
-          &   CALL ctl_stop('rdttideramp must be positive')
-       !
-       !
-       ALLOCATE( ntide(nb_harmo) )
-       DO jk = 1, nb_harmo
-          DO ji = 1, jpmax_harmo
-             IF( TRIM(clname(jk)) .eq. Wave(ji)%cname_tide ) THEN
-                ntide(jk) = ji
-                EXIT
-             END IF
-          END DO
-       END DO
-       !
-       ALLOCATE( omega_tide(nb_harmo), v0tide    (nb_harmo),   &
-          &      utide     (nb_harmo), ftide     (nb_harmo)  )
-       kt_tide = kt
-       !
+   SUBROUTINE tide_init
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE tide_init  ***
+      !!----------------------------------------------------------------------      
+      INTEGER  :: ji, jk
+      CHARACTER(LEN=4), DIMENSION(jpmax_harmo) :: clname
+      INTEGER  ::   ios                 ! Local integer output status for namelist read
+      !
+      NAMELIST/nam_tide/ln_tide_pot, ln_tide_ramp, rdttideramp, clname
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'tide_init : Initialization of the tidal components'
+         WRITE(numout,*) '~~~~~~~~~ '
       ENDIF
+      !
+      CALL tide_init_Wave
+      !
+      ! Read Namelist nam_tide
+      REWIND( numnam_ref )              ! Namelist nam_tide in reference namelist : Tides
+      READ  ( numnam_ref, nam_tide, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_tide in reference namelist', lwp )
+      !
+      REWIND( numnam_cfg )              ! Namelist nam_tide in configuration namelist : Tides
+      READ  ( numnam_cfg, nam_tide, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_tide in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, nam_tide )
+      !
+      nb_harmo=0
+      DO jk = 1, jpmax_harmo
+         DO ji = 1,jpmax_harmo
+            IF( TRIM(clname(jk)) == Wave(ji)%cname_tide )   nb_harmo = nb_harmo + 1
+         END DO
+      END DO
+      !       
+      ! Ensure that tidal components have been set in namelist_cfg
+      IF( nb_harmo == 0 )   CALL ctl_stop( 'tide_init : No tidal components set in nam_tide' )
+      !
+      IF(lwp) THEN
+         WRITE(numout,*) '   Namelist nam_tide'
+         WRITE(numout,*) '      Apply astronomical potential : ln_tide_pot  =', ln_tide_pot
+         WRITE(numout,*) '                                     nb_harmo     = ', nb_harmo
+         WRITE(numout,*) '                                     ln_tide_ramp = ', ln_tide_ramp 
+         WRITE(numout,*) '                                     rdttideramp  = ', rdttideramp
+      ENDIF
+      IF( ln_tide_ramp.AND.((nitend-nit000+1)*rdt/rday < rdttideramp) )   &
+         &   CALL ctl_stop('rdttideramp must be lower than run duration')
+      IF( ln_tide_ramp.AND.(rdttideramp<0.) ) &
+         &   CALL ctl_stop('rdttideramp must be positive')
+      !
+      ALLOCATE( ntide(nb_harmo) )
+      DO jk = 1, nb_harmo
+         DO ji = 1, jpmax_harmo
+            IF( TRIM(clname(jk)) == Wave(ji)%cname_tide ) THEN
+               ntide(jk) = ji
+               EXIT
+            ENDIF
+         END DO
+      END DO
+      !
+      ALLOCATE( omega_tide(nb_harmo), v0tide    (nb_harmo),   &
+         &      utide     (nb_harmo), ftide     (nb_harmo)  )
+      kt_tide = nit000
       !
    END SUBROUTINE tide_init
