MODULE isfdynnxt
   !!=========================================================================
   !!                       ***  MODULE  isfnxt  ***
   !! Ice shelf update: compute the dynnxt ice shelf contribution
   !!=========================================================================
   !! History :  OPA  !  2019-09  (P. Mathiot)  Original code
   !!-------------------------------------------------------------------------
  
   !!-------------------------------------------------------------------------
   !!   isfnxt       : aplly correction needed for the ice shelf to ensure conservation
   !!-------------------------------------------------------------------------

   USE isf_oce

   USE phycst , ONLY: r1_rau0                ! physical constant
   USE dom_oce, ONLY: e3t_b, e3t_n, r1_e1e2t ! time and space domain

   USE in_out_manager

   IMPLICIT NONE

   PRIVATE

   PUBLIC isf_dynnxt

CONTAINS

   SUBROUTINE isf_dynnxt ( kt, pcoef )
      !!--------------------------------------------------------------------
      !!                  ***  ROUTINE isf_dynnxt  ***
      !!
      !! ** Purpose : compute the ice shelf volume filter correction for cavity, param, ice sheet coupling case
      !!
      !!-------------------------- OUT -------------------------------------
      INTEGER ,                     INTENT(in   ) :: kt
      !
      REAL(wp),                     INTENT(in   ) :: pcoef           ! atfp * rdt * r1_rau0
      !!--------------------------------------------------------------------
      INTEGER :: jk  ! loop index
      !!--------------------------------------------------------------------
      !
      ! ice shelf cavity
      IF ( ln_isfcav_mlt ) CALL isf_dynnxt_mlt(misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, fwfisf_cav, fwfisf_cav_b, pcoef)
      !
      ! ice shelf parametrised
      IF ( ln_isfpar_mlt ) CALL isf_dynnxt_mlt(misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, fwfisf_par, fwfisf_par_b, pcoef)
      !
      IF ( ln_isfcpl .AND. ln_rstart .AND. kt == nit000+1 ) THEN
         DO jk = 1, jpkm1
            e3t_b(:,:,jk) =   e3t_b(:,:,jk) - pcoef * risfcpl_vol(:,:,jk) * r1_e1e2t(:,:)
         END DO
      END IF
      !
   END SUBROUTINE isf_dynnxt

   SUBROUTINE isf_dynnxt_mlt ( ktop, kbot, phtbl, pfrac, pfwf, pfwf_b, pcoef )
      !!--------------------------------------------------------------------
      !!                  ***  ROUTINE isf_dynnxt_mlt  ***
      !!
      !! ** Purpose : compute the ice shelf volume filter correction for cavity or param
      !!
      !!-------------------------- IN  -------------------------------------
      INTEGER , DIMENSION(jpi,jpj), INTENT(in   ) :: ktop , kbot     ! top and bottom level of tbl
      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pfrac, phtbl    ! fraction of bottom cell included in tbl, tbl thickness
      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pfwf , pfwf_b   ! now/before fwf
      REAL(wp),                     INTENT(in   ) :: pcoef           ! atfp * rdt * r1_rau0
      !!----------------------------------------------------------------------
      INTEGER :: ji,jj,jk
      REAL(wp), DIMENSION(jpi,jpj) :: zfwfinc
      !!----------------------------------------------------------------------
      !
      ! compute fwf conservation correction
      zfwfinc(:,:) = pcoef * ( pfwf_b(:,:) - pfwf(:,:) ) / phtbl(:,:) * r1_rau0
      !
      ! add the increment in the tbl
      DO jk = 1, jpkm1
         DO jj = 1, jpj
            DO ji = 1, jpi
               IF( ktop(ji,jj) <= jk .AND. jk < kbot(ji,jj)  ) THEN
                  e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zfwfinc(ji,jj) * e3t_n(ji,jj,jk)
               ELSEIF ( jk == kbot(ji,jj) ) THEN
                  e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zfwfinc(ji,jj) * e3t_n(ji,jj,jk) * pfrac(ji,jj)
               ENDIF
            END DO
         END DO
      END DO
      !
   END SUBROUTINE isf_dynnxt_mlt

END MODULE isfdynnxt