MODULE evol_h2o_ice_s_mod

IMPLICIT NONE

CONTAINS

SUBROUTINE evol_h2o_ice_s(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,qsurf,tendencies_h2o_ice_phys,STOPPING)

  use time_evol_mod,          only: dt_pem
  use comslope_mod,           only: subslope_dist, def_slope_mean
  use criterion_pem_stop_mod, only: criterion_waterice_stop
#ifndef CPP_STD
    use comcstfi_h,   only: pi
#else
    use comcstfi_mod, only: pi
#endif

  IMPLICIT NONE

!=======================================================================
!
!  Routine that compute the evolution of the water ice
!
!=======================================================================

!   arguments:
!   ----------

!   INPUT

  INTEGER, intent(in) :: ngrid                                  ! # of grid points along longitude/latitude grid;
  INTEGER, intent(in) :: nslope                                 ! # of subslope
  REAL, intent(in) ::  cell_area(ngrid)                         ! Area of each mesh grid (m^2)
  REAL, intent(in) :: delta_h2o_adsorbded(ngrid)                ! Mass of H2O adsorbded/desorbded in the soil (kg/m^2) 
  REAL, intent(in) :: delta_h2o_icetablesublim(ngrid)                ! Mass of H2O that have condensed/sublimated at the ice table (kg/m^2) 

!   OUTPUT
  REAL, INTENT(INOUT) ::  qsurf(ngrid,nslope)                   ! physical point field : Previous and actual density of water ice (kg/m^2)
  REAL, intent(inout) ::  tendencies_h2o_ice_phys(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year (kg/m^2/year)
  LOGICAL, INTENT(INOUT) :: STOPPING                            ! Stopping criterion

!   local:
!   ----

  INTEGER :: i,j,islope                                         ! loop variable
  REAL :: pos_tend, neg_tend, real_coefficient,negative_part    ! Variable to conserve water
  REAL ::  new_tendencies(ngrid,nslope)                         ! Tendencies computed in order to conserve water ice on the surface, only exchange between surface are done
  
!=======================================================================

  STOPPING=.false.

  pos_tend=0.
  neg_tend=0.
if (ngrid.NE.1) then ! to make sure we are not in 1D 
  ! We compute the amount of water accumulating and sublimating
  do i=1,ngrid
     if(delta_h2o_adsorbded(i).GT.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).GT.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 (qsurf(i,islope).GT.0) then
         if (tendencies_h2o_ice_phys(i,islope).GT.0) then
            pos_tend=pos_tend+tendencies_h2o_ice_phys(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
         else
            neg_tend=neg_tend-tendencies_h2o_ice_phys(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
         endif
     endif
     enddo
  enddo
  ! We adapt the tendencies to conserve water and do only exchange between grid points
   if(neg_tend.GT.pos_tend .and. pos_tend.GT.0) then ! We are sublimating more in the planet than condensing 
     do i=1,ngrid
       do islope=1,nslope
       if(tendencies_h2o_ice_phys(i,islope).LT.0) then ! We lower the sublimating rate by a coefficient
          new_tendencies(i,islope)=tendencies_h2o_ice_phys(i,islope)*(pos_tend/neg_tend)
       else                                            ! We dont't change the accumulating rate
          new_tendencies(i,islope)=tendencies_h2o_ice_phys(i,islope)
       endif
       enddo
     enddo
   elseif(neg_tend.LT.pos_tend .and. neg_tend.GT.0) then ! We are condensing more in the planet than sublimating
     do i=1,ngrid
       do islope=1,nslope
       if(tendencies_h2o_ice_phys(i,islope).LT.0) then ! We dont't change the sublimating rate
          new_tendencies(i,islope)=tendencies_h2o_ice_phys(i,islope)
       else                                            ! We lower the condensing rate by a coefficient
          new_tendencies(i,islope)=tendencies_h2o_ice_phys(i,islope)*(neg_tend/pos_tend)
       endif
       enddo
     enddo
   elseif(pos_tend.EQ.0 .OR. neg_tend.EQ.0) then 
    write(*,*) "Reason of stopping : There is either no water ice sublimating or no water ice increasing !!"
    write(*,*) "Tendencies on ice sublimating=", neg_tend
    write(*,*) "Tendencies on ice increasing=", pos_tend
    write(*,*) "This can be due to the absence of water ice in the PCM run!!"
    call criterion_waterice_stop(cell_area,1.,qsurf(:,:)*0.,STOPPING,ngrid,qsurf(:,:)*0.)
    do i=1,ngrid
       do islope=1,nslope
         new_tendencies(i,islope)=0
       enddo
    enddo
   endif
  negative_part = 0.

! Evolution of the water ice for each physical point
  do i=1,ngrid
    do islope=1, nslope
      qsurf(i,islope)=qsurf(i,islope)+new_tendencies(i,islope)*dt_pem
      ! We compute the amount of water that is sublimated in excess
      if (qsurf(i,islope).lt.0) then
        negative_part=negative_part-qsurf(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
        qsurf(i,islope)=0.
        tendencies_h2o_ice_phys(i,islope)=0.
      endif
    enddo
  enddo
  
  
  if(pos_tend.eq.0) then
   real_coefficient = 0.
  else 
   real_coefficient = negative_part/pos_tend ! 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
  endif
  do i=1,ngrid
    do islope=1, nslope
     if(new_tendencies(i,islope).GT.0) then  ! In the place of accumulation of ice, we remove a bit of ice in order to conserve water
         qsurf(i,islope)=qsurf(i,islope)-new_tendencies(i,islope)*real_coefficient*dt_pem*cos(def_slope_mean(islope)*pi/180.)
     endif
    enddo
  enddo
else ! ngrid==1;
  do islope=1, nslope
    qsurf(1,islope)=qsurf(1,islope)+tendencies_h2o_ice_phys(1,islope)*dt_pem
  enddo
endif

END SUBROUTINE evol_h2o_ice_s

END MODULE evol_h2o_ice_s_mod