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                                      ! Number of slopes for the statistic  (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, INTENT(IN) :: albedo_slope(ngrid,2,nslope) !  albedo on each sub-grid slope (1)
  REAL, INTENT(IN) :: emis_slope(ngrid,nslope)     !  emissivity on each sub-grid slope (1)
  REAL, INTENT(IN) :: tsurf_slope(ngrid,nslope)    !  Surface Temperature on each sub-grid slope (K)
  REAL, INTENT(IN) :: qsurf_slope(ngrid,nq,nslope) !  Surface Tracers on each sub-grid slope (kg/m^2 sloped)
  REAL, INTENT(OUT) :: albedo_meshavg(ngrid,2)     !  grid box average of the albedo (1)
  REAL, INTENT(OUT) :: emis_meshavg(ngrid)         !  grid box average of the emissivity (1)
  REAL, INTENT(OUT) :: tsurf_meshavg(ngrid)        !  grid box average of the surface temperature (K)
  REAL, INTENT(OUT) :: qsurf_meshavg(ngrid,nq)     !  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.eq.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