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