SICOPOLIS V5-dev  Revision 1368
sico_init_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : s i c o _ i n i t _ m
4 !
5 !> @file
6 !!
7 !! Initialisations for SICOPOLIS.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2018 Ralf Greve, Eduardo Flandez, Matthias Scheiter
12 !!
13 !! @section License
14 !!
15 !! This file is part of SICOPOLIS.
16 !!
17 !! SICOPOLIS is free software: you can redistribute it and/or modify
18 !! it under the terms of the GNU General Public License as published by
19 !! the Free Software Foundation, either version 3 of the License, or
20 !! (at your option) any later version.
21 !!
22 !! SICOPOLIS is distributed in the hope that it will be useful,
23 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
24 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 !! GNU General Public License for more details.
26 !!
27 !! You should have received a copy of the GNU General Public License
28 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
29 !<
30 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
31 
32 !-------------------------------------------------------------------------------
33 !> Initialisations for SICOPOLIS.
34 !<------------------------------------------------------------------------------
35 module sico_init_m
36 
37  use sico_types_m
39  use sico_vars_m
40  use error_m
41 
42  implicit none
43 
44  public
45 
46 contains
47 
48 !-------------------------------------------------------------------------------
49 !> Main routine of sico_init_m: Initialisations for SICOPOLIS.
50 !<------------------------------------------------------------------------------
51 subroutine sico_init(delta_ts, glac_index, &
52  mean_accum, &
53  dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
54  time, time_init, time_end, time_output, &
55  dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
56  z_sl, dzsl_dtau, z_mar, &
57  ii, jj, nn, &
58  ndat2d, ndat3d, n_output, &
59  runname)
60 
61  use compare_float_m
65 
66 #if (NETCDF > 1)
67  use netcdf
68  use nc_check_m
69 #endif
70 
72 
73  use boundary_m
75  use calc_enhance_m
77  use calc_vxy_m
78  use calc_vz_m
79  use calc_dxyz_m
81 
82  use output_m
83 
84 implicit none
85 
86 integer(i4b), intent(out) :: ii((imax+1)*(jmax+1)), &
87  jj((imax+1)*(jmax+1)), &
88  nn(0:jmax,0:imax)
89 integer(i4b), intent(out) :: ndat2d, ndat3d
90 integer(i4b), intent(out) :: n_output
91 real(dp), intent(out) :: delta_ts, glac_index
92 real(dp), intent(out) :: mean_accum
93 real(dp), intent(out) :: dtime, dtime_temp, dtime_wss, &
94  dtime_out, dtime_ser
95 real(dp), intent(out) :: time, time_init, time_end, time_output(100)
96 real(dp), intent(out) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
97 real(dp), intent(out) :: z_sl, dzsl_dtau, z_mar
98 character(len=100), intent(out) :: runname
99 
100 integer(i4b) :: i, j, kc, kt, kr, m, n, ir, jr
101 integer(i4b) :: ios, ios1, ios2, ios3, ios4
102 integer(i4b) :: ierr
103 real(dp) :: dtime0, dtime_temp0, dtime_wss0, dtime_out0, dtime_ser0
104 real(dp) :: time_init0, time_end0, time_output0(100)
105 real(dp) :: d_dummy
106 character(len=100) :: anfdatname, target_topo_dat_name
107 character(len=256) :: filename_with_path
108 character(len=256) :: shell_command
109 character(len=3) :: ch_nc_test
110 character :: ch_dummy
111 logical :: flag_init_output, flag_3d_output
112 
113 #if (defined(SMB_CORR_FILE))
114 real(sp), dimension(0:IMAX,0:JMAX) :: smb_corr_in_conv
115 #endif
116 
117 #if (NETCDF > 1)
118 integer(i4b) :: ncid, ncv
119 ! ncid: File ID
120 ! ncv: Variable ID
121 #endif
122 
123 character(len=64), parameter :: thisroutine = 'sico_init'
124 
125 character(len=64), parameter :: fmt1 = '(a)', &
126  fmt2 = '(a,i4)', &
127  fmt2a = '(a,i0)', &
128  fmt3 = '(a,es12.4)'
129 
130 character(len= 8) :: ch_imax
131 character(len=128) :: fmt4
132 
133 write(ch_imax, fmt='(i8)') imax
134 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
135 
136 write(unit=6, fmt='(a)') ' '
137 write(unit=6, fmt='(a)') ' -------- sico_init --------'
138 write(unit=6, fmt='(a)') ' '
139 
140 !-------- Name of the computational domain --------
141 
142 #if (defined(ANT))
143 ch_domain_long = 'Antarctica'
144 ch_domain_short = 'ant'
145 
146 #elif (defined(ASF))
147 ch_domain_long = 'Austfonna'
148 ch_domain_short = 'asf'
149 
150 #elif (defined(EMTP2SGE))
151 ch_domain_long = 'EISMINT Phase 2 Simplified Geometry Experiment'
152 ch_domain_short = 'emtp2sge'
153 
154 #elif (defined(GRL))
155 ch_domain_long = 'Greenland'
156 ch_domain_short = 'grl'
157 
158 #elif (defined(NHEM))
159 ch_domain_long = 'Northern hemisphere'
160 ch_domain_short = 'nhem'
161 
162 #elif (defined(SCAND))
163 ch_domain_long = 'Scandinavia and Eurasia'
164 ch_domain_short = 'scand'
165 
166 #elif (defined(TIBET))
167 ch_domain_long = 'Tibet'
168 ch_domain_short = 'tibet'
169 
170 #elif (defined(NMARS))
171 ch_domain_long = 'North polar cap of Mars'
172 ch_domain_short = 'nmars'
173 
174 #elif (defined(SMARS))
175 ch_domain_long = 'South polar cap of Mars'
176 ch_domain_short = 'smars'
177 
178 #elif (defined(XYZ))
179 ch_domain_long = 'XYZ'
180 ch_domain_short = 'xyz'
181 #if (defined(HEINO))
182 ch_domain_long = trim(ch_domain_long)//'/ISMIP HEINO'
183 #elif (defined(MOCHO))
184 ch_domain_long = trim(ch_domain_long)//'/MOCHO'
185 #endif
186 
187 #else
188 errormsg = ' >>> sico_init: No valid domain specified!'
189 call error(errormsg)
190 #endif
191 
192 !-------- Some initial values --------
193 
194 n_output = 0
195 
196 dtime = 0.0_dp
197 dtime_temp = 0.0_dp
198 dtime_wss = 0.0_dp
199 dtime_out = 0.0_dp
200 dtime_ser = 0.0_dp
201 
202 time = 0.0_dp
203 time_init = 0.0_dp
204 time_end = 0.0_dp
205 time_output = 0.0_dp
206 
207 !-------- Initialisation of the Library of Iterative Solvers Lis,
208 ! if required --------
209 
210 #if (CALCTHK==3 || CALCTHK==6 || MARGIN==3 || DYNAMICS==2)
211  call lis_initialize(ierr)
212 #endif
213 
214 !-------- Read physical parameters --------
215 
216 call read_phys_para()
217 
218 call ice_mat_eqs_pars(rf, r_t, kappa, c, -190, 10)
219 
220 #if (!defined(SURFACE_FORCING) || SURFACE_FORCING==1)
221 
222 errormsg = ' >>> sico_init: Option SURFACE_FORCING=1 not supported any more!'
223 call error(errormsg)
224 
225 #elif (SURFACE_FORCING==2)
226 
227 temp_0 = temp_0 ! deg C
228 gamma_t = gamma_t *1.0e-03_dp ! K/km -> K/m
229 s_0 = s_0 /year_sec ! m/a -> m/s
230 m_0 = m_0 *1.0e-03_dp/year_sec ! m/(a*km) -> 1/s
231 ela = ela *1.0e+03_dp ! km -> m
232 #if (defined(DELA_DTS))
233  dela_dts = dela_dts ! m/degC
234 #else
235  dela_dts = 0.0_dp ! m/degC
236 #endif
237 
238 #elif (SURFACE_FORCING==3)
239 
240 temp_0 = temp_0 ! deg C
241 gamma_t = gamma_t *1.0e-03_dp ! K/km -> K/m
242 s_0 = s_0 /year_sec ! m/a -> m/s
243 m_0 = m_0 *1.0e-03_dp/year_sec ! m/(a*km) -> 1/s
244 ela = ela *1.0e+03_dp ! km -> m
245 ela_amp = ela_amp ! m
246 #if (defined(DELA_DTS))
247  dela_dts = dela_dts ! m/degC
248 #else
249  dela_dts = 0.0_dp ! m/degC
250 #endif
251 phi_0 = phi_0 *pi_180 ! deg -> rad
252 x_gip = x_gip ! m
253 y_gip = y_gip ! m
254 
255 #elif (SURFACE_FORCING==4)
256 
257 temp_0 = temp_0 ! deg C
258 gamma_t = gamma_t *1.0e-03_dp ! K/km -> K/m
259 s_0 = s_0 /year_sec ! m/a -> m/s
260 m_0 = m_0 *1.0e-03_dp/year_sec ! m/(a*km) -> 1/s
261 ela = ela *1.0e+03_dp ! km -> m
262 ela_amp = ela_amp ! m
263 #if (defined(DELA_DTS))
264  dela_dts = dela_dts ! m/degC
265 #else
266  dela_dts = 0.0_dp ! m/degC
267 #endif
268 phi_0 = phi_0 *pi_180 ! deg -> rad
269 x_gip = x_gip ! m
270 y_gip = y_gip ! m
271 m_1 = m_1 *1.0e-03_dp/year_sec ! m/(a*km) -> 1/s
272 z_gc = z_gc *1.0e+03_dp ! km -> m
273 
274 #elif (SURFACE_FORCING==5)
275 
276 temp_0 = temp_0 ! deg C
277 gamma_t = gamma_t *1.0e-03_dp ! K/km -> K/m
278 s_0 = s_0 /year_sec ! m/a -> m/s
279 m_0 = m_0 *1.0e-03_dp/year_sec ! m/(a*km) -> 1/s
280 ela = ela *1.0e+03_dp ! km -> m
281 ela_amp = ela_amp ! m
282 #if (defined(DELA_DTS))
283  dela_dts = dela_dts ! m/degC
284 #else
285  dela_dts = 0.0_dp ! m/degC
286 #endif
287 phi_0 = phi_0 *pi_180 ! deg -> rad
288 m_1 = m_1 *1.0e-03_dp/year_sec ! m/(a*km) -> 1/s
289 z_gc = z_gc *1.0e+03_dp ! km -> m
290 
291 #elif (SURFACE_FORCING==6)
292 
293 temp_0 = temp_0 ! deg C
294 gamma_t = gamma_t *1.0e-03_dp ! K/km -> K/m
295 s_0 = s_0 /year_sec ! m/a -> m/s
296 m_0 = m_0 *1.0e-03_dp/year_sec ! m/(a*km) -> 1/s
297 ela = ela *1.0e+03_dp ! km -> m
298 ela_amp = ela_amp ! m
299 #if (defined(DELA_DTS))
300  dela_dts = dela_dts ! m/degC
301 #else
302  dela_dts = 0.0_dp ! m/degC
303 #endif
304 phi_0 = phi_0 *pi_180 ! deg -> rad
305 tgt = tgt ! deg
306 
307 #elif (SURFACE_FORCING==7)
308 
309 temp_0 = temp_0 ! deg C
310 gamma_t = gamma_t *1.0e-03_dp ! K/km -> K/m
311 s_0 = s_0 /year_sec ! m/a -> m/s
312 m_0 = m_0 *1.0e-03_dp/year_sec ! m/(a*km) -> 1/s
313 ela = ela *1.0e+03_dp ! km -> m
314 ela_amp = ela_amp ! m
315 #if (defined(DELA_DTS))
316  dela_dts = dela_dts ! m/degC
317 #else
318  dela_dts = 0.0_dp ! m/degC
319 #endif
320 phi_0 = phi_0 *pi_180 ! deg -> rad
321 
322 #elif (SURFACE_FORCING==8)
323 
324 temp_0 = temp_0 ! deg C
325 gamma_t = gamma_t *1.0e-03_dp ! K/km -> K/m
326 s_0 = s_0 /year_sec ! m/a -> m/s
327 m_0 = m_0 *1.0e-03_dp/year_sec ! m/(a*km) -> 1/s
328 ela = ela *1.0e+03_dp ! km -> m
329 ela_amp = ela_amp ! m
330 #if (defined(DELA_DTS))
331  dela_dts = dela_dts ! m/degC
332 #else
333  dela_dts = 0.0_dp ! m/degC
334 #endif
335 phi_0 = phi_0 *pi_180 ! deg -> rad
336 x_gip = x_gip ! m
337 y_gip = y_gip ! m
338 x_gip2 = x_gip2 ! m
339 y_gip2 = y_gip2 ! m
340 ela_amp2 = ela_amp2 ! m
341 
342 
343 #else
344 
345 errormsg = ' >>> sico_init: SURFACE_FORCING must be between 1 and 8!'
346 call error(errormsg)
347 
348 #endif
349 
350 ! ------ Some auxiliary quantities required for the enthalpy method
351 
352 call calc_c_int_table(c, -190, 10, l)
354 
355 !-------- Compatibility check of the SICOPOLIS version with the header file
356 
357 if ( trim(version) /= trim(sico_version) ) then
358  errormsg = ' >>> sico_init: ' &
359  //'SICOPOLIS version not compatible with header file!'
360  call error(errormsg)
361 end if
362 
363 !-------- Check whether the dynamics and thermodynamics modes are defined
364 
365 #if (!defined(DYNAMICS))
366 errormsg = ' >>> sico_init: DYNAMICS not defined in the header file!'
367 call error(errormsg)
368 #endif
369 
370 #if (!defined(CALCMOD))
371 errormsg = ' >>> sico_init: CALCMOD not defined in the header file!'
372 call error(errormsg)
373 #endif
374 
375 #if (defined(ENTHMOD))
376 errormsg = ' >>> sico_init: ENTHMOD must not be defined any more.' &
377  // end_of_line &
378  //' Please update your header file!'
379 call error(errormsg)
380 #endif
381 
382 !-------- Compatibility check of the horizontal resolution with the
383 ! number of grid points and the position of the origin point --------
384 
385 #if (GRID==0)
386 
387 if (approx_equal(dx, 0.2_dp, eps_sp_dp)) then
388 
389  if ((imax /= 56).or.(jmax /= 48)) then
390  errormsg = ' >>> sico_init: IMAX and/or JMAX wrong!'
391  call error(errormsg)
392  end if
393 
394  if ( (.not.(approx_equal(x0, 748.147d0, eps_sp_dp))) &
395  .or. &
396  (.not.(approx_equal(y0, 5570.528d0, eps_sp_dp))) ) then
397  errormsg = ' >>> sico_init: X0 and/or Y0 wrong!'
398  call error(errormsg)
399  end if
400 
401 else if (approx_equal(dx, 0.1_dp, eps_sp_dp)) then
402 
403  if ((imax /= 114).or.(jmax /= 98)) then
404  errormsg = ' >>> sico_init: IMAX and/or JMAX wrong!'
405  call error(errormsg)
406  end if
407 
408  if ( (.not.(approx_equal(x0, 748.225d0, eps_sp_dp))) &
409  .or. &
410  (.not.(approx_equal(y0, 5570.327d0, eps_sp_dp))) ) then
411  errormsg = ' >>> sico_init: X0 and/or Y0 wrong!'
412  call error(errormsg)
413  end if
414 
415 else if (approx_equal(dx, 0.05_dp, eps_sp_dp)) then
416 
417  if ((imax /= 231).or.(jmax /= 199)) then
418  errormsg = ' >>> sico_init: IMAX and/or JMAX wrong!'
419  call error(errormsg)
420  end if
421 
422  if ( (.not.(approx_equal(x0, 748.086d0, eps_sp_dp))) &
423  .or. &
424  (.not.(approx_equal(y0, 5570.378d0, eps_sp_dp))) ) then
425  errormsg = ' >>> sico_init: X0 and/or Y0 wrong!'
426  call error(errormsg)
427  end if
428 
429 else
430 
431  errormsg = ' >>> sico_init: DX wrong!'
432  call error(errormsg)
433 
434 end if
435 
436 #elif (GRID==1)
437 
438 errormsg = ' >>> sico_init: GRID==1 not allowed for this application!'
439 call error(errormsg)
440 
441 #elif (GRID==2)
442 
443 errormsg = ' >>> sico_init: GRID==2 not allowed for this application!'
444 call error(errormsg)
445 
446 #endif
447 
448 !-------- Compatibility check of the thermodynamics mode
449 ! (cold vs. polythermal vs. enthalpy method)
450 ! and the number of grid points in the lower (kt) ice domain --------
451 
452 #if (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1)
453 
454 if (ktmax > 2) then
455  write(6, fmt='(a)') ' >>> sico_init: For options CALCMOD==0, 2, 3 or -1,'
456  write(6, fmt='(a)') ' the separate kt domain is redundant.'
457  write(6, fmt='(a)') ' Therefore, consider setting KTMAX to 2.'
458  write(6, fmt='(a)') ' '
459 end if
460 
461 #endif
462 
463 !-------- Compatibility check of discretization schemes for the horizontal and
464 ! vertical advection terms in the temperature and age equations --------
465 
466 #if (ADV_HOR==1)
467 errormsg = ' >>> sico_init: ' &
468  //'Option ADV_HOR==1 (central differences) not defined!'
469 call error(errormsg)
470 #endif
471 
472 !-------- Check whether for the shallow shelf
473 ! or shelfy stream approximation
474 ! the chosen grid is Cartesian coordinates
475 ! without distortion correction (GRID==0) --------
476 
477 #if ((MARGIN==3 && DYNAMICS==1) || DYNAMICS==2) /* requires SSA or SStA */
478 #if (GRID != 0)
479 errormsg = ' >>> sico_init: Distortion correction not implemented' &
480  // end_of_line &
481  //' for the shallow shelf approximation (SSA)' &
482  // end_of_line &
483  //' or the shelfy stream approximation (SStA)' &
484  // end_of_line &
485  //' -> GRID==0 required!'
486 call error(errormsg)
487 #endif
488 #endif
489 
490 !-------- Setting of forcing flag --------
491 
492 forcing_flag = 1 ! forcing by delta_ts
493 
494 !-------- Initialization of numerical time steps --------
495 
496 dtime0 = dtime0
497 dtime_temp0 = dtime_temp0
498 
499 !-------- Further initializations --------
500 
501 dzeta_c = 1.0_dp/real(kcmax,dp)
502 dzeta_t = 1.0_dp/real(ktmax,dp)
503 dzeta_r = 1.0_dp/real(krmax,dp)
504 
505 ndat2d = 1
506 ndat3d = 1
507 
508 !-------- General abbreviations --------
509 
510 ! ------ kc domain
511 
512 if (deform >= eps) then
513 
514  flag_aa_nonzero = .true. ! non-equidistant grid
515 
516  aa = deform
517  ea = exp(aa)
518 
519  kc=0
520  zeta_c(kc) = 0.0_dp
521  eaz_c(kc) = 1.0_dp
522  eaz_c_quotient(kc) = 0.0_dp
523 
524  do kc=1, kcmax-1
525  zeta_c(kc) = kc*dzeta_c
526  eaz_c(kc) = exp(aa*zeta_c(kc))
527  eaz_c_quotient(kc) = (eaz_c(kc)-1.0_dp)/(ea-1.0_dp)
528  end do
529 
530  kc=kcmax
531  zeta_c(kc) = 1.0_dp
532  eaz_c(kc) = exp(aa)
533  eaz_c_quotient(kc) = 1.0_dp
534 
535 else
536 
537  flag_aa_nonzero = .false. ! equidistant grid
538 
539  aa = 0.0_dp
540  ea = 1.0_dp
541 
542  kc=0
543  zeta_c(kc) = 0.0_dp
544  eaz_c(kc) = 1.0_dp
545  eaz_c_quotient(kc) = 0.0_dp
546 
547  do kc=1, kcmax-1
548  zeta_c(kc) = kc*dzeta_c
549  eaz_c(kc) = 1.0_dp
550  eaz_c_quotient(kc) = zeta_c(kc)
551  end do
552 
553  kc=kcmax
554  zeta_c(kc) = 1.0_dp
555  eaz_c(kc) = 1.0_dp
556  eaz_c_quotient(kc) = 1.0_dp
557 
558 end if
559 
560 ! ------ kt domain
561 
562 kt=0
563 zeta_t(kt) = 0.0_dp
564 
565 do kt=1, ktmax-1
566  zeta_t(kt) = kt*dzeta_t
567 end do
568 
569 kt=ktmax
570 zeta_t(kt) = 1.0_dp
571 
572 ! ------ kr domain
573 
574 kr=0
575 zeta_r(kr) = 0.0_dp
576 
577 do kr=1, krmax-1
578  zeta_r(kr) = kr*dzeta_r
579 end do
580 
581 kr=krmax
582 zeta_r(kr) = 1.0_dp
583 
584 !-------- Construction of a vector (with index n) from a 2-d array
585 ! (with indices i, j) by diagonal numbering --------
586 
587 n=1
588 
589 do m=0, imax+jmax
590  do i=m, 0, -1
591  j = m-i
592  if ((i <= imax).and.(j <= jmax)) then
593  ii(n) = i
594  jj(n) = j
595  nn(j,i) = n
596  n=n+1
597  end if
598  end do
599 end do
600 
601 !-------- Specification of current simulation --------
602 
603 runname = runname
604 anfdatname = anfdatname
605 
606 #if (defined(YEAR_ZERO))
608 #else
609 year_zero = 2000.0_dp ! default value 2000 CE
610 #endif
611 
612 time_init0 = time_init0
613 time_end0 = time_end0
614 dtime_ser0 = dtime_ser0
615 
616 #if (OUTPUT==1 || OUTPUT==3)
617 dtime_out0 = dtime_out0
618 #endif
619 
620 #if (OUTPUT==2 || OUTPUT==3)
621 n_output = n_output
622 time_output0( 1) = time_out0_01
623 time_output0( 2) = time_out0_02
624 time_output0( 3) = time_out0_03
625 time_output0( 4) = time_out0_04
626 time_output0( 5) = time_out0_05
627 time_output0( 6) = time_out0_06
628 time_output0( 7) = time_out0_07
629 time_output0( 8) = time_out0_08
630 time_output0( 9) = time_out0_09
631 time_output0(10) = time_out0_10
632 time_output0(11) = time_out0_11
633 time_output0(12) = time_out0_12
634 time_output0(13) = time_out0_13
635 time_output0(14) = time_out0_14
636 time_output0(15) = time_out0_15
637 time_output0(16) = time_out0_16
638 time_output0(17) = time_out0_17
639 time_output0(18) = time_out0_18
640 time_output0(19) = time_out0_19
641 time_output0(20) = time_out0_20
642 #endif
643 
644 !-------- Write log file --------
645 
646 shell_command = 'if [ ! -d'
647 shell_command = trim(shell_command)//' '//outpath
648 shell_command = trim(shell_command)//' '//'] ; then mkdir'
649 shell_command = trim(shell_command)//' '//outpath
650 shell_command = trim(shell_command)//' '//'; fi'
651 call system(trim(shell_command))
652  ! Check whether directory OUTPATH exists. If not, it is created.
653 
654 filename_with_path = trim(outpath)//'/'//trim(runname)//'.log'
655 
656 open(10, iostat=ios, file=trim(filename_with_path), status='new')
657 
658 if (ios /= 0) then
659  errormsg = ' >>> sico_init: Error when opening the log file!'
660  call error(errormsg)
661 end if
662 
663 write(10, fmt=trim(fmt1)) 'Computational domain:'
664 write(10, fmt=trim(fmt1)) trim(ch_domain_long)
665 write(10, fmt=trim(fmt1)) ' '
666 
667 write(10, fmt=trim(fmt2a)) 'GRID = ', grid
668 write(10, fmt=trim(fmt1)) ' '
669 
670 write(10, fmt=trim(fmt2)) 'imax =', imax
671 write(10, fmt=trim(fmt2)) 'jmax =', jmax
672 write(10, fmt=trim(fmt2)) 'kcmax =', kcmax
673 write(10, fmt=trim(fmt2)) 'ktmax =', ktmax
674 write(10, fmt=trim(fmt2)) 'krmax =', krmax
675 write(10, fmt=trim(fmt1)) ' '
676 
677 write(10, fmt=trim(fmt3)) 'a =', aa
678 write(10, fmt=trim(fmt1)) ' '
679 
680 #if (GRID==0 || GRID==1)
681 write(10, fmt=trim(fmt3)) 'x0 =', x0
682 write(10, fmt=trim(fmt3)) 'y0 =', y0
683 write(10, fmt=trim(fmt3)) 'dx =', dx
684 #elif (GRID==2)
685 write(10, fmt=trim(fmt3)) 'dlambda =', dlambda
686 write(10, fmt=trim(fmt3)) 'dphi =', dphi
687 #endif
688 write(10, fmt=trim(fmt1)) ' '
689 
690 write(10, fmt=trim(fmt3)) 'year_zero =', year_zero
691 write(10, fmt=trim(fmt3)) 'time_init =', time_init0
692 write(10, fmt=trim(fmt3)) 'time_end =', time_end0
693 write(10, fmt=trim(fmt3)) 'dtime =', dtime0
694 write(10, fmt=trim(fmt3)) 'dtime_temp =', dtime_temp0
695 write(10, fmt=trim(fmt1)) ' '
696 
697 write(10, fmt=trim(fmt2)) 'ANF_DAT =', anf_dat
698 write(10, fmt=trim(fmt1)) 'zs_present file = '//zs_present_file
699 #if (ANF_DAT==1)
700 #if (defined(ZB_PRESENT_FILE))
701 write(10, fmt=trim(fmt1)) 'zb_present file = '//zb_present_file
702 #endif
703 write(10, fmt=trim(fmt1)) 'zl_present file = '//zl_present_file
704 #endif
705 write(10, fmt=trim(fmt1)) 'zl0 file = '//zl0_file
706 #if (defined(MASK_PRESENT_FILE))
707 write(10, fmt=trim(fmt1)) 'mask_present file = '//mask_present_file
708 #endif
709 #if (ANF_DAT==1 && defined(TEMP_INIT))
710 write(10, fmt=trim(fmt2)) 'TEMP_INIT =', temp_init
711 #endif
712 #if (ANF_DAT==3)
713 write(10, fmt=trim(fmt1)) 'Initial-value file = '//anfdatname
714 write(10, fmt=trim(fmt1)) 'Path to initial-value file = '//anfdatpath
715 #endif
716 write(10, fmt=trim(fmt1)) ' '
717 
718 #if (THK_EVOL==2)
719 write(10, fmt=trim(fmt3)) 'time_target_topo_init =', time_target_topo_init0
720 write(10, fmt=trim(fmt3)) 'time_target_topo_final =', time_target_topo_final0
721 write(10, fmt=trim(fmt1)) 'Target-topography file = '//target_topo_dat_name
722 write(10, fmt=trim(fmt1)) 'Path to target-topography file = '//target_topo_dat_path
723 write(10, fmt=trim(fmt1)) ' '
724 #endif
725 
726 #if (THK_EVOL==3)
727 write(10, fmt=trim(fmt3)) 'target_topo_tau =', target_topo_tau0
728 write(10, fmt=trim(fmt1)) 'Target-topography file = '//target_topo_dat_name
729 write(10, fmt=trim(fmt1)) 'Path to target-topography file = '//target_topo_dat_path
730 write(10, fmt=trim(fmt1)) ' '
731 #endif
732 
733 #if (THK_EVOL==4)
734 write(10, fmt=trim(fmt1)) 'Maximum ice extent mask file = '//mask_maxextent_file
735 write(10, fmt=trim(fmt1)) ' '
736 #endif
737 
738 write(10, fmt=trim(fmt1)) 'Physical-parameter file = '//phys_para_file
739 write(10, fmt=trim(fmt1)) ' '
740 
741 #if (CALCTHK==2 || CALCTHK==3 || CALCTHK==5 || CALCTHK==6)
742 write(10, fmt=trim(fmt3)) 'ovi_weight =', ovi_weight
743 #if (CALCTHK==2 || CALCTHK==5)
744 write(10, fmt=trim(fmt3)) 'omega_sor =', omega_sor
745 #if (ITER_MAX_SOR>0)
746 write(10, fmt=trim(fmt2a)) 'iter_max_sor = ', iter_max_sor
747 #endif
748 #endif
749 write(10, fmt=trim(fmt1)) ' '
750 #endif
751 
752 #if (defined(SURFACE_FORCING))
753 write(10, fmt=trim(fmt2)) 'SURFACE_FORCING =', surface_forcing
754 #endif
755 #if (!defined(SURFACE_FORCING) || SURFACE_FORCING==1)
756 errormsg = ' >>> sico_init: Option SURFACE_FORCING=1 not supported any more!'
757 call error(errormsg)
758 #elif (SURFACE_FORCING==2)
759 write(10, fmt=trim(fmt3)) 'temp_0 =', temp_0
760 write(10, fmt=trim(fmt3)) 'gamma_t =', gamma_t
761 write(10, fmt=trim(fmt3)) 's_0 =', s_0
762 write(10, fmt=trim(fmt3)) 'm_0 =', m_0
763 write(10, fmt=trim(fmt3)) 'ela =', ela
764 #endif
765 write(10, fmt=trim(fmt1)) ' '
766 
767 #if (defined(SMB_CORR_FILE))
768 if ( trim(adjustl(smb_corr_file)) /= 'none' ) then
769  write(10, fmt=trim(fmt1)) 'smb_corr_file = '//smb_corr_file
770  write(10, fmt=trim(fmt1)) ' '
771 end if
772 #endif
773 
774 #if (TSURFACE==1)
775 write(10, fmt=trim(fmt3)) 'delta_ts0 =', delta_ts0
776 #elif (TSURFACE==3)
777 write(10, fmt=trim(fmt3)) 'sine_amplit =', sine_amplit
778 write(10, fmt=trim(fmt3)) 'sine_period =', sine_period
779 #elif (TSURFACE==4)
780 write(10, fmt=trim(fmt1)) 'GRIP file = '//grip_temp_file
781 write(10, fmt=trim(fmt3)) 'grip_temp_fact =', grip_temp_fact
782 #endif
783 
784 write(10, fmt=trim(fmt2a)) 'SEA_LEVEL = ', sea_level
785 #if (SEA_LEVEL==1)
786 write(10, fmt=trim(fmt3)) 'z_sl0 =', z_sl0
787 #elif (SEA_LEVEL==3)
788 write(10, fmt=trim(fmt1)) 'sea-level file = '//sea_level_file
789 #endif
790 write(10, fmt=trim(fmt1)) ' '
791 
792 #if (MARGIN==2)
793 #if (MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3)
794 write(10, fmt=trim(fmt3)) 'z_mar =', z_mar
795 write(10, fmt=trim(fmt1)) ' '
796 #elif (MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \
797  || marine_ice_calving==6 || marine_ice_calving==7)
798 write(10, fmt=trim(fmt3)) 'fact_z_mar =', fact_z_mar
799 write(10, fmt=trim(fmt1)) ' '
800 #elif (MARINE_ICE_FORMATION==2 && MARINE_ICE_CALVING==9)
801 write(10, fmt=trim(fmt3)) 'calv_uw_coeff =', calv_uw_coeff
802 write(10, fmt=trim(fmt3)) 'r1_calv_uw =', r1_calv_uw
803 write(10, fmt=trim(fmt3)) 'r2_calv_uw =', r2_calv_uw
804 write(10, fmt=trim(fmt1)) ' '
805 #endif
806 #elif (MARGIN==3)
807 #if (ICE_SHELF_CALVING==2)
808 write(10, fmt=trim(fmt3)) 'H_calv =', h_calv
809 write(10, fmt=trim(fmt1)) ' '
810 #endif
811 #endif
812 
813 #if (defined(BASAL_HYDROLOGY))
814 write(10, fmt=trim(fmt2a)) 'BASAL_HYDROLOGY = ', basal_hydrology
815 #endif
816 write(10, fmt=trim(fmt2a)) 'SLIDE_LAW = ', slide_law
817 write(10, fmt=trim(fmt3)) 'c_slide =', c_slide
818 write(10, fmt=trim(fmt3)) 'gamma_slide =', gamma_slide
819 write(10, fmt=trim(fmt2a)) 'p_weert = ', p_weert
820 write(10, fmt=trim(fmt2a)) 'q_weert = ', q_weert
821 #if (defined(TIME_RAMP_UP_SLIDE))
822 write(10, fmt=trim(fmt3)) 'time_ramp_up_slide =', time_ramp_up_slide
823 #endif
824 #if (SLIDE_LAW==2 || SLIDE_LAW==3)
825 write(10, fmt=trim(fmt3)) 'red_pres_limit_fact =', red_pres_limit_fact
826 #endif
827 #if (BASAL_HYDROLOGY==1 \
828  && defined(hydro_slide_sat_fct) \
829  && defined(c_hw_slide) && defined(hw0_slide))
830 write(10, fmt=trim(fmt2a)) 'HYDRO_SLIDE_SAT_FCT = ', hydro_slide_sat_fct
831 write(10, fmt=trim(fmt3)) 'c_Hw_slide =', c_hw_slide
832 write(10, fmt=trim(fmt3)) 'Hw0_slide =', hw0_slide
833 #endif
834 write(10, fmt=trim(fmt1)) ' '
835 
836 write(10, fmt=trim(fmt2a)) 'Q_GEO_MOD = ', q_geo_mod
837 #if (Q_GEO_MOD==1)
838 write(10, fmt=trim(fmt3)) 'q_geo =', q_geo
839 #endif
840 write(10, fmt=trim(fmt1)) ' '
841 
842 write(10, fmt=trim(fmt2)) 'REBOUND =', rebound
843 #if (REBOUND==1)
844 write(10, fmt=trim(fmt3)) 'frac_llra =', frac_llra
845 #endif
846 #if (REBOUND==1 || REBOUND==2)
847 write(10, fmt=trim(fmt2)) 'TIME_LAG_MOD =', time_lag_mod
848 #if (TIME_LAG_MOD==1)
849 write(10, fmt=trim(fmt3)) 'time_lag =', time_lag
850 #elif (TIME_LAG_MOD==2)
851 write(10, fmt=trim(fmt1)) 'time_lag_file = '//time_lag_file
852 #else
853 errormsg = ' >>> sico_init: TIME_LAG_MOD must be either 1 or 2!'
854 call error(errormsg)
855 #endif
856 #endif
857 #if (REBOUND==2)
858 write(10, fmt=trim(fmt2)) 'FLEX_RIG_MOD =', flex_rig_mod
859 #if (FLEX_RIG_MOD==1)
860 write(10, fmt=trim(fmt3)) 'flex_rig =', flex_rig
861 #elif (FLEX_RIG_MOD==2)
862 write(10, fmt=trim(fmt1)) 'flex_rig_file = '//flex_rig_file
863 #else
864 errormsg = ' >>> sico_init: FLEX_RIG_MOD must be either 1 or 2!'
865 call error(errormsg)
866 #endif
867 #endif
868 write(10, fmt=trim(fmt2)) 'Q_LITHO =', q_litho
869 write(10, fmt=trim(fmt1)) ' '
870 
871 #if (FLOW_LAW==2)
872 write(10, fmt=trim(fmt3)) 'gr_size =', gr_size
873 write(10, fmt=trim(fmt1)) ' '
874 #endif
875 #if (FIN_VISC==2)
876 write(10, fmt=trim(fmt3)) 'sigma_res =', sigma_res
877 write(10, fmt=trim(fmt1)) ' '
878 #endif
879 
880 write(10, fmt=trim(fmt2)) 'ENHMOD =', enhmod
881 #if (ENHMOD==1 || ENHMOD==2 || ENHMOD==3)
882 write(10, fmt=trim(fmt3)) 'enh_fact =', enh_fact
883 #endif
884 #if (ENHMOD==2 || ENHMOD==3)
885 write(10, fmt=trim(fmt3)) 'enh_intg =', enh_intg
886 #endif
887 #if (ENHMOD==2)
888 write(10, fmt=trim(fmt3)) 'age_trans =', age_trans_0
889 #endif
890 #if (ENHMOD==3)
891 write(10, fmt=trim(fmt3)) 'date_trans1 =', date_trans1_0
892 write(10, fmt=trim(fmt3)) 'date_trans2 =', date_trans2_0
893 write(10, fmt=trim(fmt3)) 'date_trans3 =', date_trans3_0
894 #endif
895 #if (ENHMOD==4 || ENHMOD==5)
896 write(10, fmt=trim(fmt3)) 'enh_compr =', enh_compr
897 write(10, fmt=trim(fmt3)) 'enh_shear =', enh_shear
898 #endif
899 #if ((ENHMOD==1 || ENHMOD==2 || ENHMOD==3 || ENHMOD==4) && MARGIN==3)
900 write(10, fmt=trim(fmt3)) 'enh_shelf =', enh_shelf
901 #endif
902 write(10, fmt=trim(fmt1)) ' '
903 
904 write(10, fmt=trim(fmt3)) 'numdiff_H_t =', numdiff_h_t
905 write(10, fmt=trim(fmt3)) 'tau_cts =', tau_cts
906 write(10, fmt=trim(fmt3)) 'vh_max =', vh_max
907 write(10, fmt=trim(fmt3)) 'hd_min =', hd_min
908 write(10, fmt=trim(fmt3)) 'hd_max =', hd_max
909 #if (defined(QBM_MIN))
910 write(10, fmt=trim(fmt3)) 'qbm_min =', qbm_min
911 #elif (defined(QB_MIN)) /* obsolete */
912 write(10, fmt=trim(fmt3)) 'qb_min =', qb_min
913 #endif
914 #if (defined(QBM_MAX))
915 write(10, fmt=trim(fmt3)) 'qbm_max =', qbm_max
916 #elif (defined(QB_MAX)) /* obsolete */
917 write(10, fmt=trim(fmt3)) 'qb_max =', qb_max
918 #endif
919 write(10, fmt=trim(fmt3)) 'age_min =', age_min
920 write(10, fmt=trim(fmt3)) 'age_max =', age_max
921 write(10, fmt=trim(fmt3)) 'mean_accum =', mean_accum
922 #if (ADV_VERT==1)
923 write(10, fmt=trim(fmt3)) 'age_diff =', agediff
924 #endif
925 write(10, fmt=trim(fmt1)) ' '
926 
927 write(10, fmt=trim(fmt2a)) 'DYNAMICS = ', dynamics
928 #if ((DYNAMICS==1 && MARGIN==3) || DYNAMICS==2)
929 #if (defined(LIS_OPTS))
930 write(10, fmt=trim(fmt1)) 'lis_opts = '//lis_opts
931 #endif
932 #if (defined(N_ITER_SSA))
933 write(10, fmt=trim(fmt2a)) 'n_iter_ssa = ', n_iter_ssa
934 #endif
935 #if (defined(RELAX_FACT_SSA))
936 write(10, fmt=trim(fmt3)) 'relax_fact_ssa =', relax_fact_ssa
937 #endif
938 #endif
939 #if (DYNAMICS==2 && defined(RATIO_SL_THRESH))
940 write(10, fmt=trim(fmt3)) 'ratio_sl_thresh =', ratio_sl_thresh
941 #endif
942 #if (DYNAMICS==2 && defined(SSTA_SIA_WEIGH_FCT))
943 write(10, fmt=trim(fmt2a)) 'SSTA_SIA_WEIGH_FCT = ', ssta_sia_weigh_fct
944 #endif
945 write(10, fmt=trim(fmt1)) ' '
946 
947 write(10, fmt=trim(fmt2)) 'CALCMOD =', calcmod
948 #if (CALCMOD==-1 && defined(TEMP_CONST))
949 write(10, fmt=trim(fmt3)) 'TEMP_CONST =', temp_const
950 #endif
951 #if (CALCMOD==-1 && defined(AGE_CONST))
952 write(10, fmt=trim(fmt3)) 'AGE_CONST =', age_const
953 #endif
954 #if (CALCMOD==1 && defined(CTS_MELTING_FREEZING))
955 write(10, fmt=trim(fmt2)) 'CTS_MELTING_FREEZING =', cts_melting_freezing
956 #endif
957 write(10, fmt=trim(fmt2)) 'FLOW_LAW =', flow_law
958 write(10, fmt=trim(fmt2)) 'FIN_VISC =', fin_visc
959 write(10, fmt=trim(fmt2)) 'MARGIN =', margin
960 #if (MARGIN==2)
961 write(10, fmt=trim(fmt2)) 'MARINE_ICE_FORMATION =', marine_ice_formation
962 write(10, fmt=trim(fmt2)) 'MARINE_ICE_CALVING =', marine_ice_calving
963 #elif (MARGIN==3)
964 write(10, fmt=trim(fmt2)) 'ICE_SHELF_CALVING =', ice_shelf_calving
965 #endif
966 #if (defined(THK_EVOL))
967 write(10, fmt=trim(fmt2)) 'THK_EVOL =', thk_evol
968 #else
969 errormsg = ' >>> sico_init: Define THK_EVOL in header file!'
970 call error(errormsg)
971 #endif
972 #if (defined(CALCTHK))
973 write(10, fmt=trim(fmt2)) 'CALCTHK =', calcthk
974 #else
975 errormsg = ' >>> sico_init: Define CALCTHK in header file!'
976 call error(errormsg)
977 #endif
978 write(10, fmt=trim(fmt2)) 'ADV_HOR =', adv_hor
979 write(10, fmt=trim(fmt2)) 'ADV_VERT =', adv_vert
980 write(10, fmt=trim(fmt2)) 'TOPOGRAD =', topograd
981 write(10, fmt=trim(fmt2)) 'TSURFACE =', tsurface
982 #if (defined(MB_ACCOUNT))
983 write(10, fmt=trim(fmt2a)) 'MB_ACCOUNT = ', mb_account
984 #endif
985 write(10, fmt=trim(fmt1)) ' '
986 
987 #if (defined(OUT_TIMES))
988 write(10, fmt=trim(fmt2a)) 'OUT_TIMES = ', out_times
989 #endif
990 #if (defined(OUTPUT_INIT))
991 write(10, fmt=trim(fmt2a)) 'OUTPUT_INIT = ', output_init
992 #endif
993 write(10, fmt=trim(fmt2a)) 'OUTPUT = ', output
994 #if (OUTPUT==1 || OUTPUT==3)
995 write(10, fmt=trim(fmt3)) 'dtime_out =' , dtime_out0
996 #endif
997 write(10, fmt=trim(fmt3)) 'dtime_ser =' , dtime_ser0
998 #if (OUTPUT==1 || OUTPUT==2)
999 write(10, fmt=trim(fmt2a)) 'ERGDAT = ', ergdat
1000 #endif
1001 write(10, fmt=trim(fmt2a)) 'NETCDF = ', netcdf
1002 #if (OUTPUT==2 || OUTPUT==3)
1003 write(10, fmt=trim(fmt2a)) 'n_output = ', n_output
1004 do n=1, n_output
1005  write(10, fmt=trim(fmt3)) 'time_output =' , time_output0(n)
1006 end do
1007 #endif
1008 write(10, fmt=trim(fmt1)) ' '
1009 
1010 write(10, fmt=trim(fmt1)) 'Program version and date: '//version//' / '//date
1011 
1012 close(10, status='keep')
1013 
1014 !-------- Conversion of time quantities --------
1015 
1016 year_zero = year_zero*year_sec ! a --> s
1017 time_init = time_init0*year_sec ! a --> s
1018 time_end = time_end0*year_sec ! a --> s
1019 dtime = dtime0*year_sec ! a --> s
1020 dtime_temp = dtime_temp0*year_sec ! a --> s
1021 dtime_ser = dtime_ser0*year_sec ! a --> s
1022 #if (OUTPUT==1 || OUTPUT==3)
1023 dtime_out = dtime_out0*year_sec ! a --> s
1024 #endif
1025 #if (OUTPUT==2 || OUTPUT==3)
1026 do n=1, n_output
1027  time_output(n) = time_output0(n)*year_sec ! a --> s
1028 end do
1029 #endif
1030 
1031 if (.not.approx_integer_multiple(dtime_temp, dtime, eps_sp_dp)) then
1032  errormsg = ' >>> sico_init: dtime_temp must be a multiple of dtime!'
1033  call error(errormsg)
1034 end if
1035 
1036 if (.not.approx_integer_multiple(dtime_ser, dtime, eps_sp_dp)) then
1037  errormsg = ' >>> sico_init: dtime_ser must be a multiple of dtime!'
1038  call error(errormsg)
1039 end if
1040 
1041 #if (OUTPUT==1 || OUTPUT==3)
1042 if (.not.approx_integer_multiple(dtime_out, dtime, eps_sp_dp)) then
1043  errormsg = ' >>> sico_init: dtime_out must be a multiple of dtime!'
1044  call error(errormsg)
1045 end if
1046 #endif
1047 
1048 #if (THK_EVOL==2)
1049 time_target_topo_init = time_target_topo_init0 *year_sec ! a --> s
1050 time_target_topo_final = time_target_topo_final0*year_sec ! a --> s
1051 #endif
1052 
1053 #if (THK_EVOL==3)
1054 target_topo_tau = target_topo_tau0 *year_sec ! a --> s
1055 #endif
1056 
1057 time = time_init
1058 
1059 !-------- Mean accumulation --------
1060 
1061 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho)
1062 ! ! mm/a water equiv. --> m/s ice equiv.
1063 
1064 !-------- Reading of the prescribed target topography --------
1065 
1066 #if ((THK_EVOL==2) || (THK_EVOL==3))
1067 
1068 target_topo_dat_name = trim(target_topo_dat_name)
1069 
1070 #if (NETCDF==1)
1071 errormsg = ' >>> sico_init: Reading of target topography' &
1072  // end_of_line &
1073  //' only implemented for NetCDF files (NETCDF==2)!'
1074 call error(errormsg)
1075 #elif (NETCDF==2)
1076 call read_target_topo_nc(target_topo_dat_name)
1077 #else
1078 errormsg = ' >>> sico_init: Parameter NETCDF must be either 1 or 2!'
1079 call error(errormsg)
1080 #endif
1081 
1082 #endif
1083 
1084 !-------- Reading of the maximum ice extent mask --------
1085 
1086 #if (THK_EVOL==4)
1087 
1088 #if (GRID==0 || GRID==1)
1089 
1090 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1091  trim(mask_maxextent_file)
1092 
1093 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
1094 
1095 #elif (GRID==2)
1096 
1097 errormsg = ' >>> sico_init: GRID==2 not allowed for this application!'
1098 call error(errormsg)
1099 
1100 #endif
1101 
1102 if (ios /= 0) then
1103  errormsg = ' >>> sico_init: ' &
1104  //'Error when opening the maximum ice extent mask file!'
1105  call error(errormsg)
1106 end if
1107 
1108 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1109 
1110 do j=jmax, 0, -1
1111  read(24, fmt=trim(fmt4)) (maske_maxextent(j,i), i=0,imax)
1112 end do
1113 
1114 close(24, status='keep')
1115 
1116 #endif
1117 
1118 !-------- Read data for delta_ts --------
1119 
1120 #if (TSURFACE==4)
1121 
1122 filename_with_path = trim(inpath)//'/general/'//trim(grip_temp_file)
1123 
1124 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1125 
1126 if (ios /= 0) then
1127  errormsg = ' >>> sico_init: Error when opening the data file for delta_ts!'
1128  call error(errormsg)
1129 end if
1130 
1131 read(21, fmt=*) ch_dummy, grip_time_min, grip_time_stp, grip_time_max
1132 
1133 if (ch_dummy /= '#') then
1134  errormsg = ' >>> sico_init: grip_time_min, grip_time_stp, grip_time_max' &
1135  // end_of_line &
1136  //' not defined in data file for delta_ts!'
1137  call error(errormsg)
1138 end if
1139 
1141 
1142 allocate(griptemp(0:ndata_grip))
1143 
1144 do n=0, ndata_grip
1145  read(21, fmt=*) d_dummy, griptemp(n)
1146 end do
1147 
1148 close(21, status='keep')
1149 
1150 #endif
1151 
1152 !-------- Prescribed surface mass balance correction --------
1153 
1154 smb_corr_in = 0.0_dp
1155 
1156 #if (defined(SMB_CORR_FILE))
1157 
1158 if (trim(adjustl(smb_corr_file)) /= 'none') then
1159 
1160  filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1161  trim(smb_corr_file)
1162 
1163  n = len_trim(filename_with_path)
1164  ch_nc_test = filename_with_path(n-2:n)
1165 
1166  if (ch_nc_test == '.nc') then ! NetCDF file
1167 
1168 #if (NETCDF==1)
1169 
1170  errormsg = ' >>> sico_init: Reading of NetCDF SMB_CORR_FILE' &
1171  // end_of_line &
1172  //' requires NETCDF==2!'
1173  call error(errormsg)
1174 
1175 #elif (NETCDF==2)
1176 
1177  call check( nf90_open(trim(filename_with_path), nf90_nowrite, ncid), &
1178  thisroutine )
1179 
1180  call check( nf90_inq_varid(ncid, 'DSMB', ncv) )
1181  call check( nf90_get_var(ncid, ncv, smb_corr_in_conv) )
1182 
1183  call check( nf90_close(ncid) )
1184 
1185  do i=0, imax
1186  do j=0, jmax
1187  smb_corr_in(j,i) = smb_corr_in_conv(i,j) /year_sec
1188  ! m/a ice equiv. -> m/s ice equiv.
1189  end do
1190  end do
1191 
1192 #else
1193  errormsg = ' >>> sico_init: Parameter NETCDF must be either 1 or 2!'
1194  call error(errormsg)
1195 #endif
1196 
1197  else ! ASCII file
1198 
1199  open(21, iostat=ios, &
1200  file=trim(filename_with_path), recl=rcl1, status='old')
1201 
1202  if (ios /= 0) then
1203  errormsg = ' >>> sico_init: ' &
1204  //'Error when opening the SMB correction file!'
1205  call error(errormsg)
1206  end if
1207 
1208  do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1209 
1210  do j=jmax, 0, -1
1211  read(21, fmt=*) (smb_corr_in(j,i), i=0,imax)
1212  end do
1213 
1214  close(21, status='keep')
1215 
1216  smb_corr_in = smb_corr_in /year_sec
1217  ! m/a ice equiv. -> m/s ice equiv.
1218 
1219  end if
1220 
1221 end if
1222 
1223 #endif
1224 
1225 !-------- Read data for z_sl --------
1226 
1227 #if (SEA_LEVEL==3)
1228 
1229 filename_with_path = trim(inpath)//'/general/'//trim(sea_level_file)
1230 
1231 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1232 
1233 if (ios /= 0) then
1234  errormsg = ' >>> sico_init: Error when opening the data file for z_sl!'
1235  call error(errormsg)
1236 end if
1237 
1238 read(21, fmt=*) ch_dummy, specmap_time_min, specmap_time_stp, specmap_time_max
1239 
1240 if (ch_dummy /= '#') then
1241  errormsg = ' >>> sico_init:' &
1242  // end_of_line &
1243  //' specmap_time_min, specmap_time_stp, specmap_time_max' &
1244  // end_of_line &
1245  //' not defined in data file for z_sl!'
1246  call error(errormsg)
1247 end if
1248 
1250 
1251 allocate(specmap_zsl(0:ndata_specmap))
1252 
1253 do n=0, ndata_specmap
1254  read(21, fmt=*) d_dummy, specmap_zsl(n)
1255 end do
1256 
1257 close(21, status='keep')
1258 
1259 #endif
1260 
1261 !-------- Determination of the geothermal heat flux --------
1262 
1263 #if (Q_GEO_MOD==1)
1264 
1265 ! ------ Constant value
1266 
1267 do i=0, imax
1268 do j=0, jmax
1269  q_geo(j,i) = q_geo *1.0e-03_dp ! mW/m2 -> W/m2
1270 end do
1271 end do
1272 
1273 #elif (Q_GEO_MOD==2)
1274 
1275 errormsg = ' >>> sico_init: ' &
1276  //'Option Q_GEO_MOD==2 not available for this application!'
1277 call error(errormsg)
1278 
1279 #endif
1280 
1281 !-------- Determination of the time lag
1282 ! of the relaxing asthenosphere --------
1283 
1284 #if (REBOUND==1 || REBOUND==2)
1285 
1286 #if (TIME_LAG_MOD==1)
1287 
1288 time_lag_asth = time_lag*year_sec ! a -> s
1289 
1290 #elif (TIME_LAG_MOD==2)
1291 
1292 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1293  trim(time_lag_file)
1294 
1295 open(29, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1296 
1297 if (ios /= 0) then
1298  errormsg = ' >>> sico_init: Error when opening the time-lag file!'
1299  call error(errormsg)
1300 end if
1301 
1302 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1303 
1304 do j=jmax, 0, -1
1305  read(29, fmt=*) (time_lag_asth(j,i), i=0,imax)
1306 end do
1307 
1308 close(29, status='keep')
1309 
1310 time_lag_asth = time_lag_asth*year_sec ! a -> s
1311 
1312 #endif
1313 
1314 #elif (REBOUND==0)
1315 
1316 time_lag_asth = 0.0_dp ! dummy values
1317 
1318 #endif
1319 
1320 !-------- Determination of the flexural rigidity of the lithosphere --------
1321 
1322 #if (REBOUND==2)
1323 
1324 #if (FLEX_RIG_MOD==1)
1325 
1326 flex_rig_lith = flex_rig
1327 
1328 #elif (FLEX_RIG_MOD==2)
1329 
1330 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1331  trim(flex_rig_file)
1332 
1333 open(29, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1334 
1335 if (ios /= 0) then
1336  errormsg = ' >>> sico_init: Error when opening the flex-rig file!'
1337  call error(errormsg)
1338 end if
1339 
1340 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1341 
1342 do j=jmax, 0, -1
1343  read(29, fmt=*) (flex_rig_lith(j,i), i=0,imax)
1344 end do
1345 
1346 close(29, status='keep')
1347 
1348 #endif
1349 
1350 #elif (REBOUND==0 || REBOUND==1)
1351 
1352 flex_rig_lith = 0.0_dp ! dummy values
1353 
1354 #endif
1355 
1356 !-------- Definition of initial values --------
1357 
1358 ! ------ Present topography
1359 
1360 #if (ANF_DAT==1)
1361 
1362 call topography1(dxi, deta)
1363 
1364 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1365 
1366 call boundary(time_init, dtime, dxi, deta, &
1367  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1368 
1369 where ((maske==0_i1b).or.(maske==3_i1b))
1370  ! grounded or floating ice
1372 elsewhere ! maske==1_i1b or 2_i1b, ice-free land or sea
1373  as_perp_apl = 0.0_dp
1374 end where
1375 
1376 smb_corr = 0.0_dp
1377 
1378 q_bm = 0.0_dp
1379 q_tld = 0.0_dp
1380 q_b_tot = 0.0_dp
1381 
1382 p_b_w = 0.0_dp
1383 h_w = 0.0_dp
1384 
1385 #if (!defined(TEMP_INIT) || TEMP_INIT==1)
1387 #elif (TEMP_INIT==2)
1389 #elif (TEMP_INIT==3)
1391 #elif (TEMP_INIT==4)
1393 #else
1394  errormsg = ' >>> sico_init: TEMP_INIT must be set to either 1, 2, 3 or 4!'
1395  call error(errormsg)
1396 #endif
1397 
1398 #if (ENHMOD==1)
1399  call calc_enhance_1()
1400 #elif (ENHMOD==2)
1401  call calc_enhance_2()
1402 #elif (ENHMOD==3)
1403  call calc_enhance_3(time_init)
1404 #elif (ENHMOD==4)
1405  call calc_enhance_4()
1406 #elif (ENHMOD==5)
1407  call calc_enhance_5()
1408 #else
1409  errormsg = ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
1410  call error(errormsg)
1411 #endif
1412 
1413 ! ------ Ice-free, relaxed bedrock
1414 
1415 #elif (ANF_DAT==2)
1416 
1417 call topography2(dxi, deta)
1418 
1419 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1420 
1421 call boundary(time_init, dtime, dxi, deta, &
1422  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1423 
1424 as_perp_apl = 0.0_dp
1425 
1426 smb_corr = 0.0_dp
1427 
1428 q_bm = 0.0_dp
1429 q_tld = 0.0_dp
1430 q_b_tot = 0.0_dp
1431 
1432 p_b_w = 0.0_dp
1433 h_w = 0.0_dp
1434 
1435 call init_temp_water_age_2()
1436 
1437 #if (ENHMOD==1)
1438  call calc_enhance_1()
1439 #elif (ENHMOD==2)
1440  call calc_enhance_2()
1441 #elif (ENHMOD==3)
1442  call calc_enhance_3(time_init)
1443 #elif (ENHMOD==4)
1444  call calc_enhance_4()
1445 #elif (ENHMOD==5)
1446  call calc_enhance_5()
1447 #else
1448  errormsg = ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
1449  call error(errormsg)
1450 #endif
1451 
1452 ! ------ Read initial state from output of previous simulation
1453 
1454 #elif (ANF_DAT==3)
1455 
1456 call topography3(dxi, deta, z_sl, anfdatname)
1457 
1458 call boundary(time_init, dtime, dxi, deta, &
1459  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1460 
1461 where ((maske==0_i1b).or.(maske==3_i1b))
1462  ! grounded or floating ice
1464 elsewhere ! maske==1_i1b or 2_i1b, ice-free land or sea
1465  as_perp_apl = 0.0_dp
1466 end where
1467 
1468 smb_corr = 0.0_dp
1469 
1470 q_b_tot = q_bm + q_tld
1471 
1472 #endif
1473 
1474 !-------- Inner-point flag --------
1475 
1476 flag_inner_point = .true.
1477 
1478 flag_inner_point(0,:) = .false.
1479 flag_inner_point(jmax,:) = .false.
1480 
1481 flag_inner_point(:,0) = .false.
1482 flag_inner_point(:,imax) = .false.
1483 
1484 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
1485 
1486 #if (GRID==0 || GRID==1)
1487 
1488 do ir=-imax, imax
1489 do jr=-jmax, jmax
1490  dist_dxdy(jr,ir) = sqrt( (real(ir,dp)*dxi)**2 + (real(jr,dp)*deta)**2 )
1491  ! distortion due to stereographic projection not accounted for
1492 end do
1493 end do
1494 
1495 #endif
1496 
1497 !-------- Initial velocities --------
1498 
1499 call calc_temp_melt()
1500 call flag_update_gf_gl_cf()
1501 
1502 #if (DYNAMICS==1 || DYNAMICS==2)
1503 
1504 call calc_vxy_b_sia(time, z_sl)
1505 call calc_vxy_sia(dzeta_c, dzeta_t)
1506 
1507 #if (MARGIN==3 || DYNAMICS==2)
1508 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
1509 #endif
1510 
1511 call calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
1512 
1513 #if (MARGIN==3)
1514 call calc_vz_floating(z_sl, dxi, deta, dzeta_c)
1515 #endif
1516 
1517 #elif (DYNAMICS==0)
1518 
1519 call calc_vxy_static()
1520 call calc_vz_static()
1521 
1522 #else
1523 
1524 errormsg = ' >>> sico_init: DYNAMICS must be either 0, 1 or 2!'
1525 call error(errormsg)
1526 
1527 #endif
1528 
1529 call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
1530 
1531 !-------- Initial enthalpies --------
1532 
1533 #if (CALCMOD==0 || CALCMOD==-1)
1534 
1535 do i=0, imax
1536 do j=0, jmax
1537 
1538  do kc=0, kcmax
1539  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1540  end do
1541 
1542  do kt=0, ktmax
1543  enth_t(kt,j,i) = enth_c(0,j,i)
1544  end do
1545 
1546 end do
1547 end do
1548 
1549 #elif (CALCMOD==1)
1550 
1551 do i=0, imax
1552 do j=0, jmax
1553 
1554  do kc=0, kcmax
1555  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1556  end do
1557 
1558  if ( (maske(j,i) == 0_i1b).and.(n_cts(j,i) == 1_i1b) ) then
1559  do kt=0, ktmax
1560  enth_t(kt,j,i) = enth_fct_temp_omega(temp_t_m(kt,j,i), omega_t(kt,j,i))
1561  end do
1562  else
1563  do kt=0, ktmax
1564  enth_t(kt,j,i) = enth_c(0,j,i)
1565  end do
1566  end if
1567 
1568 end do
1569 end do
1570 
1571 #elif (CALCMOD==2 || CALCMOD==3)
1572 
1573 do i=0, imax
1574 do j=0, jmax
1575 
1576  do kc=0, kcmax
1577  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), omega_c(kc,j,i))
1578  end do
1579 
1580  do kt=0, ktmax
1581  enth_t(kt,j,i) = enth_c(0,j,i)
1582  end do
1583 
1584 end do
1585 end do
1586 
1587 #else
1588 
1589 errormsg = ' >>> sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!'
1590 call error(errormsg)
1591 
1592 #endif
1593 
1594 !-------- Initialize time-series files --------
1595 
1596 ! ------ Time-series file for the ice sheet on the whole
1597 
1598 filename_with_path = trim(outpath)//'/'//trim(runname)//'.ser'
1599 
1600 open(12, iostat=ios, file=trim(filename_with_path), status='new')
1601 
1602 if (ios /= 0) then
1603  errormsg = ' >>> sico_init: Error when opening the ser file!'
1604  call error(errormsg)
1605 end if
1606 
1607 write(12,1102)
1608 write(12,1103)
1609 
1610  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1611  ' V(m^3) V_g(m^3) V_f(m^3)', &
1612  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1613  ' V_sle(m) V_t(m^3)', &
1614  ' A_t(m^2)',/, &
1615  ' H_max(m) H_t_max(m)', &
1616  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
1617  1103 format('----------------------------------------------------', &
1618  '---------------------------------------')
1619 
1620 ! ------ Time-series file for deep boreholes
1621 
1622 n_core = 0 ! No boreholes defined
1623 
1624 filename_with_path = trim(outpath)//'/'//trim(runname)//'.core'
1625 
1626 open(14, iostat=ios, file=trim(filename_with_path), status='new')
1627 
1628 if (ios /= 0) then
1629  errormsg = ' >>> sico_init: Error when opening the core file!'
1630  call error(errormsg)
1631 end if
1632 
1633 write(14,'(1x,a)') '---------------------'
1634 write(14,'(1x,a)') 'No boreholes defined.'
1635 write(14,'(1x,a)') '---------------------'
1636 
1637 !-------- Output of the initial state --------
1638 
1639 #if (defined(OUTPUT_INIT))
1640 
1641 #if (OUTPUT_INIT==0)
1642  flag_init_output = .false.
1643 #elif (OUTPUT_INIT==1)
1644  flag_init_output = .true.
1645 #else
1646  errormsg = ' >>> sico_init: OUTPUT_INIT must be either 0 or 1!'
1647  call error(errormsg)
1648 #endif
1649 
1650 #else
1651  flag_init_output = .true. ! default for undefined parameter OUTPUT_INIT
1652 #endif
1653 
1654 #if (OUTPUT==1)
1655 
1656 #if (ERGDAT==0)
1657  flag_3d_output = .false.
1658 #elif (ERGDAT==1)
1659  flag_3d_output = .true.
1660 #else
1661  errormsg = ' >>> sico_init: ERGDAT must be either 0 or 1!'
1662  call error(errormsg)
1663 #endif
1664 
1665  if (flag_init_output) &
1666  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1667  flag_3d_output, ndat2d, ndat3d)
1668 
1669 #elif (OUTPUT==2)
1670 
1671 if (time_output(1) <= time_init+eps) then
1672 
1673 #if (ERGDAT==0)
1674  flag_3d_output = .false.
1675 #elif (ERGDAT==1)
1676  flag_3d_output = .true.
1677 #else
1678  errormsg = ' >>> sico_init: ERGDAT must be either 0 or 1!'
1679  call error(errormsg)
1680 #endif
1681 
1682  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1683  flag_3d_output, ndat2d, ndat3d)
1684 
1685 end if
1686 
1687 #elif (OUTPUT==3)
1688 
1689  flag_3d_output = .false.
1690 
1691  if (flag_init_output) &
1692  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1693  flag_3d_output, ndat2d, ndat3d)
1694 
1695 if (time_output(1) <= time_init+eps) then
1696 
1697  flag_3d_output = .true.
1698 
1699  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1700  flag_3d_output, ndat2d, ndat3d)
1701 
1702 end if
1703 
1704 #else
1705  errormsg = ' >>> sico_init: OUTPUT must be either 1, 2 or 3!'
1706  call error(errormsg)
1707 #endif
1708 
1709 if (flag_init_output) then
1710  call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1711  call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1712 end if
1713 
1714 end subroutine sico_init
1715 
1716 !-------------------------------------------------------------------------------
1717 !> Definition of the initial surface and bedrock topography
1718 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1719 !! For ice-free initial topography with relaxed lithosphere
1720 !! (not defined for this domain, thus routine stops execution of SICOPOLIS).
1721 !<------------------------------------------------------------------------------
1722 subroutine topography1(dxi, deta)
1723 
1724 #if (GRID==0 || GRID==1)
1726 #endif
1727 
1728  use metric_m
1729  use topograd_m
1730 
1731 implicit none
1732 
1733 real(dp), intent(out) :: dxi, deta
1734 
1735 integer(i4b) :: i, j, m, n
1736 integer(i4b) :: ios
1737 real(dp) :: xi0, eta0
1738 character(len=256) :: filename_with_path
1739 character :: ch_dummy
1740 
1741 !-------- Set topography --------
1742 
1743 #if (GRID==0 || GRID==1)
1744 
1745 dxi = dx *1000.0_dp ! km -> m
1746 deta = dx *1000.0_dp ! km -> m
1747 
1748 xi0 = x0 *1000.0_dp ! km -> m
1749 eta0 = y0 *1000.0_dp ! km -> m
1750 
1751 do i=0, imax
1752  xi(i) = xi0 + real(i,dp)*dxi
1753 end do
1754 
1755 do j=0, jmax
1756  eta(j) = eta0 + real(j,dp)*deta
1757 end do
1758 
1759 #elif (GRID==2)
1760 
1761 errormsg = ' >>> topography1: GRID==2 not allowed for this application!'
1762 call error(errormsg)
1763 
1764 #endif
1765 
1766 #if (defined(ZS_PRESENT_FILE))
1767 
1768 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1769  trim(zs_present_file)
1770 
1771 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1772 
1773 if (ios /= 0) then
1774  errormsg = ' >>> topography1: Error when opening the zs file!'
1775  call error(errormsg)
1776 end if
1777 
1778 m = len_trim(filename_with_path)
1779 
1780 if (filename_with_path(m-2:m) /= 'txt') then
1781  do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1782 end if
1783 
1784 do j=jmax, 0, -1
1785  read(21, fmt=*) (zs(j,i), i=0,imax)
1786 end do
1787 
1788 close(21, status='keep')
1789 
1790 #else
1791 
1792 errormsg = ' >>> topography1: ZS_PRESENT_FILE must be defined!'
1793 call error(errormsg)
1794 
1795 #endif
1796 
1797 #if (defined(ZL_PRESENT_FILE))
1798 
1799 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1800  trim(zl_present_file)
1801 
1802 open(22, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1803 
1804 if (ios /= 0) then
1805  errormsg = ' >>> topography1: Error when opening the zl file!'
1806  call error(errormsg)
1807 end if
1808 
1809 m = len_trim(filename_with_path)
1810 
1811 if (filename_with_path(m-2:m) /= 'txt') then
1812  do n=1, 6; read(22, fmt='(a)') ch_dummy; end do
1813 end if
1814 
1815 do j=jmax, 0, -1
1816  read(22, fmt=*) (zl(j,i), i=0,imax)
1817 end do
1818 
1819 close(22, status='keep')
1820 
1821 #else
1822 
1823 errormsg = ' >>> topography1: ZL_PRESENT_FILE must be defined!'
1824 call error(errormsg)
1825 
1826 #endif
1827 
1828 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1829  trim(zl0_file)
1830 
1831 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1832 
1833 if (ios /= 0) then
1834  errormsg = ' >>> topography1: Error when opening the zl0 file!'
1835  call error(errormsg)
1836 end if
1837 
1838 m = len_trim(filename_with_path)
1839 
1840 if (filename_with_path(m-2:m) /= 'txt') then
1841  do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
1842 end if
1843 
1844 do j=jmax, 0, -1
1845  read(23, fmt=*) (zl0(j,i), i=0,imax)
1846 end do
1847 
1848 close(23, status='keep')
1849 
1850 #if (defined(ZB_PRESENT_FILE))
1851 
1852 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1853  trim(zb_present_file)
1854 
1855 open(25, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1856 
1857 if (ios /= 0) then
1858  errormsg = ' >>> topography1: Error when opening the zb file!'
1859  call error(errormsg)
1860 end if
1861 
1862 m = len_trim(filename_with_path)
1863 
1864 if (filename_with_path(m-2:m) /= 'txt') then
1865  do n=1, 6; read(25, fmt='(a)') ch_dummy; end do
1866 end if
1867 
1868 do j=jmax, 0, -1
1869  read(25, fmt=*) (zb(j,i), i=0,imax)
1870 end do
1871 
1872 close(25, status='keep')
1873 
1874 #else
1875 
1876 write(6, fmt='(a)') ' >>> topography1: ZB_PRESENT_FILE not defined,'
1877 write(6, fmt='(a)') ' thus zb = zl assumed.'
1878 
1879 zb = zl
1880 
1881 #endif
1882 
1883 do i=1, imax-1
1884 do j=1, jmax-1
1885  if ((zs(j,i)-zb(j,i)) > eps_dp) then
1886  maske(j,i) = 0_i1b
1887  else
1888  zs(j,i) = zb(j,i)
1889  maske(j,i) = 1_i1b
1890  end if
1891 end do
1892 end do
1893 
1894 maske(0,:) = 1_i1b
1895 maske(jmax,:) = 1_i1b
1896 
1897 maske(:,0) = 1_i1b
1898 maske(:,imax) = 1_i1b
1899 
1900 !-------- Further stuff --------
1901 
1902 do i=0, imax
1903 do j=0, jmax
1904 
1905  zm(j,i) = zb(j,i)
1906  n_cts(j,i) = -1_i1b
1907  kc_cts(j,i) = 0
1908 
1909  h_c(j,i) = zs(j,i)-zm(j,i)
1910  h_t(j,i) = 0.0_dp
1911 
1912  dzs_dtau(j,i) = 0.0_dp
1913  dzm_dtau(j,i) = 0.0_dp
1914  dzb_dtau(j,i) = 0.0_dp
1915  dzl_dtau(j,i) = 0.0_dp
1916  dh_c_dtau(j,i) = 0.0_dp
1917  dh_t_dtau(j,i) = 0.0_dp
1918 
1919 end do
1920 end do
1921 
1922 maske_old = maske
1923 
1924 !-------- Geographic coordinates, metric tensor,
1925 ! gradients of the topography --------
1926 
1927 do i=0, imax
1928 do j=0, jmax
1929 
1930 #if (GRID==0 || GRID==1) /* Stereographic projection */
1931 
1932  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1933  lambda0, phi0, lambda(j,i), phi(j,i))
1934 
1935 #elif (GRID==2) /* Geographic coordinates */
1936 
1937  lambda(j,i) = xi(i)
1938  phi(j,i) = eta(j)
1939 
1940 #endif
1941 
1942 end do
1943 end do
1944 
1945 call metric()
1946 
1947 #if (TOPOGRAD==0)
1948 call topograd_1(dxi, deta, 1)
1949 #elif (TOPOGRAD==1)
1950 call topograd_2(dxi, deta, 1)
1951 #endif
1952 
1953 !-------- Corresponding area of grid points --------
1954 
1955 do i=0, imax
1956 do j=0, jmax
1957  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1958 end do
1959 end do
1960 
1961 end subroutine topography1
1962 
1963 !-------------------------------------------------------------------------------
1964 !> Definition of the initial surface and bedrock topography
1965 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1966 !! For ice-free initial topography with relaxed lithosphere.
1967 !<------------------------------------------------------------------------------
1968 subroutine topography2(dxi, deta)
1969 
1970 #if (GRID==0 || GRID==1)
1972 #endif
1973 
1974  use metric_m
1975  use topograd_m
1976 
1977 implicit none
1978 
1979 real(dp), intent(out) :: dxi, deta
1980 
1981 integer(i4b) :: i, j, m, n
1982 integer(i4b) :: ios
1983 real(dp) :: xi0, eta0
1984 character(len=256) :: filename_with_path
1985 character :: ch_dummy
1986 
1987 !-------- Set topography --------
1988 
1989 #if (GRID==0 || GRID==1)
1990 
1991 dxi = dx *1000.0_dp ! km -> m
1992 deta = dx *1000.0_dp ! km -> m
1993 
1994 xi0 = x0 *1000.0_dp ! km -> m
1995 eta0 = y0 *1000.0_dp ! km -> m
1996 
1997 do i=0, imax
1998  xi(i) = xi0 + real(i,dp)*dxi
1999 end do
2000 
2001 do j=0, jmax
2002  eta(j) = eta0 + real(j,dp)*deta
2003 end do
2004 
2005 #elif (GRID==2)
2006 
2007 errormsg = ' >>> topography2: GRID==2 not allowed for this application!'
2008 call error(errormsg)
2009 
2010 #endif
2011 
2012 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2013  trim(zl0_file)
2014 
2015 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
2016 
2017 if (ios /= 0) then
2018  errormsg = ' >>> topography2: Error when opening the zl0 file!'
2019  call error(errormsg)
2020 end if
2021 
2022 m = len_trim(filename_with_path)
2023 
2024 if (filename_with_path(m-2:m) /= 'txt') then
2025  do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
2026 end if
2027 
2028 do j=jmax, 0, -1
2029  read(23, fmt=*) (zl0(j,i), i=0,imax)
2030 end do
2031 
2032 close(23, status='keep')
2033 
2034 maske = 1_i1b
2035 
2036 !-------- Further stuff --------
2037 
2038 do i=0, imax
2039 do j=0, jmax
2040 
2041  zs(j,i) = zl0(j,i)
2042  zb(j,i) = zl0(j,i)
2043  zl(j,i) = zl0(j,i)
2044 
2045  zm(j,i) = zb(j,i)
2046  n_cts(j,i) = -1_i1b
2047  kc_cts(j,i) = 0
2048 
2049  h_c(j,i) = 0.0_dp
2050  h_t(j,i) = 0.0_dp
2051 
2052  dzs_dtau(j,i) = 0.0_dp
2053  dzm_dtau(j,i) = 0.0_dp
2054  dzb_dtau(j,i) = 0.0_dp
2055  dzl_dtau(j,i) = 0.0_dp
2056  dh_c_dtau(j,i) = 0.0_dp
2057  dh_t_dtau(j,i) = 0.0_dp
2058 
2059 end do
2060 end do
2061 
2062 maske_old = maske
2063 
2064 !-------- Geographic coordinates, metric tensor,
2065 ! gradients of the topography --------
2066 
2067 do i=0, imax
2068 do j=0, jmax
2069 
2070 #if (GRID==0 || GRID==1) /* Stereographic projection */
2071 
2072  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
2073  lambda0, phi0, lambda(j,i), phi(j,i))
2074 
2075 #elif (GRID==2) /* Geographic coordinates */
2076 
2077  lambda(j,i) = xi(i)
2078  phi(j,i) = eta(j)
2079 
2080 #endif
2081 
2082 end do
2083 end do
2084 
2085 call metric()
2086 
2087 #if (TOPOGRAD==0)
2088 call topograd_1(dxi, deta, 1)
2089 #elif (TOPOGRAD==1)
2090 call topograd_2(dxi, deta, 1)
2091 #endif
2092 
2093 !-------- Corresponding area of grid points --------
2094 
2095 do i=0, imax
2096 do j=0, jmax
2097  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
2098 end do
2099 end do
2100 
2101 end subroutine topography2
2102 
2103 !-------------------------------------------------------------------------------
2104 !> Definition of the initial surface and bedrock topography
2105 !! (including gradients) and of the horizontal grid spacings dxi, deta.
2106 !! For initial topography from previous simulation.
2107 !<------------------------------------------------------------------------------
2108 subroutine topography3(dxi, deta, z_sl, anfdatname)
2109 
2110  use read_m, only : read_erg_nc
2111 
2112 #if (GRID==0 || GRID==1)
2114 #endif
2115 
2116  use metric_m
2117  use topograd_m
2118 
2119 implicit none
2120 
2121 character(len=100), intent(in) :: anfdatname
2122 
2123 real(dp), intent(out) :: dxi, deta, z_sl
2124 
2125 integer(i4b) :: i, j, m, n
2126 integer(i4b) :: ios
2127 character(len=256) :: filename_with_path
2128 character :: ch_dummy
2129 
2130 !-------- Read data from time-slice file of previous simulation --------
2131 
2132 call read_erg_nc(z_sl, anfdatname)
2133 
2134 !-------- Set topography of the relaxed bedrock --------
2135 
2136 #if (GRID==0 || GRID==1)
2137 
2138 dxi = dx *1000.0_dp ! km -> m
2139 deta = dx *1000.0_dp ! km -> m
2140 
2141 #elif (GRID==2)
2142 
2143 errormsg = ' >>> topography3: GRID==2 not allowed for this application!'
2144 call error(errormsg)
2145 
2146 #endif
2147 
2148 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
2149  trim(zl0_file)
2150 
2151 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
2152 
2153 if (ios /= 0) then
2154  errormsg = ' >>> topography3: Error when opening the zl0 file!'
2155  call error(errormsg)
2156 end if
2157 
2158 m = len_trim(filename_with_path)
2159 
2160 if (filename_with_path(m-2:m) /= 'txt') then
2161  do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
2162 end if
2163 
2164 do j=jmax, 0, -1
2165  read(23, fmt=*) (zl0(j,i), i=0,imax)
2166 end do
2167 
2168 close(23, status='keep')
2169 
2170 !-------- Geographic coordinates, metric tensor,
2171 ! gradients of the topography --------
2172 
2173 do i=0, imax
2174 do j=0, jmax
2175 
2176 #if (GRID==0 || GRID==1) /* Stereographic projection */
2177 
2178  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
2179  lambda0, phi0, lambda(j,i), phi(j,i))
2180 
2181 #elif (GRID==2) /* Geographic coordinates */
2182 
2183  lambda(j,i) = xi(i)
2184  phi(j,i) = eta(j)
2185 
2186 #endif
2187 
2188 end do
2189 end do
2190 
2191 call metric()
2192 
2193 #if (TOPOGRAD==0)
2194 call topograd_1(dxi, deta, 1)
2195 #elif (TOPOGRAD==1)
2196 call topograd_2(dxi, deta, 1)
2197 #endif
2198 
2199 !-------- Corresponding area of grid points --------
2200 
2201 do i=0, imax
2202 do j=0, jmax
2203  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
2204 end do
2205 end do
2206 
2207 end subroutine topography3
2208 
2209 !-------------------------------------------------------------------------------
2210 
2211 end module sico_init_m
2212 !
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...
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_erg_nc(z_sl, filename)
Reading of time-slice files in native binary or NetCDF format.
Definition: read_m.F90:55
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...
subroutine sico_init(delta_ts, glac_index, mean_accum, dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, time, time_init, time_end, time_output, dxi, deta, dzeta_c, dzeta_t, dzeta_r, z_sl, dzsl_dtau, z_mar, ii, jj, nn, ndat2d, ndat3d, n_output, runname)
Main routine of sico_init_m: Initialisations for SICOPOLIS.
Definition: sico_init_m.F90:60
real(dp) rho_w
RHO_W: Density of pure water.
subroutine, public calc_enhance_2()
Computation of the flow enhancement factor. Case ENHMOD==2: two different values depending on age for...
real(dp) a
A: Semi-major axis of the planet.
subroutine, public calc_vz_floating(z_sl, dxi, deta, dzeta_c)
Computation of the vertical velocity vz for floating ice.
Definition: calc_vz_m.F90:259
subroutine error(error_message)
Main routine of error_m: Writing of error messages and stopping execution.
Definition: error_m.F90:47
subroutine, public output1(runname, time, delta_ts, glac_index, z_sl, flag_3d_output, ndat2d, ndat3d)
Writing of time-slice files in native binary or NetCDF format.
Definition: output_m.F90:60
integer(i4b) specmap_time_max
specmap_time_max: Maximum time of the data values for the sea level
Material properties of ice: Rate factor, heat conductivity, specific heat (heat capacity), creep function, viscosity.
real(dp), parameter eps
eps: Small number
subroutine, public read_phys_para()
Reading of physical parameters.
Definition: read_m.F90:903
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)
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) 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
real(dp), dimension(0:jmax, 0:imax) smb_corr_in
smb_corr_in(j,i): Prescribed SMB correction
Definition: sico_vars_m.F90:43
real(dp) target_topo_tau
target_topo_tau: Relaxation time for target-topography adjustment
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, 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:5214
integer(i4b) grip_time_max
grip_time_max: Maximum time of the data values for the surface temperature anomaly ...
real(dp), dimension(0:jmax, 0:imax) q_tld
Q_tld(j,i): Water drainage rate from the temperate layer.
subroutine, public calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy and the flux across ...
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars_m.F90:35
integer(i4b), parameter rcl1
rcl1: Maximum length of record for input files (with factor 3 safety margin)
Initialisations for SICOPOLIS.
Definition: sico_init_m.F90:35
Computation of the flow enhancement factor.
integer(i4b) ndata_grip
ndata_grip: Number of data values for the surface temperature anomaly
Computation of the horizontal velocity vx, vy.
Definition: calc_vxy_m.F90:36
integer(i4b), dimension(0:jmax, 0:imax) q_weert
q_weert(j,i): Weertman exponent for the basal pressure
subroutine, public calc_temp_melt()
Computation of the melting temperatures.
real(dp) ea
ea: Abbreviation for exp(aa)
real(dp), 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), 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), dimension(-190:10) kappa
KAPPA(n): Tabulated values for the heat conductivity of ice.
Computation of the vertical velocity vz.
Definition: calc_vz_m.F90:35
Writing of error messages and stopping execution.
Definition: error_m.F90:35
real(dp), dimension(0:jmax, 0:imax) dzm_dtau
dzm_dtau(j,i): Derivative of zm by tau (time)
real(dp), dimension(0:krmax) zeta_r
zeta_r(kr): Sigma coordinate zeta_r of grid point kr
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) omega_c
omega_c(kc,j,i): Water content in the upper (kc) ice domain
character(len=64) ch_domain_long
ch_domain_long: Long name of the computational domain
real(dp), parameter eps_dp
eps_dp: Small number to double-precision accuracy
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.
integer(i1b), dimension(0:jmax, 0:imax) maske_maxextent
maske_maxextent(j,i): Maximum ice extent mask. 0: not allowed to glaciate, 1: allowed to glaciate...
NetCDF error capturing.
Definition: nc_check_m.F90:35
real(dp) time_target_topo_final
time_target_topo_final: Final time for target-topography adjustment
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:561
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, 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) time_target_topo_init
time_target_topo_init: Initial time for target-topography adjustment
integer(i4b) specmap_time_stp
specmap_time_stp: Time step of the data values for the sea level
integer(i4b), parameter rcl2
rcl2: Maximum length of record for input mask files (with factor 3 safety margin) ...
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:3869
integer(i1b) forcing_flag
forcing_flag: Flag for the forcing type. 1: forcing by a spatially constant surface temperature anoma...
subroutine check(status, ch_calling_routine)
NetCDF error capturing.
Definition: nc_check_m.F90:48
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
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(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
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
Writing of output data on files.
Definition: output_m.F90:36
real(dp), dimension(0:jmax, 0:imax) dzb_dtau
dzb_dtau(j,i): Derivative of zb by tau (time)
real(dp), dimension(0:kcmax) zeta_c
zeta_c(kc): Sigma coordinate zeta_c of grid point kc
subroutine, public read_target_topo_nc(target_topo_dat_name)
Reading of the target-topography file (in NetCDF format).
Definition: read_m.F90:734
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...