SICOPOLIS V5-dev  Revision 1368
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_i1b) then ! glaciated land
173 
174 ! ------ Old vertical column cold
175 
176  if (n_cts(j,i) == -1_i1b) then
177 
178  n_cts_neu(j,i) = -1_i1b
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_i1b
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_i1b) then
207 
208  n_cts_neu(j,i) = 0_i1b
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_i1b
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_i1b
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_i1b) then ! floating ice
251 
252  n_cts_neu(j,i) = -1_i1b
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_i1b, 2_i1b (ice-free land or sea point)
277 
278  n_cts_neu(j,i) = -1_i1b
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_i1b).or.(maske(j,i) == 3_i1b) ) 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_i1b, 2_i1b (ice-free land or sea point)
328 
329  n_cts_neu(j,i) = -1_i1b
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_i1b).or.(maske(j,i) == 3_i1b) ) 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_i1b, 2_i1b (ice-free land or sea point)
374 
375  n_cts_neu(j,i) = -1_i1b
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_i1b).or.(maske(j,i) == 3_i1b) ) 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_i1b, 2_i1b (ice-free land or sea point)
420 
421  n_cts_neu(j,i) = -1_i1b
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_i1b).or.(maske(j,i) == 3_i1b) ) 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_i1b, 2_i1b (ice-free land or sea point)
466 
467  n_cts_neu(j,i) = -1_i1b
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_i1b).or.(maske(j,i) == 3_i1b) ) 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_i1b, 2_i1b (ice-free land or sea point)
515 
516  n_cts_neu(j,i) = -1_i1b
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_i1b).or.(maske(j,i) == 3_i1b) ) 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_i1b, 2_i1b (ice-free land or sea point)
560 
561  n_cts_neu(j,i) = -1_i1b
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_i1b).or.(maske(j,i) == 3_i1b) ) 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_i1b, 2_i1b (ice-free land or sea point)
611 
612  n_cts_neu(j,i) = -1_i1b
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_i1b).or.(maske(j,i) == 3_i1b) ) 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_i1b, 2_i1b (ice-free land or sea point)
656 
657  n_cts_neu(j,i) = -1_i1b
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 #ifndef ALLOW_OPENAD /* Normal */
924  ccb3 = -c_drag(j,i) &
925  * sqrt(vx_b_g(j,i)**2 &
926  +vy_b_g(j,i)**2) &
927  **(1.0_dp+p_weert_inv(j,i))
928 #else /* OpenAD: guarding against undifferentiable sqrt(0) */
929  if ((vx_b_g(j,i)**2+vy_b_g(j,i)**2) == 0) then
930  ccb3 = 0.0_dp
931  else
932  ccb3 = -c_drag(j,i) &
933  * sqrt(vx_b_g(j,i)**2 &
934  +vy_b_g(j,i)**2) &
935  **(1.0_dp+p_weert_inv(j,i))
936  end if
937 #endif /* Normal vs. OpenAD */
938  ccb4 = 0.0_dp
939 
940 end if
941 #endif
942 
943 clb1 = alb1*q_geo(j,i)
944 
945 #if (ADV_VERT==1)
946 
947 do kc=1, kcmax-1
948  ct1(kc) = at1(kc)/h_c(j,i)*0.5_dp*(vz_c(kc,j,i)+vz_c(kc-1,j,i))
949 end do
950 
951 kc=0
952 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
953  ! only needed for kc=0 ...
954 kc=kcmax-1
955 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
956  ! ... and kc=KCMAX-1
957 
958 #elif (ADV_VERT==2 || ADV_VERT==3)
959 
960 do kc=0, kcmax-1
961  ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
962 end do
963 
964 #endif
965 
966 do kc=0, kcmax
967 
968  ct2(kc) = ( at2_1(kc)*dzm_dtau(j,i) &
969  +at2_2(kc)*dh_c_dtau(j,i) )/h_c(j,i)
970  ct3(kc) = ( at3_1(kc)*dzm_dxi_g(j,i) &
971  +at3_2(kc)*dh_c_dxi_g(j,i) )/h_c(j,i) &
972  *0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1)) *insq_g11_g(j,i)
973  ct4(kc) = ( at4_1(kc)*dzm_deta_g(j,i) &
974  +at4_2(kc)*dh_c_deta_g(j,i) )/h_c(j,i) &
975  *0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i)) *insq_g22_g(j,i)
976  ce5(kc) = at5(kc)/h_c(j,i)
977 
978 #if (DYNAMICS==2)
979  if (.not.flag_shelfy_stream(j,i)) then
980 #endif
981  ce7(kc) = at7 &
982  *enh_c(kc,j,i) &
983  *ratefac_c_t(temp_c(kc,j,i), omega_c(kc,j,i), temp_c_m(kc,j,i)) &
984  *creep(sigma_c(kc,j,i)) &
985  *sigma_c(kc,j,i)*sigma_c(kc,j,i)
986 #if (DYNAMICS==2)
987  else
988  ce7(kc) = 2.0_dp*at7 &
989  *viscosity(de_c(kc,j,i), &
990  temp_c(kc,j,i), temp_c_m(kc,j,i), 0.0_dp, &
991  enh_c(kc,j,i), 0_i1b) &
992  *de_c(kc,j,i)**2
993  end if
994 #endif
995 
996  ci1(kc) = ai1(kc)/h_c(j,i)
997 
998 end do
999 
1000 #if (ADV_VERT==1)
1001 
1002 kc=0
1003 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
1004 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
1005 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
1006 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
1007 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! only needed for kc=0 ...
1008 kc=kcmax-1
1009 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
1010 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
1011 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
1012 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
1013 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! ... and kc=KCMAX-1
1014 
1015 #elif (ADV_VERT==2 || ADV_VERT==3)
1016 
1017 do kc=0, kcmax-1
1018  ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
1019  ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
1020  ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
1021  adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
1022  abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
1023 end do
1024 
1025 #endif
1026 
1027 do kc=0, kcmax-1
1028  temp_c_help(kc) = 0.5_dp*(temp_c(kc,j,i)+temp_c(kc+1,j,i))
1029  ce6(kc) = at6(kc) &
1030  *kappa_val(temp_c_help(kc))/(c_val(temp_c_help(kc))*h_c(j,i))
1031  ci2(kc) = ai2(kc)/h_c(j,i)
1032 end do
1033 
1034 #if (ADV_HOR==3)
1035 dtt_dxi = 2.0_dp*dtt_2dxi
1036 dtt_deta = 2.0_dp*dtt_2deta
1037 #endif
1038 
1039 end subroutine calc_temp_enth_1_a
1040 
1041 !-------------------------------------------------------------------------------
1042 !> Computation of temperature and age for a cold ice column with the
1043 !! enthalpy method:
1044 !! Set-up of the equations for the bedrock temperature.
1045 !<------------------------------------------------------------------------------
1046 subroutine calc_temp_enth_1_b(ctr1, clb1, i, j, &
1047  lgs_a0, lgs_a1, lgs_a2, lgs_b)
1049 implicit none
1050 
1051 #ifndef ALLOW_OPENAD /* Normal */
1052 integer(i4b), intent(in) :: i, j
1053 #else /* OpenAD */
1054 integer(i4b), intent(inout) :: i, j
1055 #endif /* Normal vs. OpenAD */
1056 
1057 real(dp), intent(in) :: ctr1, clb1
1058 
1059 real(dp), intent(out) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
1060  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
1061  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
1062  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
1063 
1064 integer(i4b) :: kr
1065 
1066 !-------- Initialisation --------
1067 
1068 lgs_a0 = 0.0_dp
1069 lgs_a1 = 0.0_dp
1070 lgs_a2 = 0.0_dp
1071 lgs_b = 0.0_dp
1072 
1073 !-------- Actual computation --------
1074 
1075 kr=0
1076 lgs_a1(kr) = 1.0_dp
1077 lgs_a2(kr) = -1.0_dp
1078 lgs_b(kr) = clb1
1079 
1080 #if (Q_LITHO==1)
1081 ! (coupled heat-conducting bedrock)
1082 
1083 do kr=1, krmax-1
1084  lgs_a0(kr) = -ctr1
1085  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
1086  lgs_a2(kr) = -ctr1
1087  lgs_b(kr) = temp_r(kr,j,i)
1088 end do
1089 
1090 #elif (Q_LITHO==0)
1091 ! (no coupled heat-conducting bedrock)
1092 
1093 do kr=1, krmax-1
1094  lgs_a0(kr) = 1.0_dp
1095  lgs_a1(kr) = 0.0_dp
1096  lgs_a2(kr) = -1.0_dp
1097  lgs_b(kr) = 2.0_dp*clb1
1098 end do
1099 
1100 #endif
1101 
1102 kr=krmax
1103 lgs_a0(kr) = 0.0_dp
1104 lgs_a1(kr) = 1.0_dp
1105 lgs_b(kr) = temp_c(0,j,i)
1106 
1107 end subroutine calc_temp_enth_1_b
1108 
1109 !-------------------------------------------------------------------------------
1110 !> Computation of temperature and age for a cold ice column with the
1111 !! enthalpy method:
1112 !! Set-up of the equations for the ice enthalpy.
1113 !<------------------------------------------------------------------------------
1114 subroutine calc_temp_enth_1_c(ct1, ct2, ct3, ct4, ce5, ce6, ce7, &
1115  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
1116  ccbe1, ccb2, ccb3, ccb4, &
1117  adv_vert_sg, abs_adv_vert_sg, &
1118  dtime_temp, dtt_dxi, dtt_deta, &
1119  dtt_2dxi, dtt_2deta, &
1120  i, j, &
1121  lgs_a0, lgs_a1, lgs_a2, lgs_b)
1122 
1124 
1125 implicit none
1126 
1127 #ifndef ALLOW_OPENAD /* Normal */
1128 integer(i4b), intent(in) :: i, j
1129 #else /* OpenAD */
1130 integer(i4b), intent(inout) :: i, j
1131 #endif /* Normal vs. OpenAD */
1132 
1133 real(dp), intent(in) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), &
1134  ct4(0:kcmax), ce5(0:kcmax), ce6(0:kcmax), &
1135  ce7(0:kcmax)
1136 real(dp), intent(in) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), &
1137  ct3_sg(0:kcmax), ct4_sg(0:kcmax), &
1138  ccbe1, ccb2, ccb3, ccb4, &
1139  adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
1140 real(dp), intent(in) :: dtime_temp, dtt_dxi, dtt_deta, dtt_2dxi, dtt_2deta
1141 
1142 real(dp), intent(out) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
1143  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
1144  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
1145  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
1146 
1147 integer(i4b) :: kc, kr
1148 real(dp) :: vx_c_help, vy_c_help
1149 real(dp) :: adv_vert_help
1150 
1151 !-------- Initialisation --------
1152 
1153 lgs_a0 = 0.0_dp
1154 lgs_a1 = 0.0_dp
1155 lgs_a2 = 0.0_dp
1156 lgs_b = 0.0_dp
1157 
1158 !-------- Actual computation --------
1159 
1160 kr=krmax
1161 kc=0
1162 lgs_a1(kc) = -ccbe1
1163 lgs_a2(kc) = ccbe1
1164 lgs_b(kc) = ccb2*(temp_r_neu(kr,j,i)-temp_r_neu(kr-1,j,i)) + ccb3 + ccb4
1165 
1166 do kc=1, kcmax-1
1167 
1168 #if (ADV_VERT==1)
1169 
1170  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
1171  -ce5(kc)*ce6(kc-1)
1172  lgs_a1(kc) = 1.0_dp+ce5(kc)*(ce6(kc)+ce6(kc-1))
1173  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
1174  -ce5(kc)*ce6(kc)
1175 
1176 #elif (ADV_VERT==2)
1177 
1178  lgs_a0(kc) &
1179  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
1180  -ce5(kc)*ce6(kc-1)
1181  lgs_a1(kc) &
1182  = 1.0_dp &
1183  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
1184  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
1185  +ce5(kc)*(ce6(kc)+ce6(kc-1))
1186  lgs_a2(kc) &
1187  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
1188  -ce5(kc)*ce6(kc)
1189 
1190 #elif (ADV_VERT==3)
1191 
1192  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
1193 
1194  lgs_a0(kc) &
1195  = -max(adv_vert_help, 0.0_dp) &
1196  -ce5(kc)*ce6(kc-1)
1197  lgs_a1(kc) &
1198  = 1.0_dp &
1199  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
1200  +ce5(kc)*(ce6(kc)+ce6(kc-1))
1201  lgs_a2(kc) &
1202  = min(adv_vert_help, 0.0_dp) &
1203  -ce5(kc)*ce6(kc)
1204 
1205 #endif
1206 
1207 #if (ADV_HOR==2)
1208 
1209  lgs_b(kc) = enth_c(kc,j,i) + ce7(kc) &
1210  -dtt_2dxi* &
1211  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
1212  *(enth_c(kc,j,i+1)-enth_c(kc,j,i)) &
1213  *insq_g11_sgx(j,i) &
1214  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
1215  *(enth_c(kc,j,i)-enth_c(kc,j,i-1)) &
1216  *insq_g11_sgx(j,i-1) ) &
1217  -dtt_2deta* &
1218  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
1219  *(enth_c(kc,j+1,i)-enth_c(kc,j,i)) &
1220  *insq_g22_sgy(j,i) &
1221  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
1222  *(enth_c(kc,j,i)-enth_c(kc,j-1,i)) &
1223  *insq_g22_sgy(j-1,i) )
1224 
1225 #elif (ADV_HOR==3)
1226 
1227  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
1228  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
1229 
1230  lgs_b(kc) = enth_c(kc,j,i) + ce7(kc) &
1231  -dtt_dxi* &
1232  ( min(vx_c_help, 0.0_dp) &
1233  *(enth_c(kc,j,i+1)-enth_c(kc,j,i)) &
1234  *insq_g11_sgx(j,i) &
1235  +max(vx_c_help, 0.0_dp) &
1236  *(enth_c(kc,j,i)-enth_c(kc,j,i-1)) &
1237  *insq_g11_sgx(j,i-1) ) &
1238  -dtt_deta* &
1239  ( min(vy_c_help, 0.0_dp) &
1240  *(enth_c(kc,j+1,i)-enth_c(kc,j,i)) &
1241  *insq_g22_sgy(j,i) &
1242  +max(vy_c_help, 0.0_dp) &
1243  *(enth_c(kc,j,i)-enth_c(kc,j-1,i)) &
1244  *insq_g22_sgy(j-1,i) )
1245 
1246 #endif
1247 
1248 end do
1249 
1250 kc=kcmax
1251 lgs_a0(kc) = 0.0_dp
1252 lgs_a1(kc) = 1.0_dp
1253 lgs_b(kc) = enth_fct_temp_omega(temp_s(j,i), 0.0_dp)
1254  ! zero water content at the ice surface
1255 
1256 end subroutine calc_temp_enth_1_c
1257 
1258 !-------------------------------------------------------------------------------
1259 !> Computation of temperature and age for a cold ice column with the
1260 !! enthalpy method:
1261 !! Set-up of the equations for the age of ice.
1262 !<------------------------------------------------------------------------------
1263 subroutine calc_temp_enth_1_d(ct1, ct2, ct3, ct4, ci1, ci2, &
1264  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
1265  adv_vert_sg, abs_adv_vert_sg, &
1266  dtime_temp, dtt_dxi, dtt_deta, &
1267  dtt_2dxi, dtt_2deta, &
1268  i, j, &
1269  lgs_a0, lgs_a1, lgs_a2, lgs_b)
1271 implicit none
1272 
1273 #ifndef ALLOW_OPENAD /* Normal */
1274 integer(i4b), intent(in) :: i, j
1275 #else /* OpenAD */
1276 integer(i4b), intent(inout) :: i, j
1277 #endif /* Normal vs. OpenAD */
1278 
1279 real(dp), intent(in) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), &
1280  ct4(0:kcmax), ci1(0:kcmax), ci2(0:kcmax)
1281 real(dp), intent(in) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), &
1282  ct3_sg(0:kcmax), ct4_sg(0:kcmax), &
1283  adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
1284 real(dp), intent(in) :: dtime_temp, dtt_dxi, dtt_deta, dtt_2dxi, dtt_2deta
1285 
1286 real(dp), intent(out) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
1287  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
1288  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
1289  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
1290 
1291 integer(i4b) :: kc
1292 real(dp) :: vx_c_help, vy_c_help
1293 real(dp) :: adv_vert_help
1294 
1295 !-------- Initialisation --------
1296 
1297 lgs_a0 = 0.0_dp
1298 lgs_a1 = 0.0_dp
1299 lgs_a2 = 0.0_dp
1300 lgs_b = 0.0_dp
1301 
1302 !-------- Actual computation --------
1303 
1304 kc=0 ! adv_vert_sg(0) <= 0
1305 lgs_a1(kc) = 1.0_dp - min(adv_vert_sg(kc), 0.0_dp) ! (directed downward)
1306 lgs_a2(kc) = min(adv_vert_sg(kc), 0.0_dp) ! assumed/enforced
1307 
1308 #if (ADV_HOR==2)
1309 
1310 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1311  -dtt_2dxi* &
1312  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
1313  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1314  *insq_g11_sgx(j,i) &
1315  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
1316  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1317  *insq_g11_sgx(j,i-1) ) &
1318  -dtt_2deta* &
1319  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
1320  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1321  *insq_g22_sgy(j,i) &
1322  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
1323  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1324  *insq_g22_sgy(j-1,i) )
1325 
1326 #elif (ADV_HOR==3)
1327 
1328 vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
1329 vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
1330 
1331 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1332  -dtt_dxi* &
1333  ( min(vx_c_help, 0.0_dp) &
1334  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1335  *insq_g11_sgx(j,i) &
1336  +max(vx_c_help, 0.0_dp) &
1337  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1338  *insq_g11_sgx(j,i-1) ) &
1339  -dtt_deta* &
1340  ( min(vy_c_help, 0.0_dp) &
1341  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1342  *insq_g22_sgy(j,i) &
1343  +max(vy_c_help, 0.0_dp) &
1344  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1345  *insq_g22_sgy(j-1,i) )
1346 
1347 #endif
1348 
1349 do kc=1, kcmax-1
1350 
1351 #if (ADV_VERT==1)
1352 
1353  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
1354  -ci1(kc)*ci2(kc-1)
1355  lgs_a1(kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
1356  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
1357  -ci1(kc)*ci2(kc)
1358 
1359 #elif (ADV_VERT==2)
1360 
1361  lgs_a0(kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
1362  lgs_a1(kc) = 1.0_dp &
1363  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
1364  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
1365  lgs_a2(kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
1366 
1367 #elif (ADV_VERT==3)
1368 
1369  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
1370 
1371  lgs_a0(kc) = -max(adv_vert_help, 0.0_dp)
1372  lgs_a1(kc) = 1.0_dp &
1373  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp)
1374  lgs_a2(kc) = min(adv_vert_help, 0.0_dp)
1375 
1376 #endif
1377 
1378 #if (ADV_HOR==2)
1379 
1380  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1381  -dtt_2dxi* &
1382  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
1383  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1384  *insq_g11_sgx(j,i) &
1385  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
1386  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1387  *insq_g11_sgx(j,i-1) ) &
1388  -dtt_2deta* &
1389  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
1390  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1391  *insq_g22_sgy(j,i) &
1392  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
1393  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1394  *insq_g22_sgy(j-1,i) )
1395 
1396 #elif (ADV_HOR==3)
1397 
1398  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
1399  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
1400 
1401  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1402  -dtt_dxi* &
1403  ( min(vx_c_help, 0.0_dp) &
1404  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1405  *insq_g11_sgx(j,i) &
1406  +max(vx_c_help, 0.0_dp) &
1407  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1408  *insq_g11_sgx(j,i-1) ) &
1409  -dtt_deta* &
1410  ( min(vy_c_help, 0.0_dp) &
1411  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1412  *insq_g22_sgy(j,i) &
1413  +max(vy_c_help, 0.0_dp) &
1414  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1415  *insq_g22_sgy(j-1,i) )
1416 
1417 #endif
1418 
1419 end do
1420 
1421 kc=kcmax
1422 if (as_perp(j,i) >= 0.0_dp) then
1423  lgs_a0(kc) = 0.0_dp
1424  lgs_a1(kc) = 1.0_dp
1425  lgs_b(kc) = 0.0_dp
1426 else
1427  lgs_a0(kc) = -max(adv_vert_sg(kc-1), 0.0_dp)
1428  lgs_a1(kc) = 1.0_dp + max(adv_vert_sg(kc-1), 0.0_dp)
1429  ! adv_vert_sg(KCMAX-1) >= 0 (directed upward)
1430  ! assumed/enforced
1431 #if (ADV_HOR==2)
1432 
1433  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1434  -dtt_2dxi* &
1435  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
1436  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1437  *insq_g11_sgx(j,i) &
1438  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
1439  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1440  *insq_g11_sgx(j,i-1) ) &
1441  -dtt_2deta* &
1442  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
1443  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1444  *insq_g22_sgy(j,i) &
1445  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
1446  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1447  *insq_g22_sgy(j-1,i) )
1448 
1449 #elif (ADV_HOR==3)
1450 
1451  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
1452  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
1453 
1454  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
1455  -dtt_dxi* &
1456  ( min(vx_c_help, 0.0_dp) &
1457  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
1458  *insq_g11_sgx(j,i) &
1459  +max(vx_c_help, 0.0_dp) &
1460  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
1461  *insq_g11_sgx(j,i-1) ) &
1462  -dtt_deta* &
1463  ( min(vy_c_help, 0.0_dp) &
1464  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
1465  *insq_g22_sgy(j,i) &
1466  +max(vy_c_help, 0.0_dp) &
1467  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
1468  *insq_g22_sgy(j-1,i) )
1469 
1470 #endif
1471 
1472 end if
1473 
1474 end subroutine calc_temp_enth_1_d
1475 
1476 !-------------------------------------------------------------------------------
1477 !> Computation of temperature and age for an ice column with a temperate base
1478 !! with the enthalpy method.
1479 !<------------------------------------------------------------------------------
1480 subroutine calc_temp_enth_2(at1, at2_1, at2_2, at3_1, at3_2, &
1481  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
1482  ai1, ai2, aqtlde, am3, &
1483  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
1484 
1485 #ifndef ALLOW_OPENAD /* Normal */
1486 use sico_maths_m, only : tri_sle
1487 #else /* OpenAD */
1488 use sico_maths_m
1489 #endif /* Normal vs. OpenAD */
1490 
1493 
1494 implicit none
1495 
1496 #ifndef ALLOW_OPENAD /* Normal */
1497 integer(i4b), intent(in) :: i, j
1498 #else /* OpenAD */
1499 integer(i4b), intent(inout) :: i, j
1500 #endif /* Normal vs. OpenAD */
1501 
1502 real(dp), intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
1503  at3_1(0:kcmax), at3_2(0:kcmax), &
1504  at4_1(0:kcmax), at4_2(0:kcmax), &
1505  at5(0:kcmax), at6(0:kcmax), at7, &
1506  ai1(0:kcmax), ai2(0:kcmax), &
1507  atr1, alb1, aqtlde(0:kcmax), am3(0:kcmax)
1508 real(dp), intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
1509 
1510 integer(i4b) :: kc, kt, kr
1511 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
1512  ce5(0:kcmax), ce6(0:kcmax), ce7(0:kcmax), ctr1, clb1
1513 real(dp) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), ct3_sg(0:kcmax), &
1514  ct4_sg(0:kcmax), adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
1515 real(dp) :: ci1(0:kcmax), ci2(0:kcmax)
1516 real(dp) :: cqtlde(0:kcmax), cm3(0:kcmax)
1517 real(dp) :: dtt_dxi, dtt_deta
1518 real(dp) :: temp_c_val(0:kcmax), omega_c_val(0:kcmax)
1519 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
1520  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
1521  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
1522  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
1523  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
1524 
1525 real(dp), parameter :: eps_omega=1.0e-12_dp
1526 
1527 #ifdef ALLOW_OPENAD /* OpenAD */
1528 logical :: kcdone
1529 #endif /* OpenAD */
1530 
1531 !-------- Check for boundary points --------
1532 
1533 if ((i == 0).or.(i == imax).or.(j == 0).or.(j == jmax)) then
1534  errormsg = ' >>> calc_temp_enth_2: Boundary points not allowed.'
1535  call error(errormsg)
1536 end if
1537 
1538 !-------- Abbreviations --------
1539 
1540 call calc_temp_enth_2_a1(at1, at2_1, at2_2, at3_1, at3_2, &
1541  at4_1, at4_2, at5, atr1, alb1, &
1542  ai1, aqtlde, &
1543  dtime_temp, dtt_2dxi, dtt_2deta, i, j, &
1544  ct1, ct2, ct3, ct4, ce5, ctr1, clb1, &
1545  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
1546  adv_vert_sg, abs_adv_vert_sg, &
1547  ci1, cqtlde, dtt_dxi, dtt_deta)
1548 
1549 do kc=0, kcmax
1550  temp_c_val(kc) = temp_c(kc,j,i)
1551  omega_c_val(kc) = omega_c(kc,j,i)
1552 end do
1553 
1554 call calc_temp_enth_2_a2(at6, at7, ai2, am3, temp_c_val, omega_c_val, &
1555  i, j, ce6, ce7, ci2, cm3)
1556 
1557 !-------- Computation of the bedrock temperature --------
1558 
1559 ! ------ Set-up of the the equations
1560 
1561 call calc_temp_enth_2_b(ctr1, clb1, i, j, &
1562  lgs_a0, lgs_a1, lgs_a2, lgs_b)
1563 
1564 ! ------ Solution of the system of linear equations
1565 
1566 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
1567 
1568 ! ------ Assignment of the result
1569 
1570 do kr=0, krmax
1571  temp_r_neu(kr,j,i) = lgs_x(kr)
1572 end do
1573 
1574 !-------- Computation of the ice enthalpy (predictor step) --------
1575 
1576 ! ------ Set-up of the equations
1577 
1578 call calc_temp_enth_2_c(ct1, ct2, ct3, ct4, ce5, ce6, ce7, &
1579  ct1_sg, ct2_sg, ct3_sg, ct4_sg, cm3, &
1580  adv_vert_sg, abs_adv_vert_sg, &
1581  dtime_temp, dtt_dxi, dtt_deta, &
1582  dtt_2dxi, dtt_2deta, &
1583  i, j, 0, &
1584  lgs_a0, lgs_a1, lgs_a2, lgs_b)
1585 
1586 ! ------ Solution of the system of linear equations
1587 
1588 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
1589 
1590 ! ------ Assignment of the result
1591 
1592 do kc=0, kcmax
1593  enth_c_neu(kc,j,i) = lgs_x(kc)
1594  temp_c_neu(kc,j,i) = temp_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
1595  omega_c_neu(kc,j,i) = omega_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
1596 end do
1597 
1598 ! ------ Determination of the CTS
1599 
1600 kc_cts_neu(j,i) = 0
1601 
1602 #ifndef ALLOW_OPENAD /* Normal */
1603 
1604 do kc=1, kcmax-1
1605  if (omega_c_neu(kc,j,i) > eps_omega) then
1606  kc_cts_neu(j,i) = kc
1607  else
1608  exit
1609  end if
1610 end do
1611 
1612 #else /* OpenAD */
1613 
1614 kcdone = .false.
1615 
1616 do kc=1, kcmax-1
1617  if (kcdone.eqv..false.) then
1618  if (omega_c_neu(kc,j,i) > eps_omega) then
1619  kc_cts_neu(j,i) = kc
1620  else
1621  kcdone = .true.
1622  end if
1623  end if
1624 end do
1625 
1626 #endif /* Normal vs. OpenAD */
1627 
1628 !-------- Computation of the ice enthalpy
1629 ! (corrector step for the cold-ice domain only
1630 ! in order to fulfull the transition condition at the CTS) --------
1631 
1632 #if (CALCMOD==3) /* ENTM scheme */
1633 
1634 if (kc_cts_neu(j,i) > 0) then
1635 
1636 ! ------ Update of the abbreviations where needed
1637 
1638  do kc=0, kcmax
1639  temp_c_val(kc) = temp_c_neu(kc,j,i)
1640  omega_c_val(kc) = omega_c_neu(kc,j,i)
1641  end do
1642 
1643  call calc_temp_enth_2_a2(at6, at7, ai2, am3, temp_c_val, omega_c_val, &
1644  i, j, ce6, ce7, ci2, cm3)
1645 
1646 ! ------ Set-up of the equations
1647 
1648  call calc_temp_enth_2_c(ct1, ct2, ct3, ct4, ce5, ce6, ce7, &
1649  ct1_sg, ct2_sg, ct3_sg, ct4_sg, cm3, &
1650  adv_vert_sg, abs_adv_vert_sg, &
1651  dtime_temp, dtt_dxi, dtt_deta, &
1652  dtt_2dxi, dtt_2deta, &
1653  i, j, kc_cts_neu(j,i), &
1654  lgs_a0, lgs_a1, lgs_a2, lgs_b)
1655 
1656 ! ------ Solution of the system of linear equations
1657 
1658  call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
1659 
1660 ! ------ Assignment of the result
1661 
1662  kc=kc_cts_neu(j,i)
1663  enth_c_neu(kc,j,i) = lgs_x(kc)
1664  temp_c_neu(kc,j,i) = temp_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
1665  omega_c_neu(kc,j,i) = omega_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
1666 
1667  do kc=kc_cts_neu(j,i)+1, kcmax
1668  enth_c_neu(kc,j,i) = lgs_x(kc)
1669  temp_c_neu(kc,j,i) = temp_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
1670  omega_c_neu(kc,j,i) = 0.0_dp ! cold-ice domain
1671  end do
1672 
1673 end if
1674 
1675 #elif (CALCMOD==2) /* ENTC scheme */
1676 
1677 !!! continue ! no corrector step
1678 
1679 #else
1680 
1681 errormsg = ' >>> calc_temp_enth_2: CALCMOD must be either 2 or 3!'
1682 call error(errormsg)
1683 
1684 #endif
1685 
1686 !-------- Water drainage from temperate ice (if existing) --------
1687 
1688 q_tld(j,i) = 0.0_dp
1689 
1690 do kc=0, kc_cts_neu(j,i)
1691 
1692  if (omega_c_neu(kc,j,i) > omega_max) then
1693 
1694  q_tld(j,i) = q_tld(j,i) + cqtlde(kc)*(omega_c_neu(kc,j,i)-omega_max)
1695 
1696  omega_c_neu(kc,j,i) = omega_max
1698 
1699  end if
1700 
1701 end do
1702 
1703 !-------- Set enthalpy and water content in the redundant,
1704 ! lower (kt) ice layer to the value at the ice base --------
1705 
1706 do kt=0, ktmax
1707  enth_t_neu(kt,j,i) = enth_c_neu(0,j,i)
1708  omega_t_neu(kt,j,i) = omega_c_neu(0,j,i)
1709 end do
1710 
1711 !-------- Computation of the age of ice --------
1712 
1713 ! ------ Set-up of the the equations
1714 
1715 call calc_temp_enth_2_d(ct1, ct2, ct3, ct4, ci1, ci2, &
1716  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
1717  adv_vert_sg, abs_adv_vert_sg, &
1718  dtime_temp, dtt_dxi, dtt_deta, &
1719  dtt_2dxi, dtt_2deta, &
1720  i, j, &
1721  lgs_a0, lgs_a1, lgs_a2, lgs_b)
1722 
1723 ! ------ Solution of the system of linear equations
1724 
1725 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
1726 
1727 ! ------ Assignment of the result
1728 ! restriction to interval [0, AGE_MAX yr]
1729 
1730 do kc=0, kcmax
1731 
1732  age_c_neu(kc,j,i) = lgs_x(kc)
1733 
1734  if (age_c_neu(kc,j,i) < (age_min*year_sec)) &
1735  age_c_neu(kc,j,i) = 0.0_dp
1736  if (age_c_neu(kc,j,i) > (age_max*year_sec)) &
1737  age_c_neu(kc,j,i) = age_max*year_sec
1738 
1739 end do
1740 
1741 !-------- Age of the ice in the redundant, lower (kt) ice layer --------
1742 
1743 do kt=0, ktmax
1744  age_t_neu(kt,j,i) = age_c_neu(0,j,i)
1745 end do
1746 
1747 end subroutine calc_temp_enth_2
1748 
1749 !-------------------------------------------------------------------------------
1750 !> Computation of temperature and age for an ice column with a temperate base
1751 !! with the enthalpy method: Abbreviations I.
1752 !<------------------------------------------------------------------------------
1753 subroutine calc_temp_enth_2_a1(at1, at2_1, at2_2, at3_1, at3_2, &
1754  at4_1, at4_2, at5, atr1, alb1, &
1755  ai1, aqtlde, &
1756  dtime_temp, dtt_2dxi, dtt_2deta, i, j, &
1757  ct1, ct2, ct3, ct4, ce5, ctr1, clb1, &
1758  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
1759  adv_vert_sg, abs_adv_vert_sg, &
1760  ci1, cqtlde, dtt_dxi, dtt_deta)
1762 implicit none
1763 
1764 #ifndef ALLOW_OPENAD /* Normal */
1765 integer(i4b), intent(in) :: i, j
1766 #else /* OpenAD */
1767 integer(i4b), intent(inout) :: i, j
1768 #endif /* Normal vs. OpenAD */
1769 
1770 real(dp), intent(in) :: at1(0:kcmax), &
1771  at2_1(0:kcmax), at2_2(0:kcmax), &
1772  at3_1(0:kcmax), at3_2(0:kcmax), &
1773  at4_1(0:kcmax), at4_2(0:kcmax), &
1774  at5(0:kcmax), ai1(0:kcmax), &
1775  atr1, alb1, aqtlde(0:kcmax)
1776 real(dp), intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
1777 
1778 real(dp), intent(out) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), &
1779  ct4(0:kcmax), ce5(0:kcmax), &
1780  ctr1, clb1
1781 real(dp), intent(out) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), &
1782  ct3_sg(0:kcmax), ct4_sg(0:kcmax), &
1783  adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
1784 real(dp), intent(out) :: ci1(0:kcmax), cqtlde(0:kcmax)
1785 real(dp), intent(out) :: dtt_dxi, dtt_deta
1786 
1787 integer(i4b) :: kc
1788 
1789 !-------- Initialisation --------
1790 
1791 ct1 = 0.0_dp
1792 ct2 = 0.0_dp
1793 ct3 = 0.0_dp
1794 ct4 = 0.0_dp
1795 ce5 = 0.0_dp
1796 ctr1 = 0.0_dp
1797 clb1 = 0.0_dp
1798 ct1_sg = 0.0_dp
1799 ct2_sg = 0.0_dp
1800 ct3_sg = 0.0_dp
1801 ct4_sg = 0.0_dp
1802 adv_vert_sg = 0.0_dp
1803 abs_adv_vert_sg = 0.0_dp
1804 ci1 = 0.0_dp
1805 cqtlde = 0.0_dp
1806 dtt_dxi = 0.0_dp
1807 dtt_deta = 0.0_dp
1808 
1809 !-------- Actual computation --------
1810 
1811 ctr1 = atr1
1812 clb1 = alb1*q_geo(j,i)
1813 
1814 #if (ADV_VERT==1)
1815 
1816 do kc=1, kcmax-1
1817  ct1(kc) = at1(kc)/h_c(j,i)*0.5_dp*(vz_c(kc,j,i)+vz_c(kc-1,j,i))
1818 end do
1819 
1820 kc=0
1821 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
1822  ! only needed for kc=0 ...
1823 kc=kcmax-1
1824 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
1825  ! ... and kc=KCMAX-1
1826 
1827 #elif (ADV_VERT==2 || ADV_VERT==3)
1828 
1829 do kc=0, kcmax-1
1830  ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
1831 end do
1832 
1833 #endif
1834 
1835 do kc=0, kcmax
1836 
1837  ct2(kc) = ( at2_1(kc)*dzm_dtau(j,i) &
1838  +at2_2(kc)*dh_c_dtau(j,i) )/h_c(j,i)
1839  ct3(kc) = ( at3_1(kc)*dzm_dxi_g(j,i) &
1840  +at3_2(kc)*dh_c_dxi_g(j,i) )/h_c(j,i) &
1841  *0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1)) *insq_g11_g(j,i)
1842  ct4(kc) = ( at4_1(kc)*dzm_deta_g(j,i) &
1843  +at4_2(kc)*dh_c_deta_g(j,i) )/h_c(j,i) &
1844  *0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i)) *insq_g22_g(j,i)
1845  ce5(kc) = at5(kc)/h_c(j,i)
1846  ci1(kc) = ai1(kc)/h_c(j,i)
1847  cqtlde(kc) = aqtlde(kc)*h_c(j,i)
1848 
1849 end do
1850 
1851 #if (ADV_VERT==1)
1852 
1853 kc=0
1854 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
1855 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
1856 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
1857 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
1858 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! only needed for kc=0 ...
1859 kc=kcmax-1
1860 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
1861 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
1862 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
1863 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
1864 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! ... and kc=KCMAX-1
1865 
1866 #elif (ADV_VERT==2 || ADV_VERT==3)
1867 
1868 do kc=0, kcmax-1
1869  ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
1870  ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
1871  ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
1872  adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
1873  abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
1874 end do
1875 
1876 #endif
1877 
1878 #if (ADV_HOR==3)
1879 dtt_dxi = 2.0_dp*dtt_2dxi
1880 dtt_deta = 2.0_dp*dtt_2deta
1881 #endif
1882 
1883 end subroutine calc_temp_enth_2_a1
1884 
1885 !-------------------------------------------------------------------------------
1886 !> Computation of temperature and age for an ice column with a temperate base
1887 !! with the enthalpy method: Abbreviations II.
1888 !<------------------------------------------------------------------------------
1889 subroutine calc_temp_enth_2_a2(at6, at7, ai2, am3, temp_c_val, omega_c_val, &
1890  i, j, ce6, ce7, ci2, cm3)
1891 
1893  creep, viscosity
1894 
1895 implicit none
1896 
1897 #ifndef ALLOW_OPENAD /* Normal */
1898 integer(i4b), intent(in) :: i, j
1899 #else /* OpenAD */
1900 integer(i4b), intent(inout) :: i, j
1901 #endif /* Normal vs. OpenAD */
1902 
1903 real(dp), intent(in) :: at6(0:kcmax), at7, ai2(0:kcmax), am3(0:kcmax)
1904 real(dp), intent(in) :: temp_c_val(0:kcmax), omega_c_val(0:kcmax)
1905 
1906 real(dp), intent(out) :: ce6(0:kcmax), ce7(0:kcmax), ci2(0:kcmax), &
1907  cm3(0:kcmax)
1908 
1909 integer(i4b) :: kc
1910 real(dp) :: temp_c_help(0:kcmax)
1911 
1912 !-------- Initialisation --------
1913 
1914 ce6 = 0.0_dp
1915 ce7 = 0.0_dp
1916 ci2 = 0.0_dp
1917 cm3 = 0.0_dp
1918 
1919 !-------- Actual computation --------
1920 
1921 do kc=0, kcmax
1922 
1923 #if (DYNAMICS==2)
1924  if (.not.flag_shelfy_stream(j,i)) then
1925 #endif
1926  ce7(kc) = at7 &
1927  *enh_c(kc,j,i) &
1928  *ratefac_c_t(temp_c_val(kc), omega_c_val(kc), temp_c_m(kc,j,i)) &
1929  *creep(sigma_c(kc,j,i)) &
1930  *sigma_c(kc,j,i)*sigma_c(kc,j,i)
1931 #if (DYNAMICS==2)
1932  else
1933  ce7(kc) = 2.0_dp*at7 &
1934  *viscosity(de_c(kc,j,i), &
1935  temp_c_val(kc), temp_c_m(kc,j,i), omega_c_val(kc), &
1936  enh_c(kc,j,i), 2_i1b) &
1937  *de_c(kc,j,i)**2
1938  end if
1939 #endif
1940 
1941  cm3(kc) = am3(kc)*h_c(j,i)*c_val(temp_c_val(kc))
1942 
1943 end do
1944 
1945 do kc=0, kc_cts_neu(j,i)-1 ! temperate layer
1946  ce6(kc) = at6(kc) &
1947  *nue/h_c(j,i)
1948  ci2(kc) = ai2(kc)/h_c(j,i)
1949 end do
1950 
1951 do kc=kc_cts_neu(j,i), kcmax-1 ! cold layer
1952  temp_c_help(kc) = 0.5_dp*(temp_c_val(kc)+temp_c_val(kc+1))
1953  ce6(kc) = at6(kc) &
1954  *kappa_val(temp_c_help(kc))/(c_val(temp_c_help(kc))*h_c(j,i))
1955  ci2(kc) = ai2(kc)/h_c(j,i)
1956 end do
1957 
1958 end subroutine calc_temp_enth_2_a2
1959 
1960 !-------------------------------------------------------------------------------
1961 !> Computation of temperature and age for an ice column with a temperate base
1962 !! with the enthalpy method:
1963 !! Set-up of the equations for the bedrock temperature.
1964 !<------------------------------------------------------------------------------
1965 subroutine calc_temp_enth_2_b(ctr1, clb1, i, j, &
1966  lgs_a0, lgs_a1, lgs_a2, lgs_b)
1968 implicit none
1969 
1970 #ifndef ALLOW_OPENAD /* Normal */
1971 integer(i4b), intent(in) :: i, j
1972 #else /* OpenAD */
1973 integer(i4b), intent(inout) :: i, j
1974 #endif /* Normal vs. OpenAD */
1975 
1976 real(dp), intent(in) :: ctr1, clb1
1977 
1978 real(dp), intent(out) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
1979  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
1980  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
1981  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
1982 
1983 integer(i4b) :: kr
1984 
1985 !-------- Initialisation --------
1986 
1987 lgs_a0 = 0.0_dp
1988 lgs_a1 = 0.0_dp
1989 lgs_a2 = 0.0_dp
1990 lgs_b = 0.0_dp
1991 
1992 !-------- Actual computation --------
1993 
1994 kr=0
1995 lgs_a1(kr) = 1.0_dp
1996 lgs_a2(kr) = -1.0_dp
1997 lgs_b(kr) = clb1
1998 
1999 #if (Q_LITHO==1)
2000 ! (coupled heat-conducting bedrock)
2001 
2002 do kr=1, krmax-1
2003  lgs_a0(kr) = - ctr1
2004  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
2005  lgs_a2(kr) = - ctr1
2006  lgs_b(kr) = temp_r(kr,j,i)
2007 end do
2008 
2009 #elif (Q_LITHO==0)
2010 ! (no coupled heat-conducting bedrock)
2011 
2012 do kr=1, krmax-1
2013  lgs_a0(kr) = 1.0_dp
2014  lgs_a1(kr) = 0.0_dp
2015  lgs_a2(kr) = -1.0_dp
2016  lgs_b(kr) = 2.0_dp*clb1
2017 end do
2018 
2019 #endif
2020 
2021 kr=krmax
2022 lgs_a0(kr) = 0.0_dp
2023 lgs_a1(kr) = 1.0_dp
2024 lgs_b(kr) = temp_t_m(0,j,i)
2025 
2026 end subroutine calc_temp_enth_2_b
2027 
2028 !-------------------------------------------------------------------------------
2029 !> Computation of temperature and age for an ice column with a temperate base
2030 !! with the enthalpy method:
2031 !! Set-up of the equations for the ice enthalpy.
2032 !<------------------------------------------------------------------------------
2033 subroutine calc_temp_enth_2_c(ct1, ct2, ct3, ct4, ce5, ce6, ce7, &
2034  ct1_sg, ct2_sg, ct3_sg, ct4_sg, cm3, &
2035  adv_vert_sg, abs_adv_vert_sg, &
2036  dtime_temp, dtt_dxi, dtt_deta, &
2037  dtt_2dxi, dtt_2deta, &
2038  i, j, kcmin, &
2039  lgs_a0, lgs_a1, lgs_a2, lgs_b)
2040 
2042 
2043 implicit none
2044 
2045 #ifndef ALLOW_OPENAD /* Normal */
2046 integer(i4b), intent(in) :: i, j
2047 #else /* OpenAD */
2048 integer(i4b), intent(inout) :: i, j
2049 #endif /* Normal vs. OpenAD */
2050 
2051 integer(i4b), intent(in) :: kcmin
2052 real(dp), intent(in) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), &
2053  ct4(0:kcmax), ce5(0:kcmax), ce6(0:kcmax), &
2054  ce7(0:kcmax)
2055 real(dp), intent(in) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), &
2056  ct3_sg(0:kcmax), ct4_sg(0:kcmax), &
2057  adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
2058 real(dp), intent(in) :: cm3(0:kcmax)
2059 real(dp), intent(in) :: dtime_temp, dtt_dxi, dtt_deta, dtt_2dxi, dtt_2deta
2060 
2061 real(dp), intent(out) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
2062  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
2063  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
2064  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
2065 
2066 integer(i4b) :: kc
2067 real(dp) :: vx_c_help, vy_c_help
2068 real(dp) :: adv_vert_help
2069 
2070 !-------- Initialisation --------
2071 
2072 lgs_a0 = 0.0_dp
2073 lgs_a1 = 0.0_dp
2074 lgs_a2 = 0.0_dp
2075 lgs_b = 0.0_dp
2076 
2077 !-------- Actual computation --------
2078 
2079 if (kcmin == 0) then ! predictor step
2080 
2081  kc=0
2082 
2083  if (kc_cts_neu(j,i) == 0) then ! temperate base without temperate layer
2084 
2085  lgs_a1(kc) = 1.0_dp
2086  lgs_a2(kc) = 0.0_dp
2087  lgs_b(kc) = enth_fct_temp_omega(temp_c_m(kc,j,i), 0.0_dp)
2088 
2089  else ! kc_cts_neu(j,i) > 0, temperate base with temperate layer
2090 
2091  lgs_a1(kc) = 1.0_dp
2092  lgs_a2(kc) = -1.0_dp
2093  lgs_b(kc) = 0.0_dp
2094 
2095  end if
2096 
2097 else ! kcmin > 0, corrector step
2098 
2099  kc=0
2100 
2101  lgs_a1(kc) = 1.0_dp ! dummy
2102  lgs_a2(kc) = 0.0_dp ! setting,
2103  lgs_b(kc) = 0.0_dp ! not needed
2104 
2105  do kc=1, kcmin-1
2106 
2107  lgs_a0(kc) = 0.0_dp ! dummy
2108  lgs_a1(kc) = 1.0_dp ! setting,
2109  lgs_a2(kc) = 0.0_dp ! not
2110  lgs_b(kc) = 0.0_dp ! needed
2111 
2112  end do
2113 
2114  kc=kcmin
2115 
2116  lgs_a0(kc) = 0.0_dp
2117  lgs_a1(kc) = 1.0_dp
2118  lgs_a2(kc) = -1.0_dp
2119  lgs_b(kc) = -cm3(kc)
2120 
2121 end if
2122 
2123 do kc=kcmin+1, kcmax-1
2124 
2125 #if (ADV_VERT==1)
2126 
2127  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2128  -ce5(kc)*ce6(kc-1)
2129  lgs_a1(kc) = 1.0_dp+ce5(kc)*(ce6(kc)+ce6(kc-1))
2130  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2131  -ce5(kc)*ce6(kc)
2132 
2133 #elif (ADV_VERT==2)
2134 
2135  lgs_a0(kc) &
2136  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2137  -ce5(kc)*ce6(kc-1)
2138  lgs_a1(kc) &
2139  = 1.0_dp &
2140  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2141  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
2142  +ce5(kc)*(ce6(kc)+ce6(kc-1))
2143  lgs_a2(kc) &
2144  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
2145  -ce5(kc)*ce6(kc)
2146 
2147 #elif (ADV_VERT==3)
2148 
2149  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
2150 
2151  lgs_a0(kc) &
2152  = -max(adv_vert_help, 0.0_dp) &
2153  -ce5(kc)*ce6(kc-1)
2154  lgs_a1(kc) &
2155  = 1.0_dp &
2156  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
2157  +ce5(kc)*(ce6(kc)+ce6(kc-1))
2158  lgs_a2(kc) &
2159  = min(adv_vert_help, 0.0_dp) &
2160  -ce5(kc)*ce6(kc)
2161 
2162 #endif
2163 
2164 #if (ADV_HOR==2)
2165 
2166  lgs_b(kc) = enth_c(kc,j,i) + ce7(kc) &
2167  -dtt_2dxi* &
2168  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2169  *(enth_c(kc,j,i+1)-enth_c(kc,j,i)) &
2170  *insq_g11_sgx(j,i) &
2171  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2172  *(enth_c(kc,j,i)-enth_c(kc,j,i-1)) &
2173  *insq_g11_sgx(j,i-1) ) &
2174  -dtt_2deta* &
2175  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2176  *(enth_c(kc,j+1,i)-enth_c(kc,j,i)) &
2177  *insq_g22_sgy(j,i) &
2178  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2179  *(enth_c(kc,j,i)-enth_c(kc,j-1,i)) &
2180  *insq_g22_sgy(j-1,i) )
2181 
2182 #elif (ADV_HOR==3)
2183 
2184  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2185  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2186 
2187  lgs_b(kc) = enth_c(kc,j,i) + ce7(kc) &
2188  -dtt_dxi* &
2189  ( min(vx_c_help, 0.0_dp) &
2190  *(enth_c(kc,j,i+1)-enth_c(kc,j,i)) &
2191  *insq_g11_sgx(j,i) &
2192  +max(vx_c_help, 0.0_dp) &
2193  *(enth_c(kc,j,i)-enth_c(kc,j,i-1)) &
2194  *insq_g11_sgx(j,i-1) ) &
2195  -dtt_deta* &
2196  ( min(vy_c_help, 0.0_dp) &
2197  *(enth_c(kc,j+1,i)-enth_c(kc,j,i)) &
2198  *insq_g22_sgy(j,i) &
2199  +max(vy_c_help, 0.0_dp) &
2200  *(enth_c(kc,j,i)-enth_c(kc,j-1,i)) &
2201  *insq_g22_sgy(j-1,i) )
2202 
2203 #endif
2204 
2205 end do
2206 
2207 kc=kcmax
2208 lgs_a0(kc) = 0.0_dp
2209 lgs_a1(kc) = 1.0_dp
2210 lgs_b(kc) = enth_fct_temp_omega(temp_s(j,i), 0.0_dp)
2211  ! zero water content at the ice surface
2212 
2213 end subroutine calc_temp_enth_2_c
2214 
2215 !-------------------------------------------------------------------------------
2216 !> Computation of temperature and age for an ice column with a temperate base
2217 !! with the enthalpy method:
2218 !! Set-up of the equations for the age of ice.
2219 !<------------------------------------------------------------------------------
2220 subroutine calc_temp_enth_2_d(ct1, ct2, ct3, ct4, ci1, ci2, &
2221  ct1_sg, ct2_sg, ct3_sg, ct4_sg, &
2222  adv_vert_sg, abs_adv_vert_sg, &
2223  dtime_temp, dtt_dxi, dtt_deta, &
2224  dtt_2dxi, dtt_2deta, &
2225  i, j, &
2226  lgs_a0, lgs_a1, lgs_a2, lgs_b)
2228 implicit none
2229 
2230 #ifndef ALLOW_OPENAD /* Normal */
2231 integer(i4b), intent(in) :: i, j
2232 #else /* OpenAD */
2233 integer(i4b), intent(inout) :: i, j
2234 #endif /* Normal vs. OpenAD */
2235 
2236 real(dp), intent(in) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), &
2237  ct4(0:kcmax), ci1(0:kcmax), ci2(0:kcmax)
2238 real(dp), intent(in) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), &
2239  ct3_sg(0:kcmax), ct4_sg(0:kcmax), &
2240  adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
2241 real(dp), intent(in) :: dtime_temp, dtt_dxi, dtt_deta, dtt_2dxi, dtt_2deta
2242 
2243 real(dp), intent(out) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
2244  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
2245  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
2246  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
2247 
2248 integer(i4b) :: kc
2249 real(dp) :: vx_c_help, vy_c_help
2250 real(dp) :: adv_vert_help
2251 
2252 !-------- Initialisation --------
2253 
2254 lgs_a0 = 0.0_dp
2255 lgs_a1 = 0.0_dp
2256 lgs_a2 = 0.0_dp
2257 lgs_b = 0.0_dp
2258 
2259 !-------- Actual computation --------
2260 
2261 kc=0 ! adv_vert_sg(0) <= 0
2262 lgs_a1(kc) = 1.0_dp - min(adv_vert_sg(kc), 0.0_dp) ! (directed downward)
2263 lgs_a2(kc) = min(adv_vert_sg(kc), 0.0_dp) ! assumed/enforced
2264 
2265 #if (ADV_HOR==2)
2266 
2267 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2268  -dtt_2dxi* &
2269  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2270  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2271  *insq_g11_sgx(j,i) &
2272  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2273  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2274  *insq_g11_sgx(j,i-1) ) &
2275  -dtt_2deta* &
2276  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2277  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2278  *insq_g22_sgy(j,i) &
2279  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2280  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2281  *insq_g22_sgy(j-1,i) )
2282 
2283 #elif (ADV_HOR==3)
2284 
2285 vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2286 vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2287 
2288 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2289  -dtt_dxi* &
2290  ( min(vx_c_help, 0.0_dp) &
2291  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2292  *insq_g11_sgx(j,i) &
2293  +max(vx_c_help, 0.0_dp) &
2294  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2295  *insq_g11_sgx(j,i-1) ) &
2296  -dtt_deta* &
2297  ( min(vy_c_help, 0.0_dp) &
2298  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2299  *insq_g22_sgy(j,i) &
2300  +max(vy_c_help, 0.0_dp) &
2301  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2302  *insq_g22_sgy(j-1,i) )
2303 
2304 #endif
2305 
2306 do kc=1, kcmax-1
2307 
2308 #if (ADV_VERT==1)
2309 
2310  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2311  -ci1(kc)*ci2(kc-1)
2312  lgs_a1(kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
2313  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2314  -ci1(kc)*ci2(kc)
2315 
2316 #elif (ADV_VERT==2)
2317 
2318  lgs_a0(kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
2319  lgs_a1(kc) = 1.0_dp &
2320  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2321  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
2322  lgs_a2(kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
2323 
2324 #elif (ADV_VERT==3)
2325 
2326  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
2327 
2328  lgs_a0(kc) = -max(adv_vert_help, 0.0_dp)
2329  lgs_a1(kc) = 1.0_dp &
2330  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp)
2331  lgs_a2(kc) = min(adv_vert_help, 0.0_dp)
2332 
2333 #endif
2334 
2335 #if (ADV_HOR==2)
2336 
2337  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2338  -dtt_2dxi* &
2339  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2340  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2341  *insq_g11_sgx(j,i) &
2342  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2343  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2344  *insq_g11_sgx(j,i-1) ) &
2345  -dtt_2deta* &
2346  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2347  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2348  *insq_g22_sgy(j,i) &
2349  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2350  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2351  *insq_g22_sgy(j-1,i) )
2352 
2353 #elif (ADV_HOR==3)
2354 
2355  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2356  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2357 
2358  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2359  -dtt_dxi* &
2360  ( min(vx_c_help, 0.0_dp) &
2361  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2362  *insq_g11_sgx(j,i) &
2363  +max(vx_c_help, 0.0_dp) &
2364  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2365  *insq_g11_sgx(j,i-1) ) &
2366  -dtt_deta* &
2367  ( min(vy_c_help, 0.0_dp) &
2368  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2369  *insq_g22_sgy(j,i) &
2370  +max(vy_c_help, 0.0_dp) &
2371  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2372  *insq_g22_sgy(j-1,i) )
2373 
2374 #endif
2375 
2376 end do
2377 
2378 kc=kcmax
2379 if (as_perp(j,i) >= 0.0_dp) then
2380  lgs_a0(kc) = 0.0_dp
2381  lgs_a1(kc) = 1.0_dp
2382  lgs_b(kc) = 0.0_dp
2383 else
2384  lgs_a0(kc) = -max(adv_vert_sg(kc-1), 0.0_dp)
2385  lgs_a1(kc) = 1.0_dp + max(adv_vert_sg(kc-1), 0.0_dp)
2386  ! adv_vert_sg(KCMAX-1) >= 0 (directed upward)
2387  ! assumed/enforced
2388 #if (ADV_HOR==2)
2389 
2390  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2391  -dtt_2dxi* &
2392  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2393  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2394  *insq_g11_sgx(j,i) &
2395  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2396  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2397  *insq_g11_sgx(j,i-1) ) &
2398  -dtt_2deta* &
2399  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2400  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2401  *insq_g22_sgy(j,i) &
2402  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2403  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2404  *insq_g22_sgy(j-1,i) )
2405 
2406 #elif (ADV_HOR==3)
2407 
2408  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2409  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2410 
2411  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2412  -dtt_dxi* &
2413  ( min(vx_c_help, 0.0_dp) &
2414  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2415  *insq_g11_sgx(j,i) &
2416  +max(vx_c_help, 0.0_dp) &
2417  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2418  *insq_g11_sgx(j,i-1) ) &
2419  -dtt_deta* &
2420  ( min(vy_c_help, 0.0_dp) &
2421  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2422  *insq_g22_sgy(j,i) &
2423  +max(vy_c_help, 0.0_dp) &
2424  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2425  *insq_g22_sgy(j-1,i) )
2426 
2427 #endif
2428 
2429 end if
2430 
2431 end subroutine calc_temp_enth_2_d
2432 
2433 !-------------------------------------------------------------------------------
2434 !> Computation of temperature, age, water content and enthalpy for an
2435 !! ice-free column.
2436 !<------------------------------------------------------------------------------
2437 subroutine calc_temp_enth_r(atr1, alb1, i, j)
2438 
2439 #ifndef ALLOW_OPENAD /* Normal */
2440 use sico_maths_m, only : tri_sle
2441 #else /* OpenAD */
2442 use sico_maths_m
2443 #endif /* Normal vs. OpenAD */
2444 
2446 
2447 implicit none
2448 
2449 #ifndef ALLOW_OPENAD /* Normal */
2450 integer(i4b), intent(in) :: i, j
2451 #else /* OpenAD */
2452 integer(i4b), intent(inout) :: i, j
2453 #endif /* Normal vs. OpenAD */
2454 
2455 real(dp), intent(in) :: atr1, alb1
2456 
2457 integer(i4b) :: kc, kt, kr
2458 real(dp) :: ctr1, clb1
2459 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
2460  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
2461  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
2462  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
2463  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
2464 real(dp) :: enth_val
2465 
2466 !-------- Abbreviations --------
2467 
2468 ctr1 = atr1
2469 clb1 = alb1*q_geo(j,i)
2470 
2471 !-------- Set up the equations for the bedrock temperature --------
2472 
2473 kr=0
2474 lgs_a1(kr) = 1.0_dp
2475 lgs_a2(kr) = -1.0_dp
2476 lgs_b(kr) = clb1
2477 
2478 #if (Q_LITHO==1)
2479 ! (coupled heat-conducting bedrock)
2480 
2481 do kr=1, krmax-1
2482  lgs_a0(kr) = - ctr1
2483  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
2484  lgs_a2(kr) = - ctr1
2485  lgs_b(kr) = temp_r(kr,j,i)
2486 end do
2487 
2488 #elif (Q_LITHO==0)
2489 ! (no coupled heat-conducting bedrock)
2490 
2491 do kr=1, krmax-1
2492  lgs_a0(kr) = 1.0_dp
2493  lgs_a1(kr) = 0.0_dp
2494  lgs_a2(kr) = -1.0_dp
2495  lgs_b(kr) = 2.0_dp*clb1
2496 end do
2497 
2498 #endif
2499 
2500 kr=krmax
2501 lgs_a0(kr) = 0.0_dp
2502 lgs_a1(kr) = 1.0_dp
2503 lgs_b(kr) = temp_s(j,i)
2504 
2505 !-------- Solve system of linear equations --------
2506 
2507 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
2508 
2509 !-------- Assign the result --------
2510 
2511 do kr=0, krmax
2512  temp_r_neu(kr,j,i) = lgs_x(kr)
2513 end do
2514 
2515 !-------- Water content, age and enthalpy
2516 ! in the non-existing lower (kt) ice layer --------
2517 
2518 enth_val = enth_fct_temp_omega(temp_s(j,i), 0.0_dp)
2519 
2520 do kt=0, ktmax
2521  omega_t_neu(kt,j,i) = 0.0_dp
2522  age_t_neu(kt,j,i) = 0.0_dp
2523  enth_t_neu(kt,j,i) = enth_val
2524 end do
2525 
2526 !-------- Temperature, age, water content and enthalpy
2527 ! in the non-existing upper (kc) ice layer --------
2528 
2529 do kc=0, kcmax
2530  temp_c_neu(kc,j,i) = temp_s(j,i)
2531  age_c_neu(kc,j,i) = 0.0_dp
2532  omega_c_neu(kc,j,i) = 0.0_dp
2533  enth_c_neu(kc,j,i) = enth_val
2534 end do
2535 
2536 end subroutine calc_temp_enth_r
2537 
2538 !-------------------------------------------------------------------------------
2539 !> Computation of temperature and age for ice shelves (floating ice)
2540 !! with the enthalpy method.
2541 !<------------------------------------------------------------------------------
2542 subroutine calc_temp_enth_ssa(at1, at2_1, at2_2, at3_1, at3_2, &
2543  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
2544  ai1, ai2, &
2545  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
2548 
2549 #ifndef ALLOW_OPENAD /* Normal */
2550 use sico_maths_m, only : tri_sle
2551 #else /* OpenAD */
2552 use sico_maths_m
2553 #endif /* Normal vs. OpenAD */
2554 
2557 
2558 implicit none
2559 
2560 #ifndef ALLOW_OPENAD /* Normal */
2561 integer(i4b), intent(in) :: i, j
2562 #else /* OpenAD */
2563 integer(i4b), intent(inout) :: i, j
2564 #endif /* Normal vs. OpenAD */
2565 
2566 real(dp), intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
2567  at3_1(0:kcmax), at3_2(0:kcmax), &
2568  at4_1(0:kcmax), at4_2(0:kcmax), &
2569  at5(0:kcmax), at6(0:kcmax), at7, &
2570  ai1(0:kcmax), ai2(0:kcmax), &
2571  atr1, alb1
2572 real(dp), intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
2573 
2574 integer(i4b) :: kc, kt, kr
2575 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
2576  ce5(0:kcmax), ce6(0:kcmax), ce7(0:kcmax), ctr1, clb1
2577 real(dp) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), ct3_sg(0:kcmax), &
2578  ct4_sg(0:kcmax), adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
2579 real(dp) :: ci1(0:kcmax), ci2(0:kcmax)
2580 real(dp) :: temp_c_help(0:kcmax)
2581 real(dp) :: vx_c_help, vy_c_help
2582 real(dp) :: adv_vert_help
2583 real(dp) :: dtt_dxi, dtt_deta
2584 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
2585  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
2586  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
2587  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
2588  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
2589 real(dp), parameter :: zero=0.0_dp
2590 
2591 !-------- Check for boundary points --------
2592 
2593 if ((i == 0).or.(i == imax).or.(j == 0).or.(j == jmax)) then
2594  errormsg = ' >>> calc_temp_enth_ssa: Boundary points not allowed.'
2595  call error(errormsg)
2596 end if
2597 
2598 !-------- Abbreviations --------
2599 
2600 ctr1 = atr1
2601 clb1 = alb1*q_geo(j,i)
2602 
2603 #if (ADV_VERT==1)
2604 
2605 do kc=1, kcmax-1
2606  ct1(kc) = at1(kc)/h_c(j,i)*0.5_dp*(vz_c(kc,j,i)+vz_c(kc-1,j,i))
2607 end do
2608 
2609 kc=0
2610 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
2611  ! only needed for kc=0 ...
2612 kc=kcmax-1
2613 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
2614  ! ... and kc=KCMAX-1
2615 
2616 #elif (ADV_VERT==2 || ADV_VERT==3)
2617 
2618 do kc=0, kcmax-1
2619  ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
2620 end do
2621 
2622 #endif
2623 
2624 do kc=0, kcmax
2625 
2626  ct2(kc) = ( at2_1(kc)*dzm_dtau(j,i) &
2627  +at2_2(kc)*dh_c_dtau(j,i) )/h_c(j,i)
2628  ct3(kc) = ( at3_1(kc)*dzm_dxi_g(j,i) &
2629  +at3_2(kc)*dh_c_dxi_g(j,i) )/h_c(j,i) &
2630  *0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1)) *insq_g11_g(j,i)
2631  ct4(kc) = ( at4_1(kc)*dzm_deta_g(j,i) &
2632  +at4_2(kc)*dh_c_deta_g(j,i) )/h_c(j,i) &
2633  *0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i)) *insq_g22_g(j,i)
2634  ce5(kc) = at5(kc)/h_c(j,i)
2635  ce7(kc) = 2.0_dp*at7 &
2636  *viscosity(de_ssa(j,i), &
2637  temp_c(kc,j,i), temp_c_m(kc,j,i), 0.0_dp, &
2638 #if (!defined(ALLOW_GRDCHK) && !defined(ALLOW_OPENAD))
2639  enh_c(kc,j,i), 0_i1b) &
2640 #else
2641  enh_c(kc,j,i), 0_i4b) &
2642 #endif
2643  *de_ssa(j,i)**2
2644  ci1(kc) = ai1(kc)/h_c(j,i)
2645 
2646 end do
2647 
2648 #if (ADV_VERT==1)
2649 
2650 kc=0
2651 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2652 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2653 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2654 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2655 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! only needed for kc=0 ...
2656 kc=kcmax-1
2657 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2658 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2659 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2660 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2661 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc)) ! ... and kc=KCMAX-1
2662 
2663 #elif (ADV_VERT==2 || ADV_VERT==3)
2664 
2665 do kc=0, kcmax-1
2666  ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2667  ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2668  ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2669  adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2670  abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
2671 end do
2672 
2673 #endif
2674 
2675 do kc=0, kcmax-1
2676  temp_c_help(kc) = 0.5_dp*(temp_c(kc,j,i)+temp_c(kc+1,j,i))
2677  ce6(kc) = at6(kc) &
2678  *kappa_val(temp_c_help(kc))/(c_val(temp_c_help(kc))*h_c(j,i))
2679  ci2(kc) = ai2(kc)/h_c(j,i)
2680 end do
2681 
2682 #if (ADV_HOR==3)
2683 dtt_dxi = 2.0_dp*dtt_2dxi
2684 dtt_deta = 2.0_dp*dtt_2deta
2685 #endif
2686 
2687 !-------- Set up the equations for the bedrock temperature --------
2688 
2689 kr=0
2690 lgs_a1(kr) = 1.0_dp
2691 lgs_a2(kr) = -1.0_dp
2692 lgs_b(kr) = clb1
2693 
2694 #if (Q_LITHO==1)
2695 ! (coupled heat-conducting bedrock)
2696 
2697 do kr=1, krmax-1
2698  lgs_a0(kr) = - ctr1
2699  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
2700  lgs_a2(kr) = - ctr1
2701  lgs_b(kr) = temp_r(kr,j,i)
2702 end do
2703 
2704 #elif (Q_LITHO==0)
2705 ! (no coupled heat-conducting bedrock)
2706 
2707 do kr=1, krmax-1
2708  lgs_a0(kr) = 1.0_dp
2709  lgs_a1(kr) = 0.0_dp
2710  lgs_a2(kr) = -1.0_dp
2711  lgs_b(kr) = 2.0_dp*clb1
2712 end do
2713 
2714 #endif
2715 
2716 kr=krmax
2717 lgs_a0(kr) = 0.0_dp
2718 lgs_a1(kr) = 1.0_dp
2719 lgs_b(kr) = temp_c_m(0,j,i)-delta_tm_sw
2720 
2721 !-------- Solve system of linear equations --------
2722 
2723 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
2724 
2725 !-------- Assign the result --------
2726 
2727 do kr=0, krmax
2728  temp_r_neu(kr,j,i) = lgs_x(kr)
2729 end do
2730 
2731 !-------- Set up the equations for the ice temperature --------
2732 
2733 kc=0
2734 lgs_a1(kc) = 1.0_dp
2735 lgs_a2(kc) = 0.0_dp
2736 lgs_b(kc) = enth_fct_temp_omega(temp_c_m(kc,j,i)-delta_tm_sw, 0.0_dp)
2737  ! zero water content assumed
2738 
2739 do kc=1, kcmax-1
2740 
2741 #if (ADV_VERT==1)
2742 
2743  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2744  -ce5(kc)*ce6(kc-1)
2745  lgs_a1(kc) = 1.0_dp+ce5(kc)*(ce6(kc)+ce6(kc-1))
2746  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2747  -ce5(kc)*ce6(kc)
2748 
2749 #elif (ADV_VERT==2)
2750 
2751  lgs_a0(kc) &
2752  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2753  -ce5(kc)*ce6(kc-1)
2754  lgs_a1(kc) &
2755  = 1.0_dp &
2756  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2757  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
2758  +ce5(kc)*(ce6(kc)+ce6(kc-1))
2759  lgs_a2(kc) &
2760  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
2761  -ce5(kc)*ce6(kc)
2762 
2763 #elif (ADV_VERT==3)
2764 
2765  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
2766 
2767  lgs_a0(kc) &
2768  = -max(adv_vert_help, 0.0_dp) &
2769  -ce5(kc)*ce6(kc-1)
2770  lgs_a1(kc) &
2771  = 1.0_dp &
2772  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
2773  +ce5(kc)*(ce6(kc)+ce6(kc-1))
2774  lgs_a2(kc) &
2775  = min(adv_vert_help, 0.0_dp) &
2776  -ce5(kc)*ce6(kc)
2777 
2778 #endif
2779 
2780 #if (ADV_HOR==2)
2781 
2782  lgs_b(kc) = enth_c(kc,j,i) + ce7(kc) &
2783  -dtt_2dxi* &
2784  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2785  *(enth_c(kc,j,i+1)-enth_c(kc,j,i)) &
2786  *insq_g11_sgx(j,i) &
2787  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2788  *(enth_c(kc,j,i)-enth_c(kc,j,i-1)) &
2789  *insq_g11_sgx(j,i-1) ) &
2790  -dtt_2deta* &
2791  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2792  *(enth_c(kc,j+1,i)-enth_c(kc,j,i)) &
2793  *insq_g22_sgy(j,i) &
2794  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2795  *(enth_c(kc,j,i)-enth_c(kc,j-1,i)) &
2796  *insq_g22_sgy(j-1,i) )
2797 
2798 #elif (ADV_HOR==3)
2799 
2800  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2801  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2802 
2803  lgs_b(kc) = enth_c(kc,j,i) + ce7(kc) &
2804  -dtt_dxi* &
2805  ( min(vx_c_help, 0.0_dp) &
2806  *(enth_c(kc,j,i+1)-enth_c(kc,j,i)) &
2807  *insq_g11_sgx(j,i) &
2808  +max(vx_c_help, 0.0_dp) &
2809  *(enth_c(kc,j,i)-enth_c(kc,j,i-1)) &
2810  *insq_g11_sgx(j,i-1) ) &
2811  -dtt_deta* &
2812  ( min(vy_c_help, 0.0_dp) &
2813  *(enth_c(kc,j+1,i)-enth_c(kc,j,i)) &
2814  *insq_g22_sgy(j,i) &
2815  +max(vy_c_help, 0.0_dp) &
2816  *(enth_c(kc,j,i)-enth_c(kc,j-1,i)) &
2817  *insq_g22_sgy(j-1,i) )
2818 
2819 #endif
2820 
2821 end do
2822 
2823 kc=kcmax
2824 lgs_a0(kc) = 0.0_dp
2825 lgs_a1(kc) = 1.0_dp
2826 lgs_b(kc) = enth_fct_temp_omega(temp_s(j,i), 0.0_dp)
2827  ! zero water content assumed
2828 
2829 !-------- Solve system of linear equations --------
2830 
2831 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
2832 
2833 !-------- Assign the result --------
2834 
2835 do kc=0, kcmax
2836  enth_c_neu(kc,j,i) = lgs_x(kc)
2837  temp_c_neu(kc,j,i) = temp_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
2838  omega_c_neu(kc,j,i) = omega_fct_enth(enth_c_neu(kc,j,i), temp_c_m(kc,j,i))
2839 end do
2840 
2841 !-------- Set enthalpy and water content in the redundant,
2842 ! lower (kt) ice layer to the value at the ice base --------
2843 
2844 do kt=0, ktmax
2845  enth_t_neu(kt,j,i) = enth_c_neu(0,j,i)
2846  omega_t_neu(kt,j,i) = omega_c_neu(0,j,i)
2847 end do
2848 
2849 !-------- Water drainage from the non-existing temperate ice --------
2850 
2851 q_tld(j,i) = 0.0_dp
2852 
2853 !-------- Set up the equations for the age of cold ice --------
2854 
2855 kc=0 ! adv_vert_sg(0) <= 0
2856 lgs_a1(kc) = 1.0_dp - min(adv_vert_sg(kc), 0.0_dp) ! (directed downward)
2857 lgs_a2(kc) = min(adv_vert_sg(kc), 0.0_dp) ! assumed/enforced
2858 
2859 #if (ADV_HOR==2)
2860 
2861 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2862  -dtt_2dxi* &
2863  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2864  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2865  *insq_g11_sgx(j,i) &
2866  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2867  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2868  *insq_g11_sgx(j,i-1) ) &
2869  -dtt_2deta* &
2870  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2871  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2872  *insq_g22_sgy(j,i) &
2873  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2874  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2875  *insq_g22_sgy(j-1,i) )
2876 
2877 #elif (ADV_HOR==3)
2878 
2879 vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2880 vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2881 
2882 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2883  -dtt_dxi* &
2884  ( min(vx_c_help, 0.0_dp) &
2885  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2886  *insq_g11_sgx(j,i) &
2887  +max(vx_c_help, 0.0_dp) &
2888  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2889  *insq_g11_sgx(j,i-1) ) &
2890  -dtt_deta* &
2891  ( min(vy_c_help, 0.0_dp) &
2892  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2893  *insq_g22_sgy(j,i) &
2894  +max(vy_c_help, 0.0_dp) &
2895  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2896  *insq_g22_sgy(j-1,i) )
2897 
2898 #endif
2899 
2900 do kc=1, kcmax-1
2901 
2902 #if (ADV_VERT==1)
2903 
2904  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2905  -ci1(kc)*ci2(kc-1)
2906  lgs_a1(kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
2907  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2908  -ci1(kc)*ci2(kc)
2909 
2910 #elif (ADV_VERT==2)
2911 
2912  lgs_a0(kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
2913  lgs_a1(kc) = 1.0_dp &
2914  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2915  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
2916  lgs_a2(kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
2917 
2918 #elif (ADV_VERT==3)
2919 
2920  adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
2921 
2922  lgs_a0(kc) = -max(adv_vert_help, 0.0_dp)
2923  lgs_a1(kc) = 1.0_dp &
2924  +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp)
2925  lgs_a2(kc) = min(adv_vert_help, 0.0_dp)
2926 
2927 #endif
2928 
2929 #if (ADV_HOR==2)
2930 
2931  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2932  -dtt_2dxi* &
2933  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2934  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2935  *insq_g11_sgx(j,i) &
2936  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2937  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2938  *insq_g11_sgx(j,i-1) ) &
2939  -dtt_2deta* &
2940  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2941  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2942  *insq_g22_sgy(j,i) &
2943  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2944  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2945  *insq_g22_sgy(j-1,i) )
2946 
2947 #elif (ADV_HOR==3)
2948 
2949  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
2950  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
2951 
2952  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2953  -dtt_dxi* &
2954  ( min(vx_c_help, 0.0_dp) &
2955  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2956  *insq_g11_sgx(j,i) &
2957  +max(vx_c_help, 0.0_dp) &
2958  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2959  *insq_g11_sgx(j,i-1) ) &
2960  -dtt_deta* &
2961  ( min(vy_c_help, 0.0_dp) &
2962  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2963  *insq_g22_sgy(j,i) &
2964  +max(vy_c_help, 0.0_dp) &
2965  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2966  *insq_g22_sgy(j-1,i) )
2967 
2968 #endif
2969 
2970 end do
2971 
2972 kc=kcmax
2973 if (as_perp(j,i) >= zero) then
2974  lgs_a0(kc) = 0.0_dp
2975  lgs_a1(kc) = 1.0_dp
2976  lgs_b(kc) = 0.0_dp
2977 else
2978  lgs_a0(kc) = -max(adv_vert_sg(kc-1), 0.0_dp)
2979  lgs_a1(kc) = 1.0_dp + max(adv_vert_sg(kc-1), 0.0_dp)
2980  ! adv_vert_sg(KCMAX-1) >= 0 (directed upward)
2981  ! assumed/enforced
2982 #if (ADV_HOR==2)
2983 
2984  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
2985  -dtt_2dxi* &
2986  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
2987  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
2988  *insq_g11_sgx(j,i) &
2989  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
2990  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
2991  *insq_g11_sgx(j,i-1) ) &
2992  -dtt_2deta* &
2993  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
2994  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
2995  *insq_g22_sgy(j,i) &
2996  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
2997  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
2998  *insq_g22_sgy(j-1,i) )
2999 
3000 #elif (ADV_HOR==3)
3001 
3002  vx_c_help = 0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1))
3003  vy_c_help = 0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i))
3004 
3005  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
3006  -dtt_dxi* &
3007  ( min(vx_c_help, 0.0_dp) &
3008  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
3009  *insq_g11_sgx(j,i) &
3010  +max(vx_c_help, 0.0_dp) &
3011  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
3012  *insq_g11_sgx(j,i-1) ) &
3013  -dtt_deta* &
3014  ( min(vy_c_help, 0.0_dp) &
3015  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
3016  *insq_g22_sgy(j,i) &
3017  +max(vy_c_help, 0.0_dp) &
3018  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
3019  *insq_g22_sgy(j-1,i) )
3020 
3021 #endif
3022 
3023 end if
3024 
3025 !-------- Solve system of linear equations --------
3026 
3027 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
3028 
3029 !-------- Assign the result,
3030 ! restriction to interval [0, AGE_MAX yr] --------
3031 
3032 do kc=0, kcmax
3033 
3034  age_c_neu(kc,j,i) = lgs_x(kc)
3035 
3036  if (age_c_neu(kc,j,i) < (age_min*year_sec)) &
3037  age_c_neu(kc,j,i) = 0.0_dp
3038  if (age_c_neu(kc,j,i) > (age_max*year_sec)) &
3039  age_c_neu(kc,j,i) = age_max*year_sec
3040 
3041 end do
3042 
3043 !-------- Age of the ice in the redundant, lower (kt) ice layer --------
3044 
3045 do kt=0, ktmax
3046  age_t_neu(kt,j,i) = age_c_neu(0,j,i)
3047 end do
3048 
3049 end subroutine calc_temp_enth_ssa
3050 
3051 !-------------------------------------------------------------------------------
3052 
3053 end module calc_temp_enth_m
3054 !
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.
integer(i1b), dimension(0:jmax, 0:imax) n_cts
n_cts(j,i): Mask for thermal conditions. -1: cold ice base, 0: temperate ice base with cold ice above...
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
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(.).
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.
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: ...
integer(i4b), 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) 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: ...
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)
integer(i1b), dimension(0:jmax, 0:imax) n_cts_neu
(.)_neu: New value of quantity (.) computed during an integration step
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...
integer(i4b), dimension(0:jmax, 0:imax) kc_cts_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp), dimension(0:jmax, 0:imax) h_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
integer(i1b), dimension(0:jmax, 0:imax) maske
maske(j,i): Ice-land-ocean mask. 0: grounded ice, 1: ice-free land, 2: ocean, 3: floating ice ...
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.
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)) ...