source: branches/GRISLIv3/SOURCES/slope_surf.f90 @ 483

!> \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, only: sdx,sdy,S,slope,slope2mx,slope2my,debug_3D
  use geography, only: nx,ny,dx,dy
  !$ USE OMP_LIB
  
  implicit none

  real,dimension(nx,ny) :: SDXMY        !< slope selon x moy selon y '^' remplace SDMX
  real,dimension(nx,ny) :: SDYMX        !< slope selon y moy selon x '>' remplace SDMY

  real :: inv_4dx        ! inverse de dx pour eviter les divisions =1/(4*dx)
  integer :: i,j

  inv_4dx = 1./(4.*dx)


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

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

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

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

  !$OMP WORKSHARE
  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(:,:)
  !$OMP END WORKSHARE
  !$OMP END PARALLEL
  
end subroutine slope_surf