!> \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_sedim

! 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) :: h_sedimmx   ! sediment thickness (or rebounded bedrock)
real,dimension(nx,ny) :: h_sedimmy   ! sediment thickness (or rebounded bedrock)
real :: seuil_sedim   ! drag reduced beyond this threshold
real :: coef_sedim    ! drag reduction factor
real :: neffmin = 1.e5 ! as in neffect, but this should be the exact same variable

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

real,dimension(nx,ny) :: h_sedim   ! sediment thickness (or rebounded bedrock)
character(len=150) :: file_sedim   ! file for sediment thickness for HN or rebounded bsoc for Antar
character(len=150) :: file_ncdf    ! fichier netcdf issue des fichiers .dat

namelist/drag_param_beta_sedim/beta_slope,beta_expo,betamax,betamin,coef_ile,file_sedim,seuil_sedim,coef_sedim

if (itracebug.eq.1)  call tracebug(' dragging param avec sedim')


! 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_sedim)

write(num_rep_42,428)'!___________________________________________________________' 
write(num_rep_42,428) '&drag_param_beta_sedim          ! 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,*) 'file_sedim     = ', file_sedim
write(num_rep_42,*) 'seuil_sedim    = ', seuil_sedim
write(num_rep_42,*) 'coef_sedim     = ', coef_sedim
write(num_rep_42,*)'/'                            
write(num_rep_42,428) '! Slope & expo of beta = - slope x Neff ** expo'
write(num_rep_42,428) '! Where h_sedim > seuil_sedim, beta*coef_sedim'

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


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

! modification of basal drag depends on where we have sediments
! for Antarctica, we assume that what is present-day below sea level has sediment
      file_sedim=trim(dirnameinp)//trim(file_sedim)
      call lect_input(1,'h_sedim',1,h_sedim(:,:),file_sedim,file_ncdf)
 
      h_sedimmx(1,:)=h_sedim(1,:)
      h_sedimmy(:,1)=h_sedim(:,1)
      do i=2,nx
         h_sedimmx(i,:)=(h_sedim(i-1,:)+h_sedim(i,:))/2.
      enddo
      do j=2,ny
         h_sedimmy(:,j)=(h_sedim(:,j-1)+h_sedim(:,j))/2.
      enddo
 
      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)=neffmin
neff(nx,:)=neffmin
! 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) !/neffmin
betamy(:,:)= beta_slope*(neffmy(:,:)**beta_expo) !/neffmin

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

where (h_sedimmx(:,:).gt.seuil_sedim) betamx(:,:) = betamx(:,:) * coef_sedim
where (h_sedimmy(:,:).gt.seuil_sedim) betamy(:,:) = betamy(:,:) * coef_sedim
!where (h_sedimmx(:,:).gt.seuil_sedim) betamx(:,:) = beta_slope*(neffmin**beta_expo)/neffmin
!where (h_sedimmy(:,:).gt.seuil_sedim) betamy(:,:) = beta_slope*(neffmin**beta_expo)/neffmin

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_sedim