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