source: branches/iLoveclim/SOURCES/Old-sources/calceps2-0.2.f @ 77

!> \file calceps2-0.2.f
!! Calcule les taux de deformation et le second invariant du tenseur des deformations
!<

!> SUBROUTINE: Calceps
!! \author Vince 
!! \date 6 jan 95 Modifie 20 jan 95
!! @note Cette routine permet de calculer les taux de deformation (epsxx epsyy epsxy)
!! @note et le second invariant du tenseur des deformations (epsilon)
!! @note qui sera reinjecte dans la boucle de resolution de l'equation diagnostique.
!! @note Used modules:
!! @note    - use module3D_phy
!!
!<
C ========================================================================
C
       SUBROUTINE CALCEPS(EPSILON)
C
C                                       Vince 6 jan 95
C                                       Modifie 20 jan 95
C                                       -----------------
C
C
C      Cette routine permet de calculer les taux de deformation
C      (epsxx epsyy epsxy)
C      et le second invariant du tenseur des deformations (epsilon)
C  qui sera reinjecte dans la boucle de resolution de l'equation diagnostique.
C
C =======================================================================

      USE module3D_phy


      IMPLICIT NONE

      REAL DUDX,DUDY,DVDX,DVDY,EPSILON(NX,NY)

      DUDX = 0.
      DUDY = 0.
      DVDX = 0.
      DVDY = 0.

      DO i=1,NX
         DO j=1,NY
!         IF (.NOT.FLOT(I,J).OR.H(i,j).eq.1.) GOTO 11
!         IF (.NOT.FRONT(I,J)) THEN
            IF (FRONT(I,J).eq.4) THEN
               DUDX = (UXBAR(I+1,J) - UXBAR(I,J))/DX

               DUDY = (UXBAR(I+1,J+1)+UXBAR(I,J+1)
     &               - UXBAR(I,J-1) - UXBAR(I+1,J-1))/(4.*DY)
               DVDX = (UYBAR(I+1,J+1)+UYBAR(I+1,J)
     &               - UYBAR(I-1,J) - UYBAR(I-1,J+1))/(4.*DX)
               DVDY = (UYBAR(I,J+1)-UYBAR(i,J))/DY


C      Cas particulier: on se trouve sur le front:
C      Si le front n'est pas sur les bords de la grille,
C      on a pas besoin de sa valeur de epsilon
C      (dans ce cas de toutes facons eps est voisin de 0)
C      ----------------------------------------------------

  
            ELSE IF (I.EQ.1.AND.J.NE.1.AND.J.NE.NY) THEN

               DUDX = (UXBAR(I+2,J) - UXBAR(I+1,J))/DX
               DVDX = 0.
               DUDY = 0.
               DVDY = (UYBAR(I,J+1) - UYBAR(I,J))/DY

            ELSE IF (I.EQ.NX.AND.J.NE.1.AND.J.NE.NY) THEN

               DUDX = (UXBAR(I,J) - UXBAR(I-1,J))/DX
               DVDX = 0.
               DUDY = 0.
               DVDY = (UYBAR(I,J+1) - UYBAR(I,J))/DY

            ELSE IF (J.EQ.1.AND.I.NE.1.AND.I.NE.NX) THEN
                
               DUDX = (UXBAR(I+1,J) - UXBAR(I,J))/DX
               DVDX = 0.
               DUDY = 0.
               DVDY = (UYBAR(I,J+2) - UYBAR(I,J+1))/DY

            ELSE IF (J.EQ.NY.AND.I.NE.1.AND.I.NE.NX) THEN

               DUDX = (UXBAR(I+1,J) - UXBAR(I,J))/DX
               DVDX = 0.
               DUDY = 0.
               DVDY = (UYBAR(I,J) - UYBAR(I,J-1))/DY

            ENDIF

            EPSILON(i,j)=(DUDX**2 + DVDY**2
     &        +(1./4.)*(DUDY+DVDX)**2+DUDX*DVDY)
            IF (EPSILON(I,J).gt.1.E-10) then 
               EPSILON(I,J)=EPSILON(I,J)**(0.5)
            else
               EPSILON(I,J)=1.E-5 
            endif
! mise en memoire de epsxx epsyy et epsxy :
! utilise dans icetemp pour chaleur de deformation          
            epsxx(I,J)=DUDX
            epsyy(I,J)=DVDY
            epsxy(I,J)=0.5*(DUDY+DVDX)

! 11       END DO
         END DO
      END DO   
c       open(66,file='M3-eps')
c         do i=1,nx
c  write(66,*)
c  write(66,*)'i=',i
c  write(66,'(5(e10.3,3x))') (epsilon(i,j),j=1,ny)
c         end do

      RETURN
      END