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