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