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