!> \file velocities-polyn-0.3.f90
!! Calcule des vitesses. 
!<

!> SUBROUTINE: VELOCITIES()
!! \author Catherine RITZ  
!! \date mars 2007
!! @note Cette routine permet de calculer les vitesses.
!! @note Used modules:
!! @note    - use module3D_phy
!! @note    - use module_choix
!!
!<

! ==========================================================================
!
      subroutine SIA_velocities()
!
!
!                                   Cat passage en f90 en mars 2007
!                                   moyenne sa_mx et sa_my
!
!     Attention la vitesse verticale est "super reduite"
! ==========================================================================
!$ USE OMP_LIB
use module3d_phy
!use deform_declar 
USE module_choix
implicit none

REAL hdd
if (itracebug.eq.1)  call tracebug(' Entree dans routine SIA_velocity')

!$OMP PARALLEL
!$OMP DO COLLAPSE(2)    
do k=1,nz
   do j=2,ny-1
      do i=2,nx-1

!      flgzmx pour tous les points qui sont traites par l'equation elliptique
!      diagnoshelf

!     Selon X
!      -------------------------------------
!         if (.not.flgzmx(i,j)) then      ! utilise ubx au lieu de ddbx (compatibilite avec mix et spinup)
!            ux(i,j,k)=ddbx(i,j)
!            sux(i,j,k)=ddbx(i,j)*cde(k)

        if (.not.flotmx(i,j)) then        ! version aout 2010 (Cat) 
                                          ! utilise ubx au lieu de ddbx (compatibilite avec mix et spinup)
                                          ! et agit sur tous les points poses

            ux(i,j,k) = 0.
            sux(i,j,k)= 0.

            do iglen=n1poly,n2poly

               ux(i,j,k)  = ux(i,j,k)+ddx(i,j,iglen)*sa_mx(i,j,k,iglen)
               sux(i,j,k) = sux(i,j,k)+ddx(i,j,iglen)*s2a_mx(i,j,k,iglen)
            end do
            ux(i,j,k)  = ux(i,j,k)*(-sdx(i,j))
            ux(i,j,k)  = ux(i,j,k) + ubx(i,j)    ! ajoute le glissement

            sux(i,j,k) = sux(i,j,k)*sdx(i,j)*hmx(i,j)
            sux(i,j,k) = sux(i,j,k) - ubx(i,j) * cde(k)

         else
          sux(i,j,k) = (nz-k)*1./(nz-1.)*uxbar(i,j)*hmx(i,j)
          ux(i,j,k)  = uxbar(i,j)
         end if

!     Selon Y
!      -------------------------------------
!         if (.not.flgzmy(i,j)) then        
!            uy(i,j,k)=ddby(i,j)
!            suy(i,j,k)=ddby(i,j)*cde(k)

        if (.not.flotmy(i,j)) then        ! version aout 2010 (Cat) 
                                          ! utilise uby au lieu de ddby (compatibilite avec mix et spinup)
                                          ! et agit sur tous les points poses

            uy(i,j,k) = 0.
            suy(i,j,k)= 0.


            do iglen=n1poly,n2poly
               uy(i,j,k)=uy(i,j,k)+ddy(i,j,iglen)*sa_my(i,j,k,iglen)
               suy(i,j,k)=suy(i,j,k)+ddy(i,j,iglen)*s2a_my(i,j,k,iglen)
            end do
            uy(i,j,k)  = uy(i,j,k)*(-sdy(i,j))
            uy(i,j,k)  = uy(i,j,k) + uby(i,j)    ! ajoute le glissement

            suy(i,j,k) = suy(i,j,k)*sdy(i,j)*hmy(i,j)
            suy(i,j,k) = suy(i,j,k) - uby(i,j) * cde(k)

         else
          suy(i,j,k) = (nz-k)*1./(nz-1)*uybar(i,j)*hmy(i,j)
          uy(i,j,k)=uybar(i,j)
         endif

      end do
   end do
end do
!$OMP END DO

!       *************************** Z VELOCITIES ******************

!$OMP DO
do j=2,ny-1
   do i=2,nx-1

      divu(i,j)=((uxbar(i+1,j)*hmx(i+1,j)-uxbar(i,j)*hmx(i,j))  &          ! formulation explicite
           + (uybar(i,j+1)*hmy(i,j+1)-uybar(i,j)*hmy(i,j)))/dx

      uzr(i,j,nz)=bmelt(i,j)
      xx(i,j)=uzr(i,j,1)
      hdd=bm(i,j)-bmelt(i,j)-divu(i,j)                                     ! formulation explicite de dH/dt

!     -------------------------------------------------------------------
!      attention uzr contient maintenant :
!      uzr=uz +u ((1-e)*dH/dx+dB/dx)+ ((1-e)*dH/dt+dB/dt)

         do k=1,nz-1
            if (.not.flot(i,j).and.(front(i,j).eq.4.).and..not.flgzmx(i,j) &
                 .and..not.flgzmx(i+1,j).and..not.flgzmy(i,j) &
                 .and..not.flgzmy(i,j+1)) then
               uzr(i,j,k)=uzr(i,j,nz)-((sux(i+1,j,k)-sux(i,j,k)) &
                +(suy(i,j+1,k)-suy(i,j,k)))/dx
               uzr(i,j,k)=uzr(i,j,k)+hdd*cde(k)

            else

!     dans l'ice-shelf la vitesse uzr est lineaire ...
!    -------------------------------------------------

!     attention, dans les zones stream, le vitesse verticale n'est pas
!     lineaire si elle est voisine d'un point .not.flgzmx
!     le calcul se fait alors par differences finies

               uzr(i,j,k)= bm(i,j)+(k-1.)/(nz-1.)*(bmelt(i,j)-bm(i,j))
           end if

        end do
        
     end do
 end do
!$OMP END DO

!$OMP DO
  do j=2,ny-1
     do i=2,nx-1  
        uzsdot(i,j)=(uzr(i,j,1)-xx(i,j))/dtt
        uzk(i,j)=(-bdot(i,j)-hdot(i,j)+bm(i,j)-bmelt(i,j))
     end do
  end do
!$OMP END DO
!$OMP END PARALLEL

end subroutine SIA_velocities