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