MODULE isfhdiv

   USE dom_oce
   USE isf
   USE phycst
   USE in_out_manager

   IMPLICIT NONE

   PRIVATE

   PUBLIC isf_hdiv

CONTAINS

   SUBROUTINE isf_hdiv( kt, phdiv )
      !!----------------------------------------------------------------------
      !!                  ***  SUBROUTINE isf_hdiv  ***
      !!       
      !! ** Purpose :   
      !!
      !! ** Method  :   
      !!
      !! ** Action  :   phdiv   decreased by the fwf inflow (isf melt in this case)
      !!----------------------------------------------------------------------
      REAL(wp), DIMENSION(:,:,:), INTENT( inout ) ::   phdiv   ! horizontal divergence
      !!----------------------------------------------------------------------
      INTEGER, INTENT(in) :: kt
      !
      ! ice shelf cavity contribution
      IF ( ln_isfcav_mlt ) CALL isf_hdiv_mlt(misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, fwfisf_cav, fwfisf_cav_b, phdiv)
      !
      ! ice shelf parametrisation contribution
      IF ( ln_isfpar_mlt ) CALL isf_hdiv_mlt(misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, fwfisf_par, fwfisf_par_b, phdiv)
      !
      ! ice sheet coupling contribution (if conservation needed)
      IF ( ll_isfcpl ) THEN
         !
         ! correct divergence only for the first time step
         IF ( kt == nit000 ) CALL isf_hdiv_cpl(risfcpl_vol, phdiv)
         !
         ! correct divergence every time step to remove any trend due to coupling
         IF ( ll_isfcpl_cons ) CALL isf_hdiv_cpl(risfcpl_cons_vol, phdiv)
         !
      END IF
      !
   END SUBROUTINE isf_hdiv

   SUBROUTINE isf_hdiv_mlt(ktop, kbot, phtbl, pfrac, pfwf, pfwf_b, phdiv)
      !!----------------------------------------------------------------------
      !!                  ***  SUBROUTINE sbc_isf_div  ***
      !!       
      !! ** Purpose :   update the horizontal divergence with the runoff inflow
      !!
      !! ** Method  :   
      !!                CAUTION : risf_tsc(:,:,jp_sal) is negative (outflow) increase the 
      !!                          divergence and expressed in m/s
      !!
      !! ** Action  :   phdivn   decreased by the runoff inflow
      !!----------------------------------------------------------------------
      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: phdiv
      !!----------------------------------------------------------------------
      INTEGER , DIMENSION(jpi,jpj), INTENT(in   ) :: ktop , kbot
      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pfrac, phtbl
      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pfwf , pfwf_b
      !!----------------------------------------------------------------------
      INTEGER  ::   ji, jj, jk   ! dummy loop indices
      INTEGER  ::   ikt, ikb 
      REAL(wp), DIMENSION(jpi,jpj) :: zqvol,ztmp
      !!----------------------------------------------------------------------
      !
      !==   fwf distributed over several levels   ==!
      !
      ! compute integrated divergence correction
      zqvol(:,:) = 0.5_wp * ( pfwf(:,:) + pfwf_b(:,:) ) * r1_rau0 / phtbl(:,:)
      !
      ! update divergence at each level affected by ice shelf top boundary layer
      DO jj = 1,jpj
         DO ji = 1,jpi
               ikt = ktop(ji,jj)
               ikb = kbot(ji,jj)
               ! level fully include in the ice shelf boundary layer
               DO jk = ikt, ikb - 1
                  phdiv(ji,jj,jk) = phdiv(ji,jj,jk) + zqvol(ji,jj)
               END DO
               ! level partially include in ice shelf boundary layer 
               phdiv(ji,jj,ikb) = phdiv(ji,jj,ikb) + zqvol(ji,jj) * pfrac(ji,jj)
         END DO
      END DO
      !
   END SUBROUTINE isf_hdiv_mlt

   SUBROUTINE isf_hdiv_cpl(pqvol, phdiv)
      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: phdiv
      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pqvol
   
   INTEGER :: jk
      
      DO jk=1,jpk 
         phdiv(:,:,jk) =  phdiv(:,:,jk) + pqvol(:,:,jk) * r1_e1e2t(:,:) / e3t_n(:,:,jk)
      END DO

   END SUBROUTINE

END MODULE isfhdiv