SICOPOLIS V3.2
 All Classes Files Functions Variables Macros
calc_temp_cold.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : c a l c _ t e m p _ c o l d
4 !
5 !> @file
6 !!
7 !! Computation of temperature and age in cold-ice 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 and age in cold-ice mode.
34 !<------------------------------------------------------------------------------
35 subroutine calc_temp_cold(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, 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 
59 !-------- Term abbreviations
60 
61 at7 = 2.0_dp/rho*dtime_temp
62 
63 aw1 = dtime_temp/dzeta_t
64 aw2 = dtime_temp/dzeta_t
65 aw3 = dtime_temp/dzeta_t
66 aw4 = dtime_temp/dzeta_t
67 aw5 = nue/rho*dtime_temp/(dzeta_t**2)
68 aw7 = 2.0_dp/(rho*l)*dtime_temp
69 aw8 = beta**2/(rho*l) &
70  *(kappa_val(0.0_dp)-kappa_val(-1.0_dp))*dtime_temp
71 aw9 = beta/l*dtime_temp
72 
73 ai3 = agediff*dtime_temp/(dzeta_t**2)
74 
75 atr1 = kappa_r/(rho_c_r*h_r**2)*dtime_temp/(dzeta_r**2)
76 
77 if (flag_aa_nonzero) then
78  am1 = aa*beta*dzeta_c/(ea-1.0_dp)
79  am2 = aa*l*rho*dzeta_c/(ea-1.0_dp)
80 else
81  am1 = beta*dzeta_c
82  am2 = l*rho*dzeta_c
83 end if
84 
85 if (flag_aa_nonzero) then
86  acb1 = (ea-1.0_dp)/aa/dzeta_c
87 else
88  acb1 = 1.0_dp/dzeta_c
89 end if
90 
91 acb2 = kappa_r/h_r/dzeta_r
92 acb3 = rho*g
93 acb4 = rho*g
94 
95 alb1 = h_r/kappa_r*dzeta_r
96 
97 aqtld = dzeta_t/dtime_temp
98 
99 dtt_2dxi = 0.5_dp*dtime_temp/dxi
100 dtt_2deta = 0.5_dp*dtime_temp/deta
101 
102 dtime_temp_inv = 1.0_dp/dtime_temp
103 
104 do kc=0, kcmax
105 
106  if (flag_aa_nonzero) then
107 
108  at1(kc) = (ea-1.0_dp)/(aa*eaz_c(kc))*dtime_temp/dzeta_c
109  at2_1(kc) = (ea-1.0_dp)/(aa*eaz_c(kc))*dtime_temp/dzeta_c
110  at2_2(kc) = (eaz_c(kc)-1.0_dp)/(aa*eaz_c(kc)) &
111  *dtime_temp/dzeta_c
112  at3_1(kc) = (ea-1.0_dp)/(aa*eaz_c(kc))*dtime_temp/dzeta_c
113  at3_2(kc) = (eaz_c(kc)-1.0_dp)/(aa*eaz_c(kc)) &
114  *dtime_temp/dzeta_c
115  at4_1(kc) = (ea-1.0_dp)/(aa*eaz_c(kc))*dtime_temp/dzeta_c
116  at4_2(kc) = (eaz_c(kc)-1.0_dp)/(aa*eaz_c(kc)) &
117  *dtime_temp/dzeta_c
118  at5(kc) = (ea-1.0_dp)/(rho*aa*eaz_c(kc)) &
119  *dtime_temp/dzeta_c
120  if (kc /= kcmax) then
121  at6(kc) = (ea-1.0_dp) &
122  /(aa*exp(aa*0.5_dp*(zeta_c(kc)+zeta_c(kc+1)))) &
123  /dzeta_c
124  else
125  at6(kc) = 0.0_dp
126  end if
127  ai1(kc) = agediff*(ea-1.0_dp)/(aa*eaz_c(kc)) &
128  *dtime_temp/dzeta_c
129  if (kc /= kcmax) then
130  ai2(kc) = (ea-1.0_dp) &
131  /(aa*exp(aa*0.5_dp*(zeta_c(kc)+zeta_c(kc+1)))) &
132  /dzeta_c
133  else
134  ai2(kc) = 0.0_dp
135  end if
136 
137  else
138 
139  at1(kc) = dtime_temp/dzeta_c
140  at2_1(kc) = dtime_temp/dzeta_c
141  at2_2(kc) = zeta_c(kc) &
142  *dtime_temp/dzeta_c
143  at3_1(kc) = dtime_temp/dzeta_c
144  at3_2(kc) = zeta_c(kc) &
145  *dtime_temp/dzeta_c
146  at4_1(kc) = dtime_temp/dzeta_c
147  at4_2(kc) = zeta_c(kc) &
148  *dtime_temp/dzeta_c
149  at5(kc) = 1.0_dp/rho &
150  *dtime_temp/dzeta_c
151  if (kc /= kcmax) then
152  at6(kc) = 1.0_dp &
153  /dzeta_c
154  else
155  at6(kc) = 0.0_dp
156  end if
157  ai1(kc) = agediff &
158  *dtime_temp/dzeta_c
159  if (kc /= kcmax) then
160  ai2(kc) = 1.0_dp &
161  /dzeta_c
162  else
163  ai2(kc) = 0.0_dp
164  end if
165 
166  end if
167 
168 end do
169 
170 !-------- Computation loop for temperature and age --------
171 
172 do i=1, imax-1 ! skipping domain margins
173 do j=1, jmax-1 ! skipping domain margins
174 
175  if (maske(j,i)==0) then ! glaciated land
176 
177  n_cts_neu(j,i) = -1
178  zm_neu(j,i) = zb(j,i)
179  h_c_neu(j,i) = h_c(j,i)
180  h_t_neu(j,i) = h_t(j,i)
181 
182  call calc_temp1(at1, at2_1, at2_2, at3_1, at3_2, &
183  at4_1, at4_2, at5, at6, at7, atr1, acb1, acb2, &
184  acb3, acb4, alb1, ai1, ai2, &
185  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
186 
187 ! ------ Reset temperatures above melting to the melting point,
188 ! look for the CTS
189 
190  kc_cts_neu(j,i) = 0
191 
192  if (temp_c_neu(0,j,i).gt.temp_c_m(0,j,i)) then
193  n_cts_neu(j,i) = 0
194  kc_cts_neu(j,i) = 0
195  temp_c_neu(0,j,i) = temp_c_m(0,j,i)
196  temp_r_neu(krmax,j,i) = temp_c_m(0,j,i)
197  end if
198 
199  do kc=1, kcmax
200  if (temp_c_neu(kc,j,i).gt.temp_c_m(kc,j,i)) then
201  kc_cts_neu(j,i) = kc
202  temp_c_neu(kc,j,i) = temp_c_m(kc,j,i)
203  end if
204  end do
205 
206 #if (MARGIN==3)
207 
208  else if (maske(j,i)==3) then ! floating ice
209 
210  n_cts_neu(j,i) = -1
211  kc_cts_neu(j,i) = 0
212  zm_neu(j,i) = zb(j,i)
213  h_c_neu(j,i) = h_c(j,i)
214  h_t_neu(j,i) = 0.0_dp
215 
216  call calc_temp_ssa(at1, at2_1, at2_2, at3_1, at3_2, &
217  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
218  ai1, ai2, &
219  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
220 
221 ! ------ Reset temperatures above melting to the melting point
222 ! (should not occur, but just in case)
223 
224  do kc=0, kcmax
225  if (temp_c_neu(kc,j,i) > temp_c_m(kc,j,i)) &
226  temp_c_neu(kc,j,i) = temp_c_m(kc,j,i)
227  end do
228 
229 #endif
230 
231  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
232 
233  n_cts_neu(j,i) = -1
234  kc_cts_neu(j,i) = 0
235  zm_neu(j,i) = zb(j,i)
236  h_c_neu(j,i) = h_c(j,i)
237  h_t_neu(j,i) = h_t(j,i)
238 
239  call calc_temp_r(atr1, alb1, i, j)
240 
241  end if
242 
243 end do
244 end do
245 
246 !-------- Extrapolate values on margins --------
247 
248 ! ------ Lower left corner
249 
250 i=0
251 j=0
252 
253 if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
254  ! glaciated land or floating ice
255  ii=i+1
256  jj=j+1
257 
258  do kc=0,kcmax
259  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
260  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
261  end do
262 
263  do kr=0,krmax
264  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
265  end do
266 
267  n_cts_neu(j,i) = n_cts_neu(jj,ii)
268  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
269  zm_neu(j,i) = zb(j,i)
270  h_c_neu(j,i) = h_c(j,i)
271  h_t_neu(j,i) = h_t(j,i)
272 
273 else ! maske(j,i).eq.1,2 (ice-free land or sea point)
274 
275  n_cts_neu(j,i) = -1
276  kc_cts_neu(j,i) = 0
277  zm_neu(j,i) = zb(j,i)
278  h_c_neu(j,i) = h_c(j,i)
279  h_t_neu(j,i) = h_t(j,i)
280 
281  call calc_temp_r(atr1, alb1, i, j)
282 
283 end if
284 
285 ! ------ Lower right corner
286 
287 i=imax
288 j=0
289 
290 if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
291  ! glaciated land or floating ice
292  ii=i-1
293  jj=j+1
294 
295  do kc=0,kcmax
296  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
297  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
298  end do
299 
300  do kr=0,krmax
301  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
302  end do
303 
304  n_cts_neu(j,i) = n_cts_neu(jj,ii)
305  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
306  zm_neu(j,i) = zb(j,i)
307  h_c_neu(j,i) = h_c(j,i)
308  h_t_neu(j,i) = h_t(j,i)
309 
310 else ! maske(j,i).eq.1,2 (ice-free land or sea point)
311 
312  n_cts_neu(j,i) = -1
313  kc_cts_neu(j,i) = 0
314  zm_neu(j,i) = zb(j,i)
315  h_c_neu(j,i) = h_c(j,i)
316  h_t_neu(j,i) = h_t(j,i)
317 
318  call calc_temp_r(atr1, alb1, i, j)
319 
320 end if
321 
322 ! ------ Upper left corner
323 
324 i=0
325 j=jmax
326 
327 if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
328  ! glaciated land or floating ice
329  ii=i+1
330  jj=j-1
331 
332  do kc=0,kcmax
333  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
334  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
335  end do
336 
337  do kr=0,krmax
338  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
339  end do
340 
341  n_cts_neu(j,i) = n_cts_neu(jj,ii)
342  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
343  zm_neu(j,i) = zb(j,i)
344  h_c_neu(j,i) = h_c(j,i)
345  h_t_neu(j,i) = h_t(j,i)
346 
347 else ! maske(j,i).eq.1,2 (ice-free land or sea point)
348 
349  n_cts_neu(j,i) = -1
350  kc_cts_neu(j,i) = 0
351  zm_neu(j,i) = zb(j,i)
352  h_c_neu(j,i) = h_c(j,i)
353  h_t_neu(j,i) = h_t(j,i)
354 
355  call calc_temp_r(atr1, alb1, i, j)
356 
357 end if
358 
359 ! ------ Upper right corner
360 
361 i=imax
362 j=jmax
363 
364 if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
365  ! glaciated land or floating ice
366  ii=i-1
367  jj=j-1
368 
369  do kc=0,kcmax
370  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
371  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
372  end do
373 
374  do kr=0,krmax
375  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
376  end do
377 
378  n_cts_neu(j,i) = n_cts_neu(jj,ii)
379  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
380  zm_neu(j,i) = zb(j,i)
381  h_c_neu(j,i) = h_c(j,i)
382  h_t_neu(j,i) = h_t(j,i)
383 
384 else ! maske(j,i).eq.1,2 (ice-free land or sea point)
385 
386  n_cts_neu(j,i) = -1
387  kc_cts_neu(j,i) = 0
388  zm_neu(j,i) = zb(j,i)
389  h_c_neu(j,i) = h_c(j,i)
390  h_t_neu(j,i) = h_t(j,i)
391 
392  call calc_temp_r(atr1, alb1, i, j)
393 
394 end if
395 
396 ! ------ Lower and upper margins
397 
398 do i=1, imax-1
399 
400 ! ---- Lower margin
401 
402  j=0
403 
404  if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
405  ! glaciated land or floating ice
406  ii=i
407  jj=j+1
408 
409  do kc=0,kcmax
410  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
411  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
412  end do
413 
414  do kr=0,krmax
415  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
416  end do
417 
418  n_cts_neu(j,i) = n_cts_neu(jj,ii)
419  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
420  zm_neu(j,i) = zb(j,i)
421  h_c_neu(j,i) = h_c(j,i)
422  h_t_neu(j,i) = h_t(j,i)
423 
424  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
425 
426  n_cts_neu(j,i) = -1
427  kc_cts_neu(j,i) = 0
428  zm_neu(j,i) = zb(j,i)
429  h_c_neu(j,i) = h_c(j,i)
430  h_t_neu(j,i) = h_t(j,i)
431 
432  call calc_temp_r(atr1, alb1, i, j)
433 
434  end if
435 
436 ! ---- Upper margin
437 
438  j=jmax
439 
440  if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
441  ! glaciated land or floating ice
442  ii=i
443  jj=j-1
444 
445  do kc=0,kcmax
446  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
447  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
448  end do
449 
450  do kr=0,krmax
451  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
452  end do
453 
454  n_cts_neu(j,i) = n_cts_neu(jj,ii)
455  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
456  zm_neu(j,i) = zb(j,i)
457  h_c_neu(j,i) = h_c(j,i)
458  h_t_neu(j,i) = h_t(j,i)
459 
460  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
461 
462  n_cts_neu(j,i) = -1
463  kc_cts_neu(j,i) = 0
464  zm_neu(j,i) = zb(j,i)
465  h_c_neu(j,i) = h_c(j,i)
466  h_t_neu(j,i) = h_t(j,i)
467 
468  call calc_temp_r(atr1, alb1, i, j)
469 
470  end if
471 
472 end do
473 
474 ! ------ Left and right margins
475 
476 do j=1, jmax-1
477 
478 ! ---- Left margin
479 
480  i=0
481 
482  if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
483  ! glaciated land or floating ice
484  ii=i+1
485  jj=j
486 
487  do kc=0,kcmax
488  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
489  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
490  end do
491 
492  do kr=0,krmax
493  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
494  end do
495 
496  n_cts_neu(j,i) = n_cts_neu(jj,ii)
497  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
498  zm_neu(j,i) = zb(j,i)
499  h_c_neu(j,i) = h_c(j,i)
500  h_t_neu(j,i) = h_t(j,i)
501 
502  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
503 
504  n_cts_neu(j,i) = -1
505  kc_cts_neu(j,i) = 0
506  zm_neu(j,i) = zb(j,i)
507  h_c_neu(j,i) = h_c(j,i)
508  h_t_neu(j,i) = h_t(j,i)
509 
510  call calc_temp_r(atr1, alb1, i, j)
511 
512  end if
513 
514 ! ---- Right margin
515 
516  i=imax
517 
518  if ( (maske(j,i).eq.0).or.(maske(j,i).eq.3) ) then
519  ! glaciated land or floating ice
520  ii=i-1
521  jj=j
522 
523  do kc=0,kcmax
524  temp_c_neu(kc,j,i) = temp_c_neu(kc,jj,ii) ! set cold-ice temperature
525  age_c_neu(kc,j,i) = age_c_neu(kc,jj,ii) ! set cold-ice age
526  end do
527 
528  do kr=0,krmax
529  temp_r_neu(kr,j,i) = temp_r_neu(kr,jj,ii) ! set bedrock temperature
530  end do
531 
532  n_cts_neu(j,i) = n_cts_neu(jj,ii)
533  kc_cts_neu(j,i) = kc_cts_neu(jj,ii)
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  else ! maske(j,i).eq.1,2 (ice-free land or sea point)
539 
540  n_cts_neu(j,i) = -1
541  kc_cts_neu(j,i) = 0
542  zm_neu(j,i) = zb(j,i)
543  h_c_neu(j,i) = h_c(j,i)
544  h_t_neu(j,i) = h_t(j,i)
545 
546  call calc_temp_r(atr1, alb1, i, j)
547 
548  end if
549 
550 end do
551 
552 !-------- Dummy values for omega_c_neu --------
553 
554 omega_c_neu = 0.0_dp ! not computed in the cold-ice mode
555 
556 end subroutine calc_temp_cold
557 !
sico_types
Declarations of kind types for SICOPOLIS.
Definition: sico_types.F90:35
sico_vars
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars.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
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_temp_cold
subroutine calc_temp_cold(dxi, deta, dzeta_c, dzeta_t, dzeta_r, dtime_temp)
Computation of temperature and age in cold-ice mode.
Definition: calc_temp_cold.F90:35
Generated on Wed Mar 2 2016 11:26:54 for SICOPOLIS V3.2 by   doxygen 1.8.5