SICOPOLIS V5-dev  Revision 1105
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-2017 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 
41  implicit none
42 
43  public
44 
45 contains
46 
47 !-------------------------------------------------------------------------------
48 !> Main routine of sico_init_m: Initialisations for SICOPOLIS.
49 !<------------------------------------------------------------------------------
50 subroutine sico_init(delta_ts, glac_index, &
51  mean_accum, &
52  dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
53  time, time_init, time_end, time_output, &
54  dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
55  z_sl, dzsl_dtau, z_mar, &
56  ii, jj, nn, &
57  ndat2d, ndat3d, n_output, &
58  runname)
59 
60  use compare_float_m
64 
65 #if (NETCDF > 1)
66  use netcdf
67  use nc_check_m
68 #endif
69 
70 #if (DISC>0 && !defined(EXEC_MAKE_C_DIS_0))
72 #endif
73 #if (DISC>0 && defined(EXEC_MAKE_C_DIS_0))
75 #endif
76 
77 #if (GRID==0 || GRID==1)
79 #endif
80 
82 
83  use boundary_m
85  use calc_enhance_m
86  use calc_vxy_m
87  use calc_vz_m
88  use calc_dxyz_m
90 
91  use output_m
92 
93 implicit none
94 
95 integer(i4b), intent(out) :: ii((imax+1)*(jmax+1)), &
96  jj((imax+1)*(jmax+1)), &
97  nn(0:jmax,0:imax)
98 integer(i4b), intent(out) :: ndat2d, ndat3d
99 integer(i4b), intent(out) :: n_output
100 real(dp), intent(out) :: delta_ts, glac_index
101 real(dp), intent(out) :: mean_accum
102 real(dp), intent(out) :: dtime, dtime_temp, dtime_wss, &
103  dtime_out, dtime_ser
104 real(dp), intent(out) :: time, time_init, time_end, time_output(100)
105 real(dp), intent(out) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
106 real(dp), intent(out) :: z_sl, dzsl_dtau, z_mar
107 
108 character(len=100), intent(out) :: runname
109 
110 integer(i4b) :: i, j, kc, kt, kr, m, n, ir, jr
111 integer(i4b) :: ios, ios1, ios2, ios3, ios4
112 integer(i4b) :: ierr
113 real(dp) :: dtime0, dtime_temp0, dtime_wss0, dtime_out0, dtime_ser0
114 real(dp) :: time_init0, time_end0, time_output0(100)
115 real(dp), dimension(0:JMAX,0:IMAX) :: precip_factor
116 real(dp) :: d_dummy
117 character(len=100) :: anfdatname, target_topo_dat_name
118 character(len=256) :: filename_with_path
119 character(len=256) :: shell_command
120 character :: ch_dummy
121 logical :: flag_precip_monthly_mean
122 logical :: flag_init_output, flag_3d_output
123 
124 #if (defined(INITMIP_SMB_ANOM_FILE))
125 real(dp), dimension(0:IMAX,0:JMAX) :: as_anom_initmip_conv
126 #endif
127 
128 #if (NETCDF > 1)
129 integer(i4b) :: dimid, ncid, ncv
130 ! dimid: Dimension ID
131 ! ncid: File ID
132 ! ncv: Variable ID
133 #endif
134 
135 real(dp), parameter :: inv_twelve = 1.0_dp/12.0_dp
136 
137 character(len=64), parameter :: thisroutine = 'sico_init'
138 
139 character(len=64), parameter :: fmt1 = '(a)', &
140  fmt2 = '(a,i4)', &
141  fmt2a = '(a,i0)', &
142  fmt3 = '(a,es12.4)'
143 
144 character(len= 8) :: ch_imax
145 character(len=128) :: fmt4
146 
147 write(ch_imax, fmt='(i8)') imax
148 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
149 
150 write(unit=6, fmt='(a)') ' '
151 write(unit=6, fmt='(a)') ' -------- sico_init --------'
152 write(unit=6, fmt='(a)') ' '
153 
154 !-------- Name of the computational domain --------
155 
156 #if (defined(ANT))
157 ch_domain_long = 'Antarctica'
158 ch_domain_short = 'ant'
159 
160 #elif (defined(ASF))
161 ch_domain_long = 'Austfonna'
162 ch_domain_short = 'asf'
163 
164 #elif (defined(EMTP2SGE))
165 ch_domain_long = 'EISMINT Phase 2 Simplified Geometry Experiment'
166 ch_domain_short = 'emtp2sge'
167 
168 #elif (defined(GRL))
169 ch_domain_long = 'Greenland'
170 ch_domain_short = 'grl'
171 
172 #elif (defined(NHEM))
173 ch_domain_long = 'Northern hemisphere'
174 ch_domain_short = 'nhem'
175 
176 #elif (defined(SCAND))
177 ch_domain_long = 'Scandinavia and Eurasia'
178 ch_domain_short = 'scand'
179 
180 #elif (defined(TIBET))
181 ch_domain_long = 'Tibet'
182 ch_domain_short = 'tibet'
183 
184 #elif (defined(NMARS))
185 ch_domain_long = 'North polar cap of Mars'
186 ch_domain_short = 'nmars'
187 
188 #elif (defined(SMARS))
189 ch_domain_long = 'South polar cap of Mars'
190 ch_domain_short = 'smars'
191 
192 #elif (defined(XYZ))
193 ch_domain_long = 'XYZ'
194 ch_domain_short = 'xyz'
195 #if (defined(HEINO))
196 ch_domain_long = trim(ch_domain_long)//'/ISMIP HEINO'
197 #endif
198 
199 #else
200 stop ' >>> sico_init: No valid domain specified!'
201 #endif
202 
203 !-------- Some initial values --------
204 
205 n_output = 0
206 
207 dtime = 0.0_dp
208 dtime_temp = 0.0_dp
209 dtime_wss = 0.0_dp
210 dtime_out = 0.0_dp
211 dtime_ser = 0.0_dp
212 
213 time = 0.0_dp
214 time_init = 0.0_dp
215 time_end = 0.0_dp
216 time_output = 0.0_dp
217 
218 !-------- Initialisation of the Library of Iterative Solvers Lis,
219 ! if required --------
220 
221 #if (CALCTHK==3 || CALCTHK==6 || MARGIN==3 || DYNAMICS==2)
222  call lis_initialize(ierr)
223 #endif
224 
225 !-------- Read physical parameters --------
226 
227 call read_phys_para()
228 
229 call ice_mat_eqs_pars(rf, r_t, kappa, c, -190, 10)
230 
231 ! ------ Some auxiliary quantities required for the enthalpy method
232 
233 call calc_c_int_table(c, -190, 10, l)
235 
236 !-------- Compatibility check of the SICOPOLIS version with the header file
237 
238 if ( trim(version) /= trim(sico_version) ) &
239  stop ' >>> sico_init: SICOPOLIS version not compatible with header file!'
240 
241 !-------- Check whether the dynamics and thermodynamics modes are defined
242 
243 #if (!defined(DYNAMICS))
244 stop ' >>> sico_init: DYNAMICS not defined in the header file!'
245 #endif
246 
247 #if (!defined(CALCMOD))
248 stop ' >>> sico_init: CALCMOD not defined in the header file!'
249 #endif
250 
251 #if (defined(ENTHMOD))
252 write(6, fmt='(a)') ' >>> sico_init: ENTHMOD must not be defined any more.'
253 write(6, fmt='(a)') ' Please update your header file!'
254 stop
255 #endif
256 
257 !-------- Compatibility check of the horizontal resolution with the
258 ! number of grid points --------
259 
260 #if (GRID==0 || GRID==1)
261 
262 if (approx_equal(dx, 20.0_dp, eps_sp_dp)) then
263 
264  if ( (imax /= 75).or.(jmax /= 140) ) then
265  stop ' >>> sico_init: IMAX and/or JMAX wrong!'
266  end if
267 
268 else if (approx_equal(dx, 10.0_dp, eps_sp_dp)) then
269 
270  if ( (imax /= 150).or.(jmax /= 280) ) then
271  stop ' >>> sico_init: IMAX and/or JMAX wrong!'
272  end if
273 
274 else if (approx_equal(dx, 5.0_dp, eps_sp_dp)) then
275 
276  if ( (imax /= 300).or.(jmax /= 560) ) then
277  stop ' >>> sico_init: IMAX and/or JMAX wrong!'
278  end if
279 
280 else
281  stop ' >>> sico_init: DX wrong!'
282 end if
283 
284 #elif (GRID==2)
285 
286 stop ' >>> sico_init: GRID==2 not allowed for this application!'
287 
288 #endif
289 
290 !-------- Compatibility check of the thermodynamics mode
291 ! (cold vs. polythermal vs. enthalpy method)
292 ! and the number of grid points in the lower (kt) ice domain --------
293 
294 #if (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1)
295 
296 if (ktmax > 2) then
297  write(6, fmt='(a)') ' >>> sico_init: For options CALCMOD==0, 2, 3 or -1,'
298  write(6, fmt='(a)') ' the separate kt domain is redundant.'
299  write(6, fmt='(a)') ' Therefore, consider setting KTMAX to 2.'
300  write(6, fmt='(a)') ' '
301 end if
302 
303 #endif
304 
305 !-------- Compatibility check of surface-temperature
306 ! and precipitation switches --------
307 
308 #if (TSURFACE == 5 && ACCSURFACE != 5)
309 stop ' >>> sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
310 #endif
311 
312 #if (TSURFACE != 5 && ACCSURFACE == 5)
313 stop ' >>> sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
314 #endif
315 
316 #if (TSURFACE == 6)
317 #if (ACCSURFACE != 6)
318 write(6, fmt='(a)') ' >>> sico_init: Options TSURFACE==6, ACCSURFACE==6'
319 write(6, fmt='(a)') ' and NETCDF==2 must be used together!'
320 stop
321 #elif (NETCDF != 2)
322 write(6, fmt='(a)') ' >>> sico_init: Options TSURFACE==6, ACCSURFACE==6'
323 write(6, fmt='(a)') ' and NETCDF==2 must be used together!'
324 stop
325 #endif
326 #endif
327 
328 #if (ACCSURFACE == 6)
329 #if (TSURFACE != 6)
330 write(6, fmt='(a)') ' >>> sico_init: Options TSURFACE==6, ACCSURFACE==6'
331 write(6, fmt='(a)') ' and NETCDF==2 must be used together!'
332 stop
333 #elif (NETCDF != 2)
334 write(6, fmt='(a)') ' >>> sico_init: Options TSURFACE==6, ACCSURFACE==6'
335 write(6, fmt='(a)') ' and NETCDF==2 must be used together!'
336 stop
337 #endif
338 #endif
339 
340 #if (defined(INITMIP_SMB_ANOM_FILE) && NETCDF != 2)
341 if ( trim(adjustl(initmip_smb_anom_file)) /= 'none' ) then
342  write(6, fmt='(a)') ' >>> sico_init: Defining INITMIP_SMB_ANOM_FILE'
343  write(6, fmt='(a)') ' must be used together with NETCDF==2!'
344  stop
345 end if
346 #endif
347 
348 !-------- Compatibility check of discretization schemes for the horizontal and
349 ! vertical advection terms in the temperature and age equations --------
350 
351 #if (ADV_HOR==1)
352 stop ' >>> sico_init: Option ADV_HOR==1 (central differences) not defined!'
353 #endif
354 
355 !-------- Check whether for the shallow shelf
356 ! or shelfy stream approximation
357 ! the chosen grid is Cartesian coordinates
358 ! without distortion correction (GRID==0) --------
359 
360 #if ((MARGIN==3 && DYNAMICS==1) || DYNAMICS==2) /* requires SSA or SStA */
361 #if (GRID != 0)
362 write(6, fmt='(a)') ' >>> sico_init: Distortion correction not implemented'
363 write(6, fmt='(a)') ' for the shallow shelf approximation (SSA)'
364 write(6, fmt='(a)') ' or the shelfy stream approximation (SStA)'
365 write(6, fmt='(a)') ' -> GRID==0 required!'
366 stop
367 #endif
368 #endif
369 
370 
371 !-------- Check whether for the hidden ablation scheme
372 ! can be applied --------
373 
374 #if (MARGIN==3 && MB_ACCOUNT==1)
375 write(6, fmt='(a)') &
376  ' >>> sico_init: MB_ACCOUNT==1 does not yet work for ice shelves!'
377 stop
378 #endif
379 
380 !-------- Setting of forcing flag --------
381 
382 #if (TSURFACE <= 4)
383 
384 forcing_flag = 1 ! forcing by delta_ts
385 
386 #elif (TSURFACE == 5)
387 
388 forcing_flag = 2 ! forcing by glac_index
389 
390 #elif (TSURFACE == 6)
391 
392 forcing_flag = 3 ! forcing by time-dependent surface temperature
393  ! and precipitation data
394 
395 #endif
396 
397 !-------- Initialization of numerical time steps --------
398 
399 dtime0 = dtime0
400 dtime_temp0 = dtime_temp0
401 #if (REBOUND==2)
402 dtime_wss0 = dtime_wss0
403 #endif
404 
405 !-------- Further initializations --------
406 
407 dzeta_c = 1.0_dp/real(kcmax,dp)
408 dzeta_t = 1.0_dp/real(ktmax,dp)
409 dzeta_r = 1.0_dp/real(krmax,dp)
410 
411 ndat2d = 1
412 ndat3d = 1
413 
414 !-------- General abbreviations --------
415 
416 ! ------ kc domain
417 
418 if (deform >= eps) then
419 
420  flag_aa_nonzero = .true. ! non-equidistant grid
421 
422  aa = deform
423  ea = exp(aa)
424 
425  kc=0
426  zeta_c(kc) = 0.0_dp
427  eaz_c(kc) = 1.0_dp
428  eaz_c_quotient(kc) = 0.0_dp
429 
430  do kc=1, kcmax-1
431  zeta_c(kc) = kc*dzeta_c
432  eaz_c(kc) = exp(aa*zeta_c(kc))
433  eaz_c_quotient(kc) = (eaz_c(kc)-1.0_dp)/(ea-1.0_dp)
434  end do
435 
436  kc=kcmax
437  zeta_c(kc) = 1.0_dp
438  eaz_c(kc) = exp(aa)
439  eaz_c_quotient(kc) = 1.0_dp
440 
441 else
442 
443  flag_aa_nonzero = .false. ! equidistant grid
444 
445  aa = 0.0_dp
446  ea = 1.0_dp
447 
448  kc=0
449  zeta_c(kc) = 0.0_dp
450  eaz_c(kc) = 1.0_dp
451  eaz_c_quotient(kc) = 0.0_dp
452 
453  do kc=1, kcmax-1
454  zeta_c(kc) = kc*dzeta_c
455  eaz_c(kc) = 1.0_dp
456  eaz_c_quotient(kc) = zeta_c(kc)
457  end do
458 
459  kc=kcmax
460  zeta_c(kc) = 1.0_dp
461  eaz_c(kc) = 1.0_dp
462  eaz_c_quotient(kc) = 1.0_dp
463 
464 end if
465 
466 ! ------ kt domain
467 
468 kt=0
469 zeta_t(kt) = 0.0_dp
470 
471 do kt=1, ktmax-1
472  zeta_t(kt) = kt*dzeta_t
473 end do
474 
475 kt=ktmax
476 zeta_t(kt) = 1.0_dp
477 
478 ! ------ kr domain
479 
480 kr=0
481 zeta_r(kr) = 0.0_dp
482 
483 do kr=1, krmax-1
484  zeta_r(kr) = kr*dzeta_r
485 end do
486 
487 kr=krmax
488 zeta_r(kr) = 1.0_dp
489 
490 !-------- Construction of a vector (with index n) from a 2-d array
491 ! (with indices i, j) by diagonal numbering --------
492 
493 n=1
494 
495 do m=0, imax+jmax
496  do i=m, 0, -1
497  j = m-i
498  if ((i <= imax).and.(j <= jmax)) then
499  ii(n) = i
500  jj(n) = j
501  nn(j,i) = n
502  n=n+1
503  end if
504  end do
505 end do
506 
507 !-------- Specification of current simulation --------
508 
509 runname = runname
510 anfdatname = anfdatname
511 
512 #if (defined(YEAR_ZERO))
514 #else
515 year_zero = 2000.0_dp ! default value 2000 CE
516 #endif
517 
518 time_init0 = time_init0
519 time_end0 = time_end0
520 dtime_ser0 = dtime_ser0
521 
522 #if (OUTPUT==1 || OUTPUT==3)
523 dtime_out0 = dtime_out0
524 #endif
525 
526 #if (OUTPUT==2 || OUTPUT==3)
527 n_output = n_output
528 time_output0( 1) = time_out0_01
529 time_output0( 2) = time_out0_02
530 time_output0( 3) = time_out0_03
531 time_output0( 4) = time_out0_04
532 time_output0( 5) = time_out0_05
533 time_output0( 6) = time_out0_06
534 time_output0( 7) = time_out0_07
535 time_output0( 8) = time_out0_08
536 time_output0( 9) = time_out0_09
537 time_output0(10) = time_out0_10
538 time_output0(11) = time_out0_11
539 time_output0(12) = time_out0_12
540 time_output0(13) = time_out0_13
541 time_output0(14) = time_out0_14
542 time_output0(15) = time_out0_15
543 time_output0(16) = time_out0_16
544 time_output0(17) = time_out0_17
545 time_output0(18) = time_out0_18
546 time_output0(19) = time_out0_19
547 time_output0(20) = time_out0_20
548 #endif
549 
550 !-------- Write log file --------
551 
552 shell_command = 'if [ ! -d'
553 shell_command = trim(shell_command)//' '//outpath
554 shell_command = trim(shell_command)//' '//'] ; then mkdir'
555 shell_command = trim(shell_command)//' '//outpath
556 shell_command = trim(shell_command)//' '//'; fi'
557 call system(trim(shell_command))
558  ! Check whether directory OUTPATH exists. If not, it is created.
559 
560 filename_with_path = trim(outpath)//'/'//trim(runname)//'.log'
561 
562 open(10, iostat=ios, file=trim(filename_with_path), status='new')
563 
564 if (ios /= 0) stop ' >>> sico_init: Error when opening the log file!'
565 
566 write(10, fmt=trim(fmt1)) 'Computational domain:'
567 write(10, fmt=trim(fmt1)) trim(ch_domain_long)
568 write(10, fmt=trim(fmt1)) ' '
569 
570 write(10, fmt=trim(fmt2a)) 'GRID = ', grid
571 write(10, fmt=trim(fmt1)) ' '
572 
573 write(10, fmt=trim(fmt2)) 'imax =', imax
574 write(10, fmt=trim(fmt2)) 'jmax =', jmax
575 write(10, fmt=trim(fmt2)) 'kcmax =', kcmax
576 write(10, fmt=trim(fmt2)) 'ktmax =', ktmax
577 write(10, fmt=trim(fmt2)) 'krmax =', krmax
578 write(10, fmt=trim(fmt1)) ' '
579 
580 write(10, fmt=trim(fmt3)) 'a =', aa
581 write(10, fmt=trim(fmt1)) ' '
582 
583 #if (GRID==0 || GRID==1)
584 write(10, fmt=trim(fmt3)) 'x0 =', x0
585 write(10, fmt=trim(fmt3)) 'y0 =', y0
586 write(10, fmt=trim(fmt3)) 'dx =', dx
587 #elif GRID==2
588 stop ' >>> sico_init: GRID==2 not allowed for this application!'
589 #endif
590 write(10, fmt=trim(fmt1)) ' '
591 
592 write(10, fmt=trim(fmt3)) 'year_zero =', year_zero
593 write(10, fmt=trim(fmt3)) 'time_init =', time_init0
594 write(10, fmt=trim(fmt3)) 'time_end =', time_end0
595 write(10, fmt=trim(fmt3)) 'dtime =', dtime0
596 write(10, fmt=trim(fmt3)) 'dtime_temp =', dtime_temp0
597 #if (REBOUND==2)
598 write(10, fmt=trim(fmt3)) 'dtime_wss =', dtime_wss0
599 #endif
600 #if (DISC>0)
601 write(10, fmt=trim(fmt3)) 'dtime_mar_coa =', dtime_mar_coa0
602 #endif
603 write(10, fmt=trim(fmt1)) ' '
604 
605 write(10, fmt=trim(fmt2)) 'ANF_DAT =', anf_dat
606 write(10, fmt=trim(fmt1)) 'zs_present file = '//zs_present_file
607 #if (ANF_DAT==1)
608 #if (defined(ZB_PRESENT_FILE))
609 write(10, fmt=trim(fmt1)) 'zb_present file = '//zb_present_file
610 #endif
611 write(10, fmt=trim(fmt1)) 'zl_present file = '//zl_present_file
612 #endif
613 write(10, fmt=trim(fmt1)) 'zl0 file = '//zl0_file
614 write(10, fmt=trim(fmt1)) 'mask_present file = '//mask_present_file
615 #if (ANF_DAT==1 && defined(TEMP_INIT))
616 write(10, fmt=trim(fmt2)) 'TEMP_INIT =', temp_init
617 #endif
618 #if (ANF_DAT==3)
619 write(10, fmt=trim(fmt1)) 'Initial-value file = '//anfdatname
620 write(10, fmt=trim(fmt1)) 'Path to initial-value file = '//anfdatpath
621 #endif
622 write(10, fmt=trim(fmt1)) ' '
623 
624 
625 #if (THK_EVOL==2)
626 write(10, fmt=trim(fmt3)) 'time_target_topo_init =', time_target_topo_init0
627 write(10, fmt=trim(fmt3)) 'time_target_topo_final =', time_target_topo_final0
628 write(10, fmt=trim(fmt1)) 'Target-topography file = '//target_topo_dat_name
629 write(10, fmt=trim(fmt1)) 'Path to target-topography file = '//target_topo_dat_path
630 write(10, fmt=trim(fmt1)) ' '
631 #endif
632 
633 #if (THK_EVOL==3)
634 write(10, fmt=trim(fmt3)) 'target_topo_tau =', target_topo_tau0
635 write(10, fmt=trim(fmt1)) 'Target-topography file = '//target_topo_dat_name
636 write(10, fmt=trim(fmt1)) 'Path to target-topography file = '//target_topo_dat_path
637 write(10, fmt=trim(fmt1)) ' '
638 #endif
639 
640 #if (THK_EVOL==4)
641 write(10, fmt=trim(fmt1)) 'Maximum ice extent mask file = '//mask_maxextent_file
642 write(10, fmt=trim(fmt1)) ' '
643 #endif
644 
645 write(10, fmt=trim(fmt1)) 'Physical-parameter file = '//phys_para_file
646 write(10, fmt=trim(fmt1)) ' '
647 
648 #if (CALCTHK==2 || CALCTHK==3 || CALCTHK==5 || CALCTHK==6)
649 write(10, fmt=trim(fmt3)) 'ovi_weight =', ovi_weight
650 #if (CALCTHK==2 || CALCTHK==5)
651 write(10, fmt=trim(fmt3)) 'omega_sor =', omega_sor
652 #if (ITER_MAX_SOR>0)
653 write(10, fmt=trim(fmt2a)) 'iter_max_sor = ', iter_max_sor
654 #endif
655 #endif
656 write(10, fmt=trim(fmt1)) ' '
657 #endif
658 
659 #if (TSURFACE==1)
660 write(10, fmt=trim(fmt3)) 'delta_ts0 =', delta_ts0
661 #elif (TSURFACE==3)
662 write(10, fmt=trim(fmt3)) 'sine_amplit =', sine_amplit
663 write(10, fmt=trim(fmt3)) 'sine_period =', sine_period
664 #elif (TSURFACE==4)
665 write(10, fmt=trim(fmt1)) 'GRIP file = '//grip_temp_file
666 write(10, fmt=trim(fmt3)) 'grip_temp_fact =', grip_temp_fact
667 #elif (TSURFACE==5)
668 write(10, fmt=trim(fmt1)) 'Glacial-index file = '//glac_ind_file
669 write(10, fmt=trim(fmt1)) 'temp_ma_anom file = '//temp_ma_anom_file
670 write(10, fmt=trim(fmt3)) 'temp_ma_anom fact = ', temp_ma_anom_fact
671 write(10, fmt=trim(fmt1)) 'temp_mj_anom file = '//temp_mj_anom_file
672 write(10, fmt=trim(fmt3)) 'temp_mj_anom fact = ', temp_mj_anom_fact
673 #elif (TSURFACE==6)
674 write(10, fmt=trim(fmt1)) 'temp_precip_anom file = '//temp_precip_anom_file
675 write(10, fmt=trim(fmt3)) 'temp_ma_anom fact = ', temp_ma_anom_fact
676 write(10, fmt=trim(fmt3)) 'temp_mj_anom fact = ', temp_mj_anom_fact
677 #endif
678 
679 #if (!defined(PRECIP_PRESENT_FILE))
680 stop ' >>> sico_init: PRECIP_PRESENT_FILE not defined in the header file!'
681 #endif
682 #if (!defined(PRECIP_MA_PRESENT_FILE))
683 stop ' >>> sico_init: PRECIP_MA_PRESENT_FILE not defined in the header file!'
684 #endif
685 #if (defined(PRECIP_PRESENT_FILE) && defined(PRECIP_MA_PRESENT_FILE))
686 if ( trim(adjustl(precip_present_file)) /= 'none' ) then
687  write(10, fmt=trim(fmt1)) 'precip_present_file = '//precip_present_file
688  flag_precip_monthly_mean = .true.
689 else if ( trim(adjustl(precip_ma_present_file)) /= 'none' ) then
690  write(10, fmt=trim(fmt1)) 'precip_ma_present_file = '//precip_ma_present_file
691  flag_precip_monthly_mean = .false.
692 else
693  write(6, fmt='(a)') ' >>> sico_init: Neither PRECIP_PRESENT_FILE'
694  write(6, fmt='(a)') ' nor PRECIP_MA_PRESENT_FILE'
695  write(6, fmt='(a)') ' specified in the header file!'
696  stop
697 end if
698 #endif
699 #if (!defined(PRECIP_ZS_REF_FILE))
700 stop ' >>> sico_init: PRECIP_ZS_REF_FILE not defined in the header file!'
701 #else
702 write(10, fmt=trim(fmt1)) 'precip_zs_ref_file = '//precip_zs_ref_file
703 #endif
704 #if (defined(PRECIP_FACTOR_FILE))
705 if ( trim(adjustl(precip_factor_file)) /= 'none' ) &
706  write(10, fmt=trim(fmt1)) 'precip_factor_file = '//precip_factor_file
707 #endif
708 #if (ACCSURFACE==1)
709 write(10, fmt=trim(fmt3)) 'accfact =', accfact
710 #elif (ACCSURFACE==2 || ACCSURFACE==3)
711 write(10, fmt=trim(fmt3)) 'gamma_s =', gamma_s
712 #elif (ACCSURFACE==5)
713 write(10, fmt=trim(fmt1)) 'precip_anom file = '//precip_anom_file
714 write(10, fmt=trim(fmt3)) 'precip_anom fact = ', precip_anom_fact
715 #elif (ACCSURFACE==6)
716 write(10, fmt=trim(fmt1)) 'temp_precip_anom file = '//temp_precip_anom_file
717 write(10, fmt=trim(fmt3)) 'precip_anom fact = ', precip_anom_fact
718 #endif
719 #if (ACCSURFACE <= 3)
720 write(10, fmt=trim(fmt2)) 'ELEV_DESERT =', elev_desert
721 #if (ELEV_DESERT == 1)
722 write(10, fmt=trim(fmt3)) 'gamma_p =', gamma_p
723 write(10, fmt=trim(fmt3)) 'zs_thresh =', zs_thresh
724 #endif
725 #endif
726 
727 #if (ABLSURFACE==3)
728 write(10, fmt=trim(fmt3)) 'lambda_lti =', lambda_lti
729 write(10, fmt=trim(fmt3)) 'temp_lti =', temp_lti
730 #endif
731 
732 write(10, fmt=trim(fmt1)) ' '
733 
734 #if (defined(INITMIP_CONST_SMB))
735 write(10, fmt=trim(fmt2a)) 'INITMIP_CONST_SMB = ', initmip_const_smb
736 #endif
737 #if (defined(INITMIP_SMB_ANOM_FILE))
738 if ( trim(adjustl(initmip_smb_anom_file)) /= 'none' ) then
739  write(10, fmt=trim(fmt1)) 'initmip_smb_anom file = '//initmip_smb_anom_file
740 end if
741 #endif
742 #if (defined(INITMIP_CONST_SMB) || defined(INITMIP_SMB_ANOM_FILE))
743 write(10, fmt=trim(fmt1)) ' '
744 #endif
745 
746 write(10, fmt=trim(fmt2a)) 'DISC = ', disc
747 #if (DISC>0)
748 write(10, fmt=trim(fmt3)) 'c_dis_0 =', c_dis_0
749 write(10, fmt=trim(fmt3)) 'c_dis_fac =', c_dis_fac
750 write(10, fmt=trim(fmt3)) 'm_H =', m_h
751 write(10, fmt=trim(fmt3)) 'm_D =', m_d
752 write(10, fmt=trim(fmt3)) 'r_mar_eff =', r_mar_eff
753 #if (defined(S_DIS))
754 write(10, fmt=trim(fmt3)) 's_dis =', s_dis
755 #endif
756 #if (defined(ALPHA_SUB))
757 write(10, fmt=trim(fmt3)) 'alpha_sub =', alpha_sub
758 #endif
759 #if (defined(ALPHA_O))
760 write(10, fmt=trim(fmt3)) 'alpha_o =', alpha_o
761 #endif
762 #endif
763 write(10, fmt=trim(fmt1)) ' '
764 
765 write(10, fmt=trim(fmt2a)) 'SEA_LEVEL = ', sea_level
766 #if (SEA_LEVEL==1)
767 write(10, fmt=trim(fmt3)) 'z_sl0 =', z_sl0
768 #elif (SEA_LEVEL==3)
769 write(10, fmt=trim(fmt1)) 'sea-level file = '//sea_level_file
770 #endif
771 write(10, fmt=trim(fmt1)) ' '
772 
773 #if (MARGIN==2)
774 #if (MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3)
775 write(10, fmt=trim(fmt3)) 'z_mar =', z_mar
776 write(10, fmt=trim(fmt1)) ' '
777 #elif (MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \
778  || marine_ice_calving==6 || marine_ice_calving==7)
779 write(10, fmt=trim(fmt3)) 'fact_z_mar =', fact_z_mar
780 write(10, fmt=trim(fmt1)) ' '
781 #elif (MARINE_ICE_FORMATION==2 && MARINE_ICE_CALVING==9)
782 write(10, fmt=trim(fmt3)) 'calv_uw_coeff =', calv_uw_coeff
783 write(10, fmt=trim(fmt3)) 'r1_calv_uw =', r1_calv_uw
784 write(10, fmt=trim(fmt3)) 'r2_calv_uw =', r2_calv_uw
785 write(10, fmt=trim(fmt1)) ' '
786 #endif
787 #elif (MARGIN==3)
788 #if (ICE_SHELF_CALVING==2)
789 write(10, fmt=trim(fmt3)) 'H_calv =', h_calv
790 write(10, fmt=trim(fmt1)) ' '
791 #endif
792 #endif
793 
794 #if (defined(BASAL_HYDROLOGY))
795 write(10, fmt=trim(fmt2a)) 'BASAL_HYDROLOGY = ', basal_hydrology
796 #endif
797 write(10, fmt=trim(fmt2a)) 'SLIDE_LAW = ', slide_law
798 write(10, fmt=trim(fmt3)) 'c_slide =', c_slide
799 write(10, fmt=trim(fmt3)) 'smw_coeff =', smw_coeff
800 write(10, fmt=trim(fmt2a)) 'r_smw = ', r_smw
801 write(10, fmt=trim(fmt2a)) 's_smw = ', s_smw
802 write(10, fmt=trim(fmt3)) 'gamma_slide =', gamma_slide
803 write(10, fmt=trim(fmt2a)) 'p_weert = ', p_weert
804 write(10, fmt=trim(fmt2a)) 'q_weert = ', q_weert
805 #if (defined(TIME_RAMP_UP_SLIDE))
806 write(10, fmt=trim(fmt3)) 'time_ramp_up_slide =', time_ramp_up_slide
807 #endif
808 #if (SLIDE_LAW==2)
809 write(10, fmt=trim(fmt3)) 'red_pres_limit_fact =', red_pres_limit_fact
810 #endif
811 #if (BASAL_HYDROLOGY==1)
812 write(10, fmt=trim(fmt3)) 'c_Hw_slide =', c_hw_slide
813 write(10, fmt=trim(fmt3)) 'Hw0_slide =', hw0_slide
814 #endif
815 write(10, fmt=trim(fmt2a)) 'ICE_STREAM = ', ice_stream
816 #if (ICE_STREAM==2)
817 write(10, fmt=trim(fmt1)) 'Sediment-mask file = '//mask_sedi_file
818 write(10, fmt=trim(fmt3)) 'c_slide_sedi =', c_slide_sedi
819 write(10, fmt=trim(fmt3)) 'smw_coeff_sedi =', smw_coeff_sedi
820 write(10, fmt=trim(fmt2a)) 'r_smw_sedi = ', r_smw_sedi
821 write(10, fmt=trim(fmt2a)) 's_smw_sedi = ', s_smw_sedi
822 write(10, fmt=trim(fmt3)) 'gamma_slide_sedi =', gamma_slide_sedi
823 write(10, fmt=trim(fmt2a)) 'p_weert_sedi = ', p_weert_sedi
824 write(10, fmt=trim(fmt2a)) 'q_weert_sedi = ', q_weert_sedi
825 #endif
826 write(10, fmt=trim(fmt1)) ' '
827 
828 write(10, fmt=trim(fmt2a)) 'Q_GEO_MOD = ', q_geo_mod
829 #if (Q_GEO_MOD==1)
830 write(10, fmt=trim(fmt3)) 'q_geo =', q_geo
831 #elif (Q_GEO_MOD==2)
832 write(10, fmt=trim(fmt1)) 'q_geo file = '//q_geo_file
833 #endif
834 write(10, fmt=trim(fmt1)) ' '
835 
836 #if (defined(MARINE_ICE_BASAL_MELTING))
837 write(10, fmt=trim(fmt2)) 'MARINE_ICE_BASAL_MELTING =', marine_ice_basal_melting
838 #if (MARINE_ICE_BASAL_MELTING==2 || MARINE_ICE_BASAL_MELTING==3)
839 write(10, fmt=trim(fmt3)) 'qbm_marine =', qbm_marine
840 #endif
841 write(10, fmt=trim(fmt1)) ' '
842 #endif
843 
844 #if (MARGIN==3)
845 write(10, fmt=trim(fmt2)) 'FLOATING_ICE_BASAL_MELTING =', floating_ice_basal_melting
846 #if (FLOATING_ICE_BASAL_MELTING<=3)
847 write(10, fmt=trim(fmt3)) 'qbm_float_1 =', qbm_float_1
848 write(10, fmt=trim(fmt3)) 'qbm_float_2 =', qbm_float_2
849 #endif
850 write(10, fmt=trim(fmt3)) 'qbm_float_3 =', qbm_float_3
851 write(10, fmt=trim(fmt3)) 'z_abyss =', z_abyss
852 #if (FLOATING_ICE_BASAL_MELTING==4)
853 write(10, fmt=trim(fmt3)) 'temp_ocean =', temp_ocean
854 write(10, fmt=trim(fmt3)) 'Omega_qbm =', omega_qbm
855 write(10, fmt=trim(fmt3)) 'alpha_qbm =', alpha_qbm
856 #endif
857 #if (FLOATING_ICE_BASAL_MELTING==4 || FLOATING_ICE_BASAL_MELTING==5)
858 write(10, fmt=trim(fmt3)) 'H_w_0 =', h_w_0
859 #endif
860 write(10, fmt=trim(fmt1)) ' '
861 #endif
862 
863 write(10, fmt=trim(fmt2)) 'REBOUND =', rebound
864 #if (REBOUND==1)
865 write(10, fmt=trim(fmt3)) 'frac_llra =', frac_llra
866 #endif
867 #if (REBOUND==1 || REBOUND==2)
868 write(10, fmt=trim(fmt2)) 'TIME_LAG_MOD =', time_lag_mod
869 #if (TIME_LAG_MOD==1)
870 write(10, fmt=trim(fmt3)) 'time_lag =', time_lag
871 #elif (TIME_LAG_MOD==2)
872 write(10, fmt=trim(fmt1)) 'time_lag_file = '//time_lag_file
873 #else
874 stop ' >>> sico_init: TIME_LAG_MOD must be either 1 or 2!'
875 #endif
876 #endif
877 #if (REBOUND==2)
878 write(10, fmt=trim(fmt2)) 'FLEX_RIG_MOD =', flex_rig_mod
879 #if (FLEX_RIG_MOD==1)
880 write(10, fmt=trim(fmt3)) 'flex_rig =', flex_rig
881 #elif (FLEX_RIG_MOD==2)
882 write(10, fmt=trim(fmt1)) 'flex_rig_file = '//flex_rig_file
883 #else
884 stop ' >>> sico_init: FLEX_RIG_MOD must be either 1 or 2!'
885 #endif
886 #endif
887 write(10, fmt=trim(fmt2)) 'Q_LITHO =', q_litho
888 write(10, fmt=trim(fmt1)) ' '
889 
890 #if (FLOW_LAW==2)
891 write(10, fmt=trim(fmt3)) 'gr_size =', gr_size
892 write(10, fmt=trim(fmt1)) ' '
893 #endif
894 #if (FIN_VISC==2)
895 write(10, fmt=trim(fmt3)) 'sigma_res =', sigma_res
896 write(10, fmt=trim(fmt1)) ' '
897 #endif
898 
899 write(10, fmt=trim(fmt2)) 'ENHMOD =', enhmod
900 #if (ENHMOD==1 || ENHMOD==2 || ENHMOD==3)
901 write(10, fmt=trim(fmt3)) 'enh_fact =', enh_fact
902 #endif
903 #if (ENHMOD==2 || ENHMOD==3)
904 write(10, fmt=trim(fmt3)) 'enh_intg =', enh_intg
905 #endif
906 #if (ENHMOD==2)
907 write(10, fmt=trim(fmt3)) 'age_trans =', age_trans_0
908 #endif
909 #if (ENHMOD==3)
910 write(10, fmt=trim(fmt3)) 'date_trans1 =', date_trans1_0
911 write(10, fmt=trim(fmt3)) 'date_trans2 =', date_trans2_0
912 write(10, fmt=trim(fmt3)) 'date_trans3 =', date_trans3_0
913 #endif
914 #if (ENHMOD==4 || ENHMOD==5)
915 write(10, fmt=trim(fmt3)) 'enh_compr =', enh_compr
916 write(10, fmt=trim(fmt3)) 'enh_shear =', enh_shear
917 #endif
918 #if ((ENHMOD==1 || ENHMOD==2 || ENHMOD==3 || ENHMOD==4) && MARGIN==3)
919 write(10, fmt=trim(fmt3)) 'enh_shelf =', enh_shelf
920 #endif
921 write(10, fmt=trim(fmt1)) ' '
922 
923 write(10, fmt=trim(fmt3)) 'numdiff_H_t =', numdiff_h_t
924 write(10, fmt=trim(fmt3)) 'tau_cts =', tau_cts
925 write(10, fmt=trim(fmt3)) 'vh_max =', vh_max
926 write(10, fmt=trim(fmt3)) 'hd_min =', hd_min
927 write(10, fmt=trim(fmt3)) 'hd_max =', hd_max
928 #if (defined(QBM_MIN))
929 write(10, fmt=trim(fmt3)) 'qbm_min =', qbm_min
930 #elif (defined(QB_MIN)) /* obsolete */
931 write(10, fmt=trim(fmt3)) 'qb_min =', qb_min
932 #endif
933 #if (defined(QBM_MAX))
934 write(10, fmt=trim(fmt3)) 'qbm_max =', qbm_max
935 #elif (defined(QB_MAX)) /* obsolete */
936 write(10, fmt=trim(fmt3)) 'qb_max =', qb_max
937 #endif
938 write(10, fmt=trim(fmt3)) 'age_min =', age_min
939 write(10, fmt=trim(fmt3)) 'age_max =', age_max
940 write(10, fmt=trim(fmt3)) 'mean_accum =', mean_accum
941 #if (ADV_VERT==1)
942 write(10, fmt=trim(fmt3)) 'age_diff =', agediff
943 #endif
944 write(10, fmt=trim(fmt1)) ' '
945 
946 write(10, fmt=trim(fmt2a)) 'DYNAMICS = ', dynamics
947 #if ((DYNAMICS==1 && MARGIN==3) || DYNAMICS==2)
948 #if (defined(LIS_OPTS))
949 write(10, fmt=trim(fmt1)) 'lis_opts = '//lis_opts
950 #endif
951 #if (defined(N_ITER_SSA))
952 write(10, fmt=trim(fmt2a)) 'n_iter_ssa = ', n_iter_ssa
953 #endif
954 #if (defined(RELAX_FACT_SSA))
955 write(10, fmt=trim(fmt3)) 'relax_fact_ssa =', relax_fact_ssa
956 #endif
957 #endif
958 #if (DYNAMICS==2 && defined(RATIO_SL_THRESH))
959 write(10, fmt=trim(fmt3)) 'ratio_sl_thresh =', ratio_sl_thresh
960 #endif
961 #if (DYNAMICS==2 && defined(SSTA_SIA_WEIGH_FCT))
962 write(10, fmt=trim(fmt2a)) 'SSTA_SIA_WEIGH_FCT = ', ssta_sia_weigh_fct
963 #endif
964 write(10, fmt=trim(fmt1)) ' '
965 
966 write(10, fmt=trim(fmt2)) 'CALCMOD =', calcmod
967 #if (CALCMOD==-1 && defined(TEMP_CONST))
968 write(10, fmt=trim(fmt3)) 'TEMP_CONST =', temp_const
969 #endif
970 #if (CALCMOD==-1 && defined(AGE_CONST))
971 write(10, fmt=trim(fmt3)) 'AGE_CONST =', age_const
972 #endif
973 #if (CALCMOD==1 && defined(CTS_MELTING_FREEZING))
974 write(10, fmt=trim(fmt2)) 'CTS_MELTING_FREEZING =', cts_melting_freezing
975 #endif
976 write(10, fmt=trim(fmt2)) 'FLOW_LAW =', flow_law
977 write(10, fmt=trim(fmt2)) 'FIN_VISC =', fin_visc
978 write(10, fmt=trim(fmt2)) 'MARGIN =', margin
979 #if (MARGIN==2)
980 write(10, fmt=trim(fmt2)) 'MARINE_ICE_FORMATION =', marine_ice_formation
981 write(10, fmt=trim(fmt2)) 'MARINE_ICE_CALVING =', marine_ice_calving
982 #elif (MARGIN==3)
983 write(10, fmt=trim(fmt2)) 'ICE_SHELF_CALVING =', ice_shelf_calving
984 #endif
985 #if (defined(THK_EVOL))
986 write(10, fmt=trim(fmt2)) 'THK_EVOL =', thk_evol
987 #else
988 stop ' >>> sico_init: Define THK_EVOL in header file!'
989 #endif
990 #if (defined(CALCTHK))
991 write(10, fmt=trim(fmt2)) 'CALCTHK =', calcthk
992 #else
993 stop ' >>> sico_init: Define CALCTHK in header file!'
994 #endif
995 write(10, fmt=trim(fmt2)) 'ADV_HOR =', adv_hor
996 write(10, fmt=trim(fmt2)) 'ADV_VERT =', adv_vert
997 write(10, fmt=trim(fmt2)) 'TOPOGRAD =', topograd
998 write(10, fmt=trim(fmt1)) ' '
999 write(10, fmt=trim(fmt2)) 'TEMP_PRESENT_PARA =', temp_present_para
1000 write(10, fmt=trim(fmt2)) 'TSURFACE =', tsurface
1001 write(10, fmt=trim(fmt2)) 'ACCSURFACE =', accsurface
1002 #if (ACCSURFACE==5)
1003 write(10, fmt=trim(fmt2)) 'PRECIP_ANOM_INTERPOL =', precip_anom_interpol
1004 #endif
1005 write(10, fmt=trim(fmt2)) 'SOLID_PRECIP =', solid_precip
1006 write(10, fmt=trim(fmt2)) 'ABLSURFACE =', ablsurface
1007 #if (defined(MB_ACCOUNT))
1008 write(10, fmt=trim(fmt2a)) 'MB_ACCOUNT = ', mb_account
1009 #endif
1010 write(10, fmt=trim(fmt1)) ' '
1011 
1012 #if (defined(OUT_TIMES))
1013 write(10, fmt=trim(fmt2a)) 'OUT_TIMES = ', out_times
1014 #endif
1015 #if (defined(OUTPUT_INIT))
1016 write(10, fmt=trim(fmt2a)) 'OUTPUT_INIT = ', output_init
1017 #endif
1018 write(10, fmt=trim(fmt2a)) 'OUTPUT = ', output
1019 #if (OUTPUT==1 || OUTPUT==3)
1020 write(10, fmt=trim(fmt3)) 'dtime_out =' , dtime_out0
1021 #endif
1022 write(10, fmt=trim(fmt3)) 'dtime_ser =' , dtime_ser0
1023 #if (OUTPUT==1 || OUTPUT==2)
1024 write(10, fmt=trim(fmt2a)) 'ERGDAT = ', ergdat
1025 #endif
1026 write(10, fmt=trim(fmt2a)) 'NETCDF = ', netcdf
1027 #if (OUTPUT==2 || OUTPUT==3)
1028 write(10, fmt=trim(fmt2a)) 'n_output = ', n_output
1029 do n=1, n_output
1030  write(10, fmt=trim(fmt3)) 'time_output =' , time_output0(n)
1031 end do
1032 #endif
1033 write(10, fmt=trim(fmt1)) ' '
1034 
1035 write(10, fmt=trim(fmt1)) 'Program version and date: '//version//' / '//date
1036 
1037 close(10, status='keep')
1038 
1039 !-------- Conversion of time quantities --------
1040 
1041 year_zero = year_zero*year_sec ! a --> s
1042 time_init = time_init0*year_sec ! a --> s
1043 time_end = time_end0*year_sec ! a --> s
1044 dtime = dtime0*year_sec ! a --> s
1045 dtime_temp = dtime_temp0*year_sec ! a --> s
1046 #if (REBOUND==2)
1047 dtime_wss = dtime_wss0*year_sec ! a --> s
1048 #endif
1049 dtime_ser = dtime_ser0*year_sec ! a --> s
1050 #if (OUTPUT==1 || OUTPUT==3)
1051 dtime_out = dtime_out0*year_sec ! a --> s
1052 #endif
1053 #if (OUTPUT==2 || OUTPUT==3)
1054 do n=1, n_output
1055  time_output(n) = time_output0(n)*year_sec ! a --> s
1056 end do
1057 #endif
1058 
1059 if (.not.approx_integer_multiple(dtime_temp, dtime, eps_sp_dp)) &
1060  stop ' >>> sico_init: dtime_temp must be a multiple of dtime!'
1061 
1062 #if (REBOUND==2)
1063 if (.not.approx_integer_multiple(dtime_wss, dtime, eps_sp_dp)) &
1064  stop ' >>> sico_init: dtime_wss must be a multiple of dtime!'
1065 #endif
1066 
1067 if (.not.approx_integer_multiple(dtime_ser, dtime, eps_sp_dp)) &
1068  stop ' >>> sico_init: dtime_ser must be a multiple of dtime!'
1069 
1070 #if (OUTPUT==1 || OUTPUT==3)
1071 if (.not.approx_integer_multiple(dtime_out, dtime, eps_sp_dp)) &
1072  stop ' >>> sico_init: dtime_out must be a multiple of dtime!'
1073 #endif
1074 
1075 #if (THK_EVOL==2)
1076 time_target_topo_init = time_target_topo_init0 *year_sec ! a --> s
1077 time_target_topo_final = time_target_topo_final0*year_sec ! a --> s
1078 #endif
1079 
1080 #if (THK_EVOL==3)
1081 target_topo_tau = target_topo_tau0 *year_sec ! a --> s
1082 #endif
1083 
1084 time = time_init
1085 
1086 !-------- Reading of present-day
1087 ! monthly mean or mean annual precipitation rate --------
1088 
1089 #if (GRID==0 || GRID==1)
1090 
1091 #if (defined(PRECIP_PRESENT_FILE) && defined(PRECIP_MA_PRESENT_FILE))
1092 
1093 if (flag_precip_monthly_mean) then
1094  filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1095  trim(precip_present_file)
1096 else
1097  filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1098  trim(precip_ma_present_file)
1099 end if
1100 
1101 #endif
1102 
1103 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1104 
1105 #elif (GRID==2)
1106 
1107 stop ' >>> sico_init: GRID==2 not allowed for this application!'
1108 
1109 #endif
1110 
1111 if (ios /= 0) stop ' >>> sico_init: Error when opening the precip file!'
1112 
1113 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
1114 
1115 if (flag_precip_monthly_mean) then
1116 
1117  do n=1, 12 ! month counter
1118  do m=1, 3; read(21, fmt='(a)') ch_dummy; end do
1119  do j=jmax, 0, -1
1120  read(21, fmt=*) (precip_present(j,i,n), i=0,imax)
1121  end do
1122  end do
1123 
1124  precip_present = precip_present *(1.0e-03_dp/year_sec)*(rho_w/rho)
1125  ! mm/a water equiv. -> m/s ice equiv.
1126 else
1127 
1128  do j=jmax, 0, -1
1129  read(21, fmt=*) (precip_ma_present(j,i), i=0,imax)
1130  end do
1131 
1132  precip_ma_present = precip_ma_present *(1.0e-03_dp/year_sec)*(rho_w/rho)
1133  ! mm/a water equiv. -> m/s ice equiv.
1134 end if
1135 
1136 close(21, status='keep')
1137 
1138 ! ------ Spatially dependent modification factor
1139 
1140 precip_factor = 1.0_dp
1141 
1142 #if (defined(PRECIP_FACTOR_FILE))
1143 
1144 if ( trim(adjustl(precip_factor_file)) /= 'none' ) then
1145 
1146  filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1147  trim(precip_factor_file)
1148 
1149  open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1150 
1151  if (ios /= 0) stop ' >>> sico_init: Error when opening the precip factor file!'
1152 
1153  do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1154 
1155  do j=jmax, 0, -1
1156  read(21, fmt=*) (precip_factor(j,i), i=0,imax)
1157  end do
1158 
1159  close(21, status='keep')
1160 
1161 end if
1162 
1163 #endif
1164 
1165 if (flag_precip_monthly_mean) then
1166 
1167  do n=1, 12 ! month counter
1168  do i=0, imax
1169  do j=0, jmax
1170  precip_present(j,i,n) = precip_present(j,i,n) * precip_factor(j,i)
1171  end do
1172  end do
1173  end do
1174 
1175 else
1176 
1177  precip_ma_present = precip_ma_present * precip_factor
1178 
1179 end if
1180 
1181 ! ------ Computation of the still undefined present-day
1182 ! mean annual or monthly mean precipitation rate
1183 
1184 if (flag_precip_monthly_mean) then
1185 
1186  precip_ma_present = 0.0_dp ! initialisation
1187 
1188  do n=1, 12 ! month counter
1189  do i=0, imax
1190  do j=0, jmax
1191  precip_ma_present(j,i) = precip_ma_present(j,i) &
1192  + precip_present(j,i,n)*inv_twelve
1193  end do
1194  end do
1195  end do
1196 
1197 else
1198 
1199  do n=1, 12 ! month counter
1200  do i=0, imax
1201  do j=0, jmax
1202  precip_present(j,i,n) = precip_ma_present(j,i)
1203  ! monthly means assumed to be equal
1204  ! to the mean annual precipitation rate
1205  end do
1206  end do
1207  end do
1208 
1209 end if
1210 
1211 ! ------ Reference topography for present-day precipitation rate
1212 
1213 #if (GRID==0 || GRID==1)
1214 
1215 #if (defined(PRECIP_ZS_REF_FILE))
1216 
1217 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1218  trim(precip_zs_ref_file)
1219 
1220 #endif
1221 
1222 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1223 
1224 #elif (GRID==2)
1225 
1226 stop ' >>> sico_init: GRID==2 not allowed for this application!'
1227 
1228 #endif
1229 
1230 if (ios /= 0) stop ' >>> sico_init: Error when opening the zs_ref file!'
1231 
1232 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1233 
1234 do j=jmax, 0, -1
1235  read(21, fmt=*) (zs_ref(j,i), i=0,imax)
1236 end do
1237 
1238 close(21, status='keep')
1239 
1240 do i=0, imax
1241 do j=0, jmax
1242  zs_ref(j,i) = max(zs_ref(j,i), 0.0_dp)
1243  ! resetting negative elevations (bathymetry data)
1244  ! to the present-day sea surface
1245 end do
1246 end do
1247 
1248 !-------- Reading of LGM mean-annual precipitation-rate anomaly --------
1249 
1250 #if (ACCSURFACE==5)
1251 
1252 #if (GRID==0 || GRID==1)
1253 
1254 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1255  trim(precip_anom_file)
1256 
1257 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1258 
1259 #endif
1260 
1261 if (ios /= 0) stop ' >>> sico_init: Error when opening the precip anomaly file!'
1262 
1263 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1264 
1265 do j=jmax, 0, -1
1266  read(21, fmt=*) (precip_ma_lgm_anom(j,i), i=0,imax)
1267 end do
1268 
1269 close(21, status='keep')
1270 
1271 precip_ma_lgm_anom = precip_ma_lgm_anom * precip_anom_fact
1272 
1273 !-------- LGM monthly precipitation-rate anomalies (assumed to be
1274 ! equal to the mean annual precipitation-rate anomaly) --------
1275 
1276 do n=1, 12 ! month counter
1277 do i=0, imax
1278 do j=0, jmax
1279 
1280  precip_lgm_anom(j,i,n) = max(precip_ma_lgm_anom(j,i), eps)
1281  ! positive values ensured
1282 
1283 #if (PRECIP_ANOM_INTERPOL==1)
1284  gamma_precip_lgm_anom(j,i,n) = 0.0_dp ! dummy values
1285 #elif (PRECIP_ANOM_INTERPOL==2)
1286  gamma_precip_lgm_anom(j,i,n) = -log(precip_lgm_anom(j,i,n))
1287 #else
1288  stop ' >>> sico_init: Wrong value of switch PRECIP_ANOM_INTERPOL!'
1289 #endif
1290 
1291 end do
1292 end do
1293 end do
1294 
1295 #endif
1296 
1297 !-------- Mean accumulation --------
1298 
1299 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho)
1300  ! mm/a water equiv. -> m/s ice equiv.
1301 
1302 !-------- Read rock/sediment mask --------
1303 
1304 #if (ICE_STREAM==2)
1305 
1306 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1307  trim(mask_sedi_file)
1308 
1309 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
1310 
1311 if (ios /= 0) stop ' >>> sico_init: Error when opening the rock/sediment mask file!'
1312 
1313 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1314 
1315 do j=jmax, 0, -1
1316  read(24, fmt=trim(fmt4)) (maske_sedi(j,i), i=0,imax)
1317 end do
1318 
1319 close(24, status='keep')
1320 
1321 #endif
1322 
1323 !-------- Reading of the prescribed target topography --------
1324 
1325 #if ((THK_EVOL==2) || (THK_EVOL==3))
1326 
1327 target_topo_dat_name = trim(target_topo_dat_name)
1328 
1329 #if (NETCDF==1)
1330 write(6, fmt='(a)') ' >>> sico_init: Reading of target topography'
1331 write(6, fmt='(a)') ' only implemented for NetCDF files (NETCDF==2)!'
1332 stop
1333 #elif (NETCDF==2)
1334 call read_target_topo_nc(target_topo_dat_name)
1335 #else
1336 stop ' >>> sico_init: Parameter NETCDF must be either 1 or 2!'
1337 #endif
1338 
1339 #endif
1340 
1341 !-------- Reading of the maximum ice extent mask --------
1342 
1343 #if (THK_EVOL==4)
1344 
1345 #if (GRID==0 || GRID==1)
1346 
1347 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1348  trim(mask_maxextent_file)
1349 
1350 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
1351 
1352 #elif (GRID==2)
1353 
1354 stop ' >>> sico_init: GRID==2 not allowed for this application!'
1355 
1356 #endif
1357 
1358 if (ios /= 0) &
1359  stop ' >>> sico_init: Error when opening the maximum ice extent mask file!'
1360 
1361 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1362 
1363 do j=jmax, 0, -1
1364  read(24, fmt=trim(fmt4)) (maske_maxextent(j,i), i=0,imax)
1365 end do
1366 
1367 close(24, status='keep')
1368 
1369 #endif
1370 
1371 !-------- Reading of LGM mean-annual and mean-July surface-temperature
1372 ! anomalies --------
1373 
1374 #if (TSURFACE==5)
1375 
1376 #if (GRID==0 || GRID==1)
1377 
1378 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1379  trim(temp_ma_anom_file)
1380 
1381 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1382 
1383 #endif
1384 
1385 if (ios /= 0) stop ' >>> sico_init: Error when opening the temp_ma anomaly file!'
1386 
1387 #if (GRID==0 || GRID==1)
1388 
1389 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1390  trim(temp_mj_anom_file)
1391 
1392 open(22, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1393 
1394 #endif
1395 
1396 if (ios /= 0) stop ' >>> sico_init: Error when opening the temp_mj anomaly file!'
1397 
1398 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1399 do n=1, 6; read(22, fmt='(a)') ch_dummy; end do
1400 
1401 do j=jmax, 0, -1
1402  read(21, fmt=*) (temp_ma_lgm_anom(j,i), i=0,imax)
1403  read(22, fmt=*) (temp_mj_lgm_anom(j,i), i=0,imax)
1404 end do
1405 
1406 close(21, status='keep')
1407 close(22, status='keep')
1408 
1409 temp_ma_lgm_anom = temp_ma_lgm_anom * temp_ma_anom_fact
1410 temp_mj_lgm_anom = temp_mj_lgm_anom * temp_mj_anom_fact
1411 
1412 #endif
1413 
1414 !-------- Read data for delta_ts --------
1415 
1416 #if (TSURFACE==4)
1417 
1418 filename_with_path = trim(inpath)//'/general/'//trim(grip_temp_file)
1419 
1420 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1421 
1422 if (ios /= 0) stop ' >>> sico_init: Error when opening the data file for delta_ts!'
1423 
1424 read(21, fmt=*) ch_dummy, grip_time_min, grip_time_stp, grip_time_max
1425 
1426 if (ch_dummy /= '#') then
1427  write(6, fmt=*) ' >>> sico_init: grip_time_min, grip_time_stp, grip_time_max'
1428  write(6, fmt=*) ' not defined in data file for delta_ts!'
1429  stop
1430 end if
1431 
1433 
1434 allocate(griptemp(0:ndata_grip))
1435 
1436 do n=0, ndata_grip
1437  read(21, fmt=*) d_dummy, griptemp(n)
1438 end do
1439 
1440 close(21, status='keep')
1441 
1442 #endif
1443 
1444 !-------- Read data for the glacial index --------
1445 
1446 #if (TSURFACE==5)
1447 
1448 filename_with_path = trim(inpath)//'/general/'//trim(glac_ind_file)
1449 
1450 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1451 
1452 if (ios /= 0) stop ' >>> sico_init: Error when opening the glacial-index file!'
1453 
1454 read(21, fmt=*) ch_dummy, gi_time_min, gi_time_stp, gi_time_max
1455 
1456 if (ch_dummy /= '#') then
1457  write(6, fmt=*) ' >>> sico_init: gi_time_min, gi_time_stp, gi_time_max'
1458  write(6, fmt=*) ' not defined in glacial-index file!'
1459  stop
1460 end if
1461 
1463 
1464 allocate(glacial_index(0:ndata_gi))
1465 
1466 do n=0, ndata_gi
1467  read(21, fmt=*) d_dummy, glacial_index(n)
1468 end do
1469 
1470 close(21, status='keep')
1471 
1472 #endif
1473 
1474 !-------- Reading of time-dependent surface-temperature
1475 ! and precipitation anomaly data --------
1476 
1477 #if (TSURFACE==6 && ACCSURFACE==6)
1478 
1479 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)// &
1480  '/'//trim(temp_precip_anom_file)
1481 
1482 call check( nf90_open(trim(filename_with_path), nf90_nowrite, &
1483  ncid_temp_precip), thisroutine )
1484 
1485 call check( nf90_inquire(ncid_temp_precip, unlimiteddimid = dimid), &
1486  thisroutine )
1487 
1488 call check( nf90_inquire_dimension(ncid_temp_precip, dimid, &
1489  len = ndata_temp_precip), thisroutine )
1490 
1492 
1494 
1495 call check( nf90_inq_varid(ncid_temp_precip, 'time', ncv), thisroutine )
1496 call check( nf90_get_var(ncid_temp_precip, ncv, temp_precip_time), thisroutine )
1497  ! in a
1498 
1501  ! in a; constant time step assumed
1503 
1504 #endif
1505 
1506 !-------- Reading of ISMIP6 InitMIP SMB anomaly data --------
1507 
1508 #if (defined(INITMIP_SMB_ANOM_FILE))
1509 
1510 if ( trim(adjustl(initmip_smb_anom_file)) /= 'none' ) then
1511 
1512  filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)// &
1513  '/'//trim(initmip_smb_anom_file)
1514 
1515  call check( nf90_open(trim(filename_with_path), nf90_nowrite, ncid), &
1516  thisroutine )
1517 
1518  call check( nf90_inq_varid(ncid, 'DSMB', ncv) )
1519  call check( nf90_get_var(ncid, ncv, as_anom_initmip_conv) )
1520 
1521  call check( nf90_close(ncid) )
1522 
1523  do i=0, imax
1524  do j=0, jmax
1525  as_anom_initmip(j,i) = as_anom_initmip_conv(i,j) /year_sec
1526  ! m/a ice equiv. -> m/s ice equiv.
1527  end do
1528  end do
1529 
1530 else
1531  as_anom_initmip = 0.0_dp
1532 end if
1533 
1534 #endif
1535 
1536 !-------- Reading of the global annual temperature anomaly
1537 ! (for parameterising the sub-ocean temperature anomaly
1538 ! for the ice discharge parameterisation)
1539 
1540 #if (DISC==2)
1541 
1542 filename_with_path = trim(inpath)//'/general/dTg_paleo.dat'
1543 
1544 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1545 
1546 if (ios /= 0) &
1547  stop ' >>> sico_init: Error when opening the data file for dT_glann_CLIMBER!'
1548 
1549 read(21, fmt=*) ch_dummy, glann_time_min, glann_time_stp, glann_time_max
1550 
1551 if (ch_dummy /= '#') then
1552  write(6, fmt=*) ' >>> sico_init: glann_time_min, glann_time_stp, glann_time_max'
1553  write(6, fmt=*) ' not defined in data file for dT_glann_CLIMBER!'
1554  stop
1555 end if
1556 
1557 ndata_glann = (glann_time_max-glann_time_min)/glann_time_stp
1558 
1559 allocate(dt_glann_climber(0:ndata_glann))
1560 
1561 do n=0, ndata_glann
1562  read(21, fmt=*) d_dummy, dt_glann_climber(n)
1563 end do
1564 
1565 close(21, status='keep')
1566 
1567 #endif
1568 
1569 !-------- Read data for z_sl --------
1570 
1571 #if (SEA_LEVEL==3)
1572 
1573 filename_with_path = trim(inpath)//'/general/'//trim(sea_level_file)
1574 
1575 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1576 
1577 if (ios /= 0) stop ' >>> sico_init: Error when opening the data file for z_sl!'
1578 
1579 read(21, fmt=*) ch_dummy, specmap_time_min, specmap_time_stp, specmap_time_max
1580 
1581 if (ch_dummy /= '#') then
1582  write(6, fmt=*) ' >>> sico_init: specmap_time_min, specmap_time_stp, specmap_time_max'
1583  write(6, fmt=*) ' not defined in data file for z_sl!'
1584  stop
1585 end if
1586 
1588 
1589 allocate(specmap_zsl(0:ndata_specmap))
1590 
1591 do n=0, ndata_specmap
1592  read(21, fmt=*) d_dummy, specmap_zsl(n)
1593 end do
1594 
1595 close(21, status='keep')
1596 
1597 #endif
1598 
1599 !-------- Determination of the geothermal heat flux --------
1600 
1601 #if (Q_GEO_MOD==1)
1602 
1603 ! ------ Constant value
1604 
1605 do i=0, imax
1606 do j=0, jmax
1607  q_geo(j,i) = q_geo *1.0e-03_dp ! mW/m2 -> W/m2
1608 end do
1609 end do
1610 
1611 #elif (Q_GEO_MOD==2)
1612 
1613 ! ------ Read data from file
1614 
1615 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1616  trim(q_geo_file)
1617 
1618 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1619 
1620 if (ios /= 0) stop ' >>> sico_init: Error when opening the qgeo file!'
1621 
1622 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1623 
1624 do j=jmax, 0, -1
1625  read(21, fmt=*) (q_geo(j,i), i=0,imax)
1626 end do
1627 
1628 close(21, status='keep')
1629 
1630 do i=0, imax
1631 do j=0, jmax
1632  q_geo(j,i) = q_geo(j,i) *1.0e-03_dp ! mW/m2 -> W/m2
1633 end do
1634 end do
1635 
1636 #endif
1637 
1638 !-------- Reading of tabulated kei function--------
1639 
1640 #if (REBOUND==0 || REBOUND==1)
1641 
1642 kei = 0.0_dp; n_data_kei = 0; kei_r_max = 0.0_dp; kei_r_incr = 0.0_dp
1643  ! dummy values
1644 #elif (REBOUND==2)
1645 
1646 call read_kei()
1647 
1648 #endif
1649 
1650 !-------- Determination of the time lag
1651 ! of the relaxing asthenosphere --------
1652 
1653 #if (REBOUND==1 || REBOUND==2)
1654 
1655 #if (TIME_LAG_MOD==1)
1656 
1657 time_lag_asth = time_lag*year_sec ! a -> s
1658 
1659 #elif (TIME_LAG_MOD==2)
1660 
1661 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1662  trim(time_lag_file)
1663 
1664 open(29, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1665 
1666 if (ios /= 0) stop ' >>> sico_init: Error when opening the time-lag file!'
1667 
1668 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1669 
1670 do j=jmax, 0, -1
1671  read(29, fmt=*) (time_lag_asth(j,i), i=0,imax)
1672 end do
1673 
1674 close(29, status='keep')
1675 
1676 time_lag_asth = time_lag_asth*year_sec ! a -> s
1677 
1678 #endif
1679 
1680 #elif (REBOUND==0)
1681 
1682 time_lag_asth = 0.0_dp ! dummy values
1683 
1684 #endif
1685 
1686 !-------- Determination of the flexural rigidity of the lithosphere --------
1687 
1688 #if (REBOUND==2)
1689 
1690 #if (FLEX_RIG_MOD==1)
1691 
1692 flex_rig_lith = flex_rig
1693 
1694 #elif (FLEX_RIG_MOD==2)
1695 
1696 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1697  trim(flex_rig_file)
1698 
1699 open(29, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1700 
1701 if (ios /= 0) stop ' >>> sico_init: Error when opening the flex-rig file!'
1702 
1703 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1704 
1705 do j=jmax, 0, -1
1706  read(29, fmt=*) (flex_rig_lith(j,i), i=0,imax)
1707 end do
1708 
1709 close(29, status='keep')
1710 
1711 #endif
1712 
1713 #elif (REBOUND==0 || REBOUND==1)
1714 
1715 flex_rig_lith = 0.0_dp ! dummy values
1716 
1717 #endif
1718 
1719 !-------- Definition of initial values --------
1720 
1721 ! ------ Present topography
1722 
1723 #if (ANF_DAT==1)
1724 
1725 call topography1(dxi, deta)
1726 
1727 #if (DISC>0) /* Ice discharge parameterization */
1728 call disc_param(dtime)
1729 call disc_fields()
1730 #endif
1731 
1732 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1733 
1734 call boundary(time_init, dtime, dxi, deta, &
1735  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1736 
1737 where ((maske==0_i2b).or.(maske==3_i2b))
1738  ! grounded or floating ice
1740 elsewhere ! maske==1_i2b or 2_i2b, ice-free land or sea
1741  as_perp_apl = 0.0_dp
1742 end where
1743 
1744 q_bm = 0.0_dp
1745 q_tld = 0.0_dp
1746 q_b_tot = 0.0_dp
1747 
1748 p_b_w = 0.0_dp
1749 h_w = 0.0_dp
1750 
1751 #if (!defined(TEMP_INIT) || TEMP_INIT==1)
1753 #elif (TEMP_INIT==2)
1755 #elif (TEMP_INIT==3)
1757 #elif (TEMP_INIT==4)
1759 #else
1760  write(6, fmt='(a)') ' >>> sico_init:'
1761  write(6, fmt='(a)') ' TEMP_INIT must be set to either 1, 2, 3 or 4!'
1762  stop
1763 #endif
1764 
1765 #if (ENHMOD==1)
1766  call calc_enhance_1()
1767 #elif (ENHMOD==2)
1768  call calc_enhance_2()
1769 #elif (ENHMOD==3)
1770  call calc_enhance_3(time_init)
1771 #elif (ENHMOD==4)
1772  call calc_enhance_4()
1773 #elif (ENHMOD==5)
1774  call calc_enhance_5()
1775 #else
1776  stop ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
1777 #endif
1778 
1779 ! ------ Ice-free, relaxed bedrock
1780 
1781 #elif (ANF_DAT==2)
1782 
1783 call topography2(dxi, deta)
1784 
1785 #if (DISC>0) /* Ice discharge parameterization */
1786 call disc_param(dtime)
1787 call disc_fields()
1788 #endif
1789 
1790 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1791 
1792 call boundary(time_init, dtime, dxi, deta, &
1793  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1794 
1795 as_perp_apl = 0.0_dp
1796 
1797 q_bm = 0.0_dp
1798 q_tld = 0.0_dp
1799 q_b_tot = 0.0_dp
1800 
1801 p_b_w = 0.0_dp
1802 h_w = 0.0_dp
1803 
1804 call init_temp_water_age_2()
1805 
1806 #if (ENHMOD==1)
1807  call calc_enhance_1()
1808 #elif (ENHMOD==2)
1809  call calc_enhance_2()
1810 #elif (ENHMOD==3)
1811  call calc_enhance_3(time_init)
1812 #elif (ENHMOD==4)
1813  call calc_enhance_4()
1814 #elif (ENHMOD==5)
1815  call calc_enhance_5()
1816 #else
1817  stop ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
1818 #endif
1819 
1820 ! ------ Read initial state from output of previous simulation
1821 
1822 #elif (ANF_DAT==3)
1823 
1824 call topography3(dxi, deta, z_sl, anfdatname)
1825 
1826 #if (DISC>0) /* Ice discharge parameterization */
1827 call disc_param(dtime)
1828 call disc_fields()
1829 #endif
1830 
1831 call boundary(time_init, dtime, dxi, deta, &
1832  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1833 
1834 where ((maske==0_i2b).or.(maske==3_i2b))
1835  ! grounded or floating ice
1837 elsewhere ! maske==1_i2b or 2_i2b, ice-free land or sea
1838  as_perp_apl = 0.0_dp
1839 end where
1840 
1841 q_b_tot = q_bm + q_tld
1842 
1843 #endif
1844 
1845 !-------- Inner-point flag --------
1846 
1847 flag_inner_point = .true.
1848 
1849 flag_inner_point(0,:) = .false.
1850 flag_inner_point(jmax,:) = .false.
1851 
1852 flag_inner_point(:,0) = .false.
1853 flag_inner_point(:,imax) = .false.
1854 
1855 !-------- Grounding line and calving front flags --------
1856 
1857 flag_grounding_line_1 = .false.
1858 flag_grounding_line_2 = .false.
1859 
1860 flag_calving_front_1 = .false.
1861 flag_calving_front_2 = .false.
1862 
1863 #if (MARGIN==1 || MARGIN==2)
1864 
1865 continue
1866 
1867 #elif (MARGIN==3)
1868 
1869 do i=1, imax-1
1870 do j=1, jmax-1
1871 
1872  if ( (maske(j,i)==0_i2b) & ! grounded ice
1873  .and. &
1874  ( (maske(j,i+1)==3_i2b) & ! with
1875  .or.(maske(j,i-1)==3_i2b) & ! one
1876  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1877  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1878  ) &
1879  flag_grounding_line_1(j,i) = .true.
1880 
1881  if ( (maske(j,i)==3_i2b) & ! floating ice
1882  .and. &
1883  ( (maske(j,i+1)==0_i2b) & ! with
1884  .or.(maske(j,i-1)==0_i2b) & ! one
1885  .or.(maske(j+1,i)==0_i2b) & ! neighbouring
1886  .or.(maske(j-1,i)==0_i2b) ) & ! grounded ice point
1887  ) &
1888  flag_grounding_line_2(j,i) = .true.
1889 
1890  if ( (maske(j,i)==3_i2b) & ! floating ice
1891  .and. &
1892  ( (maske(j,i+1)==2_i2b) & ! with
1893  .or.(maske(j,i-1)==2_i2b) & ! one
1894  .or.(maske(j+1,i)==2_i2b) & ! neighbouring
1895  .or.(maske(j-1,i)==2_i2b) ) & ! ocean point
1896  ) &
1897  flag_calving_front_1(j,i) = .true.
1898 
1899  if ( (maske(j,i)==2_i2b) & ! ocean
1900  .and. &
1901  ( (maske(j,i+1)==3_i2b) & ! with
1902  .or.(maske(j,i-1)==3_i2b) & ! one
1903  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1904  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1905  ) &
1906  flag_calving_front_2(j,i) = .true.
1907 
1908 end do
1909 end do
1910 
1911 do i=1, imax-1
1912 
1913  j=0
1914  if ( (maske(j,i)==2_i2b) & ! ocean
1915  .and. (maske(j+1,i)==3_i2b) & ! with one neighbouring
1916  ) & ! floating ice point
1917  flag_calving_front_2(j,i) = .true.
1918 
1919  j=jmax
1920  if ( (maske(j,i)==2_i2b) & ! ocean
1921  .and. (maske(j-1,i)==3_i2b) & ! with one neighbouring
1922  ) & ! floating ice point
1923  flag_calving_front_2(j,i) = .true.
1924 
1925 end do
1926 
1927 do j=1, jmax-1
1928 
1929  i=0
1930  if ( (maske(j,i)==2_i2b) & ! ocean
1931  .and. (maske(j,i+1)==3_i2b) & ! with one neighbouring
1932  ) & ! floating ice point
1933  flag_calving_front_2(j,i) = .true.
1934 
1935  i=imax
1936  if ( (maske(j,i)==2_i2b) & ! ocean
1937  .and. (maske(j,i-1)==3_i2b) & ! with one neighbouring
1938  ) & ! floating ice point
1939  flag_calving_front_2(j,i) = .true.
1940 
1941 end do
1942 
1943 #else
1944 stop ' >>> sico_init: MARGIN must be either 1, 2 or 3!'
1945 #endif
1946 
1947 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
1948 
1949 #if (GRID==0 || GRID==1)
1950 
1951 do ir=-imax, imax
1952 do jr=-jmax, jmax
1953  dist_dxdy(jr,ir) = sqrt( (real(ir,dp)*dxi)**2 + (real(jr,dp)*deta)**2 )
1954  ! distortion due to stereographic projection not accounted for
1955 end do
1956 end do
1957 
1958 #endif
1959 
1960 !-------- Initial velocities --------
1961 
1962 call calc_temp_melt()
1963 
1964 #if (DYNAMICS==1 || DYNAMICS==2)
1965 
1966 call calc_vxy_b_sia(time, z_sl)
1967 call calc_vxy_sia(dzeta_c, dzeta_t)
1968 
1969 #if (MARGIN==3 || DYNAMICS==2)
1970 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
1971 #endif
1972 
1973 call calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
1974 
1975 #if (MARGIN==3)
1976 call calc_vz_floating(z_sl, dxi, deta, dzeta_c)
1977 #endif
1978 
1979 #elif (DYNAMICS==0)
1980 
1981 call calc_vxy_static()
1982 call calc_vz_static()
1983 
1984 #else
1985  stop ' >>> sico_init: DYNAMICS must be either 0, 1 or 2!'
1986 #endif
1987 
1988 call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
1989 
1990 !-------- Initial enthalpies --------
1991 
1992 #if (CALCMOD==0 || CALCMOD==-1)
1993 
1994 do i=0, imax
1995 do j=0, jmax
1996 
1997  do kc=0, kcmax
1998  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1999  end do
2000 
2001  do kt=0, ktmax
2002  enth_t(kt,j,i) = enth_c(0,j,i)
2003  end do
2004 
2005 end do
2006 end do
2007 
2008 #elif (CALCMOD==1)
2009 
2010 do i=0, imax
2011 do j=0, jmax
2012 
2013  do kc=0, kcmax
2014  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
2015  end do
2016 
2017  if ( (maske(j,i) == 0_i2b).and.(n_cts(j,i) == 1_i2b) ) then
2018  do kt=0, ktmax
2019  enth_t(kt,j,i) = enth_fct_temp_omega(temp_t_m(kt,j,i), omega_t(kt,j,i))
2020  end do
2021  else
2022  do kt=0, ktmax
2023  enth_t(kt,j,i) = enth_c(0,j,i)
2024  end do
2025  end if
2026 
2027 end do
2028 end do
2029 
2030 #elif (CALCMOD==2 || CALCMOD==3)
2031 
2032 do i=0, imax
2033 do j=0, jmax
2034 
2035  do kc=0, kcmax
2036  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), omega_c(kc,j,i))
2037  end do
2038 
2039  do kt=0, ktmax
2040  enth_t(kt,j,i) = enth_c(0,j,i)
2041  end do
2042 
2043 end do
2044 end do
2045 
2046 #else
2047 
2048 stop ' >>> sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!'
2049 
2050 #endif
2051 
2052 !-------- Initialize time-series files --------
2053 
2054 ! ------ Time-series file for the ice sheet on the whole
2055 
2056 filename_with_path = trim(outpath)//'/'//trim(runname)//'.ser'
2057 
2058 open(12, iostat=ios, file=trim(filename_with_path), status='new')
2059 
2060 if (ios /= 0) stop ' >>> sico_init: Error when opening the ser file!'
2061 
2062 if (forcing_flag == 1) then
2063 
2064  write(12,1102)
2065  write(12,1103)
2066 
2067  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
2068  ' V(m^3) V_g(m^3) V_f(m^3)', &
2069  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
2070  ' V_sle(m) V_t(m^3)', &
2071  ' A_t(m^2)',/, &
2072  ' H_max(m) H_t_max(m)', &
2073  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
2074  1103 format('----------------------------------------------------', &
2075  '---------------------------------------')
2076 
2077 else if (forcing_flag == 2) then
2078 
2079  write(12,1112)
2080  write(12,1113)
2081 
2082  1112 format(' t(a) glac_ind(1) z_sl(m)',/, &
2083  ' V(m^3) V_g(m^3) V_f(m^3)', &
2084  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
2085  ' V_sle(m) V_t(m^3)', &
2086  ' A_t(m^2)',/, &
2087  ' H_max(m) H_t_max(m)', &
2088  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
2089  1113 format('----------------------------------------------------', &
2090  '---------------------------------------')
2091 
2092 else if (forcing_flag == 3) then
2093 
2094  write(12,1122)
2095  write(12,1123)
2096 
2097  1122 format(' t(a) (Dummy)(1) z_sl(m)',/, &
2098  ' V(m^3) V_g(m^3) V_f(m^3)', &
2099  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
2100  ' V_sle(m) V_t(m^3)', &
2101  ' A_t(m^2)',/, &
2102  ' H_max(m) H_t_max(m)', &
2103  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
2104  1123 format('----------------------------------------------------', &
2105  '---------------------------------------')
2106 
2107 end if
2108 
2109 ! ------ Time-series file for deep boreholes
2110 
2111 n_core = 7 ! GRIP, GISP2, Dye3, Camp Century (CC),
2112  ! NorthGRIP (NGRIP), NEEM, EastGRIP (EGRIP)
2113 
2114 allocate(lambda_core(n_core), phi_core(n_core), &
2116 
2117 ch_core(1) = 'GRIP'
2118 phi_core(1) = 72.58722_dp *pi_180 ! Geographical position of GRIP,
2119 lambda_core(1) = -37.64222_dp *pi_180 ! conversion deg -> rad
2120 
2121 ch_core(2) = 'GISP2'
2122 phi_core(2) = 72.58833_dp *pi_180 ! Geographical position of GISP2
2123 lambda_core(2) = -38.45750_dp *pi_180 ! conversion deg -> rad
2124 
2125 ch_core(3) = 'Dye3'
2126 phi_core(3) = 65.15139_dp *pi_180 ! Geographical position of Dye3,
2127 lambda_core(3) = -43.81722_dp *pi_180 ! conversion deg -> rad
2128 
2129 ch_core(4) = 'Camp Century'
2130 phi_core(4) = 77.17970_dp *pi_180 ! Geographical position of CC,
2131 lambda_core(4) = -61.10975_dp *pi_180 ! conversion deg -> rad
2132 
2133 ch_core(5) = 'NGRIP'
2134 phi_core(5) = 75.09694_dp *pi_180 ! Geographical position of NGRIP,
2135 lambda_core(5) = -42.31956_dp *pi_180 ! conversion deg -> rad
2136 
2137 ch_core(6) = 'NEEM'
2138 phi_core(6) = 77.5_dp *pi_180 ! Geographical position of NEEM,
2139 lambda_core(6) = -50.9_dp *pi_180 ! conversion deg -> rad
2140 
2141 ch_core(7) = 'EGRIP'
2142 phi_core(7) = 75.6299_dp *pi_180 ! Geographical position of EGRIP,
2143 lambda_core(7) = -35.9867_dp *pi_180 ! conversion deg -> rad
2144 
2145 #if (GRID==0 || GRID==1) /* Stereographic projection */
2146 
2147 do n=1, n_core
2149  a, b, lambda0, phi0, x_core(n), y_core(n))
2150 end do
2151 
2152 #elif (GRID==2) /* Geographical coordinates */
2153 
2155 y_core = phi_core
2156 
2157 #endif
2158 
2159 filename_with_path = trim(outpath)//'/'//trim(runname)//'.core'
2160 
2161 open(14, iostat=ios, file=trim(filename_with_path), status='new')
2162 
2163 if (ios /= 0) stop ' >>> sico_init: Error when opening the core file!'
2164 
2165 if (forcing_flag == 1) then
2166 
2167  write(14,1106)
2168  write(14,1107)
2169 
2170  1106 format(' t(a) D_Ts(C) z_sl(m)',/, &
2171  ' H_GR(m) H_G2(m) H_D3(m)', &
2172  ' H_CC(m) H_NG(m) H_NE(m) H_EG(m)',/, &
2173  ' v_GR(m/a) v_G2(m/a) v_D3(m/a)', &
2174  ' v_CC(m/a) v_NG(m/a) v_NE(m/a) v_EG(m/a)',/, &
2175  ' T_GR(C) T_G2(C) T_D3(C)', &
2176  ' T_CC(C) T_NG(C) T_NE(C) T_EG(C)')
2177  1107 format('----------------------------------------------------', &
2178  '----------------------------------------------------')
2179 
2180 else if (forcing_flag == 2) then
2181 
2182  write(14,1116)
2183  write(14,1117)
2184 
2185  1116 format(' t(a) glac_ind(1) z_sl(m)',/, &
2186  ' H_GR(m) H_G2(m) H_D3(m)', &
2187  ' H_CC(m) H_NG(m) H_NE(m) H_EG(m)',/, &
2188  ' v_GR(m/a) v_G2(m/a) v_D3(m/a)', &
2189  ' v_CC(m/a) v_NG(m/a) v_NE(m/a) v_EG(m/a)',/, &
2190  ' T_GR(C) T_G2(C) T_D3(C)', &
2191  ' T_CC(C) T_NG(C) T_NE(C) T_EG(C)')
2192  1117 format('----------------------------------------------------', &
2193  '----------------------------------------------------')
2194 
2195 else if (forcing_flag == 3) then
2196 
2197  write(14,1126)
2198  write(14,1127)
2199 
2200  1126 format(' t(a) (Dummy)(1) z_sl(m)',/, &
2201  ' H_GR(m) H_G2(m) H_D3(m)', &
2202  ' H_CC(m) H_NG(m) H_NE(m) H_EG(m)',/, &
2203  ' v_GR(m/a) v_G2(m/a) v_D3(m/a)', &
2204  ' v_CC(m/a) v_NG(m/a) v_NE(m/a) v_EG(m/a)',/, &
2205  ' T_GR(C) T_G2(C) T_D3(C)', &
2206  ' T_CC(C) T_NG(C) T_NE(C) T_EG(C)')
2207  1127 format('----------------------------------------------------', &
2208  '----------------------------------------------------')
2209 
2210 end if
2211 
2212 !-------- Output of the initial state --------
2213 
2214 #if (defined(OUTPUT_INIT))
2215 
2216 #if (OUTPUT_INIT==0)
2217  flag_init_output = .false.
2218 #elif (OUTPUT_INIT==1)
2219  flag_init_output = .true.
2220 #else
2221  stop ' >>> sico_init: OUTPUT_INIT must be either 0 or 1!'
2222 #endif
2223 
2224 #else
2225  flag_init_output = .true. ! default for undefined parameter OUTPUT_INIT
2226 #endif
2227 
2228 #if (OUTPUT==1)
2229 
2230 #if (ERGDAT==0)
2231  flag_3d_output = .false.
2232 #elif (ERGDAT==1)
2233  flag_3d_output = .true.
2234 #else
2235  stop ' >>> sico_init: ERGDAT must be either 0 or 1!'
2236 #endif
2237 
2238  if (flag_init_output) &
2239  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
2240  flag_3d_output, ndat2d, ndat3d)
2241 
2242 #elif (OUTPUT==2)
2243 
2244 if (time_output(1) <= time_init+eps) then
2245 
2246 #if (ERGDAT==0)
2247  flag_3d_output = .false.
2248 #elif (ERGDAT==1)
2249  flag_3d_output = .true.
2250 #else
2251  stop ' >>> sico_init: ERGDAT must be either 0 or 1!'
2252 #endif
2253 
2254  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
2255  flag_3d_output, ndat2d, ndat3d)
2256 
2257 end if
2258 
2259 #elif (OUTPUT==3)
2260 
2261  flag_3d_output = .false.
2262 
2263  if (flag_init_output) &
2264  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
2265  flag_3d_output, ndat2d, ndat3d)
2266 
2267 if (time_output(1) <= time_init+eps) then
2268 
2269  flag_3d_output = .true.
2270 
2271  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
2272  flag_3d_output, ndat2d, ndat3d)
2273 
2274 end if
2275 
2276 #else
2277  stop ' >>> sico_init: OUTPUT must be either 1, 2 or 3!'
2278 #endif
2279 
2280 if (flag_init_output) then
2281  call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
2282  call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
2283 end if
2284 
2285 #if (defined(EXEC_MAKE_C_DIS_0))
2286 
2287 #if (DISC>0)
2288 
2289 call calc_c_dis_0(dxi, deta)
2290 
2291 write(6, fmt='(a)') ' '
2292 write(6, fmt='(a)') ' * * * * * * * * * * * * * * * * * * * * * * * * * * * * '
2293 write(6, fmt='(a)') ' '
2294 write(6, fmt='(a)') ' sico_init: Routine calc_c_dis_0 successfully completed, '
2295 write(6, fmt='(a)') ' c_dis_0 written on file out_runname.dat '
2296 write(6, fmt='(a)') ' (in directory specified by OUTPATH). '
2297 write(6, fmt='(a)') ' Execution of SICOPOLIS stopped. '
2298 write(6, fmt='(a)') ' '
2299 write(6, fmt='(a)') ' * * * * * * * * * * * * * * * * * * * * * * * * * * * * '
2300 write(6, fmt='(a)') ' '
2301 
2302 stop
2303 
2304 #else
2305  stop ' >>> sico_init: EXEC_MAKE_C_DIS_0 requires DISC>0!'
2306 #endif
2307 
2308 #endif
2309 
2310 end subroutine sico_init
2311 
2312 !-------------------------------------------------------------------------------
2313 !> Definition of the initial surface and bedrock topography
2314 !! (including gradients) and of the horizontal grid spacings dxi, deta.
2315 !! For present-day initial topography.
2316 !<------------------------------------------------------------------------------
2317 subroutine topography1(dxi, deta)
2318 
2319 #if (GRID==0 || GRID==1)
2321 #endif
2322 
2323  use metric_m
2324  use topograd_m
2325 
2326 implicit none
2327 
2328 real(dp), intent(out) :: dxi, deta
2329 
2330 integer(i4b) :: i, j, n
2331 integer(i4b) :: ios, n_dummy
2332 real(dp) :: d_dummy
2333 real(dp) :: xi0, eta0
2334 real(dp) :: h_ice, freeboard_ratio
2335 character :: ch_dummy
2336 
2337 character(len= 8) :: ch_imax
2338 character(len=128) :: fmt4
2339 character(len=256) :: filename_with_path
2340 
2341 write(ch_imax, fmt='(i8)') imax
2342 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
2343 
2344 !-------- Read topography --------
2345 
2346 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2347  trim(zs_present_file)
2348 
2349 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
2350 
2351 if (ios /= 0) stop ' >>> topography1: Error when opening the zs file!'
2352 
2353 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2354  trim(zl_present_file)
2355 
2356 open(22, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
2357 
2358 if (ios /= 0) stop ' >>> topography1: Error when opening the zl file!'
2359 
2360 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2361  trim(zl0_file)
2362 
2363 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
2364 
2365 if (ios /= 0) stop ' >>> topography1: Error when opening the zl0 file!'
2366 
2367 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2368  trim(mask_present_file)
2369 
2370 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
2371 
2372 if (ios /= 0) stop ' >>> topography1: Error when opening the mask file!'
2373 
2374 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
2375 do n=1, 6; read(22, fmt='(a)') ch_dummy; end do
2376 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
2377 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
2378 
2379 do j=jmax, 0, -1
2380  read(21, fmt=*) (zs(j,i), i=0,imax)
2381  read(22, fmt=*) (zl(j,i), i=0,imax)
2382  read(23, fmt=*) (zl0(j,i), i=0,imax)
2383  read(24, fmt=trim(fmt4)) (maske(j,i), i=0,imax)
2384 end do
2385 
2386 close(21, status='keep')
2387 close(22, status='keep')
2388 close(23, status='keep')
2389 close(24, status='keep')
2390 
2391 #if (defined(ZB_PRESENT_FILE))
2392 
2393 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2394  trim(zb_present_file)
2395 
2396 open(25, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
2397 
2398 if (ios /= 0) stop ' >>> topography1: Error when opening the zb file!'
2399 
2400 do n=1, 6; read(25, fmt='(a)') ch_dummy; end do
2401 
2402 do j=jmax, 0, -1
2403  read(25, fmt=*) (zb(j,i), i=0,imax)
2404 end do
2405 
2406 close(25, status='keep')
2407 
2408 #else
2409 
2410 write(6, fmt='(a)') ' >>> topography1: ZB_PRESENT_FILE not defined,'
2411 write(6, fmt='(a)') ' thus zb = zl assumed.'
2412 
2413 zb = zl
2414 
2415 #endif
2416 
2417 !-------- Further stuff --------
2418 
2419 dxi = dx *1000.0_dp ! km -> m
2420 deta = dx *1000.0_dp ! km -> m
2421 
2422 xi0 = x0 *1000.0_dp ! km -> m
2423 eta0 = y0 *1000.0_dp ! km -> m
2424 
2425 freeboard_ratio = (rho_sw-rho)/rho_sw
2426 
2427 do i=0, imax
2428 do j=0, jmax
2429 
2430  if (maske(j,i) <= 1) then
2431 
2432  zb(j,i) = zl(j,i) ! ensure consistency
2433 
2434  else if (maske(j,i) == 2) then
2435 
2436 #if (MARGIN==1 || MARGIN==2)
2437  zs(j,i) = zl(j,i) ! ensure
2438  zb(j,i) = zl(j,i) ! consistency
2439 #elif (MARGIN==3)
2440  zs(j,i) = 0.0_dp ! present-day
2441  zb(j,i) = 0.0_dp ! sea level
2442 #endif
2443 
2444  else if (maske(j,i) == 3) then
2445 
2446 #if (MARGIN==1 || (MARGIN==2 && MARINE_ICE_FORMATION==1))
2447  maske(j,i) = 2 ! floating ice cut off
2448  zs(j,i) = zl(j,i)
2449  zb(j,i) = zl(j,i)
2450 #elif (MARGIN==2 && MARINE_ICE_FORMATION==2)
2451  maske(j,i) = 0 ! floating ice becomes "underwater ice"
2452  h_ice = zs(j,i)-zb(j,i) ! ice thickness
2453  zs(j,i) = zl(j,i)+h_ice
2454  zb(j,i) = zl(j,i)
2455 #elif (MARGIN==3)
2456  h_ice = zs(j,i)-zb(j,i) ! ice thickness
2457  zs(j,i) = freeboard_ratio*h_ice ! ensure properly
2458  zb(j,i) = zs(j,i)-h_ice ! floating ice
2459 #endif
2460 
2461  end if
2462 
2463  xi(i) = xi0 + real(i,dp)*dxi
2464  eta(j) = eta0 + real(j,dp)*deta
2465 
2466  zm(j,i) = zb(j,i)
2467  n_cts(j,i) = -1
2468  kc_cts(j,i) = 0
2469 
2470  h_c(j,i) = zs(j,i)-zm(j,i)
2471  h_t(j,i) = 0.0_dp
2472 
2473  dzs_dtau(j,i) = 0.0_dp
2474  dzm_dtau(j,i) = 0.0_dp
2475  dzb_dtau(j,i) = 0.0_dp
2476  dzl_dtau(j,i) = 0.0_dp
2477  dh_c_dtau(j,i) = 0.0_dp
2478  dh_t_dtau(j,i) = 0.0_dp
2479 
2480 end do
2481 end do
2482 
2483 maske_old = maske
2484 
2485 !-------- Geographic coordinates, metric tensor,
2486 ! gradients of the topography --------
2487 
2488 do i=0, imax
2489 do j=0, jmax
2490 
2491 #if (GRID==0 || GRID==1) /* Stereographic projection */
2492 
2493  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
2494  lambda0, phi0, lambda(j,i), phi(j,i))
2495 
2496 #elif (GRID==2) /* Geographic coordinates */
2497 
2498  lambda(j,i) = xi(i)
2499  phi(j,i) = eta(j)
2500 
2501 #endif
2502 
2503 end do
2504 end do
2505 
2506 call metric()
2507 
2508 #if (TOPOGRAD==0)
2509 call topograd_1(dxi, deta, 1)
2510 #elif (TOPOGRAD==1)
2511 call topograd_2(dxi, deta, 1)
2512 #endif
2513 
2514 !-------- Corresponding area of grid points --------
2515 
2516 do i=0, imax
2517 do j=0, jmax
2518  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
2519 end do
2520 end do
2521 
2522 end subroutine topography1
2523 
2524 !-------------------------------------------------------------------------------
2525 !> Definition of the initial surface and bedrock topography
2526 !! (including gradients) and of the horizontal grid spacings dxi, deta.
2527 !! For ice-free initial topography with relaxed lithosphere.
2528 !<------------------------------------------------------------------------------
2529 subroutine topography2(dxi, deta)
2530 
2531 #if (GRID==0 || GRID==1)
2533 #endif
2534 
2535  use metric_m
2536  use topograd_m
2537 
2538 implicit none
2539 
2540 real(dp), intent(out) :: dxi, deta
2541 
2542 integer(i4b) :: i, j, n
2543 integer(i4b) :: ios
2544 real(dp) :: xi0, eta0
2545 character :: ch_dummy
2546 
2547 character(len= 8) :: ch_imax
2548 character(len=128) :: fmt4
2549 character(len=256) :: filename_with_path
2550 
2551 write(ch_imax, fmt='(i8)') imax
2552 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
2553 
2554 !-------- Read topography --------
2555 
2556 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2557  trim(zl0_file)
2558 
2559 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
2560 
2561 if (ios /= 0) stop ' >>> topography2: Error when opening the zl0 file!'
2562 
2563 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2564  trim(mask_present_file)
2565 
2566 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
2567 
2568 if (ios /= 0) stop ' >>> topography2: Error when opening the mask file!'
2569 
2570 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
2571 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
2572 
2573 do j=jmax, 0, -1
2574  read(23, fmt=*) (zl0(j,i), i=0,imax)
2575  read(24, fmt=trim(fmt4)) (maske(j,i), i=0,imax)
2576 end do
2577 
2578 close(23, status='keep')
2579 close(24, status='keep')
2580 
2581 !-------- Further stuff --------
2582 
2583 dxi = dx *1000.0_dp ! km -> m
2584 deta = dx *1000.0_dp ! km -> m
2585 
2586 xi0 = x0 *1000.0_dp ! km -> m
2587 eta0 = y0 *1000.0_dp ! km -> m
2588 
2589 do i=0, imax
2590 do j=0, jmax
2591 
2592  if (maske(j,i) <= 1) then
2593  maske(j,i) = 1
2594  zs(j,i) = zl0(j,i)
2595  zb(j,i) = zl0(j,i)
2596  zl(j,i) = zl0(j,i)
2597  else ! (maske(j,i) >= 2)
2598  maske(j,i) = 2
2599 #if (MARGIN==1 || MARGIN==2)
2600  zs(j,i) = zl0(j,i)
2601  zb(j,i) = zl0(j,i)
2602 #elif (MARGIN==3)
2603  zs(j,i) = 0.0_dp ! present-day
2604  zb(j,i) = 0.0_dp ! sea level
2605 #endif
2606  zl(j,i) = zl0(j,i)
2607  end if
2608 
2609  xi(i) = xi0 + real(i,dp)*dxi
2610  eta(j) = eta0 + real(j,dp)*deta
2611 
2612  zm(j,i) = zb(j,i)
2613  n_cts(j,i) = -1
2614  kc_cts(j,i) = 0
2615 
2616  h_c(j,i) = 0.0_dp
2617  h_t(j,i) = 0.0_dp
2618 
2619  dzs_dtau(j,i) = 0.0_dp
2620  dzm_dtau(j,i) = 0.0_dp
2621  dzb_dtau(j,i) = 0.0_dp
2622  dzl_dtau(j,i) = 0.0_dp
2623  dh_c_dtau(j,i) = 0.0_dp
2624  dh_t_dtau(j,i) = 0.0_dp
2625 
2626 end do
2627 end do
2628 
2629 maske_old = maske
2630 
2631 !-------- Geographic coordinates, metric tensor,
2632 ! gradients of the topography --------
2633 
2634 do i=0, imax
2635 do j=0, jmax
2636 
2637 #if (GRID==0 || GRID==1) /* Stereographic projection */
2638 
2639  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
2640  lambda0, phi0, lambda(j,i), phi(j,i))
2641 
2642 #elif (GRID==2) /* Geographic coordinates */
2643 
2644  lambda(j,i) = xi(i)
2645  phi(j,i) = eta(j)
2646 
2647 #endif
2648 
2649 end do
2650 end do
2651 
2652 call metric()
2653 
2654 #if (TOPOGRAD==0)
2655 call topograd_1(dxi, deta, 1)
2656 #elif (TOPOGRAD==1)
2657 call topograd_2(dxi, deta, 1)
2658 #endif
2659 
2660 !-------- Corresponding area of grid points --------
2661 
2662 do i=0, imax
2663 do j=0, jmax
2664  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
2665 end do
2666 end do
2667 
2668 end subroutine topography2
2669 
2670 !-------------------------------------------------------------------------------
2671 !> Definition of the initial surface and bedrock topography
2672 !! (including gradients) and of the horizontal grid spacings dxi, deta.
2673 !! For initial topography from previous simulation.
2674 !<------------------------------------------------------------------------------
2675 subroutine topography3(dxi, deta, z_sl, anfdatname)
2676 
2677  use read_m, only : read_erg_nc
2678 
2679 #if (GRID==0 || GRID==1)
2681 #endif
2682 
2683  use metric_m
2684  use topograd_m
2685 
2686 implicit none
2687 
2688 character(len=100), intent(in) :: anfdatname
2689 
2690 real(dp), intent(out) :: dxi, deta, z_sl
2691 
2692 integer(i4b) :: i, j, n
2693 integer(i4b) :: ios
2694 character(len=256) :: filename_with_path
2695 character :: ch_dummy
2696 
2697 !-------- Read data from time-slice file of previous simulation --------
2698 
2699 call read_erg_nc(z_sl, anfdatname)
2700 
2701 !-------- Read topography of the relaxed bedrock --------
2702 
2703 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2704  trim(zl0_file)
2705 
2706 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
2707 
2708 if (ios /= 0) stop ' >>> topography3: Error when opening the zl0 file!'
2709 
2710 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
2711 
2712 do j=jmax, 0, -1
2713  read(23, fmt=*) (zl0(j,i), i=0,imax)
2714 end do
2715 
2716 close(23, status='keep')
2717 
2718 !-------- Further stuff --------
2719 
2720 dxi = dx *1000.0_dp ! km -> m
2721 deta = dx *1000.0_dp ! km -> m
2722 
2723 !-------- Geographic coordinates, metric tensor,
2724 ! gradients of the topography --------
2725 
2726 do i=0, imax
2727 do j=0, jmax
2728 
2729 #if (GRID==0 || GRID==1) /* Stereographic projection */
2730 
2731  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
2732  lambda0, phi0, lambda(j,i), phi(j,i))
2733 
2734 #elif (GRID==2) /* Geographic coordinates */
2735 
2736  lambda(j,i) = xi(i)
2737  phi(j,i) = eta(j)
2738 
2739 #endif
2740 
2741 end do
2742 end do
2743 
2744 call metric()
2745 
2746 #if (TOPOGRAD==0)
2747 call topograd_1(dxi, deta, 1)
2748 #elif (TOPOGRAD==1)
2749 call topograd_2(dxi, deta, 1)
2750 #endif
2751 
2752 !-------- Corresponding area of grid points --------
2753 
2754 do i=0, imax
2755 do j=0, jmax
2756  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
2757 end do
2758 end do
2759 
2760 end subroutine topography3
2761 
2762 !-------------------------------------------------------------------------------
2763 
2764 end module sico_init_m
2765 !
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c
temp_c(kc,j,i): Temperature in the upper (kc) ice domain
real(dp), dimension(0:jmax, 0:imax) c_slide
c_slide(j,i): Basal sliding coefficient
real(dp), dimension(0:jmax, 0:imax) phi
phi(j,i): Geographic latitude of grid point (i,j)
subroutine, public read_erg_nc(z_sl, filename)
Reading of time-slice files in native binary or NetCDF format.
Definition: read_m.F90:50
real(dp), dimension(:), allocatable phi_core
phi_core(n): Geographical latitude of the prescribed borehole positions
subroutine, public read_kei()
Reading of the tabulated kei function.
Definition: read_m.F90:765
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(i2b), dimension(0:jmax, 0:imax) n_cts
n_cts(j,i): Mask for thermal conditions. -1: cold ice base, 0: temperate ice base with cold ice above...
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 sico_init(delta_ts, glac_index, mean_accum, dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, time, time_init, time_end, time_output, dxi, deta, dzeta_c, dzeta_t, dzeta_r, z_sl, dzsl_dtau, z_mar, ii, jj, nn, ndat2d, ndat3d, n_output, runname)
Main routine of sico_init_m: Initialisations for SICOPOLIS.
Definition: sico_init_m.F90:59
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, 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:59
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) temp_precip_time_stp
temp_precip_time_stp: Time step of the surface-temperature and precipitation data ...
real(dp), parameter eps
eps: Small number
subroutine, public read_phys_para()
Reading of physical parameters.
Definition: read_m.F90:813
integer(i2b), dimension(0:jmax, 0:imax) maske_maxextent
maske_maxextent(j,i): Maximum ice extent mask. 0: not allowed to glaciate, 1: allowed to glaciate...
integer(i2b) forcing_flag
forcing_flag: Flag for the forcing type. 1: forcing by a spatially constant surface temperature anoma...
real(dp) temp_precip_time_max
temp_precip_time_max: Maximum time of the surface-temperature and precipitation data ...
logical, dimension(0:jmax, 0:imax) flag_grounding_line_1
flag_grounding_line_1(j,i): Grounding line flag. .true.: grounding line point (grounded ice point wit...
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) enth_t
enth_t(kt,j,i): Enthalpy in the lower (kt) ice domain
Conversion from temperature (temp) and water content (omega) to enthalpy (enth) and vice versa...
subroutine, public init_temp_water_age_2()
Initial temperature, water content and age (case ANF_DAT==2: ice-free conditions with relaxed bedrock...
real(dp), dimension(0:jmax, 0:imax) dh_t_dtau
dH_t_dtau(j,i): Derivative of H_t by tau (time)
real(dp), dimension(0:jmax, 0:imax) lambda
lambda(j,i): Geographic longitude of grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) dh_c_dtau
dH_c_dtau(j,i): Derivative of H_c by tau (time)
integer(i2b), dimension(0:jmax, 0:imax) maske
maske(j,i): Ice-land-ocean mask. 0: grounded ice, 1: ice-free land, 2: ocean, 3: floating ice ...
real(dp), dimension(0:jmax, 0:imax) p_b_w
p_b_w(j,i): Basal water pressure
integer(i4b) specmap_time_min
specmap_time_min: Minimum time of the data values for the sea level
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)
real(dp), dimension(:), allocatable glacial_index
glacial_index(n): Data values for the glacial index
subroutine, public calc_c_dis_0(dxi, deta)
Constant in ice discharge parameterization (Greenland). [Determine (amount of magnitude of) constant ...
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
integer(i2b), 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
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:56
real(dp), dimension(0:jmax, 0:imax) zs_ref
zs_ref(j,i): Reference elevation for precip_present, temp_ma_present and temp_mj_present ...
real(dp) target_topo_tau
target_topo_tau: Relaxation time for target-topography adjustment
subroutine, public calc_enhance_5()
Computation of the flow enhancement factor. Case ENHMOD==5: minimal anisotropic enhancement factor fo...
real(dp) b
B: Semi-minor axis of the planet.
subroutine, public init_temp_water_age_1_2()
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==2: present-day initial topogr...
real(dp), dimension(0:jmax, 0:imax) zl
zl(j,i): Coordinate z of the lithosphere surface
subroutine, public calc_vxy_static()
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy etc. for static ice.
subroutine, public output4(time, dxi, deta, delta_ts, glac_index, z_sl)
Writing of time-series data of the deep ice cores on file in ASCII format (and optionally in NetCDF f...
Definition: output_m.F90:4721
integer(i2b), dimension(0:jmax, 0:imax) kc_cts
kc_cts(j,i): Position kc of the CTS (for COLD, ENTC, ENTM)
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.
subroutine, public calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy and the flux across ...
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars_m.F90:35
integer(i4b), parameter rcl1
rcl1: Maximum length of record for input files (with factor 3 safety margin)
Initialisations for SICOPOLIS.
Definition: sico_init_m.F90:35
integer(i4b) ndata_gi
ndata_gi: Number of data values for the glacial index
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:35
integer(i4b), dimension(0:jmax, 0:imax) q_weert
q_weert(j,i): Weertman exponent for the basal pressure
subroutine, public calc_temp_melt()
Computation of the melting temperatures.
subroutine, public disc_param(dtime)
Ice discharge parameters (Greenland). [Assign ice discharge parameters.].
real(dp) ea
ea: Abbreviation for exp(aa)
real(dp), parameter pi_180
pi_180: pi divided by 180 (-> deg to rad)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) enth_c
enth_c(kc,j,i): Enthalpy in the upper (kc) ice domain
logical flag_aa_nonzero
flag_aa_nonzero: Flag for the exponential stretch parameter aa. .true.: aa greater than zero (non-equ...
Initial temperature, water content and age.
subroutine, public calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
Computation of the vertical velocity vz for grounded ice.
Definition: calc_vz_m.F90:52
real(dp) aa
aa: Exponential stretch parameter of the non-equidistant vertical grid in the upper (kc) ice domain ...
Reading of several input data.
Definition: read_m.F90:35
Declarations of kind types for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) zs
zs(j,i): Coordinate z of the surface topography
real(dp), dimension(0:jmax, 0:imax) precip_ma_present
precip_ma_present(j,i): Present-day mean annual precipitation rate at the ice surface ...
integer(i4b) gi_time_min
gi_time_min: Minimum time of the data values for the glacial index
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) temp_precip_time_min
temp_precip_time_min: Minimum time of the surface-temperature and precipitation data ...
subroutine, public stereo_forw_ellipsoid(lambda_val, phi_val, A, B, lambda0, phi0, x_val, y_val)
Forward stereographic projection for an ellipsoidal planet.
real(dp) kei_r_max
kei_r_max: Maximum value of the argument r of the tabulated kei function
real(dp) year_zero
year_zero: SICOPOLIS year zero in astronomical year numbering [ = signed year CE (AD) ] ...
real(dp), dimension(:), allocatable specmap_zsl
specmap_zsl(n): Data values for the sea level
Computation of the melting and basal temperatures.
subroutine, public disc_fields()
Dependence of ice discharge coefficient on latitude (Greenland). [Determine dependence of ice dischar...
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:201
real(dp), dimension(-190:10) kappa
KAPPA(n): Tabulated values for the heat conductivity of ice.
Computation of the vertical velocity vz.
Definition: calc_vz_m.F90:35
integer(i4b) gi_time_stp
gi_time_stp: Time step of the data values for the glacial index
integer(i4b) ndata_temp_precip
ndata_temp_precip: Number of surface-temperature and precipitation data
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
real(dp), dimension(:), allocatable temp_precip_time
temp_precip_time(n): Times of the surface-temperature and precipitation data
integer(i2b), dimension(0:jmax, 0:imax) maske_old
maske_old(j,i): Old value of maske (at the previous time step)
character(len=64) ch_domain_long
ch_domain_long: Long name of the computational domain
Definition of the components g11 and g22 of the metric tensor of the applied numerical coordinates...
Definition: metric_m.F90:37
real(dp), dimension(0:jmax, 0:imax) precip_ma_lgm_anom
precip_ma_lgm_anom(j,i): LGM anomaly (ratio LGM/present) of the mean annual precipitation rate at the...
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
integer(i4b) gi_time_max
gi_time_max: Maximum time of the data values for the glacial index
real(dp) phi0
PHI0: Standard parallel of the stereographic projection.
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.
real(dp), dimension(0:jmax, 0:imax, 12) gamma_precip_lgm_anom
gamma_precip_lgm_anom(j,i,n): negative natural logarithm of precip_lgm_anom(j,i,n) ...
Comparison of single- or double-precision floating-point numbers.
NetCDF error capturing.
Definition: nc_check_m.F90:35
real(dp) time_target_topo_final
time_target_topo_final: Final time for target-topography adjustment
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(0:jmax, 0:imax) temp_mj_lgm_anom
temp_mj_lgm_anom(j,i): LGM anomaly (difference LGM - present) of the mean summer (northern hemisphere...
real(dp), dimension(-190:10) rf
RF(n): Tabulated values for the rate factor of cold ice.
real(dp), dimension(0:jmax, 0:imax) q_b_tot
Q_b_tot(j,i): Sum of Q_bm and Q_tld.
real(dp), dimension(0:jmax, 0:imax) q_bm
Q_bm(j,i): Basal melting rate.
real(dp) function, public enth_fct_temp_omega(temp_val, omega_val)
Enthalpy as a function of temperature and water content.
subroutine, public calc_vxy_sia(dzeta_c, dzeta_t)
Computation of the shear stresses txz, tyz, the effective shear stress sigma, the depth-averaged flui...
Definition: calc_vxy_m.F90:478
real(dp), dimension(-190:10) c
C(n): Tabulated values for the specific heat of ice.
real(dp), dimension(0:kcmax) eaz_c_quotient
eaz_c_quotient(kc): Abbreviation for (eaz_c(kc)-1.0)/(ea-1.0)
real(dp), dimension(-10000:10000) kei
kei(n): Tabulated values of the kei function (Kelvin function of zero order)
subroutine, public init_temp_water_age_1_3()
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==3: present-day initial topogr...
subroutine boundary(time, dtime, dxi, deta, delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
Main routine of boundary_m: Computation of the surface temperature (must be less than 0 deg C!) and o...
Definition: boundary_m.F90:56
real(dp) time_target_topo_init
time_target_topo_init: Initial time for target-topography adjustment
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) ...
real(dp), dimension(0:jmax, 0:imax) temp_ma_lgm_anom
temp_ma_lgm_anom(j,i): LGM anomaly (difference LGM - present) of the mean annual surface temperature ...
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:3377
subroutine check(status, ch_calling_routine)
NetCDF error capturing.
Definition: nc_check_m.F90:46
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
integer(i4b) ncid_temp_precip
ncid_temp_precip: ID of the NetCDF file containing the surface-temperature and precipitation data as ...
character(len=16), dimension(:), allocatable ch_core
ch_core(n): Names of the prescribed borehole positions
real(dp), dimension(0:jmax, 0:imax) sq_g11_g
sq_g11_g(j,i): Square root of the coefficient g11 of the metric tensor on grid point (i...
real(dp), dimension(0:jmax, 0:imax) as_perp_apl
as_perp_apl(j,i): Applied accumulation-ablation function (SMB)
logical, dimension(0:jmax, 0:imax) flag_calving_front_1
flag_calving_front_1(j,i): Calving front flag. .true.: calving front point (floating ice point with a...
subroutine, public metric()
Main routine of module metric_m: Definition of the components g11 and g22 of the metric tensor of the...
Definition: metric_m.F90:54
subroutine, public calc_enhance_3(time)
Computation of the flow enhancement factor. Case ENHMOD==3: two different values depending on time of...
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
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:56
subroutine, public calc_enhance_4()
Computation of the flow enhancement factor. Case ENHMOD==4: minimal anisotropic enhancement factor fo...
real(dp) kei_r_incr
kei_r_incr: Increment of the argument r of the tabulated kei function
real(dp), dimension(0:ktmax) zeta_t
zeta_t(kt): Sigma coordinate zeta_t of grid point kt
real(dp) rho_sw
RHO_SW: Density of sea water.
logical, dimension(0:jmax, 0:imax) flag_calving_front_2
flag_calving_front_2(j,i): Calving front flag. .true.: calving front point (ocean point with at least...
Computation of the surface temperature (must be less than 0 deg C!) and of the accumulation-ablation ...
Definition: boundary_m.F90:37
real(dp) l
L: Latent heat of ice.
real(dp), parameter eps_sp_dp
eps_sp_dp: Small number to single-precision accuracy in double precision
integer(i4b) n_data_kei
n_data_kei: Number of tabulated values of the kei function
Ice discharge parameterization for the Greenland ice sheet.
real(dp) rho
RHO: Density of ice.
subroutine topography3(dxi, deta, z_sl, anfdatname)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
integer(i4b), dimension(0:jmax, 0:imax) p_weert
p_weert(j,i): Weertman exponent for the basal shear stress
subroutine, public init_temp_water_age_1_4()
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==4: present-day initial topogr...
subroutine, public init_temp_water_age_1_1()
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==1: present-day initial topogr...
real(dp), dimension(:), allocatable lambda_core
lambda_core(n): Geographical longitude of the prescribed borehole positions
integer(i4b) n_core
n_core: Number of positions to be considered in the time-series file for deep boreholes ...
integer, parameter dp
Double-precision reals.
subroutine topography2(dxi, deta)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
real(dp), dimension(0:jmax, 0:imax) zb
zb(j,i): Coordinate z of the ice base
logical, dimension(0:jmax, 0:imax) flag_grounding_line_2
flag_grounding_line_2(j,i): Grounding line flag. .true.: grounding line point (floating ice point wit...
real(dp), dimension(0:jmax, 0:imax, 12) precip_present
precip_present(j,i,n): Present-day mean monthly precipitation rate at the ice surface ...
Writing of output data on files.
Definition: output_m.F90:36
real(dp), dimension(0:jmax, 0:imax) dzb_dtau
dzb_dtau(j,i): Derivative of zb by tau (time)
real(dp), dimension(0:kcmax) zeta_c
zeta_c(kc): Sigma coordinate zeta_c of grid point kc
subroutine, public read_target_topo_nc(target_topo_dat_name)
Reading of the target-topography file (in NetCDF format).
Definition: read_m.F90:672
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...
real(dp), dimension(0:jmax, 0:imax, 12) precip_lgm_anom
precip_lgm_anom(j,i,n): LGM anomaly (ratio LGM/present) of the mean monthly precipitation rate at the...
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:53
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...