MODULE evol_ice_mod

implicit none

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

SUBROUTINE evol_co2_ice(ngrid,nslope,co2_ice,tend_co2_ice)

use time_evol_mod, only: dt_pem

implicit none

!=======================================================================
!
! Routine to compute the evolution of CO2 ice
!
!=======================================================================
!   arguments:
!   ----------
!   INPUT
integer, intent(in) :: ngrid, nslope ! # of grid points, # of subslopes

!   OUTPUT
real, dimension(ngrid,nslope), intent(inout) :: co2_ice      ! Previous and actual density of CO2 ice
real, dimension(ngrid,nslope), intent(inout) :: tend_co2_ice ! Evolution of perennial ice over one year

!   local:
!   ------
real, dimension(ngrid,nslope) :: co2_ice_old ! Old density of CO2 ice
!=======================================================================
! Evolution of CO2 ice for each physical point
write(*,*) 'Evolution of co2 ice'

co2_ice_old = co2_ice
co2_ice = co2_ice + tend_co2_ice*dt_pem
where (co2_ice < 0.)
    co2_ice = 0.
    tend_co2_ice = -co2_ice_old/dt_pem
end where

END SUBROUTINE evol_co2_ice

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

SUBROUTINE evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,tend_h2o_ice,stopPEM)

use time_evol_mod,          only: dt_pem
use comslope_mod,           only: subslope_dist, def_slope_mean

#ifndef CPP_STD
    use comcstfi_h,   only: pi
#else
    use comcstfi_mod, only: pi
#endif

implicit none

!=======================================================================
!
! Routine to compute the evolution of h2o ice
!
!=======================================================================
!   arguments:
!   ----------
!   INPUT
integer,                intent(in) :: ngrid, nslope            ! # of grid points, # of subslopes
real, dimension(ngrid), intent(in) :: cell_area                ! Area of each mesh grid (m^2)
real, dimension(ngrid), intent(in) :: delta_h2o_adsorbded      ! Mass of H2O adsorbded/desorbded in the soil (kg/m^2)
real, dimension(ngrid), intent(in) :: delta_h2o_icetablesublim ! Mass of H2O that have condensed/sublimated at the ice table (kg/m^2)

!   OUTPUT
real, dimension(ngrid,nslope), intent(inout) :: h2o_ice      ! Previous and actual density of h2o ice (kg/m^2)
real, dimension(ngrid,nslope), intent(inout) :: tend_h2o_ice ! Evolution of perennial ice over one year (kg/m^2/year)
integer,                       intent(inout) :: stopPEM      ! Stopping criterion code

!   local:
!   ------
integer                       :: i, islope                                           ! Loop variables
real                          :: pos_tend, neg_tend, real_coefficient, negative_part ! Variable to conserve h2o
real, dimension(ngrid,nslope) :: new_tendencies                                      ! Tendencies computed in order to conserve h2o ice on the surface, only exchange between surface are done
!=======================================================================
if (ngrid /= 1) then ! Not in 1D
    ! We compute the amount of condensing and sublimating h2o ice
    pos_tend = 0.
    neg_tend = 0.
    do i = 1,ngrid
        if (delta_h2o_adsorbded(i) > 0.) then
            pos_tend = pos_tend + delta_h2o_adsorbded(i)*cell_area(i)
        else
            neg_tend = neg_tend + delta_h2o_adsorbded(i)*cell_area(i)
        endif
        if (delta_h2o_icetablesublim(i) > 0.) then
            pos_tend = pos_tend + delta_h2o_icetablesublim(i)*cell_area(i)
        else
            neg_tend = neg_tend + delta_h2o_icetablesublim(i)*cell_area(i)
        endif
        do islope = 1,nslope
            if (h2o_ice(i,islope) > 0.) then
                if (tend_h2o_ice(i,islope) > 0.) then
                    pos_tend = pos_tend + tend_h2o_ice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
                else
                    neg_tend = neg_tend - tend_h2o_ice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
                endif
            endif
        enddo
    enddo

    new_tendencies = 0.
    if (abs(pos_tend) < 1.e-10 .or. abs(neg_tend) < 1.e-10) then
        write(*,*) "Reason of stopping: there is no sublimating or condensing h2o ice!"
        write(*,*) "Tendencies on ice sublimating =", neg_tend
        write(*,*) "Tendencies on ice increasing =", pos_tend
        write(*,*) "This can be due to the absence of h2o ice in the PCM run!"
        stopPEM = 2
    else
        ! We adapt the tendencies to conserve h2o and do only exchange between grid points
        if (neg_tend > pos_tend .and. pos_tend > 0.) then ! More sublimation on the planet than condensation
            where (tend_h2o_ice < 0.) ! We lower the sublimating rate by a coefficient
                new_tendencies = tend_h2o_ice*pos_tend/neg_tend
            elsewhere ! We don't change the condensing rate
                new_tendencies = tend_h2o_ice
            end where
        else if (neg_tend < pos_tend .and. neg_tend > 0.) then ! More condensation on the planet than sublimation
            where (tend_h2o_ice < 0.) ! We don't change the sublimating rate
                new_tendencies = tend_h2o_ice
            elsewhere ! We lower the condensing rate by a coefficient
                new_tendencies = tend_h2o_ice*neg_tend/pos_tend
            end where
        endif
    endif

    ! Evolution of the h2o ice for each physical point
    h2o_ice = h2o_ice + new_tendencies*dt_pem

    ! We compute the amount of h2o that is sublimated in excess
    negative_part = 0.
    do i = 1,ngrid
        do islope = 1,nslope
            if (h2o_ice(i,islope) < 0.) then
                negative_part = negative_part - h2o_ice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
                h2o_ice(i,islope) = 0.
                tend_h2o_ice(i,islope) = 0.
            endif
        enddo
    enddo

    ! We compute a coefficient by which we should remove the ice that has been added to places
    ! even if this ice was contributing to an unphysical negative amount of ice at other places
    if (abs(pos_tend) < 1.e-10) then
        real_coefficient = 0.
    else
        real_coefficient = negative_part/pos_tend
    endif
    ! In the place of accumulation of ice, we remove a bit of ice in order to conserve h2o
    do islope = 1,nslope
        do i = 1,ngrid
             if (new_tendencies(i,islope) > 0.) h2o_ice(i,islope) = h2o_ice(i,islope) - new_tendencies(i,islope)*real_coefficient*dt_pem*cos(def_slope_mean(islope)*pi/180.)
        enddo
     enddo
else ! ngrid == 1, i.e. in 1D
    h2o_ice = h2o_ice + tend_h2o_ice*dt_pem
endif

END SUBROUTINE evol_h2o_ice

END MODULE evol_ice_mod