source: branches/GRISLIv3/SOURCES/velocities-polyn-0.3.f90 @ 443

!> \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 runparam, only : itracebug      
  use module3d_phy, only: nx,ny,nz,dx,flotmx,flotmy,ux,uy,sux,suy,sdx,sdy,ubx,uby,&
       hmx,hmy,uxbar,uybar,iglen,cde,divu,uzr,bmelt,bm,flot,front,flgzmx,flgzmy,&
       uzsdot,dtt,uzk,bdot,hdot
  use deformation_mod_2lois, only: n1poly,n2poly,ddx,ddy,sa_mx,sa_my,s2a_mx,s2a_my

  implicit none

  integer :: i,j,k
  real, dimension(nx,ny) :: hdd
  real, dimension(nx,ny) :: xx

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

  !$OMP PARALLEL
  xx(:,:)=0.
  !$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(i,j)=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(i,j)*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