1 | !> \file velocities-polyn-0.3.f90 |
---|
2 | !! Calcule des vitesses. |
---|
3 | !< |
---|
4 | |
---|
5 | !> SUBROUTINE: VELOCITIES() |
---|
6 | !! \author Catherine RITZ |
---|
7 | !! \date mars 2007 |
---|
8 | !! @note Cette routine permet de calculer les vitesses. |
---|
9 | !! @note Used modules: |
---|
10 | !! @note - use module3D_phy |
---|
11 | !! @note - use module_choix |
---|
12 | !! |
---|
13 | !< |
---|
14 | |
---|
15 | ! ========================================================================== |
---|
16 | ! |
---|
17 | subroutine SIA_velocities() |
---|
18 | ! |
---|
19 | ! |
---|
20 | ! Cat passage en f90 en mars 2007 |
---|
21 | ! moyenne sa_mx et sa_my |
---|
22 | ! |
---|
23 | ! Attention la vitesse verticale est "super reduite" |
---|
24 | ! ========================================================================== |
---|
25 | |
---|
26 | use module3d_phy |
---|
27 | !use deform_declar |
---|
28 | USE module_choix |
---|
29 | implicit none |
---|
30 | |
---|
31 | REAL hdd |
---|
32 | if (itracebug.eq.1) call tracebug(' Entree dans routine SIA_velocity') |
---|
33 | |
---|
34 | |
---|
35 | do k=1,nz |
---|
36 | do j=2,ny-1 |
---|
37 | do i=2,nx-1 |
---|
38 | |
---|
39 | ! flgzmx pour tous les points qui sont traites par l'equation elliptique |
---|
40 | ! diagnoshelf |
---|
41 | |
---|
42 | ! Selon X |
---|
43 | ! ------------------------------------- |
---|
44 | ! if (.not.flgzmx(i,j)) then ! utilise ubx au lieu de ddbx (compatibilite avec mix et spinup) |
---|
45 | ! ux(i,j,k)=ddbx(i,j) |
---|
46 | ! sux(i,j,k)=ddbx(i,j)*cde(k) |
---|
47 | |
---|
48 | if (.not.flotmx(i,j)) then ! version aout 2010 (Cat) |
---|
49 | ! utilise ubx au lieu de ddbx (compatibilite avec mix et spinup) |
---|
50 | ! et agit sur tous les points poses |
---|
51 | |
---|
52 | ux(i,j,k) = 0. |
---|
53 | sux(i,j,k)= 0. |
---|
54 | |
---|
55 | do iglen=n1poly,n2poly |
---|
56 | |
---|
57 | ux(i,j,k) = ux(i,j,k)+ddx(i,j,iglen)*sa_mx(i,j,k,iglen) |
---|
58 | sux(i,j,k) = sux(i,j,k)+ddx(i,j,iglen)*s2a_mx(i,j,k,iglen) |
---|
59 | end do |
---|
60 | ux(i,j,k) = ux(i,j,k)*(-sdx(i,j)) |
---|
61 | ux(i,j,k) = ux(i,j,k) + ubx(i,j) ! ajoute le glissement |
---|
62 | |
---|
63 | sux(i,j,k) = sux(i,j,k)*sdx(i,j)*hmx(i,j) |
---|
64 | sux(i,j,k) = sux(i,j,k) - ubx(i,j) * cde(k) |
---|
65 | |
---|
66 | else |
---|
67 | sux(i,j,k) = (nz-k)*1./(nz-1.)*uxbar(i,j)*hmx(i,j) |
---|
68 | ux(i,j,k) = uxbar(i,j) |
---|
69 | end if |
---|
70 | |
---|
71 | ! Selon Y |
---|
72 | ! ------------------------------------- |
---|
73 | ! if (.not.flgzmy(i,j)) then |
---|
74 | ! uy(i,j,k)=ddby(i,j) |
---|
75 | ! suy(i,j,k)=ddby(i,j)*cde(k) |
---|
76 | |
---|
77 | if (.not.flotmy(i,j)) then ! version aout 2010 (Cat) |
---|
78 | ! utilise uby au lieu de ddby (compatibilite avec mix et spinup) |
---|
79 | ! et agit sur tous les points poses |
---|
80 | |
---|
81 | uy(i,j,k) = 0. |
---|
82 | suy(i,j,k)= 0. |
---|
83 | |
---|
84 | |
---|
85 | do iglen=n1poly,n2poly |
---|
86 | uy(i,j,k)=uy(i,j,k)+ddy(i,j,iglen)*sa_my(i,j,k,iglen) |
---|
87 | suy(i,j,k)=suy(i,j,k)+ddy(i,j,iglen)*s2a_my(i,j,k,iglen) |
---|
88 | end do |
---|
89 | uy(i,j,k) = uy(i,j,k)*(-sdy(i,j)) |
---|
90 | uy(i,j,k) = uy(i,j,k) + uby(i,j) ! ajoute le glissement |
---|
91 | |
---|
92 | suy(i,j,k) = suy(i,j,k)*sdy(i,j)*hmy(i,j) |
---|
93 | suy(i,j,k) = suy(i,j,k) - uby(i,j) * cde(k) |
---|
94 | |
---|
95 | else |
---|
96 | suy(i,j,k) = (nz-k)*1./(nz-1)*uybar(i,j)*hmy(i,j) |
---|
97 | uy(i,j,k)=uybar(i,j) |
---|
98 | endif |
---|
99 | |
---|
100 | end do |
---|
101 | end do |
---|
102 | end do |
---|
103 | |
---|
104 | |
---|
105 | ! *************************** Z VELOCITIES ****************** |
---|
106 | |
---|
107 | |
---|
108 | do j=2,ny-1 |
---|
109 | do i=2,nx-1 |
---|
110 | |
---|
111 | divu(i,j)=((uxbar(i+1,j)*hmx(i+1,j)-uxbar(i,j)*hmx(i,j)) & ! formulation explicite |
---|
112 | + (uybar(i,j+1)*hmy(i,j+1)-uybar(i,j)*hmy(i,j)))/dx |
---|
113 | |
---|
114 | uzr(i,j,nz)=bmelt(i,j) |
---|
115 | xx(i,j)=uzr(i,j,1) |
---|
116 | hdd=bm(i,j)-bmelt(i,j)-divu(i,j) ! formulation explicite de dH/dt |
---|
117 | |
---|
118 | ! ------------------------------------------------------------------- |
---|
119 | ! attention uzr contient maintenant : |
---|
120 | ! uzr=uz +u ((1-e)*dH/dx+dB/dx)+ ((1-e)*dH/dt+dB/dt) |
---|
121 | |
---|
122 | do k=1,nz-1 |
---|
123 | if (.not.flot(i,j).and.(front(i,j).eq.4.).and..not.flgzmx(i,j) & |
---|
124 | .and..not.flgzmx(i+1,j).and..not.flgzmy(i,j) & |
---|
125 | .and..not.flgzmy(i,j+1)) then |
---|
126 | uzr(i,j,k)=uzr(i,j,nz)-((sux(i+1,j,k)-sux(i,j,k)) & |
---|
127 | +(suy(i,j+1,k)-suy(i,j,k)))/dx |
---|
128 | uzr(i,j,k)=uzr(i,j,k)+hdd*cde(k) |
---|
129 | |
---|
130 | else |
---|
131 | |
---|
132 | ! dans l'ice-shelf la vitesse uzr est lineaire ... |
---|
133 | ! ------------------------------------------------- |
---|
134 | |
---|
135 | ! attention, dans les zones stream, le vitesse verticale n'est pas |
---|
136 | ! lineaire si elle est voisine d'un point .not.flgzmx |
---|
137 | ! le calcul se fait alors par differences finies |
---|
138 | |
---|
139 | uzr(i,j,k)= bm(i,j)+(k-1.)/(nz-1.)*(bmelt(i,j)-bm(i,j)) |
---|
140 | end if |
---|
141 | |
---|
142 | end do |
---|
143 | |
---|
144 | end do |
---|
145 | end do |
---|
146 | |
---|
147 | do j=2,ny-1 |
---|
148 | do i=2,nx-1 |
---|
149 | uzsdot(i,j)=(uzr(i,j,1)-xx(i,j))/dtt |
---|
150 | uzk(i,j)=(-bdot(i,j)-hdot(i,j)+bm(i,j)-bmelt(i,j)) |
---|
151 | end do |
---|
152 | end do |
---|
153 | |
---|
154 | |
---|
155 | end subroutine SIA_velocities |
---|