SICOPOLIS V3.3
calc_temp_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : c a l c _ t e m p _ m
4 !
5 !> @file
6 !!
7 !! Computation of temperature, water content and age.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2017 Ralf Greve
12 !!
13 !! @section License
14 !!
15 !! This file is part of SICOPOLIS.
16 !!
17 !! SICOPOLIS is free software: you can redistribute it and/or modify
18 !! it under the terms of the GNU General Public License as published by
19 !! the Free Software Foundation, either version 3 of the License, or
20 !! (at your option) any later version.
21 !!
22 !! SICOPOLIS is distributed in the hope that it will be useful,
23 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
24 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 !! GNU General Public License for more details.
26 !!
27 !! You should have received a copy of the GNU General Public License
28 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
29 !<
30 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
31 
32 !-------------------------------------------------------------------------------
33 !> Computation of temperature, water content and age.
34 !<------------------------------------------------------------------------------
35 module calc_temp_m
36 
37  use sico_types_m
38  use sico_variables_m
39  use sico_vars_m
40 
41  implicit none
42 
43  private
44  public :: calc_temp_poly, calc_temp_cold, calc_temp_const
45 
46 contains
47 
48 !-------------------------------------------------------------------------------
49 !> Computation of temperature, water content and age in polythermal mode.
50 !<------------------------------------------------------------------------------
51 subroutine calc_temp_poly(dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
52  dtime_temp)
53 
54 use ice_material_properties_m, only : kappa_val
55 
56 implicit none
57 
58 real(dp), intent(in) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
59 real(dp), intent(in) :: dtime_temp
60 
61 integer(i4b) :: i, j, kc, kt, kr, ii, jj
62 real(dp) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
63  at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
64  at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
65  acb1, acb2, acb3, acb4, &
66  ai1(0:kcmax), ai2(0:kcmax), ai3, &
67  atr1, am1, am2, alb1
68 real(dp) :: aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld
69 real(dp) :: dtime_temp_inv, dtt_2dxi, dtt_2deta
70 real(dp) :: temp_c_help(0:kcmax)
71 real(dp) :: fact_thick
72 real(dp) :: time_lag_cts
73 real(dp) :: Vol_t, Vol_t_smooth, korrfakt_t, &
74  dH_t_smooth(0:jmax,0:imax)
75 
76 !-------- Term abbreviations
77 
78 at7 = 2.0_dp/rho*dtime_temp
79 
80 aw1 = dtime_temp/dzeta_t
81 aw2 = dtime_temp/dzeta_t
82 aw3 = dtime_temp/dzeta_t
83 aw4 = dtime_temp/dzeta_t
84 aw5 = nue/rho*dtime_temp/(dzeta_t**2)
85 aw7 = 2.0_dp/(rho*l)*dtime_temp
86 aw8 = beta**2/(rho*l) &
87  *(kappa_val(0.0_dp)-kappa_val(-1.0_dp))*dtime_temp
88 aw9 = beta/l*dtime_temp
89 
90 ai3 = agediff*dtime_temp/(dzeta_t**2)
91 
92 atr1 = kappa_r/(rho_c_r*h_r**2)*dtime_temp/(dzeta_r**2)
93 
94 if (flag_aa_nonzero) then
95  am1 = aa*beta*dzeta_c/(ea-1.0_dp)
96  am2 = aa*l*rho*dzeta_c/(ea-1.0_dp)
97 else
98  am1 = beta*dzeta_c
99  am2 = l*rho*dzeta_c
100 end if
101 
102 if (flag_aa_nonzero) then
103  acb1 = (ea-1.0_dp)/aa/dzeta_c
104 else
105  acb1 = 1.0_dp/dzeta_c
106 end if
107 
108 acb2 = kappa_r/h_r/dzeta_r
109 acb3 = rho*g
110 acb4 = rho*g
111 
112 alb1 = h_r/kappa_r*dzeta_r
113 
114 aqtld = dzeta_t/dtime_temp
115 
116 dtt_2dxi = 0.5_dp*dtime_temp/dxi
117 dtt_2deta = 0.5_dp*dtime_temp/deta
118 
119 dtime_temp_inv = 1.0_dp/dtime_temp
120 
121 do kc=0, kcmax
122 
123  if (flag_aa_nonzero) then
124 
125  at1(kc) = (ea-1.0_dp)/(aa*eaz_c(kc))*dtime_temp/dzeta_c
126  at2_1(kc) = (ea-1.0_dp)/(aa*eaz_c(kc))*dtime_temp/dzeta_c
127  at2_2(kc) = (eaz_c(kc)-1.0_dp)/(aa*eaz_c(kc)) &
128  *dtime_temp/dzeta_c
129  at3_1(kc) = (ea-1.0_dp)/(aa*eaz_c(kc))*dtime_temp/dzeta_c
130  at3_2(kc) = (eaz_c(kc)-1.0_dp)/(aa*eaz_c(kc)) &
131  *dtime_temp/dzeta_c
132  at4_1(kc) = (ea-1.0_dp)/(aa*eaz_c(kc))*dtime_temp/dzeta_c
133  at4_2(kc) = (eaz_c(kc)-1.0_dp)/(aa*eaz_c(kc)) &
134  *dtime_temp/dzeta_c
135  at5(kc) = (ea-1.0_dp)/(rho*aa*eaz_c(kc)) &
136  *dtime_temp/dzeta_c
137  if (kc /= kcmax) then
138  at6(kc) = (ea-1.0_dp) &
139  /(aa*exp(aa*0.5_dp*(zeta_c(kc)+zeta_c(kc+1)))) &
140  /dzeta_c
141  else
142  at6(kc) = 0.0_dp
143  end if
144  ai1(kc) = agediff*(ea-1.0_dp)/(aa*eaz_c(kc)) &
145  *dtime_temp/dzeta_c
146  if (kc /= kcmax) then
147  ai2(kc) = (ea-1.0_dp) &
148  /(aa*exp(aa*0.5_dp*(zeta_c(kc)+zeta_c(kc+1)))) &
149  /dzeta_c
150  else
151  ai2(kc) = 0.0_dp
152  end if
153 
154  else
155 
156  at1(kc) = dtime_temp/dzeta_c
157  at2_1(kc) = dtime_temp/dzeta_c
158  at2_2(kc) = zeta_c(kc) &
159  *dtime_temp/dzeta_c
160  at3_1(kc) = dtime_temp/dzeta_c
161  at3_2(kc) = zeta_c(kc) &
162  *dtime_temp/dzeta_c
163  at4_1(kc) = dtime_temp/dzeta_c
164  at4_2(kc) = zeta_c(kc) &
165  *dtime_temp/dzeta_c
166  at5(kc) = 1.0_dp/rho &
167  *dtime_temp/dzeta_c
168  if (kc /= kcmax) then
169  at6(kc) = 1.0_dp &
170  /dzeta_c
171  else
172  at6(kc) = 0.0_dp
173  end if
174  ai1(kc) = agediff &
175  *dtime_temp/dzeta_c
176  if (kc /= kcmax) then
177  ai2(kc) = 1.0_dp &
178  /dzeta_c
179  else
180  ai2(kc) = 0.0_dp
181  end if
182 
183  end if
184 
185 end do
186 
187 !-------- Computation loop for temperature, water content and age --------
188 
189 do i=1, imax-1 ! skipping domain margins
190 do j=1, jmax-1 ! skipping domain margins
191 
192  if (maske(j,i)==0) then ! glaciated land
193 
194 ! ------ Old vertical column cold
195 
196  if (n_cts(j,i).eq.-1) then
197 
198  n_cts_neu(j,i) = n_cts(j,i)
199  zm_neu(j,i) = zb(j,i)
200  h_c_neu(j,i) = h_c(j,i)
201  h_t_neu(j,i) = h_t(j,i)
202 
203  call calc_temp1(at1, at2_1, at2_2, at3_1, at3_2, &
204  at4_1, at4_2, at5, at6, at7, atr1, acb1, acb2, &
205  acb3, acb4, alb1, ai1, ai2, &
206  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
207 
208 ! ---- Check whether base has become temperate
209 
210  if (temp_c_neu(0,j,i).gt.temp_c_m(0,j,i)) then
211 
212  n_cts_neu(j,i) = 0
213 
214  call calc_temp2(at1, at2_1, at2_2, at3_1, at3_2, &
215  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
216  ai1, ai2, &
217  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
218 
219  end if
220 
221 ! ---- Check whether even temperate layer has formed
222 
223  if ( &
224  ( n_cts_neu(j,i).eq.0 ).and. &
225  ( (temp_c_neu(1,j,i)-temp_c_neu(0,j,i)) &
226  .gt.(am1*h_c_neu(j,i)) ) &
227  ) then
228 
229  n_cts_neu(j,i) = 1
230  zm_neu(j,i) = zb(j,i)+0.001_dp
231  h_c_neu(j,i) = h_c(j,i)-0.001_dp
232  h_t_neu(j,i) = h_t(j,i)+0.001_dp
233 ! ! CTS preliminarily positioned 1 mm above ice base --------
234 
235  call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
236  at4_1, at4_2, at5, at6, at7, atr1, &
237  am1, am2, alb1, &
238  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
239  ai1, ai2, ai3, dzeta_t, &
240  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
241 
242  call shift_cts_upward(at1, at2_1, at2_2, at3_1, at3_2, &
243  at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
244  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
245  ai1, ai2, ai3, dzeta_t, &
246  dtime_temp, dtt_2dxi, dtt_2deta, dtime_temp_inv, &
247  i, j)
248 
249  end if
250 
251 ! ------ Old vertical column with temperate base
252 
253  else if (n_cts(j,i).eq.0) then
254 
255  n_cts_neu(j,i) = n_cts(j,i)
256  zm_neu(j,i) = zb(j,i)
257  h_c_neu(j,i) = h_c(j,i)
258  h_t_neu(j,i) = h_t(j,i)
259 
260  call calc_temp2(at1, at2_1, at2_2, at3_1, at3_2, &
261  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
262  ai1, ai2, &
263  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
264 
265 ! ---- Check whether temperate base becomes cold
266 
267  if ( (temp_c_neu(1,j,i)-temp_c_neu(0,j,i)) &
268  .lt. (am1*h_c(j,i)) ) then
269 
270  n_cts_neu(j,i) = -1
271 
272  call calc_temp1(at1, at2_1, at2_2, at3_1, at3_2, &
273  at4_1, at4_2, at5, at6, at7, atr1, acb1, acb2, &
274  acb3, acb4, alb1, ai1, ai2, &
275  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
276 
277  if (temp_c_neu(0,j,i).ge.temp_c_m(0,j,i)) then
278 
279  n_cts_neu(j,i) = 0
280 
281  call calc_temp2(at1, at2_1, at2_2, at3_1, at3_2, &
282  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
283  ai1, ai2, &
284  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
285 
286  end if
287 
288  end if
289 
290 ! ---- Check whether temperate layer has formed
291 
292  if ( &
293  ( n_cts_neu(j,i).eq.0 ).and. &
294  ( (temp_c_neu(1,j,i)-temp_c_neu(0,j,i)) &
295  .gt.(am1*h_c_neu(j,i)) ) &
296  ) then
297 
298  n_cts_neu(j,i) = 1
299  zm_neu(j,i) = zb(j,i)+0.001_dp
300  h_c_neu(j,i) = h_c(j,i)-0.001_dp
301  h_t_neu(j,i) = h_t(j,i)+0.001_dp
302 ! ! CTS preliminarily positioned 1 mm above ice base --------
303 
304  call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
305  at4_1, at4_2, at5, at6, at7, atr1, &
306  am1, am2, alb1, &
307  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
308  ai1, ai2, ai3, dzeta_t, &
309  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
310 
311  call shift_cts_upward(at1, at2_1, at2_2, at3_1, at3_2, &
312  at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
313  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
314  ai1, ai2, ai3, dzeta_t, &
315  dtime_temp, dtt_2dxi, dtt_2deta, dtime_temp_inv, &
316  i, j)
317 
318  end if
319 
320 ! ------ Old vertical column with temperate base and temperate layer
321 
322  else ! n_cts(j,i).eq.1
323 
324  n_cts_neu(j,i) = n_cts(j,i)
325  zm_neu(j,i) = zm(j,i)
326  h_c_neu(j,i) = h_c(j,i)
327  h_t_neu(j,i) = h_t(j,i)
328 
329  call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
330  at4_1, at4_2, at5, at6, at7, atr1, &
331  am1, am2, alb1, &
332  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
333  ai1, ai2, ai3, dzeta_t, &
334  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
335 
336  if ( (temp_c_neu(0,j,i)-(-beta*h_c_neu(j,i))).gt.0.0_dp ) &
337  then
338  call shift_cts_upward(at1, at2_1, at2_2, at3_1, at3_2, &
339  at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
340  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
341  ai1, ai2, ai3, dzeta_t, &
342  dtime_temp, dtt_2dxi, dtt_2deta, dtime_temp_inv, &
343  i, j)
344  else
345  call shift_cts_downward(at1, at2_1, at2_2, at3_1, at3_2, &
346  at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
347  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
348  ai1, ai2, ai3, dzeta_t, &
349  dtime_temp, dtt_2dxi, dtt_2deta, dtime_temp_inv, &
350  i, j)
351  end if
352 
353  end if
354 
355 #if (MARGIN==3)
356 
357  else if (maske(j,i)==3) then ! floating ice
358 
359  n_cts_neu(j,i) = -1
360  zm_neu(j,i) = zb(j,i)
361  h_c_neu(j,i) = h_c(j,i) + h_t(j,i)
362  h_t_neu(j,i) = 0.0_dp
363 
364  call calc_temp_ssa(at1, at2_1, at2_2, at3_1, at3_2, &
365  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
366  ai1, ai2, &
367  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
368 
369 ! ------ Reset temperatures above melting to the melting point
370 ! (should not occur, but just in case)
371 
372  do kc=0, kcmax
373  if (temp_c_neu(kc,j,i) > temp_c_m(kc,j,i)) &
374  temp_c_neu(kc,j,i) = temp_c_m(kc,j,i)
375  end do
376 
377 #endif
378 
379  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
380 
381  n_cts_neu(j,i) = -1
382  zm_neu(j,i) = zb(j,i)
383  h_c_neu(j,i) = h_c(j,i)
384  h_t_neu(j,i) = h_t(j,i)
385 
386  call calc_temp_r(atr1, alb1, i, j)
387 
388 endif
389 
390 end do
391 end do ! End of computation loop
392 
393 !-------- Extrapolate values on margins --------
394 
395 ! ------ Lower left corner
396 
397 i=0
398 j=0
399 
400 if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
401  ! glaciated land or floating ice
402  ii=i+1
403  jj=j+1
404 
405  do kc=0,kcmax
406  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
407  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
408  end do
409 
410  do kt=0,ktmax
411  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii) ! set temp.-ice water content
412  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii) ! set temp.-ice age
413  end do
414 
415  do kr=0,krmax
416  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
417  end do
418 
419  n_cts_neu(j,i) = min(n_cts_neu(jj,ii),0) ! temperate layer excluded
420  h_c_neu(j,i) = h_c(j,i)
421  h_t_neu(j,i) = h_t(j,i)
422  zm_neu(j,i) = zb(j,i)
423 
424 else ! maske(j,i).eq.1,2 (ice-free land or sea point)
425 
426  n_cts_neu(j,i) = -1
427  zm_neu(j,i) = zb(j,i)
428  h_c_neu(j,i) = h_c(j,i)
429  h_t_neu(j,i) = h_t(j,i)
430 
431  call calc_temp_r(atr1, alb1, i, j)
432 
433 end if
434 
435 ! ------ Lower right corner
436 
437 i=imax
438 j=0
439 
440 if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
441  ! glaciated land or floating ice
442  ii=i-1
443  jj=j+1
444 
445  do kc=0,kcmax
446  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
447  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
448  end do
449 
450  do kt=0,ktmax
451  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii) ! set temp.-ice water content
452  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii) ! set temp.-ice age
453  end do
454 
455  do kr=0,krmax
456  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
457  end do
458 
459  n_cts_neu(j,i) = min(n_cts_neu(jj,ii),0) ! temperate layer excluded
460  h_c_neu(j,i) = h_c(j,i)
461  h_t_neu(j,i) = h_t(j,i)
462  zm_neu(j,i) = zb(j,i)
463 
464 else ! maske(j,i).eq.1,2 (ice-free land or sea point)
465 
466  n_cts_neu(j,i) = -1
467  zm_neu(j,i) = zb(j,i)
468  h_c_neu(j,i) = h_c(j,i)
469  h_t_neu(j,i) = h_t(j,i)
470 
471  call calc_temp_r(atr1, alb1, i, j)
472 
473 end if
474 
475 ! ------ Upper left corner
476 
477 i=0
478 j=jmax
479 
480 if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
481  ! glaciated land or floating ice
482  ii=i+1
483  jj=j-1
484 
485  do kc=0,kcmax
486  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
487  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
488  end do
489 
490  do kt=0,ktmax
491  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii) ! set temp.-ice water content
492  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii) ! set temp.-ice age
493  end do
494 
495  do kr=0,krmax
496  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
497  end do
498 
499  n_cts_neu(j,i) = min(n_cts_neu(jj,ii),0) ! temperate layer excluded
500  h_c_neu(j,i) = h_c(j,i)
501  h_t_neu(j,i) = h_t(j,i)
502  zm_neu(j,i) = zb(j,i)
503 
504 else ! maske(j,i).eq.1,2 (ice-free land or sea point)
505 
506  n_cts_neu(j,i) = -1
507  zm_neu(j,i) = zb(j,i)
508  h_c_neu(j,i) = h_c(j,i)
509  h_t_neu(j,i) = h_t(j,i)
510 
511  call calc_temp_r(atr1, alb1, i, j)
512 
513 end if
514 
515 ! ------ Upper right corner
516 
517 i=imax
518 j=jmax
519 
520 if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
521  ! glaciated land or floating ice
522  ii=i-1
523  jj=j-1
524 
525  do kc=0,kcmax
526  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
527  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
528  end do
529 
530  do kt=0,ktmax
531  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii) ! set temp.-ice water content
532  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii) ! set temp.-ice age
533  end do
534 
535  do kr=0,krmax
536  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
537  end do
538 
539  n_cts_neu(j,i) = min(n_cts_neu(jj,ii),0) ! temperate layer excluded
540  h_c_neu(j,i) = h_c(j,i)
541  h_t_neu(j,i) = h_t(j,i)
542  zm_neu(j,i) = zb(j,i)
543 
544 else ! maske(j,i).eq.1,2 (ice-free land or sea point)
545 
546  n_cts_neu(j,i) = -1
547  zm_neu(j,i) = zb(j,i)
548  h_c_neu(j,i) = h_c(j,i)
549  h_t_neu(j,i) = h_t(j,i)
550 
551  call calc_temp_r(atr1, alb1, i, j)
552 
553 end if
554 
555 ! ------ Lower and upper margins
556 
557 do i=1, imax-1
558 
559 ! ---- Lower margin
560 
561  j=0
562 
563  if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
564  ! glaciated land or floating ice
565  ii=i
566  jj=j+1
567 
568  do kc=0,kcmax
569  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
570  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
571  end do
572 
573  do kt=0,ktmax
574  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii) ! set temp.-ice water content
575  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii) ! set temp.-ice age
576  end do
577 
578  do kr=0,krmax
579  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
580  end do
581 
582  n_cts_neu(j,i) = min(n_cts_neu(jj,ii),0) ! temperate layer excluded
583  h_c_neu(j,i) = h_c(j,i)
584  h_t_neu(j,i) = h_t(j,i)
585  zm_neu(j,i) = zb(j,i)
586 
587  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
588 
589  n_cts_neu(j,i) = -1
590  zm_neu(j,i) = zb(j,i)
591  h_c_neu(j,i) = h_c(j,i)
592  h_t_neu(j,i) = h_t(j,i)
593 
594  call calc_temp_r(atr1, alb1, i, j)
595 
596  end if
597 
598 ! ---- Upper margin
599 
600  j=jmax
601 
602  if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
603  ! glaciated land or floating ice
604  ii=i
605  jj=j-1
606 
607  do kc=0,kcmax
608  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
609  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
610  end do
611 
612  do kt=0,ktmax
613  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii) ! set temp.-ice water content
614  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii) ! set temp.-ice age
615  end do
616 
617  do kr=0,krmax
618  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
619  end do
620 
621  n_cts_neu(j,i) = min(n_cts_neu(jj,ii),0) ! temperate layer excluded
622  h_c_neu(j,i) = h_c(j,i)
623  h_t_neu(j,i) = h_t(j,i)
624  zm_neu(j,i) = zb(j,i)
625 
626  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
627 
628  n_cts_neu(j,i) = -1
629  zm_neu(j,i) = zb(j,i)
630  h_c_neu(j,i) = h_c(j,i)
631  h_t_neu(j,i) = h_t(j,i)
632 
633  call calc_temp_r(atr1, alb1, i, j)
634 
635  end if
636 
637 end do
638 
639 ! ------ Left and right margins
640 
641 do j=1, jmax-1
642 
643 ! ---- Left margin
644 
645  i=0
646 
647  if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
648  ! glaciated land or floating ice
649  ii=i+1
650  jj=j
651 
652  do kc=0,kcmax
653  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
654  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
655  end do
656 
657  do kt=0,ktmax
658  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii) ! set temp.-ice water content
659  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii) ! set temp.-ice age
660  end do
661 
662  do kr=0,krmax
663  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
664  end do
665 
666  n_cts_neu(j,i) = min(n_cts_neu(jj,ii),0) ! temperate layer excluded
667  h_c_neu(j,i) = h_c(j,i)
668  h_t_neu(j,i) = h_t(j,i)
669  zm_neu(j,i) = zb(j,i)
670 
671  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
672 
673  n_cts_neu(j,i) = -1
674  zm_neu(j,i) = zb(j,i)
675  h_c_neu(j,i) = h_c(j,i)
676  h_t_neu(j,i) = h_t(j,i)
677 
678  call calc_temp_r(atr1, alb1, i, j)
679 
680  end if
681 
682 ! ---- Right margin
683 
684  i=imax
685 
686  if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
687  ! glaciated land or floating ice
688  ii=i-1
689  jj=j
690 
691  do kc=0,kcmax
692  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
693  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
694  end do
695 
696  do kt=0,ktmax
697  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii) ! set temp.-ice water content
698  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii) ! set temp.-ice age
699  end do
700 
701  do kr=0,krmax
702  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
703  end do
704 
705  n_cts_neu(j,i) = min(n_cts_neu(jj,ii),0) ! temperate layer excluded
706  h_c_neu(j,i) = h_c(j,i)
707  h_t_neu(j,i) = h_t(j,i)
708  zm_neu(j,i) = zb(j,i)
709 
710  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
711 
712  n_cts_neu(j,i) = -1
713  zm_neu(j,i) = zb(j,i)
714  h_c_neu(j,i) = h_c(j,i)
715  h_t_neu(j,i) = h_t(j,i)
716 
717  call calc_temp_r(atr1, alb1, i, j)
718 
719  end if
720 
721 end do
722 
723 !-------- Dummy values for omega_c_neu and kc_cts_neu --------
724 
725 omega_c_neu = 0.0_dp ! not needed for
726 kc_cts_neu = 0 ! the polythermal mode
727 
728 !-------- Smoothing of H_t_neu with numerical diffusion --------
729 
730 ! ------ Volume of temperate ice without smoothing
731 
732 vol_t = 0.0_dp
733 do i=0, imax ! extended to domain margins (22.1.02 -> V1.1)
734 do j=0, jmax ! extended to domain margins (22.1.02 -> V1.1)
735  if (n_cts_neu(j,i).eq.1) then
736  vol_t = vol_t + h_t_neu(j,i)*area(j,i)
737  end if
738 end do
739 end do
740 
741 ! ------ Smoothing
742 
743 do i=1, imax-1
744 do j=1, jmax-1
745  if (n_cts_neu(j,i).ne.-1) then
746 
747  dh_t_smooth(j,i) = numdiff_h_t* ( -4.0_dp*h_t_neu(j,i) &
748  +h_t_neu(j,i+1)+h_t_neu(j,i-1) &
749  +h_t_neu(j+1,i)+h_t_neu(j-1,i) )
750  if (dh_t_smooth(j,i).gt.0.001_dp) n_cts_neu(j,i) = 1
751 
752  end if
753 end do
754 end do
755 
756 do i=1, imax-1
757 do j=1, jmax-1
758  if (n_cts_neu(j,i).eq.1) then
759  h_t_neu(j,i) = h_t_neu(j,i) + dh_t_smooth(j,i)
760  end if
761 end do
762 end do
763 
764 ! ------ Volume of temperate ice with smoothing
765 
766 vol_t_smooth = 0.0_dp
767 do i=0, imax ! extended to domain margins (22.1.02 -> V1.1)
768 do j=0, jmax ! extended to domain margins (22.1.02 -> V1.1)
769  if (n_cts_neu(j,i).eq.1) then
770  vol_t_smooth = vol_t_smooth + h_t_neu(j,i)*area(j,i)
771  end if
772 end do
773 end do
774 
775 ! ------ Correction so that volume is not changed by the smoothing
776 
777 if (vol_t_smooth.gt.0.0_dp) then
778 
779  korrfakt_t = vol_t/vol_t_smooth
780  do i=0, imax ! extended to domain margins (22.1.02 -> V1.1)
781  do j=0, jmax ! extended to domain margins (22.1.02 -> V1.1)
782  if (n_cts_neu(j,i).eq.1) then
783  h_t_neu(j,i) = h_t_neu(j,i)*korrfakt_t
784 ! zm_neu(j,i) = zb(j,i) + H_t_neu(j,i)
785 ! H_c_neu(j,i) = zs(j,i) - zm_neu(j,i)
786  end if
787  end do
788  end do
789 
790 end if
791 
792 !-------- Numerical time lag for evolution of H_t_neu --------
793 
794 time_lag_cts = tau_cts*year_sec ! yr --> s
795 
796 do i=0, imax ! extended to domain margins (22.1.02 -> V1.1)
797 do j=0, jmax ! extended to domain margins (22.1.02 -> V1.1)
798 
799  if (n_cts_neu(j,i).eq.1) then
800 
801  h_t_neu(j,i) = ( time_lag_cts*h_t(j,i) &
802  + dtime_temp*h_t_neu(j,i) ) &
803  /(time_lag_cts+dtime_temp)
804 
805  zm_neu(j,i) = zb(j,i) + h_t_neu(j,i)
806  h_c_neu(j,i) = zs(j,i) - zm_neu(j,i)
807 
808  end if
809 
810 end do
811 end do
812 
813 end subroutine calc_temp_poly
814 
815 !-------------------------------------------------------------------------------
816 !> Computation of temperature and age in cold-ice mode.
817 !<------------------------------------------------------------------------------
818 subroutine calc_temp_cold(dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
819  dtime_temp)
820 
821 use ice_material_properties_m, only : kappa_val
822 
823 implicit none
824 
825 real(dp), intent(in) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
826 real(dp), intent(in) :: dtime_temp
827 
828 integer(i4b) :: i, j, kc, kr, ii, jj
829 real(dp) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
830  at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
831  at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
832  acb1, acb2, acb3, acb4, &
833  ai1(0:kcmax), ai2(0:kcmax), ai3, &
834  atr1, am1, am2, alb1
835 real(dp) :: aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld
836 real(dp) :: dtime_temp_inv, dtt_2dxi, dtt_2deta
837 real(dp) :: temp_c_help(0:kcmax)
838 
839 !-------- Term abbreviations
840 
841 at7 = 2.0_dp/rho*dtime_temp
842 
843 aw1 = dtime_temp/dzeta_t
844 aw2 = dtime_temp/dzeta_t
845 aw3 = dtime_temp/dzeta_t
846 aw4 = dtime_temp/dzeta_t
847 aw5 = nue/rho*dtime_temp/(dzeta_t**2)
848 aw7 = 2.0_dp/(rho*l)*dtime_temp
849 aw8 = beta**2/(rho*l) &
850  *(kappa_val(0.0_dp)-kappa_val(-1.0_dp))*dtime_temp
851 aw9 = beta/l*dtime_temp
852 
853 ai3 = agediff*dtime_temp/(dzeta_t**2)
854 
855 atr1 = kappa_r/(rho_c_r*h_r**2)*dtime_temp/(dzeta_r**2)
856 
857 if (flag_aa_nonzero) then
858  am1 = aa*beta*dzeta_c/(ea-1.0_dp)
859  am2 = aa*l*rho*dzeta_c/(ea-1.0_dp)
860 else
861  am1 = beta*dzeta_c
862  am2 = l*rho*dzeta_c
863 end if
864 
865 if (flag_aa_nonzero) then
866  acb1 = (ea-1.0_dp)/aa/dzeta_c
867 else
868  acb1 = 1.0_dp/dzeta_c
869 end if
870 
871 acb2 = kappa_r/h_r/dzeta_r
872 acb3 = rho*g
873 acb4 = rho*g
874 
875 alb1 = h_r/kappa_r*dzeta_r
876 
877 aqtld = dzeta_t/dtime_temp
878 
879 dtt_2dxi = 0.5_dp*dtime_temp/dxi
880 dtt_2deta = 0.5_dp*dtime_temp/deta
881 
882 dtime_temp_inv = 1.0_dp/dtime_temp
883 
884 do kc=0, kcmax
885 
886  if (flag_aa_nonzero) then
887 
888  at1(kc) = (ea-1.0_dp)/(aa*eaz_c(kc))*dtime_temp/dzeta_c
889  at2_1(kc) = (ea-1.0_dp)/(aa*eaz_c(kc))*dtime_temp/dzeta_c
890  at2_2(kc) = (eaz_c(kc)-1.0_dp)/(aa*eaz_c(kc)) &
891  *dtime_temp/dzeta_c
892  at3_1(kc) = (ea-1.0_dp)/(aa*eaz_c(kc))*dtime_temp/dzeta_c
893  at3_2(kc) = (eaz_c(kc)-1.0_dp)/(aa*eaz_c(kc)) &
894  *dtime_temp/dzeta_c
895  at4_1(kc) = (ea-1.0_dp)/(aa*eaz_c(kc))*dtime_temp/dzeta_c
896  at4_2(kc) = (eaz_c(kc)-1.0_dp)/(aa*eaz_c(kc)) &
897  *dtime_temp/dzeta_c
898  at5(kc) = (ea-1.0_dp)/(rho*aa*eaz_c(kc)) &
899  *dtime_temp/dzeta_c
900  if (kc /= kcmax) then
901  at6(kc) = (ea-1.0_dp) &
902  /(aa*exp(aa*0.5_dp*(zeta_c(kc)+zeta_c(kc+1)))) &
903  /dzeta_c
904  else
905  at6(kc) = 0.0_dp
906  end if
907  ai1(kc) = agediff*(ea-1.0_dp)/(aa*eaz_c(kc)) &
908  *dtime_temp/dzeta_c
909  if (kc /= kcmax) then
910  ai2(kc) = (ea-1.0_dp) &
911  /(aa*exp(aa*0.5_dp*(zeta_c(kc)+zeta_c(kc+1)))) &
912  /dzeta_c
913  else
914  ai2(kc) = 0.0_dp
915  end if
916 
917  else
918 
919  at1(kc) = dtime_temp/dzeta_c
920  at2_1(kc) = dtime_temp/dzeta_c
921  at2_2(kc) = zeta_c(kc) &
922  *dtime_temp/dzeta_c
923  at3_1(kc) = dtime_temp/dzeta_c
924  at3_2(kc) = zeta_c(kc) &
925  *dtime_temp/dzeta_c
926  at4_1(kc) = dtime_temp/dzeta_c
927  at4_2(kc) = zeta_c(kc) &
928  *dtime_temp/dzeta_c
929  at5(kc) = 1.0_dp/rho &
930  *dtime_temp/dzeta_c
931  if (kc /= kcmax) then
932  at6(kc) = 1.0_dp &
933  /dzeta_c
934  else
935  at6(kc) = 0.0_dp
936  end if
937  ai1(kc) = agediff &
938  *dtime_temp/dzeta_c
939  if (kc /= kcmax) then
940  ai2(kc) = 1.0_dp &
941  /dzeta_c
942  else
943  ai2(kc) = 0.0_dp
944  end if
945 
946  end if
947 
948 end do
949 
950 !-------- Computation loop for temperature and age --------
951 
952 do i=1, imax-1 ! skipping domain margins
953 do j=1, jmax-1 ! skipping domain margins
954 
955  if (maske(j,i)==0) then ! glaciated land
956 
957  n_cts_neu(j,i) = -1
958  zm_neu(j,i) = zb(j,i)
959  h_c_neu(j,i) = h_c(j,i)
960  h_t_neu(j,i) = h_t(j,i)
961 
962  call calc_temp1(at1, at2_1, at2_2, at3_1, at3_2, &
963  at4_1, at4_2, at5, at6, at7, atr1, acb1, acb2, &
964  acb3, acb4, alb1, ai1, ai2, &
965  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
966 
967 ! ------ Reset temperatures above melting to the melting point,
968 ! look for the CTS
969 
970  kc_cts_neu(j,i) = 0
971 
972  if (temp_c_neu(0,j,i).gt.temp_c_m(0,j,i)) then
973  n_cts_neu(j,i) = 0
974  kc_cts_neu(j,i) = 0
975  temp_c_neu(0,j,i) = temp_c_m(0,j,i)
976  temp_r_neu(krmax,j,i) = temp_c_m(0,j,i)
977  end if
978 
979  do kc=1, kcmax
980  if (temp_c_neu(kc,j,i).gt.temp_c_m(kc,j,i)) then
981  kc_cts_neu(j,i) = kc
982  temp_c_neu(kc,j,i) = temp_c_m(kc,j,i)
983  end if
984  end do
985 
986 #if (MARGIN==3)
987 
988  else if (maske(j,i)==3) then ! floating ice
989 
990  n_cts_neu(j,i) = -1
991  kc_cts_neu(j,i) = 0
992  zm_neu(j,i) = zb(j,i)
993  h_c_neu(j,i) = h_c(j,i)
994  h_t_neu(j,i) = 0.0_dp
995 
996  call calc_temp_ssa(at1, at2_1, at2_2, at3_1, at3_2, &
997  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
998  ai1, ai2, &
999  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
1000 
1001 ! ------ Reset temperatures above melting to the melting point
1002 ! (should not occur, but just in case)
1003 
1004  do kc=0, kcmax
1005  if (temp_c_neu(kc,j,i) > temp_c_m(kc,j,i)) &
1006  temp_c_neu(kc,j,i) = temp_c_m(kc,j,i)
1007  end do
1008 
1009 #endif
1010 
1011  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
1012 
1013  n_cts_neu(j,i) = -1
1014  kc_cts_neu(j,i) = 0
1015  zm_neu(j,i) = zb(j,i)
1016  h_c_neu(j,i) = h_c(j,i)
1017  h_t_neu(j,i) = h_t(j,i)
1018 
1019  call calc_temp_r(atr1, alb1, i, j)
1020 
1021  end if
1022 
1023 end do
1024 end do
1025 
1026 !-------- Extrapolate values on margins --------
1027 
1028 ! ------ Lower left corner
1029 
1030 i=0
1031 j=0
1032 
1033 if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
1034  ! glaciated land or floating ice
1035  ii=i+1
1036  jj=j+1
1037 
1038  do kc=0,kcmax
1039  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
1040  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
1041  end do
1042 
1043  do kr=0,krmax
1044  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
1045  end do
1046 
1047  n_cts_neu(j,i) = n_cts_neu(jj,ii)
1048  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
1049  zm_neu(j,i) = zb(j,i)
1050  h_c_neu(j,i) = h_c(j,i)
1051  h_t_neu(j,i) = h_t(j,i)
1052 
1053 else ! maske(j,i).eq.1,2 (ice-free land or sea point)
1054 
1055  n_cts_neu(j,i) = -1
1056  kc_cts_neu(j,i) = 0
1057  zm_neu(j,i) = zb(j,i)
1058  h_c_neu(j,i) = h_c(j,i)
1059  h_t_neu(j,i) = h_t(j,i)
1060 
1061  call calc_temp_r(atr1, alb1, i, j)
1062 
1063 end if
1064 
1065 ! ------ Lower right corner
1066 
1067 i=imax
1068 j=0
1069 
1070 if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
1071  ! glaciated land or floating ice
1072  ii=i-1
1073  jj=j+1
1074 
1075  do kc=0,kcmax
1076  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
1077  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
1078  end do
1079 
1080  do kr=0,krmax
1081  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
1082  end do
1083 
1084  n_cts_neu(j,i) = n_cts_neu(jj,ii)
1085  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
1086  zm_neu(j,i) = zb(j,i)
1087  h_c_neu(j,i) = h_c(j,i)
1088  h_t_neu(j,i) = h_t(j,i)
1089 
1090 else ! maske(j,i).eq.1,2 (ice-free land or sea point)
1091 
1092  n_cts_neu(j,i) = -1
1093  kc_cts_neu(j,i) = 0
1094  zm_neu(j,i) = zb(j,i)
1095  h_c_neu(j,i) = h_c(j,i)
1096  h_t_neu(j,i) = h_t(j,i)
1097 
1098  call calc_temp_r(atr1, alb1, i, j)
1099 
1100 end if
1101 
1102 ! ------ Upper left corner
1103 
1104 i=0
1105 j=jmax
1106 
1107 if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
1108  ! glaciated land or floating ice
1109  ii=i+1
1110  jj=j-1
1111 
1112  do kc=0,kcmax
1113  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
1114  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
1115  end do
1116 
1117  do kr=0,krmax
1118  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
1119  end do
1120 
1121  n_cts_neu(j,i) = n_cts_neu(jj,ii)
1122  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
1123  zm_neu(j,i) = zb(j,i)
1124  h_c_neu(j,i) = h_c(j,i)
1125  h_t_neu(j,i) = h_t(j,i)
1126 
1127 else ! maske(j,i).eq.1,2 (ice-free land or sea point)
1128 
1129  n_cts_neu(j,i) = -1
1130  kc_cts_neu(j,i) = 0
1131  zm_neu(j,i) = zb(j,i)
1132  h_c_neu(j,i) = h_c(j,i)
1133  h_t_neu(j,i) = h_t(j,i)
1134 
1135  call calc_temp_r(atr1, alb1, i, j)
1136 
1137 end if
1138 
1139 ! ------ Upper right corner
1140 
1141 i=imax
1142 j=jmax
1143 
1144 if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
1145  ! glaciated land or floating ice
1146  ii=i-1
1147  jj=j-1
1148 
1149  do kc=0,kcmax
1150  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
1151  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
1152  end do
1153 
1154  do kr=0,krmax
1155  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
1156  end do
1157 
1158  n_cts_neu(j,i) = n_cts_neu(jj,ii)
1159  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
1160  zm_neu(j,i) = zb(j,i)
1161  h_c_neu(j,i) = h_c(j,i)
1162  h_t_neu(j,i) = h_t(j,i)
1163 
1164 else ! maske(j,i).eq.1,2 (ice-free land or sea point)
1165 
1166  n_cts_neu(j,i) = -1
1167  kc_cts_neu(j,i) = 0
1168  zm_neu(j,i) = zb(j,i)
1169  h_c_neu(j,i) = h_c(j,i)
1170  h_t_neu(j,i) = h_t(j,i)
1171 
1172  call calc_temp_r(atr1, alb1, i, j)
1173 
1174 end if
1175 
1176 ! ------ Lower and upper margins
1177 
1178 do i=1, imax-1
1179 
1180 ! ---- Lower margin
1181 
1182  j=0
1183 
1184  if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
1185  ! glaciated land or floating ice
1186  ii=i
1187  jj=j+1
1188 
1189  do kc=0,kcmax
1190  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
1191  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
1192  end do
1193 
1194  do kr=0,krmax
1195  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
1196  end do
1197 
1198  n_cts_neu(j,i) = n_cts_neu(jj,ii)
1199  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
1200  zm_neu(j,i) = zb(j,i)
1201  h_c_neu(j,i) = h_c(j,i)
1202  h_t_neu(j,i) = h_t(j,i)
1203 
1204  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
1205 
1206  n_cts_neu(j,i) = -1
1207  kc_cts_neu(j,i) = 0
1208  zm_neu(j,i) = zb(j,i)
1209  h_c_neu(j,i) = h_c(j,i)
1210  h_t_neu(j,i) = h_t(j,i)
1211 
1212  call calc_temp_r(atr1, alb1, i, j)
1213 
1214  end if
1215 
1216 ! ---- Upper margin
1217 
1218  j=jmax
1219 
1220  if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
1221  ! glaciated land or floating ice
1222  ii=i
1223  jj=j-1
1224 
1225  do kc=0,kcmax
1226  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
1227  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
1228  end do
1229 
1230  do kr=0,krmax
1231  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
1232  end do
1233 
1234  n_cts_neu(j,i) = n_cts_neu(jj,ii)
1235  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
1236  zm_neu(j,i) = zb(j,i)
1237  h_c_neu(j,i) = h_c(j,i)
1238  h_t_neu(j,i) = h_t(j,i)
1239 
1240  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
1241 
1242  n_cts_neu(j,i) = -1
1243  kc_cts_neu(j,i) = 0
1244  zm_neu(j,i) = zb(j,i)
1245  h_c_neu(j,i) = h_c(j,i)
1246  h_t_neu(j,i) = h_t(j,i)
1247 
1248  call calc_temp_r(atr1, alb1, i, j)
1249 
1250  end if
1251 
1252 end do
1253 
1254 ! ------ Left and right margins
1255 
1256 do j=1, jmax-1
1257 
1258 ! ---- Left margin
1259 
1260  i=0
1261 
1262  if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
1263  ! glaciated land or floating ice
1264  ii=i+1
1265  jj=j
1266 
1267  do kc=0,kcmax
1268  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
1269  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
1270  end do
1271 
1272  do kr=0,krmax
1273  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
1274  end do
1275 
1276  n_cts_neu(j,i) = n_cts_neu(jj,ii)
1277  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
1278  zm_neu(j,i) = zb(j,i)
1279  h_c_neu(j,i) = h_c(j,i)
1280  h_t_neu(j,i) = h_t(j,i)
1281 
1282  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
1283 
1284  n_cts_neu(j,i) = -1
1285  kc_cts_neu(j,i) = 0
1286  zm_neu(j,i) = zb(j,i)
1287  h_c_neu(j,i) = h_c(j,i)
1288  h_t_neu(j,i) = h_t(j,i)
1289 
1290  call calc_temp_r(atr1, alb1, i, j)
1291 
1292  end if
1293 
1294 ! ---- Right margin
1295 
1296  i=imax
1297 
1298  if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
1299  ! glaciated land or floating ice
1300  ii=i-1
1301  jj=j
1302 
1303  do kc=0,kcmax
1304  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
1305  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
1306  end do
1307 
1308  do kr=0,krmax
1309  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
1310  end do
1311 
1312  n_cts_neu(j,i) = n_cts_neu(jj,ii)
1313  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
1314  zm_neu(j,i) = zb(j,i)
1315  h_c_neu(j,i) = h_c(j,i)
1316  h_t_neu(j,i) = h_t(j,i)
1317 
1318  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
1319 
1320  n_cts_neu(j,i) = -1
1321  kc_cts_neu(j,i) = 0
1322  zm_neu(j,i) = zb(j,i)
1323  h_c_neu(j,i) = h_c(j,i)
1324  h_t_neu(j,i) = h_t(j,i)
1325 
1326  call calc_temp_r(atr1, alb1, i, j)
1327 
1328  end if
1329 
1330 end do
1331 
1332 !-------- Dummy values for omega_c_neu --------
1333 
1334 omega_c_neu = 0.0_dp ! not computed in the cold-ice mode
1335 
1336 end subroutine calc_temp_cold
1337 
1338 !-------------------------------------------------------------------------------
1339 !> Isothermal mode: Setting of the temperature and age to constant values.
1340 !<------------------------------------------------------------------------------
1341 subroutine calc_temp_const()
1342 
1343 implicit none
1344 
1345 #if (defined(TEMP_CONST))
1346  if ( temp_const > -eps ) &
1347  stop ' >>> calc_temp_const: TEMP_CONST must be negative!'
1348  temp_c_neu = temp_const
1349  temp_r_neu = temp_const
1350 #else
1351  temp_c_neu = -10.0_dp ! default value -10 C
1352  temp_r_neu = -10.0_dp ! default value -10 C
1353 #endif
1354 
1355 temp_c_neu = min(temp_c_neu, temp_c_m-eps)
1356  ! keep temperatures below the pressure melting point
1357 
1358 omega_t_neu = 0.0_dp
1359 omega_c_neu = 0.0_dp
1360 
1361 q_tld = 0.0_dp
1362 
1363 #if (defined(AGE_CONST))
1364  age_c_neu = age_const *year_sec ! a --> s
1365  age_t_neu = age_const *year_sec ! a --> s
1366 #else
1367  age_c_neu = 0.0_dp ! default value 0
1368  age_t_neu = 0.0_dp ! default value 0
1369 #endif
1370 
1371 n_cts_neu = -1
1372 kc_cts_neu = 0
1373 zm_neu = zb
1374 h_c_neu = h_c
1375 h_t_neu = 0.0_dp
1376 
1377 end subroutine calc_temp_const
1378 
1379 !-------------------------------------------------------------------------------
1380 !> Computation of temperature and age for a cold ice column.
1381 !<------------------------------------------------------------------------------
1382 subroutine calc_temp1(at1, at2_1, at2_2, at3_1, at3_2, &
1383  at4_1, at4_2, at5, at6, at7, atr1, acb1, acb2, &
1384  acb3, acb4, alb1, ai1, ai2, &
1385  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
1386 
1387 use ice_material_properties_m, only : ratefac_c, kappa_val, c_val, &
1388  creep, viscosity
1389 use sico_maths_m, only : tri_sle
1390 
1391 implicit none
1392 
1393 integer(i4b), intent(in) :: i, j
1394 real(dp), intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
1395  at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
1396  at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
1397  ai1(0:kcmax), ai2(0:kcmax), &
1398  atr1, acb1, acb2, acb3, acb4, alb1
1399 real(dp), intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
1400 
1401 integer(i4b) :: kc, kt, kr
1402 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
1403  ct5(0:kcmax), ct6(0:kcmax), ct7(0:kcmax), ctr1, &
1404  ccb1, ccb2, ccb3, ccb4, clb1
1405 real(dp) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), ct3_sg(0:kcmax), &
1406  ct4_sg(0:kcmax), adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
1407 real(dp) :: ci1(0:kcmax), ci2(0:kcmax)
1408 real(dp) :: temp_c_help(0:kcmax)
1409 real(dp) :: vx_c_help, vy_c_help
1410 real(dp) :: adv_vert_help
1411 real(dp) :: dtt_dxi, dtt_deta
1412 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
1413  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
1414  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
1415  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
1416  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
1417 real(dp), parameter :: zero=0.0_dp
1418 
1419 !-------- Check for boundary points --------
1420 
1421 if ((i == 0).or.(i == imax).or.(j == 0).or.(j == jmax)) &
1422  stop ' >>> calc_temp1: Boundary points not allowed.'
1423 
1424 !-------- Abbreviations --------
1425 
1426 ctr1 = atr1
1427 
1428 ccb1 = acb1 &
1429  *kappa_val(temp_c(0,j,i)) &
1430  /h_c(j,i)
1431 ccb2 = acb2
1432 
1433 #if (DYNAMICS==2)
1434 if (.not.flag_shelfy_stream(j,i)) then
1435 #endif
1436 
1437  ccb3 = acb3*0.5_dp*(vx_t(0,j,i)+vx_t(0,j,i-1)) &
1438  *h_c(j,i)*dzs_dxi_g(j,i)
1439  ccb4 = acb4*0.5_dp*(vy_t(0,j,i)+vy_t(0,j-1,i)) &
1440  *h_c(j,i)*dzs_deta_g(j,i)
1441 
1442 #if (DYNAMICS==2)
1443 else ! flag_shelfy_stream(j,i) == .true.
1444 
1445  ccb3 = -c_drag(j,i) &
1446  * sqrt(vx_b_g(j,i)**2 &
1447  +vy_b_g(j,i)**2) &
1448  **(1.0_dp+p_weert_inv(j,i))
1449  ccb4 = 0.0_dp
1450 
1451 end if
1452 #endif
1453 
1454 clb1 = alb1*q_geo(j,i)
1455 
1456 #if (ADV_VERT==1)
1457 
1458 do kc=1, kcmax-1
1459  ct1(kc) = at1(kc)/h_c(j,i)*0.5_dp*(vz_c(kc,j,i)+vz_c(kc-1,j,i))
1460 end do
1461 
1462 kc=0
1463 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
1464  ! only needed for kc=0 ...
1465 kc=kcmax-1
1466 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
1467  ! ... and kc=KCMAX-1
1468 
1469 #elif (ADV_VERT==2 || ADV_VERT==3)
1470 
1471 do kc=0, kcmax-1
1472  ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
1473 end do
1474 
1475 #endif
1476 
1477 do kc=0, kcmax
1478 
1479  ct2(kc) = ( at2_1(kc)*dzm_dtau(j,i) &
1480  +at2_2(kc)*dh_c_dtau(j,i) )/h_c(j,i)
1481  ct3(kc) = ( at3_1(kc)*dzm_dxi_g(j,i) &
1482  +at3_2(kc)*dh_c_dxi_g(j,i) )/h_c(j,i) &
1483  *0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1)) *insq_g11_g(j,i)
1484  ct4(kc) = ( at4_1(kc)*dzm_deta_g(j,i) &
1485  +at4_2(kc)*dh_c_deta_g(j,i) )/h_c(j,i) &
1486  *0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i)) *insq_g22_g(j,i)
1487  ct5(kc) = at5(kc) &
1488  /c_val(temp_c(kc,j,i)) &
1489  /h_c(j,i)
1490 
1491 #if (DYNAMICS==2)
1492  if (.not.flag_shelfy_stream(j,i)) then
1493 #endif
1494  ct7(kc) = at7 &
1495  /c_val(temp_c(kc,j,i)) &
1496  *enh_c(kc,j,i) &
1497  *ratefac_c(temp_c(kc,j,i), temp_c_m(kc,j,i)) &
1498  *creep(sigma_c(kc,j,i)) &
1499  *sigma_c(kc,j,i)*sigma_c(kc,j,i)
1500 #if (DYNAMICS==2)
1501  else
1502  ct7(kc) = 2.0_dp*at7 &
1503  /c_val(temp_c(kc,j,i)) &
1504  *viscosity(de_c(kc,j,i), &
1505  temp_c(kc,j,i), temp_c_m(kc,j,i), 0.0_dp, &
1506  enh_c(kc,j,i), 0_i2b) &
1507  *de_c(kc,j,i)**2
1508  end if
1509 #endif
1510 
1511  ci1(kc) = ai1(kc)/h_c(j,i)
1512 
1513 end do
1514 
1515 #if (ADV_VERT==1)
1516 
1517 kc=0
1518 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
1519 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
1520 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
1521 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
1522 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! only needed for kc=0 ...
1523 kc=kcmax-1
1524 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
1525 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
1526 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
1527 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
1528 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! ... and kc=KCMAX-1
1529 
1530 #elif (ADV_VERT==2 || ADV_VERT==3)
1531 
1532 do kc=0, kcmax-1
1533  ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
1534  ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
1535  ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
1536  adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
1537  abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
1538 end do
1539 
1540 #endif
1541 
1542 do kc=0, kcmax-1
1543  temp_c_help(kc) = 0.5_dp*(temp_c(kc,j,i)+temp_c(kc+1,j,i))
1544  ct6(kc) = at6(kc) &
1545  *kappa_val(temp_c_help(kc)) &
1546  /h_c(j,i)
1547  ci2(kc) = ai2(kc)/h_c(j,i)
1548 end do
1549 
1550 #if (ADV_HOR==3)
1551 dtt_dxi = 2.0_dp*dtt_2dxi
1552 dtt_deta = 2.0_dp*dtt_2deta
1553 #endif
1554 
1555 !-------- Set up the temperature equations (ice and bedrock
1556 ! simultaneously) --------
1557 
1558 kr=0
1559 lgs_a1(kr) = 1.0_dp
1560 lgs_a2(kr) = -1.0_dp
1561 lgs_b(kr) = clb1
1562 
1563 #if (Q_LITHO==1)
1564 ! (coupled heat-conducting bedrock)
1565 
1566 do kr=1, krmax-1
1567  lgs_a0(kr) = - ctr1
1568  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
1569  lgs_a2(kr) = - ctr1
1570  lgs_b(kr) = temp_r(kr,j,i)
1571 end do
1572 
1573 #elif (Q_LITHO==0)
1574 ! (no coupled heat-conducting bedrock)
1575 
1576 do kr=1, krmax-1
1577  lgs_a0(kr) = 1.0_dp
1578  lgs_a1(kr) = 0.0_dp
1579  lgs_a2(kr) = -1.0_dp
1580  lgs_b(kr) = 2.0_dp*clb1
1581 end do
1582 
1583 #endif
1584 
1585 kr=krmax
1586 kc=0
1587 lgs_a0(kr) = ccb2
1588 lgs_a1(kr) = -(ccb1+ccb2)
1589 lgs_a2(kr) = ccb1
1590 lgs_b(kr) = ccb3+ccb4
1591 
1592 do kc=1, kcmax-1
1593 
1594 #if (ADV_VERT==1)
1595 
1596  lgs_a0(krmax+kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
1597  -ct5(kc)*ct6(kc-1)
1598  lgs_a1(krmax+kc) = 1.0_dp+ct5(kc)*(ct6(kc)+ct6(kc-1))
1599  lgs_a2(krmax+kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
1600  -ct5(kc)*ct6(kc)
1601 
1602 #elif (ADV_VERT==2)
1603 
1604  lgs_a0(krmax+kc) &
1605  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
1606  -ct5(kc)*ct6(kc-1)
1607  lgs_a1(krmax+kc) &
1608  = 1.0_dp &
1609  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
1610  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
1611  +ct5(kc)*(ct6(kc)+ct6(kc-1))
1612  lgs_a2(krmax+kc) &
1613  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
1614  -ct5(kc)*ct6(kc)
1615 
1616 #elif (ADV_VERT==3)
1617 
1618  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
1619 
1620  lgs_a0(krmax+kc) &
1621  = -max(adv_vert_help, 0.0_dp) &
1622  -ct5(kc)*ct6(kc-1)
1623  lgs_a1(krmax+kc) &
1624  = 1.0_dp &
1625  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
1626  +ct5(kc)*(ct6(kc)+ct6(kc-1))
1627  lgs_a2(krmax+kc) &
1628  = min(adv_vert_help, 0.0_dp) &
1629  -ct5(kc)*ct6(kc)
1630 
1631 #endif
1632 
1633 #if (ADV_HOR==2)
1634 
1635  lgs_b(krmax+kc) = temp_c(kc,j,i) + ct7(kc) &
1636  -dtt_2dxi* &
1637  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
1638  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
1639  *insq_g11_sgx(j,i) &
1640  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
1641  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
1642  *insq_g11_sgx(j,i-1) ) &
1643  -dtt_2deta* &
1644  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
1645  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
1646  *insq_g22_sgy(j,i) &
1647  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
1648  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
1649  *insq_g22_sgy(j-1,i) )
1650 
1651 #elif (ADV_HOR==3)
1652 
1653  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
1654  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
1655 
1656  lgs_b(krmax+kc) = temp_c(kc,j,i) + ct7(kc) &
1657  -dtt_dxi* &
1658  ( min(vx_c_help, 0.0_dp) &
1659  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
1660  *insq_g11_sgx(j,i) &
1661  +max(vx_c_help, 0.0_dp) &
1662  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
1663  *insq_g11_sgx(j,i-1) ) &
1664  -dtt_deta* &
1665  ( min(vy_c_help, 0.0_dp) &
1666  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
1667  *insq_g22_sgy(j,i) &
1668  +max(vy_c_help, 0.0_dp) &
1669  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
1670  *insq_g22_sgy(j-1,i) )
1671 
1672 #endif
1673 
1674 end do
1675 
1676 kc=kcmax
1677 lgs_a0(krmax+kc) = 0.0_dp
1678 lgs_a1(krmax+kc) = 1.0_dp
1679 lgs_b(krmax+kc) = temp_s(j,i)
1680 
1681 !-------- Solve system of linear equations --------
1682 
1683 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax+krmax)
1684 
1685 !-------- Assign the result --------
1686 
1687 do kr=0, krmax
1688  temp_r_neu(kr,j,i) = lgs_x(kr)
1689 end do
1690 
1691 do kc=0, kcmax
1692  temp_c_neu(kc,j,i) = lgs_x(krmax+kc)
1693 end do
1694 
1695 !-------- Set water content in the non-existing temperate layer
1696 ! to zero --------
1697 
1698 do kt=0, ktmax
1699  omega_t_neu(kt,j,i) = 0.0_dp
1700 end do
1701 
1702 !-------- Water drainage from the non-existing temperate layer --------
1703 
1704 q_tld(j,i) = 0.0_dp
1705 
1706 !-------- Set up the equations for the age of cold ice --------
1707 
1708 kc=0 ! adv_vert_sg(0) <= 0
1709 lgs_a1(kc) = 1.0_dp - min(adv_vert_sg(kc), 0.0_dp) ! (directed downward)
1710 lgs_a2(kc) = min(adv_vert_sg(kc), 0.0_dp) ! assumed/enforced
1711 
1712 #if (ADV_HOR==2)
1713 
1714 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1715  -dtt_2dxi* &
1716  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
1717  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1718  *insq_g11_sgx(j,i) &
1719  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
1720  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1721  *insq_g11_sgx(j,i-1) ) &
1722  -dtt_2deta* &
1723  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
1724  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1725  *insq_g22_sgy(j,i) &
1726  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
1727  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1728  *insq_g22_sgy(j-1,i) )
1729 
1730 #elif (ADV_HOR==3)
1731 
1732 vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
1733 vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
1734 
1735 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1736  -dtt_dxi* &
1737  ( min(vx_c_help, 0.0_dp) &
1738  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1739  *insq_g11_sgx(j,i) &
1740  +max(vx_c_help, 0.0_dp) &
1741  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1742  *insq_g11_sgx(j,i-1) ) &
1743  -dtt_deta* &
1744  ( min(vy_c_help, 0.0_dp) &
1745  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1746  *insq_g22_sgy(j,i) &
1747  +max(vy_c_help, 0.0_dp) &
1748  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1749  *insq_g22_sgy(j-1,i) )
1750 
1751 #endif
1752 
1753 do kc=1, kcmax-1
1754 
1755 #if (ADV_VERT==1)
1756 
1757  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
1758  -ci1(kc)*ci2(kc-1)
1759  lgs_a1(kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
1760  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
1761  -ci1(kc)*ci2(kc)
1762 
1763 #elif (ADV_VERT==2)
1764 
1765  lgs_a0(kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
1766  lgs_a1(kc) = 1.0_dp &
1767  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
1768  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
1769  lgs_a2(kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
1770 
1771 #elif (ADV_VERT==3)
1772 
1773  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
1774 
1775  lgs_a0(kc) = -max(adv_vert_help, 0.0_dp)
1776  lgs_a1(kc) = 1.0_dp &
1777  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp)
1778  lgs_a2(kc) = min(adv_vert_help, 0.0_dp)
1779 
1780 #endif
1781 
1782 #if (ADV_HOR==2)
1783 
1784  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1785  -dtt_2dxi* &
1786  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
1787  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1788  *insq_g11_sgx(j,i) &
1789  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
1790  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1791  *insq_g11_sgx(j,i-1) ) &
1792  -dtt_2deta* &
1793  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
1794  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1795  *insq_g22_sgy(j,i) &
1796  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
1797  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1798  *insq_g22_sgy(j-1,i) )
1799 
1800 #elif (ADV_HOR==3)
1801 
1802  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
1803  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
1804 
1805  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1806  -dtt_dxi* &
1807  ( min(vx_c_help, 0.0_dp) &
1808  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1809  *insq_g11_sgx(j,i) &
1810  +max(vx_c_help, 0.0_dp) &
1811  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1812  *insq_g11_sgx(j,i-1) ) &
1813  -dtt_deta* &
1814  ( min(vy_c_help, 0.0_dp) &
1815  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1816  *insq_g22_sgy(j,i) &
1817  +max(vy_c_help, 0.0_dp) &
1818  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1819  *insq_g22_sgy(j-1,i) )
1820 
1821 #endif
1822 
1823 end do
1824 
1825 kc=kcmax
1826 if (as_perp(j,i) >= zero) then
1827  lgs_a0(kc) = 0.0_dp
1828  lgs_a1(kc) = 1.0_dp
1829  lgs_b(kc) = 0.0_dp
1830 else
1831  lgs_a0(kc) = -max(adv_vert_sg(kc-1), 0.0_dp)
1832  lgs_a1(kc) = 1.0_dp + max(adv_vert_sg(kc-1), 0.0_dp)
1833  ! adv_vert_sg(KCMAX-1) >= 0 (directed upward)
1834  ! assumed/enforced
1835 #if (ADV_HOR==2)
1836 
1837  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1838  -dtt_2dxi* &
1839  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
1840  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1841  *insq_g11_sgx(j,i) &
1842  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
1843  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1844  *insq_g11_sgx(j,i-1) ) &
1845  -dtt_2deta* &
1846  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
1847  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1848  *insq_g22_sgy(j,i) &
1849  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
1850  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1851  *insq_g22_sgy(j-1,i) )
1852 
1853 #elif (ADV_HOR==3)
1854 
1855  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
1856  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
1857 
1858  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1859  -dtt_dxi* &
1860  ( min(vx_c_help, 0.0_dp) &
1861  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1862  *insq_g11_sgx(j,i) &
1863  +max(vx_c_help, 0.0_dp) &
1864  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1865  *insq_g11_sgx(j,i-1) ) &
1866  -dtt_deta* &
1867  ( min(vy_c_help, 0.0_dp) &
1868  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1869  *insq_g22_sgy(j,i) &
1870  +max(vy_c_help, 0.0_dp) &
1871  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1872  *insq_g22_sgy(j-1,i) )
1873 
1874 #endif
1875 
1876 end if
1877 
1878 !-------- Solve system of linear equations --------
1879 
1880 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
1881 
1882 !-------- Assign the result,
1883 ! restriction to interval [0, AGE_MAX yr] --------
1884 
1885 do kc=0, kcmax
1886 
1887  age_c_neu(kc,j,i) = lgs_x(kc)
1888 
1889  if (age_c_neu(kc,j,i) < (age_min*year_sec)) &
1890  age_c_neu(kc,j,i) = 0.0_dp
1891  if (age_c_neu(kc,j,i) > (age_max*year_sec)) &
1892  age_c_neu(kc,j,i) = age_max*year_sec
1893 
1894 end do
1895 
1896 !-------- Age of the ice in the non-existing temperate layer --------
1897 
1898 do kt=0, ktmax
1899  age_t_neu(kt,j,i) = age_c_neu(0,j,i)
1900 end do
1901 
1902 end subroutine calc_temp1
1903 
1904 !-------------------------------------------------------------------------------
1905 !> Computation of temperature and age for an ice column with a temperate base
1906 !! overlain by cold ice.
1907 !<------------------------------------------------------------------------------
1908 subroutine calc_temp2(at1, at2_1, at2_2, at3_1, at3_2, &
1909  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
1910  ai1, ai2, &
1911  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
1912 
1913 use ice_material_properties_m, only : ratefac_c, kappa_val, c_val, &
1914  creep, viscosity
1915 use sico_maths_m, only : tri_sle
1916 
1917 implicit none
1918 
1919 integer(i4b), intent(in) :: i, j
1920 real(dp), intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
1921  at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
1922  at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
1923  ai1(0:kcmax), ai2(0:kcmax), &
1924  atr1, alb1
1925 real(dp), intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
1926 
1927 integer(i4b) :: kc, kt, kr
1928 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
1929  ct5(0:kcmax), ct6(0:kcmax), ct7(0:kcmax), ctr1, clb1
1930 real(dp) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), ct3_sg(0:kcmax), &
1931  ct4_sg(0:kcmax), adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
1932 real(dp) :: ci1(0:kcmax), ci2(0:kcmax)
1933 real(dp) :: temp_c_help(0:kcmax)
1934 real(dp) :: vx_c_help, vy_c_help
1935 real(dp) :: adv_vert_help
1936 real(dp) :: dtt_dxi, dtt_deta
1937 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
1938  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
1939  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
1940  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
1941  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
1942 real(dp), parameter :: zero=0.0_dp
1943 
1944 !-------- Check for boundary points --------
1945 
1946 if ((i == 0).or.(i == imax).or.(j == 0).or.(j == jmax)) &
1947  stop ' >>> calc_temp2: Boundary points not allowed.'
1948 
1949 !-------- Abbreviations --------
1950 
1951 ctr1 = atr1
1952 clb1 = alb1*q_geo(j,i)
1953 
1954 #if (ADV_VERT==1)
1955 
1956 do kc=1, kcmax-1
1957  ct1(kc) = at1(kc)/h_c(j,i)*0.5_dp*(vz_c(kc,j,i)+vz_c(kc-1,j,i))
1958 end do
1959 
1960 kc=0
1961 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
1962  ! only needed for kc=0 ...
1963 kc=kcmax-1
1964 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
1965  ! ... and kc=KCMAX-1
1966 
1967 #elif (ADV_VERT==2 || ADV_VERT==3)
1968 
1969 do kc=0, kcmax-1
1970  ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
1971 end do
1972 
1973 #endif
1974 
1975 do kc=0, kcmax
1976 
1977  ct2(kc) = ( at2_1(kc)*dzm_dtau(j,i) &
1978  +at2_2(kc)*dh_c_dtau(j,i) )/h_c(j,i)
1979  ct3(kc) = ( at3_1(kc)*dzm_dxi_g(j,i) &
1980  +at3_2(kc)*dh_c_dxi_g(j,i) )/h_c(j,i) &
1981  *0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1)) *insq_g11_g(j,i)
1982  ct4(kc) = ( at4_1(kc)*dzm_deta_g(j,i) &
1983  +at4_2(kc)*dh_c_deta_g(j,i) )/h_c(j,i) &
1984  *0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i)) *insq_g22_g(j,i)
1985  ct5(kc) = at5(kc) &
1986  /c_val(temp_c(kc,j,i)) &
1987  /h_c(j,i)
1988 
1989 #if (DYNAMICS==2)
1990  if (.not.flag_shelfy_stream(j,i)) then
1991 #endif
1992  ct7(kc) = at7 &
1993  /c_val(temp_c(kc,j,i)) &
1994  *enh_c(kc,j,i) &
1995  *ratefac_c(temp_c(kc,j,i), temp_c_m(kc,j,i)) &
1996  *creep(sigma_c(kc,j,i)) &
1997  *sigma_c(kc,j,i)*sigma_c(kc,j,i)
1998 #if (DYNAMICS==2)
1999  else
2000  ct7(kc) = 2.0_dp*at7 &
2001  /c_val(temp_c(kc,j,i)) &
2002  *viscosity(de_c(kc,j,i), &
2003  temp_c(kc,j,i), temp_c_m(kc,j,i), 0.0_dp, &
2004  enh_c(kc,j,i), 0_i2b) &
2005  *de_c(kc,j,i)**2
2006  end if
2007 #endif
2008 
2009  ci1(kc) = ai1(kc)/h_c(j,i)
2010 
2011 end do
2012 
2013 #if (ADV_VERT==1)
2014 
2015 kc=0
2016 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2017 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2018 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2019 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2020 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! only needed for kc=0 ...
2021 kc=kcmax-1
2022 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2023 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2024 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2025 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2026 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! ... and kc=KCMAX-1
2027 
2028 #elif (ADV_VERT==2 || ADV_VERT==3)
2029 
2030 do kc=0, kcmax-1
2031  ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2032  ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2033  ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2034  adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2035  abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
2036 end do
2037 
2038 #endif
2039 
2040 do kc=0, kcmax-1
2041  temp_c_help(kc) = 0.5_dp*(temp_c(kc,j,i)+temp_c(kc+1,j,i))
2042  ct6(kc) = at6(kc) &
2043  *kappa_val(temp_c_help(kc)) &
2044  /h_c(j,i)
2045  ci2(kc) = ai2(kc)/h_c(j,i)
2046 end do
2047 
2048 #if (ADV_HOR==3)
2049 dtt_dxi = 2.0_dp*dtt_2dxi
2050 dtt_deta = 2.0_dp*dtt_2deta
2051 #endif
2052 
2053 !-------- Set up the equations for the bedrock temperature --------
2054 
2055 kr=0
2056 lgs_a1(kr) = 1.0_dp
2057 lgs_a2(kr) = -1.0_dp
2058 lgs_b(kr) = clb1
2059 
2060 #if (Q_LITHO==1)
2061 ! (coupled heat-conducting bedrock)
2062 
2063 do kr=1, krmax-1
2064  lgs_a0(kr) = - ctr1
2065  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
2066  lgs_a2(kr) = - ctr1
2067  lgs_b(kr) = temp_r(kr,j,i)
2068 end do
2069 
2070 #elif (Q_LITHO==0)
2071 ! (no coupled heat-conducting bedrock)
2072 
2073 do kr=1, krmax-1
2074  lgs_a0(kr) = 1.0_dp
2075  lgs_a1(kr) = 0.0_dp
2076  lgs_a2(kr) = -1.0_dp
2077  lgs_b(kr) = 2.0_dp*clb1
2078 end do
2079 
2080 #endif
2081 
2082 kr=krmax
2083 lgs_a0(kr) = 0.0_dp
2084 lgs_a1(kr) = 1.0_dp
2085 lgs_b(kr) = temp_t_m(0,j,i)
2086 
2087 !-------- Solve system of linear equations --------
2088 
2089 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
2090 
2091 !-------- Assign the result --------
2092 
2093 do kr=0, krmax
2094  temp_r_neu(kr,j,i) = lgs_x(kr)
2095 end do
2096 
2097 !-------- Set up the equations for the ice temperature --------
2098 
2099 kc=0
2100 lgs_a1(kc) = 1.0_dp
2101 lgs_a2(kc) = 0.0_dp
2102 lgs_b(kc) = temp_c_m(0,j,i)
2103 
2104 do kc=1, kcmax-1
2105 
2106 #if (ADV_VERT==1)
2107 
2108  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2109  -ct5(kc)*ct6(kc-1)
2110  lgs_a1(kc) = 1.0_dp+ct5(kc)*(ct6(kc)+ct6(kc-1))
2111  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2112  -ct5(kc)*ct6(kc)
2113 
2114 #elif (ADV_VERT==2)
2115 
2116  lgs_a0(kc) &
2117  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2118  -ct5(kc)*ct6(kc-1)
2119  lgs_a1(kc) &
2120  = 1.0_dp &
2121  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2122  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
2123  +ct5(kc)*(ct6(kc)+ct6(kc-1))
2124  lgs_a2(kc) &
2125  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
2126  -ct5(kc)*ct6(kc)
2127 
2128 #elif (ADV_VERT==3)
2129 
2130  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
2131 
2132  lgs_a0(kc) &
2133  = -max(adv_vert_help, 0.0_dp) &
2134  -ct5(kc)*ct6(kc-1)
2135  lgs_a1(kc) &
2136  = 1.0_dp &
2137  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
2138  +ct5(kc)*(ct6(kc)+ct6(kc-1))
2139  lgs_a2(kc) &
2140  = min(adv_vert_help, 0.0_dp) &
2141  -ct5(kc)*ct6(kc)
2142 
2143 #endif
2144 
2145 #if (ADV_HOR==2)
2146 
2147  lgs_b(kc) = temp_c(kc,j,i) + ct7(kc) &
2148  -dtt_2dxi* &
2149  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2150  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
2151  *insq_g11_sgx(j,i) &
2152  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2153  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
2154  *insq_g11_sgx(j,i-1) ) &
2155  -dtt_2deta* &
2156  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2157  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
2158  *insq_g22_sgy(j,i) &
2159  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2160  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
2161  *insq_g22_sgy(j-1,i) )
2162 
2163 #elif (ADV_HOR==3)
2164 
2165  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2166  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2167 
2168  lgs_b(kc) = temp_c(kc,j,i) + ct7(kc) &
2169  -dtt_dxi* &
2170  ( min(vx_c_help, 0.0_dp) &
2171  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
2172  *insq_g11_sgx(j,i) &
2173  +max(vx_c_help, 0.0_dp) &
2174  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
2175  *insq_g11_sgx(j,i-1) ) &
2176  -dtt_deta* &
2177  ( min(vy_c_help, 0.0_dp) &
2178  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
2179  *insq_g22_sgy(j,i) &
2180  +max(vy_c_help, 0.0_dp) &
2181  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
2182  *insq_g22_sgy(j-1,i) )
2183 
2184 #endif
2185 
2186 end do
2187 
2188 kc=kcmax
2189 lgs_a0(kc) = 0.0_dp
2190 lgs_a1(kc) = 1.0_dp
2191 lgs_b(kc) = temp_s(j,i)
2192 
2193 !-------- Solve system of linear equations --------
2194 
2195 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
2196 
2197 !-------- Assign the result --------
2198 
2199 do kc=0, kcmax
2200  temp_c_neu(kc,j,i) = lgs_x(kc)
2201 end do
2202 
2203 !-------- Set water content in the non-existing temperate layer
2204 ! to zero --------
2205 
2206 do kt=0, ktmax
2207  omega_t_neu(kt,j,i) = 0.0_dp
2208 end do
2209 
2210 !-------- Water drainage from the non-existing temperate layer --------
2211 
2212 q_tld(j,i) = 0.0_dp
2213 
2214 !-------- Set up the equations for the age of cold ice --------
2215 
2216 kc=0 ! adv_vert_sg(0) <= 0
2217 lgs_a1(kc) = 1.0_dp - min(adv_vert_sg(kc), 0.0_dp) ! (directed downward)
2218 lgs_a2(kc) = min(adv_vert_sg(kc), 0.0_dp) ! assumed/enforced
2219 
2220 #if (ADV_HOR==2)
2221 
2222 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2223  -dtt_2dxi* &
2224  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2225  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2226  *insq_g11_sgx(j,i) &
2227  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2228  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2229  *insq_g11_sgx(j,i-1) ) &
2230  -dtt_2deta* &
2231  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2232  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2233  *insq_g22_sgy(j,i) &
2234  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2235  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2236  *insq_g22_sgy(j-1,i) )
2237 
2238 #elif (ADV_HOR==3)
2239 
2240 vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2241 vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2242 
2243 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2244  -dtt_dxi* &
2245  ( min(vx_c_help, 0.0_dp) &
2246  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2247  *insq_g11_sgx(j,i) &
2248  +max(vx_c_help, 0.0_dp) &
2249  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2250  *insq_g11_sgx(j,i-1) ) &
2251  -dtt_deta* &
2252  ( min(vy_c_help, 0.0_dp) &
2253  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2254  *insq_g22_sgy(j,i) &
2255  +max(vy_c_help, 0.0_dp) &
2256  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2257  *insq_g22_sgy(j-1,i) )
2258 
2259 #endif
2260 
2261 do kc=1, kcmax-1
2262 
2263 #if (ADV_VERT==1)
2264 
2265  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2266  -ci1(kc)*ci2(kc-1)
2267  lgs_a1(kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
2268  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2269  -ci1(kc)*ci2(kc)
2270 
2271 #elif (ADV_VERT==2)
2272 
2273  lgs_a0(kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
2274  lgs_a1(kc) = 1.0_dp &
2275  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2276  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
2277  lgs_a2(kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
2278 
2279 #elif (ADV_VERT==3)
2280 
2281  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
2282 
2283  lgs_a0(kc) = -max(adv_vert_help, 0.0_dp)
2284  lgs_a1(kc) = 1.0_dp &
2285  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp)
2286  lgs_a2(kc) = min(adv_vert_help, 0.0_dp)
2287 
2288 #endif
2289 
2290 #if (ADV_HOR==2)
2291 
2292  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2293  -dtt_2dxi* &
2294  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2295  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2296  *insq_g11_sgx(j,i) &
2297  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2298  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2299  *insq_g11_sgx(j,i-1) ) &
2300  -dtt_2deta* &
2301  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2302  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2303  *insq_g22_sgy(j,i) &
2304  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2305  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2306  *insq_g22_sgy(j-1,i) )
2307 
2308 #elif (ADV_HOR==3)
2309 
2310  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2311  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2312 
2313  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2314  -dtt_dxi* &
2315  ( min(vx_c_help, 0.0_dp) &
2316  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2317  *insq_g11_sgx(j,i) &
2318  +max(vx_c_help, 0.0_dp) &
2319  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2320  *insq_g11_sgx(j,i-1) ) &
2321  -dtt_deta* &
2322  ( min(vy_c_help, 0.0_dp) &
2323  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2324  *insq_g22_sgy(j,i) &
2325  +max(vy_c_help, 0.0_dp) &
2326  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2327  *insq_g22_sgy(j-1,i) )
2328 
2329 #endif
2330 
2331 end do
2332 
2333 kc=kcmax
2334 if (as_perp(j,i) >= zero) then
2335  lgs_a0(kc) = 0.0_dp
2336  lgs_a1(kc) = 1.0_dp
2337  lgs_b(kc) = 0.0_dp
2338 else
2339  lgs_a0(kc) = -max(adv_vert_sg(kc-1), 0.0_dp)
2340  lgs_a1(kc) = 1.0_dp + max(adv_vert_sg(kc-1), 0.0_dp)
2341  ! adv_vert_sg(KCMAX-1) >= 0 (directed upward)
2342  ! assumed/enforced
2343 #if (ADV_HOR==2)
2344 
2345  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2346  -dtt_2dxi* &
2347  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2348  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2349  *insq_g11_sgx(j,i) &
2350  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2351  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2352  *insq_g11_sgx(j,i-1) ) &
2353  -dtt_2deta* &
2354  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2355  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2356  *insq_g22_sgy(j,i) &
2357  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2358  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2359  *insq_g22_sgy(j-1,i) )
2360 
2361 #elif (ADV_HOR==3)
2362 
2363  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2364  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2365 
2366  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2367  -dtt_dxi* &
2368  ( min(vx_c_help, 0.0_dp) &
2369  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2370  *insq_g11_sgx(j,i) &
2371  +max(vx_c_help, 0.0_dp) &
2372  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2373  *insq_g11_sgx(j,i-1) ) &
2374  -dtt_deta* &
2375  ( min(vy_c_help, 0.0_dp) &
2376  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2377  *insq_g22_sgy(j,i) &
2378  +max(vy_c_help, 0.0_dp) &
2379  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2380  *insq_g22_sgy(j-1,i) )
2381 
2382 #endif
2383 
2384 end if
2385 
2386 !-------- Solve system of linear equations --------
2387 
2388 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
2389 
2390 !-------- Assign the result,
2391 ! restriction to interval [0, AGE_MAX yr] --------
2392 
2393 do kc=0, kcmax
2394 
2395  age_c_neu(kc,j,i) = lgs_x(kc)
2396 
2397  if (age_c_neu(kc,j,i) < (age_min*year_sec)) &
2398  age_c_neu(kc,j,i) = 0.0_dp
2399  if (age_c_neu(kc,j,i) > (age_max*year_sec)) &
2400  age_c_neu(kc,j,i) = age_max*year_sec
2401 
2402 end do
2403 
2404 !-------- Age of the ice in the non-existing temperate layer --------
2405 
2406 do kt=0, ktmax
2407  age_t_neu(kt,j,i) = age_c_neu(0,j,i)
2408 end do
2409 
2410 end subroutine calc_temp2
2411 
2412 !-------------------------------------------------------------------------------
2413 !> Computation of temperature, water content and age for an ice column with a
2414 !! temperate base overlain by a temperate-ice layer.
2415 !<------------------------------------------------------------------------------
2416 subroutine calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
2417  at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
2418  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
2419  ai1, ai2, ai3, dzeta_t, &
2420  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
2421 
2422 use ice_material_properties_m, only : ratefac_c, ratefac_t, kappa_val, c_val, &
2423  creep, viscosity
2424 use sico_maths_m, only : tri_sle
2425 
2426 implicit none
2427 
2428 integer(i4b), intent(in) :: i, j
2429 real(dp), intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
2430  at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
2431  at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
2432  ai1(0:kcmax), ai2(0:kcmax), ai3, &
2433  atr1, am1, am2, alb1
2434 real(dp), intent(in) :: aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld
2435 real(dp), intent(in) :: dzeta_t, dtime_temp, dtt_2dxi, dtt_2deta
2436 
2437 integer(i4b) :: kc, kt, kr
2438 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
2439  ct5(0:kcmax), ct6(0:kcmax), ct7(0:kcmax), ctr1, cm1, cm2, &
2440  clb1
2441 real(dp) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), ct3_sg(0:kcmax), &
2442  ct4_sg(0:kcmax), adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
2443 real(dp) :: ci1(0:kcmax), ci2(0:kcmax), ci3
2444 real(dp) :: cw1(0:ktmax), cw2(0:ktmax), cw3(0:ktmax), cw4(0:ktmax), &
2445  cw5, cw7(0:ktmax), cw8, cw9(0:ktmax)
2446 real(dp) :: cw1_sg(0:ktmax), cw2_sg(0:ktmax), cw3_sg(0:ktmax), &
2447  cw4_sg(0:ktmax), adv_vert_w_sg(0:ktmax), abs_adv_vert_w_sg(0:ktmax)
2448 real(dp) :: sigma_c_help(0:kcmax), sigma_t_help(0:ktmax), &
2449  temp_c_help(0:kcmax)
2450 real(dp) :: vx_c_help, vy_c_help, vx_t_help, vy_t_help
2451 real(dp) :: adv_vert_help, adv_vert_w_help
2452 real(dp) :: dtt_dxi, dtt_deta
2453 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
2454  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
2455  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
2456  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
2457  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
2458 real(dp), parameter :: zero=0.0_dp
2459 
2460 !-------- Check for boundary points --------
2461 
2462 if ((i == 0).or.(i == imax).or.(j == 0).or.(j == jmax)) &
2463  stop ' >>> calc_temp3: Boundary points not allowed.'
2464 
2465 !-------- Abbreviations --------
2466 
2467 ctr1 = atr1
2468 cm1 = am1*h_c_neu(j,i)
2469 clb1 = alb1*q_geo(j,i)
2470 
2471 #if (ADV_VERT==1)
2472 
2473 do kc=1, kcmax-1
2474  ct1(kc) = at1(kc)/h_c_neu(j,i)*0.5_dp*(vz_c(kc,j,i)+vz_c(kc-1,j,i))
2475 end do
2476 
2477 kc=0
2478 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c_neu(j,i)*vz_c(kc,j,i)
2479  ! only needed for kc=0 ...
2480 kc=kcmax-1
2481 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c_neu(j,i)*vz_c(kc,j,i)
2482  ! ... and kc=KCMAX-1
2483 
2484 #elif (ADV_VERT==2 || ADV_VERT==3)
2485 
2486 do kc=0, kcmax-1
2487  ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c_neu(j,i)*vz_c(kc,j,i)
2488 end do
2489 
2490 #endif
2491 
2492 do kc=0, kcmax
2493 
2494  ct2(kc) = ( at2_1(kc)*dzm_dtau(j,i) &
2495  +at2_2(kc)*dh_c_dtau(j,i) )/h_c_neu(j,i)
2496  ct3(kc) = ( at3_1(kc)*dzm_dxi_g(j,i) &
2497  +at3_2(kc)*dh_c_dxi_g(j,i) )/h_c_neu(j,i) &
2498  *0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1)) *insq_g11_g(j,i)
2499  ct4(kc) = ( at4_1(kc)*dzm_deta_g(j,i) &
2500  +at4_2(kc)*dh_c_deta_g(j,i) )/h_c_neu(j,i) &
2501  *0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i)) *insq_g22_g(j,i)
2502  ct5(kc) = at5(kc) &
2503  /c_val(temp_c(kc,j,i)) &
2504  /h_c_neu(j,i)
2505 
2506  sigma_c_help(kc) &
2507  = rho*g*h_c_neu(j,i)*(1.0_dp-eaz_c_quotient(kc)) &
2508  *(dzs_dxi_g(j,i)**2+dzs_deta_g(j,i)**2)**0.5_dp
2509 
2510 #if (DYNAMICS==2)
2511  if (.not.flag_shelfy_stream(j,i)) then
2512 #endif
2513  ct7(kc) = at7 &
2514  /c_val(temp_c(kc,j,i)) &
2515  *enh_c(kc,j,i) &
2516  *ratefac_c(temp_c(kc,j,i), temp_c_m(kc,j,i)) &
2517  *creep(sigma_c_help(kc)) &
2518  *sigma_c_help(kc)*sigma_c_help(kc)
2519 #if (DYNAMICS==2)
2520  else
2521  ct7(kc) = 2.0_dp*at7 &
2522  /c_val(temp_c(kc,j,i)) &
2523  *viscosity(de_c(kc,j,i), &
2524  temp_c(kc,j,i), temp_c_m(kc,j,i), 0.0_dp, &
2525  enh_c(kc,j,i), 0_i2b) &
2526  *de_c(kc,j,i)**2
2527  end if
2528 #endif
2529 
2530  ci1(kc) = ai1(kc)/h_c_neu(j,i)
2531 
2532 end do
2533 
2534 #if (ADV_VERT==1)
2535 
2536 kc=0
2537 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2538 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2539 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2540 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2541 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! only needed for kc=0 ...
2542 kc=kcmax-1
2543 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2544 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2545 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2546 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2547 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! ... and kc=KCMAX-1
2548 
2549 #elif (ADV_VERT==2 || ADV_VERT==3)
2550 
2551 do kc=0, kcmax-1
2552  ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2553  ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2554  ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2555  adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2556  abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
2557 end do
2558 
2559 #endif
2560 
2561 do kc=0, kcmax-1
2562  temp_c_help(kc) = 0.5_dp*(temp_c(kc,j,i)+temp_c(kc+1,j,i))
2563  ct6(kc) = at6(kc) &
2564  *kappa_val(temp_c_help(kc)) &
2565  /h_c_neu(j,i)
2566  ci2(kc) = ai2(kc)/h_c_neu(j,i)
2567 end do
2568 
2569 cw5 = aw5/(h_t_neu(j,i)**2)
2570 cw8 = aw8
2571 ci3 = ai3/(h_t_neu(j,i)**2)
2572 
2573 #if (ADV_VERT==1)
2574 
2575 do kt=1, ktmax-1
2576  cw1(kt) = aw1/h_t_neu(j,i)*0.5_dp*(vz_t(kt,j,i)+vz_t(kt-1,j,i))
2577 end do
2578 
2579 kt=ktmax
2580 cw1(kt) = aw1/h_t_neu(j,i)*0.5_dp*(vz_t(kt-1,j,i)+vz_c(0,j,i))
2581 
2582 kt=0
2583 cw1_sg(kt) = aw1/h_t_neu(j,i)*vz_t(kt,j,i)
2584  ! only needed for kt=0 ...
2585 kt=ktmax-1
2586 cw1_sg(kt) = aw1/h_t_neu(j,i)*vz_t(kt,j,i)
2587  ! ... and kt=KTMAX-1
2588 
2589 #elif (ADV_VERT==2 || ADV_VERT==3)
2590 
2591 do kt=0, ktmax-1
2592  cw1_sg(kt) = aw1/h_t_neu(j,i)*vz_t(kt,j,i)
2593 end do
2594 
2595 #endif
2596 
2597 do kt=1, ktmax-1
2598  cw9(kt) = aw9 &
2599  *c_val(temp_t_m(kt,j,i)) &
2600  *( dzs_dtau(j,i) &
2601  +0.5_dp*(vx_t(kt,j,i)+vx_t(kt,j,i-1))*dzs_dxi_g(j,i) &
2602  +0.5_dp*(vy_t(kt,j,i)+vy_t(kt,j-1,i))*dzs_deta_g(j,i) &
2603  -0.5_dp*(vz_t(kt,j,i)+vz_t(kt-1,j,i)) )
2604 end do
2605 
2606 do kt=0, ktmax
2607 
2608  cw2(kt) = aw2*(dzb_dtau(j,i)+zeta_t(kt)*dh_t_dtau(j,i)) &
2609  /h_t_neu(j,i)
2610  cw3(kt) = aw3*(dzb_dxi_g(j,i)+zeta_t(kt)*dh_t_dxi_g(j,i)) &
2611  /h_t_neu(j,i) &
2612  *0.5_dp*(vx_t(kt,j,i)+vx_t(kt,j,i-1)) *insq_g11_g(j,i)
2613  cw4(kt) = aw4*(dzb_deta_g(j,i)+zeta_t(kt)*dh_t_deta_g(j,i)) &
2614  /h_t_neu(j,i) &
2615  *0.5_dp*(vy_t(kt,j,i)+vy_t(kt,j-1,i)) *insq_g22_g(j,i)
2616  sigma_t_help(kt) &
2617  = sigma_c_help(0) &
2618  + rho*g*h_t_neu(j,i)*(1.0_dp-zeta_t(kt)) &
2619  *(dzs_dxi_g(j,i)**2+dzs_deta_g(j,i)**2)**0.5_dp
2620 
2621 #if (DYNAMICS==2)
2622  if (.not.flag_shelfy_stream(j,i)) then
2623 #endif
2624  cw7(kt) = aw7 &
2625  *enh_t(kt,j,i) &
2626  *ratefac_t(omega_t(kt,j,i)) &
2627  *creep(sigma_t_help(kt)) &
2628  *sigma_t_help(kt)*sigma_t_help(kt)
2629 #if (DYNAMICS==2)
2630  else
2631  cw7(kt) = 2.0_dp*aw7 &
2632  *viscosity(de_t(kt,j,i), &
2633  temp_t_m(kt,j,i), temp_t_m(kt,j,i), &
2634  omega_t(kt,j,i), &
2635  enh_t(kt,j,i), 1_i2b) &
2636  *de_t(kt,j,i)**2
2637  end if
2638 #endif
2639 
2640 end do
2641 
2642 #if (ADV_VERT==1)
2643 
2644 kt=0
2645 cw2_sg(kt) = 0.5_dp*(cw2(kt)+cw2(kt+1))
2646 cw3_sg(kt) = 0.5_dp*(cw3(kt)+cw3(kt+1))
2647 cw4_sg(kt) = 0.5_dp*(cw4(kt)+cw4(kt+1))
2648 adv_vert_w_sg(kt) = cw1_sg(kt)-cw2_sg(kt)-cw3_sg(kt)-cw4_sg(kt)
2649 abs_adv_vert_w_sg(kt) = abs(adv_vert_w_sg(kt)) ! only needed for kt=0 ...
2650 kt=ktmax-1
2651 cw2_sg(kt) = 0.5_dp*(cw2(kt)+cw2(kt+1))
2652 cw3_sg(kt) = 0.5_dp*(cw3(kt)+cw3(kt+1))
2653 cw4_sg(kt) = 0.5_dp*(cw4(kt)+cw4(kt+1))
2654 adv_vert_w_sg(kt) = cw1_sg(kt)-cw2_sg(kt)-cw3_sg(kt)-cw4_sg(kt)
2655 abs_adv_vert_w_sg(kt) = abs(adv_vert_w_sg(kt)) ! ... and kt=KTMAX-1
2656 
2657 #elif (ADV_VERT==2 || ADV_VERT==3)
2658 
2659 do kt=0, ktmax-1
2660  cw2_sg(kt) = 0.5_dp*(cw2(kt)+cw2(kt+1))
2661  cw3_sg(kt) = 0.5_dp*(cw3(kt)+cw3(kt+1))
2662  cw4_sg(kt) = 0.5_dp*(cw4(kt)+cw4(kt+1))
2663  adv_vert_w_sg(kt) = cw1_sg(kt)-cw2_sg(kt)-cw3_sg(kt)-cw4_sg(kt)
2664  abs_adv_vert_w_sg(kt) = abs(adv_vert_w_sg(kt))
2665 end do
2666 
2667 #endif
2668 
2669 #if (ADV_HOR==3)
2670 dtt_dxi = 2.0_dp*dtt_2dxi
2671 dtt_deta = 2.0_dp*dtt_2deta
2672 #endif
2673 
2674 !-------- Set up the equations for the bedrock temperature --------
2675 
2676 kr=0
2677 lgs_a1(kr) = 1.0_dp
2678 lgs_a2(kr) = -1.0_dp
2679 lgs_b(kr) = clb1
2680 
2681 #if (Q_LITHO==1)
2682 ! (coupled heat-conducting bedrock)
2683 
2684 do kr=1, krmax-1
2685  lgs_a0(kr) = - ctr1
2686  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
2687  lgs_a2(kr) = - ctr1
2688  lgs_b(kr) = temp_r(kr,j,i)
2689 end do
2690 
2691 #elif (Q_LITHO==0)
2692 ! (no coupled heat-conducting bedrock)
2693 
2694 do kr=1, krmax-1
2695  lgs_a0(kr) = 1.0_dp
2696  lgs_a1(kr) = 0.0_dp
2697  lgs_a2(kr) = -1.0_dp
2698  lgs_b(kr) = 2.0_dp*clb1
2699 end do
2700 
2701 #endif
2702 
2703 kr=krmax
2704 lgs_a0(kr) = 0.0_dp
2705 lgs_a1(kr) = 1.0_dp
2706 lgs_b(kr) = temp_t_m(0,j,i)
2707 
2708 !-------- Solve system of linear equations --------
2709 
2710 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
2711 
2712 !-------- Assign the result --------
2713 
2714 do kr=0, krmax
2715  temp_r_neu(kr,j,i) = lgs_x(kr)
2716 end do
2717 
2718 !-------- Set up the equations for the water content in
2719 ! temperate ice --------
2720 
2721 kt=0
2722 lgs_a1(kt) = 1.0_dp
2723 lgs_a2(kt) = -1.0_dp
2724 lgs_b(kt) = 0.0_dp
2725 
2726 do kt=1, ktmax-1
2727 
2728 #if (ADV_VERT==1)
2729 
2730  lgs_a0(kt) = -0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - cw5
2731  lgs_a1(kt) = 1.0_dp + 2.0_dp*cw5
2732  lgs_a2(kt) = 0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - cw5
2733 
2734 #elif (ADV_VERT==2)
2735 
2736  lgs_a0(kt) &
2737  = -0.5_dp*(adv_vert_w_sg(kt-1)+abs_adv_vert_w_sg(kt-1)) &
2738  -cw5
2739  lgs_a1(kt) &
2740  = 1.0_dp &
2741  +0.5_dp*(adv_vert_w_sg(kt-1)+abs_adv_vert_w_sg(kt-1)) &
2742  -0.5_dp*(adv_vert_w_sg(kt) -abs_adv_vert_w_sg(kt) ) &
2743  +2.0_dp*cw5
2744  lgs_a2(kt) &
2745  = 0.5_dp*(adv_vert_w_sg(kt) -abs_adv_vert_w_sg(kt) ) &
2746  -cw5
2747 
2748 #elif (ADV_VERT==3)
2749 
2750  adv_vert_w_help = 0.5_dp*(adv_vert_w_sg(kt)+adv_vert_w_sg(kt-1))
2751 
2752  lgs_a0(kt) &
2753  = -max(adv_vert_w_help, 0.0_dp) &
2754  -cw5
2755  lgs_a1(kt) &
2756  = 1.0_dp &
2757  +max(adv_vert_w_help, 0.0_dp)-min(adv_vert_w_help, 0.0_dp) &
2758  +2.0_dp*cw5
2759  lgs_a2(kt) &
2760  = min(adv_vert_w_help, 0.0_dp) &
2761  -cw5
2762 
2763 #endif
2764 
2765 #if (ADV_HOR==2)
2766 
2767  lgs_b(kt) = omega_t(kt,j,i) + cw7(kt) + cw8 + cw9(kt) &
2768  -dtt_2dxi* &
2769  ( (vx_t(kt,j,i)-abs(vx_t(kt,j,i))) &
2770  *(omega_t(kt,j,i+1)-omega_t(kt,j,i)) &
2771  *insq_g11_sgx(j,i) &
2772  +(vx_t(kt,j,i-1)+abs(vx_t(kt,j,i-1))) &
2773  *(omega_t(kt,j,i)-omega_t(kt,j,i-1)) &
2774  *insq_g11_sgx(j,i-1) ) &
2775  -dtt_2deta* &
2776  ( (vy_t(kt,j,i)-abs(vy_t(kt,j,i))) &
2777  *(omega_t(kt,j+1,i)-omega_t(kt,j,i)) &
2778  *insq_g22_sgy(j,i) &
2779  +(vy_t(kt,j-1,i)+abs(vy_t(kt,j-1,i))) &
2780  *(omega_t(kt,j,i)-omega_t(kt,j-1,i)) &
2781  *insq_g22_sgy(j-1,i) )
2782 
2783 #elif (ADV_HOR==3)
2784 
2785  vx_t_help = 0.5_dp*(vx_t(kt,j,i)+vx_t(kt,j,i-1))
2786  vy_t_help = 0.5_dp*(vy_t(kt,j,i)+vy_t(kt,j-1,i))
2787 
2788  lgs_b(kt) = omega_t(kt,j,i) + cw7(kt) + cw8 + cw9(kt) &
2789  -dtt_dxi* &
2790  ( min(vx_t_help, 0.0_dp) &
2791  *(omega_t(kt,j,i+1)-omega_t(kt,j,i)) &
2792  *insq_g11_sgx(j,i) &
2793  +max(vx_t_help, 0.0_dp) &
2794  *(omega_t(kt,j,i)-omega_t(kt,j,i-1)) &
2795  *insq_g11_sgx(j,i-1) ) &
2796  -dtt_deta* &
2797  ( min(vy_t_help, 0.0_dp) &
2798  *(omega_t(kt,j+1,i)-omega_t(kt,j,i)) &
2799  *insq_g22_sgy(j,i) &
2800  +max(vy_t_help, 0.0_dp) &
2801  *(omega_t(kt,j,i)-omega_t(kt,j-1,i)) &
2802  *insq_g22_sgy(j-1,i) )
2803 
2804 #endif
2805 
2806 end do
2807 
2808 kt=ktmax
2809 
2810 #if (!defined(CTS_MELTING_FREEZING) || CTS_MELTING_FREEZING==1)
2811 
2812 if (am_perp(j,i) >= zero) then ! melting condition
2813  lgs_a0(kt) = 0.0_dp
2814  lgs_a1(kt) = 1.0_dp
2815  lgs_b(kt) = 0.0_dp
2816 else ! am_perp(j,i) < 0.0, freezing condition
2817  lgs_a0(kt) = -1.0_dp
2818  lgs_a1(kt) = 1.0_dp
2819  lgs_b(kt) = 0.0_dp
2820 end if
2821 
2822 #elif (CTS_MELTING_FREEZING==2)
2823 
2824 lgs_a0(kt) = 0.0_dp
2825 lgs_a1(kt) = 1.0_dp ! melting condition assumed
2826 lgs_b(kt) = 0.0_dp
2827 
2828 #else
2829 stop ' >>> calc_temp3: CTS_MELTING_FREEZING must be either 1 or 2!'
2830 #endif
2831 
2832 !-------- Solve system of linear equations --------
2833 
2834 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, ktmax)
2835 
2836 !-------- Assign the result, compute the water drainage --------
2837 
2838 q_tld(j,i) = 0.0_dp
2839 
2840 do kt=0, ktmax
2841 
2842  if (lgs_x(kt) < zero) then
2843  omega_t_neu(kt,j,i) = 0.0_dp ! (as a precaution)
2844  else if (lgs_x(kt) < omega_max) then
2845  omega_t_neu(kt,j,i) = lgs_x(kt)
2846  else
2847  omega_t_neu(kt,j,i) = omega_max
2848  q_tld(j,i) = q_tld(j,i) &
2849  +aqtld*h_t_neu(j,i)*(lgs_x(kt)-omega_max)
2850  end if
2851 
2852 end do
2853 
2854 !-------- Set up the equations for the ice temperature --------
2855 
2856 ! ------ Abbreviation for the jump of the temperature gradient with
2857 ! the new omega
2858 
2859 #if (!defined(CTS_MELTING_FREEZING) || CTS_MELTING_FREEZING==1)
2860 
2861 if (am_perp(j,i) >= zero) then ! melting condition
2862  cm2 = 0.0_dp
2863 else ! am_perp(j,i) < 0.0, freezing condition
2864  cm2 = am2*h_c_neu(j,i)*omega_t_neu(ktmax,j,i)*am_perp(j,i) &
2865  /kappa_val(temp_c(0,j,i))
2866 end if
2867 
2868 #elif (CTS_MELTING_FREEZING==2)
2869 
2870 cm2 = 0.0_dp ! melting condition assumed
2871 
2872 #else
2873 stop ' >>> calc_temp3: CTS_MELTING_FREEZING must be either 1 or 2!'
2874 #endif
2875 
2876 kc=0
2877 lgs_a1(kc) = 1.0_dp
2878 lgs_a2(kc) = -1.0_dp
2879 lgs_b(kc) = -cm1-cm2
2880 
2881 do kc=1, kcmax-1
2882 
2883 #if (ADV_VERT==1)
2884 
2885  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2886  -ct5(kc)*ct6(kc-1)
2887  lgs_a1(kc) = 1.0_dp+ct5(kc)*(ct6(kc)+ct6(kc-1))
2888  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2889  -ct5(kc)*ct6(kc)
2890 
2891 #elif (ADV_VERT==2)
2892 
2893  lgs_a0(kc) &
2894  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2895  -ct5(kc)*ct6(kc-1)
2896  lgs_a1(kc) &
2897  = 1.0_dp &
2898  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2899  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
2900  +ct5(kc)*(ct6(kc)+ct6(kc-1))
2901  lgs_a2(kc) &
2902  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
2903  -ct5(kc)*ct6(kc)
2904 
2905 #elif (ADV_VERT==3)
2906 
2907  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
2908 
2909  lgs_a0(kc) &
2910  = -max(adv_vert_help, 0.0_dp) &
2911  -ct5(kc)*ct6(kc-1)
2912  lgs_a1(kc) &
2913  = 1.0_dp &
2914  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
2915  +ct5(kc)*(ct6(kc)+ct6(kc-1))
2916  lgs_a2(kc) &
2917  = min(adv_vert_help, 0.0_dp) &
2918  -ct5(kc)*ct6(kc)
2919 
2920 #endif
2921 
2922 #if (ADV_HOR==2)
2923 
2924  lgs_b(kc) = temp_c(kc,j,i) + ct7(kc) &
2925  -dtt_2dxi* &
2926  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2927  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
2928  *insq_g11_sgx(j,i) &
2929  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2930  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
2931  *insq_g11_sgx(j,i-1) ) &
2932  -dtt_2deta* &
2933  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2934  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
2935  *insq_g22_sgy(j,i) &
2936  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2937  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
2938  *insq_g22_sgy(j-1,i) )
2939 
2940 #elif (ADV_HOR==3)
2941 
2942  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2943  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2944 
2945  lgs_b(kc) = temp_c(kc,j,i) + ct7(kc) &
2946  -dtt_dxi* &
2947  ( min(vx_c_help, 0.0_dp) &
2948  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
2949  *insq_g11_sgx(j,i) &
2950  +max(vx_c_help, 0.0_dp) &
2951  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
2952  *insq_g11_sgx(j,i-1) ) &
2953  -dtt_deta* &
2954  ( min(vy_c_help, 0.0_dp) &
2955  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
2956  *insq_g22_sgy(j,i) &
2957  +max(vy_c_help, 0.0_dp) &
2958  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
2959  *insq_g22_sgy(j-1,i) )
2960 
2961 #endif
2962 
2963 end do
2964 
2965 kc=kcmax
2966 lgs_a0(kc) = 0.0_dp
2967 lgs_a1(kc) = 1.0_dp
2968 lgs_b(kc) = temp_s(j,i)
2969 
2970 !-------- Solve system of linear equations --------
2971 
2972 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
2973 
2974 !-------- Assign the result --------
2975 
2976 do kc=0, kcmax
2977  temp_c_neu(kc,j,i) = lgs_x(kc)
2978 end do
2979 
2980 !-------- Set up the equations for the age (cold and temperate ice
2981 ! simultaneously) --------
2982 
2983 kt=0 ! adv_vert_w_sg(0) <= 0
2984 lgs_a1(kt) = 1.0_dp - min(adv_vert_w_sg(kt), 0.0_dp) ! (directed downward)
2985 lgs_a2(kt) = min(adv_vert_w_sg(kt), 0.0_dp) ! assumed/enforced
2986 
2987 #if (ADV_HOR==2)
2988 
2989 lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
2990  -dtt_2dxi* &
2991  ( (vx_t(kt,j,i)-abs(vx_t(kt,j,i))) &
2992  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
2993  *insq_g11_sgx(j,i) &
2994  +(vx_t(kt,j,i-1)+abs(vx_t(kt,j,i-1))) &
2995  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
2996  *insq_g11_sgx(j,i-1) ) &
2997  -dtt_2deta* &
2998  ( (vy_t(kt,j,i)-abs(vy_t(kt,j,i))) &
2999  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
3000  *insq_g22_sgy(j,i) &
3001  +(vy_t(kt,j-1,i)+abs(vy_t(kt,j-1,i))) &
3002  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
3003  *insq_g22_sgy(j-1,i) )
3004 
3005 #elif (ADV_HOR==3)
3006 
3007 vx_t_help = 0.5_dp*(vx_t(kt,j,i)+vx_t(kt,j,i-1))
3008 vy_t_help = 0.5_dp*(vy_t(kt,j,i)+vy_t(kt,j-1,i))
3009 
3010 lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
3011  -dtt_dxi* &
3012  ( min(vx_t_help, 0.0_dp) &
3013  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
3014  *insq_g11_sgx(j,i) &
3015  +max(vx_t_help, 0.0_dp) &
3016  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
3017  *insq_g11_sgx(j,i-1) ) &
3018  -dtt_deta* &
3019  ( min(vy_t_help, 0.0_dp) &
3020  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
3021  *insq_g22_sgy(j,i) &
3022  +max(vy_t_help, 0.0_dp) &
3023  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
3024  *insq_g22_sgy(j-1,i) )
3025 
3026 #endif
3027 
3028 do kt=1, ktmax-1
3029 
3030 #if (ADV_VERT==1)
3031 
3032  lgs_a0(kt) = -0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - ci3
3033  lgs_a1(kt) = 1.0_dp + 2.0_dp*ci3
3034  lgs_a2(kt) = 0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - ci3
3035 
3036 #elif (ADV_VERT==2)
3037 
3038  lgs_a0(kt) = -0.5_dp*(adv_vert_w_sg(kt-1)+abs_adv_vert_w_sg(kt-1))
3039  lgs_a1(kt) = 1.0_dp &
3040  +0.5_dp*(adv_vert_w_sg(kt-1)+abs_adv_vert_w_sg(kt-1)) &
3041  -0.5_dp*(adv_vert_w_sg(kt) -abs_adv_vert_w_sg(kt) )
3042  lgs_a2(kt) = 0.5_dp*(adv_vert_w_sg(kt) -abs_adv_vert_w_sg(kt) )
3043 
3044 #elif (ADV_VERT==3)
3045 
3046  adv_vert_w_help = 0.5_dp*(adv_vert_w_sg(kt)+adv_vert_w_sg(kt-1))
3047 
3048  lgs_a0(kt) = -max(adv_vert_w_help, 0.0_dp)
3049  lgs_a1(kt) = 1.0_dp &
3050  +max(adv_vert_w_help, 0.0_dp)-min(adv_vert_w_help, 0.0_dp)
3051  lgs_a2(kt) = min(adv_vert_w_help, 0.0_dp)
3052 
3053 #endif
3054 
3055 #if (ADV_HOR==2)
3056 
3057  lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
3058  -dtt_2dxi* &
3059  ( (vx_t(kt,j,i)-abs(vx_t(kt,j,i))) &
3060  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
3061  *insq_g11_sgx(j,i) &
3062  +(vx_t(kt,j,i-1)+abs(vx_t(kt,j,i-1))) &
3063  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
3064  *insq_g11_sgx(j,i-1) ) &
3065  -dtt_2deta* &
3066  ( (vy_t(kt,j,i)-abs(vy_t(kt,j,i))) &
3067  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
3068  *insq_g22_sgy(j,i) &
3069  +(vy_t(kt,j-1,i)+abs(vy_t(kt,j-1,i))) &
3070  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
3071  *insq_g22_sgy(j-1,i) )
3072 
3073 #elif (ADV_HOR==3)
3074 
3075  vx_t_help = 0.5_dp*(vx_t(kt,j,i)+vx_t(kt,j,i-1))
3076  vy_t_help = 0.5_dp*(vy_t(kt,j,i)+vy_t(kt,j-1,i))
3077 
3078  lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
3079  -dtt_dxi* &
3080  ( min(vx_t_help, 0.0_dp) &
3081  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
3082  *insq_g11_sgx(j,i) &
3083  +max(vx_t_help, 0.0_dp) &
3084  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
3085  *insq_g11_sgx(j,i-1) ) &
3086  -dtt_deta* &
3087  ( min(vy_t_help, 0.0_dp) &
3088  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
3089  *insq_g22_sgy(j,i) &
3090  +max(vy_t_help, 0.0_dp) &
3091  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
3092  *insq_g22_sgy(j-1,i) )
3093 
3094 #endif
3095 
3096 end do
3097 
3098 #if (ADV_VERT==1)
3099 
3100 kt=ktmax
3101 kc=0
3102 
3103 lgs_a0(kt) = -0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - ci3
3104 lgs_a1(kt) = 1.0_dp + 2.0_dp*ci3
3105 lgs_a2(kt) = 0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - ci3
3106 
3107 #if (ADV_HOR==2)
3108 
3109 lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
3110  -dtt_2dxi* &
3111  ( (vx_t(kt,j,i)-abs(vx_t(kt,j,i))) &
3112  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
3113  *insq_g11_sgx(j,i) &
3114  +(vx_t(kt,j,i-1)+abs(vx_t(kt,j,i-1))) &
3115  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
3116  *insq_g11_sgx(j,i-1) ) &
3117  -dtt_2deta* &
3118  ( (vy_t(kt,j,i)-abs(vy_t(kt,j,i))) &
3119  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
3120  *insq_g22_sgy(j,i) &
3121  +(vy_t(kt,j-1,i)+abs(vy_t(kt,j-1,i))) &
3122  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
3123  *insq_g22_sgy(j-1,i) )
3124 
3125 #elif (ADV_HOR==3)
3126 
3127 vx_t_help = 0.5_dp*(vx_t(kt,j,i)+vx_t(kt,j,i-1))
3128 vy_t_help = 0.5_dp*(vy_t(kt,j,i)+vy_t(kt,j-1,i))
3129 
3130 lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
3131  -dtt_dxi* &
3132  ( min(vx_t_help, 0.0_dp) &
3133  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
3134  *insq_g11_sgx(j,i) &
3135  +max(vx_t_help, 0.0_dp) &
3136  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
3137  *insq_g11_sgx(j,i-1) ) &
3138  -dtt_deta* &
3139  ( min(vy_t_help, 0.0_dp) &
3140  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
3141  *insq_g22_sgy(j,i) &
3142  +max(vy_t_help, 0.0_dp) &
3143  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
3144  *insq_g22_sgy(j-1,i) )
3145 
3146 #endif
3147 
3148 #elif (ADV_VERT==2 || ADV_VERT==3)
3149 
3150 kt=ktmax
3151 kc=0
3152 
3153 if (adv_vert_sg(kc) <= zero) then
3154 
3155  lgs_a0(ktmax+kc) = 0.0_dp
3156  lgs_a1(ktmax+kc) = 1.0_dp - adv_vert_sg(kc)
3157  lgs_a2(ktmax+kc) = adv_vert_sg(kc)
3158 
3159 #if (ADV_HOR==2)
3160 
3161  lgs_b(ktmax+kc) = age_c(kc,j,i) + dtime_temp &
3162  -dtt_2dxi* &
3163  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
3164  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
3165  *insq_g11_sgx(j,i) &
3166  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
3167  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
3168  *insq_g11_sgx(j,i-1) ) &
3169  -dtt_2deta* &
3170  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
3171  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
3172  *insq_g22_sgy(j,i) &
3173  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
3174  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
3175  *insq_g22_sgy(j-1,i) )
3176 
3177 #elif (ADV_HOR==3)
3178 
3179  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
3180  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
3181 
3182  lgs_b(ktmax+kc) = age_c(kc,j,i) + dtime_temp &
3183  -dtt_dxi* &
3184  ( min(vx_c_help, 0.0_dp) &
3185  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
3186  *insq_g11_sgx(j,i) &
3187  +max(vx_c_help, 0.0_dp) &
3188  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
3189  *insq_g11_sgx(j,i-1) ) &
3190  -dtt_deta* &
3191  ( min(vy_c_help, 0.0_dp) &
3192  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
3193  *insq_g22_sgy(j,i) &
3194  +max(vy_c_help, 0.0_dp) &
3195  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
3196  *insq_g22_sgy(j-1,i) )
3197 
3198 #endif
3199 
3200 else if (adv_vert_w_sg(kt-1) >= zero) then
3201 
3202  lgs_a0(kt) = -adv_vert_w_sg(kt-1)
3203  lgs_a1(kt) = 1.0_dp + adv_vert_w_sg(kt-1)
3204  lgs_a2(kt) = 0.0_dp
3205 
3206 #if (ADV_HOR==2)
3207 
3208  lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
3209  -dtt_2dxi* &
3210  ( (vx_t(kt,j,i)-abs(vx_t(kt,j,i))) &
3211  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
3212  *insq_g11_sgx(j,i) &
3213  +(vx_t(kt,j,i-1)+abs(vx_t(kt,j,i-1))) &
3214  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
3215  *insq_g11_sgx(j,i-1) ) &
3216  -dtt_2deta* &
3217  ( (vy_t(kt,j,i)-abs(vy_t(kt,j,i))) &
3218  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
3219  *insq_g22_sgy(j,i) &
3220  +(vy_t(kt,j-1,i)+abs(vy_t(kt,j-1,i))) &
3221  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
3222  *insq_g22_sgy(j-1,i) )
3223 
3224 #elif (ADV_HOR==3)
3225 
3226  vx_t_help = 0.5_dp*(vx_t(kt,j,i)+vx_t(kt,j,i-1))
3227  vy_t_help = 0.5_dp*(vy_t(kt,j,i)+vy_t(kt,j-1,i))
3228 
3229  lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
3230  -dtt_dxi* &
3231  ( min(vx_t_help, 0.0_dp) &
3232  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
3233  *insq_g11_sgx(j,i) &
3234  +max(vx_t_help, 0.0_dp) &
3235  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
3236  *insq_g11_sgx(j,i-1) ) &
3237  -dtt_deta* &
3238  ( min(vy_t_help, 0.0_dp) &
3239  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
3240  *insq_g22_sgy(j,i) &
3241  +max(vy_t_help, 0.0_dp) &
3242  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
3243  *insq_g22_sgy(j-1,i) )
3244 
3245 #endif
3246 
3247 else
3248 
3249  lgs_a0(kt) = -0.5_dp
3250  lgs_a1(kt) = 1.0_dp
3251  lgs_a2(kt) = -0.5_dp
3252  lgs_b(kt) = 0.0_dp
3253  ! Makeshift: Average of age_c(kc=1) and age_t(kt=KTMAX-1)
3254 
3255 end if
3256 
3257 #endif
3258 
3259 do kc=1, kcmax-1
3260 
3261 #if (ADV_VERT==1)
3262 
3263  lgs_a0(ktmax+kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
3264  -ci1(kc)*ci2(kc-1)
3265  lgs_a1(ktmax+kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
3266  lgs_a2(ktmax+kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
3267  -ci1(kc)*ci2(kc)
3268 
3269 #elif (ADV_VERT==2)
3270 
3271  lgs_a0(ktmax+kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
3272  lgs_a1(ktmax+kc) = 1.0_dp &
3273  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
3274  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
3275  lgs_a2(ktmax+kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
3276 
3277 #elif (ADV_VERT==3)
3278 
3279  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
3280 
3281  lgs_a0(ktmax+kc) = -max(adv_vert_help, 0.0_dp)
3282  lgs_a1(ktmax+kc) = 1.0_dp &
3283  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp)
3284  lgs_a2(ktmax+kc) = min(adv_vert_help, 0.0_dp)
3285 
3286 #endif
3287 
3288 #if (ADV_HOR==2)
3289 
3290  lgs_b(ktmax+kc) = age_c(kc,j,i) + dtime_temp &
3291  -dtt_2dxi* &
3292  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
3293  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
3294  *insq_g11_sgx(j,i) &
3295  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
3296  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
3297  *insq_g11_sgx(j,i-1) ) &
3298  -dtt_2deta* &
3299  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
3300  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
3301  *insq_g22_sgy(j,i) &
3302  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
3303  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
3304  *insq_g22_sgy(j-1,i) )
3305 
3306 #elif (ADV_HOR==3)
3307 
3308  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
3309  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
3310 
3311  lgs_b(ktmax+kc) = age_c(kc,j,i) + dtime_temp &
3312  -dtt_dxi* &
3313  ( min(vx_c_help, 0.0_dp) &
3314  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
3315  *insq_g11_sgx(j,i) &
3316  +max(vx_c_help, 0.0_dp) &
3317  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
3318  *insq_g11_sgx(j,i-1) ) &
3319  -dtt_deta* &
3320  ( min(vy_c_help, 0.0_dp) &
3321  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
3322  *insq_g22_sgy(j,i) &
3323  +max(vy_c_help, 0.0_dp) &
3324  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
3325  *insq_g22_sgy(j-1,i) )
3326 
3327 #endif
3328 
3329 end do
3330 
3331 kc=kcmax
3332 if (as_perp(j,i) >= zero) then
3333  lgs_a0(ktmax+kc) = 0.0_dp
3334  lgs_a1(ktmax+kc) = 1.0_dp
3335  lgs_b(ktmax+kc) = 0.0_dp
3336 else
3337  lgs_a0(ktmax+kc) = -max(adv_vert_sg(kc-1), 0.0_dp)
3338  lgs_a1(ktmax+kc) = 1.0_dp + max(adv_vert_sg(kc-1), 0.0_dp)
3339  ! adv_vert_sg(KCMAX-1) >= 0 (directed upward)
3340  ! assumed/enforced
3341 #if (ADV_HOR==2)
3342 
3343  lgs_b(ktmax+kc) = age_c(kc,j,i) + dtime_temp &
3344  -dtt_2dxi* &
3345  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
3346  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
3347  *insq_g11_sgx(j,i) &
3348  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
3349  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
3350  *insq_g11_sgx(j,i-1) ) &
3351  -dtt_2deta* &
3352  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
3353  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
3354  *insq_g22_sgy(j,i) &
3355  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
3356  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
3357  *insq_g22_sgy(j-1,i) )
3358 
3359 #elif (ADV_HOR==3)
3360 
3361  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
3362  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
3363 
3364  lgs_b(ktmax+kc) = age_c(kc,j,i) + dtime_temp &
3365  -dtt_dxi* &
3366  ( min(vx_c_help, 0.0_dp) &
3367  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
3368  *insq_g11_sgx(j,i) &
3369  +max(vx_c_help, 0.0_dp) &
3370  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
3371  *insq_g11_sgx(j,i-1) ) &
3372  -dtt_deta* &
3373  ( min(vy_c_help, 0.0_dp) &
3374  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
3375  *insq_g22_sgy(j,i) &
3376  +max(vy_c_help, 0.0_dp) &
3377  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
3378  *insq_g22_sgy(j-1,i) )
3379 
3380 #endif
3381 
3382 end if
3383 
3384 !-------- Solve system of linear equations --------
3385 
3386 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax+ktmax)
3387 
3388 !-------- Assign the result,
3389 ! restriction to interval [0, AGE_MAX yr] --------
3390 
3391 do kt=0, ktmax
3392 
3393  age_t_neu(kt,j,i) = lgs_x(kt)
3394 
3395  if (age_t_neu(kt,j,i) < (age_min*year_sec)) &
3396  age_t_neu(kt,j,i) = 0.0_dp
3397  if (age_t_neu(kt,j,i) > (age_max*year_sec)) &
3398  age_t_neu(kt,j,i) = age_max*year_sec
3399 
3400 end do
3401 
3402 do kc=0, kcmax
3403 
3404  age_c_neu(kc,j,i) = lgs_x(ktmax+kc)
3405 
3406  if (age_c_neu(kc,j,i) < (age_min*year_sec)) &
3407  age_c_neu(kc,j,i) = 0.0_dp
3408  if (age_c_neu(kc,j,i) > (age_max*year_sec)) &
3409  age_c_neu(kc,j,i) = age_max*year_sec
3410 
3411 end do
3412 
3413 end subroutine calc_temp3
3414 
3415 !-------------------------------------------------------------------------------
3416 !> Computation of temperature and age for an ice-free column.
3417 !<------------------------------------------------------------------------------
3418 subroutine calc_temp_r(atr1, alb1, i, j)
3419 
3420 use sico_maths_m, only : tri_sle
3421 
3422 implicit none
3423 
3424 integer(i4b), intent(in) :: i, j
3425 real(dp), intent(in) :: atr1, alb1
3426 
3427 integer(i4b) :: kc, kt, kr
3428 real(dp) :: ctr1, clb1
3429 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
3430  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
3431  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
3432  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
3433  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
3434 
3435 !-------- Abbreviations --------
3436 
3437 ctr1 = atr1
3438 clb1 = alb1*q_geo(j,i)
3439 
3440 !-------- Set up the equations for the bedrock temperature --------
3441 
3442 kr=0
3443 lgs_a1(kr) = 1.0_dp
3444 lgs_a2(kr) = -1.0_dp
3445 lgs_b(kr) = clb1
3446 
3447 #if (Q_LITHO==1)
3448 ! (coupled heat-conducting bedrock)
3449 
3450 do kr=1, krmax-1
3451  lgs_a0(kr) = - ctr1
3452  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
3453  lgs_a2(kr) = - ctr1
3454  lgs_b(kr) = temp_r(kr,j,i)
3455 end do
3456 
3457 #elif (Q_LITHO==0)
3458 ! (no coupled heat-conducting bedrock)
3459 
3460 do kr=1, krmax-1
3461  lgs_a0(kr) = 1.0_dp
3462  lgs_a1(kr) = 0.0_dp
3463  lgs_a2(kr) = -1.0_dp
3464  lgs_b(kr) = 2.0_dp*clb1
3465 end do
3466 
3467 #endif
3468 
3469 kr=krmax
3470 lgs_a0(kr) = 0.0_dp
3471 lgs_a1(kr) = 1.0_dp
3472 lgs_b(kr) = temp_s(j,i)
3473 
3474 !-------- Solve system of linear equations --------
3475 
3476 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
3477 
3478 !-------- Assign the result --------
3479 
3480 do kr=0, krmax
3481  temp_r_neu(kr,j,i) = lgs_x(kr)
3482 end do
3483 
3484 !-------- Water content and age
3485 ! in the non-existing lower (kt) ice layer --------
3486 
3487 do kt=0, ktmax
3488  omega_t_neu(kt,j,i) = 0.0_dp
3489  age_t_neu(kt,j,i) = 0.0_dp
3490 end do
3491 
3492 !-------- Temperature and age
3493 ! in the non-existing upper (kc) ice layer --------
3494 
3495 do kc=0, kcmax
3496  temp_c_neu(kc,j,i) = temp_s(j,i)
3497  age_c_neu(kc,j,i) = 0.0_dp
3498 end do
3499 
3500 end subroutine calc_temp_r
3501 
3502 !-------------------------------------------------------------------------------
3503 !> Upward shifting of the CTS.
3504 !<------------------------------------------------------------------------------
3505 subroutine shift_cts_upward(at1, at2_1, at2_2, at3_1, at3_2, &
3506  at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
3507  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
3508  ai1, ai2, ai3, dzeta_t, &
3509  dtime_temp, dtt_2dxi, dtt_2deta, dtime_temp_inv, &
3510  i, j)
3511 
3512 implicit none
3513 
3514 integer(i4b), intent(in) :: i, j
3515 real(dp), intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
3516  at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
3517  at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
3518  ai1(0:kcmax), ai2(0:kcmax), ai3, &
3519  atr1, am1, am2, alb1
3520 real(dp), intent(in) :: aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld
3521 real(dp), intent(in) :: dzeta_t
3522 real(dp), intent(in) :: dtime_temp, dtime_temp_inv, dtt_2dxi, dtt_2deta
3523 
3524 real(dp) :: zm_shift
3525 real(dp) :: difftemp_a, difftemp_b, interpol
3526 
3527 zm_shift = 1.0_dp ! CTS shift in intervals of 1 m
3528 
3529 !-------- Temperature discrepancy from the computation of the main
3530 ! program --------
3531 
3532 difftemp_a = temp_c_neu(0,j,i)-(-beta*h_c_neu(j,i))
3533 if (difftemp_a <= 0.0_dp) return
3534 
3535 !-------- Shift CTS upward until it is too high --------
3536 
3537  10 continue
3538 
3539  zm_neu(j,i) = zm_neu(j,i) + zm_shift
3540  if (zm_neu(j,i) >= zs(j,i)) then
3541  zm_neu(j,i) = zm_neu(j,i) - zm_shift
3542  return
3543  end if
3544  h_c_neu(j,i) = h_c_neu(j,i) - zm_shift
3545  h_t_neu(j,i) = h_t_neu(j,i) + zm_shift
3546 
3547  dh_t_dtau(j,i) = (zm_neu(j,i)-zm(j,i))*dtime_temp_inv
3548  dzm_dtau(j,i) = dzb_dtau(j,i)+dh_t_dtau(j,i)
3549  dh_c_dtau(j,i) = dzs_dtau(j,i)-dzm_dtau(j,i)
3550 
3551  am_perp(j,i) = am_perp_st(j,i) + dzm_dtau(j,i)
3552 
3553  call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
3554  at4_1, at4_2, at5, at6, at7, atr1, &
3555  am1, am2, alb1, &
3556  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
3557  ai1, ai2, ai3, dzeta_t, &
3558  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
3559 
3560  difftemp_b = difftemp_a
3561  difftemp_a = temp_c_neu(0,j,i)-(-beta*h_c_neu(j,i))
3562 
3563 if (difftemp_a > 0.0_dp) go to 10
3564 
3565 !-------- Interpolate the CTS position from the last (_a) and the
3566 ! last but one (_b) value, weighed with the temperature
3567 ! discrepancies at the CTS --------
3568 
3569 interpol = difftemp_a/(difftemp_b-difftemp_a)*zm_shift
3570 
3571 zm_neu(j,i) = zm_neu(j,i) + interpol
3572 h_c_neu(j,i) = h_c_neu(j,i) - interpol
3573 h_t_neu(j,i) = h_t_neu(j,i) + interpol
3574 
3575 dh_t_dtau(j,i) = (zm_neu(j,i)-zm(j,i))*dtime_temp_inv
3576 dzm_dtau(j,i) = dzb_dtau(j,i)+dh_t_dtau(j,i)
3577 dh_c_dtau(j,i) = dzs_dtau(j,i)-dzm_dtau(j,i)
3578 
3579 am_perp(j,i) = am_perp_st(j,i) + dzm_dtau(j,i)
3580 
3581 call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
3582  at4_1, at4_2, at5, at6, at7, atr1, &
3583  am1, am2, alb1, &
3584  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
3585  ai1, ai2, ai3, dzeta_t, &
3586  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
3587 
3588 end subroutine shift_cts_upward
3589 
3590 !-------------------------------------------------------------------------------
3591 !> Downward shifting of the CTS.
3592 !<------------------------------------------------------------------------------
3593 subroutine shift_cts_downward(at1, at2_1, at2_2, at3_1, at3_2, &
3594  at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
3595  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
3596  ai1, ai2, ai3, dzeta_t, &
3597  dtime_temp, dtt_2dxi, dtt_2deta, dtime_temp_inv, &
3598  i, j)
3599 
3600 implicit none
3601 
3602 integer(i4b), intent(in) :: i, j
3603 real(dp), intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
3604  at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
3605  at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
3606  ai1(0:kcmax), ai2(0:kcmax), ai3, &
3607  atr1, am1, am2, alb1
3608 real(dp), intent(in) :: aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld
3609 real(dp), intent(in) :: dzeta_t
3610 real(dp), intent(in) :: dtt_2dxi, dtt_2deta, dtime_temp, dtime_temp_inv
3611 
3612 real(dp) :: zm_shift
3613 real(dp) :: difftemp_a, difftemp_b, interpol
3614 
3615 zm_shift = 1.0_dp ! CTS shift in intervals of 1 m
3616 
3617 !-------- Temperature discrepancy from the computation of the main
3618 ! program --------
3619 
3620 difftemp_a = temp_c_neu(0,j,i)-(-beta*h_c_neu(j,i))
3621 if (difftemp_a >= 0.0_dp) return
3622 
3623 !-------- Shift CTS downward until it is too low --------
3624 
3625  10 continue
3626 
3627  zm_neu(j,i) = zm_neu(j,i) - zm_shift
3628 
3629 ! ------ Special case: CTS too close to the base
3630 
3631  if (zm_neu(j,i) <= zb(j,i)) then
3632 
3633  zm_shift = (zm_neu(j,i)+zm_shift)-(zb(j,i)+0.001_dp)
3634  zm_neu(j,i) = zb(j,i)+0.001_dp
3635  h_c_neu(j,i) = h_c_neu(j,i)+h_t_neu(j,i)-0.001_dp
3636  h_t_neu(j,i) = 0.001_dp
3637 ! ! CTS positioned 1 mm above ice base --------
3638 
3639  dh_t_dtau(j,i) = (zm_neu(j,i)-zm(j,i))*dtime_temp_inv
3640  dzm_dtau(j,i) = dzb_dtau(j,i)+dh_t_dtau(j,i)
3641  dh_c_dtau(j,i) = dzs_dtau(j,i)-dzm_dtau(j,i)
3642 
3643  am_perp(j,i) = am_perp_st(j,i) + dzm_dtau(j,i)
3644 
3645  call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
3646  at4_1, at4_2, at5, at6, at7, atr1, &
3647  am1, am2, alb1, &
3648  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
3649  ai1, ai2, ai3, dzeta_t, &
3650  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
3651 
3652  difftemp_b = difftemp_a
3653  difftemp_a = temp_c_neu(0,j,i)-(-beta*h_c_neu(j,i))
3654 
3655  if (difftemp_a >= 0.0_dp) then ! CTS remains above the base
3656 
3657 ! ---- Interpolate the CTS position from the last (_a) and the
3658 ! last but one (_b) value, weighed with the temperature
3659 ! discrepancies at the CTS --------
3660 
3661  interpol = difftemp_a/(difftemp_a-difftemp_b)*zm_shift
3662 
3663  zm_neu(j,i) = zm_neu(j,i) + interpol
3664  h_c_neu(j,i) = h_c_neu(j,i) - interpol
3665  h_t_neu(j,i) = h_t_neu(j,i) + interpol
3666 
3667  dh_t_dtau(j,i) = (zm_neu(j,i)-zm(j,i))*dtime_temp_inv
3668  dzm_dtau(j,i) = dzb_dtau(j,i)+dh_t_dtau(j,i)
3669  dh_c_dtau(j,i) = dzs_dtau(j,i)-dzm_dtau(j,i)
3670 
3671  am_perp(j,i) = am_perp_st(j,i) + dzm_dtau(j,i)
3672 
3673  call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
3674  at4_1, at4_2, at5, at6, at7, atr1, &
3675  am1, am2, alb1, &
3676  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
3677  ai1, ai2, ai3, dzeta_t, &
3678  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
3679 
3680  else ! CTS disappears
3681 
3682  n_cts_neu(j,i) = 0
3683  zm_neu(j,i) = zb(j,i)
3684  h_c_neu(j,i) = h_c_neu(j,i)+h_t_neu(j,i)
3685  h_t_neu(j,i) = 0.0_dp
3686 
3687  dh_t_dtau(j,i) = (zm_neu(j,i)-zm(j,i))*dtime_temp_inv
3688  dzm_dtau(j,i) = dzb_dtau(j,i)+dh_t_dtau(j,i)
3689  dh_c_dtau(j,i) = dzs_dtau(j,i)-dzm_dtau(j,i)
3690 
3691  am_perp(j,i) = am_perp_st(j,i) + dzm_dtau(j,i)
3692 
3693  call calc_temp2(at1, at2_1, at2_2, at3_1, at3_2, &
3694  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
3695  ai1, ai2, &
3696  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
3697 
3698  end if
3699 
3700  return
3701  end if
3702 
3703 ! ------ End of treatment of special case
3704 
3705  h_c_neu(j,i) = h_c_neu(j,i) + zm_shift
3706  h_t_neu(j,i) = h_t_neu(j,i) - zm_shift
3707 
3708  dh_t_dtau(j,i) = (zm_neu(j,i)-zm(j,i))*dtime_temp_inv
3709  dzm_dtau(j,i) = dzb_dtau(j,i)+dh_t_dtau(j,i)
3710  dh_c_dtau(j,i) = dzs_dtau(j,i)-dzm_dtau(j,i)
3711 
3712  am_perp(j,i) = am_perp_st(j,i) + dzm_dtau(j,i)
3713 
3714  call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
3715  at4_1, at4_2, at5, at6, at7, atr1, &
3716  am1, am2, alb1, &
3717  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
3718  ai1, ai2, ai3, dzeta_t, &
3719  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
3720 
3721  difftemp_b = difftemp_a
3722  difftemp_a = temp_c_neu(0,j,i)-(-beta*h_c_neu(j,i))
3723 
3724 if (difftemp_a < 0.0_dp) go to 10
3725 
3726 !-------- Interpolate the CTS position from the last (_a) and the
3727 ! last but one (_b) value, weighed with the temperature
3728 ! discrepancies at the CTS --------
3729 
3730 interpol = difftemp_a/(difftemp_a-difftemp_b)*zm_shift
3731 
3732 zm_neu(j,i) = zm_neu(j,i) + interpol
3733 h_c_neu(j,i) = h_c_neu(j,i) - interpol
3734 h_t_neu(j,i) = h_t_neu(j,i) + interpol
3735 
3736 dh_t_dtau(j,i) = (zm_neu(j,i)-zm(j,i))*dtime_temp_inv
3737 dzm_dtau(j,i) = dzb_dtau(j,i)+dh_t_dtau(j,i)
3738 dh_c_dtau(j,i) = dzs_dtau(j,i)-dzm_dtau(j,i)
3739 
3740 am_perp(j,i) = am_perp_st(j,i) + dzm_dtau(j,i)
3741 
3742 call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
3743  at4_1, at4_2, at5, at6, at7, atr1, &
3744  am1, am2, alb1, &
3745  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
3746  ai1, ai2, ai3, dzeta_t, &
3747  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
3748 
3749 end subroutine shift_cts_downward
3750 
3751 !-------------------------------------------------------------------------------
3752 !> Computation of temperature and age for ice shelves (floating ice).
3753 !<------------------------------------------------------------------------------
3754 subroutine calc_temp_ssa(at1, at2_1, at2_2, at3_1, at3_2, &
3755  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
3756  ai1, ai2, &
3757  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
3758 
3759 use ice_material_properties_m, only : kappa_val, c_val, viscosity
3760 use sico_maths_m, only : tri_sle
3761 
3762 implicit none
3763 
3764 integer(i4b), intent(in) :: i, j
3765 real(dp), intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
3766  at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
3767  at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
3768  ai1(0:kcmax), ai2(0:kcmax), &
3769  atr1, alb1
3770 real(dp), intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
3771 
3772 integer(i4b) :: kc, kt, kr
3773 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
3774  ct5(0:kcmax), ct6(0:kcmax), ct7(0:kcmax), ctr1, clb1
3775 real(dp) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), ct3_sg(0:kcmax), &
3776  ct4_sg(0:kcmax), adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
3777 real(dp) :: ci1(0:kcmax), ci2(0:kcmax)
3778 real(dp) :: temp_c_help(0:kcmax)
3779 real(dp) :: vx_c_help, vy_c_help
3780 real(dp) :: adv_vert_help
3781 real(dp) :: dtt_dxi, dtt_deta
3782 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
3783  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
3784  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
3785  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
3786  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
3787 real(dp), parameter :: zero=0.0_dp
3788 
3789 !-------- Check for boundary points --------
3790 
3791 if ((i == 0).or.(i == imax).or.(j == 0).or.(j == jmax)) &
3792  stop ' >>> calc_temp_ssa: Boundary points not allowed.'
3793 
3794 !-------- Abbreviations --------
3795 
3796 ctr1 = atr1
3797 clb1 = alb1*q_geo(j,i)
3798 
3799 #if (ADV_VERT==1)
3800 
3801 do kc=1, kcmax-1
3802  ct1(kc) = at1(kc)/h_c(j,i)*0.5_dp*(vz_c(kc,j,i)+vz_c(kc-1,j,i))
3803 end do
3804 
3805 kc=0
3806 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
3807  ! only needed for kc=0 ...
3808 kc=kcmax-1
3809 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
3810  ! ... and kc=KCMAX-1
3811 
3812 #elif (ADV_VERT==2 || ADV_VERT==3)
3813 
3814 do kc=0, kcmax-1
3815  ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
3816 end do
3817 
3818 #endif
3819 
3820 do kc=0, kcmax
3821 
3822  ct2(kc) = ( at2_1(kc)*dzm_dtau(j,i) &
3823  +at2_2(kc)*dh_c_dtau(j,i) )/h_c(j,i)
3824  ct3(kc) = ( at3_1(kc)*dzm_dxi_g(j,i) &
3825  +at3_2(kc)*dh_c_dxi_g(j,i) )/h_c(j,i) &
3826  *0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1)) *insq_g11_g(j,i)
3827  ct4(kc) = ( at4_1(kc)*dzm_deta_g(j,i) &
3828  +at4_2(kc)*dh_c_deta_g(j,i) )/h_c(j,i) &
3829  *0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i)) *insq_g22_g(j,i)
3830  ct5(kc) = at5(kc) &
3831  /c_val(temp_c(kc,j,i)) &
3832  /h_c(j,i)
3833  ct7(kc) = 2.0_dp*at7 &
3834  /c_val(temp_c(kc,j,i)) &
3835  *viscosity(de_ssa(j,i), &
3836  temp_c(kc,j,i), temp_c_m(kc,j,i), 0.0_dp, &
3837  enh_c(kc,j,i), 0_i2b) &
3838  *de_ssa(j,i)**2
3839  ci1(kc) = ai1(kc)/h_c(j,i)
3840 
3841 end do
3842 
3843 #if (ADV_VERT==1)
3844 
3845 kc=0
3846 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
3847 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
3848 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
3849 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
3850 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! only needed for kc=0 ...
3851 kc=kcmax-1
3852 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
3853 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
3854 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
3855 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
3856 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! ... and kc=KCMAX-1
3857 
3858 #elif (ADV_VERT==2 || ADV_VERT==3)
3859 
3860 do kc=0, kcmax-1
3861  ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
3862  ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
3863  ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
3864  adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
3865  abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
3866 end do
3867 
3868 #endif
3869 
3870 do kc=0, kcmax-1
3871  temp_c_help(kc) = 0.5_dp*(temp_c(kc,j,i)+temp_c(kc+1,j,i))
3872  ct6(kc) = at6(kc) &
3873  *kappa_val(temp_c_help(kc)) &
3874  /h_c(j,i)
3875  ci2(kc) = ai2(kc)/h_c(j,i)
3876 end do
3877 
3878 #if (ADV_HOR==3)
3879 dtt_dxi = 2.0_dp*dtt_2dxi
3880 dtt_deta = 2.0_dp*dtt_2deta
3881 #endif
3882 
3883 !-------- Set up the equations for the bedrock temperature --------
3884 
3885 kr=0
3886 lgs_a1(kr) = 1.0_dp
3887 lgs_a2(kr) = -1.0_dp
3888 lgs_b(kr) = clb1
3889 
3890 #if (Q_LITHO==1)
3891 ! (coupled heat-conducting bedrock)
3892 
3893 do kr=1, krmax-1
3894  lgs_a0(kr) = - ctr1
3895  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
3896  lgs_a2(kr) = - ctr1
3897  lgs_b(kr) = temp_r(kr,j,i)
3898 end do
3899 
3900 #elif (Q_LITHO==0)
3901 ! (no coupled heat-conducting bedrock)
3902 
3903 do kr=1, krmax-1
3904  lgs_a0(kr) = 1.0_dp
3905  lgs_a1(kr) = 0.0_dp
3906  lgs_a2(kr) = -1.0_dp
3907  lgs_b(kr) = 2.0_dp*clb1
3908 end do
3909 
3910 #endif
3911 
3912 kr=krmax
3913 lgs_a0(kr) = 0.0_dp
3914 lgs_a1(kr) = 1.0_dp
3915 lgs_b(kr) = temp_c_m(0,j,i)-delta_tm_sw
3916 
3917 !-------- Solve system of linear equations --------
3918 
3919 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
3920 
3921 !-------- Assign the result --------
3922 
3923 do kr=0, krmax
3924  temp_r_neu(kr,j,i) = lgs_x(kr)
3925 end do
3926 
3927 !-------- Set up the equations for the ice temperature --------
3928 
3929 kc=0
3930 lgs_a1(kc) = 1.0_dp
3931 lgs_a2(kc) = 0.0_dp
3932 lgs_b(kc) = temp_c_m(0,j,i)-delta_tm_sw
3933 
3934 do kc=1, kcmax-1
3935 
3936 #if (ADV_VERT==1)
3937 
3938  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
3939  -ct5(kc)*ct6(kc-1)
3940  lgs_a1(kc) = 1.0_dp+ct5(kc)*(ct6(kc)+ct6(kc-1))
3941  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
3942  -ct5(kc)*ct6(kc)
3943 
3944 #elif (ADV_VERT==2)
3945 
3946  lgs_a0(kc) &
3947  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
3948  -ct5(kc)*ct6(kc-1)
3949  lgs_a1(kc) &
3950  = 1.0_dp &
3951  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
3952  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
3953  +ct5(kc)*(ct6(kc)+ct6(kc-1))
3954  lgs_a2(kc) &
3955  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
3956  -ct5(kc)*ct6(kc)
3957 
3958 #elif (ADV_VERT==3)
3959 
3960  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
3961 
3962  lgs_a0(kc) &
3963  = -max(adv_vert_help, 0.0_dp) &
3964  -ct5(kc)*ct6(kc-1)
3965  lgs_a1(kc) &
3966  = 1.0_dp &
3967  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
3968  +ct5(kc)*(ct6(kc)+ct6(kc-1))
3969  lgs_a2(kc) &
3970  = min(adv_vert_help, 0.0_dp) &
3971  -ct5(kc)*ct6(kc)
3972 
3973 #endif
3974 
3975 #if (ADV_HOR==2)
3976 
3977  lgs_b(kc) = temp_c(kc,j,i) + ct7(kc) &
3978  -dtt_2dxi* &
3979  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
3980  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
3981  *insq_g11_sgx(j,i) &
3982  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
3983  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
3984  *insq_g11_sgx(j,i-1) ) &
3985  -dtt_2deta* &
3986  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
3987  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
3988  *insq_g22_sgy(j,i) &
3989  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
3990  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
3991  *insq_g22_sgy(j-1,i) )
3992 
3993 #elif (ADV_HOR==3)
3994 
3995  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
3996  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
3997 
3998  lgs_b(kc) = temp_c(kc,j,i) + ct7(kc) &
3999  -dtt_dxi* &
4000  ( min(vx_c_help, 0.0_dp) &
4001  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
4002  *insq_g11_sgx(j,i) &
4003  +max(vx_c_help, 0.0_dp) &
4004  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
4005  *insq_g11_sgx(j,i-1) ) &
4006  -dtt_deta* &
4007  ( min(vy_c_help, 0.0_dp) &
4008  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
4009  *insq_g22_sgy(j,i) &
4010  +max(vy_c_help, 0.0_dp) &
4011  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
4012  *insq_g22_sgy(j-1,i) )
4013 
4014 #endif
4015 
4016 end do
4017 
4018 kc=kcmax
4019 lgs_a0(kc) = 0.0_dp
4020 lgs_a1(kc) = 1.0_dp
4021 lgs_b(kc) = temp_s(j,i)
4022 
4023 !-------- Solve system of linear equations --------
4024 
4025 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
4026 
4027 !-------- Assign the result --------
4028 
4029 do kc=0, kcmax
4030  temp_c_neu(kc,j,i) = lgs_x(kc)
4031 end do
4032 
4033 !-------- Set water content in the non-existing temperate layer
4034 ! to zero --------
4035 
4036 do kt=0, ktmax
4037  omega_t_neu(kt,j,i) = 0.0_dp
4038 end do
4039 
4040 !-------- Water drainage from the non-existing temperate layer --------
4041 
4042 q_tld(j,i) = 0.0_dp
4043 
4044 !-------- Set up the equations for the age of cold ice --------
4045 
4046 kc=0 ! adv_vert_sg(0) <= 0
4047 lgs_a1(kc) = 1.0_dp - min(adv_vert_sg(kc), 0.0_dp) ! (directed downward)
4048 lgs_a2(kc) = min(adv_vert_sg(kc), 0.0_dp) ! assumed/enforced
4049 
4050 #if (ADV_HOR==2)
4051 
4052 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
4053  -dtt_2dxi* &
4054  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
4055  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
4056  *insq_g11_sgx(j,i) &
4057  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
4058  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
4059  *insq_g11_sgx(j,i-1) ) &
4060  -dtt_2deta* &
4061  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
4062  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
4063  *insq_g22_sgy(j,i) &
4064  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
4065  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
4066  *insq_g22_sgy(j-1,i) )
4067 
4068 #elif (ADV_HOR==3)
4069 
4070 vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
4071 vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
4072 
4073 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
4074  -dtt_dxi* &
4075  ( min(vx_c_help, 0.0_dp) &
4076  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
4077  *insq_g11_sgx(j,i) &
4078  +max(vx_c_help, 0.0_dp) &
4079  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
4080  *insq_g11_sgx(j,i-1) ) &
4081  -dtt_deta* &
4082  ( min(vy_c_help, 0.0_dp) &
4083  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
4084  *insq_g22_sgy(j,i) &
4085  +max(vy_c_help, 0.0_dp) &
4086  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
4087  *insq_g22_sgy(j-1,i) )
4088 
4089 #endif
4090 
4091 do kc=1, kcmax-1
4092 
4093 #if (ADV_VERT==1)
4094 
4095  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
4096  -ci1(kc)*ci2(kc-1)
4097  lgs_a1(kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
4098  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
4099  -ci1(kc)*ci2(kc)
4100 
4101 #elif (ADV_VERT==2)
4102 
4103  lgs_a0(kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
4104  lgs_a1(kc) = 1.0_dp &
4105  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
4106  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
4107  lgs_a2(kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
4108 
4109 #elif (ADV_VERT==3)
4110 
4111  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
4112 
4113  lgs_a0(kc) = -max(adv_vert_help, 0.0_dp)
4114  lgs_a1(kc) = 1.0_dp &
4115  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp)
4116  lgs_a2(kc) = min(adv_vert_help, 0.0_dp)
4117 
4118 #endif
4119 
4120 #if (ADV_HOR==2)
4121 
4122  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
4123  -dtt_2dxi* &
4124  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
4125  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
4126  *insq_g11_sgx(j,i) &
4127  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
4128  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
4129  *insq_g11_sgx(j,i-1) ) &
4130  -dtt_2deta* &
4131  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
4132  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
4133  *insq_g22_sgy(j,i) &
4134  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
4135  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
4136  *insq_g22_sgy(j-1,i) )
4137 
4138 #elif (ADV_HOR==3)
4139 
4140  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
4141  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
4142 
4143  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
4144  -dtt_dxi* &
4145  ( min(vx_c_help, 0.0_dp) &
4146  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
4147  *insq_g11_sgx(j,i) &
4148  +max(vx_c_help, 0.0_dp) &
4149  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
4150  *insq_g11_sgx(j,i-1) ) &
4151  -dtt_deta* &
4152  ( min(vy_c_help, 0.0_dp) &
4153  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
4154  *insq_g22_sgy(j,i) &
4155  +max(vy_c_help, 0.0_dp) &
4156  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
4157  *insq_g22_sgy(j-1,i) )
4158 
4159 #endif
4160 
4161 end do
4162 
4163 kc=kcmax
4164 if (as_perp(j,i) >= zero) then
4165  lgs_a0(kc) = 0.0_dp
4166  lgs_a1(kc) = 1.0_dp
4167  lgs_b(kc) = 0.0_dp
4168 else
4169  lgs_a0(kc) = -max(adv_vert_sg(kc-1), 0.0_dp)
4170  lgs_a1(kc) = 1.0_dp + max(adv_vert_sg(kc-1), 0.0_dp)
4171  ! adv_vert_sg(KCMAX-1) >= 0 (directed upward)
4172  ! assumed/enforced
4173 #if (ADV_HOR==2)
4174 
4175  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
4176  -dtt_2dxi* &
4177  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
4178  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
4179  *insq_g11_sgx(j,i) &
4180  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
4181  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
4182  *insq_g11_sgx(j,i-1) ) &
4183  -dtt_2deta* &
4184  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
4185  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
4186  *insq_g22_sgy(j,i) &
4187  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
4188  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
4189  *insq_g22_sgy(j-1,i) )
4190 
4191 #elif (ADV_HOR==3)
4192 
4193  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
4194  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
4195 
4196  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
4197  -dtt_dxi* &
4198  ( min(vx_c_help, 0.0_dp) &
4199  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
4200  *insq_g11_sgx(j,i) &
4201  +max(vx_c_help, 0.0_dp) &
4202  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
4203  *insq_g11_sgx(j,i-1) ) &
4204  -dtt_deta* &
4205  ( min(vy_c_help, 0.0_dp) &
4206  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
4207  *insq_g22_sgy(j,i) &
4208  +max(vy_c_help, 0.0_dp) &
4209  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
4210  *insq_g22_sgy(j-1,i) )
4211 
4212 #endif
4213 
4214 end if
4215 
4216 !-------- Solve system of linear equations --------
4217 
4218 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
4219 
4220 !-------- Assign the result,
4221 ! restriction to interval [0, AGE_MAX yr] --------
4222 
4223 do kc=0, kcmax
4224 
4225  age_c_neu(kc,j,i) = lgs_x(kc)
4226 
4227  if (age_c_neu(kc,j,i) < (age_min*year_sec)) &
4228  age_c_neu(kc,j,i) = 0.0_dp
4229  if (age_c_neu(kc,j,i) > (age_max*year_sec)) &
4230  age_c_neu(kc,j,i) = age_max*year_sec
4231 
4232 end do
4233 
4234 !-------- Age of the ice in the non-existing temperate layer --------
4235 
4236 do kt=0, ktmax
4237  age_t_neu(kt,j,i) = age_c_neu(0,j,i)
4238 end do
4239 
4240 end subroutine calc_temp_ssa
4241 
4242 !-------------------------------------------------------------------------------
4243 
4244 end module calc_temp_m
4245 !
sico_variables_m::vx_c
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vx_c
vx_c(kc,j,i): Velocity in x-direction in the upper (kc) ice domain (at (i+1/2,j,kc)) ...
Definition: sico_variables_m.F90:402
sico_variables_m::temp_c
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c
temp_c(kc,j,i): Temperature in the upper (kc) ice domain
Definition: sico_variables_m.F90:410
sico_variables_m::kappa_r
real(dp) kappa_r
KAPPA_R: Heat conductivity of the lithosphere.
Definition: sico_variables_m.F90:546
sico_variables_m::delta_tm_sw
real(dp) delta_tm_sw
DELTA_TM_SW: Melting point depression of sea water due to its average salinity.
Definition: sico_variables_m.F90:538
sico_variables_m::temp_r_neu
real(dp), dimension(0:krmax, 0:jmax, 0:imax) temp_r_neu
(.)_neu: New value of quantity (.) computed during an integration step
Definition: sico_variables_m.F90:474
sico_variables_m::de_c
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) de_c
de_c(kc,j,i): Full effective strain rate in the upper (kc) ice domain
Definition: sico_variables_m.F90:501
sico_variables_m::age_c
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) age_c
age_c(kc,j,i): Age in the upper (kc) ice domain
Definition: sico_variables_m.F90:416
sico_variables_m::n_cts
integer(i2b), dimension(0:jmax, 0:imax) n_cts
n_cts(j,i): Mask for thermal conditions. -1: cold ice base, 0: temperate ice base with cold ice above...
Definition: sico_variables_m.F90:65
sico_variables_m::insq_g22_sgy
real(dp), dimension(0:jmax, 0:imax) insq_g22_sgy
insq_g22_sgy(j,i): Inverse square root of g22, at (i,j+1/2)
Definition: sico_variables_m.F90:165
sico_variables_m::de_t
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) de_t
de_t(kt,j,i): Full effective strain rate in the lower (kt) ice domain
Definition: sico_variables_m.F90:516
sico_variables_m::eaz_c
real(dp), dimension(0:kcmax) eaz_c
eaz_c(kc): Abbreviation for exp(aa*zeta(kc))
Definition: sico_variables_m.F90:134
sico_variables_m::q_geo
real(dp), dimension(0:jmax, 0:imax) q_geo
q_geo(j,i): Geothermal heat flux
Definition: sico_variables_m.F90:299
sico_variables_m::vx_b_g
real(dp), dimension(0:jmax, 0:imax) vx_b_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
Definition: sico_variables_m.F90:274
sico_variables_m::age_t
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) age_t
age_t(kt,j,i): Age in the lower (kt) ice domain
Definition: sico_variables_m.F90:457
sico_variables_m::enh_c
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) enh_c
enh_c(kc,j,i): Flow enhancement factor in the upper (kc) ice domain
Definition: sico_variables_m.F90:428
sico_variables_m::temp_r
real(dp), dimension(0:krmax, 0:jmax, 0:imax) temp_r
temp_r(kr,j,i): Temperature in the bedrock
Definition: sico_variables_m.F90:472
sico_variables_m::age_t_neu
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) age_t_neu
(.)_neu: New value of quantity (.) computed during an integration step
Definition: sico_variables_m.F90:459
sico_variables_m::h_r
real(dp) h_r
H_R: Thickness of the modelled lithosphere layer.
Definition: sico_variables_m.F90:542
ice_material_properties_m
Material properties of ice: Rate factor, heat conductivity, specific heat (heat capacity), creep function, viscosity.
Definition: ice_material_properties_m.F90:39
calc_temp_m::calc_temp_r
subroutine calc_temp_r(atr1, alb1, i, j)
Computation of temperature and age for an ice-free column.
Definition: calc_temp_m.F90:3419
sico_variables_m::eps
real(dp), parameter eps
eps: Small number
Definition: sico_variables_m.F90:610
sico_variables_m::omega_max
real(dp) omega_max
OMEGA_MAX: Threshold value for the water content of temperate ice.
Definition: sico_variables_m.F90:540
sico_variables_m::insq_g11_sgx
real(dp), dimension(0:jmax, 0:imax) insq_g11_sgx
insq_g11_sgx(j,i): Inverse square root of g11, at (i+1/2,j)
Definition: sico_variables_m.F90:163
sico_maths_m
Several mathematical tools used by SICOPOLIS.
Definition: sico_maths_m.F90:35
sico_variables_m::insq_g22_g
real(dp), dimension(0:jmax, 0:imax) insq_g22_g
insq_g22_g(j,i): Inverse square root of g22 on grid point (i,j)
Definition: sico_variables_m.F90:153
sico_variables_m::dh_t_dtau
real(dp), dimension(0:jmax, 0:imax) dh_t_dtau
dH_t_dtau(j,i): Derivative of H_t by tau (time)
Definition: sico_variables_m.F90:244
sico_variables_m::dh_c_dtau
real(dp), dimension(0:jmax, 0:imax) dh_c_dtau
dH_c_dtau(j,i): Derivative of H_c by tau (time)
Definition: sico_variables_m.F90:242
ice_material_properties_m::kappa_val
real(dp) function, public kappa_val(temp_val)
Heat conductivity of ice: Linear interpolation of tabulated values in KAPPA(.).
Definition: ice_material_properties_m.F90:247
sico_variables_m::maske
integer(i2b), dimension(0:jmax, 0:imax) maske
maske(j,i): Ice-land-ocean mask. 0: grounded ice, 1: ice-free land, 2: ocean, 3: floating ice ...
Definition: sico_variables_m.F90:49
sico_variables_m::dh_c_dxi_g
real(dp), dimension(0:jmax, 0:imax) dh_c_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
Definition: sico_variables_m.F90:220
calc_temp_m::calc_temp_cold
subroutine, public calc_temp_cold(dxi, deta, dzeta_c, dzeta_t, dzeta_r, dtime_temp)
Computation of temperature and age in cold-ice mode.
Definition: calc_temp_m.F90:820
calc_temp_m::calc_temp2
subroutine calc_temp2(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 an ice column with a temperate base overlain by cold ice...
Definition: calc_temp_m.F90:1912
sico_variables_m::as_perp
real(dp), dimension(0:jmax, 0:imax) as_perp
as_perp(j,i): Accumulation-ablation function at the ice surface (SMB)
Definition: sico_variables_m.F90:320
sico_variables_m::area
real(dp), dimension(0:jmax, 0:imax) area
area(j,i): Area of grid cell associated with grid point (i,j)
Definition: sico_variables_m.F90:143
sico_variables_m::dzs_deta_g
real(dp), dimension(0:jmax, 0:imax) dzs_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
Definition: sico_variables_m.F90:224
sico_variables_m::temp_t_m
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) temp_t_m
temp_t_m(kt,j,i): Melting temperature in the lower (kt) ice domain
Definition: sico_variables_m.F90:455
sico_variables_m::de_ssa
real(dp), dimension(0:jmax, 0:imax) de_ssa
de_ssa(j,i): Effective strain rate of the SSA, at (i,j)
Definition: sico_variables_m.F90:431
sico_variables_m::h_t_neu
real(dp), dimension(0:jmax, 0:imax) h_t_neu
(.)_neu: New value of quantity (.) computed during an integration step
Definition: sico_variables_m.F90:351
calc_temp_m::calc_temp_const
subroutine, public calc_temp_const()
Isothermal mode: Setting of the temperature and age to constant values.
Definition: calc_temp_m.F90:1342
sico_maths_m::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_maths_m.F90:134
ice_material_properties_m::creep
real(dp) function, public creep(sigma_val)
Creep response function for ice.
Definition: ice_material_properties_m.F90:330
sico_variables_m::temp_s
real(dp), dimension(0:jmax, 0:imax) temp_s
temp_s(j,i): Ice surface temperature
Definition: sico_variables_m.F90:334
sico_variables_m::vz_t
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) vz_t
vz_t(kt,j,i): Velocity in z-direction in the lower (kt) ice domain (at (i,j,kt+1/2)) ...
Definition: sico_variables_m.F90:449
sico_variables_m::g
real(dp) g
G: Acceleration due to gravity.
Definition: sico_variables_m.F90:531
sico_variables_m::q_tld
real(dp), dimension(0:jmax, 0:imax) q_tld
Q_tld(j,i): Water drainage rate from the temperate layer.
Definition: sico_variables_m.F90:307
sico_vars_m
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars_m.F90:35
calc_temp_m::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_temp_m.F90:2421
ice_material_properties_m::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_properties_m.F90:171
sico_variables_m::ea
real(dp) ea
ea: Abbreviation for exp(aa)
Definition: sico_variables_m.F90:132
sico_variables_m::flag_aa_nonzero
logical flag_aa_nonzero
flag_aa_nonzero: Flag for the exponential stretch parameter aa. .true.: aa greater than zero (non-equ...
Definition: sico_variables_m.F90:130
sico_variables_m::age_c_neu
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) age_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
Definition: sico_variables_m.F90:418
sico_variables_m::aa
real(dp) aa
aa: Exponential stretch parameter of the non-equidistant vertical grid in the upper (kc) ice domain ...
Definition: sico_variables_m.F90:126
sico_variables_m::rho_c_r
real(dp) rho_c_r
RHO_C_R: Density times specific heat of the lithosphere.
Definition: sico_variables_m.F90:544
sico_types_m
Declarations of kind types for SICOPOLIS.
Definition: sico_types_m.F90:35
sico_variables_m::zs
real(dp), dimension(0:jmax, 0:imax) zs
zs(j,i): Coordinate z of the surface topography
Definition: sico_variables_m.F90:167
sico_variables_m::dh_c_deta_g
real(dp), dimension(0:jmax, 0:imax) dh_c_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
Definition: sico_variables_m.F90:230
sico_variables_m::enh_t
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) enh_t
enh_t(kt,j,i): Flow enhancement factor in the lower (kt) ice domain
Definition: sico_variables_m.F90:469
sico_variables_m::kc_cts_neu
integer(i2b), dimension(0:jmax, 0:imax) kc_cts_neu
(.)_neu: New value of quantity (.) computed during an integration step
Definition: sico_variables_m.F90:71
ice_material_properties_m::ratefac_t
real(dp) function, public ratefac_t(omega_val)
Rate factor for temperate ice.
Definition: ice_material_properties_m.F90:199
sico_variables_m::temp_c_neu
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
Definition: sico_variables_m.F90:412
sico_variables_m::h_c
real(dp), dimension(0:jmax, 0:imax) h_c
H_c(j,i): Thickness of ice in the upper (kc) domain (thickness of the cold-ice layer for POLY...
Definition: sico_variables_m.F90:188
calc_temp_m
Computation of temperature, water content and age.
Definition: calc_temp_m.F90:35
sico_variables_m::vy_b_g
real(dp), dimension(0:jmax, 0:imax) vy_b_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
Definition: sico_variables_m.F90:276
sico_variables_m::omega_t_neu
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) omega_t_neu
(.)_neu: New value of quantity (.) computed during an integration step
Definition: sico_variables_m.F90:453
sico_variables_m::dzm_dtau
real(dp), dimension(0:jmax, 0:imax) dzm_dtau
dzm_dtau(j,i): Derivative of zm by tau (time)
Definition: sico_variables_m.F90:236
sico_variables_m::temp_c_m
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c_m
temp_c_m(kc,j,i): Melting temperature in the upper (kc) ice domain
Definition: sico_variables_m.F90:414
sico_variables_m::dzb_deta_g
real(dp), dimension(0:jmax, 0:imax) dzb_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
Definition: sico_variables_m.F90:228
sico_variables_m::nue
real(dp) nue
NUE: Water diffusivity in ice.
Definition: sico_variables_m.F90:533
ice_material_properties_m::c_val
real(dp) function, public c_val(temp_val)
Specific heat of ice: Linear interpolation of tabulated values in C(.).
Definition: ice_material_properties_m.F90:289
sico_variables_m::p_weert_inv
real(dp), dimension(0:jmax, 0:imax) p_weert_inv
p_weert_inv(j,i): Inverse of p_weert
Definition: sico_variables_m.F90:250
sico_variables_m::flag_shelfy_stream
logical, dimension(0:jmax, 0:imax) flag_shelfy_stream
flag_shelfy_stream(j,i): Shelfy stream flag on the main grid. .true.: grounded ice, and at least one neighbour on the staggered grid is a shelfy stream point .false.: otherwise
Definition: sico_variables_m.F90:113
sico_variables_m::vz_c
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vz_c
vz_c(kc,j,i): Velocity in z-direction in the upper (kc) ice domain (at (i,j,kc+1/2)) ...
Definition: sico_variables_m.F90:408
sico_variables_m::h_c_neu
real(dp), dimension(0:jmax, 0:imax) h_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
Definition: sico_variables_m.F90:349
sico_variables_m::zm
real(dp), dimension(0:jmax, 0:imax) zm
zm(j,i): Coordinate z of the bottom of the upper (kc) ice domain = top of the lower (kt) ice domain (...
Definition: sico_variables_m.F90:172
sico_variables_m::eaz_c_quotient
real(dp), dimension(0:kcmax) eaz_c_quotient
eaz_c_quotient(kc): Abbreviation for (eaz_c(kc)-1.0)/(ea-1.0)
Definition: sico_variables_m.F90:136
sico_variables_m::dh_t_deta_g
real(dp), dimension(0:jmax, 0:imax) dh_t_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
Definition: sico_variables_m.F90:232
sico_variables_m::dh_t_dxi_g
real(dp), dimension(0:jmax, 0:imax) dh_t_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
Definition: sico_variables_m.F90:222
calc_temp_m::shift_cts_upward
subroutine shift_cts_upward(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, dtime_temp_inv, i, j)
Upward shifting of the CTS.
Definition: calc_temp_m.F90:3511
sico_variables_m::dzm_dxi_g
real(dp), dimension(0:jmax, 0:imax) dzm_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
Definition: sico_variables_m.F90:216
sico_variables_m::dzs_dxi_g
real(dp), dimension(0:jmax, 0:imax) dzs_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
Definition: sico_variables_m.F90:214
sico_variables_m::vy_t
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) vy_t
vy_t(kt,j,i): Velocity in y-direction in the lower (kt) ice domain (at (i,j+1/2,kt)) ...
Definition: sico_variables_m.F90:446
sico_variables_m::am_perp
real(dp), dimension(0:jmax, 0:imax) am_perp
am_perp(j,i): Ice volume flux across the z=zm interface
Definition: sico_variables_m.F90:336
sico_variables_m::beta
real(dp) beta
BETA: Clausius-Clapeyron gradient of ice.
Definition: sico_variables_m.F90:535
sico_variables_m::dzm_deta_g
real(dp), dimension(0:jmax, 0:imax) dzm_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
Definition: sico_variables_m.F90:226
sico_variables_m::insq_g11_g
real(dp), dimension(0:jmax, 0:imax) insq_g11_g
insq_g11_g(j,i): Inverse square root of g11 on grid point (i,j)
Definition: sico_variables_m.F90:151
sico_variables_m::c_drag
real(dp), dimension(0:jmax, 0:imax) c_drag
c_drag(j,i): Auxiliary quantity for the computation of the basal drag
Definition: sico_variables_m.F90:256
sico_variables_m::vy_c
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vy_c
vy_c(kc,j,i): Velocity in y-direction in the upper (kc) ice domain (at (i,j+1/2,kc)) ...
Definition: sico_variables_m.F90:405
sico_variables_m::omega_t
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) omega_t
omega_t(kt,j,i): Water content in the lower (kt) ice domain
Definition: sico_variables_m.F90:451
sico_variables_m::zeta_t
real(dp), dimension(0:ktmax) zeta_t
zeta_t(kt): Sigma coordinate zeta_t of grid point kt
Definition: sico_variables_m.F90:121
ice_material_properties_m::viscosity
real(dp) function, public 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: ice_material_properties_m.F90:406
sico_variables_m::l
real(dp) l
L: Latent heat of ice.
Definition: sico_variables_m.F90:529
sico_variables_m::zm_neu
real(dp), dimension(0:jmax, 0:imax) zm_neu
(.)_neu: New value of quantity (.) computed during an integration step
Definition: sico_variables_m.F90:343
sico_variables_m::rho
real(dp) rho
RHO: Density of ice.
Definition: sico_variables_m.F90:523
sico_variables_m::n_cts_neu
integer(i2b), dimension(0:jmax, 0:imax) n_cts_neu
(.)_neu: New value of quantity (.) computed during an integration step
Definition: sico_variables_m.F90:67
sico_variables_m::sigma_c
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) sigma_c
sigma_c(kc,j,i): Effective stress in the upper (kc) ice domain
Definition: sico_variables_m.F90:426
sico_variables_m::zb
real(dp), dimension(0:jmax, 0:imax) zb
zb(j,i): Coordinate z of the ice base
Definition: sico_variables_m.F90:174
sico_variables_m::omega_c_neu
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) omega_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
Definition: sico_variables_m.F90:483
calc_temp_m::calc_temp_poly
subroutine, public calc_temp_poly(dxi, deta, dzeta_c, dzeta_t, dzeta_r, dtime_temp)
Computation of temperature, water content and age in polythermal mode.
Definition: calc_temp_m.F90:53
sico_variables_m::am_perp_st
real(dp), dimension(0:jmax, 0:imax) am_perp_st
am_perp_st(j,i): Steady-state part of am_perp (without contribution of dzm_dtau)
Definition: sico_variables_m.F90:339
sico_variables_m::dzb_dtau
real(dp), dimension(0:jmax, 0:imax) dzb_dtau
dzb_dtau(j,i): Derivative of zb by tau (time)
Definition: sico_variables_m.F90:238
sico_variables_m::zeta_c
real(dp), dimension(0:kcmax) zeta_c
zeta_c(kc): Sigma coordinate zeta_c of grid point kc
Definition: sico_variables_m.F90:119
sico_variables_m::h_t
real(dp), dimension(0:jmax, 0:imax) h_t
H_t(j,i): Thickness of ice in the lower (kt) domain (thickness of the temperate layer for POLY...
Definition: sico_variables_m.F90:192
sico_variables_m
Declarations of global variables for SICOPOLIS.
Definition: sico_variables_m.F90:35
sico_variables_m::dzs_dtau
real(dp), dimension(0:jmax, 0:imax) dzs_dtau
dzs_dtau(j,i): Derivative of zs by tau (time)
Definition: sico_variables_m.F90:234
sico_variables_m::vx_t
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) vx_t
vx_t(kt,j,i): Velocity in x-direction in the lower (kt) ice domain (at (i+1/2,j,kt)) ...
Definition: sico_variables_m.F90:443
sico_variables_m::dzb_dxi_g
real(dp), dimension(0:jmax, 0:imax) dzb_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
Definition: sico_variables_m.F90:218
Generated on Thu Aug 3 2017 21:26:01 for SICOPOLIS V3.3 by   doxygen 1.8.11