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