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