!> \file surface-0.3.f90
!! Calcule de l'altitude de la surface tous les DTT (pas de temps)
!<

!> SUBROUTINE:  SURFACE()
!! \author Christophe Dumas
!! \date Mai 2000
!! @note  Routine qui calcul l'altitude de la surface tous les DTT (pas de temps)
!! @note  Iterations jusqu'a stabilisation de la surface.
!! @note  Les changements d'altitude de la surface sont limites a 1 m a chaque
!! iteration.
!! @note Used modules:
!! @note    - use module3D_phy
!! @note    - use RELAXATION_MOD
!!
!<


!  *****************************************************************
!
!
!
!                Christophe Dumas, Mai 2000
! 
!   routine qui calcul l'altitude de la surface tous les DTT (pas de temps)
!   iterations jusqu'a stabilisation de la surface.
!   les changements d'altitude de la surface sont limites a 1 m a chaque
!   iteration.
!
!  *****************************************************************

subroutine SURFACE()

       USE module3D_phy

  USE RELAXATION_MOD ! module contenant la routine relaxation
                     ! avec interface explicite

  IMPLICIT NONE

  INTEGER :: ntour  
  INTEGER :: compteur_limitation !< compteur du nbr de points ou on a limite 
  INTEGER :: NXlocal,NYlocal !< pour pouvoir passer NX et NY en parametre a la routine relaxation
  REAL :: dtsdx
  REAL :: difH
  REAL,dimension(NX,NY) :: vieuxH,vieuxH1 !< ancienne valeur de H
  REAL :: SURNET !< permet le calcul de l'altitude de la surface qd la glace flotte

! initialisation des variables
!------------------------------
  NXlocal=NX
  NYlocal=NY
  ntour=0
  DT= DTT
  dtsdx= dt/dx
  compteur_limitation= 51

! Copie de H sur vieuxH
! --------------------


  vieuxH(:,:)= H(:,:)
  vieuxH1(:,:)=H(:,:)
  write(6,*) H(70,60)
! resolution de l'equation de relaxation avec le pas de temps DTT


!  DO WHILE ((compteur_limitation.GT.50))
  DO WHILE ((compteur_limitation.GT.50).AND.(ntour.LT.20))

  compteur_limitation = 0
  ntour = ntour + 1


  CALL RELAXATION(NXlocal,NYlocal,MK,MARINE,FLOT,DTT,DX,vieuxH,BM,BMELT,UXBAR,UYBAR,H,HDOT)



! pour empecher de partir vers des valeurs trop elevees on limite a 1 m les variations
  write(6,*) H(70,60)
  do i=1,NX
     do j=1,NY
        difH= H(i,j) - vieuxH1(i,j)

        IF ((ABS(difH)).GT.1.) THEN
           compteur_limitation = compteur_limitation + 1
           IF (difH.GT.0.) THEN
              H(I,J)= vieuxH1(I,J) + 1.
           ELSE
              H(I,J)= vieuxH1(I,J) - 1.
           ENDIF
        ENDIF
     enddo
  enddo
  write(6,*) H(70,60)
  H(:,:)=MAX(H(:,:),1.)
  vieuxH1(:,:)=H(:,:)


!*************************************
! passage par diffusiv et Neffect pour avoir les vitesses diagnostiques
  call Neffect()
  call DIFFUSIV()


  do i=1,nx
     do j=1,ny
        IF (.NOT.FLOT(I,J)) then
           S(I,J)=BSOC(I,J)+H(I,J)
        ELSE
           SURNET = H(I,J) * (1.-RO/ROW)
           S(I,J) = SURNET + SEALEVEL
           B(I,J) = S(I,J) - H(I,J)
        END IF
     end do
  end do



debug_3D(:,:,87)=S(:,:)-debug_3D(:,:,88)

  ENDDO


  TIME= TIME + DTT
  ISYNCHRO=1

  write(6,*) 'nbr d iterations dans surface', ntour
  
END SUBROUTINE SURFACE