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 
287  errormsg = ' >>> sico_init: DX wrong!'
288  call error(errormsg)
289 
290 end if
291 
292 #elif (GRID==1)
293 
294 errormsg = ' >>> sico_init: GRID==1 not allowed for this application!'
295 call error(errormsg)
296 
297 #elif (GRID==2)
298 
299 errormsg = ' >>> sico_init: GRID==2 not allowed for this application!'
300 call error(errormsg)
301 
302 #endif
303 
304 !-------- Compatibility check of the thermodynamics mode
305 ! (cold vs. polythermal vs. enthalpy method)
306 ! and the number of grid points in the lower (kt) ice domain --------
307 
308 #if (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1)
309 
310 if (ktmax > 2) then
311  write(6, fmt='(a)') ' >>> sico_init: For options CALCMOD==0, 2, 3 or -1,'
312  write(6, fmt='(a)') ' the separate kt domain is redundant.'
313  write(6, fmt='(a)') ' Therefore, consider setting KTMAX to 2.'
314  write(6, fmt='(a)') ' '
315 end if
316 
317 #endif
318 
319 !-------- Compatibility check of discretization schemes for the horizontal and
320 ! vertical advection terms in the temperature and age equations --------
321 
322 #if (ADV_HOR==1)
323 errormsg = ' >>> sico_init: ' &
324  //'Option ADV_HOR==1 (central differences) not defined!'
325 call error(errormsg)
326 #endif
327 
328 !-------- Check whether for the shallow shelf
329 ! or shelfy stream approximation
330 ! the chosen grid is Cartesian coordinates
331 ! without distortion correction (GRID==0) --------
332 
333 #if ((MARGIN==3 && DYNAMICS==1) || DYNAMICS==2) /* requires SSA or SStA */
334 #if (GRID != 0)
335 errormsg = ' >>> sico_init: Distortion correction not implemented' &
336  // end_of_line &
337  //' for the shallow shelf approximation (SSA)' &
338  // end_of_line &
339  //' or the shelfy stream approximation (SStA)' &
340  // end_of_line &
341  //' -> GRID==0 required!'
342 call error(errormsg)
343 #endif
344 #endif
345 
346 !-------- Setting of forcing flag --------
347 
348 forcing_flag = 1 ! forcing by delta_ts
349 
350 !-------- Initialization of numerical time steps --------
351 
352 dtime0 = dtime0
353 dtime_temp0 = dtime_temp0
354 
355 !-------- Further initializations --------
356 
357 dzeta_c = 1.0_dp/real(kcmax,dp)
358 dzeta_t = 1.0_dp/real(ktmax,dp)
359 dzeta_r = 1.0_dp/real(krmax,dp)
360 
361 ndat2d = 1
362 ndat3d = 1
363 
364 !-------- General abbreviations --------
365 
366 ! ------ kc domain
367 
368 if (deform >= eps) then
369 
370  flag_aa_nonzero = .true. ! non-equidistant grid
371 
372  aa = deform
373  ea = exp(aa)
374 
375  kc=0
376  zeta_c(kc) = 0.0_dp
377  eaz_c(kc) = 1.0_dp
378  eaz_c_quotient(kc) = 0.0_dp
379 
380  do kc=1, kcmax-1
381  zeta_c(kc) = kc*dzeta_c
382  eaz_c(kc) = exp(aa*zeta_c(kc))
383  eaz_c_quotient(kc) = (eaz_c(kc)-1.0_dp)/(ea-1.0_dp)
384  end do
385 
386  kc=kcmax
387  zeta_c(kc) = 1.0_dp
388  eaz_c(kc) = exp(aa)
389  eaz_c_quotient(kc) = 1.0_dp
390 
391 else
392 
393  flag_aa_nonzero = .false. ! equidistant grid
394 
395  aa = 0.0_dp
396  ea = 1.0_dp
397 
398  kc=0
399  zeta_c(kc) = 0.0_dp
400  eaz_c(kc) = 1.0_dp
401  eaz_c_quotient(kc) = 0.0_dp
402 
403  do kc=1, kcmax-1
404  zeta_c(kc) = kc*dzeta_c
405  eaz_c(kc) = 1.0_dp
406  eaz_c_quotient(kc) = zeta_c(kc)
407  end do
408 
409  kc=kcmax
410  zeta_c(kc) = 1.0_dp
411  eaz_c(kc) = 1.0_dp
412  eaz_c_quotient(kc) = 1.0_dp
413 
414 end if
415 
416 ! ------ kt domain
417 
418 kt=0
419 zeta_t(kt) = 0.0_dp
420 
421 do kt=1, ktmax-1
422  zeta_t(kt) = kt*dzeta_t
423 end do
424 
425 kt=ktmax
426 zeta_t(kt) = 1.0_dp
427 
428 ! ------ kr domain
429 
430 kr=0
431 zeta_r(kr) = 0.0_dp
432 
433 do kr=1, krmax-1
434  zeta_r(kr) = kr*dzeta_r
435 end do
436 
437 kr=krmax
438 zeta_r(kr) = 1.0_dp
439 
440 !-------- Construction of a vector (with index n) from a 2-d array
441 ! (with indices i, j) by diagonal numbering --------
442 
443 n=1
444 
445 do m=0, imax+jmax
446  do i=m, 0, -1
447  j = m-i
448  if ((i <= imax).and.(j <= jmax)) then
449  ii(n) = i
450  jj(n) = j
451  nn(j,i) = n
452  n=n+1
453  end if
454  end do
455 end do
456 
457 !-------- Specification of current simulation --------
458 
459 runname = runname
460 anfdatname = anfdatname
461 
462 #if (defined(YEAR_ZERO))
464 #else
465 year_zero = 2000.0_dp ! default value 2000 CE
466 #endif
467 
468 time_init0 = time_init0
469 time_end0 = time_end0
470 dtime_ser0 = dtime_ser0
471 
472 #if (OUTPUT==1 || OUTPUT==3)
473 dtime_out0 = dtime_out0
474 #endif
475 
476 #if (OUTPUT==2 || OUTPUT==3)
477 n_output = n_output
478 time_output0( 1) = time_out0_01
479 time_output0( 2) = time_out0_02
480 time_output0( 3) = time_out0_03
481 time_output0( 4) = time_out0_04
482 time_output0( 5) = time_out0_05
483 time_output0( 6) = time_out0_06
484 time_output0( 7) = time_out0_07
485 time_output0( 8) = time_out0_08
486 time_output0( 9) = time_out0_09
487 time_output0(10) = time_out0_10
488 time_output0(11) = time_out0_11
489 time_output0(12) = time_out0_12
490 time_output0(13) = time_out0_13
491 time_output0(14) = time_out0_14
492 time_output0(15) = time_out0_15
493 time_output0(16) = time_out0_16
494 time_output0(17) = time_out0_17
495 time_output0(18) = time_out0_18
496 time_output0(19) = time_out0_19
497 time_output0(20) = time_out0_20
498 #endif
499 
500 !-------- Write log file --------
501 
502 shell_command = 'if [ ! -d'
503 shell_command = trim(shell_command)//' '//outpath
504 shell_command = trim(shell_command)//' '//'] ; then mkdir'
505 shell_command = trim(shell_command)//' '//outpath
506 shell_command = trim(shell_command)//' '//'; fi'
507 call system(trim(shell_command))
508  ! Check whether directory OUTPATH exists. If not, it is created.
509 
510 filename_with_path = trim(outpath)//'/'//trim(runname)//'.log'
511 
512 open(10, iostat=ios, file=trim(filename_with_path), status='new')
513 
514 if (ios /= 0) then
515  errormsg = ' >>> sico_init: Error when opening the log file!'
516  call error(errormsg)
517 end if
518 
519 write(10, fmt=trim(fmt1)) 'Computational domain:'
520 write(10, fmt=trim(fmt1)) trim(ch_domain_long)
521 write(10, fmt=trim(fmt1)) ' '
522 
523 write(10, fmt=trim(fmt2a)) 'GRID = ', grid
524 write(10, fmt=trim(fmt1)) ' '
525 
526 write(10, fmt=trim(fmt2)) 'imax =', imax
527 write(10, fmt=trim(fmt2)) 'jmax =', jmax
528 write(10, fmt=trim(fmt2)) 'kcmax =', kcmax
529 write(10, fmt=trim(fmt2)) 'ktmax =', ktmax
530 write(10, fmt=trim(fmt2)) 'krmax =', krmax
531 write(10, fmt=trim(fmt1)) ' '
532 
533 write(10, fmt=trim(fmt3)) 'a =', aa
534 write(10, fmt=trim(fmt1)) ' '
535 
536 #if (GRID==0 || GRID==1)
537 write(10, fmt=trim(fmt3)) 'x0 =', x0
538 write(10, fmt=trim(fmt3)) 'y0 =', y0
539 write(10, fmt=trim(fmt3)) 'dx =', dx
540 #elif (GRID==2)
541 write(10, fmt=trim(fmt3)) 'dlambda =', dlambda
542 write(10, fmt=trim(fmt3)) 'dphi =', dphi
543 #endif
544 write(10, fmt=trim(fmt1)) ' '
545 
546 write(10, fmt=trim(fmt3)) 'year_zero =', year_zero
547 write(10, fmt=trim(fmt3)) 'time_init =', time_init0
548 write(10, fmt=trim(fmt3)) 'time_end =', time_end0
549 write(10, fmt=trim(fmt3)) 'dtime =', dtime0
550 write(10, fmt=trim(fmt3)) 'dtime_temp =', dtime_temp0
551 write(10, fmt=trim(fmt1)) ' '
552 
553 write(10, fmt=trim(fmt2)) 'ANF_DAT =', anf_dat
554 #if (ANF_DAT==1 && defined(TEMP_INIT))
555 write(10, fmt=trim(fmt2)) 'TEMP_INIT =', temp_init
556 #endif
557 #if (ANF_DAT==3)
558 write(10, fmt=trim(fmt1)) 'Initial-value file = '//anfdatname
559 write(10, fmt=trim(fmt1)) 'Path to initial-value file = '//anfdatpath
560 #endif
561 write(10, fmt=trim(fmt1)) ' '
562 
563 write(10, fmt=trim(fmt1)) 'Physical-parameter file = '//phys_para_file
564 write(10, fmt=trim(fmt1)) ' '
565 
566 #if (CALCTHK==2 || CALCTHK==3 || CALCTHK==5 || CALCTHK==6)
567 write(10, fmt=trim(fmt3)) 'ovi_weight =', ovi_weight
568 #if (CALCTHK==2 || CALCTHK==5)
569 write(10, fmt=trim(fmt3)) 'omega_sor =', omega_sor
570 #if (ITER_MAX_SOR>0)
571 write(10, fmt=trim(fmt2a)) 'iter_max_sor = ', iter_max_sor
572 #endif
573 #endif
574 write(10, fmt=trim(fmt1)) ' '
575 #endif
576 
577 #if (defined(SURFACE_FORCING))
578 write(10, fmt=trim(fmt2)) 'SURFACE_FORCING =', surface_forcing
579 #endif
580 #if (!defined(SURFACE_FORCING) || SURFACE_FORCING==1)
581 write(10, fmt=trim(fmt3)) 'temp_min =', temp_min
582 write(10, fmt=trim(fmt3)) 's_t =', s_t
583 write(10, fmt=trim(fmt3)) 'x_hat =', x_hat
584 write(10, fmt=trim(fmt3)) 'y_hat =', y_hat
585 write(10, fmt=trim(fmt3)) 'b_max =', b_max
586 write(10, fmt=trim(fmt3)) 's_b =', s_b
587 write(10, fmt=trim(fmt3)) 'eld =', eld
588 #elif (SURFACE_FORCING==2)
589 write(10, fmt=trim(fmt3)) 'temp_0 =', temp_0
590 write(10, fmt=trim(fmt3)) 'gamma_t =', gamma_t
591 write(10, fmt=trim(fmt3)) 's_0 =', s_0
592 write(10, fmt=trim(fmt3)) 'm_0 =', m_0
593 write(10, fmt=trim(fmt3)) 'ela =', ela
594 #endif
595 write(10, fmt=trim(fmt1)) ' '
596 
597 #if (TSURFACE==1)
598 write(10, fmt=trim(fmt3)) 'delta_ts0 =', delta_ts0
599 #elif (TSURFACE==3)
600 write(10, fmt=trim(fmt3)) 'sine_amplit =', sine_amplit
601 write(10, fmt=trim(fmt3)) 'sine_period =', sine_period
602 #elif (TSURFACE==4)
603 write(10, fmt=trim(fmt1)) 'GRIP file = '//grip_temp_file
604 write(10, fmt=trim(fmt3)) 'grip_temp_fact =', grip_temp_fact
605 #endif
606 
607 write(10, fmt=trim(fmt2a)) 'SEA_LEVEL = ', sea_level
608 #if (SEA_LEVEL==1)
609 write(10, fmt=trim(fmt3)) 'z_sl0 =', z_sl0
610 #elif (SEA_LEVEL==3)
611 write(10, fmt=trim(fmt1)) 'sea-level file = '//sea_level_file
612 #endif
613 write(10, fmt=trim(fmt1)) ' '
614 
615 #if (MARGIN==2)
616 #if (MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3)
617 write(10, fmt=trim(fmt3)) 'z_mar =', z_mar
618 write(10, fmt=trim(fmt1)) ' '
619 #elif (MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \
620  || marine_ice_calving==6 || marine_ice_calving==7)
621 write(10, fmt=trim(fmt3)) 'fact_z_mar =', fact_z_mar
622 write(10, fmt=trim(fmt1)) ' '
623 #elif (MARINE_ICE_FORMATION==2 && MARINE_ICE_CALVING==9)
624 write(10, fmt=trim(fmt3)) 'calv_uw_coeff =', calv_uw_coeff
625 write(10, fmt=trim(fmt3)) 'r1_calv_uw =', r1_calv_uw
626 write(10, fmt=trim(fmt3)) 'r2_calv_uw =', r2_calv_uw
627 write(10, fmt=trim(fmt1)) ' '
628 #endif
629 #elif (MARGIN==3)
630 #if (ICE_SHELF_CALVING==2)
631 write(10, fmt=trim(fmt3)) 'H_calv =', h_calv
632 write(10, fmt=trim(fmt1)) ' '
633 #endif
634 #endif
635 
636 #if (defined(BASAL_HYDROLOGY))
637 write(10, fmt=trim(fmt2a)) 'BASAL_HYDROLOGY = ', basal_hydrology
638 #endif
639 write(10, fmt=trim(fmt2a)) 'SLIDE_LAW = ', slide_law
640 write(10, fmt=trim(fmt3)) 'c_slide =', c_slide
641 write(10, fmt=trim(fmt3)) 'gamma_slide =', gamma_slide
642 write(10, fmt=trim(fmt2a)) 'p_weert = ', p_weert
643 write(10, fmt=trim(fmt2a)) 'q_weert = ', q_weert
644 #if (defined(TIME_RAMP_UP_SLIDE))
645 write(10, fmt=trim(fmt3)) 'time_ramp_up_slide =', time_ramp_up_slide
646 #endif
647 #if (SLIDE_LAW==2)
648 write(10, fmt=trim(fmt3)) 'red_pres_limit_fact =', red_pres_limit_fact
649 #endif
650 #if (BASAL_HYDROLOGY==1 \
651  && defined(hydro_slide_sat_fct) \
652  && defined(c_hw_slide) && defined(hw0_slide))
653 write(10, fmt=trim(fmt2a)) 'HYDRO_SLIDE_SAT_FCT = ', hydro_slide_sat_fct
654 write(10, fmt=trim(fmt3)) 'c_Hw_slide =', c_hw_slide
655 write(10, fmt=trim(fmt3)) 'Hw0_slide =', hw0_slide
656 #endif
657 write(10, fmt=trim(fmt1)) ' '
658 
659 write(10, fmt=trim(fmt2a)) 'Q_GEO_MOD = ', q_geo_mod
660 #if (Q_GEO_MOD==1)
661 write(10, fmt=trim(fmt3)) 'q_geo =', q_geo
662 #endif
663 write(10, fmt=trim(fmt1)) ' '
664 
665 write(10, fmt=trim(fmt2)) 'REBOUND =', rebound
666 #if (REBOUND==1)
667 write(10, fmt=trim(fmt3)) 'frac_llra =', frac_llra
668 #endif
669 #if (REBOUND==1 || REBOUND==2)
670 write(10, fmt=trim(fmt2)) 'TIME_LAG_MOD =', time_lag_mod
671 #if (TIME_LAG_MOD==1)
672 write(10, fmt=trim(fmt3)) 'time_lag =', time_lag
673 #elif (TIME_LAG_MOD==2)
674 write(10, fmt=trim(fmt1)) 'time_lag_file = '//time_lag_file
675 #else
676 errormsg = ' >>> sico_init: TIME_LAG_MOD must be either 1 or 2!'
677 call error(errormsg)
678 #endif
679 #endif
680 #if (REBOUND==2)
681 write(10, fmt=trim(fmt2)) 'FLEX_RIG_MOD =', flex_rig_mod
682 #if (FLEX_RIG_MOD==1)
683 write(10, fmt=trim(fmt3)) 'flex_rig =', flex_rig
684 #elif (FLEX_RIG_MOD==2)
685 write(10, fmt=trim(fmt1)) 'flex_rig_file = '//flex_rig_file
686 #else
687 errormsg = ' >>> sico_init: FLEX_RIG_MOD must be either 1 or 2!'
688 call error(errormsg)
689 #endif
690 #endif
691 write(10, fmt=trim(fmt2)) 'Q_LITHO =', q_litho
692 write(10, fmt=trim(fmt1)) ' '
693 
694 #if (FLOW_LAW==2)
695 write(10, fmt=trim(fmt3)) 'gr_size =', gr_size
696 write(10, fmt=trim(fmt1)) ' '
697 #endif
698 #if (FIN_VISC==2)
699 write(10, fmt=trim(fmt3)) 'sigma_res =', sigma_res
700 write(10, fmt=trim(fmt1)) ' '
701 #endif
702 
703 write(10, fmt=trim(fmt2)) 'ENHMOD =', enhmod
704 #if (ENHMOD==1 || ENHMOD==2 || ENHMOD==3)
705 write(10, fmt=trim(fmt3)) 'enh_fact =', enh_fact
706 #endif
707 #if (ENHMOD==2 || ENHMOD==3)
708 write(10, fmt=trim(fmt3)) 'enh_intg =', enh_intg
709 #endif
710 #if (ENHMOD==2)
711 write(10, fmt=trim(fmt3)) 'age_trans =', age_trans_0
712 #endif
713 #if (ENHMOD==3)
714 write(10, fmt=trim(fmt3)) 'date_trans1 =', date_trans1_0
715 write(10, fmt=trim(fmt3)) 'date_trans2 =', date_trans2_0
716 write(10, fmt=trim(fmt3)) 'date_trans3 =', date_trans3_0
717 #endif
718 #if (ENHMOD==4 || ENHMOD==5)
719 write(10, fmt=trim(fmt3)) 'enh_compr =', enh_compr
720 write(10, fmt=trim(fmt3)) 'enh_shear =', enh_shear
721 #endif
722 #if ((ENHMOD==1 || ENHMOD==2 || ENHMOD==3 || ENHMOD==4) && MARGIN==3)
723 write(10, fmt=trim(fmt3)) 'enh_shelf =', enh_shelf
724 #endif
725 write(10, fmt=trim(fmt1)) ' '
726 
727 write(10, fmt=trim(fmt3)) 'numdiff_H_t =', numdiff_h_t
728 write(10, fmt=trim(fmt3)) 'tau_cts =', tau_cts
729 write(10, fmt=trim(fmt3)) 'vh_max =', vh_max
730 write(10, fmt=trim(fmt3)) 'hd_min =', hd_min
731 write(10, fmt=trim(fmt3)) 'hd_max =', hd_max
732 #if (defined(QBM_MIN))
733 write(10, fmt=trim(fmt3)) 'qbm_min =', qbm_min
734 #elif (defined(QB_MIN)) /* obsolete */
735 write(10, fmt=trim(fmt3)) 'qb_min =', qb_min
736 #endif
737 #if (defined(QBM_MAX))
738 write(10, fmt=trim(fmt3)) 'qbm_max =', qbm_max
739 #elif (defined(QB_MAX)) /* obsolete */
740 write(10, fmt=trim(fmt3)) 'qb_max =', qb_max
741 #endif
742 write(10, fmt=trim(fmt3)) 'age_min =', age_min
743 write(10, fmt=trim(fmt3)) 'age_max =', age_max
744 write(10, fmt=trim(fmt3)) 'mean_accum =', mean_accum
745 #if (ADV_VERT==1)
746 write(10, fmt=trim(fmt3)) 'age_diff =', agediff
747 #endif
748 write(10, fmt=trim(fmt1)) ' '
749 
750 write(10, fmt=trim(fmt2a)) 'DYNAMICS = ', dynamics
751 #if ((DYNAMICS==1 && MARGIN==3) || DYNAMICS==2)
752 #if (defined(LIS_OPTS))
753 write(10, fmt=trim(fmt1)) 'lis_opts = '//lis_opts
754 #endif
755 #if (defined(N_ITER_SSA))
756 write(10, fmt=trim(fmt2a)) 'n_iter_ssa = ', n_iter_ssa
757 #endif
758 #if (defined(RELAX_FACT_SSA))
759 write(10, fmt=trim(fmt3)) 'relax_fact_ssa =', relax_fact_ssa
760 #endif
761 #endif
762 #if (DYNAMICS==2 && defined(RATIO_SL_THRESH))
763 write(10, fmt=trim(fmt3)) 'ratio_sl_thresh =', ratio_sl_thresh
764 #endif
765 #if (DYNAMICS==2 && defined(SSTA_SIA_WEIGH_FCT))
766 write(10, fmt=trim(fmt2a)) 'SSTA_SIA_WEIGH_FCT = ', ssta_sia_weigh_fct
767 #endif
768 write(10, fmt=trim(fmt1)) ' '
769 
770 write(10, fmt=trim(fmt2)) 'CALCMOD =', calcmod
771 #if (CALCMOD==-1 && defined(TEMP_CONST))
772 write(10, fmt=trim(fmt3)) 'TEMP_CONST =', temp_const
773 #endif
774 #if (CALCMOD==-1 && defined(AGE_CONST))
775 write(10, fmt=trim(fmt3)) 'AGE_CONST =', age_const
776 #endif
777 #if (CALCMOD==1 && defined(CTS_MELTING_FREEZING))
778 write(10, fmt=trim(fmt2)) 'CTS_MELTING_FREEZING =', cts_melting_freezing
779 #endif
780 write(10, fmt=trim(fmt2)) 'FLOW_LAW =', flow_law
781 write(10, fmt=trim(fmt2)) 'FIN_VISC =', fin_visc
782 write(10, fmt=trim(fmt2)) 'MARGIN =', margin
783 #if (MARGIN==2)
784 write(10, fmt=trim(fmt2)) 'MARINE_ICE_FORMATION =', marine_ice_formation
785 write(10, fmt=trim(fmt2)) 'MARINE_ICE_CALVING =', marine_ice_calving
786 #elif (MARGIN==3)
787 write(10, fmt=trim(fmt2)) 'ICE_SHELF_CALVING =', ice_shelf_calving
788 #endif
789 #if (defined(THK_EVOL))
790 write(10, fmt=trim(fmt2)) 'THK_EVOL =', thk_evol
791 #else
792 errormsg = ' >>> sico_init: Define THK_EVOL in header file!'
793 call error(errormsg)
794 #endif
795 #if (defined(CALCTHK))
796 write(10, fmt=trim(fmt2)) 'CALCTHK =', calcthk
797 #else
798 errormsg = ' >>> sico_init: Define CALCTHK in header file!'
799 call error(errormsg)
800 #endif
801 write(10, fmt=trim(fmt2)) 'ADV_HOR =', adv_hor
802 write(10, fmt=trim(fmt2)) 'ADV_VERT =', adv_vert
803 write(10, fmt=trim(fmt2)) 'TOPOGRAD =', topograd
804 write(10, fmt=trim(fmt2)) 'TSURFACE =', tsurface
805 #if (defined(MB_ACCOUNT))
806 write(10, fmt=trim(fmt2a)) 'MB_ACCOUNT = ', mb_account
807 #endif
808 write(10, fmt=trim(fmt1)) ' '
809 
810 #if (defined(OUT_TIMES))
811 write(10, fmt=trim(fmt2a)) 'OUT_TIMES = ', out_times
812 #endif
813 #if (defined(OUTPUT_INIT))
814 write(10, fmt=trim(fmt2a)) 'OUTPUT_INIT = ', output_init
815 #endif
816 write(10, fmt=trim(fmt2a)) 'OUTPUT = ', output
817 #if (OUTPUT==1 || OUTPUT==3)
818 write(10, fmt=trim(fmt3)) 'dtime_out =' , dtime_out0
819 #endif
820 write(10, fmt=trim(fmt3)) 'dtime_ser =' , dtime_ser0
821 #if (OUTPUT==1 || OUTPUT==2)
822 write(10, fmt=trim(fmt2a)) 'ERGDAT = ', ergdat
823 #endif
824 write(10, fmt=trim(fmt2a)) 'NETCDF = ', netcdf
825 #if (OUTPUT==2 || OUTPUT==3)
826 write(10, fmt=trim(fmt2a)) 'n_output = ', n_output
827 do n=1, n_output
828  write(10, fmt=trim(fmt3)) 'time_output =' , time_output0(n)
829 end do
830 #endif
831 write(10, fmt=trim(fmt1)) ' '
832 
833 write(10, fmt=trim(fmt1)) 'Program version and date: '//version//' / '//date
834 
835 close(10, status='keep')
836 
837 !-------- Conversion of time quantities --------
838 
839 year_zero = year_zero*year_sec ! a --> s
840 time_init = time_init0*year_sec ! a --> s
841 time_end = time_end0*year_sec ! a --> s
842 dtime = dtime0*year_sec ! a --> s
843 dtime_temp = dtime_temp0*year_sec ! a --> s
844 dtime_ser = dtime_ser0*year_sec ! a --> s
845 #if (OUTPUT==1 || OUTPUT==3)
846 dtime_out = dtime_out0*year_sec ! a --> s
847 #endif
848 #if (OUTPUT==2 || OUTPUT==3)
849 do n=1, n_output
850  time_output(n) = time_output0(n)*year_sec ! a --> s
851 end do
852 #endif
853 
854 if (.not.approx_integer_multiple(dtime_temp, dtime, eps_sp_dp)) then
855  errormsg = ' >>> sico_init: dtime_temp must be a multiple of dtime!'
856  call error(errormsg)
857 end if
858 
859 if (.not.approx_integer_multiple(dtime_ser, dtime, eps_sp_dp)) then
860  errormsg = ' >>> sico_init: dtime_ser must be a multiple of dtime!'
861  call error(errormsg)
862 end if
863 
864 #if (OUTPUT==1 || OUTPUT==3)
865 if (.not.approx_integer_multiple(dtime_out, dtime, eps_sp_dp)) then
866  errormsg = ' >>> sico_init: dtime_out must be a multiple of dtime!'
867  call error(errormsg)
868 end if
869 #endif
870 
871 time = time_init
872 
873 !-------- Mean accumulation --------
874 
875 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho)
876 ! ! mm/a water equiv. --> m/s ice equiv.
877 
878 !-------- Read data for delta_ts --------
879 
880 #if (TSURFACE==4)
881 
882 filename_with_path = trim(inpath)//'/general/'//trim(grip_temp_file)
883 
884 open(21, iostat=ios, file=trim(filename_with_path), status='old')
885 
886 if (ios /= 0) then
887  errormsg = ' >>> sico_init: Error when opening the data file for delta_ts!'
888  call error(errormsg)
889 end if
890 
891 read(21, fmt=*) ch_dummy, grip_time_min, grip_time_stp, grip_time_max
892 
893 if (ch_dummy /= '#') then
894  errormsg = ' >>> sico_init: grip_time_min, grip_time_stp, grip_time_max' &
895  // end_of_line &
896  //' not defined in data file for delta_ts!'
897  call error(errormsg)
898 end if
899 
901 
902 allocate(griptemp(0:ndata_grip))
903 
904 do n=0, ndata_grip
905  read(21, fmt=*) d_dummy, griptemp(n)
906 end do
907 
908 close(21, status='keep')
909 
910 #endif
911 
912 !-------- Read data for z_sl --------
913 
914 #if (SEA_LEVEL==3)
915 
916 filename_with_path = trim(inpath)//'/general/'//trim(sea_level_file)
917 
918 open(21, iostat=ios, file=trim(filename_with_path), status='old')
919 
920 if (ios /= 0) then
921  errormsg = ' >>> sico_init: Error when opening the data file for z_sl!'
922  call error(errormsg)
923 end if
924 
925 read(21, fmt=*) ch_dummy, specmap_time_min, specmap_time_stp, specmap_time_max
926 
927 if (ch_dummy /= '#') then
928  errormsg = ' >>> sico_init:' &
929  // end_of_line &
930  //' specmap_time_min, specmap_time_stp, specmap_time_max' &
931  // end_of_line &
932  //' not defined in data file for z_sl!'
933  call error(errormsg)
934 end if
935 
937 
938 allocate(specmap_zsl(0:ndata_specmap))
939 
940 do n=0, ndata_specmap
941  read(21, fmt=*) d_dummy, specmap_zsl(n)
942 end do
943 
944 close(21, status='keep')
945 
946 #endif
947 
948 !-------- Determination of the geothermal heat flux --------
949 
950 #if (Q_GEO_MOD==1)
951 
952 ! ------ Constant value
953 
954 do i=0, imax
955 do j=0, jmax
956  q_geo(j,i) = q_geo *1.0e-03_dp ! mW/m2 -> W/m2
957 end do
958 end do
959 
960 #elif (Q_GEO_MOD==2)
961 
962 errormsg = ' >>> sico_init: ' &
963  //'Option Q_GEO_MOD==2 not available for this application!'
964 call error(errormsg)
965 
966 #endif
967 
968 !-------- Determination of the time lag
969 ! of the relaxing asthenosphere --------
970 
971 #if (REBOUND==1 || REBOUND==2)
972 
973 #if (TIME_LAG_MOD==1)
974 
975 time_lag_asth = time_lag*year_sec ! a -> s
976 
977 #elif (TIME_LAG_MOD==2)
978 
979 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
980  trim(time_lag_file)
981 
982 open(29, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
983 
984 if (ios /= 0) then
985  errormsg = ' >>> sico_init: Error when opening the time-lag file!'
986  call error(errormsg)
987 end if
988 
989 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
990 
991 do j=jmax, 0, -1
992  read(29, fmt=*) (time_lag_asth(j,i), i=0,imax)
993 end do
994 
995 close(29, status='keep')
996 
997 time_lag_asth = time_lag_asth*year_sec ! a -> s
998 
999 #endif
1000 
1001 #elif (REBOUND==0)
1002 
1003 time_lag_asth = 0.0_dp ! dummy values
1004 
1005 #endif
1006 
1007 !-------- Determination of the flexural rigidity of the lithosphere --------
1008 
1009 #if (REBOUND==2)
1010 
1011 #if (FLEX_RIG_MOD==1)
1012 
1013 flex_rig_lith = flex_rig
1014 
1015 #elif (FLEX_RIG_MOD==2)
1016 
1017 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1018  trim(flex_rig_file)
1019 
1020 open(29, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1021 
1022 if (ios /= 0) then
1023  errormsg = ' >>> sico_init: Error when opening the flex-rig file!'
1024  call error(errormsg)
1025 end if
1026 
1027 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1028 
1029 do j=jmax, 0, -1
1030  read(29, fmt=*) (flex_rig_lith(j,i), i=0,imax)
1031 end do
1032 
1033 close(29, status='keep')
1034 
1035 #endif
1036 
1037 #elif (REBOUND==0 || REBOUND==1)
1038 
1039 flex_rig_lith = 0.0_dp ! dummy values
1040 
1041 #endif
1042 
1043 !-------- Definition of initial values --------
1044 
1045 ! ------ Present topography
1046 
1047 #if (ANF_DAT==1)
1048 
1049 errormsg = ' >>> sico_init: ANF_DAT==1 not allowed for this application!'
1050 call error(errormsg)
1051 
1052 ! ------ Ice-free, relaxed bedrock
1053 
1054 #elif (ANF_DAT==2)
1055 
1056 call topography2(dxi, deta)
1057 
1058 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1059 
1060 call boundary(time_init, dtime, dxi, deta, &
1061  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1062 
1063 as_perp_apl = 0.0_dp
1064 
1065 smb_corr = 0.0_dp
1066 
1067 q_bm = 0.0_dp
1068 q_tld = 0.0_dp
1069 q_b_tot = 0.0_dp
1070 
1071 p_b_w = 0.0_dp
1072 h_w = 0.0_dp
1073 
1074 call init_temp_water_age_2()
1075 
1076 #if (ENHMOD==1)
1077  call calc_enhance_1()
1078 #elif (ENHMOD==2)
1079  call calc_enhance_2()
1080 #elif (ENHMOD==3)
1081  call calc_enhance_3(time_init)
1082 #elif (ENHMOD==4)
1083  call calc_enhance_4()
1084 #elif (ENHMOD==5)
1085  call calc_enhance_5()
1086 #else
1087  errormsg = ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
1088  call error(errormsg)
1089 #endif
1090 
1091 ! ------ Read initial state from output of previous simulation
1092 
1093 #elif (ANF_DAT==3)
1094 
1095 call topography3(dxi, deta, z_sl, anfdatname)
1096 
1097 call boundary(time_init, dtime, dxi, deta, &
1098  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1099 
1100 where ((maske==0_i2b).or.(maske==3_i2b))
1101  ! grounded or floating ice
1103 elsewhere ! maske==1_i2b or 2_i2b, ice-free land or sea
1104  as_perp_apl = 0.0_dp
1105 end where
1106 
1107 smb_corr = 0.0_dp
1108 
1109 q_b_tot = q_bm + q_tld
1110 
1111 #endif
1112 
1113 !-------- Inner-point flag --------
1114 
1115 flag_inner_point = .true.
1116 
1117 flag_inner_point(0,:) = .false.
1118 flag_inner_point(jmax,:) = .false.
1119 
1120 flag_inner_point(:,0) = .false.
1121 flag_inner_point(:,imax) = .false.
1122 
1123 !-------- Grounding line and calving front flags --------
1124 
1125 flag_grounding_line_1 = .false.
1126 flag_grounding_line_2 = .false.
1127 
1128 flag_calving_front_1 = .false.
1129 flag_calving_front_2 = .false.
1130 
1131 #if (MARGIN==1 || MARGIN==2)
1132 
1133 continue
1134 
1135 #elif (MARGIN==3)
1136 
1137 do i=1, imax-1
1138 do j=1, jmax-1
1139 
1140  if ( (maske(j,i)==0_i2b) & ! grounded ice
1141  .and. &
1142  ( (maske(j,i+1)==3_i2b) & ! with
1143  .or.(maske(j,i-1)==3_i2b) & ! one
1144  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1145  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1146  ) &
1147  flag_grounding_line_1(j,i) = .true.
1148 
1149  if ( (maske(j,i)==3_i2b) & ! floating ice
1150  .and. &
1151  ( (maske(j,i+1)==0_i2b) & ! with
1152  .or.(maske(j,i-1)==0_i2b) & ! one
1153  .or.(maske(j+1,i)==0_i2b) & ! neighbouring
1154  .or.(maske(j-1,i)==0_i2b) ) & ! grounded ice point
1155  ) &
1156  flag_grounding_line_2(j,i) = .true.
1157 
1158  if ( (maske(j,i)==3_i2b) & ! floating ice
1159  .and. &
1160  ( (maske(j,i+1)==2_i2b) & ! with
1161  .or.(maske(j,i-1)==2_i2b) & ! one
1162  .or.(maske(j+1,i)==2_i2b) & ! neighbouring
1163  .or.(maske(j-1,i)==2_i2b) ) & ! ocean point
1164  ) &
1165  flag_calving_front_1(j,i) = .true.
1166 
1167  if ( (maske(j,i)==2_i2b) & ! ocean
1168  .and. &
1169  ( (maske(j,i+1)==3_i2b) & ! with
1170  .or.(maske(j,i-1)==3_i2b) & ! one
1171  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1172  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1173  ) &
1174  flag_calving_front_2(j,i) = .true.
1175 
1176 end do
1177 end do
1178 
1179 do i=1, imax-1
1180 
1181  j=0
1182  if ( (maske(j,i)==2_i2b) & ! ocean
1183  .and. (maske(j+1,i)==3_i2b) & ! with one neighbouring
1184  ) & ! floating ice point
1185  flag_calving_front_2(j,i) = .true.
1186 
1187  j=jmax
1188  if ( (maske(j,i)==2_i2b) & ! ocean
1189  .and. (maske(j-1,i)==3_i2b) & ! with one neighbouring
1190  ) & ! floating ice point
1191  flag_calving_front_2(j,i) = .true.
1192 
1193 end do
1194 
1195 do j=1, jmax-1
1196 
1197  i=0
1198  if ( (maske(j,i)==2_i2b) & ! ocean
1199  .and. (maske(j,i+1)==3_i2b) & ! with one neighbouring
1200  ) & ! floating ice point
1201  flag_calving_front_2(j,i) = .true.
1202 
1203  i=imax
1204  if ( (maske(j,i)==2_i2b) & ! ocean
1205  .and. (maske(j,i-1)==3_i2b) & ! with one neighbouring
1206  ) & ! floating ice point
1207  flag_calving_front_2(j,i) = .true.
1208 
1209 end do
1210 
1211 #else
1212 
1213 errormsg = ' >>> sico_init: MARGIN must be either 1, 2 or 3!'
1214 call error(errormsg)
1215 
1216 #endif
1217 
1218 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
1219 
1220 #if (GRID==0 || GRID==1)
1221 
1222 do ir=-imax, imax
1223 do jr=-jmax, jmax
1224  dist_dxdy(jr,ir) = sqrt( (real(ir,dp)*dxi)**2 + (real(jr,dp)*deta)**2 )
1225  ! distortion due to stereographic projection not accounted for
1226 end do
1227 end do
1228 
1229 #endif
1230 
1231 !-------- Initial velocities --------
1232 
1233 call calc_temp_melt()
1234 
1235 #if (DYNAMICS==1 || DYNAMICS==2)
1236 
1237 call calc_vxy_b_sia(time, z_sl)
1238 call calc_vxy_sia(dzeta_c, dzeta_t)
1239 
1240 #if (MARGIN==3 || DYNAMICS==2)
1241 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
1242 #endif
1243 
1244 call calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
1245 
1246 #if (MARGIN==3)
1247 call calc_vz_floating(z_sl, dxi, deta, dzeta_c)
1248 #endif
1249 
1250 #elif (DYNAMICS==0)
1251 
1252 call calc_vxy_static()
1253 call calc_vz_static()
1254 
1255 #else
1256 
1257 errormsg = ' >>> sico_init: DYNAMICS must be either 0, 1 or 2!'
1258 call error(errormsg)
1259 
1260 #endif
1261 
1262 call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
1263 
1264 !-------- Initial enthalpies --------
1265 
1266 #if (CALCMOD==0 || CALCMOD==-1)
1267 
1268 do i=0, imax
1269 do j=0, jmax
1270 
1271  do kc=0, kcmax
1272  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1273  end do
1274 
1275  do kt=0, ktmax
1276  enth_t(kt,j,i) = enth_c(0,j,i)
1277  end do
1278 
1279 end do
1280 end do
1281 
1282 #elif (CALCMOD==1)
1283 
1284 do i=0, imax
1285 do j=0, jmax
1286 
1287  do kc=0, kcmax
1288  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1289  end do
1290 
1291  if ( (maske(j,i) == 0_i2b).and.(n_cts(j,i) == 1_i2b) ) then
1292  do kt=0, ktmax
1293  enth_t(kt,j,i) = enth_fct_temp_omega(temp_t_m(kt,j,i), omega_t(kt,j,i))
1294  end do
1295  else
1296  do kt=0, ktmax
1297  enth_t(kt,j,i) = enth_c(0,j,i)
1298  end do
1299  end if
1300 
1301 end do
1302 end do
1303 
1304 #elif (CALCMOD==2 || CALCMOD==3)
1305 
1306 do i=0, imax
1307 do j=0, jmax
1308 
1309  do kc=0, kcmax
1310  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), omega_c(kc,j,i))
1311  end do
1312 
1313  do kt=0, ktmax
1314  enth_t(kt,j,i) = enth_c(0,j,i)
1315  end do
1316 
1317 end do
1318 end do
1319 
1320 #else
1321 
1322 errormsg = ' >>> sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!'
1323 call error(errormsg)
1324 
1325 #endif
1326 
1327 !-------- Initialize time-series files --------
1328 
1329 ! ------ Time-series file for the ice sheet on the whole
1330 
1331 filename_with_path = trim(outpath)//'/'//trim(runname)//'.ser'
1332 
1333 open(12, iostat=ios, file=trim(filename_with_path), status='new')
1334 
1335 if (ios /= 0) then
1336  errormsg = ' >>> sico_init: Error when opening the ser file!'
1337  call error(errormsg)
1338 end if
1339 
1340 write(12,1102)
1341 write(12,1103)
1342 
1343  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1344  ' V(m^3) V_g(m^3) V_f(m^3)', &
1345  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1346  ' V_sle(m) V_t(m^3)', &
1347  ' A_t(m^2)',/, &
1348  ' H_max(m) H_t_max(m)', &
1349  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
1350  1103 format('----------------------------------------------------', &
1351  '---------------------------------------')
1352 
1353 ! ------ Time-series file for selected positions ("deep boreholes")
1354 
1355 n_core = 2 ! central dome, position halfway to coast
1356 
1357 allocate(lambda_core(n_core), phi_core(n_core), &
1359 
1360 ch_core(1) = 'P1'
1361 lambda_core(1) = 0.0_dp ! dummy
1362 phi_core(1) = 0.0_dp ! dummy
1363 x_core(1) = 750.0_dp *1.0e+03_dp ! Position of the central dome
1364 y_core(1) = 750.0_dp *1.0e+03_dp ! (750 km, 750 km),
1365  ! conversion km -> m
1366 
1367 ch_core(2) = 'P2'
1368 lambda_core(2) = 0.0_dp ! dummy
1369 phi_core(2) = 0.0_dp ! dummy
1370 x_core(2) = 750.0_dp *1.0e+03_dp ! Position halfway to the coast
1371 y_core(2) = 1125.0_dp *1.0e+03_dp ! (750 km, 1125 km),
1372  ! conversion km -> m
1373 
1374 filename_with_path = trim(outpath)//'/'//trim(runname)//'.core'
1375 
1376 open(14, iostat=ios, file=trim(filename_with_path), status='new')
1377 
1378 if (ios /= 0) then
1379  errormsg = ' >>> sico_init: Error when opening the core file!'
1380  call error(errormsg)
1381 end if
1382 
1383 if (forcing_flag == 1) then
1384 
1385  write(14,1106)
1386  write(14,1107)
1387 
1388  1106 format(' t(a) D_Ts(C) z_sl(m)',/, &
1389  ' H_P1(m) H_P2(m)',/, &
1390  ' v_P1(m/a) v_P2(m/a)',/, &
1391  ' T_P1(C) T_P2(C)')
1392  1107 format('---------------------------------------')
1393 
1394 end if
1395 
1396 !-------- Output of the initial state --------
1397 
1398 #if (defined(OUTPUT_INIT))
1399 
1400 #if (OUTPUT_INIT==0)
1401  flag_init_output = .false.
1402 #elif (OUTPUT_INIT==1)
1403  flag_init_output = .true.
1404 #else
1405  errormsg = ' >>> sico_init: OUTPUT_INIT must be either 0 or 1!'
1406  call error(errormsg)
1407 #endif
1408 
1409 #else
1410  flag_init_output = .true. ! default for undefined parameter OUTPUT_INIT
1411 #endif
1412 
1413 #if (OUTPUT==1)
1414 
1415 #if (ERGDAT==0)
1416  flag_3d_output = .false.
1417 #elif (ERGDAT==1)
1418  flag_3d_output = .true.
1419 #else
1420  errormsg = ' >>> sico_init: ERGDAT must be either 0 or 1!'
1421  call error(errormsg)
1422 #endif
1423 
1424  if (flag_init_output) &
1425  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1426  flag_3d_output, ndat2d, ndat3d)
1427 
1428 #elif (OUTPUT==2)
1429 
1430 if (time_output(1) <= time_init+eps) then
1431 
1432 #if (ERGDAT==0)
1433  flag_3d_output = .false.
1434 #elif (ERGDAT==1)
1435  flag_3d_output = .true.
1436 #else
1437  errormsg = ' >>> sico_init: ERGDAT must be either 0 or 1!'
1438  call error(errormsg)
1439 #endif
1440 
1441  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1442  flag_3d_output, ndat2d, ndat3d)
1443 
1444 end if
1445 
1446 #elif (OUTPUT==3)
1447 
1448  flag_3d_output = .false.
1449 
1450  if (flag_init_output) &
1451  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1452  flag_3d_output, ndat2d, ndat3d)
1453 
1454 if (time_output(1) <= time_init+eps) then
1455 
1456  flag_3d_output = .true.
1457 
1458  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1459  flag_3d_output, ndat2d, ndat3d)
1460 
1461 end if
1462 
1463 #else
1464 
1465  errormsg = ' >>> sico_init: OUTPUT must be either 1, 2 or 3!'
1466  call error(errormsg)
1467 
1468 #endif
1469 
1470 if (flag_init_output) then
1471  call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1472  call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1473 end if
1474 
1475 end subroutine sico_init
1476 
1477 !-------------------------------------------------------------------------------
1478 !> Definition of the initial surface and bedrock topography
1479 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1480 !! For ice-free initial topography with relaxed lithosphere
1481 !! (not defined for this domain, thus routine stops execution of SICOPOLIS).
1482 !<------------------------------------------------------------------------------
1483 subroutine topography1(dxi, deta)
1484 
1485 implicit none
1486 
1487 real(dp), intent(out) :: dxi, deta
1488 
1489 ! integer(i4b) :: i, j
1490 ! integer(i4b) :: ios
1491 ! real(dp) :: xi0, eta0
1492 
1493 dxi=0.0_dp; deta=0.0_dp ! dummy values
1494 
1495 errormsg = ' >>> topography1: ANF_DAT==1 not defined for this domain!'
1496 call error(errormsg)
1497 
1498 end subroutine topography1
1499 
1500 !-------------------------------------------------------------------------------
1501 !> Definition of the initial surface and bedrock topography
1502 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1503 !! For ice-free initial topography with relaxed lithosphere.
1504 !<------------------------------------------------------------------------------
1505 subroutine topography2(dxi, deta)
1506 
1507 #if (GRID==0 || GRID==1)
1509 #endif
1510 
1511  use metric_m
1512  use topograd_m
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 !-------- Set topography --------
1523 
1524 zl0 = 0.0_dp
1525 
1526 maske = 1
1527 
1528 maske(0,:) = 2
1529 maske(jmax,:) = 2
1530 
1531 maske(:,0) = 2
1532 maske(:,imax) = 2
1533 
1534 !-------- Further stuff --------
1535 
1536 #if (GRID==0 || GRID==1)
1537 
1538 dxi = dx *1000.0_dp ! km -> m
1539 deta = dx *1000.0_dp ! km -> m
1540 
1541 xi0 = x0 *1000.0_dp ! km -> m
1542 eta0 = y0 *1000.0_dp ! km -> m
1543 
1544 #elif (GRID==2)
1545 
1546 errormsg = ' >>> topography2: GRID==2 not allowed for this application!'
1547 call error(errormsg)
1548 
1549 #endif
1550 
1551 do i=0, imax
1552 do j=0, jmax
1553 
1554  zs(j,i) = zl0(j,i)
1555  zb(j,i) = zl0(j,i)
1556  zl(j,i) = zl0(j,i)
1557 
1558  xi(i) = xi0 + real(i,dp)*dxi
1559  eta(j) = eta0 + real(j,dp)*deta
1560 
1561  zm(j,i) = zb(j,i)
1562  n_cts(j,i) = -1
1563  kc_cts(j,i) = 0
1564 
1565  h_c(j,i) = 0.0_dp
1566  h_t(j,i) = 0.0_dp
1567 
1568  dzs_dtau(j,i) = 0.0_dp
1569  dzm_dtau(j,i) = 0.0_dp
1570  dzb_dtau(j,i) = 0.0_dp
1571  dzl_dtau(j,i) = 0.0_dp
1572  dh_c_dtau(j,i) = 0.0_dp
1573  dh_t_dtau(j,i) = 0.0_dp
1574 
1575 end do
1576 end do
1577 
1578 maske_old = maske
1579 
1580 !-------- Geographic coordinates, metric tensor,
1581 ! gradients of the topography --------
1582 
1583 do i=0, imax
1584 do j=0, jmax
1585 
1586 #if (GRID==0 || GRID==1) /* Stereographic projection */
1587 
1588  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1589  lambda0, phi0, lambda(j,i), phi(j,i))
1590 
1591 #elif (GRID==2) /* Geographic coordinates */
1592 
1593  lambda(j,i) = xi(i)
1594  phi(j,i) = eta(j)
1595 
1596 #endif
1597 
1598 end do
1599 end do
1600 
1601 call metric()
1602 
1603 #if (TOPOGRAD==0)
1604 call topograd_1(dxi, deta, 1)
1605 #elif (TOPOGRAD==1)
1606 call topograd_2(dxi, deta, 1)
1607 #endif
1608 
1609 !-------- Corresponding area of grid points --------
1610 
1611 do i=0, imax
1612 do j=0, jmax
1613  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1614 end do
1615 end do
1616 
1617 end subroutine topography2
1618 
1619 !-------------------------------------------------------------------------------
1620 !> Definition of the initial surface and bedrock topography
1621 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1622 !! For initial topography from previous simulation.
1623 !<------------------------------------------------------------------------------
1624 subroutine topography3(dxi, deta, z_sl, anfdatname)
1625 
1626  use read_m, only : read_erg_nc
1627 
1628 #if (GRID==0 || GRID==1)
1630 #endif
1631 
1632  use metric_m
1633  use topograd_m
1634 
1635 implicit none
1636 
1637 character(len=100), intent(in) :: anfdatname
1638 
1639 real(dp), intent(out) :: dxi, deta, z_sl
1640 
1641 integer(i4b) :: i, j
1642 
1643 !-------- Read data from time-slice file of previous simulation --------
1644 
1645 call read_erg_nc(z_sl, anfdatname)
1646 
1647 !-------- Set topography of the relaxed bedrock --------
1648 
1649 zl0 = 0.0_dp
1650 
1651 !-------- Further stuff --------
1652 
1653 #if (GRID==0 || GRID==1)
1654 
1655 dxi = dx *1000.0_dp ! km -> m
1656 deta = dx *1000.0_dp ! km -> m
1657 
1658 #elif (GRID==2)
1659 
1660 errormsg = ' >>> topography3: GRID==2 not allowed for this application!'
1661 call error(errormsg)
1662 
1663 #endif
1664 
1665 !-------- Geographic coordinates, metric tensor,
1666 ! gradients of the topography --------
1667 
1668 do i=0, imax
1669 do j=0, jmax
1670 
1671 #if (GRID==0 || GRID==1) /* Stereographic projection */
1672 
1673  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1674  lambda0, phi0, lambda(j,i), phi(j,i))
1675 
1676 #elif (GRID==2) /* Geographic coordinates */
1677 
1678  lambda(j,i) = xi(i)
1679  phi(j,i) = eta(j)
1680 
1681 #endif
1682 
1683 end do
1684 end do
1685 
1686 call metric()
1687 
1688 #if (TOPOGRAD==0)
1689 call topograd_1(dxi, deta, 1)
1690 #elif (TOPOGRAD==1)
1691 call topograd_2(dxi, deta, 1)
1692 #endif
1693 
1694 !-------- Corresponding area of grid points --------
1695 
1696 do i=0, imax
1697 do j=0, jmax
1698  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1699 end do
1700 end do
1701 
1702 end subroutine topography3
1703 
1704 !-------------------------------------------------------------------------------
1705 
1706 end module sico_init_m
1707 !
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...