SICOPOLIS V3.2
 All Classes Files Functions Variables Macros
calc_temp_enth_ssa.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : c a l c _ t e m p _ e n t h _ s s a
4 !
5 !> @file
6 !!
7 !! Computation of temperature and age for ice shelves (floating ice)
8 !! with the enthalpy method.
9 !!
10 !! @section Copyright
11 !!
12 !! Copyright 2013-2016 Ralf Greve, Heinz Blatter
13 !!
14 !! @section License
15 !!
16 !! This file is part of SICOPOLIS.
17 !!
18 !! SICOPOLIS is free software: you can redistribute it and/or modify
19 !! it under the terms of the GNU General Public License as published by
20 !! the Free Software Foundation, either version 3 of the License, or
21 !! (at your option) any later version.
22 !!
23 !! SICOPOLIS is distributed in the hope that it will be useful,
24 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
25 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26 !! GNU General Public License for more details.
27 !!
28 !! You should have received a copy of the GNU General Public License
29 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
30 !<
31 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
32 
33 !-------------------------------------------------------------------------------
34 !> Computation of temperature and age for ice shelves (floating ice)
35 !! with the enthalpy method.
36 !<------------------------------------------------------------------------------
37 subroutine calc_temp_enth_ssa(at1, at2_1, at2_2, at3_1, at3_2, &
38  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
39  ai1, ai2, &
40  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
41 
42 use sico_types
43 use sico_variables
44 use sico_vars
45 use ice_material_quantities, only : kappa_val, c_val
46 use sico_sle_solvers, only : tri_sle
47 use enth_temp_omega, only : enth_fct_temp_omega, temp_fct_enth, omega_fct_enth
48 
49 implicit none
50 
51 integer(i4b), intent(in) :: i, j
52 real(dp), intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
53  at3_1(0:kcmax), at3_2(0:kcmax), &
54  at4_1(0:kcmax), at4_2(0:kcmax), &
55  at5(0:kcmax), at6(0:kcmax), at7, &
56  ai1(0:kcmax), ai2(0:kcmax), &
57  atr1, alb1
58 real(dp), intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
59 
60 integer(i4b) :: kc, kt, kr
61 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
62  ce5(0:kcmax), ce6(0:kcmax), ce7(0:kcmax), ctr1, clb1
63 real(dp) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), ct3_sg(0:kcmax), &
64  ct4_sg(0:kcmax), adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
65 real(dp) :: ci1(0:kcmax), ci2(0:kcmax)
66 real(dp) :: temp_c_help(0:kcmax)
67 real(dp) :: vx_c_help, vy_c_help
68 real(dp) :: adv_vert_help
69 real(dp) :: dtt_dxi, dtt_deta
70 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
71  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
72  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
73  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
74  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
75 real(dp), parameter :: zero=0.0_dp
76 
77 !-------- Check for boundary points --------
78 
79 if ((i == 0).or.(i == imax).or.(j == 0).or.(j == jmax)) &
80  stop ' calc_temp_enth_ssa: Boundary points not allowed.'
81 
82 !-------- Abbreviations --------
83 
84 ctr1 = atr1
85 clb1 = alb1*q_geo(j,i)
86 
87 #if (ADV_VERT==1)
88 
89 do kc=1, kcmax-1
90  ct1(kc) = at1(kc)/h_c(j,i)*0.5_dp*(vz_c(kc,j,i)+vz_c(kc-1,j,i))
91 end do
92 
93 kc=0
94 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
95  ! only needed for kc=0 ...
96 kc=kcmax-1
97 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
98  ! ... and kc=KCMAX-1
99 
100 #elif (ADV_VERT==2 || ADV_VERT==3)
101 
102 do kc=0, kcmax-1
103  ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
104 end do
105 
106 #endif
107 
108 do kc=0, kcmax
109 
110  ct2(kc) = ( at2_1(kc)*dzm_dtau(j,i) &
111  +at2_2(kc)*dh_c_dtau(j,i) )/h_c(j,i)
112  ct3(kc) = ( at3_1(kc)*dzm_dxi_g(j,i) &
113  +at3_2(kc)*dh_c_dxi_g(j,i) )/h_c(j,i) &
114  *0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1)) *insq_g11_g(j,i)
115  ct4(kc) = ( at4_1(kc)*dzm_deta_g(j,i) &
116  +at4_2(kc)*dh_c_deta_g(j,i) )/h_c(j,i) &
117  *0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i)) *insq_g22_g(j,i)
118  ce5(kc) = at5(kc)/h_c(j,i)
119  ce7(kc) = 2.0_dp*at7 &
120  *viscosity(de_ssa(j,i), &
121  temp_c(kc,j,i), temp_c_m(kc,j,i), 0.0_dp, &
122  enh_c(kc,j,i), 0_i2b) &
123  *de_ssa(j,i)**2
124  ci1(kc) = ai1(kc)/h_c(j,i)
125 
126 end do
127 
128 #if (ADV_VERT==1)
129 
130 kc=0
131 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
132 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
133 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
134 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
135 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! only needed for kc=0 ...
136 kc=kcmax-1
137 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
138 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
139 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
140 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
141 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! ... and kc=KCMAX-1
142 
143 #elif (ADV_VERT==2 || ADV_VERT==3)
144 
145 do kc=0, kcmax-1
146  ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
147  ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
148  ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
149  adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
150  abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
151 end do
152 
153 #endif
154 
155 do kc=0, kcmax-1
156  temp_c_help(kc) = 0.5_dp*(temp_c(kc,j,i)+temp_c(kc+1,j,i))
157  ce6(kc) = at6(kc) &
158  *kappa_val(temp_c_help(kc))/(c_val(temp_c_help(kc))*h_c(j,i))
159  ci2(kc) = ai2(kc)/h_c(j,i)
160 end do
161 
162 #if (ADV_HOR==3)
163 dtt_dxi = 2.0_dp*dtt_2dxi
164 dtt_deta = 2.0_dp*dtt_2deta
165 #endif
166 
167 !-------- Set up the equations for the bedrock temperature --------
168 
169 kr=0
170 lgs_a1(kr) = 1.0_dp
171 lgs_a2(kr) = -1.0_dp
172 lgs_b(kr) = clb1
173 
174 #if (Q_LITHO==1)
175 ! (coupled heat-conducting bedrock)
176 
177 do kr=1, krmax-1
178  lgs_a0(kr) = - ctr1
179  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
180  lgs_a2(kr) = - ctr1
181  lgs_b(kr) = temp_r(kr,j,i)
182 end do
183 
184 #elif (Q_LITHO==0)
185 ! (no coupled heat-conducting bedrock)
186 
187 do kr=1, krmax-1
188  lgs_a0(kr) = 1.0_dp
189  lgs_a1(kr) = 0.0_dp
190  lgs_a2(kr) = -1.0_dp
191  lgs_b(kr) = 2.0_dp*clb1
192 end do
193 
194 #endif
195 
196 kr=krmax
197 lgs_a0(kr) = 0.0_dp
198 lgs_a1(kr) = 1.0_dp
199 lgs_b(kr) = temp_c_m(0,j,i)-delta_tm_sw
200 
201 !-------- Solve system of linear equations --------
202 
203 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
204 
205 !-------- Assign the result --------
206 
207 do kr=0, krmax
208  temp_r_neu(kr,j,i) = lgs_x(kr)
209 end do
210 
211 !-------- Set up the equations for the ice temperature --------
212 
213 kc=0
214 lgs_a1(kc) = 1.0_dp
215 lgs_a2(kc) = 0.0_dp
216 lgs_b(kc) = enth_fct_temp_omega(temp_c_m(kc,j,i)-delta_tm_sw, 0.0_dp)
217  ! zero water content assumed
218 
219 do kc=1, kcmax-1
220 
221 #if (ADV_VERT==1)
222 
223  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
224  -ce5(kc)*ce6(kc-1)
225  lgs_a1(kc) = 1.0_dp+ce5(kc)*(ce6(kc)+ce6(kc-1))
226  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
227  -ce5(kc)*ce6(kc)
228 
229 #elif (ADV_VERT==2)
230 
231  lgs_a0(kc) &
232  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
233  -ce5(kc)*ce6(kc-1)
234  lgs_a1(kc) &
235  = 1.0_dp &
236  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
237  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
238  +ce5(kc)*(ce6(kc)+ce6(kc-1))
239  lgs_a2(kc) &
240  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
241  -ce5(kc)*ce6(kc)
242 
243 #elif (ADV_VERT==3)
244 
245  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
246 
247  lgs_a0(kc) &
248  = -max(adv_vert_help, 0.0_dp) &
249  -ce5(kc)*ce6(kc-1)
250  lgs_a1(kc) &
251  = 1.0_dp &
252  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
253  +ce5(kc)*(ce6(kc)+ce6(kc-1))
254  lgs_a2(kc) &
255  = min(adv_vert_help, 0.0_dp) &
256  -ce5(kc)*ce6(kc)
257 
258 #endif
259 
260 #if (ADV_HOR==2)
261 
262  lgs_b(kc) = enth_c(kc,j,i) + ce7(kc) &
263  -dtt_2dxi* &
264  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
265  *(enth_c(kc,j,i+1)-enth_c(kc,j,i)) &
266  *insq_g11_sgx(j,i) &
267  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
268  *(enth_c(kc,j,i)-enth_c(kc,j,i-1)) &
269  *insq_g11_sgx(j,i-1) ) &
270  -dtt_2deta* &
271  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
272  *(enth_c(kc,j+1,i)-enth_c(kc,j,i)) &
273  *insq_g22_sgy(j,i) &
274  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
275  *(enth_c(kc,j,i)-enth_c(kc,j-1,i)) &
276  *insq_g22_sgy(j-1,i) )
277 
278 #elif (ADV_HOR==3)
279 
280  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
281  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
282 
283  lgs_b(kc) = enth_c(kc,j,i) + ce7(kc) &
284  -dtt_dxi* &
285  ( min(vx_c_help, 0.0_dp) &
286  *(enth_c(kc,j,i+1)-enth_c(kc,j,i)) &
287  *insq_g11_sgx(j,i) &
288  +max(vx_c_help, 0.0_dp) &
289  *(enth_c(kc,j,i)-enth_c(kc,j,i-1)) &
290  *insq_g11_sgx(j,i-1) ) &
291  -dtt_deta* &
292  ( min(vy_c_help, 0.0_dp) &
293  *(enth_c(kc,j+1,i)-enth_c(kc,j,i)) &
294  *insq_g22_sgy(j,i) &
295  +max(vy_c_help, 0.0_dp) &
296  *(enth_c(kc,j,i)-enth_c(kc,j-1,i)) &
297  *insq_g22_sgy(j-1,i) )
298 
299 #endif
300 
301 end do
302 
303 kc=kcmax
304 lgs_a0(kc) = 0.0_dp
305 lgs_a1(kc) = 1.0_dp
306 lgs_b(kc) = enth_fct_temp_omega(temp_s(j,i), 0.0_dp)
307  ! zero water content assumed
308 
309 !-------- Solve system of linear equations --------
310 
311 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
312 
313 !-------- Assign the result --------
314 
315 do kc=0, kcmax
316  enth_c_neu(kc,j,i) = lgs_x(kc)
317  temp_c_neu(kc,j,i) = temp_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
318  omega_c_neu(kc,j,i) = omega_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
319 end do
320 
321 !-------- Set enthalpy and water content in the redundant,
322 ! lower (kt) ice layer to the value at the ice base --------
323 
324 do kt=0, ktmax
325  enth_t_neu(kt,j,i) = enth_c_neu(0,j,i)
326  omega_t_neu(kt,j,i) = omega_c_neu(0,j,i)
327 end do
328 
329 !-------- Water drainage from the non-existing temperate ice --------
330 
331 q_tld(j,i) = 0.0_dp
332 
333 !-------- Set up the equations for the age of cold ice --------
334 
335 kc=0 ! adv_vert_sg(0) <= 0
336 lgs_a1(kc) = 1.0_dp - min(adv_vert_sg(kc), 0.0_dp) ! (directed downward)
337 lgs_a2(kc) = min(adv_vert_sg(kc), 0.0_dp) ! assumed/enforced
338 
339 #if (ADV_HOR==2)
340 
341 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
342  -dtt_2dxi* &
343  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
344  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
345  *insq_g11_sgx(j,i) &
346  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
347  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
348  *insq_g11_sgx(j,i-1) ) &
349  -dtt_2deta* &
350  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
351  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
352  *insq_g22_sgy(j,i) &
353  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
354  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
355  *insq_g22_sgy(j-1,i) )
356 
357 #elif (ADV_HOR==3)
358 
359 vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
360 vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
361 
362 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
363  -dtt_dxi* &
364  ( min(vx_c_help, 0.0_dp) &
365  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
366  *insq_g11_sgx(j,i) &
367  +max(vx_c_help, 0.0_dp) &
368  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
369  *insq_g11_sgx(j,i-1) ) &
370  -dtt_deta* &
371  ( min(vy_c_help, 0.0_dp) &
372  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
373  *insq_g22_sgy(j,i) &
374  +max(vy_c_help, 0.0_dp) &
375  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
376  *insq_g22_sgy(j-1,i) )
377 
378 #endif
379 
380 do kc=1, kcmax-1
381 
382 #if (ADV_VERT==1)
383 
384  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
385  -ci1(kc)*ci2(kc-1)
386  lgs_a1(kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
387  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
388  -ci1(kc)*ci2(kc)
389 
390 #elif (ADV_VERT==2)
391 
392  lgs_a0(kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
393  lgs_a1(kc) = 1.0_dp &
394  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
395  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
396  lgs_a2(kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
397 
398 #elif (ADV_VERT==3)
399 
400  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
401 
402  lgs_a0(kc) = -max(adv_vert_help, 0.0_dp)
403  lgs_a1(kc) = 1.0_dp &
404  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp)
405  lgs_a2(kc) = min(adv_vert_help, 0.0_dp)
406 
407 #endif
408 
409 #if (ADV_HOR==2)
410 
411  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
412  -dtt_2dxi* &
413  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
414  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
415  *insq_g11_sgx(j,i) &
416  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
417  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
418  *insq_g11_sgx(j,i-1) ) &
419  -dtt_2deta* &
420  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
421  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
422  *insq_g22_sgy(j,i) &
423  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
424  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
425  *insq_g22_sgy(j-1,i) )
426 
427 #elif (ADV_HOR==3)
428 
429  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
430  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
431 
432  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
433  -dtt_dxi* &
434  ( min(vx_c_help, 0.0_dp) &
435  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
436  *insq_g11_sgx(j,i) &
437  +max(vx_c_help, 0.0_dp) &
438  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
439  *insq_g11_sgx(j,i-1) ) &
440  -dtt_deta* &
441  ( min(vy_c_help, 0.0_dp) &
442  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
443  *insq_g22_sgy(j,i) &
444  +max(vy_c_help, 0.0_dp) &
445  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
446  *insq_g22_sgy(j-1,i) )
447 
448 #endif
449 
450 end do
451 
452 kc=kcmax
453 if (as_perp(j,i) >= zero) then
454  lgs_a0(kc) = 0.0_dp
455  lgs_a1(kc) = 1.0_dp
456  lgs_b(kc) = 0.0_dp
457 else
458  lgs_a0(kc) = -max(adv_vert_sg(kc-1), 0.0_dp)
459  lgs_a1(kc) = 1.0_dp + max(adv_vert_sg(kc-1), 0.0_dp)
460  ! adv_vert_sg(KCMAX-1) >= 0 (directed upward)
461  ! assumed/enforced
462 #if (ADV_HOR==2)
463 
464  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
465  -dtt_2dxi* &
466  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
467  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
468  *insq_g11_sgx(j,i) &
469  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
470  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
471  *insq_g11_sgx(j,i-1) ) &
472  -dtt_2deta* &
473  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
474  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
475  *insq_g22_sgy(j,i) &
476  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
477  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
478  *insq_g22_sgy(j-1,i) )
479 
480 #elif (ADV_HOR==3)
481 
482  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
483  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
484 
485  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
486  -dtt_dxi* &
487  ( min(vx_c_help, 0.0_dp) &
488  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
489  *insq_g11_sgx(j,i) &
490  +max(vx_c_help, 0.0_dp) &
491  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
492  *insq_g11_sgx(j,i-1) ) &
493  -dtt_deta* &
494  ( min(vy_c_help, 0.0_dp) &
495  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
496  *insq_g22_sgy(j,i) &
497  +max(vy_c_help, 0.0_dp) &
498  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
499  *insq_g22_sgy(j-1,i) )
500 
501 #endif
502 
503 end if
504 
505 !-------- Solve system of linear equations --------
506 
507 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
508 
509 !-------- Assign the result,
510 ! restriction to interval [0, AGE_MAX yr] --------
511 
512 do kc=0, kcmax
513 
514  age_c_neu(kc,j,i) = lgs_x(kc)
515 
516  if (age_c_neu(kc,j,i) < (age_min*year_sec)) &
517  age_c_neu(kc,j,i) = 0.0_dp
518  if (age_c_neu(kc,j,i) > (age_max*year_sec)) &
519  age_c_neu(kc,j,i) = age_max*year_sec
520 
521 end do
522 
523 !-------- Age of the ice in the redundant, lower (kt) ice layer --------
524 
525 do kt=0, ktmax
526  age_t_neu(kt,j,i) = age_c_neu(0,j,i)
527 end do
528 
529 end subroutine calc_temp_enth_ssa
530 !
sico_sle_solvers::tri_sle
subroutine tri_sle(a0, a1, a2, x, b, nrows)
Solution of a system of linear equations Ax=b with tridiagonal matrix A.
Definition: sico_sle_solvers.F90:122
sico_types
Declarations of kind types for SICOPOLIS.
Definition: sico_types.F90:35
viscosity
real(dp) function viscosity(de_val, temp_val, temp_m_val, omega_val, enh_val, i_flag_cold_temp)
Ice viscosity as a function of the effective strain rate and the temperature (in cold ice) or the wat...
Definition: viscosity.F90:39
enth_temp_omega::temp_fct_enth
real(dp) function, public temp_fct_enth(enth_val, temp_m_val)
Temperature as a function of enthalpy.
Definition: enth_temp_omega.F90:274
sico_vars
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars.F90:35
enth_temp_omega::omega_fct_enth
real(dp) function, public omega_fct_enth(enth_val, temp_m_val)
Water content as a function of enthalpy.
Definition: enth_temp_omega.F90:298
sico_sle_solvers
Solvers for systems of linear equations used by SICOPOLIS.
Definition: sico_sle_solvers.F90:35
ice_material_quantities::kappa_val
real(dp) function, public kappa_val(temp_val)
Heat conductivity of ice: Linear interpolation of tabulated values in KAPPA(.).
Definition: ice_material_quantities.F90:217
ice_material_quantities::c_val
real(dp) function, public c_val(temp_val)
Specific heat of ice: Linear interpolation of tabulated values in C(.).
Definition: ice_material_quantities.F90:254
ice_material_quantities
Material quantities of ice: Rate factor, heat conductivity, specific heat (heat capacity).
Definition: ice_material_quantities.F90:37
calc_temp_enth_ssa
subroutine calc_temp_enth_ssa(at1, at2_1, at2_2, at3_1, at3_2, at4_1, at4_2, at5, at6, at7, atr1, alb1, ai1, ai2, dtime_temp, dtt_2dxi, dtt_2deta, i, j)
Computation of temperature and age for ice shelves (floating ice) with the enthalpy method...
Definition: calc_temp_enth_ssa.F90:37
sico_variables
Declarations of global variables for SICOPOLIS.
Definition: sico_variables.F90:35
enth_temp_omega
Conversion from temperature (temp) and water content (omega) to enthalpy (enth) and vice versa...
Definition: enth_temp_omega.F90:37
enth_temp_omega::enth_fct_temp_omega
real(dp) function, public enth_fct_temp_omega(temp_val, omega_val)
Enthalpy as a function of temperature and water content.
Definition: enth_temp_omega.F90:259
Generated on Wed Mar 2 2016 11:26:54 for SICOPOLIS V3.2 by   doxygen 1.8.5