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 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 smb_corr = 0.0_dp
999 
1000 q_bm = 0.0_dp
1001 q_tld = 0.0_dp
1002 q_b_tot = 0.0_dp
1003 
1004 p_b_w = 0.0_dp
1005 h_w = 0.0_dp
1006 
1007 #if (!defined(TEMP_INIT) || TEMP_INIT==1)
1009 #elif (TEMP_INIT==2)
1011 #elif (TEMP_INIT==3)
1013 #elif (TEMP_INIT==4)
1015 #else
1016  write(6, fmt='(a)') ' >>> sico_init:'
1017  write(6, fmt='(a)') ' TEMP_INIT must be set to either 1, 2, 3 or 4!'
1018  stop
1019 #endif
1020 
1021 #if (ENHMOD==1)
1022  call calc_enhance_1()
1023 #elif (ENHMOD==2)
1024  call calc_enhance_2()
1025 #elif (ENHMOD==3)
1026  call calc_enhance_3(time_init)
1027 #elif (ENHMOD==4)
1028  call calc_enhance_4()
1029 #elif (ENHMOD==5)
1030  call calc_enhance_5()
1031 #else
1032  stop ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
1033 #endif
1034 
1035 ! ------ Ice-free, relaxed bedrock
1036 
1037 #elif (ANF_DAT==2)
1038 
1039 call topography2(dxi, deta)
1040 
1041 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1042 
1043 call boundary(time_init, dtime, dxi, deta, &
1044  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1045 
1046 as_perp_apl = 0.0_dp
1047 
1048 smb_corr = 0.0_dp
1049 
1050 q_bm = 0.0_dp
1051 q_tld = 0.0_dp
1052 q_b_tot = 0.0_dp
1053 
1054 p_b_w = 0.0_dp
1055 h_w = 0.0_dp
1056 
1057 call init_temp_water_age_2()
1058 
1059 #if (ENHMOD==1)
1060  call calc_enhance_1()
1061 #elif (ENHMOD==2)
1062  call calc_enhance_2()
1063 #elif (ENHMOD==3)
1064  call calc_enhance_3(time_init)
1065 #elif (ENHMOD==4)
1066  call calc_enhance_4()
1067 #elif (ENHMOD==5)
1068  call calc_enhance_5()
1069 #else
1070  stop ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
1071 #endif
1072 
1073 ! ------ Read initial state from output of previous simulation
1074 
1075 #elif (ANF_DAT==3)
1076 
1077 call topography3(dxi, deta, z_sl, anfdatname)
1078 
1079 call boundary(time_init, dtime, dxi, deta, &
1080  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1081 
1082 where ((maske==0_i2b).or.(maske==3_i2b))
1083  ! grounded or floating ice
1085 elsewhere ! maske==1_i2b or 2_i2b, ice-free land or sea
1086  as_perp_apl = 0.0_dp
1087 end where
1088 
1089 smb_corr = 0.0_dp
1090 
1091 q_b_tot = q_bm + q_tld
1092 
1093 #endif
1094 
1095 !-------- Inner-point flag --------
1096 
1097 flag_inner_point = .true.
1098 
1099 flag_inner_point(0,:) = .false.
1100 flag_inner_point(jmax,:) = .false.
1101 
1102 flag_inner_point(:,0) = .false.
1103 flag_inner_point(:,imax) = .false.
1104 
1105 !-------- Grounding line and calving front flags --------
1106 
1107 flag_grounding_line_1 = .false.
1108 flag_grounding_line_2 = .false.
1109 
1110 flag_calving_front_1 = .false.
1111 flag_calving_front_2 = .false.
1112 
1113 #if (MARGIN==1 || MARGIN==2)
1114 
1115 continue
1116 
1117 #elif (MARGIN==3)
1118 
1119 do i=1, imax-1
1120 do j=1, jmax-1
1121 
1122  if ( (maske(j,i)==0_i2b) & ! grounded ice
1123  .and. &
1124  ( (maske(j,i+1)==3_i2b) & ! with
1125  .or.(maske(j,i-1)==3_i2b) & ! one
1126  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1127  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1128  ) &
1129  flag_grounding_line_1(j,i) = .true.
1130 
1131  if ( (maske(j,i)==3_i2b) & ! floating ice
1132  .and. &
1133  ( (maske(j,i+1)==0_i2b) & ! with
1134  .or.(maske(j,i-1)==0_i2b) & ! one
1135  .or.(maske(j+1,i)==0_i2b) & ! neighbouring
1136  .or.(maske(j-1,i)==0_i2b) ) & ! grounded ice point
1137  ) &
1138  flag_grounding_line_2(j,i) = .true.
1139 
1140  if ( (maske(j,i)==3_i2b) & ! floating ice
1141  .and. &
1142  ( (maske(j,i+1)==2_i2b) & ! with
1143  .or.(maske(j,i-1)==2_i2b) & ! one
1144  .or.(maske(j+1,i)==2_i2b) & ! neighbouring
1145  .or.(maske(j-1,i)==2_i2b) ) & ! ocean point
1146  ) &
1147  flag_calving_front_1(j,i) = .true.
1148 
1149  if ( (maske(j,i)==2_i2b) & ! ocean
1150  .and. &
1151  ( (maske(j,i+1)==3_i2b) & ! with
1152  .or.(maske(j,i-1)==3_i2b) & ! one
1153  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1154  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1155  ) &
1156  flag_calving_front_2(j,i) = .true.
1157 
1158 end do
1159 end do
1160 
1161 do i=1, imax-1
1162 
1163  j=0
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  j=jmax
1170  if ( (maske(j,i)==2_i2b) & ! ocean
1171  .and. (maske(j-1,i)==3_i2b) & ! with one neighbouring
1172  ) & ! floating ice point
1173  flag_calving_front_2(j,i) = .true.
1174 
1175 end do
1176 
1177 do j=1, jmax-1
1178 
1179  i=0
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  i=imax
1186  if ( (maske(j,i)==2_i2b) & ! ocean
1187  .and. (maske(j,i-1)==3_i2b) & ! with one neighbouring
1188  ) & ! floating ice point
1189  flag_calving_front_2(j,i) = .true.
1190 
1191 end do
1192 
1193 #else
1194 stop ' >>> sico_init: MARGIN must be either 1, 2 or 3!'
1195 #endif
1196 
1197 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
1198 
1199 #if (GRID==0 || GRID==1)
1200 
1201 do ir=-imax, imax
1202 do jr=-jmax, jmax
1203  dist_dxdy(jr,ir) = sqrt( (real(ir,dp)*dxi)**2 + (real(jr,dp)*deta)**2 )
1204  ! distortion due to stereographic projection not accounted for
1205 end do
1206 end do
1207 
1208 #endif
1209 
1210 !-------- Initial velocities --------
1211 
1212 call calc_temp_melt()
1213 
1214 #if (DYNAMICS==1 || DYNAMICS==2)
1215 
1216 call calc_vxy_b_sia(time, z_sl)
1217 call calc_vxy_sia(dzeta_c, dzeta_t)
1218 
1219 #if (MARGIN==3 || DYNAMICS==2)
1220 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
1221 #endif
1222 
1223 call calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
1224 
1225 #if (MARGIN==3)
1226 call calc_vz_floating(z_sl, dxi, deta, dzeta_c)
1227 #endif
1228 
1229 #elif (DYNAMICS==0)
1230 
1231 call calc_vxy_static()
1232 call calc_vz_static()
1233 
1234 #else
1235  stop ' >>> sico_init: DYNAMICS must be either 0, 1 or 2!'
1236 #endif
1237 
1238 call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
1239 
1240 !-------- Initial enthalpies --------
1241 
1242 #if (CALCMOD==0 || CALCMOD==-1)
1243 
1244 do i=0, imax
1245 do j=0, jmax
1246 
1247  do kc=0, kcmax
1248  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1249  end do
1250 
1251  do kt=0, ktmax
1252  enth_t(kt,j,i) = enth_c(0,j,i)
1253  end do
1254 
1255 end do
1256 end do
1257 
1258 #elif (CALCMOD==1)
1259 
1260 do i=0, imax
1261 do j=0, jmax
1262 
1263  do kc=0, kcmax
1264  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1265  end do
1266 
1267  if ( (maske(j,i) == 0_i2b).and.(n_cts(j,i) == 1_i2b) ) then
1268  do kt=0, ktmax
1269  enth_t(kt,j,i) = enth_fct_temp_omega(temp_t_m(kt,j,i), omega_t(kt,j,i))
1270  end do
1271  else
1272  do kt=0, ktmax
1273  enth_t(kt,j,i) = enth_c(0,j,i)
1274  end do
1275  end if
1276 
1277 end do
1278 end do
1279 
1280 #elif (CALCMOD==2 || CALCMOD==3)
1281 
1282 do i=0, imax
1283 do j=0, jmax
1284 
1285  do kc=0, kcmax
1286  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), omega_c(kc,j,i))
1287  end do
1288 
1289  do kt=0, ktmax
1290  enth_t(kt,j,i) = enth_c(0,j,i)
1291  end do
1292 
1293 end do
1294 end do
1295 
1296 #else
1297 
1298 stop ' >>> sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!'
1299 
1300 #endif
1301 
1302 !-------- Initialize time-series files --------
1303 
1304 ! ------ Time-series file for the ice sheet on the whole
1305 
1306 filename_with_path = trim(outpath)//'/'//trim(runname)//'.ser'
1307 
1308 open(12, iostat=ios, file=trim(filename_with_path), status='new')
1309 
1310 if (ios /= 0) stop ' >>> sico_init: Error when opening the ser file!'
1311 
1312 write(12,1102)
1313 write(12,1103)
1314 
1315  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1316  ' V(m^3) V_g(m^3) V_f(m^3)', &
1317  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1318  ' V_sle(m) V_t(m^3)', &
1319  ' A_t(m^2)',/, &
1320  ' H_max(m) H_t_max(m)', &
1321  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
1322  1103 format('----------------------------------------------------', &
1323  '---------------------------------------')
1324 
1325 ! ------ Time-series file for the sediment area ("Hudson Bay, Hudson Strait")
1326 
1327 filename_with_path = trim(outpath)//'/'//trim(runname)//'.sed'
1328 
1329 open(15, iostat=ios, file=trim(filename_with_path), status='new')
1330 
1331 if (ios /= 0) stop ' >>> sico_init: Error when opening the sed file!'
1332 
1333 write(15,1108)
1334 write(15,1109)
1335 
1336 1108 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1337  ' H_ave(m) Tbh_ave(C) Atb(m^2)')
1338 1109 format('----------------------------------------------------')
1339 
1340 ! ------ Time-series file for selected positions ("deep boreholes")
1341 
1342 n_core = 7 ! Points P1 - P7
1343 
1344 allocate(lambda_core(n_core), phi_core(n_core), &
1346 
1347 ch_core(1) = 'P1'
1348 lambda_core(1) = 0.0_dp ! dummy
1349 phi_core(1) = 0.0_dp ! dummy
1350 x_core(1) = 3900.0_dp *1.0e+03_dp ! Point P1,
1351 y_core(1) = 2000.0_dp *1.0e+03_dp ! conversion km -> m
1352 
1353 ch_core(2) = 'P2'
1354 lambda_core(2) = 0.0_dp ! dummy
1355 phi_core(2) = 0.0_dp ! dummy
1356 x_core(2) = 3800.0_dp *1.0e+03_dp ! Point P2,
1357 y_core(2) = 2000.0_dp *1.0e+03_dp ! conversion km -> m
1358 
1359 ch_core(3) = 'P3'
1360 lambda_core(3) = 0.0_dp ! dummy
1361 phi_core(3) = 0.0_dp ! dummy
1362 x_core(3) = 3700.0_dp *1.0e+03_dp ! Point P3,
1363 y_core(3) = 2000.0_dp *1.0e+03_dp ! conversion km -> m
1364 
1365 ch_core(4) = 'P4'
1366 lambda_core(4) = 0.0_dp ! dummy
1367 phi_core(4) = 0.0_dp ! dummy
1368 x_core(4) = 3500.0_dp *1.0e+03_dp ! Point P4,
1369 y_core(4) = 2000.0_dp *1.0e+03_dp ! conversion km -> m
1370 
1371 ch_core(5) = 'P5'
1372 lambda_core(5) = 0.0_dp ! dummy
1373 phi_core(5) = 0.0_dp ! dummy
1374 x_core(5) = 3200.0_dp *1.0e+03_dp ! Point P5,
1375 y_core(5) = 2000.0_dp *1.0e+03_dp ! conversion km -> m
1376 
1377 ch_core(6) = 'P6'
1378 lambda_core(6) = 0.0_dp ! dummy
1379 phi_core(6) = 0.0_dp ! dummy
1380 x_core(6) = 2900.0_dp *1.0e+03_dp ! Point P6,
1381 y_core(6) = 2000.0_dp *1.0e+03_dp ! conversion km -> m
1382 
1383 ch_core(7) = 'P7'
1384 lambda_core(7) = 0.0_dp ! dummy
1385 phi_core(7) = 0.0_dp ! dummy
1386 x_core(7) = 2600.0_dp *1.0e+03_dp ! Point P7,
1387 y_core(7) = 2000.0_dp *1.0e+03_dp ! conversion km -> m
1388 
1389 filename_with_path = trim(outpath)//'/'//trim(runname)//'.core'
1390 
1391 open(14, iostat=ios, file=trim(filename_with_path), status='new')
1392 
1393 if (ios /= 0) stop ' >>> sico_init: Error when opening the core file!'
1394 
1395 if (forcing_flag == 1) then
1396 
1397  write(14,1106)
1398  write(14,1107)
1399 
1400  1106 format(' t(a) D_Ts(C) z_sl(m)',/, &
1401  ' H_P1(m) H_P2(m) H_P3(m) H_P4(m)', &
1402  ' H_P5(m) H_P6(m) H_P7(m)',/, &
1403  ' v_P1(m/a) v_P2(m/a) v_P3(m/a) v_P4(m/a)', &
1404  ' v_P5(m/a) v_P6(m/a) v_P7(m/a)',/, &
1405  ' T_P1(C) T_P2(C) T_P3(C) T_P4(C)', &
1406  ' T_P5(C) T_P6(C) T_P7(C)')
1407  1107 format('-----------------------------------------------------------------', &
1408  '---------------------------------------')
1409 
1410 end if
1411 
1412 !-------- Output of the initial state --------
1413 
1414 #if (defined(OUTPUT_INIT))
1415 
1416 #if (OUTPUT_INIT==0)
1417  flag_init_output = .false.
1418 #elif (OUTPUT_INIT==1)
1419  flag_init_output = .true.
1420 #else
1421  stop ' >>> sico_init: OUTPUT_INIT must be either 0 or 1!'
1422 #endif
1423 
1424 #else
1425  flag_init_output = .true. ! default for undefined parameter OUTPUT_INIT
1426 #endif
1427 
1428 #if (OUTPUT==1)
1429 
1430 #if (ERGDAT==0)
1431  flag_3d_output = .false.
1432 #elif (ERGDAT==1)
1433  flag_3d_output = .true.
1434 #else
1435  stop ' >>> sico_init: ERGDAT must be either 0 or 1!'
1436 #endif
1437 
1438  if (flag_init_output) &
1439  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1440  flag_3d_output, ndat2d, ndat3d)
1441 
1442 #elif (OUTPUT==2)
1443 
1444 if (time_output(1) <= time_init+eps) then
1445 
1446 #if (ERGDAT==0)
1447  flag_3d_output = .false.
1448 #elif (ERGDAT==1)
1449  flag_3d_output = .true.
1450 #else
1451  stop ' >>> sico_init: ERGDAT must be either 0 or 1!'
1452 #endif
1453 
1454  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1455  flag_3d_output, ndat2d, ndat3d)
1456 
1457 end if
1458 
1459 #elif (OUTPUT==3)
1460 
1461  flag_3d_output = .false.
1462 
1463  if (flag_init_output) &
1464  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1465  flag_3d_output, ndat2d, ndat3d)
1466 
1467 if (time_output(1) <= time_init+eps) then
1468 
1469  flag_3d_output = .true.
1470 
1471  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1472  flag_3d_output, ndat2d, ndat3d)
1473 
1474 end if
1475 
1476 #else
1477  stop ' >>> sico_init: OUTPUT must be either 1, 2 or 3!'
1478 #endif
1479 
1480 if (flag_init_output) then
1481  call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1482  call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1483 end if
1484 
1485 end subroutine sico_init
1486 
1487 !-------------------------------------------------------------------------------
1488 !> Definition of the initial surface and bedrock topography
1489 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1490 !! For an initial topography with a very thin ice layer everywhere on the
1491 !! land area.
1492 !<------------------------------------------------------------------------------
1493 subroutine topography1(dxi, deta)
1494 
1495 #if (GRID==0 || GRID==1)
1497 #endif
1498 
1499  use metric_m
1500  use topograd_m
1501 
1502 implicit none
1503 
1504 real(dp), intent(out) :: dxi, deta
1505 
1506 integer(i4b) :: i, j, n
1507 integer(i4b) :: ios, n_dummy
1508 real(dp) :: d_dummy
1509 real(dp) :: xi0, eta0
1510 character :: ch_dummy
1511 
1512 real(dp), parameter :: h_init = 1.0e-03_dp ! initial ice thickness
1513 
1514 character(len= 8) :: ch_imax
1515 character(len=128) :: fmt4
1516 character(len=256) :: filename_with_path
1517 
1518 write(ch_imax, fmt='(i8)') imax
1519 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
1520 
1521 !-------- Set topography --------
1522 
1523 zl0 = 0.0_dp
1524 
1525 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1526  trim(mask_present_file)
1527 
1528 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
1529 
1530 if (ios /= 0) stop ' >>> topography1: Error when opening the mask file!'
1531 
1532 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1533 
1534 do j=jmax, 0, -1
1535  read(24, fmt=trim(fmt4)) (maske(j,i), i=0,imax)
1536 end do
1537 
1538 close(24, status='keep')
1539 
1540 !-------- Further stuff --------
1541 
1542 #if (GRID==0 || GRID==1)
1543 
1544 dxi = dx *1000.0_dp ! km -> m
1545 deta = dx *1000.0_dp ! km -> m
1546 
1547 xi0 = x0 *1000.0_dp ! km -> m
1548 eta0 = y0 *1000.0_dp ! km -> m
1549 
1550 #elif (GRID==2)
1551 
1552 stop ' >>> topography1: GRID==2 not allowed for this application!'
1553 
1554 #endif
1555 
1556 do i=0, imax
1557 do j=0, jmax
1558 
1559  zs(j,i) = zl0(j,i)
1560  zb(j,i) = zl0(j,i)
1561  zl(j,i) = zl0(j,i)
1562 
1563  if (maske(j,i) <= 1) then
1564  maske(j,i) = 0
1565  zs(j,i) = zs(j,i) + h_init
1566  end if
1567 
1568  xi(i) = xi0 + real(i,dp)*dxi
1569  eta(j) = eta0 + real(j,dp)*deta
1570 
1571  zm(j,i) = zb(j,i)
1572  n_cts(j,i) = -1
1573  kc_cts(j,i) = 0
1574 
1575  h_c(j,i) = zs(j,i)-zm(j,i)
1576  h_t(j,i) = 0.0_dp
1577 
1578  dzs_dtau(j,i) = 0.0_dp
1579  dzm_dtau(j,i) = 0.0_dp
1580  dzb_dtau(j,i) = 0.0_dp
1581  dzl_dtau(j,i) = 0.0_dp
1582  dh_c_dtau(j,i) = 0.0_dp
1583  dh_t_dtau(j,i) = 0.0_dp
1584 
1585 end do
1586 end do
1587 
1588 maske_old = maske
1589 
1590 !-------- Geographic coordinates, metric tensor,
1591 ! gradients of the topography --------
1592 
1593 do i=0, imax
1594 do j=0, jmax
1595 
1596 #if (GRID==0 || GRID==1) /* Stereographic projection */
1597 
1598  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1599  lambda0, phi0, lambda(j,i), phi(j,i))
1600 
1601 #elif (GRID==2) /* Geographic coordinates */
1602 
1603  lambda(j,i) = xi(i)
1604  phi(j,i) = eta(j)
1605 
1606 #endif
1607 
1608 end do
1609 end do
1610 
1611 call metric()
1612 
1613 #if (TOPOGRAD==0)
1614 call topograd_1(dxi, deta, 1)
1615 #elif (TOPOGRAD==1)
1616 call topograd_2(dxi, deta, 1)
1617 #endif
1618 
1619 !-------- Corresponding area of grid points --------
1620 
1621 do i=0, imax
1622 do j=0, jmax
1623  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1624 end do
1625 end do
1626 
1627 end subroutine topography1
1628 
1629 !-------------------------------------------------------------------------------
1630 !> Definition of the initial surface and bedrock topography
1631 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1632 !! For ice-free initial topography with relaxed lithosphere.
1633 !<------------------------------------------------------------------------------
1634 subroutine topography2(dxi, deta)
1635 
1636 #if (GRID==0 || GRID==1)
1638 #endif
1639 
1640  use metric_m
1641  use topograd_m
1642 
1643 implicit none
1644 
1645 real(dp), intent(out) :: dxi, deta
1646 
1647 integer(i4b) :: i, j, n
1648 integer(i4b) :: ios
1649 real(dp) :: xi0, eta0
1650 character :: ch_dummy
1651 
1652 character(len= 8) :: ch_imax
1653 character(len=128) :: fmt4
1654 character(len=256) :: filename_with_path
1655 
1656 write(ch_imax, fmt='(i8)') imax
1657 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
1658 
1659 !-------- Set topography --------
1660 
1661 zl0 = 0.0_dp
1662 
1663 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1664  trim(mask_present_file)
1665 
1666 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
1667 
1668 if (ios /= 0) stop ' >>> topography2: Error when opening the mask file!'
1669 
1670 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1671 
1672 do j=jmax, 0, -1
1673  read(24, fmt=trim(fmt4)) (maske(j,i), i=0,imax)
1674 end do
1675 
1676 close(24, status='keep')
1677 
1678 !-------- Further stuff --------
1679 
1680 #if (GRID==0 || GRID==1)
1681 
1682 dxi = dx *1000.0_dp ! km -> m
1683 deta = dx *1000.0_dp ! km -> m
1684 
1685 xi0 = x0 *1000.0_dp ! km -> m
1686 eta0 = y0 *1000.0_dp ! km -> m
1687 
1688 #elif (GRID==2)
1689 
1690 stop ' >>> topography2: GRID==2 not allowed for this application!'
1691 
1692 #endif
1693 
1694 do i=0, imax
1695 do j=0, jmax
1696 
1697  zs(j,i) = zl0(j,i)
1698  zb(j,i) = zl0(j,i)
1699  zl(j,i) = zl0(j,i)
1700 
1701  xi(i) = xi0 + real(i,dp)*dxi
1702  eta(j) = eta0 + real(j,dp)*deta
1703 
1704  zm(j,i) = zb(j,i)
1705  n_cts(j,i) = -1
1706  kc_cts(j,i) = 0
1707 
1708  h_c(j,i) = 0.0_dp
1709  h_t(j,i) = 0.0_dp
1710 
1711  dzs_dtau(j,i) = 0.0_dp
1712  dzm_dtau(j,i) = 0.0_dp
1713  dzb_dtau(j,i) = 0.0_dp
1714  dzl_dtau(j,i) = 0.0_dp
1715  dh_c_dtau(j,i) = 0.0_dp
1716  dh_t_dtau(j,i) = 0.0_dp
1717 
1718 end do
1719 end do
1720 
1721 maske_old = maske
1722 
1723 !-------- Geographic coordinates, metric tensor,
1724 ! gradients of the topography --------
1725 
1726 do i=0, imax
1727 do j=0, jmax
1728 
1729 #if (GRID==0 || GRID==1) /* Stereographic projection */
1730 
1731  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1732  lambda0, phi0, lambda(j,i), phi(j,i))
1733 
1734 #elif (GRID==2) /* Geographic coordinates */
1735 
1736  lambda(j,i) = xi(i)
1737  phi(j,i) = eta(j)
1738 
1739 #endif
1740 
1741 end do
1742 end do
1743 
1744 call metric()
1745 
1746 #if (TOPOGRAD==0)
1747 call topograd_1(dxi, deta, 1)
1748 #elif (TOPOGRAD==1)
1749 call topograd_2(dxi, deta, 1)
1750 #endif
1751 
1752 !-------- Corresponding area of grid points --------
1753 
1754 do i=0, imax
1755 do j=0, jmax
1756  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1757 end do
1758 end do
1759 
1760 end subroutine topography2
1761 
1762 !-------------------------------------------------------------------------------
1763 !> Definition of the initial surface and bedrock topography
1764 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1765 !! For initial topography from previous simulation.
1766 !<------------------------------------------------------------------------------
1767 subroutine topography3(dxi, deta, z_sl, anfdatname)
1768 
1769  use read_m, only : read_erg_nc
1770 
1771 #if (GRID==0 || GRID==1)
1773 #endif
1774 
1775  use metric_m
1776  use topograd_m
1777 
1778 implicit none
1779 
1780 character(len=100), intent(in) :: anfdatname
1781 
1782 real(dp), intent(out) :: dxi, deta, z_sl
1783 
1784 integer(i4b) :: i, j
1785 
1786 !-------- Read data from time-slice file of previous simulation --------
1787 
1788 call read_erg_nc(z_sl, anfdatname)
1789 
1790 !-------- Set topography of the relaxed bedrock --------
1791 
1792 zl0 = 0.0_dp
1793 
1794 !-------- Further stuff --------
1795 
1796 #if (GRID==0 || GRID==1)
1797 
1798 dxi = dx *1000.0_dp ! km -> m
1799 deta = dx *1000.0_dp ! km -> m
1800 
1801 #elif (GRID==2)
1802 
1803 stop ' >>> topography3: GRID==2 not allowed for this application!'
1804 
1805 #endif
1806 
1807 !-------- Geographic coordinates, metric tensor,
1808 ! gradients of the topography --------
1809 
1810 do i=0, imax
1811 do j=0, jmax
1812 
1813 #if (GRID==0 || GRID==1) /* Stereographic projection */
1814 
1815  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1816  lambda0, phi0, lambda(j,i), phi(j,i))
1817 
1818 #elif (GRID==2) /* Geographic coordinates */
1819 
1820  lambda(j,i) = xi(i)
1821  phi(j,i) = eta(j)
1822 
1823 #endif
1824 
1825 end do
1826 end do
1827 
1828 call metric()
1829 
1830 #if (TOPOGRAD==0)
1831 call topograd_1(dxi, deta, 1)
1832 #elif (TOPOGRAD==1)
1833 call topograd_2(dxi, deta, 1)
1834 #endif
1835 
1836 !-------- Corresponding area of grid points --------
1837 
1838 do i=0, imax
1839 do j=0, jmax
1840  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1841 end do
1842 end do
1843 
1844 end subroutine topography3
1845 
1846 !-------------------------------------------------------------------------------
1847 
1848 end module sico_init_m
1849 !
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: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
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: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) 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), 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), 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: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) 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...