source: trunk/SOURCES/new-flot-0.3.f90 @ 23

Last change on this file since 23 was 4, checked in by dumas, 10 years ago

initial import GRISLI trunk
File size: 6.6 KB
Line 
1!> \file new-flot-0.3.f90
2!!  La routine calcule la vitesse des nouveau points flottant entre 2 pas de temps dtt.
3!<
4
5!> SUBROUTINE: new_flot
6!! \author ...
7!! \date 17 decembre 2001
8!! @note La routine calcul la vitesse des nouveau points flottant.
9!! @note Subroutine appelee lorsqu'un point devient flottant entre 2 pas
10!! de temps dtt.
11!!  @note Les nouveaux points flottants sont detectes par flottab
12!! @note  Pour le calcul de la vitesse on utilise l'equation de vitesse du shelf
13!! @note Les points que deviennent flottants sont .true. dans la variable
14!! tableau new_flot_point(nx,ny) new_flotmy(nx,ny) et new_flotmx(nx,ny)
15!! @note Les nouvelles vitesses sont mises en memoires dans une variable
16!! locale new_Uxbar(nx,ny) et new_Uybar(nx,ny)
17!! @note A la fin de la routine on remplace les anciennes valeurs par les nouvelles
18!! @note Used modules:
19!! @note    - use module3D_phy
20!! @note    - use param_phy_mod
21!<
22
23subroutine new_flot
24
25  USE module3D_phy
26  USE param_phy_mod
27
28  implicit none
29
30  REAL,dimension(nx,ny) :: new_Uxbar, new_Uybar !< nouvelles valeurs de Uxbar et Uybar
31
32  ! variables de calcul de la vitesse :
33  REAL :: TUIMJ, TUMIJ, TUIJ, TUPIJ, TUIPJ, TVMIJ, TVIJ, TVIPJ, TVMIPJ
34  REAL :: TVIMJ, TVPIJ, TUPIMJ
35  REAL :: beta !<      beta est le terme de frottement basal
36  REAL :: SCAL !< tous les termes de l'equation sont divise par SCAL
37  REAL :: moteur !< entre dans le calcul de BDR
38  REAL :: unorm !< pour le calcul de BDR
39  REAL :: BDR !< terme de droite de l'equation
40
41
42  DO j=3,ny-2
43     DO i=3,nx-2
44
45        IF (new_flotmx(i,j)) then
46
47           if (gzmx(i,j).and..not.ilemx(i,j)) then
48              beta= min(abs(tobmx(i,j)),betamax)  ! 12 juin 2000
49              FROTMX(I,J)=beta
50              unorm=((uybar(i,j)+uybar(i-1,j))+(uybar(i,j+1)+uybar(i-1,j+1)))/4.
51              unorm=unorm*unorm+uxbar(i,j)*uxbar(i,j)
52              unorm=max(1.,unorm)
53              unorm=sqrt(unorm)
54              beta=beta/unorm*dx*dx
55           else
56              beta=0.
57              FROTMX(I,J)=0.
58           endif
59
60
61           !   Terme en u(i,j):
62           !   _______________
63
64           TUIJ = -4.*PVI(I,J) - 4.*PVI(I-1,J) - PVM(i,j+1)-PVM(I,J)-beta
65
66           !   Terme en u(i-1,j):
67           !   _________________
68
69           TUMIJ = 4.*PVI(I-1,J)
70
71
72           !   Terme en u(i+1,j):
73           !   _________________
74
75           TUPIJ = 4.*PVI(I,J)
76
77
78           !   Terme en u(i,j-1):
79           !   _________________
80
81           TUIMJ = PVM(I,J)
82
83
84           !   Terme en u(i,j+1):
85           !   _________________
86
87           TUIPJ = PVM(I,J+1)
88
89
90           !   Terme en v(i,j):
91           !   _______________
92
93           TVIJ = -2.*PVI(I,J)-PVM(I,J)
94
95           !   Terme en v(i-1,j):
96           !   _________________
97
98           TVMIJ = 2.*PVI(I-1,J)+PVM(I,J)
99
100           !   Terme en v(i-1,j+1):
101           !   ___________________
102
103           TVMIPJ = -2.*PVI(I-1,J)-PVM(I,J+1)
104
105           !   Terme en v(i,j+1):
106           !   _________________
107
108           TVIPJ = 2.*PVI(I,J)+PVM(I,J+1)
109
110           ! on divise tous les termes de l'equation par SCAL
111           SCAL = TUIJ
112
113           ! terme de droite de l'equation : bdr dans remplimat
114           !      limitation de la force motrice
115           moteur= RO*G*HMX(I,J)*(S(I,J)-S(I-1,J))/dx
116           moteur=min(moteur,moteurmax)
117           moteur=max(moteur,-moteurmax)
118
119           ! modif tof le 3-09-2002
120           if (SCAL.NE.0.) then
121              BDR = moteur*dx*dx/SCAL 
122
123              ! calcul de U(i,j) :
124
125              new_Uxbar(i,j)= BDR - (1/SCAL)* (Uxbar(i-1,j)*TUMIJ &
126                   + Uxbar(i+1,j)*TUPIJ &
127                   + Uxbar(i,j-1)*TUIMJ + Uxbar(i,j+1)*TUIPJ + Uybar(i,j)*TVIJ &
128                   + Uybar(i-1,j)*TVMIJ + Uybar(i-1,j+1)*TVMIPJ + Uybar(i,j+1)*TVIPJ)
129           else
130              new_Uxbar(i,j)= Uxbar(i,j)
131           endif
132        endif
133
134
135
136        !-------------------------------------------------
137        ! Boucle pour les vitesses selon y
138        !-------------------------------------------------
139        !
140
141        IF (new_flotmy(i,j)) then
142
143           !      beta est le terme de frottement basal
144           if (gzmy(i,j).and..not.ilemy(i,j)) then
145              beta= min(abs(tobmy(i,j)),betamax)   ! 12 juin 2000
146              FROTMY(I,J)=beta
147              unorm=((uxbar(i,j)+uxbar(i,j-1))+(uxbar(i+1,j)+uxbar(i+1,j-1)))/4.
148              unorm=unorm*unorm+uybar(i,j)*uybar(i,j)
149              unorm=max(1.,unorm)
150              unorm=sqrt(unorm)
151              beta=beta/unorm*dx*dx
152           else
153              beta=0.
154              FROTMY(I,J)=0.
155           endif
156
157
158           ! Terme en v(i,j):
159           ! ________________
160
161           TVIJ = -4.*PVI(I,J)-4.*PVI(I,J-1) - PVM(I+1,J)-PVM(I,J)-beta
162
163
164           ! Terme en v(i,j-1):
165           ! __________________
166
167           TVIMJ = 4.*PVI(I,J-1)
168
169           ! Terme en v(i,j+1):
170           ! __________________
171
172           TVIPJ = 4.*PVI(I,J)
173
174           ! Terme en v(i-1,j):
175           ! __________________
176
177           TVMIJ = PVM(I,J) 
178
179           ! Terme en v(i+1,j):
180           ! __________________
181
182           TVPIJ = PVM(I+1,J)
183
184           ! Terme en u(i,j):
185           ! ________________
186
187           TUIJ = -2.*PVI(I,J)-PVM(I,J)
188
189           ! Terme en u(i,j-1):
190           ! __________________
191
192           TUIMJ = 2.*PVI(I,J-1)+PVM(I,J)
193
194
195           ! Terme en u(i+1,j-1):
196           ! ___________________
197
198           TUPIMJ = -2.*PVI(I,J-1)-PVM(I+1,J)
199
200
201           ! Terme en u(i+1,j):
202           ! ___________________
203
204           TUPIJ = 2.*PVI(I,J)+PVM(I+1,J)
205
206
207           ! on divise tous les termes de l'equation par SCAL
208           SCAL = TVIJ
209
210           ! terme de droite de l'equation : bdr dans remplimat
211           !      limitation de la force motrice
212           moteur= RO*G*HMY(I,J)*(S(I,J)-S(I,J-1))/dx
213           moteur=min(moteur,moteurmax)
214           moteur=max(moteur,-moteurmax)
215
216           ! modif tof 3-09-2002
217           if (SCAL.NE.0.) then
218              BDR= moteur*dx*dx/SCAL 
219
220              ! calcul de V(i,j) :
221
222              new_Uybar(i,j)= BDR - (1/SCAL)* (Uybar(i,j-1)*TVIMJ &
223                   + Uybar(i,j+1)*TVIPJ &
224                   + Uybar(i-1,j)*TVMIJ + Uybar(i+1,j)*TVPIJ + Uxbar(i,j)*TUIJ &
225                   + Uxbar(i,j-1)*TUIMJ + Uxbar(i+1,j-1)*TUPIMJ &
226                   + Uxbar(i+1,j)*TUPIJ)
227           else
228              new_Uybar(i,j)= Uybar(i,j)
229           endif
230        endif
231     enddo
232  enddo
233
234  ! remplacement des anciennes vitesses par les nouvelles
235  do j=3,ny-2
236     do i=3,nx-2
237        if (new_flotmx(i,j)) then
238           UXBAR(i,j)=new_uxbar(i,j)
239        endif
240        if (new_flotmy(i,j)) then
241           UYBAR(i,j)=new_uybar(i,j)
242        endif
243     enddo
244  enddo
245
246
247END SUBROUTINE new_flot
Note: See TracBrowser for help on using the repository browser.