SICOPOLIS V3.1
 All Classes Files Functions Variables Macros
calc_temp3.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : c a l c _ t e m p 3
4 !
5 !> @file
6 !!
7 !! Computation of temperature, water content and age for an ice column with a
8 !! temperate base overlain by a temperate-ice layer.
9 !!
10 !! @section Copyright
11 !!
12 !! Copyright 2009-2013 Ralf Greve
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, water content and age for an ice column with a
35 !! temperate base overlain by a temperate-ice layer.
36 !<------------------------------------------------------------------------------
37 subroutine calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
38  at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
39  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
40  ai1, ai2, ai3, mean_accum_inv, dzeta_t, &
41  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
42 
43 use sico_types
44 use sico_variables
45 use sico_sle_solvers, only : tri_sle
46 
47 implicit none
48 
49 integer(i4b), intent(in) :: i, j
50 real(dp), intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
51  at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
52  at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
53  ai1(0:kcmax), ai2(0:kcmax), ai3, &
54  atr1, am1, am2, alb1
55 real(dp), intent(in) :: aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld
56 real(dp), intent(in) :: mean_accum_inv
57 real(dp), intent(in) :: dzeta_t, dtime_temp, dtt_2dxi, dtt_2deta
58 
59 integer(i4b) :: kc, kt, kr
60 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
61  ct5(0:kcmax), ct6(0:kcmax), ct7(0:kcmax), ctr1, cm1, cm2, &
62  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), ci3
66 real(dp) :: cw1(0:ktmax), cw2(0:ktmax), cw3(0:ktmax), cw4(0:ktmax), &
67  cw5, cw7(0:ktmax), cw8, cw9(0:ktmax)
68 real(dp) :: cw1_sg(0:ktmax), cw2_sg(0:ktmax), cw3_sg(0:ktmax), &
69  cw4_sg(0:ktmax), adv_vert_w_sg(0:ktmax), abs_adv_vert_w_sg(0:ktmax)
70 real(dp) :: ftx_c_l(0:kcmax), ftx_c_r(0:kcmax), &
71  fty_c_l(0:kcmax), fty_c_r(0:kcmax), &
72  fax_c_l(0:kcmax), fax_c_r(0:kcmax), &
73  fay_c_l(0:kcmax), fay_c_r(0:kcmax)
74 real(dp) :: fwx_t_l(0:ktmax), fwx_t_r(0:ktmax), &
75  fwy_t_l(0:ktmax), fwy_t_r(0:ktmax), &
76  fax_t_l(0:ktmax), fax_t_r(0:ktmax), &
77  fay_t_l(0:ktmax), fay_t_r(0:ktmax)
78 real(dp) :: sigma_c_help(0:kcmax), sigma_t_help(0:ktmax), &
79  temp_c_help(0:kcmax)
80 real(dp) :: vx_c_help, vy_c_help, vx_t_help, vy_t_help
81 real(dp) :: adv_vert_help, adv_vert_w_help
82 real(dp) :: dtt_dxi, dtt_deta
83 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
84  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
85  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
86  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
87  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
88 real(dp), parameter :: zero=0.0_dp
89 
90 !-------- Check for boundary points --------
91 
92 if ((i == 0).or.(i == imax).or.(j == 0).or.(j == jmax)) &
93  stop ' calc_temp3: Boundary points not allowed.'
94 
95 !-------- Abbreviations --------
96 
97 ctr1 = atr1
98 cm1 = am1*h_c_neu(j,i)
99 clb1 = alb1*q_geo(j,i)
100 
101 #if ADV_VERT==1
102 
103 do kc=1, kcmax-1
104  ct1(kc) = at1(kc)/h_c_neu(j,i)*0.5_dp*(vz_c(kc,j,i)+vz_c(kc-1,j,i))
105 end do
106 
107 kc=0
108 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c_neu(j,i)*vz_c(kc,j,i)
109  ! only needed for kc=0 ...
110 kc=kcmax-1
111 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c_neu(j,i)*vz_c(kc,j,i)
112  ! ... and kc=KCMAX-1
113 
114 #elif ( ADV_VERT==2 || ADV_VERT==3 )
115 
116 do kc=0, kcmax-1
117  ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c_neu(j,i)*vz_c(kc,j,i)
118 end do
119 
120 #endif
121 
122 do kc=0, kcmax
123 
124  ct2(kc) = ( at2_1(kc)*dzm_dtau(j,i) &
125  +at2_2(kc)*dh_c_dtau(j,i) )/h_c_neu(j,i)
126  ct3(kc) = ( at3_1(kc)*dzm_dxi_g(j,i) &
127  +at3_2(kc)*dh_c_dxi_g(j,i) )/h_c_neu(j,i) &
128  *0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1)) *insq_g11_g(j,i)
129  ct4(kc) = ( at4_1(kc)*dzm_deta_g(j,i) &
130  +at4_2(kc)*dh_c_deta_g(j,i) )/h_c_neu(j,i) &
131  *0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i)) *insq_g22_g(j,i)
132  ct5(kc) = at5(kc) &
133  /c_val(temp_c(kc,j,i)) &
134  /h_c_neu(j,i)
135 
136  sigma_c_help(kc) &
137  = rho*g*h_c_neu(j,i)*(1.0_dp-eaz_c_quotient(kc)) &
138  *(dzs_dxi_g(j,i)**2+dzs_deta_g(j,i)**2)**0.5_dp
139  ct7(kc) = at7 &
140  /c_val(temp_c(kc,j,i)) &
141  *enh_c(kc,j,i) &
142  *ratefac(temp_c(kc,j,i), temp_c_m(kc,j,i)) &
143  *creep(sigma_c_help(kc)) &
144  *sigma_c_help(kc)*sigma_c_help(kc)
145 
146  ci1(kc) = ai1(kc)/h_c_neu(j,i)
147 
148 ! ------ Fluxes for horizontal advection of temperature and age
149 
150 #if ADV_HOR==4
151  if (vx_c(kc,j,i-1) >= zero) then
152  ftx_c_l(kc) = temp_c(kc,j,i-1)*vx_c(kc,j,i-1)
153  fax_c_l(kc) = age_c(kc,j,i-1)*vx_c(kc,j,i-1)
154  else
155  ftx_c_l(kc) = temp_c(kc,j,i)*vx_c(kc,j,i-1)
156  fax_c_l(kc) = age_c(kc,j,i)*vx_c(kc,j,i-1)
157  end if
158 
159  if (vx_c(kc,j,i) >= zero) then
160  ftx_c_r(kc) = temp_c(kc,j,i)*vx_c(kc,j,i)
161  fax_c_r(kc) = age_c(kc,j,i)*vx_c(kc,j,i)
162  else
163  ftx_c_r(kc) = temp_c(kc,j,i+1)*vx_c(kc,j,i)
164  fax_c_r(kc) = age_c(kc,j,i+1)*vx_c(kc,j,i)
165  end if
166 
167  if (vy_c(kc,j-1,i) >= zero) then
168  fty_c_l(kc) = temp_c(kc,j-1,i)*vy_c(kc,j-1,i)
169  fay_c_l(kc) = age_c(kc,j-1,i)*vy_c(kc,j-1,i)
170  else
171  fty_c_l(kc) = temp_c(kc,j,i)*vy_c(kc,j-1,i)
172  fay_c_l(kc) = age_c(kc,j,i)*vy_c(kc,j-1,i)
173  end if
174 
175  if (vy_c(kc,j,i) >= zero) then
176  fty_c_r(kc) = temp_c(kc,j,i)*vy_c(kc,j,i)
177  fay_c_r(kc) = age_c(kc,j,i)*vy_c(kc,j,i)
178  else
179  fty_c_r(kc) = temp_c(kc,j+1,i)*vy_c(kc,j,i)
180  fay_c_r(kc) = age_c(kc,j+1,i)*vy_c(kc,j,i)
181  end if
182 #endif
183 
184 end do
185 
186 #if ADV_VERT==1
187 
188 kc=0
189 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
190 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
191 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
192 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
193 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! only needed for kc=0 ...
194 kc=kcmax-1
195 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
196 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
197 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
198 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
199 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! ... and kc=KCMAX-1
200 
201 #elif ( ADV_VERT==2 || ADV_VERT==3 )
202 
203 do kc=0, kcmax-1
204  ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
205  ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
206  ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
207  adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
208  abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
209 end do
210 
211 #endif
212 
213 do kc=0, kcmax-1
214  temp_c_help(kc) = 0.5_dp*(temp_c(kc,j,i)+temp_c(kc+1,j,i))
215  ct6(kc) = at6(kc) &
216  *kappa_val(temp_c_help(kc)) &
217  /h_c_neu(j,i)
218  ci2(kc) = ai2(kc)/h_c_neu(j,i)
219 end do
220 
221 cw5 = aw5/(h_t_neu(j,i)**2)
222 cw8 = aw8
223 ci3 = ai3/(h_t_neu(j,i)**2)
224 
225 #if ADV_VERT==1
226 
227 do kt=1, ktmax-1
228  cw1(kt) = aw1/h_t_neu(j,i)*0.5_dp*(vz_t(kt,j,i)+vz_t(kt-1,j,i))
229 end do
230 
231 kt=ktmax
232 cw1(kt) = aw1/h_t_neu(j,i)*0.5_dp*(vz_t(kt-1,j,i)+vz_c(0,j,i))
233 
234 kt=0
235 cw1_sg(kt) = aw1/h_t_neu(j,i)*vz_t(kt,j,i)
236  ! only needed for kt=0 ...
237 kt=ktmax-1
238 cw1_sg(kt) = aw1/h_t_neu(j,i)*vz_t(kt,j,i)
239  ! ... and kt=KTMAX-1
240 
241 #elif ( ADV_VERT==2 || ADV_VERT==3 )
242 
243 do kt=0, ktmax-1
244  cw1_sg(kt) = aw1/h_t_neu(j,i)*vz_t(kt,j,i)
245 end do
246 
247 #endif
248 
249 do kt=1, ktmax-1
250  cw9(kt) = aw9 &
251  *c_val(temp_t_m(kt,j,i)) &
252  *( dzs_dtau(j,i) &
253  +0.5_dp*(vx_t(kt,j,i)+vx_t(kt,j,i-1))*dzs_dxi_g(j,i) &
254  +0.5_dp*(vy_t(kt,j,i)+vy_t(kt,j-1,i))*dzs_deta_g(j,i) &
255  -0.5_dp*(vz_t(kt,j,i)+vz_t(kt-1,j,i)) )
256 end do
257 
258 do kt=0, ktmax
259 
260  cw2(kt) = aw2*(dzb_dtau(j,i)+zeta_t(kt)*dh_t_dtau(j,i)) &
261  /h_t_neu(j,i)
262  cw3(kt) = aw3*(dzb_dxi_g(j,i)+zeta_t(kt)*dh_t_dxi_g(j,i)) &
263  /h_t_neu(j,i) &
264  *0.5_dp*(vx_t(kt,j,i)+vx_t(kt,j,i-1)) *insq_g11_g(j,i)
265  cw4(kt) = aw4*(dzb_deta_g(j,i)+zeta_t(kt)*dh_t_deta_g(j,i)) &
266  /h_t_neu(j,i) &
267  *0.5_dp*(vy_t(kt,j,i)+vy_t(kt,j-1,i)) *insq_g22_g(j,i)
268  sigma_t_help(kt) &
269  = sigma_c_help(0) &
270  + rho*g*h_t_neu(j,i)*(1.0_dp-zeta_t(kt)) &
271  *(dzs_dxi_g(j,i)**2+dzs_deta_g(j,i)**2)**0.5_dp
272  cw7(kt) = aw7 &
273  *enh_t(kt,j,i) &
274  *ratefac_t(omega_t(kt,j,i)) &
275  *creep(sigma_t_help(kt)) &
276  *sigma_t_help(kt)*sigma_t_help(kt)
277 
278 ! ------ Fluxes for horizontal advection of water content and age
279 
280 #if ADV_HOR==4
281  if (vx_t(kt,j,i-1) >= zero) then
282  fwx_t_l(kt) = omega_t(kt,j,i-1)*vx_t(kt,j,i-1)
283  fax_t_l(kt) = age_t(kt,j,i-1)*vx_t(kt,j,i-1)
284  else
285  fwx_t_l(kt) = omega_t(kt,j,i)*vx_t(kt,j,i-1)
286  fax_t_l(kt) = age_t(kt,j,i)*vx_t(kt,j,i-1)
287  end if
288 
289  if (vx_t(kt,j,i) >= zero) then
290  fwx_t_r(kt) = omega_t(kt,j,i)*vx_t(kt,j,i)
291  fax_t_r(kt) = age_t(kt,j,i)*vx_t(kt,j,i)
292  else
293  fwx_t_r(kt) = omega_t(kt,j,i+1)*vx_t(kt,j,i)
294  fax_t_r(kt) = age_t(kt,j,i+1)*vx_t(kt,j,i)
295  end if
296 
297  if (vy_t(kt,j-1,i) >= zero) then
298  fwy_t_l(kt) = omega_t(kt,j-1,i)*vy_t(kt,j-1,i)
299  fay_t_l(kt) = age_t(kt,j-1,i)*vy_t(kt,j-1,i)
300  else
301  fwy_t_l(kt) = omega_t(kt,j,i)*vy_t(kt,j-1,i)
302  fay_t_l(kt) = age_t(kt,j,i)*vy_t(kt,j-1,i)
303  end if
304 
305  if (vy_t(kt,j,i) >= zero) then
306  fwy_t_r(kt) = omega_t(kt,j,i)*vy_t(kt,j,i)
307  fay_t_r(kt) = age_t(kt,j,i)*vy_t(kt,j,i)
308  else
309  fwy_t_r(kt) = omega_t(kt,j+1,i)*vy_t(kt,j,i)
310  fay_t_r(kt) = age_t(kt,j+1,i)*vy_t(kt,j,i)
311  end if
312 #endif
313 
314 end do
315 
316 #if ADV_VERT==1
317 
318 kt=0
319 cw2_sg(kt) = 0.5_dp*(cw2(kt)+cw2(kt+1))
320 cw3_sg(kt) = 0.5_dp*(cw3(kt)+cw3(kt+1))
321 cw4_sg(kt) = 0.5_dp*(cw4(kt)+cw4(kt+1))
322 adv_vert_w_sg(kt) = cw1_sg(kt)-cw2_sg(kt)-cw3_sg(kt)-cw4_sg(kt)
323 abs_adv_vert_w_sg(kt) = abs(adv_vert_w_sg(kt)) ! only needed for kt=0 ...
324 kt=ktmax-1
325 cw2_sg(kt) = 0.5_dp*(cw2(kt)+cw2(kt+1))
326 cw3_sg(kt) = 0.5_dp*(cw3(kt)+cw3(kt+1))
327 cw4_sg(kt) = 0.5_dp*(cw4(kt)+cw4(kt+1))
328 adv_vert_w_sg(kt) = cw1_sg(kt)-cw2_sg(kt)-cw3_sg(kt)-cw4_sg(kt)
329 abs_adv_vert_w_sg(kt) = abs(adv_vert_w_sg(kt)) ! ... and kt=KTMAX-1
330 
331 #elif ( ADV_VERT==2 || ADV_VERT==3 )
332 
333 do kt=0, ktmax-1
334  cw2_sg(kt) = 0.5_dp*(cw2(kt)+cw2(kt+1))
335  cw3_sg(kt) = 0.5_dp*(cw3(kt)+cw3(kt+1))
336  cw4_sg(kt) = 0.5_dp*(cw4(kt)+cw4(kt+1))
337  adv_vert_w_sg(kt) = cw1_sg(kt)-cw2_sg(kt)-cw3_sg(kt)-cw4_sg(kt)
338  abs_adv_vert_w_sg(kt) = abs(adv_vert_w_sg(kt))
339 end do
340 
341 #endif
342 
343 #if ( ADV_HOR==3 || ADV_HOR==4 )
344 dtt_dxi = 2.0_dp*dtt_2dxi
345 dtt_deta = 2.0_dp*dtt_2deta
346 #endif
347 
348 !-------- Set up the equations for the bedrock temperature --------
349 
350 kr=0
351 lgs_a1(kr) = 1.0_dp
352 lgs_a2(kr) = -1.0_dp
353 lgs_b(kr) = clb1
354 
355 #if Q_LITHO==1
356 ! (coupled heat-conducting bedrock)
357 
358 do kr=1, krmax-1
359  lgs_a0(kr) = - ctr1
360  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
361  lgs_a2(kr) = - ctr1
362  lgs_b(kr) = temp_r(kr,j,i)
363 end do
364 
365 #elif Q_LITHO==0
366 ! (no coupled heat-conducting bedrock)
367 
368 do kr=1, krmax-1
369  lgs_a0(kr) = 1.0_dp
370  lgs_a1(kr) = 0.0_dp
371  lgs_a2(kr) = -1.0_dp
372  lgs_b(kr) = 2.0_dp*clb1
373 end do
374 
375 #endif
376 
377 kr=krmax
378 lgs_a0(kr) = 0.0_dp
379 lgs_a1(kr) = 1.0_dp
380 lgs_b(kr) = temp_t_m(0,j,i)
381 
382 !-------- Solve system of linear equations --------
383 
384 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
385 
386 !-------- Assign the result --------
387 
388 do kr=0, krmax
389  temp_r_neu(kr,j,i) = lgs_x(kr)
390 end do
391 
392 !-------- Set up the equations for the water content in
393 ! temperate ice --------
394 
395 kt=0
396 lgs_a1(kt) = 1.0_dp
397 lgs_a2(kt) = -1.0_dp
398 lgs_b(kt) = 0.0_dp
399 
400 do kt=1, ktmax-1
401 
402 #if ADV_VERT==1
403 
404  lgs_a0(kt) = -0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - cw5
405  lgs_a1(kt) = 1.0_dp + 2.0_dp*cw5
406  lgs_a2(kt) = 0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - cw5
407 
408 #elif ADV_VERT==2
409 
410  lgs_a0(kt) &
411  = -0.5_dp*(adv_vert_w_sg(kt-1)+abs_adv_vert_w_sg(kt-1)) &
412  -cw5
413  lgs_a1(kt) &
414  = 1.0_dp &
415  +0.5_dp*(adv_vert_w_sg(kt-1)+abs_adv_vert_w_sg(kt-1)) &
416  -0.5_dp*(adv_vert_w_sg(kt) -abs_adv_vert_w_sg(kt) ) &
417  +2.0_dp*cw5
418  lgs_a2(kt) &
419  = 0.5_dp*(adv_vert_w_sg(kt) -abs_adv_vert_w_sg(kt) ) &
420  -cw5
421 
422 #elif ADV_VERT==3
423 
424  adv_vert_w_help = 0.5_dp*(adv_vert_w_sg(kt)+adv_vert_w_sg(kt-1))
425 
426  lgs_a0(kt) &
427  = -max(adv_vert_w_help, 0.0_dp) &
428  -cw5
429  lgs_a1(kt) &
430  = 1.0_dp &
431  +max(adv_vert_w_help, 0.0_dp)-min(adv_vert_w_help, 0.0_dp) &
432  +2.0_dp*cw5
433  lgs_a2(kt) &
434  = min(adv_vert_w_help, 0.0_dp) &
435  -cw5
436 
437 #endif
438 
439 #if ADV_HOR==2
440 
441  lgs_b(kt) = omega_t(kt,j,i) + cw7(kt) + cw8 + cw9(kt) &
442  -dtt_2dxi* &
443  ( (vx_t(kt,j,i)-abs(vx_t(kt,j,i))) &
444  *(omega_t(kt,j,i+1)-omega_t(kt,j,i)) &
445  *insq_g11_sgx(j,i) &
446  +(vx_t(kt,j,i-1)+abs(vx_t(kt,j,i-1))) &
447  *(omega_t(kt,j,i)-omega_t(kt,j,i-1)) &
448  *insq_g11_sgx(j,i-1) ) &
449  -dtt_2deta* &
450  ( (vy_t(kt,j,i)-abs(vy_t(kt,j,i))) &
451  *(omega_t(kt,j+1,i)-omega_t(kt,j,i)) &
452  *insq_g22_sgy(j,i) &
453  +(vy_t(kt,j-1,i)+abs(vy_t(kt,j-1,i))) &
454  *(omega_t(kt,j,i)-omega_t(kt,j-1,i)) &
455  *insq_g22_sgy(j-1,i) )
456 
457 #elif ADV_HOR==3
458 
459  vx_t_help = 0.5_dp*(vx_t(kt,j,i)+vx_t(kt,j,i-1))
460  vy_t_help = 0.5_dp*(vy_t(kt,j,i)+vy_t(kt,j-1,i))
461 
462  lgs_b(kt) = omega_t(kt,j,i) + cw7(kt) + cw8 + cw9(kt) &
463  -dtt_dxi* &
464  ( min(vx_t_help, 0.0_dp) &
465  *(omega_t(kt,j,i+1)-omega_t(kt,j,i)) &
466  *insq_g11_sgx(j,i) &
467  +max(vx_t_help, 0.0_dp) &
468  *(omega_t(kt,j,i)-omega_t(kt,j,i-1)) &
469  *insq_g11_sgx(j,i-1) ) &
470  -dtt_deta* &
471  ( min(vy_t_help, 0.0_dp) &
472  *(omega_t(kt,j+1,i)-omega_t(kt,j,i)) &
473  *insq_g22_sgy(j,i) &
474  +max(vy_t_help, 0.0_dp) &
475  *(omega_t(kt,j,i)-omega_t(kt,j-1,i)) &
476  *insq_g22_sgy(j-1,i) )
477 
478 #elif ADV_HOR==4
479 
480  lgs_b(kt) = omega_t(kt,j,i) + cw7(kt) + cw8 + cw9(kt) &
481  -dtt_dxi *(fwx_t_r(kt)-fwx_t_l(kt)) &
482  -dtt_deta*(fwy_t_r(kt)-fwy_t_l(kt))
483 
484 #endif
485 
486 end do
487 
488 kt=ktmax
489 
490 if (am_perp(j,i) >= zero) then ! melting condition
491  lgs_a0(kt) = 0.0_dp
492  lgs_a1(kt) = 1.0_dp
493  lgs_b(kt) = 0.0_dp
494 else ! am_perp(j,i) < 0.0, freezing condition
495  lgs_a0(kt) = -1.0_dp
496  lgs_a1(kt) = 1.0_dp
497  lgs_b(kt) = 0.0_dp
498 end if
499 
500 !-------- Solve system of linear equations --------
501 
502 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, ktmax)
503 
504 !-------- Assign the result, compute the water drainage --------
505 
506 q_tld(j,i) = 0.0_dp
507 
508 do kt=0, ktmax
509 
510  if (lgs_x(kt) < zero) then
511  omega_t_neu(kt,j,i) = 0.0_dp ! (as a precaution)
512  else if (lgs_x(kt) < omega_max) then
513  omega_t_neu(kt,j,i) = lgs_x(kt)
514  else
515  omega_t_neu(kt,j,i) = omega_max
516  q_tld(j,i) = q_tld(j,i) &
517  +aqtld*h_t_neu(j,i)*(lgs_x(kt)-omega_max)
518  end if
519 
520 end do
521 
522 !-------- Set up the equations for the ice temperature --------
523 
524 ! ------ Abbreviation for the jump of the temperature gradient with
525 ! the new omega
526 
527 if (am_perp(j,i) >= zero) then ! melting condition
528  cm2 = 0.0_dp
529 else ! am_perp(j,i) < 0.0, freezing condition
530  cm2 = am2*h_c_neu(j,i)*omega_t_neu(ktmax,j,i)*am_perp(j,i) &
531  /kappa_val(temp_c(0,j,i))
532 end if
533 
534 kc=0
535 lgs_a1(kc) = 1.0_dp
536 lgs_a2(kc) = -1.0_dp
537 lgs_b(kc) = -cm1-cm2
538 
539 do kc=1, kcmax-1
540 
541 #if ADV_VERT==1
542 
543  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
544  -ct5(kc)*ct6(kc-1)
545  lgs_a1(kc) = 1.0_dp+ct5(kc)*(ct6(kc)+ct6(kc-1))
546  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
547  -ct5(kc)*ct6(kc)
548 
549 #elif ADV_VERT==2
550 
551  lgs_a0(kc) &
552  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
553  -ct5(kc)*ct6(kc-1)
554  lgs_a1(kc) &
555  = 1.0_dp &
556  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
557  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
558  +ct5(kc)*(ct6(kc)+ct6(kc-1))
559  lgs_a2(kc) &
560  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
561  -ct5(kc)*ct6(kc)
562 
563 #elif ADV_VERT==3
564 
565  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
566 
567  lgs_a0(kc) &
568  = -max(adv_vert_help, 0.0_dp) &
569  -ct5(kc)*ct6(kc-1)
570  lgs_a1(kc) &
571  = 1.0_dp &
572  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
573  +ct5(kc)*(ct6(kc)+ct6(kc-1))
574  lgs_a2(kc) &
575  = min(adv_vert_help, 0.0_dp) &
576  -ct5(kc)*ct6(kc)
577 
578 #endif
579 
580 #if ADV_HOR==2
581 
582  lgs_b(kc) = temp_c(kc,j,i) + ct7(kc) &
583  -dtt_2dxi* &
584  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
585  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
586  *insq_g11_sgx(j,i) &
587  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
588  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
589  *insq_g11_sgx(j,i-1) ) &
590  -dtt_2deta* &
591  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
592  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
593  *insq_g22_sgy(j,i) &
594  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
595  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
596  *insq_g22_sgy(j-1,i) )
597 
598 #elif ADV_HOR==3
599 
600  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
601  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
602 
603  lgs_b(kc) = temp_c(kc,j,i) + ct7(kc) &
604  -dtt_dxi* &
605  ( min(vx_c_help, 0.0_dp) &
606  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
607  *insq_g11_sgx(j,i) &
608  +max(vx_c_help, 0.0_dp) &
609  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
610  *insq_g11_sgx(j,i-1) ) &
611  -dtt_deta* &
612  ( min(vy_c_help, 0.0_dp) &
613  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
614  *insq_g22_sgy(j,i) &
615  +max(vy_c_help, 0.0_dp) &
616  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
617  *insq_g22_sgy(j-1,i) )
618 
619 #elif ADV_HOR==4
620 
621  lgs_b(kc) = temp_c(kc,j,i) + ct7(kc) &
622  -dtt_dxi *(ftx_c_r(kc)-ftx_c_l(kc)) &
623  -dtt_deta*(fty_c_r(kc)-fty_c_l(kc))
624 
625 #endif
626 
627 end do
628 
629 kc=kcmax
630 lgs_a0(kc) = 0.0_dp
631 lgs_a1(kc) = 1.0_dp
632 lgs_b(kc) = temp_s(j,i)
633 
634 !-------- Solve system of linear equations --------
635 
636 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
637 
638 !-------- Assign the result --------
639 
640 do kc=0, kcmax
641  temp_c_neu(kc,j,i) = lgs_x(kc)
642 end do
643 
644 !-------- Set up the equations for the age (cold and temperate ice
645 ! simultaneously) --------
646 
647 kt=0 ! adv_vert_w_sg(0) <= 0
648 lgs_a1(kt) = 1.0_dp - min(adv_vert_w_sg(kt), 0.0_dp) ! (directed downward)
649 lgs_a2(kt) = min(adv_vert_w_sg(kt), 0.0_dp) ! assumed/enforced
650 
651 #if ADV_HOR==2
652 
653 lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
654  -dtt_2dxi* &
655  ( (vx_t(kt,j,i)-abs(vx_t(kt,j,i))) &
656  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
657  *insq_g11_sgx(j,i) &
658  +(vx_t(kt,j,i-1)+abs(vx_t(kt,j,i-1))) &
659  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
660  *insq_g11_sgx(j,i-1) ) &
661  -dtt_2deta* &
662  ( (vy_t(kt,j,i)-abs(vy_t(kt,j,i))) &
663  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
664  *insq_g22_sgy(j,i) &
665  +(vy_t(kt,j-1,i)+abs(vy_t(kt,j-1,i))) &
666  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
667  *insq_g22_sgy(j-1,i) )
668 
669 #elif ADV_HOR==3
670 
671 vx_t_help = 0.5_dp*(vx_t(kt,j,i)+vx_t(kt,j,i-1))
672 vy_t_help = 0.5_dp*(vy_t(kt,j,i)+vy_t(kt,j-1,i))
673 
674 lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
675  -dtt_dxi* &
676  ( min(vx_t_help, 0.0_dp) &
677  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
678  *insq_g11_sgx(j,i) &
679  +max(vx_t_help, 0.0_dp) &
680  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
681  *insq_g11_sgx(j,i-1) ) &
682  -dtt_deta* &
683  ( min(vy_t_help, 0.0_dp) &
684  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
685  *insq_g22_sgy(j,i) &
686  +max(vy_t_help, 0.0_dp) &
687  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
688  *insq_g22_sgy(j-1,i) )
689 
690 #elif ADV_HOR==4
691 
692 lgs_b(kt) = ...
693 
694 #endif
695 
696 do kt=1, ktmax-1
697 
698 #if ADV_VERT==1
699 
700  lgs_a0(kt) = -0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - ci3
701  lgs_a1(kt) = 1.0_dp + 2.0_dp*ci3
702  lgs_a2(kt) = 0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - ci3
703 
704 #elif ADV_VERT==2
705 
706  lgs_a0(kt) = -0.5_dp*(adv_vert_w_sg(kt-1)+abs_adv_vert_w_sg(kt-1))
707  lgs_a1(kt) = 1.0_dp &
708  +0.5_dp*(adv_vert_w_sg(kt-1)+abs_adv_vert_w_sg(kt-1)) &
709  -0.5_dp*(adv_vert_w_sg(kt) -abs_adv_vert_w_sg(kt) )
710  lgs_a2(kt) = 0.5_dp*(adv_vert_w_sg(kt) -abs_adv_vert_w_sg(kt) )
711 
712 #elif ADV_VERT==3
713 
714  adv_vert_w_help = 0.5_dp*(adv_vert_w_sg(kt)+adv_vert_w_sg(kt-1))
715 
716  lgs_a0(kt) = -max(adv_vert_w_help, 0.0_dp)
717  lgs_a1(kt) = 1.0_dp &
718  +max(adv_vert_w_help, 0.0_dp)-min(adv_vert_w_help, 0.0_dp)
719  lgs_a2(kt) = min(adv_vert_w_help, 0.0_dp)
720 
721 #endif
722 
723 #if ADV_HOR==2
724 
725  lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
726  -dtt_2dxi* &
727  ( (vx_t(kt,j,i)-abs(vx_t(kt,j,i))) &
728  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
729  *insq_g11_sgx(j,i) &
730  +(vx_t(kt,j,i-1)+abs(vx_t(kt,j,i-1))) &
731  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
732  *insq_g11_sgx(j,i-1) ) &
733  -dtt_2deta* &
734  ( (vy_t(kt,j,i)-abs(vy_t(kt,j,i))) &
735  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
736  *insq_g22_sgy(j,i) &
737  +(vy_t(kt,j-1,i)+abs(vy_t(kt,j-1,i))) &
738  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
739  *insq_g22_sgy(j-1,i) )
740 
741 #elif ADV_HOR==3
742 
743  vx_t_help = 0.5_dp*(vx_t(kt,j,i)+vx_t(kt,j,i-1))
744  vy_t_help = 0.5_dp*(vy_t(kt,j,i)+vy_t(kt,j-1,i))
745 
746  lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
747  -dtt_dxi* &
748  ( min(vx_t_help, 0.0_dp) &
749  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
750  *insq_g11_sgx(j,i) &
751  +max(vx_t_help, 0.0_dp) &
752  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
753  *insq_g11_sgx(j,i-1) ) &
754  -dtt_deta* &
755  ( min(vy_t_help, 0.0_dp) &
756  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
757  *insq_g22_sgy(j,i) &
758  +max(vy_t_help, 0.0_dp) &
759  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
760  *insq_g22_sgy(j-1,i) )
761 
762 #elif ADV_HOR==4
763 
764  lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
765  -dtt_dxi *(fax_t_r(kt)-fax_t_l(kt)) &
766  -dtt_deta*(fay_t_r(kt)-fay_t_l(kt))
767 
768 #endif
769 
770 end do
771 
772 #if ADV_VERT==1
773 
774 kt=ktmax
775 kc=0
776 
777 lgs_a0(kt) = -0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - ci3
778 lgs_a1(kt) = 1.0_dp + 2.0_dp*ci3
779 lgs_a2(kt) = 0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - ci3
780 
781 #if ADV_HOR==2
782 
783 lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
784  -dtt_2dxi* &
785  ( (vx_t(kt,j,i)-abs(vx_t(kt,j,i))) &
786  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
787  *insq_g11_sgx(j,i) &
788  +(vx_t(kt,j,i-1)+abs(vx_t(kt,j,i-1))) &
789  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
790  *insq_g11_sgx(j,i-1) ) &
791  -dtt_2deta* &
792  ( (vy_t(kt,j,i)-abs(vy_t(kt,j,i))) &
793  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
794  *insq_g22_sgy(j,i) &
795  +(vy_t(kt,j-1,i)+abs(vy_t(kt,j-1,i))) &
796  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
797  *insq_g22_sgy(j-1,i) )
798 
799 #elif ADV_HOR==3
800 
801 vx_t_help = 0.5_dp*(vx_t(kt,j,i)+vx_t(kt,j,i-1))
802 vy_t_help = 0.5_dp*(vy_t(kt,j,i)+vy_t(kt,j-1,i))
803 
804 lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
805  -dtt_dxi* &
806  ( min(vx_t_help, 0.0_dp) &
807  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
808  *insq_g11_sgx(j,i) &
809  +max(vx_t_help, 0.0_dp) &
810  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
811  *insq_g11_sgx(j,i-1) ) &
812  -dtt_deta* &
813  ( min(vy_t_help, 0.0_dp) &
814  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
815  *insq_g22_sgy(j,i) &
816  +max(vy_t_help, 0.0_dp) &
817  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
818  *insq_g22_sgy(j-1,i) )
819 
820 #elif ADV_HOR==4
821 
822 lgs_b(kt) = ...
823 
824 #endif
825 
826 #elif ( ADV_VERT==2 || ADV_VERT==3 )
827 
828 kt=ktmax
829 kc=0
830 
831 if (adv_vert_sg(kc) <= zero) then
832 
833  lgs_a0(ktmax+kc) = 0.0_dp
834  lgs_a1(ktmax+kc) = 1.0_dp - adv_vert_sg(kc)
835  lgs_a2(ktmax+kc) = adv_vert_sg(kc)
836 
837 #if ADV_HOR==2
838 
839  lgs_b(ktmax+kc) = age_c(kc,j,i) + dtime_temp &
840  -dtt_2dxi* &
841  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
842  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
843  *insq_g11_sgx(j,i) &
844  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
845  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
846  *insq_g11_sgx(j,i-1) ) &
847  -dtt_2deta* &
848  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
849  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
850  *insq_g22_sgy(j,i) &
851  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
852  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
853  *insq_g22_sgy(j-1,i) )
854 
855 #elif ADV_HOR==3
856 
857  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
858  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
859 
860  lgs_b(ktmax+kc) = age_c(kc,j,i) + dtime_temp &
861  -dtt_dxi* &
862  ( min(vx_c_help, 0.0_dp) &
863  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
864  *insq_g11_sgx(j,i) &
865  +max(vx_c_help, 0.0_dp) &
866  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
867  *insq_g11_sgx(j,i-1) ) &
868  -dtt_deta* &
869  ( min(vy_c_help, 0.0_dp) &
870  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
871  *insq_g22_sgy(j,i) &
872  +max(vy_c_help, 0.0_dp) &
873  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
874  *insq_g22_sgy(j-1,i) )
875 
876 #elif ADV_HOR==4
877 
878  lgs_b(ktmax+kc) = ...
879 
880 #endif
881 
882 else if (adv_vert_w_sg(kt-1) >= zero) then
883 
884  lgs_a0(kt) = -adv_vert_w_sg(kt-1)
885  lgs_a1(kt) = 1.0_dp + adv_vert_w_sg(kt-1)
886  lgs_a2(kt) = 0.0_dp
887 
888 #if ADV_HOR==2
889 
890  lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
891  -dtt_2dxi* &
892  ( (vx_t(kt,j,i)-abs(vx_t(kt,j,i))) &
893  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
894  *insq_g11_sgx(j,i) &
895  +(vx_t(kt,j,i-1)+abs(vx_t(kt,j,i-1))) &
896  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
897  *insq_g11_sgx(j,i-1) ) &
898  -dtt_2deta* &
899  ( (vy_t(kt,j,i)-abs(vy_t(kt,j,i))) &
900  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
901  *insq_g22_sgy(j,i) &
902  +(vy_t(kt,j-1,i)+abs(vy_t(kt,j-1,i))) &
903  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
904  *insq_g22_sgy(j-1,i) )
905 
906 #elif ADV_HOR==3
907 
908  vx_t_help = 0.5_dp*(vx_t(kt,j,i)+vx_t(kt,j,i-1))
909  vy_t_help = 0.5_dp*(vy_t(kt,j,i)+vy_t(kt,j-1,i))
910 
911  lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
912  -dtt_dxi* &
913  ( min(vx_t_help, 0.0_dp) &
914  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
915  *insq_g11_sgx(j,i) &
916  +max(vx_t_help, 0.0_dp) &
917  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
918  *insq_g11_sgx(j,i-1) ) &
919  -dtt_deta* &
920  ( min(vy_t_help, 0.0_dp) &
921  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
922  *insq_g22_sgy(j,i) &
923  +max(vy_t_help, 0.0_dp) &
924  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
925  *insq_g22_sgy(j-1,i) )
926 
927 #elif ADV_HOR==4
928 
929  lgs_b(kt) = ...
930 
931 #endif
932 
933 else
934 
935  lgs_a0(kt) = -0.5_dp
936  lgs_a1(kt) = 1.0_dp
937  lgs_a2(kt) = -0.5_dp
938  lgs_b(kt) = 0.0_dp
939  ! Makeshift: Average of age_c(kc=1) and age_t(kt=KTMAX-1)
940 
941 end if
942 
943 #endif
944 
945 do kc=1, kcmax-1
946 
947 #if ADV_VERT==1
948 
949  lgs_a0(ktmax+kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
950  -ci1(kc)*ci2(kc-1)
951  lgs_a1(ktmax+kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
952  lgs_a2(ktmax+kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
953  -ci1(kc)*ci2(kc)
954 
955 #elif ADV_VERT==2
956 
957  lgs_a0(ktmax+kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
958  lgs_a1(ktmax+kc) = 1.0_dp &
959  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
960  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
961  lgs_a2(ktmax+kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
962 
963 #elif ADV_VERT==3
964 
965  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
966 
967  lgs_a0(ktmax+kc) = -max(adv_vert_help, 0.0_dp)
968  lgs_a1(ktmax+kc) = 1.0_dp &
969  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp)
970  lgs_a2(ktmax+kc) = min(adv_vert_help, 0.0_dp)
971 
972 #endif
973 
974 #if ADV_HOR==2
975 
976  lgs_b(ktmax+kc) = age_c(kc,j,i) + dtime_temp &
977  -dtt_2dxi* &
978  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
979  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
980  *insq_g11_sgx(j,i) &
981  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
982  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
983  *insq_g11_sgx(j,i-1) ) &
984  -dtt_2deta* &
985  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
986  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
987  *insq_g22_sgy(j,i) &
988  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
989  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
990  *insq_g22_sgy(j-1,i) )
991 
992 #elif ADV_HOR==3
993 
994  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
995  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
996 
997  lgs_b(ktmax+kc) = age_c(kc,j,i) + dtime_temp &
998  -dtt_dxi* &
999  ( min(vx_c_help, 0.0_dp) &
1000  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1001  *insq_g11_sgx(j,i) &
1002  +max(vx_c_help, 0.0_dp) &
1003  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1004  *insq_g11_sgx(j,i-1) ) &
1005  -dtt_deta* &
1006  ( min(vy_c_help, 0.0_dp) &
1007  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1008  *insq_g22_sgy(j,i) &
1009  +max(vy_c_help, 0.0_dp) &
1010  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1011  *insq_g22_sgy(j-1,i) )
1012 
1013 #elif ADV_HOR==4
1014 
1015  lgs_b(ktmax+kc) = age_c(kc,j,i) + dtime_temp &
1016  -dtt_dxi *(fax_c_r(kc)-fax_c_l(kc)) &
1017  -dtt_deta*(fay_c_r(kc)-fay_c_l(kc))
1018 
1019 #endif
1020 
1021 end do
1022 
1023 kc=kcmax
1024 if (as_perp(j,i) >= zero) then
1025  lgs_a0(ktmax+kc) = 0.0_dp
1026  lgs_a1(ktmax+kc) = 1.0_dp
1027  lgs_b(ktmax+kc) = 0.0_dp
1028 else
1029  lgs_a0(ktmax+kc) = -max(adv_vert_sg(kc-1), 0.0_dp)
1030  lgs_a1(ktmax+kc) = 1.0_dp + max(adv_vert_sg(kc-1), 0.0_dp)
1031  ! adv_vert_sg(KCMAX-1) >= 0 (directed upward)
1032  ! assumed/enforced
1033 #if ADV_HOR==2
1034 
1035  lgs_b(ktmax+kc) = age_c(kc,j,i) + dtime_temp &
1036  -dtt_2dxi* &
1037  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
1038  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1039  *insq_g11_sgx(j,i) &
1040  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
1041  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1042  *insq_g11_sgx(j,i-1) ) &
1043  -dtt_2deta* &
1044  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
1045  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1046  *insq_g22_sgy(j,i) &
1047  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
1048  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1049  *insq_g22_sgy(j-1,i) )
1050 
1051 #elif ADV_HOR==3
1052 
1053  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
1054  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
1055 
1056  lgs_b(ktmax+kc) = age_c(kc,j,i) + dtime_temp &
1057  -dtt_dxi* &
1058  ( min(vx_c_help, 0.0_dp) &
1059  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1060  *insq_g11_sgx(j,i) &
1061  +max(vx_c_help, 0.0_dp) &
1062  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1063  *insq_g11_sgx(j,i-1) ) &
1064  -dtt_deta* &
1065  ( min(vy_c_help, 0.0_dp) &
1066  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1067  *insq_g22_sgy(j,i) &
1068  +max(vy_c_help, 0.0_dp) &
1069  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1070  *insq_g22_sgy(j-1,i) )
1071 
1072 #elif ADV_HOR==4
1073 
1074  lgs_b(ktmax+kc) = ...
1075 
1076 #endif
1077 
1078 end if
1079 
1080 !-------- Solve system of linear equations --------
1081 
1082 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax+ktmax)
1083 
1084 !-------- Assign the result,
1085 ! restriction to interval [0, AGE_MAX yr] --------
1086 
1087 do kt=0, ktmax
1088 
1089  age_t_neu(kt,j,i) = lgs_x(kt)
1090 
1091  if (age_t_neu(kt,j,i) < (age_min*year_sec)) &
1092  age_t_neu(kt,j,i) = 0.0_dp
1093  if (age_t_neu(kt,j,i) > (age_max*year_sec)) &
1094  age_t_neu(kt,j,i) = age_max*year_sec
1095 
1096 end do
1097 
1098 do kc=0, kcmax
1099 
1100  age_c_neu(kc,j,i) = lgs_x(ktmax+kc)
1101 
1102  if (age_c_neu(kc,j,i) < (age_min*year_sec)) &
1103  age_c_neu(kc,j,i) = 0.0_dp
1104  if (age_c_neu(kc,j,i) > (age_max*year_sec)) &
1105  age_c_neu(kc,j,i) = age_max*year_sec
1106 
1107 end do
1108 
1109 end subroutine calc_temp3
1110 !
Generated on Thu Sep 12 2013 09:45:10 for SICOPOLIS V3.1 by   doxygen 1.8.2