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