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