!> \file prescribe-H-i2s_mod.f90
!!
!!  Determine mask and thickness for prescription of H 
!!  include ice2sea prescribed retreat  (Cat april 2012)
!<

!> \namespace prescribe_H
!! 
!! Determine mask and thickness for prescription of H 
!!  
!! \author CatRitz
!! \date june 2009
!!
!! @note use  module3D_phy
!! @note use param_phy_mod
!<

module prescribe_H

  use module3D_phy
  use param_phy_mod
  use interface_input

! pour les reculs de type ice2sea
!  use toy_retreat_mod

  implicit none
!  real,dimension(nx,ny)    :: Hp           !< H value if prescribed
!  real,dimension(nx,ny)    :: Hp0          !< H value if prescribed (reference value)
!  real,dimension(nx,ny)    :: Delta_H      !< Delta_H value if prescribed
!  integer,dimension(nx,ny) :: i_delta_H    !< 1 if Delta_H is prescribed on this node, else 0
!  integer,dimension(nx,ny) :: i_Hp         !< 1 if H is prescribed on this node, else 0
!  integer,dimension(nx,ny) :: i_Hp0        !< i_hp mask reference value does not change with time
!  integer,dimension(nx,ny) :: MK_gl0       !< mask grounding line initial
!  integer,dimension(nx,ny) :: MK_flot0     !< mask float initial

! pour grounding line retreat, ice2sea
! nouvelle version interactive pour imposer des contraintes variees.
! sanity check, temps de depart en fonction des bassins, ....


! nouvelle version recul grounding line pour ice2sea
! ------------------------------------------------------



contains

  !______________________________________________________________________________________
  !> initialise prescribe H
  !! calculate the initial grounding line position
  !! @return MK_flot0 and Mk_gl0

  subroutine init_prescribe_H

    implicit none
    integer    :: voisin   !< pour test sur des masques

    if (itracebug.eq.1)  call tracebug(' Entree dans routine init_prescribe_H')

    ! keep the mask of floating points
    where (flot(:,:))
       MK_flot0(:,:)= 1
    elsewhere 
       MK_flot0(:,:)= 0
    end where


    ! determine the initial grounding line

    MK_gl0(:,:)=0
    HP0(:,:)=H0(:,:)


    do j=2,ny-1
       do i=2,nx-1
          voisin = MK_flot0(i-1,j) + MK_flot0(i+1,j) + MK_flot0(i,j-1) + MK_flot0(i,j+1)
          if ((MK_flot0(i,j).eq.0).and.(voisin.gt.0))  then  ! grounded and at least one neighbour floating
             MK_gl0(i,j) = 1                                ! tagged mask grounding line

          end if
       end do
    end do


    call prescribe_grid_boundary_0


! Dans le cas du Groenland ice2sea, mk_init = 3 -> noeuds qui ne bougent pas

    if (geoplace(1:4).eq.'GI2S') then

       do j = 1,ny
          do i = 1,nx
             if (mk_init(i,j).eq.3) then
                i_hp0 (i,j) = 1
                Hp0  (i,j)  = H0(i,j)
             endif
          end do
       end do
    end if

  end subroutine init_prescribe_H
  !______________________________________________________________________________________
  !>  function prescribe_present_H_gl
  !! calculate the initial grounding line position
  !! @return i_hp(:,:) and hp(:,:)

  subroutine prescribe_present_H_gl

    implicit none

    if (itracebug.eq.1)  call tracebug(' Entree dans routine  prescribe_present_H_gl')

    where ((MK_flot0(:,:).eq. 1).or.(MK_gl0(:,:) .eq. 1)) ! floating or grounding line
       i_hp(:,:) = 1                                      ! thickness prescribed to present value
       hp(:,:)   = Hp0(:,:)
    end where
if (itracebug.eq.1)  call tracebug(' fin prescribe_present_H_gl')


  end subroutine prescribe_present_H_gl
  
   !______________________________________________________________________________________
  !>  function prescribe_present_H_gl_bmelt
  !! calculate the initial grounding line position
  !! only grounding line points are prescribed to present value
  !! @return i_hp(:,:) and hp(:,:)

  subroutine prescribe_present_H_gl_bmelt

    implicit none

    if (itracebug.eq.1)  call tracebug(' Entree dans routine  prescribe_present_H_gl_bmelt')

!    where (MK_gl0(:,:) .eq. 1) ! grounding line only !cdc pour calcule fonte basale
!       i_hp(:,:) = 1                                 ! thickness prescribed to present value
!       hp(:,:)   = Hp0(:,:)
!    end where
    i_hp(:,:) = 0
if (itracebug.eq.1)  call tracebug(' fin prescribe_present_H_gl_bmelt')


  end subroutine prescribe_present_H_gl_bmelt

  !______________________________________________________________________________________
  !>  function prescribe_fixed_points
  !!  fixes  thickness for some points
  !! @return i_hp(:,:) and hp(:,:)

  subroutine prescribe_fixed_points

if (itracebug.eq.1)  call tracebug(' Entree dans routine  prescribe_fixed_points')
    where (i_hp0(:,:).eq.1)             ! les points i_hp0 le sont pour tout le run
       i_hp(:,:) = i_hp0(:,:)
       Hp(:,:)   = Hp0(:,:)
    end where

if (itracebug.eq.1)  call tracebug(' fin prescribe_fixed_points')


  end subroutine prescribe_fixed_points

  !______________________________________________________________________________________
  !>  function prescribe_paleo_gl_shelf
  !! calculate the  grounding line thickness
  !! no ice shelves
  !! @return i_hp(:,:) and hp(:,:)

  subroutine prescribe_paleo_gl_shelf

    implicit none

    if (itracebug.eq.1)  call tracebug(' Entree dans routine prescribe_paleo_gl_shelf')


!   noeuds qui doivent être imposés

    where ((MK_flot0(:,:).eq. 1).or.(MK_gl0(:,:) .eq. 1)) ! floating or grounding line
       i_hp(:,:) = 1                                      ! thickness prescribed
    end where

! valeur imposee
    where (MK_flot0(:,:).eq. 1)                          ! paleo shelf epaisseur a 1
       hp(:,:)   = 1.
    end where


    where (MK_gl0(:,:) .eq. 1)
       hp(:,:) = (sealevel_2d(:,:)-Bsoc(:,:))*row/ro +1.
    end where


  end subroutine prescribe_paleo_gl_shelf

  !______________________________________________________________________________________
  !>  function prescribe_grid_boundary
  !! prescribe the grid boundary thickness to 0
  !! @return i_hp(:,:) and hp(:,:)

  subroutine prescribe_grid_boundary_0

    implicit none
    if (itracebug.eq.1)  call tracebug(' Entree dans routine prescribe_grid_boundary_0')

    ! prescribe thickness on the grid boundaries
    i_hp(:,:) = 0

    i_hp(1,:) = 1
    i_hp(nx,:)= 1
    hp(1,:)   = 0
    hp(nx,:)  = 0

    i_hp(:,1) = 1
    i_hp(:,ny)= 1
    hp(:,1)   = 0
    hp(:,ny)  = 0

! valeurs de reference
    i_hp0(1,:)  = 1
    i_hp0(nx,:) = 1
    hp0(1,:)    = 0
    hp0(nx,:)   = 0

    i_hp0(:,1)  = 1
    i_hp0(:,ny) = 1
    hp0(:,1)    = 0
    hp0(:,ny)   = 0


  end subroutine prescribe_grid_boundary_0
  !______________________________________________________________________________________

  !______________________________________________________________________________________
  !> break ice shelves
  !! 
  !!
  subroutine break_all_ice_shelves

    implicit none
    if (itracebug.eq.1)  call tracebug(' Entree dans routine break_all_ice_shelves ')


    debug_3D(:,:,89)=0.

    where (flot(:,:))                           ! floating from the beginning and eventually from grounding line retreat
       i_hp(:,:) = 1                                      ! thickness prescribed to 1 m
       hp (:,:)  = 1.      
       hp0 (:,:)  = 1.    
       H(:,:)     =1.
       debug_3D(:,:,89)=1.

    end where
    where (MK_flot0(:,:).eq. 1)
       i_hp(:,:) = 1                                      ! thickness prescribed to 1 m
       hp (:,:)  = 1.  
       hp0 (:,:)  = 1.   
       H(:,:)     =1.
    end where


  end subroutine break_all_ice_shelves

  !
  !______________________________________________________________________________________
  !> break ice shelves
  !! 
  !!
  subroutine melt_ice_shelves

    implicit none
    integer :: nbvoisins,nbdeglac,iv,jv
    real    :: hmoy,hmin,hmax

    if (itracebug.eq.1)  call tracebug(' Entree dans routine melt_ice_shelves')

    do j = 2,ny-1
       do i=2, nx-1
          if  (flot(i,j)) then                                 ! floating 
             i_hp(i,j) = 1                                      ! thickness prescribed 
             hmoy=0.
             hmin=0.
             hmax=0.
             nbvoisins=0
             nbdeglac=0
             do jv=-1,1
                do iv=1,1
                   if (flot(i+iv,j+jv)) then
                      nbvoisins=nbvoisins+1
                      hmoy=hmoy+h(i+iv,j+jv)
                      hmin=min(hmin,h(i+iv,j+jv))
                      hmax=max(hmax,h(i+iv,j+jv))
                      if (h(i+jv,j+jv).lt.10.) nbdeglac=nbdeglac+1
                   end if
                end do
             end do
             hmoy=hmoy/(max(nbvoisins,1))

             if(H(i,j).gt.hmoy) then  ! melt only if smooth
                hp(i,j)=H(i,j)-1.*dt
             else if (hmax-hmin.lt.10.) then 
                hp(i,j)=H(i,j)-1.*dt
             else if ((H(i,j).lt.50.).and.(Hmin.lt.30.)) then
                hp(i,j)=1.
             else if ((H(i,j).lt.50.).and.(nbdeglac.ge.4)) then
                hp(i,j)=1.
             endif
             hp(i,j)   = max(hp(i,j),1.)
          end if
       end do
    end do

  end subroutine melt_ice_shelves


   !______________________________________________________________________________________
  !>  function prescribe_present_H_gl
  !! calculate the initial grounding line position
  !! @return i_hp(:,:) and hp(:,:)

  subroutine prescribe_present_H_gl_copy

    implicit none
    if (itracebug.eq.1)  call tracebug(' Entree dans routine  prescribe_present_H_gl-copy')

    do j=1,ny
       do i=1,nx
          if ((MK_flot0(i,j).eq. 1).or.(MK_gl0(i,j) .eq. 1)) then
             i_hp(i,j) = 1                                      ! thickness prescribed to present value
             hp(i,j)   = Hp0(i,j)
          end if
       end do
    end do

!!$where ((MK_flot0(:,:).eq. 1).or.(MK_gl0(:,:) .eq. 1)) ! floating or grounding line
!!$   i_hp(:,:) = 1                                      ! thickness prescribed to present value
!!$!   hp(:,:)   = 1   !Hp0(:,:)
!!$end where


  end subroutine prescribe_present_H_gl_copy
!______________________________________________________________________________________


end module prescribe_H