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