!> \file new-flot-0.3.f90
!!  La routine calcule la vitesse des nouveau points flottant entre 2 pas de temps dtt.
!<

!> SUBROUTINE: new_flot
!! \author ...
!! \date 17 decembre 2001
!! @note La routine calcul la vitesse des nouveau points flottant.
!! @note Subroutine appelee lorsqu'un point devient flottant entre 2 pas
!! de temps dtt.
!!  @note Les nouveaux points flottants sont detectes par flottab
!! @note  Pour le calcul de la vitesse on utilise l'equation de vitesse du shelf
!! @note Les points que deviennent flottants sont .true. dans la variable
!! tableau new_flot_point(nx,ny) new_flotmy(nx,ny) et new_flotmx(nx,ny)
!! @note Les nouvelles vitesses sont mises en memoires dans une variable 
!! locale new_Uxbar(nx,ny) et new_Uybar(nx,ny)
!! @note A la fin de la routine on remplace les anciennes valeurs par les nouvelles
!! @note Used modules:
!! @note    - use module3D_phy
!! @note    - use param_phy_mod
!<

subroutine new_flot

  USE module3D_phy
  USE param_phy_mod

  implicit none

  REAL,dimension(nx,ny) :: new_Uxbar, new_Uybar !< nouvelles valeurs de Uxbar et Uybar

  ! variables de calcul de la vitesse :
  REAL :: TUIMJ, TUMIJ, TUIJ, TUPIJ, TUIPJ, TVMIJ, TVIJ, TVIPJ, TVMIPJ
  REAL :: TVIMJ, TVPIJ, TUPIMJ
  REAL :: beta !<      beta est le terme de frottement basal
  REAL :: SCAL !< tous les termes de l'equation sont divise par SCAL
  REAL :: moteur !< entre dans le calcul de BDR
  REAL :: unorm !< pour le calcul de BDR
  REAL :: BDR !< terme de droite de l'equation


  DO j=3,ny-2
     DO i=3,nx-2

        IF (new_flotmx(i,j)) then

           if (gzmx(i,j).and..not.ilemx(i,j)) then 
              beta= min(abs(tobmx(i,j)),betamax)  ! 12 juin 2000
              FROTMX(I,J)=beta
              unorm=((uybar(i,j)+uybar(i-1,j))+(uybar(i,j+1)+uybar(i-1,j+1)))/4.
              unorm=unorm*unorm+uxbar(i,j)*uxbar(i,j)
              unorm=max(1.,unorm)
              unorm=sqrt(unorm)
              beta=beta/unorm*dx*dx
           else
              beta=0.
              FROTMX(I,J)=0.
           endif


           !   Terme en u(i,j):
           !   _______________

           TUIJ = -4.*PVI(I,J) - 4.*PVI(I-1,J) - PVM(i,j+1)-PVM(I,J)-beta

           !   Terme en u(i-1,j):
           !   _________________

           TUMIJ = 4.*PVI(I-1,J)


           !   Terme en u(i+1,j):
           !   _________________

           TUPIJ = 4.*PVI(I,J)


           !   Terme en u(i,j-1):
           !   _________________

           TUIMJ = PVM(I,J)


           !   Terme en u(i,j+1):
           !   _________________

           TUIPJ = PVM(I,J+1)


           !   Terme en v(i,j):
           !   _______________

           TVIJ = -2.*PVI(I,J)-PVM(I,J)

           !   Terme en v(i-1,j):
           !   _________________

           TVMIJ = 2.*PVI(I-1,J)+PVM(I,J)

           !   Terme en v(i-1,j+1):
           !   ___________________

           TVMIPJ = -2.*PVI(I-1,J)-PVM(I,J+1)

           !   Terme en v(i,j+1):
           !   _________________

           TVIPJ = 2.*PVI(I,J)+PVM(I,J+1)

           ! on divise tous les termes de l'equation par SCAL
           SCAL = TUIJ

           ! terme de droite de l'equation : bdr dans remplimat
           !      limitation de la force motrice
           moteur= RO*G*HMX(I,J)*(S(I,J)-S(I-1,J))/dx
           moteur=min(moteur,moteurmax)
           moteur=max(moteur,-moteurmax)

           ! modif tof le 3-09-2002
           if (SCAL.NE.0.) then
              BDR = moteur*dx*dx/SCAL 

              ! calcul de U(i,j) :

              new_Uxbar(i,j)= BDR - (1/SCAL)* (Uxbar(i-1,j)*TUMIJ &
                   + Uxbar(i+1,j)*TUPIJ &
                   + Uxbar(i,j-1)*TUIMJ + Uxbar(i,j+1)*TUIPJ + Uybar(i,j)*TVIJ &
                   + Uybar(i-1,j)*TVMIJ + Uybar(i-1,j+1)*TVMIPJ + Uybar(i,j+1)*TVIPJ)
           else
              new_Uxbar(i,j)= Uxbar(i,j)
           endif
        endif



        !-------------------------------------------------
        ! Boucle pour les vitesses selon y
        !-------------------------------------------------
        !

        IF (new_flotmy(i,j)) then

           !      beta est le terme de frottement basal
           if (gzmy(i,j).and..not.ilemy(i,j)) then 
              beta= min(abs(tobmy(i,j)),betamax)   ! 12 juin 2000
              FROTMY(I,J)=beta
              unorm=((uxbar(i,j)+uxbar(i,j-1))+(uxbar(i+1,j)+uxbar(i+1,j-1)))/4.
              unorm=unorm*unorm+uybar(i,j)*uybar(i,j)
              unorm=max(1.,unorm)
              unorm=sqrt(unorm)
              beta=beta/unorm*dx*dx
           else
              beta=0.
              FROTMY(I,J)=0.
           endif


           ! Terme en v(i,j):
           ! ________________

           TVIJ = -4.*PVI(I,J)-4.*PVI(I,J-1) - PVM(I+1,J)-PVM(I,J)-beta


           ! Terme en v(i,j-1):
           ! __________________

           TVIMJ = 4.*PVI(I,J-1)

           ! Terme en v(i,j+1):
           ! __________________

           TVIPJ = 4.*PVI(I,J)

           ! Terme en v(i-1,j):
           ! __________________

           TVMIJ = PVM(I,J) 

           ! Terme en v(i+1,j):
           ! __________________

           TVPIJ = PVM(I+1,J)

           ! Terme en u(i,j):
           ! ________________

           TUIJ = -2.*PVI(I,J)-PVM(I,J)

           ! Terme en u(i,j-1):
           ! __________________

           TUIMJ = 2.*PVI(I,J-1)+PVM(I,J)


           ! Terme en u(i+1,j-1):
           ! ___________________

           TUPIMJ = -2.*PVI(I,J-1)-PVM(I+1,J)


           ! Terme en u(i+1,j):
           ! ___________________

           TUPIJ = 2.*PVI(I,J)+PVM(I+1,J)


           ! on divise tous les termes de l'equation par SCAL
           SCAL = TVIJ

           ! terme de droite de l'equation : bdr dans remplimat
           !      limitation de la force motrice
           moteur= RO*G*HMY(I,J)*(S(I,J)-S(I,J-1))/dx
           moteur=min(moteur,moteurmax)
           moteur=max(moteur,-moteurmax)

           ! modif tof 3-09-2002
           if (SCAL.NE.0.) then
              BDR= moteur*dx*dx/SCAL 

              ! calcul de V(i,j) :

              new_Uybar(i,j)= BDR - (1/SCAL)* (Uybar(i,j-1)*TVIMJ &
                   + Uybar(i,j+1)*TVIPJ &
                   + Uybar(i-1,j)*TVMIJ + Uybar(i+1,j)*TVPIJ + Uxbar(i,j)*TUIJ &
                   + Uxbar(i,j-1)*TUIMJ + Uxbar(i+1,j-1)*TUPIMJ &
                   + Uxbar(i+1,j)*TUPIJ)
           else
              new_Uybar(i,j)= Uybar(i,j)
           endif
        endif
     enddo
  enddo

  ! remplacement des anciennes vitesses par les nouvelles
  do j=3,ny-2
     do i=3,nx-2
        if (new_flotmx(i,j)) then
           UXBAR(i,j)=new_uxbar(i,j)
        endif
        if (new_flotmy(i,j)) then
           UYBAR(i,j)=new_uybar(i,j)
        endif
     enddo
  enddo


END SUBROUTINE new_flot