Index: branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
===================================================================
--- branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 5947)
+++ branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 5948)
@@ -18,19 +18,18 @@
    USE eosbn2          ! equation of state
    USE sbc_oce         ! surface boundary condition: ocean fields
+   USE zdfbfr          !
+   !
+   USE in_out_manager  ! I/O manager
+   USE iom             ! I/O manager library
+   USE fldread         ! read input field at current time step
    USE lbclnk          !
-   USE iom             ! I/O manager library
-   USE in_out_manager  ! I/O manager
    USE wrk_nemo        ! Memory allocation
    USE timing          ! Timing
    USE lib_fortran     ! glob_sum
-   USE zdfbfr
-   USE fldread         ! read input field at current time step
-
-
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   sbc_isf, sbc_isf_div, sbc_isf_alloc  ! routine called in sbcmod and divcur
+   PUBLIC   sbc_isf, sbc_isf_div, sbc_isf_alloc  ! routine called in sbcmod and divhor
 
    ! public in order to be able to output then 
@@ -53,12 +52,5 @@
    REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  risfLeff               !:effective length (Leff) BG03 nn_isf==2
    REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  ttbl, stbl, utbl, vtbl !:top boundary layer variable at T point
-#if defined key_agrif
-   ! AGRIF can not handle these arrays as integers. The reason is a mystery but problems avoided by declaring them as reals
-   REAL(wp),    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  misfkt, misfkb         !:Level of ice shelf base
-                                                                                          !: (first wet level and last level include in the tbl)
-#else
    INTEGER,    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  misfkt, misfkb         !:Level of ice shelf base
-#endif
-
 
    REAL(wp), PUBLIC, SAVE ::   rcpi   = 2000.0_wp     ! phycst ?
@@ -79,26 +71,29 @@
 #  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)
+   !! NEMO/OPA 3.7 , LOCEAN-IPSL (2015)
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
  
-  SUBROUTINE sbc_isf(kt)
-    INTEGER, INTENT(in)          ::   kt         ! ocean time step
-    INTEGER                      ::   ji, jj, jk, ijkmin, inum, ierror
-    INTEGER                      ::   ikt, ikb   ! top and bottom level of the isf boundary layer
-    REAL(wp)                     ::   rmin
-    REAL(wp)                     ::   zhk
-    CHARACTER(len=256)           ::   cfisf, cvarzisf, cvarhisf   ! name for isf file
-    CHARACTER(LEN=256)           :: cnameis                     ! name of iceshelf file
-    CHARACTER (LEN=32)           :: cvarLeff                    ! variable name for efficient Length scale
-    INTEGER           ::   ios           ! Local integer output status for namelist read
-      !
+   SUBROUTINE sbc_isf(kt)
       !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE sbc_isf  ***
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in)          ::   kt         ! ocean time step
+      !
+      INTEGER                      ::   ji, jj, jk, ijkmin, inum, ierror
+      INTEGER                      ::   ikt, ikb   ! top and bottom level of the isf boundary layer
+      REAL(wp)                     ::   rmin
+      REAL(wp)                     ::   zhk
+      REAL(wp)                     ::   zt_frz, zpress
+      CHARACTER(len=256)           ::   cfisf , cvarzisf, cvarhisf   ! name for isf file
+      CHARACTER(LEN=256)           :: cnameis                     ! name of iceshelf file
+      CHARACTER (LEN=32)           :: cvarLeff                    ! variable name for efficient Length scale
+      INTEGER           ::   ios           ! Local integer output status for namelist read
+      !!
       NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, ln_divisf, ln_conserve, rn_gammat0, rn_gammas0, nn_gammablk, &
-                         & sn_fwfisf, sn_qisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf
-      !
+         &                 sn_fwfisf, sn_qisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf
+      !!---------------------------------------------------------------------
       !
       !                                         ! ====================== !
@@ -113,5 +108,4 @@
 902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp )
          IF(lwm) WRITE ( numond, namsbc_isf )
-
 
          IF ( lwp ) WRITE(numout,*)
@@ -194,7 +188,6 @@
          END IF
          
+         ! compute bottom level of isf tbl and thickness of tbl below the ice shelf
          rhisf_tbl_0(:,:) = rhisf_tbl(:,:)
-
-         ! compute bottom level of isf tbl and thickness of tbl below the ice shelf
          DO jj = 1,jpj
             DO ji = 1,jpi
@@ -217,5 +210,5 @@
             END DO
          END DO
-         
+         !
       END IF
 
@@ -270,10 +263,19 @@
          END IF
          ! compute tsc due to isf
-         ! WARNING water add at temp = 0C, correction term is added in trasbc, maybe better here but need a 3D variable).
-         risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp !
+         ! WARNING water add at temp = 0C, correction term is added, maybe better here but need a 3D variable).
+!         zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04
+         zt_frz = -1.9 !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress )
+         risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - rdivisf * fwfisf(:,:) * zt_frz * r1_rau0 !
          
          ! salt effect already take into account in vertical advection
          risf_tsc(:,:,jp_sal) = (1.0_wp-rdivisf) * fwfisf(:,:) * stbl(:,:) * r1_rau0
-          
+
+         ! output
+         IF( iom_use('qisf'  ) )   CALL iom_put('qisf'  , qisf)
+         IF( iom_use('fwfisf') )   CALL iom_put('fwfisf', fwfisf * stbl(:,:) / soce )
+
+         ! if apply only on the trend and not as a volume flux (rdivisf = 0), fwfisf have to be set to 0 now
+         fwfisf(:,:) = rdivisf * fwfisf(:,:)         
+ 
          ! lbclnk
          CALL lbc_lnk(risf_tsc(:,:,jp_tem),'T',1.)
@@ -295,10 +297,8 @@
          ENDIF
          ! 
-         ! output
-         CALL iom_put('qisf'  , qisf)
-         IF( iom_use('fwfisf') )   CALL iom_put('fwfisf', fwfisf * stbl(:,:) / soce )
       END IF
-  
+      !  
   END SUBROUTINE sbc_isf
+
 
   INTEGER FUNCTION sbc_isf_alloc()
@@ -321,26 +321,26 @@
   END FUNCTION
 
-  SUBROUTINE sbc_isf_bg03(kt)
-   !!==========================================================================
-   !!                 *** SUBROUTINE sbcisf_bg03  ***
-   !! add net heat and fresh water flux from ice shelf melting
-   !! into the adjacent ocean using the parameterisation by
-   !! Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean
-   !!     interaction for climate models", Ocean Modelling 5(2003) 157-170.
-   !!  (hereafter BG)
-   !!==========================================================================
-   !!----------------------------------------------------------------------
-   !!   sbc_isf_bg03      : routine called from sbcmod
-   !!----------------------------------------------------------------------
-   !!
-   !! ** Purpose   :   Add heat and fresh water fluxes due to ice shelf melting
-   !! ** Reference :   Beckmann et Goosse, 2003, Ocean Modelling
-   !!
-   !! History :
-   !!      !  06-02  (C. Wang) Original code
-   !!----------------------------------------------------------------------
-
-    INTEGER, INTENT ( in ) :: kt
-
+
+   SUBROUTINE sbc_isf_bg03(kt)
+      !!==========================================================================
+      !!                 *** SUBROUTINE sbcisf_bg03  ***
+      !! add net heat and fresh water flux from ice shelf melting
+      !! into the adjacent ocean using the parameterisation by
+      !! Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean
+      !!     interaction for climate models", Ocean Modelling 5(2003) 157-170.
+      !!  (hereafter BG)
+      !!==========================================================================
+      !!----------------------------------------------------------------------
+      !!   sbc_isf_bg03      : routine called from sbcmod
+      !!----------------------------------------------------------------------
+      !!
+      !! ** Purpose   :   Add heat and fresh water fluxes due to ice shelf melting
+      !! ** Reference :   Beckmann et Goosse, 2003, Ocean Modelling
+      !!
+      !! History :
+      !!      !  06-02  (C. Wang) Original code
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT ( in ) :: kt
+      !
     INTEGER :: ji, jj, jk, jish  !temporary integer
     INTEGER :: ijkmin
@@ -370,5 +370,5 @@
              ! Calculate freezing temperature
                 zpress = grav*rau0*fsdept(ji,jj,ik)*1.e-04 
-                zt_frz = eos_fzp(tsb(ji,jj,ik,jp_sal), zpress) 
+                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
@@ -386,9 +386,11 @@
              qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp
           END IF
-       ENDDO
-    ENDDO
+       END DO
+    END DO
     !
     IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_bg03')
+      !
   END SUBROUTINE sbc_isf_bg03
+
 
    SUBROUTINE sbc_isf_cav( kt )
@@ -439,5 +441,4 @@
       !
       !
-!CDIR COLLAPSE
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -452,5 +453,5 @@
       zti(:,:)=tinsitu( ttbl, stbl, zpress )
 ! Calculate freezing temperature
-      zfrz(:,:)=eos_fzp( sss_m(:,:), zpress )
+      CALL eos_fzp( sss_m(:,:), zfrz(:,:), zpress )
 
       
@@ -472,9 +473,6 @@
 
                      nit = nit + 1
-                     IF (nit .GE. 100) THEN
-                        !WRITE(numout,*) "sbcisf : too many iteration ... ", zhtflx, zhtflx_b,zgammat, rn_gammat0, rn_tfri2, nn_gammablk, ji,jj
-                        !WRITE(numout,*) "sbcisf : too many iteration ... ", (zhtflx - zhtflx_b)/zhtflx
-                        CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
-                     END IF
+                     IF (nit .GE. 100) CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
+
 ! save gammat and compute zhtflx_b
                      zgammat2d(ji,jj)=zgammat
@@ -496,5 +494,4 @@
 
 ! More complicated 3 equation thermodynamics as in MITgcm
-!CDIR COLLAPSE
          DO jj = 2, jpj
             DO ji = 2, jpi
@@ -565,6 +562,7 @@
       !
       IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_cav')
-
+      !
    END SUBROUTINE sbc_isf_cav
+
 
    SUBROUTINE sbc_isf_gammats(gt, gs, zqhisf, zqwisf, ji, jj, lit )
@@ -693,6 +691,7 @@
                END IF
       END IF
-
+      !
    END SUBROUTINE
+
 
    SUBROUTINE sbc_isf_tbl( varin, varout, cptin )
@@ -756,5 +755,5 @@
       IF (cptin == 'T') CALL lbc_lnk(varout,'T',1.)
       IF (cptin == 'U' .OR. cptin == 'V') CALL lbc_lnk(varout,'T',-1.)
-
+      !
    END SUBROUTINE sbc_isf_tbl
       
@@ -794,5 +793,5 @@
                ! test on tmask useless ?????
                DO jk = ikt, mbkt(ji,jj)
-!                  IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
+                  IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
                END DO
                rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
@@ -823,5 +822,6 @@
       !
    END SUBROUTINE sbc_isf_div
-                        
+
+
    FUNCTION tinsitu( ptem, psal, ppress ) RESULT( pti )
       !!----------------------------------------------------------------------
@@ -874,5 +874,6 @@
       !
    END FUNCTION tinsitu
-   !
+
+
    FUNCTION fsatg( pfps, pfpt, pfphp )
       !!----------------------------------------------------------------------
