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