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