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