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, 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 smb_corr = 0.0_dp
1118 
1119 q_bm = 0.0_dp
1120 q_tld = 0.0_dp
1121 q_b_tot = 0.0_dp
1122 
1123 p_b_w = 0.0_dp
1124 h_w = 0.0_dp
1125 
1126 #if (!defined(TEMP_INIT) || TEMP_INIT==1)
1128 #elif (TEMP_INIT==2)
1130 #elif (TEMP_INIT==3)
1132 #elif (TEMP_INIT==4)
1134 #else
1135  write(6, fmt='(a)') ' >>> sico_init:'
1136  write(6, fmt='(a)') ' TEMP_INIT must be set to either 1, 2, 3 or 4!'
1137  stop
1138 #endif
1139 
1140 #if (ENHMOD==1)
1141  call calc_enhance_1()
1142 #elif (ENHMOD==2)
1143  call calc_enhance_2()
1144 #elif (ENHMOD==3)
1145  call calc_enhance_3(time_init)
1146 #elif (ENHMOD==4)
1147  call calc_enhance_4()
1148 #elif (ENHMOD==5)
1149  call calc_enhance_5()
1150 #else
1151  stop ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
1152 #endif
1153 
1154 ! ------ Ice-free, relaxed bedrock
1155 
1156 #elif (ANF_DAT==2)
1157 
1158 call topography2(dxi, deta)
1159 
1160 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1161 
1162 call boundary(time_init, dtime, dxi, deta, &
1163  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1164 
1165 as_perp_apl = 0.0_dp
1166 
1167 smb_corr = 0.0_dp
1168 
1169 q_bm = 0.0_dp
1170 q_tld = 0.0_dp
1171 q_b_tot = 0.0_dp
1172 
1173 p_b_w = 0.0_dp
1174 h_w = 0.0_dp
1175 
1176 call init_temp_water_age_2()
1177 
1178 #if (ENHMOD==1)
1179  call calc_enhance_1()
1180 #elif (ENHMOD==2)
1181  call calc_enhance_2()
1182 #elif (ENHMOD==3)
1183  call calc_enhance_3(time_init)
1184 #elif (ENHMOD==4)
1185  call calc_enhance_4()
1186 #elif (ENHMOD==5)
1187  call calc_enhance_5()
1188 #else
1189  stop ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
1190 #endif
1191 
1192 ! ------ Read initial state from output of previous simulation
1193 
1194 #elif (ANF_DAT==3)
1195 
1196 call topography3(dxi, deta, z_sl, anfdatname)
1197 
1198 call boundary(time_init, dtime, dxi, deta, &
1199  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1200 
1201 where ((maske==0_i2b).or.(maske==3_i2b))
1202  ! grounded or floating ice
1204 elsewhere ! maske==1_i2b or 2_i2b, ice-free land or sea
1205  as_perp_apl = 0.0_dp
1206 end where
1207 
1208 smb_corr = 0.0_dp
1209 
1210 q_b_tot = q_bm + q_tld
1211 
1212 #endif
1213 
1214 !-------- Inner-point flag --------
1215 
1216 flag_inner_point = .true.
1217 
1218 flag_inner_point(0,:) = .false.
1219 flag_inner_point(jmax,:) = .false.
1220 
1221 flag_inner_point(:,0) = .false.
1222 flag_inner_point(:,imax) = .false.
1223 
1224 !-------- Grounding line and calving front flags --------
1225 
1226 flag_grounding_line_1 = .false.
1227 flag_grounding_line_2 = .false.
1228 
1229 flag_calving_front_1 = .false.
1230 flag_calving_front_2 = .false.
1231 
1232 #if (MARGIN==1 || MARGIN==2)
1233 
1234 continue
1235 
1236 #elif (MARGIN==3)
1237 
1238 do i=1, imax-1
1239 do j=1, jmax-1
1240 
1241  if ( (maske(j,i)==0_i2b) & ! grounded ice
1242  .and. &
1243  ( (maske(j,i+1)==3_i2b) & ! with
1244  .or.(maske(j,i-1)==3_i2b) & ! one
1245  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1246  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1247  ) &
1248  flag_grounding_line_1(j,i) = .true.
1249 
1250  if ( (maske(j,i)==3_i2b) & ! floating ice
1251  .and. &
1252  ( (maske(j,i+1)==0_i2b) & ! with
1253  .or.(maske(j,i-1)==0_i2b) & ! one
1254  .or.(maske(j+1,i)==0_i2b) & ! neighbouring
1255  .or.(maske(j-1,i)==0_i2b) ) & ! grounded ice point
1256  ) &
1257  flag_grounding_line_2(j,i) = .true.
1258 
1259  if ( (maske(j,i)==3_i2b) & ! floating ice
1260  .and. &
1261  ( (maske(j,i+1)==2_i2b) & ! with
1262  .or.(maske(j,i-1)==2_i2b) & ! one
1263  .or.(maske(j+1,i)==2_i2b) & ! neighbouring
1264  .or.(maske(j-1,i)==2_i2b) ) & ! ocean point
1265  ) &
1266  flag_calving_front_1(j,i) = .true.
1267 
1268  if ( (maske(j,i)==2_i2b) & ! ocean
1269  .and. &
1270  ( (maske(j,i+1)==3_i2b) & ! with
1271  .or.(maske(j,i-1)==3_i2b) & ! one
1272  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1273  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1274  ) &
1275  flag_calving_front_2(j,i) = .true.
1276 
1277 end do
1278 end do
1279 
1280 do i=1, imax-1
1281 
1282  j=0
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  j=jmax
1289  if ( (maske(j,i)==2_i2b) & ! ocean
1290  .and. (maske(j-1,i)==3_i2b) & ! with one neighbouring
1291  ) & ! floating ice point
1292  flag_calving_front_2(j,i) = .true.
1293 
1294 end do
1295 
1296 do j=1, jmax-1
1297 
1298  i=0
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  i=imax
1305  if ( (maske(j,i)==2_i2b) & ! ocean
1306  .and. (maske(j,i-1)==3_i2b) & ! with one neighbouring
1307  ) & ! floating ice point
1308  flag_calving_front_2(j,i) = .true.
1309 
1310 end do
1311 
1312 #else
1313 stop ' >>> sico_init: MARGIN must be either 1, 2 or 3!'
1314 #endif
1315 
1316 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
1317 
1318 #if (GRID==0 || GRID==1)
1319 
1320 do ir=-imax, imax
1321 do jr=-jmax, jmax
1322  dist_dxdy(jr,ir) = sqrt( (real(ir,dp)*dxi)**2 + (real(jr,dp)*deta)**2 )
1323  ! distortion due to stereographic projection not accounted for
1324 end do
1325 end do
1326 
1327 #endif
1328 
1329 !-------- Initial velocities --------
1330 
1331 call calc_temp_melt()
1332 
1333 #if (DYNAMICS==1 || DYNAMICS==2)
1334 
1335 call calc_vxy_b_sia(time, z_sl)
1336 call calc_vxy_sia(dzeta_c, dzeta_t)
1337 
1338 #if (MARGIN==3 || DYNAMICS==2)
1339 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
1340 #endif
1341 
1342 call calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
1343 
1344 #if (MARGIN==3)
1345 call calc_vz_floating(z_sl, dxi, deta, dzeta_c)
1346 #endif
1347 
1348 #elif (DYNAMICS==0)
1349 
1350 call calc_vxy_static()
1351 call calc_vz_static()
1352 
1353 #else
1354  stop ' >>> sico_init: DYNAMICS must be either 0, 1 or 2!'
1355 #endif
1356 
1357 call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
1358 
1359 !-------- Initial enthalpies --------
1360 
1361 #if (CALCMOD==0 || CALCMOD==-1)
1362 
1363 do i=0, imax
1364 do j=0, jmax
1365 
1366  do kc=0, kcmax
1367  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1368  end do
1369 
1370  do kt=0, ktmax
1371  enth_t(kt,j,i) = enth_c(0,j,i)
1372  end do
1373 
1374 end do
1375 end do
1376 
1377 #elif (CALCMOD==1)
1378 
1379 do i=0, imax
1380 do j=0, jmax
1381 
1382  do kc=0, kcmax
1383  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1384  end do
1385 
1386  if ( (maske(j,i) == 0_i2b).and.(n_cts(j,i) == 1_i2b) ) then
1387  do kt=0, ktmax
1388  enth_t(kt,j,i) = enth_fct_temp_omega(temp_t_m(kt,j,i), omega_t(kt,j,i))
1389  end do
1390  else
1391  do kt=0, ktmax
1392  enth_t(kt,j,i) = enth_c(0,j,i)
1393  end do
1394  end if
1395 
1396 end do
1397 end do
1398 
1399 #elif (CALCMOD==2 || CALCMOD==3)
1400 
1401 do i=0, imax
1402 do j=0, jmax
1403 
1404  do kc=0, kcmax
1405  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), omega_c(kc,j,i))
1406  end do
1407 
1408  do kt=0, ktmax
1409  enth_t(kt,j,i) = enth_c(0,j,i)
1410  end do
1411 
1412 end do
1413 end do
1414 
1415 #else
1416 
1417 stop ' >>> sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!'
1418 
1419 #endif
1420 
1421 !-------- Initialize time-series files --------
1422 
1423 ! ------ Time-series file for the ice sheet on the whole
1424 
1425 filename_with_path = trim(outpath)//'/'//trim(runname)//'.ser'
1426 
1427 open(12, iostat=ios, file=trim(filename_with_path), status='new')
1428 
1429 if (ios /= 0) stop ' >>> sico_init: Error when opening the ser file!'
1430 
1431 write(12,1102)
1432 write(12,1103)
1433 
1434 1102 format(' t(a) D_Ts(C) z_sl(m)',/, &
1435  ' V(m^3) V_g(m^3) V_f(m^3)', &
1436  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1437  ' V_sle(m) V_t(m^3)', &
1438  ' A_t(m^2)',/, &
1439  ' H_max(m) H_t_max(m)', &
1440  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
1441 1103 format('----------------------------------------------------', &
1442  '---------------------------------------')
1443 
1444 ! ------ Time-series file for deep boreholes
1445 
1446 n_core = 0 ! No boreholes defined
1447 
1448 filename_with_path = trim(outpath)//'/'//trim(runname)//'.core'
1449 
1450 open(14, iostat=ios, file=trim(filename_with_path), status='new')
1451 
1452 write(14,'(1x,a)') '---------------------'
1453 write(14,'(1x,a)') 'No boreholes defined.'
1454 write(14,'(1x,a)') '---------------------'
1455 
1456 !-------- Output of the initial state --------
1457 
1458 #if (defined(OUTPUT_INIT))
1459 
1460 #if (OUTPUT_INIT==0)
1461  flag_init_output = .false.
1462 #elif (OUTPUT_INIT==1)
1463  flag_init_output = .true.
1464 #else
1465  stop ' >>> sico_init: OUTPUT_INIT must be either 0 or 1!'
1466 #endif
1467 
1468 #else
1469  flag_init_output = .true. ! default for undefined parameter OUTPUT_INIT
1470 #endif
1471 
1472 #if (OUTPUT==1)
1473 
1474 #if (ERGDAT==0)
1475  flag_3d_output = .false.
1476 #elif (ERGDAT==1)
1477  flag_3d_output = .true.
1478 #else
1479  stop ' >>> sico_init: ERGDAT must be either 0 or 1!'
1480 #endif
1481 
1482  if (flag_init_output) &
1483  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1484  flag_3d_output, ndat2d, ndat3d)
1485 
1486 #elif (OUTPUT==2)
1487 
1488 if (time_output(1) <= time_init+eps) then
1489 
1490 #if (ERGDAT==0)
1491  flag_3d_output = .false.
1492 #elif (ERGDAT==1)
1493  flag_3d_output = .true.
1494 #else
1495  stop ' >>> sico_init: ERGDAT must be either 0 or 1!'
1496 #endif
1497 
1498  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1499  flag_3d_output, ndat2d, ndat3d)
1500 
1501 end if
1502 
1503 #elif (OUTPUT==3)
1504 
1505  flag_3d_output = .false.
1506 
1507  if (flag_init_output) &
1508  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1509  flag_3d_output, ndat2d, ndat3d)
1510 
1511 if (time_output(1) <= time_init+eps) then
1512 
1513  flag_3d_output = .true.
1514 
1515  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1516  flag_3d_output, ndat2d, ndat3d)
1517 
1518 end if
1519 
1520 #else
1521  stop ' >>> sico_init: OUTPUT must be either 1, 2 or 3!'
1522 #endif
1523 
1524 if (flag_init_output) then
1525  call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1526  call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1527 end if
1528 
1529 end subroutine sico_init
1530 
1531 !-------------------------------------------------------------------------------
1532 !> Definition of the initial surface and bedrock topography
1533 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1534 !! For present-day initial topography.
1535 !<------------------------------------------------------------------------------
1536 subroutine topography1(dxi, deta)
1537 
1538 #if (GRID==0 || GRID==1)
1540 #endif
1541 
1542  use metric_m
1543  use topograd_m
1544 
1545 implicit none
1546 
1547 real(dp), intent(out) :: dxi, deta
1548 
1549 integer(i4b) :: i, j, n
1550 integer(i4b) :: ios, n_dummy
1551 real(dp) :: d_dummy
1552 real(dp) :: xi0, eta0
1553 character :: ch_dummy
1554 
1555 character(len= 8) :: ch_imax
1556 character(len=128) :: fmt4
1557 character(len=256) :: filename_with_path
1558 
1559 write(ch_imax, fmt='(i8)') imax
1560 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
1561 
1562 !-------- Read topography --------
1563 
1564 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1565  trim(zs_present_file)
1566 
1567 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1568 
1569 if (ios /= 0) stop ' >>> topography1: Error when opening the zs file!'
1570 
1571 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1572 
1573 do j=jmax, 0, -1
1574  read(21, fmt=*) (zs(j,i), i=0,imax)
1575 end do
1576 
1577 close(21, status='keep')
1578 
1579 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1580  trim(zl_present_file)
1581 
1582 open(22, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1583 
1584 if (ios /= 0) stop ' >>> topography1: Error when opening the zl file!'
1585 
1586 do n=1, 6; read(22, fmt='(a)') ch_dummy; end do
1587 
1588 do j=jmax, 0, -1
1589  read(22, fmt=*) (zl(j,i), i=0,imax)
1590 end do
1591 
1592 close(22, status='keep')
1593 
1594 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1595  trim(zl0_file)
1596 
1597 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1598 
1599 if (ios /= 0) stop ' >>> topography1: Error when opening the zl0 file!'
1600 
1601 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
1602 
1603 do j=jmax, 0, -1
1604  read(23, fmt=*) (zl0(j,i), i=0,imax)
1605 end do
1606 
1607 close(23, status='keep')
1608 
1609 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1610  trim(mask_present_file)
1611 
1612 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
1613 
1614 if (ios /= 0) stop ' >>> topography1: Error when opening the mask file!'
1615 
1616 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1617 
1618 do j=jmax, 0, -1
1619  read(24, fmt=trim(fmt4)) (maske(j,i), i=0,imax)
1620 end do
1621 
1622 close(24, status='keep')
1623 
1624 !-------- Further stuff --------
1625 
1626 dxi = dx *1000.0_dp ! km -> m
1627 deta = dx *1000.0_dp ! km -> m
1628 
1629 xi0 = x0 *1000.0_dp ! km -> m
1630 eta0 = y0 *1000.0_dp ! km -> m
1631 
1632 do i=0, imax
1633 do j=0, jmax
1634 
1635  if (maske(j,i) <= 1) then
1636  zb(j,i) = zl(j,i)
1637  else ! (maske(j,i>=2)
1638  stop ' >>> topography1: maske(j,i)>=2 not allowed for initial topography!'
1639  end if
1640 
1641  zs(j,i) = zs(j,i) *1000.0_dp
1642  zb(j,i) = zb(j,i) *1000.0_dp ! km --> m
1643  zl(j,i) = zl(j,i) *1000.0_dp
1644  zl0(j,i) = zl0(j,i)*1000.0_dp
1645 
1646  xi(i) = xi0 + real(i,dp)*dxi
1647  eta(j) = eta0 + real(j,dp)*deta
1648 
1649  zm(j,i) = zb(j,i)
1650  n_cts(j,i) = -1
1651  kc_cts(j,i) = 0
1652 
1653  h_c(j,i) = zs(j,i)-zm(j,i)
1654  h_t(j,i) = 0.0_dp
1655 
1656  dzs_dtau(j,i) = 0.0_dp
1657  dzm_dtau(j,i) = 0.0_dp
1658  dzb_dtau(j,i) = 0.0_dp
1659  dzl_dtau(j,i) = 0.0_dp
1660  dh_c_dtau(j,i) = 0.0_dp
1661  dh_t_dtau(j,i) = 0.0_dp
1662 
1663 end do
1664 end do
1665 
1666 maske_old = maske
1667 
1668 !-------- Geographic coordinates, metric tensor,
1669 ! gradients of the topography --------
1670 
1671 do i=0, imax
1672 do j=0, jmax
1673 
1674 #if (GRID==0 || GRID==1) /* Stereographic projection */
1675 
1676  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1677  lambda0, phi0, lambda(j,i), phi(j,i))
1678 
1679 #elif (GRID==2) /* Geographic coordinates */
1680 
1681  lambda(j,i) = xi(i)
1682  phi(j,i) = eta(j)
1683 
1684 #endif
1685 
1686 end do
1687 end do
1688 
1689 call metric()
1690 
1691 #if (TOPOGRAD==0)
1692 call topograd_1(dxi, deta, 1)
1693 #elif (TOPOGRAD==1)
1694 call topograd_2(dxi, deta, 1)
1695 #endif
1696 
1697 !-------- Corresponding area of grid points --------
1698 
1699 do i=0, imax
1700 do j=0, jmax
1701  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1702 end do
1703 end do
1704 
1705 end subroutine topography1
1706 
1707 !-------------------------------------------------------------------------------
1708 !> Definition of the initial surface and bedrock topography
1709 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1710 !! For ice-free initial topography with relaxed lithosphere.
1711 !<------------------------------------------------------------------------------
1712 subroutine topography2(dxi, deta)
1713 
1714 #if (GRID==0 || GRID==1)
1716 #endif
1717 
1718  use metric_m
1719  use topograd_m
1720 
1721 implicit none
1722 
1723 real(dp), intent(out) :: dxi, deta
1724 
1725 integer(i4b) :: i, j, n
1726 integer(i4b) :: ios
1727 real(dp) :: xi0, eta0
1728 character :: ch_dummy
1729 
1730 character(len= 8) :: ch_imax
1731 character(len=128) :: fmt4
1732 character(len=256) :: filename_with_path
1733 
1734 write(ch_imax, fmt='(i8)') imax
1735 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
1736 
1737 !-------- Read topography --------
1738 
1739 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1740  trim(zl0_file)
1741 
1742 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1743 
1744 if (ios /= 0) stop ' >>> topography2: Error when opening the zl0 file!'
1745 
1746 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
1747 
1748 do j=jmax, 0, -1
1749  read(23, fmt=*) (zl0(j,i), i=0,imax)
1750 end do
1751 
1752 close(23, status='keep')
1753 
1754 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1755  trim(mask_present_file)
1756 
1757 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
1758 
1759 if (ios /= 0) stop ' >>> topography2: Error when opening the mask file!'
1760 
1761 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1762 
1763 do j=jmax, 0, -1
1764  read(24, fmt=trim(fmt4)) (maske(j,i), i=0,imax)
1765 end do
1766 
1767 close(24, status='keep')
1768 
1769 !-------- Further stuff --------
1770 
1771 dxi = dx *1000.0_dp ! km -> m
1772 deta = dx *1000.0_dp ! km -> m
1773 
1774 xi0 = x0 *1000.0_dp ! km -> m
1775 eta0 = y0 *1000.0_dp ! km -> m
1776 
1777 do i=0, imax
1778 do j=0, jmax
1779 
1780  if (maske(j,i) <= 1) then
1781  maske(j,i) = 1
1782  zs(j,i) = zl0(j,i)
1783  zb(j,i) = zl0(j,i)
1784  zl(j,i) = zl0(j,i)
1785  else ! (maske(j,i)>=2)
1786  stop ' >>> topography2: maske(j,i)>=2 not allowed for initial topography!'
1787  end if
1788 
1789  zs(j,i) = zs(j,i) *1000.0_dp
1790  zb(j,i) = zb(j,i) *1000.0_dp ! km --> m
1791  zl(j,i) = zl(j,i) *1000.0_dp
1792  zl0(j,i) = zl0(j,i)*1000.0_dp
1793 
1794  xi(i) = xi0 + real(i,dp)*dxi
1795  eta(j) = eta0 + real(j,dp)*deta
1796 
1797  zm(j,i) = zb(j,i)
1798  n_cts(j,i) = -1
1799  kc_cts(j,i) = 0
1800 
1801  h_c(j,i) = 0.0_dp
1802  h_t(j,i) = 0.0_dp
1803 
1804  dzs_dtau(j,i) = 0.0_dp
1805  dzm_dtau(j,i) = 0.0_dp
1806  dzb_dtau(j,i) = 0.0_dp
1807  dzl_dtau(j,i) = 0.0_dp
1808  dh_c_dtau(j,i) = 0.0_dp
1809  dh_t_dtau(j,i) = 0.0_dp
1810 
1811 end do
1812 end do
1813 
1814 maske_old = maske
1815 
1816 !-------- Geographic coordinates, metric tensor,
1817 ! gradients of the topography --------
1818 
1819 do i=0, imax
1820 do j=0, jmax
1821 
1822 #if (GRID==0 || GRID==1) /* Stereographic projection */
1823 
1824  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1825  lambda0, phi0, lambda(j,i), phi(j,i))
1826 
1827 #elif (GRID==2) /* Geographic coordinates */
1828 
1829  lambda(j,i) = xi(i)
1830  phi(j,i) = eta(j)
1831 
1832 #endif
1833 
1834 end do
1835 end do
1836 
1837 call metric()
1838 
1839 #if (TOPOGRAD==0)
1840 call topograd_1(dxi, deta, 1)
1841 #elif (TOPOGRAD==1)
1842 call topograd_2(dxi, deta, 1)
1843 #endif
1844 
1845 !-------- Corresponding area of grid points --------
1846 
1847 do i=0, imax
1848 do j=0, jmax
1849  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1850 end do
1851 end do
1852 
1853 end subroutine topography2
1854 
1855 !-------------------------------------------------------------------------------
1856 !> Definition of the initial surface and bedrock topography
1857 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1858 !! For initial topography from previous simulation.
1859 !<------------------------------------------------------------------------------
1860 subroutine topography3(dxi, deta, z_sl, anfdatname)
1861 
1862  use read_m, only : read_erg_nc
1863 
1864 #if (GRID==0 || GRID==1)
1866 #endif
1867 
1868  use metric_m
1869  use topograd_m
1870 
1871 implicit none
1872 
1873 character(len=100), intent(in) :: anfdatname
1874 
1875 real(dp), intent(out) :: dxi, deta, z_sl
1876 
1877 integer(i4b) :: i, j, n
1878 integer(i4b) :: ios
1879 character(len=256) :: filename_with_path
1880 character :: ch_dummy
1881 
1882 !-------- Read data from time-slice file of previous simulation --------
1883 
1884 call read_erg_nc(z_sl, anfdatname)
1885 
1886 !-------- Read topography of the relaxed bedrock --------
1887 
1888 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1889  trim(zl0_file)
1890 
1891 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1892 
1893 if (ios /= 0) stop ' >>> topography3: Error when opening the zl0 file!'
1894 
1895 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
1896 
1897 do j=jmax, 0, -1
1898  read(23, fmt=*) (zl0(j,i), i=0,imax)
1899 end do
1900 
1901 close(23, status='keep')
1902 
1903 !-------- Further stuff --------
1904 
1905 dxi = dx *1000.0_dp ! km -> m
1906 deta = dx *1000.0_dp ! km -> m
1907 
1908 zl0 = zl0 *1000.0_dp ! km -> m
1909 
1910 !-------- Geographic coordinates, metric tensor,
1911 ! gradients of the topography --------
1912 
1913 do i=0, imax
1914 do j=0, jmax
1915 
1916 #if (GRID==0 || GRID==1) /* Stereographic projection */
1917 
1918  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1919  lambda0, phi0, lambda(j,i), phi(j,i))
1920 
1921 #elif (GRID==2) /* Geographic coordinates */
1922 
1923  lambda(j,i) = xi(i)
1924  phi(j,i) = eta(j)
1925 
1926 #endif
1927 
1928 end do
1929 end do
1930 
1931 call metric()
1932 
1933 #if (TOPOGRAD==0)
1934 call topograd_1(dxi, deta, 1)
1935 #elif (TOPOGRAD==1)
1936 call topograd_2(dxi, deta, 1)
1937 #endif
1938 
1939 !-------- Corresponding area of grid points --------
1940 
1941 do i=0, imax
1942 do j=0, jmax
1943  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1944 end do
1945 end do
1946 
1947 end subroutine topography3
1948 
1949 !-------------------------------------------------------------------------------
1950 
1951 end module sico_init_m
1952 !
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:775
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:823
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:4751
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), dimension(0:jmax, 0:imax) smb_corr
smb_corr(j,i): Diagnosed SMB correction
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: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
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...