SICOPOLIS V3.0
 All Classes Files Functions Variables Macros
calc_temp_ssa.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : c a l c _ t e m p _ s s a
4 !
5 !> @file
6 !!
7 !! Computation of temperature and age for ice shelves (floating ice).
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2013 Ralf Greve, Tatsuru Sato
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 for ice shelves (floating ice).
34 !<------------------------------------------------------------------------------
35 subroutine calc_temp_ssa(at1, at2_1, at2_2, at3_1, at3_2, &
36  at4_1, at4_2, at5, at6, at7, atr1, alb1, &
37  ai1, ai2, ai4, mean_accum_inv, &
38  dtime_temp, dtt_2dxi, dtt_2deta, i, j)
39 
40 use sico_types
41 use sico_variables
42 
43 implicit none
44 integer(i4b) :: i, j, kc, kt, kr
45 real(dp) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
46  at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
47  at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
48  ai1(0:kcmax), ai2(0:kcmax), ai4, &
49  atr1, alb1
50 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
51  ct5(0:kcmax), ct6(0:kcmax), ct7(0:kcmax), ctr1, clb1
52 real(dp) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), ct3_sg(0:kcmax), &
53  ct4_sg(0:kcmax), adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
54 real(dp) :: ci1(0:kcmax), ci2(0:kcmax)
55 real(dp) :: ftx_c_l(0:kcmax), ftx_c_r(0:kcmax), &
56  fty_c_l(0:kcmax), fty_c_r(0:kcmax), &
57  fax_c_l(0:kcmax), fax_c_r(0:kcmax), &
58  fay_c_l(0:kcmax), fay_c_r(0:kcmax)
59 real(dp) :: temp_c_help(0:kcmax)
60 real(dp) :: mean_accum_inv
61 real(dp) :: dtt_2dxi, dtt_2deta, dtime_temp
62 real(dp) :: dtt_dxi, dtt_deta
63 real(dp) :: viscosity, kappa_val, c_val
64 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
65  lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
66  lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
67  lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
68  lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
69 real(dp), parameter :: zero=0.0_dp
70 
71 !-------- Check for boundary points --------
72 
73 if ((i.eq.0).or.(i.eq.imax).or.(j.eq.0).or.(j.eq.jmax)) &
74  stop ' calc_temp_ssa: Boundary points not allowed.'
75 
76 !-------- Abbreviations --------
77 
78 ctr1 = atr1
79 clb1 = alb1*q_geo(j,i)
80 
81 #if ADV_VERT==1
82 
83 do kc=1, kcmax
84  ct1(kc) = at1(kc)/h_c(j,i)*0.5_dp*(vz_c(kc,j,i)+vz_c(kc-1,j,i))
85 end do
86 
87 #elif (ADV_VERT==2 || ADV_VERT==3)
88 
89 do kc=0, kcmax-1
90  ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/h_c(j,i)*vz_c(kc,j,i)
91 end do
92 
93 #endif
94 
95 do kc=0, kcmax
96 
97  ct2(kc) = ( at2_1(kc)*dzm_dtau(j,i) &
98  +at2_2(kc)*dh_c_dtau(j,i) )/h_c(j,i)
99  ct3(kc) = ( at3_1(kc)*dzm_dxi_g(j,i) &
100  +at3_2(kc)*dh_c_dxi_g(j,i) )/h_c(j,i) &
101  *0.5_dp*(vx_c(kc,j,i)+vx_c(kc,j,i-1)) *insq_g11_g(j,i)
102  ct4(kc) = ( at4_1(kc)*dzm_deta_g(j,i) &
103  +at4_2(kc)*dh_c_deta_g(j,i) )/h_c(j,i) &
104  *0.5_dp*(vy_c(kc,j,i)+vy_c(kc,j-1,i)) *insq_g22_g(j,i)
105  ct5(kc) = at5(kc) &
106  /c_val(temp_c(kc,j,i)) &
107  /h_c(j,i)
108  ct7(kc) = 2.0_dp*at7 &
109  /c_val(temp_c(kc,j,i)) &
110  *viscosity(de_ssa(j,i), &
111  temp_c(kc,j,i), temp_c_m(kc,j,i), 0.0_dp, &
112  enh_c(kc,j,i), .true.) &
113  *de_ssa(j,i)**2
114  ci1(kc) = ai1(kc)/h_c(j,i)
115 
116 ! ------ Fluxes for horizontal advection of temperature and age
117 
118 #if ADV_HOR==4
119  if (vx_c(kc,j,i-1).ge.zero) then
120  ftx_c_l(kc) = temp_c(kc,j,i-1)*vx_c(kc,j,i-1)
121  fax_c_l(kc) = age_c(kc,j,i-1)*vx_c(kc,j,i-1)
122  else
123  ftx_c_l(kc) = temp_c(kc,j,i)*vx_c(kc,j,i-1)
124  fax_c_l(kc) = age_c(kc,j,i)*vx_c(kc,j,i-1)
125  end if
126 
127  if (vx_c(kc,j,i).ge.zero) then
128  ftx_c_r(kc) = temp_c(kc,j,i)*vx_c(kc,j,i)
129  fax_c_r(kc) = age_c(kc,j,i)*vx_c(kc,j,i)
130  else
131  ftx_c_r(kc) = temp_c(kc,j,i+1)*vx_c(kc,j,i)
132  fax_c_r(kc) = age_c(kc,j,i+1)*vx_c(kc,j,i)
133  end if
134 
135  if (vy_c(kc,j-1,i).ge.zero) then
136  fty_c_l(kc) = temp_c(kc,j-1,i)*vy_c(kc,j-1,i)
137  fay_c_l(kc) = age_c(kc,j-1,i)*vy_c(kc,j-1,i)
138  else
139  fty_c_l(kc) = temp_c(kc,j,i)*vy_c(kc,j-1,i)
140  fay_c_l(kc) = age_c(kc,j,i)*vy_c(kc,j-1,i)
141  end if
142 
143  if (vy_c(kc,j,i).ge.zero) then
144  fty_c_r(kc) = temp_c(kc,j,i)*vy_c(kc,j,i)
145  fay_c_r(kc) = age_c(kc,j,i)*vy_c(kc,j,i)
146  else
147  fty_c_r(kc) = temp_c(kc,j+1,i)*vy_c(kc,j,i)
148  fay_c_r(kc) = age_c(kc,j+1,i)*vy_c(kc,j,i)
149  end if
150 #endif
151 
152 end do
153 
154 #if (ADV_VERT==2 || ADV_VERT==3)
155 
156 do kc=0, kcmax-1
157  ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
158  ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
159  ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
160  adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
161  abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
162 end do
163 
164 #endif
165 
166 do kc=0, kcmax-1
167  temp_c_help(kc) = 0.5_dp*(temp_c(kc,j,i)+temp_c(kc+1,j,i))
168  ct6(kc) = at6(kc) &
169  *kappa_val(temp_c_help(kc)) &
170  /h_c(j,i)
171  ci2(kc) = ai2(kc)/h_c(j,i)
172 end do
173 
174 #if ADV_HOR==4
175 dtt_dxi = 2.0_dp*dtt_2dxi
176 dtt_deta = 2.0_dp*dtt_2deta
177 #endif
178 
179 !-------- Set up the equations for the bedrock temperature --------
180 
181 kr=0
182 lgs_a1(kr) = 1.0_dp
183 lgs_a2(kr) = -1.0_dp
184 lgs_b(kr) = clb1
185 
186 #if Q_LITHO==1
187 ! (coupled heat-conducting bedrock)
188 
189 do kr=1, krmax-1
190  lgs_a0(kr) = - ctr1
191  lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
192  lgs_a2(kr) = - ctr1
193  lgs_b(kr) = temp_r(kr,j,i)
194 end do
195 
196 #elif Q_LITHO==0
197 ! (no coupled heat-conducting bedrock)
198 
199 do kr=1, krmax-1
200  lgs_a0(kr) = 1.0_dp
201  lgs_a1(kr) = 0.0_dp
202  lgs_a2(kr) = -1.0_dp
203  lgs_b(kr) = 2.0_dp*clb1
204 end do
205 
206 #endif
207 
208 kr=krmax
209 lgs_a0(kr) = 0.0_dp
210 lgs_a1(kr) = 1.0_dp
211 lgs_b(kr) = temp_c_m(0,j,i)-delta_tm_sw
212 
213 !-------- Solve system of linear equations --------
214 
215 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
216 
217 !-------- Assign the result --------
218 
219 do kr=0, krmax
220  temp_r_neu(kr,j,i) = lgs_x(kr)
221 end do
222 
223 !-------- Set up the equations for the ice temperature --------
224 
225 kc=0
226 lgs_a1(kc) = 1.0_dp
227 lgs_a2(kc) = 0.0_dp
228 lgs_b(kc) = temp_c_m(0,j,i)-delta_tm_sw
229 
230 do kc=1, kcmax-1
231 
232 #if ADV_VERT==1
233 
234  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
235  -ct5(kc)*ct6(kc-1)
236  lgs_a1(kc) = 1.0_dp+ct5(kc)*(ct6(kc)+ct6(kc-1))
237  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
238  -ct5(kc)*ct6(kc)
239 
240 #elif (ADV_VERT==2 || ADV_VERT==3)
241 
242  lgs_a0(kc) &
243  = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
244  -ct5(kc)*ct6(kc-1)
245  lgs_a1(kc) &
246  = 1.0_dp &
247  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
248  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
249  +ct5(kc)*(ct6(kc)+ct6(kc-1))
250  lgs_a2(kc) &
251  = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
252  -ct5(kc)*ct6(kc)
253 
254 #endif
255 
256 #if (ADV_HOR==2 || ADV_HOR==3)
257 
258  lgs_b(kc) = temp_c(kc,j,i) + ct7(kc) &
259  -dtt_2dxi* &
260  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
261  *(temp_c(kc,j,i+1)-temp_c(kc,j,i)) &
262  *insq_g11_sgx(j,i) &
263  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
264  *(temp_c(kc,j,i)-temp_c(kc,j,i-1)) &
265  *insq_g11_sgx(j,i-1) ) &
266  -dtt_2deta* &
267  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
268  *(temp_c(kc,j+1,i)-temp_c(kc,j,i)) &
269  *insq_g22_sgy(j,i) &
270  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
271  *(temp_c(kc,j,i)-temp_c(kc,j-1,i)) &
272  *insq_g22_sgy(j-1,i) )
273 
274 #elif ADV_HOR==4
275 
276  lgs_b(kc) = temp_c(kc,j,i) + ct7(kc) &
277  -dtt_dxi *(ftx_c_r(kc)-ftx_c_l(kc)) &
278  -dtt_deta*(fty_c_r(kc)-fty_c_l(kc))
279 
280 #endif
281 
282 end do
283 
284 kc=kcmax
285 lgs_a0(kc) = 0.0_dp
286 lgs_a1(kc) = 1.0_dp
287 lgs_b(kc) = temp_s(j,i)
288 
289 !-------- Solve system of linear equations --------
290 
291 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
292 
293 !-------- Assign the result --------
294 
295 do kc=0, kcmax
296  temp_c_neu(kc,j,i) = lgs_x(kc)
297 end do
298 
299 !-------- Set water content in the non-existing temperate layer
300 ! to zero --------
301 
302 do kt=0, ktmax
303  omega_t_neu(kt,j,i) = 0.0_dp
304 end do
305 
306 !-------- Water drainage from the non-existing temperate layer --------
307 
308 q_tld(j,i) = 0.0_dp
309 
310 !-------- Set up the equations for the age of cold ice --------
311 
312 #if (ADV_HOR==3 && ADV_VERT==3)
313 
314 call calc_age_c(dtime_temp, i, j) ! Algorithm by B. Muegge
315 
316 #else
317 
318 #if ADV_VERT==1
319 
320 kc=0
321 lgs_a1(kc) = 1.0_dp
322 lgs_a2(kc) = -1.0_dp
323 lgs_b(kc) = m_age*h_c(j,i)*ai4*mean_accum_inv
324 
325 #elif ADV_VERT==2
326 
327 kc=0
328 lgs_a1(kc) = 1.0_dp - adv_vert_sg(kc)
329 lgs_a2(kc) = adv_vert_sg(kc)
330 lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
331  -dtt_2dxi* &
332  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
333  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
334  *insq_g11_sgx(j,i) &
335  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
336  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
337  *insq_g11_sgx(j,i-1) ) &
338  -dtt_2deta* &
339  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
340  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
341  *insq_g22_sgy(j,i) &
342  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
343  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
344  *insq_g22_sgy(j-1,i) )
345 
346 #endif
347 
348 do kc=1, kcmax-1
349 
350 #if ADV_VERT==1
351 
352  lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
353  -ci1(kc)*ci2(kc-1)
354  lgs_a1(kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
355  lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
356  -ci1(kc)*ci2(kc)
357 
358 #elif ADV_VERT==2
359 
360  lgs_a0(kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
361  lgs_a1(kc) = 1.0_dp &
362  +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
363  -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
364  lgs_a2(kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
365 
366 #endif
367 
368 #if ADV_HOR==2
369 
370  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
371  -dtt_2dxi* &
372  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
373  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
374  *insq_g11_sgx(j,i) &
375  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
376  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
377  *insq_g11_sgx(j,i-1) ) &
378  -dtt_2deta* &
379  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
380  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
381  *insq_g22_sgy(j,i) &
382  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
383  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
384  *insq_g22_sgy(j-1,i) )
385 
386 #elif ADV_HOR==4
387 
388  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
389  -dtt_dxi *(fax_c_r(kc)-fax_c_l(kc)) &
390  -dtt_deta*(fay_c_r(kc)-fay_c_l(kc))
391 
392 #endif
393 
394 end do
395 
396 #if ADV_VERT==1
397 
398 kc=kcmax
399 if (as_perp(j,i).ge.zero) then
400  lgs_a0(kc) = 0.0_dp
401  lgs_a1(kc) = 1.0_dp
402  lgs_b(kc) = 0.0_dp
403 else
404  lgs_a0(kc) = -1.0_dp
405  lgs_a1(kc) = 1.0_dp
406  lgs_b(kc) = 0.0_dp
407 end if
408 
409 #elif ADV_VERT==2
410 
411 kc=kcmax
412 if (as_perp(j,i).ge.zero) then
413  lgs_a0(kc) = 0.0_dp
414  lgs_a1(kc) = 1.0_dp
415  lgs_b(kc) = 0.0_dp
416 else
417  lgs_a0(kc) = -adv_vert_sg(kc-1)
418  lgs_a1(kc) = 1.0_dp + adv_vert_sg(kc-1)
419  lgs_b(kc) = age_c(kc,j,i) + dtime_temp &
420  -dtt_2dxi* &
421  ( (vx_c(kc,j,i)-abs(vx_c(kc,j,i))) &
422  *(age_c(kc,j,i+1)-age_c(kc,j,i)) &
423  *insq_g11_sgx(j,i) &
424  +(vx_c(kc,j,i-1)+abs(vx_c(kc,j,i-1))) &
425  *(age_c(kc,j,i)-age_c(kc,j,i-1)) &
426  *insq_g11_sgx(j,i-1) ) &
427  -dtt_2deta* &
428  ( (vy_c(kc,j,i)-abs(vy_c(kc,j,i))) &
429  *(age_c(kc,j+1,i)-age_c(kc,j,i)) &
430  *insq_g22_sgy(j,i) &
431  +(vy_c(kc,j-1,i)+abs(vy_c(kc,j-1,i))) &
432  *(age_c(kc,j,i)-age_c(kc,j-1,i)) &
433  *insq_g22_sgy(j-1,i) )
434 end if
435 
436 #endif
437 
438 !-------- Solve system of linear equations --------
439 
440 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
441 
442 !-------- Assign the result,
443 ! restriction to interval [0, AGE_MAX yr] --------
444 
445 do kc=0, kcmax
446 
447  age_c_neu(kc,j,i) = lgs_x(kc)
448 
449  if (age_c_neu(kc,j,i).lt.(age_min*year_sec)) &
450  age_c_neu(kc,j,i) = 0.0_dp
451  if (age_c_neu(kc,j,i).gt.(age_max*year_sec)) &
452  age_c_neu(kc,j,i) = age_max*year_sec
453 
454 end do
455 
456 !-------- Age of the ice in the non-existing temperate layer --------
457 
458 do kt=0, ktmax
459  age_t_neu(kt,j,i) = age_c_neu(0,j,i)
460 end do
461 
462 #endif
463 
464 end subroutine calc_temp_ssa
465 !
Generated on Mon May 20 2013 10:18:41 for SICOPOLIS V3.0 by   doxygen 1.8.1