source: trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod.F90 @ 3136

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, 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) :: qsurf                   ! physical point field: Previous and actual density of water ice (kg/m^2)
real, dimension(ngrid,nslope), intent(inout) :: tendencies_h2o_ice_phys ! physical point field: Evolution of perennial 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