SICOPOLIS V3.2
 All Classes Files Functions Variables Macros
calc_temp_poly.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : c a l c _ t e m p _ p o l y
4 !
5 !> @file
6 !!
7 !! Computation of temperature, water content and age in polythermal mode.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2016 Ralf Greve
12 !!
13 !! @section License
14 !!
15 !! This file is part of SICOPOLIS.
16 !!
17 !! SICOPOLIS is free software: you can redistribute it and/or modify
18 !! it under the terms of the GNU General Public License as published by
19 !! the Free Software Foundation, either version 3 of the License, or
20 !! (at your option) any later version.
21 !!
22 !! SICOPOLIS is distributed in the hope that it will be useful,
23 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
24 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 !! GNU General Public License for more details.
26 !!
27 !! You should have received a copy of the GNU General Public License
28 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
29 !<
30 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
31 
32 !-------------------------------------------------------------------------------
33 !> Computation of temperature, water content and age in polythermal mode.
34 !<------------------------------------------------------------------------------
35 subroutine calc_temp_poly(dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
36  dtime_temp)
37 
38 use sico_types
39 use sico_variables
40 use sico_vars
41 use ice_material_quantities, only : kappa_val
42 
43 implicit none
44 
45 real(dp), intent(in) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
46 real(dp), intent(in) :: dtime_temp
47 
48 integer(i4b) :: i, j, kc, kt, kr, ii, jj
49 real(dp) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
50  at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
51  at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
52  acb1, acb2, acb3, acb4, &
53  ai1(0:kcmax), ai2(0:kcmax), ai3, &
54  atr1, am1, am2, alb1
55 real(dp) :: aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld
56 real(dp) :: dtime_temp_inv, dtt_2dxi, dtt_2deta
57 real(dp) :: temp_c_help(0:kcmax)
58 real(dp) :: fact_thick
59 real(dp) :: time_lag_cts
60 real(dp) :: vol_t, vol_t_smooth, korrfakt_t, &
61  dh_t_smooth(0:jmax,0:imax)
62 
63 !-------- Term abbreviations
64 
65 at7 = 2.0_dp/rho*dtime_temp
66 
67 aw1 = dtime_temp/dzeta_t
68 aw2 = dtime_temp/dzeta_t
69 aw3 = dtime_temp/dzeta_t
70 aw4 = dtime_temp/dzeta_t
71 aw5 = nue/rho*dtime_temp/(dzeta_t**2)
72 aw7 = 2.0_dp/(rho*l)*dtime_temp
73 aw8 = beta**2/(rho*l) &
74  *(kappa_val(0.0_dp)-kappa_val(-1.0_dp))*dtime_temp
75 aw9 = beta/l*dtime_temp
76 
77 ai3 = agediff*dtime_temp/(dzeta_t**2)
78 
79 atr1 = kappa_r/(rho_c_r*h_r**2)*dtime_temp/(dzeta_r**2)
80 
81 if (flag_aa_nonzero) then
82  am1 = aa*beta*dzeta_c/(ea-1.0_dp)
83  am2 = aa*l*rho*dzeta_c/(ea-1.0_dp)
84 else
85  am1 = beta*dzeta_c
86  am2 = l*rho*dzeta_c
87 end if
88 
89 if (flag_aa_nonzero) then
90  acb1 = (ea-1.0_dp)/aa/dzeta_c
91 else
92  acb1 = 1.0_dp/dzeta_c
93 end if
94 
95 acb2 = kappa_r/h_r/dzeta_r
96 acb3 = rho*g
97 acb4 = rho*g
98 
99 alb1 = h_r/kappa_r*dzeta_r
100 
101 aqtld = dzeta_t/dtime_temp
102 
103 dtt_2dxi = 0.5_dp*dtime_temp/dxi
104 dtt_2deta = 0.5_dp*dtime_temp/deta
105 
106 dtime_temp_inv = 1.0_dp/dtime_temp
107 
108 do kc=0, kcmax
109 
110  if (flag_aa_nonzero) then
111 
112  at1(kc) = (ea-1.0_dp)/(aa*eaz_c(kc))*dtime_temp/dzeta_c
113  at2_1(kc) = (ea-1.0_dp)/(aa*eaz_c(kc))*dtime_temp/dzeta_c
114  at2_2(kc) = (eaz_c(kc)-1.0_dp)/(aa*eaz_c(kc)) &
115  *dtime_temp/dzeta_c
116  at3_1(kc) = (ea-1.0_dp)/(aa*eaz_c(kc))*dtime_temp/dzeta_c
117  at3_2(kc) = (eaz_c(kc)-1.0_dp)/(aa*eaz_c(kc)) &
118  *dtime_temp/dzeta_c
119  at4_1(kc) = (ea-1.0_dp)/(aa*eaz_c(kc))*dtime_temp/dzeta_c
120  at4_2(kc) = (eaz_c(kc)-1.0_dp)/(aa*eaz_c(kc)) &
121  *dtime_temp/dzeta_c
122  at5(kc) = (ea-1.0_dp)/(rho*aa*eaz_c(kc)) &
123  *dtime_temp/dzeta_c
124  if (kc /= kcmax) then
125  at6(kc) = (ea-1.0_dp) &
126  /(aa*exp(aa*0.5_dp*(zeta_c(kc)+zeta_c(kc+1)))) &
127  /dzeta_c
128  else
129  at6(kc) = 0.0_dp
130  end if
131  ai1(kc) = agediff*(ea-1.0_dp)/(aa*eaz_c(kc)) &
132  *dtime_temp/dzeta_c
133  if (kc /= kcmax) then
134  ai2(kc) = (ea-1.0_dp) &
135  /(aa*exp(aa*0.5_dp*(zeta_c(kc)+zeta_c(kc+1)))) &
136  /dzeta_c
137  else
138  ai2(kc) = 0.0_dp
139  end if
140 
141  else
142 
143  at1(kc) = dtime_temp/dzeta_c
144  at2_1(kc) = dtime_temp/dzeta_c
145  at2_2(kc) = zeta_c(kc) &
146  *dtime_temp/dzeta_c
147  at3_1(kc) = dtime_temp/dzeta_c
148  at3_2(kc) = zeta_c(kc) &
149  *dtime_temp/dzeta_c
150  at4_1(kc) = dtime_temp/dzeta_c
151  at4_2(kc) = zeta_c(kc) &
152  *dtime_temp/dzeta_c
153  at5(kc) = 1.0_dp/rho &
154  *dtime_temp/dzeta_c
155  if (kc /= kcmax) then
156  at6(kc) = 1.0_dp &
157  /dzeta_c
158  else
159  at6(kc) = 0.0_dp
160  end if
161  ai1(kc) = agediff &
162  *dtime_temp/dzeta_c
163  if (kc /= kcmax) then
164  ai2(kc) = 1.0_dp &
165  /dzeta_c
166  else
167  ai2(kc) = 0.0_dp
168  end if
169 
170  end if
171 
172 end do
173 
174 !-------- Computation loop for temperature, water content and age --------
175 
176 do i=1, imax-1 ! skipping domain margins
177 do j=1, jmax-1 ! skipping domain margins
178 
179  if (maske(j,i)==0) then ! glaciated land
180 
181 ! ------ Old vertical column cold
182 
183  if (n_cts(j,i).eq.-1) then
184 
185  n_cts_neu(j,i) = n_cts(j,i)
186  zm_neu(j,i) = zb(j,i)
187  h_c_neu(j,i) = h_c(j,i)
188  h_t_neu(j,i) = h_t(j,i)
189 
190  call calc_temp1(at1, at2_1, at2_2, at3_1, at3_2, &
191  at4_1, at4_2, at5, at6, at7, atr1, acb1, acb2, &
192  acb3, acb4, alb1, ai1, ai2, &
193  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
194 
195 ! ---- Check whether base has become temperate
196 
197  if (temp_c_neu(0,j,i).gt.temp_c_m(0,j,i)) then
198 
199  n_cts_neu(j,i) = 0
200 
201  call calc_temp2(at1, at2_1, at2_2, at3_1, at3_2, &
202  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
203  ai1, ai2, &
204  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
205 
206  end if
207 
208 ! ---- Check whether even temperate layer has formed
209 
210  if ( &
211  ( n_cts_neu(j,i).eq.0 ).and. &
212  ( (temp_c_neu(1,j,i)-temp_c_neu(0,j,i)) &
213  .gt.(am1*h_c_neu(j,i)) ) &
214  ) then
215 
216  n_cts_neu(j,i) = 1
217  zm_neu(j,i) = zb(j,i)+0.001_dp
218  h_c_neu(j,i) = h_c(j,i)-0.001_dp
219  h_t_neu(j,i) = h_t(j,i)+0.001_dp
220 ! ! CTS preliminarily positioned 1 mm above ice base --------
221 
222  call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
223  at4_1, at4_2, at5, at6, at7, atr1, &
224  am1, am2, alb1, &
225  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
226  ai1, ai2, ai3, dzeta_t, &
227  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
228 
229  call shift_cts_upward(at1, at2_1, at2_2, at3_1, at3_2, &
230  at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
231  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
232  ai1, ai2, ai3, dzeta_t, &
233  dtime_temp, dtt_2dxi, dtt_2deta, dtime_temp_inv, &
234  i, j)
235 
236  end if
237 
238 ! ------ Old vertical column with temperate base
239 
240  else if (n_cts(j,i).eq.0) then
241 
242  n_cts_neu(j,i) = n_cts(j,i)
243  zm_neu(j,i) = zb(j,i)
244  h_c_neu(j,i) = h_c(j,i)
245  h_t_neu(j,i) = h_t(j,i)
246 
247  call calc_temp2(at1, at2_1, at2_2, at3_1, at3_2, &
248  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
249  ai1, ai2, &
250  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
251 
252 ! ---- Check whether temperate base becomes cold
253 
254  if ( (temp_c_neu(1,j,i)-temp_c_neu(0,j,i)) &
255  .lt. (am1*h_c(j,i)) ) then
256 
257  n_cts_neu(j,i) = -1
258 
259  call calc_temp1(at1, at2_1, at2_2, at3_1, at3_2, &
260  at4_1, at4_2, at5, at6, at7, atr1, acb1, acb2, &
261  acb3, acb4, alb1, ai1, ai2, &
262  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
263 
264  if (temp_c_neu(0,j,i).ge.temp_c_m(0,j,i)) then
265 
266  n_cts_neu(j,i) = 0
267 
268  call calc_temp2(at1, at2_1, at2_2, at3_1, at3_2, &
269  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
270  ai1, ai2, &
271  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
272 
273  end if
274 
275  end if
276 
277 ! ---- Check whether temperate layer has formed
278 
279  if ( &
280  ( n_cts_neu(j,i).eq.0 ).and. &
281  ( (temp_c_neu(1,j,i)-temp_c_neu(0,j,i)) &
282  .gt.(am1*h_c_neu(j,i)) ) &
283  ) then
284 
285  n_cts_neu(j,i) = 1
286  zm_neu(j,i) = zb(j,i)+0.001_dp
287  h_c_neu(j,i) = h_c(j,i)-0.001_dp
288  h_t_neu(j,i) = h_t(j,i)+0.001_dp
289 ! ! CTS preliminarily positioned 1 mm above ice base --------
290 
291  call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
292  at4_1, at4_2, at5, at6, at7, atr1, &
293  am1, am2, alb1, &
294  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
295  ai1, ai2, ai3, dzeta_t, &
296  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
297 
298  call shift_cts_upward(at1, at2_1, at2_2, at3_1, at3_2, &
299  at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
300  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
301  ai1, ai2, ai3, dzeta_t, &
302  dtime_temp, dtt_2dxi, dtt_2deta, dtime_temp_inv, &
303  i, j)
304 
305  end if
306 
307 ! ------ Old vertical column with temperate base and temperate layer
308 
309  else ! n_cts(j,i).eq.1
310 
311  n_cts_neu(j,i) = n_cts(j,i)
312  zm_neu(j,i) = zm(j,i)
313  h_c_neu(j,i) = h_c(j,i)
314  h_t_neu(j,i) = h_t(j,i)
315 
316  call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
317  at4_1, at4_2, at5, at6, at7, atr1, &
318  am1, am2, alb1, &
319  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
320  ai1, ai2, ai3, dzeta_t, &
321  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
322 
323  if ( (temp_c_neu(0,j,i)-(-beta*h_c_neu(j,i))).gt.0.0_dp ) &
324  then
325  call shift_cts_upward(at1, at2_1, at2_2, at3_1, at3_2, &
326  at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
327  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
328  ai1, ai2, ai3, dzeta_t, &
329  dtime_temp, dtt_2dxi, dtt_2deta, dtime_temp_inv, &
330  i, j)
331  else
332  call shift_cts_downward(at1, at2_1, at2_2, at3_1, at3_2, &
333  at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
334  aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
335  ai1, ai2, ai3, dzeta_t, &
336  dtime_temp, dtt_2dxi, dtt_2deta, dtime_temp_inv, &
337  i, j)
338  end if
339 
340  end if
341 
342 #if (MARGIN==3)
343 
344  else if (maske(j,i)==3) then ! floating ice
345 
346  n_cts_neu(j,i) = -1
347  zm_neu(j,i) = zb(j,i)
348  h_c_neu(j,i) = h_c(j,i) + h_t(j,i)
349  h_t_neu(j,i) = 0.0_dp
350 
351  call calc_temp_ssa(at1, at2_1, at2_2, at3_1, at3_2, &
352  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
353  ai1, ai2, &
354  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
355 
356 ! ------ Reset temperatures above melting to the melting point
357 ! (should not occur, but just in case)
358 
359  do kc=0, kcmax
360  if (temp_c_neu(kc,j,i) > temp_c_m(kc,j,i)) &
361  temp_c_neu(kc,j,i) = temp_c_m(kc,j,i)
362  end do
363 
364 #endif
365 
366  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
367 
368  n_cts_neu(j,i) = -1
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  call calc_temp_r(atr1, alb1, i, j)
374 
375 endif
376 
377 end do
378 end do ! End of computation loop
379 
380 !-------- Extrapolate values on margins --------
381 
382 ! ------ Lower left corner
383 
384 i=0
385 j=0
386 
387 if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
388  ! glaciated land or floating ice
389  ii=i+1
390  jj=j+1
391 
392  do kc=0,kcmax
393  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
394  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
395  end do
396 
397  do kt=0,ktmax
398  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii) ! set temp.-ice water content
399  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii) ! set temp.-ice age
400  end do
401 
402  do kr=0,krmax
403  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
404  end do
405 
406  n_cts_neu(j,i) = min(n_cts_neu(jj,ii),0) ! temperate layer excluded
407  h_c_neu(j,i) = h_c(j,i)
408  h_t_neu(j,i) = h_t(j,i)
409  zm_neu(j,i) = zb(j,i)
410 
411 else ! maske(j,i).eq.1,2 (ice-free land or sea point)
412 
413  n_cts_neu(j,i) = -1
414  zm_neu(j,i) = zb(j,i)
415  h_c_neu(j,i) = h_c(j,i)
416  h_t_neu(j,i) = h_t(j,i)
417 
418  call calc_temp_r(atr1, alb1, i, j)
419 
420 end if
421 
422 ! ------ Lower right corner
423 
424 i=imax
425 j=0
426 
427 if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
428  ! glaciated land or floating ice
429  ii=i-1
430  jj=j+1
431 
432  do kc=0,kcmax
433  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
434  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
435  end do
436 
437  do kt=0,ktmax
438  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii) ! set temp.-ice water content
439  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii) ! set temp.-ice age
440  end do
441 
442  do kr=0,krmax
443  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
444  end do
445 
446  n_cts_neu(j,i) = min(n_cts_neu(jj,ii),0) ! temperate layer excluded
447  h_c_neu(j,i) = h_c(j,i)
448  h_t_neu(j,i) = h_t(j,i)
449  zm_neu(j,i) = zb(j,i)
450 
451 else ! maske(j,i).eq.1,2 (ice-free land or sea point)
452 
453  n_cts_neu(j,i) = -1
454  zm_neu(j,i) = zb(j,i)
455  h_c_neu(j,i) = h_c(j,i)
456  h_t_neu(j,i) = h_t(j,i)
457 
458  call calc_temp_r(atr1, alb1, i, j)
459 
460 end if
461 
462 ! ------ Upper left corner
463 
464 i=0
465 j=jmax
466 
467 if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
468  ! glaciated land or floating ice
469  ii=i+1
470  jj=j-1
471 
472  do kc=0,kcmax
473  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
474  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
475  end do
476 
477  do kt=0,ktmax
478  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii) ! set temp.-ice water content
479  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii) ! set temp.-ice age
480  end do
481 
482  do kr=0,krmax
483  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
484  end do
485 
486  n_cts_neu(j,i) = min(n_cts_neu(jj,ii),0) ! temperate layer excluded
487  h_c_neu(j,i) = h_c(j,i)
488  h_t_neu(j,i) = h_t(j,i)
489  zm_neu(j,i) = zb(j,i)
490 
491 else ! maske(j,i).eq.1,2 (ice-free land or sea point)
492 
493  n_cts_neu(j,i) = -1
494  zm_neu(j,i) = zb(j,i)
495  h_c_neu(j,i) = h_c(j,i)
496  h_t_neu(j,i) = h_t(j,i)
497 
498  call calc_temp_r(atr1, alb1, i, j)
499 
500 end if
501 
502 ! ------ Upper right corner
503 
504 i=imax
505 j=jmax
506 
507 if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
508  ! glaciated land or floating ice
509  ii=i-1
510  jj=j-1
511 
512  do kc=0,kcmax
513  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
514  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
515  end do
516 
517  do kt=0,ktmax
518  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii) ! set temp.-ice water content
519  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii) ! set temp.-ice age
520  end do
521 
522  do kr=0,krmax
523  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
524  end do
525 
526  n_cts_neu(j,i) = min(n_cts_neu(jj,ii),0) ! temperate layer excluded
527  h_c_neu(j,i) = h_c(j,i)
528  h_t_neu(j,i) = h_t(j,i)
529  zm_neu(j,i) = zb(j,i)
530 
531 else ! maske(j,i).eq.1,2 (ice-free land or sea point)
532 
533  n_cts_neu(j,i) = -1
534  zm_neu(j,i) = zb(j,i)
535  h_c_neu(j,i) = h_c(j,i)
536  h_t_neu(j,i) = h_t(j,i)
537 
538  call calc_temp_r(atr1, alb1, i, j)
539 
540 end if
541 
542 ! ------ Lower and upper margins
543 
544 do i=1, imax-1
545 
546 ! ---- Lower margin
547 
548  j=0
549 
550  if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
551  ! glaciated land or floating ice
552  ii=i
553  jj=j+1
554 
555  do kc=0,kcmax
556  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
557  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
558  end do
559 
560  do kt=0,ktmax
561  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii) ! set temp.-ice water content
562  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii) ! set temp.-ice age
563  end do
564 
565  do kr=0,krmax
566  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
567  end do
568 
569  n_cts_neu(j,i) = min(n_cts_neu(jj,ii),0) ! temperate layer excluded
570  h_c_neu(j,i) = h_c(j,i)
571  h_t_neu(j,i) = h_t(j,i)
572  zm_neu(j,i) = zb(j,i)
573 
574  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
575 
576  n_cts_neu(j,i) = -1
577  zm_neu(j,i) = zb(j,i)
578  h_c_neu(j,i) = h_c(j,i)
579  h_t_neu(j,i) = h_t(j,i)
580 
581  call calc_temp_r(atr1, alb1, i, j)
582 
583  end if
584 
585 ! ---- Upper margin
586 
587  j=jmax
588 
589  if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
590  ! glaciated land or floating ice
591  ii=i
592  jj=j-1
593 
594  do kc=0,kcmax
595  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
596  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
597  end do
598 
599  do kt=0,ktmax
600  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii) ! set temp.-ice water content
601  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii) ! set temp.-ice age
602  end do
603 
604  do kr=0,krmax
605  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
606  end do
607 
608  n_cts_neu(j,i) = min(n_cts_neu(jj,ii),0) ! temperate layer excluded
609  h_c_neu(j,i) = h_c(j,i)
610  h_t_neu(j,i) = h_t(j,i)
611  zm_neu(j,i) = zb(j,i)
612 
613  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
614 
615  n_cts_neu(j,i) = -1
616  zm_neu(j,i) = zb(j,i)
617  h_c_neu(j,i) = h_c(j,i)
618  h_t_neu(j,i) = h_t(j,i)
619 
620  call calc_temp_r(atr1, alb1, i, j)
621 
622  end if
623 
624 end do
625 
626 ! ------ Left and right margins
627 
628 do j=1, jmax-1
629 
630 ! ---- Left margin
631 
632  i=0
633 
634  if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
635  ! glaciated land or floating ice
636  ii=i+1
637  jj=j
638 
639  do kc=0,kcmax
640  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
641  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
642  end do
643 
644  do kt=0,ktmax
645  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii) ! set temp.-ice water content
646  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii) ! set temp.-ice age
647  end do
648 
649  do kr=0,krmax
650  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
651  end do
652 
653  n_cts_neu(j,i) = min(n_cts_neu(jj,ii),0) ! temperate layer excluded
654  h_c_neu(j,i) = h_c(j,i)
655  h_t_neu(j,i) = h_t(j,i)
656  zm_neu(j,i) = zb(j,i)
657 
658  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
659 
660  n_cts_neu(j,i) = -1
661  zm_neu(j,i) = zb(j,i)
662  h_c_neu(j,i) = h_c(j,i)
663  h_t_neu(j,i) = h_t(j,i)
664 
665  call calc_temp_r(atr1, alb1, i, j)
666 
667  end if
668 
669 ! ---- Right margin
670 
671  i=imax
672 
673  if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
674  ! glaciated land or floating ice
675  ii=i-1
676  jj=j
677 
678  do kc=0,kcmax
679  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
680  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
681  end do
682 
683  do kt=0,ktmax
684  omega_t_neu(kt,j,i) = omega_t_neu(kt,jj,ii) ! set temp.-ice water content
685  age_t_neu(kt,j,i) = age_t_neu(kt,jj,ii) ! set temp.-ice age
686  end do
687 
688  do kr=0,krmax
689  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
690  end do
691 
692  n_cts_neu(j,i) = min(n_cts_neu(jj,ii),0) ! temperate layer excluded
693  h_c_neu(j,i) = h_c(j,i)
694  h_t_neu(j,i) = h_t(j,i)
695  zm_neu(j,i) = zb(j,i)
696 
697  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
698 
699  n_cts_neu(j,i) = -1
700  zm_neu(j,i) = zb(j,i)
701  h_c_neu(j,i) = h_c(j,i)
702  h_t_neu(j,i) = h_t(j,i)
703 
704  call calc_temp_r(atr1, alb1, i, j)
705 
706  end if
707 
708 end do
709 
710 !-------- Dummy values for omega_c_neu and kc_cts_neu --------
711 
712 omega_c_neu = 0.0_dp ! not needed for
713 kc_cts_neu = 0 ! the polythermal mode
714 
715 !-------- Smoothing of H_t_neu with numerical diffusion --------
716 
717 ! ------ Volume of temperate ice without smoothing
718 
719 vol_t = 0.0_dp
720 do i=0, imax ! extended to domain margins (22.1.02 -> V1.1)
721 do j=0, jmax ! extended to domain margins (22.1.02 -> V1.1)
722  if (n_cts_neu(j,i).eq.1) then
723  vol_t = vol_t + h_t_neu(j,i)*area(j,i)
724  end if
725 end do
726 end do
727 
728 ! ------ Smoothing
729 
730 do i=1, imax-1
731 do j=1, jmax-1
732  if (n_cts_neu(j,i).ne.-1) then
733 
734  dh_t_smooth(j,i) = numdiff_h_t* ( -4.0_dp*h_t_neu(j,i) &
735  +h_t_neu(j,i+1)+h_t_neu(j,i-1) &
736  +h_t_neu(j+1,i)+h_t_neu(j-1,i) )
737  if (dh_t_smooth(j,i).gt.0.001_dp) n_cts_neu(j,i) = 1
738 
739  end if
740 end do
741 end do
742 
743 do i=1, imax-1
744 do j=1, jmax-1
745  if (n_cts_neu(j,i).eq.1) then
746  h_t_neu(j,i) = h_t_neu(j,i) + dh_t_smooth(j,i)
747  end if
748 end do
749 end do
750 
751 ! ------ Volume of temperate ice with smoothing
752 
753 vol_t_smooth = 0.0_dp
754 do i=0, imax ! extended to domain margins (22.1.02 -> V1.1)
755 do j=0, jmax ! extended to domain margins (22.1.02 -> V1.1)
756  if (n_cts_neu(j,i).eq.1) then
757  vol_t_smooth = vol_t_smooth + h_t_neu(j,i)*area(j,i)
758  end if
759 end do
760 end do
761 
762 ! ------ Correction so that volume is not changed by the smoothing
763 
764 if (vol_t_smooth.gt.0.0_dp) then
765 
766  korrfakt_t = vol_t/vol_t_smooth
767  do i=0, imax ! extended to domain margins (22.1.02 -> V1.1)
768  do j=0, jmax ! extended to domain margins (22.1.02 -> V1.1)
769  if (n_cts_neu(j,i).eq.1) then
770  h_t_neu(j,i) = h_t_neu(j,i)*korrfakt_t
771 ! zm_neu(j,i) = zb(j,i) + H_t_neu(j,i)
772 ! H_c_neu(j,i) = zs(j,i) - zm_neu(j,i)
773  end if
774  end do
775  end do
776 
777 end if
778 
779 !-------- Numerical time lag for evolution of H_t_neu --------
780 
781 time_lag_cts = tau_cts*year_sec ! yr --> s
782 
783 do i=0, imax ! extended to domain margins (22.1.02 -> V1.1)
784 do j=0, jmax ! extended to domain margins (22.1.02 -> V1.1)
785 
786  if (n_cts_neu(j,i).eq.1) then
787 
788  h_t_neu(j,i) = ( time_lag_cts*h_t(j,i) &
789  + dtime_temp*h_t_neu(j,i) ) &
790  /(time_lag_cts+dtime_temp)
791 
792  zm_neu(j,i) = zb(j,i) + h_t_neu(j,i)
793  h_c_neu(j,i) = zs(j,i) - zm_neu(j,i)
794 
795  end if
796 
797 end do
798 end do
799 
800 end subroutine calc_temp_poly
801 !
calc_temp3
subroutine calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, ai1, ai2, ai3, dzeta_t, dtime_temp, dtt_2dxi, dtt_2deta, i, j)
Computation of temperature, water content and age for an ice column with a temperate base overlain by...
Definition: calc_temp3.F90:37
sico_types
Declarations of kind types for SICOPOLIS.
Definition: sico_types.F90:35
shift_cts_downward
subroutine shift_cts_downward(at1, at2_1, at2_2, at3_1, at3_2, at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, ai1, ai2, ai3, dzeta_t, dtime_temp, dtt_2dxi, dtt_2deta, dtime_temp_inv, i, j)
Downward shifting of the CTS.
Definition: shift_cts_downward.F90:35
sico_vars
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars.F90:35
calc_temp_poly
subroutine calc_temp_poly(dxi, deta, dzeta_c, dzeta_t, dzeta_r, dtime_temp)
Computation of temperature, water content and age in polythermal mode.
Definition: calc_temp_poly.F90:35
ice_material_quantities::kappa_val
real(dp) function, public kappa_val(temp_val)
Heat conductivity of ice: Linear interpolation of tabulated values in KAPPA(.).
Definition: ice_material_quantities.F90:217
calc_temp1
subroutine calc_temp1(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)
Computation of temperature and age for a cold ice column.
Definition: calc_temp1.F90:35
ice_material_quantities
Material quantities of ice: Rate factor, heat conductivity, specific heat (heat capacity).
Definition: ice_material_quantities.F90:37
calc_temp_ssa
subroutine calc_temp_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).
Definition: calc_temp_ssa.F90:35
shift_cts_upward
subroutine shift_cts_upward(at1, at2_1, at2_2, at3_1, at3_2, at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, ai1, ai2, ai3, dzeta_t, dtime_temp, dtt_2dxi, dtt_2deta, dtime_temp_inv, i, j)
Upward shifting of the CTS.
Definition: shift_cts_upward.F90:35
sico_variables
Declarations of global variables for SICOPOLIS.
Definition: sico_variables.F90:35
calc_temp_r
subroutine calc_temp_r(atr1, alb1, i, j)
Computation of temperature for an ice-free column.
Definition: calc_temp_r.F90:35
calc_temp2
subroutine calc_temp2(at1, at2_1, at2_2, at3_1, at3_2, at4_1, at4_2, at5, at6, at7, atr1, alb1, ai1, ai2, dtime_temp, dtt_2dxi, dtt_2deta, i, j)
Computation of temperature and age for an ice column with a temperate base overlain by cold ice...
Definition: calc_temp2.F90:37
Generated on Wed Mar 2 2016 11:26:55 for SICOPOLIS V3.2 by   doxygen 1.8.5