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