SICOPOLIS V5-dev  Revision 1420
sico_init_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : s i c o _ i n i t _ m
4 !
5 !> @file
6 !!
7 !! Initialisations for SICOPOLIS.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2019 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 !> Initialisations for SICOPOLIS.
34 !<------------------------------------------------------------------------------
35 module sico_init_m
36 
37  use sico_types_m
39  use sico_vars_m
40  use error_m
41 
42  implicit none
43 
44  public
45 
46 contains
47 
48 !-------------------------------------------------------------------------------
49 !> Main routine of sico_init_m: Initialisations for SICOPOLIS.
50 !<------------------------------------------------------------------------------
51 subroutine sico_init(delta_ts, glac_index, &
52  mean_accum, &
53  dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
54  time, time_init, time_end, time_output, &
55  dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
56  z_sl, dzsl_dtau, z_mar, &
57  ndat2d, ndat3d, n_output, &
58  runname)
59 
60  use compare_float_m
64 
65  use read_m, only : read_phys_para, read_kei
66 
67  use boundary_m
69  use calc_enhance_m
71  use calc_vxy_m
72  use calc_vz_m
73  use calc_dxyz_m
75 
76  use output_m
77 
78 implicit none
79 
80 integer(i4b), intent(out) :: ndat2d, ndat3d
81 integer(i4b), intent(out) :: n_output
82 real(dp), intent(out) :: delta_ts, glac_index
83 real(dp), intent(out) :: mean_accum
84 real(dp), intent(out) :: dtime, dtime_temp, dtime_wss, &
85  dtime_out, dtime_ser
86 real(dp), intent(out) :: time, time_init, time_end, time_output(100)
87 real(dp), intent(out) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
88 real(dp), intent(out) :: z_sl, dzsl_dtau, z_mar
89 character(len=100), intent(out) :: runname
90 
91 integer(i4b) :: i, j, kc, kt, kr, m, n, ir, jr
92 integer(i4b) :: ios, ios1, ios2, ios3, ios4
93 integer(i4b) :: ierr
94 integer(i4b) :: ndata_insol
95 real(dp) :: dtime0, dtime_temp0, dtime_wss0, dtime_out0, dtime_ser0
96 real(dp) :: time_init0, time_end0, time_output0(100)
97 real(dp) :: time_chasm_init0, time_chasm_end0
98 real(dp) :: d_dummy
99 character(len=100) :: anfdatname
100 character(len=256) :: filename_with_path
101 character(len=256) :: shell_command
102 character :: ch_dummy
103 logical :: flag_init_output, flag_3d_output
104 
105 character(len=64), parameter :: fmt1 = '(a)', &
106  fmt2 = '(a,i4)', &
107  fmt2a = '(a,i0)', &
108  fmt3 = '(a,es12.4)'
109 
110 character(len= 8) :: ch_imax
111 character(len=128) :: fmt4
112 
113 write(ch_imax, fmt='(i8)') imax
114 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
115 
116 write(unit=6, fmt='(a)') ' '
117 write(unit=6, fmt='(a)') ' -------- sico_init --------'
118 write(unit=6, fmt='(a)') ' '
119 
120 !-------- Name of the computational domain --------
121 
122 #if (defined(ANT))
123 ch_domain_long = 'Antarctica'
124 ch_domain_short = 'ant'
125 
126 #elif (defined(ASF))
127 ch_domain_long = 'Austfonna'
128 ch_domain_short = 'asf'
129 
130 #elif (defined(EMTP2SGE))
131 ch_domain_long = 'EISMINT Phase 2 Simplified Geometry Experiment'
132 ch_domain_short = 'emtp2sge'
133 
134 #elif (defined(GRL))
135 ch_domain_long = 'Greenland'
136 ch_domain_short = 'grl'
137 
138 #elif (defined(NHEM))
139 ch_domain_long = 'Northern hemisphere'
140 ch_domain_short = 'nhem'
141 
142 #elif (defined(SCAND))
143 ch_domain_long = 'Scandinavia and Eurasia'
144 ch_domain_short = 'scand'
145 
146 #elif (defined(TIBET))
147 ch_domain_long = 'Tibet'
148 ch_domain_short = 'tibet'
149 
150 #elif (defined(NMARS))
151 ch_domain_long = 'North polar cap of Mars'
152 ch_domain_short = 'nmars'
153 
154 #elif (defined(SMARS))
155 ch_domain_long = 'South polar cap of Mars'
156 ch_domain_short = 'smars'
157 
158 #elif (defined(XYZ))
159 ch_domain_long = 'XYZ'
160 ch_domain_short = 'xyz'
161 #if (defined(HEINO))
162 ch_domain_long = trim(ch_domain_long)//'/ISMIP HEINO'
163 #endif
164 
165 #else
166 
167 errormsg = ' >>> sico_init: No valid domain specified!'
168 call error(errormsg)
169 
170 #endif
171 
172 !-------- Some initial values --------
173 
174 n_output = 0
175 
176 dtime = 0.0_dp
177 dtime_temp = 0.0_dp
178 dtime_wss = 0.0_dp
179 dtime_out = 0.0_dp
180 dtime_ser = 0.0_dp
181 
182 time = 0.0_dp
183 time_init = 0.0_dp
184 time_end = 0.0_dp
185 time_output = 0.0_dp
186 
187 !-------- Initialisation of the Library of Iterative Solvers Lis,
188 ! if required --------
189 
190 #if (CALCTHK==3 || CALCTHK==6 || MARGIN==3 || DYNAMICS==2)
191  call lis_initialize(ierr)
192 #endif
193 
194 !-------- Read physical parameters --------
195 
196 call read_phys_para()
197 
198 call ice_mat_eqs_pars(rf, r_t, kappa, c, -190, 10, rho_i, rho_c, kappa_c, c_c)
199 
200 ! ------ Density of ice-dust mixture
201 
202 rho = (1.0_dp-frac_dust)*rho_i + frac_dust*rho_c
203 
204 rho_inv = 1.0_dp/rho
205 
206 ! ------ Some auxiliary quantities required for the enthalpy method
207 
208 call calc_c_int_table(c, -190, 10, l)
210 
211 !-------- Compatibility check of the SICOPOLIS version with the header file
212 
213 if ( trim(version) /= trim(sico_version) ) then
214  errormsg = ' >>> sico_init: ' &
215  //'SICOPOLIS version not compatible with header file!'
216  call error(errormsg)
217 end if
218 
219 !-------- Check whether the dynamics and thermodynamics modes are defined
220 
221 #if (!defined(DYNAMICS))
222 errormsg = ' >>> sico_init: DYNAMICS not defined in the header file!'
223 call error(errormsg)
224 #endif
225 
226 #if (!defined(CALCMOD))
227 errormsg = ' >>> sico_init: CALCMOD not defined in the header file!'
228 call error(errormsg)
229 #endif
230 
231 #if (defined(ENTHMOD))
232 errormsg = ' >>> sico_init: ENTHMOD must not be defined any more.' &
233  // end_of_line &
234  //' Please update your header file!'
235 call error(errormsg)
236 #endif
237 
238 !-------- Compatibility check of the horizontal resolution with the
239 ! number of grid points --------
240 
241 #if (GRID==0 || GRID==1)
242 
243 if (approx_equal(dx, 20.0_dp, eps_sp_dp)) then
244 
245  if ((imax /= 90).or.(jmax /= 90)) then
246  errormsg = ' >>> sico_init: Wrong values for IMAX and JMAX!'
247  call error(errormsg)
248  end if
249 
250 else if (approx_equal(dx, 10.0_dp, eps_sp_dp)) then
251 
252  if ((imax /= 180).or.(jmax /= 180)) then
253  errormsg = ' >>> sico_init: Wrong values for IMAX and JMAX!'
254  call error(errormsg)
255  end if
256 
257 else if (approx_equal(dx, 5.0_dp, eps_sp_dp)) then
258 
259  if ((imax /= 360).or.(jmax /= 360)) then
260  errormsg = ' >>> sico_init: Wrong values for IMAX and JMAX!'
261  call error(errormsg)
262  end if
263 
264 else
265 
266  errormsg = ' >>> sico_init: Wrong value for DX!'
267  call error(errormsg)
268 
269 end if
270 
271 #elif (GRID==2)
272 
273 errormsg = ' >>> sico_init: GRID==2 not allowed for nmars application!'
274 call error(errormsg)
275 
276 #endif
277 
278 !-------- Compatibility check of the thermodynamics mode
279 ! (cold vs. polythermal vs. enthalpy method)
280 ! and the number of grid points in the lower (kt) ice domain --------
281 
282 #if (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1)
283 
284 if (ktmax > 2) then
285  write(6, fmt='(a)') ' >>> sico_init: For options CALCMOD==0, 2, 3 or -1,'
286  write(6, fmt='(a)') ' the separate kt domain is redundant.'
287  write(6, fmt='(a)') ' Therefore, consider setting KTMAX to 2.'
288  write(6, fmt='(a)') ' '
289 end if
290 
291 #endif
292 
293 !-------- Compatibility check of discretization schemes for the horizontal and
294 ! vertical advection terms in the temperature and age equations --------
295 
296 #if (ADV_HOR==1)
297 errormsg = ' >>> sico_init: ' &
298  //'Option ADV_HOR==1 (central differences) not defined!'
299 call error(errormsg)
300 #endif
301 
302 !-------- Check whether for the shallow shelf
303 ! or shelfy stream approximation
304 ! the chosen grid is Cartesian coordinates
305 ! without distortion correction (GRID==0) --------
306 
307 #if ((MARGIN==3 && DYNAMICS==1) || DYNAMICS==2) /* requires SSA or SStA */
308 #if (GRID != 0)
309 errormsg = ' >>> sico_init: Distortion correction not implemented' &
310  // end_of_line &
311  //' for the shallow shelf approximation (SSA)' &
312  // end_of_line &
313  //' or the shelfy stream approximation (SStA)' &
314  // end_of_line &
315  //' -> GRID==0 required!'
316 call error(errormsg)
317 #endif
318 #endif
319 
320 !-------- Setting of forcing flag --------
321 
322 forcing_flag = 1 ! forcing by delta_ts
323 
324 !-------- Initialization of numerical time steps --------
325 
326 dtime0 = dtime0
327 dtime_temp0 = dtime_temp0
328 #if (REBOUND==2)
329 dtime_wss0 = dtime_wss0
330 #endif
331 
332 !-------- Further initializations --------
333 
334 dzeta_c = 1.0_dp/real(kcmax,dp)
335 dzeta_t = 1.0_dp/real(ktmax,dp)
336 dzeta_r = 1.0_dp/real(krmax,dp)
337 
338 ndat2d = 1
339 ndat3d = 1
340 
341 !-------- General abbreviations --------
342 
343 ! ------ kc domain
344 
345 if (deform >= eps) then
346 
347  flag_aa_nonzero = .true. ! non-equidistant grid
348 
349  aa = deform
350  ea = exp(aa)
351 
352  kc=0
353  zeta_c(kc) = 0.0_dp
354  eaz_c(kc) = 1.0_dp
355  eaz_c_quotient(kc) = 0.0_dp
356 
357  do kc=1, kcmax-1
358  zeta_c(kc) = kc*dzeta_c
359  eaz_c(kc) = exp(aa*zeta_c(kc))
360  eaz_c_quotient(kc) = (eaz_c(kc)-1.0_dp)/(ea-1.0_dp)
361  end do
362 
363  kc=kcmax
364  zeta_c(kc) = 1.0_dp
365  eaz_c(kc) = exp(aa)
366  eaz_c_quotient(kc) = 1.0_dp
367 
368 else
369 
370  flag_aa_nonzero = .false. ! equidistant grid
371 
372  aa = 0.0_dp
373  ea = 1.0_dp
374 
375  kc=0
376  zeta_c(kc) = 0.0_dp
377  eaz_c(kc) = 1.0_dp
378  eaz_c_quotient(kc) = 0.0_dp
379 
380  do kc=1, kcmax-1
381  zeta_c(kc) = kc*dzeta_c
382  eaz_c(kc) = 1.0_dp
383  eaz_c_quotient(kc) = zeta_c(kc)
384  end do
385 
386  kc=kcmax
387  zeta_c(kc) = 1.0_dp
388  eaz_c(kc) = 1.0_dp
389  eaz_c_quotient(kc) = 1.0_dp
390 
391 end if
392 
393 ! ------ kt domain
394 
395 kt=0
396 zeta_t(kt) = 0.0_dp
397 
398 do kt=1, ktmax-1
399  zeta_t(kt) = kt*dzeta_t
400 end do
401 
402 kt=ktmax
403 zeta_t(kt) = 1.0_dp
404 
405 ! ------ kr domain
406 
407 kr=0
408 zeta_r(kr) = 0.0_dp
409 
410 do kr=1, krmax-1
411  zeta_r(kr) = kr*dzeta_r
412 end do
413 
414 kr=krmax
415 zeta_r(kr) = 1.0_dp
416 
417 !-------- Reshaping of a 2-d array (with indices i, j)
418 ! to a vector (with index n) --------
419 
420 n=1
421 
422 do i=0, imax
423 do j=0, jmax
424  n2i(n) = i
425  n2j(n) = j
426  ij2n(j,i) = n
427  n=n+1
428 end do
429 end do
430 
431 !-------- Specification of current simulation --------
432 
433 runname = runname
434 anfdatname = anfdatname
435 
436 #if (defined(YEAR_ZERO))
438 #else
439 year_zero = 2000.0_dp ! default value 2000 CE
440 #endif
441 
442 time_init0 = time_init0
443 time_end0 = time_end0
444 dtime_ser0 = dtime_ser0
445 
446 #if (CHASM==2)
447 time_chasm_init0 = time_chasm_init0
448 time_chasm_end0 = time_chasm_end0
449 #endif
450 
451 #if (OUTPUT==1 || OUTPUT==3)
452 dtime_out0 = dtime_out0
453 #endif
454 
455 #if (OUTPUT==2 || OUTPUT==3)
456 n_output = n_output
457 time_output0( 1) = time_out0_01
458 time_output0( 2) = time_out0_02
459 time_output0( 3) = time_out0_03
460 time_output0( 4) = time_out0_04
461 time_output0( 5) = time_out0_05
462 time_output0( 6) = time_out0_06
463 time_output0( 7) = time_out0_07
464 time_output0( 8) = time_out0_08
465 time_output0( 9) = time_out0_09
466 time_output0(10) = time_out0_10
467 time_output0(11) = time_out0_11
468 time_output0(12) = time_out0_12
469 time_output0(13) = time_out0_13
470 time_output0(14) = time_out0_14
471 time_output0(15) = time_out0_15
472 time_output0(16) = time_out0_16
473 time_output0(17) = time_out0_17
474 time_output0(18) = time_out0_18
475 time_output0(19) = time_out0_19
476 time_output0(20) = time_out0_20
477 #endif
478 
479 !-------- Write log file --------
480 
481 shell_command = 'if [ ! -d'
482 shell_command = trim(shell_command)//' '//outpath
483 shell_command = trim(shell_command)//' '//'] ; then mkdir'
484 shell_command = trim(shell_command)//' '//outpath
485 shell_command = trim(shell_command)//' '//'; fi'
486 call system(trim(shell_command))
487  ! Check whether directory OUTPATH exists. If not, it is created.
488 
489 filename_with_path = trim(outpath)//'/'//trim(runname)//'.log'
490 
491 open(10, iostat=ios, file=trim(filename_with_path), status='new')
492 
493 if (ios /= 0) then
494  errormsg = ' >>> sico_init: Error when opening the log file!'
495  call error(errormsg)
496 end if
497 
498 write(10, fmt=trim(fmt1)) 'Computational domain:'
499 write(10, fmt=trim(fmt1)) trim(ch_domain_long)
500 write(10, fmt=trim(fmt1)) ' '
501 
502 write(10, fmt=trim(fmt2a)) 'GRID = ', grid
503 write(10, fmt=trim(fmt1)) ' '
504 
505 write(10, fmt=trim(fmt2)) 'imax =', imax
506 write(10, fmt=trim(fmt2)) 'jmax =', jmax
507 write(10, fmt=trim(fmt2)) 'kcmax =', kcmax
508 write(10, fmt=trim(fmt2)) 'ktmax =', ktmax
509 write(10, fmt=trim(fmt2)) 'krmax =', krmax
510 write(10, fmt=trim(fmt1)) ' '
511 
512 write(10, fmt=trim(fmt3)) 'a =', aa
513 write(10, fmt=trim(fmt1)) ' '
514 
515 write(10, fmt=trim(fmt3)) 'x0 =', x0
516 write(10, fmt=trim(fmt3)) 'y0 =', y0
517 write(10, fmt=trim(fmt3)) 'dx =', dx
518 write(10, fmt=trim(fmt1)) ' '
519 
520 write(10, fmt=trim(fmt3)) 'year_zero =', year_zero
521 write(10, fmt=trim(fmt3)) 'time_init =', time_init0
522 write(10, fmt=trim(fmt3)) 'time_end =', time_end0
523 write(10, fmt=trim(fmt3)) 'dtime =', dtime0
524 write(10, fmt=trim(fmt3)) 'dtime_temp =', dtime_temp0
525 #if (REBOUND==2)
526 write(10, fmt=trim(fmt3)) 'dtime_wss =', dtime_wss0
527 #endif
528 write(10, fmt=trim(fmt1)) ' '
529 
530 write(10, fmt=trim(fmt2)) 'ANF_DAT =', anf_dat
531 write(10, fmt=trim(fmt1)) 'zs_present file = '//zs_present_file
532 #if (ANF_DAT==1)
533 write(10, fmt=trim(fmt1)) 'zl_present file = '//zl_present_file
534 #endif
535 write(10, fmt=trim(fmt1)) 'zl0 file = '//zl0_file
536 write(10, fmt=trim(fmt1)) 'mask_present file = '//mask_present_file
537 #if (ANF_DAT==1 && defined(TEMP_INIT))
538 write(10, fmt=trim(fmt2)) 'TEMP_INIT =', temp_init
539 #endif
540 #if (ANF_DAT==3 || (ANF_DAT==1 && TEMP_INIT==5))
541 write(10, fmt=trim(fmt1)) 'Initial-value file = '//anfdatname
542 write(10, fmt=trim(fmt1)) 'Path to initial-value file = '//anfdatpath
543 #endif
544 write(10, fmt=trim(fmt1)) ' '
545 
546 write(10, fmt=trim(fmt1)) 'Physical-parameter file = '//phys_para_file
547 write(10, fmt=trim(fmt1)) ' '
548 
549 #if (defined(THK_EVOL))
550 write(10, fmt=trim(fmt2a)) 'THK_EVOL = ', thk_evol
551 #else
552 errormsg = ' >>> sico_init: Define THK_EVOL in header file!'
553 call error(errormsg)
554 #endif
555 #if (defined(CALCTHK))
556 write(10, fmt=trim(fmt2a)) 'CALCTHK = ', calcthk
557 #else
558 errormsg = ' >>> sico_init: Define CALCTHK in header file!'
559 call error(errormsg)
560 #endif
561 
562 #if (defined(H_ISOL_MAX))
563 write(10, fmt=trim(fmt3)) 'H_isol_max =', h_isol_max
564 #endif
565 
566 #if (CALCTHK==2 || CALCTHK==3 || CALCTHK==5 || CALCTHK==6)
567 write(10, fmt=trim(fmt3)) 'ovi_weight =', ovi_weight
568 #if (CALCTHK==2 || CALCTHK==5)
569 write(10, fmt=trim(fmt3)) 'omega_sor =', omega_sor
570 #if (ITER_MAX_SOR>0)
571 write(10, fmt=trim(fmt2a)) 'iter_max_sor = ', iter_max_sor
572 #endif
573 #endif
574 #endif
575 
576 write(10, fmt=trim(fmt1)) ' '
577 
578 #if (TSURFACE==1 || TSURFACE==6)
579 write(10, fmt=trim(fmt3)) 'delta_ts0 =', delta_ts0
580 #elif (TSURFACE==3)
581 write(10, fmt=trim(fmt3)) 'sine_amplit =', sine_amplit
582 write(10, fmt=trim(fmt3)) 'sine_period =', sine_period
583 #endif
584 #if (TSURFACE==1 || TSURFACE==2 || TSURFACE==3)
585 write(10, fmt=trim(fmt3)) 'temp0_ma_90N =', temp0_ma_90n
586 #endif
587 #if (TSURFACE==1 || TSURFACE==2 || TSURFACE==3 || TSURFACE==4 || TSURFACE==5)
588 write(10, fmt=trim(fmt3)) 'c_ma =', c_ma
589 write(10, fmt=trim(fmt3)) 'gamma_ma =', gamma_ma
590 #endif
591 #if (TSURFACE==5)
592 write(10, fmt=trim(fmt1)) 'Insolation file = '//insol_ma_90n_file
593 #endif
594 #if (TSURFACE==4 || TSURFACE==5 || TSURFACE==6)
595 write(10, fmt=trim(fmt3)) 'albedo =', albedo
596 #endif
597 
598 #if (ACC_PRESENT==1)
599 write(10, fmt=trim(fmt3)) 'acc_present_val =', acc_present_val
600 #endif
601 
602 #if (ACCSURFACE==1 || ACCSURFACE==2)
603 write(10, fmt=trim(fmt3)) 'gamma_s =', gamma_s
604 #endif
605 
606 #if (ABLSURFACE==1 || ABLSURFACE==2)
607 write(10, fmt=trim(fmt3)) 'eld_0 =', eld_0
608 write(10, fmt=trim(fmt3)) 'g_mb_0 =', g_0
609 write(10, fmt=trim(fmt3)) 'gamma_eld =', gamma_eld
610 write(10, fmt=trim(fmt3)) 'gamma_g =', gamma_g
611 #endif
612 
613 write(10, fmt=trim(fmt2a)) 'SEA_LEVEL = ', sea_level
614 #if (SEA_LEVEL==1)
615 write(10, fmt=trim(fmt3)) 'z_sl0 =', z_sl0
616 #elif (SEA_LEVEL==3)
617 write(10, fmt=trim(fmt1)) 'sea-level file = '//sea_level_file
618 #endif
619 write(10, fmt=trim(fmt1)) ' '
620 
621 #if (MARGIN==2)
622 #if ( MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3 )
623 write(10, fmt=trim(fmt3)) 'z_mar =', z_mar
624 write(10, fmt=trim(fmt1)) ' '
625 #elif (MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \
626  || marine_ice_calving==6 || marine_ice_calving==7)
627 write(10, fmt=trim(fmt3)) 'fact_z_mar =', fact_z_mar
628 write(10, fmt=trim(fmt1)) ' '
629 #elif (MARINE_ICE_FORMATION==2 && MARINE_ICE_CALVING==9)
630 write(10, fmt=trim(fmt3)) 'calv_uw_coeff =', calv_uw_coeff
631 write(10, fmt=trim(fmt3)) 'r1_calv_uw =', r1_calv_uw
632 write(10, fmt=trim(fmt3)) 'r2_calv_uw =', r2_calv_uw
633 write(10, fmt=trim(fmt1)) ' '
634 #endif
635 #elif (MARGIN==3)
636 #if (ICE_SHELF_CALVING==2)
637 write(10, fmt=trim(fmt3)) 'H_calv =', h_calv
638 write(10, fmt=trim(fmt1)) ' '
639 #endif
640 #endif
641 
642 #if (defined(BASAL_HYDROLOGY))
643 write(10, fmt=trim(fmt2a)) 'BASAL_HYDROLOGY = ', basal_hydrology
644 #endif
645 write(10, fmt=trim(fmt2a)) 'SLIDE_LAW = ', slide_law
646 write(10, fmt=trim(fmt3)) 'c_slide =', c_slide
647 write(10, fmt=trim(fmt3)) 'gamma_slide =', gamma_slide
648 write(10, fmt=trim(fmt2a)) 'p_weert = ', p_weert
649 write(10, fmt=trim(fmt2a)) 'q_weert = ', q_weert
650 #if (defined(TIME_RAMP_UP_SLIDE))
651 write(10, fmt=trim(fmt3)) 'time_ramp_up_slide =', time_ramp_up_slide
652 #endif
653 #if (SLIDE_LAW==2 || SLIDE_LAW==3)
654 write(10, fmt=trim(fmt3)) 'red_pres_limit_fact =', red_pres_limit_fact
655 #endif
656 #if (BASAL_HYDROLOGY==1 \
657  && defined(hydro_slide_sat_fct) \
658  && defined(c_hw_slide) && defined(hw0_slide))
659 write(10, fmt=trim(fmt2a)) 'HYDRO_SLIDE_SAT_FCT = ', hydro_slide_sat_fct
660 write(10, fmt=trim(fmt3)) 'c_Hw_slide =', c_hw_slide
661 write(10, fmt=trim(fmt3)) 'Hw0_slide =', hw0_slide
662 #endif
663 write(10, fmt=trim(fmt1)) ' '
664 
665 write(10, fmt=trim(fmt2a)) 'Q_GEO_MOD = ', q_geo_mod
666 #if (Q_GEO_MOD==1)
667 write(10, fmt=trim(fmt3)) 'q_geo =', q_geo
668 #elif (Q_GEO_MOD==2)
669 write(10, fmt=trim(fmt1)) 'q_geo file = '//q_geo_file
670 #endif
671 write(10, fmt=trim(fmt1)) ' '
672 
673 write(10, fmt=trim(fmt2)) 'REBOUND =', rebound
674 #if (REBOUND==1)
675 write(10, fmt=trim(fmt3)) 'frac_llra =', frac_llra
676 #endif
677 #if (REBOUND==1 || REBOUND==2)
678 write(10, fmt=trim(fmt2)) 'TIME_LAG_MOD =', time_lag_mod
679 #if (TIME_LAG_MOD==1)
680 write(10, fmt=trim(fmt3)) 'time_lag =', time_lag
681 #elif (TIME_LAG_MOD==2)
682 write(10, fmt=trim(fmt1)) 'time_lag_file = '//time_lag_file
683 #else
684 errormsg = ' >>> sico_init: TIME_LAG_MOD must be either 1 or 2!'
685 call error(errormsg)
686 #endif
687 #endif
688 #if (REBOUND==2)
689 write(10, fmt=trim(fmt2)) 'FLEX_RIG_MOD =', flex_rig_mod
690 #if (FLEX_RIG_MOD==1)
691 write(10, fmt=trim(fmt3)) 'flex_rig =', flex_rig
692 #elif (FLEX_RIG_MOD==2)
693 write(10, fmt=trim(fmt1)) 'flex_rig_file = '//flex_rig_file
694 #else
695 errormsg = ' >>> sico_init: FLEX_RIG_MOD must be either 1 or 2!'
696 call error(errormsg)
697 #endif
698 #endif
699 write(10, fmt=trim(fmt2)) 'Q_LITHO =', q_litho
700 write(10, fmt=trim(fmt1)) ' '
701 
702 #if (FLOW_LAW==2)
703 write(10, fmt=trim(fmt3)) 'gr_size =', gr_size
704 write(10, fmt=trim(fmt1)) ' '
705 #endif
706 #if (FIN_VISC==2)
707 write(10, fmt=trim(fmt3)) 'sigma_res =', sigma_res
708 write(10, fmt=trim(fmt1)) ' '
709 #endif
710 write(10, fmt=trim(fmt3)) 'frac_dust =', frac_dust
711 write(10, fmt=trim(fmt1)) ' '
712 
713 write(10, fmt=trim(fmt2)) 'ENHMOD =', enhmod
714 #if (ENHMOD==1 || ENHMOD==2 || ENHMOD==3)
715 write(10, fmt=trim(fmt3)) 'enh_fact =', enh_fact
716 #endif
717 #if (ENHMOD==2 || ENHMOD==3)
718 write(10, fmt=trim(fmt3)) 'enh_intg =', enh_intg
719 #endif
720 #if (ENHMOD==2)
721 write(10, fmt=trim(fmt3)) 'age_trans =', age_trans_0
722 #endif
723 #if (ENHMOD==3)
724 write(10, fmt=trim(fmt3)) 'date_trans1 =', date_trans1_0
725 write(10, fmt=trim(fmt3)) 'date_trans2 =', date_trans2_0
726 write(10, fmt=trim(fmt3)) 'date_trans3 =', date_trans3_0
727 #endif
728 #if (ENHMOD==4 || ENHMOD==5)
729 write(10, fmt=trim(fmt3)) 'enh_compr =', enh_compr
730 write(10, fmt=trim(fmt3)) 'enh_shear =', enh_shear
731 #endif
732 #if ((ENHMOD==1 || ENHMOD==2 || ENHMOD==3 || ENHMOD==4) && MARGIN==3)
733 write(10, fmt=trim(fmt3)) 'enh_shelf =', enh_shelf
734 #endif
735 write(10, fmt=trim(fmt1)) ' '
736 
737 write(10, fmt=trim(fmt2a)) 'CHASM = ', chasm
738 #if (CHASM==2)
739 write(10, fmt=trim(fmt1)) 'mask_chasm file = '//mask_chasm_file
740 write(10, fmt=trim(fmt3)) 'time_chasm_init =', time_chasm_init0
741 write(10, fmt=trim(fmt3)) 'time_chasm_end =', time_chasm_end0
742 write(10, fmt=trim(fmt3)) 'q_geo_chasm =', q_geo_chasm
743 write(10, fmt=trim(fmt3)) 'erosion_chasm =', erosion_chasm
744 write(10, fmt=trim(fmt1)) ' '
745 #endif
746 
747 write(10, fmt=trim(fmt2a)) 'DYNAMICS = ', dynamics
748 #if ((DYNAMICS==1 && MARGIN==3) || DYNAMICS==2)
749 #if (defined(LIS_OPTS))
750 write(10, fmt=trim(fmt1)) 'lis_opts = '//lis_opts
751 #endif
752 #if (defined(TOL_ITER_SSA))
753 write(10, fmt=trim(fmt3)) 'tol_iter_ssa =', tol_iter_ssa
754 #endif
755 #if (defined(N_ITER_SSA))
756 write(10, fmt=trim(fmt2a)) 'n_iter_ssa = ', n_iter_ssa
757 #endif
758 #if (defined(ITER_INIT_SSA))
759 write(10, fmt=trim(fmt2a)) 'iter_init_ssa = ', iter_init_ssa
760 #endif
761 #if (defined(VISC_INIT_SSA))
762 write(10, fmt=trim(fmt3)) 'visc_init_ssa =', visc_init_ssa
763 #endif
764 #if (defined(N_VISC_SMOOTH))
765 write(10, fmt=trim(fmt2a)) 'n_visc_smooth = ', n_visc_smooth
766 #endif
767 #if (defined(VISC_SMOOTH_DIFF))
768 write(10, fmt=trim(fmt3)) 'visc_smooth_diff =', visc_smooth_diff
769 #endif
770 #if (defined(RELAX_FACT_SSA))
771 write(10, fmt=trim(fmt3)) 'relax_fact_ssa =', relax_fact_ssa
772 #endif
773 #endif
774 #if (DYNAMICS==2 && defined(RATIO_SL_THRESH))
775 write(10, fmt=trim(fmt3)) 'ratio_sl_thresh =', ratio_sl_thresh
776 #endif
777 #if (DYNAMICS==2 && defined(SSTA_SIA_WEIGH_FCT))
778 write(10, fmt=trim(fmt2a)) 'SSTA_SIA_WEIGH_FCT = ', ssta_sia_weigh_fct
779 #endif
780 write(10, fmt=trim(fmt1)) ' '
781 
782 write(10, fmt=trim(fmt2)) 'CALCMOD =', calcmod
783 #if (CALCMOD==-1 && defined(TEMP_CONST))
784 write(10, fmt=trim(fmt3)) 'TEMP_CONST =', temp_const
785 #endif
786 #if (CALCMOD==-1 && defined(AGE_CONST))
787 write(10, fmt=trim(fmt3)) 'AGE_CONST =', age_const
788 #endif
789 #if (CALCMOD==1 && defined(CTS_MELTING_FREEZING))
790 write(10, fmt=trim(fmt2)) 'CTS_MELTING_FREEZING =', cts_melting_freezing
791 #endif
792 write(10, fmt=trim(fmt2)) 'FLOW_LAW =', flow_law
793 write(10, fmt=trim(fmt2)) 'FIN_VISC =', fin_visc
794 write(10, fmt=trim(fmt2)) 'MARGIN =', margin
795 #if (MARGIN==2)
796 write(10, fmt=trim(fmt2)) 'MARINE_ICE_FORMATION =', marine_ice_formation
797 write(10, fmt=trim(fmt2)) 'MARINE_ICE_CALVING =', marine_ice_calving
798 #elif (MARGIN==3)
799 write(10, fmt=trim(fmt2)) 'ICE_SHELF_CALVING =', ice_shelf_calving
800 #endif
801 write(10, fmt=trim(fmt2)) 'ADV_HOR =', adv_hor
802 write(10, fmt=trim(fmt2)) 'ADV_VERT =', adv_vert
803 write(10, fmt=trim(fmt2)) 'TOPOGRAD =', topograd
804 write(10, fmt=trim(fmt2)) 'TSURFACE =', tsurface
805 write(10, fmt=trim(fmt2)) 'ACC_UNIT =', acc_unit
806 write(10, fmt=trim(fmt2)) 'ACC_PRESENT=', acc_present
807 write(10, fmt=trim(fmt2)) 'ACCSURFACE =', accsurface
808 write(10, fmt=trim(fmt2)) 'ABLSURFACE =', ablsurface
809 #if (defined(MB_ACCOUNT))
810 write(10, fmt=trim(fmt2a)) 'MB_ACCOUNT = ', mb_account
811 #endif
812 write(10, fmt=trim(fmt1)) ' '
813 
814 write(10, fmt=trim(fmt3)) 'numdiff_H_t =', numdiff_h_t
815 write(10, fmt=trim(fmt3)) 'tau_cts =', tau_cts
816 write(10, fmt=trim(fmt3)) 'vh_max =', vh_max
817 write(10, fmt=trim(fmt3)) 'hd_min =', hd_min
818 write(10, fmt=trim(fmt3)) 'hd_max =', hd_max
819 #if (defined(VISC_MIN) && defined(VISC_MAX))
820 write(10, fmt=trim(fmt3)) 'visc_min =', visc_min
821 write(10, fmt=trim(fmt3)) 'visc_max =', visc_max
822 #endif
823 write(10, fmt=trim(fmt3)) 'qbm_min =', qbm_min
824 write(10, fmt=trim(fmt3)) 'qbm_max =', qbm_max
825 write(10, fmt=trim(fmt3)) 'age_min =', age_min
826 write(10, fmt=trim(fmt3)) 'age_max =', age_max
827 write(10, fmt=trim(fmt3)) 'mean_accum =', mean_accum
828 #if (ADV_VERT==1)
829 write(10, fmt=trim(fmt3)) 'age_diff =', agediff
830 #endif
831 write(10, fmt=trim(fmt1)) ' '
832 
833 #if (defined(OUT_TIMES))
834 write(10, fmt=trim(fmt2a)) 'OUT_TIMES = ', out_times
835 #endif
836 #if (defined(OUTPUT_INIT))
837 write(10, fmt=trim(fmt2a)) 'OUTPUT_INIT = ', output_init
838 #endif
839 write(10, fmt=trim(fmt2a)) 'OUTPUT = ', output
840 #if (OUTPUT==1 || OUTPUT==3)
841 write(10, fmt=trim(fmt3)) 'dtime_out =' , dtime_out0
842 #endif
843 write(10, fmt=trim(fmt3)) 'dtime_ser =' , dtime_ser0
844 #if (OUTPUT==1 || OUTPUT==2)
845 write(10, fmt=trim(fmt2a)) 'ERGDAT = ', ergdat
846 #endif
847 write(10, fmt=trim(fmt2a)) 'NETCDF = ', netcdf
848 #if (OUTPUT==2 || OUTPUT==3)
849 write(10, fmt=trim(fmt2a)) 'n_output = ', n_output
850 do n=1, n_output
851  write(10, fmt=trim(fmt3)) 'time_output =' , time_output0(n)
852 end do
853 #endif
854 write(10, fmt=trim(fmt1)) ' '
855 
856 write(10, fmt=trim(fmt1)) 'Program version and date: '//version//' / '//date
857 
858 close(10, status='keep')
859 
860 !-------- Conversion of time quantities --------
861 
862 year_zero = year_zero*year_sec ! a --> s
863 time_init = time_init0*year_sec ! a --> s
864 time_end = time_end0*year_sec ! a --> s
865 dtime = dtime0*year_sec ! a --> s
866 dtime_temp = dtime_temp0*year_sec ! a --> s
867 #if (REBOUND==2)
868 dtime_wss = dtime_wss0*year_sec ! a --> s
869 #endif
870 #if (CHASM==2)
871 time_chasm_init = time_chasm_init0 *year_sec ! a --> s
872 time_chasm_end = time_chasm_end0 *year_sec ! a --> s
873 #endif
874 dtime_ser = dtime_ser0*year_sec ! a --> s
875 #if (OUTPUT==1 || OUTPUT==3)
876 dtime_out = dtime_out0*year_sec ! a --> s
877 #endif
878 #if (OUTPUT==2 || OUTPUT==3)
879 do n=1, n_output
880  time_output(n) = time_output0(n)*year_sec ! a --> s
881 end do
882 #endif
883 
884 if (.not.approx_integer_multiple(dtime_temp, dtime, eps_sp_dp)) then
885  errormsg = ' >>> sico_init: dtime_temp must be a multiple of dtime!'
886  call error(errormsg)
887 end if
888 
889 #if (REBOUND==2)
890 if (.not.approx_integer_multiple(dtime_wss, dtime, eps_sp_dp)) then
891  errormsg = ' >>> sico_init: dtime_wss must be a multiple of dtime!'
892  call error(errormsg)
893 end if
894 #endif
895 
896 if (.not.approx_integer_multiple(dtime_ser, dtime, eps_sp_dp)) then
897  errormsg = ' >>> sico_init: dtime_ser must be a multiple of dtime!'
898  call error(errormsg)
899 end if
900 
901 #if (OUTPUT==1 || OUTPUT==3)
902 if (.not.approx_integer_multiple(dtime_out, dtime, eps_sp_dp)) then
903  errormsg = ' >>> sico_init: dtime_out must be a multiple of dtime!'
904  call error(errormsg)
905 end if
906 #endif
907 
908 time = time_init
909 
910 !-------- Read accumulation measurements --------
911 
912 #if (ACC_PRESENT==1)
913 
914 do i=0, imax
915 do j=0, jmax
916  accum_present(j,i) = acc_present_val
917 end do
918 end do
919 
920 #endif
921 
922 ! ------ Conversion mm/a water or (ice+dust) equivalent
923 ! --> m/s (ice+dust) equivalent
924 
925 do i=0, imax
926 do j=0, jmax
927 
928 #if (ACC_UNIT==1)
929 
930  accum_present(j,i) = accum_present(j,i) &
931  *(1.0e-03_dp/year_sec)*(rho_w/rho_i) &
932  *(1.0_dp/(1.0_dp-frac_dust))
933 
934 #elif (ACC_UNIT==2)
935 
936  accum_present(j,i) = accum_present(j,i)*(1.0e-03_dp/year_sec)
937 
938 #endif
939 
940 end do
941 end do
942 
943 !-------- Mean accumulation --------
944 
945 #if (ACC_UNIT==1)
946 
947 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho_i) &
948  *(1.0_dp/(1.0_dp-frac_dust))
949 ! ! mm/a water equiv. --> m/s (ice+dust) equiv.
950 
951 #elif (ACC_UNIT==2)
952 
953 mean_accum = mean_accum*(1.0e-03_dp/year_sec)
954 ! ! mm/a (ice+dust) equiv. --> m/s (ice+dust) equiv.
955 
956 #endif
957 
958 !-------- Read chasm mask --------
959 
960 #if (CHASM==2)
961 
962 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
963  trim(mask_chasm_file)
964 
965 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
966 
967 if (ios /= 0) then
968  errormsg = ' >>> sico_init: Error when opening the chasm mask file!'
969  call error(errormsg)
970 end if
971 
972 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
973 
974 do j=jmax, 0, -1
975  read(24, fmt=trim(fmt4)) (maske_chasm(j,i), i=0,imax)
976 end do
977 
978 close(24, status='keep')
979 
980 #endif
981 
982 !-------- Read data for z_sl --------
983 
984 #if (SEA_LEVEL==3)
985 
986 errormsg = ' >>> sico_init: SEA_LEVEL==3 not allowed for nmars application!'
987 call error(errormsg)
988 
989 #endif
990 
991 !-------- Reading of insolation data --------
992 
993 insol_ma_90 = 0.0_dp ! Assignment of dummy values
994 obl_data = 0.0_dp
995 ecc_data = 0.0_dp
996 ave_data = 0.0_dp
997 cp_data = 0.0_dp
998 
999 #if (TSURFACE==5 || TSURFACE==6)
1000 
1001 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1002  trim(insol_ma_90n_file)
1003 
1004 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1005 
1006 if (ios /= 0) then
1007  errormsg = ' >>> sico_init: Error when opening the insolation-data file!'
1008  call error(errormsg)
1009 end if
1010 
1011 read(21, fmt=*) ch_dummy, insol_time_min, insol_time_stp, insol_time_max
1012 
1013 if (ch_dummy /= '#') then
1014  errormsg = ' >>> sico_init: insol_time_min, insol_time_stp, insol_time_max' &
1015  // end_of_line &
1016  //' not defined in insolation-data file!'
1017  call error(errormsg)
1018 end if
1019 
1021 
1022 if (ndata_insol > 100000) then
1023  errormsg = ' >>> sico_init: Too many data in insolation-data file!'
1024  call error(errormsg)
1025 end if
1026 
1027 do n=0, ndata_insol
1028  read(21, fmt=*) d_dummy, ecc_data(n), obl_data(n), cp_data(n), ave_data(n), insol_ma_90(n)
1029  obl_data(n) = obl_data(n) *pi_180
1030  ave_data(n) = ave_data(n) *pi_180
1031 end do
1032 
1033 close(21, status='keep')
1034 
1035 #endif
1036 
1037 !-------- Determination of the normal geothermal heat flux
1038 ! (without the possible contribution from an active chasm area) --------
1039 
1040 #if (Q_GEO_MOD==1)
1041 
1042 ! ------ Constant value
1043 
1044 do i=0, imax
1045 do j=0, jmax
1046  q_geo_normal(j,i) = q_geo *1.0e-03_dp ! mW/m2 -> W/m2
1047 end do
1048 end do
1049 
1050 #elif (Q_GEO_MOD==2)
1051 
1052 ! ------ Read data from file
1053 
1054 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1055  trim(q_geo_file)
1056 
1057 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1058 
1059 if (ios /= 0) then
1060  errormsg = ' >>> sico_init: Error when opening the qgeo file!'
1061  call error(errormsg)
1062 end if
1063 
1064 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1065 
1066 do j=jmax, 0, -1
1067  read(21, fmt=*) (q_geo_normal(j,i), i=0,imax)
1068 end do
1069 
1070 close(21, status='keep')
1071 
1072 do i=0, imax
1073 do j=0, jmax
1074  q_geo_normal(j,i) = q_geo_normal(j,i) *1.0e-03_dp ! mW/m2 -> W/m2
1075 end do
1076 end do
1077 
1078 #endif
1079 
1080 !-------- Reading of tabulated kei function--------
1081 
1082 #if (REBOUND==0 || REBOUND==1)
1083 
1084 kei = 0.0_dp; n_data_kei = 0; kei_r_max = 0.0_dp; kei_r_incr = 0.0_dp
1085  ! dummy values
1086 #elif (REBOUND==2)
1087 
1088 call read_kei()
1089 
1090 #endif
1091 
1092 !-------- Determination of the time lag
1093 ! of the relaxing asthenosphere --------
1094 
1095 #if (REBOUND==1 || REBOUND==2)
1096 
1097 #if (TIME_LAG_MOD==1)
1098 
1099 time_lag_asth = time_lag*year_sec ! a -> s
1100 
1101 #elif (TIME_LAG_MOD==2)
1102 
1103 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1104  trim(time_lag_file)
1105 
1106 open(29, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1107 
1108 if (ios /= 0) then
1109  errormsg = ' >>> sico_init: Error when opening the time-lag file!'
1110  call error(errormsg)
1111 end if
1112 
1113 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1114 
1115 do j=jmax, 0, -1
1116  read(29, fmt=*) (time_lag_asth(j,i), i=0,imax)
1117 end do
1118 
1119 close(29, status='keep')
1120 
1121 time_lag_asth = time_lag_asth*year_sec ! a -> s
1122 
1123 #endif
1124 
1125 #elif (REBOUND==0)
1126 
1127 time_lag_asth = 0.0_dp ! dummy values
1128 
1129 #endif
1130 
1131 !-------- Determination of the flexural rigidity of the lithosphere --------
1132 
1133 #if (REBOUND==2)
1134 
1135 #if (FLEX_RIG_MOD==1)
1136 
1137 flex_rig_lith = flex_rig
1138 
1139 #elif (FLEX_RIG_MOD==2)
1140 
1141 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1142  trim(flex_rig_file)
1143 
1144 open(29, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1145 
1146 if (ios /= 0) then
1147  errormsg = ' >>> sico_init: Error when opening the flex-rig file!'
1148  call error(errormsg)
1149 end if
1150 
1151 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1152 
1153 do j=jmax, 0, -1
1154  read(29, fmt=*) (flex_rig_lith(j,i), i=0,imax)
1155 end do
1156 
1157 close(29, status='keep')
1158 
1159 #endif
1160 
1161 #elif (REBOUND==0 || REBOUND==1)
1162 
1163 flex_rig_lith = 0.0_dp ! dummy values
1164 
1165 #endif
1166 
1167 !-------- Definition of initial values --------
1168 
1169 ! ------ Present topography
1170 
1171 #if (ANF_DAT==1)
1172 
1173 call topography1(dxi, deta)
1174 
1175 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1176 
1177 call boundary(time_init, dtime, dxi, deta, &
1178  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1179 
1180 where ((maske==0_i1b).or.(maske==3_i1b))
1181  ! grounded or floating ice
1183 elsewhere ! maske==1_i1b or 2_i1b, ice-free land or sea
1184  as_perp_apl = 0.0_dp
1185 end where
1186 
1187 smb_corr = 0.0_dp
1188 
1189 q_bm = 0.0_dp
1190 q_tld = 0.0_dp
1191 q_b_tot = 0.0_dp
1192 
1193 p_b_w = 0.0_dp
1194 h_w = 0.0_dp
1195 
1196 #if (!defined(TEMP_INIT) || TEMP_INIT==1)
1198 #elif (TEMP_INIT==2)
1200 #elif (TEMP_INIT==3)
1202 #elif (TEMP_INIT==4)
1204 #elif (TEMP_INIT==5)
1205  call init_temp_water_age_1_5(z_sl, anfdatname)
1206 #else
1207  errormsg = ' >>> sico_init: TEMP_INIT must be between 1 and 5!'
1208  call error(errormsg)
1209 #endif
1210 
1211 #if (ENHMOD==1)
1212  call calc_enhance_1()
1213 #elif (ENHMOD==2)
1214  call calc_enhance_2()
1215 #elif (ENHMOD==3)
1216  call calc_enhance_3(time_init)
1217 #elif (ENHMOD==4)
1218  call calc_enhance_4()
1219 #elif (ENHMOD==5)
1220  call calc_enhance_5()
1221 #else
1222  errormsg = ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
1223  call error(errormsg)
1224 #endif
1225 
1226 vx_m_sia = 0.0_dp
1227 vy_m_sia = 0.0_dp
1228 vx_m_ssa = 0.0_dp
1229 vy_m_ssa = 0.0_dp
1230 
1231 #if (defined(VISC_INIT_SSA))
1232  vis_ave_g = visc_init_ssa
1233 #else
1234  vis_ave_g = 1.0e+15_dp ! Pa s
1235 #endif
1236 
1237 vis_int_g = 0.0_dp
1238 
1239 ! ------ Ice-free, relaxed bedrock
1240 
1241 #elif (ANF_DAT==2)
1242 
1243 call topography2(dxi, deta)
1244 
1245 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1246 
1247 call boundary(time_init, dtime, dxi, deta, &
1248  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1249 
1250 as_perp_apl = 0.0_dp
1251 
1252 smb_corr = 0.0_dp
1253 
1254 q_bm = 0.0_dp
1255 q_tld = 0.0_dp
1256 q_b_tot = 0.0_dp
1257 
1258 p_b_w = 0.0_dp
1259 h_w = 0.0_dp
1260 
1261 call init_temp_water_age_2()
1262 
1263 #if (ENHMOD==1)
1264  call calc_enhance_1()
1265 #elif (ENHMOD==2)
1266  call calc_enhance_2()
1267 #elif (ENHMOD==3)
1268  call calc_enhance_3(time_init)
1269 #elif (ENHMOD==4)
1270  call calc_enhance_4()
1271 #elif (ENHMOD==5)
1272  call calc_enhance_5()
1273 #else
1274  errormsg = ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
1275  call error(errormsg)
1276 #endif
1277 
1278 vx_m_sia = 0.0_dp
1279 vy_m_sia = 0.0_dp
1280 vx_m_ssa = 0.0_dp
1281 vy_m_ssa = 0.0_dp
1282 
1283 #if (defined(VISC_MIN) && defined(VISC_MAX))
1284  vis_ave_g = visc_max
1285 #else
1286  vis_ave_g = 1.0e+25_dp ! Pa s
1287 #endif
1288 
1289 vis_int_g = 0.0_dp
1290 
1291 ! ------ Read initial state from output of previous simulation
1292 
1293 #elif (ANF_DAT==3)
1294 
1295 call topography3(dxi, deta, z_sl, anfdatname)
1296 
1297 call boundary(time_init, dtime, dxi, deta, &
1298  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1299 
1300 where ((maske==0_i1b).or.(maske==3_i1b))
1301  ! grounded or floating ice
1303 elsewhere ! maske==1_i1b or 2_i1b, ice-free land or sea
1304  as_perp_apl = 0.0_dp
1305 end where
1306 
1307 smb_corr = 0.0_dp
1308 
1309 q_b_tot = q_bm + q_tld
1310 
1311 #endif
1312 
1313 !-------- Inner-point flag --------
1314 
1315 flag_inner_point = .true.
1316 
1317 flag_inner_point(0,:) = .false.
1318 flag_inner_point(jmax,:) = .false.
1319 
1320 flag_inner_point(:,0) = .false.
1321 flag_inner_point(:,imax) = .false.
1322 
1323 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
1324 
1325 #if (GRID==0 || GRID==1)
1326 
1327 do ir=-imax, imax
1328 do jr=-jmax, jmax
1329  dist_dxdy(jr,ir) = sqrt( (real(ir,dp)*dxi)**2 + (real(jr,dp)*deta)**2 )
1330  ! distortion due to stereographic projection not accounted for
1331 end do
1332 end do
1333 
1334 #endif
1335 
1336 !-------- Initial velocities --------
1337 
1338 call calc_temp_melt()
1339 call flag_update_gf_gl_cf()
1340 
1341 #if (DYNAMICS==1 || DYNAMICS==2)
1342 
1343 call calc_vxy_b_sia(time, z_sl)
1344 call calc_vxy_sia(dzeta_c, dzeta_t)
1345 
1346 #if (MARGIN==3 || DYNAMICS==2)
1347 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t)
1348 #endif
1349 
1350 call calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
1351 
1352 #if (MARGIN==3)
1353 call calc_vz_floating(z_sl, dxi, deta, dzeta_c)
1354 #endif
1355 
1356 #elif (DYNAMICS==0)
1357 
1358 call calc_vxy_static()
1359 call calc_vz_static()
1360 
1361 #else
1362 
1363 errormsg = ' >>> sico_init: DYNAMICS must be either 0, 1 or 2!'
1364 call error(errormsg)
1365 
1366 #endif
1367 
1368 call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
1369 
1370 !-------- Initial enthalpies --------
1371 
1372 #if (CALCMOD==0 || CALCMOD==-1)
1373 
1374 do i=0, imax
1375 do j=0, jmax
1376 
1377  do kc=0, kcmax
1378  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1379  end do
1380 
1381  do kt=0, ktmax
1382  enth_t(kt,j,i) = enth_c(0,j,i)
1383  end do
1384 
1385 end do
1386 end do
1387 
1388 #elif (CALCMOD==1)
1389 
1390 do i=0, imax
1391 do j=0, jmax
1392 
1393  do kc=0, kcmax
1394  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1395  end do
1396 
1397  if ( (maske(j,i) == 0_i1b).and.(n_cts(j,i) == 1_i1b) ) then
1398  do kt=0, ktmax
1399  enth_t(kt,j,i) = enth_fct_temp_omega(temp_t_m(kt,j,i), omega_t(kt,j,i))
1400  end do
1401  else
1402  do kt=0, ktmax
1403  enth_t(kt,j,i) = enth_c(0,j,i)
1404  end do
1405  end if
1406 
1407 end do
1408 end do
1409 
1410 #elif (CALCMOD==2 || CALCMOD==3)
1411 
1412 do i=0, imax
1413 do j=0, jmax
1414 
1415  do kc=0, kcmax
1416  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), omega_c(kc,j,i))
1417  end do
1418 
1419  do kt=0, ktmax
1420  enth_t(kt,j,i) = enth_c(0,j,i)
1421  end do
1422 
1423 end do
1424 end do
1425 
1426 #else
1427 
1428 errormsg = ' >>> sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!'
1429 call error(errormsg)
1430 
1431 #endif
1432 
1433 !-------- Initialize time-series files --------
1434 
1435 ! ------ Time-series file for the ice sheet on the whole
1436 
1437 filename_with_path = trim(outpath)//'/'//trim(runname)//'.ser'
1438 
1439 open(12, iostat=ios, file=trim(filename_with_path), status='new')
1440 
1441 if (ios /= 0) then
1442  errormsg = ' >>> sico_init: Error when opening the ser file!'
1443  call error(errormsg)
1444 end if
1445 
1446 write(12,1102)
1447 write(12,1103)
1448 
1449 1102 format(' t(a) D_Ts(C) z_sl(m)',/, &
1450  ' V(m^3) V_g(m^3) V_f(m^3)', &
1451  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1452  ' V_sle(m) V_t(m^3)', &
1453  ' A_t(m^2)',/, &
1454  ' H_max(m) H_t_max(m)', &
1455  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
1456 1103 format('----------------------------------------------------', &
1457  '---------------------------------------')
1458 
1459 ! ------ Time-series file for selected positions ("deep boreholes")
1460 
1461 n_core = 3 ! Points NP (north pole), C1, C2 (in Chasma Borealis)
1462 
1463 allocate(lambda_core(n_core), phi_core(n_core), &
1465 
1466 ch_core(1) = 'North Pole'
1467 lambda_core(1) = 0.0_dp ! dummy
1468 phi_core(1) = 0.0_dp ! dummy
1469 x_core(1) = 0.0_dp *1.0e+03_dp ! North pole,
1470 y_core(1) = 0.0_dp *1.0e+03_dp ! conversion km -> m
1471 
1472 ch_core(2) = 'C1'
1473 lambda_core(2) = 0.0_dp ! dummy
1474 phi_core(2) = 0.0_dp ! dummy
1475 x_core(2) = -150.0_dp *1.0e+03_dp ! Point C1,
1476 y_core(2) = -290.0_dp *1.0e+03_dp ! conversion km -> m
1477 
1478 ch_core(3) = 'C2'
1479 lambda_core(3) = 0.0_dp ! dummy
1480 phi_core(3) = 0.0_dp ! dummy
1481 x_core(3) = -300.0_dp *1.0e+03_dp ! Point C2,
1482 y_core(3) = -280.0_dp *1.0e+03_dp ! conversion km -> m
1483 
1484 filename_with_path = trim(outpath)//'/'//trim(runname)//'.core'
1485 
1486 open(14, iostat=ios, file=trim(filename_with_path), status='new')
1487 
1488 if (ios /= 0) then
1489  errormsg = ' >>> sico_init: Error when opening the core file!'
1490  call error(errormsg)
1491 end if
1492 
1493 if (forcing_flag == 1) then
1494 
1495  write(14,1106)
1496  write(14,1107)
1497 
1498  1106 format(' t(a) D_Ts(C) z_sl(m)',/, &
1499  ' H_NP(m) H_C1(m) H_C2(m)',/, &
1500  ' v_NP(m/a) v_C1(m/a) v_C2(m/a)')
1501  1107 format('----------------------------------------------------')
1502 
1503 end if
1504 
1505 !-------- Output of the initial state --------
1506 
1507 #if (defined(OUTPUT_INIT))
1508 
1509 #if (OUTPUT_INIT==0)
1510  flag_init_output = .false.
1511 #elif (OUTPUT_INIT==1)
1512  flag_init_output = .true.
1513 #else
1514  errormsg = ' >>> sico_init: OUTPUT_INIT must be either 0 or 1!'
1515  call error(errormsg)
1516 #endif
1517 
1518 #else
1519  flag_init_output = .true. ! default for undefined parameter OUTPUT_INIT
1520 #endif
1521 
1522 #if (OUTPUT==1)
1523 
1524 #if (ERGDAT==0)
1525  flag_3d_output = .false.
1526 #elif (ERGDAT==1)
1527  flag_3d_output = .true.
1528 #else
1529  errormsg = ' >>> sico_init: ERGDAT must be either 0 or 1!'
1530  call error(errormsg)
1531 #endif
1532 
1533  if (flag_init_output) &
1534  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1535  flag_3d_output, ndat2d, ndat3d)
1536 
1537 #elif (OUTPUT==2)
1538 
1539 if (time_output(1) <= time_init+eps) then
1540 
1541 #if (ERGDAT==0)
1542  flag_3d_output = .false.
1543 #elif (ERGDAT==1)
1544  flag_3d_output = .true.
1545 #else
1546  errormsg = ' >>> sico_init: ERGDAT must be either 0 or 1!'
1547  call error(errormsg)
1548 #endif
1549 
1550  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1551  flag_3d_output, ndat2d, ndat3d)
1552 
1553 end if
1554 
1555 #elif (OUTPUT==3)
1556 
1557  flag_3d_output = .false.
1558 
1559  if (flag_init_output) &
1560  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1561  flag_3d_output, ndat2d, ndat3d)
1562 
1563 if (time_output(1) <= time_init+eps) then
1564 
1565  flag_3d_output = .true.
1566 
1567  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1568  flag_3d_output, ndat2d, ndat3d)
1569 
1570 end if
1571 
1572 #else
1573 
1574  errormsg = ' >>> sico_init: OUTPUT must be either 1, 2 or 3!'
1575  call error(errormsg)
1576 
1577 #endif
1578 
1579 if (flag_init_output) then
1580  call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1581  call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1582 end if
1583 
1584 end subroutine sico_init
1585 
1586 !-------------------------------------------------------------------------------
1587 !> Definition of the initial surface and bedrock topography
1588 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1589 !! For present-day initial topography.
1590 !<------------------------------------------------------------------------------
1591 subroutine topography1(dxi, deta)
1592 
1593 #if (GRID==0 || GRID==1)
1595 #endif
1596 
1597  use metric_m
1598  use topograd_m
1599 
1600 implicit none
1601 
1602 real(dp), intent(out) :: dxi, deta
1603 
1604 integer(i4b) :: i, j, n
1605 integer(i4b) :: ios, n_dummy
1606 real(dp) :: d_dummy
1607 real(dp) :: xi0, eta0
1608 character :: ch_dummy
1609 
1610 character(len= 8) :: ch_imax
1611 character(len=128) :: fmt4
1612 character(len=256) :: filename_with_path
1613 
1614 write(ch_imax, fmt='(i8)') imax
1615 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
1616 
1617 !-------- Read topography --------
1618 
1619 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1620  trim(zs_present_file)
1621 
1622 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1623 
1624 if (ios /= 0) then
1625  errormsg = ' >>> topography1: Error when opening the zs file!'
1626  call error(errormsg)
1627 end if
1628 
1629 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1630 
1631 do j=jmax, 0, -1
1632  read(21, fmt=*) (zs(j,i), i=0,imax)
1633 end do
1634 
1635 close(21, status='keep')
1636 
1637 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1638  trim(zl_present_file)
1639 
1640 open(22, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1641 
1642 if (ios /= 0) then
1643  errormsg = ' >>> topography1: Error when opening the zl file!'
1644  call error(errormsg)
1645 end if
1646 
1647 do n=1, 6; read(22, fmt='(a)') ch_dummy; end do
1648 
1649 do j=jmax, 0, -1
1650  read(22, fmt=*) (zl(j,i), i=0,imax)
1651 end do
1652 
1653 close(22, status='keep')
1654 
1655 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1656  trim(zl0_file)
1657 
1658 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1659 
1660 if (ios /= 0) then
1661  errormsg = ' >>> topography1: Error when opening the zl0 file!'
1662  call error(errormsg)
1663 end if
1664 
1665 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
1666 
1667 do j=jmax, 0, -1
1668  read(23, fmt=*) (zl0(j,i), i=0,imax)
1669 end do
1670 
1671 close(23, status='keep')
1672 
1673 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1674  trim(mask_present_file)
1675 
1676 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
1677 
1678 if (ios /= 0) then
1679  errormsg = ' >>> topography1: Error when opening the mask file!'
1680  call error(errormsg)
1681 end if
1682 
1683 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1684 
1685 do j=jmax, 0, -1
1686  read(24, fmt=trim(fmt4)) (maske(j,i), i=0,imax)
1687 end do
1688 
1689 close(24, status='keep')
1690 
1691 !-------- Further stuff --------
1692 
1693 dxi = dx *1000.0_dp ! km -> m
1694 deta = dx *1000.0_dp ! km -> m
1695 
1696 xi0 = x0 *1000.0_dp ! km -> m
1697 eta0 = y0 *1000.0_dp ! km -> m
1698 
1699 do i=0, imax
1700 do j=0, jmax
1701 
1702  if (maske(j,i) <= 1_i1b) then
1703  zb(j,i) = zl(j,i)
1704  else ! (maske(j,i)>=2_i1b)
1705  errormsg = ' >>> topography1: ' &
1706  //'maske(j,i)>=2 not allowed for initial topography!'
1707  call error(errormsg)
1708  end if
1709 
1710  zs(j,i) = zs(j,i) *1000.0_dp
1711  zb(j,i) = zb(j,i) *1000.0_dp ! km --> m
1712  zl(j,i) = zl(j,i) *1000.0_dp
1713  zl0(j,i) = zl0(j,i)*1000.0_dp
1714 
1715  xi(i) = xi0 + real(i,dp)*dxi
1716  eta(j) = eta0 + real(j,dp)*deta
1717 
1718  zm(j,i) = zb(j,i)
1719  n_cts(j,i) = -1_i1b
1720  kc_cts(j,i) = 0
1721 
1722  h_c(j,i) = zs(j,i)-zm(j,i)
1723  h_t(j,i) = 0.0_dp
1724 
1725  dzs_dtau(j,i) = 0.0_dp
1726  dzm_dtau(j,i) = 0.0_dp
1727  dzb_dtau(j,i) = 0.0_dp
1728  dzl_dtau(j,i) = 0.0_dp
1729  dh_c_dtau(j,i) = 0.0_dp
1730  dh_t_dtau(j,i) = 0.0_dp
1731 
1732 end do
1733 end do
1734 
1735 maske_old = maske
1736 
1737 !-------- Geographic coordinates, metric tensor,
1738 ! gradients of the topography --------
1739 
1740 do i=0, imax
1741 do j=0, jmax
1742 
1743 #if (GRID==0 || GRID==1) /* Stereographic projection */
1744 
1745  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1746  lambda0, phi0, lambda(j,i), phi(j,i))
1747 
1748 #elif (GRID==2) /* Geographic coordinates */
1749 
1750  lambda(j,i) = xi(i)
1751  phi(j,i) = eta(j)
1752 
1753 #endif
1754 
1755 end do
1756 end do
1757 
1758 call metric()
1759 
1760 #if (TOPOGRAD==0)
1761 call topograd_1(dxi, deta, 1)
1762 #elif (TOPOGRAD==1)
1763 call topograd_2(dxi, deta, 1)
1764 #endif
1765 
1766 !-------- Corresponding area of grid points --------
1767 
1768 do i=0, imax
1769 do j=0, jmax
1770  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1771 end do
1772 end do
1773 
1774 end subroutine topography1
1775 
1776 !-------------------------------------------------------------------------------
1777 !> Definition of the initial surface and bedrock topography
1778 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1779 !! For ice-free initial topography with relaxed lithosphere.
1780 !<------------------------------------------------------------------------------
1781 subroutine topography2(dxi, deta)
1782 
1783 #if (GRID==0 || GRID==1)
1785 #endif
1786 
1787  use metric_m
1788  use topograd_m
1789 
1790 implicit none
1791 
1792 real(dp), intent(out) :: dxi, deta
1793 
1794 integer(i4b) :: i, j, n
1795 integer(i4b) :: ios
1796 real(dp) :: xi0, eta0
1797 character :: ch_dummy
1798 
1799 character(len= 8) :: ch_imax
1800 character(len=128) :: fmt4
1801 character(len=256) :: filename_with_path
1802 
1803 write(ch_imax, fmt='(i8)') imax
1804 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
1805 
1806 !-------- Read topography --------
1807 
1808 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1809  trim(zl0_file)
1810 
1811 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1812 
1813 if (ios /= 0) then
1814  errormsg = ' >>> topography2: Error when opening the zl0 file!'
1815  call error(errormsg)
1816 end if
1817 
1818 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
1819 
1820 do j=jmax, 0, -1
1821  read(23, fmt=*) (zl0(j,i), i=0,imax)
1822 end do
1823 
1824 close(23, status='keep')
1825 
1826 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1827  trim(mask_present_file)
1828 
1829 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
1830 
1831 if (ios /= 0) then
1832  errormsg = ' >>> topography2: Error when opening the mask file!'
1833  call error(errormsg)
1834 end if
1835 
1836 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1837 
1838 do j=jmax, 0, -1
1839  read(24, fmt=trim(fmt4)) (maske(j,i), i=0,imax)
1840 end do
1841 
1842 close(24, status='keep')
1843 
1844 !-------- Further stuff --------
1845 
1846 dxi = dx *1000.0_dp ! km -> m
1847 deta = dx *1000.0_dp ! km -> m
1848 
1849 xi0 = x0 *1000.0_dp ! km -> m
1850 eta0 = y0 *1000.0_dp ! km -> m
1851 
1852 do i=0, imax
1853 do j=0, jmax
1854 
1855  if (maske(j,i) <= 1_i1b) then
1856  maske(j,i) = 1_i1b
1857  zs(j,i) = zl0(j,i)
1858  zb(j,i) = zl0(j,i)
1859  zl(j,i) = zl0(j,i)
1860  else ! (maske(j,i)>=2_i1b)
1861  errormsg = ' >>> topography2: ' &
1862  //'maske(j,i)>=2 not allowed for initial topography!'
1863  call error(errormsg)
1864  end if
1865 
1866  zs(j,i) = zs(j,i) *1000.0_dp
1867  zb(j,i) = zb(j,i) *1000.0_dp ! km --> m
1868  zl(j,i) = zl(j,i) *1000.0_dp
1869  zl0(j,i) = zl0(j,i)*1000.0_dp
1870 
1871  xi(i) = xi0 + real(i,dp)*dxi
1872  eta(j) = eta0 + real(j,dp)*deta
1873 
1874  zm(j,i) = zb(j,i)
1875  n_cts(j,i) = -1_i1b
1876  kc_cts(j,i) = 0
1877 
1878  h_c(j,i) = 0.0_dp
1879  h_t(j,i) = 0.0_dp
1880 
1881  dzs_dtau(j,i) = 0.0_dp
1882  dzm_dtau(j,i) = 0.0_dp
1883  dzb_dtau(j,i) = 0.0_dp
1884  dzl_dtau(j,i) = 0.0_dp
1885  dh_c_dtau(j,i) = 0.0_dp
1886  dh_t_dtau(j,i) = 0.0_dp
1887 
1888 end do
1889 end do
1890 
1891 maske_old = maske
1892 
1893 !-------- Geographic coordinates, metric tensor,
1894 ! gradients of the topography --------
1895 
1896 do i=0, imax
1897 do j=0, jmax
1898 
1899 #if (GRID==0 || GRID==1) /* Stereographic projection */
1900 
1901  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1902  lambda0, phi0, lambda(j,i), phi(j,i))
1903 
1904 #elif (GRID==2) /* Geographic coordinates */
1905 
1906  lambda(j,i) = xi(i)
1907  phi(j,i) = eta(j)
1908 
1909 #endif
1910 
1911 end do
1912 end do
1913 
1914 call metric()
1915 
1916 #if (TOPOGRAD==0)
1917 call topograd_1(dxi, deta, 1)
1918 #elif (TOPOGRAD==1)
1919 call topograd_2(dxi, deta, 1)
1920 #endif
1921 
1922 !-------- Corresponding area of grid points --------
1923 
1924 do i=0, imax
1925 do j=0, jmax
1926  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1927 end do
1928 end do
1929 
1930 end subroutine topography2
1931 
1932 !-------------------------------------------------------------------------------
1933 !> Definition of the initial surface and bedrock topography
1934 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1935 !! For initial topography from previous simulation.
1936 !<------------------------------------------------------------------------------
1937 subroutine topography3(dxi, deta, z_sl, anfdatname)
1938 
1939  use read_m, only : read_erg_nc
1940 
1941 #if (GRID==0 || GRID==1)
1943 #endif
1944 
1945  use metric_m
1946  use topograd_m
1947 
1948 implicit none
1949 
1950 character(len=100), intent(in) :: anfdatname
1951 
1952 real(dp), intent(out) :: dxi, deta, z_sl
1953 
1954 integer(i4b) :: i, j, n
1955 integer(i4b) :: ios
1956 character(len=256) :: filename_with_path
1957 character :: ch_dummy
1958 
1959 !-------- Read data from time-slice file of previous simulation --------
1960 
1961 call read_erg_nc(z_sl, anfdatname)
1962 
1963 !-------- Read topography of the relaxed bedrock --------
1964 
1965 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1966  trim(zl0_file)
1967 
1968 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1969 
1970 if (ios /= 0) then
1971  errormsg = ' >>> topography3: Error when opening the zl0 file!'
1972  call error(errormsg)
1973 end if
1974 
1975 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
1976 
1977 do j=jmax, 0, -1
1978  read(23, fmt=*) (zl0(j,i), i=0,imax)
1979 end do
1980 
1981 close(23, status='keep')
1982 
1983 !-------- Further stuff --------
1984 
1985 dxi = dx *1000.0_dp ! km -> m
1986 deta = dx *1000.0_dp ! km -> m
1987 
1988 zl0 = zl0 *1000.0_dp ! km -> m
1989 
1990 !-------- Geographic coordinates, metric tensor,
1991 ! gradients of the topography --------
1992 
1993 do i=0, imax
1994 do j=0, jmax
1995 
1996 #if (GRID==0 || GRID==1) /* Stereographic projection */
1997 
1998  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1999  lambda0, phi0, lambda(j,i), phi(j,i))
2000 
2001 #elif (GRID==2) /* Geographic coordinates */
2002 
2003  lambda(j,i) = xi(i)
2004  phi(j,i) = eta(j)
2005 
2006 #endif
2007 
2008 end do
2009 end do
2010 
2011 call metric()
2012 
2013 #if (TOPOGRAD==0)
2014 call topograd_1(dxi, deta, 1)
2015 #elif (TOPOGRAD==1)
2016 call topograd_2(dxi, deta, 1)
2017 #endif
2018 
2019 !-------- Corresponding area of grid points --------
2020 
2021 do i=0, imax
2022 do j=0, jmax
2023  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
2024 end do
2025 end do
2026 
2027 end subroutine topography3
2028 
2029 !-------------------------------------------------------------------------------
2030 
2031 end module sico_init_m
2032 !
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c
temp_c(kc,j,i): Temperature in the upper (kc) ice domain
real(dp), dimension(0:jmax, 0:imax) c_slide
c_slide(j,i): Basal sliding coefficient
integer(i1b), dimension(0:jmax, 0:imax) n_cts
n_cts(j,i): Mask for thermal conditions. -1: cold ice base, 0: temperate ice base with cold ice above...
integer(i4b), dimension((imax+1)*(jmax+1)) n2j
n2j: Conversion from linear index n to 2d index j
Update of the flags for the land-terminating grounded front, marine-terminating grounded front...
real(dp), dimension(0:jmax, 0:imax) phi
phi(j,i): Geographic latitude of grid point (i,j)
real(dp), dimension(:), allocatable phi_core
phi_core(n): Geographical latitude of the prescribed borehole positions
subroutine, public read_kei()
Reading of the tabulated kei function.
Definition: read_m.F90:1061
logical, dimension(0:jmax, 0:imax) flag_inner_point
flag_inner_point(j,i): Inner-point flag. .true.: inner point, .false.: margin point ...
subroutine, public calc_c_int_inv_table()
Computation of the inverse of the temperature integral of the specific heat of ice as a table (c_int_...
real(dp), dimension(0:kcmax) eaz_c
eaz_c(kc): Abbreviation for exp(aa*zeta(kc))
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
subroutine, public flag_update_gf_gl_cf()
Main subroutine of flag_update_gf_gl_cf_m: Update of the flags for the land-terminating grounded fron...
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.
subroutine, public calc_enhance_2()
Computation of the flow enhancement factor. Case ENHMOD==2: two different values depending on age for...
real(dp) a
A: Semi-major axis of the planet.
subroutine, public calc_vz_floating(z_sl, dxi, deta, dzeta_c)
Computation of the vertical velocity vz for floating ice.
Definition: calc_vz_m.F90:259
subroutine error(error_message)
Main routine of error_m: Writing of error messages and stopping execution.
Definition: error_m.F90:47
subroutine, public output1(runname, time, delta_ts, glac_index, z_sl, flag_3d_output, ndat2d, ndat3d)
Writing of time-slice files in native binary or NetCDF format.
Definition: output_m.F90:60
Material properties of ice: Rate factor, heat conductivity, specific heat (heat capacity), creep function, viscosity.
real(dp), parameter eps
eps: Small number
subroutine, public read_phys_para()
Reading of physical parameters.
Definition: read_m.F90:1119
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) vis_ave_g
vis_ave_g(j,i): Depth-averaged viscosity of the SIA/SSA, at (i,j)
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) enth_t
enth_t(kt,j,i): Enthalpy in the lower (kt) ice domain
Conversion from temperature (temp) and water content (omega) to enthalpy (enth) and vice versa...
subroutine, public init_temp_water_age_2()
Initial temperature, water content and age (case ANF_DAT==2: ice-free conditions with relaxed bedrock...
real(dp), dimension(0:jmax, 0:imax) dh_t_dtau
dH_t_dtau(j,i): Derivative of H_t by tau (time)
real(dp), dimension(0:jmax, 0:imax) lambda
lambda(j,i): Geographic longitude of grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) dh_c_dtau
dH_c_dtau(j,i): Derivative of H_c by tau (time)
subroutine, public read_erg_nc(z_sl, filename, opt_maske, opt_n_cts, opt_kc_cts, opt_H_cold, opt_H_temp, opt_H, opt_temp_r, opt_omega_t, opt_age_t, opt_temp_c, opt_age_c, opt_omega_c, opt_flag_temp_age_only)
Reading of time-slice files in native binary or NetCDF format.
Definition: read_m.F90:60
real(dp), dimension(0:jmax, 0:imax) p_b_w
p_b_w(j,i): Basal water pressure
integer(i1b), dimension(0:jmax, 0:imax) maske_old
maske_old(j,i): Old value of maske (at the previous time step)
real(dp), dimension(0:jmax, 0:imax) vx_m_ssa
vx_m_ssa(j,i): Mean (depth-averaged) SSA velocity in x-direction, at (i+1/2,j)
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) area
area(j,i): Area of grid cell associated with grid point (i,j)
Computation of all components of the strain-rate tensor, the full effective strain rate and the shear...
Definition: calc_dxyz_m.F90:37
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) temp_t_m
temp_t_m(kt,j,i): Melting temperature in the lower (kt) ice domain
subroutine, public stereo_inv_ellipsoid(x_val, y_val, A, B, lambda0, phi0, lambda_val, phi_val)
Inverse stereographic projection for an ellipsoidal planet.
subroutine topograd_1(dxi, deta, n_switch)
Calculation of topography gradients on the staggered grid and on the grid points (the latter by secon...
Definition: topograd_m.F90:57
subroutine, public calc_enhance_5()
Computation of the flow enhancement factor. Case ENHMOD==5: minimal anisotropic enhancement factor fo...
real(dp) b
B: Semi-minor axis of the planet.
subroutine, public init_temp_water_age_1_2()
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==2: present-day initial topogr...
real(dp), dimension(0:jmax, 0:imax) zl
zl(j,i): Coordinate z of the lithosphere surface
subroutine sico_init(delta_ts, glac_index, mean_accum, dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, time, time_init, time_end, time_output, dxi, deta, dzeta_c, dzeta_t, dzeta_r, z_sl, dzsl_dtau, z_mar, ndat2d, ndat3d, n_output, runname)
Main routine of sico_init_m: Initialisations for SICOPOLIS.
Definition: sico_init_m.F90:59
subroutine, public calc_vxy_static()
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy etc. for static ice.
subroutine, public output4(time, dxi, deta, delta_ts, glac_index, z_sl)
Writing of time-series data of the deep ice cores on file in ASCII format (and optionally in NetCDF f...
Definition: output_m.F90:5384
real(dp), dimension(0:jmax, 0:imax) q_tld
Q_tld(j,i): Water drainage rate from the temperate layer.
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
integer(i4b), parameter rcl1
rcl1: Maximum length of record for input files (with factor 3 safety margin)
real(dp), dimension(0:jmax, 0:imax) vis_int_g
vis_int_g(j,i): Depth-integrated viscosity of the SIA/SSA, at (i,j)
Initialisations for SICOPOLIS.
Definition: sico_init_m.F90:35
Computation of the flow enhancement factor.
Computation of the horizontal velocity vx, vy.
Definition: calc_vxy_m.F90:36
integer(i4b), dimension(0:jmax, 0:imax) q_weert
q_weert(j,i): Weertman exponent for the basal pressure
real(dp), dimension(0:jmax, 0:imax) vy_m_ssa
vy_m_ssa(j,i): Mean (depth-averaged) SSA velocity in y-direction, at (i,j+1/2)
subroutine, public calc_temp_melt()
Computation of the melting temperatures.
real(dp) ea
ea: Abbreviation for exp(aa)
real(dp), parameter pi_180
pi_180: pi divided by 180 (-> deg to rad)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) enth_c
enth_c(kc,j,i): Enthalpy in the upper (kc) ice domain
logical flag_aa_nonzero
flag_aa_nonzero: Flag for the exponential stretch parameter aa. .true.: aa greater than zero (non-equ...
Initial temperature, water content and age.
subroutine, public calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
Computation of the vertical velocity vz for grounded ice.
Definition: calc_vz_m.F90:52
real(dp) aa
aa: Exponential stretch parameter of the non-equidistant vertical grid in the upper (kc) ice domain ...
Reading of several input data.
Definition: read_m.F90:35
Declarations of kind types for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) zs
zs(j,i): Coordinate z of the surface topography
integer(i4b), dimension(0:jmax, 0:imax) kc_cts
kc_cts(j,i): Position kc of the CTS (for COLD, ENTC, ENTM)
real(dp) rho_c
RHO_C: Density of crustal material (dust)
Definition: sico_vars_m.F90:79
real(dp), dimension(0:jmax, 0:imax) h_c
H_c(j,i): Thickness of ice in the upper (kc) domain (thickness of the cold-ice layer for POLY...
real(dp), dimension(0:jmax, 0:imax) smb_corr
smb_corr(j,i): Diagnosed SMB correction
real(dp) time_chasm_end
time_chasm_end: Final time for active chasma area
Definition: sico_vars_m.F90:72
real(dp) kei_r_max
kei_r_max: Maximum value of the argument r of the tabulated kei function
real(dp) year_zero
year_zero: SICOPOLIS year zero in astronomical year numbering [ = signed year CE (AD) ] ...
Computation of the melting and basal temperatures.
subroutine topograd_2(dxi, deta, n_switch)
Calculation of topography gradients on the staggered grid and on the grid points (the latter by fourt...
Definition: topograd_m.F90:203
real(dp), dimension(-190:10) kappa
KAPPA(n): Tabulated values for the heat conductivity of ice.
subroutine, public calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t)
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy and the flux across ...
Computation of the vertical velocity vz.
Definition: calc_vz_m.F90:35
Writing of error messages and stopping execution.
Definition: error_m.F90:35
real(dp), dimension(0:jmax, 0:imax) dzm_dtau
dzm_dtau(j,i): Derivative of zm by tau (time)
real(dp), dimension(0:krmax) zeta_r
zeta_r(kr): Sigma coordinate zeta_r of grid point kr
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) omega_c
omega_c(kc,j,i): Water content in the upper (kc) ice domain
real(dp), dimension(:), allocatable y_core
y_core(n): Coordinate eta (= y) of the prescribed borehole positions
character(len=64) ch_domain_long
ch_domain_long: Long name of the computational domain
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) kappa_c
KAPPA_C: Heat conductivity of crustal material (dust)
Definition: sico_vars_m.F90:81
real(dp) rho_i
RHO_I: Density of ice.
Definition: sico_vars_m.F90:77
Definition of the components g11 and g22 of the metric tensor of the applied numerical coordinates...
Definition: metric_m.F90:37
real(dp) c_c
C_C: Specific heat of crustal material (dust)
Definition: sico_vars_m.F90:83
Computation of the forward or inverse stereographic projection, alternatively for a spherical or an e...
real(dp) rho_inv
rho_inv: Inverse of the density of ice-dust mixture
Definition: sico_vars_m.F90:85
real(dp), dimension(0:jmax, 0:imax) zl0
zl0(j,i): zl for isostatically relaxed ice-free conditions
subroutine, public calc_vz_static()
Computation of the vertical velocity vz for static ice.
Definition: calc_vz_m.F90:374
real(dp) time_chasm_init
time_chasm_init: Initial time for active chasma area
Definition: sico_vars_m.F90:70
real(dp) phi0
PHI0: Standard parallel of the stereographic projection.
character, parameter end_of_line
end_of_line: End-of-line string
subroutine, public ice_mat_eqs_pars(RF_table, R_T_val, KAPPA_table, C_table, n_tmp_min, n_tmp_max, RHO_I_val, RHO_C_val, KAPPA_C_val, C_C_val)
Setting of required physical parameters.
Comparison of single- or double-precision floating-point numbers.
real(dp), dimension(0:jmax, 0:imax) vy_m_sia
vy_m_sia(j,i): Mean (depth-averaged) SIA velocity in y-direction, at (i,j+1/2)
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
subroutine, public calc_enhance_1()
Computation of the flow enhancement factor. Case ENHMOD==1: constant for grounded ice...
real(dp), dimension(0:jmax, 0:imax) zm
zm(j,i): Coordinate z of the bottom of the upper (kc) ice domain = top of the lower (kt) ice domain (...
real(dp), dimension(-190:10) rf
RF(n): Tabulated values for the rate factor of cold ice.
integer(i1b), 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 ...
real(dp), dimension(0:jmax, 0:imax) q_b_tot
Q_b_tot(j,i): Sum of Q_bm and Q_tld.
real(dp), dimension(0:jmax, 0:imax) q_bm
Q_bm(j,i): Basal melting rate.
real(dp) function, public enth_fct_temp_omega(temp_val, omega_val)
Enthalpy as a function of temperature and water content.
subroutine, public calc_vxy_sia(dzeta_c, dzeta_t)
Computation of the shear stresses txz, tyz, the effective shear stress sigma, the depth-averaged flui...
Definition: calc_vxy_m.F90:573
real(dp), dimension(-190:10) c
C(n): Tabulated values for the specific heat of ice.
real(dp), dimension(0:kcmax) eaz_c_quotient
eaz_c_quotient(kc): Abbreviation for (eaz_c(kc)-1.0)/(ea-1.0)
real(dp), dimension(-10000:10000) kei
kei(n): Tabulated values of the kei function (Kelvin function of zero order)
subroutine, public init_temp_water_age_1_3()
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==3: present-day initial topogr...
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:61
real(dp), dimension(0:100000) obl_data
obl_data(n): Data values for the obliquity
Definition: sico_vars_m.F90:52
integer(i4b), parameter rcl2
rcl2: Maximum length of record for input mask files (with factor 3 safety margin) ...
integer(i4b), dimension(0:jmax, 0:imax) ij2n
ij2n: Conversion from 2d index pair (i,j) to linear index n
subroutine, public output2(time, dxi, deta, delta_ts, glac_index, z_sl)
Writing of time-series data on file in ASCII format (and optionally in NetCDF format).
Definition: output_m.F90:4039
integer(i1b) forcing_flag
forcing_flag: Flag for the forcing type. 1: forcing by a spatially constant surface temperature anoma...
real(dp), dimension(-jmax:jmax,-imax:imax) dist_dxdy
dist_dxdy(jr,ir): Distance between grid points with delta_i=ir, delta_j=jr
Calculation of topography gradients on the staggered grid and on the grid points (including length re...
Definition: topograd_m.F90:39
character(len=16), dimension(:), allocatable ch_core
ch_core(n): Names of the prescribed borehole positions
real(dp), dimension(0:jmax, 0:imax) sq_g11_g
sq_g11_g(j,i): Square root of the coefficient g11 of the metric tensor on grid point (i...
real(dp), dimension(0:jmax, 0:imax) as_perp_apl
as_perp_apl(j,i): Applied accumulation-ablation function (SMB)
subroutine, public metric()
Main routine of module metric_m: Definition of the components g11 and g22 of the metric tensor of the...
Definition: metric_m.F90:54
integer(i1b), 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
subroutine, public calc_enhance_3(time)
Computation of the flow enhancement factor. Case ENHMOD==3: two different values depending on time of...
real(dp), dimension(0:jmax, 0:imax) time_lag_asth
time_lag_asth(j,i): Time lag of the relaxing asthenosphere
real(dp), dimension(0:jmax) eta
eta(j): Coordinate eta (= y) of grid point j
real(dp), dimension(:), allocatable x_core
x_core(n): Coordinate xi (= x) of the prescribed borehole positions
integer(i4b), dimension((imax+1)*(jmax+1)) n2i
n2i: Conversion from linear index n to 2d index i
real(dp), dimension(0:imax) xi
xi(i): Coordinate xi (= x) of grid point i
character(len=16) ch_domain_short
ch_domain_short: Short name of the computational domain
real(dp) r_t
R_T: Coefficient of the water-content dependence in the rate factor for temperate ice...
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) omega_t
omega_t(kt,j,i): Water content in the lower (kt) ice domain
subroutine, public calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
Main subroutine of calc_dxyz_m: Computation of all components of the strain-rate tensor, the full effective strain rate and the shear fraction.
Definition: calc_dxyz_m.F90:57
subroutine, public calc_enhance_4()
Computation of the flow enhancement factor. Case ENHMOD==4: minimal anisotropic enhancement factor fo...
real(dp) kei_r_incr
kei_r_incr: Increment of the argument r of the tabulated kei function
real(dp), dimension(0:ktmax) zeta_t
zeta_t(kt): Sigma coordinate zeta_t of grid point kt
subroutine, public init_temp_water_age_1_5(z_sl, filename)
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==5: present-day initial topogr...
Computation of the surface temperature (must be less than 0 deg C!) and of the accumulation-ablation ...
Definition: boundary_m.F90:37
real(dp) l
L: Latent heat of ice.
real(dp), parameter eps_sp_dp
eps_sp_dp: Small number to single-precision accuracy in double precision
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
integer(i4b) n_data_kei
n_data_kei: Number of tabulated values of the kei function
character(len=256) errormsg
errormsg: Error-message string
real(dp) rho
RHO: Density of ice.
subroutine topography3(dxi, deta, z_sl, anfdatname)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
integer(i4b), dimension(0:jmax, 0:imax) p_weert
p_weert(j,i): Weertman exponent for the basal shear stress
subroutine, public init_temp_water_age_1_4()
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==4: present-day initial topogr...
subroutine, public init_temp_water_age_1_1()
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==1: present-day initial topogr...
real(dp), dimension(:), allocatable lambda_core
lambda_core(n): Geographical longitude of the prescribed borehole positions
integer(i4b) n_core
n_core: Number of positions to be considered in the time-series file for deep boreholes ...
integer, parameter dp
Double-precision reals.
subroutine topography2(dxi, deta)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
real(dp), dimension(0:jmax, 0:imax) zb
zb(j,i): Coordinate z of the ice base
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
real(dp), dimension(0:jmax, 0:imax) vx_m_sia
vx_m_sia(j,i): Mean (depth-averaged) SIA velocity in x-direction, at (i+1/2,j)
real(dp), dimension(0:jmax, 0:imax) dzb_dtau
dzb_dtau(j,i): Derivative of zb by tau (time)
real(dp), dimension(0:kcmax) zeta_c
zeta_c(kc): Sigma coordinate zeta_c of grid point kc
real(dp) lambda0
LAMBDA0: Reference longitude (central meridian) of the stereographic projection.
subroutine, public calc_c_int_table(c_table, n_tmp_min, n_tmp_max, L_val)
Computation of the temperature integral of the specific heat of ice as a table (c_int_table). Further, definition of the latent heat of ice.
real(dp), dimension(0:jmax, 0:imax) h_t
H_t(j,i): Thickness of ice in the lower (kt) domain (thickness of the temperate layer for POLY...
subroutine topography1(dxi, deta)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
real(dp), dimension(0:jmax, 0:imax) h_w
H_w(j,i): Thickness of the water column under the ice base.
subroutine, public calc_vxy_b_sia(time, z_sl)
Computation of the basal horizontal velocity vx_b, vy_b in the shallow ice approximation.
Definition: calc_vxy_m.F90:55
Declarations of global variables for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) dzs_dtau
dzs_dtau(j,i): Derivative of zs by tau (time)
real(dp), dimension(0:jmax, 0:imax) flex_rig_lith
flex_rig_lith(j,i): Flexural rigidity of the lithosphere
real(dp), dimension(0:jmax, 0:imax) dzl_dtau
dzl_dtau(j,i): Derivative of zl by tau (time)
real(dp), dimension(0:jmax, 0:imax) sq_g22_g
sq_g22_g(j,i): Square root of the coefficient g22 of the metric tensor on grid point (i...