MODULE isfparmlt
   !!======================================================================
   !!                       ***  MODULE  sbcisf  ***
   !! Surface module :  update surface ocean boundary condition under ice
   !!                   shelf
   !!======================================================================
   !! History :  4.0  ! original code
   !!----------------------------------------------------------------------

   USE oce            ! ocean dynamics and tracers
   USE isf
   USE isftbl
   USE dom_oce        ! ocean space and time domain
   USE phycst         ! physical constants
   USE eosbn2         ! equation of state

   USE in_out_manager ! I/O manager
   USE iom            ! I/O library
   USE fldread

   IMPLICIT NONE

   PRIVATE

   PUBLIC  isfpar_mlt 
   
   !!----------------------------------------------------------------------
   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
   !! $Id: sbcisf.F90 10536 2019-01-16 19:21:09Z mathiot $
   !! Software governed by the CeCILL license (see ./LICENSE)
   !!----------------------------------------------------------------------
CONTAINS

! -------------------------------------------------------------------------------------------------------
! -------------------------------- PUBLIC SUBROUTINE ----------------------------------------------------
! -------------------------------------------------------------------------------------------------------

  SUBROUTINE isfpar_mlt( kt, pqfwf, pqoce, pqhc )
      !!---------------------------------------------------------------------
      !!                  ***  ROUTINE sbc_isf  ***
      !!
      !! ** Purpose : Compute Salt and Heat fluxes related to ice_shelf 
      !!              melting and freezing 
      !!
      !! ** Method  :  2 parameterizations are available according to XXXXX
      !!                        2 : Beckmann & Goose parameterization
      !!                        3 : Specified runoff in deptht (Mathiot & al. 2017)
      !!----------------------------------------------------------------------
      !!-------------------------- OUT -------------------------------------
      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pqfwf,pqoce, pqhc  ! fresh water, ice-ocean heat and heat content fluxes
      !!-------------------------- IN  -------------------------------------
      INTEGER, INTENT(in) ::   kt   ! ocean time step
      !!---------------------------------------------------------------------
      !
      ! Choose among the available ice shelf parametrisation
      SELECT CASE ( cn_isfpar_mlt )
      CASE ( 'spe' )    ! specified runoff in depth (Mathiot et al., 2017 in preparation)
         CALL isfpar_mlt_spe(kt, pqhc, pqoce, pqfwf)
      CASE ( 'bg03' )    ! Beckmann and Goosse parametrisation 
         CALL isfpar_mlt_bg03(kt, pqhc, pqoce, pqfwf)
      CASE ( 'oasis' )
         !CALL isfpar_mlt_oasis
      CASE DEFAULT
         CALL ctl_stop('STOP', 'unknown isf melt formulation : cn_isfpar (should not see this)')
      END SELECT
      !
   END SUBROUTINE isfpar_mlt

! -------------------------------------------------------------------------------------------------------
! -------------------------------- PRIVATE SUBROUTINE ---------------------------------------------------
! -------------------------------------------------------------------------------------------------------

   SUBROUTINE isfpar_mlt_spe(kt, pqhc, pqfwf, pqoce)
      !!---------------------------------------------------------------------
      !!                  ***  ROUTINE sbc_isf_bg03  ***
      !!
      !! ** Purpose : prescribed ice shelf melting in case ice shelf cavities are closed.
      !!              data read into a forcing files.
      !!
      !! ** Reference : Mathiot et al. (2017)
      !!----------------------------------------------------------------------
      !!-------------------------- OUT -------------------------------------
      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pqhc, pqfwf,pqoce  ! fresh water and ice-ocean heat fluxes
      !!-------------------------- IN  -------------------------------------
      INTEGER,  INTENT(in) :: kt
      !!--------------------------------------------------------------------
      INTEGER :: jk
      REAL(wp), DIMENSION(jpi,jpj,jpk)  :: ztfrz3d
      REAL(wp), DIMENSION(jpi,jpj)      :: ztfrz
      !!--------------------------------------------------------------------
      !
      ! specified runoff in depth (Mathiot et al., 2017)
      CALL fld_read ( kt, nn_fsbc, sf_isfpar_fwf   )
      !
      ! compute ptfrz
      ! 0. ------------Mean freezing point
      DO jk = 1,jpk
         CALL eos_fzp(tsn(:,:,jk,jp_sal), ztfrz3d(:,:,jk), gdept_n(:,:,jk))
      END DO
      CALL isf_tbl(ztfrz3d, ztfrz, 'T', misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par )
      !
      pqfwf(:,:) = - sf_isfpar_fwf(1)%fnow(:,:,1)      ! fresh water flux from the isf (fwfisf <0 mean melting) 
      pqoce(:,:) =   pqfwf(:,:) * rLfusisf             ! ocean/ice shelf flux assume to be equal to latent heat flux
      pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp     ! heat content flux 
      !
   END SUBROUTINE isfpar_mlt_spe

   SUBROUTINE isfpar_mlt_bg03(kt, pqhc, pqoce, pqfwf)
      !!---------------------------------------------------------------------
      !!                  ***  ROUTINE sbc_isf_bg03  ***
      !!
      !! ** Purpose : compute an estimate of ice shelf melting in case cavities are closed
      !!              based on the far fields T and S properties. 
      !!
      !! ** 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)
      !!----------------------------------------------------------------------
      !!-------------------------- OUT -------------------------------------
      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pqhc, pqfwf,pqoce  ! fresh water and ice-ocean heat fluxes
      !!-------------------------- IN  -------------------------------------
      INTEGER,  INTENT(in) :: kt
      !!--------------------------------------------------------------------
      INTEGER :: jk
      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztfrz3d        ! freezing point
      REAL(wp), DIMENSION(jpi,jpj)     :: ztfrz          ! freezing point
      REAL(wp), DIMENSION(jpi,jpj)     :: ztavg          ! temperature avg
      !!----------------------------------------------------------------------
      !
      ! 0. ------------Mean freezing point
      DO jk = 1,jpk
         CALL eos_fzp(tsn(:,:,jk,jp_sal), ztfrz3d(:,:,jk), gdept_n(:,:,jk))
      END DO
      CALL isf_tbl(ztfrz3d, ztfrz, 'T', misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par )
      !
      ! 1. ------------Mean temperature
      CALL isf_tbl(tsn(:,:,jk,jp_tem), ztavg, 'T', misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par )
      !
      ! 2. ------------Net heat flux and fresh water flux due to the ice shelf
      pqoce(:,:) = - rau0 * rcp * rn_gammat0 * risfLeff(:,:) * e1t(:,:) * ( ztfrz(:,:) - ztavg(:,:) ) * r1_e1e2t(:,:)
      pqfwf(:,:) =   pqoce(:,:) / rLfusisf
      pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp     ! heat content flux 
      !
      ! output ttbl
      CALL iom_put('ttbl_par', ztavg(:,:) )
      !
      ! output thermal driving
      CALL iom_put('isfthermald_par',( ztfrz(:,:) - ztavg(:,:) ))
      !
      !
   END SUBROUTINE isfpar_mlt_bg03

   SUBROUTINE isfpar_mlt_oasis
   !TODO
   END SUBROUTINE isfpar_mlt_oasis

END MODULE isfparmlt