MODULE recomp_tend_co2_slope_mod

implicit none

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

SUBROUTINE recomp_tend_co2_slope(ngrid,nslope,timelen,d_co2ice_phys,d_co2ice_phys_ini,co2ice_slope,emissivity_slope, &
                                 vmr_co2_PCM,vmr_co2_pem,ps_PCM_2,global_avg_press_PCM,global_avg_press_new)

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
!
!=======================================================================

!   arguments:
!   ----------
!   INPUT
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_2                    ! physical point field: Surface pressure in the PCM
real,                           intent(in) :: global_avg_press_PCM        ! global averaged pressure at previous timestep
real,                           intent(in) :: global_avg_press_new        ! global averaged pressure at current timestep
real, dimension(ngrid,nslope),  intent(in) :: d_co2ice_phys_ini ! physical point field: Evolution of perennial ice over one year
real, dimension(ngrid,nslope),  intent(in) :: co2ice_slope                ! CO2 ice per mesh and sub-grid slope (kg/m^2)
real, dimension(ngrid,nslope),  intent(in) :: emissivity_slope            ! Emissivity per mesh and sub-grid slope(1)
!   OUTPUT
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, ave

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_slope(i,islope)*sigmaB/Lco2
        ave = 0.
        if (co2ice_slope(i,islope) > 1.e-4 .and. abs(d_co2ice_phys(i,islope)) > 1.e-5) then
            do t=1,timelen
                ave = ave + (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PCM(i,t)*ps_PCM_2(i,t)/100.)))**4 &
                      - (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem(i,t)*ps_PCM_2(i,t)*(global_avg_press_new/global_avg_press_PCM)/100.)))**4
            enddo
            if (ave < 1e-4) ave = 0.
            d_co2ice_phys(i,islope) = d_co2ice_phys_ini(i,islope) - coef*ave/timelen
        endif
    enddo
enddo

END SUBROUTINE recomp_tend_co2_slope

END MODULE recomp_tend_co2_slope_mod