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