!> \file dragging_map_Ant_mod.f90
!! Definition des zones de streams apres lecture d'une carte evaluee offline
!<

!> \namespace dragging_map_Ant
!! Definition zone stream. A faire off line
!! \author ...
!! \date ...
!! @note active les zones de stream avec :
!! @note   - un critere sur la hauteur d'eau 
!! @note   - un critere de fleuves actuels
!! @note Attention il faut fonctionner en additionnant SIA et Stream pour éviter les pb
!! de fleuves en pointille
!! @note Used module
!! @note   - use module3D_phy
!! @note   - use param_phy_mod
!! @note peut consdidérer des lois differentes
!<


module dragging_map_allow

! Defini les zones de stream avec :
! * un critere sur la hauteur d'eau 
! * un critere de fleuves actuels (determine off line)
! attention il faut fonctionner en additionnant SIA et Stream pour éviter les pb
! de fleuves en pointille


use module3d_phy
use param_phy_mod


implicit none
logical,dimension(nx,ny) :: fleuvemx         !< fleuves sont les tableaux courants (dep. time)   
logical,dimension(nx,ny) :: fleuvemy
logical,dimension(nx,ny) :: fleuve      
!logical,dimension(nx,ny) :: cote    
!real, dimension(nx,ny) :: Vcol_x              !< uniquement pour compatibilite dans le spinup cat (aurel)
!real, dimension(nx,ny) :: Vcol_y              !< uniquement pour compatibilite dans le spinup cat (aurel)
!real, dimension(nx,ny) :: Vsl_x               !< pour compatibilite dragging beta (aurel)
!real, dimension(nx,ny) :: Vsl_y               !< pour compatibilite dragging beta (aurel)

!logical :: corr_def = .false.                 !< for deformation correction, pour compatibilite beta (aurel)

character(len=100) :: masque_stream        !<  file masque streams autorises (sur noeuds majeurs)

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


real :: tostick      ! pour la glace posee
real :: tob_ile      ! pour les iles

contains
!-------------------------------------------------------------------------------


!> SUBROUTINE: init_dragging
!! Cette routine fait l'initialisation du dragging.
!<   

subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.

implicit none

integer :: bidon

 
namelist/drag_map_allow/hwatstream,cf,betamax,toblim,masque_stream

if (itracebug.eq.1) write(num_tracebug,*)' init_dragging avec dragging_neem'


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

! lecture des parametres du run                    
!--------------------------------------------------------------------

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

write(num_rep_42,428) '!___________________________________________________________' 
write(num_rep_42,428) '! module dragging_neem '
write(num_rep_42,drag_neem)                    
write(num_rep_42,428) '! hwatstream (m) :  critere de passage en stream vitesses de bilan'
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) '! seulement pour les points cotiers'
write(num_rep_42,428) '! betamax : (Pa) frottement maxi sous les streams '
write(num_rep_42,428) '! toblim : (Pa)  pour les iles '
write(num_rep_42,428) '! masque_stream : nom du fichier qui contient le masque de streams autorises'
write(num_rep_42,*)

tostick=1.e5   ! valeurs pour les points non flgzmx
tob_ile=betamax/2.
moteurmax=toblim
!moteurmax=0.

!-------------------------------------------------------------------
! masque stream : lecture du masque
! attention on prend la colonne 1 qui correspond aux noeuds majeurs
masque_stream  = trim(dirnameinp)//trim(masque_stream)
print*, 'le fichier des streams ', masque_stream
call lect_datfile(nx,ny,mstream,1,masque_stream)     
open (20,file=masque_stream)
read(20,*)
do j=1,ny
   do i=1,nx
      read(20,*) mstream(i,j),bidon, bidon! mstream_mx(i:i+1,j), mstream_my(i,j:j+1)
   end do
end do
close(20)

do j=1,ny
   do i=1,nx
      if (mstream(i,j).eq.1) then ! allowed ice streams
         mstream_mx(i:i+1,j) = 1
         mstream_my(i,j:j+1) = 1
      end if
   end do
end do


debug_3D(:,:,25) = mstream(:,:)


! 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(rog*Hmx(i,j)*sdx(i,j)).lt.toblim)
      ilemy(i,j)=ilemy(i,j).and.(abs(rog*Hmy(i,j)*sdy(i,j)).lt.toblim)
   end do
end do

fleuvemx(:,:)=.false.
fleuvemy(:,:)=.false.
fleuve(:,:)=.false.
cote(:,:)=.false.
cotemx(:,:)=.false.
cotemy(:,:)=.false.

! detection des cotes
do  j=2,ny-1  
   do i=2,nx-1  
      if ((.not.flot(i,j)).and. & 
      ((flot(i+1,j)).or.(flot(i,j+1)).or.(flot(i-1,j)).or.(flot(i,j-1)))) then
         cote(i,j)=.true.
         cotemx(i:i+1,j) = .true.
         cotemy(i,j:j+1) = .true.
      endif
   end do
end do

! le calcul des fleuves se fait sur les mailles majeures
! normalement, les points cotiers sont deja tagges gzmx dans la routine flottab

do j=2,ny-1
   do i=2,nx-1
      
      if ((hwater(i,j).gt.hwatstream).and.(mstream(i,j).eq.1).and.(.not.flot(i,j))) then
         fleuve(i,j)=.true.
         fleuvemx(i:i+1,j)=.true.
         fleuvemy(i,j:j+1)=.true.
         gzmx(i:i+1,j)=.true.
         gzmy(i,j:j+1)=.true.
      endif
   end do
end do


! calcul du frottement basal (ce bloc etait avant dans neffect)

drag_mx : do j=2,ny-1
   do i=2,nx-1

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

      if (cotemx(i,j)) then                        ! point cotier peut etre diminue si pression d'eau
         betamx(i,j)=cf*neffmx(i,j) 
         betamx(i,j)=min(betamx(i,j),betamax)
      else if ((gzmx(i,j)).and.(.not.cotemx(i,j))) then  ! stream -> betamax
         betamx(i,j)=betamax
         betamx(i,j)=min(betamx(i,j),betamax)
         betamx(i,j)=max(betamx(i,j),10.)

!         if (cf*neffmx(i,j).gt.1500.) then
!            gzmx(i,j)=.false.
!            betamx(i,j)=tostick
!         endif

      else if (ilemx(i,j)) then                     ! ile : tob_ile=betamax/2 + peut etre diminue si pression d'eau
         betamx(i,j)=cf*neffmx(i,j) 
         betamx(i,j)=min(betamx(i,j),tob_ile)       
      else if (flotmx(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 drag_mx

drag_my : do j=2,ny-1
   do i=2,nx-1

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

      if (cotemy(i,j)) then                        ! point cotier peut etre diminue si pression d'eau
         betamy(i,j)=cf*neffmy(i,j) 
         betamy(i,j)=min(betamy(i,j),betamax)
      else if ((gzmy(i,j)).and.(.not.cotemy(i,j))) then  ! stream -> betamax
         betamy(i,j)=betamax
         betamy(i,j)=min(betamy(i,j),betamax)
         betamy(i,j)=max(betamy(i,j),10.)

!         if (cf*neffmy(i,j).gt.1500.) then
!            gzmy(i,j)=.false.
!            betamy(i,j)=tostick
!         endif

      else if (ilemy(i,j)) then                     ! ile : tob_ile=betamax/2 + peut etre diminue si pression d'eau
         betamy(i,j)=cf*neffmy(i,j) 
         betamy(i,j)=min(betamy(i,j),tob_ile)       
      else if (flotmy(i,j)) then                    ! flottant ou ile
         betamy(i,j)=0.
      else                                         ! grounded, SIA
         betamy(i,j)=tostick                       ! frottement glace posee (1 bar)
      endif

   end do
end do drag_my


where (fleuve(:,:))
   debug_3D(:,:,1)=1
elsewhere (flot(:,:))
   debug_3D(:,:,1)=2
elsewhere
   debug_3D(:,:,1)=0
endwhere

where (cote(:,:))
   debug_3D(:,:,2)=1
elsewhere
   debug_3D(:,:,2)=0
endwhere

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

debug_3D(:,:,23)= abs(RO*G*HMX(:,:)*sdx(:,:)*1.e-5) 
debug_3D(:,:,24)= abs(RO*G*HMY(:,:)*sdy(:,:)*1.e-5) 

return
end subroutine dragging

end module dragging_neem