MODULE evol_h2o_ice_s_mod

IMPLICIT NONE

CONTAINS

SUBROUTINE evol_h2o_ice_s(qsurf,tendencies_h2o_ice_phys,&
                             iim_input,jjm_input,ngrid,cell_area,STOPPING,nslope)

  USE temps_mod_evol, ONLY: dt_pem
  use comslope_mod, ONLY: subslope_dist
  use criterion_pem_stop_mod, only: criterion_waterice_stop

  IMPLICIT NONE

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

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

!   INPUT

  INTEGER, intent(in) :: iim_input,jjm_input, ngrid,nslope      ! # of grid points along longitude/latitude/ total
  REAL, intent(in) ::  cell_area(ngrid)                         ! Area of each mesh grid

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

!   local:
!   ----

  INTEGER :: i,j,ig0,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
     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)
         else
            neg_tend=neg_tend-tendencies_h2o_ice_phys(i,islope)*cell_area(i)*subslope_dist(i,islope)
         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 
    print *, "Reason of stopping : There is either no water ice sublimating or no water ice increasing !!"
    print *, "Tendencies on ice sublimating=", neg_tend
    print *, "Tendencies on ice increasing=", pos_tend
    print *, "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)
        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
     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