SICOPOLIS V5-dev  Revision 1368
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-2018 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 #ifndef 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 #ifndef 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 #ifndef 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  enh_c(kc,j,i), 0_i1b) &
1523  *de_c(kc,j,i)**2
1524  end if
1525 #endif
1526 
1527  ci1(kc) = ai1(kc)/h_c(j,i)
1528 
1529 end do
1530 
1531 #if (ADV_VERT==1)
1532 
1533 kc=0
1534 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
1535 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
1536 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
1537 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
1538 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! only needed for kc=0 ...
1539 kc=kcmax-1
1540 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
1541 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
1542 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
1543 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
1544 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! ... and kc=KCMAX-1
1545 
1546 #elif (ADV_VERT==2 || ADV_VERT==3)
1547 
1548 do kc=0, kcmax-1
1549  ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
1550  ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
1551  ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
1552  adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
1553  abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
1554 end do
1555 
1556 #endif
1557 
1558 do kc=0, kcmax-1
1559  temp_c_help(kc) = 0.5_dp*(temp_c(kc,j,i)+temp_c(kc+1,j,i))
1560  ct6(kc) = at6(kc) &
1561  *kappa_val(temp_c_help(kc)) &
1562  /h_c(j,i)
1563  ci2(kc) = ai2(kc)/h_c(j,i)
1564 end do
1565 
1566 #if (ADV_HOR==3)
1567 dtt_dxi = 2.0_dp*dtt_2dxi
1568 dtt_deta = 2.0_dp*dtt_2deta
1569 #endif
1570 
1571 !-------- Set up the temperature equations (ice and bedrock
1572 ! simultaneously) --------
1573 
1574 kr=0
1575 lgs_a1(kr) = 1.0_dp
1576 lgs_a2(kr) = -1.0_dp
1577 lgs_b(kr) = clb1
1578 
1579 #if (Q_LITHO==1)
1580 ! (coupled heat-conducting bedrock)
1581 
1582 do kr=1, krmax-1
1583  lgs_a0(kr) = - ctr1
1584  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
1585  lgs_a2(kr) = - ctr1
1586  lgs_b(kr) = temp_r(kr,j,i)
1587 end do
1588 
1589 #elif (Q_LITHO==0)
1590 ! (no coupled heat-conducting bedrock)
1591 
1592 do kr=1, krmax-1
1593  lgs_a0(kr) = 1.0_dp
1594  lgs_a1(kr) = 0.0_dp
1595  lgs_a2(kr) = -1.0_dp
1596  lgs_b(kr) = 2.0_dp*clb1
1597 end do
1598 
1599 #endif
1600 
1601 kr=krmax
1602 kc=0
1603 lgs_a0(kr) = ccb2
1604 lgs_a1(kr) = -(ccb1+ccb2)
1605 lgs_a2(kr) = ccb1
1606 lgs_b(kr) = ccb3+ccb4
1607 
1608 do kc=1, kcmax-1
1609 
1610 #if (ADV_VERT==1)
1611 
1612  lgs_a0(krmax+kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
1613  -ct5(kc)*ct6(kc-1)
1614  lgs_a1(krmax+kc) = 1.0_dp+ct5(kc)*(ct6(kc)+ct6(kc-1))
1615  lgs_a2(krmax+kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
1616  -ct5(kc)*ct6(kc)
1617 
1618 #elif (ADV_VERT==2)
1619 
1620  lgs_a0(krmax+kc) &
1621  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
1622  -ct5(kc)*ct6(kc-1)
1623  lgs_a1(krmax+kc) &
1624  = 1.0_dp &
1625  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
1626  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
1627  +ct5(kc)*(ct6(kc)+ct6(kc-1))
1628  lgs_a2(krmax+kc) &
1629  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
1630  -ct5(kc)*ct6(kc)
1631 
1632 #elif (ADV_VERT==3)
1633 
1634  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
1635 
1636  lgs_a0(krmax+kc) &
1637  = -max(adv_vert_help, 0.0_dp) &
1638  -ct5(kc)*ct6(kc-1)
1639  lgs_a1(krmax+kc) &
1640  = 1.0_dp &
1641  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
1642  +ct5(kc)*(ct6(kc)+ct6(kc-1))
1643  lgs_a2(krmax+kc) &
1644  = min(adv_vert_help, 0.0_dp) &
1645  -ct5(kc)*ct6(kc)
1646 
1647 #endif
1648 
1649 #if (ADV_HOR==2)
1650 
1651  lgs_b(krmax+kc) = temp_c(kc,j,i) + ct7(kc) &
1652  -dtt_2dxi* &
1653  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
1654  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
1655  *insq_g11_sgx(j,i) &
1656  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
1657  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
1658  *insq_g11_sgx(j,i-1) ) &
1659  -dtt_2deta* &
1660  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
1661  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
1662  *insq_g22_sgy(j,i) &
1663  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
1664  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
1665  *insq_g22_sgy(j-1,i) )
1666 
1667 #elif (ADV_HOR==3)
1668 
1669  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
1670  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
1671 
1672  lgs_b(krmax+kc) = temp_c(kc,j,i) + ct7(kc) &
1673  -dtt_dxi* &
1674  ( min(vx_c_help, 0.0_dp) &
1675  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
1676  *insq_g11_sgx(j,i) &
1677  +max(vx_c_help, 0.0_dp) &
1678  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
1679  *insq_g11_sgx(j,i-1) ) &
1680  -dtt_deta* &
1681  ( min(vy_c_help, 0.0_dp) &
1682  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
1683  *insq_g22_sgy(j,i) &
1684  +max(vy_c_help, 0.0_dp) &
1685  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
1686  *insq_g22_sgy(j-1,i) )
1687 
1688 #endif
1689 
1690 end do
1691 
1692 kc=kcmax
1693 lgs_a0(krmax+kc) = 0.0_dp
1694 lgs_a1(krmax+kc) = 1.0_dp
1695 lgs_b(krmax+kc) = temp_s(j,i)
1696 
1697 !-------- Solve system of linear equations --------
1698 
1699 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax+krmax)
1700 
1701 !-------- Assign the result --------
1702 
1703 do kr=0, krmax
1704  temp_r_neu(kr,j,i) = lgs_x(kr)
1705 end do
1706 
1707 do kc=0, kcmax
1708  temp_c_neu(kc,j,i) = lgs_x(krmax+kc)
1709 end do
1710 
1711 !-------- Set water content in the non-existing temperate layer
1712 ! to zero --------
1713 
1714 do kt=0, ktmax
1715  omega_t_neu(kt,j,i) = 0.0_dp
1716 end do
1717 
1718 !-------- Water drainage from the non-existing temperate layer --------
1719 
1720 q_tld(j,i) = 0.0_dp
1721 
1722 !-------- Set up the equations for the age of cold ice --------
1723 
1724 kc=0 ! adv_vert_sg(0) <= 0
1725 lgs_a1(kc) = 1.0_dp - min(adv_vert_sg(kc), 0.0_dp) ! (directed downward)
1726 lgs_a2(kc) = min(adv_vert_sg(kc), 0.0_dp) ! assumed/enforced
1727 
1728 #if (ADV_HOR==2)
1729 
1730 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1731  -dtt_2dxi* &
1732  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
1733  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1734  *insq_g11_sgx(j,i) &
1735  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
1736  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1737  *insq_g11_sgx(j,i-1) ) &
1738  -dtt_2deta* &
1739  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
1740  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1741  *insq_g22_sgy(j,i) &
1742  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
1743  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1744  *insq_g22_sgy(j-1,i) )
1745 
1746 #elif (ADV_HOR==3)
1747 
1748 vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
1749 vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
1750 
1751 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1752  -dtt_dxi* &
1753  ( min(vx_c_help, 0.0_dp) &
1754  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1755  *insq_g11_sgx(j,i) &
1756  +max(vx_c_help, 0.0_dp) &
1757  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1758  *insq_g11_sgx(j,i-1) ) &
1759  -dtt_deta* &
1760  ( min(vy_c_help, 0.0_dp) &
1761  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1762  *insq_g22_sgy(j,i) &
1763  +max(vy_c_help, 0.0_dp) &
1764  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1765  *insq_g22_sgy(j-1,i) )
1766 
1767 #endif
1768 
1769 do kc=1, kcmax-1
1770 
1771 #if (ADV_VERT==1)
1772 
1773  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
1774  -ci1(kc)*ci2(kc-1)
1775  lgs_a1(kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
1776  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
1777  -ci1(kc)*ci2(kc)
1778 
1779 #elif (ADV_VERT==2)
1780 
1781  lgs_a0(kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
1782  lgs_a1(kc) = 1.0_dp &
1783  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
1784  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
1785  lgs_a2(kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
1786 
1787 #elif (ADV_VERT==3)
1788 
1789  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
1790 
1791  lgs_a0(kc) = -max(adv_vert_help, 0.0_dp)
1792  lgs_a1(kc) = 1.0_dp &
1793  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp)
1794  lgs_a2(kc) = min(adv_vert_help, 0.0_dp)
1795 
1796 #endif
1797 
1798 #if (ADV_HOR==2)
1799 
1800  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1801  -dtt_2dxi* &
1802  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
1803  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1804  *insq_g11_sgx(j,i) &
1805  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
1806  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1807  *insq_g11_sgx(j,i-1) ) &
1808  -dtt_2deta* &
1809  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
1810  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1811  *insq_g22_sgy(j,i) &
1812  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
1813  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1814  *insq_g22_sgy(j-1,i) )
1815 
1816 #elif (ADV_HOR==3)
1817 
1818  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
1819  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
1820 
1821  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1822  -dtt_dxi* &
1823  ( min(vx_c_help, 0.0_dp) &
1824  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1825  *insq_g11_sgx(j,i) &
1826  +max(vx_c_help, 0.0_dp) &
1827  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1828  *insq_g11_sgx(j,i-1) ) &
1829  -dtt_deta* &
1830  ( min(vy_c_help, 0.0_dp) &
1831  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1832  *insq_g22_sgy(j,i) &
1833  +max(vy_c_help, 0.0_dp) &
1834  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1835  *insq_g22_sgy(j-1,i) )
1836 
1837 #endif
1838 
1839 end do
1840 
1841 kc=kcmax
1842 if (as_perp(j,i) >= zero) then
1843  lgs_a0(kc) = 0.0_dp
1844  lgs_a1(kc) = 1.0_dp
1845  lgs_b(kc) = 0.0_dp
1846 else
1847  lgs_a0(kc) = -max(adv_vert_sg(kc-1), 0.0_dp)
1848  lgs_a1(kc) = 1.0_dp + max(adv_vert_sg(kc-1), 0.0_dp)
1849  ! adv_vert_sg(KCMAX-1) >= 0 (directed upward)
1850  ! assumed/enforced
1851 #if (ADV_HOR==2)
1852 
1853  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1854  -dtt_2dxi* &
1855  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
1856  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1857  *insq_g11_sgx(j,i) &
1858  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
1859  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1860  *insq_g11_sgx(j,i-1) ) &
1861  -dtt_2deta* &
1862  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
1863  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1864  *insq_g22_sgy(j,i) &
1865  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
1866  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1867  *insq_g22_sgy(j-1,i) )
1868 
1869 #elif (ADV_HOR==3)
1870 
1871  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
1872  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
1873 
1874  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1875  -dtt_dxi* &
1876  ( min(vx_c_help, 0.0_dp) &
1877  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1878  *insq_g11_sgx(j,i) &
1879  +max(vx_c_help, 0.0_dp) &
1880  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1881  *insq_g11_sgx(j,i-1) ) &
1882  -dtt_deta* &
1883  ( min(vy_c_help, 0.0_dp) &
1884  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1885  *insq_g22_sgy(j,i) &
1886  +max(vy_c_help, 0.0_dp) &
1887  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1888  *insq_g22_sgy(j-1,i) )
1889 
1890 #endif
1891 
1892 end if
1893 
1894 !-------- Solve system of linear equations --------
1895 
1896 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
1897 
1898 !-------- Assign the result,
1899 ! restriction to interval [0, AGE_MAX yr] --------
1900 
1901 do kc=0, kcmax
1902 
1903  age_c_neu(kc,j,i) = lgs_x(kc)
1904 
1905  if (age_c_neu(kc,j,i) < (age_min*year_sec)) &
1906  age_c_neu(kc,j,i) = 0.0_dp
1907  if (age_c_neu(kc,j,i) > (age_max*year_sec)) &
1908  age_c_neu(kc,j,i) = age_max*year_sec
1909 
1910 end do
1911 
1912 !-------- Age of the ice in the non-existing temperate layer --------
1913 
1914 do kt=0, ktmax
1915  age_t_neu(kt,j,i) = age_c_neu(0,j,i)
1916 end do
1917 
1918 end subroutine calc_temp1
1919 
1920 !-------------------------------------------------------------------------------
1921 !> Computation of temperature and age for an ice column with a temperate base
1922 !! overlain by cold ice.
1923 !<------------------------------------------------------------------------------
1924 subroutine calc_temp2(at1, at2_1, at2_2, at3_1, at3_2, &
1925  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
1926  ai1, ai2, &
1927  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
1928 
1929 use ice_material_properties_m, only : ratefac_c, kappa_val, c_val, &
1930  creep, viscosity
1931 #ifndef ALLOW_OPENAD /* Normal */
1932 use sico_maths_m, only : tri_sle
1933 #else /* OpenAD */
1934 use sico_maths_m
1935 #endif /* Normal vs. OpenAD */
1936 
1937 implicit none
1938 
1939 #ifndef ALLOW_OPENAD /* Normal */
1940 integer(i4b), intent(in) :: i, j
1941 #else /* OpenAD */
1942 integer(i4b), intent(inout) :: i, j
1943 #endif /* Normal vs. OpenAD */
1944 real(dp), intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
1945  at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
1946  at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
1947  ai1(0:kcmax), ai2(0:kcmax), &
1948  atr1, alb1
1949 real(dp), intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
1950 
1951 integer(i4b) :: kc, kt, kr
1952 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
1953  ct5(0:kcmax), ct6(0:kcmax), ct7(0:kcmax), ctr1, clb1
1954 real(dp) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), ct3_sg(0:kcmax), &
1955  ct4_sg(0:kcmax), adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
1956 real(dp) :: ci1(0:kcmax), ci2(0:kcmax)
1957 real(dp) :: temp_c_help(0:kcmax)
1958 real(dp) :: vx_c_help, vy_c_help
1959 real(dp) :: adv_vert_help
1960 real(dp) :: dtt_dxi, dtt_deta
1961 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
1962  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
1963  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
1964  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
1965  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
1966 real(dp), parameter :: zero=0.0_dp
1967 
1968 !-------- Check for boundary points --------
1969 
1970 if ((i == 0).or.(i == imax).or.(j == 0).or.(j == jmax)) then
1971  errormsg = ' >>> calc_temp2: Boundary points not allowed!'
1972  call error(errormsg)
1973 end if
1974 
1975 !-------- Abbreviations --------
1976 
1977 ctr1 = atr1
1978 clb1 = alb1*q_geo(j,i)
1979 
1980 #if (ADV_VERT==1)
1981 
1982 do kc=1, kcmax-1
1983  ct1(kc) = at1(kc)/h_c(j,i)*0.5_dp*(vz_c(kc,j,i)+vz_c(kc-1,j,i))
1984 end do
1985 
1986 kc=0
1987 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
1988  ! only needed for kc=0 ...
1989 kc=kcmax-1
1990 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
1991  ! ... and kc=KCMAX-1
1992 
1993 #elif (ADV_VERT==2 || ADV_VERT==3)
1994 
1995 do kc=0, kcmax-1
1996  ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
1997 end do
1998 
1999 #endif
2000 
2001 do kc=0, kcmax
2002 
2003  ct2(kc) = ( at2_1(kc)*dzm_dtau(j,i) &
2004  +at2_2(kc)*dh_c_dtau(j,i) )/h_c(j,i)
2005  ct3(kc) = ( at3_1(kc)*dzm_dxi_g(j,i) &
2006  +at3_2(kc)*dh_c_dxi_g(j,i) )/h_c(j,i) &
2007  *0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1)) *insq_g11_g(j,i)
2008  ct4(kc) = ( at4_1(kc)*dzm_deta_g(j,i) &
2009  +at4_2(kc)*dh_c_deta_g(j,i) )/h_c(j,i) &
2010  *0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i)) *insq_g22_g(j,i)
2011  ct5(kc) = at5(kc) &
2012  /c_val(temp_c(kc,j,i)) &
2013  /h_c(j,i)
2014 
2015 #if (DYNAMICS==2)
2016  if (.not.flag_shelfy_stream(j,i)) then
2017 #endif
2018  ct7(kc) = at7 &
2019  /c_val(temp_c(kc,j,i)) &
2020  *enh_c(kc,j,i) &
2021  *ratefac_c(temp_c(kc,j,i), temp_c_m(kc,j,i)) &
2022  *creep(sigma_c(kc,j,i)) &
2023  *sigma_c(kc,j,i)*sigma_c(kc,j,i)
2024 #if (DYNAMICS==2)
2025  else
2026  ct7(kc) = 2.0_dp*at7 &
2027  /c_val(temp_c(kc,j,i)) &
2028  *viscosity(de_c(kc,j,i), &
2029  temp_c(kc,j,i), temp_c_m(kc,j,i), 0.0_dp, &
2030  enh_c(kc,j,i), 0_i1b) &
2031  *de_c(kc,j,i)**2
2032  end if
2033 #endif
2034 
2035  ci1(kc) = ai1(kc)/h_c(j,i)
2036 
2037 end do
2038 
2039 #if (ADV_VERT==1)
2040 
2041 kc=0
2042 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2043 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2044 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2045 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2046 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! only needed for kc=0 ...
2047 kc=kcmax-1
2048 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2049 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2050 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2051 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2052 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! ... and kc=KCMAX-1
2053 
2054 #elif (ADV_VERT==2 || ADV_VERT==3)
2055 
2056 do kc=0, kcmax-1
2057  ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2058  ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2059  ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2060  adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2061  abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
2062 end do
2063 
2064 #endif
2065 
2066 do kc=0, kcmax-1
2067  temp_c_help(kc) = 0.5_dp*(temp_c(kc,j,i)+temp_c(kc+1,j,i))
2068  ct6(kc) = at6(kc) &
2069  *kappa_val(temp_c_help(kc)) &
2070  /h_c(j,i)
2071  ci2(kc) = ai2(kc)/h_c(j,i)
2072 end do
2073 
2074 #if (ADV_HOR==3)
2075 dtt_dxi = 2.0_dp*dtt_2dxi
2076 dtt_deta = 2.0_dp*dtt_2deta
2077 #endif
2078 
2079 !-------- Set up the equations for the bedrock temperature --------
2080 
2081 kr=0
2082 lgs_a1(kr) = 1.0_dp
2083 lgs_a2(kr) = -1.0_dp
2084 lgs_b(kr) = clb1
2085 
2086 #if (Q_LITHO==1)
2087 ! (coupled heat-conducting bedrock)
2088 
2089 do kr=1, krmax-1
2090  lgs_a0(kr) = - ctr1
2091  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
2092  lgs_a2(kr) = - ctr1
2093  lgs_b(kr) = temp_r(kr,j,i)
2094 end do
2095 
2096 #elif (Q_LITHO==0)
2097 ! (no coupled heat-conducting bedrock)
2098 
2099 do kr=1, krmax-1
2100  lgs_a0(kr) = 1.0_dp
2101  lgs_a1(kr) = 0.0_dp
2102  lgs_a2(kr) = -1.0_dp
2103  lgs_b(kr) = 2.0_dp*clb1
2104 end do
2105 
2106 #endif
2107 
2108 kr=krmax
2109 lgs_a0(kr) = 0.0_dp
2110 lgs_a1(kr) = 1.0_dp
2111 lgs_b(kr) = temp_t_m(0,j,i)
2112 
2113 !-------- Solve system of linear equations --------
2114 
2115 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
2116 
2117 !-------- Assign the result --------
2118 
2119 do kr=0, krmax
2120  temp_r_neu(kr,j,i) = lgs_x(kr)
2121 end do
2122 
2123 !-------- Set up the equations for the ice temperature --------
2124 
2125 kc=0
2126 lgs_a1(kc) = 1.0_dp
2127 lgs_a2(kc) = 0.0_dp
2128 lgs_b(kc) = temp_c_m(0,j,i)
2129 
2130 do kc=1, kcmax-1
2131 
2132 #if (ADV_VERT==1)
2133 
2134  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2135  -ct5(kc)*ct6(kc-1)
2136  lgs_a1(kc) = 1.0_dp+ct5(kc)*(ct6(kc)+ct6(kc-1))
2137  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2138  -ct5(kc)*ct6(kc)
2139 
2140 #elif (ADV_VERT==2)
2141 
2142  lgs_a0(kc) &
2143  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2144  -ct5(kc)*ct6(kc-1)
2145  lgs_a1(kc) &
2146  = 1.0_dp &
2147  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2148  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
2149  +ct5(kc)*(ct6(kc)+ct6(kc-1))
2150  lgs_a2(kc) &
2151  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
2152  -ct5(kc)*ct6(kc)
2153 
2154 #elif (ADV_VERT==3)
2155 
2156  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
2157 
2158  lgs_a0(kc) &
2159  = -max(adv_vert_help, 0.0_dp) &
2160  -ct5(kc)*ct6(kc-1)
2161  lgs_a1(kc) &
2162  = 1.0_dp &
2163  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
2164  +ct5(kc)*(ct6(kc)+ct6(kc-1))
2165  lgs_a2(kc) &
2166  = min(adv_vert_help, 0.0_dp) &
2167  -ct5(kc)*ct6(kc)
2168 
2169 #endif
2170 
2171 #if (ADV_HOR==2)
2172 
2173  lgs_b(kc) = temp_c(kc,j,i) + ct7(kc) &
2174  -dtt_2dxi* &
2175  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2176  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
2177  *insq_g11_sgx(j,i) &
2178  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2179  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
2180  *insq_g11_sgx(j,i-1) ) &
2181  -dtt_2deta* &
2182  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2183  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
2184  *insq_g22_sgy(j,i) &
2185  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2186  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
2187  *insq_g22_sgy(j-1,i) )
2188 
2189 #elif (ADV_HOR==3)
2190 
2191  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2192  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2193 
2194  lgs_b(kc) = temp_c(kc,j,i) + ct7(kc) &
2195  -dtt_dxi* &
2196  ( min(vx_c_help, 0.0_dp) &
2197  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
2198  *insq_g11_sgx(j,i) &
2199  +max(vx_c_help, 0.0_dp) &
2200  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
2201  *insq_g11_sgx(j,i-1) ) &
2202  -dtt_deta* &
2203  ( min(vy_c_help, 0.0_dp) &
2204  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
2205  *insq_g22_sgy(j,i) &
2206  +max(vy_c_help, 0.0_dp) &
2207  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
2208  *insq_g22_sgy(j-1,i) )
2209 
2210 #endif
2211 
2212 end do
2213 
2214 kc=kcmax
2215 lgs_a0(kc) = 0.0_dp
2216 lgs_a1(kc) = 1.0_dp
2217 lgs_b(kc) = temp_s(j,i)
2218 
2219 !-------- Solve system of linear equations --------
2220 
2221 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
2222 
2223 !-------- Assign the result --------
2224 
2225 do kc=0, kcmax
2226  temp_c_neu(kc,j,i) = lgs_x(kc)
2227 end do
2228 
2229 !-------- Set water content in the non-existing temperate layer
2230 ! to zero --------
2231 
2232 do kt=0, ktmax
2233  omega_t_neu(kt,j,i) = 0.0_dp
2234 end do
2235 
2236 !-------- Water drainage from the non-existing temperate layer --------
2237 
2238 q_tld(j,i) = 0.0_dp
2239 
2240 !-------- Set up the equations for the age of cold ice --------
2241 
2242 kc=0 ! adv_vert_sg(0) <= 0
2243 lgs_a1(kc) = 1.0_dp - min(adv_vert_sg(kc), 0.0_dp) ! (directed downward)
2244 lgs_a2(kc) = min(adv_vert_sg(kc), 0.0_dp) ! assumed/enforced
2245 
2246 #if (ADV_HOR==2)
2247 
2248 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2249  -dtt_2dxi* &
2250  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2251  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2252  *insq_g11_sgx(j,i) &
2253  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2254  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2255  *insq_g11_sgx(j,i-1) ) &
2256  -dtt_2deta* &
2257  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2258  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2259  *insq_g22_sgy(j,i) &
2260  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2261  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2262  *insq_g22_sgy(j-1,i) )
2263 
2264 #elif (ADV_HOR==3)
2265 
2266 vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2267 vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2268 
2269 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2270  -dtt_dxi* &
2271  ( min(vx_c_help, 0.0_dp) &
2272  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2273  *insq_g11_sgx(j,i) &
2274  +max(vx_c_help, 0.0_dp) &
2275  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2276  *insq_g11_sgx(j,i-1) ) &
2277  -dtt_deta* &
2278  ( min(vy_c_help, 0.0_dp) &
2279  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2280  *insq_g22_sgy(j,i) &
2281  +max(vy_c_help, 0.0_dp) &
2282  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2283  *insq_g22_sgy(j-1,i) )
2284 
2285 #endif
2286 
2287 do kc=1, kcmax-1
2288 
2289 #if (ADV_VERT==1)
2290 
2291  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2292  -ci1(kc)*ci2(kc-1)
2293  lgs_a1(kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
2294  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2295  -ci1(kc)*ci2(kc)
2296 
2297 #elif (ADV_VERT==2)
2298 
2299  lgs_a0(kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
2300  lgs_a1(kc) = 1.0_dp &
2301  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2302  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
2303  lgs_a2(kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
2304 
2305 #elif (ADV_VERT==3)
2306 
2307  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
2308 
2309  lgs_a0(kc) = -max(adv_vert_help, 0.0_dp)
2310  lgs_a1(kc) = 1.0_dp &
2311  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp)
2312  lgs_a2(kc) = min(adv_vert_help, 0.0_dp)
2313 
2314 #endif
2315 
2316 #if (ADV_HOR==2)
2317 
2318  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2319  -dtt_2dxi* &
2320  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2321  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2322  *insq_g11_sgx(j,i) &
2323  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2324  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2325  *insq_g11_sgx(j,i-1) ) &
2326  -dtt_2deta* &
2327  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2328  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2329  *insq_g22_sgy(j,i) &
2330  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2331  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2332  *insq_g22_sgy(j-1,i) )
2333 
2334 #elif (ADV_HOR==3)
2335 
2336  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2337  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2338 
2339  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2340  -dtt_dxi* &
2341  ( min(vx_c_help, 0.0_dp) &
2342  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2343  *insq_g11_sgx(j,i) &
2344  +max(vx_c_help, 0.0_dp) &
2345  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2346  *insq_g11_sgx(j,i-1) ) &
2347  -dtt_deta* &
2348  ( min(vy_c_help, 0.0_dp) &
2349  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2350  *insq_g22_sgy(j,i) &
2351  +max(vy_c_help, 0.0_dp) &
2352  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2353  *insq_g22_sgy(j-1,i) )
2354 
2355 #endif
2356 
2357 end do
2358 
2359 kc=kcmax
2360 if (as_perp(j,i) >= zero) then
2361  lgs_a0(kc) = 0.0_dp
2362  lgs_a1(kc) = 1.0_dp
2363  lgs_b(kc) = 0.0_dp
2364 else
2365  lgs_a0(kc) = -max(adv_vert_sg(kc-1), 0.0_dp)
2366  lgs_a1(kc) = 1.0_dp + max(adv_vert_sg(kc-1), 0.0_dp)
2367  ! adv_vert_sg(KCMAX-1) >= 0 (directed upward)
2368  ! assumed/enforced
2369 #if (ADV_HOR==2)
2370 
2371  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2372  -dtt_2dxi* &
2373  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2374  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2375  *insq_g11_sgx(j,i) &
2376  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2377  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2378  *insq_g11_sgx(j,i-1) ) &
2379  -dtt_2deta* &
2380  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2381  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2382  *insq_g22_sgy(j,i) &
2383  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2384  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2385  *insq_g22_sgy(j-1,i) )
2386 
2387 #elif (ADV_HOR==3)
2388 
2389  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2390  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2391 
2392  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2393  -dtt_dxi* &
2394  ( min(vx_c_help, 0.0_dp) &
2395  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2396  *insq_g11_sgx(j,i) &
2397  +max(vx_c_help, 0.0_dp) &
2398  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2399  *insq_g11_sgx(j,i-1) ) &
2400  -dtt_deta* &
2401  ( min(vy_c_help, 0.0_dp) &
2402  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2403  *insq_g22_sgy(j,i) &
2404  +max(vy_c_help, 0.0_dp) &
2405  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2406  *insq_g22_sgy(j-1,i) )
2407 
2408 #endif
2409 
2410 end if
2411 
2412 !-------- Solve system of linear equations --------
2413 
2414 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
2415 
2416 !-------- Assign the result,
2417 ! restriction to interval [0, AGE_MAX yr] --------
2418 
2419 do kc=0, kcmax
2420 
2421  age_c_neu(kc,j,i) = lgs_x(kc)
2422 
2423  if (age_c_neu(kc,j,i) < (age_min*year_sec)) &
2424  age_c_neu(kc,j,i) = 0.0_dp
2425  if (age_c_neu(kc,j,i) > (age_max*year_sec)) &
2426  age_c_neu(kc,j,i) = age_max*year_sec
2427 
2428 end do
2429 
2430 !-------- Age of the ice in the non-existing temperate layer --------
2431 
2432 do kt=0, ktmax
2433  age_t_neu(kt,j,i) = age_c_neu(0,j,i)
2434 end do
2435 
2436 end subroutine calc_temp2
2437 
2438 !-------------------------------------------------------------------------------
2439 !> Computation of temperature, water content and age for an ice column with a
2440 !! temperate base overlain by a temperate-ice layer.
2441 !<------------------------------------------------------------------------------
2442 subroutine calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
2443  at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
2444  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
2445  ai1, ai2, ai3, dzeta_t, &
2446  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
2447 
2448 use ice_material_properties_m, only : ratefac_c, ratefac_t, kappa_val, c_val, &
2449  creep, viscosity
2450 #ifndef ALLOW_OPENAD /* Normal */
2451 use sico_maths_m, only : tri_sle
2452 #else /* OpenAD */
2453 use sico_maths_m
2454 #endif /* Normal vs. OpenAD */
2455 
2456 implicit none
2457 
2458 #ifndef ALLOW_OPENAD /* Normal */
2459 integer(i4b), intent(in) :: i, j
2460 #else /* OpenAD */
2461 integer(i4b), intent(inout) :: i, j
2462 #endif /* Normal vs. OpenAD */
2463 real(dp), intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
2464  at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
2465  at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
2466  ai1(0:kcmax), ai2(0:kcmax), ai3, &
2467  atr1, am1, am2, alb1
2468 real(dp), intent(in) :: aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld
2469 real(dp), intent(in) :: dzeta_t, dtime_temp, dtt_2dxi, dtt_2deta
2470 
2471 integer(i4b) :: kc, kt, kr
2472 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
2473  ct5(0:kcmax), ct6(0:kcmax), ct7(0:kcmax), ctr1, cm1, cm2, &
2474  clb1
2475 real(dp) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), ct3_sg(0:kcmax), &
2476  ct4_sg(0:kcmax), adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
2477 real(dp) :: ci1(0:kcmax), ci2(0:kcmax), ci3
2478 real(dp) :: cw1(0:ktmax), cw2(0:ktmax), cw3(0:ktmax), cw4(0:ktmax), &
2479  cw5, cw7(0:ktmax), cw8, cw9(0:ktmax)
2480 real(dp) :: cw1_sg(0:ktmax), cw2_sg(0:ktmax), cw3_sg(0:ktmax), &
2481  cw4_sg(0:ktmax), adv_vert_w_sg(0:ktmax), abs_adv_vert_w_sg(0:ktmax)
2482 real(dp) :: sigma_c_help(0:kcmax), sigma_t_help(0:ktmax), &
2483  temp_c_help(0:kcmax)
2484 real(dp) :: vx_c_help, vy_c_help, vx_t_help, vy_t_help
2485 real(dp) :: adv_vert_help, adv_vert_w_help
2486 real(dp) :: dtt_dxi, dtt_deta
2487 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
2488  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
2489  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
2490  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
2491  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
2492 real(dp), parameter :: zero=0.0_dp
2493 
2494 !-------- Check for boundary points --------
2495 
2496 if ((i == 0).or.(i == imax).or.(j == 0).or.(j == jmax)) then
2497  errormsg = ' >>> calc_temp3: Boundary points not allowed!'
2498  call error(errormsg)
2499 end if
2500 
2501 !-------- Abbreviations --------
2502 
2503 ctr1 = atr1
2504 cm1 = am1*h_c_neu(j,i)
2505 clb1 = alb1*q_geo(j,i)
2506 
2507 #if (ADV_VERT==1)
2508 
2509 do kc=1, kcmax-1
2510  ct1(kc) = at1(kc)/h_c_neu(j,i)*0.5_dp*(vz_c(kc,j,i)+vz_c(kc-1,j,i))
2511 end do
2512 
2513 kc=0
2514 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c_neu(j,i)*vz_c(kc,j,i)
2515  ! only needed for kc=0 ...
2516 kc=kcmax-1
2517 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c_neu(j,i)*vz_c(kc,j,i)
2518  ! ... and kc=KCMAX-1
2519 
2520 #elif (ADV_VERT==2 || ADV_VERT==3)
2521 
2522 do kc=0, kcmax-1
2523  ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c_neu(j,i)*vz_c(kc,j,i)
2524 end do
2525 
2526 #endif
2527 
2528 do kc=0, kcmax
2529 
2530  ct2(kc) = ( at2_1(kc)*dzm_dtau(j,i) &
2531  +at2_2(kc)*dh_c_dtau(j,i) )/h_c_neu(j,i)
2532  ct3(kc) = ( at3_1(kc)*dzm_dxi_g(j,i) &
2533  +at3_2(kc)*dh_c_dxi_g(j,i) )/h_c_neu(j,i) &
2534  *0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1)) *insq_g11_g(j,i)
2535  ct4(kc) = ( at4_1(kc)*dzm_deta_g(j,i) &
2536  +at4_2(kc)*dh_c_deta_g(j,i) )/h_c_neu(j,i) &
2537  *0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i)) *insq_g22_g(j,i)
2538  ct5(kc) = at5(kc) &
2539  /c_val(temp_c(kc,j,i)) &
2540  /h_c_neu(j,i)
2541 
2542  sigma_c_help(kc) &
2543  = rho*g*h_c_neu(j,i)*(1.0_dp-eaz_c_quotient(kc)) &
2544  *(dzs_dxi_g(j,i)**2+dzs_deta_g(j,i)**2)**0.5_dp
2545 
2546 #if (DYNAMICS==2)
2547  if (.not.flag_shelfy_stream(j,i)) then
2548 #endif
2549  ct7(kc) = at7 &
2550  /c_val(temp_c(kc,j,i)) &
2551  *enh_c(kc,j,i) &
2552  *ratefac_c(temp_c(kc,j,i), temp_c_m(kc,j,i)) &
2553  *creep(sigma_c_help(kc)) &
2554  *sigma_c_help(kc)*sigma_c_help(kc)
2555 #if (DYNAMICS==2)
2556  else
2557  ct7(kc) = 2.0_dp*at7 &
2558  /c_val(temp_c(kc,j,i)) &
2559  *viscosity(de_c(kc,j,i), &
2560  temp_c(kc,j,i), temp_c_m(kc,j,i), 0.0_dp, &
2561  enh_c(kc,j,i), 0_i1b) &
2562  *de_c(kc,j,i)**2
2563  end if
2564 #endif
2565 
2566  ci1(kc) = ai1(kc)/h_c_neu(j,i)
2567 
2568 end do
2569 
2570 #if (ADV_VERT==1)
2571 
2572 kc=0
2573 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2574 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2575 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2576 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2577 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! only needed for kc=0 ...
2578 kc=kcmax-1
2579 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2580 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2581 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2582 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2583 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! ... and kc=KCMAX-1
2584 
2585 #elif (ADV_VERT==2 || ADV_VERT==3)
2586 
2587 do kc=0, kcmax-1
2588  ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2589  ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2590  ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2591  adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2592  abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
2593 end do
2594 
2595 #endif
2596 
2597 do kc=0, kcmax-1
2598  temp_c_help(kc) = 0.5_dp*(temp_c(kc,j,i)+temp_c(kc+1,j,i))
2599  ct6(kc) = at6(kc) &
2600  *kappa_val(temp_c_help(kc)) &
2601  /h_c_neu(j,i)
2602  ci2(kc) = ai2(kc)/h_c_neu(j,i)
2603 end do
2604 
2605 cw5 = aw5/(h_t_neu(j,i)**2)
2606 cw8 = aw8
2607 ci3 = ai3/(h_t_neu(j,i)**2)
2608 
2609 #if (ADV_VERT==1)
2610 
2611 do kt=1, ktmax-1
2612  cw1(kt) = aw1/h_t_neu(j,i)*0.5_dp*(vz_t(kt,j,i)+vz_t(kt-1,j,i))
2613 end do
2614 
2615 kt=ktmax
2616 cw1(kt) = aw1/h_t_neu(j,i)*0.5_dp*(vz_t(kt-1,j,i)+vz_c(0,j,i))
2617 
2618 kt=0
2619 cw1_sg(kt) = aw1/h_t_neu(j,i)*vz_t(kt,j,i)
2620  ! only needed for kt=0 ...
2621 kt=ktmax-1
2622 cw1_sg(kt) = aw1/h_t_neu(j,i)*vz_t(kt,j,i)
2623  ! ... and kt=KTMAX-1
2624 
2625 #elif (ADV_VERT==2 || ADV_VERT==3)
2626 
2627 do kt=0, ktmax-1
2628  cw1_sg(kt) = aw1/h_t_neu(j,i)*vz_t(kt,j,i)
2629 end do
2630 
2631 #endif
2632 
2633 do kt=1, ktmax-1
2634  cw9(kt) = aw9 &
2635  *c_val(temp_t_m(kt,j,i)) &
2636  *( dzs_dtau(j,i) &
2637  +0.5_dp*(vx_t(kt,j,i)+vx_t(kt,j,i-1))*dzs_dxi_g(j,i) &
2638  +0.5_dp*(vy_t(kt,j,i)+vy_t(kt,j-1,i))*dzs_deta_g(j,i) &
2639  -0.5_dp*(vz_t(kt,j,i)+vz_t(kt-1,j,i)) )
2640 end do
2641 
2642 do kt=0, ktmax
2643 
2644  cw2(kt) = aw2*(dzb_dtau(j,i)+zeta_t(kt)*dh_t_dtau(j,i)) &
2645  /h_t_neu(j,i)
2646  cw3(kt) = aw3*(dzb_dxi_g(j,i)+zeta_t(kt)*dh_t_dxi_g(j,i)) &
2647  /h_t_neu(j,i) &
2648  *0.5_dp*(vx_t(kt,j,i)+vx_t(kt,j,i-1)) *insq_g11_g(j,i)
2649  cw4(kt) = aw4*(dzb_deta_g(j,i)+zeta_t(kt)*dh_t_deta_g(j,i)) &
2650  /h_t_neu(j,i) &
2651  *0.5_dp*(vy_t(kt,j,i)+vy_t(kt,j-1,i)) *insq_g22_g(j,i)
2652  sigma_t_help(kt) &
2653  = sigma_c_help(0) &
2654  + rho*g*h_t_neu(j,i)*(1.0_dp-zeta_t(kt)) &
2655  *(dzs_dxi_g(j,i)**2+dzs_deta_g(j,i)**2)**0.5_dp
2656 
2657 #if (DYNAMICS==2)
2658  if (.not.flag_shelfy_stream(j,i)) then
2659 #endif
2660  cw7(kt) = aw7 &
2661  *enh_t(kt,j,i) &
2662  *ratefac_t(omega_t(kt,j,i)) &
2663  *creep(sigma_t_help(kt)) &
2664  *sigma_t_help(kt)*sigma_t_help(kt)
2665 #if (DYNAMICS==2)
2666  else
2667  cw7(kt) = 2.0_dp*aw7 &
2668  *viscosity(de_t(kt,j,i), &
2669  temp_t_m(kt,j,i), temp_t_m(kt,j,i), &
2670  omega_t(kt,j,i), &
2671  enh_t(kt,j,i), 1_i1b) &
2672  *de_t(kt,j,i)**2
2673  end if
2674 #endif
2675 
2676 end do
2677 
2678 #if (ADV_VERT==1)
2679 
2680 kt=0
2681 cw2_sg(kt) = 0.5_dp*(cw2(kt)+cw2(kt+1))
2682 cw3_sg(kt) = 0.5_dp*(cw3(kt)+cw3(kt+1))
2683 cw4_sg(kt) = 0.5_dp*(cw4(kt)+cw4(kt+1))
2684 adv_vert_w_sg(kt) = cw1_sg(kt)-cw2_sg(kt)-cw3_sg(kt)-cw4_sg(kt)
2685 abs_adv_vert_w_sg(kt) = abs(adv_vert_w_sg(kt)) ! only needed for kt=0 ...
2686 kt=ktmax-1
2687 cw2_sg(kt) = 0.5_dp*(cw2(kt)+cw2(kt+1))
2688 cw3_sg(kt) = 0.5_dp*(cw3(kt)+cw3(kt+1))
2689 cw4_sg(kt) = 0.5_dp*(cw4(kt)+cw4(kt+1))
2690 adv_vert_w_sg(kt) = cw1_sg(kt)-cw2_sg(kt)-cw3_sg(kt)-cw4_sg(kt)
2691 abs_adv_vert_w_sg(kt) = abs(adv_vert_w_sg(kt)) ! ... and kt=KTMAX-1
2692 
2693 #elif (ADV_VERT==2 || ADV_VERT==3)
2694 
2695 do kt=0, ktmax-1
2696  cw2_sg(kt) = 0.5_dp*(cw2(kt)+cw2(kt+1))
2697  cw3_sg(kt) = 0.5_dp*(cw3(kt)+cw3(kt+1))
2698  cw4_sg(kt) = 0.5_dp*(cw4(kt)+cw4(kt+1))
2699  adv_vert_w_sg(kt) = cw1_sg(kt)-cw2_sg(kt)-cw3_sg(kt)-cw4_sg(kt)
2700  abs_adv_vert_w_sg(kt) = abs(adv_vert_w_sg(kt))
2701 end do
2702 
2703 #endif
2704 
2705 #if (ADV_HOR==3)
2706 dtt_dxi = 2.0_dp*dtt_2dxi
2707 dtt_deta = 2.0_dp*dtt_2deta
2708 #endif
2709 
2710 !-------- Set up the equations for the bedrock temperature --------
2711 
2712 kr=0
2713 lgs_a1(kr) = 1.0_dp
2714 lgs_a2(kr) = -1.0_dp
2715 lgs_b(kr) = clb1
2716 
2717 #if (Q_LITHO==1)
2718 ! (coupled heat-conducting bedrock)
2719 
2720 do kr=1, krmax-1
2721  lgs_a0(kr) = - ctr1
2722  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
2723  lgs_a2(kr) = - ctr1
2724  lgs_b(kr) = temp_r(kr,j,i)
2725 end do
2726 
2727 #elif (Q_LITHO==0)
2728 ! (no coupled heat-conducting bedrock)
2729 
2730 do kr=1, krmax-1
2731  lgs_a0(kr) = 1.0_dp
2732  lgs_a1(kr) = 0.0_dp
2733  lgs_a2(kr) = -1.0_dp
2734  lgs_b(kr) = 2.0_dp*clb1
2735 end do
2736 
2737 #endif
2738 
2739 kr=krmax
2740 lgs_a0(kr) = 0.0_dp
2741 lgs_a1(kr) = 1.0_dp
2742 lgs_b(kr) = temp_t_m(0,j,i)
2743 
2744 !-------- Solve system of linear equations --------
2745 
2746 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
2747 
2748 !-------- Assign the result --------
2749 
2750 do kr=0, krmax
2751  temp_r_neu(kr,j,i) = lgs_x(kr)
2752 end do
2753 
2754 !-------- Set up the equations for the water content in
2755 ! temperate ice --------
2756 
2757 kt=0
2758 lgs_a1(kt) = 1.0_dp
2759 lgs_a2(kt) = -1.0_dp
2760 lgs_b(kt) = 0.0_dp
2761 
2762 do kt=1, ktmax-1
2763 
2764 #if (ADV_VERT==1)
2765 
2766  lgs_a0(kt) = -0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - cw5
2767  lgs_a1(kt) = 1.0_dp + 2.0_dp*cw5
2768  lgs_a2(kt) = 0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - cw5
2769 
2770 #elif (ADV_VERT==2)
2771 
2772  lgs_a0(kt) &
2773  = -0.5_dp*(adv_vert_w_sg(kt-1)+abs_adv_vert_w_sg(kt-1)) &
2774  -cw5
2775  lgs_a1(kt) &
2776  = 1.0_dp &
2777  +0.5_dp*(adv_vert_w_sg(kt-1)+abs_adv_vert_w_sg(kt-1)) &
2778  -0.5_dp*(adv_vert_w_sg(kt) -abs_adv_vert_w_sg(kt) ) &
2779  +2.0_dp*cw5
2780  lgs_a2(kt) &
2781  = 0.5_dp*(adv_vert_w_sg(kt) -abs_adv_vert_w_sg(kt) ) &
2782  -cw5
2783 
2784 #elif (ADV_VERT==3)
2785 
2786  adv_vert_w_help = 0.5_dp*(adv_vert_w_sg(kt)+adv_vert_w_sg(kt-1))
2787 
2788  lgs_a0(kt) &
2789  = -max(adv_vert_w_help, 0.0_dp) &
2790  -cw5
2791  lgs_a1(kt) &
2792  = 1.0_dp &
2793  +max(adv_vert_w_help, 0.0_dp)-min(adv_vert_w_help, 0.0_dp) &
2794  +2.0_dp*cw5
2795  lgs_a2(kt) &
2796  = min(adv_vert_w_help, 0.0_dp) &
2797  -cw5
2798 
2799 #endif
2800 
2801 #if (ADV_HOR==2)
2802 
2803  lgs_b(kt) = omega_t(kt,j,i) + cw7(kt) + cw8 + cw9(kt) &
2804  -dtt_2dxi* &
2805  ( (vx_t(kt,j,i)-abs(vx_t(kt,j,i))) &
2806  *(omega_t(kt,j,i+1)-omega_t(kt,j,i)) &
2807  *insq_g11_sgx(j,i) &
2808  +(vx_t(kt,j,i-1)+abs(vx_t(kt,j,i-1))) &
2809  *(omega_t(kt,j,i)-omega_t(kt,j,i-1)) &
2810  *insq_g11_sgx(j,i-1) ) &
2811  -dtt_2deta* &
2812  ( (vy_t(kt,j,i)-abs(vy_t(kt,j,i))) &
2813  *(omega_t(kt,j+1,i)-omega_t(kt,j,i)) &
2814  *insq_g22_sgy(j,i) &
2815  +(vy_t(kt,j-1,i)+abs(vy_t(kt,j-1,i))) &
2816  *(omega_t(kt,j,i)-omega_t(kt,j-1,i)) &
2817  *insq_g22_sgy(j-1,i) )
2818 
2819 #elif (ADV_HOR==3)
2820 
2821  vx_t_help = 0.5_dp*(vx_t(kt,j,i)+vx_t(kt,j,i-1))
2822  vy_t_help = 0.5_dp*(vy_t(kt,j,i)+vy_t(kt,j-1,i))
2823 
2824  lgs_b(kt) = omega_t(kt,j,i) + cw7(kt) + cw8 + cw9(kt) &
2825  -dtt_dxi* &
2826  ( min(vx_t_help, 0.0_dp) &
2827  *(omega_t(kt,j,i+1)-omega_t(kt,j,i)) &
2828  *insq_g11_sgx(j,i) &
2829  +max(vx_t_help, 0.0_dp) &
2830  *(omega_t(kt,j,i)-omega_t(kt,j,i-1)) &
2831  *insq_g11_sgx(j,i-1) ) &
2832  -dtt_deta* &
2833  ( min(vy_t_help, 0.0_dp) &
2834  *(omega_t(kt,j+1,i)-omega_t(kt,j,i)) &
2835  *insq_g22_sgy(j,i) &
2836  +max(vy_t_help, 0.0_dp) &
2837  *(omega_t(kt,j,i)-omega_t(kt,j-1,i)) &
2838  *insq_g22_sgy(j-1,i) )
2839 
2840 #endif
2841 
2842 end do
2843 
2844 kt=ktmax
2845 
2846 #if (!defined(CTS_MELTING_FREEZING) || CTS_MELTING_FREEZING==1)
2847 
2848 if (am_perp(j,i) >= zero) then ! melting condition
2849  lgs_a0(kt) = 0.0_dp
2850  lgs_a1(kt) = 1.0_dp
2851  lgs_b(kt) = 0.0_dp
2852 else ! am_perp(j,i) < 0.0, freezing condition
2853  lgs_a0(kt) = -1.0_dp
2854  lgs_a1(kt) = 1.0_dp
2855  lgs_b(kt) = 0.0_dp
2856 end if
2857 
2858 #elif (CTS_MELTING_FREEZING==2)
2859 
2860 lgs_a0(kt) = 0.0_dp
2861 lgs_a1(kt) = 1.0_dp ! melting condition assumed
2862 lgs_b(kt) = 0.0_dp
2863 
2864 #else
2865 
2866 errormsg = ' >>> calc_temp3: CTS_MELTING_FREEZING must be either 1 or 2!'
2867 call error(errormsg)
2868 
2869 #endif
2870 
2871 !-------- Solve system of linear equations --------
2872 
2873 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, ktmax)
2874 
2875 !-------- Assign the result, compute the water drainage --------
2876 
2877 q_tld(j,i) = 0.0_dp
2878 
2879 do kt=0, ktmax
2880 
2881  if (lgs_x(kt) < zero) then
2882  omega_t_neu(kt,j,i) = 0.0_dp ! (as a precaution)
2883  else if (lgs_x(kt) < omega_max) then
2884  omega_t_neu(kt,j,i) = lgs_x(kt)
2885  else
2886  omega_t_neu(kt,j,i) = omega_max
2887  q_tld(j,i) = q_tld(j,i) &
2888  +aqtld*h_t_neu(j,i)*(lgs_x(kt)-omega_max)
2889  end if
2890 
2891 end do
2892 
2893 !-------- Set up the equations for the ice temperature --------
2894 
2895 ! ------ Abbreviation for the jump of the temperature gradient with
2896 ! the new omega
2897 
2898 #if (!defined(CTS_MELTING_FREEZING) || CTS_MELTING_FREEZING==1)
2899 
2900 if (am_perp(j,i) >= zero) then ! melting condition
2901  cm2 = 0.0_dp
2902 else ! am_perp(j,i) < 0.0, freezing condition
2903  cm2 = am2*h_c_neu(j,i)*omega_t_neu(ktmax,j,i)*am_perp(j,i) &
2904  /kappa_val(temp_c(0,j,i))
2905 end if
2906 
2907 #elif (CTS_MELTING_FREEZING==2)
2908 
2909 cm2 = 0.0_dp ! melting condition assumed
2910 
2911 #else
2912 
2913 errormsg = ' >>> calc_temp3: CTS_MELTING_FREEZING must be either 1 or 2!'
2914 call error(errormsg)
2915 
2916 #endif
2917 
2918 kc=0
2919 lgs_a1(kc) = 1.0_dp
2920 lgs_a2(kc) = -1.0_dp
2921 lgs_b(kc) = -cm1-cm2
2922 
2923 do kc=1, kcmax-1
2924 
2925 #if (ADV_VERT==1)
2926 
2927  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2928  -ct5(kc)*ct6(kc-1)
2929  lgs_a1(kc) = 1.0_dp+ct5(kc)*(ct6(kc)+ct6(kc-1))
2930  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2931  -ct5(kc)*ct6(kc)
2932 
2933 #elif (ADV_VERT==2)
2934 
2935  lgs_a0(kc) &
2936  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2937  -ct5(kc)*ct6(kc-1)
2938  lgs_a1(kc) &
2939  = 1.0_dp &
2940  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2941  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
2942  +ct5(kc)*(ct6(kc)+ct6(kc-1))
2943  lgs_a2(kc) &
2944  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
2945  -ct5(kc)*ct6(kc)
2946 
2947 #elif (ADV_VERT==3)
2948 
2949  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
2950 
2951  lgs_a0(kc) &
2952  = -max(adv_vert_help, 0.0_dp) &
2953  -ct5(kc)*ct6(kc-1)
2954  lgs_a1(kc) &
2955  = 1.0_dp &
2956  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
2957  +ct5(kc)*(ct6(kc)+ct6(kc-1))
2958  lgs_a2(kc) &
2959  = min(adv_vert_help, 0.0_dp) &
2960  -ct5(kc)*ct6(kc)
2961 
2962 #endif
2963 
2964 #if (ADV_HOR==2)
2965 
2966  lgs_b(kc) = temp_c(kc,j,i) + ct7(kc) &
2967  -dtt_2dxi* &
2968  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2969  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
2970  *insq_g11_sgx(j,i) &
2971  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2972  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
2973  *insq_g11_sgx(j,i-1) ) &
2974  -dtt_2deta* &
2975  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2976  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
2977  *insq_g22_sgy(j,i) &
2978  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2979  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
2980  *insq_g22_sgy(j-1,i) )
2981 
2982 #elif (ADV_HOR==3)
2983 
2984  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2985  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2986 
2987  lgs_b(kc) = temp_c(kc,j,i) + ct7(kc) &
2988  -dtt_dxi* &
2989  ( min(vx_c_help, 0.0_dp) &
2990  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
2991  *insq_g11_sgx(j,i) &
2992  +max(vx_c_help, 0.0_dp) &
2993  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
2994  *insq_g11_sgx(j,i-1) ) &
2995  -dtt_deta* &
2996  ( min(vy_c_help, 0.0_dp) &
2997  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
2998  *insq_g22_sgy(j,i) &
2999  +max(vy_c_help, 0.0_dp) &
3000  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
3001  *insq_g22_sgy(j-1,i) )
3002 
3003 #endif
3004 
3005 end do
3006 
3007 kc=kcmax
3008 lgs_a0(kc) = 0.0_dp
3009 lgs_a1(kc) = 1.0_dp
3010 lgs_b(kc) = temp_s(j,i)
3011 
3012 !-------- Solve system of linear equations --------
3013 
3014 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
3015 
3016 !-------- Assign the result --------
3017 
3018 do kc=0, kcmax
3019  temp_c_neu(kc,j,i) = lgs_x(kc)
3020 end do
3021 
3022 !-------- Set up the equations for the age (cold and temperate ice
3023 ! simultaneously) --------
3024 
3025 kt=0 ! adv_vert_w_sg(0) <= 0
3026 lgs_a1(kt) = 1.0_dp - min(adv_vert_w_sg(kt), 0.0_dp) ! (directed downward)
3027 lgs_a2(kt) = min(adv_vert_w_sg(kt), 0.0_dp) ! assumed/enforced
3028 
3029 #if (ADV_HOR==2)
3030 
3031 lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
3032  -dtt_2dxi* &
3033  ( (vx_t(kt,j,i)-abs(vx_t(kt,j,i))) &
3034  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
3035  *insq_g11_sgx(j,i) &
3036  +(vx_t(kt,j,i-1)+abs(vx_t(kt,j,i-1))) &
3037  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
3038  *insq_g11_sgx(j,i-1) ) &
3039  -dtt_2deta* &
3040  ( (vy_t(kt,j,i)-abs(vy_t(kt,j,i))) &
3041  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
3042  *insq_g22_sgy(j,i) &
3043  +(vy_t(kt,j-1,i)+abs(vy_t(kt,j-1,i))) &
3044  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
3045  *insq_g22_sgy(j-1,i) )
3046 
3047 #elif (ADV_HOR==3)
3048 
3049 vx_t_help = 0.5_dp*(vx_t(kt,j,i)+vx_t(kt,j,i-1))
3050 vy_t_help = 0.5_dp*(vy_t(kt,j,i)+vy_t(kt,j-1,i))
3051 
3052 lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
3053  -dtt_dxi* &
3054  ( min(vx_t_help, 0.0_dp) &
3055  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
3056  *insq_g11_sgx(j,i) &
3057  +max(vx_t_help, 0.0_dp) &
3058  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
3059  *insq_g11_sgx(j,i-1) ) &
3060  -dtt_deta* &
3061  ( min(vy_t_help, 0.0_dp) &
3062  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
3063  *insq_g22_sgy(j,i) &
3064  +max(vy_t_help, 0.0_dp) &
3065  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
3066  *insq_g22_sgy(j-1,i) )
3067 
3068 #endif
3069 
3070 do kt=1, ktmax-1
3071 
3072 #if (ADV_VERT==1)
3073 
3074  lgs_a0(kt) = -0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - ci3
3075  lgs_a1(kt) = 1.0_dp + 2.0_dp*ci3
3076  lgs_a2(kt) = 0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - ci3
3077 
3078 #elif (ADV_VERT==2)
3079 
3080  lgs_a0(kt) = -0.5_dp*(adv_vert_w_sg(kt-1)+abs_adv_vert_w_sg(kt-1))
3081  lgs_a1(kt) = 1.0_dp &
3082  +0.5_dp*(adv_vert_w_sg(kt-1)+abs_adv_vert_w_sg(kt-1)) &
3083  -0.5_dp*(adv_vert_w_sg(kt) -abs_adv_vert_w_sg(kt) )
3084  lgs_a2(kt) = 0.5_dp*(adv_vert_w_sg(kt) -abs_adv_vert_w_sg(kt) )
3085 
3086 #elif (ADV_VERT==3)
3087 
3088  adv_vert_w_help = 0.5_dp*(adv_vert_w_sg(kt)+adv_vert_w_sg(kt-1))
3089 
3090  lgs_a0(kt) = -max(adv_vert_w_help, 0.0_dp)
3091  lgs_a1(kt) = 1.0_dp &
3092  +max(adv_vert_w_help, 0.0_dp)-min(adv_vert_w_help, 0.0_dp)
3093  lgs_a2(kt) = min(adv_vert_w_help, 0.0_dp)
3094 
3095 #endif
3096 
3097 #if (ADV_HOR==2)
3098 
3099  lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
3100  -dtt_2dxi* &
3101  ( (vx_t(kt,j,i)-abs(vx_t(kt,j,i))) &
3102  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
3103  *insq_g11_sgx(j,i) &
3104  +(vx_t(kt,j,i-1)+abs(vx_t(kt,j,i-1))) &
3105  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
3106  *insq_g11_sgx(j,i-1) ) &
3107  -dtt_2deta* &
3108  ( (vy_t(kt,j,i)-abs(vy_t(kt,j,i))) &
3109  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
3110  *insq_g22_sgy(j,i) &
3111  +(vy_t(kt,j-1,i)+abs(vy_t(kt,j-1,i))) &
3112  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
3113  *insq_g22_sgy(j-1,i) )
3114 
3115 #elif (ADV_HOR==3)
3116 
3117  vx_t_help = 0.5_dp*(vx_t(kt,j,i)+vx_t(kt,j,i-1))
3118  vy_t_help = 0.5_dp*(vy_t(kt,j,i)+vy_t(kt,j-1,i))
3119 
3120  lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
3121  -dtt_dxi* &
3122  ( min(vx_t_help, 0.0_dp) &
3123  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
3124  *insq_g11_sgx(j,i) &
3125  +max(vx_t_help, 0.0_dp) &
3126  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
3127  *insq_g11_sgx(j,i-1) ) &
3128  -dtt_deta* &
3129  ( min(vy_t_help, 0.0_dp) &
3130  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
3131  *insq_g22_sgy(j,i) &
3132  +max(vy_t_help, 0.0_dp) &
3133  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
3134  *insq_g22_sgy(j-1,i) )
3135 
3136 #endif
3137 
3138 end do
3139 
3140 #if (ADV_VERT==1)
3141 
3142 kt=ktmax
3143 kc=0
3144 
3145 lgs_a0(kt) = -0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - ci3
3146 lgs_a1(kt) = 1.0_dp + 2.0_dp*ci3
3147 lgs_a2(kt) = 0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - ci3
3148 
3149 #if (ADV_HOR==2)
3150 
3151 lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
3152  -dtt_2dxi* &
3153  ( (vx_t(kt,j,i)-abs(vx_t(kt,j,i))) &
3154  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
3155  *insq_g11_sgx(j,i) &
3156  +(vx_t(kt,j,i-1)+abs(vx_t(kt,j,i-1))) &
3157  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
3158  *insq_g11_sgx(j,i-1) ) &
3159  -dtt_2deta* &
3160  ( (vy_t(kt,j,i)-abs(vy_t(kt,j,i))) &
3161  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
3162  *insq_g22_sgy(j,i) &
3163  +(vy_t(kt,j-1,i)+abs(vy_t(kt,j-1,i))) &
3164  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
3165  *insq_g22_sgy(j-1,i) )
3166 
3167 #elif (ADV_HOR==3)
3168 
3169 vx_t_help = 0.5_dp*(vx_t(kt,j,i)+vx_t(kt,j,i-1))
3170 vy_t_help = 0.5_dp*(vy_t(kt,j,i)+vy_t(kt,j-1,i))
3171 
3172 lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
3173  -dtt_dxi* &
3174  ( min(vx_t_help, 0.0_dp) &
3175  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
3176  *insq_g11_sgx(j,i) &
3177  +max(vx_t_help, 0.0_dp) &
3178  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
3179  *insq_g11_sgx(j,i-1) ) &
3180  -dtt_deta* &
3181  ( min(vy_t_help, 0.0_dp) &
3182  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
3183  *insq_g22_sgy(j,i) &
3184  +max(vy_t_help, 0.0_dp) &
3185  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
3186  *insq_g22_sgy(j-1,i) )
3187 
3188 #endif
3189 
3190 #elif (ADV_VERT==2 || ADV_VERT==3)
3191 
3192 kt=ktmax
3193 kc=0
3194 
3195 if (adv_vert_sg(kc) <= zero) then
3196 
3197  lgs_a0(ktmax+kc) = 0.0_dp
3198  lgs_a1(ktmax+kc) = 1.0_dp - adv_vert_sg(kc)
3199  lgs_a2(ktmax+kc) = adv_vert_sg(kc)
3200 
3201 #if (ADV_HOR==2)
3202 
3203  lgs_b(ktmax+kc) = age_c(kc,j,i) + dtime_temp &
3204  -dtt_2dxi* &
3205  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
3206  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
3207  *insq_g11_sgx(j,i) &
3208  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
3209  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
3210  *insq_g11_sgx(j,i-1) ) &
3211  -dtt_2deta* &
3212  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
3213  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
3214  *insq_g22_sgy(j,i) &
3215  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
3216  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
3217  *insq_g22_sgy(j-1,i) )
3218 
3219 #elif (ADV_HOR==3)
3220 
3221  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
3222  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
3223 
3224  lgs_b(ktmax+kc) = age_c(kc,j,i) + dtime_temp &
3225  -dtt_dxi* &
3226  ( min(vx_c_help, 0.0_dp) &
3227  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
3228  *insq_g11_sgx(j,i) &
3229  +max(vx_c_help, 0.0_dp) &
3230  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
3231  *insq_g11_sgx(j,i-1) ) &
3232  -dtt_deta* &
3233  ( min(vy_c_help, 0.0_dp) &
3234  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
3235  *insq_g22_sgy(j,i) &
3236  +max(vy_c_help, 0.0_dp) &
3237  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
3238  *insq_g22_sgy(j-1,i) )
3239 
3240 #endif
3241 
3242 else if (adv_vert_w_sg(kt-1) >= zero) then
3243 
3244  lgs_a0(kt) = -adv_vert_w_sg(kt-1)
3245  lgs_a1(kt) = 1.0_dp + adv_vert_w_sg(kt-1)
3246  lgs_a2(kt) = 0.0_dp
3247 
3248 #if (ADV_HOR==2)
3249 
3250  lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
3251  -dtt_2dxi* &
3252  ( (vx_t(kt,j,i)-abs(vx_t(kt,j,i))) &
3253  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
3254  *insq_g11_sgx(j,i) &
3255  +(vx_t(kt,j,i-1)+abs(vx_t(kt,j,i-1))) &
3256  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
3257  *insq_g11_sgx(j,i-1) ) &
3258  -dtt_2deta* &
3259  ( (vy_t(kt,j,i)-abs(vy_t(kt,j,i))) &
3260  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
3261  *insq_g22_sgy(j,i) &
3262  +(vy_t(kt,j-1,i)+abs(vy_t(kt,j-1,i))) &
3263  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
3264  *insq_g22_sgy(j-1,i) )
3265 
3266 #elif (ADV_HOR==3)
3267 
3268  vx_t_help = 0.5_dp*(vx_t(kt,j,i)+vx_t(kt,j,i-1))
3269  vy_t_help = 0.5_dp*(vy_t(kt,j,i)+vy_t(kt,j-1,i))
3270 
3271  lgs_b(kt) = age_t(kt,j,i) + dtime_temp &
3272  -dtt_dxi* &
3273  ( min(vx_t_help, 0.0_dp) &
3274  *(age_t(kt,j,i+1)-age_t(kt,j,i)) &
3275  *insq_g11_sgx(j,i) &
3276  +max(vx_t_help, 0.0_dp) &
3277  *(age_t(kt,j,i)-age_t(kt,j,i-1)) &
3278  *insq_g11_sgx(j,i-1) ) &
3279  -dtt_deta* &
3280  ( min(vy_t_help, 0.0_dp) &
3281  *(age_t(kt,j+1,i)-age_t(kt,j,i)) &
3282  *insq_g22_sgy(j,i) &
3283  +max(vy_t_help, 0.0_dp) &
3284  *(age_t(kt,j,i)-age_t(kt,j-1,i)) &
3285  *insq_g22_sgy(j-1,i) )
3286 
3287 #endif
3288 
3289 else
3290 
3291  lgs_a0(kt) = -0.5_dp
3292  lgs_a1(kt) = 1.0_dp
3293  lgs_a2(kt) = -0.5_dp
3294  lgs_b(kt) = 0.0_dp
3295  ! Makeshift: Average of age_c(kc=1) and age_t(kt=KTMAX-1)
3296 
3297 end if
3298 
3299 #endif
3300 
3301 do kc=1, kcmax-1
3302 
3303 #if (ADV_VERT==1)
3304 
3305  lgs_a0(ktmax+kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
3306  -ci1(kc)*ci2(kc-1)
3307  lgs_a1(ktmax+kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
3308  lgs_a2(ktmax+kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
3309  -ci1(kc)*ci2(kc)
3310 
3311 #elif (ADV_VERT==2)
3312 
3313  lgs_a0(ktmax+kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
3314  lgs_a1(ktmax+kc) = 1.0_dp &
3315  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
3316  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
3317  lgs_a2(ktmax+kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
3318 
3319 #elif (ADV_VERT==3)
3320 
3321  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
3322 
3323  lgs_a0(ktmax+kc) = -max(adv_vert_help, 0.0_dp)
3324  lgs_a1(ktmax+kc) = 1.0_dp &
3325  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp)
3326  lgs_a2(ktmax+kc) = min(adv_vert_help, 0.0_dp)
3327 
3328 #endif
3329 
3330 #if (ADV_HOR==2)
3331 
3332  lgs_b(ktmax+kc) = age_c(kc,j,i) + dtime_temp &
3333  -dtt_2dxi* &
3334  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
3335  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
3336  *insq_g11_sgx(j,i) &
3337  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
3338  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
3339  *insq_g11_sgx(j,i-1) ) &
3340  -dtt_2deta* &
3341  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
3342  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
3343  *insq_g22_sgy(j,i) &
3344  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
3345  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
3346  *insq_g22_sgy(j-1,i) )
3347 
3348 #elif (ADV_HOR==3)
3349 
3350  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
3351  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
3352 
3353  lgs_b(ktmax+kc) = age_c(kc,j,i) + dtime_temp &
3354  -dtt_dxi* &
3355  ( min(vx_c_help, 0.0_dp) &
3356  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
3357  *insq_g11_sgx(j,i) &
3358  +max(vx_c_help, 0.0_dp) &
3359  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
3360  *insq_g11_sgx(j,i-1) ) &
3361  -dtt_deta* &
3362  ( min(vy_c_help, 0.0_dp) &
3363  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
3364  *insq_g22_sgy(j,i) &
3365  +max(vy_c_help, 0.0_dp) &
3366  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
3367  *insq_g22_sgy(j-1,i) )
3368 
3369 #endif
3370 
3371 end do
3372 
3373 kc=kcmax
3374 if (as_perp(j,i) >= zero) then
3375  lgs_a0(ktmax+kc) = 0.0_dp
3376  lgs_a1(ktmax+kc) = 1.0_dp
3377  lgs_b(ktmax+kc) = 0.0_dp
3378 else
3379  lgs_a0(ktmax+kc) = -max(adv_vert_sg(kc-1), 0.0_dp)
3380  lgs_a1(ktmax+kc) = 1.0_dp + max(adv_vert_sg(kc-1), 0.0_dp)
3381  ! adv_vert_sg(KCMAX-1) >= 0 (directed upward)
3382  ! assumed/enforced
3383 #if (ADV_HOR==2)
3384 
3385  lgs_b(ktmax+kc) = age_c(kc,j,i) + dtime_temp &
3386  -dtt_2dxi* &
3387  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
3388  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
3389  *insq_g11_sgx(j,i) &
3390  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
3391  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
3392  *insq_g11_sgx(j,i-1) ) &
3393  -dtt_2deta* &
3394  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
3395  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
3396  *insq_g22_sgy(j,i) &
3397  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
3398  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
3399  *insq_g22_sgy(j-1,i) )
3400 
3401 #elif (ADV_HOR==3)
3402 
3403  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
3404  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
3405 
3406  lgs_b(ktmax+kc) = age_c(kc,j,i) + dtime_temp &
3407  -dtt_dxi* &
3408  ( min(vx_c_help, 0.0_dp) &
3409  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
3410  *insq_g11_sgx(j,i) &
3411  +max(vx_c_help, 0.0_dp) &
3412  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
3413  *insq_g11_sgx(j,i-1) ) &
3414  -dtt_deta* &
3415  ( min(vy_c_help, 0.0_dp) &
3416  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
3417  *insq_g22_sgy(j,i) &
3418  +max(vy_c_help, 0.0_dp) &
3419  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
3420  *insq_g22_sgy(j-1,i) )
3421 
3422 #endif
3423 
3424 end if
3425 
3426 !-------- Solve system of linear equations --------
3427 
3428 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax+ktmax)
3429 
3430 !-------- Assign the result,
3431 ! restriction to interval [0, AGE_MAX yr] --------
3432 
3433 do kt=0, ktmax
3434 
3435  age_t_neu(kt,j,i) = lgs_x(kt)
3436 
3437  if (age_t_neu(kt,j,i) < (age_min*year_sec)) &
3438  age_t_neu(kt,j,i) = 0.0_dp
3439  if (age_t_neu(kt,j,i) > (age_max*year_sec)) &
3440  age_t_neu(kt,j,i) = age_max*year_sec
3441 
3442 end do
3443 
3444 do kc=0, kcmax
3445 
3446  age_c_neu(kc,j,i) = lgs_x(ktmax+kc)
3447 
3448  if (age_c_neu(kc,j,i) < (age_min*year_sec)) &
3449  age_c_neu(kc,j,i) = 0.0_dp
3450  if (age_c_neu(kc,j,i) > (age_max*year_sec)) &
3451  age_c_neu(kc,j,i) = age_max*year_sec
3452 
3453 end do
3454 
3455 end subroutine calc_temp3
3456 
3457 !-------------------------------------------------------------------------------
3458 !> Computation of temperature and age for an ice-free column.
3459 !<------------------------------------------------------------------------------
3460 subroutine calc_temp_r(atr1, alb1, i, j)
3461 
3462 #ifndef ALLOW_OPENAD /* Normal */
3463 use sico_maths_m, only : tri_sle
3464 #else /* OpenAD */
3465 use sico_maths_m
3466 #endif /* Normal vs. OpenAD */
3467 
3468 implicit none
3469 
3470 #ifndef ALLOW_OPENAD /* Normal */
3471 integer(i4b), intent(in) :: i, j
3472 #else /* OpenAD */
3473 integer(i4b), intent(inout) :: i, j
3474 #endif /* Normal vs. OpenAD */
3475 real(dp), intent(in) :: atr1, alb1
3476 
3477 integer(i4b) :: kc, kt, kr
3478 real(dp) :: ctr1, clb1
3479 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
3480  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
3481  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
3482  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
3483  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
3484 
3485 !-------- Abbreviations --------
3486 
3487 ctr1 = atr1
3488 clb1 = alb1*q_geo(j,i)
3489 
3490 !-------- Set up the equations for the bedrock temperature --------
3491 
3492 kr=0
3493 #ifdef ALLOW_OPENAD /* OpenAD */
3494 lgs_a0(kr) = 0.0_dp
3495 #endif /* OpenAD */
3496 lgs_a1(kr) = 1.0_dp
3497 lgs_a2(kr) = -1.0_dp
3498 lgs_b(kr) = clb1
3499 
3500 #if (Q_LITHO==1)
3501 ! (coupled heat-conducting bedrock)
3502 
3503 do kr=1, krmax-1
3504  lgs_a0(kr) = - ctr1
3505  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
3506  lgs_a2(kr) = - ctr1
3507  lgs_b(kr) = temp_r(kr,j,i)
3508 end do
3509 
3510 #elif (Q_LITHO==0)
3511 ! (no coupled heat-conducting bedrock)
3512 
3513 do kr=1, krmax-1
3514  lgs_a0(kr) = 1.0_dp
3515  lgs_a1(kr) = 0.0_dp
3516  lgs_a2(kr) = -1.0_dp
3517  lgs_b(kr) = 2.0_dp*clb1
3518 end do
3519 
3520 #endif
3521 
3522 kr=krmax
3523 lgs_a0(kr) = 0.0_dp
3524 lgs_a1(kr) = 1.0_dp
3525 lgs_b(kr) = temp_s(j,i)
3526 
3527 !-------- Solve system of linear equations --------
3528 
3529 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
3530 
3531 !-------- Assign the result --------
3532 
3533 do kr=0, krmax
3534  temp_r_neu(kr,j,i) = lgs_x(kr)
3535 end do
3536 
3537 !-------- Water content and age
3538 ! in the non-existing lower (kt) ice layer --------
3539 
3540 do kt=0, ktmax
3541  omega_t_neu(kt,j,i) = 0.0_dp
3542  age_t_neu(kt,j,i) = 0.0_dp
3543 end do
3544 
3545 !-------- Temperature and age
3546 ! in the non-existing upper (kc) ice layer --------
3547 
3548 do kc=0, kcmax
3549  temp_c_neu(kc,j,i) = temp_s(j,i)
3550  age_c_neu(kc,j,i) = 0.0_dp
3551 end do
3552 
3553 end subroutine calc_temp_r
3554 
3555 !-------------------------------------------------------------------------------
3556 !> Upward shifting of the CTS.
3557 !<------------------------------------------------------------------------------
3558 subroutine shift_cts_upward(at1, at2_1, at2_2, at3_1, at3_2, &
3559  at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
3560  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
3561  ai1, ai2, ai3, dzeta_t, &
3562  dtime_temp, dtt_2dxi, dtt_2deta, dtime_temp_inv, &
3563  i, j)
3564 
3565 implicit none
3566 
3567 #ifndef ALLOW_OPENAD /* Normal */
3568 integer(i4b), intent(in) :: i, j
3569 #else /* OpenAD */
3570 integer(i4b), intent(inout) :: i, j
3571 #endif /* Normal vs. OpenAD */
3572 real(dp), intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
3573  at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
3574  at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
3575  ai1(0:kcmax), ai2(0:kcmax), ai3, &
3576  atr1, am1, am2, alb1
3577 real(dp), intent(in) :: aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld
3578 real(dp), intent(in) :: dzeta_t
3579 real(dp), intent(in) :: dtime_temp, dtime_temp_inv, dtt_2dxi, dtt_2deta
3580 
3581 real(dp) :: zm_shift
3582 real(dp) :: difftemp_a, difftemp_b, interpol
3583 
3584 zm_shift = 1.0_dp ! CTS shift in intervals of 1 m
3585 
3586 !-------- Temperature discrepancy from the computation of the main
3587 ! program --------
3588 
3589 difftemp_a = temp_c_neu(0,j,i)-(-beta*h_c_neu(j,i))
3590 if (difftemp_a <= 0.0_dp) return
3591 
3592 !-------- Shift CTS upward until it is too high --------
3593 
3594 do while (difftemp_a > 0.0_dp)
3595 
3596  zm_neu(j,i) = zm_neu(j,i) + zm_shift
3597  if (zm_neu(j,i) >= zs(j,i)) then
3598  zm_neu(j,i) = zm_neu(j,i) - zm_shift
3599  return
3600  end if
3601  h_c_neu(j,i) = h_c_neu(j,i) - zm_shift
3602  h_t_neu(j,i) = h_t_neu(j,i) + zm_shift
3603 
3604  dh_t_dtau(j,i) = (zm_neu(j,i)-zm(j,i))*dtime_temp_inv
3605  dzm_dtau(j,i) = dzb_dtau(j,i)+dh_t_dtau(j,i)
3606  dh_c_dtau(j,i) = dzs_dtau(j,i)-dzm_dtau(j,i)
3607 
3608  am_perp(j,i) = am_perp_st(j,i) + dzm_dtau(j,i)
3609 
3610  call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
3611  at4_1, at4_2, at5, at6, at7, atr1, &
3612  am1, am2, alb1, &
3613  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
3614  ai1, ai2, ai3, dzeta_t, &
3615  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
3616 
3617  difftemp_b = difftemp_a
3618  difftemp_a = temp_c_neu(0,j,i)-(-beta*h_c_neu(j,i))
3619 
3620 end do
3621 
3622 !-------- Interpolate the CTS position from the last (_a) and the
3623 ! last but one (_b) value, weighed with the temperature
3624 ! discrepancies at the CTS --------
3625 
3626 interpol = difftemp_a/(difftemp_b-difftemp_a)*zm_shift
3627 
3628 zm_neu(j,i) = zm_neu(j,i) + interpol
3629 h_c_neu(j,i) = h_c_neu(j,i) - interpol
3630 h_t_neu(j,i) = h_t_neu(j,i) + interpol
3631 
3632 dh_t_dtau(j,i) = (zm_neu(j,i)-zm(j,i))*dtime_temp_inv
3633 dzm_dtau(j,i) = dzb_dtau(j,i)+dh_t_dtau(j,i)
3634 dh_c_dtau(j,i) = dzs_dtau(j,i)-dzm_dtau(j,i)
3635 
3636 am_perp(j,i) = am_perp_st(j,i) + dzm_dtau(j,i)
3637 
3638 call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
3639  at4_1, at4_2, at5, at6, at7, atr1, &
3640  am1, am2, alb1, &
3641  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
3642  ai1, ai2, ai3, dzeta_t, &
3643  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
3644 
3645 end subroutine shift_cts_upward
3646 
3647 !-------------------------------------------------------------------------------
3648 !> Downward shifting of the CTS.
3649 !<------------------------------------------------------------------------------
3650 subroutine shift_cts_downward(at1, at2_1, at2_2, at3_1, at3_2, &
3651  at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
3652  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
3653  ai1, ai2, ai3, dzeta_t, &
3654  dtime_temp, dtt_2dxi, dtt_2deta, dtime_temp_inv, &
3655  i, j)
3656 
3657 implicit none
3658 
3659 #ifndef ALLOW_OPENAD /* Normal */
3660 integer(i4b), intent(in) :: i, j
3661 #else /* OpenAD */
3662 integer(i4b), intent(inout) :: i, j
3663 #endif /* Normal vs. OpenAD */
3664 real(dp), intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
3665  at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
3666  at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
3667  ai1(0:kcmax), ai2(0:kcmax), ai3, &
3668  atr1, am1, am2, alb1
3669 real(dp), intent(in) :: aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld
3670 real(dp), intent(in) :: dzeta_t
3671 real(dp), intent(in) :: dtt_2dxi, dtt_2deta, dtime_temp, dtime_temp_inv
3672 
3673 real(dp) :: zm_shift
3674 real(dp) :: difftemp_a, difftemp_b, interpol
3675 
3676 zm_shift = 1.0_dp ! CTS shift in intervals of 1 m
3677 
3678 !-------- Temperature discrepancy from the computation of the main
3679 ! program --------
3680 
3681 difftemp_a = temp_c_neu(0,j,i)-(-beta*h_c_neu(j,i))
3682 if (difftemp_a >= 0.0_dp) return
3683 
3684 !-------- Shift CTS downward until it is too low --------
3685 
3686 do while (difftemp_a < 0.0_dp)
3687 
3688  zm_neu(j,i) = zm_neu(j,i) - zm_shift
3689 
3690 ! ------ Special case: CTS too close to the base
3691 
3692  if (zm_neu(j,i) <= zb(j,i)) then
3693 
3694  zm_shift = (zm_neu(j,i)+zm_shift)-(zb(j,i)+0.001_dp)
3695  zm_neu(j,i) = zb(j,i)+0.001_dp
3696  h_c_neu(j,i) = h_c_neu(j,i)+h_t_neu(j,i)-0.001_dp
3697  h_t_neu(j,i) = 0.001_dp
3698 ! ! CTS positioned 1 mm above ice base --------
3699 
3700  dh_t_dtau(j,i) = (zm_neu(j,i)-zm(j,i))*dtime_temp_inv
3701  dzm_dtau(j,i) = dzb_dtau(j,i)+dh_t_dtau(j,i)
3702  dh_c_dtau(j,i) = dzs_dtau(j,i)-dzm_dtau(j,i)
3703 
3704  am_perp(j,i) = am_perp_st(j,i) + dzm_dtau(j,i)
3705 
3706  call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
3707  at4_1, at4_2, at5, at6, at7, atr1, &
3708  am1, am2, alb1, &
3709  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
3710  ai1, ai2, ai3, dzeta_t, &
3711  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
3712 
3713  difftemp_b = difftemp_a
3714  difftemp_a = temp_c_neu(0,j,i)-(-beta*h_c_neu(j,i))
3715 
3716  if (difftemp_a >= 0.0_dp) then ! CTS remains above the base
3717 
3718 ! ---- Interpolate the CTS position from the last (_a) and the
3719 ! last but one (_b) value, weighed with the temperature
3720 ! discrepancies at the CTS --------
3721 
3722  interpol = difftemp_a/(difftemp_a-difftemp_b)*zm_shift
3723 
3724  zm_neu(j,i) = zm_neu(j,i) + interpol
3725  h_c_neu(j,i) = h_c_neu(j,i) - interpol
3726  h_t_neu(j,i) = h_t_neu(j,i) + interpol
3727 
3728  dh_t_dtau(j,i) = (zm_neu(j,i)-zm(j,i))*dtime_temp_inv
3729  dzm_dtau(j,i) = dzb_dtau(j,i)+dh_t_dtau(j,i)
3730  dh_c_dtau(j,i) = dzs_dtau(j,i)-dzm_dtau(j,i)
3731 
3732  am_perp(j,i) = am_perp_st(j,i) + dzm_dtau(j,i)
3733 
3734  call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
3735  at4_1, at4_2, at5, at6, at7, atr1, &
3736  am1, am2, alb1, &
3737  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
3738  ai1, ai2, ai3, dzeta_t, &
3739  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
3740 
3741  else ! CTS disappears
3742 
3743  n_cts_neu(j,i) = 0_i1b
3744  zm_neu(j,i) = zb(j,i)
3745  h_c_neu(j,i) = h_c_neu(j,i)+h_t_neu(j,i)
3746  h_t_neu(j,i) = 0.0_dp
3747 
3748  dh_t_dtau(j,i) = (zm_neu(j,i)-zm(j,i))*dtime_temp_inv
3749  dzm_dtau(j,i) = dzb_dtau(j,i)+dh_t_dtau(j,i)
3750  dh_c_dtau(j,i) = dzs_dtau(j,i)-dzm_dtau(j,i)
3751 
3752  am_perp(j,i) = am_perp_st(j,i) + dzm_dtau(j,i)
3753 
3754  call calc_temp2(at1, at2_1, at2_2, at3_1, at3_2, &
3755  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
3756  ai1, ai2, &
3757  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
3758 
3759  end if
3760 
3761  return
3762  end if
3763 
3764 ! ------ End of treatment of special case
3765 
3766  h_c_neu(j,i) = h_c_neu(j,i) + zm_shift
3767  h_t_neu(j,i) = h_t_neu(j,i) - zm_shift
3768 
3769  dh_t_dtau(j,i) = (zm_neu(j,i)-zm(j,i))*dtime_temp_inv
3770  dzm_dtau(j,i) = dzb_dtau(j,i)+dh_t_dtau(j,i)
3771  dh_c_dtau(j,i) = dzs_dtau(j,i)-dzm_dtau(j,i)
3772 
3773  am_perp(j,i) = am_perp_st(j,i) + dzm_dtau(j,i)
3774 
3775  call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
3776  at4_1, at4_2, at5, at6, at7, atr1, &
3777  am1, am2, alb1, &
3778  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
3779  ai1, ai2, ai3, dzeta_t, &
3780  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
3781 
3782  difftemp_b = difftemp_a
3783  difftemp_a = temp_c_neu(0,j,i)-(-beta*h_c_neu(j,i))
3784 
3785 end do
3786 
3787 !-------- Interpolate the CTS position from the last (_a) and the
3788 ! last but one (_b) value, weighed with the temperature
3789 ! discrepancies at the CTS --------
3790 
3791 interpol = difftemp_a/(difftemp_a-difftemp_b)*zm_shift
3792 
3793 zm_neu(j,i) = zm_neu(j,i) + interpol
3794 h_c_neu(j,i) = h_c_neu(j,i) - interpol
3795 h_t_neu(j,i) = h_t_neu(j,i) + interpol
3796 
3797 dh_t_dtau(j,i) = (zm_neu(j,i)-zm(j,i))*dtime_temp_inv
3798 dzm_dtau(j,i) = dzb_dtau(j,i)+dh_t_dtau(j,i)
3799 dh_c_dtau(j,i) = dzs_dtau(j,i)-dzm_dtau(j,i)
3800 
3801 am_perp(j,i) = am_perp_st(j,i) + dzm_dtau(j,i)
3802 
3803 call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
3804  at4_1, at4_2, at5, at6, at7, atr1, &
3805  am1, am2, alb1, &
3806  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
3807  ai1, ai2, ai3, dzeta_t, &
3808  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
3809 
3810 end subroutine shift_cts_downward
3811 
3812 !-------------------------------------------------------------------------------
3813 !> Computation of temperature and age for ice shelves (floating ice).
3814 !<------------------------------------------------------------------------------
3815 subroutine calc_temp_ssa(at1, at2_1, at2_2, at3_1, at3_2, &
3816  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
3817  ai1, ai2, &
3818  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
3819 
3820 use ice_material_properties_m, only : kappa_val, c_val, viscosity
3821 
3822 #ifndef ALLOW_OPENAD /* Normal */
3823 use sico_maths_m, only : tri_sle
3824 #else /* OpenAD */
3825 use sico_maths_m
3826 #endif /* Normal vs. OpenAD */
3827 
3828 implicit none
3829 
3830 #ifndef ALLOW_OPENAD /* Normal */
3831 integer(i4b), intent(in) :: i, j
3832 #else /* OpenAD */
3833 integer(i4b), intent(inout) :: i, j
3834 #endif /* Normal vs. OpenAD */
3835 real(dp), intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
3836  at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
3837  at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
3838  ai1(0:kcmax), ai2(0:kcmax), &
3839  atr1, alb1
3840 real(dp), intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
3841 
3842 integer(i4b) :: kc, kt, kr
3843 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
3844  ct5(0:kcmax), ct6(0:kcmax), ct7(0:kcmax), ctr1, clb1
3845 real(dp) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), ct3_sg(0:kcmax), &
3846  ct4_sg(0:kcmax), adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
3847 real(dp) :: ci1(0:kcmax), ci2(0:kcmax)
3848 real(dp) :: temp_c_help(0:kcmax)
3849 real(dp) :: vx_c_help, vy_c_help
3850 real(dp) :: adv_vert_help
3851 real(dp) :: dtt_dxi, dtt_deta
3852 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
3853  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
3854  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
3855  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
3856  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
3857 real(dp), parameter :: zero=0.0_dp
3858 
3859 !-------- Check for boundary points --------
3860 
3861 if ((i == 0).or.(i == imax).or.(j == 0).or.(j == jmax)) then
3862  errormsg = ' >>> calc_temp_ssa: Boundary points not allowed!'
3863  call error(errormsg)
3864 end if
3865 
3866 !-------- Abbreviations --------
3867 
3868 ctr1 = atr1
3869 clb1 = alb1*q_geo(j,i)
3870 
3871 #if (ADV_VERT==1)
3872 
3873 do kc=1, kcmax-1
3874  ct1(kc) = at1(kc)/h_c(j,i)*0.5_dp*(vz_c(kc,j,i)+vz_c(kc-1,j,i))
3875 end do
3876 
3877 kc=0
3878 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
3879  ! only needed for kc=0 ...
3880 kc=kcmax-1
3881 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
3882  ! ... and kc=KCMAX-1
3883 
3884 #elif (ADV_VERT==2 || ADV_VERT==3)
3885 
3886 do kc=0, kcmax-1
3887  ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
3888 end do
3889 
3890 #endif
3891 
3892 do kc=0, kcmax
3893 
3894  ct2(kc) = ( at2_1(kc)*dzm_dtau(j,i) &
3895  +at2_2(kc)*dh_c_dtau(j,i) )/h_c(j,i)
3896  ct3(kc) = ( at3_1(kc)*dzm_dxi_g(j,i) &
3897  +at3_2(kc)*dh_c_dxi_g(j,i) )/h_c(j,i) &
3898  *0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1)) *insq_g11_g(j,i)
3899  ct4(kc) = ( at4_1(kc)*dzm_deta_g(j,i) &
3900  +at4_2(kc)*dh_c_deta_g(j,i) )/h_c(j,i) &
3901  *0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i)) *insq_g22_g(j,i)
3902  ct5(kc) = at5(kc) &
3903  /c_val(temp_c(kc,j,i)) &
3904  /h_c(j,i)
3905  ct7(kc) = 2.0_dp*at7 &
3906  /c_val(temp_c(kc,j,i)) &
3907  *viscosity(de_ssa(j,i), &
3908  temp_c(kc,j,i), temp_c_m(kc,j,i), 0.0_dp, &
3909  enh_c(kc,j,i), 0_i1b) &
3910  *de_ssa(j,i)**2
3911  ci1(kc) = ai1(kc)/h_c(j,i)
3912 
3913 end do
3914 
3915 #if (ADV_VERT==1)
3916 
3917 kc=0
3918 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
3919 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
3920 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
3921 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
3922 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! only needed for kc=0 ...
3923 kc=kcmax-1
3924 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
3925 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
3926 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
3927 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
3928 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! ... and kc=KCMAX-1
3929 
3930 #elif (ADV_VERT==2 || ADV_VERT==3)
3931 
3932 do kc=0, kcmax-1
3933  ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
3934  ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
3935  ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
3936  adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
3937  abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
3938 end do
3939 
3940 #endif
3941 
3942 do kc=0, kcmax-1
3943  temp_c_help(kc) = 0.5_dp*(temp_c(kc,j,i)+temp_c(kc+1,j,i))
3944  ct6(kc) = at6(kc) &
3945  *kappa_val(temp_c_help(kc)) &
3946  /h_c(j,i)
3947  ci2(kc) = ai2(kc)/h_c(j,i)
3948 end do
3949 
3950 #if (ADV_HOR==3)
3951 dtt_dxi = 2.0_dp*dtt_2dxi
3952 dtt_deta = 2.0_dp*dtt_2deta
3953 #endif
3954 
3955 !-------- Set up the equations for the bedrock temperature --------
3956 
3957 kr=0
3958 lgs_a1(kr) = 1.0_dp
3959 lgs_a2(kr) = -1.0_dp
3960 lgs_b(kr) = clb1
3961 
3962 #if (Q_LITHO==1)
3963 ! (coupled heat-conducting bedrock)
3964 
3965 do kr=1, krmax-1
3966  lgs_a0(kr) = - ctr1
3967  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
3968  lgs_a2(kr) = - ctr1
3969  lgs_b(kr) = temp_r(kr,j,i)
3970 end do
3971 
3972 #elif (Q_LITHO==0)
3973 ! (no coupled heat-conducting bedrock)
3974 
3975 do kr=1, krmax-1
3976  lgs_a0(kr) = 1.0_dp
3977  lgs_a1(kr) = 0.0_dp
3978  lgs_a2(kr) = -1.0_dp
3979  lgs_b(kr) = 2.0_dp*clb1
3980 end do
3981 
3982 #endif
3983 
3984 kr=krmax
3985 lgs_a0(kr) = 0.0_dp
3986 lgs_a1(kr) = 1.0_dp
3987 lgs_b(kr) = temp_c_m(0,j,i)-delta_tm_sw
3988 
3989 !-------- Solve system of linear equations --------
3990 
3991 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
3992 
3993 !-------- Assign the result --------
3994 
3995 do kr=0, krmax
3996  temp_r_neu(kr,j,i) = lgs_x(kr)
3997 end do
3998 
3999 !-------- Set up the equations for the ice temperature --------
4000 
4001 kc=0
4002 lgs_a1(kc) = 1.0_dp
4003 lgs_a2(kc) = 0.0_dp
4004 lgs_b(kc) = temp_c_m(0,j,i)-delta_tm_sw
4005 
4006 do kc=1, kcmax-1
4007 
4008 #if (ADV_VERT==1)
4009 
4010  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
4011  -ct5(kc)*ct6(kc-1)
4012  lgs_a1(kc) = 1.0_dp+ct5(kc)*(ct6(kc)+ct6(kc-1))
4013  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
4014  -ct5(kc)*ct6(kc)
4015 
4016 #elif (ADV_VERT==2)
4017 
4018  lgs_a0(kc) &
4019  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
4020  -ct5(kc)*ct6(kc-1)
4021  lgs_a1(kc) &
4022  = 1.0_dp &
4023  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
4024  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
4025  +ct5(kc)*(ct6(kc)+ct6(kc-1))
4026  lgs_a2(kc) &
4027  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
4028  -ct5(kc)*ct6(kc)
4029 
4030 #elif (ADV_VERT==3)
4031 
4032  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
4033 
4034  lgs_a0(kc) &
4035  = -max(adv_vert_help, 0.0_dp) &
4036  -ct5(kc)*ct6(kc-1)
4037  lgs_a1(kc) &
4038  = 1.0_dp &
4039  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
4040  +ct5(kc)*(ct6(kc)+ct6(kc-1))
4041  lgs_a2(kc) &
4042  = min(adv_vert_help, 0.0_dp) &
4043  -ct5(kc)*ct6(kc)
4044 
4045 #endif
4046 
4047 #if (ADV_HOR==2)
4048 
4049  lgs_b(kc) = temp_c(kc,j,i) + ct7(kc) &
4050  -dtt_2dxi* &
4051  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
4052  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
4053  *insq_g11_sgx(j,i) &
4054  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
4055  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
4056  *insq_g11_sgx(j,i-1) ) &
4057  -dtt_2deta* &
4058  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
4059  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
4060  *insq_g22_sgy(j,i) &
4061  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
4062  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
4063  *insq_g22_sgy(j-1,i) )
4064 
4065 #elif (ADV_HOR==3)
4066 
4067  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
4068  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
4069 
4070  lgs_b(kc) = temp_c(kc,j,i) + ct7(kc) &
4071  -dtt_dxi* &
4072  ( min(vx_c_help, 0.0_dp) &
4073  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
4074  *insq_g11_sgx(j,i) &
4075  +max(vx_c_help, 0.0_dp) &
4076  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
4077  *insq_g11_sgx(j,i-1) ) &
4078  -dtt_deta* &
4079  ( min(vy_c_help, 0.0_dp) &
4080  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
4081  *insq_g22_sgy(j,i) &
4082  +max(vy_c_help, 0.0_dp) &
4083  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
4084  *insq_g22_sgy(j-1,i) )
4085 
4086 #endif
4087 
4088 end do
4089 
4090 kc=kcmax
4091 lgs_a0(kc) = 0.0_dp
4092 lgs_a1(kc) = 1.0_dp
4093 lgs_b(kc) = temp_s(j,i)
4094 
4095 !-------- Solve system of linear equations --------
4096 
4097 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
4098 
4099 !-------- Assign the result --------
4100 
4101 do kc=0, kcmax
4102  temp_c_neu(kc,j,i) = lgs_x(kc)
4103 end do
4104 
4105 !-------- Set water content in the non-existing temperate layer
4106 ! to zero --------
4107 
4108 do kt=0, ktmax
4109  omega_t_neu(kt,j,i) = 0.0_dp
4110 end do
4111 
4112 !-------- Water drainage from the non-existing temperate layer --------
4113 
4114 q_tld(j,i) = 0.0_dp
4115 
4116 !-------- Set up the equations for the age of cold ice --------
4117 
4118 kc=0 ! adv_vert_sg(0) <= 0
4119 lgs_a1(kc) = 1.0_dp - min(adv_vert_sg(kc), 0.0_dp) ! (directed downward)
4120 lgs_a2(kc) = min(adv_vert_sg(kc), 0.0_dp) ! assumed/enforced
4121 
4122 #if (ADV_HOR==2)
4123 
4124 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
4125  -dtt_2dxi* &
4126  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
4127  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
4128  *insq_g11_sgx(j,i) &
4129  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
4130  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
4131  *insq_g11_sgx(j,i-1) ) &
4132  -dtt_2deta* &
4133  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
4134  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
4135  *insq_g22_sgy(j,i) &
4136  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
4137  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
4138  *insq_g22_sgy(j-1,i) )
4139 
4140 #elif (ADV_HOR==3)
4141 
4142 vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
4143 vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
4144 
4145 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
4146  -dtt_dxi* &
4147  ( min(vx_c_help, 0.0_dp) &
4148  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
4149  *insq_g11_sgx(j,i) &
4150  +max(vx_c_help, 0.0_dp) &
4151  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
4152  *insq_g11_sgx(j,i-1) ) &
4153  -dtt_deta* &
4154  ( min(vy_c_help, 0.0_dp) &
4155  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
4156  *insq_g22_sgy(j,i) &
4157  +max(vy_c_help, 0.0_dp) &
4158  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
4159  *insq_g22_sgy(j-1,i) )
4160 
4161 #endif
4162 
4163 do kc=1, kcmax-1
4164 
4165 #if (ADV_VERT==1)
4166 
4167  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
4168  -ci1(kc)*ci2(kc-1)
4169  lgs_a1(kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
4170  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
4171  -ci1(kc)*ci2(kc)
4172 
4173 #elif (ADV_VERT==2)
4174 
4175  lgs_a0(kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
4176  lgs_a1(kc) = 1.0_dp &
4177  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
4178  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
4179  lgs_a2(kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
4180 
4181 #elif (ADV_VERT==3)
4182 
4183  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
4184 
4185  lgs_a0(kc) = -max(adv_vert_help, 0.0_dp)
4186  lgs_a1(kc) = 1.0_dp &
4187  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp)
4188  lgs_a2(kc) = min(adv_vert_help, 0.0_dp)
4189 
4190 #endif
4191 
4192 #if (ADV_HOR==2)
4193 
4194  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
4195  -dtt_2dxi* &
4196  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
4197  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
4198  *insq_g11_sgx(j,i) &
4199  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
4200  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
4201  *insq_g11_sgx(j,i-1) ) &
4202  -dtt_2deta* &
4203  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
4204  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
4205  *insq_g22_sgy(j,i) &
4206  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
4207  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
4208  *insq_g22_sgy(j-1,i) )
4209 
4210 #elif (ADV_HOR==3)
4211 
4212  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
4213  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
4214 
4215  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
4216  -dtt_dxi* &
4217  ( min(vx_c_help, 0.0_dp) &
4218  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
4219  *insq_g11_sgx(j,i) &
4220  +max(vx_c_help, 0.0_dp) &
4221  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
4222  *insq_g11_sgx(j,i-1) ) &
4223  -dtt_deta* &
4224  ( min(vy_c_help, 0.0_dp) &
4225  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
4226  *insq_g22_sgy(j,i) &
4227  +max(vy_c_help, 0.0_dp) &
4228  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
4229  *insq_g22_sgy(j-1,i) )
4230 
4231 #endif
4232 
4233 end do
4234 
4235 kc=kcmax
4236 if (as_perp(j,i) >= zero) then
4237  lgs_a0(kc) = 0.0_dp
4238  lgs_a1(kc) = 1.0_dp
4239  lgs_b(kc) = 0.0_dp
4240 else
4241  lgs_a0(kc) = -max(adv_vert_sg(kc-1), 0.0_dp)
4242  lgs_a1(kc) = 1.0_dp + max(adv_vert_sg(kc-1), 0.0_dp)
4243  ! adv_vert_sg(KCMAX-1) >= 0 (directed upward)
4244  ! assumed/enforced
4245 #if (ADV_HOR==2)
4246 
4247  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
4248  -dtt_2dxi* &
4249  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
4250  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
4251  *insq_g11_sgx(j,i) &
4252  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
4253  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
4254  *insq_g11_sgx(j,i-1) ) &
4255  -dtt_2deta* &
4256  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
4257  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
4258  *insq_g22_sgy(j,i) &
4259  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
4260  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
4261  *insq_g22_sgy(j-1,i) )
4262 
4263 #elif (ADV_HOR==3)
4264 
4265  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
4266  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
4267 
4268  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
4269  -dtt_dxi* &
4270  ( min(vx_c_help, 0.0_dp) &
4271  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
4272  *insq_g11_sgx(j,i) &
4273  +max(vx_c_help, 0.0_dp) &
4274  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
4275  *insq_g11_sgx(j,i-1) ) &
4276  -dtt_deta* &
4277  ( min(vy_c_help, 0.0_dp) &
4278  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
4279  *insq_g22_sgy(j,i) &
4280  +max(vy_c_help, 0.0_dp) &
4281  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
4282  *insq_g22_sgy(j-1,i) )
4283 
4284 #endif
4285 
4286 end if
4287 
4288 !-------- Solve system of linear equations --------
4289 
4290 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
4291 
4292 !-------- Assign the result,
4293 ! restriction to interval [0, AGE_MAX yr] --------
4294 
4295 do kc=0, kcmax
4296 
4297  age_c_neu(kc,j,i) = lgs_x(kc)
4298 
4299  if (age_c_neu(kc,j,i) < (age_min*year_sec)) &
4300  age_c_neu(kc,j,i) = 0.0_dp
4301  if (age_c_neu(kc,j,i) > (age_max*year_sec)) &
4302  age_c_neu(kc,j,i) = age_max*year_sec
4303 
4304 end do
4305 
4306 !-------- Age of the ice in the non-existing temperate layer --------
4307 
4308 do kt=0, ktmax
4309  age_t_neu(kt,j,i) = age_c_neu(0,j,i)
4310 end do
4311 
4312 end subroutine calc_temp_ssa
4313 
4314 !-------------------------------------------------------------------------------
4315 
4316 end module calc_temp_m
4317 !
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:448
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:456
sico_variables_m::kappa_r
real(dp) kappa_r
KAPPA_R: Heat conductivity of the lithosphere.
Definition: sico_variables_m.F90:592
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:584
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:520
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:547
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:462
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:562
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:317
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:292
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:503
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:474
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:518
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:505
sico_variables_m::h_r
real(dp) h_r
H_R: Thickness of the modelled lithosphere layer.
Definition: sico_variables_m.F90:588
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:3461
sico_variables_m::eps
real(dp), parameter eps
eps: Small number
Definition: sico_variables_m.F90:656
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:586
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:1928
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:355
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:501
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:477
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:397
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:380
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:495
sico_variables_m::g
real(dp) g
G: Acceleration due to gravity.
Definition: sico_variables_m.F90:577
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:325
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:2447
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:464
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:590
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:515
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:458
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:294
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:499
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:460
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:579
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:454
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:395
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:3564
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:492
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:382
sico_variables_m::beta
real(dp) beta
BETA: Clausius-Clapeyron gradient of ice.
Definition: sico_variables_m.F90:581
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:451
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:497
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:575
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:389
sico_variables_m::errormsg
character(len=256) errormsg
errormsg: Error-message string
Definition: sico_variables_m.F90:832
sico_variables_m::rho
real(dp) rho
RHO: Density of ice.
Definition: sico_variables_m.F90:569
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:472
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:529
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:385
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:489
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 Wed Oct 31 2018 17:45:43 for SICOPOLIS V5-dev by   doxygen 1.8.11