SICOPOLIS V3.1
 All Classes Files Functions Variables Macros
boundary.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : b o u n d a r y
4 !
5 !> @file
6 !!
7 !! Mars Atmosphere-Ice Coupler MAIC-1.5:
8 !! Computation of the surface temperature (must be less than 0 deg C)
9 !! and of the accumulation-ablation rate for the south polar cap of Mars.
10 !! Computation of the geothermal heat flux.
11 !!
12 !! @section Copyright
13 !!
14 !! Copyright 2009-2013 Ralf Greve
15 !!
16 !! @section License
17 !!
18 !! This file is part of SICOPOLIS.
19 !!
20 !! SICOPOLIS is free software: you can redistribute it and/or modify
21 !! it under the terms of the GNU General Public License as published by
22 !! the Free Software Foundation, either version 3 of the License, or
23 !! (at your option) any later version.
24 !!
25 !! SICOPOLIS is distributed in the hope that it will be useful,
26 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
27 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
28 !! GNU General Public License for more details.
29 !!
30 !! You should have received a copy of the GNU General Public License
31 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
32 !<
33 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
34 
35 !-------------------------------------------------------------------------------
36 !> Mars Atmosphere-Ice Coupler MAIC-1.5:
37 !! Computation of the surface temperature (must be less than 0 deg C)
38 !! and of the accumulation-ablation rate for the south polar cap of Mars.
39 !! Computation of the geothermal heat flux.
40 !<------------------------------------------------------------------------------
41 subroutine boundary(time, dtime, dxi, deta, &
42  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
43 
44 use sico_types
45 use sico_variables
46 
47 #if TSURFACE==6
48 use instemp
49 #endif
50 
51 implicit none
52 
53 real(dp), intent(in) :: time, dtime, dxi, deta
54 
55 real(dp), intent(out) :: delta_ts, glac_index, dzsl_dtau, z_mar
56 real(dp), intent(inout) :: z_sl
57 
58 ! Further return variables
59 ! (defined as global variables in module sico_variables):
60 !
61 ! accum(j,i), evap(j,i), runoff(j,i), as_perp(j,i), temp_s(j,i)
62 
63 integer(i4b) :: i, j
64 integer(i4b) :: i_gr, i_kl
65 integer(i4b) :: ndata_insol
66 real(dp) :: z_sl_old
67 real(dp) :: z_sl_min, t1, t2, t3, t4, t5, t6
68 real(dp) :: time_gr, time_kl
69 real(dp) :: z_sle_present, z_sle_help
70 real(dp), dimension(0:JMAX,0:IMAX) :: temp_ma
71 real(dp) :: eld, g_mb, accum_factor
72 real(dp) :: erosion_chasm
73 real(dp), dimension(0:JMAX,0:IMAX) :: et
74 real(dp) :: t_obliq_main, t_obliq_mod, &
75  obliq_ampl_max, obliq_ampl_min, obliq_ampl, obliq0, obliq, &
76  ecc, ecc0
77 real(dp) :: insol_ma_90_now, insol_ma_90_present, &
78  obl_now, obl_present, ecc_now, ecc_present, &
79  ave_now, ave_present, cp_now, cp_present, &
80  temp_ma_90, temp_ma_90_present, zs_90
81 real(dp), dimension(0:JMAX,0:IMAX) :: dist
82 real(dp) :: ave_data_i_kl, ave_data_i_gr
83 real(dp) :: q_geo_chasm
84 logical, dimension(0:JMAX,0:IMAX) :: check_point
85 
86 #if TSURFACE==6
87 type (ins) :: temp_now, temp_present
88 #endif
89 
90 real(dp), parameter :: &
91  time_present = 0.0_dp, & ! Present time [s]
92  zs_90_present = 3.85e+03_dp ! Present elevation of the
93  ! south pole [m]
94 real(dp), parameter :: &
95  sol = 590.0_dp, & ! Solar constant [W/m2]
96  epsilon = 1.0_dp, & ! Emissivity
97  sigma = 5.67e-08_dp ! Stefan-Boltzmann constant [W/(m2*K4)]
98 
99 real(dp), parameter :: &
100  lambda_h2o = 2.86e+06_dp, & ! Latent heat [J/kg]
101  r_h2o = 461.5_dp, & ! Gas constant [J/(kg*K)]
102  temp0_ref = 173.0_dp ! Atmopheric reference temperature [K]
103 
104 real(dp), parameter :: &
105  temp_s_min = -128.0_dp ! Minimum ice-surface temperature
106  ! (sublimation temperature of CO2) [C]
107 
108 !-------- Initialization of intent(out) variables --------
109 
110 z_sl_old = z_sl
111 
112 delta_ts = 0.0_dp
113 glac_index = 0.0_dp
114 z_sl = 0.0_dp
115 dzsl_dtau = 0.0_dp
116 z_mar = 0.0_dp
117 
118 !-------- Surface-temperature deviation from present values --------
119 
120 #if TSURFACE==1
121 delta_ts = delta_ts0
122 ! ! Steady state with prescribed constant
123 ! ! air-temperature deviation
124 #elif TSURFACE==3
125 delta_ts = -sine_amplit &
126  *cos(2.0_dp*pi*time/(sine_period*year_sec)) &
127  +sine_amplit
128 ! ! Sinusoidal air-temperature forcing
129 #elif TSURFACE==4
130 
131 ! ------ Obliquity (main cycle and first modulation) and eccentricity
132 
133 t_obliq_main = 1.25e+05_dp *year_sec
134 t_obliq_mod = 1.3e+06_dp *year_sec
135 
136 obliq0 = 25.2_dp *pi_180
137 obliq_ampl_max = 10.0_dp *pi_180
138 obliq_ampl_min = 2.5_dp *pi_180
139 
140 obliq_ampl = 0.5_dp*(obliq_ampl_max+obliq_ampl_min) &
141  -0.5_dp*(obliq_ampl_max-obliq_ampl_min) &
142  *cos(2.0_dp*pi*time/t_obliq_mod)
143 
144 obliq = obliq0 + obliq_ampl*sin(2.0_dp*pi*time/t_obliq_main)
145 
146 ecc = 0.0_dp ! Values not correct, but influence
147 ecc0 = 0.0_dp ! on mean-annual south-polar insolation is negligible
148 
149 ! ------ Mean annual insolation at the south pole
150 
151 insol_ma_90_now = (sol/pi)*sin(obliq)/sqrt(1.0_dp-ecc**2)
152 
153 insol_ma_90_present = (sol/pi)*sin(obliq0)/sqrt(1.0_dp-ecc0**2)
154  ! present value
155 
156 #elif ( TSURFACE==5 || TSURFACE==6 )
157 
158 ! ------ Mean annual insolation at the south pole
159 
160 ndata_insol = (insol_time_max-insol_time_min)/insol_time_stp
161 
162 if (time/year_sec.lt.real(insol_time_min,dp)) then
163 
164  insol_ma_90_now = insol_ma_90(0)
165  obl_now = obl_data(0)
166  ecc_now = ecc_data(0)
167  ave_now = ave_data(0)
168  cp_now = cp_data(0)
169 
170 else if (time/year_sec.lt.real(insol_time_max,dp)) then
171 
172  i_kl = floor(((time/year_sec) &
173  -real(insol_time_min,dp))/real(insol_time_stp,dp))
174  i_kl = max(i_kl, 0)
175 
176  i_gr = ceiling(((time/year_sec) &
177  -real(insol_time_min,dp))/real(insol_time_stp,dp))
178  i_gr = min(i_gr, ndata_insol)
179 
180  if (i_kl.eq.i_gr) then
181 
182  insol_ma_90_now = insol_ma_90(i_kl)
183  obl_now = obl_data(i_kl)
184  ecc_now = ecc_data(i_kl)
185  ave_now = ave_data(i_kl)
186  cp_now = cp_data(i_kl)
187 
188  else
189 
190  time_kl = (insol_time_min + i_kl*insol_time_stp) *year_sec
191  time_gr = (insol_time_min + i_gr*insol_time_stp) *year_sec
192 
193  insol_ma_90_now = insol_ma_90(i_kl) &
194  +(insol_ma_90(i_gr)-insol_ma_90(i_kl)) &
195  *(time-time_kl)/(time_gr-time_kl)
196  obl_now = obl_data(i_kl) &
197  +(obl_data(i_gr)-obl_data(i_kl)) &
198  *(time-time_kl)/(time_gr-time_kl)
199  ecc_now = ecc_data(i_kl) &
200  +(ecc_data(i_gr)-ecc_data(i_kl)) &
201  *(time-time_kl)/(time_gr-time_kl)
202 
203  if ( abs(ave_data(i_gr)-ave_data(i_kl)) < pi ) then ! regular case
204  ave_data_i_kl = ave_data(i_kl)
205  ave_data_i_gr = ave_data(i_gr)
206  else
207  if ( ave_data(i_gr) > ave_data(i_kl) ) then
208  ave_data_i_kl = ave_data(i_kl) + 2.0_dp*pi
209  ave_data_i_gr = ave_data(i_gr)
210  else
211  ave_data_i_kl = ave_data(i_kl)
212  ave_data_i_gr = ave_data(i_gr) + 2.0_dp*pi
213  end if
214  end if
215 
216  ave_now = ave_data_i_kl &
217  +(ave_data_i_gr-ave_data_i_kl) &
218  *(time-time_kl)/(time_gr-time_kl)
219 
220  cp_now = cp_data(i_kl) &
221  +(cp_data(i_gr)-cp_data(i_kl)) &
222  *(time-time_kl)/(time_gr-time_kl)
223  ! linear interpolation of the data
224  end if
225 
226 else
227 
228  insol_ma_90_now = insol_ma_90(ndata_insol)
229  obl_now = obl_data(ndata_insol)
230  ecc_now = ecc_data(ndata_insol)
231  ave_now = ave_data(ndata_insol)
232  cp_now = cp_data(ndata_insol)
233 
234 end if
235 
236 ! ---- Present value
237 
238 if (time_present/year_sec.lt.real(insol_time_min,dp)) then
239 
240  insol_ma_90_present = insol_ma_90(0)
241  obl_present = obl_data(0)
242  ecc_present = ecc_data(0)
243  ave_present = ave_data(0)
244  cp_present = cp_data(0)
245 
246 else if (time_present/year_sec.lt.real(insol_time_max,dp)) then
247 
248  i_kl = floor(((time_present/year_sec) &
249  -real(insol_time_min,dp))/real(insol_time_stp,dp))
250  i_kl = max(i_kl, 0)
251 
252  i_gr = ceiling(((time_present/year_sec) &
253  -real(insol_time_min,dp))/real(insol_time_stp,dp))
254  i_gr = min(i_gr, ndata_insol)
255 
256  if (i_kl.eq.i_gr) then
257 
258  insol_ma_90_present = insol_ma_90(i_kl)
259  obl_present = obl_data(i_kl)
260  ecc_present = ecc_data(i_kl)
261  ave_present = ave_data(i_kl)
262  cp_present = cp_data(i_kl)
263 
264  else
265 
266  time_kl = (insol_time_min + i_kl*insol_time_stp) *year_sec
267  time_gr = (insol_time_min + i_gr*insol_time_stp) *year_sec
268 
269  insol_ma_90_present = insol_ma_90(i_kl) &
270  +(insol_ma_90(i_gr)-insol_ma_90(i_kl)) &
271  *(time_present-time_kl)/(time_gr-time_kl)
272  obl_present = obl_data(i_kl) &
273  +(obl_data(i_gr)-obl_data(i_kl)) &
274  *(time_present-time_kl)/(time_gr-time_kl)
275  ecc_present = ecc_data(i_kl) &
276  +(ecc_data(i_gr)-ecc_data(i_kl)) &
277  *(time_present-time_kl)/(time_gr-time_kl)
278 
279  if ( abs(ave_data(i_gr)-ave_data(i_kl)) < pi ) then ! regular case
280  ave_data_i_kl = ave_data(i_kl)
281  ave_data_i_gr = ave_data(i_gr)
282  else
283  if ( ave_data(i_gr) > ave_data(i_kl) ) then
284  ave_data_i_kl = ave_data(i_kl) + 2.0_dp*pi
285  ave_data_i_gr = ave_data(i_gr)
286  else
287  ave_data_i_kl = ave_data(i_kl)
288  ave_data_i_gr = ave_data(i_gr) + 2.0_dp*pi
289  end if
290  end if
291 
292  ave_present = ave_data_i_kl &
293  +(ave_data_i_gr-ave_data_i_kl) &
294  *(time_present-time_kl)/(time_gr-time_kl)
295 
296  cp_present = cp_data(i_kl) &
297  +(cp_data(i_gr)-cp_data(i_kl)) &
298  *(time_present-time_kl)/(time_gr-time_kl)
299  ! linear interpolation of the data
300  end if
301 
302 else
303 
304  insol_ma_90_present = insol_ma_90(ndata_insol)
305  obl_present = obl_data(ndata_insol)
306  ecc_present = ecc_data(ndata_insol)
307  ave_present = ave_data(ndata_insol)
308  cp_present = cp_data(ndata_insol)
309 
310 end if
311 
312 #endif
313 
314 #if (TSURFACE==4 || TSURFACE==5)
315 
316 ! ------ Mean-annual surface temperature at the south pole
317 
318 temp_ma_90 = sqrt(sqrt(insol_ma_90_now*(1.0_dp-albedo)/(epsilon*sigma))) &
319  -273.15_dp ! K -> C
320 
321 ! ---- Present value
322 
323 temp_ma_90_present &
324  = sqrt(sqrt(insol_ma_90_present*(1.0_dp-albedo)/(epsilon*sigma))) &
325  -273.15_dp ! K -> C
326 
327 ! ------ Surface-temperature deviation at the south pole
328 
329 delta_ts = temp_ma_90 - temp_ma_90_present
330 
331 #elif TSURFACE==6
332 
333 ! ------ Mean-annual surface temperature at the south pole
334 
335 call setinstemp(temp_now, &
336  ecc = ecc_now, ave = ave_now*pi_180_inv, &
337  obl = obl_now*pi_180_inv, sa = albedo, ct = 148.7_dp)
338 
339 temp_ma_90 = instam(temp_now, -90.0_dp) + delta_ts0 &
340  - 273.15_dp ! K -> C
341 
342 ! ---- Present value
343 
344 call setinstemp(temp_present, &
345  ecc = ecc_present, ave = ave_present*pi_180_inv, &
346  obl = obl_present*pi_180_inv, sa = albedo, ct = 148.7_dp)
347 
348 temp_ma_90_present = instam(temp_present, -90.0_dp) + delta_ts0 &
349  - 273.15_dp ! K -> C
350 
351 ! ------ Surface-temperature deviation at the south pole
352 
353 delta_ts = temp_ma_90 - temp_ma_90_present
354 
355 #endif
356 
357 !-------- Sea level --------
358 
359 #if SEA_LEVEL==1
360 ! ------ constant sea level
361 z_sl = z_sl0
362 
363 #elif SEA_LEVEL==2
364 
365 stop ' boundary: SEA_LEVEL==2 not allowed for smars application!'
366 
367 #elif SEA_LEVEL==3
368 
369 stop ' boundary: SEA_LEVEL==3 not allowed for smars application!'
370 
371 #endif
372 
373 ! ------ Time derivative of the sea level
374 
375 if ( z_sl_old > -999999.9_dp ) then
376  dzsl_dtau = (z_sl-z_sl_old)/dtime
377 else ! only dummy value for z_sl_old available
378  dzsl_dtau = 0.0_dp
379 end if
380 
381 ! ------ Minimum bedrock elevation for extent of marine ice
382 
383 #if MARGIN==2
384 
385 #if ( MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3 )
386 z_mar = z_mar
387 #elif ( MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 )
388 z_mar = fact_z_mar*z_sl
389 #elif ( MARINE_ICE_CALVING==6 || MARINE_ICE_CALVING==7 )
390 if (z_sl >= -80.0_dp) then
391  z_mar = 2.5_dp*z_sl
392 else
393  z_mar = 10.25_dp*(z_sl+80.0_dp)-200.0_dp
394 end if
395 z_mar = fact_z_mar*z_mar
396 #endif
397 
398 #endif
399 
400 ! ------ Update of the mask according to the sea level
401 
402 ! ---- Check all sea and floating-ice points and their direct
403 ! neighbours
404 
405 do i=0, imax
406 do j=0, jmax
407  check_point(j,i) = .false.
408 end do
409 end do
410 
411 do i=1, imax-1
412 do j=1, jmax-1
413  if (maske(j,i).ge.2) then
414  check_point(j ,i ) = .true.
415  check_point(j ,i+1) = .true.
416  check_point(j ,i-1) = .true.
417  check_point(j+1,i ) = .true.
418  check_point(j-1,i ) = .true.
419  end if
420 end do
421 end do
422 
423 do i=1, imax-1
424 do j=1, jmax-1
425  if (check_point(j,i)) then
426  maske_neu(j,i) = mask_update(z_sl, i, j)
427  end if
428 end do
429 end do
430 
431 ! ---- Assign new values of the mask
432 
433 do i=1, imax-1
434 do j=1, jmax-1
435  if (check_point(j,i)) then
436  maske(j,i) = maske_neu(j,i)
437  end if
438 end do
439 end do
440 
441 !-------- Surface air temperature !
442 
443 call borehole(zs, 0.0_dp, 0.0_dp, dxi, deta, 'grid', zs_90)
444  ! zs_90: elevation of the south pole
445 
446 do i=0, imax
447 do j=0, jmax
448 
449 ! ------ Present mean annual air temperature temp_ma
450 
451 #if (TSURFACE==1 || TSURFACE==2 || TSURFACE==3)
452 
453  temp_ma_90_present = temp0_ma_90s
454 
455  temp_ma(j,i) = temp_ma_90_present &
456  + (gamma_ma)*(zs(j,i)-zs_90_present) &
457  + (c_ma*pi_180_inv)*((90.0_dp*pi_180)-abs(phi(j,i)))
458 
459 #elif (TSURFACE==4 || TSURFACE==5)
460 
461  temp_ma(j,i) = temp_ma_90_present &
462  + (gamma_ma)*(zs(j,i)-zs_90_present) &
463  + (c_ma*pi_180_inv)*((90.0_dp*pi_180)-abs(phi(j,i)))
464 
465 #endif
466 
467 ! ------ Correction with deviation delta_ts
468 
469 #if (TSURFACE==1 || TSURFACE==2 || TSURFACE==3)
470 
471  temp_ma(j,i) = temp_ma(j,i) + delta_ts
472 
473 #elif (TSURFACE==4 || TSURFACE==5)
474 
475  temp_ma(j,i) = temp_ma(j,i) + delta_ts &
476  - (gamma_ma)*(zs_90-zs_90_present)
477 
478 #endif
479 
480 ! ------ Direct computation of the mean annual air temperature temp_ma
481 
482 #if TSURFACE==6
483 
484  temp_ma(j,i) = instam(temp_now, phi(j,i)*pi_180_inv) + delta_ts0 &
485  - 273.15_dp ! K -> C
486 
487 #endif
488 
489 end do
490 end do
491 
492 !-------- Accumulation-ablation function as_perp --------
493 
494 ! ------ Accumulation
495 
496 #if ACCSURFACE==1
497  accum_factor = 1.0_dp+gamma_s*delta_ts
498 #elif ACCSURFACE==2
499  accum_factor = exp(gamma_s*delta_ts)
500 #elif ACCSURFACE==3
501  accum_factor = exp(lambda_h2o/(r_h2o*temp0_ref) &
502  -lambda_h2o/(r_h2o*(temp0_ref+delta_ts)))
503 #endif
504 
505 do i=0, imax
506 do j=0, jmax
507  accum(j,i) = accum_present(j,i)*accum_factor
508  if (accum(j,i).lt.0.0_dp) stop ' boundary: Negative accumulation rate!'
509 end do
510 end do
511 
512 ! ------ Ablation
513 
514 #if ( ABLSURFACE==1 || ABLSURFACE==2)
515 
516 eld = eld_0 *1.0e+03_dp ! km -> m
517 
518 #if ACC_UNIT==1
519 g_mb = g_0 *(1.0e-06_dp/year_sec)*(rho_w/rho_i) &
520  *(1.0_dp/(1.0_dp-frac_dust))
521  ! [mm/a water equiv.]/km -> [m/s (ice+dust) equiv.]/m
522 #elif ACC_UNIT==2
523 g_mb = g_0 *(1.0e-06_dp/year_sec)
524  ! [mm/a (ice+dust) equiv.]/km -> [m/s (ice+dust) equiv.]/m
525 #endif
526 
527 #if ABLSURFACE==1
528 eld = eld*(1.0_dp+gamma_eld*delta_ts)
529 g_mb = g_mb*(gamma_g*accum_factor)
530 #elif ABLSURFACE==2
531 eld = eld*exp(gamma_eld*delta_ts)
532 g_mb = g_mb*(gamma_g*accum_factor)
533 #endif
534 
535 if (eld.lt.eps) stop ' boundary: Negative equilibrium line distance eld!'
536 
537 do i=0, imax
538 do j=0, jmax
539  dist(j,i) = sqrt( xi(i)**2 + eta(j)**2 )
540 end do
541 end do
542 
543 #endif
544 
545 #if CHASM==2
546 
547 #if ACC_UNIT==1
548 erosion_chasm = erosion_chasm *(1.0e-03_dp/year_sec)*(rho_w/rho_i) &
549  *(1.0_dp/(1.0_dp-frac_dust))
550  ! [mm/a water equiv.] -> [m/s (ice+dust) equiv.]
551 #elif ACC_UNIT==2
552 erosion_chasm = erosion_chasm *(1.0e-03_dp/year_sec)
553  ! [mm/a (ice+dust) equiv.] -> [m/s (ice+dust) equiv.]
554 #endif
555 
556 #endif
557 
558 do i=0, imax
559 do j=0, jmax
560 
561 #if ( ABLSURFACE==1 || ABLSURFACE==2 )
562 
563  as_perp(j,i) = g_mb*(eld-dist(j,i))
564  if (as_perp(j,i).gt.accum(j,i)) as_perp(j,i) = accum(j,i)
565 
566 #endif
567 
568  evap(j,i) = accum(j,i) - as_perp(j,i)
569  runoff(j,i) = 0.0_dp
570 
571 #if CHASM==2
572 
573 ! ---- Correction for additional wind erosion in the chasm area
574 
575  if ( (maske_chasm(j,i) == 7) &
576  .and.(time >= time_chasm_init) &
577  .and.(time <= time_chasm_end) ) then ! active chasm area
578  evap(j,i) = evap(j,i) + erosion_chasm
579  as_perp(j,i) = accum(j,i) - evap(j,i)
580  end if
581 
582 #endif
583 
584 end do
585 end do
586 
587 !-------- Ice-surface temperature (10-m firn temperature) temp_s --------
588 
589 temp_s = min(temp_ma, -eps) ! Cut-off of positive air temperatures
590 
591 temp_s = max(temp_s, temp_s_min) ! Cut-off of air temperatures below the
592  ! sublimation temperature of CO2
593 
594 !-------- Geothermal heat flux --------
595 
596 #if CHASM==1
597 
598 q_geo = q_geo_normal
599 
600 #elif CHASM==2
601 
602 q_geo_chasm = q_geo_chasm *1.0e-03_dp ! mW/m2 -> W/m2
603 
604 do i=0, imax
605 do j=0, jmax
606 
607  if ( (maske_chasm(j,i) == 7) &
608  .and.(time >= time_chasm_init) &
609  .and.(time <= time_chasm_end) ) then ! active chasm area
610  q_geo(j,i) = q_geo_chasm
611  else
612  q_geo(j,i) = q_geo_normal(j,i)
613  end if
614 
615 end do
616 end do
617 
618 #endif
619 
620 end subroutine boundary
621 !
Generated on Thu Sep 12 2013 09:44:55 for SICOPOLIS V3.1 by   doxygen 1.8.2