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