MODULE slope_mod

implicit none

real, save, dimension(:), allocatable :: theta_sl ! slope angle versus horizontal (deg)
real, save, dimension(:), allocatable :: psi_sl   ! slope orientation (deg)

!$OMP THREADPRIVATE(theta_sl,psi_sl)

!=======================================================================
contains
!=======================================================================

SUBROUTINE getslopes(ngrid,geopot)

use geometry_mod,       only: longitude, latitude ! in radians
use comcstfi_h,         only: g, rad, pi
use mod_phys_lmdz_para, only: is_parallel
use mod_grid_phy_lmdz,  only: nbp_lon, nbp_lat

implicit none

! 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

! Input arguments
integer,                intent(in) :: ngrid  ! nnumber of atmospheric columns
real, dimension(ngrid), intent(in) :: geopot ! geopotential on phy grid

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

if (is_parallel) then
    ! This routine only works in serial mode so stop now.
    write(*,*) "getslopes Error: this routine is not designed to run in parallel"
    call abort_physic("getslopes",'cannot be run in parallel',1)
endif

id2  = 2
idm1 = nbp_lon-1

! rearrange topography on a 2d array
do j = 2,nbp_lat-1
    ig0 = 1 + (j - 2)*nbp_lon
    do i = 1,nbp_lon
        topogrid(i,j) = geopot(ig0 + i)/g
        latigrid(i,j) = latitude(ig0 + i)
        longgrid(i,j) = longitude(ig0 + i)
    enddo
enddo

! poles:
topogrid(:,1) = geopot(1)/g
latigrid(:,1) = latitude(1)
longgrid(:,1) = longitude(1)
topogrid(:,nbp_lat) = geopot(ngrid)/g
latigrid(:,nbp_lat) = latitude(ngrid)
longgrid(:,nbp_lat) = longitude(ngrid)

! compute topography gradient
! topogrid and rad are both in meters
do j = 2,nbp_lat - 1
    do i=1,nbp_lon
        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(nbp_lon,j))/(2*pi + longgrid(id2,j) - longgrid(nbp_lon,j))
    grady(1,j) = grady(1,j) / rad
    grady(nbp_lon,j) = (topogrid(1,j) - topogrid(idm1,j))/(2*pi + longgrid(1,j) - longgrid(idm1,j))
    grady(nbp_lon,j) = grady(nbp_lon,j)/rad
    do i = 2,nbp_lon - 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(:,nbp_lat) = 0.
grady(:,nbp_lat) = 0.

! compute slope inclination and orientation:
theta_sl = 0.
psi_sl  = 0.
do j = 2,nbp_lat - 1
    do i = 1,nbp_lon
        ig0 = 1 + (j - 2)*nbp_lon

        theta_val = atan(sqrt((gradx(i,j))**2 + (grady(i,j))**2))

        psi_val = 0.
        if (gradx(i,j) /= 0.) psi_val = -pi/2. - atan(grady(i,j)/gradx(i,j))
        if (gradx(i,j) >= 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

!=======================================================================

SUBROUTINE param_slope(csza,declin,rho,latitude,taudust,albedo,theta_s,psi_s,fdir_0,ftot_0,ftot)
!***********************************************************************
!
! SUBROUTINE:
! param_slope
! 
!
! PURPOSE: 
! computes total solar irradiance on a given Martian slope 
!
!
! INPUTS:
! csza		cosine solar zenith angle
! declin	sun declination (rad)
! rho		sun right ascension (rad)
! latitude	latitude (deg)
! taudust	dust optical depth at reference wavelength 0.67 mic. 
! albedo	spectrally integrated surface Lambertian reflection albedo 
! theta_s	slope inclination angle (deg)
!			0 is horizontal, 90 is vertical
! phi_s		slope azimuth (deg) 
!       		0 >> Northward
!       		90 >> Eastward
!       		180 >> Southward
!       		270 >> Westward
! ftot_0	spectrally integrated total irradiance on an horizontal surface (W/m2)
! fdir_0	spectrally integrated direct irradiance on an horizontal surface (W/m2)
!
!
! OUTPUTS:
! ftot		spectrally integrated total irradiance on the slope (W/m2)
!
! REFERENCE: 
! "Fast and accurate estimation of irradiance on Martian slopes"
! A. Spiga & F. Forget
! .....
!
! AUTHOR: 
! A. Spiga (spiga@lmd.jussieu.fr)
! March 2008
!
!***********************************************************************

use comcstfi_h, only: pi

implicit none

! Input arguments
real, intent(in) :: csza, declin, rho, latitude, taudust, theta_s, psi_s, albedo, ftot_0 , fdir_0

! Output arguments
real, intent(out) :: ftot

! Local variables
real                 :: deg2rad, a, mu_s, sigma_s, fdir, fscat, fscat_0, fref, ratio
real, dimension(4,2) :: mat_M, mat_N, mat_T
real, dimension(2)   :: g_vector
real, dimension(4)   :: s_vector
!***********************************************************************
! Prerequisite
deg2rad = pi/180.
if ((theta_s > 90.) .or. (theta_s < 0.)) then
    write(*,*) 'please set theta_s between 0 and 90', theta_s
    call abort_physic("param_slopes","invalid theta_s",1) 
endif

! Solar Zenith angle (radian)
if (csza < 0.01) then
    !print *, 'sun below horizon' 
    !fdir_0=0.
    fdir    = 0.
    fscat_0 = 0.
    fscat   = 0.   
    fref    = 0.
else
! Low incidence fix 
!    if (csza < 0.15) csza = 0.15

! 'Slope vs Sun' azimuth (radian)
    if (cos(declin)*sin(rho) == 0. .and. &
        sin(deg2rad*latitude)*cos(declin)*cos(rho) - cos(deg2rad*latitude)*sin(declin) == 0.) then
        a = deg2rad*psi_s ! some compilator need specfying value for atan2(0,0)  
    else
        a = deg2rad*psi_s + atan2(cos(declin)*sin(rho),sin(deg2rad*latitude)*cos(declin)*cos(rho)-cos(deg2rad*latitude)*sin(declin))
    endif

! Cosine of slope-sun phase angle 
    mu_s = csza*cos(deg2rad*theta_s) - cos(a)*sin(deg2rad*theta_s)*sqrt(1-csza**2)
    if (mu_s <= 0.) mu_s = 0.

! Sky-view factor
    sigma_s=0.5*(1. + cos(deg2rad*theta_s))

! Direct flux on the slope
    fdir = fdir_0*mu_s/csza

! Reflected flux on the slope
    fref = albedo*(1 - sigma_s)*ftot_0

! Scattered flux on a flat surface
    fscat_0 = ftot_0 - fdir_0

! Scattering vector (slope vs sky)
    s_vector = (/ 1., exp(-taudust), sin(deg2rad*theta_s), sin(deg2rad*theta_s)*exp(-taudust) /)

! Geometry vector (slope vs sun)
    g_vector = (/ mu_s/csza, 1. /)

! Coupling matrix
    if (csza >= 0.5) then
        mat_M(:,1) = (/ -0.264,  1.309,  0.208, -0.828 /)
        mat_M(:,2) = (/  1.291*sigma_s, -1.371*sigma_s, -0.581,  1.641 /)
        mat_N(:,1) = (/  0.911, -0.777, -0.223,  0.623 /)
        mat_N(:,2) = (/ -0.933*sigma_s,  0.822*sigma_s,  0.514, -1.195 /)
    else
        mat_M(:,1) = (/ -0.373,  0.792, -0.095,  0.398 /)
        mat_M(:,2) = (/  1.389*sigma_s, -0.794*sigma_s, -0.325,  0.183 /)
        mat_N(:,1) = (/  1.079,  0.275,  0.419, -1.855 /)
        mat_N(:,2) = (/ -1.076*sigma_s, -0.357*sigma_s, -0.075,  1.844 /)
    endif
    mat_T = mat_M + csza*mat_N

! Scattered flux slope ratio
    if (deg2rad*theta_s <= 0.0872664626) then
    ! low angles
        s_vector = (/ 1., exp(-taudust) , sin(0.0872664626), sin(0.0872664626)*exp(-taudust) /)
        ratio = dot_product(matmul(s_vector, mat_T),g_vector)
        ratio = 1. + (ratio - 1.)*deg2rad*theta_s/0.0872664626
    else
    ! general case
    ratio = dot_product(matmul(s_vector,mat_T),g_vector)  
    ! NB: ratio = dot_product(s_vector,matmul(mat_T,g_vector)) is equivalent
    endif

! Scattered flux on the slope
    fscat = ratio*fscat_0
endif ! if (csza < 0.01)

! Total flux on the slope
ftot = fdir + fref + fscat

! Display results
!  print *, 'sca component 0 ', ftot_0-fdir_0
!  print *, 'dir component 0 ', fdir_0
!  print *, 'scattered component ', fscat
!  print *, 'direct component ', fdir
!  print *, 'reflected component ', fref

END SUBROUTINE param_slope

!=======================================================================

SUBROUTINE ini_slope_mod(ngrid)

implicit none

integer, intent(in) :: ngrid ! number of atmospheric columns

allocate(theta_sl(ngrid))
allocate(psi_sl(ngrid))

END SUBROUTINE ini_slope_mod

!=======================================================================

SUBROUTINE end_slope_mod

implicit none

if (allocated(theta_sl)) deallocate(theta_sl)
if (allocated(psi_sl)) deallocate(psi_sl)

END SUBROUTINE end_slope_mod

END MODULE slope_mod