SICOPOLIS V3.0
 All Classes Files Functions Variables Macros
calc_vz_sia.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : c a l c _ v z _ s i a
4 !
5 !> @file
6 !!
7 !! Computation of the vertical velocity vz in the shallow ice approximation.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2013 Ralf Greve
12 !!
13 !! @section License
14 !!
15 !! This file is part of SICOPOLIS.
16 !!
17 !! SICOPOLIS is free software: you can redistribute it and/or modify
18 !! it under the terms of the GNU General Public License as published by
19 !! the Free Software Foundation, either version 3 of the License, or
20 !! (at your option) any later version.
21 !!
22 !! SICOPOLIS is distributed in the hope that it will be useful,
23 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
24 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 !! GNU General Public License for more details.
26 !!
27 !! You should have received a copy of the GNU General Public License
28 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
29 !<
30 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
31 
32 !-------------------------------------------------------------------------------
33 !> Computation of the vertical velocity vz in the shallow ice approximation.
34 !<------------------------------------------------------------------------------
35 subroutine calc_vz_sia(dxi, deta, dzeta_c, dzeta_t)
36 
37 use sico_types
38 use sico_variables
39 
40 implicit none
41 integer(i4b) :: i, j, kc, kt
42 real(dp) :: dxi, deta, dzeta_c, dzeta_t
43 real(dp) :: vz_m(0:jmax,0:imax)
44 real(dp) :: avz3(0:kcmax)
45 real(dp) :: cvz00, cvz0(0:ktmax), cvz1(0:ktmax), cvz2(0:ktmax), &
46  cvz3(0:kcmax), cvz4(0:kcmax), cvz5(0:kcmax)
47 real(dp) :: dxi_inv, deta_inv
48 
49 !-------- Term abbreviations --------
50 
51 dxi_inv = 1.0_dp/dxi
52 deta_inv = 1.0_dp/deta
53 
54 do kc=0, kcmax
55  avz3(kc) = deform*eaz_c(kc)/(ea-1.0_dp)*dzeta_c
56 end do
57 
58 do i=1, imax-1
59 do j=1, jmax-1
60 
61  if (maske(j,i) == 0) then ! glaciated land
62 
63 !-------- Abbreviations --------
64 
65  cvz00 = 0.5_dp*(vx_t(0,j,i)+vx_t(0,j,i-1))*dzb_dxi_g(j,i) &
66  +0.5_dp*(vy_t(0,j,i)+vy_t(0,j-1,i))*dzb_deta_g(j,i) &
67  +dzb_dtau(j,i)-q_b_tot(j,i)
68 
69  kt=0
70  cvz0(kt) = h_t(j,i)*insq_g11_g(j,i)*insq_g22_g(j,i) &
71  *( (vx_t(kt,j,i)*sq_g22_sgx(j,i) &
72  -vx_t(kt,j,i-1)*sq_g22_sgx(j,i-1))*dxi_inv &
73  +(vy_t(kt,j,i)*sq_g11_sgy(j,i) &
74  -vy_t(kt,j-1,i)*sq_g11_sgy(j-1,i))*deta_inv ) &
75  *dzeta_t
76  cvz1(kt) = (dzb_dxi_g(j,i)+zeta_t(kt)*dh_t_dxi_g(j,i)) &
77  *(vx_t(kt+1,j,i)+vx_t(kt+1,j,i-1) &
78  -vx_t(kt,j,i) -vx_t(kt,j,i-1))*0.5_dp
79  cvz2(kt) = (dzb_deta_g(j,i)+zeta_t(kt)*dh_t_deta_g(j,i)) &
80  *(vy_t(kt+1,j,i)+vy_t(kt+1,j-1,i) &
81  -vy_t(kt,j,i) -vy_t(kt,j-1,i))*0.5_dp
82 
83  do kt=1, ktmax-1
84  cvz0(kt) = h_t(j,i)*insq_g11_g(j,i)*insq_g22_g(j,i) &
85  *( (vx_t(kt,j,i)*sq_g22_sgx(j,i) &
86  -vx_t(kt,j,i-1)*sq_g22_sgx(j,i-1))*dxi_inv &
87  +(vy_t(kt,j,i)*sq_g11_sgy(j,i) &
88  -vy_t(kt,j-1,i)*sq_g11_sgy(j-1,i))*deta_inv ) &
89  *dzeta_t
90  cvz1(kt) = (dzb_dxi_g(j,i)+zeta_t(kt)*dh_t_dxi_g(j,i)) &
91  *(vx_t(kt+1,j,i)+vx_t(kt+1,j,i-1) &
92  -vx_t(kt-1,j,i)-vx_t(kt-1,j,i-1))*0.25_dp
93  cvz2(kt) = (dzb_deta_g(j,i)+zeta_t(kt)*dh_t_deta_g(j,i)) &
94  *(vy_t(kt+1,j,i)+vy_t(kt+1,j-1,i) &
95  -vy_t(kt-1,j,i)-vy_t(kt-1,j-1,i))*0.25_dp
96  end do
97 
98  kt=ktmax
99  cvz0(kt) = h_t(j,i)*insq_g11_g(j,i)*insq_g22_g(j,i) &
100  *( (vx_t(kt,j,i)*sq_g22_sgx(j,i) &
101  -vx_t(kt,j,i-1)*sq_g22_sgx(j,i-1))*dxi_inv &
102  +(vy_t(kt,j,i)*sq_g11_sgy(j,i) &
103  -vy_t(kt,j-1,i)*sq_g11_sgy(j-1,i))*deta_inv ) &
104  *dzeta_t
105  cvz1(kt) = (dzb_dxi_g(j,i)+zeta_t(kt)*dh_t_dxi_g(j,i)) &
106  *(vx_t(kt,j,i) +vx_t(kt,j,i-1) &
107  -vx_t(kt-1,j,i)-vx_t(kt-1,j,i-1))*0.5_dp
108  cvz2(kt) = (dzb_deta_g(j,i)+zeta_t(kt)*dh_t_deta_g(j,i)) &
109  *(vy_t(kt,j,i) +vy_t(kt,j-1,i) &
110  -vy_t(kt-1,j,i)-vy_t(kt-1,j-1,i))*0.5_dp
111 
112  kc=0
113  cvz3(kc) = avz3(kc)*h_c(j,i) &
114  *insq_g11_g(j,i)*insq_g22_g(j,i) &
115  *( (vx_c(kc,j,i)*sq_g22_sgx(j,i) &
116  -vx_c(kc,j,i-1)*sq_g22_sgx(j,i-1))*dxi_inv &
117  +(vy_c(kc,j,i)*sq_g11_sgy(j,i) &
118  -vy_c(kc,j-1,i)*sq_g11_sgy(j-1,i))*deta_inv )
119  cvz4(kc) = (dzm_dxi_g(j,i) &
120  +eaz_c_quotient(kc)*dh_c_dxi_g(j,i)) &
121  *(vx_c(kc+1,j,i)+vx_c(kc+1,j,i-1) &
122  -vx_c(kc,j,i) -vx_c(kc,j,i-1))*0.5_dp
123  cvz5(kc) = (dzm_deta_g(j,i) &
124  +eaz_c_quotient(kc)*dh_c_deta_g(j,i)) &
125  *(vy_c(kc+1,j,i)+vy_c(kc+1,j-1,i) &
126  -vy_c(kc,j,i) -vy_c(kc,j-1,i))*0.5_dp
127 
128  do kc=1, kcmax-1
129  cvz3(kc) = avz3(kc)*h_c(j,i) &
130  *insq_g11_g(j,i)*insq_g22_g(j,i) &
131  *( (vx_c(kc,j,i)*sq_g22_sgx(j,i) &
132  -vx_c(kc,j,i-1)*sq_g22_sgx(j,i-1))*dxi_inv &
133  +(vy_c(kc,j,i)*sq_g11_sgy(j,i) &
134  -vy_c(kc,j-1,i)*sq_g11_sgy(j-1,i))*deta_inv )
135  cvz4(kc) = (dzm_dxi_g(j,i) &
136  +eaz_c_quotient(kc)*dh_c_dxi_g(j,i)) &
137  *(vx_c(kc+1,j,i)+vx_c(kc+1,j,i-1) &
138  -vx_c(kc-1,j,i)-vx_c(kc-1,j,i-1))*0.25_dp
139  cvz5(kc) = (dzm_deta_g(j,i) &
140  +eaz_c_quotient(kc)*dh_c_deta_g(j,i)) &
141  *(vy_c(kc+1,j,i)+vy_c(kc+1,j-1,i) &
142  -vy_c(kc-1,j,i)-vy_c(kc-1,j-1,i))*0.25_dp
143  end do
144 
145  kc=kcmax
146  cvz3(kc) = avz3(kc)*h_c(j,i) &
147  *insq_g11_g(j,i)*insq_g22_g(j,i) &
148  *( (vx_c(kc,j,i)*sq_g22_sgx(j,i) &
149  -vx_c(kc,j,i-1)*sq_g22_sgx(j,i-1))*dxi_inv &
150  +(vy_c(kc,j,i)*sq_g11_sgy(j,i) &
151  -vy_c(kc,j-1,i)*sq_g11_sgy(j-1,i))*deta_inv )
152  cvz4(kc) = (dzm_dxi_g(j,i) &
153  +eaz_c_quotient(kc)*dh_c_dxi_g(j,i)) &
154  *(vx_c(kc,j,i) +vx_c(kc,j,i-1) &
155  -vx_c(kc-1,j,i)-vx_c(kc-1,j,i-1))*0.5_dp
156  cvz5(kc) = (dzm_deta_g(j,i) &
157  +eaz_c_quotient(kc)*dh_c_deta_g(j,i)) &
158  *(vy_c(kc,j,i) +vy_c(kc,j-1,i) &
159  -vy_c(kc-1,j,i)-vy_c(kc-1,j-1,i))*0.5_dp
160 
161 !-------- Computation of vz_b --------
162 
163  vz_b(j,i) = cvz00
164 
165 !-------- Computation of vz --------
166 
167  if ((n_cts(j,i) == -1).or.(n_cts(j,i) == 0)) then
168  ! cold ice base, temperate ice base
169 
170  do kt=0, ktmax-1
171  vz_t(kt,j,i) = vz_b(j,i)
172  end do
173 
174  vz_m(j,i) = vz_b(j,i)
175 
176  vz_c(0,j,i) = vz_m(j,i) &
177  +0.5_dp*(-cvz3(0)+cvz4(0)+cvz5(0))
178 
179  do kc=1, kcmax-1
180  vz_c(kc,j,i) = vz_c(kc-1,j,i) &
181  -cvz3(kc)+cvz4(kc)+cvz5(kc)
182  end do
183 
184  else ! n_cts(j,i) == 1, temperate ice layer
185 
186  vz_t(0,j,i) = vz_b(j,i) &
187  +0.5_dp*(-cvz0(0)+cvz1(0)+cvz2(0))
188 
189  do kt=1, ktmax-1
190  vz_t(kt,j,i) = vz_t(kt-1,j,i) &
191  -cvz0(kt)+cvz1(kt)+cvz2(kt)
192  end do
193 
194  vz_m(j,i) = vz_t(ktmax-1,j,i) &
195  +0.5_dp*(-cvz0(ktmax)+cvz1(ktmax)+cvz2(ktmax))
196 
197  vz_c(0,j,i) = vz_m(j,i) &
198  +0.5_dp*(-cvz3(0)+cvz4(0)+cvz5(0))
199 
200  do kc=1, kcmax-1
201  vz_c(kc,j,i) = vz_c(kc-1,j,i) &
202  -cvz3(kc)+cvz4(kc)+cvz5(kc)
203  end do
204 
205  end if
206 
207 !-------- Computation of vz_s --------
208 
209  kc=kcmax
210 
211  vz_s(j,i) = vz_c(kc-1,j,i) &
212  +0.5_dp*(-cvz3(kc)+cvz4(kc)+cvz5(kc))
213 
214  else ! maske(j,i) == (1 or 2)
215 
216  vz_b(j,i) = 0.0_dp
217 
218  do kt=0, ktmax-1
219  vz_t(kt,j,i) = 0.0_dp
220  end do
221 
222  do kc=0, kcmax-1
223  vz_c(kc,j,i) = 0.0_dp
224  end do
225 
226  vz_s(j,i) = 0.0_dp
227 
228  end if
229 
230 end do
231 end do
232 
233 end subroutine calc_vz_sia
234 !
Generated on Mon May 20 2013 10:18:41 for SICOPOLIS V3.0 by   doxygen 1.8.1