SICOPOLIS V5-dev  Revision 1173
boundary_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : b o u n d a r y _ m
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-2017 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 module boundary_m
42 
43  use sico_types_m
45  use sico_vars_m
46 
47  implicit none
48 
49  public
50 
51 contains
52 
53 !-------------------------------------------------------------------------------
54 !> Main routine of boundary_m:
55 !! Mars Atmosphere-Ice Coupler MAIC-1.5:
56 !! Computation of the surface temperature (must be less than 0 deg C)
57 !! and of the accumulation-ablation rate for the south polar cap of Mars.
58 !! Computation of the geothermal heat flux.
59 !<------------------------------------------------------------------------------
60 subroutine boundary(time, dtime, dxi, deta, &
61  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
62 
64 
65  use output_m, only : borehole
66 
67 #if (TSURFACE==6)
68  use mars_instemp_m
69 #endif
70 
71 #if ((MARGIN==2) \
72  && (marine_ice_formation==2) \
73  && (marine_ice_calving==9))
75 #endif
76 
77 implicit none
78 
79 real(dp), intent(in) :: time, dtime, dxi, deta
80 
81 real(dp), intent(out) :: delta_ts, glac_index, dzsl_dtau, z_mar
82 real(dp), intent(inout) :: z_sl
83 
84 ! Further return variables
85 ! (defined as global variables in module sico_variables_m):
86 !
87 ! accum(j,i), evap(j,i), runoff(j,i), as_perp(j,i),
88 ! calv_grounded(j,i), temp_s(j,i)
89 
90 integer(i4b) :: i, j
91 integer(i4b) :: i_gr, i_kl
92 integer(i4b) :: ndata_insol
93 real(dp) :: z_sl_old
94 real(dp) :: z_sl_min, t1, t2, t3, t4, t5, t6
95 real(dp) :: time_gr, time_kl
96 real(dp) :: z_sle_present, z_sle_help
97 real(dp), dimension(0:JMAX,0:IMAX) :: temp_ma
98 real(dp) :: eld, g_mb, accum_factor
99 real(dp) :: erosion_chasm
100 real(dp), dimension(0:JMAX,0:IMAX) :: et
101 real(dp) :: t_obliq_main, t_obliq_mod, &
102  obliq_ampl_max, obliq_ampl_min, obliq_ampl, obliq0, obliq, &
103  ecc, ecc0
104 real(dp) :: insol_ma_90_now, insol_ma_90_present, &
105  obl_now, obl_present, ecc_now, ecc_present, &
106  ave_now, ave_present, cp_now, cp_present, &
107  temp_ma_90, temp_ma_90_present, zs_90
108 real(dp), dimension(0:JMAX,0:IMAX) :: dist
109 real(dp) :: ave_data_i_kl, ave_data_i_gr
110 real(dp) :: q_geo_chasm
111 logical, dimension(0:JMAX,0:IMAX) :: check_point
112 logical, save :: firstcall = .true.
113 
114 #if (TSURFACE==6)
115 type(ins) :: temp_now, temp_present
116 #endif
117 
118 real(dp), parameter :: &
119  time_present = 0.0_dp, & ! Present time [s]
120  zs_90_present = 3.85e+03_dp ! Present elevation of the
121  ! south pole [m]
122 real(dp), parameter :: &
123  sol = 590.0_dp, & ! Solar constant [W/m2]
124  epsilon = 1.0_dp, & ! Emissivity
125  sigma = 5.67e-08_dp ! Stefan-Boltzmann constant [W/(m2*K4)]
126 
127 real(dp), parameter :: &
128  lambda_h2o = 2.86e+06_dp, & ! Latent heat [J/kg]
129  r_h2o = 461.5_dp, & ! Gas constant [J/(kg*K)]
130  temp0_ref = 173.0_dp ! Atmopheric reference temperature [K]
131 
132 real(dp), parameter :: &
133  temp_s_min = -128.0_dp ! Minimum ice-surface temperature
134  ! (sublimation temperature of CO2) [C]
135 
136 !-------- Initialization of intent(out) variables --------
137 
138 z_sl_old = z_sl
139 
140 delta_ts = 0.0_dp
141 glac_index = 0.0_dp
142 z_sl = 0.0_dp
143 dzsl_dtau = 0.0_dp
144 z_mar = 0.0_dp
145 
146 !-------- Surface-temperature deviation from present values --------
147 
148 #if (TSURFACE==1)
149 delta_ts = delta_ts0
150 ! ! Steady state with prescribed constant
151 ! ! air-temperature deviation
152 #elif (TSURFACE==3)
153 delta_ts = -sine_amplit &
154  *cos(2.0_dp*pi*time/(sine_period*year_sec)) &
155  +sine_amplit
156 ! ! Sinusoidal air-temperature forcing
157 #elif (TSURFACE==4)
158 
159 ! ------ Obliquity (main cycle and first modulation) and eccentricity
160 
161 t_obliq_main = 1.25e+05_dp *year_sec
162 t_obliq_mod = 1.3e+06_dp *year_sec
163 
164 obliq0 = 25.2_dp *pi_180
165 obliq_ampl_max = 10.0_dp *pi_180
166 obliq_ampl_min = 2.5_dp *pi_180
167 
168 obliq_ampl = 0.5_dp*(obliq_ampl_max+obliq_ampl_min) &
169  -0.5_dp*(obliq_ampl_max-obliq_ampl_min) &
170  *cos(2.0_dp*pi*time/t_obliq_mod)
171 
172 obliq = obliq0 + obliq_ampl*sin(2.0_dp*pi*time/t_obliq_main)
173 
174 ecc = 0.0_dp ! Values not correct, but influence
175 ecc0 = 0.0_dp ! on mean-annual south-polar insolation is negligible
176 
177 ! ------ Mean annual insolation at the south pole
178 
179 insol_ma_90_now = (sol/pi)*sin(obliq)/sqrt(1.0_dp-ecc**2)
180 
181 insol_ma_90_present = (sol/pi)*sin(obliq0)/sqrt(1.0_dp-ecc0**2)
182  ! present value
183 
184 #elif (TSURFACE==5 || TSURFACE==6)
185 
186 ! ------ Mean annual insolation at the south pole
187 
189 
190 if (time/year_sec.lt.real(insol_time_min,dp)) then
191 
192  insol_ma_90_now = insol_ma_90(0)
193  obl_now = obl_data(0)
194  ecc_now = ecc_data(0)
195  ave_now = ave_data(0)
196  cp_now = cp_data(0)
197 
198 else if (time/year_sec.lt.real(insol_time_max,dp)) then
199 
200  i_kl = floor(((time/year_sec) &
201  -real(insol_time_min,dp))/real(insol_time_stp,dp))
202  i_kl = max(i_kl, 0)
203 
204  i_gr = ceiling(((time/year_sec) &
205  -real(insol_time_min,dp))/real(insol_time_stp,dp))
206  i_gr = min(i_gr, ndata_insol)
207 
208  if (i_kl.eq.i_gr) then
209 
210  insol_ma_90_now = insol_ma_90(i_kl)
211  obl_now = obl_data(i_kl)
212  ecc_now = ecc_data(i_kl)
213  ave_now = ave_data(i_kl)
214  cp_now = cp_data(i_kl)
215 
216  else
217 
218  time_kl = (insol_time_min + i_kl*insol_time_stp) *year_sec
219  time_gr = (insol_time_min + i_gr*insol_time_stp) *year_sec
220 
221  insol_ma_90_now = insol_ma_90(i_kl) &
222  +(insol_ma_90(i_gr)-insol_ma_90(i_kl)) &
223  *(time-time_kl)/(time_gr-time_kl)
224  obl_now = obl_data(i_kl) &
225  +(obl_data(i_gr)-obl_data(i_kl)) &
226  *(time-time_kl)/(time_gr-time_kl)
227  ecc_now = ecc_data(i_kl) &
228  +(ecc_data(i_gr)-ecc_data(i_kl)) &
229  *(time-time_kl)/(time_gr-time_kl)
230 
231  if ( abs(ave_data(i_gr)-ave_data(i_kl)) < pi ) then ! regular case
232  ave_data_i_kl = ave_data(i_kl)
233  ave_data_i_gr = ave_data(i_gr)
234  else
235  if ( ave_data(i_gr) > ave_data(i_kl) ) then
236  ave_data_i_kl = ave_data(i_kl) + 2.0_dp*pi
237  ave_data_i_gr = ave_data(i_gr)
238  else
239  ave_data_i_kl = ave_data(i_kl)
240  ave_data_i_gr = ave_data(i_gr) + 2.0_dp*pi
241  end if
242  end if
243 
244  ave_now = ave_data_i_kl &
245  +(ave_data_i_gr-ave_data_i_kl) &
246  *(time-time_kl)/(time_gr-time_kl)
247 
248  cp_now = cp_data(i_kl) &
249  +(cp_data(i_gr)-cp_data(i_kl)) &
250  *(time-time_kl)/(time_gr-time_kl)
251  ! linear interpolation of the data
252  end if
253 
254 else
255 
256  insol_ma_90_now = insol_ma_90(ndata_insol)
257  obl_now = obl_data(ndata_insol)
258  ecc_now = ecc_data(ndata_insol)
259  ave_now = ave_data(ndata_insol)
260  cp_now = cp_data(ndata_insol)
261 
262 end if
263 
264 ! ---- Present value
265 
266 if (time_present/year_sec.lt.real(insol_time_min,dp)) then
267 
268  insol_ma_90_present = insol_ma_90(0)
269  obl_present = obl_data(0)
270  ecc_present = ecc_data(0)
271  ave_present = ave_data(0)
272  cp_present = cp_data(0)
273 
274 else if (time_present/year_sec.lt.real(insol_time_max,dp)) then
275 
276  i_kl = floor(((time_present/year_sec) &
277  -real(insol_time_min,dp))/real(insol_time_stp,dp))
278  i_kl = max(i_kl, 0)
279 
280  i_gr = ceiling(((time_present/year_sec) &
281  -real(insol_time_min,dp))/real(insol_time_stp,dp))
282  i_gr = min(i_gr, ndata_insol)
283 
284  if (i_kl.eq.i_gr) then
285 
286  insol_ma_90_present = insol_ma_90(i_kl)
287  obl_present = obl_data(i_kl)
288  ecc_present = ecc_data(i_kl)
289  ave_present = ave_data(i_kl)
290  cp_present = cp_data(i_kl)
291 
292  else
293 
294  time_kl = (insol_time_min + i_kl*insol_time_stp) *year_sec
295  time_gr = (insol_time_min + i_gr*insol_time_stp) *year_sec
296 
297  insol_ma_90_present = insol_ma_90(i_kl) &
298  +(insol_ma_90(i_gr)-insol_ma_90(i_kl)) &
299  *(time_present-time_kl)/(time_gr-time_kl)
300  obl_present = obl_data(i_kl) &
301  +(obl_data(i_gr)-obl_data(i_kl)) &
302  *(time_present-time_kl)/(time_gr-time_kl)
303  ecc_present = ecc_data(i_kl) &
304  +(ecc_data(i_gr)-ecc_data(i_kl)) &
305  *(time_present-time_kl)/(time_gr-time_kl)
306 
307  if ( abs(ave_data(i_gr)-ave_data(i_kl)) < pi ) then ! regular case
308  ave_data_i_kl = ave_data(i_kl)
309  ave_data_i_gr = ave_data(i_gr)
310  else
311  if ( ave_data(i_gr) > ave_data(i_kl) ) then
312  ave_data_i_kl = ave_data(i_kl) + 2.0_dp*pi
313  ave_data_i_gr = ave_data(i_gr)
314  else
315  ave_data_i_kl = ave_data(i_kl)
316  ave_data_i_gr = ave_data(i_gr) + 2.0_dp*pi
317  end if
318  end if
319 
320  ave_present = ave_data_i_kl &
321  +(ave_data_i_gr-ave_data_i_kl) &
322  *(time_present-time_kl)/(time_gr-time_kl)
323 
324  cp_present = cp_data(i_kl) &
325  +(cp_data(i_gr)-cp_data(i_kl)) &
326  *(time_present-time_kl)/(time_gr-time_kl)
327  ! linear interpolation of the data
328  end if
329 
330 else
331 
332  insol_ma_90_present = insol_ma_90(ndata_insol)
333  obl_present = obl_data(ndata_insol)
334  ecc_present = ecc_data(ndata_insol)
335  ave_present = ave_data(ndata_insol)
336  cp_present = cp_data(ndata_insol)
337 
338 end if
339 
340 #endif
341 
342 #if (TSURFACE==4 || TSURFACE==5)
343 
344 ! ------ Mean-annual surface temperature at the south pole
345 
346 temp_ma_90 = sqrt(sqrt(insol_ma_90_now*(1.0_dp-albedo)/(epsilon*sigma))) &
347  -273.15_dp ! K -> C
348 
349 ! ---- Present value
350 
351 temp_ma_90_present &
352  = sqrt(sqrt(insol_ma_90_present*(1.0_dp-albedo)/(epsilon*sigma))) &
353  -273.15_dp ! K -> C
354 
355 ! ------ Surface-temperature deviation at the south pole
356 
357 delta_ts = temp_ma_90 - temp_ma_90_present
358 
359 #elif (TSURFACE==6)
360 
361 ! ------ Mean-annual surface temperature at the south pole
362 
363 call setinstemp(temp_now, &
364  ecc = ecc_now, ave = ave_now*pi_180_inv, &
365  obl = obl_now*pi_180_inv, sa = albedo, ct = 148.7_dp)
366 
367 temp_ma_90 = instam(temp_now, -90.0_dp) + (delta_ts0) &
368  - 273.15_dp ! K -> C
369 
370 ! ---- Present value
371 
372 call setinstemp(temp_present, &
373  ecc = ecc_present, ave = ave_present*pi_180_inv, &
374  obl = obl_present*pi_180_inv, sa = albedo, ct = 148.7_dp)
375 
376 temp_ma_90_present = instam(temp_present, -90.0_dp) + (delta_ts0) &
377  - 273.15_dp ! K -> C
378 
379 ! ------ Surface-temperature deviation at the south pole
380 
381 delta_ts = temp_ma_90 - temp_ma_90_present
382 
383 #endif
384 
385 !-------- Sea level --------
386 
387 #if (SEA_LEVEL==1)
388 ! ------ constant sea level
389 z_sl = z_sl0
390 
391 #elif (SEA_LEVEL==2)
392 
393 stop ' >>> boundary: SEA_LEVEL==2 not allowed for smars application!'
394 
395 #elif (SEA_LEVEL==3)
396 
397 stop ' >>> boundary: SEA_LEVEL==3 not allowed for smars application!'
398 
399 #endif
400 
401 ! ------ Time derivative of the sea level
402 
403 if ( z_sl_old > -999999.9_dp ) then
404  dzsl_dtau = (z_sl-z_sl_old)/dtime
405 else ! only dummy value for z_sl_old available
406  dzsl_dtau = 0.0_dp
407 end if
408 
409 ! ------ Minimum bedrock elevation for extent of marine ice
410 
411 #if (MARGIN==2)
412 
413 #if (MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3)
414 z_mar = z_mar
415 #elif (MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5)
416 z_mar = fact_z_mar*z_sl
417 #elif (MARINE_ICE_CALVING==6 || MARINE_ICE_CALVING==7)
418 if (z_sl >= -80.0_dp) then
419  z_mar = 2.5_dp*z_sl
420 else
421  z_mar = 10.25_dp*(z_sl+80.0_dp)-200.0_dp
422 end if
423 z_mar = fact_z_mar*z_mar
424 #endif
425 
426 #endif
427 
428 ! ------ Update of the mask according to the sea level
429 
430 ! ---- Check all sea and floating-ice points and their direct
431 ! neighbours
432 
433 do i=0, imax
434 do j=0, jmax
435  check_point(j,i) = .false.
436 end do
437 end do
438 
439 do i=1, imax-1
440 do j=1, jmax-1
441  if (maske(j,i).ge.2) then
442  check_point(j ,i ) = .true.
443  check_point(j ,i+1) = .true.
444  check_point(j ,i-1) = .true.
445  check_point(j+1,i ) = .true.
446  check_point(j-1,i ) = .true.
447  end if
448 end do
449 end do
450 
451 do i=1, imax-1
452 do j=1, jmax-1
453  if (check_point(j,i)) then
454  maske_neu(j,i) = mask_update_sea_level(z_sl, i, j)
455  end if
456 end do
457 end do
458 
459 ! ---- Assign new values of the mask
460 
461 do i=1, imax-1
462 do j=1, jmax-1
463  if (check_point(j,i)) then
464  maske(j,i) = maske_neu(j,i)
465  end if
466 end do
467 end do
468 
469 !-------- Surface air temperature !
470 
471 call borehole(zs, 0.0_dp, 0.0_dp, dxi, deta, 'grid', zs_90)
472  ! zs_90: elevation of the south pole
473 
474 do i=0, imax
475 do j=0, jmax
476 
477 ! ------ Present mean annual air temperature temp_ma
478 
479 #if (TSURFACE==1 || TSURFACE==2 || TSURFACE==3)
480 
481  temp_ma_90_present = temp0_ma_90s
482 
483  temp_ma(j,i) = temp_ma_90_present &
484  + (gamma_ma)*(zs(j,i)-zs_90_present) &
485  + (c_ma*pi_180_inv)*((90.0_dp*pi_180)-abs(phi(j,i)))
486 
487 #elif (TSURFACE==4 || TSURFACE==5)
488 
489  temp_ma(j,i) = temp_ma_90_present &
490  + (gamma_ma)*(zs(j,i)-zs_90_present) &
491  + (c_ma*pi_180_inv)*((90.0_dp*pi_180)-abs(phi(j,i)))
492 
493 #endif
494 
495 ! ------ Correction with deviation delta_ts
496 
497 #if (TSURFACE==1 || TSURFACE==2 || TSURFACE==3)
498 
499  temp_ma(j,i) = temp_ma(j,i) + delta_ts
500 
501 #elif (TSURFACE==4 || TSURFACE==5)
502 
503  temp_ma(j,i) = temp_ma(j,i) + delta_ts &
504  - (gamma_ma)*(zs_90-zs_90_present)
505 
506 #endif
507 
508 ! ------ Direct computation of the mean annual air temperature temp_ma
509 
510 #if (TSURFACE==6)
511 
512  temp_ma(j,i) = instam(temp_now, phi(j,i)*pi_180_inv) + (delta_ts0) &
513  - 273.15_dp ! K -> C
514 
515 #endif
516 
517 end do
518 end do
519 
520 !-------- Accumulation-ablation function as_perp --------
521 
522 ! ------ Accumulation
523 
524 #if (ACCSURFACE==1)
525  accum_factor = 1.0_dp+gamma_s*delta_ts
526 #elif (ACCSURFACE==2)
527  accum_factor = exp(gamma_s*delta_ts)
528 #elif (ACCSURFACE==3)
529  accum_factor = exp(lambda_h2o/(r_h2o*temp0_ref) &
530  -lambda_h2o/(r_h2o*(temp0_ref+delta_ts)))
531 #endif
532 
533 do i=0, imax
534 do j=0, jmax
535  accum(j,i) = accum_present(j,i)*accum_factor
536  if (accum(j,i).lt.0.0_dp) stop ' >>> boundary: Negative accumulation rate!'
537 end do
538 end do
539 
540 ! ------ Ablation
541 
542 #if (ABLSURFACE==1 || ABLSURFACE==2)
543 
544 eld = eld_0 *1.0e+03_dp ! km -> m
545 
546 #if (ACC_UNIT==1)
547 g_mb = g_0 *(1.0e-06_dp/year_sec)*(rho_w/rho_i) &
548  *(1.0_dp/(1.0_dp-frac_dust))
549  ! [mm/a water equiv.]/km -> [m/s (ice+dust) equiv.]/m
550 #elif (ACC_UNIT==2)
551 g_mb = g_0 *(1.0e-06_dp/year_sec)
552  ! [mm/a (ice+dust) equiv.]/km -> [m/s (ice+dust) equiv.]/m
553 #endif
554 
555 #if (ABLSURFACE==1)
556 eld = eld*(1.0_dp+gamma_eld*delta_ts)
557 g_mb = g_mb*(gamma_g*accum_factor)
558 #elif (ABLSURFACE==2)
559 eld = eld*exp(gamma_eld*delta_ts)
560 g_mb = g_mb*(gamma_g*accum_factor)
561 #endif
562 
563 if (eld.lt.eps) stop ' >>> boundary: Negative equilibrium line distance eld!'
564 
565 do i=0, imax
566 do j=0, jmax
567  dist(j,i) = sqrt( xi(i)**2 + eta(j)**2 )
568 end do
569 end do
570 
571 #endif
572 
573 #if (CHASM==2)
574 
575 #if (ACC_UNIT==1)
576 erosion_chasm = erosion_chasm *(1.0e-03_dp/year_sec)*(rho_w/rho_i) &
577  *(1.0_dp/(1.0_dp-frac_dust))
578  ! [mm/a water equiv.] -> [m/s (ice+dust) equiv.]
579 #elif (ACC_UNIT==2)
580 erosion_chasm = erosion_chasm *(1.0e-03_dp/year_sec)
581  ! [mm/a (ice+dust) equiv.] -> [m/s (ice+dust) equiv.]
582 #endif
583 
584 #endif
585 
586 do i=0, imax
587 do j=0, jmax
588 
589 #if (ABLSURFACE==1 || ABLSURFACE==2)
590 
591  as_perp(j,i) = g_mb*(eld-dist(j,i))
592  if (as_perp(j,i).gt.accum(j,i)) as_perp(j,i) = accum(j,i)
593 
594 #endif
595 
596  evap(j,i) = accum(j,i) - as_perp(j,i)
597  runoff(j,i) = 0.0_dp
598 
599 #if (CHASM==2)
600 
601 ! ---- Correction for additional wind erosion in the chasm area
602 
603  if ( (maske_chasm(j,i) == 7) &
604  .and.(time >= time_chasm_init) &
605  .and.(time <= time_chasm_end) ) then ! active chasm area
606  evap(j,i) = evap(j,i) + erosion_chasm
607  as_perp(j,i) = accum(j,i) - evap(j,i)
608  end if
609 
610 #endif
611 
612 end do
613 end do
614 
615 !-------- Ice-surface temperature (10-m firn temperature) temp_s --------
616 
617 temp_s = min(temp_ma, -eps) ! Cut-off of positive air temperatures
618 
619 temp_s = max(temp_s, temp_s_min) ! Cut-off of air temperatures below the
620  ! sublimation temperature of CO2
621 
622 !-------- Calving rate of grounded ice --------
623 
624 calv_grounded = 0.0_dp
625 
626 #if ((MARGIN==2) \
627  && (marine_ice_formation==2) \
628  && (marine_ice_calving==9))
629 
630 call calving_underwater_ice(z_sl)
632 
633 #endif
634 
635 !-------- Geothermal heat flux --------
636 
637 #if (CHASM==1)
638 
640 
641 #elif (CHASM==2)
642 
643 q_geo_chasm = q_geo_chasm *1.0e-03_dp ! mW/m2 -> W/m2
644 
645 do i=0, imax
646 do j=0, jmax
647 
648  if ( (maske_chasm(j,i) == 7) &
649  .and.(time >= time_chasm_init) &
650  .and.(time <= time_chasm_end) ) then ! active chasm area
651  q_geo(j,i) = q_geo_chasm
652  else
653  q_geo(j,i) = q_geo_normal(j,i)
654  end if
655 
656 end do
657 end do
658 
659 #endif
660 
661 if (firstcall) firstcall = .false.
662 
663 end subroutine boundary
664 
665 !-------------------------------------------------------------------------------
666 
667 end module boundary_m
668 !
real(dp), dimension(0:jmax, 0:imax) phi
phi(j,i): Geographic latitude of grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) q_geo
q_geo(j,i): Geothermal heat flux
integer(i4b) insol_time_stp
insol_time_stp: Time step of the data values for the insolation etc.
Definition: sico_vars_m.F90:45
real(dp), dimension(0:100000) ecc_data
ecc_data(n): Data values for the eccentricity
Definition: sico_vars_m.F90:54
real(dp) rho_w
RHO_W: Density of pure water.
real(dp), parameter eps
eps: Small number
real(dp), dimension(0:jmax, 0:imax) q_geo_normal
q_geo_normal(j,i): Geothermal heat flux for normal (non-chasma) areas
Definition: sico_vars_m.F90:74
real(dp), dimension(0:jmax, 0:imax) calv_grounded
calv_grounded(j,i): Calving rate of grounded ice
integer(i2b), dimension(0:jmax, 0:imax) maske
maske(j,i): Ice-land-ocean mask. 0: grounded ice, 1: ice-free land, 2: ocean, 3: floating ice ...
integer(i4b) insol_time_min
insol_time_min: Minimum time of the data values for the insolation etc.
Definition: sico_vars_m.F90:43
real(dp), dimension(0:jmax, 0:imax) as_perp
as_perp(j,i): Accumulation-ablation function at the ice surface (SMB)
real(dp), dimension(0:jmax, 0:imax) runoff
runoff(j,i): Runoff rate at the ice surface
Martian surface temperatures.
integer(i2b) function, public mask_update_sea_level(z_sl, i, j)
Main function of mask_update_m: Update of the ice-land-ocean mask due to changes of the sea level...
real(dp), parameter pi_180_inv
pi_180_inv: 180 divided by pi (-> rad to deg)
real(dp), dimension(0:jmax, 0:imax) temp_s
temp_s(j,i): Ice surface temperature
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars_m.F90:35
integer(i4b) insol_time_max
insol_time_max: Maximum time of the data values for the insolation etc.
Definition: sico_vars_m.F90:47
subroutine, public borehole(field, x_pos, y_pos, dxi, deta, ch_grid, field_val)
Computation of an arbitrary field quantity for a given borehole position x_pos, y_pos by weighed aver...
Definition: output_m.F90:5586
real(dp), dimension(0:jmax, 0:imax) calv_uw_ice
real(dp), parameter pi_180
pi_180: pi divided by 180 (-> deg to rad)
Calving of "underwater ice".
integer(i2b), dimension(0:jmax, 0:imax) maske_chasm
maske_chasm(j,i): Chasma mask. 0: grounded ice, 1: ice-free land (normal area), 7: chasma area ...
Definition: sico_vars_m.F90:68
Declarations of kind types for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) zs
zs(j,i): Coordinate z of the surface topography
real(dp), dimension(0:jmax, 0:imax) evap
evap(j,i): Evaporation rate at the ice surface
real(dp), dimension(0:jmax, 0:imax) accum
accum(j,i): Accumulation rate at the ice surface (includes liquid precipitation = rainfall!) ...
real(dp) time_chasm_end
time_chasm_end: Final time for active chasma area
Definition: sico_vars_m.F90:72
real(dp) eld
eld: Equilibrium line distance
Definition: sico_vars_m.F90:56
real(dp), dimension(0:jmax, 0:imax) accum_present
accum_present(j,i): Present-day accumulation rate at the ice surface (for EISMINT, ISMIP HEINO and the north and south polar caps of Mars)
real(dp) rho_i
RHO_I: Density of ice.
Definition: sico_vars_m.F90:77
real(dp) function instam(o, phi)
Annual mean temperature at latitude phi.
real(dp) time_chasm_init
time_chasm_init: Initial time for active chasma area
Definition: sico_vars_m.F90:70
Computation of the daily mean surface temperature of Mars based on obliquity, eccentricity and the an...
real(dp), dimension(0:100000) cp_data
cp_data(n): Data values for Laskar&#39;s climate parameter = eccentricity *sin(Laskar&#39;s longitude of peri...
Definition: sico_vars_m.F90:63
integer(i2b), dimension(0:jmax, 0:imax) maske_neu
maske_neu(j,i): New value of maske computed during an integration step
subroutine setinstemp(o, ecc, ave, obl, sma, sa, sac, op, ct)
Main subroutine of module mars_instemp_m.
subroutine boundary(time, dtime, dxi, deta, delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
Main routine of boundary_m: Computation of the surface temperature (must be less than 0 deg C!) and o...
Definition: boundary_m.F90:56
real(dp), dimension(0:100000) obl_data
obl_data(n): Data values for the obliquity
Definition: sico_vars_m.F90:52
real(dp), parameter pi
pi: Constant pi
real(dp), dimension(0:jmax) eta
eta(j): Coordinate eta (= y) of grid point j
real(dp), dimension(0:imax) xi
xi(i): Coordinate xi (= x) of grid point i
subroutine calving_underwater_ice(z_sl)
Main routine: Calving of "underwater ice".
Computation of the surface temperature (must be less than 0 deg C!) and of the accumulation-ablation ...
Definition: boundary_m.F90:37
real(dp), dimension(0:100000) insol_ma_90
insol_ma_90(n): Data values for the mean-annual north- or south-polar insolation
Definition: sico_vars_m.F90:50
Update of the ice-land-ocean mask due to changes of the sea level.
integer, parameter dp
Double-precision reals.
real(dp), dimension(0:100000) ave_data
ave_data(n): Data values for the anomaly of vernal equinox (= 360 deg - Ls of perihelion ) ...
Definition: sico_vars_m.F90:57
Writing of output data on files.
Definition: output_m.F90:36
Declarations of global variables for SICOPOLIS.