SICOPOLIS V3.1
 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 
42 real(dp), intent(in) :: dxi, deta, dzeta_c, dzeta_t
43 
44 integer(i4b) :: i, j, kc, kt
45 real(dp) :: vz_m(0:jmax,0:imax)
46 real(dp) :: avz3(0:kcmax)
47 real(dp) :: cvz00, cvz0(0:ktmax), cvz1(0:ktmax), cvz2(0:ktmax), &
48  cvz3(0:kcmax), cvz4(0:kcmax), cvz5(0:kcmax)
49 real(dp) :: dxi_inv, deta_inv
50 
51 !-------- Term abbreviations --------
52 
53 dxi_inv = 1.0_dp/dxi
54 deta_inv = 1.0_dp/deta
55 
56 do kc=0, kcmax
57  avz3(kc) = deform*eaz_c(kc)/(ea-1.0_dp)*dzeta_c
58 end do
59 
60 do i=1, imax-1
61 do j=1, jmax-1
62 
63  if (maske(j,i) == 0) then ! glaciated land
64 
65 !-------- Abbreviations --------
66 
67  cvz00 = 0.5_dp*(vx_t(0,j,i)+vx_t(0,j,i-1))*dzb_dxi_g(j,i) &
68  +0.5_dp*(vy_t(0,j,i)+vy_t(0,j-1,i))*dzb_deta_g(j,i) &
69  +dzb_dtau(j,i)-q_b_tot(j,i)
70 
71  kt=0
72  cvz0(kt) = h_t(j,i)*insq_g11_g(j,i)*insq_g22_g(j,i) &
73  *( (vx_t(kt,j,i)*sq_g22_sgx(j,i) &
74  -vx_t(kt,j,i-1)*sq_g22_sgx(j,i-1))*dxi_inv &
75  +(vy_t(kt,j,i)*sq_g11_sgy(j,i) &
76  -vy_t(kt,j-1,i)*sq_g11_sgy(j-1,i))*deta_inv ) &
77  *dzeta_t
78  cvz1(kt) = (dzb_dxi_g(j,i)+zeta_t(kt)*dh_t_dxi_g(j,i)) &
79  *(vx_t(kt+1,j,i)+vx_t(kt+1,j,i-1) &
80  -vx_t(kt,j,i) -vx_t(kt,j,i-1))*0.5_dp
81  cvz2(kt) = (dzb_deta_g(j,i)+zeta_t(kt)*dh_t_deta_g(j,i)) &
82  *(vy_t(kt+1,j,i)+vy_t(kt+1,j-1,i) &
83  -vy_t(kt,j,i) -vy_t(kt,j-1,i))*0.5_dp
84 
85  do kt=1, ktmax-1
86  cvz0(kt) = h_t(j,i)*insq_g11_g(j,i)*insq_g22_g(j,i) &
87  *( (vx_t(kt,j,i)*sq_g22_sgx(j,i) &
88  -vx_t(kt,j,i-1)*sq_g22_sgx(j,i-1))*dxi_inv &
89  +(vy_t(kt,j,i)*sq_g11_sgy(j,i) &
90  -vy_t(kt,j-1,i)*sq_g11_sgy(j-1,i))*deta_inv ) &
91  *dzeta_t
92  cvz1(kt) = (dzb_dxi_g(j,i)+zeta_t(kt)*dh_t_dxi_g(j,i)) &
93  *(vx_t(kt+1,j,i)+vx_t(kt+1,j,i-1) &
94  -vx_t(kt-1,j,i)-vx_t(kt-1,j,i-1))*0.25_dp
95  cvz2(kt) = (dzb_deta_g(j,i)+zeta_t(kt)*dh_t_deta_g(j,i)) &
96  *(vy_t(kt+1,j,i)+vy_t(kt+1,j-1,i) &
97  -vy_t(kt-1,j,i)-vy_t(kt-1,j-1,i))*0.25_dp
98  end do
99 
100  kt=ktmax
101  cvz0(kt) = h_t(j,i)*insq_g11_g(j,i)*insq_g22_g(j,i) &
102  *( (vx_t(kt,j,i)*sq_g22_sgx(j,i) &
103  -vx_t(kt,j,i-1)*sq_g22_sgx(j,i-1))*dxi_inv &
104  +(vy_t(kt,j,i)*sq_g11_sgy(j,i) &
105  -vy_t(kt,j-1,i)*sq_g11_sgy(j-1,i))*deta_inv ) &
106  *dzeta_t
107  cvz1(kt) = (dzb_dxi_g(j,i)+zeta_t(kt)*dh_t_dxi_g(j,i)) &
108  *(vx_t(kt,j,i) +vx_t(kt,j,i-1) &
109  -vx_t(kt-1,j,i)-vx_t(kt-1,j,i-1))*0.5_dp
110  cvz2(kt) = (dzb_deta_g(j,i)+zeta_t(kt)*dh_t_deta_g(j,i)) &
111  *(vy_t(kt,j,i) +vy_t(kt,j-1,i) &
112  -vy_t(kt-1,j,i)-vy_t(kt-1,j-1,i))*0.5_dp
113 
114  kc=0
115  cvz3(kc) = avz3(kc)*h_c(j,i) &
116  *insq_g11_g(j,i)*insq_g22_g(j,i) &
117  *( (vx_c(kc,j,i)*sq_g22_sgx(j,i) &
118  -vx_c(kc,j,i-1)*sq_g22_sgx(j,i-1))*dxi_inv &
119  +(vy_c(kc,j,i)*sq_g11_sgy(j,i) &
120  -vy_c(kc,j-1,i)*sq_g11_sgy(j-1,i))*deta_inv )
121  cvz4(kc) = (dzm_dxi_g(j,i) &
122  +eaz_c_quotient(kc)*dh_c_dxi_g(j,i)) &
123  *(vx_c(kc+1,j,i)+vx_c(kc+1,j,i-1) &
124  -vx_c(kc,j,i) -vx_c(kc,j,i-1))*0.5_dp
125  cvz5(kc) = (dzm_deta_g(j,i) &
126  +eaz_c_quotient(kc)*dh_c_deta_g(j,i)) &
127  *(vy_c(kc+1,j,i)+vy_c(kc+1,j-1,i) &
128  -vy_c(kc,j,i) -vy_c(kc,j-1,i))*0.5_dp
129 
130  do kc=1, kcmax-1
131  cvz3(kc) = avz3(kc)*h_c(j,i) &
132  *insq_g11_g(j,i)*insq_g22_g(j,i) &
133  *( (vx_c(kc,j,i)*sq_g22_sgx(j,i) &
134  -vx_c(kc,j,i-1)*sq_g22_sgx(j,i-1))*dxi_inv &
135  +(vy_c(kc,j,i)*sq_g11_sgy(j,i) &
136  -vy_c(kc,j-1,i)*sq_g11_sgy(j-1,i))*deta_inv )
137  cvz4(kc) = (dzm_dxi_g(j,i) &
138  +eaz_c_quotient(kc)*dh_c_dxi_g(j,i)) &
139  *(vx_c(kc+1,j,i)+vx_c(kc+1,j,i-1) &
140  -vx_c(kc-1,j,i)-vx_c(kc-1,j,i-1))*0.25_dp
141  cvz5(kc) = (dzm_deta_g(j,i) &
142  +eaz_c_quotient(kc)*dh_c_deta_g(j,i)) &
143  *(vy_c(kc+1,j,i)+vy_c(kc+1,j-1,i) &
144  -vy_c(kc-1,j,i)-vy_c(kc-1,j-1,i))*0.25_dp
145  end do
146 
147  kc=kcmax
148  cvz3(kc) = avz3(kc)*h_c(j,i) &
149  *insq_g11_g(j,i)*insq_g22_g(j,i) &
150  *( (vx_c(kc,j,i)*sq_g22_sgx(j,i) &
151  -vx_c(kc,j,i-1)*sq_g22_sgx(j,i-1))*dxi_inv &
152  +(vy_c(kc,j,i)*sq_g11_sgy(j,i) &
153  -vy_c(kc,j-1,i)*sq_g11_sgy(j-1,i))*deta_inv )
154  cvz4(kc) = (dzm_dxi_g(j,i) &
155  +eaz_c_quotient(kc)*dh_c_dxi_g(j,i)) &
156  *(vx_c(kc,j,i) +vx_c(kc,j,i-1) &
157  -vx_c(kc-1,j,i)-vx_c(kc-1,j,i-1))*0.5_dp
158  cvz5(kc) = (dzm_deta_g(j,i) &
159  +eaz_c_quotient(kc)*dh_c_deta_g(j,i)) &
160  *(vy_c(kc,j,i) +vy_c(kc,j-1,i) &
161  -vy_c(kc-1,j,i)-vy_c(kc-1,j-1,i))*0.5_dp
162 
163 !-------- Computation of vz_b --------
164 
165  vz_b(j,i) = cvz00
166 
167 !-------- Computation of vz --------
168 
169  if ((n_cts(j,i) == -1).or.(n_cts(j,i) == 0)) then
170  ! cold ice base, temperate ice base
171 
172  do kt=0, ktmax-1
173  vz_t(kt,j,i) = vz_b(j,i)
174  end do
175 
176  vz_m(j,i) = vz_b(j,i)
177 
178  vz_c(0,j,i) = vz_m(j,i) &
179  +0.5_dp*(-cvz3(0)+cvz4(0)+cvz5(0))
180 
181  do kc=1, kcmax-1
182  vz_c(kc,j,i) = vz_c(kc-1,j,i) &
183  -cvz3(kc)+cvz4(kc)+cvz5(kc)
184  end do
185 
186  else ! n_cts(j,i) == 1, temperate ice layer
187 
188  vz_t(0,j,i) = vz_b(j,i) &
189  +0.5_dp*(-cvz0(0)+cvz1(0)+cvz2(0))
190 
191  do kt=1, ktmax-1
192  vz_t(kt,j,i) = vz_t(kt-1,j,i) &
193  -cvz0(kt)+cvz1(kt)+cvz2(kt)
194  end do
195 
196  vz_m(j,i) = vz_t(ktmax-1,j,i) &
197  +0.5_dp*(-cvz0(ktmax)+cvz1(ktmax)+cvz2(ktmax))
198 
199  vz_c(0,j,i) = vz_m(j,i) &
200  +0.5_dp*(-cvz3(0)+cvz4(0)+cvz5(0))
201 
202  do kc=1, kcmax-1
203  vz_c(kc,j,i) = vz_c(kc-1,j,i) &
204  -cvz3(kc)+cvz4(kc)+cvz5(kc)
205  end do
206 
207  end if
208 
209 !-------- Computation of vz_s --------
210 
211  kc=kcmax
212 
213  vz_s(j,i) = vz_c(kc-1,j,i) &
214  +0.5_dp*(-cvz3(kc)+cvz4(kc)+cvz5(kc))
215 
216  else ! maske(j,i) == (1 or 2)
217 
218  vz_b(j,i) = 0.0_dp
219 
220  do kt=0, ktmax-1
221  vz_t(kt,j,i) = 0.0_dp
222  end do
223 
224  do kc=0, kcmax-1
225  vz_c(kc,j,i) = 0.0_dp
226  end do
227 
228  vz_s(j,i) = 0.0_dp
229 
230  end if
231 
232 end do
233 end do
234 
235 end subroutine calc_vz_sia
236 !
Generated on Thu Sep 12 2013 09:45:13 for SICOPOLIS V3.1 by   doxygen 1.8.2