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