SICOPOLIS V5-dev  Revision 1288
calc_temp_enth_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : c a l c _ t e m p _ e n t h _ m
4 !
5 !> @file
6 !!
7 !! Computation of temperature, water content and age with the enthalpy method.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2013-2018 Ralf Greve, Heinz Blatter
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 with the enthalpy method.
34 !<------------------------------------------------------------------------------
36 
37  use sico_types_m
39  use sico_vars_m
40  use error_m
41 
42  implicit none
43 
44  private
45  public :: calc_temp_enth
46 
47 contains
48 
49 !-------------------------------------------------------------------------------
50 !> Main subroutine of calc_temp_enth_m:
51 !! Computation of temperature, water content and age with the enthalpy method.
52 !<------------------------------------------------------------------------------
53 subroutine calc_temp_enth(dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
54  dtime_temp)
55 
56 implicit none
57 
58 real(dp), intent(in) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
59 real(dp), intent(in) :: dtime_temp
60 
61 integer(i4b) :: i, j, kc, kt, kr, ii, jj
62 real(dp) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
63  at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
64  at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
65  ai1(0:kcmax), ai2(0:kcmax), &
66  atr1, acb1, acb2, acb3, acb4, alb1, aqtlde(0:kcmax), &
67  am1, am3(0:kcmax)
68 real(dp) :: dtt_2dxi, dtt_2deta
69 
70 !-------- Term abbreviations
71 
72 at7 = 2.0_dp/rho*dtime_temp
73 
74 atr1 = kappa_r/(rho_c_r*h_r**2)*dtime_temp/(dzeta_r**2)
75 
76 if (flag_aa_nonzero) then
77  am1 = aa*beta*dzeta_c/(ea-1.0_dp)
78 else
79  am1 = beta*dzeta_c
80 end if
81 
82 if (flag_aa_nonzero) then
83  acb1 = (ea-1.0_dp)/aa/dzeta_c
84 else
85  acb1 = 1.0_dp/dzeta_c
86 end if
87 
88 acb2 = kappa_r/h_r/dzeta_r
89 acb3 = rho*g
90 acb4 = rho*g
91 
92 alb1 = h_r/kappa_r*dzeta_r
93 
94 dtt_2dxi = 0.5_dp*dtime_temp/dxi
95 dtt_2deta = 0.5_dp*dtime_temp/deta
96 
97 do kc=0, kcmax
98 
99  if (flag_aa_nonzero) then
100 
101  at1(kc) = (ea-1.0_dp)/(aa*eaz_c(kc))*dtime_temp/dzeta_c
102  at2_1(kc) = (ea-1.0_dp)/(aa*eaz_c(kc))*dtime_temp/dzeta_c
103  at2_2(kc) = (eaz_c(kc)-1.0_dp)/(aa*eaz_c(kc)) &
104  *dtime_temp/dzeta_c
105  at3_1(kc) = (ea-1.0_dp)/(aa*eaz_c(kc))*dtime_temp/dzeta_c
106  at3_2(kc) = (eaz_c(kc)-1.0_dp)/(aa*eaz_c(kc)) &
107  *dtime_temp/dzeta_c
108  at4_1(kc) = (ea-1.0_dp)/(aa*eaz_c(kc))*dtime_temp/dzeta_c
109  at4_2(kc) = (eaz_c(kc)-1.0_dp)/(aa*eaz_c(kc)) &
110  *dtime_temp/dzeta_c
111  at5(kc) = (ea-1.0_dp)/(rho*aa*eaz_c(kc)) &
112  *dtime_temp/dzeta_c
113  if (kc /= kcmax) then
114  at6(kc) = (ea-1.0_dp) &
115  /(aa*exp(aa*0.5_dp*(zeta_c(kc)+zeta_c(kc+1)))) &
116  /dzeta_c
117  else
118  at6(kc) = 0.0_dp
119  end if
120  ai1(kc) = agediff*(ea-1.0_dp)/(aa*eaz_c(kc)) &
121  *dtime_temp/dzeta_c
122  if (kc /= kcmax) then
123  ai2(kc) = (ea-1.0_dp) &
124  /(aa*exp(aa*0.5_dp*(zeta_c(kc)+zeta_c(kc+1)))) &
125  /dzeta_c
126  else
127  ai2(kc) = 0.0_dp
128  end if
129  aqtlde(kc) = (aa*eaz_c(kc))/(ea-1.0_dp)*dzeta_c/dtime_temp
130  am3(kc) = (aa*eaz_c(kc))/(ea-1.0_dp)*dzeta_c*beta
131 
132  else
133 
134  at1(kc) = dtime_temp/dzeta_c
135  at2_1(kc) = dtime_temp/dzeta_c
136  at2_2(kc) = zeta_c(kc) &
137  *dtime_temp/dzeta_c
138  at3_1(kc) = dtime_temp/dzeta_c
139  at3_2(kc) = zeta_c(kc) &
140  *dtime_temp/dzeta_c
141  at4_1(kc) = dtime_temp/dzeta_c
142  at4_2(kc) = zeta_c(kc) &
143  *dtime_temp/dzeta_c
144  at5(kc) = 1.0_dp/rho &
145  *dtime_temp/dzeta_c
146  if (kc /= kcmax) then
147  at6(kc) = 1.0_dp &
148  /dzeta_c
149  else
150  at6(kc) = 0.0_dp
151  end if
152  ai1(kc) = agediff &
153  *dtime_temp/dzeta_c
154  if (kc /= kcmax) then
155  ai2(kc) = 1.0_dp &
156  /dzeta_c
157  else
158  ai2(kc) = 0.0_dp
159  end if
160  aqtlde(kc) = dzeta_c/dtime_temp
161  am3(kc) = dzeta_c*beta
162 
163  end if
164 
165 end do
166 
167 !-------- Computation loop --------
168 
169 do i=1, imax-1 ! skipping domain margins
170 do j=1, jmax-1 ! skipping domain margins
171 
172  if (maske(j,i)==0) then ! glaciated land
173 
174 ! ------ Old vertical column cold
175 
176  if (n_cts(j,i) == -1) then
177 
178  n_cts_neu(j,i) = -1
179  kc_cts_neu(j,i) = 0
180  zm_neu(j,i) = zb(j,i)
181  h_c_neu(j,i) = h_c(j,i)
182  h_t_neu(j,i) = 0.0_dp
183 
184  call calc_temp_enth_1(at1, at2_1, at2_2, at3_1, at3_2, &
185  at4_1, at4_2, at5, at6, at7, &
186  atr1, acb1, acb2, acb3, acb4, alb1, &
187  ai1, ai2, &
188  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
189 
190 ! ---- Check whether base has become temperate
191 
192  if (temp_c_neu(0,j,i) > temp_c_m(0,j,i)-eps) then
193 
194  n_cts_neu(j,i) = 0
195  kc_cts_neu(j,i) = 0
196 
197  call calc_temp_enth_2(at1, at2_1, at2_2, at3_1, at3_2, &
198  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
199  ai1, ai2, aqtlde, am3, &
200  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
201 
202  end if
203 
204 ! ------ Old vertical column with temperate base
205 
206  else if (n_cts(j,i) == 0) then
207 
208  n_cts_neu(j,i) = 0
209  kc_cts_neu(j,i) = kc_cts(j,i)
210  zm_neu(j,i) = zb(j,i)
211  h_c_neu(j,i) = h_c(j,i)
212  h_t_neu(j,i) = h_t(j,i)
213 
214  call calc_temp_enth_2(at1, at2_1, at2_2, at3_1, at3_2, &
215  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
216  ai1, ai2, aqtlde, am3, &
217  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
218 
219 ! ---- Check whether temperate base becomes cold
220 
221  if ( (temp_c_neu(1,j,i)-temp_c_neu(0,j,i)) < (am1*h_c(j,i)) ) then
222 
223  n_cts_neu(j,i) = -1
224  kc_cts_neu(j,i) = 0
225 
226  call calc_temp_enth_1(at1, at2_1, at2_2, at3_1, at3_2, &
227  at4_1, at4_2, at5, at6, at7, &
228  atr1, acb1, acb2, acb3, acb4, alb1, &
229  ai1, ai2, &
230  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
231 
232  if (temp_c_neu(0,j,i) > temp_c_m(0,j,i)-eps) then
233 
234  n_cts_neu(j,i) = 0
235  kc_cts_neu(j,i) = 0
236 
237  call calc_temp_enth_2(at1, at2_1, at2_2, at3_1, at3_2, &
238  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
239  ai1, ai2, aqtlde, am3, &
240  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
241 
242  end if
243 
244  end if
245 
246  end if
247 
248 #if (MARGIN==3)
249 
250  else if (maske(j,i)==3) then ! floating ice
251 
252  n_cts_neu(j,i) = -1
253  kc_cts_neu(j,i) = 0
254  zm_neu(j,i) = zb(j,i)
255  h_c_neu(j,i) = h_c(j,i)
256  h_t_neu(j,i) = 0.0_dp
257 
258  call calc_temp_enth_ssa(at1, at2_1, at2_2, at3_1, at3_2, &
259  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
260  ai1, ai2, &
261  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
262 
263 ! ------ Reset temperatures above melting to the melting point
264 ! and water contents above zero to zero
265 ! (should not occur, but just in case)
266 
267  do kc=0, kcmax
268  if (temp_c_neu(kc,j,i) > temp_c_m(kc,j,i)) &
269  temp_c_neu(kc,j,i) = temp_c_m(kc,j,i)
270  if (omega_c_neu(kc,j,i) > 0.0_dp) &
271  omega_c_neu(kc,j,i) = 0.0_dp
272  end do
273 
274 #endif
275 
276  else ! maske(j,i) == 1,2 (ice-free land or sea point)
277 
278  n_cts_neu(j,i) = -1
279  kc_cts_neu(j,i) = 0
280  zm_neu(j,i) = zb(j,i)
281  h_c_neu(j,i) = h_c(j,i)
282  h_t_neu(j,i) = 0.0_dp
283 
284  call calc_temp_enth_r(atr1, alb1, i, j)
285 
286  end if
287 
288 end do
289 end do
290 
291 !-------- Extrapolate values on margins --------
292 
293 ! ------ Lower left corner
294 
295 i=0
296 j=0
297 
298 if ( (maske(j,i) == 0).or.(maske(j,i) == 3) ) then
299  ! glaciated land or floating ice
300  ii=i+1
301  jj=j+1
302 
303  do kc=0, kcmax
304  enth_c_neu(kc,j,i) = enth_c_neu(kc,jj,ii)
305  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii)
306  omega_c_neu(kc,j,i) = omega_c_neu(kc,jj,ii)
307  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii)
308  end do
309 
310  do kt=0, ktmax
311  ! redundant, lower (kt) ice layer
312  enth_t_neu(kt,j,i) = enth_t_neu(kt,jj,ii)
313  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii)
314  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii)
315  end do
316 
317  do kr=0, krmax
318  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii)
319  end do
320 
321  n_cts_neu(j,i) = n_cts_neu(jj,ii)
322  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
323  zm_neu(j,i) = zb(j,i)
324  h_c_neu(j,i) = h_c(j,i)
325  h_t_neu(j,i) = h_t(j,i)
326 
327 else ! maske(j,i) == 1,2 (ice-free land or sea point)
328 
329  n_cts_neu(j,i) = -1
330  kc_cts_neu(j,i) = 0
331  zm_neu(j,i) = zb(j,i)
332  h_c_neu(j,i) = h_c(j,i)
333  h_t_neu(j,i) = h_t(j,i)
334 
335  call calc_temp_enth_r(atr1, alb1, i, j)
336 
337 end if
338 
339 ! ------ Lower right corner
340 
341 i=imax
342 j=0
343 
344 if ( (maske(j,i) == 0).or.(maske(j,i) == 3) ) then
345  ! glaciated land or floating ice
346  ii=i-1
347  jj=j+1
348 
349  do kc=0, kcmax
350  enth_c_neu(kc,j,i) = enth_c_neu(kc,jj,ii)
351  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii)
352  omega_c_neu(kc,j,i) = omega_c_neu(kc,jj,ii)
353  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii)
354  end do
355 
356  do kt=0, ktmax
357  ! redundant, lower (kt) ice layer
358  enth_t_neu(kt,j,i) = enth_t_neu(kt,jj,ii)
359  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii)
360  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii)
361  end do
362 
363  do kr=0, krmax
364  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii)
365  end do
366 
367  n_cts_neu(j,i) = n_cts_neu(jj,ii)
368  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
369  zm_neu(j,i) = zb(j,i)
370  h_c_neu(j,i) = h_c(j,i)
371  h_t_neu(j,i) = h_t(j,i)
372 
373 else ! maske(j,i) == 1,2 (ice-free land or sea point)
374 
375  n_cts_neu(j,i) = -1
376  kc_cts_neu(j,i) = 0
377  zm_neu(j,i) = zb(j,i)
378  h_c_neu(j,i) = h_c(j,i)
379  h_t_neu(j,i) = h_t(j,i)
380 
381  call calc_temp_enth_r(atr1, alb1, i, j)
382 
383 end if
384 
385 ! ------ Upper left corner
386 
387 i=0
388 j=jmax
389 
390 if ( (maske(j,i) == 0).or.(maske(j,i) == 3) ) then
391  ! glaciated land or floating ice
392  ii=i+1
393  jj=j-1
394 
395  do kc=0, kcmax
396  enth_c_neu(kc,j,i) = enth_c_neu(kc,jj,ii)
397  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii)
398  omega_c_neu(kc,j,i) = omega_c_neu(kc,jj,ii)
399  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii)
400  end do
401 
402  do kt=0, ktmax
403  ! redundant, lower (kt) ice layer
404  enth_t_neu(kt,j,i) = enth_t_neu(kt,jj,ii)
405  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii)
406  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii)
407  end do
408 
409  do kr=0, krmax
410  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii)
411  end do
412 
413  n_cts_neu(j,i) = n_cts_neu(jj,ii)
414  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
415  zm_neu(j,i) = zb(j,i)
416  h_c_neu(j,i) = h_c(j,i)
417  h_t_neu(j,i) = h_t(j,i)
418 
419 else ! maske(j,i) == 1,2 (ice-free land or sea point)
420 
421  n_cts_neu(j,i) = -1
422  kc_cts_neu(j,i) = 0
423  zm_neu(j,i) = zb(j,i)
424  h_c_neu(j,i) = h_c(j,i)
425  h_t_neu(j,i) = h_t(j,i)
426 
427  call calc_temp_enth_r(atr1, alb1, i, j)
428 
429 end if
430 
431 ! ------ Upper right corner
432 
433 i=imax
434 j=jmax
435 
436 if ( (maske(j,i) == 0).or.(maske(j,i) == 3) ) then
437  ! glaciated land or floating ice
438  ii=i-1
439  jj=j-1
440 
441  do kc=0, kcmax
442  enth_c_neu(kc,j,i) = enth_c_neu(kc,jj,ii)
443  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii)
444  omega_c_neu(kc,j,i) = omega_c_neu(kc,jj,ii)
445  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii)
446  end do
447 
448  do kt=0, ktmax
449  ! redundant, lower (kt) ice layer
450  enth_t_neu(kt,j,i) = enth_t_neu(kt,jj,ii)
451  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii)
452  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii)
453  end do
454 
455  do kr=0, krmax
456  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii)
457  end do
458 
459  n_cts_neu(j,i) = n_cts_neu(jj,ii)
460  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
461  zm_neu(j,i) = zb(j,i)
462  h_c_neu(j,i) = h_c(j,i)
463  h_t_neu(j,i) = h_t(j,i)
464 
465 else ! maske(j,i) == 1,2 (ice-free land or sea point)
466 
467  n_cts_neu(j,i) = -1
468  kc_cts_neu(j,i) = 0
469  zm_neu(j,i) = zb(j,i)
470  h_c_neu(j,i) = h_c(j,i)
471  h_t_neu(j,i) = h_t(j,i)
472 
473  call calc_temp_enth_r(atr1, alb1, i, j)
474 
475 end if
476 
477 ! ------ Lower and upper margins
478 
479 do i=1, imax-1
480 
481 ! ---- Lower margin
482 
483  j=0
484 
485  if ( (maske(j,i) == 0).or.(maske(j,i) == 3) ) then
486  ! glaciated land or floating ice
487  ii=i
488  jj=j+1
489 
490  do kc=0, kcmax
491  enth_c_neu(kc,j,i) = enth_c_neu(kc,jj,ii)
492  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii)
493  omega_c_neu(kc,j,i) = omega_c_neu(kc,jj,ii)
494  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii)
495  end do
496 
497  do kt=0, ktmax
498  ! redundant, lower (kt) ice layer
499  enth_t_neu(kt,j,i) = enth_t_neu(kt,jj,ii)
500  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii)
501  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii)
502  end do
503 
504  do kr=0, krmax
505  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii)
506  end do
507 
508  n_cts_neu(j,i) = n_cts_neu(jj,ii)
509  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
510  zm_neu(j,i) = zb(j,i)
511  h_c_neu(j,i) = h_c(j,i)
512  h_t_neu(j,i) = h_t(j,i)
513 
514  else ! maske(j,i) == 1,2 (ice-free land or sea point)
515 
516  n_cts_neu(j,i) = -1
517  kc_cts_neu(j,i) = 0
518  zm_neu(j,i) = zb(j,i)
519  h_c_neu(j,i) = h_c(j,i)
520  h_t_neu(j,i) = h_t(j,i)
521 
522  call calc_temp_enth_r(atr1, alb1, i, j)
523 
524  end if
525 
526 ! ---- Upper margin
527 
528  j=jmax
529 
530  if ( (maske(j,i) == 0).or.(maske(j,i) == 3) ) then
531  ! glaciated land or floating ice
532  ii=i
533  jj=j-1
534 
535  do kc=0, kcmax
536  enth_c_neu(kc,j,i) = enth_c_neu(kc,jj,ii)
537  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii)
538  omega_c_neu(kc,j,i) = omega_c_neu(kc,jj,ii)
539  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii)
540  end do
541 
542  do kt=0, ktmax
543  ! redundant, lower (kt) ice layer
544  enth_t_neu(kt,j,i) = enth_t_neu(kt,jj,ii)
545  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii)
546  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii)
547  end do
548 
549  do kr=0, krmax
550  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii)
551  end do
552 
553  n_cts_neu(j,i) = n_cts_neu(jj,ii)
554  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
555  zm_neu(j,i) = zb(j,i)
556  h_c_neu(j,i) = h_c(j,i)
557  h_t_neu(j,i) = h_t(j,i)
558 
559  else ! maske(j,i) == 1,2 (ice-free land or sea point)
560 
561  n_cts_neu(j,i) = -1
562  kc_cts_neu(j,i) = 0
563  zm_neu(j,i) = zb(j,i)
564  h_c_neu(j,i) = h_c(j,i)
565  h_t_neu(j,i) = h_t(j,i)
566 
567  call calc_temp_enth_r(atr1, alb1, i, j)
568 
569  end if
570 
571 end do
572 
573 ! ------ Left and right margins
574 
575 do j=1, jmax-1
576 
577 ! ---- Left margin
578 
579  i=0
580 
581  if ( (maske(j,i) == 0).or.(maske(j,i) == 3) ) then
582  ! glaciated land or floating ice
583  ii=i+1
584  jj=j
585 
586  do kc=0, kcmax
587  enth_c_neu(kc,j,i) = enth_c_neu(kc,jj,ii)
588  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii)
589  omega_c_neu(kc,j,i) = omega_c_neu(kc,jj,ii)
590  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii)
591  end do
592 
593  do kt=0, ktmax
594  ! redundant, lower (kt) ice layer
595  enth_t_neu(kt,j,i) = enth_t_neu(kt,jj,ii)
596  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii)
597  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii)
598  end do
599 
600  do kr=0, krmax
601  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii)
602  end do
603 
604  n_cts_neu(j,i) = n_cts_neu(jj,ii)
605  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
606  zm_neu(j,i) = zb(j,i)
607  h_c_neu(j,i) = h_c(j,i)
608  h_t_neu(j,i) = h_t(j,i)
609 
610  else ! maske(j,i) == 1,2 (ice-free land or sea point)
611 
612  n_cts_neu(j,i) = -1
613  kc_cts_neu(j,i) = 0
614  zm_neu(j,i) = zb(j,i)
615  h_c_neu(j,i) = h_c(j,i)
616  h_t_neu(j,i) = h_t(j,i)
617 
618  call calc_temp_enth_r(atr1, alb1, i, j)
619 
620  end if
621 
622 ! ---- Right margin
623 
624  i=imax
625 
626  if ( (maske(j,i) == 0).or.(maske(j,i) == 3) ) then
627  ! glaciated land or floating ice
628  ii=i-1
629  jj=j
630 
631  do kc=0, kcmax
632  enth_c_neu(kc,j,i) = enth_c_neu(kc,jj,ii)
633  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii)
634  omega_c_neu(kc,j,i) = omega_c_neu(kc,jj,ii)
635  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii)
636  end do
637 
638  do kt=0, ktmax
639  ! redundant, lower (kt) ice layer
640  enth_t_neu(kt,j,i) = enth_t_neu(kt,jj,ii)
641  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii)
642  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii)
643  end do
644 
645  do kr=0, krmax
646  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii)
647  end do
648 
649  n_cts_neu(j,i) = n_cts_neu(jj,ii)
650  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
651  zm_neu(j,i) = zb(j,i)
652  h_c_neu(j,i) = h_c(j,i)
653  h_t_neu(j,i) = h_t(j,i)
654 
655  else ! maske(j,i) == 1,2 (ice-free land or sea point)
656 
657  n_cts_neu(j,i) = -1
658  kc_cts_neu(j,i) = 0
659  zm_neu(j,i) = zb(j,i)
660  h_c_neu(j,i) = h_c(j,i)
661  h_t_neu(j,i) = h_t(j,i)
662 
663  call calc_temp_enth_r(atr1, alb1, i, j)
664 
665  end if
666 
667 end do
668 
669 end subroutine calc_temp_enth
670 
671 !-------------------------------------------------------------------------------
672 !> Computation of temperature and age for a cold ice column with the
673 !! enthalpy method.
674 !<------------------------------------------------------------------------------
675 subroutine calc_temp_enth_1(at1, at2_1, at2_2, at3_1, at3_2, &
676  at4_1, at4_2, at5, at6, at7, &
677  atr1, acb1, acb2, acb3, acb4, alb1, &
678  ai1, ai2, &
679  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
680 
681 #ifndef ALLOW_OPENAD /* Normal */
682 use sico_maths_m, only : tri_sle
683 #else /* OpenAD */
684 use sico_maths_m
685 #endif /* Normal vs. OpenAD */
686 
688 
689 implicit none
690 
691 #ifndef ALLOW_OPENAD /* Normal */
692 integer(i4b), intent(in) :: i, j
693 #else /* OpenAD */
694 integer(i4b), intent(inout) :: i, j
695 #endif /* Normal vs. OpenAD */
696 
697 real(dp), intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
698  at3_1(0:kcmax), at3_2(0:kcmax), &
699  at4_1(0:kcmax), at4_2(0:kcmax), &
700  at5(0:kcmax), at6(0:kcmax), at7, &
701  ai1(0:kcmax), ai2(0:kcmax), &
702  atr1, acb1, acb2, acb3, acb4, alb1
703 real(dp), intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
704 
705 integer(i4b) :: kc, kt, kr
706 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
707  ce5(0:kcmax), ce6(0:kcmax), ce7(0:kcmax), &
708  ctr1, ccbe1, ccb2, ccb3, ccb4, clb1
709 real(dp) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), ct3_sg(0:kcmax), &
710  ct4_sg(0:kcmax), adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
711 real(dp) :: ci1(0:kcmax), ci2(0:kcmax)
712 real(dp) :: dtt_dxi, dtt_deta
713 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
714  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
715  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
716  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
717  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
718 
719 !-------- Check for boundary points --------
720 
721 if ((i == 0).or.(i == imax).or.(j == 0).or.(j == jmax)) then
722  errormsg = ' >>> calc_temp_enth_1: Boundary points not allowed.'
723  call error(errormsg)
724 end if
725 
726 !-------- Abbreviations --------
727 
728 call calc_temp_enth_1_a(at1, at2_1, at2_2, at3_1, at3_2, &
729  at4_1, at4_2, at5, at6, at7, &
730  atr1, acb1, acb2, acb3, acb4, alb1, &
731  ai1, ai2, &
732  dtime_temp, dtt_2dxi, dtt_2deta, i, j, &
733  ct1, ct2, ct3, ct4, ce5, ce6, ce7, &
734  ctr1, ccbe1, ccb2, ccb3, ccb4, clb1, &
735  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
736  adv_vert_sg, abs_adv_vert_sg, &
737  ci1, ci2, dtt_dxi, dtt_deta)
738 
739 !-------- Computation of the bedrock temperature
740 ! (upper boundary condition: old temperature at the ice base) --------
741 
742 ! ------ Set-up of the the equations
743 
744 call calc_temp_enth_1_b(ctr1, clb1, i, j, &
745  lgs_a0, lgs_a1, lgs_a2, lgs_b)
746 
747 ! ------ Solution of the system of linear equations
748 
749 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
750 
751 ! ------ Assignment of the result (predictor values)
752 
753 do kr=0, krmax
754  temp_r_neu(kr,j,i) = lgs_x(kr)
755 end do
756 
757 !-------- Computation of the ice enthalpy --------
758 
759 ! ------ Set-up of the the equations
760 
761 call calc_temp_enth_1_c(ct1, ct2, ct3, ct4, ce5, ce6, ce7, &
762  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
763  ccbe1, ccb2, ccb3, ccb4, &
764  adv_vert_sg, abs_adv_vert_sg, &
765  dtime_temp, dtt_dxi, dtt_deta, &
766  dtt_2dxi, dtt_2deta, &
767  i, j, &
768  lgs_a0, lgs_a1, lgs_a2, lgs_b)
769 
770 ! ------ Solution of the system of linear equations
771 
772 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
773 
774 ! ------ Assignment of the result
775 
776 do kc=0, kcmax
777  enth_c_neu(kc,j,i) = lgs_x(kc)
778  temp_c_neu(kc,j,i) = temp_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
779  omega_c_neu(kc,j,i) = 0.0_dp ! solution is supposed to be for cold ice
780 end do
781 
782 !-------- Water drainage from the non-existing temperate ice --------
783 
784 q_tld(j,i) = 0.0_dp
785 
786 !-------- Set enthalpy and water content in the redundant,
787 ! lower (kt) ice layer to the value at the ice base --------
788 
789 do kt=0, ktmax
790  enth_t_neu(kt,j,i) = enth_c_neu(0,j,i)
791  omega_t_neu(kt,j,i) = omega_c_neu(0,j,i)
792 end do
793 
794 !-------- Computation of the age of ice --------
795 
796 ! ------ Set-up of the the equations
797 
798 call calc_temp_enth_1_d(ct1, ct2, ct3, ct4, ci1, ci2, &
799  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
800  adv_vert_sg, abs_adv_vert_sg, &
801  dtime_temp, dtt_dxi, dtt_deta, &
802  dtt_2dxi, dtt_2deta, &
803  i, j, &
804  lgs_a0, lgs_a1, lgs_a2, lgs_b)
805 
806 ! ------ Solution of the system of linear equations
807 
808 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
809 
810 ! ------ Assignment of the result,
811 ! restriction to interval [0, AGE_MAX yr] --------
812 
813 do kc=0, kcmax
814 
815  age_c_neu(kc,j,i) = lgs_x(kc)
816 
817  if (age_c_neu(kc,j,i) < (age_min*year_sec)) &
818  age_c_neu(kc,j,i) = 0.0_dp
819  if (age_c_neu(kc,j,i) > (age_max*year_sec)) &
820  age_c_neu(kc,j,i) = age_max*year_sec
821 
822 end do
823 
824 !-------- Age of the ice in the redundant, lower (kt) ice layer --------
825 
826 do kt=0, ktmax
827  age_t_neu(kt,j,i) = age_c_neu(0,j,i)
828 end do
829 
830 end subroutine calc_temp_enth_1
831 
832 !-------------------------------------------------------------------------------
833 !> Computation of temperature and age for a cold ice column with the
834 !! enthalpy method: Abbreviations.
835 !<------------------------------------------------------------------------------
836 subroutine calc_temp_enth_1_a(at1, at2_1, at2_2, at3_1, at3_2, &
837  at4_1, at4_2, at5, at6, at7, &
838  atr1, acb1, acb2, acb3, acb4, alb1, &
839  ai1, ai2, &
840  dtime_temp, dtt_2dxi, dtt_2deta, i, j, &
841  ct1, ct2, ct3, ct4, ce5, ce6, ce7, &
842  ctr1, ccbe1, ccb2, ccb3, ccb4, clb1, &
843  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
844  adv_vert_sg, abs_adv_vert_sg, &
845  ci1, ci2, dtt_dxi, dtt_deta)
849 
850 implicit none
851 
852 #ifndef ALLOW_OPENAD /* Normal */
853 integer(i4b), intent(in) :: i, j
854 #else /* OpenAD */
855 integer(i4b), intent(inout) :: i, j
856 #endif /* Normal vs. OpenAD */
857 
858 real(dp), intent(in) :: at1(0:kcmax), &
859  at2_1(0:kcmax), at2_2(0:kcmax), &
860  at3_1(0:kcmax), at3_2(0:kcmax), &
861  at4_1(0:kcmax), at4_2(0:kcmax), &
862  at5(0:kcmax), at6(0:kcmax), at7, &
863  ai1(0:kcmax), ai2(0:kcmax), &
864  atr1, acb1, acb2, acb3, acb4, alb1
865 real(dp), intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
866 
867 real(dp), intent(out) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), &
868  ct4(0:kcmax), ce5(0:kcmax), ce6(0:kcmax), &
869  ce7(0:kcmax), &
870  ctr1, ccbe1, ccb2, ccb3, ccb4, clb1
871 real(dp), intent(out) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), &
872  ct3_sg(0:kcmax), ct4_sg(0:kcmax), &
873  adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
874 real(dp), intent(out) :: ci1(0:kcmax), ci2(0:kcmax)
875 real(dp), intent(out) :: dtt_dxi, dtt_deta
876 
877 integer(i4b) :: kc
878 real(dp) :: temp_c_help(0:kcmax)
879 
880 !-------- Initialisation --------
881 
882 ct1 = 0.0_dp
883 ct2 = 0.0_dp
884 ct3 = 0.0_dp
885 ct4 = 0.0_dp
886 ce5 = 0.0_dp
887 ce6 = 0.0_dp
888 ce7 = 0.0_dp
889 ctr1 = 0.0_dp
890 clb1 = 0.0_dp
891 ct1_sg = 0.0_dp
892 ct2_sg = 0.0_dp
893 ct3_sg = 0.0_dp
894 ct4_sg = 0.0_dp
895 adv_vert_sg = 0.0_dp
896 abs_adv_vert_sg = 0.0_dp
897 ci1 = 0.0_dp
898 ci2 = 0.0_dp
899 dtt_dxi = 0.0_dp
900 dtt_deta = 0.0_dp
901 
902 !-------- Actual computation --------
903 
904 ctr1 = atr1
905 
906 ccbe1 = acb1 &
907  *kappa_val(temp_c(0,j,i)) &
908  /(c_val(temp_c(0,j,i))*h_c(j,i))
909 ccb2 = acb2
910 
911 #if (DYNAMICS==2)
912 if (.not.flag_shelfy_stream(j,i)) then
913 #endif
914 
915  ccb3 = acb3*0.5_dp*(vx_t(0,j,i)+vx_t(0,j,i-1)) &
916  *h_c(j,i)*dzs_dxi_g(j,i)
917  ccb4 = acb4*0.5_dp*(vy_t(0,j,i)+vy_t(0,j-1,i)) &
918  *h_c(j,i)*dzs_deta_g(j,i)
919 
920 #if (DYNAMICS==2)
921 else ! flag_shelfy_stream(j,i) == .true.
922 
923  ccb3 = -c_drag(j,i) &
924  * sqrt(vx_b_g(j,i)**2 &
925  +vy_b_g(j,i)**2) &
926  **(1.0_dp+p_weert_inv(j,i))
927  ccb4 = 0.0_dp
928 
929 end if
930 #endif
931 
932 clb1 = alb1*q_geo(j,i)
933 
934 #if (ADV_VERT==1)
935 
936 do kc=1, kcmax-1
937  ct1(kc) = at1(kc)/h_c(j,i)*0.5_dp*(vz_c(kc,j,i)+vz_c(kc-1,j,i))
938 end do
939 
940 kc=0
941 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
942  ! only needed for kc=0 ...
943 kc=kcmax-1
944 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
945  ! ... and kc=KCMAX-1
946 
947 #elif (ADV_VERT==2 || ADV_VERT==3)
948 
949 do kc=0, kcmax-1
950  ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
951 end do
952 
953 #endif
954 
955 do kc=0, kcmax
956 
957  ct2(kc) = ( at2_1(kc)*dzm_dtau(j,i) &
958  +at2_2(kc)*dh_c_dtau(j,i) )/h_c(j,i)
959  ct3(kc) = ( at3_1(kc)*dzm_dxi_g(j,i) &
960  +at3_2(kc)*dh_c_dxi_g(j,i) )/h_c(j,i) &
961  *0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1)) *insq_g11_g(j,i)
962  ct4(kc) = ( at4_1(kc)*dzm_deta_g(j,i) &
963  +at4_2(kc)*dh_c_deta_g(j,i) )/h_c(j,i) &
964  *0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i)) *insq_g22_g(j,i)
965  ce5(kc) = at5(kc)/h_c(j,i)
966 
967 #if (DYNAMICS==2)
968  if (.not.flag_shelfy_stream(j,i)) then
969 #endif
970  ce7(kc) = at7 &
971  *enh_c(kc,j,i) &
972  *ratefac_c_t(temp_c(kc,j,i), omega_c(kc,j,i), temp_c_m(kc,j,i)) &
973  *creep(sigma_c(kc,j,i)) &
974  *sigma_c(kc,j,i)*sigma_c(kc,j,i)
975 #if (DYNAMICS==2)
976  else
977  ce7(kc) = 2.0_dp*at7 &
978  *viscosity(de_c(kc,j,i), &
979  temp_c(kc,j,i), temp_c_m(kc,j,i), 0.0_dp, &
980  enh_c(kc,j,i), 0_i2b) &
981  *de_c(kc,j,i)**2
982  end if
983 #endif
984 
985  ci1(kc) = ai1(kc)/h_c(j,i)
986 
987 end do
988 
989 #if (ADV_VERT==1)
990 
991 kc=0
992 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
993 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
994 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
995 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
996 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! only needed for kc=0 ...
997 kc=kcmax-1
998 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
999 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
1000 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
1001 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
1002 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! ... and kc=KCMAX-1
1003 
1004 #elif (ADV_VERT==2 || ADV_VERT==3)
1005 
1006 do kc=0, kcmax-1
1007  ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
1008  ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
1009  ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
1010  adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
1011  abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
1012 end do
1013 
1014 #endif
1015 
1016 do kc=0, kcmax-1
1017  temp_c_help(kc) = 0.5_dp*(temp_c(kc,j,i)+temp_c(kc+1,j,i))
1018  ce6(kc) = at6(kc) &
1019  *kappa_val(temp_c_help(kc))/(c_val(temp_c_help(kc))*h_c(j,i))
1020  ci2(kc) = ai2(kc)/h_c(j,i)
1021 end do
1022 
1023 #if (ADV_HOR==3)
1024 dtt_dxi = 2.0_dp*dtt_2dxi
1025 dtt_deta = 2.0_dp*dtt_2deta
1026 #endif
1027 
1028 end subroutine calc_temp_enth_1_a
1029 
1030 !-------------------------------------------------------------------------------
1031 !> Computation of temperature and age for a cold ice column with the
1032 !! enthalpy method:
1033 !! Set-up of the equations for the bedrock temperature.
1034 !<------------------------------------------------------------------------------
1035 subroutine calc_temp_enth_1_b(ctr1, clb1, i, j, &
1036  lgs_a0, lgs_a1, lgs_a2, lgs_b)
1038 implicit none
1039 
1040 #ifndef ALLOW_OPENAD /* Normal */
1041 integer(i4b), intent(in) :: i, j
1042 #else /* OpenAD */
1043 integer(i4b), intent(inout) :: i, j
1044 #endif /* Normal vs. OpenAD */
1045 
1046 real(dp), intent(in) :: ctr1, clb1
1047 
1048 real(dp), intent(out) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
1049  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
1050  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
1051  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
1052 
1053 integer(i4b) :: kr
1054 
1055 !-------- Initialisation --------
1056 
1057 lgs_a0 = 0.0_dp
1058 lgs_a1 = 0.0_dp
1059 lgs_a2 = 0.0_dp
1060 lgs_b = 0.0_dp
1061 
1062 !-------- Actual computation --------
1063 
1064 kr=0
1065 lgs_a1(kr) = 1.0_dp
1066 lgs_a2(kr) = -1.0_dp
1067 lgs_b(kr) = clb1
1068 
1069 #if (Q_LITHO==1)
1070 ! (coupled heat-conducting bedrock)
1071 
1072 do kr=1, krmax-1
1073  lgs_a0(kr) = -ctr1
1074  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
1075  lgs_a2(kr) = -ctr1
1076  lgs_b(kr) = temp_r(kr,j,i)
1077 end do
1078 
1079 #elif (Q_LITHO==0)
1080 ! (no coupled heat-conducting bedrock)
1081 
1082 do kr=1, krmax-1
1083  lgs_a0(kr) = 1.0_dp
1084  lgs_a1(kr) = 0.0_dp
1085  lgs_a2(kr) = -1.0_dp
1086  lgs_b(kr) = 2.0_dp*clb1
1087 end do
1088 
1089 #endif
1090 
1091 kr=krmax
1092 lgs_a0(kr) = 0.0_dp
1093 lgs_a1(kr) = 1.0_dp
1094 lgs_b(kr) = temp_c(0,j,i)
1095 
1096 end subroutine calc_temp_enth_1_b
1097 
1098 !-------------------------------------------------------------------------------
1099 !> Computation of temperature and age for a cold ice column with the
1100 !! enthalpy method:
1101 !! Set-up of the equations for the ice enthalpy.
1102 !<------------------------------------------------------------------------------
1103 subroutine calc_temp_enth_1_c(ct1, ct2, ct3, ct4, ce5, ce6, ce7, &
1104  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
1105  ccbe1, ccb2, ccb3, ccb4, &
1106  adv_vert_sg, abs_adv_vert_sg, &
1107  dtime_temp, dtt_dxi, dtt_deta, &
1108  dtt_2dxi, dtt_2deta, &
1109  i, j, &
1110  lgs_a0, lgs_a1, lgs_a2, lgs_b)
1111 
1113 
1114 implicit none
1115 
1116 #ifndef ALLOW_OPENAD /* Normal */
1117 integer(i4b), intent(in) :: i, j
1118 #else /* OpenAD */
1119 integer(i4b), intent(inout) :: i, j
1120 #endif /* Normal vs. OpenAD */
1121 
1122 real(dp), intent(in) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), &
1123  ct4(0:kcmax), ce5(0:kcmax), ce6(0:kcmax), &
1124  ce7(0:kcmax)
1125 real(dp), intent(in) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), &
1126  ct3_sg(0:kcmax), ct4_sg(0:kcmax), &
1127  ccbe1, ccb2, ccb3, ccb4, &
1128  adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
1129 real(dp), intent(in) :: dtime_temp, dtt_dxi, dtt_deta, dtt_2dxi, dtt_2deta
1130 
1131 real(dp), intent(out) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
1132  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
1133  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
1134  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
1135 
1136 integer(i4b) :: kc, kr
1137 real(dp) :: vx_c_help, vy_c_help
1138 real(dp) :: adv_vert_help
1139 
1140 !-------- Initialisation --------
1141 
1142 lgs_a0 = 0.0_dp
1143 lgs_a1 = 0.0_dp
1144 lgs_a2 = 0.0_dp
1145 lgs_b = 0.0_dp
1146 
1147 !-------- Actual computation --------
1148 
1149 kr=krmax
1150 kc=0
1151 lgs_a1(kc) = -ccbe1
1152 lgs_a2(kc) = ccbe1
1153 lgs_b(kc) = ccb2*(temp_r_neu(kr,j,i)-temp_r_neu(kr-1,j,i)) + ccb3 + ccb4
1154 
1155 do kc=1, kcmax-1
1156 
1157 #if (ADV_VERT==1)
1158 
1159  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
1160  -ce5(kc)*ce6(kc-1)
1161  lgs_a1(kc) = 1.0_dp+ce5(kc)*(ce6(kc)+ce6(kc-1))
1162  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
1163  -ce5(kc)*ce6(kc)
1164 
1165 #elif (ADV_VERT==2)
1166 
1167  lgs_a0(kc) &
1168  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
1169  -ce5(kc)*ce6(kc-1)
1170  lgs_a1(kc) &
1171  = 1.0_dp &
1172  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
1173  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
1174  +ce5(kc)*(ce6(kc)+ce6(kc-1))
1175  lgs_a2(kc) &
1176  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
1177  -ce5(kc)*ce6(kc)
1178 
1179 #elif (ADV_VERT==3)
1180 
1181  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
1182 
1183  lgs_a0(kc) &
1184  = -max(adv_vert_help, 0.0_dp) &
1185  -ce5(kc)*ce6(kc-1)
1186  lgs_a1(kc) &
1187  = 1.0_dp &
1188  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
1189  +ce5(kc)*(ce6(kc)+ce6(kc-1))
1190  lgs_a2(kc) &
1191  = min(adv_vert_help, 0.0_dp) &
1192  -ce5(kc)*ce6(kc)
1193 
1194 #endif
1195 
1196 #if (ADV_HOR==2)
1197 
1198  lgs_b(kc) = enth_c(kc,j,i) + ce7(kc) &
1199  -dtt_2dxi* &
1200  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
1201  *(enth_c(kc,j,i+1)-enth_c(kc,j,i)) &
1202  *insq_g11_sgx(j,i) &
1203  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
1204  *(enth_c(kc,j,i)-enth_c(kc,j,i-1)) &
1205  *insq_g11_sgx(j,i-1) ) &
1206  -dtt_2deta* &
1207  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
1208  *(enth_c(kc,j+1,i)-enth_c(kc,j,i)) &
1209  *insq_g22_sgy(j,i) &
1210  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
1211  *(enth_c(kc,j,i)-enth_c(kc,j-1,i)) &
1212  *insq_g22_sgy(j-1,i) )
1213 
1214 #elif (ADV_HOR==3)
1215 
1216  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
1217  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
1218 
1219  lgs_b(kc) = enth_c(kc,j,i) + ce7(kc) &
1220  -dtt_dxi* &
1221  ( min(vx_c_help, 0.0_dp) &
1222  *(enth_c(kc,j,i+1)-enth_c(kc,j,i)) &
1223  *insq_g11_sgx(j,i) &
1224  +max(vx_c_help, 0.0_dp) &
1225  *(enth_c(kc,j,i)-enth_c(kc,j,i-1)) &
1226  *insq_g11_sgx(j,i-1) ) &
1227  -dtt_deta* &
1228  ( min(vy_c_help, 0.0_dp) &
1229  *(enth_c(kc,j+1,i)-enth_c(kc,j,i)) &
1230  *insq_g22_sgy(j,i) &
1231  +max(vy_c_help, 0.0_dp) &
1232  *(enth_c(kc,j,i)-enth_c(kc,j-1,i)) &
1233  *insq_g22_sgy(j-1,i) )
1234 
1235 #endif
1236 
1237 end do
1238 
1239 kc=kcmax
1240 lgs_a0(kc) = 0.0_dp
1241 lgs_a1(kc) = 1.0_dp
1242 lgs_b(kc) = enth_fct_temp_omega(temp_s(j,i), 0.0_dp)
1243  ! zero water content at the ice surface
1244 
1245 end subroutine calc_temp_enth_1_c
1246 
1247 !-------------------------------------------------------------------------------
1248 !> Computation of temperature and age for a cold ice column with the
1249 !! enthalpy method:
1250 !! Set-up of the equations for the age of ice.
1251 !<------------------------------------------------------------------------------
1252 subroutine calc_temp_enth_1_d(ct1, ct2, ct3, ct4, ci1, ci2, &
1253  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
1254  adv_vert_sg, abs_adv_vert_sg, &
1255  dtime_temp, dtt_dxi, dtt_deta, &
1256  dtt_2dxi, dtt_2deta, &
1257  i, j, &
1258  lgs_a0, lgs_a1, lgs_a2, lgs_b)
1260 implicit none
1261 
1262 #ifndef ALLOW_OPENAD /* Normal */
1263 integer(i4b), intent(in) :: i, j
1264 #else /* OpenAD */
1265 integer(i4b), intent(inout) :: i, j
1266 #endif /* Normal vs. OpenAD */
1267 
1268 real(dp), intent(in) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), &
1269  ct4(0:kcmax), ci1(0:kcmax), ci2(0:kcmax)
1270 real(dp), intent(in) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), &
1271  ct3_sg(0:kcmax), ct4_sg(0:kcmax), &
1272  adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
1273 real(dp), intent(in) :: dtime_temp, dtt_dxi, dtt_deta, dtt_2dxi, dtt_2deta
1274 
1275 real(dp), intent(out) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
1276  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
1277  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
1278  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
1279 
1280 integer(i4b) :: kc
1281 real(dp) :: vx_c_help, vy_c_help
1282 real(dp) :: adv_vert_help
1283 
1284 !-------- Initialisation --------
1285 
1286 lgs_a0 = 0.0_dp
1287 lgs_a1 = 0.0_dp
1288 lgs_a2 = 0.0_dp
1289 lgs_b = 0.0_dp
1290 
1291 !-------- Actual computation --------
1292 
1293 kc=0 ! adv_vert_sg(0) <= 0
1294 lgs_a1(kc) = 1.0_dp - min(adv_vert_sg(kc), 0.0_dp) ! (directed downward)
1295 lgs_a2(kc) = min(adv_vert_sg(kc), 0.0_dp) ! assumed/enforced
1296 
1297 #if (ADV_HOR==2)
1298 
1299 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1300  -dtt_2dxi* &
1301  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
1302  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1303  *insq_g11_sgx(j,i) &
1304  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
1305  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1306  *insq_g11_sgx(j,i-1) ) &
1307  -dtt_2deta* &
1308  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
1309  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1310  *insq_g22_sgy(j,i) &
1311  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
1312  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1313  *insq_g22_sgy(j-1,i) )
1314 
1315 #elif (ADV_HOR==3)
1316 
1317 vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
1318 vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
1319 
1320 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1321  -dtt_dxi* &
1322  ( min(vx_c_help, 0.0_dp) &
1323  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1324  *insq_g11_sgx(j,i) &
1325  +max(vx_c_help, 0.0_dp) &
1326  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1327  *insq_g11_sgx(j,i-1) ) &
1328  -dtt_deta* &
1329  ( min(vy_c_help, 0.0_dp) &
1330  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1331  *insq_g22_sgy(j,i) &
1332  +max(vy_c_help, 0.0_dp) &
1333  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1334  *insq_g22_sgy(j-1,i) )
1335 
1336 #endif
1337 
1338 do kc=1, kcmax-1
1339 
1340 #if (ADV_VERT==1)
1341 
1342  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
1343  -ci1(kc)*ci2(kc-1)
1344  lgs_a1(kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
1345  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
1346  -ci1(kc)*ci2(kc)
1347 
1348 #elif (ADV_VERT==2)
1349 
1350  lgs_a0(kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
1351  lgs_a1(kc) = 1.0_dp &
1352  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
1353  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
1354  lgs_a2(kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
1355 
1356 #elif (ADV_VERT==3)
1357 
1358  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
1359 
1360  lgs_a0(kc) = -max(adv_vert_help, 0.0_dp)
1361  lgs_a1(kc) = 1.0_dp &
1362  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp)
1363  lgs_a2(kc) = min(adv_vert_help, 0.0_dp)
1364 
1365 #endif
1366 
1367 #if (ADV_HOR==2)
1368 
1369  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1370  -dtt_2dxi* &
1371  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
1372  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1373  *insq_g11_sgx(j,i) &
1374  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
1375  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1376  *insq_g11_sgx(j,i-1) ) &
1377  -dtt_2deta* &
1378  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
1379  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1380  *insq_g22_sgy(j,i) &
1381  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
1382  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1383  *insq_g22_sgy(j-1,i) )
1384 
1385 #elif (ADV_HOR==3)
1386 
1387  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
1388  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
1389 
1390  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1391  -dtt_dxi* &
1392  ( min(vx_c_help, 0.0_dp) &
1393  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1394  *insq_g11_sgx(j,i) &
1395  +max(vx_c_help, 0.0_dp) &
1396  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1397  *insq_g11_sgx(j,i-1) ) &
1398  -dtt_deta* &
1399  ( min(vy_c_help, 0.0_dp) &
1400  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1401  *insq_g22_sgy(j,i) &
1402  +max(vy_c_help, 0.0_dp) &
1403  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1404  *insq_g22_sgy(j-1,i) )
1405 
1406 #endif
1407 
1408 end do
1409 
1410 kc=kcmax
1411 if (as_perp(j,i) >= 0.0_dp) then
1412  lgs_a0(kc) = 0.0_dp
1413  lgs_a1(kc) = 1.0_dp
1414  lgs_b(kc) = 0.0_dp
1415 else
1416  lgs_a0(kc) = -max(adv_vert_sg(kc-1), 0.0_dp)
1417  lgs_a1(kc) = 1.0_dp + max(adv_vert_sg(kc-1), 0.0_dp)
1418  ! adv_vert_sg(KCMAX-1) >= 0 (directed upward)
1419  ! assumed/enforced
1420 #if (ADV_HOR==2)
1421 
1422  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1423  -dtt_2dxi* &
1424  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
1425  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1426  *insq_g11_sgx(j,i) &
1427  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
1428  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1429  *insq_g11_sgx(j,i-1) ) &
1430  -dtt_2deta* &
1431  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
1432  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1433  *insq_g22_sgy(j,i) &
1434  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
1435  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1436  *insq_g22_sgy(j-1,i) )
1437 
1438 #elif (ADV_HOR==3)
1439 
1440  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
1441  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
1442 
1443  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1444  -dtt_dxi* &
1445  ( min(vx_c_help, 0.0_dp) &
1446  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1447  *insq_g11_sgx(j,i) &
1448  +max(vx_c_help, 0.0_dp) &
1449  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1450  *insq_g11_sgx(j,i-1) ) &
1451  -dtt_deta* &
1452  ( min(vy_c_help, 0.0_dp) &
1453  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1454  *insq_g22_sgy(j,i) &
1455  +max(vy_c_help, 0.0_dp) &
1456  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1457  *insq_g22_sgy(j-1,i) )
1458 
1459 #endif
1460 
1461 end if
1462 
1463 end subroutine calc_temp_enth_1_d
1464 
1465 !-------------------------------------------------------------------------------
1466 !> Computation of temperature and age for an ice column with a temperate base
1467 !! with the enthalpy method.
1468 !<------------------------------------------------------------------------------
1469 subroutine calc_temp_enth_2(at1, at2_1, at2_2, at3_1, at3_2, &
1470  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
1471  ai1, ai2, aqtlde, am3, &
1472  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
1473 
1474 #ifndef ALLOW_OPENAD /* Normal */
1475 use sico_maths_m, only : tri_sle
1476 #else /* OpenAD */
1477 use sico_maths_m
1478 #endif /* Normal vs. OpenAD */
1479 
1482 
1483 implicit none
1484 
1485 #ifndef ALLOW_OPENAD /* Normal */
1486 integer(i4b), intent(in) :: i, j
1487 #else /* OpenAD */
1488 integer(i4b), intent(inout) :: i, j
1489 #endif /* Normal vs. OpenAD */
1490 
1491 real(dp), intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
1492  at3_1(0:kcmax), at3_2(0:kcmax), &
1493  at4_1(0:kcmax), at4_2(0:kcmax), &
1494  at5(0:kcmax), at6(0:kcmax), at7, &
1495  ai1(0:kcmax), ai2(0:kcmax), &
1496  atr1, alb1, aqtlde(0:kcmax), am3(0:kcmax)
1497 real(dp), intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
1498 
1499 integer(i4b) :: kc, kt, kr
1500 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
1501  ce5(0:kcmax), ce6(0:kcmax), ce7(0:kcmax), ctr1, clb1
1502 real(dp) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), ct3_sg(0:kcmax), &
1503  ct4_sg(0:kcmax), adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
1504 real(dp) :: ci1(0:kcmax), ci2(0:kcmax)
1505 real(dp) :: cqtlde(0:kcmax), cm3(0:kcmax)
1506 real(dp) :: dtt_dxi, dtt_deta
1507 real(dp) :: temp_c_val(0:kcmax), omega_c_val(0:kcmax)
1508 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
1509  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
1510  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
1511  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
1512  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
1513 
1514 real(dp), parameter :: eps_omega=1.0e-12_dp
1515 
1516 #ifdef ALLOW_OPENAD /* OpenAD */
1517 logical :: kcdone
1518 #endif /* OpenAD */
1519 
1520 !-------- Check for boundary points --------
1521 
1522 if ((i == 0).or.(i == imax).or.(j == 0).or.(j == jmax)) then
1523  errormsg = ' >>> calc_temp_enth_2: Boundary points not allowed.'
1524  call error(errormsg)
1525 end if
1526 
1527 !-------- Abbreviations --------
1528 
1529 call calc_temp_enth_2_a1(at1, at2_1, at2_2, at3_1, at3_2, &
1530  at4_1, at4_2, at5, atr1, alb1, &
1531  ai1, aqtlde, &
1532  dtime_temp, dtt_2dxi, dtt_2deta, i, j, &
1533  ct1, ct2, ct3, ct4, ce5, ctr1, clb1, &
1534  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
1535  adv_vert_sg, abs_adv_vert_sg, &
1536  ci1, cqtlde, dtt_dxi, dtt_deta)
1537 
1538 do kc=0, kcmax
1539  temp_c_val(kc) = temp_c(kc,j,i)
1540  omega_c_val(kc) = omega_c(kc,j,i)
1541 end do
1542 
1543 call calc_temp_enth_2_a2(at6, at7, ai2, am3, temp_c_val, omega_c_val, &
1544  i, j, ce6, ce7, ci2, cm3)
1545 
1546 !-------- Computation of the bedrock temperature --------
1547 
1548 ! ------ Set-up of the the equations
1549 
1550 call calc_temp_enth_2_b(ctr1, clb1, i, j, &
1551  lgs_a0, lgs_a1, lgs_a2, lgs_b)
1552 
1553 ! ------ Solution of the system of linear equations
1554 
1555 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
1556 
1557 ! ------ Assignment of the result
1558 
1559 do kr=0, krmax
1560  temp_r_neu(kr,j,i) = lgs_x(kr)
1561 end do
1562 
1563 !-------- Computation of the ice enthalpy (predictor step) --------
1564 
1565 ! ------ Set-up of the equations
1566 
1567 call calc_temp_enth_2_c(ct1, ct2, ct3, ct4, ce5, ce6, ce7, &
1568  ct1_sg, ct2_sg, ct3_sg, ct4_sg, cm3, &
1569  adv_vert_sg, abs_adv_vert_sg, &
1570  dtime_temp, dtt_dxi, dtt_deta, &
1571  dtt_2dxi, dtt_2deta, &
1572 #ifndef ALLOW_OPENAD /* Normal */
1573  i, j, 0_i2b, &
1574 #else /* OpenAD */
1575  i, j, 0_i4b, &
1576 #endif /* Normal vs. OpenAD */
1577  lgs_a0, lgs_a1, lgs_a2, lgs_b)
1578 
1579 ! ------ Solution of the system of linear equations
1580 
1581 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
1582 
1583 ! ------ Assignment of the result
1584 
1585 do kc=0, kcmax
1586  enth_c_neu(kc,j,i) = lgs_x(kc)
1587  temp_c_neu(kc,j,i) = temp_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
1588  omega_c_neu(kc,j,i) = omega_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
1589 end do
1590 
1591 ! ------ Determination of the CTS
1592 
1593 kc_cts_neu(j,i) = 0
1594 
1595 #ifndef ALLOW_OPENAD /* Normal */
1596 
1597 do kc=1, kcmax-1
1598  if (omega_c_neu(kc,j,i) > eps_omega) then
1599  kc_cts_neu(j,i) = kc
1600  else
1601  exit
1602  end if
1603 end do
1604 
1605 #else /* OpenAD */
1606 
1607 kcdone = .false.
1608 
1609 do kc=1, kcmax-1
1610  if (kcdone.eqv..false.) then
1611  if (omega_c_neu(kc,j,i) > eps_omega) then
1612  kc_cts_neu(j,i) = kc
1613  else
1614  kcdone = .true.
1615  end if
1616  end if
1617 end do
1618 
1619 #endif /* Normal vs. OpenAD */
1620 
1621 !-------- Computation of the ice enthalpy
1622 ! (corrector step for the cold-ice domain only
1623 ! in order to fulfull the transition condition at the CTS) --------
1624 
1625 #if (CALCMOD==3) /* ENTM scheme */
1626 
1627 if (kc_cts_neu(j,i) > 0) then
1628 
1629 ! ------ Update of the abbreviations where needed
1630 
1631  do kc=0, kcmax
1632  temp_c_val(kc) = temp_c_neu(kc,j,i)
1633  omega_c_val(kc) = omega_c_neu(kc,j,i)
1634  end do
1635 
1636  call calc_temp_enth_2_a2(at6, at7, ai2, am3, temp_c_val, omega_c_val, &
1637  i, j, ce6, ce7, ci2, cm3)
1638 
1639 ! ------ Set-up of the equations
1640 
1641  call calc_temp_enth_2_c(ct1, ct2, ct3, ct4, ce5, ce6, ce7, &
1642  ct1_sg, ct2_sg, ct3_sg, ct4_sg, cm3, &
1643  adv_vert_sg, abs_adv_vert_sg, &
1644  dtime_temp, dtt_dxi, dtt_deta, &
1645  dtt_2dxi, dtt_2deta, &
1646  i, j, kc_cts_neu(j,i), &
1647  lgs_a0, lgs_a1, lgs_a2, lgs_b)
1648 
1649 ! ------ Solution of the system of linear equations
1650 
1651  call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
1652 
1653 ! ------ Assignment of the result
1654 
1655  kc=kc_cts_neu(j,i)
1656  enth_c_neu(kc,j,i) = lgs_x(kc)
1657  temp_c_neu(kc,j,i) = temp_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
1658  omega_c_neu(kc,j,i) = omega_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
1659 
1660  do kc=kc_cts_neu(j,i)+1, kcmax
1661  enth_c_neu(kc,j,i) = lgs_x(kc)
1662  temp_c_neu(kc,j,i) = temp_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
1663  omega_c_neu(kc,j,i) = 0.0_dp ! cold-ice domain
1664  end do
1665 
1666 end if
1667 
1668 #elif (CALCMOD==2) /* ENTC scheme */
1669 
1670 !!! continue ! no corrector step
1671 
1672 #else
1673 
1674 errormsg = ' >>> calc_temp_enth_2: CALCMOD must be either 2 or 3!'
1675 call error(errormsg)
1676 
1677 #endif
1678 
1679 !-------- Water drainage from temperate ice (if existing) --------
1680 
1681 q_tld(j,i) = 0.0_dp
1682 
1683 do kc=0, kc_cts_neu(j,i)
1684 
1685  if (omega_c_neu(kc,j,i) > omega_max) then
1686 
1687  q_tld(j,i) = q_tld(j,i) + cqtlde(kc)*(omega_c_neu(kc,j,i)-omega_max)
1688 
1689  omega_c_neu(kc,j,i) = omega_max
1691 
1692  end if
1693 
1694 end do
1695 
1696 !-------- Set enthalpy and water content in the redundant,
1697 ! lower (kt) ice layer to the value at the ice base --------
1698 
1699 do kt=0, ktmax
1700  enth_t_neu(kt,j,i) = enth_c_neu(0,j,i)
1701  omega_t_neu(kt,j,i) = omega_c_neu(0,j,i)
1702 end do
1703 
1704 !-------- Computation of the age of ice --------
1705 
1706 ! ------ Set-up of the the equations
1707 
1708 call calc_temp_enth_2_d(ct1, ct2, ct3, ct4, ci1, ci2, &
1709  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
1710  adv_vert_sg, abs_adv_vert_sg, &
1711  dtime_temp, dtt_dxi, dtt_deta, &
1712  dtt_2dxi, dtt_2deta, &
1713  i, j, &
1714  lgs_a0, lgs_a1, lgs_a2, lgs_b)
1715 
1716 ! ------ Solution of the system of linear equations
1717 
1718 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
1719 
1720 ! ------ Assignment of the result
1721 ! restriction to interval [0, AGE_MAX yr]
1722 
1723 do kc=0, kcmax
1724 
1725  age_c_neu(kc,j,i) = lgs_x(kc)
1726 
1727  if (age_c_neu(kc,j,i) < (age_min*year_sec)) &
1728  age_c_neu(kc,j,i) = 0.0_dp
1729  if (age_c_neu(kc,j,i) > (age_max*year_sec)) &
1730  age_c_neu(kc,j,i) = age_max*year_sec
1731 
1732 end do
1733 
1734 !-------- Age of the ice in the redundant, lower (kt) ice layer --------
1735 
1736 do kt=0, ktmax
1737  age_t_neu(kt,j,i) = age_c_neu(0,j,i)
1738 end do
1739 
1740 end subroutine calc_temp_enth_2
1741 
1742 !-------------------------------------------------------------------------------
1743 !> Computation of temperature and age for an ice column with a temperate base
1744 !! with the enthalpy method: Abbreviations I.
1745 !<------------------------------------------------------------------------------
1746 subroutine calc_temp_enth_2_a1(at1, at2_1, at2_2, at3_1, at3_2, &
1747  at4_1, at4_2, at5, atr1, alb1, &
1748  ai1, aqtlde, &
1749  dtime_temp, dtt_2dxi, dtt_2deta, i, j, &
1750  ct1, ct2, ct3, ct4, ce5, ctr1, clb1, &
1751  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
1752  adv_vert_sg, abs_adv_vert_sg, &
1753  ci1, cqtlde, dtt_dxi, dtt_deta)
1755 implicit none
1756 
1757 #ifndef ALLOW_OPENAD /* Normal */
1758 integer(i4b), intent(in) :: i, j
1759 #else /* OpenAD */
1760 integer(i4b), intent(inout) :: i, j
1761 #endif /* Normal vs. OpenAD */
1762 
1763 real(dp), intent(in) :: at1(0:kcmax), &
1764  at2_1(0:kcmax), at2_2(0:kcmax), &
1765  at3_1(0:kcmax), at3_2(0:kcmax), &
1766  at4_1(0:kcmax), at4_2(0:kcmax), &
1767  at5(0:kcmax), ai1(0:kcmax), &
1768  atr1, alb1, aqtlde(0:kcmax)
1769 real(dp), intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
1770 
1771 real(dp), intent(out) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), &
1772  ct4(0:kcmax), ce5(0:kcmax), &
1773  ctr1, clb1
1774 real(dp), intent(out) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), &
1775  ct3_sg(0:kcmax), ct4_sg(0:kcmax), &
1776  adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
1777 real(dp), intent(out) :: ci1(0:kcmax), cqtlde(0:kcmax)
1778 real(dp), intent(out) :: dtt_dxi, dtt_deta
1779 
1780 integer(i4b) :: kc
1781 
1782 !-------- Initialisation --------
1783 
1784 ct1 = 0.0_dp
1785 ct2 = 0.0_dp
1786 ct3 = 0.0_dp
1787 ct4 = 0.0_dp
1788 ce5 = 0.0_dp
1789 ctr1 = 0.0_dp
1790 clb1 = 0.0_dp
1791 ct1_sg = 0.0_dp
1792 ct2_sg = 0.0_dp
1793 ct3_sg = 0.0_dp
1794 ct4_sg = 0.0_dp
1795 adv_vert_sg = 0.0_dp
1796 abs_adv_vert_sg = 0.0_dp
1797 ci1 = 0.0_dp
1798 cqtlde = 0.0_dp
1799 dtt_dxi = 0.0_dp
1800 dtt_deta = 0.0_dp
1801 
1802 !-------- Actual computation --------
1803 
1804 ctr1 = atr1
1805 clb1 = alb1*q_geo(j,i)
1806 
1807 #if (ADV_VERT==1)
1808 
1809 do kc=1, kcmax-1
1810  ct1(kc) = at1(kc)/h_c(j,i)*0.5_dp*(vz_c(kc,j,i)+vz_c(kc-1,j,i))
1811 end do
1812 
1813 kc=0
1814 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
1815  ! only needed for kc=0 ...
1816 kc=kcmax-1
1817 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
1818  ! ... and kc=KCMAX-1
1819 
1820 #elif (ADV_VERT==2 || ADV_VERT==3)
1821 
1822 do kc=0, kcmax-1
1823  ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
1824 end do
1825 
1826 #endif
1827 
1828 do kc=0, kcmax
1829 
1830  ct2(kc) = ( at2_1(kc)*dzm_dtau(j,i) &
1831  +at2_2(kc)*dh_c_dtau(j,i) )/h_c(j,i)
1832  ct3(kc) = ( at3_1(kc)*dzm_dxi_g(j,i) &
1833  +at3_2(kc)*dh_c_dxi_g(j,i) )/h_c(j,i) &
1834  *0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1)) *insq_g11_g(j,i)
1835  ct4(kc) = ( at4_1(kc)*dzm_deta_g(j,i) &
1836  +at4_2(kc)*dh_c_deta_g(j,i) )/h_c(j,i) &
1837  *0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i)) *insq_g22_g(j,i)
1838  ce5(kc) = at5(kc)/h_c(j,i)
1839  ci1(kc) = ai1(kc)/h_c(j,i)
1840  cqtlde(kc) = aqtlde(kc)*h_c(j,i)
1841 
1842 end do
1843 
1844 #if (ADV_VERT==1)
1845 
1846 kc=0
1847 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
1848 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
1849 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
1850 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
1851 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! only needed for kc=0 ...
1852 kc=kcmax-1
1853 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
1854 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
1855 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
1856 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
1857 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! ... and kc=KCMAX-1
1858 
1859 #elif (ADV_VERT==2 || ADV_VERT==3)
1860 
1861 do kc=0, kcmax-1
1862  ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
1863  ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
1864  ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
1865  adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
1866  abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
1867 end do
1868 
1869 #endif
1870 
1871 #if (ADV_HOR==3)
1872 dtt_dxi = 2.0_dp*dtt_2dxi
1873 dtt_deta = 2.0_dp*dtt_2deta
1874 #endif
1875 
1876 end subroutine calc_temp_enth_2_a1
1877 
1878 !-------------------------------------------------------------------------------
1879 !> Computation of temperature and age for an ice column with a temperate base
1880 !! with the enthalpy method: Abbreviations II.
1881 !<------------------------------------------------------------------------------
1882 subroutine calc_temp_enth_2_a2(at6, at7, ai2, am3, temp_c_val, omega_c_val, &
1883  i, j, ce6, ce7, ci2, cm3)
1884 
1886  creep, viscosity
1887 
1888 implicit none
1889 
1890 #ifndef ALLOW_OPENAD /* Normal */
1891 integer(i4b), intent(in) :: i, j
1892 #else /* OpenAD */
1893 integer(i4b), intent(inout) :: i, j
1894 #endif /* Normal vs. OpenAD */
1895 
1896 real(dp), intent(in) :: at6(0:kcmax), at7, ai2(0:kcmax), am3(0:kcmax)
1897 real(dp), intent(in) :: temp_c_val(0:kcmax), omega_c_val(0:kcmax)
1898 
1899 real(dp), intent(out) :: ce6(0:kcmax), ce7(0:kcmax), ci2(0:kcmax), &
1900  cm3(0:kcmax)
1901 
1902 integer(i4b) :: kc
1903 real(dp) :: temp_c_help(0:kcmax)
1904 
1905 !-------- Initialisation --------
1906 
1907 ce6 = 0.0_dp
1908 ce7 = 0.0_dp
1909 ci2 = 0.0_dp
1910 cm3 = 0.0_dp
1911 
1912 !-------- Actual computation --------
1913 
1914 do kc=0, kcmax
1915 
1916 #if (DYNAMICS==2)
1917  if (.not.flag_shelfy_stream(j,i)) then
1918 #endif
1919  ce7(kc) = at7 &
1920  *enh_c(kc,j,i) &
1921  *ratefac_c_t(temp_c_val(kc), omega_c_val(kc), temp_c_m(kc,j,i)) &
1922  *creep(sigma_c(kc,j,i)) &
1923  *sigma_c(kc,j,i)*sigma_c(kc,j,i)
1924 #if (DYNAMICS==2)
1925  else
1926  ce7(kc) = 2.0_dp*at7 &
1927  *viscosity(de_c(kc,j,i), &
1928  temp_c_val(kc), temp_c_m(kc,j,i), omega_c_val(kc), &
1929  enh_c(kc,j,i), 2_i2b) &
1930  *de_c(kc,j,i)**2
1931  end if
1932 #endif
1933 
1934  cm3(kc) = am3(kc)*h_c(j,i)*c_val(temp_c_val(kc))
1935 
1936 end do
1937 
1938 do kc=0, kc_cts_neu(j,i)-1 ! temperate layer
1939  ce6(kc) = at6(kc) &
1940  *nue/h_c(j,i)
1941  ci2(kc) = ai2(kc)/h_c(j,i)
1942 end do
1943 
1944 do kc=kc_cts_neu(j,i), kcmax-1 ! cold layer
1945  temp_c_help(kc) = 0.5_dp*(temp_c_val(kc)+temp_c_val(kc+1))
1946  ce6(kc) = at6(kc) &
1947  *kappa_val(temp_c_help(kc))/(c_val(temp_c_help(kc))*h_c(j,i))
1948  ci2(kc) = ai2(kc)/h_c(j,i)
1949 end do
1950 
1951 end subroutine calc_temp_enth_2_a2
1952 
1953 !-------------------------------------------------------------------------------
1954 !> Computation of temperature and age for an ice column with a temperate base
1955 !! with the enthalpy method:
1956 !! Set-up of the equations for the bedrock temperature.
1957 !<------------------------------------------------------------------------------
1958 subroutine calc_temp_enth_2_b(ctr1, clb1, i, j, &
1959  lgs_a0, lgs_a1, lgs_a2, lgs_b)
1961 implicit none
1962 
1963 #ifndef ALLOW_OPENAD /* Normal */
1964 integer(i4b), intent(in) :: i, j
1965 #else /* OpenAD */
1966 integer(i4b), intent(inout) :: i, j
1967 #endif /* Normal vs. OpenAD */
1968 
1969 real(dp), intent(in) :: ctr1, clb1
1970 
1971 real(dp), intent(out) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
1972  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
1973  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
1974  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
1975 
1976 integer(i4b) :: kr
1977 
1978 !-------- Initialisation --------
1979 
1980 lgs_a0 = 0.0_dp
1981 lgs_a1 = 0.0_dp
1982 lgs_a2 = 0.0_dp
1983 lgs_b = 0.0_dp
1984 
1985 !-------- Actual computation --------
1986 
1987 kr=0
1988 lgs_a1(kr) = 1.0_dp
1989 lgs_a2(kr) = -1.0_dp
1990 lgs_b(kr) = clb1
1991 
1992 #if (Q_LITHO==1)
1993 ! (coupled heat-conducting bedrock)
1994 
1995 do kr=1, krmax-1
1996  lgs_a0(kr) = - ctr1
1997  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
1998  lgs_a2(kr) = - ctr1
1999  lgs_b(kr) = temp_r(kr,j,i)
2000 end do
2001 
2002 #elif (Q_LITHO==0)
2003 ! (no coupled heat-conducting bedrock)
2004 
2005 do kr=1, krmax-1
2006  lgs_a0(kr) = 1.0_dp
2007  lgs_a1(kr) = 0.0_dp
2008  lgs_a2(kr) = -1.0_dp
2009  lgs_b(kr) = 2.0_dp*clb1
2010 end do
2011 
2012 #endif
2013 
2014 kr=krmax
2015 lgs_a0(kr) = 0.0_dp
2016 lgs_a1(kr) = 1.0_dp
2017 lgs_b(kr) = temp_t_m(0,j,i)
2018 
2019 end subroutine calc_temp_enth_2_b
2020 
2021 !-------------------------------------------------------------------------------
2022 !> Computation of temperature and age for an ice column with a temperate base
2023 !! with the enthalpy method:
2024 !! Set-up of the equations for the ice enthalpy.
2025 !<------------------------------------------------------------------------------
2026 subroutine calc_temp_enth_2_c(ct1, ct2, ct3, ct4, ce5, ce6, ce7, &
2027  ct1_sg, ct2_sg, ct3_sg, ct4_sg, cm3, &
2028  adv_vert_sg, abs_adv_vert_sg, &
2029  dtime_temp, dtt_dxi, dtt_deta, &
2030  dtt_2dxi, dtt_2deta, &
2031  i, j, kcmin, &
2032  lgs_a0, lgs_a1, lgs_a2, lgs_b)
2033 
2035 
2036 implicit none
2037 
2038 #ifndef ALLOW_OPENAD /* Normal */
2039 integer(i4b), intent(in) :: i, j
2040 #else /* OpenAD */
2041 integer(i4b), intent(inout) :: i, j
2042 #endif /* Normal vs. OpenAD */
2043 
2044 integer(i2b), intent(in) :: kcmin
2045 real(dp), intent(in) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), &
2046  ct4(0:kcmax), ce5(0:kcmax), ce6(0:kcmax), &
2047  ce7(0:kcmax)
2048 real(dp), intent(in) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), &
2049  ct3_sg(0:kcmax), ct4_sg(0:kcmax), &
2050  adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
2051 real(dp), intent(in) :: cm3(0:kcmax)
2052 real(dp), intent(in) :: dtime_temp, dtt_dxi, dtt_deta, dtt_2dxi, dtt_2deta
2053 
2054 real(dp), intent(out) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
2055  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
2056  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
2057  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
2058 
2059 integer(i4b) :: kc
2060 real(dp) :: vx_c_help, vy_c_help
2061 real(dp) :: adv_vert_help
2062 
2063 !-------- Initialisation --------
2064 
2065 lgs_a0 = 0.0_dp
2066 lgs_a1 = 0.0_dp
2067 lgs_a2 = 0.0_dp
2068 lgs_b = 0.0_dp
2069 
2070 !-------- Actual computation --------
2071 
2072 if (kcmin == 0) then ! predictor step
2073 
2074  kc=0
2075 
2076  if (kc_cts_neu(j,i) == 0) then ! temperate base without temperate layer
2077 
2078  lgs_a1(kc) = 1.0_dp
2079  lgs_a2(kc) = 0.0_dp
2080  lgs_b(kc) = enth_fct_temp_omega(temp_c_m(kc,j,i), 0.0_dp)
2081 
2082  else ! kc_cts_neu(j,i) > 0, temperate base with temperate layer
2083 
2084  lgs_a1(kc) = 1.0_dp
2085  lgs_a2(kc) = -1.0_dp
2086  lgs_b(kc) = 0.0_dp
2087 
2088  end if
2089 
2090 else ! kcmin > 0, corrector step
2091 
2092  kc=0
2093 
2094  lgs_a1(kc) = 1.0_dp ! dummy
2095  lgs_a2(kc) = 0.0_dp ! setting,
2096  lgs_b(kc) = 0.0_dp ! not needed
2097 
2098  do kc=1, kcmin-1
2099 
2100  lgs_a0(kc) = 0.0_dp ! dummy
2101  lgs_a1(kc) = 1.0_dp ! setting,
2102  lgs_a2(kc) = 0.0_dp ! not
2103  lgs_b(kc) = 0.0_dp ! needed
2104 
2105  end do
2106 
2107  kc=kcmin
2108 
2109  lgs_a0(kc) = 0.0_dp
2110  lgs_a1(kc) = 1.0_dp
2111  lgs_a2(kc) = -1.0_dp
2112  lgs_b(kc) = -cm3(kc)
2113 
2114 end if
2115 
2116 do kc=kcmin+1, kcmax-1
2117 
2118 #if (ADV_VERT==1)
2119 
2120  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2121  -ce5(kc)*ce6(kc-1)
2122  lgs_a1(kc) = 1.0_dp+ce5(kc)*(ce6(kc)+ce6(kc-1))
2123  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2124  -ce5(kc)*ce6(kc)
2125 
2126 #elif (ADV_VERT==2)
2127 
2128  lgs_a0(kc) &
2129  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2130  -ce5(kc)*ce6(kc-1)
2131  lgs_a1(kc) &
2132  = 1.0_dp &
2133  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2134  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
2135  +ce5(kc)*(ce6(kc)+ce6(kc-1))
2136  lgs_a2(kc) &
2137  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
2138  -ce5(kc)*ce6(kc)
2139 
2140 #elif (ADV_VERT==3)
2141 
2142  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
2143 
2144  lgs_a0(kc) &
2145  = -max(adv_vert_help, 0.0_dp) &
2146  -ce5(kc)*ce6(kc-1)
2147  lgs_a1(kc) &
2148  = 1.0_dp &
2149  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
2150  +ce5(kc)*(ce6(kc)+ce6(kc-1))
2151  lgs_a2(kc) &
2152  = min(adv_vert_help, 0.0_dp) &
2153  -ce5(kc)*ce6(kc)
2154 
2155 #endif
2156 
2157 #if (ADV_HOR==2)
2158 
2159  lgs_b(kc) = enth_c(kc,j,i) + ce7(kc) &
2160  -dtt_2dxi* &
2161  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2162  *(enth_c(kc,j,i+1)-enth_c(kc,j,i)) &
2163  *insq_g11_sgx(j,i) &
2164  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2165  *(enth_c(kc,j,i)-enth_c(kc,j,i-1)) &
2166  *insq_g11_sgx(j,i-1) ) &
2167  -dtt_2deta* &
2168  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2169  *(enth_c(kc,j+1,i)-enth_c(kc,j,i)) &
2170  *insq_g22_sgy(j,i) &
2171  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2172  *(enth_c(kc,j,i)-enth_c(kc,j-1,i)) &
2173  *insq_g22_sgy(j-1,i) )
2174 
2175 #elif (ADV_HOR==3)
2176 
2177  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2178  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2179 
2180  lgs_b(kc) = enth_c(kc,j,i) + ce7(kc) &
2181  -dtt_dxi* &
2182  ( min(vx_c_help, 0.0_dp) &
2183  *(enth_c(kc,j,i+1)-enth_c(kc,j,i)) &
2184  *insq_g11_sgx(j,i) &
2185  +max(vx_c_help, 0.0_dp) &
2186  *(enth_c(kc,j,i)-enth_c(kc,j,i-1)) &
2187  *insq_g11_sgx(j,i-1) ) &
2188  -dtt_deta* &
2189  ( min(vy_c_help, 0.0_dp) &
2190  *(enth_c(kc,j+1,i)-enth_c(kc,j,i)) &
2191  *insq_g22_sgy(j,i) &
2192  +max(vy_c_help, 0.0_dp) &
2193  *(enth_c(kc,j,i)-enth_c(kc,j-1,i)) &
2194  *insq_g22_sgy(j-1,i) )
2195 
2196 #endif
2197 
2198 end do
2199 
2200 kc=kcmax
2201 lgs_a0(kc) = 0.0_dp
2202 lgs_a1(kc) = 1.0_dp
2203 lgs_b(kc) = enth_fct_temp_omega(temp_s(j,i), 0.0_dp)
2204  ! zero water content at the ice surface
2205 
2206 end subroutine calc_temp_enth_2_c
2207 
2208 !-------------------------------------------------------------------------------
2209 !> Computation of temperature and age for an ice column with a temperate base
2210 !! with the enthalpy method:
2211 !! Set-up of the equations for the age of ice.
2212 !<------------------------------------------------------------------------------
2213 subroutine calc_temp_enth_2_d(ct1, ct2, ct3, ct4, ci1, ci2, &
2214  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
2215  adv_vert_sg, abs_adv_vert_sg, &
2216  dtime_temp, dtt_dxi, dtt_deta, &
2217  dtt_2dxi, dtt_2deta, &
2218  i, j, &
2219  lgs_a0, lgs_a1, lgs_a2, lgs_b)
2221 implicit none
2222 
2223 #ifndef ALLOW_OPENAD /* Normal */
2224 integer(i4b), intent(in) :: i, j
2225 #else /* OpenAD */
2226 integer(i4b), intent(inout) :: i, j
2227 #endif /* Normal vs. OpenAD */
2228 
2229 real(dp), intent(in) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), &
2230  ct4(0:kcmax), ci1(0:kcmax), ci2(0:kcmax)
2231 real(dp), intent(in) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), &
2232  ct3_sg(0:kcmax), ct4_sg(0:kcmax), &
2233  adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
2234 real(dp), intent(in) :: dtime_temp, dtt_dxi, dtt_deta, dtt_2dxi, dtt_2deta
2235 
2236 real(dp), intent(out) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
2237  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
2238  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
2239  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
2240 
2241 integer(i4b) :: kc
2242 real(dp) :: vx_c_help, vy_c_help
2243 real(dp) :: adv_vert_help
2244 
2245 !-------- Initialisation --------
2246 
2247 lgs_a0 = 0.0_dp
2248 lgs_a1 = 0.0_dp
2249 lgs_a2 = 0.0_dp
2250 lgs_b = 0.0_dp
2251 
2252 !-------- Actual computation --------
2253 
2254 kc=0 ! adv_vert_sg(0) <= 0
2255 lgs_a1(kc) = 1.0_dp - min(adv_vert_sg(kc), 0.0_dp) ! (directed downward)
2256 lgs_a2(kc) = min(adv_vert_sg(kc), 0.0_dp) ! assumed/enforced
2257 
2258 #if (ADV_HOR==2)
2259 
2260 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2261  -dtt_2dxi* &
2262  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2263  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2264  *insq_g11_sgx(j,i) &
2265  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2266  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2267  *insq_g11_sgx(j,i-1) ) &
2268  -dtt_2deta* &
2269  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2270  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2271  *insq_g22_sgy(j,i) &
2272  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2273  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2274  *insq_g22_sgy(j-1,i) )
2275 
2276 #elif (ADV_HOR==3)
2277 
2278 vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2279 vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2280 
2281 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2282  -dtt_dxi* &
2283  ( min(vx_c_help, 0.0_dp) &
2284  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2285  *insq_g11_sgx(j,i) &
2286  +max(vx_c_help, 0.0_dp) &
2287  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2288  *insq_g11_sgx(j,i-1) ) &
2289  -dtt_deta* &
2290  ( min(vy_c_help, 0.0_dp) &
2291  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2292  *insq_g22_sgy(j,i) &
2293  +max(vy_c_help, 0.0_dp) &
2294  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2295  *insq_g22_sgy(j-1,i) )
2296 
2297 #endif
2298 
2299 do kc=1, kcmax-1
2300 
2301 #if (ADV_VERT==1)
2302 
2303  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2304  -ci1(kc)*ci2(kc-1)
2305  lgs_a1(kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
2306  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2307  -ci1(kc)*ci2(kc)
2308 
2309 #elif (ADV_VERT==2)
2310 
2311  lgs_a0(kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
2312  lgs_a1(kc) = 1.0_dp &
2313  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2314  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
2315  lgs_a2(kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
2316 
2317 #elif (ADV_VERT==3)
2318 
2319  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
2320 
2321  lgs_a0(kc) = -max(adv_vert_help, 0.0_dp)
2322  lgs_a1(kc) = 1.0_dp &
2323  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp)
2324  lgs_a2(kc) = min(adv_vert_help, 0.0_dp)
2325 
2326 #endif
2327 
2328 #if (ADV_HOR==2)
2329 
2330  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2331  -dtt_2dxi* &
2332  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2333  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2334  *insq_g11_sgx(j,i) &
2335  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2336  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2337  *insq_g11_sgx(j,i-1) ) &
2338  -dtt_2deta* &
2339  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2340  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2341  *insq_g22_sgy(j,i) &
2342  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2343  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2344  *insq_g22_sgy(j-1,i) )
2345 
2346 #elif (ADV_HOR==3)
2347 
2348  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2349  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2350 
2351  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2352  -dtt_dxi* &
2353  ( min(vx_c_help, 0.0_dp) &
2354  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2355  *insq_g11_sgx(j,i) &
2356  +max(vx_c_help, 0.0_dp) &
2357  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2358  *insq_g11_sgx(j,i-1) ) &
2359  -dtt_deta* &
2360  ( min(vy_c_help, 0.0_dp) &
2361  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2362  *insq_g22_sgy(j,i) &
2363  +max(vy_c_help, 0.0_dp) &
2364  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2365  *insq_g22_sgy(j-1,i) )
2366 
2367 #endif
2368 
2369 end do
2370 
2371 kc=kcmax
2372 if (as_perp(j,i) >= 0.0_dp) then
2373  lgs_a0(kc) = 0.0_dp
2374  lgs_a1(kc) = 1.0_dp
2375  lgs_b(kc) = 0.0_dp
2376 else
2377  lgs_a0(kc) = -max(adv_vert_sg(kc-1), 0.0_dp)
2378  lgs_a1(kc) = 1.0_dp + max(adv_vert_sg(kc-1), 0.0_dp)
2379  ! adv_vert_sg(KCMAX-1) >= 0 (directed upward)
2380  ! assumed/enforced
2381 #if (ADV_HOR==2)
2382 
2383  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2384  -dtt_2dxi* &
2385  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2386  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2387  *insq_g11_sgx(j,i) &
2388  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2389  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2390  *insq_g11_sgx(j,i-1) ) &
2391  -dtt_2deta* &
2392  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2393  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2394  *insq_g22_sgy(j,i) &
2395  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2396  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2397  *insq_g22_sgy(j-1,i) )
2398 
2399 #elif (ADV_HOR==3)
2400 
2401  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2402  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2403 
2404  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2405  -dtt_dxi* &
2406  ( min(vx_c_help, 0.0_dp) &
2407  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2408  *insq_g11_sgx(j,i) &
2409  +max(vx_c_help, 0.0_dp) &
2410  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2411  *insq_g11_sgx(j,i-1) ) &
2412  -dtt_deta* &
2413  ( min(vy_c_help, 0.0_dp) &
2414  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2415  *insq_g22_sgy(j,i) &
2416  +max(vy_c_help, 0.0_dp) &
2417  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2418  *insq_g22_sgy(j-1,i) )
2419 
2420 #endif
2421 
2422 end if
2423 
2424 end subroutine calc_temp_enth_2_d
2425 
2426 !-------------------------------------------------------------------------------
2427 !> Computation of temperature, age, water content and enthalpy for an
2428 !! ice-free column.
2429 !<------------------------------------------------------------------------------
2430 subroutine calc_temp_enth_r(atr1, alb1, i, j)
2431 
2432 #ifndef ALLOW_OPENAD /* Normal */
2433 use sico_maths_m, only : tri_sle
2434 #else /* OpenAD */
2435 use sico_maths_m
2436 #endif /* Normal vs. OpenAD */
2437 
2439 
2440 implicit none
2441 
2442 #ifndef ALLOW_OPENAD /* Normal */
2443 integer(i4b), intent(in) :: i, j
2444 #else /* OpenAD */
2445 integer(i4b), intent(inout) :: i, j
2446 #endif /* Normal vs. OpenAD */
2447 
2448 real(dp), intent(in) :: atr1, alb1
2449 
2450 integer(i4b) :: kc, kt, kr
2451 real(dp) :: ctr1, clb1
2452 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
2453  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
2454  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
2455  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
2456  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
2457 real(dp) :: enth_val
2458 
2459 !-------- Abbreviations --------
2460 
2461 ctr1 = atr1
2462 clb1 = alb1*q_geo(j,i)
2463 
2464 !-------- Set up the equations for the bedrock temperature --------
2465 
2466 kr=0
2467 lgs_a1(kr) = 1.0_dp
2468 lgs_a2(kr) = -1.0_dp
2469 lgs_b(kr) = clb1
2470 
2471 #if (Q_LITHO==1)
2472 ! (coupled heat-conducting bedrock)
2473 
2474 do kr=1, krmax-1
2475  lgs_a0(kr) = - ctr1
2476  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
2477  lgs_a2(kr) = - ctr1
2478  lgs_b(kr) = temp_r(kr,j,i)
2479 end do
2480 
2481 #elif (Q_LITHO==0)
2482 ! (no coupled heat-conducting bedrock)
2483 
2484 do kr=1, krmax-1
2485  lgs_a0(kr) = 1.0_dp
2486  lgs_a1(kr) = 0.0_dp
2487  lgs_a2(kr) = -1.0_dp
2488  lgs_b(kr) = 2.0_dp*clb1
2489 end do
2490 
2491 #endif
2492 
2493 kr=krmax
2494 lgs_a0(kr) = 0.0_dp
2495 lgs_a1(kr) = 1.0_dp
2496 lgs_b(kr) = temp_s(j,i)
2497 
2498 !-------- Solve system of linear equations --------
2499 
2500 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
2501 
2502 !-------- Assign the result --------
2503 
2504 do kr=0, krmax
2505  temp_r_neu(kr,j,i) = lgs_x(kr)
2506 end do
2507 
2508 !-------- Water content, age and enthalpy
2509 ! in the non-existing lower (kt) ice layer --------
2510 
2511 enth_val = enth_fct_temp_omega(temp_s(j,i), 0.0_dp)
2512 
2513 do kt=0, ktmax
2514  omega_t_neu(kt,j,i) = 0.0_dp
2515  age_t_neu(kt,j,i) = 0.0_dp
2516  enth_t_neu(kt,j,i) = enth_val
2517 end do
2518 
2519 !-------- Temperature, age, water content and enthalpy
2520 ! in the non-existing upper (kc) ice layer --------
2521 
2522 do kc=0, kcmax
2523  temp_c_neu(kc,j,i) = temp_s(j,i)
2524  age_c_neu(kc,j,i) = 0.0_dp
2525  omega_c_neu(kc,j,i) = 0.0_dp
2526  enth_c_neu(kc,j,i) = enth_val
2527 end do
2528 
2529 end subroutine calc_temp_enth_r
2530 
2531 !-------------------------------------------------------------------------------
2532 !> Computation of temperature and age for ice shelves (floating ice)
2533 !! with the enthalpy method.
2534 !<------------------------------------------------------------------------------
2535 subroutine calc_temp_enth_ssa(at1, at2_1, at2_2, at3_1, at3_2, &
2536  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
2537  ai1, ai2, &
2538  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
2541 
2542 #ifndef ALLOW_OPENAD /* Normal */
2543 use sico_maths_m, only : tri_sle
2544 #else /* OpenAD */
2545 use sico_maths_m
2546 #endif /* Normal vs. OpenAD */
2547 
2550 
2551 implicit none
2552 
2553 #ifndef ALLOW_OPENAD /* Normal */
2554 integer(i4b), intent(in) :: i, j
2555 #else /* OpenAD */
2556 integer(i4b), intent(inout) :: i, j
2557 #endif /* Normal vs. OpenAD */
2558 
2559 real(dp), intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
2560  at3_1(0:kcmax), at3_2(0:kcmax), &
2561  at4_1(0:kcmax), at4_2(0:kcmax), &
2562  at5(0:kcmax), at6(0:kcmax), at7, &
2563  ai1(0:kcmax), ai2(0:kcmax), &
2564  atr1, alb1
2565 real(dp), intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
2566 
2567 integer(i4b) :: kc, kt, kr
2568 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
2569  ce5(0:kcmax), ce6(0:kcmax), ce7(0:kcmax), ctr1, clb1
2570 real(dp) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), ct3_sg(0:kcmax), &
2571  ct4_sg(0:kcmax), adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
2572 real(dp) :: ci1(0:kcmax), ci2(0:kcmax)
2573 real(dp) :: temp_c_help(0:kcmax)
2574 real(dp) :: vx_c_help, vy_c_help
2575 real(dp) :: adv_vert_help
2576 real(dp) :: dtt_dxi, dtt_deta
2577 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
2578  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
2579  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
2580  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
2581  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
2582 real(dp), parameter :: zero=0.0_dp
2583 
2584 !-------- Check for boundary points --------
2585 
2586 if ((i == 0).or.(i == imax).or.(j == 0).or.(j == jmax)) then
2587  errormsg = ' >>> calc_temp_enth_ssa: Boundary points not allowed.'
2588  call error(errormsg)
2589 end if
2590 
2591 !-------- Abbreviations --------
2592 
2593 ctr1 = atr1
2594 clb1 = alb1*q_geo(j,i)
2595 
2596 #if (ADV_VERT==1)
2597 
2598 do kc=1, kcmax-1
2599  ct1(kc) = at1(kc)/h_c(j,i)*0.5_dp*(vz_c(kc,j,i)+vz_c(kc-1,j,i))
2600 end do
2601 
2602 kc=0
2603 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
2604  ! only needed for kc=0 ...
2605 kc=kcmax-1
2606 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
2607  ! ... and kc=KCMAX-1
2608 
2609 #elif (ADV_VERT==2 || ADV_VERT==3)
2610 
2611 do kc=0, kcmax-1
2612  ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
2613 end do
2614 
2615 #endif
2616 
2617 do kc=0, kcmax
2618 
2619  ct2(kc) = ( at2_1(kc)*dzm_dtau(j,i) &
2620  +at2_2(kc)*dh_c_dtau(j,i) )/h_c(j,i)
2621  ct3(kc) = ( at3_1(kc)*dzm_dxi_g(j,i) &
2622  +at3_2(kc)*dh_c_dxi_g(j,i) )/h_c(j,i) &
2623  *0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1)) *insq_g11_g(j,i)
2624  ct4(kc) = ( at4_1(kc)*dzm_deta_g(j,i) &
2625  +at4_2(kc)*dh_c_deta_g(j,i) )/h_c(j,i) &
2626  *0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i)) *insq_g22_g(j,i)
2627  ce5(kc) = at5(kc)/h_c(j,i)
2628  ce7(kc) = 2.0_dp*at7 &
2629  *viscosity(de_ssa(j,i), &
2630  temp_c(kc,j,i), temp_c_m(kc,j,i), 0.0_dp, &
2631 #ifndef ALLOW_OPENAD /* Normal */
2632  enh_c(kc,j,i), 0_i2b) &
2633 #else /* OpenAD */
2634  enh_c(kc,j,i), 0_i4b) &
2635 #endif /* Normal vs. OpenAD */
2636  *de_ssa(j,i)**2
2637  ci1(kc) = ai1(kc)/h_c(j,i)
2638 
2639 end do
2640 
2641 #if (ADV_VERT==1)
2642 
2643 kc=0
2644 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2645 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2646 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2647 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2648 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! only needed for kc=0 ...
2649 kc=kcmax-1
2650 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2651 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2652 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2653 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2654 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! ... and kc=KCMAX-1
2655 
2656 #elif (ADV_VERT==2 || ADV_VERT==3)
2657 
2658 do kc=0, kcmax-1
2659  ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2660  ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2661  ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2662  adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2663  abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
2664 end do
2665 
2666 #endif
2667 
2668 do kc=0, kcmax-1
2669  temp_c_help(kc) = 0.5_dp*(temp_c(kc,j,i)+temp_c(kc+1,j,i))
2670  ce6(kc) = at6(kc) &
2671  *kappa_val(temp_c_help(kc))/(c_val(temp_c_help(kc))*h_c(j,i))
2672  ci2(kc) = ai2(kc)/h_c(j,i)
2673 end do
2674 
2675 #if (ADV_HOR==3)
2676 dtt_dxi = 2.0_dp*dtt_2dxi
2677 dtt_deta = 2.0_dp*dtt_2deta
2678 #endif
2679 
2680 !-------- Set up the equations for the bedrock temperature --------
2681 
2682 kr=0
2683 lgs_a1(kr) = 1.0_dp
2684 lgs_a2(kr) = -1.0_dp
2685 lgs_b(kr) = clb1
2686 
2687 #if (Q_LITHO==1)
2688 ! (coupled heat-conducting bedrock)
2689 
2690 do kr=1, krmax-1
2691  lgs_a0(kr) = - ctr1
2692  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
2693  lgs_a2(kr) = - ctr1
2694  lgs_b(kr) = temp_r(kr,j,i)
2695 end do
2696 
2697 #elif (Q_LITHO==0)
2698 ! (no coupled heat-conducting bedrock)
2699 
2700 do kr=1, krmax-1
2701  lgs_a0(kr) = 1.0_dp
2702  lgs_a1(kr) = 0.0_dp
2703  lgs_a2(kr) = -1.0_dp
2704  lgs_b(kr) = 2.0_dp*clb1
2705 end do
2706 
2707 #endif
2708 
2709 kr=krmax
2710 lgs_a0(kr) = 0.0_dp
2711 lgs_a1(kr) = 1.0_dp
2712 lgs_b(kr) = temp_c_m(0,j,i)-delta_tm_sw
2713 
2714 !-------- Solve system of linear equations --------
2715 
2716 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
2717 
2718 !-------- Assign the result --------
2719 
2720 do kr=0, krmax
2721  temp_r_neu(kr,j,i) = lgs_x(kr)
2722 end do
2723 
2724 !-------- Set up the equations for the ice temperature --------
2725 
2726 kc=0
2727 lgs_a1(kc) = 1.0_dp
2728 lgs_a2(kc) = 0.0_dp
2729 lgs_b(kc) = enth_fct_temp_omega(temp_c_m(kc,j,i)-delta_tm_sw, 0.0_dp)
2730  ! zero water content assumed
2731 
2732 do kc=1, kcmax-1
2733 
2734 #if (ADV_VERT==1)
2735 
2736  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2737  -ce5(kc)*ce6(kc-1)
2738  lgs_a1(kc) = 1.0_dp+ce5(kc)*(ce6(kc)+ce6(kc-1))
2739  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2740  -ce5(kc)*ce6(kc)
2741 
2742 #elif (ADV_VERT==2)
2743 
2744  lgs_a0(kc) &
2745  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2746  -ce5(kc)*ce6(kc-1)
2747  lgs_a1(kc) &
2748  = 1.0_dp &
2749  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2750  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
2751  +ce5(kc)*(ce6(kc)+ce6(kc-1))
2752  lgs_a2(kc) &
2753  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
2754  -ce5(kc)*ce6(kc)
2755 
2756 #elif (ADV_VERT==3)
2757 
2758  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
2759 
2760  lgs_a0(kc) &
2761  = -max(adv_vert_help, 0.0_dp) &
2762  -ce5(kc)*ce6(kc-1)
2763  lgs_a1(kc) &
2764  = 1.0_dp &
2765  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
2766  +ce5(kc)*(ce6(kc)+ce6(kc-1))
2767  lgs_a2(kc) &
2768  = min(adv_vert_help, 0.0_dp) &
2769  -ce5(kc)*ce6(kc)
2770 
2771 #endif
2772 
2773 #if (ADV_HOR==2)
2774 
2775  lgs_b(kc) = enth_c(kc,j,i) + ce7(kc) &
2776  -dtt_2dxi* &
2777  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2778  *(enth_c(kc,j,i+1)-enth_c(kc,j,i)) &
2779  *insq_g11_sgx(j,i) &
2780  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2781  *(enth_c(kc,j,i)-enth_c(kc,j,i-1)) &
2782  *insq_g11_sgx(j,i-1) ) &
2783  -dtt_2deta* &
2784  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2785  *(enth_c(kc,j+1,i)-enth_c(kc,j,i)) &
2786  *insq_g22_sgy(j,i) &
2787  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2788  *(enth_c(kc,j,i)-enth_c(kc,j-1,i)) &
2789  *insq_g22_sgy(j-1,i) )
2790 
2791 #elif (ADV_HOR==3)
2792 
2793  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2794  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2795 
2796  lgs_b(kc) = enth_c(kc,j,i) + ce7(kc) &
2797  -dtt_dxi* &
2798  ( min(vx_c_help, 0.0_dp) &
2799  *(enth_c(kc,j,i+1)-enth_c(kc,j,i)) &
2800  *insq_g11_sgx(j,i) &
2801  +max(vx_c_help, 0.0_dp) &
2802  *(enth_c(kc,j,i)-enth_c(kc,j,i-1)) &
2803  *insq_g11_sgx(j,i-1) ) &
2804  -dtt_deta* &
2805  ( min(vy_c_help, 0.0_dp) &
2806  *(enth_c(kc,j+1,i)-enth_c(kc,j,i)) &
2807  *insq_g22_sgy(j,i) &
2808  +max(vy_c_help, 0.0_dp) &
2809  *(enth_c(kc,j,i)-enth_c(kc,j-1,i)) &
2810  *insq_g22_sgy(j-1,i) )
2811 
2812 #endif
2813 
2814 end do
2815 
2816 kc=kcmax
2817 lgs_a0(kc) = 0.0_dp
2818 lgs_a1(kc) = 1.0_dp
2819 lgs_b(kc) = enth_fct_temp_omega(temp_s(j,i), 0.0_dp)
2820  ! zero water content assumed
2821 
2822 !-------- Solve system of linear equations --------
2823 
2824 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
2825 
2826 !-------- Assign the result --------
2827 
2828 do kc=0, kcmax
2829  enth_c_neu(kc,j,i) = lgs_x(kc)
2830  temp_c_neu(kc,j,i) = temp_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
2831  omega_c_neu(kc,j,i) = omega_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
2832 end do
2833 
2834 !-------- Set enthalpy and water content in the redundant,
2835 ! lower (kt) ice layer to the value at the ice base --------
2836 
2837 do kt=0, ktmax
2838  enth_t_neu(kt,j,i) = enth_c_neu(0,j,i)
2839  omega_t_neu(kt,j,i) = omega_c_neu(0,j,i)
2840 end do
2841 
2842 !-------- Water drainage from the non-existing temperate ice --------
2843 
2844 q_tld(j,i) = 0.0_dp
2845 
2846 !-------- Set up the equations for the age of cold ice --------
2847 
2848 kc=0 ! adv_vert_sg(0) <= 0
2849 lgs_a1(kc) = 1.0_dp - min(adv_vert_sg(kc), 0.0_dp) ! (directed downward)
2850 lgs_a2(kc) = min(adv_vert_sg(kc), 0.0_dp) ! assumed/enforced
2851 
2852 #if (ADV_HOR==2)
2853 
2854 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2855  -dtt_2dxi* &
2856  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2857  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2858  *insq_g11_sgx(j,i) &
2859  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2860  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2861  *insq_g11_sgx(j,i-1) ) &
2862  -dtt_2deta* &
2863  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2864  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2865  *insq_g22_sgy(j,i) &
2866  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2867  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2868  *insq_g22_sgy(j-1,i) )
2869 
2870 #elif (ADV_HOR==3)
2871 
2872 vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2873 vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2874 
2875 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2876  -dtt_dxi* &
2877  ( min(vx_c_help, 0.0_dp) &
2878  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2879  *insq_g11_sgx(j,i) &
2880  +max(vx_c_help, 0.0_dp) &
2881  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2882  *insq_g11_sgx(j,i-1) ) &
2883  -dtt_deta* &
2884  ( min(vy_c_help, 0.0_dp) &
2885  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2886  *insq_g22_sgy(j,i) &
2887  +max(vy_c_help, 0.0_dp) &
2888  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2889  *insq_g22_sgy(j-1,i) )
2890 
2891 #endif
2892 
2893 do kc=1, kcmax-1
2894 
2895 #if (ADV_VERT==1)
2896 
2897  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2898  -ci1(kc)*ci2(kc-1)
2899  lgs_a1(kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
2900  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2901  -ci1(kc)*ci2(kc)
2902 
2903 #elif (ADV_VERT==2)
2904 
2905  lgs_a0(kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
2906  lgs_a1(kc) = 1.0_dp &
2907  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2908  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
2909  lgs_a2(kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
2910 
2911 #elif (ADV_VERT==3)
2912 
2913  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
2914 
2915  lgs_a0(kc) = -max(adv_vert_help, 0.0_dp)
2916  lgs_a1(kc) = 1.0_dp &
2917  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp)
2918  lgs_a2(kc) = min(adv_vert_help, 0.0_dp)
2919 
2920 #endif
2921 
2922 #if (ADV_HOR==2)
2923 
2924  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2925  -dtt_2dxi* &
2926  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2927  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2928  *insq_g11_sgx(j,i) &
2929  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2930  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2931  *insq_g11_sgx(j,i-1) ) &
2932  -dtt_2deta* &
2933  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2934  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2935  *insq_g22_sgy(j,i) &
2936  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2937  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2938  *insq_g22_sgy(j-1,i) )
2939 
2940 #elif (ADV_HOR==3)
2941 
2942  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2943  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2944 
2945  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2946  -dtt_dxi* &
2947  ( min(vx_c_help, 0.0_dp) &
2948  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2949  *insq_g11_sgx(j,i) &
2950  +max(vx_c_help, 0.0_dp) &
2951  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2952  *insq_g11_sgx(j,i-1) ) &
2953  -dtt_deta* &
2954  ( min(vy_c_help, 0.0_dp) &
2955  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2956  *insq_g22_sgy(j,i) &
2957  +max(vy_c_help, 0.0_dp) &
2958  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2959  *insq_g22_sgy(j-1,i) )
2960 
2961 #endif
2962 
2963 end do
2964 
2965 kc=kcmax
2966 if (as_perp(j,i) >= zero) then
2967  lgs_a0(kc) = 0.0_dp
2968  lgs_a1(kc) = 1.0_dp
2969  lgs_b(kc) = 0.0_dp
2970 else
2971  lgs_a0(kc) = -max(adv_vert_sg(kc-1), 0.0_dp)
2972  lgs_a1(kc) = 1.0_dp + max(adv_vert_sg(kc-1), 0.0_dp)
2973  ! adv_vert_sg(KCMAX-1) >= 0 (directed upward)
2974  ! assumed/enforced
2975 #if (ADV_HOR==2)
2976 
2977  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2978  -dtt_2dxi* &
2979  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2980  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2981  *insq_g11_sgx(j,i) &
2982  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2983  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2984  *insq_g11_sgx(j,i-1) ) &
2985  -dtt_2deta* &
2986  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2987  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2988  *insq_g22_sgy(j,i) &
2989  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2990  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2991  *insq_g22_sgy(j-1,i) )
2992 
2993 #elif (ADV_HOR==3)
2994 
2995  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2996  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2997 
2998  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2999  -dtt_dxi* &
3000  ( min(vx_c_help, 0.0_dp) &
3001  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
3002  *insq_g11_sgx(j,i) &
3003  +max(vx_c_help, 0.0_dp) &
3004  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
3005  *insq_g11_sgx(j,i-1) ) &
3006  -dtt_deta* &
3007  ( min(vy_c_help, 0.0_dp) &
3008  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
3009  *insq_g22_sgy(j,i) &
3010  +max(vy_c_help, 0.0_dp) &
3011  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
3012  *insq_g22_sgy(j-1,i) )
3013 
3014 #endif
3015 
3016 end if
3017 
3018 !-------- Solve system of linear equations --------
3019 
3020 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
3021 
3022 !-------- Assign the result,
3023 ! restriction to interval [0, AGE_MAX yr] --------
3024 
3025 do kc=0, kcmax
3026 
3027  age_c_neu(kc,j,i) = lgs_x(kc)
3028 
3029  if (age_c_neu(kc,j,i) < (age_min*year_sec)) &
3030  age_c_neu(kc,j,i) = 0.0_dp
3031  if (age_c_neu(kc,j,i) > (age_max*year_sec)) &
3032  age_c_neu(kc,j,i) = age_max*year_sec
3033 
3034 end do
3035 
3036 !-------- Age of the ice in the redundant, lower (kt) ice layer --------
3037 
3038 do kt=0, ktmax
3039  age_t_neu(kt,j,i) = age_c_neu(0,j,i)
3040 end do
3041 
3042 end subroutine calc_temp_enth_ssa
3043 
3044 !-------------------------------------------------------------------------------
3045 
3046 end module calc_temp_enth_m
3047 !
subroutine calc_temp_enth_1_a(at1, at2_1, at2_2, at3_1, at3_2, at4_1, at4_2, at5, at6, at7, atr1, acb1, acb2, acb3, acb4, alb1, ai1, ai2, dtime_temp, dtt_2dxi, dtt_2deta, i, j, ct1, ct2, ct3, ct4, ce5, ce6, ce7, ctr1, ccbe1, ccb2, ccb3, ccb4, clb1, ct1_sg, ct2_sg, ct3_sg, ct4_sg, adv_vert_sg, abs_adv_vert_sg, ci1, ci2, dtt_dxi, dtt_deta)
Computation of temperature and age for a cold ice column with the enthalpy method: Abbreviations...
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)) ...
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c
temp_c(kc,j,i): Temperature in the upper (kc) ice domain
real(dp) kappa_r
KAPPA_R: Heat conductivity of the lithosphere.
real(dp) delta_tm_sw
DELTA_TM_SW: Melting point depression of sea water due to its average salinity.
real(dp), dimension(0:krmax, 0:jmax, 0:imax) temp_r_neu
(.)_neu: New value of quantity (.) computed during an integration step
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
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) age_c
age_c(kc,j,i): Age in the upper (kc) ice domain
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...
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)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) enth_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp), dimension(0:kcmax) eaz_c
eaz_c(kc): Abbreviation for exp(aa*zeta(kc))
real(dp), dimension(0:jmax, 0:imax) q_geo
q_geo(j,i): Geothermal heat flux
real(dp), dimension(0:jmax, 0:imax) vx_b_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
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
real(dp), dimension(0:krmax, 0:jmax, 0:imax) temp_r
temp_r(kr,j,i): Temperature in the bedrock
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) age_t_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp) h_r
H_R: Thickness of the modelled lithosphere layer.
subroutine error(error_message)
Main routine of error_m: Writing of error messages and stopping execution.
Definition: error_m.F90:47
Material properties of ice: Rate factor, heat conductivity, specific heat (heat capacity), creep function, viscosity.
real(dp), parameter eps
eps: Small number
real(dp) omega_max
OMEGA_MAX: Threshold value for the water content of temperate ice.
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)
Conversion from temperature (temp) and water content (omega) to enthalpy (enth) and vice versa...
Several mathematical tools used by SICOPOLIS.
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)
real(dp), dimension(0:jmax, 0:imax) dh_c_dtau
dH_c_dtau(j,i): Derivative of H_c by tau (time)
real(dp) function, public kappa_val(temp_val)
Heat conductivity of ice: Linear interpolation of tabulated values in KAPPA(.).
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 ...
real(dp), dimension(0:jmax, 0:imax) dh_c_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp) function, public temp_fct_enth(enth_val, temp_m_val)
Temperature as a function of enthalpy.
real(dp), dimension(0:jmax, 0:imax) as_perp
as_perp(j,i): Accumulation-ablation function at the ice surface (SMB)
real(dp), dimension(0:jmax, 0:imax) dzs_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
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
real(dp), dimension(0:jmax, 0:imax) de_ssa
de_ssa(j,i): Effective strain rate of the SSA, at (i,j)
real(dp), dimension(0:jmax, 0:imax) h_t_neu
(.)_neu: New value of quantity (.) computed during an integration step
subroutine tri_sle(a0, a1, a2, x, b, nrows)
Solution of a system of linear equations Ax=b with tridiagonal matrix A.
real(dp) function, public creep(sigma_val)
Creep response function for ice.
real(dp), dimension(0:jmax, 0:imax) temp_s
temp_s(j,i): Ice surface temperature
real(dp) g
G: Acceleration due to gravity.
integer(i2b), dimension(0:jmax, 0:imax) kc_cts
kc_cts(j,i): Position kc of the CTS (for COLD, ENTC, ENTM)
real(dp), dimension(0:jmax, 0:imax) q_tld
Q_tld(j,i): Water drainage rate from the temperate layer.
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars_m.F90:35
real(dp) function, public omega_fct_enth(enth_val, temp_m_val)
Water content as a function of enthalpy.
real(dp) ea
ea: Abbreviation for exp(aa)
subroutine calc_temp_enth_1_d(ct1, ct2, ct3, ct4, ci1, ci2, ct1_sg, ct2_sg, ct3_sg, ct4_sg, adv_vert_sg, abs_adv_vert_sg, dtime_temp, dtt_dxi, dtt_deta, dtt_2dxi, dtt_2deta, i, j, lgs_a0, lgs_a1, lgs_a2, lgs_b)
Computation of temperature and age for a cold ice column with the enthalpy method: Set-up of the equa...
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) enth_c
enth_c(kc,j,i): Enthalpy in the upper (kc) ice domain
logical flag_aa_nonzero
flag_aa_nonzero: Flag for the exponential stretch parameter aa. .true.: aa greater than zero (non-equ...
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) age_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp) aa
aa: Exponential stretch parameter of the non-equidistant vertical grid in the upper (kc) ice domain ...
real(dp) rho_c_r
RHO_C_R: Density times specific heat of the lithosphere.
Declarations of kind types for SICOPOLIS.
subroutine calc_temp_enth_2_a1(at1, at2_1, at2_2, at3_1, at3_2, at4_1, at4_2, at5, atr1, alb1, ai1, aqtlde, dtime_temp, dtt_2dxi, dtt_2deta, i, j, ct1, ct2, ct3, ct4, ce5, ctr1, clb1, ct1_sg, ct2_sg, ct3_sg, ct4_sg, adv_vert_sg, abs_adv_vert_sg, ci1, cqtlde, dtt_dxi, dtt_deta)
Computation of temperature and age for an ice column with a temperate base with the enthalpy method: ...
real(dp), dimension(0:jmax, 0:imax) dh_c_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
subroutine calc_temp_enth_2_b(ctr1, clb1, i, j, lgs_a0, lgs_a1, lgs_a2, lgs_b)
Computation of temperature and age for an ice column with a temperate base with the enthalpy method: ...
integer(i2b), dimension(0:jmax, 0:imax) kc_cts_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
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...
real(dp), dimension(0:jmax, 0:imax) vy_b_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) omega_t_neu
(.)_neu: New value of quantity (.) computed during an integration step
Writing of error messages and stopping execution.
Definition: error_m.F90:35
real(dp), dimension(0:jmax, 0:imax) dzm_dtau
dzm_dtau(j,i): Derivative of zm by tau (time)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) omega_c
omega_c(kc,j,i): Water content in the upper (kc) ice domain
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
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) enth_t_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp) nue
NUE: Water diffusivity in ice.
real(dp) function, public c_val(temp_val)
Specific heat of ice: Linear interpolation of tabulated values in C(.).
real(dp), dimension(0:jmax, 0:imax) p_weert_inv
p_weert_inv(j,i): Inverse of p_weert
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
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)) ...
subroutine calc_temp_enth_1_b(ctr1, clb1, i, j, lgs_a0, lgs_a1, lgs_a2, lgs_b)
Computation of temperature and age for a cold ice column with the enthalpy method: Set-up of the equa...
real(dp), dimension(0:jmax, 0:imax) h_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp) function, public enth_fct_temp_omega(temp_val, omega_val)
Enthalpy as a function of temperature and water content.
Computation of temperature, water content and age with the enthalpy method.
subroutine calc_temp_enth_ssa(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 ice shelves (floating ice) with the enthalpy method...
real(dp), dimension(0:jmax, 0:imax) dzm_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) dzs_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
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)) ...
real(dp) beta
BETA: Clausius-Clapeyron gradient of ice.
subroutine calc_temp_enth_2_d(ct1, ct2, ct3, ct4, ci1, ci2, ct1_sg, ct2_sg, ct3_sg, ct4_sg, adv_vert_sg, abs_adv_vert_sg, dtime_temp, dtt_dxi, dtt_deta, dtt_2dxi, dtt_2deta, i, j, lgs_a0, lgs_a1, lgs_a2, lgs_b)
Computation of temperature and age for an ice column with a temperate base with the enthalpy method: ...
real(dp), dimension(0:jmax, 0:imax) dzm_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
subroutine, public calc_temp_enth(dxi, deta, dzeta_c, dzeta_t, dzeta_r, dtime_temp)
Main subroutine of calc_temp_enth_m: Computation of temperature, water content and age with the entha...
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)
real(dp), dimension(0:jmax, 0:imax) c_drag
c_drag(j,i): Auxiliary quantity for the computation of the basal drag
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)) ...
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...
real(dp), dimension(0:jmax, 0:imax) zm_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp) function, public ratefac_c_t(temp_val, omega_val, temp_m_val)
Rate factor for cold and temperate ice: Combination of ratefac_c and ratefac_t (only for the enthalpy...
character(len=256) errormsg
errormsg: Error-message string
real(dp) rho
RHO: Density of ice.
integer(i2b), dimension(0:jmax, 0:imax) n_cts_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) sigma_c
sigma_c(kc,j,i): Effective stress in the upper (kc) ice domain
real(dp), dimension(0:jmax, 0:imax) zb
zb(j,i): Coordinate z of the ice base
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) omega_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp), dimension(0:kcmax) zeta_c
zeta_c(kc): Sigma coordinate zeta_c of grid point kc
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...
Declarations of global variables for SICOPOLIS.
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)) ...