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