MODULE recomp_tend_mod

implicit none

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

SUBROUTINE recomp_tend_co2(ngrid,nslope,timelen,d_co2ice_phys,d_co2ice_ini,co2ice,emissivity, &
                           vmr_co2_PCM,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global)

use constants_marspem_mod, only : alpha_clap_co2, beta_clap_co2, sigmaB, Lco2, sols_per_my, sec_per_sol

implicit none

!=======================================================================
!
!  To compute the evolution of the tendency for co2 ice
!
!=======================================================================

! Inputs
! ------
integer,                        intent(in) :: timelen, ngrid, nslope
real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PCM       ! physical point field: Volume mixing ratio of co2 in the first layer
real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM       ! physical point field: Volume mixing ratio of co2 in the first layer
real, dimension(ngrid,timelen), intent(in) :: ps_PCM            ! physical point field: Surface pressure in the PCM
real,                           intent(in) :: ps_avg_global_ini ! global averaged pressure at previous timestep
real,                           intent(in) :: ps_avg_global     ! global averaged pressure at current timestep
real, dimension(ngrid,nslope),  intent(in) :: d_co2ice_ini      ! physical point field: Evolution of perennial ice over one year
real, dimension(ngrid,nslope),  intent(in) :: co2ice            ! CO2 ice per mesh and sub-grid slope (kg/m^2)
real, dimension(ngrid,nslope),  intent(in) :: emissivity        ! Emissivity per mesh and sub-grid slope(1)
! Outputs
! -------
real, dimension(ngrid,nslope), intent(inout) :: d_co2ice_phys ! physical point field: Evolution of perennial ice over one year

! Local:
! ------
integer :: i, t, islope
real    :: coef, avg

write(*,*) "> Updating the CO2 ice tendency for the new pressure"

! Evolution of the water ice for each physical point
do i = 1,ngrid
    do islope = 1,nslope
        coef = sols_per_my*sec_per_sol*emissivity(i,islope)*sigmaB/Lco2
        avg = 0.
        if (co2ice(i,islope) > 1.e-4 .and. abs(d_co2ice_phys(i,islope)) > 1.e-5) then
            do t = 1,timelen
                avg = avg + (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PCM(i,t)*ps_PCM(i,t)/100.)))**4 &
                      - (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PEM(i,t)*ps_PCM(i,t)*(ps_avg_global/ps_avg_global_ini)/100.)))**4
            enddo
            if (avg < 1.e-4) avg = 0.
            d_co2ice_phys(i,islope) = d_co2ice_ini(i,islope) - coef*avg/timelen
        endif
    enddo
enddo

END SUBROUTINE recomp_tend_co2
!=======================================================================

SUBROUTINE recomp_tend_h2o(icetable_depth_old,icetable_depth_new,tsurf,tsoil_PEM_timeseries_old,tsoil_PEM_timeseries_new,d_h2oice)

use compute_soiltemp_mod, only: itp_tsoil
use fast_subs_mars,       only: psv

implicit none

!=======================================================================
!
!  To compute the evolution of the tendency for h2o ice
!
!=======================================================================

! Inputs
! ------
real,                 intent(in) :: icetable_depth_old, icetable_depth_new, tsurf
real, dimension(:,:), intent(in) :: tsoil_PEM_timeseries_old, tsoil_PEM_timeseries_new
! Outputs
! -------
real, intent(inout) :: d_h2oice ! Evolution of perennial ice over one year

! Local:
! ------
real            :: Rz_old, Rz_new, R_dec, hum_dec, psv_max_old, psv_max_new
integer         :: t
real, parameter :: coef_diff = 4.e-4 ! Diffusion coefficient
real, parameter :: zcdv = 0.0325     ! Drag coefficient

write(*,*) "> Updating the H2O tendency due to lag layer"

! Higher resistance due to growing lag layer (higher depth)
Rz_old = icetable_depth_old*zcdv/coef_diff ! Old resistance from PCM
Rz_new = icetable_depth_new*zcdv/coef_diff ! New resistance based on new depth
R_dec = Rz_new/Rz_old ! Decrease because of resistance

! The maxmimum of the daily averages over one year for the saturation vapor pressure at the ice table location
psv_max_old = 0.
psv_max_new = 0.
do t = 1,size(tsoil_PEM_timeseries_old,2)
    psv_max_old = max(psv_max_old,psv(itp_tsoil(tsoil_PEM_timeseries_old(:,t),tsurf,icetable_depth_old)))
    psv_max_new = max(psv_max_new,psv(itp_tsoil(tsoil_PEM_timeseries_new(:,t),tsurf,icetable_depth_new)))
enddo

! Lower humidity due to growing lag layer (higher depth)
hum_dec = psv_max_old/psv_max_new ! Decrease because of lower water vapor pressure at the new depth

! Flux correction (decrease)
d_h2oice = d_h2oice*R_dec*hum_dec

END SUBROUTINE recomp_tend_h2o

END MODULE recomp_tend_mod