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
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 real(dp) :: dtime0, dtime_temp0, dtime_wss0, dtime_out0, dtime_ser0
95 real(dp) :: time_init0, time_end0, time_output0(100)
96 real(dp) :: d_dummy
97 character(len=100) :: anfdatname
98 character(len=256) :: filename_with_path
99 character(len=256) :: shell_command
100 character :: ch_dummy
101 logical :: flag_init_output, flag_3d_output
102 
103 character(len=64), parameter :: fmt1 = '(a)', &
104  fmt2 = '(a,i4)', &
105  fmt2a = '(a,i0)', &
106  fmt3 = '(a,es12.4)'
107 
108 write(unit=6, fmt='(a)') ' '
109 write(unit=6, fmt='(a)') ' -------- sico_init --------'
110 write(unit=6, fmt='(a)') ' '
111 
112 !-------- Name of the computational domain --------
113 
114 #if (defined(ANT))
115 ch_domain_long = 'Antarctica'
116 ch_domain_short = 'ant'
117 
118 #elif (defined(ASF))
119 ch_domain_long = 'Austfonna'
120 ch_domain_short = 'asf'
121 
122 #elif (defined(EMTP2SGE))
123 ch_domain_long = 'EISMINT Phase 2 Simplified Geometry Experiment'
124 ch_domain_short = 'emtp2sge'
125 
126 #elif (defined(GRL))
127 ch_domain_long = 'Greenland'
128 ch_domain_short = 'grl'
129 
130 #elif (defined(NHEM))
131 ch_domain_long = 'Northern hemisphere'
132 ch_domain_short = 'nhem'
133 
134 #elif (defined(SCAND))
135 ch_domain_long = 'Scandinavia and Eurasia'
136 ch_domain_short = 'scand'
137 
138 #elif (defined(TIBET))
139 ch_domain_long = 'Tibet'
140 ch_domain_short = 'tibet'
141 
142 #elif (defined(NMARS))
143 ch_domain_long = 'North polar cap of Mars'
144 ch_domain_short = 'nmars'
145 
146 #elif (defined(SMARS))
147 ch_domain_long = 'South polar cap of Mars'
148 ch_domain_short = 'smars'
149 
150 #elif (defined(XYZ))
151 ch_domain_long = 'XYZ'
152 ch_domain_short = 'xyz'
153 #if (defined(HEINO))
154 ch_domain_long = trim(ch_domain_long)//'/ISMIP HEINO'
155 #endif
156 
157 #else
158 
159 errormsg = ' >>> sico_init: No valid domain specified!'
160 call error(errormsg)
161 
162 #endif
163 
164 !-------- Some initial values --------
165 
166 n_output = 0
167 
168 dtime = 0.0_dp
169 dtime_temp = 0.0_dp
170 dtime_wss = 0.0_dp
171 dtime_out = 0.0_dp
172 dtime_ser = 0.0_dp
173 
174 time = 0.0_dp
175 time_init = 0.0_dp
176 time_end = 0.0_dp
177 time_output = 0.0_dp
178 
179 !-------- Initialisation of the Library of Iterative Solvers Lis,
180 ! if required --------
181 
182 #if (CALCTHK==3 || CALCTHK==6 || MARGIN==3 || DYNAMICS==2)
183  call lis_initialize(ierr)
184 #endif
185 
186 !-------- Read physical parameters --------
187 
188 call read_phys_para()
189 
190 call ice_mat_eqs_pars(rf, r_t, kappa, c, -190, 10)
191 
192 #if (!defined(SURFACE_FORCING) || SURFACE_FORCING==1)
193 
194 temp_min = temp_min ! deg C
195 s_t = s_t *1.0e-03_dp ! K/km -> K/m
196 x_hat = x_hat *1.0e+03_dp ! km -> m
197 y_hat = y_hat *1.0e+03_dp ! km -> m
198 b_max = b_max /year_sec ! m/a -> m/s
199 s_b = s_b *1.0e-03_dp/year_sec ! m/(a*km) -> 1/s
200 eld = eld *1.0e+03_dp ! km -> m
201 
202 #elif (SURFACE_FORCING==2)
203 
204 temp_0 = temp_0 ! deg C
205 gamma_t = gamma_t *1.0e-03_dp ! K/km -> K/m
206 s_0 = s_0 /year_sec ! m/a -> m/s
207 m_0 = m_0 *1.0e-03_dp/year_sec ! m/(a*km) -> 1/s
208 ela = ela *1.0e+03_dp ! km -> m
209 
210 #else
211 
212 errormsg = ' >>> sico_init: SURFACE_FORCING must be either 1 or 2!'
213 call error(errormsg)
214 
215 #endif
216 
217 ! ------ Some auxiliary quantities required for the enthalpy method
218 
219 call calc_c_int_table(c, -190, 10, l)
221 
222 !-------- Compatibility check of the SICOPOLIS version with the header file
223 
224 if ( trim(version) /= trim(sico_version) ) then
225  errormsg = ' >>> sico_init: ' &
226  //'SICOPOLIS version not compatible with header file!'
227  call error(errormsg)
228 end if
229 
230 !-------- Check whether the dynamics and thermodynamics modes are defined
231 
232 #if (!defined(DYNAMICS))
233 errormsg = ' >>> sico_init: DYNAMICS not defined in the header file!'
234 call error(errormsg)
235 #endif
236 
237 #if (!defined(CALCMOD))
238 errormsg = ' >>> sico_init: CALCMOD not defined in the header file!'
239 call error(errormsg)
240 #endif
241 
242 #if (defined(ENTHMOD))
243 errormsg = ' >>> sico_init: ENTHMOD must not be defined any more.' &
244  // end_of_line &
245  //' Please update your header file!'
246 call error(errormsg)
247 #endif
248 
249 !-------- Compatibility check of the horizontal resolution with the
250 ! number of grid points --------
251 
252 #if (GRID==0)
253 
254 if (approx_equal(dx, 5.0_dp, eps_sp_dp)) then
255 
256  if ((imax /= 300).or.(jmax /= 300)) then
257  errormsg = ' >>> sico_init: IMAX and/or JMAX wrong!'
258  call error(errormsg)
259  end if
260 
261 else if (approx_equal(dx, 10.0_dp, eps_sp_dp)) then
262 
263  if ((imax /= 150).or.(jmax /= 150)) then
264  errormsg = ' >>> sico_init: IMAX and/or JMAX wrong!'
265  call error(errormsg)
266  end if
267 
268 else if (approx_equal(dx, 25.0_dp, eps_sp_dp)) then
269 
270  if ((imax /= 60).or.(jmax /= 60)) then
271  errormsg = ' >>> sico_init: IMAX and/or JMAX wrong!'
272  call error(errormsg)
273  end if
274 
275 else if (approx_equal(dx, 75.0_dp, eps_sp_dp)) then
276 
277  if ((imax /= 20).or.(jmax /= 20)) then
278  errormsg = ' >>> sico_init: IMAX and/or JMAX wrong!'
279  call error(errormsg)
280  end if
281 
282 else
283 
284  errormsg = ' >>> sico_init: DX wrong!'
285  call error(errormsg)
286 
287 end if
288 
289 #elif (GRID==1)
290 
291 errormsg = ' >>> sico_init: GRID==1 not allowed for this application!'
292 call error(errormsg)
293 
294 #elif (GRID==2)
295 
296 errormsg = ' >>> sico_init: GRID==2 not allowed for this application!'
297 call error(errormsg)
298 
299 #endif
300 
301 !-------- Compatibility check of the thermodynamics mode
302 ! (cold vs. polythermal vs. enthalpy method)
303 ! and the number of grid points in the lower (kt) ice domain --------
304 
305 #if (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1)
306 
307 if (ktmax > 2) then
308  write(6, fmt='(a)') ' >>> sico_init: For options CALCMOD==0, 2, 3 or -1,'
309  write(6, fmt='(a)') ' the separate kt domain is redundant.'
310  write(6, fmt='(a)') ' Therefore, consider setting KTMAX to 2.'
311  write(6, fmt='(a)') ' '
312 end if
313 
314 #endif
315 
316 !-------- Compatibility check of discretization schemes for the horizontal and
317 ! vertical advection terms in the temperature and age equations --------
318 
319 #if (ADV_HOR==1)
320 errormsg = ' >>> sico_init: ' &
321  //'Option ADV_HOR==1 (central differences) not defined!'
322 call error(errormsg)
323 #endif
324 
325 !-------- Check whether for the shallow shelf
326 ! or shelfy stream approximation
327 ! the chosen grid is Cartesian coordinates
328 ! without distortion correction (GRID==0) --------
329 
330 #if ((MARGIN==3 && DYNAMICS==1) || DYNAMICS==2) /* requires SSA or SStA */
331 #if (GRID != 0)
332 errormsg = ' >>> sico_init: Distortion correction not implemented' &
333  // end_of_line &
334  //' for the shallow shelf approximation (SSA)' &
335  // end_of_line &
336  //' or the shelfy stream approximation (SStA)' &
337  // end_of_line &
338  //' -> GRID==0 required!'
339 call error(errormsg)
340 #endif
341 #endif
342 
343 !-------- Setting of forcing flag --------
344 
345 forcing_flag = 1 ! forcing by delta_ts
346 
347 !-------- Initialization of numerical time steps --------
348 
349 dtime0 = dtime0
350 dtime_temp0 = dtime_temp0
351 
352 !-------- Further initializations --------
353 
354 dzeta_c = 1.0_dp/real(kcmax,dp)
355 dzeta_t = 1.0_dp/real(ktmax,dp)
356 dzeta_r = 1.0_dp/real(krmax,dp)
357 
358 ndat2d = 1
359 ndat3d = 1
360 
361 !-------- General abbreviations --------
362 
363 ! ------ kc domain
364 
365 if (deform >= eps) then
366 
367  flag_aa_nonzero = .true. ! non-equidistant grid
368 
369  aa = deform
370  ea = exp(aa)
371 
372  kc=0
373  zeta_c(kc) = 0.0_dp
374  eaz_c(kc) = 1.0_dp
375  eaz_c_quotient(kc) = 0.0_dp
376 
377  do kc=1, kcmax-1
378  zeta_c(kc) = kc*dzeta_c
379  eaz_c(kc) = exp(aa*zeta_c(kc))
380  eaz_c_quotient(kc) = (eaz_c(kc)-1.0_dp)/(ea-1.0_dp)
381  end do
382 
383  kc=kcmax
384  zeta_c(kc) = 1.0_dp
385  eaz_c(kc) = exp(aa)
386  eaz_c_quotient(kc) = 1.0_dp
387 
388 else
389 
390  flag_aa_nonzero = .false. ! equidistant grid
391 
392  aa = 0.0_dp
393  ea = 1.0_dp
394 
395  kc=0
396  zeta_c(kc) = 0.0_dp
397  eaz_c(kc) = 1.0_dp
398  eaz_c_quotient(kc) = 0.0_dp
399 
400  do kc=1, kcmax-1
401  zeta_c(kc) = kc*dzeta_c
402  eaz_c(kc) = 1.0_dp
403  eaz_c_quotient(kc) = zeta_c(kc)
404  end do
405 
406  kc=kcmax
407  zeta_c(kc) = 1.0_dp
408  eaz_c(kc) = 1.0_dp
409  eaz_c_quotient(kc) = 1.0_dp
410 
411 end if
412 
413 ! ------ kt domain
414 
415 kt=0
416 zeta_t(kt) = 0.0_dp
417 
418 do kt=1, ktmax-1
419  zeta_t(kt) = kt*dzeta_t
420 end do
421 
422 kt=ktmax
423 zeta_t(kt) = 1.0_dp
424 
425 ! ------ kr domain
426 
427 kr=0
428 zeta_r(kr) = 0.0_dp
429 
430 do kr=1, krmax-1
431  zeta_r(kr) = kr*dzeta_r
432 end do
433 
434 kr=krmax
435 zeta_r(kr) = 1.0_dp
436 
437 !-------- Reshaping of a 2-d array (with indices i, j)
438 ! to a vector (with index n) --------
439 
440 n=1
441 
442 do i=0, imax
443 do j=0, jmax
444  n2i(n) = i
445  n2j(n) = j
446  ij2n(j,i) = n
447  n=n+1
448 end do
449 end do
450 
451 !-------- Specification of current simulation --------
452 
453 runname = runname
454 anfdatname = anfdatname
455 
456 #if (defined(YEAR_ZERO))
458 #else
459 year_zero = 2000.0_dp ! default value 2000 CE
460 #endif
461 
462 time_init0 = time_init0
463 time_end0 = time_end0
464 dtime_ser0 = dtime_ser0
465 
466 #if (OUTPUT==1 || OUTPUT==3)
467 dtime_out0 = dtime_out0
468 #endif
469 
470 #if (OUTPUT==2 || OUTPUT==3)
471 n_output = n_output
472 time_output0( 1) = time_out0_01
473 time_output0( 2) = time_out0_02
474 time_output0( 3) = time_out0_03
475 time_output0( 4) = time_out0_04
476 time_output0( 5) = time_out0_05
477 time_output0( 6) = time_out0_06
478 time_output0( 7) = time_out0_07
479 time_output0( 8) = time_out0_08
480 time_output0( 9) = time_out0_09
481 time_output0(10) = time_out0_10
482 time_output0(11) = time_out0_11
483 time_output0(12) = time_out0_12
484 time_output0(13) = time_out0_13
485 time_output0(14) = time_out0_14
486 time_output0(15) = time_out0_15
487 time_output0(16) = time_out0_16
488 time_output0(17) = time_out0_17
489 time_output0(18) = time_out0_18
490 time_output0(19) = time_out0_19
491 time_output0(20) = time_out0_20
492 #endif
493 
494 !-------- Write log file --------
495 
496 shell_command = 'if [ ! -d'
497 shell_command = trim(shell_command)//' '//outpath
498 shell_command = trim(shell_command)//' '//'] ; then mkdir'
499 shell_command = trim(shell_command)//' '//outpath
500 shell_command = trim(shell_command)//' '//'; fi'
501 call system(trim(shell_command))
502  ! Check whether directory OUTPATH exists. If not, it is created.
503 
504 filename_with_path = trim(outpath)//'/'//trim(runname)//'.log'
505 
506 open(10, iostat=ios, file=trim(filename_with_path), status='new')
507 
508 if (ios /= 0) then
509  errormsg = ' >>> sico_init: Error when opening the log file!'
510  call error(errormsg)
511 end if
512 
513 write(10, fmt=trim(fmt1)) 'Computational domain:'
514 write(10, fmt=trim(fmt1)) trim(ch_domain_long)
515 write(10, fmt=trim(fmt1)) ' '
516 
517 write(10, fmt=trim(fmt2a)) 'GRID = ', grid
518 write(10, fmt=trim(fmt1)) ' '
519 
520 write(10, fmt=trim(fmt2)) 'imax =', imax
521 write(10, fmt=trim(fmt2)) 'jmax =', jmax
522 write(10, fmt=trim(fmt2)) 'kcmax =', kcmax
523 write(10, fmt=trim(fmt2)) 'ktmax =', ktmax
524 write(10, fmt=trim(fmt2)) 'krmax =', krmax
525 write(10, fmt=trim(fmt1)) ' '
526 
527 write(10, fmt=trim(fmt3)) 'a =', aa
528 write(10, fmt=trim(fmt1)) ' '
529 
530 #if (GRID==0 || GRID==1)
531 write(10, fmt=trim(fmt3)) 'x0 =', x0
532 write(10, fmt=trim(fmt3)) 'y0 =', y0
533 write(10, fmt=trim(fmt3)) 'dx =', dx
534 #elif (GRID==2)
535 write(10, fmt=trim(fmt3)) 'dlambda =', dlambda
536 write(10, fmt=trim(fmt3)) 'dphi =', dphi
537 #endif
538 write(10, fmt=trim(fmt1)) ' '
539 
540 write(10, fmt=trim(fmt3)) 'year_zero =', year_zero
541 write(10, fmt=trim(fmt3)) 'time_init =', time_init0
542 write(10, fmt=trim(fmt3)) 'time_end =', time_end0
543 write(10, fmt=trim(fmt3)) 'dtime =', dtime0
544 write(10, fmt=trim(fmt3)) 'dtime_temp =', dtime_temp0
545 write(10, fmt=trim(fmt1)) ' '
546 
547 write(10, fmt=trim(fmt2)) 'ANF_DAT =', anf_dat
548 #if (ANF_DAT==1 && defined(TEMP_INIT))
549 write(10, fmt=trim(fmt2)) 'TEMP_INIT =', temp_init
550 #endif
551 #if (ANF_DAT==3 || (ANF_DAT==1 && TEMP_INIT==5))
552 write(10, fmt=trim(fmt1)) 'Initial-value file = '//anfdatname
553 write(10, fmt=trim(fmt1)) 'Path to initial-value file = '//anfdatpath
554 #endif
555 write(10, fmt=trim(fmt1)) ' '
556 
557 write(10, fmt=trim(fmt1)) 'Physical-parameter file = '//phys_para_file
558 write(10, fmt=trim(fmt1)) ' '
559 
560 #if (defined(THK_EVOL))
561 write(10, fmt=trim(fmt2a)) 'THK_EVOL = ', thk_evol
562 #else
563 errormsg = ' >>> sico_init: Define THK_EVOL in header file!'
564 call error(errormsg)
565 #endif
566 #if (defined(CALCTHK))
567 write(10, fmt=trim(fmt2a)) 'CALCTHK = ', calcthk
568 #else
569 errormsg = ' >>> sico_init: Define CALCTHK in header file!'
570 call error(errormsg)
571 #endif
572 
573 #if (defined(H_ISOL_MAX))
574 write(10, fmt=trim(fmt3)) 'H_isol_max =', h_isol_max
575 #endif
576 
577 #if (CALCTHK==2 || CALCTHK==3 || CALCTHK==5 || CALCTHK==6)
578 write(10, fmt=trim(fmt3)) 'ovi_weight =', ovi_weight
579 #if (CALCTHK==2 || CALCTHK==5)
580 write(10, fmt=trim(fmt3)) 'omega_sor =', omega_sor
581 #if (ITER_MAX_SOR>0)
582 write(10, fmt=trim(fmt2a)) 'iter_max_sor = ', iter_max_sor
583 #endif
584 #endif
585 #endif
586 
587 write(10, fmt=trim(fmt1)) ' '
588 
589 #if (defined(SURFACE_FORCING))
590 write(10, fmt=trim(fmt2)) 'SURFACE_FORCING =', surface_forcing
591 #endif
592 #if (!defined(SURFACE_FORCING) || SURFACE_FORCING==1)
593 write(10, fmt=trim(fmt3)) 'temp_min =', temp_min
594 write(10, fmt=trim(fmt3)) 's_t =', s_t
595 write(10, fmt=trim(fmt3)) 'x_hat =', x_hat
596 write(10, fmt=trim(fmt3)) 'y_hat =', y_hat
597 write(10, fmt=trim(fmt3)) 'b_max =', b_max
598 write(10, fmt=trim(fmt3)) 's_b =', s_b
599 write(10, fmt=trim(fmt3)) 'eld =', eld
600 #elif (SURFACE_FORCING==2)
601 write(10, fmt=trim(fmt3)) 'temp_0 =', temp_0
602 write(10, fmt=trim(fmt3)) 'gamma_t =', gamma_t
603 write(10, fmt=trim(fmt3)) 's_0 =', s_0
604 write(10, fmt=trim(fmt3)) 'm_0 =', m_0
605 write(10, fmt=trim(fmt3)) 'ela =', ela
606 #endif
607 write(10, fmt=trim(fmt1)) ' '
608 
609 #if (TSURFACE==1)
610 write(10, fmt=trim(fmt3)) 'delta_ts0 =', delta_ts0
611 #elif (TSURFACE==3)
612 write(10, fmt=trim(fmt3)) 'sine_amplit =', sine_amplit
613 write(10, fmt=trim(fmt3)) 'sine_period =', sine_period
614 #elif (TSURFACE==4)
615 write(10, fmt=trim(fmt1)) 'GRIP file = '//grip_temp_file
616 write(10, fmt=trim(fmt3)) 'grip_temp_fact =', grip_temp_fact
617 #endif
618 
619 write(10, fmt=trim(fmt2a)) 'SEA_LEVEL = ', sea_level
620 #if (SEA_LEVEL==1)
621 write(10, fmt=trim(fmt3)) 'z_sl0 =', z_sl0
622 #elif (SEA_LEVEL==3)
623 write(10, fmt=trim(fmt1)) 'sea-level file = '//sea_level_file
624 #endif
625 write(10, fmt=trim(fmt1)) ' '
626 
627 #if (MARGIN==2)
628 #if (MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3)
629 write(10, fmt=trim(fmt3)) 'z_mar =', z_mar
630 write(10, fmt=trim(fmt1)) ' '
631 #elif (MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \
632  || marine_ice_calving==6 || marine_ice_calving==7)
633 write(10, fmt=trim(fmt3)) 'fact_z_mar =', fact_z_mar
634 write(10, fmt=trim(fmt1)) ' '
635 #elif (MARINE_ICE_FORMATION==2 && MARINE_ICE_CALVING==9)
636 write(10, fmt=trim(fmt3)) 'calv_uw_coeff =', calv_uw_coeff
637 write(10, fmt=trim(fmt3)) 'r1_calv_uw =', r1_calv_uw
638 write(10, fmt=trim(fmt3)) 'r2_calv_uw =', r2_calv_uw
639 write(10, fmt=trim(fmt1)) ' '
640 #endif
641 #elif (MARGIN==3)
642 #if (ICE_SHELF_CALVING==2)
643 write(10, fmt=trim(fmt3)) 'H_calv =', h_calv
644 write(10, fmt=trim(fmt1)) ' '
645 #endif
646 #endif
647 
648 #if (defined(BASAL_HYDROLOGY))
649 write(10, fmt=trim(fmt2a)) 'BASAL_HYDROLOGY = ', basal_hydrology
650 #endif
651 write(10, fmt=trim(fmt2a)) 'SLIDE_LAW = ', slide_law
652 write(10, fmt=trim(fmt3)) 'c_slide =', c_slide
653 write(10, fmt=trim(fmt3)) 'gamma_slide =', gamma_slide
654 write(10, fmt=trim(fmt2a)) 'p_weert = ', p_weert
655 write(10, fmt=trim(fmt2a)) 'q_weert = ', q_weert
656 #if (defined(TIME_RAMP_UP_SLIDE))
657 write(10, fmt=trim(fmt3)) 'time_ramp_up_slide =', time_ramp_up_slide
658 #endif
659 #if (SLIDE_LAW==2 || SLIDE_LAW==3)
660 write(10, fmt=trim(fmt3)) 'red_pres_limit_fact =', red_pres_limit_fact
661 #endif
662 #if (BASAL_HYDROLOGY==1 \
663  && defined(hydro_slide_sat_fct) \
664  && defined(c_hw_slide) && defined(hw0_slide))
665 write(10, fmt=trim(fmt2a)) 'HYDRO_SLIDE_SAT_FCT = ', hydro_slide_sat_fct
666 write(10, fmt=trim(fmt3)) 'c_Hw_slide =', c_hw_slide
667 write(10, fmt=trim(fmt3)) 'Hw0_slide =', hw0_slide
668 #endif
669 write(10, fmt=trim(fmt1)) ' '
670 
671 write(10, fmt=trim(fmt2a)) 'Q_GEO_MOD = ', q_geo_mod
672 #if (Q_GEO_MOD==1)
673 write(10, fmt=trim(fmt3)) 'q_geo =', q_geo
674 #endif
675 write(10, fmt=trim(fmt1)) ' '
676 
677 write(10, fmt=trim(fmt2)) 'REBOUND =', rebound
678 #if (REBOUND==1)
679 write(10, fmt=trim(fmt3)) 'frac_llra =', frac_llra
680 #endif
681 #if (REBOUND==1 || REBOUND==2)
682 write(10, fmt=trim(fmt2)) 'TIME_LAG_MOD =', time_lag_mod
683 #if (TIME_LAG_MOD==1)
684 write(10, fmt=trim(fmt3)) 'time_lag =', time_lag
685 #elif (TIME_LAG_MOD==2)
686 write(10, fmt=trim(fmt1)) 'time_lag_file = '//time_lag_file
687 #else
688 errormsg = ' >>> sico_init: TIME_LAG_MOD must be either 1 or 2!'
689 call error(errormsg)
690 #endif
691 #endif
692 #if (REBOUND==2)
693 write(10, fmt=trim(fmt2)) 'FLEX_RIG_MOD =', flex_rig_mod
694 #if (FLEX_RIG_MOD==1)
695 write(10, fmt=trim(fmt3)) 'flex_rig =', flex_rig
696 #elif (FLEX_RIG_MOD==2)
697 write(10, fmt=trim(fmt1)) 'flex_rig_file = '//flex_rig_file
698 #else
699 errormsg = ' >>> sico_init: FLEX_RIG_MOD must be either 1 or 2!'
700 call error(errormsg)
701 #endif
702 #endif
703 write(10, fmt=trim(fmt2)) 'Q_LITHO =', q_litho
704 write(10, fmt=trim(fmt1)) ' '
705 
706 #if (FLOW_LAW==2)
707 write(10, fmt=trim(fmt3)) 'gr_size =', gr_size
708 write(10, fmt=trim(fmt1)) ' '
709 #endif
710 #if (FIN_VISC==2)
711 write(10, fmt=trim(fmt3)) 'sigma_res =', sigma_res
712 write(10, fmt=trim(fmt1)) ' '
713 #endif
714 
715 write(10, fmt=trim(fmt2)) 'ENHMOD =', enhmod
716 #if (ENHMOD==1 || ENHMOD==2 || ENHMOD==3)
717 write(10, fmt=trim(fmt3)) 'enh_fact =', enh_fact
718 #endif
719 #if (ENHMOD==2 || ENHMOD==3)
720 write(10, fmt=trim(fmt3)) 'enh_intg =', enh_intg
721 #endif
722 #if (ENHMOD==2)
723 write(10, fmt=trim(fmt3)) 'age_trans =', age_trans_0
724 #endif
725 #if (ENHMOD==3)
726 write(10, fmt=trim(fmt3)) 'date_trans1 =', date_trans1_0
727 write(10, fmt=trim(fmt3)) 'date_trans2 =', date_trans2_0
728 write(10, fmt=trim(fmt3)) 'date_trans3 =', date_trans3_0
729 #endif
730 #if (ENHMOD==4 || ENHMOD==5)
731 write(10, fmt=trim(fmt3)) 'enh_compr =', enh_compr
732 write(10, fmt=trim(fmt3)) 'enh_shear =', enh_shear
733 #endif
734 #if ((ENHMOD==1 || ENHMOD==2 || ENHMOD==3 || ENHMOD==4) && MARGIN==3)
735 write(10, fmt=trim(fmt3)) 'enh_shelf =', enh_shelf
736 #endif
737 write(10, fmt=trim(fmt1)) ' '
738 
739 write(10, fmt=trim(fmt2a)) 'DYNAMICS = ', dynamics
740 #if ((DYNAMICS==1 && MARGIN==3) || DYNAMICS==2)
741 #if (defined(LIS_OPTS))
742 write(10, fmt=trim(fmt1)) 'lis_opts = '//lis_opts
743 #endif
744 #if (defined(TOL_ITER_SSA))
745 write(10, fmt=trim(fmt3)) 'tol_iter_ssa =', tol_iter_ssa
746 #endif
747 #if (defined(N_ITER_SSA))
748 write(10, fmt=trim(fmt2a)) 'n_iter_ssa = ', n_iter_ssa
749 #endif
750 #if (defined(ITER_INIT_SSA))
751 write(10, fmt=trim(fmt2a)) 'iter_init_ssa = ', iter_init_ssa
752 #endif
753 #if (defined(VISC_INIT_SSA))
754 write(10, fmt=trim(fmt3)) 'visc_init_ssa =', visc_init_ssa
755 #endif
756 #if (defined(N_VISC_SMOOTH))
757 write(10, fmt=trim(fmt2a)) 'n_visc_smooth = ', n_visc_smooth
758 #endif
759 #if (defined(VISC_SMOOTH_DIFF))
760 write(10, fmt=trim(fmt3)) 'visc_smooth_diff =', visc_smooth_diff
761 #endif
762 #if (defined(RELAX_FACT_SSA))
763 write(10, fmt=trim(fmt3)) 'relax_fact_ssa =', relax_fact_ssa
764 #endif
765 #endif
766 #if (DYNAMICS==2 && defined(RATIO_SL_THRESH))
767 write(10, fmt=trim(fmt3)) 'ratio_sl_thresh =', ratio_sl_thresh
768 #endif
769 #if (DYNAMICS==2 && defined(SSTA_SIA_WEIGH_FCT))
770 write(10, fmt=trim(fmt2a)) 'SSTA_SIA_WEIGH_FCT = ', ssta_sia_weigh_fct
771 #endif
772 write(10, fmt=trim(fmt1)) ' '
773 
774 write(10, fmt=trim(fmt2)) 'CALCMOD =', calcmod
775 #if (CALCMOD==-1 && defined(TEMP_CONST))
776 write(10, fmt=trim(fmt3)) 'TEMP_CONST =', temp_const
777 #endif
778 #if (CALCMOD==-1 && defined(AGE_CONST))
779 write(10, fmt=trim(fmt3)) 'AGE_CONST =', age_const
780 #endif
781 #if (CALCMOD==1 && defined(CTS_MELTING_FREEZING))
782 write(10, fmt=trim(fmt2)) 'CTS_MELTING_FREEZING =', cts_melting_freezing
783 #endif
784 write(10, fmt=trim(fmt2)) 'FLOW_LAW =', flow_law
785 write(10, fmt=trim(fmt2)) 'FIN_VISC =', fin_visc
786 write(10, fmt=trim(fmt2)) 'MARGIN =', margin
787 #if (MARGIN==2)
788 write(10, fmt=trim(fmt2)) 'MARINE_ICE_FORMATION =', marine_ice_formation
789 write(10, fmt=trim(fmt2)) 'MARINE_ICE_CALVING =', marine_ice_calving
790 #elif (MARGIN==3)
791 write(10, fmt=trim(fmt2)) 'ICE_SHELF_CALVING =', ice_shelf_calving
792 #endif
793 write(10, fmt=trim(fmt2)) 'ADV_HOR =', adv_hor
794 write(10, fmt=trim(fmt2)) 'ADV_VERT =', adv_vert
795 write(10, fmt=trim(fmt2)) 'TOPOGRAD =', topograd
796 write(10, fmt=trim(fmt2)) 'TSURFACE =', tsurface
797 #if (defined(MB_ACCOUNT))
798 write(10, fmt=trim(fmt2a)) 'MB_ACCOUNT = ', mb_account
799 #endif
800 write(10, fmt=trim(fmt1)) ' '
801 
802 write(10, fmt=trim(fmt3)) 'numdiff_H_t =', numdiff_h_t
803 write(10, fmt=trim(fmt3)) 'tau_cts =', tau_cts
804 write(10, fmt=trim(fmt3)) 'vh_max =', vh_max
805 write(10, fmt=trim(fmt3)) 'hd_min =', hd_min
806 write(10, fmt=trim(fmt3)) 'hd_max =', hd_max
807 #if (defined(VISC_MIN) && defined(VISC_MAX))
808 write(10, fmt=trim(fmt3)) 'visc_min =', visc_min
809 write(10, fmt=trim(fmt3)) 'visc_max =', visc_max
810 #endif
811 write(10, fmt=trim(fmt3)) 'qbm_min =', qbm_min
812 write(10, fmt=trim(fmt3)) 'qbm_max =', qbm_max
813 write(10, fmt=trim(fmt3)) 'age_min =', age_min
814 write(10, fmt=trim(fmt3)) 'age_max =', age_max
815 write(10, fmt=trim(fmt3)) 'mean_accum =', mean_accum
816 #if (ADV_VERT==1)
817 write(10, fmt=trim(fmt3)) 'age_diff =', agediff
818 #endif
819 write(10, fmt=trim(fmt1)) ' '
820 
821 #if (defined(OUT_TIMES))
822 write(10, fmt=trim(fmt2a)) 'OUT_TIMES = ', out_times
823 #endif
824 #if (defined(OUTPUT_INIT))
825 write(10, fmt=trim(fmt2a)) 'OUTPUT_INIT = ', output_init
826 #endif
827 write(10, fmt=trim(fmt2a)) 'OUTPUT = ', output
828 #if (OUTPUT==1 || OUTPUT==3)
829 write(10, fmt=trim(fmt3)) 'dtime_out =' , dtime_out0
830 #endif
831 write(10, fmt=trim(fmt3)) 'dtime_ser =' , dtime_ser0
832 #if (OUTPUT==1 || OUTPUT==2)
833 write(10, fmt=trim(fmt2a)) 'ERGDAT = ', ergdat
834 #endif
835 write(10, fmt=trim(fmt2a)) 'NETCDF = ', netcdf
836 #if (OUTPUT==2 || OUTPUT==3)
837 write(10, fmt=trim(fmt2a)) 'n_output = ', n_output
838 do n=1, n_output
839  write(10, fmt=trim(fmt3)) 'time_output =' , time_output0(n)
840 end do
841 #endif
842 write(10, fmt=trim(fmt1)) ' '
843 
844 write(10, fmt=trim(fmt1)) 'Program version and date: '//version//' / '//date
845 
846 close(10, status='keep')
847 
848 !-------- Conversion of time quantities --------
849 
850 year_zero = year_zero*year_sec ! a --> s
851 time_init = time_init0*year_sec ! a --> s
852 time_end = time_end0*year_sec ! a --> s
853 dtime = dtime0*year_sec ! a --> s
854 dtime_temp = dtime_temp0*year_sec ! a --> s
855 dtime_ser = dtime_ser0*year_sec ! a --> s
856 #if (OUTPUT==1 || OUTPUT==3)
857 dtime_out = dtime_out0*year_sec ! a --> s
858 #endif
859 #if (OUTPUT==2 || OUTPUT==3)
860 do n=1, n_output
861  time_output(n) = time_output0(n)*year_sec ! a --> s
862 end do
863 #endif
864 
865 if (.not.approx_integer_multiple(dtime_temp, dtime, eps_sp_dp)) then
866  errormsg = ' >>> sico_init: dtime_temp must be a multiple of dtime!'
867  call error(errormsg)
868 end if
869 
870 if (.not.approx_integer_multiple(dtime_ser, dtime, eps_sp_dp)) then
871  errormsg = ' >>> sico_init: dtime_ser must be a multiple of dtime!'
872  call error(errormsg)
873 end if
874 
875 #if (OUTPUT==1 || OUTPUT==3)
876 if (.not.approx_integer_multiple(dtime_out, dtime, eps_sp_dp)) then
877  errormsg = ' >>> sico_init: dtime_out must be a multiple of dtime!'
878  call error(errormsg)
879 end if
880 #endif
881 
882 time = time_init
883 
884 !-------- Mean accumulation --------
885 
886 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho)
887 ! ! mm/a water equiv. --> m/s ice equiv.
888 
889 !-------- Read data for delta_ts --------
890 
891 #if (TSURFACE==4)
892 
893 filename_with_path = trim(inpath)//'/general/'//trim(grip_temp_file)
894 
895 open(21, iostat=ios, file=trim(filename_with_path), status='old')
896 
897 if (ios /= 0) then
898  errormsg = ' >>> sico_init: Error when opening the data file for delta_ts!'
899  call error(errormsg)
900 end if
901 
902 read(21, fmt=*) ch_dummy, grip_time_min, grip_time_stp, grip_time_max
903 
904 if (ch_dummy /= '#') then
905  errormsg = ' >>> sico_init: grip_time_min, grip_time_stp, grip_time_max' &
906  // end_of_line &
907  //' not defined in data file for delta_ts!'
908  call error(errormsg)
909 end if
910 
912 
913 allocate(griptemp(0:ndata_grip))
914 
915 do n=0, ndata_grip
916  read(21, fmt=*) d_dummy, griptemp(n)
917 end do
918 
919 close(21, status='keep')
920 
921 #endif
922 
923 !-------- Read data for z_sl --------
924 
925 #if (SEA_LEVEL==3)
926 
927 filename_with_path = trim(inpath)//'/general/'//trim(sea_level_file)
928 
929 open(21, iostat=ios, file=trim(filename_with_path), status='old')
930 
931 if (ios /= 0) then
932  errormsg = ' >>> sico_init: Error when opening the data file for z_sl!'
933  call error(errormsg)
934 end if
935 
936 read(21, fmt=*) ch_dummy, specmap_time_min, specmap_time_stp, specmap_time_max
937 
938 if (ch_dummy /= '#') then
939  errormsg = ' >>> sico_init:' &
940  // end_of_line &
941  //' specmap_time_min, specmap_time_stp, specmap_time_max' &
942  // end_of_line &
943  //' not defined in data file for z_sl!'
944  call error(errormsg)
945 end if
946 
948 
949 allocate(specmap_zsl(0:ndata_specmap))
950 
951 do n=0, ndata_specmap
952  read(21, fmt=*) d_dummy, specmap_zsl(n)
953 end do
954 
955 close(21, status='keep')
956 
957 #endif
958 
959 !-------- Determination of the geothermal heat flux --------
960 
961 #if (Q_GEO_MOD==1)
962 
963 ! ------ Constant value
964 
965 do i=0, imax
966 do j=0, jmax
967  q_geo(j,i) = q_geo *1.0e-03_dp ! mW/m2 -> W/m2
968 end do
969 end do
970 
971 #elif (Q_GEO_MOD==2)
972 
973 errormsg = ' >>> sico_init: ' &
974  //'Option Q_GEO_MOD==2 not available for this application!'
975 call error(errormsg)
976 
977 #endif
978 
979 !-------- Determination of the time lag
980 ! of the relaxing asthenosphere --------
981 
982 #if (REBOUND==1 || REBOUND==2)
983 
984 #if (TIME_LAG_MOD==1)
985 
986 time_lag_asth = time_lag*year_sec ! a -> s
987 
988 #elif (TIME_LAG_MOD==2)
989 
990 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
991  trim(time_lag_file)
992 
993 open(29, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
994 
995 if (ios /= 0) then
996  errormsg = ' >>> sico_init: Error when opening the time-lag file!'
997  call error(errormsg)
998 end if
999 
1000 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1001 
1002 do j=jmax, 0, -1
1003  read(29, fmt=*) (time_lag_asth(j,i), i=0,imax)
1004 end do
1005 
1006 close(29, status='keep')
1007 
1008 time_lag_asth = time_lag_asth*year_sec ! a -> s
1009 
1010 #endif
1011 
1012 #elif (REBOUND==0)
1013 
1014 time_lag_asth = 0.0_dp ! dummy values
1015 
1016 #endif
1017 
1018 !-------- Determination of the flexural rigidity of the lithosphere --------
1019 
1020 #if (REBOUND==2)
1021 
1022 #if (FLEX_RIG_MOD==1)
1023 
1024 flex_rig_lith = flex_rig
1025 
1026 #elif (FLEX_RIG_MOD==2)
1027 
1028 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1029  trim(flex_rig_file)
1030 
1031 open(29, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1032 
1033 if (ios /= 0) then
1034  errormsg = ' >>> sico_init: Error when opening the flex-rig file!'
1035  call error(errormsg)
1036 end if
1037 
1038 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1039 
1040 do j=jmax, 0, -1
1041  read(29, fmt=*) (flex_rig_lith(j,i), i=0,imax)
1042 end do
1043 
1044 close(29, status='keep')
1045 
1046 #endif
1047 
1048 #elif (REBOUND==0 || REBOUND==1)
1049 
1050 flex_rig_lith = 0.0_dp ! dummy values
1051 
1052 #endif
1053 
1054 !-------- Definition of initial values --------
1055 
1056 ! ------ Present topography
1057 
1058 #if (ANF_DAT==1)
1059 
1060 errormsg = ' >>> sico_init: ANF_DAT==1 not allowed for this application!'
1061 call error(errormsg)
1062 
1063 ! ------ Ice-free, relaxed bedrock
1064 
1065 #elif (ANF_DAT==2)
1066 
1067 call topography2(dxi, deta)
1068 
1069 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1070 
1071 call boundary(time_init, dtime, dxi, deta, &
1072  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1073 
1074 as_perp_apl = 0.0_dp
1075 
1076 smb_corr = 0.0_dp
1077 
1078 q_bm = 0.0_dp
1079 q_tld = 0.0_dp
1080 q_b_tot = 0.0_dp
1081 
1082 p_b_w = 0.0_dp
1083 h_w = 0.0_dp
1084 
1085 call init_temp_water_age_2()
1086 
1087 #if (ENHMOD==1)
1088  call calc_enhance_1()
1089 #elif (ENHMOD==2)
1090  call calc_enhance_2()
1091 #elif (ENHMOD==3)
1092  call calc_enhance_3(time_init)
1093 #elif (ENHMOD==4)
1094  call calc_enhance_4()
1095 #elif (ENHMOD==5)
1096  call calc_enhance_5()
1097 #else
1098  errormsg = ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
1099  call error(errormsg)
1100 #endif
1101 
1102 vx_m_sia = 0.0_dp
1103 vy_m_sia = 0.0_dp
1104 vx_m_ssa = 0.0_dp
1105 vy_m_ssa = 0.0_dp
1106 
1107 #if (defined(VISC_MIN) && defined(VISC_MAX))
1108  vis_ave_g = visc_max
1109 #else
1110  vis_ave_g = 1.0e+25_dp ! Pa s
1111 #endif
1112 
1113 vis_int_g = 0.0_dp
1114 
1115 ! ------ Read initial state from output of previous simulation
1116 
1117 #elif (ANF_DAT==3)
1118 
1119 call topography3(dxi, deta, z_sl, anfdatname)
1120 
1121 call boundary(time_init, dtime, dxi, deta, &
1122  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1123 
1124 where ((maske==0_i1b).or.(maske==3_i1b))
1125  ! grounded or floating ice
1127 elsewhere ! maske==1_i1b or 2_i1b, ice-free land or sea
1128  as_perp_apl = 0.0_dp
1129 end where
1130 
1131 smb_corr = 0.0_dp
1132 
1133 q_b_tot = q_bm + q_tld
1134 
1135 #endif
1136 
1137 !-------- Inner-point flag --------
1138 
1139 flag_inner_point = .true.
1140 
1141 flag_inner_point(0,:) = .false.
1142 flag_inner_point(jmax,:) = .false.
1143 
1144 flag_inner_point(:,0) = .false.
1145 flag_inner_point(:,imax) = .false.
1146 
1147 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
1148 
1149 #if (GRID==0 || GRID==1)
1150 
1151 do ir=-imax, imax
1152 do jr=-jmax, jmax
1153  dist_dxdy(jr,ir) = sqrt( (real(ir,dp)*dxi)**2 + (real(jr,dp)*deta)**2 )
1154  ! distortion due to stereographic projection not accounted for
1155 end do
1156 end do
1157 
1158 #endif
1159 
1160 !-------- Initial velocities --------
1161 
1162 call calc_temp_melt()
1163 call flag_update_gf_gl_cf()
1164 
1165 #if (DYNAMICS==1 || DYNAMICS==2)
1166 
1167 call calc_vxy_b_sia(time, z_sl)
1168 call calc_vxy_sia(dzeta_c, dzeta_t)
1169 
1170 #if (MARGIN==3 || DYNAMICS==2)
1171 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t)
1172 #endif
1173 
1174 call calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
1175 
1176 #if (MARGIN==3)
1177 call calc_vz_floating(z_sl, dxi, deta, dzeta_c)
1178 #endif
1179 
1180 #elif (DYNAMICS==0)
1181 
1182 call calc_vxy_static()
1183 call calc_vz_static()
1184 
1185 #else
1186 
1187 errormsg = ' >>> sico_init: DYNAMICS must be either 0, 1 or 2!'
1188 call error(errormsg)
1189 
1190 #endif
1191 
1192 call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
1193 
1194 !-------- Initial enthalpies --------
1195 
1196 #if (CALCMOD==0 || CALCMOD==-1)
1197 
1198 do i=0, imax
1199 do j=0, jmax
1200 
1201  do kc=0, kcmax
1202  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1203  end do
1204 
1205  do kt=0, ktmax
1206  enth_t(kt,j,i) = enth_c(0,j,i)
1207  end do
1208 
1209 end do
1210 end do
1211 
1212 #elif (CALCMOD==1)
1213 
1214 do i=0, imax
1215 do j=0, jmax
1216 
1217  do kc=0, kcmax
1218  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1219  end do
1220 
1221  if ( (maske(j,i) == 0_i1b).and.(n_cts(j,i) == 1_i1b) ) then
1222  do kt=0, ktmax
1223  enth_t(kt,j,i) = enth_fct_temp_omega(temp_t_m(kt,j,i), omega_t(kt,j,i))
1224  end do
1225  else
1226  do kt=0, ktmax
1227  enth_t(kt,j,i) = enth_c(0,j,i)
1228  end do
1229  end if
1230 
1231 end do
1232 end do
1233 
1234 #elif (CALCMOD==2 || CALCMOD==3)
1235 
1236 do i=0, imax
1237 do j=0, jmax
1238 
1239  do kc=0, kcmax
1240  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), omega_c(kc,j,i))
1241  end do
1242 
1243  do kt=0, ktmax
1244  enth_t(kt,j,i) = enth_c(0,j,i)
1245  end do
1246 
1247 end do
1248 end do
1249 
1250 #else
1251 
1252 errormsg = ' >>> sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!'
1253 call error(errormsg)
1254 
1255 #endif
1256 
1257 !-------- Initialize time-series files --------
1258 
1259 ! ------ Time-series file for the ice sheet on the whole
1260 
1261 filename_with_path = trim(outpath)//'/'//trim(runname)//'.ser'
1262 
1263 open(12, iostat=ios, file=trim(filename_with_path), status='new')
1264 
1265 if (ios /= 0) then
1266  errormsg = ' >>> sico_init: Error when opening the ser file!'
1267  call error(errormsg)
1268 end if
1269 
1270 write(12,1102)
1271 write(12,1103)
1272 
1273  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1274  ' V(m^3) V_g(m^3) V_f(m^3)', &
1275  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1276  ' V_sle(m) V_t(m^3)', &
1277  ' A_t(m^2)',/, &
1278  ' H_max(m) H_t_max(m)', &
1279  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
1280  1103 format('----------------------------------------------------', &
1281  '---------------------------------------')
1282 
1283 ! ------ Time-series file for selected positions ("deep boreholes")
1284 
1285 n_core = 2 ! central dome, position halfway to coast
1286 
1287 allocate(lambda_core(n_core), phi_core(n_core), &
1289 
1290 ch_core(1) = 'P1'
1291 lambda_core(1) = 0.0_dp ! dummy
1292 phi_core(1) = 0.0_dp ! dummy
1293 x_core(1) = 750.0_dp *1.0e+03_dp ! Position of the central dome
1294 y_core(1) = 750.0_dp *1.0e+03_dp ! (750 km, 750 km),
1295  ! conversion km -> m
1296 
1297 ch_core(2) = 'P2'
1298 lambda_core(2) = 0.0_dp ! dummy
1299 phi_core(2) = 0.0_dp ! dummy
1300 x_core(2) = 750.0_dp *1.0e+03_dp ! Position halfway to the coast
1301 y_core(2) = 1125.0_dp *1.0e+03_dp ! (750 km, 1125 km),
1302  ! conversion km -> m
1303 
1304 filename_with_path = trim(outpath)//'/'//trim(runname)//'.core'
1305 
1306 open(14, iostat=ios, file=trim(filename_with_path), status='new')
1307 
1308 if (ios /= 0) then
1309  errormsg = ' >>> sico_init: Error when opening the core file!'
1310  call error(errormsg)
1311 end if
1312 
1313 if (forcing_flag == 1) then
1314 
1315  write(14,1106)
1316  write(14,1107)
1317 
1318  1106 format(' t(a) D_Ts(C) z_sl(m)',/, &
1319  ' H_P1(m) H_P2(m)',/, &
1320  ' v_P1(m/a) v_P2(m/a)',/, &
1321  ' T_P1(C) T_P2(C)')
1322  1107 format('---------------------------------------')
1323 
1324 end if
1325 
1326 !-------- Output of the initial state --------
1327 
1328 #if (defined(OUTPUT_INIT))
1329 
1330 #if (OUTPUT_INIT==0)
1331  flag_init_output = .false.
1332 #elif (OUTPUT_INIT==1)
1333  flag_init_output = .true.
1334 #else
1335  errormsg = ' >>> sico_init: OUTPUT_INIT must be either 0 or 1!'
1336  call error(errormsg)
1337 #endif
1338 
1339 #else
1340  flag_init_output = .true. ! default for undefined parameter OUTPUT_INIT
1341 #endif
1342 
1343 #if (OUTPUT==1)
1344 
1345 #if (ERGDAT==0)
1346  flag_3d_output = .false.
1347 #elif (ERGDAT==1)
1348  flag_3d_output = .true.
1349 #else
1350  errormsg = ' >>> sico_init: ERGDAT must be either 0 or 1!'
1351  call error(errormsg)
1352 #endif
1353 
1354  if (flag_init_output) &
1355  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1356  flag_3d_output, ndat2d, ndat3d)
1357 
1358 #elif (OUTPUT==2)
1359 
1360 if (time_output(1) <= time_init+eps) then
1361 
1362 #if (ERGDAT==0)
1363  flag_3d_output = .false.
1364 #elif (ERGDAT==1)
1365  flag_3d_output = .true.
1366 #else
1367  errormsg = ' >>> sico_init: ERGDAT must be either 0 or 1!'
1368  call error(errormsg)
1369 #endif
1370 
1371  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1372  flag_3d_output, ndat2d, ndat3d)
1373 
1374 end if
1375 
1376 #elif (OUTPUT==3)
1377 
1378  flag_3d_output = .false.
1379 
1380  if (flag_init_output) &
1381  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1382  flag_3d_output, ndat2d, ndat3d)
1383 
1384 if (time_output(1) <= time_init+eps) then
1385 
1386  flag_3d_output = .true.
1387 
1388  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1389  flag_3d_output, ndat2d, ndat3d)
1390 
1391 end if
1392 
1393 #else
1394 
1395  errormsg = ' >>> sico_init: OUTPUT must be either 1, 2 or 3!'
1396  call error(errormsg)
1397 
1398 #endif
1399 
1400 if (flag_init_output) then
1401  call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1402  call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1403 end if
1404 
1405 end subroutine sico_init
1406 
1407 !-------------------------------------------------------------------------------
1408 !> Definition of the initial surface and bedrock topography
1409 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1410 !! For ice-free initial topography with relaxed lithosphere
1411 !! (not defined for this domain, thus routine stops execution of SICOPOLIS).
1412 !<------------------------------------------------------------------------------
1413 subroutine topography1(dxi, deta)
1414 
1415 implicit none
1416 
1417 real(dp), intent(out) :: dxi, deta
1418 
1419 ! integer(i4b) :: i, j
1420 ! integer(i4b) :: ios
1421 ! real(dp) :: xi0, eta0
1422 
1423 dxi=0.0_dp; deta=0.0_dp ! dummy values
1424 
1425 errormsg = ' >>> topography1: ANF_DAT==1 not defined for this domain!'
1426 call error(errormsg)
1427 
1428 end subroutine topography1
1429 
1430 !-------------------------------------------------------------------------------
1431 !> Definition of the initial surface and bedrock topography
1432 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1433 !! For ice-free initial topography with relaxed lithosphere.
1434 !<------------------------------------------------------------------------------
1435 subroutine topography2(dxi, deta)
1436 
1437 #if (GRID==0 || GRID==1)
1439 #endif
1440 
1441  use metric_m
1442  use topograd_m
1443 
1444 implicit none
1445 
1446 real(dp), intent(out) :: dxi, deta
1447 
1448 integer(i4b) :: i, j
1449 integer(i4b) :: ios
1450 real(dp) :: xi0, eta0
1451 
1452 !-------- Set topography --------
1453 
1454 zl0 = 0.0_dp
1455 
1456 maske = 1_i1b
1457 
1458 maske(0,:) = 2_i1b
1459 maske(jmax,:) = 2_i1b
1460 
1461 maske(:,0) = 2_i1b
1462 maske(:,imax) = 2_i1b
1463 
1464 !-------- Further stuff --------
1465 
1466 #if (GRID==0 || GRID==1)
1467 
1468 dxi = dx *1000.0_dp ! km -> m
1469 deta = dx *1000.0_dp ! km -> m
1470 
1471 xi0 = x0 *1000.0_dp ! km -> m
1472 eta0 = y0 *1000.0_dp ! km -> m
1473 
1474 #elif (GRID==2)
1475 
1476 errormsg = ' >>> topography2: GRID==2 not allowed for this application!'
1477 call error(errormsg)
1478 
1479 #endif
1480 
1481 do i=0, imax
1482 do j=0, jmax
1483 
1484  zs(j,i) = zl0(j,i)
1485  zb(j,i) = zl0(j,i)
1486  zl(j,i) = zl0(j,i)
1487 
1488  xi(i) = xi0 + real(i,dp)*dxi
1489  eta(j) = eta0 + real(j,dp)*deta
1490 
1491  zm(j,i) = zb(j,i)
1492  n_cts(j,i) = -1_i1b
1493  kc_cts(j,i) = 0
1494 
1495  h_c(j,i) = 0.0_dp
1496  h_t(j,i) = 0.0_dp
1497 
1498  dzs_dtau(j,i) = 0.0_dp
1499  dzm_dtau(j,i) = 0.0_dp
1500  dzb_dtau(j,i) = 0.0_dp
1501  dzl_dtau(j,i) = 0.0_dp
1502  dh_c_dtau(j,i) = 0.0_dp
1503  dh_t_dtau(j,i) = 0.0_dp
1504 
1505 end do
1506 end do
1507 
1508 maske_old = maske
1509 
1510 !-------- Geographic coordinates, metric tensor,
1511 ! gradients of the topography --------
1512 
1513 do i=0, imax
1514 do j=0, jmax
1515 
1516 #if (GRID==0 || GRID==1) /* Stereographic projection */
1517 
1518  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1519  lambda0, phi0, lambda(j,i), phi(j,i))
1520 
1521 #elif (GRID==2) /* Geographic coordinates */
1522 
1523  lambda(j,i) = xi(i)
1524  phi(j,i) = eta(j)
1525 
1526 #endif
1527 
1528 end do
1529 end do
1530 
1531 call metric()
1532 
1533 #if (TOPOGRAD==0)
1534 call topograd_1(dxi, deta, 1)
1535 #elif (TOPOGRAD==1)
1536 call topograd_2(dxi, deta, 1)
1537 #endif
1538 
1539 !-------- Corresponding area of grid points --------
1540 
1541 do i=0, imax
1542 do j=0, jmax
1543  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1544 end do
1545 end do
1546 
1547 end subroutine topography2
1548 
1549 !-------------------------------------------------------------------------------
1550 !> Definition of the initial surface and bedrock topography
1551 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1552 !! For initial topography from previous simulation.
1553 !<------------------------------------------------------------------------------
1554 subroutine topography3(dxi, deta, z_sl, anfdatname)
1555 
1556  use read_m, only : read_erg_nc
1557 
1558 #if (GRID==0 || GRID==1)
1560 #endif
1561 
1562  use metric_m
1563  use topograd_m
1564 
1565 implicit none
1566 
1567 character(len=100), intent(in) :: anfdatname
1568 
1569 real(dp), intent(out) :: dxi, deta, z_sl
1570 
1571 integer(i4b) :: i, j
1572 
1573 !-------- Read data from time-slice file of previous simulation --------
1574 
1575 call read_erg_nc(z_sl, anfdatname)
1576 
1577 !-------- Set topography of the relaxed bedrock --------
1578 
1579 zl0 = 0.0_dp
1580 
1581 !-------- Further stuff --------
1582 
1583 #if (GRID==0 || GRID==1)
1584 
1585 dxi = dx *1000.0_dp ! km -> m
1586 deta = dx *1000.0_dp ! km -> m
1587 
1588 #elif (GRID==2)
1589 
1590 errormsg = ' >>> topography3: GRID==2 not allowed for this application!'
1591 call error(errormsg)
1592 
1593 #endif
1594 
1595 !-------- Geographic coordinates, metric tensor,
1596 ! gradients of the topography --------
1597 
1598 do i=0, imax
1599 do j=0, jmax
1600 
1601 #if (GRID==0 || GRID==1) /* Stereographic projection */
1602 
1603  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1604  lambda0, phi0, lambda(j,i), phi(j,i))
1605 
1606 #elif (GRID==2) /* Geographic coordinates */
1607 
1608  lambda(j,i) = xi(i)
1609  phi(j,i) = eta(j)
1610 
1611 #endif
1612 
1613 end do
1614 end do
1615 
1616 call metric()
1617 
1618 #if (TOPOGRAD==0)
1619 call topograd_1(dxi, deta, 1)
1620 #elif (TOPOGRAD==1)
1621 call topograd_2(dxi, deta, 1)
1622 #endif
1623 
1624 !-------- Corresponding area of grid points --------
1625 
1626 do i=0, imax
1627 do j=0, jmax
1628  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1629 end do
1630 end do
1631 
1632 end subroutine topography3
1633 
1634 !-------------------------------------------------------------------------------
1635 
1636 end module sico_init_m
1637 !
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
logical, dimension(0:jmax, 0:imax) flag_inner_point
flag_inner_point(j,i): Inner-point flag. .true.: inner point, .false.: margin point ...
integer(i4b) grip_time_stp
grip_time_stp: Time step of the data values for the surface temperature anomaly
integer(i4b) ndata_specmap
ndata_specmap: Number of data values for the sea level
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
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) 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
integer(i4b) specmap_time_max
specmap_time_max: Maximum time of the data values for the sea level
Material properties of ice: Rate factor, heat conductivity, specific heat (heat capacity), creep function, viscosity.
real(dp) s_b
s_b: Gradient of accumulation rate change with horizontal distance
Definition: sico_vars_m.F90:54
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) 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(i4b) specmap_time_min
specmap_time_min: Minimum time of the data values for the sea level
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)
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.
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
integer(i4b) grip_time_max
grip_time_max: Maximum time of the data values for the surface temperature anomaly ...
real(dp), dimension(0:jmax, 0:imax) q_tld
Q_tld(j,i): Water drainage rate from the temperate layer.
real(dp) b_max
b_max: Maximum accumulation rate
Definition: sico_vars_m.F90:52
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars_m.F90:35
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.
integer(i4b) ndata_grip
ndata_grip: Number of data values for the surface temperature anomaly
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), 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), 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) year_zero
year_zero: SICOPOLIS year zero in astronomical year numbering [ = signed year CE (AD) ] ...
real(dp), dimension(:), allocatable specmap_zsl
specmap_zsl(n): Data values for the sea level
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) eld
eld: Equilibrium line distance
Definition: sico_vars_m.F90:56
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
Definition of the components g11 and g22 of the metric tensor of the applied numerical coordinates...
Definition: metric_m.F90:37
Computation of the forward or inverse stereographic projection, alternatively for a spherical or an e...
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) 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)
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) s_t
s_t: Gradient of surface temperature change with horizontal distance
Definition: sico_vars_m.F90:46
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)
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
integer(i4b) specmap_time_stp
specmap_time_stp: Time step of the data values for the sea level
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
subroutine, public calc_enhance_3(time)
Computation of the flow enhancement factor. Case ENHMOD==3: two different values depending on time of...
integer(i4b) grip_time_min
grip_time_min: Minimum time of the data values for the surface temperature anomaly ...
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) temp_min
temp_min: Minimum surface temperature
Definition: sico_vars_m.F90:44
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), dimension(0:ktmax) zeta_t
zeta_t(kt): Sigma coordinate zeta_t of grid point kt
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) y_hat
y_hat: Coordinate eta (= y) of the centre of the model domain
Definition: sico_vars_m.F90:50
character(len=256) errormsg
errormsg: Error-message string
real(dp) rho
RHO: Density of ice.
real(dp) x_hat
x_hat: Coordinate xi (= x) of the centre of the model domain
Definition: sico_vars_m.F90:48
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
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
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(:), allocatable griptemp
griptemp(n): Data values for the surface temperature anomaly
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...