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