SICOPOLIS V3.1
 All Classes Files Functions Variables Macros
calc_vxy_ssa.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : c a l c _ v x y _ s s a
4 !
5 !> @file
6 !!
7 !! Computation of the horizontal velocity vx, vy, the horizontal volume flux
8 !! qx, qy and the flux across the grounding line q_gl_g in the shallow shelf
9 !! approximation.
10 !!
11 !! @section Copyright
12 !!
13 !! Copyright 2009-2013 Ralf Greve, Tatsuru Sato
14 !!
15 !! @section License
16 !!
17 !! This file is part of SICOPOLIS.
18 !!
19 !! SICOPOLIS is free software: you can redistribute it and/or modify
20 !! it under the terms of the GNU General Public License as published by
21 !! the Free Software Foundation, either version 3 of the License, or
22 !! (at your option) any later version.
23 !!
24 !! SICOPOLIS is distributed in the hope that it will be useful,
25 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
26 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
27 !! GNU General Public License for more details.
28 !!
29 !! You should have received a copy of the GNU General Public License
30 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
31 !<
32 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
33 
34 !-------------------------------------------------------------------------------
35 !> Computation of the horizontal velocity vx, vy, the horizontal volume flux
36 !! qx, qy and the flux across the grounding line q_gl_g in the shallow shelf
37 !! approximation.
38 !<------------------------------------------------------------------------------
39 subroutine calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, ii, jj, nn)
40 
41 use sico_types
43 
44 implicit none
45 
46 integer(i4b), intent(in) :: ii((imax+1)*(jmax+1)), jj((imax+1)*(jmax+1))
47 integer(i4b), intent(in) :: nn(0:jmax,0:imax)
48 real(dp), intent(in) :: z_sl, dxi, deta, dzeta_c
49 
50 integer(i4b) :: i, j, kc, kt, m
51 real(dp), dimension(0:JMAX,0:IMAX) :: vx_m_prev, vy_m_prev
52 real(dp) :: dxi_inv, deta_inv
53 real(dp) :: vh_max
54 real(dp) :: qx_gl_g, qy_gl_g
55 
56 integer(i4b), parameter :: iter_ssa = 2 ! number of iterations
57 real(dp), parameter :: rel = 0.7_dp ! relaxation factor
58 
59 !-------- Detection of the grounding line and the calving front --------
60 
61 flag_grounding_line = .false.
62 flag_calving_front = .false.
63 
64 do i=1, imax-1
65 do j=1, jmax-1
66 
67  if ( (maske(j,i)==0_i2b) & ! grounded ice
68  .and. &
69  ( (maske(j,i+1)==3_i2b) & ! with
70  .or.(maske(j,i-1)==3_i2b) & ! one
71  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
72  .or.(maske(j-1,i)==3_i2b) ) & ! ice shelf point
73  ) &
74  flag_grounding_line(j,i) = .true.
75 
76  if ( (maske(j,i)==3_i2b) & ! floating ice
77  .and. &
78  ( (maske(j,i+1)==2_i2b) & ! with
79  .or.(maske(j,i-1)==2_i2b) & ! one
80  .or.(maske(j+1,i)==2_i2b) & ! neighbouring
81  .or.(maske(j-1,i)==2_i2b) ) & ! ocean point
82  ) &
83  flag_calving_front(j,i) = .true.
84 
85 end do
86 end do
87 
88 !-------- First iteration --------
89 
90 m=1
91 
92 ! ------ Depth-integrated viscosity vis_int_g
93 
94 do i=0, imax
95 do j=0, jmax
96  vis_int_g(j,i) = 1.0e+15_dp*(h_c(j,i)+h_t(j,i))
97  ! viscosity of 10^15 Pa s assumed
98 end do
99 end do
100 
101 ! ------ Horizontal velocity vx_m, vy_m
102 
103 call calc_vxy_ssa_matrix(z_sl, dxi, deta, ii, jj, nn)
104 
105 !-------- Further iterations --------
106 
107 do m=2, iter_ssa
108 
109 ! ------ Save velocities from previous iteration
110 
111  vx_m_prev = vx_m
112  vy_m_prev = vy_m
113 
114 ! ------ Depth-integrated viscosity vis_int_g
115 
116  call calc_vis_ssa(dxi, deta, dzeta_c)
117 
118 ! ------ Horizontal velocity vx_m, vy_m
119 
120  call calc_vxy_ssa_matrix(z_sl, dxi, deta, ii, jj, nn)
121 
122 ! ------ Relaxation scheme
123 
124  do i=0, imax
125  do j=0, jmax
126  vx_m(j,i) = rel*vx_m(j,i) + (1.0_dp-rel)*vx_m_prev(j,i)
127  vy_m(j,i) = rel*vy_m(j,i) + (1.0_dp-rel)*vy_m_prev(j,i)
128  end do
129  end do
130 
131 end do
132 
133 !-------- Limitation of computed vx_m and vy_m to the interval
134 ! [-VH_MAX, VH_MAX] --------
135 
136 vh_max = vh_max/year_sec
137 
138 do i=0, imax
139 do j=0, jmax
140  vx_m(j,i) = max(vx_m(j,i), -vh_max)
141  vx_m(j,i) = min(vx_m(j,i), vh_max)
142  vy_m(j,i) = max(vy_m(j,i), -vh_max)
143  vy_m(j,i) = min(vy_m(j,i), vh_max)
144 end do
145 end do
146 
147 !-------- 3-D velocities, basal velocities and volume flux --------
148 
149 ! ------ x-component
150 
151 do i=0, imax-1
152 do j=0, jmax
153 
154  if ( (maske(j,i)==3_i2b).or.(maske(j,i+1)==3_i2b) ) then
155  ! at least one neighbour point is floating ice
156  do kt=0, ktmax
157  vx_t(kt,j,i) = vx_m(j,i)
158  end do
159 
160  do kc=0, kcmax
161  vx_c(kc,j,i) = vx_m(j,i)
162  end do
163 
164  vx_b(j,i) = vx_m(j,i)
165 
166  qx(j,i) = vx_m(j,i) * 0.5_dp*(h_c(j,i)+h_t(j,i)+h_c(j,i+1)+h_t(j,i+1))
167 
168  end if
169 
170 end do
171 end do
172 
173 ! ------ y-component
174 
175 do i=0, imax
176 do j=0, jmax-1
177 
178  if ( (maske(j,i)==3_i2b).or.(maske(j+1,i)==3_i2b) ) then
179  ! at least one neighbour point is floating ice
180  do kt=0, ktmax
181  vy_t(kt,j,i) = vy_m(j,i)
182  end do
183 
184  do kc=0, kcmax
185  vy_c(kc,j,i) = vy_m(j,i)
186  end do
187 
188  vy_b(j,i) = vy_m(j,i)
189 
190  qy(j,i) = vy_m(j,i) * 0.5_dp*(h_c(j,i)+h_t(j,i)+h_c(j+1,i)+h_t(j+1,i))
191 
192  end if
193 
194 end do
195 end do
196 
197 !-------- Surface and basal velocities vx_s_g vy_s_g, vx_b_g vy_b_g
198 ! (defined at (i,j)) --------
199 
200 do i=1, imax-1
201 do j=1, jmax-1
202 
203  if (maske(j,i)==3_i2b) then
204 
205  vx_s_g(j,i) = 0.5_dp*(vx_m(j,i-1)+vx_m(j,i))
206  vx_b_g(j,i) = vx_s_g(j,i)
207 
208  vy_s_g(j,i) = 0.5_dp*(vy_m(j-1,i)+vy_m(j,i))
209  vy_b_g(j,i) = vy_s_g(j,i)
210 
211  end if
212 
213 end do
214 end do
215 
216 !-------- Computation of the flux across the grounding line q_gl_g
217 
218 do i=1, imax-1
219 do j=1, jmax-1
220 
221  if ( flag_grounding_line(j,i) ) then ! grounding line
222 
223  qx_gl_g = 0.5_dp*(qx(j,i)+qx(j,i-1))
224  qy_gl_g = 0.5_dp*(qy(j,i)+qy(j-1,i))
225 
226  q_gl_g(j,i) = sqrt(qx_gl_g*qx_gl_g+qy_gl_g*qy_gl_g)
227 
228  end if
229 
230 end do
231 end do
232 
233 end subroutine calc_vxy_ssa
234 !