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