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