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