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