!> \file dragging_neff_contmaj_mod.f90
!! Defini les zones de stream avec 3 criteres 
!<

!> \namespace  dragging_neff_contmaj
!! Defini les zones de stream avec 3 criteres
!! \author ...
!! \date ...
!! @note Les trois criteres sont:
!! @note   * un critere sur la pression effective 
!! @note   * un critere de continuite depuis la cote
!! @note   * un critere sur la courbure du socle (si negatif, vallees) 
!! @note Used module
!! @note   - use module3D_phy
!! @note   - use param_phy_mod
!<

module dragging_param_beta

! Defini les zones de stream avec :
! * un critere sur la pression effective 
! * un critere de continuite depuis la cote
! * un critere sur la courbure du socle (si negatif, vallees) 

use module3d_phy
use param_phy_mod
use interface_input
use io_netcdf_grisli
use fake_beta_iter_vitbil_mod

implicit none

logical,dimension(nx,ny) :: cote

real,dimension(nx,ny) :: beta_param ! local variable

real :: betamin   ! betamin : (Pa) frottement mini sous les streams

real :: beta_slope        ! = A in:     B = A x Neff ** K
real :: beta_expo         ! = K in:     B = A x Neff ** K

real,dimension(nx,ny) :: neff   ! pression effective noeuds majeurs
real,dimension(nx,ny) :: hwatmx ! hauteur d'eau staggered grille - afq
real,dimension(nx,ny) :: hwatmy ! hauteur d'eau staggered grille - afq

real :: coef_ile      ! coef frottement zones iles   (0.1)

real, dimension(nx,ny) :: Vcol_x           !< uniquement pour compatibilite avec spinup cat
real, dimension(nx,ny) :: Vcol_y           !< uniquement pour compatibilite avec spinup cat
real, dimension(nx,ny) :: Vsl_x            !< uniquement pour compatibilite avec spinup cat
real, dimension(nx,ny) :: Vsl_y            !< uniquement pour compatibilite avec spinup cat
logical :: corr_def = .false.              !< for deformation correction, pour compatibilite beta

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_param_beta/beta_slope,beta_expo,betamax,betamin,coef_ile

if (itracebug.eq.1)  call tracebug(' dragging avec neff et slope')


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

! lecture des parametres du run                      block drag neff
!--------------------------------------------------------------------

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

write(num_rep_42,428)'!___________________________________________________________' 
write(num_rep_42,428) '&drag_param_beta          ! nom du bloc dragging param beta'
write(num_rep_42,*)
write(num_rep_42,*) 'beta_slope            = ', beta_slope
write(num_rep_42,*) 'beta_expo             = ', beta_expo
write(num_rep_42,*) 'betamax               = ', betamax
write(num_rep_42,*) 'betamin               = ', betamin
write(num_rep_42,*)'/'                            
write(num_rep_42,428) '! Slope & expo of beta = - slope x Neff ** expo'

inv_beta=0
!-------------------------------------------------------------------
! masque stream
mstream_mx(:,:)=1
mstream_my(:,:)=1


! coefficient permettant de modifier le basal drag.
drag_mx(:,:)=1.
drag_my(:,:)=1.
      
betamax_2d(:,:) = betamax
   
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
!$ USE OMP_LIB

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

!$OMP PARALLEL
!$OMP DO
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
!$OMP END DO


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


!$OMP WORKSHARE
fleuvemx(:,:)=.false.
fleuvemy(:,:)=.false.
cote(:,:)=.false.
gzmx(:,:)=.true.
gzmy(:,:)=.true.
flgzmx(:,:)=.false.
flgzmy(:,:)=.false.
!$OMP END WORKSHARE

! detection des cotes
!$OMP DO
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.
      endif
   end do
end do
!$OMP END DO

! calcul de neff (pression effective sur noeuds majeurs)
if (sum(neffmx(:,:)).le.0.) neffmx(:,:) =1.e8
if (sum(neffmy(:,:)).le.0.) neffmy(:,:) =1.e8

!$OMP DO
do j=1,ny-1
   do i=1,nx-1
        neff(i,j)=(neffmx(i,j)+neffmx(i+1,j)+neffmy(i,j)+neffmy(i,j+1))/4
   enddo
enddo
!$OMP END DO
!aurel, for the borders:
!$OMP WORKSHARE
neff(:,ny)=1.e5
neff(nx,:)=1.e5
! calcul de hwat (staggered)
hwatmx(:,:)=0.
hwatmy(:,:)=0.
!$OMP END WORKSHARE
!$OMP DO
do i=2,nx
   hwatmx(i,:)=(hwater(i-1,:)+hwater(i,:))/2.
enddo
!$OMP END DO
!$OMP DO
do j=2,ny
   hwatmy(:,j)=(hwater(:,j-1)+hwater(:,j))/2.
enddo
!$OMP END DO

!$OMP WORKSHARE

! new parametrisation of beta on Neff:
betamx(:,:)= beta_slope*(neffmx(:,:)**beta_expo)
betamy(:,:)= beta_slope*(neffmy(:,:)**beta_expo)

where (ilemx(:,:)) betamx(:,:) = betamx(:,:) * coef_ile
where (ilemy(:,:)) betamy(:,:) = betamy(:,:) * coef_ile


betamx(:,:)=max(betamx(:,:),betamin)
betamy(:,:)=max(betamy(:,:),betamin)
betamx(:,:)=min(betamx(:,:),betamax)
betamy(:,:)=min(betamy(:,:),betamax)

where ( hwatmx(:,:) .lt. 0.5 ) betamx(:,:) = betamax
where ( hwatmy(:,:) .lt. 0.5 ) betamy(:,:) = betamax

!$OMP END WORKSHARE


! calcul de gzmx 

! points flottants : flgzmx mais pas gzmx
!$OMP DO
do j=2,ny
   do i=2,nx
      if (flot(i,j).and.(flot(i-1,j))) then
         flotmx(i,j)=.true.
         flgzmx(i,j)=.true.
         gzmx(i,j)=.false.
         betamx(i,j)=0.
      end if
      if (flot(i,j).and.(flot(i,j-1))) then
         flotmy(i,j)=.true.
         flgzmy(i,j)=.true.
         gzmy(i,j)=.false.
         betamy(i,j)=0.
      end if
   end do
end do
!$OMP END DO

!--------- autres criteres

! points poses
!$OMP WORKSHARE
gzmx(:,:)=gzmx(:,:).and.(betamx(:,:).lt.betamax)   !  Pas de calcul pour les points
gzmy(:,:)=gzmy(:,:).and.(betamy(:,:).lt.betamax)   !  au fort beta

flgzmx(:,:) = flgzmx(:,:) .or. gzmx(:,:)
flgzmy(:,:) = flgzmy(:,:) .or. gzmy(:,:)

where (hmx(:,:).eq.0.)
  flgzmx(:,:) = .false.
endwhere
where (hmy(:,:).eq.0.)
  flgzmy(:,:) = .false.
endwhere

fleuvemx(:,:)=gzmx(:,:)
fleuvemy(:,:)=gzmy(:,:)
!$OMP END WORKSHARE

!$OMP DO
do j=2,ny-1
   do i=2,nx-1
      beta_centre(i,j) = ((betamx(i,j)+betamx(i+1,j)) &
          + (betamy(i,j)+betamy(i,j+1)))*0.25
   end do
end do
!$OMP END DO
!$OMP END PARALLEL

return
end subroutine dragging

end module dragging_param_beta