SICOPOLIS V5-dev  Revision 1173
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 smb_corr = 0.0_dp
1004 
1005 q_bm = 0.0_dp
1006 q_tld = 0.0_dp
1007 q_b_tot = 0.0_dp
1008 
1009 p_b_w = 0.0_dp
1010 h_w = 0.0_dp
1011 
1012 call init_temp_water_age_2()
1013 
1014 #if (ENHMOD==1)
1015  call calc_enhance_1()
1016 #elif (ENHMOD==2)
1017  call calc_enhance_2()
1018 #elif (ENHMOD==3)
1019  call calc_enhance_3(time_init)
1020 #elif (ENHMOD==4)
1021  call calc_enhance_4()
1022 #elif (ENHMOD==5)
1023  call calc_enhance_5()
1024 #else
1025  stop ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
1026 #endif
1027 
1028 ! ------ Read initial state from output of previous simulation
1029 
1030 #elif (ANF_DAT==3)
1031 
1032 call topography3(dxi, deta, z_sl, anfdatname)
1033 
1034 call boundary(time_init, dtime, dxi, deta, &
1035  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1036 
1037 where ((maske==0_i2b).or.(maske==3_i2b))
1038  ! grounded or floating ice
1040 elsewhere ! maske==1_i2b or 2_i2b, ice-free land or sea
1041  as_perp_apl = 0.0_dp
1042 end where
1043 
1044 smb_corr = 0.0_dp
1045 
1046 q_b_tot = q_bm + q_tld
1047 
1048 #endif
1049 
1050 !-------- Inner-point flag --------
1051 
1052 flag_inner_point = .true.
1053 
1054 flag_inner_point(0,:) = .false.
1055 flag_inner_point(jmax,:) = .false.
1056 
1057 flag_inner_point(:,0) = .false.
1058 flag_inner_point(:,imax) = .false.
1059 
1060 !-------- Grounding line and calving front flags --------
1061 
1062 flag_grounding_line_1 = .false.
1063 flag_grounding_line_2 = .false.
1064 
1065 flag_calving_front_1 = .false.
1066 flag_calving_front_2 = .false.
1067 
1068 #if (MARGIN==1 || MARGIN==2)
1069 
1070 continue
1071 
1072 #elif (MARGIN==3)
1073 
1074 do i=1, imax-1
1075 do j=1, jmax-1
1076 
1077  if ( (maske(j,i)==0_i2b) & ! grounded ice
1078  .and. &
1079  ( (maske(j,i+1)==3_i2b) & ! with
1080  .or.(maske(j,i-1)==3_i2b) & ! one
1081  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1082  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1083  ) &
1084  flag_grounding_line_1(j,i) = .true.
1085 
1086  if ( (maske(j,i)==3_i2b) & ! floating ice
1087  .and. &
1088  ( (maske(j,i+1)==0_i2b) & ! with
1089  .or.(maske(j,i-1)==0_i2b) & ! one
1090  .or.(maske(j+1,i)==0_i2b) & ! neighbouring
1091  .or.(maske(j-1,i)==0_i2b) ) & ! grounded ice point
1092  ) &
1093  flag_grounding_line_2(j,i) = .true.
1094 
1095  if ( (maske(j,i)==3_i2b) & ! floating ice
1096  .and. &
1097  ( (maske(j,i+1)==2_i2b) & ! with
1098  .or.(maske(j,i-1)==2_i2b) & ! one
1099  .or.(maske(j+1,i)==2_i2b) & ! neighbouring
1100  .or.(maske(j-1,i)==2_i2b) ) & ! ocean point
1101  ) &
1102  flag_calving_front_1(j,i) = .true.
1103 
1104  if ( (maske(j,i)==2_i2b) & ! ocean
1105  .and. &
1106  ( (maske(j,i+1)==3_i2b) & ! with
1107  .or.(maske(j,i-1)==3_i2b) & ! one
1108  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1109  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1110  ) &
1111  flag_calving_front_2(j,i) = .true.
1112 
1113 end do
1114 end do
1115 
1116 do i=1, imax-1
1117 
1118  j=0
1119  if ( (maske(j,i)==2_i2b) & ! ocean
1120  .and. (maske(j+1,i)==3_i2b) & ! with one neighbouring
1121  ) & ! floating ice point
1122  flag_calving_front_2(j,i) = .true.
1123 
1124  j=jmax
1125  if ( (maske(j,i)==2_i2b) & ! ocean
1126  .and. (maske(j-1,i)==3_i2b) & ! with one neighbouring
1127  ) & ! floating ice point
1128  flag_calving_front_2(j,i) = .true.
1129 
1130 end do
1131 
1132 do j=1, jmax-1
1133 
1134  i=0
1135  if ( (maske(j,i)==2_i2b) & ! ocean
1136  .and. (maske(j,i+1)==3_i2b) & ! with one neighbouring
1137  ) & ! floating ice point
1138  flag_calving_front_2(j,i) = .true.
1139 
1140  i=imax
1141  if ( (maske(j,i)==2_i2b) & ! ocean
1142  .and. (maske(j,i-1)==3_i2b) & ! with one neighbouring
1143  ) & ! floating ice point
1144  flag_calving_front_2(j,i) = .true.
1145 
1146 end do
1147 
1148 #else
1149 stop ' >>> sico_init: MARGIN must be either 1, 2 or 3!'
1150 #endif
1151 
1152 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
1153 
1154 #if (GRID==0 || GRID==1)
1155 
1156 do ir=-imax, imax
1157 do jr=-jmax, jmax
1158  dist_dxdy(jr,ir) = sqrt( (real(ir,dp)*dxi)**2 + (real(jr,dp)*deta)**2 )
1159  ! distortion due to stereographic projection not accounted for
1160 end do
1161 end do
1162 
1163 #endif
1164 
1165 !-------- Initial velocities --------
1166 
1167 call calc_temp_melt()
1168 
1169 #if (DYNAMICS==1 || DYNAMICS==2)
1170 
1171 call calc_vxy_b_sia(time, z_sl)
1172 call calc_vxy_sia(dzeta_c, dzeta_t)
1173 
1174 #if (MARGIN==3 || DYNAMICS==2)
1175 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
1176 #endif
1177 
1178 call calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
1179 
1180 #if (MARGIN==3)
1181 call calc_vz_floating(z_sl, dxi, deta, dzeta_c)
1182 #endif
1183 
1184 #elif (DYNAMICS==0)
1185 
1186 call calc_vxy_static()
1187 call calc_vz_static()
1188 
1189 #else
1190  stop ' >>> sico_init: DYNAMICS must be either 0, 1 or 2!'
1191 #endif
1192 
1193 call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
1194 
1195 !-------- Initial enthalpies --------
1196 
1197 #if (CALCMOD==0 || CALCMOD==-1)
1198 
1199 do i=0, imax
1200 do j=0, jmax
1201 
1202  do kc=0, kcmax
1203  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1204  end do
1205 
1206  do kt=0, ktmax
1207  enth_t(kt,j,i) = enth_c(0,j,i)
1208  end do
1209 
1210 end do
1211 end do
1212 
1213 #elif (CALCMOD==1)
1214 
1215 do i=0, imax
1216 do j=0, jmax
1217 
1218  do kc=0, kcmax
1219  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1220  end do
1221 
1222  if ( (maske(j,i) == 0_i2b).and.(n_cts(j,i) == 1_i2b) ) then
1223  do kt=0, ktmax
1224  enth_t(kt,j,i) = enth_fct_temp_omega(temp_t_m(kt,j,i), omega_t(kt,j,i))
1225  end do
1226  else
1227  do kt=0, ktmax
1228  enth_t(kt,j,i) = enth_c(0,j,i)
1229  end do
1230  end if
1231 
1232 end do
1233 end do
1234 
1235 #elif (CALCMOD==2 || CALCMOD==3)
1236 
1237 do i=0, imax
1238 do j=0, jmax
1239 
1240  do kc=0, kcmax
1241  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), omega_c(kc,j,i))
1242  end do
1243 
1244  do kt=0, ktmax
1245  enth_t(kt,j,i) = enth_c(0,j,i)
1246  end do
1247 
1248 end do
1249 end do
1250 
1251 #else
1252 
1253 stop ' >>> sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!'
1254 
1255 #endif
1256 
1257 !-------- Initialize time-series files --------
1258 
1259 ! ------ Time-series file for the ice sheet on the whole
1260 
1261 filename_with_path = trim(outpath)//'/'//trim(runname)//'.ser'
1262 
1263 open(12, iostat=ios, file=trim(filename_with_path), status='new')
1264 
1265 if (ios /= 0) stop ' >>> sico_init: Error when opening the ser file!'
1266 
1267 write(12,1102)
1268 write(12,1103)
1269 
1270  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1271  ' V(m^3) V_g(m^3) V_f(m^3)', &
1272  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1273  ' V_sle(m) V_t(m^3)', &
1274  ' A_t(m^2)',/, &
1275  ' H_max(m) H_t_max(m)', &
1276  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
1277  1103 format('----------------------------------------------------', &
1278  '---------------------------------------')
1279 
1280 ! ------ Time-series file for selected positions ("deep boreholes")
1281 
1282 n_core = 2 ! central dome, position halfway to coast
1283 
1284 allocate(lambda_core(n_core), phi_core(n_core), &
1286 
1287 ch_core(1) = 'P1'
1288 lambda_core(1) = 0.0_dp ! dummy
1289 phi_core(1) = 0.0_dp ! dummy
1290 x_core(1) = 750.0_dp *1.0e+03_dp ! Position of the central dome
1291 y_core(1) = 750.0_dp *1.0e+03_dp ! (750 km, 750 km),
1292  ! conversion km -> m
1293 
1294 ch_core(2) = 'P2'
1295 lambda_core(2) = 0.0_dp ! dummy
1296 phi_core(2) = 0.0_dp ! dummy
1297 x_core(2) = 750.0_dp *1.0e+03_dp ! Position halfway to the coast
1298 y_core(2) = 1125.0_dp *1.0e+03_dp ! (750 km, 1125 km),
1299  ! conversion km -> m
1300 
1301 filename_with_path = trim(outpath)//'/'//trim(runname)//'.core'
1302 
1303 open(14, iostat=ios, file=trim(filename_with_path), status='new')
1304 
1305 if (ios /= 0) stop ' >>> sico_init: Error when opening the core file!'
1306 
1307 if (forcing_flag == 1) then
1308 
1309  write(14,1106)
1310  write(14,1107)
1311 
1312  1106 format(' t(a) D_Ts(C) z_sl(m)',/, &
1313  ' H_P1(m) H_P2(m)',/, &
1314  ' v_P1(m/a) v_P2(m/a)',/, &
1315  ' T_P1(C) T_P2(C)')
1316  1107 format('---------------------------------------')
1317 
1318 end if
1319 
1320 !-------- Output of the initial state --------
1321 
1322 #if (defined(OUTPUT_INIT))
1323 
1324 #if (OUTPUT_INIT==0)
1325  flag_init_output = .false.
1326 #elif (OUTPUT_INIT==1)
1327  flag_init_output = .true.
1328 #else
1329  stop ' >>> sico_init: OUTPUT_INIT must be either 0 or 1!'
1330 #endif
1331 
1332 #else
1333  flag_init_output = .true. ! default for undefined parameter OUTPUT_INIT
1334 #endif
1335 
1336 #if (OUTPUT==1)
1337 
1338 #if (ERGDAT==0)
1339  flag_3d_output = .false.
1340 #elif (ERGDAT==1)
1341  flag_3d_output = .true.
1342 #else
1343  stop ' >>> sico_init: ERGDAT must be either 0 or 1!'
1344 #endif
1345 
1346  if (flag_init_output) &
1347  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1348  flag_3d_output, ndat2d, ndat3d)
1349 
1350 #elif (OUTPUT==2)
1351 
1352 if (time_output(1) <= time_init+eps) then
1353 
1354 #if (ERGDAT==0)
1355  flag_3d_output = .false.
1356 #elif (ERGDAT==1)
1357  flag_3d_output = .true.
1358 #else
1359  stop ' >>> sico_init: ERGDAT must be either 0 or 1!'
1360 #endif
1361 
1362  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1363  flag_3d_output, ndat2d, ndat3d)
1364 
1365 end if
1366 
1367 #elif (OUTPUT==3)
1368 
1369  flag_3d_output = .false.
1370 
1371  if (flag_init_output) &
1372  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1373  flag_3d_output, ndat2d, ndat3d)
1374 
1375 if (time_output(1) <= time_init+eps) then
1376 
1377  flag_3d_output = .true.
1378 
1379  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1380  flag_3d_output, ndat2d, ndat3d)
1381 
1382 end if
1383 
1384 #else
1385  stop ' >>> sico_init: OUTPUT must be either 1, 2 or 3!'
1386 #endif
1387 
1388 if (flag_init_output) then
1389  call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1390  call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1391 end if
1392 
1393 end subroutine sico_init
1394 
1395 !-------------------------------------------------------------------------------
1396 !> Definition of the initial surface and bedrock topography
1397 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1398 !! For ice-free initial topography with relaxed lithosphere
1399 !! (not defined for this domain, thus routine stops execution of SICOPOLIS).
1400 !<------------------------------------------------------------------------------
1401 subroutine topography1(dxi, deta)
1402 
1403 implicit none
1404 
1405 real(dp), intent(out) :: dxi, deta
1406 
1407 ! integer(i4b) :: i, j
1408 ! integer(i4b) :: ios
1409 ! real(dp) :: xi0, eta0
1410 
1411 dxi=0.0_dp; deta=0.0_dp ! dummy values
1412 
1413 write(6,'(/1x,a)') '>>> topography1: ANF_DAT==1 not defined for this domain!'
1414 stop
1415 
1416 end subroutine topography1
1417 
1418 !-------------------------------------------------------------------------------
1419 !> Definition of the initial surface and bedrock topography
1420 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1421 !! For ice-free initial topography with relaxed lithosphere.
1422 !<------------------------------------------------------------------------------
1423 subroutine topography2(dxi, deta)
1424 
1425 #if (GRID==0 || GRID==1)
1427 #endif
1428 
1429  use metric_m
1430  use topograd_m
1431 
1432 implicit none
1433 
1434 real(dp), intent(out) :: dxi, deta
1435 
1436 integer(i4b) :: i, j
1437 integer(i4b) :: ios
1438 real(dp) :: xi0, eta0
1439 
1440 !-------- Set topography --------
1441 
1442 zl0 = 0.0_dp
1443 
1444 maske = 1
1445 
1446 maske(0,:) = 2
1447 maske(jmax,:) = 2
1448 
1449 maske(:,0) = 2
1450 maske(:,imax) = 2
1451 
1452 !-------- Further stuff --------
1453 
1454 #if (GRID==0 || GRID==1)
1455 
1456 dxi = dx *1000.0_dp ! km -> m
1457 deta = dx *1000.0_dp ! km -> m
1458 
1459 xi0 = x0 *1000.0_dp ! km -> m
1460 eta0 = y0 *1000.0_dp ! km -> m
1461 
1462 #elif (GRID==2)
1463 
1464 stop ' >>> topography2: GRID==2 not allowed for this application!'
1465 
1466 #endif
1467 
1468 do i=0, imax
1469 do j=0, jmax
1470 
1471  zs(j,i) = zl0(j,i)
1472  zb(j,i) = zl0(j,i)
1473  zl(j,i) = zl0(j,i)
1474 
1475  xi(i) = xi0 + real(i,dp)*dxi
1476  eta(j) = eta0 + real(j,dp)*deta
1477 
1478  zm(j,i) = zb(j,i)
1479  n_cts(j,i) = -1
1480  kc_cts(j,i) = 0
1481 
1482  h_c(j,i) = 0.0_dp
1483  h_t(j,i) = 0.0_dp
1484 
1485  dzs_dtau(j,i) = 0.0_dp
1486  dzm_dtau(j,i) = 0.0_dp
1487  dzb_dtau(j,i) = 0.0_dp
1488  dzl_dtau(j,i) = 0.0_dp
1489  dh_c_dtau(j,i) = 0.0_dp
1490  dh_t_dtau(j,i) = 0.0_dp
1491 
1492 end do
1493 end do
1494 
1495 maske_old = maske
1496 
1497 !-------- Geographic coordinates, metric tensor,
1498 ! gradients of the topography --------
1499 
1500 do i=0, imax
1501 do j=0, jmax
1502 
1503 #if (GRID==0 || GRID==1) /* Stereographic projection */
1504 
1505  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1506  lambda0, phi0, lambda(j,i), phi(j,i))
1507 
1508 #elif (GRID==2) /* Geographic coordinates */
1509 
1510  lambda(j,i) = xi(i)
1511  phi(j,i) = eta(j)
1512 
1513 #endif
1514 
1515 end do
1516 end do
1517 
1518 call metric()
1519 
1520 #if (TOPOGRAD==0)
1521 call topograd_1(dxi, deta, 1)
1522 #elif (TOPOGRAD==1)
1523 call topograd_2(dxi, deta, 1)
1524 #endif
1525 
1526 !-------- Corresponding area of grid points --------
1527 
1528 do i=0, imax
1529 do j=0, jmax
1530  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1531 end do
1532 end do
1533 
1534 end subroutine topography2
1535 
1536 !-------------------------------------------------------------------------------
1537 !> Definition of the initial surface and bedrock topography
1538 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1539 !! For initial topography from previous simulation.
1540 !<------------------------------------------------------------------------------
1541 subroutine topography3(dxi, deta, z_sl, anfdatname)
1542 
1543  use read_m, only : read_erg_nc
1544 
1545 #if (GRID==0 || GRID==1)
1547 #endif
1548 
1549  use metric_m
1550  use topograd_m
1551 
1552 implicit none
1553 
1554 character(len=100), intent(in) :: anfdatname
1555 
1556 real(dp), intent(out) :: dxi, deta, z_sl
1557 
1558 integer(i4b) :: i, j
1559 
1560 !-------- Read data from time-slice file of previous simulation --------
1561 
1562 call read_erg_nc(z_sl, anfdatname)
1563 
1564 !-------- Set topography of the relaxed bedrock --------
1565 
1566 zl0 = 0.0_dp
1567 
1568 !-------- Further stuff --------
1569 
1570 #if (GRID==0 || GRID==1)
1571 
1572 dxi = dx *1000.0_dp ! km -> m
1573 deta = dx *1000.0_dp ! km -> m
1574 
1575 #elif (GRID==2)
1576 
1577 stop ' >>> topography3: GRID==2 not allowed for this application!'
1578 
1579 #endif
1580 
1581 !-------- Geographic coordinates, metric tensor,
1582 ! gradients of the topography --------
1583 
1584 do i=0, imax
1585 do j=0, jmax
1586 
1587 #if (GRID==0 || GRID==1) /* Stereographic projection */
1588 
1589  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1590  lambda0, phi0, lambda(j,i), phi(j,i))
1591 
1592 #elif (GRID==2) /* Geographic coordinates */
1593 
1594  lambda(j,i) = xi(i)
1595  phi(j,i) = eta(j)
1596 
1597 #endif
1598 
1599 end do
1600 end do
1601 
1602 call metric()
1603 
1604 #if (TOPOGRAD==0)
1605 call topograd_1(dxi, deta, 1)
1606 #elif (TOPOGRAD==1)
1607 call topograd_2(dxi, deta, 1)
1608 #endif
1609 
1610 !-------- Corresponding area of grid points --------
1611 
1612 do i=0, imax
1613 do j=0, jmax
1614  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1615 end do
1616 end do
1617 
1618 end subroutine topography3
1619 
1620 !-------------------------------------------------------------------------------
1621 
1622 end module sico_init_m
1623 !
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:823
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:4751
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), 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: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:3407
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...