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