MODULE evol_ice_mod

implicit none

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

SUBROUTINE evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice,zshift_surf)

use time_evol_mod, only: dt
use glaciers_mod,  only: rho_co2ice

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) :: d_co2ice    ! Evolution of perennial ice over one year
real, dimension(ngrid,nslope), intent(out) :: zshift_surf ! Elevation shift for the surface [m]

! 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 + d_co2ice*dt
where (co2_ice < 0.)
    co2_ice = 0.
    d_co2ice = -co2_ice_old/dt
end where

zshift_surf = d_co2ice*dt/rho_co2ice

END SUBROUTINE evol_co2_ice

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

SUBROUTINE evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopPEM)

use time_evol_mod, only: dt
use comslope_mod,  only: subslope_dist, def_slope_mean
use glaciers_mod,  only: rho_h2oice

#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_adsorbed      ! Mass of H2O adsorbed/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) :: d_h2oice    ! Evolution of perennial ice over one year (kg/m^2/year)
integer,                       intent(inout) :: stopPEM     ! Stopping criterion code
real, dimension(ngrid,nslope), intent(out) :: zshift_surf ! Elevation shift for the surface [m]

! local:
! ------
integer                       :: i, islope                                           ! Loop variables
real                          :: pos_tend, neg_tend, real_coefficient, negative_part ! Variable to conserve h2o
real, dimension(ngrid,nslope) :: new_tend                                            ! 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_adsorbed(i) > 0.) then
            pos_tend = pos_tend + delta_h2o_adsorbed(i)*cell_area(i)
        else
            neg_tend = neg_tend + delta_h2o_adsorbed(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 (d_h2oice(i,islope) > 0.) then
                    pos_tend = pos_tend + d_h2oice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
                else
                    neg_tend = neg_tend - d_h2oice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
                endif
            endif
        enddo
    enddo

    new_tend = 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 (d_h2oice < 0.) ! We lower the sublimating rate by a coefficient
                new_tend = d_h2oice*pos_tend/neg_tend
            elsewhere ! We don't change the condensing rate
                new_tend = d_h2oice
            end where
        else if (neg_tend < pos_tend .and. neg_tend > 0.) then ! More condensation on the planet than sublimation
            where (d_h2oice < 0.) ! We don't change the sublimating rate
                new_tend = d_h2oice
            elsewhere ! We lower the condensing rate by a coefficient
                new_tend = d_h2oice*neg_tend/pos_tend
            end where
        endif
    endif

    ! Evolution of the h2o ice for each physical point
    h2o_ice = h2o_ice + new_tend*dt

    ! 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.
                d_h2oice(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_tend(i,islope) > 0.) h2o_ice(i,islope) = h2o_ice(i,islope) - new_tend(i,islope)*real_coefficient*dt*cos(def_slope_mean(islope)*pi/180.)
        enddo
     enddo
else ! ngrid == 1, i.e. in 1D
    h2o_ice = h2o_ice + d_h2oice*dt
endif

zshift_surf = d_h2oice*dt/rho_h2oice

END SUBROUTINE evol_h2o_ice

END MODULE evol_ice_mod