SICOPOLIS V3.1
 All Classes Files Functions Variables Macros
calc_temp_ssa.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : c a l c _ t e m p _ s s a
4 !
5 !> @file
6 !!
7 !! Computation of temperature and age for ice shelves (floating ice).
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2013 Ralf Greve, Tatsuru Sato
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 temperature and age for ice shelves (floating ice).
34 !<------------------------------------------------------------------------------
35 subroutine calc_temp_ssa(at1, at2_1, at2_2, at3_1, at3_2, &
36  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
37  ai1, ai2, ai4, mean_accum_inv, &
38  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
39 
40 use sico_types
41 use sico_variables
42 use sico_sle_solvers, only : tri_sle
43 
44 implicit none
45 
46 integer(i4b), intent(in) :: i, j
47 real(dp), intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
48  at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
49  at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
50  ai1(0:kcmax), ai2(0:kcmax), ai4, &
51  atr1, alb1
52 real(dp), intent(in) :: mean_accum_inv
53 real(dp), intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
54 
55 integer(i4b) :: kc, kt, kr
56 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
57  ct5(0:kcmax), ct6(0:kcmax), ct7(0:kcmax), ctr1, clb1
58 real(dp) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), ct3_sg(0:kcmax), &
59  ct4_sg(0:kcmax), adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
60 real(dp) :: ci1(0:kcmax), ci2(0:kcmax)
61 real(dp) :: ftx_c_l(0:kcmax), ftx_c_r(0:kcmax), &
62  fty_c_l(0:kcmax), fty_c_r(0:kcmax), &
63  fax_c_l(0:kcmax), fax_c_r(0:kcmax), &
64  fay_c_l(0:kcmax), fay_c_r(0:kcmax)
65 real(dp) :: temp_c_help(0:kcmax)
66 real(dp) :: vx_c_help, vy_c_help
67 real(dp) :: adv_vert_help
68 real(dp) :: dtt_dxi, dtt_deta
69 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
70  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
71  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
72  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
73  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
74 real(dp), parameter :: zero=0.0_dp
75 
76 !-------- Check for boundary points --------
77 
78 if ((i == 0).or.(i == imax).or.(j == 0).or.(j == jmax)) &
79  stop ' calc_temp_ssa: Boundary points not allowed.'
80 
81 !-------- Abbreviations --------
82 
83 ctr1 = atr1
84 clb1 = alb1*q_geo(j,i)
85 
86 #if ADV_VERT==1
87 
88 do kc=1, kcmax-1
89  ct1(kc) = at1(kc)/h_c(j,i)*0.5_dp*(vz_c(kc,j,i)+vz_c(kc-1,j,i))
90 end do
91 
92 kc=0
93 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
94  ! only needed for kc=0 ...
95 kc=kcmax-1
96 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
97  ! ... and kc=KCMAX-1
98 
99 #elif ( ADV_VERT==2 || ADV_VERT==3 )
100 
101 do kc=0, kcmax-1
102  ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
103 end do
104 
105 #endif
106 
107 do kc=0, kcmax
108 
109  ct2(kc) = ( at2_1(kc)*dzm_dtau(j,i) &
110  +at2_2(kc)*dh_c_dtau(j,i) )/h_c(j,i)
111  ct3(kc) = ( at3_1(kc)*dzm_dxi_g(j,i) &
112  +at3_2(kc)*dh_c_dxi_g(j,i) )/h_c(j,i) &
113  *0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1)) *insq_g11_g(j,i)
114  ct4(kc) = ( at4_1(kc)*dzm_deta_g(j,i) &
115  +at4_2(kc)*dh_c_deta_g(j,i) )/h_c(j,i) &
116  *0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i)) *insq_g22_g(j,i)
117  ct5(kc) = at5(kc) &
118  /c_val(temp_c(kc,j,i)) &
119  /h_c(j,i)
120  ct7(kc) = 2.0_dp*at7 &
121  /c_val(temp_c(kc,j,i)) &
122  *viscosity(de_ssa(j,i), &
123  temp_c(kc,j,i), temp_c_m(kc,j,i), 0.0_dp, &
124  enh_c(kc,j,i), .true.) &
125  *de_ssa(j,i)**2
126  ci1(kc) = ai1(kc)/h_c(j,i)
127 
128 ! ------ Fluxes for horizontal advection of temperature and age
129 
130 #if ADV_HOR==4
131  if (vx_c(kc,j,i-1) >= zero) then
132  ftx_c_l(kc) = temp_c(kc,j,i-1)*vx_c(kc,j,i-1)
133  fax_c_l(kc) = age_c(kc,j,i-1)*vx_c(kc,j,i-1)
134  else
135  ftx_c_l(kc) = temp_c(kc,j,i)*vx_c(kc,j,i-1)
136  fax_c_l(kc) = age_c(kc,j,i)*vx_c(kc,j,i-1)
137  end if
138 
139  if (vx_c(kc,j,i) >= zero) then
140  ftx_c_r(kc) = temp_c(kc,j,i)*vx_c(kc,j,i)
141  fax_c_r(kc) = age_c(kc,j,i)*vx_c(kc,j,i)
142  else
143  ftx_c_r(kc) = temp_c(kc,j,i+1)*vx_c(kc,j,i)
144  fax_c_r(kc) = age_c(kc,j,i+1)*vx_c(kc,j,i)
145  end if
146 
147  if (vy_c(kc,j-1,i) >= zero) then
148  fty_c_l(kc) = temp_c(kc,j-1,i)*vy_c(kc,j-1,i)
149  fay_c_l(kc) = age_c(kc,j-1,i)*vy_c(kc,j-1,i)
150  else
151  fty_c_l(kc) = temp_c(kc,j,i)*vy_c(kc,j-1,i)
152  fay_c_l(kc) = age_c(kc,j,i)*vy_c(kc,j-1,i)
153  end if
154 
155  if (vy_c(kc,j,i) >= zero) then
156  fty_c_r(kc) = temp_c(kc,j,i)*vy_c(kc,j,i)
157  fay_c_r(kc) = age_c(kc,j,i)*vy_c(kc,j,i)
158  else
159  fty_c_r(kc) = temp_c(kc,j+1,i)*vy_c(kc,j,i)
160  fay_c_r(kc) = age_c(kc,j+1,i)*vy_c(kc,j,i)
161  end if
162 #endif
163 
164 end do
165 
166 #if ADV_VERT==1
167 
168 kc=0
169 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
170 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
171 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
172 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
173 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! only needed for kc=0 ...
174 kc=kcmax-1
175 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
176 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
177 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
178 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
179 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! ... and kc=KCMAX-1
180 
181 #elif ( ADV_VERT==2 || ADV_VERT==3 )
182 
183 do kc=0, kcmax-1
184  ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
185  ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
186  ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
187  adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
188  abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
189 end do
190 
191 #endif
192 
193 do kc=0, kcmax-1
194  temp_c_help(kc) = 0.5_dp*(temp_c(kc,j,i)+temp_c(kc+1,j,i))
195  ct6(kc) = at6(kc) &
196  *kappa_val(temp_c_help(kc)) &
197  /h_c(j,i)
198  ci2(kc) = ai2(kc)/h_c(j,i)
199 end do
200 
201 #if ( ADV_HOR==3 || ADV_HOR==4 )
202 dtt_dxi = 2.0_dp*dtt_2dxi
203 dtt_deta = 2.0_dp*dtt_2deta
204 #endif
205 
206 !-------- Set up the equations for the bedrock temperature --------
207 
208 kr=0
209 lgs_a1(kr) = 1.0_dp
210 lgs_a2(kr) = -1.0_dp
211 lgs_b(kr) = clb1
212 
213 #if Q_LITHO==1
214 ! (coupled heat-conducting bedrock)
215 
216 do kr=1, krmax-1
217  lgs_a0(kr) = - ctr1
218  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
219  lgs_a2(kr) = - ctr1
220  lgs_b(kr) = temp_r(kr,j,i)
221 end do
222 
223 #elif Q_LITHO==0
224 ! (no coupled heat-conducting bedrock)
225 
226 do kr=1, krmax-1
227  lgs_a0(kr) = 1.0_dp
228  lgs_a1(kr) = 0.0_dp
229  lgs_a2(kr) = -1.0_dp
230  lgs_b(kr) = 2.0_dp*clb1
231 end do
232 
233 #endif
234 
235 kr=krmax
236 lgs_a0(kr) = 0.0_dp
237 lgs_a1(kr) = 1.0_dp
238 lgs_b(kr) = temp_c_m(0,j,i)-delta_tm_sw
239 
240 !-------- Solve system of linear equations --------
241 
242 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
243 
244 !-------- Assign the result --------
245 
246 do kr=0, krmax
247  temp_r_neu(kr,j,i) = lgs_x(kr)
248 end do
249 
250 !-------- Set up the equations for the ice temperature --------
251 
252 kc=0
253 lgs_a1(kc) = 1.0_dp
254 lgs_a2(kc) = 0.0_dp
255 lgs_b(kc) = temp_c_m(0,j,i)-delta_tm_sw
256 
257 do kc=1, kcmax-1
258 
259 #if ADV_VERT==1
260 
261  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
262  -ct5(kc)*ct6(kc-1)
263  lgs_a1(kc) = 1.0_dp+ct5(kc)*(ct6(kc)+ct6(kc-1))
264  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
265  -ct5(kc)*ct6(kc)
266 
267 #elif ADV_VERT==2
268 
269  lgs_a0(kc) &
270  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
271  -ct5(kc)*ct6(kc-1)
272  lgs_a1(kc) &
273  = 1.0_dp &
274  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
275  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
276  +ct5(kc)*(ct6(kc)+ct6(kc-1))
277  lgs_a2(kc) &
278  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
279  -ct5(kc)*ct6(kc)
280 
281 #elif ADV_VERT==3
282 
283  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
284 
285  lgs_a0(kc) &
286  = -max(adv_vert_help, 0.0_dp) &
287  -ct5(kc)*ct6(kc-1)
288  lgs_a1(kc) &
289  = 1.0_dp &
290  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
291  +ct5(kc)*(ct6(kc)+ct6(kc-1))
292  lgs_a2(kc) &
293  = min(adv_vert_help, 0.0_dp) &
294  -ct5(kc)*ct6(kc)
295 
296 #endif
297 
298 #if ADV_HOR==2
299 
300  lgs_b(kc) = temp_c(kc,j,i) + ct7(kc) &
301  -dtt_2dxi* &
302  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
303  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
304  *insq_g11_sgx(j,i) &
305  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
306  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
307  *insq_g11_sgx(j,i-1) ) &
308  -dtt_2deta* &
309  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
310  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
311  *insq_g22_sgy(j,i) &
312  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
313  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
314  *insq_g22_sgy(j-1,i) )
315 
316 #elif ADV_HOR==3
317 
318  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
319  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
320 
321  lgs_b(kc) = temp_c(kc,j,i) + ct7(kc) &
322  -dtt_dxi* &
323  ( min(vx_c_help, 0.0_dp) &
324  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
325  *insq_g11_sgx(j,i) &
326  +max(vx_c_help, 0.0_dp) &
327  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
328  *insq_g11_sgx(j,i-1) ) &
329  -dtt_deta* &
330  ( min(vy_c_help, 0.0_dp) &
331  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
332  *insq_g22_sgy(j,i) &
333  +max(vy_c_help, 0.0_dp) &
334  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
335  *insq_g22_sgy(j-1,i) )
336 
337 #elif ADV_HOR==4
338 
339  lgs_b(kc) = temp_c(kc,j,i) + ct7(kc) &
340  -dtt_dxi *(ftx_c_r(kc)-ftx_c_l(kc)) &
341  -dtt_deta*(fty_c_r(kc)-fty_c_l(kc))
342 
343 #endif
344 
345 end do
346 
347 kc=kcmax
348 lgs_a0(kc) = 0.0_dp
349 lgs_a1(kc) = 1.0_dp
350 lgs_b(kc) = temp_s(j,i)
351 
352 !-------- Solve system of linear equations --------
353 
354 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
355 
356 !-------- Assign the result --------
357 
358 do kc=0, kcmax
359  temp_c_neu(kc,j,i) = lgs_x(kc)
360 end do
361 
362 !-------- Set water content in the non-existing temperate layer
363 ! to zero --------
364 
365 do kt=0, ktmax
366  omega_t_neu(kt,j,i) = 0.0_dp
367 end do
368 
369 !-------- Water drainage from the non-existing temperate layer --------
370 
371 q_tld(j,i) = 0.0_dp
372 
373 !-------- Set up the equations for the age of cold ice --------
374 
375 kc=0 ! adv_vert_sg(0) <= 0
376 lgs_a1(kc) = 1.0_dp - min(adv_vert_sg(kc), 0.0_dp) ! (directed downward)
377 lgs_a2(kc) = min(adv_vert_sg(kc), 0.0_dp) ! assumed/enforced
378 
379 #if ADV_HOR==2
380 
381 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
382  -dtt_2dxi* &
383  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
384  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
385  *insq_g11_sgx(j,i) &
386  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
387  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
388  *insq_g11_sgx(j,i-1) ) &
389  -dtt_2deta* &
390  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
391  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
392  *insq_g22_sgy(j,i) &
393  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
394  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
395  *insq_g22_sgy(j-1,i) )
396 
397 #elif ADV_HOR==3
398 
399 vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
400 vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
401 
402 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
403  -dtt_dxi* &
404  ( min(vx_c_help, 0.0_dp) &
405  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
406  *insq_g11_sgx(j,i) &
407  +max(vx_c_help, 0.0_dp) &
408  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
409  *insq_g11_sgx(j,i-1) ) &
410  -dtt_deta* &
411  ( min(vy_c_help, 0.0_dp) &
412  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
413  *insq_g22_sgy(j,i) &
414  +max(vy_c_help, 0.0_dp) &
415  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
416  *insq_g22_sgy(j-1,i) )
417 
418 #elif ADV_HOR==4
419 
420 lgs_b(kc) = ...
421 
422 #endif
423 
424 do kc=1, kcmax-1
425 
426 #if ADV_VERT==1
427 
428  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
429  -ci1(kc)*ci2(kc-1)
430  lgs_a1(kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
431  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
432  -ci1(kc)*ci2(kc)
433 
434 #elif ADV_VERT==2
435 
436  lgs_a0(kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
437  lgs_a1(kc) = 1.0_dp &
438  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
439  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
440  lgs_a2(kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
441 
442 #elif ADV_VERT==3
443 
444  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
445 
446  lgs_a0(kc) = -max(adv_vert_help, 0.0_dp)
447  lgs_a1(kc) = 1.0_dp &
448  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp)
449  lgs_a2(kc) = min(adv_vert_help, 0.0_dp)
450 
451 #endif
452 
453 #if ADV_HOR==2
454 
455  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
456  -dtt_2dxi* &
457  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
458  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
459  *insq_g11_sgx(j,i) &
460  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
461  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
462  *insq_g11_sgx(j,i-1) ) &
463  -dtt_2deta* &
464  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
465  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
466  *insq_g22_sgy(j,i) &
467  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
468  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
469  *insq_g22_sgy(j-1,i) )
470 
471 #elif ADV_HOR==3
472 
473  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
474  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
475 
476  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
477  -dtt_dxi* &
478  ( min(vx_c_help, 0.0_dp) &
479  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
480  *insq_g11_sgx(j,i) &
481  +max(vx_c_help, 0.0_dp) &
482  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
483  *insq_g11_sgx(j,i-1) ) &
484  -dtt_deta* &
485  ( min(vy_c_help, 0.0_dp) &
486  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
487  *insq_g22_sgy(j,i) &
488  +max(vy_c_help, 0.0_dp) &
489  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
490  *insq_g22_sgy(j-1,i) )
491 
492 #elif ADV_HOR==4
493 
494  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
495  -dtt_dxi *(fax_c_r(kc)-fax_c_l(kc)) &
496  -dtt_deta*(fay_c_r(kc)-fay_c_l(kc))
497 
498 #endif
499 
500 end do
501 
502 kc=kcmax
503 if (as_perp(j,i) >= zero) then
504  lgs_a0(kc) = 0.0_dp
505  lgs_a1(kc) = 1.0_dp
506  lgs_b(kc) = 0.0_dp
507 else
508  lgs_a0(kc) = -max(adv_vert_sg(kc-1), 0.0_dp)
509  lgs_a1(kc) = 1.0_dp + max(adv_vert_sg(kc-1), 0.0_dp)
510  ! adv_vert_sg(KCMAX-1) >= 0 (directed upward)
511  ! assumed/enforced
512 #if ADV_HOR==2
513 
514  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
515  -dtt_2dxi* &
516  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
517  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
518  *insq_g11_sgx(j,i) &
519  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
520  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
521  *insq_g11_sgx(j,i-1) ) &
522  -dtt_2deta* &
523  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
524  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
525  *insq_g22_sgy(j,i) &
526  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
527  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
528  *insq_g22_sgy(j-1,i) )
529 
530 #elif ADV_HOR==3
531 
532  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
533  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
534 
535  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
536  -dtt_dxi* &
537  ( min(vx_c_help, 0.0_dp) &
538  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
539  *insq_g11_sgx(j,i) &
540  +max(vx_c_help, 0.0_dp) &
541  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
542  *insq_g11_sgx(j,i-1) ) &
543  -dtt_deta* &
544  ( min(vy_c_help, 0.0_dp) &
545  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
546  *insq_g22_sgy(j,i) &
547  +max(vy_c_help, 0.0_dp) &
548  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
549  *insq_g22_sgy(j-1,i) )
550 
551 #elif ADV_HOR==4
552 
553  lgs_b(kc) = ...
554 
555 #endif
556 
557 end if
558 
559 !-------- Solve system of linear equations --------
560 
561 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
562 
563 !-------- Assign the result,
564 ! restriction to interval [0, AGE_MAX yr] --------
565 
566 do kc=0, kcmax
567 
568  age_c_neu(kc,j,i) = lgs_x(kc)
569 
570  if (age_c_neu(kc,j,i) < (age_min*year_sec)) &
571  age_c_neu(kc,j,i) = 0.0_dp
572  if (age_c_neu(kc,j,i) > (age_max*year_sec)) &
573  age_c_neu(kc,j,i) = age_max*year_sec
574 
575 end do
576 
577 !-------- Age of the ice in the non-existing temperate layer --------
578 
579 do kt=0, ktmax
580  age_t_neu(kt,j,i) = age_c_neu(0,j,i)
581 end do
582 
583 end subroutine calc_temp_ssa
584 !
Generated on Thu Sep 12 2013 09:45:11 for SICOPOLIS V3.1 by   doxygen 1.8.2