!> \file dragging_hwat-contigu_mod.f90
!! Defini les zones de stream avec 2 criteres
!<

!> \namespace  dragging_hwat_cont
!! Defini les zones de stream avec 2 criteres
!! \author ...
!! \date ... 
!! @note Les 2 criteres sont: 
!! @note   * un critere sur la hauteur d'eau
!! @note   * un critere de continuite
!! depuis la cote
!! @note Used module
!! @note   - use module3D_phy
!<
module dragging_hwat_cont

! Defini les zones de stream avec un critere sur la hauteur d'eau + un critere de continuite
! depuis la cote

use module3d_phy
use fake_beta_iter_vitbil_mod

implicit none
!logical,dimension(nx,ny) :: fleuvemx
!logical,dimension(nx,ny) :: fleuvemy
real,dimension(nx,ny) :: hwatmx
real,dimension(nx,ny) :: hwatmY

real :: valmax
integer :: imax,jmax
integer :: i_moins1,i_plus1,j_moins1,j_plus1
integer :: lmax=20
integer :: idep,jdep,iloc,jloc

real :: tostick

contains
!-------------------------------------------------------------------------------
!> SUBROUTINE: init_dragging
!! Cette routine fait l'initialisation du dragging.
!>
subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.

implicit none
 
namelist/drag_hwat_cont/hwatstream,cf,betamax,toblim,tostick    

if (itracebug.eq.1)  call tracebug(' Dragging avec hwatermax')


! formats pour les ecritures dans 42
428 format(A)

! lecture des parametres du run                      block draghwat
!--------------------------------------------------------------------

rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
read(num_param,drag_hwat_cont)

write(num_rep_42,428)'!___________________________________________________________' 
write(num_rep_42,428) '&drag_hwat_cont              ! nom du bloc hwatermax '
write(num_rep_42,*)
write(num_rep_42,*) 'hwatsream      = ',hwatstream 
write(num_rep_42,*) 'cf             = ',cf
write(num_rep_42,*) 'betamax         = ', betamax
write(num_rep_42,*) 'toblim         = ', toblim
write(num_rep_42,*)'/'                            
write(num_rep_42,428) '! hwatstream (m) :  critere de passage en stream en partant de la cote'
write(num_rep_42,428) '!  si hwater > hwatstream '
write(num_rep_42,428) '! cf coefficient de la loi de frottement fonction Neff'
write(num_rep_42,428) '! betamax : (Pa) frottement maxi '
write(num_rep_42,428) '! toblim : (Pa) tes pour les iles '
write(num_rep_42,*)

inv_beta=0
tostick=1.e5   ! valeurs pour les points non flgzmx

!-------------------------------------------------------------------
! masque stream

      mstream_mx(:,:)=1
      mstream_my(:,:)=1


! coefficient permettant de modifier le basal drag.
      drag_mx(:,:)=1.
      drag_my(:,:)=1.

   
      return
      end subroutine init_dragging
!________________________________________________________________________________

!> SUBROUTINE: dragging
!! Defini la localisation des streams et le frottement basal
!>
!-------------------------------------------------------------------------
subroutine dragging   ! defini la localisation des streams et le frottement basal


!         les iles n'ont pas de condition neff mais ont la condition toblim
!         (ce bloc etait avant dans flottab)


do j=2,ny
   do i=2,nx
      ilemx(i,j)=ilemx(i,j).and.(abs(tobmx(i,j)).lt.toblim)
      ilemy(i,j)=ilemy(i,j).and.(abs(tobmy(i,j)).lt.toblim)
   end do
end do

! calcul de la hauteur d'eau sur les mailles mineures
do j=1,ny
   do i=2,nx
      hwatmx(i,j)=0.5*(hwater(i,j)+hwater(i-1,j))
   end do
end do

do j=2,ny
   do i=1,nx
      hwatmy(i,j)=0.5*(hwater(i,j)+hwater(i,j-1))
   end do
end do


debug_3D(:,:,1)=hwatmx(:,:)
debug_3D(:,:,2)=hwatmy(:,:)


! le gzmx initial a ete calculé dans flottab pour les points a cheval sur la grounding line

! Partant de la cote, le fleuve est determine en remontant vers le point de hauteur d'eau max
! tant que ce point est superieur a hwatermax.

! Attention : ce systeme ne permet pas d'avoir des fleuves convergents
!             et les fleuves n'ont une épaisseur que de 1 noeud.

fleuvemx(:,:)=.false.
fleuvemy(:,:)=.false.

fleuve_mx: do j=1,ny  
   do i=2,nx                      ! le noeud 1 n'est pas attribue

cote_mx :  if (gzmx(i,j)) then    ! a ce niveau les gzmx sont sur la grounding line
         idep=i
         jdep=j

suit_x : do l=1,lmax         ! debut de la boucle de suivi, lmax longueur maxi des fleuves
           i_moins1=max(idep-1,2)
           j_moins1=max(jdep-1,1)
           i_plus1=min(idep+1,nx)
           j_plus1=min(jdep+1,ny)
           
           valmax=0.

! recherche du max en excluant les points flottants
         do jloc=j_moins1,j_plus1
            do iloc=i_moins1,i_plus1
 
               if ((hwatmx(iloc,jloc).gt.valmax).and.(.not.flotmx(iloc,jloc))) then
                  imax=iloc
                  jmax=jloc
                  valmax=hwatmx(iloc,jloc)
               endif
            end do
         end do

         if (hwatmx(imax,jmax).gt.hwatstream) then
            fleuvemx(imax,jmax)=.true.
            idep=imax
            jdep=jmax
         end if
      end do suit_x

   end if cote_mx

if ((.not.ilemx(i,j)).and.(fleuvemx(i,j))) gzmx(i,j)=.true. 



! calcul du frottement basal (ce bloc etait avant dans neffect)
         if (gzmx(i,j)) then                          ! stream
            betamx(i,j)=cf*neffmx(i,j) 
         else if (flotmx(i,j).or.ilemx(i,j)) then     ! flottant ou ile
            betamx(i,j)=0.
         else                                         ! grounded, SIA
            betamx(i,j)=tostick                       ! frottement glace posee (1 bar)
         endif      

 end do
end do fleuve_mx




fleuve_my: do j=2,ny                  ! le noeud 1 n'est pas attribue
   do i=1,nx   

cote_my :  if (gzmy(i,j)) then    ! a ce niveau les gzmx sont sur la grounding line
         idep=i
         jdep=j

suit_y : do l=1,lmax         ! debut de la boucle de suivi, lmax longueur maxi des fleuves
           i_moins1=max(idep-1,1)
           j_moins1=max(jdep-1,2)
           i_plus1=min(idep+1,nx)
           j_plus1=min(jdep+1,ny)
           
           valmax=0.

! recherche du max en excluant les points flottants
         do jloc=j_moins1,j_plus1
            do iloc=i_moins1,i_plus1
 
               if ((hwatmy(iloc,jloc).gt.valmax).and.(.not.flotmy(iloc,jloc))) then
                  imax=iloc
                  jmax=jloc
                  valmax=hwatmy(iloc,jloc)
               endif
            end do
         end do

         if (hwatmy(imax,jmax).gt.hwatstream) then
            fleuvemy(imax,jmax)=.true.
            idep=imax
            jdep=jmax
         end if
      end do suit_y

   end if cote_my

if ((.not.ilemy(i,j)).and.(fleuvemy(i,j))) gzmy(i,j)=.true. 


 ! calcul du frottement basal (ce bloc etait avant dans neffect)
   if (gzmy(i,j)) then                ! stream
      betamy(i,j)=cf*neffmy(i,j) 
   else if (flotmy(i,j).or.ilemy(i,j)) then     ! flottant ou ile
      betamy(i,j)=0.
   else                               ! grounded, SIA
      betamy(i,j)=tostick             ! frottement glace posee
   endif


 end do
end do fleuve_my



where (fleuvemx(:,:))
   debug_3D(:,:,3)=1
elsewhere
   debug_3D(:,:,3)=0
endwhere

where (flotmx(:,:))
   debug_3D(:,:,3)=-1
endwhere

!_____________________


where (fleuvemy(:,:))
   debug_3D(:,:,4)=1
elsewhere
   debug_3D(:,:,4)=0
endwhere

where (flotmy(:,:))
   debug_3D(:,:,4)=-1
end where


return
end subroutine dragging

end module dragging_hwat_cont