MODULE recomp_tend_co2_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

!=======================================================================
!
!  Routine that compute the evolution of the tendencie 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(*,*) "Update of the CO2 tendency from the current 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(ngrid,nslope,timelen,d_h2oice,PCM_temp,PEM_temp)

implicit none

!=======================================================================
!
!  Routine that compute the evolution of the tendencie for h2o ice
!
!=======================================================================

! Inputs
! ------
integer,            intent(in) :: timelen, ngrid, nslope
real, dimension(:), intent(in) :: PCM_temp, PEM_temp
! Outputs
! -------
real, dimension(ngrid,nslope), intent(inout) :: d_h2oice ! physical point field: Evolution of perennial ice over one year

! Local:
! ------

write(*,*) "Update of the H2O tendency due to lag layer"

! Flux correction due to lag layer
!~ Rz_old = h2oice_depth_old*0.0325/4.e-4              ! resistance from PCM
!~ Rz_new = h2oice_depth_new*0.0325/4.e-4              ! new resistance based on new depth
!~ R_dec = (1./Rz_old)/(1./Rz_new)                     ! decrease because of resistance
!~ soil_psv_old = psv(max(PCM_temp(h2oice_depth_old))) ! the maxmimum annual mean saturation vapor pressure at the temperature of the GCM run temperature at the old ice location
!~ soil_psv_new = psv(max(PEM_temp(h2oice_depth_new))) ! the maxmimum annual mean saturation vapor pressure at the temperature of the PEM run temperature at the new ice location
!~ hum_dec = soil_psv_old/soil_psv_new                 ! decrease because of lower water vapor pressure at the new depth
!~ d_h2oice = d_h2oice*R_dec*hum_dec                   ! decrease of flux

END SUBROUTINE recomp_tend_h2o

END MODULE recomp_tend_co2_mod