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