MODULE comslope_mod
!======================================================================================================================!
! Subject:
!---------
!   Module used for the parametrization of subgrid scale slope
!----------------------------------------------------------------------------------------------------------------------!
! Reference:
!-----------
!   Lange et al. (2023 in prep.), 'Introduction of Sub-Grid Scale Slope Microclimates in the Mars Planetary Climate Model', JGR
!
!======================================================================================================================!

implicit none

integer, save                              :: nslope         ! Number of slopes for the statistic (1)
integer, save                              :: iflat          ! Flat slope for islope (1)
real,    save, allocatable, dimension(:)   :: def_slope      ! Bound of the slope statistics / repartitions (°)
real,    save, allocatable, dimension(:)   :: def_slope_mean ! Mean slope of each bin (°)
real,    save, allocatable, dimension(:)   :: sky_slope      ! Portion of the sky view by each slopes (1)
real,    save, allocatable, dimension(:,:) :: subslope_dist  ! Distribution of the slopes (1)
integer, save, allocatable, dimension(:)   :: major_slope    ! Index of the subslope that occupies most of themesh (1)
!$OMP THREADPRIVATE(nslope,def_slope,def_slope_mean,sky_slope,subslope_dist,iflat,major_slope)

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

SUBROUTINE ini_comslope_h(ngrid,nslope_in)

implicit none

integer, intent(in) :: ngrid
integer, intent(in) :: nslope_in

allocate(def_slope(nslope_in+1))
allocate(def_slope_mean(nslope_in))
allocate(sky_slope(nslope_in))
allocate(subslope_dist(ngrid,nslope_in))
allocate(major_slope(ngrid))

END SUBROUTINE ini_comslope_h

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

SUBROUTINE end_comslope_h

implicit none

if (allocated(def_slope)) deallocate(def_slope)
if (allocated(def_slope_mean)) deallocate(def_slope_mean)
if (allocated(sky_slope)) deallocate(sky_slope)
if (allocated(subslope_dist)) deallocate(subslope_dist)
if (allocated(major_slope)) deallocate(major_slope)

END SUBROUTINE end_comslope_h

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

SUBROUTINE compute_meshgridavg(ngrid,nq,albedo_slope,emis_slope,tsurf_slope,qsurf_slope,albedo_meshavg,emis_meshavg,tsurf_meshavg, qsurf_meshavg)

use comcstfi_h, only: pi

implicit none

integer, intent(in) :: ngrid, nq                             ! # points on the physical grid, tracers (1)
real, dimension(ngrid,2,nslope),  intent(in) :: albedo_slope ! albedo on each sub-grid slope (1)
real, dimension(ngrid,nslope),    intent(in) :: emis_slope   ! emissivity on each sub-grid slope (1)
real, dimension(ngrid,nslope),    intent(in) :: tsurf_slope  ! Surface Temperature on each sub-grid slope (K)
real, dimension(ngrid,nq,nslope), intent(in) :: qsurf_slope  ! Surface Tracers on each sub-grid slope (kg/m^2 sloped)
real, dimension(ngrid,2),  intent(out) :: albedo_meshavg ! grid box average of the albedo (1)
real, dimension(ngrid),    intent(out) :: emis_meshavg   ! grid box average of the emissivity (1)
real, dimension(ngrid),    intent(out) :: tsurf_meshavg  ! grid box average of the surface temperature (K)
real, dimension(ngrid,nq), intent(out) :: qsurf_meshavg  ! grid box average of the surface tracers (kg/m^2 flat)
! Local
integer :: islope, ig, l, iq

!-------------------

albedo_meshavg = 0.
emis_meshavg = 0.
tsurf_meshavg = 0.
qsurf_meshavg = 0.

if (nslope == 1) then
    albedo_meshavg = albedo_slope(:,:,1)
    emis_meshavg = emis_slope(:,1)
    tsurf_meshavg = tsurf_slope(:,1)
    qsurf_meshavg = qsurf_slope(:,:,1)
else
    do ig = 1,ngrid
        do islope = 1,nslope
            do l = 1,2
                albedo_meshavg(ig,l) = albedo_meshavg(ig,l) + albedo_slope(ig,l,islope)*subslope_dist(ig,islope)
            enddo
            do iq = 1,nq
                qsurf_meshavg(ig,iq) = qsurf_meshavg(ig,iq) + qsurf_slope(ig,iq,islope)*subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
            enddo
            emis_meshavg(ig) = emis_meshavg(ig) + emis_slope(ig,islope)*subslope_dist(ig,islope)
            tsurf_meshavg(ig) = tsurf_meshavg(ig) + (emis_slope(ig,islope)*tsurf_slope(ig,islope)**4)*subslope_dist(ig,islope)
        enddo
        tsurf_meshavg(ig) = (tsurf_meshavg(ig)/emis_meshavg(ig))**(0.25)
    enddo
endif

END SUBROUTINE compute_meshgridavg

END MODULE comslope_mod