SICOPOLIS V5-dev  Revision 1105
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 q_bm = 0.0_dp
1035 q_tld = 0.0_dp
1036 q_b_tot = 0.0_dp
1037 
1038 p_b_w = 0.0_dp
1039 h_w = 0.0_dp
1040 
1041 call init_temp_water_age_2()
1042 
1043 #if (ENHMOD==1)
1044  call calc_enhance_1()
1045 #elif (ENHMOD==2)
1046  call calc_enhance_2()
1047 #elif (ENHMOD==3)
1048  call calc_enhance_3(time_init)
1049 #elif (ENHMOD==4)
1050  call calc_enhance_4()
1051 #elif (ENHMOD==5)
1052  call calc_enhance_5()
1053 #else
1054  stop ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
1055 #endif
1056 
1057 ! ------ Read initial state from output of previous simulation
1058 
1059 #elif (ANF_DAT==3)
1060 
1061 call topography3(dxi, deta, z_sl, anfdatname)
1062 
1063 call boundary(time_init, dtime, dxi, deta, &
1064  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1065 
1066 where ((maske==0_i2b).or.(maske==3_i2b))
1067  ! grounded or floating ice
1069 elsewhere ! maske==1_i2b or 2_i2b, ice-free land or sea
1070  as_perp_apl = 0.0_dp
1071 end where
1072 
1073 q_b_tot = q_bm + q_tld
1074 
1075 #endif
1076 
1077 !-------- Inner-point flag --------
1078 
1079 flag_inner_point = .true.
1080 
1081 flag_inner_point(0,:) = .false.
1082 flag_inner_point(jmax,:) = .false.
1083 
1084 flag_inner_point(:,0) = .false.
1085 flag_inner_point(:,imax) = .false.
1086 
1087 !-------- Grounding line and calving front flags --------
1088 
1089 flag_grounding_line_1 = .false.
1090 flag_grounding_line_2 = .false.
1091 
1092 flag_calving_front_1 = .false.
1093 flag_calving_front_2 = .false.
1094 
1095 #if (MARGIN==1 || MARGIN==2)
1096 
1097 continue
1098 
1099 #elif (MARGIN==3)
1100 
1101 do i=1, imax-1
1102 do j=1, jmax-1
1103 
1104  if ( (maske(j,i)==0_i2b) & ! grounded ice
1105  .and. &
1106  ( (maske(j,i+1)==3_i2b) & ! with
1107  .or.(maske(j,i-1)==3_i2b) & ! one
1108  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1109  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1110  ) &
1111  flag_grounding_line_1(j,i) = .true.
1112 
1113  if ( (maske(j,i)==3_i2b) & ! floating ice
1114  .and. &
1115  ( (maske(j,i+1)==0_i2b) & ! with
1116  .or.(maske(j,i-1)==0_i2b) & ! one
1117  .or.(maske(j+1,i)==0_i2b) & ! neighbouring
1118  .or.(maske(j-1,i)==0_i2b) ) & ! grounded ice point
1119  ) &
1120  flag_grounding_line_2(j,i) = .true.
1121 
1122  if ( (maske(j,i)==3_i2b) & ! floating ice
1123  .and. &
1124  ( (maske(j,i+1)==2_i2b) & ! with
1125  .or.(maske(j,i-1)==2_i2b) & ! one
1126  .or.(maske(j+1,i)==2_i2b) & ! neighbouring
1127  .or.(maske(j-1,i)==2_i2b) ) & ! ocean point
1128  ) &
1129  flag_calving_front_1(j,i) = .true.
1130 
1131  if ( (maske(j,i)==2_i2b) & ! ocean
1132  .and. &
1133  ( (maske(j,i+1)==3_i2b) & ! with
1134  .or.(maske(j,i-1)==3_i2b) & ! one
1135  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1136  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1137  ) &
1138  flag_calving_front_2(j,i) = .true.
1139 
1140 end do
1141 end do
1142 
1143 do i=1, imax-1
1144 
1145  j=0
1146  if ( (maske(j,i)==2_i2b) & ! ocean
1147  .and. (maske(j+1,i)==3_i2b) & ! with one neighbouring
1148  ) & ! floating ice point
1149  flag_calving_front_2(j,i) = .true.
1150 
1151  j=jmax
1152  if ( (maske(j,i)==2_i2b) & ! ocean
1153  .and. (maske(j-1,i)==3_i2b) & ! with one neighbouring
1154  ) & ! floating ice point
1155  flag_calving_front_2(j,i) = .true.
1156 
1157 end do
1158 
1159 do j=1, jmax-1
1160 
1161  i=0
1162  if ( (maske(j,i)==2_i2b) & ! ocean
1163  .and. (maske(j,i+1)==3_i2b) & ! with one neighbouring
1164  ) & ! floating ice point
1165  flag_calving_front_2(j,i) = .true.
1166 
1167  i=imax
1168  if ( (maske(j,i)==2_i2b) & ! ocean
1169  .and. (maske(j,i-1)==3_i2b) & ! with one neighbouring
1170  ) & ! floating ice point
1171  flag_calving_front_2(j,i) = .true.
1172 
1173 end do
1174 
1175 #else
1176 stop ' >>> sico_init: MARGIN must be either 1, 2 or 3!'
1177 #endif
1178 
1179 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
1180 
1181 #if (GRID==0 || GRID==1)
1182 
1183 do ir=-imax, imax
1184 do jr=-jmax, jmax
1185  dist_dxdy(jr,ir) = sqrt( (real(ir,dp)*dxi)**2 + (real(jr,dp)*deta)**2 )
1186  ! distortion due to stereographic projection not accounted for
1187 end do
1188 end do
1189 
1190 #endif
1191 
1192 !-------- Initial velocities --------
1193 
1194 call calc_temp_melt()
1195 
1196 #if (DYNAMICS==1 || DYNAMICS==2)
1197 
1198 call calc_vxy_b_sia(time, z_sl)
1199 call calc_vxy_sia(dzeta_c, dzeta_t)
1200 
1201 #if (MARGIN==3 || DYNAMICS==2)
1202 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
1203 #endif
1204 
1205 call calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
1206 
1207 #if (MARGIN==3)
1208 call calc_vz_floating(z_sl, dxi, deta, dzeta_c)
1209 #endif
1210 
1211 #elif (DYNAMICS==0)
1212 
1213 call calc_vxy_static()
1214 call calc_vz_static()
1215 
1216 #else
1217  stop ' >>> sico_init: DYNAMICS must be either 0, 1 or 2!'
1218 #endif
1219 
1220 call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
1221 
1222 !-------- Initial enthalpies --------
1223 
1224 #if (CALCMOD==0 || CALCMOD==-1)
1225 
1226 do i=0, imax
1227 do j=0, jmax
1228 
1229  do kc=0, kcmax
1230  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1231  end do
1232 
1233  do kt=0, ktmax
1234  enth_t(kt,j,i) = enth_c(0,j,i)
1235  end do
1236 
1237 end do
1238 end do
1239 
1240 #elif (CALCMOD==1)
1241 
1242 do i=0, imax
1243 do j=0, jmax
1244 
1245  do kc=0, kcmax
1246  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1247  end do
1248 
1249  if ( (maske(j,i) == 0_i2b).and.(n_cts(j,i) == 1_i2b) ) then
1250  do kt=0, ktmax
1251  enth_t(kt,j,i) = enth_fct_temp_omega(temp_t_m(kt,j,i), omega_t(kt,j,i))
1252  end do
1253  else
1254  do kt=0, ktmax
1255  enth_t(kt,j,i) = enth_c(0,j,i)
1256  end do
1257  end if
1258 
1259 end do
1260 end do
1261 
1262 #elif (CALCMOD==2 || CALCMOD==3)
1263 
1264 do i=0, imax
1265 do j=0, jmax
1266 
1267  do kc=0, kcmax
1268  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), omega_c(kc,j,i))
1269  end do
1270 
1271  do kt=0, ktmax
1272  enth_t(kt,j,i) = enth_c(0,j,i)
1273  end do
1274 
1275 end do
1276 end do
1277 
1278 #else
1279 
1280 stop ' >>> sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!'
1281 
1282 #endif
1283 
1284 !-------- Initialize time-series files --------
1285 
1286 ! ------ Time-series file for the ice sheet on the whole
1287 
1288 filename_with_path = trim(outpath)//'/'//trim(runname)//'.ser'
1289 
1290 open(12, iostat=ios, file=trim(filename_with_path), status='new')
1291 
1292 if (ios /= 0) stop ' >>> sico_init: Error when opening the ser file!'
1293 
1294 write(12,1102)
1295 write(12,1103)
1296 
1297  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1298  ' V(m^3) V_g(m^3) V_f(m^3)', &
1299  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1300  ' V_sle(m) V_t(m^3)', &
1301  ' A_t(m^2)',/, &
1302  ' H_max(m) H_t_max(m)', &
1303  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
1304  1103 format('----------------------------------------------------', &
1305  '---------------------------------------')
1306 
1307 ! ------ Time-series file for selected positions ("deep boreholes")
1308 
1309 n_core = 2 ! central dome, position halfway to coast
1310 
1311 allocate(lambda_core(n_core), phi_core(n_core), &
1313 
1314 ch_core(1) = 'P1'
1315 lambda_core(1) = 0.0_dp ! dummy
1316 phi_core(1) = 0.0_dp ! dummy
1317 x_core(1) = 750.0_dp *1.0e+03_dp ! Position of the central dome
1318 y_core(1) = 750.0_dp *1.0e+03_dp ! (750 km, 750 km),
1319  ! conversion km -> m
1320 
1321 ch_core(2) = 'P2'
1322 lambda_core(2) = 0.0_dp ! dummy
1323 phi_core(2) = 0.0_dp ! dummy
1324 x_core(2) = 750.0_dp *1.0e+03_dp ! Position halfway to the coast
1325 y_core(2) = 1125.0_dp *1.0e+03_dp ! (750 km, 1125 km),
1326  ! conversion km -> m
1327 
1328 filename_with_path = trim(outpath)//'/'//trim(runname)//'.core'
1329 
1330 open(14, iostat=ios, file=trim(filename_with_path), status='new')
1331 
1332 if (ios /= 0) stop ' >>> sico_init: Error when opening the core file!'
1333 
1334 if (forcing_flag == 1) then
1335 
1336  write(14,1106)
1337  write(14,1107)
1338 
1339  1106 format(' t(a) D_Ts(C) z_sl(m)',/, &
1340  ' H_P1(m) H_P2(m)',/, &
1341  ' v_P1(m/a) v_P2(m/a)',/, &
1342  ' T_P1(C) T_P2(C)')
1343  1107 format('---------------------------------------')
1344 
1345 end if
1346 
1347 !-------- Output of the initial state --------
1348 
1349 #if (defined(OUTPUT_INIT))
1350 
1351 #if (OUTPUT_INIT==0)
1352  flag_init_output = .false.
1353 #elif (OUTPUT_INIT==1)
1354  flag_init_output = .true.
1355 #else
1356  stop ' >>> sico_init: OUTPUT_INIT must be either 0 or 1!'
1357 #endif
1358 
1359 #else
1360  flag_init_output = .true. ! default for undefined parameter OUTPUT_INIT
1361 #endif
1362 
1363 #if (OUTPUT==1)
1364 
1365 #if (ERGDAT==0)
1366  flag_3d_output = .false.
1367 #elif (ERGDAT==1)
1368  flag_3d_output = .true.
1369 #else
1370  stop ' >>> sico_init: ERGDAT must be either 0 or 1!'
1371 #endif
1372 
1373  if (flag_init_output) &
1374  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1375  flag_3d_output, ndat2d, ndat3d)
1376 
1377 #elif (OUTPUT==2)
1378 
1379 if (time_output(1) <= time_init+eps) then
1380 
1381 #if (ERGDAT==0)
1382  flag_3d_output = .false.
1383 #elif (ERGDAT==1)
1384  flag_3d_output = .true.
1385 #else
1386  stop ' >>> sico_init: ERGDAT must be either 0 or 1!'
1387 #endif
1388 
1389  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1390  flag_3d_output, ndat2d, ndat3d)
1391 
1392 end if
1393 
1394 #elif (OUTPUT==3)
1395 
1396  flag_3d_output = .false.
1397 
1398  if (flag_init_output) &
1399  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1400  flag_3d_output, ndat2d, ndat3d)
1401 
1402 if (time_output(1) <= time_init+eps) then
1403 
1404  flag_3d_output = .true.
1405 
1406  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1407  flag_3d_output, ndat2d, ndat3d)
1408 
1409 end if
1410 
1411 #else
1412  stop ' >>> sico_init: OUTPUT must be either 1, 2 or 3!'
1413 #endif
1414 
1415 if (flag_init_output) then
1416  call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1417  call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1418 end if
1419 
1420 end subroutine sico_init
1421 
1422 !-------------------------------------------------------------------------------
1423 !> Definition of the initial surface and bedrock topography
1424 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1425 !! For ice-free initial topography with relaxed lithosphere
1426 !! (not defined for this domain, thus routine stops execution of SICOPOLIS).
1427 !<------------------------------------------------------------------------------
1428 subroutine topography1(dxi, deta)
1429 
1430 implicit none
1431 
1432 real(dp), intent(out) :: dxi, deta
1433 
1434 ! integer(i4b) :: i, j
1435 ! integer(i4b) :: ios
1436 ! real(dp) :: xi0, eta0
1437 
1438 dxi=0.0_dp; deta=0.0_dp ! dummy values
1439 
1440 write(6,'(/1x,a)') '>>> topography1: ANF_DAT==1 not defined for this domain!'
1441 stop
1442 
1443 end subroutine topography1
1444 
1445 !-------------------------------------------------------------------------------
1446 !> Definition of the initial surface and bedrock topography
1447 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1448 !! For ice-free initial topography with relaxed lithosphere.
1449 !<------------------------------------------------------------------------------
1450 subroutine topography2(dxi, deta)
1451 
1452 #if (GRID==0 || GRID==1)
1454 #endif
1455 
1456  use metric_m
1457  use topograd_m
1458 
1459 implicit none
1460 
1461 real(dp), intent(out) :: dxi, deta
1462 
1463 integer(i4b) :: i, j, n
1464 integer(i4b) :: ios
1465 real(dp) :: xi0, eta0
1466 
1467 real(dp) :: half_width
1468 real(dp), dimension(0:JMAX,0:IMAX) :: zl0_update
1469 
1470 real(dp), parameter :: zl0_ocean = -2000.0_dp
1471 real(dp), parameter :: epss = 0.01_dp
1472 
1473 !-------- Set topography --------
1474 
1475 #if (!defined(LAND_OCEAN_TRANSITION_WIDTH))
1476 
1477 write(6, fmt='(a)') ' >>> topography2: LAND_OCEAN_TRANSITION_WIDTH not defined.'
1478 write(6, fmt='(a)') ' Value 0 assumed by default!'
1479 write(6, fmt='(a)') ' '
1480 
1481 half_width = 0.0_dp
1482 
1483 #else
1484 
1485 half_width = 0.5_dp*abs(land_ocean_transition_width) *1000.0_dp ! km -> m
1486 
1487 #endif
1488 
1489 zl0 = 10.0_dp
1490 
1491 zl0(0,:) = zl0_ocean
1492 zl0(jmax,:) = zl0_ocean
1493 
1494 zl0(:,0) = zl0_ocean
1495 zl0(:,imax) = zl0_ocean
1496 
1497 maske = 1
1498 
1499 maske(0,:) = 2
1500 maske(jmax,:) = 2
1501 
1502 maske(:,0) = 2
1503 maske(:,imax) = 2
1504 
1505 !-------- Further stuff --------
1506 
1507 #if (GRID==0 || GRID==1)
1508 
1509 dxi = dx *1000.0_dp ! km -> m
1510 deta = dx *1000.0_dp ! km -> m
1511 
1512 xi0 = x0 *1000.0_dp ! km -> m
1513 eta0 = y0 *1000.0_dp ! km -> m
1514 
1515 #elif (GRID==2)
1516 
1517 stop ' >>> topography2: GRID==2 not allowed for this application!'
1518 
1519 #endif
1520 
1521 do i=0, imax
1522  xi(i) = xi0 + real(i,dp)*dxi
1523 end do
1524 
1525 do j=0, jmax
1526  eta(j) = eta0 + real(j,dp)*deta
1527 end do
1528 
1529 do i=1, imax-1
1530 do j=1, jmax-1
1531 
1532  if ( (xi(i) >= (1150.0e+03_dp+half_width-epss)) &
1533  .and.(eta(j) >= (625.0e+03_dp+half_width-epss)) &
1534  .and.(eta(j) <= (875.0e+03_dp-half_width+epss)) ) then
1535 
1536  maske(j,i) = 2
1537  zl0(j,i) = zl0_ocean
1538 
1539  end if
1540 
1541 end do
1542 end do
1543 
1544 if (half_width > epss) then
1545 
1546  zl0_update = zl0
1547 
1548  do n=1, 1000
1549 
1550  do i=1, imax-1
1551  do j=1, jmax-1
1552 
1553  if ( ( .not.( (xi(i) >= (1150.0e+03_dp+half_width-epss)) &
1554  .and.(eta(j) >= (625.0e+03_dp+half_width-epss)) &
1555  .and.(eta(j) <= (875.0e+03_dp-half_width+epss)) ) ) &
1556  .and. &
1557  ( (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)) ) ) then
1560 
1561  maske(j,i) = 2
1562 
1563  zl0_update(j,i) = zl0(j,i) &
1564  + 0.1_dp &
1565  * (zl0(j,i+1)-2.0_dp*zl0(j,i)+zl0(j,i-1)) &
1566  + 0.1_dp &
1567  * (zl0(j+1,i)-2.0_dp*zl0(j,i)+zl0(j-1,i))
1568  ! Solving the diffusion equation
1569 
1570  end if
1571 
1572  end do
1573  end do
1574 
1575  zl0 = zl0_update
1576 
1577  end do
1578 
1579 end if
1580 
1581 do i=0, imax
1582 do j=0, jmax
1583 
1584  if (maske(j,i) <= 1) then
1585  maske(j,i) = 1
1586  zs(j,i) = zl0(j,i)
1587  zb(j,i) = zl0(j,i)
1588  zl(j,i) = zl0(j,i)
1589  else ! (maske(j,i) >= 2)
1590  maske(j,i) = 2
1591 #if (MARGIN==1 || MARGIN==2)
1592  zs(j,i) = zl0(j,i)
1593  zb(j,i) = zl0(j,i)
1594 #elif (MARGIN==3)
1595  zs(j,i) = 0.0_dp ! present-day
1596  zb(j,i) = 0.0_dp ! sea level
1597 #endif
1598  zl(j,i) = zl0(j,i)
1599  end if
1600 
1601  zm(j,i) = zb(j,i)
1602  n_cts(j,i) = -1
1603  kc_cts(j,i) = 0
1604 
1605  h_c(j,i) = 0.0_dp
1606  h_t(j,i) = 0.0_dp
1607 
1608  dzs_dtau(j,i) = 0.0_dp
1609  dzm_dtau(j,i) = 0.0_dp
1610  dzb_dtau(j,i) = 0.0_dp
1611  dzl_dtau(j,i) = 0.0_dp
1612  dh_c_dtau(j,i) = 0.0_dp
1613  dh_t_dtau(j,i) = 0.0_dp
1614 
1615 end do
1616 end do
1617 
1618 maske_old = maske
1619 
1620 !-------- Geographic coordinates, metric tensor,
1621 ! gradients of the topography --------
1622 
1623 do i=0, imax
1624 do j=0, jmax
1625 
1626 #if (GRID==0 || GRID==1) /* Stereographic projection */
1627 
1628  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1629  lambda0, phi0, lambda(j,i), phi(j,i))
1630 
1631 #elif (GRID==2) /* Geographic coordinates */
1632 
1633  lambda(j,i) = xi(i)
1634  phi(j,i) = eta(j)
1635 
1636 #endif
1637 
1638 end do
1639 end do
1640 
1641 call metric()
1642 
1643 #if (TOPOGRAD==0)
1644 call topograd_1(dxi, deta, 1)
1645 #elif (TOPOGRAD==1)
1646 call topograd_2(dxi, deta, 1)
1647 #endif
1648 
1649 !-------- Corresponding area of grid points --------
1650 
1651 do i=0, imax
1652 do j=0, jmax
1653  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1654 end do
1655 end do
1656 
1657 end subroutine topography2
1658 
1659 !-------------------------------------------------------------------------------
1660 !> Definition of the initial surface and bedrock topography
1661 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1662 !! For initial topography from previous simulation.
1663 !<------------------------------------------------------------------------------
1664 subroutine topography3(dxi, deta, z_sl, anfdatname)
1665 
1666  use read_m, only : read_erg_nc
1667 
1668 #if (GRID==0 || GRID==1)
1670 #endif
1671 
1672  use metric_m
1673  use topograd_m
1674 
1675 implicit none
1676 
1677 character(len=100), intent(in) :: anfdatname
1678 
1679 real(dp), intent(out) :: dxi, deta, z_sl
1680 
1681 integer(i4b) :: i, j, n
1682 
1683 real(dp) :: half_width
1684 real(dp), dimension(0:JMAX,0:IMAX) :: zl0_update
1685 
1686 real(dp), parameter :: zl0_ocean = -2000.0_dp
1687 real(dp), parameter :: epss = 0.01_dp
1688 
1689 !-------- Read data from time-slice file of previous simulation --------
1690 
1691 call read_erg_nc(z_sl, anfdatname)
1692 
1693 !-------- Set topography of the relaxed bedrock --------
1694 
1695 #if (!defined(LAND_OCEAN_TRANSITION_WIDTH))
1696 
1697 write(6, fmt='(a)') ' >>> topography3: LAND_OCEAN_TRANSITION_WIDTH not defined.'
1698 write(6, fmt='(a)') ' Value 0 assumed by default!'
1699 write(6, fmt='(a)') ' '
1700 
1701 half_width = 0.0_dp
1702 
1703 #else
1704 
1705 half_width = 0.5_dp*abs(land_ocean_transition_width) *1000.0_dp ! km -> m
1706 
1707 #endif
1708 
1709 zl0 = 10.0_dp
1710 
1711 zl0(0,:) = zl0_ocean
1712 zl0(jmax,:) = zl0_ocean
1713 
1714 zl0(:,0) = zl0_ocean
1715 zl0(:,imax) = zl0_ocean
1716 
1717 !-------- Further stuff --------
1718 
1719 #if (GRID==0 || GRID==1)
1720 
1721 dxi = dx *1000.0_dp ! km -> m
1722 deta = dx *1000.0_dp ! km -> m
1723 
1724 #elif (GRID==2)
1725 
1726 stop ' >>> topography3: GRID==2 not allowed for this application!'
1727 
1728 #endif
1729 
1730 do i=1, imax-1
1731 do j=1, jmax-1
1732 
1733  if ( (xi(i) >= (1150.0e+03_dp+half_width-epss)) &
1734  .and.(eta(j) >= (625.0e+03_dp+half_width-epss)) &
1735  .and.(eta(j) <= (875.0e+03_dp-half_width+epss)) ) then
1736 
1737  zl0(j,i) = zl0_ocean
1738 
1739  end if
1740 
1741 end do
1742 end do
1743 
1744 if (half_width > epss) then
1745 
1746  zl0_update = zl0
1747 
1748  do n=1, 1000
1749 
1750  do i=1, imax-1
1751  do j=1, jmax-1
1752 
1753  if ( ( .not.( (xi(i) >= (1150.0e+03_dp+half_width-epss)) &
1754  .and.(eta(j) >= (625.0e+03_dp+half_width-epss)) &
1755  .and.(eta(j) <= (875.0e+03_dp-half_width+epss)) ) ) &
1756  .and. &
1757  ( (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)) ) ) then
1760 
1761  zl0_update(j,i) = zl0(j,i) &
1762  + 0.1_dp &
1763  * (zl0(j,i+1)-2.0_dp*zl0(j,i)+zl0(j,i-1)) &
1764  + 0.1_dp &
1765  * (zl0(j+1,i)-2.0_dp*zl0(j,i)+zl0(j-1,i))
1766  ! Solving the diffusion equation
1767 
1768  end if
1769 
1770  end do
1771  end do
1772 
1773  zl0 = zl0_update
1774 
1775  end do
1776 
1777 end if
1778 
1779 !-------- Geographic coordinates, metric tensor,
1780 ! gradients of the topography --------
1781 
1782 do i=0, imax
1783 do j=0, jmax
1784 
1785 #if (GRID==0 || GRID==1) /* Stereographic projection */
1786 
1787  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1788  lambda0, phi0, lambda(j,i), phi(j,i))
1789 
1790 #elif (GRID==2) /* Geographic coordinates */
1791 
1792  lambda(j,i) = xi(i)
1793  phi(j,i) = eta(j)
1794 
1795 #endif
1796 
1797 end do
1798 end do
1799 
1800 call metric()
1801 
1802 #if (TOPOGRAD==0)
1803 call topograd_1(dxi, deta, 1)
1804 #elif (TOPOGRAD==1)
1805 call topograd_2(dxi, deta, 1)
1806 #endif
1807 
1808 !-------- Corresponding area of grid points --------
1809 
1810 do i=0, imax
1811 do j=0, jmax
1812  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1813 end do
1814 end do
1815 
1816 end subroutine topography3
1817 
1818 !-------------------------------------------------------------------------------
1819 
1820 end module sico_init_m
1821 !
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:813
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:4721
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) 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:3377
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...