!> \file slope_surf.f90
!!  File to compute the surface slope by finite difference
!<

!> SUBROUTINE:  slope_surf
!! Compute the surface slope by finite difference
!! \author Catritz
!! \date ...
!! @note Used module
!! @note   - use module3D_phy
!<

subroutine slope_surf

use module3D_phy
implicit none

real :: inv_4dx        ! inverse de dx pour eviter les divisions =1/(4*dx)
inv_4dx = 1./(4.*dx)



do j=1,ny                                 ! slope along x on node >
   do i=2,nx
      sdx(i,j)=(s(i,j)-s(i-1,j))/dx
   end do
end do

do j=2,ny                                 ! slope along y on node ^
   do i=1,nx
      sdy(i,j)=(s(i,j)-s(i,j-1))/dy
   end do
end do

do j=2,ny-1                               ! slope amplitude on node o
   do i=2,nx-1
      slope(i,j)=(s(i+1,j)-s(i-1,j))**2 + (s(i,j+1)-s(i,j-1))**2
      slope(i,j)=slope(i,j)**0.5 / dx /2.
   end do
end do


do j=2,ny                                 ! slope along x on node ^ (average)
   do i=2,nx-1
      sdxmy(i,j)= ((s(i+1,j)-s(i-1,j))+s(i+1,j-1)-s(i-1,j-1))*inv_4dx
   end do
end do

do j=2,ny-1                               ! slope along y on node > (average) 
   do i=2,nx
      sdymx(i,j)= ((s(i,j+1)-s(i,j-1))+s(i-1,j+1)-s(i-1,j-1))*inv_4dx
   enddo
enddo

slope2mx(:,:)=sdx(:,:)**2+sdymx(:,:)**2   ! slope amplitude on node  >
slope2my(:,:)=sdy(:,:)**2+sdxmy(:,:)**2   ! slope amplitude on node  ^

! conditions on the grid boundaries

slope(1,:)    =  0.
slope(:,1)    =  0.
slope(nx,:)   =  0.
slope(:,ny)   =  0.

slope2mx(1,:) =  0.
slope2mx(:,1) =  0.
slope2mx(:,ny) = 0.

slope2my(1,:) =  0.
slope2my(:,1) =  0.
slope2my(nx,:) = 0.


! limitation to 2e-2
sdx(:,:)=max(sdx(:,:),-2.e-2)
sdx(:,:)=min(sdx(:,:),2.e-2)
sdxmy(:,:)=max(sdxmy(:,:),-2.e-2)
sdxmy(:,:)=min(sdxmy(:,:),2.e-2)
slope2mx(:,:)=min(slope2mx(:,:),4.e-4)

sdy(:,:)=max(sdy(:,:),-2.e-2)
sdy(:,:)=min(sdy(:,:),2.e-2)
sdymx(:,:)=max(sdymx(:,:),-2.e-2)
sdymx(:,:)=min(sdymx(:,:),2.e-2)
slope2my(:,:)=min(slope2my(:,:),4.e-4)

slope(:,:)=min(slope(:,:),2.e-2)


debug_3D(:,:,21)=sdx(:,:)
debug_3D(:,:,22)=sdy(:,:)
end subroutine slope_surf