SICOPOLIS V5-dev  Revision 1105
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-2017 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 !<------------------------------------------------------------------------------
36 
37  use sico_types_m
39  use sico_vars_m
40 
41  implicit none
42 
43  public
44 
45 contains
46 
47 !-------------------------------------------------------------------------------
48 !> Main routine of sico_init_m: Initialisations for SICOPOLIS.
49 !<------------------------------------------------------------------------------
50 subroutine sico_init(delta_ts, glac_index, &
51  mean_accum, &
52  dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
53  time, time_init, time_end, time_output, &
54  dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
55  z_sl, dzsl_dtau, z_mar, &
56  ii, jj, nn, &
57  ndat2d, ndat3d, n_output, &
58  runname)
59 
60  use compare_float_m
64 
65 #if (NETCDF > 1)
66  use netcdf
67  use nc_check_m
68 #endif
69 
70 #if (GRID==0 || GRID==1)
72 #endif
73 
75 
76  use boundary_m
78  use calc_enhance_m
79  use calc_vxy_m
80  use calc_vz_m
81  use calc_dxyz_m
83 
84  use output_m
85 
86 implicit none
87 
88 integer(i4b), intent(out) :: ii((imax+1)*(jmax+1)), &
89  jj((imax+1)*(jmax+1)), &
90  nn(0:jmax,0:imax)
91 integer(i4b), intent(out) :: ndat2d, ndat3d
92 integer(i4b), intent(out) :: n_output
93 real(dp), intent(out) :: delta_ts, glac_index
94 real(dp), intent(out) :: mean_accum
95 real(dp), intent(out) :: dtime, dtime_temp, dtime_wss, &
96  dtime_out, dtime_ser
97 real(dp), intent(out) :: time, time_init, time_end, time_output(100)
98 real(dp), intent(out) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
99 real(dp), intent(out) :: z_sl, dzsl_dtau, z_mar
100 
101 character(len=100), intent(out) :: runname
102 
103 integer(i4b) :: i, j, kc, kt, kr, m, n, ir, jr
104 integer(i4b) :: ios, ios1, ios2, ios3, ios4
105 integer(i4b) :: ierr
106 integer(i2b), dimension(0:JMAX,0:IMAX) :: maske_ref
107 real(dp) :: dtime0, dtime_temp0, dtime_wss0, dtime_out0, dtime_ser0
108 real(dp) :: time_init0, time_end0, time_output0(100)
109 real(dp), dimension(0:JMAX,0:IMAX) :: precip_factor
110 real(dp) :: d_dummy
111 real(dp) :: transition_length_sliding, quasi_time, d_quasi_time
112 real(dp), dimension(0:JMAX,0:IMAX) :: r_mask_sedi_new
113 character(len=100) :: anfdatname, target_topo_dat_name
114 character(len=256) :: filename_with_path
115 character(len=256) :: shell_command
116 character :: ch_dummy
117 logical :: flag_init_output, flag_3d_output
118 
119 #if (defined(INITMIP_SMB_ANOM_FILE))
120 real(sp), dimension(0:IMAX,0:JMAX) :: as_anom_initmip_conv
121 #endif
122 
123 #if (defined(INITMIP_BMB_ANOM_FILE))
124 real(sp), dimension(0:IMAX,0:JMAX) :: ab_anom_initmip_conv
125 #endif
126 
127 #if (NETCDF > 1)
128 integer(i4b) :: dimid, ncid, ncv
129 ! dimid: Dimension ID
130 ! ncid: File ID
131 ! ncv: Variable ID
132 #endif
133 
134 character(len=64), parameter :: thisroutine = 'sico_init'
135 
136 character(len=64), parameter :: fmt1 = '(a)', &
137  fmt2 = '(a,i4)', &
138  fmt2a = '(a,i0)', &
139  fmt3 = '(a,es12.4)'
140 
141 character(len= 8) :: ch_imax
142 character(len=128) :: fmt4
143 
144 write(ch_imax, fmt='(i8)') imax
145 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
146 
147 write(unit=6, fmt='(a)') ' '
148 write(unit=6, fmt='(a)') ' -------- sico_init --------'
149 write(unit=6, fmt='(a)') ' '
150 
151 !-------- Name of the computational domain --------
152 
153 #if (defined(ANT))
154 ch_domain_long = 'Antarctica'
155 ch_domain_short = 'ant'
156 
157 #elif (defined(ASF))
158 ch_domain_long = 'Austfonna'
159 ch_domain_short = 'asf'
160 
161 #elif (defined(EMTP2SGE))
162 ch_domain_long = 'EISMINT Phase 2 Simplified Geometry Experiment'
163 ch_domain_short = 'emtp2sge'
164 
165 #elif (defined(GRL))
166 ch_domain_long = 'Greenland'
167 ch_domain_short = 'grl'
168 
169 #elif (defined(NHEM))
170 ch_domain_long = 'Northern hemisphere'
171 ch_domain_short = 'nhem'
172 
173 #elif (defined(SCAND))
174 ch_domain_long = 'Scandinavia and Eurasia'
175 ch_domain_short = 'scand'
176 
177 #elif (defined(TIBET))
178 ch_domain_long = 'Tibet'
179 ch_domain_short = 'tibet'
180 
181 #elif (defined(NMARS))
182 ch_domain_long = 'North polar cap of Mars'
183 ch_domain_short = 'nmars'
184 
185 #elif (defined(SMARS))
186 ch_domain_long = 'South polar cap of Mars'
187 ch_domain_short = 'smars'
188 
189 #elif (defined(XYZ))
190 ch_domain_long = 'XYZ'
191 ch_domain_short = 'xyz'
192 #if (defined(HEINO))
193 ch_domain_long = trim(ch_domain_long)//'/ISMIP HEINO'
194 #endif
195 
196 #else
197 stop ' >>> sico_init: No valid domain specified!'
198 #endif
199 
200 !-------- Some initial values --------
201 
202 n_output = 0
203 
204 dtime = 0.0_dp
205 dtime_temp = 0.0_dp
206 dtime_wss = 0.0_dp
207 dtime_out = 0.0_dp
208 dtime_ser = 0.0_dp
209 
210 time = 0.0_dp
211 time_init = 0.0_dp
212 time_end = 0.0_dp
213 time_output = 0.0_dp
214 
215 !-------- Initialisation of the Library of Iterative Solvers Lis,
216 ! if required --------
217 
218 #if (CALCTHK==3 || CALCTHK==6 || MARGIN==3 || DYNAMICS==2)
219  call lis_initialize(ierr)
220 #endif
221 
222 !-------- Read physical parameters --------
223 
224 call read_phys_para()
225 
226 call ice_mat_eqs_pars(rf, r_t, kappa, c, -190, 10)
227 
228 ! ------ Some auxiliary quantities required for the enthalpy method
229 
230 call calc_c_int_table(c, -190, 10, l)
232 
233 !-------- Compatibility check of the SICOPOLIS version with the header file
234 
235 if ( trim(version) /= trim(sico_version) ) &
236  stop ' >>> sico_init: SICOPOLIS version not compatible with header file!'
237 
238 !-------- Check whether the dynamics and thermodynamics modes are defined
239 
240 #if (!defined(DYNAMICS))
241 stop ' >>> sico_init: DYNAMICS not defined in the header file!'
242 #endif
243 
244 #if (!defined(CALCMOD))
245 stop ' >>> sico_init: CALCMOD not defined in the header file!'
246 #endif
247 
248 #if (defined(ENTHMOD))
249 write(6, fmt='(a)') ' >>> sico_init: ENTHMOD must not be defined any more.'
250 write(6, fmt='(a)') ' Please update your header file!'
251 stop
252 #endif
253 
254 !-------- Compatibility check of the horizontal resolution with the
255 ! number of grid points --------
256 
257 #if (GRID==0 || GRID==1)
258 
259 if (approx_equal(dx, 64.0_dp, eps_sp_dp)) then
260 
261  if ( (imax /= 95).or.(jmax /= 95) ) then
262  stop ' >>> sico_init: IMAX and/or JMAX wrong!'
263  end if
264 
265 else if (approx_equal(dx, 40.0_dp, eps_sp_dp)) then
266 
267  if ( (imax /= 152).or.(jmax /= 152) ) then
268  stop ' >>> sico_init: IMAX and/or JMAX wrong!'
269  end if
270 
271 else if (approx_equal(dx, 32.0_dp, eps_sp_dp)) then
272 
273  if ( (imax /= 190).or.(jmax /= 190) ) then
274  stop ' >>> sico_init: IMAX and/or JMAX wrong!'
275  end if
276 
277 else if (approx_equal(dx, 20.0_dp, eps_sp_dp)) then
278 
279  if ( (imax /= 304).or.(jmax /= 304) ) then
280  stop ' >>> sico_init: IMAX and/or JMAX wrong!'
281  end if
282 
283 else if (approx_equal(dx, 16.0_dp, eps_sp_dp)) then
284 
285  if ( (imax /= 380).or.(jmax /= 380) ) then
286  stop ' >>> sico_init: IMAX and/or JMAX wrong!'
287  end if
288 
289 else if (approx_equal(dx, 10.0_dp, eps_sp_dp)) then
290 
291  if ( (imax /= 608).or.(jmax /= 608) ) then
292  stop ' >>> sico_init: IMAX and/or JMAX wrong!'
293  end if
294 
295 else if (approx_equal(dx, 8.0_dp, eps_sp_dp)) then
296 
297  if ( (imax /= 760).or.(jmax /= 760) ) then
298  stop ' >>> sico_init: IMAX and/or JMAX wrong!'
299  end if
300 
301 else
302  stop ' >>> sico_init: DX wrong!'
303 end if
304 
305 #elif (GRID==2)
306 
307 stop ' >>> sico_init: GRID==2 not allowed for this application!'
308 
309 #endif
310 
311 !-------- Compatibility check of the thermodynamics mode
312 ! (cold vs. polythermal vs. enthalpy method)
313 ! and the number of grid points in the lower (kt) ice domain --------
314 
315 #if (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1)
316 
317 if (ktmax > 2) then
318  write(6, fmt='(a)') ' >>> sico_init: For options CALCMOD==0, 2, 3 or -1,'
319  write(6, fmt='(a)') ' the separate kt domain is redundant.'
320  write(6, fmt='(a)') ' Therefore, consider setting KTMAX to 2.'
321  write(6, fmt='(a)') ' '
322 end if
323 
324 #endif
325 
326 !-------- Compatibility check of surface-temperature
327 ! and precipitation switches --------
328 
329 #if (TSURFACE == 5 && ACCSURFACE != 5)
330 stop ' >>> sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
331 #endif
332 
333 #if (TSURFACE != 5 && ACCSURFACE == 5)
334 stop ' >>> sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
335 #endif
336 
337 #if (TSURFACE == 6)
338 #if (ACCSURFACE != 6)
339 write(6, fmt='(a)') ' >>> sico_init: Options TSURFACE==6, ACCSURFACE==6'
340 write(6, fmt='(a)') ' and NETCDF==2 must be used together!'
341 stop
342 #elif (NETCDF != 2)
343 write(6, fmt='(a)') ' >>> sico_init: Options TSURFACE==6, ACCSURFACE==6'
344 write(6, fmt='(a)') ' and NETCDF==2 must be used together!'
345 stop
346 #endif
347 #endif
348 
349 #if (ACCSURFACE == 6)
350 #if (TSURFACE != 6)
351 write(6, fmt='(a)') ' >>> sico_init: Options TSURFACE==6, ACCSURFACE==6'
352 write(6, fmt='(a)') ' and NETCDF==2 must be used together!'
353 stop
354 #elif (NETCDF != 2)
355 write(6, fmt='(a)') ' >>> sico_init: Options TSURFACE==6, ACCSURFACE==6'
356 write(6, fmt='(a)') ' and NETCDF==2 must be used together!'
357 stop
358 #endif
359 #endif
360 
361 #if (defined(INITMIP_SMB_ANOM_FILE) && NETCDF != 2)
362 if ( trim(adjustl(initmip_smb_anom_file)) /= 'none' ) then
363  write(6, fmt='(a)') ' >>> sico_init: Defining INITMIP_SMB_ANOM_FILE'
364  write(6, fmt='(a)') ' must be used together with NETCDF==2!'
365  stop
366 end if
367 #endif
368 
369 #if (defined(INITMIP_BMB_ANOM_FILE) && NETCDF != 2)
370 if ( trim(adjustl(initmip_bmb_anom_file)) /= 'none' ) then
371  write(6, fmt='(a)') ' >>> sico_init: Defining INITMIP_BMB_ANOM_FILE'
372  write(6, fmt='(a)') ' must be used together with NETCDF==2!'
373  stop
374 end if
375 #endif
376 
377 !-------- Compatibility check of discretization schemes for the horizontal and
378 ! vertical advection terms in the temperature and age equations --------
379 
380 #if (ADV_HOR==1)
381 stop ' >>> sico_init: Option ADV_HOR==1 (central differences) not defined!'
382 #endif
383 
384 !-------- Check whether for the shallow shelf
385 ! or shelfy stream approximation
386 ! the chosen grid is Cartesian coordinates
387 ! without distortion correction (GRID==0) --------
388 
389 #if ((MARGIN==3 && DYNAMICS==1) || DYNAMICS==2) /* requires SSA or SStA */
390 #if (GRID != 0)
391 write(6, fmt='(a)') ' >>> sico_init: Distortion correction not implemented'
392 write(6, fmt='(a)') ' for the shallow shelf approximation (SSA)'
393 write(6, fmt='(a)') ' or the shelfy stream approximation (SStA)'
394 write(6, fmt='(a)') ' -> GRID==0 required!'
395 stop
396 #endif
397 #endif
398 
399 !-------- Setting of forcing flag --------
400 
401 #if (TSURFACE <= 4)
402 
403 forcing_flag = 1 ! forcing by delta_ts
404 
405 #elif (TSURFACE == 5)
406 
407 forcing_flag = 2 ! forcing by glac_index
408 
409 #elif (TSURFACE == 6)
410 
411 forcing_flag = 3 ! forcing by time-dependent surface temperature
412  ! and precipitation data
413 
414 #endif
415 
416 !-------- Initialization of numerical time steps --------
417 
418 dtime0 = dtime0
419 dtime_temp0 = dtime_temp0
420 #if (REBOUND==2)
421 dtime_wss0 = dtime_wss0
422 #endif
423 
424 !-------- Further initializations --------
425 
426 dzeta_c = 1.0_dp/real(kcmax,dp)
427 dzeta_t = 1.0_dp/real(ktmax,dp)
428 dzeta_r = 1.0_dp/real(krmax,dp)
429 
430 ndat2d = 1
431 ndat3d = 1
432 
433 !-------- General abbreviations --------
434 
435 ! ------ kc domain
436 
437 if (deform >= eps) then
438 
439  flag_aa_nonzero = .true. ! non-equidistant grid
440 
441  aa = deform
442  ea = exp(aa)
443 
444  kc=0
445  zeta_c(kc) = 0.0_dp
446  eaz_c(kc) = 1.0_dp
447  eaz_c_quotient(kc) = 0.0_dp
448 
449  do kc=1, kcmax-1
450  zeta_c(kc) = kc*dzeta_c
451  eaz_c(kc) = exp(aa*zeta_c(kc))
452  eaz_c_quotient(kc) = (eaz_c(kc)-1.0_dp)/(ea-1.0_dp)
453  end do
454 
455  kc=kcmax
456  zeta_c(kc) = 1.0_dp
457  eaz_c(kc) = exp(aa)
458  eaz_c_quotient(kc) = 1.0_dp
459 
460 else
461 
462  flag_aa_nonzero = .false. ! equidistant grid
463 
464  aa = 0.0_dp
465  ea = 1.0_dp
466 
467  kc=0
468  zeta_c(kc) = 0.0_dp
469  eaz_c(kc) = 1.0_dp
470  eaz_c_quotient(kc) = 0.0_dp
471 
472  do kc=1, kcmax-1
473  zeta_c(kc) = kc*dzeta_c
474  eaz_c(kc) = 1.0_dp
475  eaz_c_quotient(kc) = zeta_c(kc)
476  end do
477 
478  kc=kcmax
479  zeta_c(kc) = 1.0_dp
480  eaz_c(kc) = 1.0_dp
481  eaz_c_quotient(kc) = 1.0_dp
482 
483 end if
484 
485 ! ------ kt domain
486 
487 kt=0
488 zeta_t(kt) = 0.0_dp
489 
490 do kt=1, ktmax-1
491  zeta_t(kt) = kt*dzeta_t
492 end do
493 
494 kt=ktmax
495 zeta_t(kt) = 1.0_dp
496 
497 ! ------ kr domain
498 
499 kr=0
500 zeta_r(kr) = 0.0_dp
501 
502 do kr=1, krmax-1
503  zeta_r(kr) = kr*dzeta_r
504 end do
505 
506 kr=krmax
507 zeta_r(kr) = 1.0_dp
508 
509 !-------- Construction of a vector (with index n) from a 2-d array
510 ! (with indices i, j) by diagonal numbering --------
511 
512 n=1
513 
514 do m=0, imax+jmax
515  do i=m, 0, -1
516  j = m-i
517  if ((i <= imax).and.(j <= jmax)) then
518  ii(n) = i
519  jj(n) = j
520  nn(j,i) = n
521  n=n+1
522  end if
523  end do
524 end do
525 
526 !-------- Specification of current simulation --------
527 
528 runname = runname
529 anfdatname = anfdatname
530 
531 #if (defined(YEAR_ZERO))
533 #else
534 year_zero = 2000.0_dp ! default value 2000 CE
535 #endif
536 
537 time_init0 = time_init0
538 time_end0 = time_end0
539 dtime_ser0 = dtime_ser0
540 
541 #if (OUTPUT==1 || OUTPUT==3)
542 dtime_out0 = dtime_out0
543 #endif
544 
545 #if (OUTPUT==2 || OUTPUT==3)
546 n_output = n_output
547 time_output0( 1) = time_out0_01
548 time_output0( 2) = time_out0_02
549 time_output0( 3) = time_out0_03
550 time_output0( 4) = time_out0_04
551 time_output0( 5) = time_out0_05
552 time_output0( 6) = time_out0_06
553 time_output0( 7) = time_out0_07
554 time_output0( 8) = time_out0_08
555 time_output0( 9) = time_out0_09
556 time_output0(10) = time_out0_10
557 time_output0(11) = time_out0_11
558 time_output0(12) = time_out0_12
559 time_output0(13) = time_out0_13
560 time_output0(14) = time_out0_14
561 time_output0(15) = time_out0_15
562 time_output0(16) = time_out0_16
563 time_output0(17) = time_out0_17
564 time_output0(18) = time_out0_18
565 time_output0(19) = time_out0_19
566 time_output0(20) = time_out0_20
567 #endif
568 
569 !-------- Write log file --------
570 
571 shell_command = 'if [ ! -d'
572 shell_command = trim(shell_command)//' '//outpath
573 shell_command = trim(shell_command)//' '//'] ; then mkdir'
574 shell_command = trim(shell_command)//' '//outpath
575 shell_command = trim(shell_command)//' '//'; fi'
576 call system(trim(shell_command))
577  ! Check whether directory OUTPATH exists. If not, it is created.
578 
579 filename_with_path = trim(outpath)//'/'//trim(runname)//'.log'
580 
581 open(10, iostat=ios, file=trim(filename_with_path), status='new')
582 
583 if (ios /= 0) stop ' >>> sico_init: Error when opening the log file!'
584 
585 write(10, fmt=trim(fmt1)) 'Computational domain:'
586 write(10, fmt=trim(fmt1)) trim(ch_domain_long)
587 write(10, fmt=trim(fmt1)) ' '
588 
589 write(10, fmt=trim(fmt2a)) 'GRID = ', grid
590 write(10, fmt=trim(fmt1)) ' '
591 
592 write(10, fmt=trim(fmt2)) 'imax =', imax
593 write(10, fmt=trim(fmt2)) 'jmax =', jmax
594 write(10, fmt=trim(fmt2)) 'kcmax =', kcmax
595 write(10, fmt=trim(fmt2)) 'ktmax =', ktmax
596 write(10, fmt=trim(fmt2)) 'krmax =', krmax
597 write(10, fmt=trim(fmt1)) ' '
598 
599 write(10, fmt=trim(fmt3)) 'a =', aa
600 write(10, fmt=trim(fmt1)) ' '
601 
602 #if (GRID==0 || GRID==1)
603 write(10, fmt=trim(fmt3)) 'x0 =', x0
604 write(10, fmt=trim(fmt3)) 'y0 =', y0
605 write(10, fmt=trim(fmt3)) 'dx =', dx
606 #elif (GRID==2)
607 stop ' >>> sico_init: GRID==2 not allowed for this application!'
608 #endif
609 write(10, fmt=trim(fmt1)) ' '
610 
611 write(10, fmt=trim(fmt3)) 'year_zero =', year_zero
612 write(10, fmt=trim(fmt3)) 'time_init =', time_init0
613 write(10, fmt=trim(fmt3)) 'time_end =', time_end0
614 write(10, fmt=trim(fmt3)) 'dtime =', dtime0
615 write(10, fmt=trim(fmt3)) 'dtime_temp =', dtime_temp0
616 #if (REBOUND==2)
617 write(10, fmt=trim(fmt3)) 'dtime_wss =', dtime_wss0
618 #endif
619 write(10, fmt=trim(fmt1)) ' '
620 
621 write(10, fmt=trim(fmt2)) 'ANF_DAT =', anf_dat
622 write(10, fmt=trim(fmt1)) 'zs_present file = '//zs_present_file
623 #if (ANF_DAT==1)
624 #if (defined(ZB_PRESENT_FILE))
625 write(10, fmt=trim(fmt1)) 'zb_present file = '//zb_present_file
626 #endif
627 write(10, fmt=trim(fmt1)) 'zl_present file = '//zl_present_file
628 #endif
629 write(10, fmt=trim(fmt1)) 'zl0 file = '//zl0_file
630 write(10, fmt=trim(fmt1)) 'mask_present file = '//mask_present_file
631 #if (ANF_DAT==1 && defined(TEMP_INIT))
632 write(10, fmt=trim(fmt2)) 'TEMP_INIT =', temp_init
633 #endif
634 #if (ANF_DAT==3)
635 write(10, fmt=trim(fmt1)) 'Initial-value file = '//anfdatname
636 write(10, fmt=trim(fmt1)) 'Path to initial-value file = '//anfdatpath
637 #endif
638 write(10, fmt=trim(fmt1)) ' '
639 
640 #if (THK_EVOL==2)
641 write(10, fmt=trim(fmt3)) 'time_target_topo_init =', time_target_topo_init0
642 write(10, fmt=trim(fmt3)) 'time_target_topo_final =', time_target_topo_final0
643 write(10, fmt=trim(fmt1)) 'Target-topography file = '//target_topo_dat_name
644 write(10, fmt=trim(fmt1)) 'Path to target-topography file = '//target_topo_dat_path
645 write(10, fmt=trim(fmt1)) ' '
646 #endif
647 
648 #if (THK_EVOL==3)
649 write(10, fmt=trim(fmt3)) 'target_topo_tau =', target_topo_tau0
650 write(10, fmt=trim(fmt1)) 'Target-topography file = '//target_topo_dat_name
651 write(10, fmt=trim(fmt1)) 'Path to target-topography file = '//target_topo_dat_path
652 write(10, fmt=trim(fmt1)) ' '
653 #endif
654 
655 #if (THK_EVOL==4)
656 write(10, fmt=trim(fmt1)) 'Maximum ice extent mask file = '//mask_maxextent_file
657 write(10, fmt=trim(fmt1)) ' '
658 #endif
659 
660 write(10, fmt=trim(fmt1)) 'Physical-parameter file = '//phys_para_file
661 write(10, fmt=trim(fmt1)) ' '
662 
663 #if (CALCTHK==2 || CALCTHK==3 || CALCTHK==5 || CALCTHK==6)
664 write(10, fmt=trim(fmt3)) 'ovi_weight =', ovi_weight
665 #if (CALCTHK==2 || CALCTHK==5)
666 write(10, fmt=trim(fmt3)) 'omega_sor =', omega_sor
667 #if (ITER_MAX_SOR>0)
668 write(10, fmt=trim(fmt2a)) 'iter_max_sor = ', iter_max_sor
669 #endif
670 #endif
671 write(10, fmt=trim(fmt1)) ' '
672 #endif
673 
674 #if (TSURFACE==1)
675 write(10, fmt=trim(fmt3)) 'delta_ts0 =', delta_ts0
676 #elif (TSURFACE==3)
677 write(10, fmt=trim(fmt3)) 'sine_amplit =', sine_amplit
678 write(10, fmt=trim(fmt3)) 'sine_period =', sine_period
679 #elif (TSURFACE==4)
680 write(10, fmt=trim(fmt1)) 'GRIP file = '//grip_temp_file
681 write(10, fmt=trim(fmt3)) 'grip_temp_fact =', grip_temp_fact
682 #elif (TSURFACE==5)
683 write(10, fmt=trim(fmt1)) 'Glacial-index file = '//glac_ind_file
684 write(10, fmt=trim(fmt1)) 'temp_ma_anom file = '//temp_ma_anom_file
685 write(10, fmt=trim(fmt3)) 'temp_ma_anom fact = ', temp_ma_anom_fact
686 write(10, fmt=trim(fmt1)) 'temp_mj_anom file = '//temp_mj_anom_file
687 write(10, fmt=trim(fmt3)) 'temp_mj_anom fact = ', temp_mj_anom_fact
688 #elif (TSURFACE==6)
689 write(10, fmt=trim(fmt1)) 'temp_precip_anom file = '//temp_precip_anom_file
690 write(10, fmt=trim(fmt3)) 'temp_ma_anom fact = ', temp_ma_anom_fact
691 write(10, fmt=trim(fmt3)) 'temp_mj_anom fact = ', temp_mj_anom_fact
692 #endif
693 
694 write(10, fmt=trim(fmt1)) 'precip_present_file = '//precip_present_file
695 #if (defined(PRECIP_FACTOR_FILE))
696 if ( trim(adjustl(precip_factor_file)) /= 'none' ) &
697  write(10, fmt=trim(fmt1)) 'precip_factor_file = '//precip_factor_file
698 #endif
699 #if (ACCSURFACE==1)
700 write(10, fmt=trim(fmt3)) 'accfact =', accfact
701 #elif (ACCSURFACE==2 || ACCSURFACE==3)
702 write(10, fmt=trim(fmt3)) 'gamma_s =', gamma_s
703 #elif (ACCSURFACE==5)
704 write(10, fmt=trim(fmt1)) 'precip_anom file = '//precip_anom_file
705 write(10, fmt=trim(fmt3)) 'precip_anom fact = ', precip_anom_fact
706 #elif (ACCSURFACE==6)
707 write(10, fmt=trim(fmt1)) 'temp_precip_anom file = '//temp_precip_anom_file
708 write(10, fmt=trim(fmt3)) 'precip_anom fact = ', precip_anom_fact
709 #endif
710 #if (ACCSURFACE <= 3)
711 write(10, fmt=trim(fmt2)) 'ELEV_DESERT =', elev_desert
712 #if (ELEV_DESERT == 1)
713 write(10, fmt=trim(fmt3)) 'gamma_p =', gamma_p
714 write(10, fmt=trim(fmt3)) 'zs_thresh =', zs_thresh
715 #endif
716 #endif
717 
718 #if (ABLSURFACE==3)
719 write(10, fmt=trim(fmt3)) 'lambda_lti =', lambda_lti
720 write(10, fmt=trim(fmt3)) 'temp_lti =', temp_lti
721 #endif
722 
723 write(10, fmt=trim(fmt1)) ' '
724 
725 #if (defined(INITMIP_CONST_SMB))
726 write(10, fmt=trim(fmt2a)) 'INITMIP_CONST_SMB = ', initmip_const_smb
727 #endif
728 #if (defined(INITMIP_SMB_ANOM_FILE))
729 if ( trim(adjustl(initmip_smb_anom_file)) /= 'none' ) then
730  write(10, fmt=trim(fmt1)) 'initmip_smb_anom file = '//initmip_smb_anom_file
731 end if
732 #endif
733 #if (defined(INITMIP_CONST_SMB) || defined(INITMIP_SMB_ANOM_FILE))
734 write(10, fmt=trim(fmt1)) ' '
735 #endif
736 
737 write(10, fmt=trim(fmt2a)) 'SEA_LEVEL = ', sea_level
738 #if (SEA_LEVEL==1)
739 write(10, fmt=trim(fmt3)) 'z_sl0 =', z_sl0
740 #elif (SEA_LEVEL==3)
741 write(10, fmt=trim(fmt1)) 'sea-level file = '//sea_level_file
742 #endif
743 write(10, fmt=trim(fmt1)) ' '
744 
745 #if (MARGIN==2)
746 #if (MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3)
747 write(10, fmt=trim(fmt3)) 'z_mar =', z_mar
748 write(10, fmt=trim(fmt1)) ' '
749 #elif (MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \
750  || marine_ice_calving==6 || marine_ice_calving==7)
751 write(10, fmt=trim(fmt3)) 'fact_z_mar =', fact_z_mar
752 write(10, fmt=trim(fmt1)) ' '
753 #elif (MARINE_ICE_FORMATION==2 && MARINE_ICE_CALVING==9)
754 write(10, fmt=trim(fmt3)) 'calv_uw_coeff =', calv_uw_coeff
755 write(10, fmt=trim(fmt3)) 'r1_calv_uw =', r1_calv_uw
756 write(10, fmt=trim(fmt3)) 'r2_calv_uw =', r2_calv_uw
757 write(10, fmt=trim(fmt1)) ' '
758 #endif
759 #elif (MARGIN==3)
760 #if (ICE_SHELF_CALVING==2)
761 write(10, fmt=trim(fmt3)) 'H_calv =', h_calv
762 write(10, fmt=trim(fmt1)) ' '
763 #endif
764 #endif
765 
766 #if (defined(BASAL_HYDROLOGY))
767 write(10, fmt=trim(fmt2a)) 'BASAL_HYDROLOGY = ', basal_hydrology
768 #endif
769 write(10, fmt=trim(fmt2a)) 'SLIDE_LAW = ', slide_law
770 write(10, fmt=trim(fmt3)) 'c_slide =', c_slide
771 write(10, fmt=trim(fmt3)) 'gamma_slide =', gamma_slide
772 write(10, fmt=trim(fmt2a)) 'p_weert = ', p_weert
773 write(10, fmt=trim(fmt2a)) 'q_weert = ', q_weert
774 #if (defined(TIME_RAMP_UP_SLIDE))
775 write(10, fmt=trim(fmt3)) 'time_ramp_up_slide =', time_ramp_up_slide
776 #endif
777 #if (SLIDE_LAW==2)
778 write(10, fmt=trim(fmt3)) 'red_pres_limit_fact =', red_pres_limit_fact
779 #endif
780 #if (BASAL_HYDROLOGY==1)
781 write(10, fmt=trim(fmt3)) 'c_Hw_slide =', c_hw_slide
782 write(10, fmt=trim(fmt3)) 'Hw0_slide =', hw0_slide
783 #endif
784 #if (defined(SEDI_SLIDE))
785 write(10, fmt=trim(fmt2a)) 'SEDI_SLIDE = ', sedi_slide
786 #endif
787 #if (defined(SEDI_SLIDE) && SEDI_SLIDE==2)
788 write(10, fmt=trim(fmt1)) 'Sediment-mask file = '//mask_sedi_file
789 write(10, fmt=trim(fmt3)) 'c_slide_sedi =', c_slide_sedi
790 write(10, fmt=trim(fmt3)) 'gamma_slide_sedi =', gamma_slide_sedi
791 write(10, fmt=trim(fmt2a)) 'p_weert_sedi = ', p_weert_sedi
792 write(10, fmt=trim(fmt2a)) 'q_weert_sedi = ', q_weert_sedi
793 #if (defined(TRANS_LENGTH_SL))
794 write(10, fmt=trim(fmt3)) 'trans_length_sl =', trans_length_sl
795 #endif
796 #endif
797 write(10, fmt=trim(fmt1)) ' '
798 
799 write(10, fmt=trim(fmt2a)) 'Q_GEO_MOD = ', q_geo_mod
800 #if (Q_GEO_MOD==1)
801 write(10, fmt=trim(fmt3)) 'q_geo =', q_geo
802 #elif (Q_GEO_MOD==2)
803 write(10, fmt=trim(fmt1)) 'q_geo file = '//q_geo_file
804 #endif
805 write(10, fmt=trim(fmt1)) ' '
806 
807 #if (defined(MARINE_ICE_BASAL_MELTING))
808 write(10, fmt=trim(fmt2)) 'MARINE_ICE_BASAL_MELTING =', marine_ice_basal_melting
809 #if (MARINE_ICE_BASAL_MELTING==2 || MARINE_ICE_BASAL_MELTING==3)
810 write(10, fmt=trim(fmt3)) 'qbm_marine =', qbm_marine
811 #endif
812 write(10, fmt=trim(fmt1)) ' '
813 #endif
814 
815 #if (MARGIN==3)
816 write(10, fmt=trim(fmt2)) 'FLOATING_ICE_BASAL_MELTING =', floating_ice_basal_melting
817 #if (FLOATING_ICE_BASAL_MELTING<=3)
818 write(10, fmt=trim(fmt3)) 'qbm_float_1 =', qbm_float_1
819 write(10, fmt=trim(fmt3)) 'qbm_float_2 =', qbm_float_2
820 #endif
821 write(10, fmt=trim(fmt3)) 'qbm_float_3 =', qbm_float_3
822 write(10, fmt=trim(fmt3)) 'z_abyss =', z_abyss
823 #if (FLOATING_ICE_BASAL_MELTING==4)
824 write(10, fmt=trim(fmt3)) 'temp_ocean =', temp_ocean
825 write(10, fmt=trim(fmt3)) 'Omega_qbm =', omega_qbm
826 write(10, fmt=trim(fmt3)) 'alpha_qbm =', alpha_qbm
827 #endif
828 #if (FLOATING_ICE_BASAL_MELTING==4 || FLOATING_ICE_BASAL_MELTING==5)
829 write(10, fmt=trim(fmt3)) 'H_w_0 =', h_w_0
830 #endif
831 #if (defined(INITMIP_BMB_ANOM_FILE))
832 if ( trim(adjustl(initmip_bmb_anom_file)) /= 'none' ) then
833  write(10, fmt=trim(fmt1)) 'initmip_bmb_anom file = '//initmip_bmb_anom_file
834 end if
835 #endif
836 write(10, fmt=trim(fmt1)) ' '
837 #endif
838 
839 write(10, fmt=trim(fmt2)) 'REBOUND =', rebound
840 #if (REBOUND==1)
841 write(10, fmt=trim(fmt3)) 'frac_llra =', frac_llra
842 #endif
843 #if (REBOUND==1 || REBOUND==2)
844 write(10, fmt=trim(fmt2)) 'TIME_LAG_MOD =', time_lag_mod
845 #if (TIME_LAG_MOD==1)
846 write(10, fmt=trim(fmt3)) 'time_lag =', time_lag
847 #elif (TIME_LAG_MOD==2)
848 write(10, fmt=trim(fmt1)) 'time_lag_file = '//time_lag_file
849 #else
850 stop ' >>> sico_init: TIME_LAG_MOD must be either 1 or 2!'
851 #endif
852 #endif
853 #if (REBOUND==2)
854 write(10, fmt=trim(fmt2)) 'FLEX_RIG_MOD =', flex_rig_mod
855 #if (FLEX_RIG_MOD==1)
856 write(10, fmt=trim(fmt3)) 'flex_rig =', flex_rig
857 #elif (FLEX_RIG_MOD==2)
858 write(10, fmt=trim(fmt1)) 'flex_rig_file = '//flex_rig_file
859 #else
860 stop ' >>> sico_init: FLEX_RIG_MOD must be either 1 or 2!'
861 #endif
862 #endif
863 write(10, fmt=trim(fmt2)) 'Q_LITHO =', q_litho
864 write(10, fmt=trim(fmt1)) ' '
865 
866 #if (FLOW_LAW==2)
867 write(10, fmt=trim(fmt3)) 'gr_size =', gr_size
868 write(10, fmt=trim(fmt1)) ' '
869 #endif
870 #if (FIN_VISC==2)
871 write(10, fmt=trim(fmt3)) 'sigma_res =', sigma_res
872 write(10, fmt=trim(fmt1)) ' '
873 #endif
874 
875 write(10, fmt=trim(fmt2)) 'ENHMOD =', enhmod
876 #if (ENHMOD==1 || ENHMOD==2 || ENHMOD==3)
877 write(10, fmt=trim(fmt3)) 'enh_fact =', enh_fact
878 #endif
879 #if (ENHMOD==2 || ENHMOD==3)
880 write(10, fmt=trim(fmt3)) 'enh_intg =', enh_intg
881 #endif
882 #if (ENHMOD==2)
883 write(10, fmt=trim(fmt3)) 'age_trans =', age_trans_0
884 #endif
885 #if (ENHMOD==3)
886 write(10, fmt=trim(fmt3)) 'date_trans1 =', date_trans1_0
887 write(10, fmt=trim(fmt3)) 'date_trans2 =', date_trans2_0
888 write(10, fmt=trim(fmt3)) 'date_trans3 =', date_trans3_0
889 #endif
890 #if (ENHMOD==4 || ENHMOD==5)
891 write(10, fmt=trim(fmt3)) 'enh_compr =', enh_compr
892 write(10, fmt=trim(fmt3)) 'enh_shear =', enh_shear
893 #endif
894 #if ((ENHMOD==1 || ENHMOD==2 || ENHMOD==3 || ENHMOD==4) && MARGIN==3)
895 write(10, fmt=trim(fmt3)) 'enh_shelf =', enh_shelf
896 #endif
897 write(10, fmt=trim(fmt1)) ' '
898 
899 write(10, fmt=trim(fmt3)) 'numdiff_H_t =', numdiff_h_t
900 write(10, fmt=trim(fmt3)) 'tau_cts =', tau_cts
901 write(10, fmt=trim(fmt3)) 'vh_max =', vh_max
902 write(10, fmt=trim(fmt3)) 'hd_min =', hd_min
903 write(10, fmt=trim(fmt3)) 'hd_max =', hd_max
904 #if (defined(QBM_MIN))
905 write(10, fmt=trim(fmt3)) 'qbm_min =', qbm_min
906 #elif (defined(QB_MIN)) /* obsolete */
907 write(10, fmt=trim(fmt3)) 'qb_min =', qb_min
908 #endif
909 #if (defined(QBM_MAX))
910 write(10, fmt=trim(fmt3)) 'qbm_max =', qbm_max
911 #elif (defined(QB_MAX)) /* obsolete */
912 write(10, fmt=trim(fmt3)) 'qb_max =', qb_max
913 #endif
914 write(10, fmt=trim(fmt3)) 'age_min =', age_min
915 write(10, fmt=trim(fmt3)) 'age_max =', age_max
916 write(10, fmt=trim(fmt3)) 'mean_accum =', mean_accum
917 #if (ADV_VERT==1)
918 write(10, fmt=trim(fmt3)) 'age_diff =', agediff
919 #endif
920 write(10, fmt=trim(fmt1)) ' '
921 
922 write(10, fmt=trim(fmt2a)) 'DYNAMICS = ', dynamics
923 #if ((DYNAMICS==1 && MARGIN==3) || DYNAMICS==2)
924 #if (defined(LIS_OPTS))
925 write(10, fmt=trim(fmt1)) 'lis_opts = '//lis_opts
926 #endif
927 #if (defined(N_ITER_SSA))
928 write(10, fmt=trim(fmt2a)) 'n_iter_ssa = ', n_iter_ssa
929 #endif
930 #if (defined(RELAX_FACT_SSA))
931 write(10, fmt=trim(fmt3)) 'relax_fact_ssa =', relax_fact_ssa
932 #endif
933 #endif
934 #if (DYNAMICS==2 && defined(RATIO_SL_THRESH))
935 write(10, fmt=trim(fmt3)) 'ratio_sl_thresh =', ratio_sl_thresh
936 #endif
937 #if (DYNAMICS==2 && defined(SSTA_SIA_WEIGH_FCT))
938 write(10, fmt=trim(fmt2a)) 'SSTA_SIA_WEIGH_FCT = ', ssta_sia_weigh_fct
939 #endif
940 write(10, fmt=trim(fmt1)) ' '
941 
942 write(10, fmt=trim(fmt2)) 'CALCMOD =', calcmod
943 #if (CALCMOD==-1 && defined(TEMP_CONST))
944 write(10, fmt=trim(fmt3)) 'TEMP_CONST =', temp_const
945 #endif
946 #if (CALCMOD==-1 && defined(AGE_CONST))
947 write(10, fmt=trim(fmt3)) 'AGE_CONST =', age_const
948 #endif
949 #if (CALCMOD==1 && defined(CTS_MELTING_FREEZING))
950 write(10, fmt=trim(fmt2)) 'CTS_MELTING_FREEZING =', cts_melting_freezing
951 #endif
952 write(10, fmt=trim(fmt2)) 'FLOW_LAW =', flow_law
953 write(10, fmt=trim(fmt2)) 'FIN_VISC =', fin_visc
954 write(10, fmt=trim(fmt2)) 'MARGIN =', margin
955 #if (MARGIN==2)
956 write(10, fmt=trim(fmt2)) 'MARINE_ICE_FORMATION =', marine_ice_formation
957 write(10, fmt=trim(fmt2)) 'MARINE_ICE_CALVING =', marine_ice_calving
958 #elif (MARGIN==3)
959 write(10, fmt=trim(fmt2)) 'ICE_SHELF_CALVING =', ice_shelf_calving
960 #endif
961 #if (defined(THK_EVOL))
962 write(10, fmt=trim(fmt2)) 'THK_EVOL =', thk_evol
963 #else
964 stop ' >>> sico_init: Define THK_EVOL in header file!'
965 #endif
966 #if (defined(CALCTHK))
967 write(10, fmt=trim(fmt2)) 'CALCTHK =', calcthk
968 #else
969 stop ' >>> sico_init: Define CALCTHK in header file!'
970 #endif
971 write(10, fmt=trim(fmt2)) 'ADV_HOR =', adv_hor
972 write(10, fmt=trim(fmt2)) 'ADV_VERT =', adv_vert
973 write(10, fmt=trim(fmt2)) 'TOPOGRAD =', topograd
974 write(10, fmt=trim(fmt2)) 'TEMP_PRESENT_PARA =', temp_present_para
975 write(10, fmt=trim(fmt2)) 'TSURFACE =', tsurface
976 write(10, fmt=trim(fmt2)) 'ACCSURFACE =', accsurface
977 #if (ACCSURFACE==5)
978 write(10, fmt=trim(fmt2)) 'PRECIP_ANOM_INTERPOL =', precip_anom_interpol
979 #endif
980 write(10, fmt=trim(fmt2)) 'SOLID_PRECIP =', solid_precip
981 write(10, fmt=trim(fmt2)) 'ABLSURFACE =', ablsurface
982 #if (defined(MB_ACCOUNT))
983 write(10, fmt=trim(fmt2a)) 'MB_ACCOUNT = ', mb_account
984 #endif
985 write(10, fmt=trim(fmt1)) ' '
986 
987 #if (defined(OUT_TIMES))
988 write(10, fmt=trim(fmt2a)) 'OUT_TIMES = ', out_times
989 #endif
990 #if (defined(OUTPUT_INIT))
991 write(10, fmt=trim(fmt2a)) 'OUTPUT_INIT = ', output_init
992 #endif
993 write(10, fmt=trim(fmt2a)) 'OUTPUT = ', output
994 #if (OUTPUT==1 || OUTPUT==3)
995 write(10, fmt=trim(fmt3)) 'dtime_out =' , dtime_out0
996 #endif
997 write(10, fmt=trim(fmt3)) 'dtime_ser =' , dtime_ser0
998 #if (OUTPUT==1 || OUTPUT==2)
999 write(10, fmt=trim(fmt2a)) 'ERGDAT = ', ergdat
1000 #endif
1001 write(10, fmt=trim(fmt2a)) 'NETCDF = ', netcdf
1002 #if (OUTPUT==2 || OUTPUT==3)
1003 write(10, fmt=trim(fmt2a)) 'n_output = ', n_output
1004 do n=1, n_output
1005  write(10, fmt=trim(fmt3)) 'time_output =' , time_output0(n)
1006 end do
1007 #endif
1008 write(10, fmt=trim(fmt1)) ' '
1009 
1010 write(10, fmt=trim(fmt1)) 'Program version and date: '//version//' / '//date
1011 
1012 close(10, status='keep')
1013 
1014 !-------- Conversion of time quantities --------
1015 
1016 year_zero = year_zero*year_sec ! a --> s
1017 time_init = time_init0*year_sec ! a --> s
1018 time_end = time_end0*year_sec ! a --> s
1019 dtime = dtime0*year_sec ! a --> s
1020 dtime_temp = dtime_temp0*year_sec ! a --> s
1021 #if (REBOUND==2)
1022 dtime_wss = dtime_wss0*year_sec ! a --> s
1023 #endif
1024 dtime_ser = dtime_ser0*year_sec ! a --> s
1025 #if (OUTPUT==1 || OUTPUT==3)
1026 dtime_out = dtime_out0*year_sec ! a --> s
1027 #endif
1028 #if (OUTPUT==2 || OUTPUT==3)
1029 do n=1, n_output
1030  time_output(n) = time_output0(n)*year_sec ! a --> s
1031 end do
1032 #endif
1033 
1034 if (.not.approx_integer_multiple(dtime_temp, dtime, eps_sp_dp)) &
1035  stop ' >>> sico_init: dtime_temp must be a multiple of dtime!'
1036 
1037 #if (REBOUND==2)
1038 if (.not.approx_integer_multiple(dtime_wss, dtime, eps_sp_dp)) &
1039  stop ' >>> sico_init: dtime_wss must be a multiple of dtime!'
1040 #endif
1041 
1042 if (.not.approx_integer_multiple(dtime_ser, dtime, eps_sp_dp)) &
1043  stop ' >>> sico_init: dtime_ser must be a multiple of dtime!'
1044 
1045 #if (OUTPUT==1 || OUTPUT==3)
1046 if (.not.approx_integer_multiple(dtime_out, dtime, eps_sp_dp)) &
1047  stop ' >>> sico_init: dtime_out must be a multiple of dtime!'
1048 #endif
1049 
1050 #if (THK_EVOL==2)
1051 time_target_topo_init = time_target_topo_init0 *year_sec ! a --> s
1052 time_target_topo_final = time_target_topo_final0*year_sec ! a --> s
1053 #endif
1054 
1055 #if (THK_EVOL==3)
1056 target_topo_tau = target_topo_tau0 *year_sec ! a --> s
1057 #endif
1058 
1059 
1060 time = time_init
1061 
1062 !-------- Reading of present mean-annual precipitation rate --------
1063 
1064 #if (GRID==0 || GRID==1)
1065 
1066 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1067  trim(precip_present_file)
1068 
1069 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1070 
1071 #elif (GRID==2)
1072 
1073 stop ' >>> sico_init: GRID==2 not allowed for this application!'
1074 
1075 #endif
1076 
1077 if (ios /= 0) stop ' >>> sico_init: Error when opening the precip file!'
1078 
1079 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1080 
1081 do j=jmax, 0, -1
1082  read(21, fmt=*) (precip_ma_present(j,i), i=0,imax)
1083 end do
1084 
1085 close(21, status='keep')
1086 
1087 precip_ma_present = precip_ma_present *(1.0e-03_dp/year_sec)*(rho_w/rho)
1088  ! mm/a water equiv. -> m/s ice equiv.
1089 
1090 precip_factor = 1.0_dp
1091 
1092 #if (defined(PRECIP_FACTOR_FILE))
1093 
1094 if ( trim(adjustl(precip_factor_file)) /= 'none' ) then
1095 
1096  filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1097  trim(precip_factor_file)
1098 
1099  open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1100 
1101  if (ios /= 0) stop ' >>> sico_init: Error when opening the precip factor file!'
1102 
1103  do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1104 
1105  do j=jmax, 0, -1
1106  read(21, fmt=*) (precip_factor(j,i), i=0,imax)
1107  end do
1108 
1109  close(21, status='keep')
1110 
1111 end if
1112 
1113 #endif
1114 
1115 precip_ma_present = precip_ma_present * precip_factor
1116 
1117 !-------- Present monthly precipitation rates
1118 ! (assumed to be equal to the mean annual precipitation rate) --------
1119 
1120 do n=1, 12 ! month counter
1121 do i=0, imax
1122 do j=0, jmax
1123  precip_present(j,i,n) = precip_ma_present(j,i)
1124 end do
1125 end do
1126 end do
1127 
1128 !-------- Reading of LGM mean-annual precipitation-rate anomaly --------
1129 
1130 #if (ACCSURFACE==5)
1131 
1132 #if (GRID==0 || GRID==1)
1133 
1134 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1135  trim(precip_anom_file)
1136 
1137 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1138 
1139 #endif
1140 
1141 if (ios /= 0) stop ' >>> sico_init: Error when opening the precip anomaly file!'
1142 
1143 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1144 
1145 do j=jmax, 0, -1
1146  read(21, fmt=*) (precip_ma_lgm_anom(j,i), i=0,imax)
1147 end do
1148 
1149 close(21, status='keep')
1150 
1151 precip_ma_lgm_anom = precip_ma_lgm_anom * precip_anom_fact
1152 
1153 !-------- LGM monthly precipitation-rate anomalies (assumed to be
1154 ! equal to the mean annual precipitation-rate anomaly) --------
1155 
1156 do n=1, 12 ! month counter
1157 do i=0, imax
1158 do j=0, jmax
1159 
1160  precip_lgm_anom(j,i,n) = max(precip_ma_lgm_anom(j,i), eps)
1161  ! positive values ensured
1162 
1163 #if (PRECIP_ANOM_INTERPOL==1)
1164  gamma_precip_lgm_anom(j,i,n) = 0.0_dp ! dummy values
1165 #elif (PRECIP_ANOM_INTERPOL==2)
1166  gamma_precip_lgm_anom(j,i,n) = -log(precip_lgm_anom(j,i,n))
1167 #else
1168  stop ' >>> sico_init: Wrong value of switch PRECIP_ANOM_INTERPOL!'
1169 #endif
1170 
1171 end do
1172 end do
1173 end do
1174 
1175 #endif
1176 
1177 !-------- Mean accumulation --------
1178 
1179 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho)
1180  ! mm/a water equiv. -> m/s ice equiv.
1181 
1182 !-------- Read present topography mask --------
1183 
1184 #if (GRID==0 || GRID==1)
1185 
1186 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1187  trim(mask_present_file)
1188 
1189 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
1190 
1191 #elif (GRID==2)
1192 
1193 stop ' >>> sico_init: GRID==2 not allowed for this application!'
1194 
1195 #endif
1196 
1197 if (ios /= 0) stop ' >>> sico_init: Error when opening the mask file!'
1198 
1199 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1200 
1201 do j=jmax, 0, -1
1202  read(24, fmt=trim(fmt4)) (maske_ref(j,i), i=0,imax)
1203 end do
1204 
1205 close(24, status='keep')
1206 
1207 !-------- Read rock/sediment mask --------
1208 
1209 #if (defined(SEDI_SLIDE) && SEDI_SLIDE==2)
1210 
1211 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1212  trim(mask_sedi_file)
1213 
1214 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
1215 
1216 if (ios /= 0) stop ' >>> sico_init: Error when opening the rock/sediment mask file!'
1217 
1218 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1219 
1220 do j=jmax, 0, -1
1221  read(24, fmt=trim(fmt4)) (maske_sedi(j,i), i=0,imax)
1222 end do
1223 
1224 close(24, status='keep')
1225 
1226 do i=0, imax
1227 do j=0, jmax
1228  if (maske_sedi(j,i) == 2) maske_sedi(j,i) = 7
1229  ! ocean to be treated like soft sediment
1230 end do
1231 end do
1232 
1233 #endif
1234 
1235 !-------- Initialisation of the marker for the different sectors for
1236 ! ice shelf basal melting --------
1237 
1238 n_sector = 0_i2b
1239 
1240 !-------- Present reference elevation (for precipitation data) --------
1241 
1242 #if (GRID==0 || GRID==1)
1243 
1244 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1245  trim(zs_present_file)
1246 
1247 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1248 
1249 #elif (GRID==2)
1250 
1251 stop ' >>> sico_init: GRID==2 not allowed for this application!'
1252 
1253 #endif
1254 
1255 if (ios /= 0) stop ' >>> sico_init: Error when opening the zs file!'
1256 
1257 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1258 
1259 do j=jmax, 0, -1
1260  read(21, fmt=*) (zs_ref(j,i), i=0,imax)
1261 end do
1262 
1263 close(21, status='keep')
1264 
1265 ! ------ Reset bathymetry data (sea floor elevation) to the
1266 ! sea surface
1267 
1268 do i=0, imax
1269 do j=0, jmax
1270  if (maske_ref(j,i) >= 2) zs_ref(j,i) = 0.0_dp
1271 end do
1272 end do
1273 
1274 !-------- Reading of the prescribed target topography --------
1275 
1276 #if ((THK_EVOL==2) || (THK_EVOL==3))
1277 
1278 target_topo_dat_name = trim(target_topo_dat_name)
1279 
1280 #if (NETCDF==1)
1281 write(6, fmt='(a)') ' >>> sico_init: Reading of target topography'
1282 write(6, fmt='(a)') ' only implemented for NetCDF files (NETCDF==2)!'
1283 stop
1284 #elif (NETCDF==2)
1285 call read_target_topo_nc(target_topo_dat_name)
1286 #else
1287 stop ' >>> sico_init: Parameter NETCDF must be either 1 or 2!'
1288 #endif
1289 
1290 #endif
1291 
1292 !-------- Reading of the maximum ice extent mask --------
1293 
1294 #if (THK_EVOL==4)
1295 
1296 #if (GRID==0 || GRID==1)
1297 
1298 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1299  trim(mask_maxextent_file)
1300 
1301 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
1302 
1303 #elif (GRID==2)
1304 
1305 stop ' >>> sico_init: GRID==2 not allowed for this application!'
1306 
1307 #endif
1308 
1309 if (ios /= 0) &
1310  stop ' >>> sico_init: Error when opening the maximum ice extent mask file!'
1311 
1312 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1313 
1314 do j=jmax, 0, -1
1315  read(24, fmt=trim(fmt4)) (maske_maxextent(j,i), i=0,imax)
1316 end do
1317 
1318 close(24, status='keep')
1319 
1320 #endif
1321 
1322 !-------- Reading of LGM mean-annual and mean-January (summer)
1323 ! surface-temperature anomalies --------
1324 
1325 #if (TSURFACE==5)
1326 
1327 #if (GRID==0 || GRID==1)
1328 
1329 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1330  trim(temp_ma_anom_file)
1331 
1332 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1333 
1334 #endif
1335 
1336 if (ios /= 0) stop ' >>> sico_init: Error when opening the temp_ma anomaly file!'
1337 
1338 #if (GRID==0 || GRID==1)
1339 
1340 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1341  trim(temp_mj_anom_file)
1342 
1343 open(22, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1344 
1345 #endif
1346 
1347 if (ios /= 0) stop ' >>> sico_init: Error when opening the temp_mj anomaly file!'
1348 
1349 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1350 do n=1, 6; read(22, fmt='(a)') ch_dummy; end do
1351 
1352 do j=jmax, 0, -1
1353  read(21, fmt=*) (temp_ma_lgm_anom(j,i), i=0,imax)
1354  read(22, fmt=*) (temp_mj_lgm_anom(j,i), i=0,imax)
1355 end do
1356 
1357 close(21, status='keep')
1358 close(22, status='keep')
1359 
1360 temp_ma_lgm_anom = temp_ma_lgm_anom * temp_ma_anom_fact
1361 temp_mj_lgm_anom = temp_mj_lgm_anom * temp_mj_anom_fact
1362 
1363 #endif
1364 
1365 !-------- Read data for delta_ts --------
1366 
1367 #if (TSURFACE==4)
1368 
1369 filename_with_path = trim(inpath)//'/general/'//trim(grip_temp_file)
1370 
1371 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1372 
1373 if (ios /= 0) stop ' >>> sico_init: Error when opening the data file for delta_ts!'
1374 
1375 read(21, fmt=*) ch_dummy, grip_time_min, grip_time_stp, grip_time_max
1376 
1377 if (ch_dummy /= '#') then
1378  write(6, fmt=*) ' >>> sico_init: grip_time_min, grip_time_stp, grip_time_max'
1379  write(6, fmt=*) ' not defined in data file for delta_ts!'
1380  stop
1381 end if
1382 
1384 
1385 allocate(griptemp(0:ndata_grip))
1386 
1387 do n=0, ndata_grip
1388  read(21, fmt=*) d_dummy, griptemp(n)
1389 end do
1390 
1391 close(21, status='keep')
1392 
1393 #endif
1394 
1395 !-------- Read data for the glacial index --------
1396 
1397 #if (TSURFACE==5)
1398 
1399 filename_with_path = trim(inpath)//'/general/'//trim(glac_ind_file)
1400 
1401 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1402 
1403 if (ios /= 0) stop ' >>> sico_init: Error when opening the glacial-index file!'
1404 
1405 read(21, fmt=*) ch_dummy, gi_time_min, gi_time_stp, gi_time_max
1406 
1407 if (ch_dummy /= '#') then
1408  write(6, fmt=*) ' >>> sico_init: gi_time_min, gi_time_stp, gi_time_max'
1409  write(6, fmt=*) ' not defined in glacial-index file!'
1410  stop
1411 end if
1412 
1414 
1415 allocate(glacial_index(0:ndata_gi))
1416 
1417 do n=0, ndata_gi
1418  read(21, fmt=*) d_dummy, glacial_index(n)
1419 end do
1420 
1421 close(21, status='keep')
1422 
1423 #endif
1424 
1425 !-------- Reading of time-dependent surface-temperature
1426 ! and precipitation anomaly data --------
1427 
1428 #if (TSURFACE==6 && ACCSURFACE==6)
1429 
1430 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)// &
1431  '/'//trim(temp_precip_anom_file)
1432 
1433 call check( nf90_open(trim(filename_with_path), nf90_nowrite, &
1434  ncid_temp_precip), thisroutine )
1435 
1436 call check( nf90_inquire(ncid_temp_precip, unlimiteddimid = dimid), &
1437  thisroutine )
1438 
1439 call check( nf90_inquire_dimension(ncid_temp_precip, dimid, &
1440  len = ndata_temp_precip), thisroutine )
1441 
1443 
1445 
1446 call check( nf90_inq_varid(ncid_temp_precip, 'time', ncv), thisroutine )
1447 call check( nf90_get_var(ncid_temp_precip, ncv, temp_precip_time), thisroutine )
1448  ! in a
1449 
1452  ! in a; constant time step assumed
1454 
1455 #endif
1456 
1457 !-------- Reading of ISMIP6 InitMIP SMB and BMB anomaly data --------
1458 
1459 ! ------ SMB
1460 
1461 #if (defined(INITMIP_SMB_ANOM_FILE))
1462 
1463 if ( trim(adjustl(initmip_smb_anom_file)) /= 'none' ) then
1464 
1465  filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)// &
1466  '/'//trim(initmip_smb_anom_file)
1467 
1468  call check( nf90_open(trim(filename_with_path), nf90_nowrite, ncid), &
1469  thisroutine )
1470 
1471  call check( nf90_inq_varid(ncid, 'asmb', ncv) )
1472  call check( nf90_get_var(ncid, ncv, as_anom_initmip_conv) )
1473 
1474  call check( nf90_close(ncid) )
1475 
1476  do i=0, imax
1477  do j=0, jmax
1478  as_anom_initmip(j,i) = real(as_anom_initmip_conv(i,j),dp) /YEAR_SEC
1479  ! m/a ice equiv. -> m/s ice equiv.
1480  end do
1481  end do
1482 
1483 else
1484  as_anom_initmip = 0.0_dp
1485 end if
1486 
1487 #endif
1488 
1489 ! ------ BMB
1490 
1491 #if (defined(INITMIP_BMB_ANOM_FILE))
1492 
1493 if ( trim(adjustl(initmip_bmb_anom_file)) /= 'none' ) then
1494 
1495  filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)// &
1496  '/'//trim(initmip_bmb_anom_file)
1497 
1498  call check( nf90_open(trim(filename_with_path), nf90_nowrite, ncid), &
1499  thisroutine )
1500 
1501  call check( nf90_inq_varid(ncid, 'abmb', ncv) )
1502  call check( nf90_get_var(ncid, ncv, ab_anom_initmip_conv) )
1503 
1504  call check( nf90_close(ncid) )
1505 
1506  do i=0, imax
1507  do j=0, jmax
1508  ab_anom_initmip(j,i) = real(ab_anom_initmip_conv(i,j),dp) /YEAR_SEC
1509  ! m/a ice equiv. -> m/s ice equiv.
1510  end do
1511  end do
1512 
1513 else
1514  ab_anom_initmip = 0.0_dp
1515 end if
1516 
1517 #endif
1518 
1519 !-------- Read data for z_sl --------
1520 
1521 #if (SEA_LEVEL==3)
1522 
1523 filename_with_path = trim(inpath)//'/general/'//trim(sea_level_file)
1524 
1525 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1526 
1527 if (ios /= 0) stop ' >>> sico_init: Error when opening the data file for z_sl!'
1528 
1529 read(21, fmt=*) ch_dummy, specmap_time_min, specmap_time_stp, specmap_time_max
1530 
1531 if (ch_dummy /= '#') then
1532  write(6, fmt=*) ' >>> sico_init: specmap_time_min, specmap_time_stp, specmap_time_max'
1533  write(6, fmt=*) ' not defined in data file for z_sl!'
1534  stop
1535 end if
1536 
1538 
1539 allocate(specmap_zsl(0:ndata_specmap))
1540 
1541 do n=0, ndata_specmap
1542  read(21, fmt=*) d_dummy, specmap_zsl(n)
1543 end do
1544 
1545 close(21, status='keep')
1546 
1547 #endif
1548 
1549 !-------- Determination of the geothermal heat flux --------
1550 
1551 #if (Q_GEO_MOD==1)
1552 
1553 ! ------ Constant value
1554 
1555 do i=0, imax
1556 do j=0, jmax
1557  q_geo(j,i) = q_geo *1.0e-03_dp ! mW/m2 -> W/m2
1558 end do
1559 end do
1560 
1561 #elif (Q_GEO_MOD==2)
1562 
1563 ! ------ Read data from file
1564 
1565 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1566  trim(q_geo_file)
1567 
1568 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1569 
1570 if (ios /= 0) stop ' >>> sico_init: Error when opening the qgeo file!'
1571 
1572 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1573 
1574 do j=jmax, 0, -1
1575  read(21, fmt=*) (q_geo(j,i), i=0,imax)
1576 end do
1577 
1578 close(21, status='keep')
1579 
1580 do i=0, imax
1581 do j=0, jmax
1582  q_geo(j,i) = q_geo(j,i) *1.0e-03_dp ! mW/m2 -> W/m2
1583 end do
1584 end do
1585 
1586 #endif
1587 
1588 !-------- Reading of tabulated kei function--------
1589 
1590 #if (REBOUND==0 || REBOUND==1)
1591 
1592 kei = 0.0_dp; n_data_kei = 0; kei_r_max = 0.0_dp; kei_r_incr = 0.0_dp
1593  ! dummy values
1594 #elif (REBOUND==2)
1595 
1596 call read_kei()
1597 
1598 #endif
1599 
1600 !-------- Determination of the time lag
1601 ! of the relaxing asthenosphere --------
1602 
1603 #if (REBOUND==1 || REBOUND==2)
1604 
1605 #if (TIME_LAG_MOD==1)
1606 
1607 time_lag_asth = time_lag*year_sec ! a -> s
1608 
1609 #elif (TIME_LAG_MOD==2)
1610 
1611 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1612  trim(time_lag_file)
1613 
1614 open(29, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1615 
1616 if (ios /= 0) stop ' >>> sico_init: Error when opening the time-lag file!'
1617 
1618 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1619 
1620 do j=jmax, 0, -1
1621  read(29, fmt=*) (time_lag_asth(j,i), i=0,imax)
1622 end do
1623 
1624 close(29, status='keep')
1625 
1626 time_lag_asth = time_lag_asth*year_sec ! a -> s
1627 
1628 #endif
1629 
1630 #elif (REBOUND==0)
1631 
1632 time_lag_asth = 0.0_dp ! dummy values
1633 
1634 #endif
1635 
1636 !-------- Determination of the flexural rigidity of the lithosphere --------
1637 
1638 #if (REBOUND==2)
1639 
1640 #if (FLEX_RIG_MOD==1)
1641 
1642 flex_rig_lith = flex_rig
1643 
1644 #elif (FLEX_RIG_MOD==2)
1645 
1646 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1647  trim(flex_rig_file)
1648 
1649 open(29, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1650 
1651 if (ios /= 0) stop ' >>> sico_init: Error when opening the flex-rig file!'
1652 
1653 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1654 
1655 do j=jmax, 0, -1
1656  read(29, fmt=*) (flex_rig_lith(j,i), i=0,imax)
1657 end do
1658 
1659 close(29, status='keep')
1660 
1661 #endif
1662 
1663 #elif (REBOUND==0 || REBOUND==1)
1664 
1665 flex_rig_lith = 0.0_dp ! dummy values
1666 
1667 #endif
1668 
1669 !-------- Definition of initial values --------
1670 
1671 ! ------ Present topography
1672 
1673 #if (ANF_DAT==1)
1674 
1675 call topography1(dxi, deta)
1676 
1677 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1678 
1679 call boundary(time_init, dtime, dxi, deta, &
1680  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1681 
1682 where ((maske==0_i2b).or.(maske==3_i2b))
1683  ! grounded or floating ice
1685 elsewhere ! maske==1_i2b or 2_i2b, ice-free land or sea
1686  as_perp_apl = 0.0_dp
1687 end where
1688 
1689 q_bm = 0.0_dp
1690 q_tld = 0.0_dp
1691 q_b_tot = 0.0_dp
1692 
1693 p_b_w = 0.0_dp
1694 h_w = 0.0_dp
1695 
1696 #if (!defined(TEMP_INIT) || TEMP_INIT==1)
1698 #elif (TEMP_INIT==2)
1700 #elif (TEMP_INIT==3)
1702 #elif (TEMP_INIT==4)
1704 #else
1705  write(6, fmt='(a)') ' >>> sico_init:'
1706  write(6, fmt='(a)') ' TEMP_INIT must be set to either 1, 2, 3 or 4!'
1707  stop
1708 #endif
1709 
1710 #if (ENHMOD==1)
1711  call calc_enhance_1()
1712 #elif (ENHMOD==2)
1713  call calc_enhance_2()
1714 #elif (ENHMOD==3)
1715  call calc_enhance_3(time_init)
1716 #elif (ENHMOD==4)
1717  call calc_enhance_4()
1718 #elif (ENHMOD==5)
1719  call calc_enhance_5()
1720 #else
1721  stop ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
1722 #endif
1723 
1724 ! ------ Ice-free, relaxed bedrock
1725 
1726 #elif (ANF_DAT==2)
1727 
1728 call topography2(dxi, deta)
1729 
1730 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1731 
1732 call boundary(time_init, dtime, dxi, deta, &
1733  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1734 
1735 as_perp_apl = 0.0_dp
1736 
1737 q_bm = 0.0_dp
1738 q_tld = 0.0_dp
1739 q_b_tot = 0.0_dp
1740 
1741 p_b_w = 0.0_dp
1742 h_w = 0.0_dp
1743 
1744 call init_temp_water_age_2()
1745 
1746 #if (ENHMOD==1)
1747  call calc_enhance_1()
1748 #elif (ENHMOD==2)
1749  call calc_enhance_2()
1750 #elif (ENHMOD==3)
1751  call calc_enhance_3(time_init)
1752 #elif (ENHMOD==4)
1753  call calc_enhance_4()
1754 #elif (ENHMOD==5)
1755  call calc_enhance_5()
1756 #else
1757  stop ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
1758 #endif
1759 
1760 ! ------ Read initial state from output of previous simulation
1761 
1762 #elif (ANF_DAT==3)
1763 
1764 call topography3(dxi, deta, z_sl, anfdatname)
1765 
1766 call boundary(time_init, dtime, dxi, deta, &
1767  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1768 
1769 where ((maske==0_i2b).or.(maske==3_i2b))
1770  ! grounded or floating ice
1772 elsewhere ! maske==1_i2b or 2_i2b, ice-free land or sea
1773  as_perp_apl = 0.0_dp
1774 end where
1775 
1776 q_b_tot = q_bm + q_tld
1777 
1778 #endif
1779 
1780 !-------- Inner-point flag --------
1781 
1782 flag_inner_point = .true.
1783 
1784 flag_inner_point(0,:) = .false.
1785 flag_inner_point(jmax,:) = .false.
1786 
1787 flag_inner_point(:,0) = .false.
1788 flag_inner_point(:,imax) = .false.
1789 
1790 !-------- Grounding line and calving front flags --------
1791 
1792 flag_grounding_line_1 = .false.
1793 flag_grounding_line_2 = .false.
1794 
1795 flag_calving_front_1 = .false.
1796 flag_calving_front_2 = .false.
1797 
1798 #if (MARGIN==1 || MARGIN==2)
1799 
1800 continue
1801 
1802 #elif (MARGIN==3)
1803 
1804 do i=1, imax-1
1805 do j=1, jmax-1
1806 
1807  if ( (maske(j,i)==0_i2b) & ! grounded ice
1808  .and. &
1809  ( (maske(j,i+1)==3_i2b) & ! with
1810  .or.(maske(j,i-1)==3_i2b) & ! one
1811  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1812  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1813  ) &
1814  flag_grounding_line_1(j,i) = .true.
1815 
1816  if ( (maske(j,i)==3_i2b) & ! floating ice
1817  .and. &
1818  ( (maske(j,i+1)==0_i2b) & ! with
1819  .or.(maske(j,i-1)==0_i2b) & ! one
1820  .or.(maske(j+1,i)==0_i2b) & ! neighbouring
1821  .or.(maske(j-1,i)==0_i2b) ) & ! grounded ice point
1822  ) &
1823  flag_grounding_line_2(j,i) = .true.
1824 
1825  if ( (maske(j,i)==3_i2b) & ! floating ice
1826  .and. &
1827  ( (maske(j,i+1)==2_i2b) & ! with
1828  .or.(maske(j,i-1)==2_i2b) & ! one
1829  .or.(maske(j+1,i)==2_i2b) & ! neighbouring
1830  .or.(maske(j-1,i)==2_i2b) ) & ! ocean point
1831  ) &
1832  flag_calving_front_1(j,i) = .true.
1833 
1834  if ( (maske(j,i)==2_i2b) & ! ocean
1835  .and. &
1836  ( (maske(j,i+1)==3_i2b) & ! with
1837  .or.(maske(j,i-1)==3_i2b) & ! one
1838  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1839  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1840  ) &
1841  flag_calving_front_2(j,i) = .true.
1842 
1843 end do
1844 end do
1845 
1846 do i=1, imax-1
1847 
1848  j=0
1849  if ( (maske(j,i)==2_i2b) & ! ocean
1850  .and. (maske(j+1,i)==3_i2b) & ! with one neighbouring
1851  ) & ! floating ice point
1852  flag_calving_front_2(j,i) = .true.
1853 
1854  j=jmax
1855  if ( (maske(j,i)==2_i2b) & ! ocean
1856  .and. (maske(j-1,i)==3_i2b) & ! with one neighbouring
1857  ) & ! floating ice point
1858  flag_calving_front_2(j,i) = .true.
1859 
1860 end do
1861 
1862 do j=1, jmax-1
1863 
1864  i=0
1865  if ( (maske(j,i)==2_i2b) & ! ocean
1866  .and. (maske(j,i+1)==3_i2b) & ! with one neighbouring
1867  ) & ! floating ice point
1868  flag_calving_front_2(j,i) = .true.
1869 
1870  i=imax
1871  if ( (maske(j,i)==2_i2b) & ! ocean
1872  .and. (maske(j,i-1)==3_i2b) & ! with one neighbouring
1873  ) & ! floating ice point
1874  flag_calving_front_2(j,i) = .true.
1875 
1876 end do
1877 
1878 #else
1879 stop ' >>> sico_init: MARGIN must be either 1, 2 or 3!'
1880 #endif
1881 
1882 !-------- Rock/sediment mask: smoothing over the transition length --------
1883 
1884 #if (defined(SEDI_SLIDE) && SEDI_SLIDE==2 && defined(TRANS_LENGTH_SL))
1885 
1886 do i=0, imax
1887 do j=0, jmax
1888  if (maske_sedi(j,i) /= 7) then ! no soft sediment
1889  r_mask_sedi(j,i) = 0.0_dp
1890  else ! maske_sedi(j,i) == 7, soft sediment
1891  r_mask_sedi(j,i) = 1.0_dp
1892  end if
1893 end do
1894 end do
1895 
1896 if (trans_length_sl > eps) then
1897 
1898  if ((p_weert /= p_weert_sedi).or.(q_weert /= q_weert_sedi)) then
1899  write(6, fmt='(a)') ' >>> sico_init: Specifying a transition length'
1900  write(6, fmt='(a)') ' TRANS_LENGTH_SL > 0 for the basal sliding'
1901  write(6, fmt='(a)') ' regimes requires P_WEERT==P_WEERT_SEDI'
1902  write(6, fmt='(a)') ' and Q_WEERT==Q_WEERT_SEDI!'
1903  stop
1904  end if
1905 
1906  transition_length_sliding = trans_length_sl *1.0e+03_dp ! km -> m
1907 
1908  quasi_time = 0.25_dp*transition_length_sliding**2
1909  ! Time scale resulting from the fundamental solution
1910  ! of the diffusion equation for diffusitivity = 1
1911  d_quasi_time = 0.1_dp*quasi_time
1912  ! First guess for the time step
1913 
1914  do
1915  if ( d_quasi_time < 0.1_dp*min(dxi,deta)**2 ) exit
1916  ! CFL condition fulfilled (with some leeway)
1917  d_quasi_time = 0.5_dp*d_quasi_time
1918  end do
1919 
1920  m = nint(quasi_time/d_quasi_time)
1921 
1922  do n=1, m
1923 
1924  r_mask_sedi_new = r_mask_sedi
1925 
1926  do i=1, imax-1
1927  do j=1, jmax-1
1928  r_mask_sedi_new(j,i) = r_mask_sedi(j,i) &
1929  + d_quasi_time/dxi**2 &
1930  * (r_mask_sedi(j,i+1)-2.0_dp*r_mask_sedi(j,i)+r_mask_sedi(j,i-1)) &
1931  + d_quasi_time/deta**2 &
1932  * (r_mask_sedi(j+1,i)-2.0_dp*r_mask_sedi(j,i)+r_mask_sedi(j-1,i))
1933  ! Solving the diffusion equation for diffusitivity = 1
1934  end do
1935  end do
1936 
1937  r_mask_sedi = r_mask_sedi_new
1938 
1939  end do
1940 
1941 end if
1942 
1943 #endif
1944 
1945 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
1946 
1947 #if (GRID==0 || GRID==1)
1948 
1949 do ir=-imax, imax
1950 do jr=-jmax, jmax
1951  dist_dxdy(jr,ir) = sqrt( (real(ir,dp)*dxi)**2 + (real(jr,dp)*deta)**2 )
1952  ! distortion due to stereographic projection not accounted for
1953 end do
1954 end do
1955 
1956 #endif
1957 
1958 !-------- Initial velocities --------
1959 
1960 call calc_temp_melt()
1961 
1962 #if (DYNAMICS==1 || DYNAMICS==2)
1963 
1964 call calc_vxy_b_sia(time, z_sl)
1965 call calc_vxy_sia(dzeta_c, dzeta_t)
1966 
1967 #if (MARGIN==3 || DYNAMICS==2)
1968 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
1969 #endif
1970 
1971 call calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
1972 
1973 #if (MARGIN==3)
1974 call calc_vz_floating(z_sl, dxi, deta, dzeta_c)
1975 #endif
1976 
1977 #elif (DYNAMICS==0)
1978 
1979 call calc_vxy_static()
1980 call calc_vz_static()
1981 
1982 #else
1983  stop ' >>> sico_init: DYNAMICS must be either 0, 1 or 2!'
1984 #endif
1985 
1986 call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
1987 
1988 !-------- Initial enthalpies --------
1989 
1990 #if (CALCMOD==0 || CALCMOD==-1)
1991 
1992 do i=0, imax
1993 do j=0, jmax
1994 
1995  do kc=0, kcmax
1996  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1997  end do
1998 
1999  do kt=0, ktmax
2000  enth_t(kt,j,i) = enth_c(0,j,i)
2001  end do
2002 
2003 end do
2004 end do
2005 
2006 #elif (CALCMOD==1)
2007 
2008 do i=0, imax
2009 do j=0, jmax
2010 
2011  do kc=0, kcmax
2012  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
2013  end do
2014 
2015  if ( (maske(j,i) == 0_i2b).and.(n_cts(j,i) == 1_i2b) ) then
2016  do kt=0, ktmax
2017  enth_t(kt,j,i) = enth_fct_temp_omega(temp_t_m(kt,j,i), omega_t(kt,j,i))
2018  end do
2019  else
2020  do kt=0, ktmax
2021  enth_t(kt,j,i) = enth_c(0,j,i)
2022  end do
2023  end if
2024 
2025 end do
2026 end do
2027 
2028 #elif (CALCMOD==2 || CALCMOD==3)
2029 
2030 do i=0, imax
2031 do j=0, jmax
2032 
2033  do kc=0, kcmax
2034  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), omega_c(kc,j,i))
2035  end do
2036 
2037  do kt=0, ktmax
2038  enth_t(kt,j,i) = enth_c(0,j,i)
2039  end do
2040 
2041 end do
2042 end do
2043 
2044 #else
2045 
2046 stop ' >>> sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!'
2047 
2048 #endif
2049 
2050 !-------- Initialize time-series files --------
2051 
2052 ! ------ Time-series file for the ice sheet on the whole
2053 
2054 filename_with_path = trim(outpath)//'/'//trim(runname)//'.ser'
2055 
2056 open(12, iostat=ios, file=trim(filename_with_path), status='new')
2057 
2058 if (ios /= 0) stop ' >>> sico_init: Error when opening the ser file!'
2059 
2060 if (forcing_flag == 1) then
2061 
2062  write(12,1102)
2063  write(12,1103)
2064 
2065  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
2066  ' V(m^3) V_g(m^3) V_f(m^3)', &
2067  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
2068  ' V_sle(m) V_t(m^3)', &
2069  ' A_t(m^2)',/, &
2070  ' H_max(m) H_t_max(m)', &
2071  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
2072  1103 format('----------------------------------------------------', &
2073  '---------------------------------------')
2074 
2075 else if (forcing_flag == 2) then
2076 
2077  write(12,1112)
2078  write(12,1113)
2079 
2080  1112 format(' t(a) glac_ind(1) z_sl(m)',/, &
2081  ' V(m^3) V_g(m^3) V_f(m^3)', &
2082  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
2083  ' V_sle(m) V_t(m^3)', &
2084  ' A_t(m^2)',/, &
2085  ' H_max(m) H_t_max(m)', &
2086  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
2087  1113 format('----------------------------------------------------', &
2088  '---------------------------------------')
2089 
2090 else if (forcing_flag == 3) then
2091 
2092  write(12,1122)
2093  write(12,1123)
2094 
2095  1122 format(' t(a) (Dummy)(1) z_sl(m)',/, &
2096  ' V(m^3) V_g(m^3) V_f(m^3)', &
2097  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
2098  ' V_sle(m) V_t(m^3)', &
2099  ' A_t(m^2)',/, &
2100  ' H_max(m) H_t_max(m)', &
2101  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
2102  1123 format('----------------------------------------------------', &
2103  '---------------------------------------')
2104 
2105 end if
2106 
2107 ! ------ Time-series file for deep boreholes
2108 
2109 n_core = 6 ! Vostok, Dome A, Dome C, Dome F, Kohnen, Byrd
2110 
2111 allocate(lambda_core(n_core), phi_core(n_core), &
2113 
2114 ch_core(1) = 'Vostok'
2115 phi_core(1) = -78.467_dp *pi_180 ! Geographical position of Vostok,
2116 lambda_core(1) = 106.8_dp *pi_180 ! conversion deg -> rad
2117 
2118 ch_core(2) = 'Dome A'
2119 phi_core(2) = -80.367_dp *pi_180 ! Geographical position of Dome A,
2120 lambda_core(2) = 77.35_dp *pi_180 ! conversion deg -> rad
2121 
2122 ch_core(3) = 'Dome C'
2123 phi_core(3) = -75.1_dp *pi_180 ! Geographical position of Dome C,
2124 lambda_core(3) = 123.4_dp *pi_180 ! conversion deg -> rad
2125 
2126 ch_core(4) = 'Dome F'
2127 phi_core(4) = -77.317_dp *pi_180 ! Geographical position of Dome F,
2128 lambda_core(4) = 39.7_dp *pi_180 ! conversion deg -> rad
2129 
2130 ch_core(5) = 'Kohnen'
2131 phi_core(5) = -75.0_dp *pi_180 ! Geographical position of Kohnen,
2132 lambda_core(5) = 0.067_dp *pi_180 ! conversion deg -> rad
2133 
2134 ch_core(6) = 'Byrd'
2135 phi_core(6) = -80.017_dp *pi_180 ! Geographical position of Byrd,
2136 lambda_core(6) = -119.517_dp *pi_180 ! conversion deg -> rad
2137 
2138 #if (GRID==0 || GRID==1) /* Stereographic projection */
2139 
2140 do n=1, n_core
2142  a, b, lambda0, phi0, x_core(n), y_core(n))
2143 end do
2144 
2145 #elif (GRID==2) /* Geographical coordinates */
2146 
2148 y_core = phi_core
2149 
2150 #endif
2151 
2152 filename_with_path = trim(outpath)//'/'//trim(runname)//'.core'
2153 
2154 open(14, iostat=ios, file=trim(filename_with_path), status='new')
2155 
2156 if (ios /= 0) stop ' >>> sico_init: Error when opening the core file!'
2157 
2158 if (forcing_flag == 1) then
2159 
2160  write(14,1106)
2161  write(14,1107)
2162 
2163  1106 format(' t(a) D_Ts(C) z_sl(m)',/, &
2164  ' H_Vo(m) H_DA(m) H_DC(m)', &
2165  ' H_DF(m) H_Ko(m) H_By(m)',/, &
2166  ' v_Vo(m/a) v_DA(m/a) v_DC(m/a)', &
2167  ' v_DF(m/a) v_Ko(m/a) v_By(m/a)',/, &
2168  ' T_Vo(C) T_DA(C) T_DC(C)', &
2169  ' T_DF(C) T_Ko(C) T_By(C)')
2170  1107 format('----------------------------------------------------', &
2171  '---------------------------------------')
2172 
2173 else if (forcing_flag == 2) then
2174 
2175  write(14,1116)
2176  write(14,1117)
2177 
2178  1116 format(' t(a) glac_ind(1) z_sl(m)',/, &
2179  ' H_Vo(m) H_DA(m) H_DC(m)', &
2180  ' H_DF(m) H_Ko(m) H_By(m)',/, &
2181  ' v_Vo(m/a) v_DA(m/a) v_DC(m/a)', &
2182  ' v_DF(m/a) v_Ko(m/a) v_By(m/a)',/, &
2183  ' T_Vo(C) T_DA(C) T_DC(C)', &
2184  ' T_DF(C) T_Ko(C) T_By(C)')
2185  1117 format('----------------------------------------------------', &
2186  '---------------------------------------')
2187 
2188 else if (forcing_flag == 3) then
2189 
2190  write(14,1126)
2191  write(14,1127)
2192 
2193  1126 format(' t(a) (Dummy)(1) z_sl(m)',/, &
2194  ' H_Vo(m) H_DA(m) H_DC(m)', &
2195  ' H_DF(m) H_Ko(m) H_By(m)',/, &
2196  ' v_Vo(m/a) v_DA(m/a) v_DC(m/a)', &
2197  ' v_DF(m/a) v_Ko(m/a) v_By(m/a)',/, &
2198  ' T_Vo(C) T_DA(C) T_DC(C)', &
2199  ' T_DF(C) T_Ko(C) T_By(C)')
2200  1127 format('----------------------------------------------------', &
2201  '---------------------------------------')
2202 
2203 end if
2204 
2205 !-------- Output of the initial state --------
2206 
2207 #if (defined(OUTPUT_INIT))
2208 
2209 #if (OUTPUT_INIT==0)
2210  flag_init_output = .false.
2211 #elif (OUTPUT_INIT==1)
2212  flag_init_output = .true.
2213 #else
2214  stop ' >>> sico_init: OUTPUT_INIT must be either 0 or 1!'
2215 #endif
2216 
2217 #else
2218  flag_init_output = .true. ! default for undefined parameter OUTPUT_INIT
2219 #endif
2220 
2221 #if (OUTPUT==1)
2222 
2223 #if (ERGDAT==0)
2224  flag_3d_output = .false.
2225 #elif (ERGDAT==1)
2226  flag_3d_output = .true.
2227 #else
2228  stop ' >>> sico_init: ERGDAT must be either 0 or 1!'
2229 #endif
2230 
2231  if (flag_init_output) &
2232  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
2233  flag_3d_output, ndat2d, ndat3d)
2234 
2235 #elif (OUTPUT==2)
2236 
2237 if (time_output(1) <= time_init+eps) then
2238 
2239 #if (ERGDAT==0)
2240  flag_3d_output = .false.
2241 #elif (ERGDAT==1)
2242  flag_3d_output = .true.
2243 #else
2244  stop ' >>> sico_init: ERGDAT must be either 0 or 1!'
2245 #endif
2246 
2247  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
2248  flag_3d_output, ndat2d, ndat3d)
2249 
2250 end if
2251 
2252 #elif (OUTPUT==3)
2253 
2254  flag_3d_output = .false.
2255 
2256  if (flag_init_output) &
2257  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
2258  flag_3d_output, ndat2d, ndat3d)
2259 
2260 if (time_output(1) <= time_init+eps) then
2261 
2262  flag_3d_output = .true.
2263 
2264  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
2265  flag_3d_output, ndat2d, ndat3d)
2266 
2267 end if
2268 
2269 #else
2270  stop ' >>> sico_init: OUTPUT must be either 1, 2 or 3!'
2271 #endif
2272 
2273 if (flag_init_output) then
2274  call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
2275  call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
2276 end if
2277 
2278 end subroutine sico_init
2279 
2280 !-------------------------------------------------------------------------------
2281 !> Definition of the initial surface and bedrock topography
2282 !! (including gradients) and of the horizontal grid spacings dxi, deta.
2283 !! For present-day initial topography.
2284 !<------------------------------------------------------------------------------
2285 subroutine topography1(dxi, deta)
2287 #if (GRID==0 || GRID==1)
2289 #endif
2290 
2291  use metric_m
2292  use topograd_m
2293 
2294 implicit none
2295 
2296 real(dp), intent(out) :: dxi, deta
2297 
2298 integer(i4b) :: i, j, n
2299 integer(i4b) :: ios, n_dummy
2300 real(dp) :: d_dummy
2301 real(dp) :: xi0, eta0
2302 real(dp) :: H_ice, freeboard_ratio
2303 character :: ch_dummy
2304 
2305 character(len= 8) :: ch_imax
2306 character(len=128) :: fmt4
2307 character(len=256) :: filename_with_path
2308 
2309 write(ch_imax, fmt='(i8)') imax
2310 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
2311 
2312 !-------- Read topography --------
2313 
2314 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2315  trim(zs_present_file)
2316 
2317 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
2318 
2319 if (ios /= 0) stop ' >>> topography1: Error when opening the zs file!'
2320 
2321 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2322  trim(zl_present_file)
2323 
2324 open(22, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
2325 
2326 if (ios /= 0) stop ' >>> topography1: Error when opening the zl file!'
2327 
2328 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2329  trim(zl0_file)
2330 
2331 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
2332 
2333 if (ios /= 0) stop ' >>> topography1: Error when opening the zl0 file!'
2334 
2335 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2336  trim(mask_present_file)
2337 
2338 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
2339 
2340 if (ios /= 0) stop ' >>> topography1: Error when opening the mask file!'
2341 
2342 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
2343 do n=1, 6; read(22, fmt='(a)') ch_dummy; end do
2344 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
2345 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
2346 
2347 do j=jmax, 0, -1
2348  read(21, fmt=*) (zs(j,i), i=0,imax)
2349  read(22, fmt=*) (zl(j,i), i=0,imax)
2350  read(23, fmt=*) (zl0(j,i), i=0,imax)
2351  read(24, fmt=trim(fmt4)) (maske(j,i), i=0,imax)
2352 end do
2353 
2354 close(21, status='keep')
2355 close(22, status='keep')
2356 close(23, status='keep')
2357 close(24, status='keep')
2358 
2359 #if (defined(ZB_PRESENT_FILE))
2360 
2361 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2362  trim(zb_present_file)
2363 
2364 open(25, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
2365 
2366 if (ios /= 0) stop ' >>> topography1: Error when opening the zb file!'
2367 
2368 do n=1, 6; read(25, fmt='(a)') ch_dummy; end do
2369 
2370 do j=jmax, 0, -1
2371  read(25, fmt=*) (zb(j,i), i=0,imax)
2372 end do
2373 
2374 close(25, status='keep')
2375 
2376 #else
2377 
2378 write(6, fmt='(a)') ' >>> topography1: ZB_PRESENT_FILE not defined,'
2379 write(6, fmt='(a)') ' thus zb = zl assumed.'
2380 
2381 zb = zl
2382 
2383 #endif
2384 
2385 !-------- Further stuff --------
2386 
2387 dxi = dx *1000.0_dp ! km -> m
2388 deta = dx *1000.0_dp ! km -> m
2389 
2390 xi0 = x0 *1000.0_dp ! km -> m
2391 eta0 = y0 *1000.0_dp ! km -> m
2392 
2393 freeboard_ratio = (rho_sw-rho)/rho_sw
2394 
2395 do i=0, imax
2396 do j=0, jmax
2397 
2398  if (maske(j,i) <= 1) then
2399 
2400  zb(j,i) = zl(j,i) ! ensure consistency
2401 
2402  else if (maske(j,i) == 2) then
2403 
2404 #if (MARGIN==1 || MARGIN==2)
2405  zs(j,i) = zl(j,i) ! ensure
2406  zb(j,i) = zl(j,i) ! consistency
2407 #elif (MARGIN==3)
2408  zs(j,i) = 0.0_dp ! present-day
2409  zb(j,i) = 0.0_dp ! sea level
2410 #endif
2411 
2412  else if (maske(j,i) == 3) then
2413 
2414 #if (MARGIN==1 || (MARGIN==2 && MARINE_ICE_FORMATION==1))
2415  maske(j,i) = 2 ! floating ice cut off
2416  zs(j,i) = zl(j,i)
2417  zb(j,i) = zl(j,i)
2418 #elif (MARGIN==2 && MARINE_ICE_FORMATION==2)
2419  maske(j,i) = 0 ! floating ice becomes "underwater ice"
2420  h_ice = zs(j,i)-zb(j,i) ! ice thickness
2421  zs(j,i) = zl(j,i)+h_ice
2422  zb(j,i) = zl(j,i)
2423 #elif (MARGIN==3)
2424  h_ice = zs(j,i)-zb(j,i) ! ice thickness
2425  zs(j,i) = freeboard_ratio*h_ice ! ensure properly
2426  zb(j,i) = zs(j,i)-h_ice ! floating ice
2427 #endif
2428 
2429  end if
2430 
2431  xi(i) = xi0 + real(i,dp)*dxi
2432  eta(j) = eta0 + real(j,dp)*deta
2433 
2434  zm(j,i) = zb(j,i)
2435  n_cts(j,i) = -1
2436  kc_cts(j,i) = 0
2437 
2438  h_c(j,i) = zs(j,i)-zm(j,i)
2439  h_t(j,i) = 0.0_dp
2440 
2441  dzs_dtau(j,i) = 0.0_dp
2442  dzm_dtau(j,i) = 0.0_dp
2443  dzb_dtau(j,i) = 0.0_dp
2444  dzl_dtau(j,i) = 0.0_dp
2445  dh_c_dtau(j,i) = 0.0_dp
2446  dh_t_dtau(j,i) = 0.0_dp
2447 
2448 end do
2449 end do
2450 
2451 maske_old = maske
2452 
2453 !-------- Geographic coordinates, metric tensor,
2454 ! gradients of the topography --------
2455 
2456 do i=0, imax
2457 do j=0, jmax
2458 
2459 #if (GRID==0 || GRID==1) /* Stereographic projection */
2460 
2461  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
2462  lambda0, phi0, lambda(j,i), phi(j,i))
2463 
2464 #elif (GRID==2) /* Geographic coordinates */
2465 
2466  lambda(j,i) = xi(i)
2467  phi(j,i) = eta(j)
2468 
2469 #endif
2470 
2471 end do
2472 end do
2473 
2474 call metric()
2475 
2476 #if (TOPOGRAD==0)
2477 call topograd_1(dxi, deta, 1)
2478 #elif (TOPOGRAD==1)
2479 call topograd_2(dxi, deta, 1)
2480 #endif
2481 
2482 !-------- Corresponding area of grid points --------
2483 
2484 do i=0, imax
2485 do j=0, jmax
2486  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
2487 end do
2488 end do
2489 
2490 end subroutine topography1
2491 
2492 !-------------------------------------------------------------------------------
2493 !> Definition of the initial surface and bedrock topography
2494 !! (including gradients) and of the horizontal grid spacings dxi, deta.
2495 !! For ice-free initial topography with relaxed lithosphere.
2496 !<------------------------------------------------------------------------------
2497 subroutine topography2(dxi, deta)
2499 #if (GRID==0 || GRID==1)
2501 #endif
2502 
2503  use metric_m
2504  use topograd_m
2505 
2506 implicit none
2507 
2508 real(dp), intent(out) :: dxi, deta
2509 
2510 integer(i4b) :: i, j, n
2511 integer(i4b) :: ios
2512 real(dp) :: xi0, eta0
2513 character :: ch_dummy
2514 
2515 character(len= 8) :: ch_imax
2516 character(len=128) :: fmt4
2517 character(len=256) :: filename_with_path
2518 
2519 write(ch_imax, fmt='(i8)') imax
2520 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
2521 
2522 !-------- Read topography --------
2523 
2524 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2525  trim(zl0_file)
2526 
2527 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
2528 
2529 if (ios /= 0) stop ' >>> topography2: Error when opening the zl0 file!'
2530 
2531 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2532  trim(mask_present_file)
2533 
2534 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
2535 
2536 if (ios /= 0) stop ' >>> topography2: Error when opening the mask file!'
2537 
2538 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
2539 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
2540 
2541 do j=jmax, 0, -1
2542  read(23, fmt=*) (zl0(j,i), i=0,imax)
2543  read(24, fmt=trim(fmt4)) (maske(j,i), i=0,imax)
2544 end do
2545 
2546 close(23, status='keep')
2547 close(24, status='keep')
2548 
2549 !-------- Further stuff --------
2550 
2551 dxi = dx *1000.0_dp ! km -> m
2552 deta = dx *1000.0_dp ! km -> m
2553 
2554 xi0 = x0 *1000.0_dp ! km -> m
2555 eta0 = y0 *1000.0_dp ! km -> m
2556 
2557 do i=0, imax
2558 do j=0, jmax
2559 
2560  if (maske(j,i) <= 1) then
2561  maske(j,i) = 1
2562  zs(j,i) = zl0(j,i)
2563  zb(j,i) = zl0(j,i)
2564  zl(j,i) = zl0(j,i)
2565  else ! (maske(j,i) >= 2)
2566  maske(j,i) = 2
2567 #if (MARGIN==1 || MARGIN==2)
2568  zs(j,i) = zl0(j,i)
2569  zb(j,i) = zl0(j,i)
2570 #elif (MARGIN==3)
2571  zs(j,i) = 0.0_dp ! present-day
2572  zb(j,i) = 0.0_dp ! sea level
2573 #endif
2574  zl(j,i) = zl0(j,i)
2575  end if
2576 
2577  xi(i) = xi0 + real(i,dp)*dxi
2578  eta(j) = eta0 + real(j,dp)*deta
2579 
2580  zm(j,i) = zb(j,i)
2581  n_cts(j,i) = -1
2582  kc_cts(j,i) = 0
2583 
2584  h_c(j,i) = 0.0_dp
2585  h_t(j,i) = 0.0_dp
2586 
2587  dzs_dtau(j,i) = 0.0_dp
2588  dzm_dtau(j,i) = 0.0_dp
2589  dzb_dtau(j,i) = 0.0_dp
2590  dzl_dtau(j,i) = 0.0_dp
2591  dh_c_dtau(j,i) = 0.0_dp
2592  dh_t_dtau(j,i) = 0.0_dp
2593 
2594 end do
2595 end do
2596 
2597 maske_old = maske
2598 
2599 !-------- Geographic coordinates, metric tensor,
2600 ! gradients of the topography --------
2601 
2602 do i=0, imax
2603 do j=0, jmax
2604 
2605 #if (GRID==0 || GRID==1) /* Stereographic projection */
2606 
2607  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
2608  lambda0, phi0, lambda(j,i), phi(j,i))
2609 
2610 #elif (GRID==2) /* Geographic coordinates */
2611 
2612  lambda(j,i) = xi(i)
2613  phi(j,i) = eta(j)
2614 
2615 #endif
2616 
2617 end do
2618 end do
2619 
2620 call metric()
2621 
2622 #if (TOPOGRAD==0)
2623 call topograd_1(dxi, deta, 1)
2624 #elif (TOPOGRAD==1)
2625 call topograd_2(dxi, deta, 1)
2626 #endif
2627 
2628 !-------- Corresponding area of grid points --------
2629 
2630 do i=0, imax
2631 do j=0, jmax
2632  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
2633 end do
2634 end do
2635 
2636 end subroutine topography2
2637 
2638 !-------------------------------------------------------------------------------
2639 !> Definition of the initial surface and bedrock topography
2640 !! (including gradients) and of the horizontal grid spacings dxi, deta.
2641 !! For initial topography from previous simulation.
2642 !<------------------------------------------------------------------------------
2643 subroutine topography3(dxi, deta, z_sl, anfdatname)
2645  use read_m, only : read_erg_nc
2646 
2647 #if (GRID==0 || GRID==1)
2649 #endif
2650 
2651  use metric_m
2652  use topograd_m
2653 
2654 implicit none
2655 
2656 character(len=100), intent(in) :: anfdatname
2657 
2658 real(dp), intent(out) :: dxi, deta, z_sl
2659 
2660 integer(i4b) :: i, j, n
2661 integer(i4b) :: ios
2662 character(len=256) :: filename_with_path
2663 character :: ch_dummy
2664 
2665 !-------- Read data from time-slice file of previous simulation --------
2666 
2667 call read_erg_nc(z_sl, anfdatname)
2668 
2669 !-------- Read topography of the relaxed bedrock --------
2670 
2671 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2672  trim(zl0_file)
2673 
2674 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
2675 
2676 if (ios /= 0) stop ' >>> topography3: Error when opening the zl0 file!'
2677 
2678 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
2679 
2680 do j=jmax, 0, -1
2681  read(23, fmt=*) (zl0(j,i), i=0,imax)
2682 end do
2683 
2684 close(23, status='keep')
2685 
2686 !-------- Further stuff --------
2687 
2688 dxi = dx *1000.0_dp ! km -> m
2689 deta = dx *1000.0_dp ! km -> m
2690 
2691 !-------- Geographic coordinates, metric tensor,
2692 ! gradients of the topography --------
2693 
2694 do i=0, imax
2695 do j=0, jmax
2696 
2697 #if (GRID==0 || GRID==1) /* Stereographic projection */
2698 
2699  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
2700  lambda0, phi0, lambda(j,i), phi(j,i))
2701 
2702 #elif (GRID==2) /* Geographic coordinates */
2703 
2704  lambda(j,i) = xi(i)
2705  phi(j,i) = eta(j)
2706 
2707 #endif
2708 
2709 end do
2710 end do
2711 
2712 call metric()
2713 
2714 #if (TOPOGRAD==0)
2715 call topograd_1(dxi, deta, 1)
2716 #elif (TOPOGRAD==1)
2717 call topograd_2(dxi, deta, 1)
2718 #endif
2719 
2720 !-------- Corresponding area of grid points --------
2721 
2722 do i=0, imax
2723 do j=0, jmax
2724  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
2725 end do
2726 end do
2727 
2728 end subroutine topography3
2729 
2730 !-------------------------------------------------------------------------------
2731 
2732 end module sico_init_m
2733 !
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:50
real(dp), dimension(:), allocatable phi_core
phi_core(n): Geographical latitude of the prescribed borehole positions
subroutine, public read_kei()
Reading of the tabulated kei function.
Definition: read_m.F90:765
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
integer(i2b), dimension(0:jmax, 0:imax) n_sector
n_sector(j,i): Marker for the different sectors for ice shelf basal melting.
Definition: sico_vars_m.F90:55
subroutine, public calc_c_int_inv_table()
Computation of the inverse of the temperature integral of the specific heat of ice as a table (c_int_...
real(dp), dimension(0:kcmax) eaz_c
eaz_c(kc): Abbreviation for exp(aa*zeta(kc))
real(dp), dimension(0:jmax, 0:imax) q_geo
q_geo(j,i): Geothermal heat flux
subroutine sico_init(delta_ts, glac_index, mean_accum, dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, time, time_init, time_end, time_output, dxi, deta, dzeta_c, dzeta_t, dzeta_r, z_sl, dzsl_dtau, z_mar, ii, jj, nn, ndat2d, ndat3d, n_output, runname)
Main routine of sico_init_m: Initialisations for SICOPOLIS.
Definition: sico_init_m.F90:59
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, 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:59
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) temp_precip_time_stp
temp_precip_time_stp: Time step of the surface-temperature and precipitation data ...
real(dp), parameter eps
eps: Small number
subroutine, public read_phys_para()
Reading of physical parameters.
Definition: read_m.F90:813
integer(i2b), dimension(0:jmax, 0:imax) maske_maxextent
maske_maxextent(j,i): Maximum ice extent mask. 0: not allowed to glaciate, 1: allowed to glaciate...
integer(i2b) forcing_flag
forcing_flag: Flag for the forcing type. 1: forcing by a spatially constant surface temperature anoma...
real(dp) temp_precip_time_max
temp_precip_time_max: Maximum time of the surface-temperature and precipitation data ...
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
integer(i2b), dimension(0:jmax, 0:imax) maske_sedi
maske_sedi(j,i): Sediment mask. 1: hard rock, 7: soft sediment, 2: ocean.
Definition: sico_vars_m.F90:46
subroutine, public stereo_inv_ellipsoid(x_val, y_val, A, B, lambda0, phi0, lambda_val, phi_val)
Inverse stereographic projection for an ellipsoidal planet.
subroutine topograd_1(dxi, deta, n_switch)
Calculation of topography gradients on the staggered grid and on the grid points (the latter by secon...
Definition: topograd_m.F90:56
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 ...
real(dp) target_topo_tau
target_topo_tau: Relaxation time for target-topography adjustment
subroutine, public calc_enhance_5()
Computation of the flow enhancement factor. Case ENHMOD==5: minimal anisotropic enhancement factor fo...
real(dp) b
B: Semi-minor axis of the planet.
subroutine, public init_temp_water_age_1_2()
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==2: present-day initial topogr...
real(dp), dimension(0:jmax, 0:imax) zl
zl(j,i): Coordinate z of the lithosphere surface
subroutine, public calc_vxy_static()
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy etc. for static ice.
subroutine, public output4(time, dxi, deta, delta_ts, glac_index, z_sl)
Writing of time-series data of the deep ice cores on file in ASCII format (and optionally in NetCDF f...
Definition: output_m.F90:4721
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:35
integer(i4b), dimension(0:jmax, 0:imax) q_weert
q_weert(j,i): Weertman exponent for the basal pressure
subroutine, public calc_temp_melt()
Computation of the melting temperatures.
real(dp) ea
ea: Abbreviation for exp(aa)
real(dp), parameter pi_180
pi_180: pi divided by 180 (-> deg to rad)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) enth_c
enth_c(kc,j,i): Enthalpy in the upper (kc) ice domain
logical flag_aa_nonzero
flag_aa_nonzero: Flag for the exponential stretch parameter aa. .true.: aa greater than zero (non-equ...
Initial temperature, water content and age.
subroutine, public calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
Computation of the vertical velocity vz for grounded ice.
Definition: calc_vz_m.F90:52
real(dp) aa
aa: Exponential stretch parameter of the non-equidistant vertical grid in the upper (kc) ice domain ...
Reading of several input data.
Definition: read_m.F90:35
Declarations of kind types for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) zs
zs(j,i): Coordinate z of the surface topography
real(dp), dimension(0:jmax, 0:imax) precip_ma_present
precip_ma_present(j,i): Present-day mean annual precipitation rate at the ice surface ...
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) temp_precip_time_min
temp_precip_time_min: Minimum time of the surface-temperature and precipitation data ...
subroutine, public stereo_forw_ellipsoid(lambda_val, phi_val, A, B, lambda0, phi0, x_val, y_val)
Forward stereographic projection for an ellipsoidal planet.
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:201
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
integer(i4b) ndata_temp_precip
ndata_temp_precip: Number of surface-temperature and precipitation data
real(dp), dimension(0:jmax, 0:imax) dzm_dtau
dzm_dtau(j,i): Derivative of zm by tau (time)
real(dp), dimension(0:krmax) zeta_r
zeta_r(kr): Sigma coordinate zeta_r of grid point kr
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) omega_c
omega_c(kc,j,i): Water content in the upper (kc) ice domain
real(dp), dimension(:), allocatable y_core
y_core(n): Coordinate eta (= y) of the prescribed borehole positions
real(dp), dimension(:), allocatable temp_precip_time
temp_precip_time(n): Times of the surface-temperature and precipitation data
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
real(dp), dimension(0:jmax, 0:imax) precip_ma_lgm_anom
precip_ma_lgm_anom(j,i): LGM anomaly (ratio LGM/present) of the mean annual precipitation rate at the...
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.
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.
NetCDF error capturing.
Definition: nc_check_m.F90:35
real(dp) time_target_topo_final
time_target_topo_final: Final time for target-topography adjustment
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(0:jmax, 0:imax) temp_mj_lgm_anom
temp_mj_lgm_anom(j,i): LGM anomaly (difference LGM - present) of the mean summer (northern hemisphere...
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:478
real(dp), dimension(-190:10) c
C(n): Tabulated values for the specific heat of ice.
real(dp), dimension(0:kcmax) eaz_c_quotient
eaz_c_quotient(kc): Abbreviation for (eaz_c(kc)-1.0)/(ea-1.0)
real(dp), dimension(-10000:10000) kei
kei(n): Tabulated values of the kei function (Kelvin function of zero order)
subroutine, public init_temp_water_age_1_3()
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==3: present-day initial topogr...
subroutine boundary(time, dtime, dxi, deta, delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
Main routine of boundary_m: Computation of the surface temperature (must be less than 0 deg C!) and o...
Definition: boundary_m.F90:56
real(dp) time_target_topo_init
time_target_topo_init: Initial time for target-topography adjustment
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) ...
real(dp), dimension(0:jmax, 0:imax) temp_ma_lgm_anom
temp_ma_lgm_anom(j,i): LGM anomaly (difference LGM - present) of the mean annual surface temperature ...
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:3377
subroutine check(status, ch_calling_routine)
NetCDF error capturing.
Definition: nc_check_m.F90:46
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
integer(i4b) ncid_temp_precip
ncid_temp_precip: ID of the NetCDF file containing the surface-temperature and precipitation data as ...
character(len=16), dimension(:), allocatable ch_core
ch_core(n): Names of the prescribed borehole positions
real(dp), dimension(0:jmax, 0:imax) sq_g11_g
sq_g11_g(j,i): Square root of the coefficient g11 of the metric tensor on grid point (i...
real(dp), dimension(0:jmax, 0:imax) as_perp_apl
as_perp_apl(j,i): Applied accumulation-ablation function (SMB)
logical, dimension(0:jmax, 0:imax) flag_calving_front_1
flag_calving_front_1(j,i): Calving front flag. .true.: calving front point (floating ice point with a...
subroutine, public metric()
Main routine of module metric_m: Definition of the components g11 and g22 of the metric tensor of the...
Definition: metric_m.F90:54
subroutine, public calc_enhance_3(time)
Computation of the flow enhancement factor. Case ENHMOD==3: two different values depending on time of...
integer(i4b) grip_time_min
grip_time_min: Minimum time of the data values for the surface temperature anomaly ...
real(dp), dimension(0:jmax, 0:imax) time_lag_asth
time_lag_asth(j,i): Time lag of the relaxing asthenosphere
real(dp), dimension(0:jmax) eta
eta(j): Coordinate eta (= y) of grid point j
real(dp), dimension(:), allocatable x_core
x_core(n): Coordinate xi (= x) of the prescribed borehole positions
real(dp), 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:56
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
real(dp) rho_sw
RHO_SW: Density of sea water.
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
real(dp) rho
RHO: Density of ice.
subroutine topography3(dxi, deta, z_sl, anfdatname)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
integer(i4b), dimension(0:jmax, 0:imax) p_weert
p_weert(j,i): Weertman exponent for the basal shear stress
subroutine, public init_temp_water_age_1_4()
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==4: present-day initial topogr...
subroutine, public init_temp_water_age_1_1()
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==1: present-day initial topogr...
real(dp), dimension(:), allocatable lambda_core
lambda_core(n): Geographical longitude of the prescribed borehole positions
integer(i4b) n_core
n_core: Number of positions to be considered in the time-series file for deep boreholes ...
integer, parameter dp
Double-precision reals.
subroutine topography2(dxi, deta)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
real(dp), dimension(0:jmax, 0:imax) zb
zb(j,i): Coordinate z of the ice base
logical, dimension(0:jmax, 0:imax) flag_grounding_line_2
flag_grounding_line_2(j,i): Grounding line flag. .true.: grounding line point (floating ice point wit...
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
subroutine, public read_target_topo_nc(target_topo_dat_name)
Reading of the target-topography file (in NetCDF format).
Definition: read_m.F90:672
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) 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:53
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...