source: trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod_slope.F90 @ 2945

!
! $Id $
!
SUBROUTINE evol_h2o_ice_s_slope(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)

!   OUTPUT
  REAL, INTENT(INOUT) ::  qsurf(ngrid,nslope)                ! physical point field : Previous and actual density of water ice
  LOGICAL, INTENT(INOUT) :: STOPPING
  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
  REAL ::  new_tendencies(ngrid,nslope)

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

  STOPPING=.false.

  pos_tend=0.
  neg_tend=0.

  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

  if(neg_tend.GT.pos_tend .and. pos_tend.GT.0) then
     do i=1,ngrid
       do islope=1,nslope
       if(tendencies_h2o_ice_phys(i,islope).LT.0) then
          new_tendencies(i,islope)=tendencies_h2o_ice_phys(i,islope)*(pos_tend/neg_tend)
       else
          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
     do i=1,ngrid
       do islope=1,nslope
       if(tendencies_h2o_ice_phys(i,islope).LT.0) then
          new_tendencies(i,islope)=tendencies_h2o_ice_phys(i,islope)
       else
          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
      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
  endif
  do i=1,ngrid
    do islope=1, nslope
     if(new_tendencies(i,islope).GT.0) then
         qsurf(i,islope)=qsurf(i,islope)-new_tendencies(i,islope)*real_coefficient*dt_pem
     endif
    enddo
  enddo



END SUBROUTINE evol_h2o_ice_s_slope