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