1 | !> \file strain_rate.f90 |
---|
2 | !! Calcule les taux de deformation et le second invariant du tenseur des deformations |
---|
3 | !< |
---|
4 | |
---|
5 | !> SUBROUTINE: STRAIN_RATE() |
---|
6 | !! \author Vincent Peyaud |
---|
7 | !! \date octobre 2003 |
---|
8 | !! @note Calcul le taux de deformation horizontaux (moyenne verticale) |
---|
9 | !! @note Cette version est adaptee a l'Antarctique |
---|
10 | !! @note Modele couple "ice sheet - ice shelves" |
---|
11 | !! @note Used modules: |
---|
12 | !! @note - use module3D_phy |
---|
13 | !< |
---|
14 | |
---|
15 | |
---|
16 | |
---|
17 | ! ***************************************************************** |
---|
18 | ! |
---|
19 | ! |
---|
20 | ! Vincent Peyaud : octobre 2003 |
---|
21 | ! Calcul le taux de deformation horizontaux (moyenne verticale) |
---|
22 | ! cette version est adaptee a l'Antarctique |
---|
23 | ! Modele couple "ice sheet - ice shelves" |
---|
24 | ! |
---|
25 | ! |
---|
26 | ! ***************************************************************** |
---|
27 | |
---|
28 | |
---|
29 | subroutine strain_rate() |
---|
30 | |
---|
31 | use module3d_phy |
---|
32 | |
---|
33 | implicit none |
---|
34 | |
---|
35 | real dvdx,dudy |
---|
36 | |
---|
37 | !time=tbegin |
---|
38 | !les eps sont initialises a 0 dans INITIAL2 |
---|
39 | |
---|
40 | !=======calcul============================================! |
---|
41 | if (itracebug.eq.1) call tracebug(' Entree dans routine strain_rate') |
---|
42 | |
---|
43 | DVDX= 0. |
---|
44 | DUDY= 0. |
---|
45 | |
---|
46 | do j=2,ny-1 |
---|
47 | do i=2,nx-1 |
---|
48 | |
---|
49 | !les epsilons xx et yy sont toujours centres (noeud majeur) |
---|
50 | epsxx(i,j) = (uxbar(i+1,j) - uxbar(i,j))/dx |
---|
51 | epsyy(i,j) = (uybar(i,j+1) - uybar(i,j))/dy |
---|
52 | |
---|
53 | ! calcul de epsxy sur les noeuds majeurs--------------------------------- |
---|
54 | ! |
---|
55 | dudy = ((uxbar(i+1,j+1) - uxbar(i+1,j-1)) & |
---|
56 | + (uxbar(i,j+1) - uxbar(i,j-1)))/(4.*dy) |
---|
57 | dvdx = ((uybar(i+1,j+1) - uybar(i-1,j+1)) & |
---|
58 | + (uybar(i+1,j) - uybar(i-1,j)))/(4.*dx) |
---|
59 | |
---|
60 | epsxy(i,j) =0.5* (dudy+dvdx) ! sur un noeud majeur |
---|
61 | eps(i,j)=epsxx(i,j)**2+epsyy(i,j)**2+epsxx(i,j)*epsyy(i,j)+epsxy(i,j)**2 |
---|
62 | ! |
---|
63 | ! fin de calcul de epsxy sur les noeuds majeurs-------------------------- |
---|
64 | |
---|
65 | ! calcul de epsxy sur les noeuds mineurs ------------------------------- |
---|
66 | ! attention dans le calcul suivant epsxy est calcule sur le noeud mineur |
---|
67 | ! -1/2,-1/2 |
---|
68 | ! dudy = (uxbar(i,j) - uxbar(i,j-1))/dy |
---|
69 | ! dvdx = (uybar(i,j) - uybar(i-1,j))/dx |
---|
70 | |
---|
71 | ! epsxy(i,j) =0.5* (dudy+dvdx) ! sur un noeud mineur |
---|
72 | |
---|
73 | ! attention dans le calcul suivant eps(i,j) ne tient pas compte |
---|
74 | ! du cisaillement horizontal. |
---|
75 | |
---|
76 | ! eps(i,j)=epsxx(i,j)**2+epsyy(i,j)**2+epsxx(i,j)*epsyy(i,j) |
---|
77 | |
---|
78 | ! si epsxy est sur les noeuds mineurs |
---|
79 | ! il faudrait tenir compte du epsxy dans le calcul de pvm (voir diagno_L2) |
---|
80 | ! fin de calcul de epsxy sur les noeuds mineurs----------------------------- |
---|
81 | |
---|
82 | eps(i,j)=eps(i,j)**0.5 |
---|
83 | |
---|
84 | end do |
---|
85 | end do |
---|
86 | |
---|
87 | ! attention : mise au point mismip |
---|
88 | !eps(:,:) = abs(epsxx(:,:)) |
---|
89 | return |
---|
90 | end subroutine |
---|
91 | |
---|