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