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, read_kei
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 integer(i4b) :: ndata_insol
95 real(dp) :: dtime0, dtime_temp0, dtime_wss0, dtime_out0, dtime_ser0
96 real(dp) :: time_init0, time_end0, time_output0(100)
97 real(dp) :: time_chasm_init0, time_chasm_end0
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 character(len= 8) :: ch_imax
111 character(len=128) :: fmt4
112 
113 write(ch_imax, fmt='(i8)') imax
114 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
115 
116 write(unit=6, fmt='(a)') ' '
117 write(unit=6, fmt='(a)') ' -------- sico_init --------'
118 write(unit=6, fmt='(a)') ' '
119 
120 !-------- Name of the computational domain --------
121 
122 #if (defined(ANT))
123 ch_domain_long = 'Antarctica'
124 ch_domain_short = 'ant'
125 
126 #elif (defined(ASF))
127 ch_domain_long = 'Austfonna'
128 ch_domain_short = 'asf'
129 
130 #elif (defined(EMTP2SGE))
131 ch_domain_long = 'EISMINT Phase 2 Simplified Geometry Experiment'
132 ch_domain_short = 'emtp2sge'
133 
134 #elif (defined(GRL))
135 ch_domain_long = 'Greenland'
136 ch_domain_short = 'grl'
137 
138 #elif (defined(NHEM))
139 ch_domain_long = 'Northern hemisphere'
140 ch_domain_short = 'nhem'
141 
142 #elif (defined(SCAND))
143 ch_domain_long = 'Scandinavia and Eurasia'
144 ch_domain_short = 'scand'
145 
146 #elif (defined(TIBET))
147 ch_domain_long = 'Tibet'
148 ch_domain_short = 'tibet'
149 
150 #elif (defined(NMARS))
151 ch_domain_long = 'North polar cap of Mars'
152 ch_domain_short = 'nmars'
153 
154 #elif (defined(SMARS))
155 ch_domain_long = 'South polar cap of Mars'
156 ch_domain_short = 'smars'
157 
158 #elif (defined(XYZ))
159 ch_domain_long = 'XYZ'
160 ch_domain_short = 'xyz'
161 #if (defined(HEINO))
162 ch_domain_long = trim(ch_domain_long)//'/ISMIP HEINO'
163 #endif
164 
165 #else
166 
167 errormsg = ' >>> sico_init: No valid domain specified!'
168 call error(errormsg)
169 
170 #endif
171 
172 !-------- Some initial values --------
173 
174 n_output = 0
175 
176 dtime = 0.0_dp
177 dtime_temp = 0.0_dp
178 dtime_wss = 0.0_dp
179 dtime_out = 0.0_dp
180 dtime_ser = 0.0_dp
181 
182 time = 0.0_dp
183 time_init = 0.0_dp
184 time_end = 0.0_dp
185 time_output = 0.0_dp
186 
187 !-------- Initialisation of the Library of Iterative Solvers Lis,
188 ! if required --------
189 
190 #if (CALCTHK==3 || CALCTHK==6 || MARGIN==3 || DYNAMICS==2)
191  call lis_initialize(ierr)
192 #endif
193 
194 !-------- Read physical parameters --------
195 
196 call read_phys_para()
197 
198 call ice_mat_eqs_pars(rf, r_t, kappa, c, -190, 10, rho_i, rho_c, kappa_c, c_c)
199 
200 ! ------ Density of ice-dust mixture
201 
202 rho = (1.0_dp-frac_dust)*rho_i + frac_dust*rho_c
203 
204 rho_inv = 1.0_dp/rho
205 
206 ! ------ Some auxiliary quantities required for the enthalpy method
207 
208 call calc_c_int_table(c, -190, 10, l)
210 
211 !-------- Compatibility check of the SICOPOLIS version with the header file
212 
213 if ( trim(version) /= trim(sico_version) ) then
214  errormsg = ' >>> sico_init: ' &
215  //'SICOPOLIS version not compatible with header file!'
216  call error(errormsg)
217 end if
218 
219 !-------- Check whether the dynamics and thermodynamics modes are defined
220 
221 #if (!defined(DYNAMICS))
222 errormsg = ' >>> sico_init: DYNAMICS not defined in the header file!'
223 call error(errormsg)
224 #endif
225 
226 #if (!defined(CALCMOD))
227 errormsg = ' >>> sico_init: CALCMOD not defined in the header file!'
228 call error(errormsg)
229 #endif
230 
231 #if (defined(ENTHMOD))
232 errormsg = ' >>> sico_init: ENTHMOD must not be defined any more.' &
233  // end_of_line &
234  //' Please update your header file!'
235 call error(errormsg)
236 #endif
237 
238 !-------- Compatibility check of the horizontal resolution with the
239 ! number of grid points --------
240 
241 #if (GRID==0 || GRID==1)
242 
243 if (approx_equal(dx, 20.0_dp, eps_sp_dp)) then
244 
245  if ((imax /= 120).or.(jmax /= 120)) then
246  errormsg = ' >>> sico_init: Wrong values for IMAX and JMAX!'
247  call error(errormsg)
248  end if
249 
250 else if (approx_equal(dx, 10.0_dp, eps_sp_dp)) then
251 
252  if ((imax /= 240).or.(jmax /= 240)) then
253  errormsg = ' >>> sico_init: Wrong values for IMAX and JMAX!'
254  call error(errormsg)
255  end if
256 
257 else if (approx_equal(dx, 5.0_dp, eps_sp_dp)) then
258 
259  if ((imax /= 480).or.(jmax /= 480)) then
260  errormsg = ' >>> sico_init: Wrong values for IMAX and JMAX!'
261  call error(errormsg)
262  end if
263 
264 else
265 
266  errormsg = ' >>> sico_init: Wrong value for DX!'
267  call error(errormsg)
268 
269 end if
270 
271 #elif (GRID==2)
272 
273 errormsg = ' >>> sico_init: GRID==2 not allowed for smars application!'
274 call error(errormsg)
275 
276 #endif
277 
278 !-------- Compatibility check of the thermodynamics mode
279 ! (cold vs. polythermal vs. enthalpy method)
280 ! and the number of grid points in the lower (kt) ice domain --------
281 
282 #if (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1)
283 
284 if (ktmax > 2) then
285  write(6, fmt='(a)') ' >>> sico_init: For options CALCMOD==0, 2, 3 or -1,'
286  write(6, fmt='(a)') ' the separate kt domain is redundant.'
287  write(6, fmt='(a)') ' Therefore, consider setting KTMAX to 2.'
288  write(6, fmt='(a)') ' '
289 end if
290 
291 #endif
292 
293 !-------- Compatibility check of discretization schemes for the horizontal and
294 ! vertical advection terms in the temperature and age equations --------
295 
296 #if (ADV_HOR==1)
297 errormsg = ' >>> sico_init: ' &
298  //'Option ADV_HOR==1 (central differences) not defined!'
299 call error(errormsg)
300 #endif
301 
302 !-------- Check whether for the shallow shelf
303 ! or shelfy stream approximation
304 ! the chosen grid is Cartesian coordinates
305 ! without distortion correction (GRID==0) --------
306 
307 #if ((MARGIN==3 && DYNAMICS==1) || DYNAMICS==2) /* requires SSA or SStA */
308 #if (GRID != 0)
309 errormsg = ' >>> sico_init: Distortion correction not implemented' &
310  // end_of_line &
311  //' for the shallow shelf approximation (SSA)' &
312  // end_of_line &
313  //' or the shelfy stream approximation (SStA)' &
314  // end_of_line &
315  //' -> GRID==0 required!'
316 call error(errormsg)
317 #endif
318 #endif
319 
320 !-------- Setting of forcing flag --------
321 
322 forcing_flag = 1 ! forcing by delta_ts
323 
324 !-------- Initialization of numerical time steps --------
325 
326 dtime0 = dtime0
327 dtime_temp0 = dtime_temp0
328 #if (REBOUND==2)
329 dtime_wss0 = dtime_wss0
330 #endif
331 
332 !-------- Further initializations --------
333 
334 dzeta_c = 1.0_dp/real(kcmax,dp)
335 dzeta_t = 1.0_dp/real(ktmax,dp)
336 dzeta_r = 1.0_dp/real(krmax,dp)
337 
338 ndat2d = 1
339 ndat3d = 1
340 
341 !-------- General abbreviations --------
342 
343 ! ------ kc domain
344 
345 if (deform >= eps) then
346 
347  flag_aa_nonzero = .true. ! non-equidistant grid
348 
349  aa = deform
350  ea = exp(aa)
351 
352  kc=0
353  zeta_c(kc) = 0.0_dp
354  eaz_c(kc) = 1.0_dp
355  eaz_c_quotient(kc) = 0.0_dp
356 
357  do kc=1, kcmax-1
358  zeta_c(kc) = kc*dzeta_c
359  eaz_c(kc) = exp(aa*zeta_c(kc))
360  eaz_c_quotient(kc) = (eaz_c(kc)-1.0_dp)/(ea-1.0_dp)
361  end do
362 
363  kc=kcmax
364  zeta_c(kc) = 1.0_dp
365  eaz_c(kc) = exp(aa)
366  eaz_c_quotient(kc) = 1.0_dp
367 
368 else
369 
370  flag_aa_nonzero = .false. ! equidistant grid
371 
372  aa = 0.0_dp
373  ea = 1.0_dp
374 
375  kc=0
376  zeta_c(kc) = 0.0_dp
377  eaz_c(kc) = 1.0_dp
378  eaz_c_quotient(kc) = 0.0_dp
379 
380  do kc=1, kcmax-1
381  zeta_c(kc) = kc*dzeta_c
382  eaz_c(kc) = 1.0_dp
383  eaz_c_quotient(kc) = zeta_c(kc)
384  end do
385 
386  kc=kcmax
387  zeta_c(kc) = 1.0_dp
388  eaz_c(kc) = 1.0_dp
389  eaz_c_quotient(kc) = 1.0_dp
390 
391 end if
392 
393 ! ------ kt domain
394 
395 kt=0
396 zeta_t(kt) = 0.0_dp
397 
398 do kt=1, ktmax-1
399  zeta_t(kt) = kt*dzeta_t
400 end do
401 
402 kt=ktmax
403 zeta_t(kt) = 1.0_dp
404 
405 ! ------ kr domain
406 
407 kr=0
408 zeta_r(kr) = 0.0_dp
409 
410 do kr=1, krmax-1
411  zeta_r(kr) = kr*dzeta_r
412 end do
413 
414 kr=krmax
415 zeta_r(kr) = 1.0_dp
416 
417 !-------- Reshaping of a 2-d array (with indices i, j)
418 ! to a vector (with index n) --------
419 
420 n=1
421 
422 do i=0, imax
423 do j=0, jmax
424  n2i(n) = i
425  n2j(n) = j
426  ij2n(j,i) = n
427  n=n+1
428 end do
429 end do
430 
431 !-------- Specification of current simulation --------
432 
433 runname = runname
434 anfdatname = anfdatname
435 
436 #if (defined(YEAR_ZERO))
438 #else
439 year_zero = 2000.0_dp ! default value 2000 CE
440 #endif
441 
442 time_init0 = time_init0
443 time_end0 = time_end0
444 dtime_ser0 = dtime_ser0
445 
446 #if (CHASM==2)
447 time_chasm_init0 = time_chasm_init0
448 time_chasm_end0 = time_chasm_end0
449 #endif
450 
451 #if (OUTPUT==1 || OUTPUT==3)
452 dtime_out0 = dtime_out0
453 #endif
454 
455 #if (OUTPUT==2 || OUTPUT==3)
456 n_output = n_output
457 time_output0( 1) = time_out0_01
458 time_output0( 2) = time_out0_02
459 time_output0( 3) = time_out0_03
460 time_output0( 4) = time_out0_04
461 time_output0( 5) = time_out0_05
462 time_output0( 6) = time_out0_06
463 time_output0( 7) = time_out0_07
464 time_output0( 8) = time_out0_08
465 time_output0( 9) = time_out0_09
466 time_output0(10) = time_out0_10
467 time_output0(11) = time_out0_11
468 time_output0(12) = time_out0_12
469 time_output0(13) = time_out0_13
470 time_output0(14) = time_out0_14
471 time_output0(15) = time_out0_15
472 time_output0(16) = time_out0_16
473 time_output0(17) = time_out0_17
474 time_output0(18) = time_out0_18
475 time_output0(19) = time_out0_19
476 time_output0(20) = time_out0_20
477 #endif
478 
479 !-------- Write log file --------
480 
481 shell_command = 'if [ ! -d'
482 shell_command = trim(shell_command)//' '//outpath
483 shell_command = trim(shell_command)//' '//'] ; then mkdir'
484 shell_command = trim(shell_command)//' '//outpath
485 shell_command = trim(shell_command)//' '//'; fi'
486 call system(trim(shell_command))
487  ! Check whether directory OUTPATH exists. If not, it is created.
488 
489 filename_with_path = trim(outpath)//'/'//trim(runname)//'.log'
490 
491 open(10, iostat=ios, file=trim(filename_with_path), status='new')
492 
493 if (ios /= 0) then
494  errormsg = ' >>> sico_init: Error when opening the log file!'
495  call error(errormsg)
496 end if
497 
498 write(10, fmt=trim(fmt1)) 'Computational domain:'
499 write(10, fmt=trim(fmt1)) trim(ch_domain_long)
500 write(10, fmt=trim(fmt1)) ' '
501 
502 write(10, fmt=trim(fmt2a)) 'GRID = ', grid
503 write(10, fmt=trim(fmt1)) ' '
504 
505 write(10, fmt=trim(fmt2)) 'imax =', imax
506 write(10, fmt=trim(fmt2)) 'jmax =', jmax
507 write(10, fmt=trim(fmt2)) 'kcmax =', kcmax
508 write(10, fmt=trim(fmt2)) 'ktmax =', ktmax
509 write(10, fmt=trim(fmt2)) 'krmax =', krmax
510 write(10, fmt=trim(fmt1)) ' '
511 
512 write(10, fmt=trim(fmt3)) 'a =', aa
513 write(10, fmt=trim(fmt1)) ' '
514 
515 write(10, fmt=trim(fmt3)) 'x0 =', x0
516 write(10, fmt=trim(fmt3)) 'y0 =', y0
517 write(10, fmt=trim(fmt3)) 'dx =', dx
518 write(10, fmt=trim(fmt1)) ' '
519 
520 write(10, fmt=trim(fmt3)) 'year_zero =', year_zero
521 write(10, fmt=trim(fmt3)) 'time_init =', time_init0
522 write(10, fmt=trim(fmt3)) 'time_end =', time_end0
523 write(10, fmt=trim(fmt3)) 'dtime =', dtime0
524 write(10, fmt=trim(fmt3)) 'dtime_temp =', dtime_temp0
525 #if (REBOUND==2)
526 write(10, fmt=trim(fmt3)) 'dtime_wss =', dtime_wss0
527 #endif
528 write(10, fmt=trim(fmt1)) ' '
529 
530 write(10, fmt=trim(fmt2)) 'ANF_DAT =', anf_dat
531 write(10, fmt=trim(fmt1)) 'zs_present file = '//zs_present_file
532 #if (ANF_DAT==1)
533 write(10, fmt=trim(fmt1)) 'zl_present file = '//zl_present_file
534 #endif
535 write(10, fmt=trim(fmt1)) 'zl0 file = '//zl0_file
536 write(10, fmt=trim(fmt1)) 'mask_present file = '//mask_present_file
537 #if (ANF_DAT==1 && defined(TEMP_INIT))
538 write(10, fmt=trim(fmt2)) 'TEMP_INIT =', temp_init
539 #endif
540 #if (ANF_DAT==3 || (ANF_DAT==1 && TEMP_INIT==5))
541 write(10, fmt=trim(fmt1)) 'Initial-value file = '//anfdatname
542 write(10, fmt=trim(fmt1)) 'Path to initial-value file = '//anfdatpath
543 #endif
544 write(10, fmt=trim(fmt1)) ' '
545 
546 write(10, fmt=trim(fmt1)) 'Physical-parameter file = '//phys_para_file
547 write(10, fmt=trim(fmt1)) ' '
548 
549 #if (defined(THK_EVOL))
550 write(10, fmt=trim(fmt2a)) 'THK_EVOL = ', thk_evol
551 #else
552 errormsg = ' >>> sico_init: Define THK_EVOL in header file!'
553 call error(errormsg)
554 #endif
555 #if (defined(CALCTHK))
556 write(10, fmt=trim(fmt2a)) 'CALCTHK = ', calcthk
557 #else
558 errormsg = ' >>> sico_init: Define CALCTHK in header file!'
559 call error(errormsg)
560 #endif
561 
562 #if (defined(H_ISOL_MAX))
563 write(10, fmt=trim(fmt3)) 'H_isol_max =', h_isol_max
564 #endif
565 
566 #if (CALCTHK==2 || CALCTHK==3 || CALCTHK==5 || CALCTHK==6)
567 write(10, fmt=trim(fmt3)) 'ovi_weight =', ovi_weight
568 #if (CALCTHK==2 || CALCTHK==5)
569 write(10, fmt=trim(fmt3)) 'omega_sor =', omega_sor
570 #if (ITER_MAX_SOR>0)
571 write(10, fmt=trim(fmt2a)) 'iter_max_sor = ', iter_max_sor
572 #endif
573 #endif
574 #endif
575 
576 write(10, fmt=trim(fmt1)) ' '
577 
578 #if (TSURFACE==1 || TSURFACE==6)
579 write(10, fmt=trim(fmt3)) 'delta_ts0 =', delta_ts0
580 #elif (TSURFACE==3)
581 write(10, fmt=trim(fmt3)) 'sine_amplit =', sine_amplit
582 write(10, fmt=trim(fmt3)) 'sine_period =', sine_period
583 #endif
584 #if (TSURFACE==1 || TSURFACE==2 || TSURFACE==3)
585 write(10, fmt=trim(fmt3)) 'temp0_ma_90S =', temp0_ma_90s
586 #endif
587 #if (TSURFACE==1 || TSURFACE==2 || TSURFACE==3 || TSURFACE==4 || TSURFACE==5)
588 write(10, fmt=trim(fmt3)) 'c_ma =', c_ma
589 write(10, fmt=trim(fmt3)) 'gamma_ma =', gamma_ma
590 #endif
591 #if (TSURFACE==5)
592 write(10, fmt=trim(fmt1)) 'Insolation file = '//insol_ma_90s_file
593 #endif
594 #if (TSURFACE==4 || TSURFACE==5 || TSURFACE==6)
595 write(10, fmt=trim(fmt3)) 'albedo =', albedo
596 #endif
597 
598 #if (ACC_PRESENT==1)
599 write(10, fmt=trim(fmt3)) 'acc_present_val =', acc_present_val
600 #endif
601 
602 #if (ACCSURFACE==1 || ACCSURFACE==2)
603 write(10, fmt=trim(fmt3)) 'gamma_s =', gamma_s
604 #endif
605 
606 #if (ABLSURFACE==1 || ABLSURFACE==2)
607 write(10, fmt=trim(fmt3)) 'eld_0 =', eld_0
608 write(10, fmt=trim(fmt3)) 'g_mb_0 =', g_0
609 write(10, fmt=trim(fmt3)) 'gamma_eld =', gamma_eld
610 write(10, fmt=trim(fmt3)) 'gamma_g =', gamma_g
611 #endif
612 
613 write(10, fmt=trim(fmt2a)) 'SEA_LEVEL = ', sea_level
614 #if (SEA_LEVEL==1)
615 write(10, fmt=trim(fmt3)) 'z_sl0 =', z_sl0
616 #elif (SEA_LEVEL==3)
617 write(10, fmt=trim(fmt1)) 'sea-level file = '//sea_level_file
618 #endif
619 write(10, fmt=trim(fmt1)) ' '
620 
621 #if (MARGIN==2)
622 #if (MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3)
623 write(10, fmt=trim(fmt3)) 'z_mar =', z_mar
624 write(10, fmt=trim(fmt1)) ' '
625 #elif (MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \
626  || marine_ice_calving==6 || marine_ice_calving==7)
627 write(10, fmt=trim(fmt3)) 'fact_z_mar =', fact_z_mar
628 write(10, fmt=trim(fmt1)) ' '
629 #elif (MARINE_ICE_FORMATION==2 && MARINE_ICE_CALVING==9)
630 write(10, fmt=trim(fmt3)) 'calv_uw_coeff =', calv_uw_coeff
631 write(10, fmt=trim(fmt3)) 'r1_calv_uw =', r1_calv_uw
632 write(10, fmt=trim(fmt3)) 'r2_calv_uw =', r2_calv_uw
633 write(10, fmt=trim(fmt1)) ' '
634 #endif
635 #elif (MARGIN==3)
636 #if (ICE_SHELF_CALVING==2)
637 write(10, fmt=trim(fmt3)) 'H_calv =', h_calv
638 write(10, fmt=trim(fmt1)) ' '
639 #endif
640 #endif
641 
642 #if (defined(BASAL_HYDROLOGY))
643 write(10, fmt=trim(fmt2a)) 'BASAL_HYDROLOGY = ', basal_hydrology
644 #endif
645 write(10, fmt=trim(fmt2a)) 'SLIDE_LAW = ', slide_law
646 write(10, fmt=trim(fmt3)) 'c_slide =', c_slide
647 write(10, fmt=trim(fmt3)) 'gamma_slide =', gamma_slide
648 write(10, fmt=trim(fmt2a)) 'p_weert = ', p_weert
649 write(10, fmt=trim(fmt2a)) 'q_weert = ', q_weert
650 #if (defined(TIME_RAMP_UP_SLIDE))
651 write(10, fmt=trim(fmt3)) 'time_ramp_up_slide =', time_ramp_up_slide
652 #endif
653 #if (SLIDE_LAW==2 || SLIDE_LAW==3)
654 write(10, fmt=trim(fmt3)) 'red_pres_limit_fact =', red_pres_limit_fact
655 #endif
656 #if (BASAL_HYDROLOGY==1 \
657  && defined(hydro_slide_sat_fct) \
658  && defined(c_hw_slide) && defined(hw0_slide))
659 write(10, fmt=trim(fmt2a)) 'HYDRO_SLIDE_SAT_FCT = ', hydro_slide_sat_fct
660 write(10, fmt=trim(fmt3)) 'c_Hw_slide =', c_hw_slide
661 write(10, fmt=trim(fmt3)) 'Hw0_slide =', hw0_slide
662 #endif
663 write(10, fmt=trim(fmt1)) ' '
664 
665 write(10, fmt=trim(fmt2a)) 'Q_GEO_MOD = ', q_geo_mod
666 #if (Q_GEO_MOD==1)
667 write(10, fmt=trim(fmt3)) 'q_geo =', q_geo
668 #elif (Q_GEO_MOD==2)
669 write(10, fmt=trim(fmt1)) 'q_geo file = '//q_geo_file
670 #endif
671 write(10, fmt=trim(fmt1)) ' '
672 
673 write(10, fmt=trim(fmt2)) 'REBOUND =', rebound
674 #if (REBOUND==1)
675 write(10, fmt=trim(fmt3)) 'frac_llra =', frac_llra
676 #endif
677 #if (REBOUND==1 || REBOUND==2)
678 write(10, fmt=trim(fmt2)) 'TIME_LAG_MOD =', time_lag_mod
679 #if (TIME_LAG_MOD==1)
680 write(10, fmt=trim(fmt3)) 'time_lag =', time_lag
681 #elif (TIME_LAG_MOD==2)
682 write(10, fmt=trim(fmt1)) 'time_lag_file = '//time_lag_file
683 #else
684 errormsg = ' >>> sico_init: TIME_LAG_MOD must be either 1 or 2!'
685 call error(errormsg)
686 #endif
687 #endif
688 #if (REBOUND==2)
689 write(10, fmt=trim(fmt2)) 'FLEX_RIG_MOD =', flex_rig_mod
690 #if (FLEX_RIG_MOD==1)
691 write(10, fmt=trim(fmt3)) 'flex_rig =', flex_rig
692 #elif (FLEX_RIG_MOD==2)
693 write(10, fmt=trim(fmt1)) 'flex_rig_file = '//flex_rig_file
694 #else
695 errormsg = ' >>> sico_init: FLEX_RIG_MOD must be either 1 or 2!'
696 call error(errormsg)
697 #endif
698 #endif
699 write(10, fmt=trim(fmt2)) 'Q_LITHO =', q_litho
700 write(10, fmt=trim(fmt1)) ' '
701 
702 #if (FLOW_LAW==2)
703 write(10, fmt=trim(fmt3)) 'gr_size =', gr_size
704 write(10, fmt=trim(fmt1)) ' '
705 #endif
706 #if (FIN_VISC==2)
707 write(10, fmt=trim(fmt3)) 'sigma_res =', sigma_res
708 write(10, fmt=trim(fmt1)) ' '
709 #endif
710 write(10, fmt=trim(fmt3)) 'frac_dust =', frac_dust
711 write(10, fmt=trim(fmt1)) ' '
712 
713 write(10, fmt=trim(fmt2)) 'ENHMOD =', enhmod
714 #if (ENHMOD==1 || ENHMOD==2 || ENHMOD==3)
715 write(10, fmt=trim(fmt3)) 'enh_fact =', enh_fact
716 #endif
717 #if (ENHMOD==2 || ENHMOD==3)
718 write(10, fmt=trim(fmt3)) 'enh_intg =', enh_intg
719 #endif
720 #if (ENHMOD==2)
721 write(10, fmt=trim(fmt3)) 'age_trans =', age_trans_0
722 #endif
723 #if (ENHMOD==3)
724 write(10, fmt=trim(fmt3)) 'date_trans1 =', date_trans1_0
725 write(10, fmt=trim(fmt3)) 'date_trans2 =', date_trans2_0
726 write(10, fmt=trim(fmt3)) 'date_trans3 =', date_trans3_0
727 #endif
728 #if (ENHMOD==4 || ENHMOD==5)
729 write(10, fmt=trim(fmt3)) 'enh_compr =', enh_compr
730 write(10, fmt=trim(fmt3)) 'enh_shear =', enh_shear
731 #endif
732 #if ((ENHMOD==1 || ENHMOD==2 || ENHMOD==3 || ENHMOD==4) && MARGIN==3)
733 write(10, fmt=trim(fmt3)) 'enh_shelf =', enh_shelf
734 #endif
735 write(10, fmt=trim(fmt1)) ' '
736 
737 write(10, fmt=trim(fmt2a)) 'CHASM =', chasm
738 #if (CHASM==2)
739 write(10, fmt=trim(fmt1)) 'mask_chasm file = '//mask_chasm_file
740 write(10, fmt=trim(fmt3)) 'time_chasm_init =', time_chasm_init0
741 write(10, fmt=trim(fmt3)) 'time_chasm_end =', time_chasm_end0
742 write(10, fmt=trim(fmt3)) 'q_geo_chasm =', q_geo_chasm
743 write(10, fmt=trim(fmt3)) 'erosion_chasm =', erosion_chasm
744 write(10, fmt=trim(fmt1)) ' '
745 #endif
746 
747 write(10, fmt=trim(fmt2a)) 'DYNAMICS = ', dynamics
748 #if ((DYNAMICS==1 && MARGIN==3) || DYNAMICS==2)
749 #if (defined(LIS_OPTS))
750 write(10, fmt=trim(fmt1)) 'lis_opts = '//lis_opts
751 #endif
752 #if (defined(TOL_ITER_SSA))
753 write(10, fmt=trim(fmt3)) 'tol_iter_ssa =', tol_iter_ssa
754 #endif
755 #if (defined(N_ITER_SSA))
756 write(10, fmt=trim(fmt2a)) 'n_iter_ssa = ', n_iter_ssa
757 #endif
758 #if (defined(ITER_INIT_SSA))
759 write(10, fmt=trim(fmt2a)) 'iter_init_ssa = ', iter_init_ssa
760 #endif
761 #if (defined(VISC_INIT_SSA))
762 write(10, fmt=trim(fmt3)) 'visc_init_ssa =', visc_init_ssa
763 #endif
764 #if (defined(N_VISC_SMOOTH))
765 write(10, fmt=trim(fmt2a)) 'n_visc_smooth = ', n_visc_smooth
766 #endif
767 #if (defined(VISC_SMOOTH_DIFF))
768 write(10, fmt=trim(fmt3)) 'visc_smooth_diff =', visc_smooth_diff
769 #endif
770 #if (defined(RELAX_FACT_SSA))
771 write(10, fmt=trim(fmt3)) 'relax_fact_ssa =', relax_fact_ssa
772 #endif
773 #endif
774 #if (DYNAMICS==2 && defined(RATIO_SL_THRESH))
775 write(10, fmt=trim(fmt3)) 'ratio_sl_thresh =', ratio_sl_thresh
776 #endif
777 #if (DYNAMICS==2 && defined(SSTA_SIA_WEIGH_FCT))
778 write(10, fmt=trim(fmt2a)) 'SSTA_SIA_WEIGH_FCT = ', ssta_sia_weigh_fct
779 #endif
780 write(10, fmt=trim(fmt1)) ' '
781 
782 write(10, fmt=trim(fmt2)) 'CALCMOD =', calcmod
783 #if (CALCMOD==-1 && defined(TEMP_CONST))
784 write(10, fmt=trim(fmt3)) 'TEMP_CONST =', temp_const
785 #endif
786 #if (CALCMOD==-1 && defined(AGE_CONST))
787 write(10, fmt=trim(fmt3)) 'AGE_CONST =', age_const
788 #endif
789 #if (CALCMOD==1 && defined(CTS_MELTING_FREEZING))
790 write(10, fmt=trim(fmt2)) 'CTS_MELTING_FREEZING =', cts_melting_freezing
791 #endif
792 write(10, fmt=trim(fmt2)) 'FLOW_LAW =', flow_law
793 write(10, fmt=trim(fmt2)) 'FIN_VISC =', fin_visc
794 write(10, fmt=trim(fmt2)) 'MARGIN =', margin
795 #if (MARGIN==2)
796 write(10, fmt=trim(fmt2)) 'MARINE_ICE_FORMATION =', marine_ice_formation
797 write(10, fmt=trim(fmt2)) 'MARINE_ICE_CALVING =', marine_ice_calving
798 #elif (MARGIN==3)
799 write(10, fmt=trim(fmt2)) 'ICE_SHELF_CALVING =', ice_shelf_calving
800 #endif
801 write(10, fmt=trim(fmt2)) 'ADV_HOR =', adv_hor
802 write(10, fmt=trim(fmt2)) 'ADV_VERT =', adv_vert
803 write(10, fmt=trim(fmt2)) 'TOPOGRAD =', topograd
804 write(10, fmt=trim(fmt2)) 'TSURFACE =', tsurface
805 write(10, fmt=trim(fmt2)) 'ACC_UNIT =', acc_unit
806 write(10, fmt=trim(fmt2)) 'ACC_PRESENT=', acc_present
807 write(10, fmt=trim(fmt2)) 'ACCSURFACE =', accsurface
808 write(10, fmt=trim(fmt2)) 'ABLSURFACE =', ablsurface
809 #if (defined(MB_ACCOUNT))
810 write(10, fmt=trim(fmt2a)) 'MB_ACCOUNT = ', mb_account
811 #endif
812 write(10, fmt=trim(fmt1)) ' '
813 
814 write(10, fmt=trim(fmt3)) 'numdiff_H_t =', numdiff_h_t
815 write(10, fmt=trim(fmt3)) 'tau_cts =', tau_cts
816 write(10, fmt=trim(fmt3)) 'vh_max =', vh_max
817 write(10, fmt=trim(fmt3)) 'hd_min =', hd_min
818 write(10, fmt=trim(fmt3)) 'hd_max =', hd_max
819 #if (defined(VISC_MIN) && defined(VISC_MAX))
820 write(10, fmt=trim(fmt3)) 'visc_min =', visc_min
821 write(10, fmt=trim(fmt3)) 'visc_max =', visc_max
822 #endif
823 write(10, fmt=trim(fmt3)) 'qbm_min =', qbm_min
824 write(10, fmt=trim(fmt3)) 'qbm_max =', qbm_max
825 write(10, fmt=trim(fmt3)) 'age_min =', age_min
826 write(10, fmt=trim(fmt3)) 'age_max =', age_max
827 write(10, fmt=trim(fmt3)) 'mean_accum =', mean_accum
828 #if (ADV_VERT==1)
829 write(10, fmt=trim(fmt3)) 'age_diff =', agediff
830 #endif
831 write(10, fmt=trim(fmt1)) ' '
832 
833 #if (defined(OUT_TIMES))
834 write(10, fmt=trim(fmt2a)) 'OUT_TIMES = ', out_times
835 #endif
836 #if (defined(OUTPUT_INIT))
837 write(10, fmt=trim(fmt2a)) 'OUTPUT_INIT = ', output_init
838 #endif
839 write(10, fmt=trim(fmt2a)) 'OUTPUT = ', output
840 #if (OUTPUT==1 || OUTPUT==3)
841 write(10, fmt=trim(fmt3)) 'dtime_out =' , dtime_out0
842 #endif
843 write(10, fmt=trim(fmt3)) 'dtime_ser =' , dtime_ser0
844 #if (OUTPUT==1 || OUTPUT==2)
845 write(10, fmt=trim(fmt2a)) 'ERGDAT = ', ergdat
846 #endif
847 write(10, fmt=trim(fmt2a)) 'NETCDF = ', netcdf
848 #if (OUTPUT==2 || OUTPUT==3)
849 write(10, fmt=trim(fmt2a)) 'n_output = ', n_output
850 do n=1, n_output
851  write(10, fmt=trim(fmt3)) 'time_output =' , time_output0(n)
852 end do
853 #endif
854 write(10, fmt=trim(fmt1)) ' '
855 
856 write(10, fmt=trim(fmt1)) 'Program version and date: '//version//' / '//date
857 
858 close(10, status='keep')
859 
860 !-------- Conversion of time quantities --------
861 
862 year_zero = year_zero*year_sec ! a --> s
863 time_init = time_init0*year_sec ! a --> s
864 time_end = time_end0*year_sec ! a --> s
865 dtime = dtime0*year_sec ! a --> s
866 dtime_temp = dtime_temp0*year_sec ! a --> s
867 #if (REBOUND==2)
868 dtime_wss = dtime_wss0*year_sec ! a --> s
869 #endif
870 #if (CHASM==2)
871 time_chasm_init = time_chasm_init0 *year_sec ! a --> s
872 time_chasm_end = time_chasm_end0 *year_sec ! a --> s
873 #endif
874 dtime_ser = dtime_ser0*year_sec ! a --> s
875 #if (OUTPUT==1 || OUTPUT==3)
876 dtime_out = dtime_out0*year_sec ! a --> s
877 #endif
878 #if (OUTPUT==2 || OUTPUT==3)
879 do n=1, n_output
880  time_output(n) = time_output0(n)*year_sec ! a --> s
881 end do
882 #endif
883 
884 if (.not.approx_integer_multiple(dtime_temp, dtime, eps_sp_dp)) then
885  errormsg = ' >>> sico_init: dtime_temp must be a multiple of dtime!'
886  call error(errormsg)
887 end if
888 
889 #if (REBOUND==2)
890 if (.not.approx_integer_multiple(dtime_wss, dtime, eps_sp_dp)) then
891  errormsg = ' >>> sico_init: dtime_wss must be a multiple of dtime!'
892  call error(errormsg)
893 end if
894 #endif
895 
896 if (.not.approx_integer_multiple(dtime_ser, dtime, eps_sp_dp)) then
897  errormsg = ' >>> sico_init: dtime_ser must be a multiple of dtime!'
898  call error(errormsg)
899 end if
900 
901 #if (OUTPUT==1 || OUTPUT==3)
902 if (.not.approx_integer_multiple(dtime_out, dtime, eps_sp_dp)) then
903  errormsg = ' >>> sico_init: dtime_out must be a multiple of dtime!'
904  call error(errormsg)
905 end if
906 #endif
907 
908 time = time_init
909 
910 !-------- Read accumulation measurements --------
911 
912 #if (ACC_PRESENT==1)
913 
914 do i=0, imax
915 do j=0, jmax
916  accum_present(j,i) = acc_present_val
917 end do
918 end do
919 
920 #endif
921 
922 ! ------ Conversion mm/a water or (ice+dust) equivalent
923 ! --> m/s (ice+dust) equivalent
924 
925 do i=0, imax
926 do j=0, jmax
927 
928 #if (ACC_UNIT==1)
929 
930  accum_present(j,i) = accum_present(j,i) &
931  *(1.0e-03_dp/year_sec)*(rho_w/rho_i) &
932  *(1.0_dp/(1.0_dp-frac_dust))
933 
934 #elif (ACC_UNIT==2)
935 
936  accum_present(j,i) = accum_present(j,i)*(1.0e-03_dp/year_sec)
937 
938 #endif
939 
940 end do
941 end do
942 
943 !-------- Mean accumulation --------
944 
945 #if (ACC_UNIT==1)
946 
947 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho_i) &
948  *(1.0_dp/(1.0_dp-frac_dust))
949 ! ! mm/a water equiv. --> m/s (ice+dust) equiv.
950 
951 #elif (ACC_UNIT==2)
952 
953 mean_accum = mean_accum*(1.0e-03_dp/year_sec)
954 ! ! mm/a (ice+dust) equiv. --> m/s (ice+dust) equiv.
955 
956 #endif
957 
958 !-------- Read chasm mask --------
959 
960 #if (CHASM==2)
961 
962 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
963  trim(mask_chasm_file)
964 
965 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
966 
967 if (ios /= 0) then
968  errormsg = ' >>> sico_init: Error when opening the chasm mask file!'
969  call error(errormsg)
970 end if
971 
972 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
973 
974 do j=jmax, 0, -1
975  read(24, fmt=trim(fmt4)) (maske_chasm(j,i), i=0,imax)
976 end do
977 
978 close(24, status='keep')
979 
980 #endif
981 
982 !-------- Read data for z_sl --------
983 
984 #if (SEA_LEVEL==3)
985 
986 errormsg = ' >>> sico_init: SEA_LEVEL==3 not allowed for smars application!'
987 call error(errormsg)
988 
989 #endif
990 
991 !-------- Reading of insolation data --------
992 
993 insol_ma_90 = 0.0_dp ! Assignment of dummy values
994 obl_data = 0.0_dp
995 ecc_data = 0.0_dp
996 ave_data = 0.0_dp
997 cp_data = 0.0_dp
998 
999 #if (TSURFACE==5 || TSURFACE==6)
1000 
1001 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1002  trim(insol_ma_90s_file)
1003 
1004 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1005 
1006 if (ios /= 0) then
1007  errormsg = ' >>> sico_init: Error when opening the insolation-data file!'
1008  call error(errormsg)
1009 end if
1010 
1011 read(21, fmt=*) ch_dummy, insol_time_min, insol_time_stp, insol_time_max
1012 
1013 if (ch_dummy /= '#') then
1014  errormsg = ' >>> sico_init: insol_time_min, insol_time_stp, insol_time_max' &
1015  // end_of_line &
1016  //' not defined in insolation-data file!'
1017  call error(errormsg)
1018 end if
1019 
1021 
1022 if (ndata_insol > 100000) then
1023  errormsg = ' >>> sico_init: Too many data in insolation-data file!'
1024  call error(errormsg)
1025 end if
1026 
1027 do n=0, ndata_insol
1028  read(21, fmt=*) d_dummy, ecc_data(n), obl_data(n), cp_data(n), ave_data(n), insol_ma_90(n)
1029  obl_data(n) = obl_data(n) *pi_180
1030  ave_data(n) = ave_data(n) *pi_180
1031 end do
1032 
1033 close(21, status='keep')
1034 
1035 #endif
1036 
1037 !-------- Determination of the normal geothermal heat flux
1038 ! (without the possible contribution from an active chasm area) --------
1039 
1040 #if (Q_GEO_MOD==1)
1041 
1042 ! ------ Constant value
1043 
1044 do i=0, imax
1045 do j=0, jmax
1046  q_geo_normal(j,i) = q_geo *1.0e-03_dp ! mW/m2 -> W/m2
1047 end do
1048 end do
1049 
1050 #elif (Q_GEO_MOD==2)
1051 
1052 ! ------ Read data from file
1053 
1054 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1055  trim(q_geo_file)
1056 
1057 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1058 
1059 if (ios /= 0) then
1060  errormsg = ' >>> sico_init: Error when opening the qgeo file!'
1061  call error(errormsg)
1062 end if
1063 
1064 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1065 
1066 do j=jmax, 0, -1
1067  read(21, fmt=*) (q_geo_normal(j,i), i=0,imax)
1068 end do
1069 
1070 close(21, status='keep')
1071 
1072 do i=0, imax
1073 do j=0, jmax
1074  q_geo_normal(j,i) = q_geo_normal(j,i) *1.0e-03_dp ! mW/m2 -> W/m2
1075 end do
1076 end do
1077 
1078 #endif
1079 
1080 !-------- Reading of tabulated kei function--------
1081 
1082 #if (REBOUND==0 || REBOUND==1)
1083 
1084 kei = 0.0_dp; n_data_kei = 0; kei_r_max = 0.0_dp; kei_r_incr = 0.0_dp
1085  ! dummy values
1086 #elif (REBOUND==2)
1087 
1088 call read_kei()
1089 
1090 #endif
1091 
1092 !-------- Determination of the time lag
1093 ! of the relaxing asthenosphere --------
1094 
1095 #if (REBOUND==1 || REBOUND==2)
1096 
1097 #if (TIME_LAG_MOD==1)
1098 
1099 time_lag_asth = time_lag*year_sec ! a -> s
1100 
1101 #elif (TIME_LAG_MOD==2)
1102 
1103 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1104  trim(time_lag_file)
1105 
1106 open(29, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1107 
1108 if (ios /= 0) then
1109  errormsg = ' >>> sico_init: Error when opening the time-lag file!'
1110  call error(errormsg)
1111 end if
1112 
1113 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1114 
1115 do j=jmax, 0, -1
1116  read(29, fmt=*) (time_lag_asth(j,i), i=0,imax)
1117 end do
1118 
1119 close(29, status='keep')
1120 
1121 time_lag_asth = time_lag_asth*year_sec ! a -> s
1122 
1123 #endif
1124 
1125 #elif (REBOUND==0)
1126 
1127 time_lag_asth = 0.0_dp ! dummy values
1128 
1129 #endif
1130 
1131 !-------- Determination of the flexural rigidity of the lithosphere --------
1132 
1133 #if (REBOUND==2)
1134 
1135 #if (FLEX_RIG_MOD==1)
1136 
1137 flex_rig_lith = flex_rig
1138 
1139 #elif (FLEX_RIG_MOD==2)
1140 
1141 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1142  trim(flex_rig_file)
1143 
1144 open(29, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1145 
1146 if (ios /= 0) then
1147  errormsg = ' >>> sico_init: Error when opening the flex-rig file!'
1148  call error(errormsg)
1149 end if
1150 
1151 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1152 
1153 do j=jmax, 0, -1
1154  read(29, fmt=*) (flex_rig_lith(j,i), i=0,imax)
1155 end do
1156 
1157 close(29, status='keep')
1158 
1159 #endif
1160 
1161 #elif (REBOUND==0 || REBOUND==1)
1162 
1163 flex_rig_lith = 0.0_dp ! dummy values
1164 
1165 #endif
1166 
1167 !-------- Definition of initial values --------
1168 
1169 ! ------ Present topography
1170 
1171 #if (ANF_DAT==1)
1172 
1173 call topography1(dxi, deta)
1174 
1175 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1176 
1177 call boundary(time_init, dtime, dxi, deta, &
1178  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1179 
1180 where ((maske==0_i1b).or.(maske==3_i1b))
1181  ! grounded or floating ice
1183 elsewhere ! maske==1_i1b or 2_i1b, ice-free land or sea
1184  as_perp_apl = 0.0_dp
1185 end where
1186 
1187 smb_corr = 0.0_dp
1188 
1189 q_bm = 0.0_dp
1190 q_tld = 0.0_dp
1191 q_b_tot = 0.0_dp
1192 
1193 p_b_w = 0.0_dp
1194 h_w = 0.0_dp
1195 
1196 #if (!defined(TEMP_INIT) || TEMP_INIT==1)
1198 #elif (TEMP_INIT==2)
1200 #elif (TEMP_INIT==3)
1202 #elif (TEMP_INIT==4)
1204 #elif (TEMP_INIT==5)
1205  call init_temp_water_age_1_5(z_sl, anfdatname)
1206 #else
1207  errormsg = ' >>> sico_init: TEMP_INIT must be between 1 and 5!'
1208  call error(errormsg)
1209 #endif
1210 
1211 #if (ENHMOD==1)
1212  call calc_enhance_1()
1213 #elif (ENHMOD==2)
1214  call calc_enhance_2()
1215 #elif (ENHMOD==3)
1216  call calc_enhance_3(time_init)
1217 #elif (ENHMOD==4)
1218  call calc_enhance_4()
1219 #elif (ENHMOD==5)
1220  call calc_enhance_5()
1221 #else
1222  errormsg = ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
1223  call error(errormsg)
1224 #endif
1225 
1226 vx_m_sia = 0.0_dp
1227 vy_m_sia = 0.0_dp
1228 vx_m_ssa = 0.0_dp
1229 vy_m_ssa = 0.0_dp
1230 
1231 #if (defined(VISC_INIT_SSA))
1232  vis_ave_g = visc_init_ssa
1233 #else
1234  vis_ave_g = 1.0e+15_dp ! Pa s
1235 #endif
1236 
1237 vis_int_g = 0.0_dp
1238 
1239 ! ------ Ice-free, relaxed bedrock
1240 
1241 #elif (ANF_DAT==2)
1242 
1243 call topography2(dxi, deta)
1244 
1245 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1246 
1247 call boundary(time_init, dtime, dxi, deta, &
1248  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1249 
1250 as_perp_apl = 0.0_dp
1251 
1252 smb_corr = 0.0_dp
1253 
1254 q_bm = 0.0_dp
1255 q_tld = 0.0_dp
1256 q_b_tot = 0.0_dp
1257 
1258 p_b_w = 0.0_dp
1259 h_w = 0.0_dp
1260 
1261 call init_temp_water_age_2()
1262 
1263 #if (ENHMOD==1)
1264  call calc_enhance_1()
1265 #elif (ENHMOD==2)
1266  call calc_enhance_2()
1267 #elif (ENHMOD==3)
1268  call calc_enhance_3(time_init)
1269 #elif (ENHMOD==4)
1270  call calc_enhance_4()
1271 #elif (ENHMOD==5)
1272  call calc_enhance_5()
1273 #else
1274  errormsg = ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
1275  call error(errormsg)
1276 #endif
1277 
1278 vx_m_sia = 0.0_dp
1279 vy_m_sia = 0.0_dp
1280 vx_m_ssa = 0.0_dp
1281 vy_m_ssa = 0.0_dp
1282 
1283 #if (defined(VISC_MIN) && defined(VISC_MAX))
1284  vis_ave_g = visc_max
1285 #else
1286  vis_ave_g = 1.0e+25_dp ! Pa s
1287 #endif
1288 
1289 vis_int_g = 0.0_dp
1290 
1291 ! ------ Read initial state from output of previous simulation
1292 
1293 #elif (ANF_DAT==3)
1294 
1295 call topography3(dxi, deta, z_sl, anfdatname)
1296 
1297 call boundary(time_init, dtime, dxi, deta, &
1298  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1299 
1300 where ((maske==0_i1b).or.(maske==3_i1b))
1301  ! grounded or floating ice
1303 elsewhere ! maske==1_i1b or 2_i1b, ice-free land or sea
1304  as_perp_apl = 0.0_dp
1305 end where
1306 
1307 smb_corr = 0.0_dp
1308 
1309 q_b_tot = q_bm + q_tld
1310 
1311 #endif
1312 
1313 !-------- Inner-point flag --------
1314 
1315 flag_inner_point = .true.
1316 
1317 flag_inner_point(0,:) = .false.
1318 flag_inner_point(jmax,:) = .false.
1319 
1320 flag_inner_point(:,0) = .false.
1321 flag_inner_point(:,imax) = .false.
1322 
1323 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
1324 
1325 #if (GRID==0 || GRID==1)
1326 
1327 do ir=-imax, imax
1328 do jr=-jmax, jmax
1329  dist_dxdy(jr,ir) = sqrt( (real(ir,dp)*dxi)**2 + (real(jr,dp)*deta)**2 )
1330  ! distortion due to stereographic projection not accounted for
1331 end do
1332 end do
1333 
1334 #endif
1335 
1336 !-------- Initial velocities --------
1337 
1338 call calc_temp_melt()
1339 call flag_update_gf_gl_cf()
1340 
1341 #if (DYNAMICS==1 || DYNAMICS==2)
1342 
1343 call calc_vxy_b_sia(time, z_sl)
1344 call calc_vxy_sia(dzeta_c, dzeta_t)
1345 
1346 #if (MARGIN==3 || DYNAMICS==2)
1347 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t)
1348 #endif
1349 
1350 call calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
1351 
1352 #if (MARGIN==3)
1353 call calc_vz_floating(z_sl, dxi, deta, dzeta_c)
1354 #endif
1355 
1356 #elif (DYNAMICS==0)
1357 
1358 call calc_vxy_static()
1359 call calc_vz_static()
1360 
1361 #else
1362 
1363 errormsg = ' >>> sico_init: DYNAMICS must be either 0, 1 or 2!'
1364 call error(errormsg)
1365 
1366 #endif
1367 
1368 call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
1369 
1370 !-------- Initial enthalpies --------
1371 
1372 #if (CALCMOD==0 || CALCMOD==-1)
1373 
1374 do i=0, imax
1375 do j=0, jmax
1376 
1377  do kc=0, kcmax
1378  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1379  end do
1380 
1381  do kt=0, ktmax
1382  enth_t(kt,j,i) = enth_c(0,j,i)
1383  end do
1384 
1385 end do
1386 end do
1387 
1388 #elif (CALCMOD==1)
1389 
1390 do i=0, imax
1391 do j=0, jmax
1392 
1393  do kc=0, kcmax
1394  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1395  end do
1396 
1397  if ( (maske(j,i) == 0_i1b).and.(n_cts(j,i) == 1_i1b) ) then
1398  do kt=0, ktmax
1399  enth_t(kt,j,i) = enth_fct_temp_omega(temp_t_m(kt,j,i), omega_t(kt,j,i))
1400  end do
1401  else
1402  do kt=0, ktmax
1403  enth_t(kt,j,i) = enth_c(0,j,i)
1404  end do
1405  end if
1406 
1407 end do
1408 end do
1409 
1410 #elif (CALCMOD==2 || CALCMOD==3)
1411 
1412 do i=0, imax
1413 do j=0, jmax
1414 
1415  do kc=0, kcmax
1416  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), omega_c(kc,j,i))
1417  end do
1418 
1419  do kt=0, ktmax
1420  enth_t(kt,j,i) = enth_c(0,j,i)
1421  end do
1422 
1423 end do
1424 end do
1425 
1426 #else
1427 
1428 errormsg = ' >>> sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!'
1429 call error(errormsg)
1430 
1431 #endif
1432 
1433 !-------- Initialize time-series files --------
1434 
1435 ! ------ Time-series file for the ice sheet on the whole
1436 
1437 filename_with_path = trim(outpath)//'/'//trim(runname)//'.ser'
1438 
1439 open(12, iostat=ios, file=trim(filename_with_path), status='new')
1440 
1441 if (ios /= 0) then
1442  errormsg = ' >>> sico_init: Error when opening the ser file!'
1443  call error(errormsg)
1444 end if
1445 
1446 write(12,1102)
1447 write(12,1103)
1448 
1449 1102 format(' t(a) D_Ts(C) z_sl(m)',/, &
1450  ' V(m^3) V_g(m^3) V_f(m^3)', &
1451  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1452  ' V_sle(m) V_t(m^3)', &
1453  ' A_t(m^2)',/, &
1454  ' H_max(m) H_t_max(m)', &
1455  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
1456 1103 format('----------------------------------------------------', &
1457  '---------------------------------------')
1458 
1459 ! ------ Time-series file for deep boreholes
1460 
1461 n_core = 0 ! No boreholes defined
1462 
1463 filename_with_path = trim(outpath)//'/'//trim(runname)//'.core'
1464 
1465 open(14, iostat=ios, file=trim(filename_with_path), status='new')
1466 
1467 write(14,'(1x,a)') '---------------------'
1468 write(14,'(1x,a)') 'No boreholes defined.'
1469 write(14,'(1x,a)') '---------------------'
1470 
1471 !-------- Output of the initial state --------
1472 
1473 #if (defined(OUTPUT_INIT))
1474 
1475 #if (OUTPUT_INIT==0)
1476  flag_init_output = .false.
1477 #elif (OUTPUT_INIT==1)
1478  flag_init_output = .true.
1479 #else
1480  errormsg = ' >>> sico_init: OUTPUT_INIT must be either 0 or 1!'
1481  call error(errormsg)
1482 #endif
1483 
1484 #else
1485  flag_init_output = .true. ! default for undefined parameter OUTPUT_INIT
1486 #endif
1487 
1488 #if (OUTPUT==1)
1489 
1490 #if (ERGDAT==0)
1491  flag_3d_output = .false.
1492 #elif (ERGDAT==1)
1493  flag_3d_output = .true.
1494 #else
1495  errormsg = ' >>> sico_init: ERGDAT must be either 0 or 1!'
1496  call error(errormsg)
1497 #endif
1498 
1499  if (flag_init_output) &
1500  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1501  flag_3d_output, ndat2d, ndat3d)
1502 
1503 #elif (OUTPUT==2)
1504 
1505 if (time_output(1) <= time_init+eps) then
1506 
1507 #if (ERGDAT==0)
1508  flag_3d_output = .false.
1509 #elif (ERGDAT==1)
1510  flag_3d_output = .true.
1511 #else
1512  errormsg = ' >>> sico_init: ERGDAT must be either 0 or 1!'
1513  call error(errormsg)
1514 #endif
1515 
1516  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1517  flag_3d_output, ndat2d, ndat3d)
1518 
1519 end if
1520 
1521 #elif (OUTPUT==3)
1522 
1523  flag_3d_output = .false.
1524 
1525  if (flag_init_output) &
1526  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1527  flag_3d_output, ndat2d, ndat3d)
1528 
1529 if (time_output(1) <= time_init+eps) then
1530 
1531  flag_3d_output = .true.
1532 
1533  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1534  flag_3d_output, ndat2d, ndat3d)
1535 
1536 end if
1537 
1538 #else
1539 
1540  errormsg = ' >>> sico_init: OUTPUT must be either 1, 2 or 3!'
1541  call error(errormsg)
1542 
1543 #endif
1544 
1545 if (flag_init_output) then
1546  call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1547  call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1548 end if
1549 
1550 end subroutine sico_init
1551 
1552 !-------------------------------------------------------------------------------
1553 !> Definition of the initial surface and bedrock topography
1554 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1555 !! For present-day initial topography.
1556 !<------------------------------------------------------------------------------
1557 subroutine topography1(dxi, deta)
1558 
1559 #if (GRID==0 || GRID==1)
1561 #endif
1562 
1563  use metric_m
1564  use topograd_m
1565 
1566 implicit none
1567 
1568 real(dp), intent(out) :: dxi, deta
1569 
1570 integer(i4b) :: i, j, n
1571 integer(i4b) :: ios, n_dummy
1572 real(dp) :: d_dummy
1573 real(dp) :: xi0, eta0
1574 character :: ch_dummy
1575 
1576 character(len= 8) :: ch_imax
1577 character(len=128) :: fmt4
1578 character(len=256) :: filename_with_path
1579 
1580 write(ch_imax, fmt='(i8)') imax
1581 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
1582 
1583 !-------- Read topography --------
1584 
1585 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1586  trim(zs_present_file)
1587 
1588 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1589 
1590 if (ios /= 0) then
1591  errormsg = ' >>> topography1: Error when opening the zs file!'
1592  call error(errormsg)
1593 end if
1594 
1595 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1596 
1597 do j=jmax, 0, -1
1598  read(21, fmt=*) (zs(j,i), i=0,imax)
1599 end do
1600 
1601 close(21, status='keep')
1602 
1603 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1604  trim(zl_present_file)
1605 
1606 open(22, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1607 
1608 if (ios /= 0) then
1609  errormsg = ' >>> topography1: Error when opening the zl file!'
1610  call error(errormsg)
1611 end if
1612 
1613 do n=1, 6; read(22, fmt='(a)') ch_dummy; end do
1614 
1615 do j=jmax, 0, -1
1616  read(22, fmt=*) (zl(j,i), i=0,imax)
1617 end do
1618 
1619 close(22, status='keep')
1620 
1621 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1622  trim(zl0_file)
1623 
1624 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1625 
1626 if (ios /= 0) then
1627  errormsg = ' >>> topography1: Error when opening the zl0 file!'
1628  call error(errormsg)
1629 end if
1630 
1631 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
1632 
1633 do j=jmax, 0, -1
1634  read(23, fmt=*) (zl0(j,i), i=0,imax)
1635 end do
1636 
1637 close(23, status='keep')
1638 
1639 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1640  trim(mask_present_file)
1641 
1642 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
1643 
1644 if (ios /= 0) then
1645  errormsg = ' >>> topography1: Error when opening the mask file!'
1646  call error(errormsg)
1647 end if
1648 
1649 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1650 
1651 do j=jmax, 0, -1
1652  read(24, fmt=trim(fmt4)) (maske(j,i), i=0,imax)
1653 end do
1654 
1655 close(24, status='keep')
1656 
1657 !-------- Further stuff --------
1658 
1659 dxi = dx *1000.0_dp ! km -> m
1660 deta = dx *1000.0_dp ! km -> m
1661 
1662 xi0 = x0 *1000.0_dp ! km -> m
1663 eta0 = y0 *1000.0_dp ! km -> m
1664 
1665 do i=0, imax
1666 do j=0, jmax
1667 
1668  if (maske(j,i) <= 1_i1b) then
1669  zb(j,i) = zl(j,i)
1670  else ! (maske(j,i>=2_i1b)
1671  errormsg = ' >>> topography1: ' &
1672  //'maske(j,i)>=2 not allowed for initial topography!'
1673  call error(errormsg)
1674  end if
1675 
1676  zs(j,i) = zs(j,i) *1000.0_dp
1677  zb(j,i) = zb(j,i) *1000.0_dp ! km --> m
1678  zl(j,i) = zl(j,i) *1000.0_dp
1679  zl0(j,i) = zl0(j,i)*1000.0_dp
1680 
1681  xi(i) = xi0 + real(i,dp)*dxi
1682  eta(j) = eta0 + real(j,dp)*deta
1683 
1684  zm(j,i) = zb(j,i)
1685  n_cts(j,i) = -1_i1b
1686  kc_cts(j,i) = 0
1687 
1688  h_c(j,i) = zs(j,i)-zm(j,i)
1689  h_t(j,i) = 0.0_dp
1690 
1691  dzs_dtau(j,i) = 0.0_dp
1692  dzm_dtau(j,i) = 0.0_dp
1693  dzb_dtau(j,i) = 0.0_dp
1694  dzl_dtau(j,i) = 0.0_dp
1695  dh_c_dtau(j,i) = 0.0_dp
1696  dh_t_dtau(j,i) = 0.0_dp
1697 
1698 end do
1699 end do
1700 
1701 maske_old = maske
1702 
1703 !-------- Geographic coordinates, metric tensor,
1704 ! gradients of the topography --------
1705 
1706 do i=0, imax
1707 do j=0, jmax
1708 
1709 #if (GRID==0 || GRID==1) /* Stereographic projection */
1710 
1711  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1712  lambda0, phi0, lambda(j,i), phi(j,i))
1713 
1714 #elif (GRID==2) /* Geographic coordinates */
1715 
1716  lambda(j,i) = xi(i)
1717  phi(j,i) = eta(j)
1718 
1719 #endif
1720 
1721 end do
1722 end do
1723 
1724 call metric()
1725 
1726 #if (TOPOGRAD==0)
1727 call topograd_1(dxi, deta, 1)
1728 #elif (TOPOGRAD==1)
1729 call topograd_2(dxi, deta, 1)
1730 #endif
1731 
1732 !-------- Corresponding area of grid points --------
1733 
1734 do i=0, imax
1735 do j=0, jmax
1736  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1737 end do
1738 end do
1739 
1740 end subroutine topography1
1741 
1742 !-------------------------------------------------------------------------------
1743 !> Definition of the initial surface and bedrock topography
1744 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1745 !! For ice-free initial topography with relaxed lithosphere.
1746 !<------------------------------------------------------------------------------
1747 subroutine topography2(dxi, deta)
1748 
1749 #if (GRID==0 || GRID==1)
1751 #endif
1752 
1753  use metric_m
1754  use topograd_m
1755 
1756 implicit none
1757 
1758 real(dp), intent(out) :: dxi, deta
1759 
1760 integer(i4b) :: i, j, n
1761 integer(i4b) :: ios
1762 real(dp) :: xi0, eta0
1763 character :: ch_dummy
1764 
1765 character(len= 8) :: ch_imax
1766 character(len=128) :: fmt4
1767 character(len=256) :: filename_with_path
1768 
1769 write(ch_imax, fmt='(i8)') imax
1770 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
1771 
1772 !-------- Read topography --------
1773 
1774 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1775  trim(zl0_file)
1776 
1777 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1778 
1779 if (ios /= 0) then
1780  errormsg = ' >>> topography2: Error when opening the zl0 file!'
1781  call error(errormsg)
1782 end if
1783 
1784 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
1785 
1786 do j=jmax, 0, -1
1787  read(23, fmt=*) (zl0(j,i), i=0,imax)
1788 end do
1789 
1790 close(23, status='keep')
1791 
1792 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1793  trim(mask_present_file)
1794 
1795 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
1796 
1797 if (ios /= 0) then
1798  errormsg = ' >>> topography2: Error when opening the mask file!'
1799  call error(errormsg)
1800 end if
1801 
1802 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1803 
1804 do j=jmax, 0, -1
1805  read(24, fmt=trim(fmt4)) (maske(j,i), i=0,imax)
1806 end do
1807 
1808 close(24, status='keep')
1809 
1810 !-------- Further stuff --------
1811 
1812 dxi = dx *1000.0_dp ! km -> m
1813 deta = dx *1000.0_dp ! km -> m
1814 
1815 xi0 = x0 *1000.0_dp ! km -> m
1816 eta0 = y0 *1000.0_dp ! km -> m
1817 
1818 do i=0, imax
1819 do j=0, jmax
1820 
1821  if (maske(j,i) <= 1_i1b) then
1822  maske(j,i) = 1_i1b
1823  zs(j,i) = zl0(j,i)
1824  zb(j,i) = zl0(j,i)
1825  zl(j,i) = zl0(j,i)
1826  else ! (maske(j,i)>=2_i1b)
1827  errormsg = ' >>> topography2: ' &
1828  //'maske(j,i)>=2 not allowed for initial topography!'
1829  call error(errormsg)
1830  end if
1831 
1832  zs(j,i) = zs(j,i) *1000.0_dp
1833  zb(j,i) = zb(j,i) *1000.0_dp ! km --> m
1834  zl(j,i) = zl(j,i) *1000.0_dp
1835  zl0(j,i) = zl0(j,i)*1000.0_dp
1836 
1837  xi(i) = xi0 + real(i,dp)*dxi
1838  eta(j) = eta0 + real(j,dp)*deta
1839 
1840  zm(j,i) = zb(j,i)
1841  n_cts(j,i) = -1_i1b
1842  kc_cts(j,i) = 0
1843 
1844  h_c(j,i) = 0.0_dp
1845  h_t(j,i) = 0.0_dp
1846 
1847  dzs_dtau(j,i) = 0.0_dp
1848  dzm_dtau(j,i) = 0.0_dp
1849  dzb_dtau(j,i) = 0.0_dp
1850  dzl_dtau(j,i) = 0.0_dp
1851  dh_c_dtau(j,i) = 0.0_dp
1852  dh_t_dtau(j,i) = 0.0_dp
1853 
1854 end do
1855 end do
1856 
1857 maske_old = maske
1858 
1859 !-------- Geographic coordinates, metric tensor,
1860 ! gradients of the topography --------
1861 
1862 do i=0, imax
1863 do j=0, jmax
1864 
1865 #if (GRID==0 || GRID==1) /* Stereographic projection */
1866 
1867  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1868  lambda0, phi0, lambda(j,i), phi(j,i))
1869 
1870 #elif (GRID==2) /* Geographic coordinates */
1871 
1872  lambda(j,i) = xi(i)
1873  phi(j,i) = eta(j)
1874 
1875 #endif
1876 
1877 end do
1878 end do
1879 
1880 call metric()
1881 
1882 #if (TOPOGRAD==0)
1883 call topograd_1(dxi, deta, 1)
1884 #elif (TOPOGRAD==1)
1885 call topograd_2(dxi, deta, 1)
1886 #endif
1887 
1888 !-------- Corresponding area of grid points --------
1889 
1890 do i=0, imax
1891 do j=0, jmax
1892  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1893 end do
1894 end do
1895 
1896 end subroutine topography2
1897 
1898 !-------------------------------------------------------------------------------
1899 !> Definition of the initial surface and bedrock topography
1900 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1901 !! For initial topography from previous simulation.
1902 !<------------------------------------------------------------------------------
1903 subroutine topography3(dxi, deta, z_sl, anfdatname)
1904 
1905  use read_m, only : read_erg_nc
1906 
1907 #if (GRID==0 || GRID==1)
1909 #endif
1910 
1911  use metric_m
1912  use topograd_m
1913 
1914 implicit none
1915 
1916 character(len=100), intent(in) :: anfdatname
1917 
1918 real(dp), intent(out) :: dxi, deta, z_sl
1919 
1920 integer(i4b) :: i, j, n
1921 integer(i4b) :: ios
1922 character(len=256) :: filename_with_path
1923 character :: ch_dummy
1924 
1925 !-------- Read data from time-slice file of previous simulation --------
1926 
1927 call read_erg_nc(z_sl, anfdatname)
1928 
1929 !-------- Read topography of the relaxed bedrock --------
1930 
1931 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1932  trim(zl0_file)
1933 
1934 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1935 
1936 if (ios /= 0) then
1937  errormsg = ' >>> topography3: Error when opening the zl0 file!'
1938  call error(errormsg)
1939 end if
1940 
1941 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
1942 
1943 do j=jmax, 0, -1
1944  read(23, fmt=*) (zl0(j,i), i=0,imax)
1945 end do
1946 
1947 close(23, status='keep')
1948 
1949 !-------- Further stuff --------
1950 
1951 dxi = dx *1000.0_dp ! km -> m
1952 deta = dx *1000.0_dp ! km -> m
1953 
1954 zl0 = zl0 *1000.0_dp ! km -> m
1955 
1956 !-------- Geographic coordinates, metric tensor,
1957 ! gradients of the topography --------
1958 
1959 do i=0, imax
1960 do j=0, jmax
1961 
1962 #if (GRID==0 || GRID==1) /* Stereographic projection */
1963 
1964  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1965  lambda0, phi0, lambda(j,i), phi(j,i))
1966 
1967 #elif (GRID==2) /* Geographic coordinates */
1968 
1969  lambda(j,i) = xi(i)
1970  phi(j,i) = eta(j)
1971 
1972 #endif
1973 
1974 end do
1975 end do
1976 
1977 call metric()
1978 
1979 #if (TOPOGRAD==0)
1980 call topograd_1(dxi, deta, 1)
1981 #elif (TOPOGRAD==1)
1982 call topograd_2(dxi, deta, 1)
1983 #endif
1984 
1985 !-------- Corresponding area of grid points --------
1986 
1987 do i=0, imax
1988 do j=0, jmax
1989  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1990 end do
1991 end do
1992 
1993 end subroutine topography3
1994 
1995 !-------------------------------------------------------------------------------
1996 
1997 end module sico_init_m
1998 !
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)
subroutine, public read_kei()
Reading of the tabulated kei function.
Definition: read_m.F90:1061
logical, dimension(0:jmax, 0:imax) flag_inner_point
flag_inner_point(j,i): Inner-point flag. .true.: inner point, .false.: margin point ...
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
integer(i4b) insol_time_stp
insol_time_stp: Time step of the data values for the insolation etc.
Definition: sico_vars_m.F90:45
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), dimension(0:100000) ecc_data
ecc_data(n): Data values for the eccentricity
Definition: sico_vars_m.F90:54
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
Material properties of ice: Rate factor, heat conductivity, specific heat (heat capacity), creep function, viscosity.
real(dp), parameter eps
eps: Small number
subroutine, public read_phys_para()
Reading of physical parameters.
Definition: read_m.F90:1119
real(dp), dimension(0:jmax, 0:imax) q_geo_normal
q_geo_normal(j,i): Geothermal heat flux for normal (non-chasma) areas
Definition: sico_vars_m.F90:74
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(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)
integer(i4b) insol_time_min
insol_time_min: Minimum time of the data values for the insolation etc.
Definition: sico_vars_m.F90:43
real(dp), dimension(0:jmax, 0:imax) as_perp
as_perp(j,i): Accumulation-ablation function at the ice surface (SMB)
real(dp), dimension(0:jmax, 0:imax) area
area(j,i): Area of grid cell associated with grid point (i,j)
Computation of all components of the strain-rate tensor, the full effective strain rate and the shear...
Definition: calc_dxyz_m.F90:37
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) temp_t_m
temp_t_m(kt,j,i): Melting temperature in the lower (kt) ice domain
subroutine, public stereo_inv_ellipsoid(x_val, y_val, A, B, lambda0, phi0, lambda_val, phi_val)
Inverse stereographic projection for an ellipsoidal planet.
subroutine topograd_1(dxi, deta, n_switch)
Calculation of topography gradients on the staggered grid and on the grid points (the latter by secon...
Definition: topograd_m.F90:57
subroutine, public calc_enhance_5()
Computation of the flow enhancement factor. Case ENHMOD==5: minimal anisotropic enhancement factor fo...
real(dp) b
B: Semi-minor axis of the planet.
subroutine, public init_temp_water_age_1_2()
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==2: present-day initial topogr...
real(dp), dimension(0:jmax, 0:imax) zl
zl(j,i): Coordinate z of the lithosphere surface
subroutine 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
real(dp), dimension(0:jmax, 0:imax) q_tld
Q_tld(j,i): Water drainage rate from the temperate layer.
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars_m.F90:35
integer(i4b) insol_time_max
insol_time_max: Maximum time of the data values for the insolation etc.
Definition: sico_vars_m.F90:47
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.
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), parameter pi_180
pi_180: pi divided by 180 (-> deg to rad)
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) rho_c
RHO_C: Density of crustal material (dust)
Definition: sico_vars_m.F90:79
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) time_chasm_end
time_chasm_end: Final time for active chasma area
Definition: sico_vars_m.F90:72
real(dp) kei_r_max
kei_r_max: Maximum value of the argument r of the tabulated kei function
real(dp) year_zero
year_zero: SICOPOLIS year zero in astronomical year numbering [ = signed year CE (AD) ] ...
Computation of the melting and basal temperatures.
subroutine topograd_2(dxi, deta, n_switch)
Calculation of topography gradients on the staggered grid and on the grid points (the latter by fourt...
Definition: topograd_m.F90:203
real(dp), dimension(-190:10) kappa
KAPPA(n): Tabulated values for the heat conductivity of ice.
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
character(len=64) ch_domain_long
ch_domain_long: Long name of the computational domain
real(dp), dimension(0:jmax, 0:imax) accum_present
accum_present(j,i): Present-day accumulation rate at the ice surface (for EISMINT, ISMIP HEINO and the north and south polar caps of Mars)
real(dp) kappa_c
KAPPA_C: Heat conductivity of crustal material (dust)
Definition: sico_vars_m.F90:81
real(dp) rho_i
RHO_I: Density of ice.
Definition: sico_vars_m.F90:77
Definition of the components g11 and g22 of the metric tensor of the applied numerical coordinates...
Definition: metric_m.F90:37
real(dp) c_c
C_C: Specific heat of crustal material (dust)
Definition: sico_vars_m.F90:83
Computation of the forward or inverse stereographic projection, alternatively for a spherical or an e...
real(dp) rho_inv
rho_inv: Inverse of the density of ice-dust mixture
Definition: sico_vars_m.F90:85
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) time_chasm_init
time_chasm_init: Initial time for active chasma area
Definition: sico_vars_m.F90:70
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)
real(dp), dimension(0:100000) cp_data
cp_data(n): Data values for Laskar&#39;s climate parameter = eccentricity *sin(Laskar&#39;s longitude of peri...
Definition: sico_vars_m.F90:63
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.
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)
real(dp), dimension(-10000:10000) kei
kei(n): Tabulated values of the kei function (Kelvin function of zero order)
subroutine, public init_temp_water_age_1_3()
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==3: present-day initial topogr...
subroutine boundary(time, dtime, dxi, deta, delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
Main routine of boundary_m: Computation of the surface temperature (must be less than 0 deg C!) and o...
Definition: boundary_m.F90:61
real(dp), dimension(0:100000) obl_data
obl_data(n): Data values for the obliquity
Definition: sico_vars_m.F90:52
integer(i4b), parameter rcl2
rcl2: Maximum length of record for input mask files (with factor 3 safety margin) ...
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
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
integer(i1b), dimension(0:jmax, 0:imax) maske_chasm
maske_chasm(j,i): Chasma mask. 0: grounded ice, 1: ice-free land (normal area), 7: chasma area ...
Definition: sico_vars_m.F90:68
subroutine, public calc_enhance_3(time)
Computation of the flow enhancement factor. Case ENHMOD==3: two different values depending on time of...
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
integer(i4b), dimension((imax+1)*(jmax+1)) n2i
n2i: Conversion from linear index n to 2d index i
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) kei_r_incr
kei_r_incr: Increment of the argument r of the tabulated kei function
real(dp), dimension(0:ktmax) zeta_t
zeta_t(kt): Sigma coordinate zeta_t of grid point kt
subroutine, public init_temp_water_age_1_5(z_sl, filename)
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==5: present-day initial topogr...
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), dimension(0:100000) insol_ma_90
insol_ma_90(n): Data values for the mean-annual north- or south-polar insolation
Definition: sico_vars_m.F90:50
integer(i4b) n_data_kei
n_data_kei: Number of tabulated values of the kei function
character(len=256) errormsg
errormsg: Error-message string
real(dp) rho
RHO: Density of ice.
subroutine topography3(dxi, deta, z_sl, anfdatname)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
integer(i4b), dimension(0:jmax, 0:imax) p_weert
p_weert(j,i): Weertman exponent for the basal shear stress
subroutine, public init_temp_water_age_1_4()
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==4: present-day initial topogr...
subroutine, public init_temp_water_age_1_1()
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==1: present-day initial topogr...
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
real(dp), dimension(0:100000) ave_data
ave_data(n): Data values for the anomaly of vernal equinox (= 360 deg - Ls of perihelion ) ...
Definition: sico_vars_m.F90:57
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(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...