subroutine getslopes(geopot)
    
implicit none

#include "dimensions.h"
#include "dimphys.h"
#include "slope.h"
#include "comgeomfi.h"
#include "comcstfi.h"


! This routine computes slope inclination and orientation for the GCM (callslope=.true. in callphys.def)
! It works fine with a non-regular grid for zoomed simulations.
! slope inclination angle (deg)  0 == horizontal, 90 == vertical
! slope orientation angle (deg)  0 == Northward,  90 == Eastward, 180 == Southward, 270 == Westward
! TN 04/1013


real geopot(ngridmx)     ! geopotential on phy grid
real topogrid(iim,jjm+1) ! topography on lat/lon grid with poles and only one -180/180 point
real latigrid(iim,jjm+1),longgrid(iim,jjm+1) ! meshgrid of latitude and longitude values (radians)
real theta_val ! slope inclination
real psi_val   ! slope orientation
real gradx(iim,jjm+1) ! x: latitude-wise topography gradient,  increasing northward
real grady(iim,jjm+1) ! y: longitude-wise topography gradient, increasing westward
integer i,j,ig0
integer id2,idm1 ! a trick to compile testphys1d with debug option

id2  = 2
idm1 = iim-1

! rearrange topography on a 2d array
do j=2,jjm
   ig0= 1+(j-2)*iim
   do i=1,iim
      topogrid(i,j)=geopot(ig0+i)/g
      latigrid(i,j)=lati(ig0+i)
      longgrid(i,j)=long(ig0+i)
   enddo
enddo
!poles :
topogrid(:,1) = geopot(1)/g
latigrid(:,1) = lati(1)
longgrid(:,1) = long(1)
topogrid(:,jjm+1) = geopot(ngridmx)/g
latigrid(:,jjm+1) = lati(ngridmx)
longgrid(:,jjm+1) = long(ngridmx)



! compute topography gradient
! topogrid and rad are both in meters
do j=2,jjm
   do i=1,iim
     gradx(i,j) = (topogrid(i,j+1) - topogrid(i,j-1)) / (latigrid(i,j+1)-latigrid(i,j-1))
     gradx(i,j) = gradx(i,j) / rad
   enddo
   grady(1,j) = (topogrid(id2,j) - topogrid(iim,j)) / (2*pi+longgrid(id2,j)-longgrid(iim,j)) 
   grady(1,j) = grady(1,j) / rad
   grady(iim,j) = (topogrid(1,j) - topogrid(idm1,j)) / (2*pi+longgrid(1,j)-longgrid(idm1,j)) 
   grady(iim,j) = grady(iim,j) / rad
   do i=2,iim-1
     grady(i,j) = (topogrid(i+1,j) - topogrid(i-1,j)) / (longgrid(i+1,j)-longgrid(i-1,j)) 
     grady(i,j) = grady(i,j) / rad
   enddo
enddo
! poles :
gradx(:,1) = 0.
grady(:,1) = 0.
gradx(:,jjm+1) = 0.
grady(:,jjm+1) = 0.



! compute slope inclination and orientation :
theta_sl(:) = 0.
psi_sl(:)   = 0.
do j=2,jjm
   do i=1,iim
   
     ig0= 1+(j-2)*iim
   
     theta_val=atan(sqrt( (gradx(i,j))**2 + (grady(i,j))**2 ))
     
     psi_val=0.
     if (gradx(i,j) .ne. 0.) psi_val= -pi/2. - atan(grady(i,j)/gradx(i,j))
     if (gradx(i,j) .ge. 0.) psi_val= psi_val - pi
     psi_val = 3*pi/2. - psi_val
     psi_val = psi_val*180./pi
     psi_val = MODULO(psi_val,360.)
     
     theta_sl(ig0+i) = theta_val
     psi_sl(ig0+i)   = psi_val
          
   enddo
enddo



end subroutine getslopes