source: trunk/SOURCES/strain_rate.f90 @ 23

!> \file strain_rate.f90
!! Calcule les taux de deformation et le second invariant du tenseur des deformations
!<

!> SUBROUTINE: STRAIN_RATE()
!! \author Vincent Peyaud 
!! \date octobre 2003
!! @note Calcul le taux de deformation horizontaux (moyenne verticale)
!! @note Cette version est adaptee a l'Antarctique 
!! @note Modele couple "ice sheet - ice shelves"
!! @note Used modules:
!! @note    - use module3D_phy
!<



!  *****************************************************************
!
!
!     Vincent Peyaud : octobre 2003
!     Calcul le taux de deformation horizontaux (moyenne verticale)
!     cette version est adaptee a l'Antarctique 
!     Modele couple "ice sheet - ice shelves"
!
!
!  *****************************************************************


      subroutine strain_rate()   

      use module3d_phy
 
      implicit none
    
      real dvdx,dudy

!time=tbegin
!les eps sont initialises a 0 dans INITIAL2

!=======calcul============================================!
      if (itracebug.eq.1)  call tracebug(' Entree dans routine strain_rate')

      DVDX= 0.
      DUDY= 0.
  
      do j=2,ny-1
      do i=2,nx-1

!les epsilons xx et yy sont toujours centres (noeud majeur)
         epsxx(i,j) = (uxbar(i+1,j) - uxbar(i,j))/dx
         epsyy(i,j)  = (uybar(i,j+1) - uybar(i,j))/dy

! calcul de epsxy sur les noeuds majeurs---------------------------------
!
                dudy = ((uxbar(i+1,j+1) - uxbar(i+1,j-1))    &
                     + (uxbar(i,j+1) - uxbar(i,j-1)))/(4.*dy)
                dvdx = ((uybar(i+1,j+1) - uybar(i-1,j+1))    &
                     + (uybar(i+1,j) - uybar(i-1,j)))/(4.*dx)

                epsxy(i,j) =0.5* (dudy+dvdx)        ! sur un noeud majeur
                eps(i,j)=epsxx(i,j)**2+epsyy(i,j)**2+epsxx(i,j)*epsyy(i,j)+epsxy(i,j)**2
!
! fin de calcul de epsxy sur les noeuds majeurs--------------------------

! calcul de epsxy sur les noeuds mineurs -------------------------------
! attention dans le calcul suivant epsxy est calcule sur le noeud mineur
! -1/2,-1/2
!         dudy = (uxbar(i,j) - uxbar(i,j-1))/dy
!         dvdx = (uybar(i,j) - uybar(i-1,j))/dx  

!         epsxy(i,j) =0.5* (dudy+dvdx)        ! sur un noeud mineur 

! attention dans le calcul suivant eps(i,j) ne tient pas compte
! du cisaillement horizontal.

!         eps(i,j)=epsxx(i,j)**2+epsyy(i,j)**2+epsxx(i,j)*epsyy(i,j)

! si epsxy est sur les noeuds mineurs
! il faudrait tenir compte du epsxy dans le calcul de pvm (voir diagno_L2)
! fin de calcul de epsxy sur les noeuds mineurs-----------------------------

         eps(i,j)=eps(i,j)**0.5

      end do
      end do

! attention : mise au point mismip
!eps(:,:) = abs(epsxx(:,:))
      return
      end subroutine