37 dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
38 time, time_init, time_end, time_output, &
39 dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
40 z_sl, dzsl_dtau, z_mar, &
42 ndat2d, ndat3d, n_output, &
54 integer(i4b),
intent(out) :: ii((imax+1)*(jmax+1)), &
55 jj((imax+1)*(jmax+1)), &
57 integer(i4b),
intent(out) :: ndat2d, ndat3d
58 integer(i4b),
intent(out) :: n_output
59 real(dp),
intent(out) :: delta_ts, glac_index
60 real(dp),
intent(out) :: mean_accum
61 real(dp),
intent(out) :: dtime, dtime_temp, dtime_wss, &
63 real(dp),
intent(out) :: time, time_init, time_end, time_output(100)
64 real(dp),
intent(out) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
65 real(dp),
intent(out) :: z_sl, dzsl_dtau, z_mar
66 character(len=100),
intent(out) :: runname
68 integer(i4b) :: i, j, kc, kt, kr, m, n, ir, jr
69 integer(i4b) :: ios, ios1, ios2, ios3, ios4
71 real(dp) :: dtime0, dtime_temp0, dtime_wss0, dtime_out0, dtime_ser0
72 real(dp) :: time_init0, time_end0, time_output0(100)
74 character (len=100) :: anfdatname
75 character (len=256) :: shell_command
77 logical :: flag_3d_output
79 character(len=64),
parameter :: fmt1 =
'(a)', &
83 character(len= 8) :: ch_imax
84 character(len=128) :: fmt4
86 write(ch_imax, fmt=
'(i8)') imax
87 write(fmt4, fmt=
'(a)')
'('//trim(adjustl(ch_imax))//
'(i1),i1)'
92 ch_domain_long =
'Antarctica'
93 ch_domain_short =
'ant'
96 ch_domain_long =
'Austfonna'
97 ch_domain_short =
'asf'
99 #elif defined(EMTP2SGE)
100 ch_domain_long =
'EISMINT Phase 2 Simplified Geometry Experiment'
101 ch_domain_short =
'emtp2sge'
104 ch_domain_long =
'Greenland'
105 ch_domain_short =
'grl'
108 ch_domain_long =
'Northern hemisphere'
109 ch_domain_short =
'nhem'
112 ch_domain_long =
'Scandinavia and Eurasia'
113 ch_domain_short =
'scand'
116 ch_domain_long =
'Tibet'
117 ch_domain_short =
'tibet'
120 ch_domain_long =
'North polar cap of Mars'
121 ch_domain_short =
'nmars'
124 ch_domain_long =
'South polar cap of Mars'
125 ch_domain_short =
'smars'
128 ch_domain_long =
'XYZ'
129 ch_domain_short =
'xyz'
131 ch_domain_long = trim(ch_domain_long)//
'/ISMIP HEINO'
135 stop
' sico_init: No valid domain specified!'
156 #if (CALCZS==4 || MARGIN==3 || DYNAMICS==2)
159 call lis_initialize(ierr)
170 s_t = s_t *1.0e-09_dp
171 x_hat = x_hat *1.0e+03_dp
172 y_hat = y_hat *1.0e+03_dp
173 rad = rad *1.0e+03_dp
174 b_min = b_min /year_sec
175 b_max = b_max /year_sec
184 if ( trim(version) /= trim(sico_version) ) &
185 stop
' sico_init: SICOPOLIS version not compatible with header file!'
189 #if (!defined(DYNAMICS))
190 stop
' sico_init: DYNAMICS not defined in the header file!'
193 #if (!defined(CALCMOD))
194 stop
' sico_init: CALCMOD not defined in the header file!'
197 #if (defined(ENTHMOD))
198 write(6, fmt=
'(a)')
' sico_init: ENTHMOD must not be defined any more.'
199 write(6, fmt=
'(a)')
' Please update your header file!'
208 if (abs(dx-50.0_dp) < eps)
then
210 if ((imax /= 80).or.(jmax /= 80))
then
211 stop
' sico_init: IMAX and/or JMAX wrong!'
214 else if (abs(dx-25.0_dp) < eps)
then
216 if ((imax /= 160).or.(jmax /= 160))
then
217 stop
' sico_init: IMAX and/or JMAX wrong!'
221 stop
' sico_init: DX wrong!'
226 stop
' sico_init: GRID==1 not allowed for this application!'
230 stop
' sico_init: GRID==2 not allowed for this application!'
238 #if (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1)
241 write(6, fmt=
'(a)')
' sico_init: For options CALCMOD==0, 2, 3 or -1,'
242 write(6, fmt=
'(a)')
' the separate kt domain is redundant.'
243 write(6, fmt=
'(a)')
' Therefore, consider setting KTMAX to 1.'
244 write(6, fmt=
'(a)')
' '
253 stop
' sico_init: Option ADV_HOR==1 (central differences) not defined!'
261 #if ((MARGIN==3 && DYNAMICS==1) || DYNAMICS==2) /* requires SSA or SStA */
263 write(6, fmt=
'(a)')
' sico_init: Distortion correction not implemented'
264 write(6, fmt=
'(a)')
' for the shallow shelf approximation (SSA)'
265 write(6, fmt=
'(a)')
' or the shelfy stream approximation (SStA)'
266 write(6, fmt=
'(a)')
' -> GRID==0 required!'
278 dtime_temp0 = dtime_temp0
282 dzeta_c = 1.0_dp/
real(kcmax,dp)
283 dzeta_t = 1.0_dp/
real(ktmax,dp)
284 dzeta_r = 1.0_dp/
real(krmax,dp)
293 if (deform >= eps)
then
295 flag_aa_nonzero = .true.
303 eaz_c_quotient(kc) = 0.0_dp
306 zeta_c(kc) = kc*dzeta_c
307 eaz_c(kc) = exp(aa*zeta_c(kc))
308 eaz_c_quotient(kc) = (eaz_c(kc)-1.0_dp)/(ea-1.0_dp)
314 eaz_c_quotient(kc) = 1.0_dp
318 flag_aa_nonzero = .false.
326 eaz_c_quotient(kc) = 0.0_dp
329 zeta_c(kc) = kc*dzeta_c
331 eaz_c_quotient(kc) = zeta_c(kc)
337 eaz_c_quotient(kc) = 1.0_dp
347 zeta_t(kt) = kt*dzeta_t
359 zeta_r(kr) = kr*dzeta_r
373 if ((i <= imax).and.(j <= jmax))
then
385 anfdatname = anfdatname
387 #if (defined(YEAR_ZERO))
388 year_zero = year_zero
390 year_zero = 2000.0_dp
393 time_init0 = time_init0
394 time_end0 = time_end0
395 dtime_ser0 = dtime_ser0
397 #if (OUTPUT==1 || OUTPUT==3)
398 dtime_out0 = dtime_out0
401 #if (OUTPUT==2 || OUTPUT==3)
403 time_output0( 1) = time_out0_01
404 time_output0( 2) = time_out0_02
405 time_output0( 3) = time_out0_03
406 time_output0( 4) = time_out0_04
407 time_output0( 5) = time_out0_05
408 time_output0( 6) = time_out0_06
409 time_output0( 7) = time_out0_07
410 time_output0( 8) = time_out0_08
411 time_output0( 9) = time_out0_09
412 time_output0(10) = time_out0_10
413 time_output0(11) = time_out0_11
414 time_output0(12) = time_out0_12
415 time_output0(13) = time_out0_13
416 time_output0(14) = time_out0_14
417 time_output0(15) = time_out0_15
418 time_output0(16) = time_out0_16
419 time_output0(17) = time_out0_17
420 time_output0(18) = time_out0_18
421 time_output0(19) = time_out0_19
422 time_output0(20) = time_out0_20
427 shell_command =
'if [ ! -d'
428 shell_command = trim(shell_command)//
' '//outpath
429 shell_command = trim(shell_command)//
' '//
'] ; then mkdir'
430 shell_command = trim(shell_command)//
' '//outpath
431 shell_command = trim(shell_command)//
' '//
'; fi'
432 call system(trim(shell_command))
435 open(10, iostat=ios, &
436 file=outpath//
'/'//trim(runname)//
'.log', &
438 if (ios /= 0) stop
' sico_init: Error when opening the log file!'
440 write(10, fmt=trim(fmt1))
'Computational domain:'
441 write(10, fmt=trim(fmt1)) trim(ch_domain_long)
442 write(10, fmt=trim(fmt1))
' '
444 write(10, fmt=trim(fmt2))
'imax =', imax
445 write(10, fmt=trim(fmt2))
'jmax =', jmax
446 write(10, fmt=trim(fmt2))
'kcmax =', kcmax
447 write(10, fmt=trim(fmt2))
'ktmax =', ktmax
448 write(10, fmt=trim(fmt2))
'krmax =', krmax
449 write(10, fmt=trim(fmt1))
' '
451 write(10, fmt=trim(fmt3))
'a =', aa
452 write(10, fmt=trim(fmt1))
' '
454 #if (GRID==0 || GRID==1)
455 write(10, fmt=trim(fmt3))
'x0 =', x0
456 write(10, fmt=trim(fmt3))
'y0 =', y0
457 write(10, fmt=trim(fmt3))
'dx =', dx
459 write(10, fmt=trim(fmt3))
'dlambda =', dlambda
460 write(10, fmt=trim(fmt3))
'dphi =', dphi
462 write(10, fmt=trim(fmt1))
' '
464 write(10, fmt=trim(fmt3))
'year_zero =', year_zero
465 write(10, fmt=trim(fmt3))
'time_init =', time_init0
466 write(10, fmt=trim(fmt3))
'time_end =', time_end0
467 write(10, fmt=trim(fmt3))
'dtime =', dtime0
468 write(10, fmt=trim(fmt3))
'dtime_temp =', dtime_temp0
469 #if (OUTPUT==1 || OUTPUT==3)
470 write(10, fmt=trim(fmt3))
'dtime_out =', dtime_out0
472 write(10, fmt=trim(fmt3))
'dtime_ser =', dtime_ser0
473 write(10, fmt=trim(fmt1))
' '
475 write(10, fmt=trim(fmt2))
'ANF_DAT =', anf_dat
476 write(10, fmt=trim(fmt1))
'mask_present file = '//mask_present_file
477 #if (ANF_DAT==1 && defined(TEMP_INIT))
478 write(10, fmt=trim(fmt2))
'TEMP_INIT =', temp_init
481 write(10, fmt=trim(fmt1))
'Initial-value file = '//anfdatname
482 write(10, fmt=trim(fmt1))
'Path to initial-value file = '//anfdatpath
484 write(10, fmt=trim(fmt1))
' '
486 write(10, fmt=trim(fmt1))
'Physical-parameter file = '//phys_para_file
487 write(10, fmt=trim(fmt1))
' '
489 #if (CALCZS==3 || CALCZS==4)
490 write(10, fmt=trim(fmt3))
'ovi_weight =', ovi_weight
492 write(10, fmt=trim(fmt3))
'omega_sor =', omega_sor
494 write(10, fmt=trim(fmt1))
' '
497 write(10, fmt=trim(fmt3))
'temp_min =', temp_min
498 write(10, fmt=trim(fmt3))
's_t =', s_t
499 write(10, fmt=trim(fmt3))
'x_hat =', x_hat
500 write(10, fmt=trim(fmt3))
'y_hat =', y_hat
501 write(10, fmt=trim(fmt3))
'rad =', rad
502 write(10, fmt=trim(fmt3))
'b_min =', b_min
503 write(10, fmt=trim(fmt3))
'b_max =', b_max
504 write(10, fmt=trim(fmt1))
' '
507 write(10, fmt=trim(fmt3))
'delta_ts0 =', delta_ts0
509 write(10, fmt=trim(fmt3))
'sine_amplit =', sine_amplit
510 write(10, fmt=trim(fmt3))
'sine_period =', sine_period
512 write(10, fmt=trim(fmt1))
'GRIP file = '//grip_temp_file
513 write(10, fmt=trim(fmt3))
'grip_temp_fact =', grip_temp_fact
517 write(10, fmt=trim(fmt3))
'z_sl0 =', z_sl0
519 write(10, fmt=trim(fmt1))
'sea-level file = '//sea_level_file
521 write(10, fmt=trim(fmt1))
' '
524 #if (MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3)
525 write(10, fmt=trim(fmt3))
'z_mar =', z_mar
526 write(10, fmt=trim(fmt1))
' '
527 #elif (MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \
528 || marine_ice_calving==6 || marine_ice_calving==7)
529 write(10, fmt=trim(fmt3))
'fact_z_mar =', fact_z_mar
530 write(10, fmt=trim(fmt1))
' '
531 #elif (MARINE_ICE_FORMATION==2 && MARINE_ICE_CALVING==9)
532 write(10, fmt=trim(fmt3))
'calv_uw_coeff =', calv_uw_coeff
533 write(10, fmt=trim(fmt3))
'r1_calv_uw =', r1_calv_uw
534 write(10, fmt=trim(fmt3))
'r2_calv_uw =', r2_calv_uw
535 write(10, fmt=trim(fmt1))
' '
538 #if (ICE_SHELF_CALVING==2)
539 write(10, fmt=trim(fmt3))
'H_calv =', h_calv
540 write(10, fmt=trim(fmt1))
' '
544 write(10, fmt=trim(fmt1))
'Sediment-mask file = '//mask_sedi_file
545 write(10, fmt=trim(fmt1))
' '
547 write(10, fmt=trim(fmt3))
'c_slide =', c_slide
548 write(10, fmt=trim(fmt3))
'gamma_slide =', gamma_slide
549 write(10, fmt=trim(fmt2))
'p_weert =', p_weert
550 write(10, fmt=trim(fmt2))
'q_weert =', q_weert
552 write(10, fmt=trim(fmt3))
'red_pres_limit_fact =', red_pres_limit_fact
554 write(10, fmt=trim(fmt1))
' '
556 write(10, fmt=trim(fmt3))
'c_slide_sedi =', c_slide_sedi
557 write(10, fmt=trim(fmt3))
'gamma_slide_sedi =', gamma_slide_sedi
558 write(10, fmt=trim(fmt2))
'p_weert_sedi =', p_weert_sedi
559 write(10, fmt=trim(fmt2))
'q_weert_sedi =', q_weert_sedi
560 write(10, fmt=trim(fmt1))
' '
563 write(10, fmt=trim(fmt3))
'q_geo =', q_geo
564 write(10, fmt=trim(fmt1))
' '
567 write(10, fmt=trim(fmt2))
'REBOUND =', rebound
569 write(10, fmt=trim(fmt3))
'frac_llra =', frac_llra
571 #if (REBOUND==1 || REBOUND==2)
572 write(10, fmt=trim(fmt2))
'TIME_LAG_MOD =', time_lag_mod
573 #if (TIME_LAG_MOD==1)
574 write(10, fmt=trim(fmt3))
'time_lag =', time_lag
575 #elif (TIME_LAG_MOD==2)
576 write(10, fmt=trim(fmt1))
'time_lag_file = '//time_lag_file
578 stop
' sico_init: TIME_LAG_MOD must be either 1 or 2!'
582 write(10, fmt=trim(fmt2))
'FLEX_RIG_MOD =', flex_rig_mod
583 #if (FLEX_RIG_MOD==1)
584 write(10, fmt=trim(fmt3))
'flex_rig =', flex_rig
585 #elif (FLEX_RIG_MOD==2)
586 write(10, fmt=trim(fmt1))
'flex_rig_file = '//flex_rig_file
588 stop
' sico_init: FLEX_RIG_MOD must be either 1 or 2!'
591 write(10, fmt=trim(fmt2))
'Q_LITHO =', q_litho
592 write(10, fmt=trim(fmt1))
' '
595 write(10, fmt=trim(fmt3))
'gr_size =', gr_size
596 write(10, fmt=trim(fmt1))
' '
599 write(10, fmt=trim(fmt3))
'sigma_res =', sigma_res
600 write(10, fmt=trim(fmt1))
' '
603 write(10, fmt=trim(fmt2))
'ENHMOD =', enhmod
604 #if (ENHMOD==1 || ENHMOD==2 || ENHMOD==3)
605 write(10, fmt=trim(fmt3))
'enh_fact =', enh_fact
607 #if (ENHMOD==2 || ENHMOD==3)
608 write(10, fmt=trim(fmt3))
'enh_intg =', enh_intg
611 write(10, fmt=trim(fmt3))
'age_trans =', age_trans_0
614 write(10, fmt=trim(fmt3))
'date_trans1 =', date_trans1_0
615 write(10, fmt=trim(fmt3))
'date_trans2 =', date_trans2_0
616 write(10, fmt=trim(fmt3))
'date_trans3 =', date_trans3_0
618 #if (ENHMOD==4 || ENHMOD==5)
619 write(10, fmt=trim(fmt3))
'enh_compr =', enh_compr
620 write(10, fmt=trim(fmt3))
'enh_shear =', enh_shear
622 #if ((ENHMOD==1 || ENHMOD==2 || ENHMOD==3 || ENHMOD==4) && MARGIN==3)
623 write(10, fmt=trim(fmt3))
'enh_shelf =', enh_shelf
625 write(10, fmt=trim(fmt1))
' '
627 write(10, fmt=trim(fmt3))
'numdiff_H_t =', numdiff_h_t
628 write(10, fmt=trim(fmt3))
'tau_cts =', tau_cts
629 write(10, fmt=trim(fmt3))
'H_min =', h_min
630 write(10, fmt=trim(fmt3))
'vh_max =', vh_max
631 write(10, fmt=trim(fmt3))
'hd_min =', hd_min
632 write(10, fmt=trim(fmt3))
'hd_max =', hd_max
633 #if (defined(QBM_MIN))
634 write(10, fmt=trim(fmt3))
'qbm_min =', qbm_min
635 #elif (defined(QB_MIN)) /* obsolete */
636 write(10, fmt=trim(fmt3))
'qb_min =', qb_min
638 #if (defined(QBM_MAX))
639 write(10, fmt=trim(fmt3))
'qbm_max =', qbm_max
640 #elif (defined(QB_MAX)) /* obsolete */
641 write(10, fmt=trim(fmt3))
'qb_max =', qb_max
643 write(10, fmt=trim(fmt3))
'age_min =', age_min
644 write(10, fmt=trim(fmt3))
'age_max =', age_max
645 write(10, fmt=trim(fmt3))
'mean_accum =', mean_accum
647 write(10, fmt=trim(fmt3))
'age_diff =', agediff
649 write(10, fmt=trim(fmt1))
' '
651 write(10, fmt=trim(fmt2))
'GRID =', grid
652 write(10, fmt=trim(fmt2))
'DYNAMICS =', dynamics
653 #if (DYNAMICS==2 && defined(RATIO_SL_THRESH))
654 write(10, fmt=trim(fmt3))
'ratio_sl_thresh =', ratio_sl_thresh
656 write(10, fmt=trim(fmt2))
'CALCMOD =', calcmod
657 #if (CALCMOD==-1 && defined(TEMP_CONST))
658 write(10, fmt=trim(fmt3))
'TEMP_CONST =', temp_const
660 #if (CALCMOD==-1 && defined(AGE_CONST))
661 write(10, fmt=trim(fmt3))
'AGE_CONST =', age_const
663 #if (CALCMOD==1 && defined(CTS_MELTING_FREEZING))
664 write(10, fmt=trim(fmt2))
'CTS_MELTING_FREEZING =', cts_melting_freezing
666 write(10, fmt=trim(fmt2))
'FLOW_LAW =', flow_law
667 write(10, fmt=trim(fmt2))
'FIN_VISC =', fin_visc
668 write(10, fmt=trim(fmt2))
'MARGIN =', margin
670 write(10, fmt=trim(fmt2))
'MARINE_ICE_FORMATION =', marine_ice_formation
671 write(10, fmt=trim(fmt2))
'MARINE_ICE_CALVING =', marine_ice_calving
673 write(10, fmt=trim(fmt2))
'ICE_SHELF_CALVING =', ice_shelf_calving
675 write(10, fmt=trim(fmt2))
'ZS_EVOL =', zs_evol
676 write(10, fmt=trim(fmt2))
'CALCZS =', calczs
677 write(10, fmt=trim(fmt2))
'ADV_HOR =', adv_hor
678 write(10, fmt=trim(fmt2))
'ADV_VERT =', adv_vert
679 write(10, fmt=trim(fmt2))
'TOPOGRAD =', topograd
680 write(10, fmt=trim(fmt2))
'TSURFACE =', tsurface
681 write(10, fmt=trim(fmt2))
'SEA_LEVEL =', sea_level
682 write(10, fmt=trim(fmt2))
'SLIDE_LAW =', slide_law
683 write(10, fmt=trim(fmt2))
'Q_GEO_MOD =', q_geo_mod
684 write(10, fmt=trim(fmt1))
' '
686 #if (defined(OUT_TIMES))
687 write(10, fmt=trim(fmt2))
'OUT_TIMES =', out_times
689 write(10, fmt=trim(fmt2))
'OUTPUT =', output
690 write(10, fmt=trim(fmt2))
'OUTSER =', outser
691 #if (defined(OUTSER_V_A))
692 write(10, fmt=trim(fmt2))
'OUTSER_V_A =', outser_v_a
694 #if (OUTPUT==1 || OUTPUT==2)
695 write(10, fmt=trim(fmt2))
'ERGDAT =', ergdat
697 write(10, fmt=trim(fmt2))
'NETCDF =', netcdf
698 write(10, fmt=trim(fmt1))
' '
699 #if (OUTPUT==2 || OUTPUT==3)
700 write(10, fmt=trim(fmt2))
'n_output =', n_output
702 write(10, fmt=trim(fmt3))
'time_output =', time_output0(n)
704 write(10, fmt=trim(fmt1))
' '
707 write(10, fmt=trim(fmt1))
'Program version and date: '//version//
' / '//date
709 close(10, status=
'keep')
713 year_zero = year_zero*year_sec
714 time_init = time_init0*year_sec
715 time_end = time_end0*year_sec
716 dtime = dtime0*year_sec
717 dtime_temp = dtime_temp0*year_sec
718 dtime_ser = dtime_ser0*year_sec
719 #if (OUTPUT==1 || OUTPUT==3)
720 dtime_out = dtime_out0*year_sec
722 #if (OUTPUT==2 || OUTPUT==3)
724 time_output(n) = time_output0(n)*year_sec
728 if (abs(dtime_temp/dtime-nint(dtime_temp/dtime)) > eps) &
729 stop
' sico_init: dtime_temp must be a multiple of dtime!'
735 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho)
740 open(24, iostat=ios, &
741 file=inpath//
'/'//trim(ch_domain_short)//
'/'//mask_present_file, &
742 recl=2048, status=
'old')
744 if (ios /= 0) stop
' sico_init: Error when opening the mask file!'
746 do n=1, 6;
read(24, fmt=
'(a)') ch_dummy;
end do
749 read(24, fmt=trim(fmt4)) (maske_help(j,i), i=0,imax)
752 close(24, status=
'keep')
756 open(24, iostat=ios, &
757 file=inpath//
'/'//trim(ch_domain_short)//
'/'//mask_sedi_file, &
758 recl=2048, status=
'old')
760 if (ios /= 0) stop
' sico_init: Error when opening the rock/sediment mask file!'
762 do n=1, 6;
read(24, fmt=
'(a)') ch_dummy;
end do
765 read(24, fmt=trim(fmt4)) (maske_sedi(j,i), i=0,imax)
768 close(24, status=
'keep')
774 open(21, iostat=ios, &
775 file=inpath//
'/general/'//grip_temp_file, &
777 if (ios /= 0) stop
' sico_init: Error when opening the data file for delta_ts!'
779 read(21, fmt=*) ch_dummy, grip_time_min, grip_time_stp, grip_time_max
781 if (ch_dummy /=
'#')
then
782 write(6, fmt=*)
' sico_init: grip_time_min, grip_time_stp, grip_time_max'
783 write(6, fmt=*)
' not defined in data file for delta_ts!'
787 ndata_grip = (grip_time_max-grip_time_min)/grip_time_stp
789 allocate(griptemp(0:ndata_grip))
792 read(21, fmt=*) d_dummy, griptemp(n)
795 close(21, status=
'keep')
803 open(21, iostat=ios, &
804 file=inpath//
'/general/'//sea_level_file, &
806 if (ios /= 0) stop
' sico_init: Error when opening the data file for z_sl!'
808 read(21, fmt=*) ch_dummy, specmap_time_min, specmap_time_stp, specmap_time_max
810 if (ch_dummy /=
'#')
then
811 write(6, fmt=*)
' sico_init: specmap_time_min, specmap_time_stp, specmap_time_max'
812 write(6, fmt=*)
' not defined in data file for z_sl!'
816 ndata_specmap = (specmap_time_max-specmap_time_min)/specmap_time_stp
818 allocate(specmap_zsl(0:ndata_specmap))
820 do n=0, ndata_specmap
821 read(21, fmt=*) d_dummy, specmap_zsl(n)
824 close(21, status=
'keep')
836 q_geo(j,i) = q_geo *1.0e-03_dp
842 stop
' sico_init: Option Q_GEO_MOD==2 not available for this application!'
849 #if (REBOUND==1 || REBOUND==2)
851 #if (TIME_LAG_MOD==1)
853 time_lag_asth = time_lag*year_sec
855 #elif (TIME_LAG_MOD==2)
857 open(29, iostat=ios, &
858 file=inpath//
'/'//trim(ch_domain_short)//
'/'//time_lag_file, &
859 recl=8192, status=
'old')
861 if (ios /= 0) stop
' sico_init: Error when opening the time-lag file!'
863 do n=1, 6;
read(29, fmt=
'(a)') ch_dummy;
end do
866 read(29, fmt=*) (time_lag_asth(j,i), i=0,imax)
869 close(29, status=
'keep')
871 time_lag_asth = time_lag_asth*year_sec
877 time_lag_asth = 0.0_dp
885 #if (FLEX_RIG_MOD==1)
887 flex_rig_lith = flex_rig
889 #elif (FLEX_RIG_MOD==2)
891 open(29, iostat=ios, &
892 file=inpath//
'/'//trim(ch_domain_short)//
'/'//flex_rig_file, &
893 recl=8192, status=
'old')
895 if (ios /= 0) stop
' sico_init: Error when opening the flex-rig file!'
897 do n=1, 6;
read(29, fmt=
'(a)') ch_dummy;
end do
900 read(29, fmt=*) (flex_rig_lith(j,i), i=0,imax)
903 close(29, status=
'keep')
907 #elif (REBOUND==0 || REBOUND==1)
909 flex_rig_lith = 0.0_dp
923 call
boundary(time_init, dtime, dxi, deta, &
924 delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
945 call
boundary(time_init, dtime, dxi, deta, &
946 delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
968 stop
' sico_init: Parameter NETCDF must be either 1 or 2!'
971 call
boundary(time_init, dtime, dxi, deta, &
972 delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
974 q_b_tot = q_bm + q_tld
980 flag_grounding_line_1 = .false.
981 flag_grounding_line_2 = .false.
983 flag_calving_front_1 = .false.
984 flag_calving_front_2 = .false.
986 #if (MARGIN==1 || MARGIN==2)
995 if ( (maske(j,i)==0_i2b) &
997 ( (maske(j,i+1)==3_i2b) &
998 .or.(maske(j,i-1)==3_i2b) &
999 .or.(maske(j+1,i)==3_i2b) &
1000 .or.(maske(j-1,i)==3_i2b) ) &
1002 flag_grounding_line_1(j,i) = .true.
1004 if ( (maske(j,i)==3_i2b) &
1006 ( (maske(j,i+1)==0_i2b) &
1007 .or.(maske(j,i-1)==0_i2b) &
1008 .or.(maske(j+1,i)==0_i2b) &
1009 .or.(maske(j-1,i)==0_i2b) ) &
1011 flag_grounding_line_2(j,i) = .true.
1013 if ( (maske(j,i)==3_i2b) &
1015 ( (maske(j,i+1)==2_i2b) &
1016 .or.(maske(j,i-1)==2_i2b) &
1017 .or.(maske(j+1,i)==2_i2b) &
1018 .or.(maske(j-1,i)==2_i2b) ) &
1020 flag_calving_front_1(j,i) = .true.
1022 if ( (maske(j,i)==2_i2b) &
1024 ( (maske(j,i+1)==3_i2b) &
1025 .or.(maske(j,i-1)==3_i2b) &
1026 .or.(maske(j+1,i)==3_i2b) &
1027 .or.(maske(j-1,i)==3_i2b) ) &
1029 flag_calving_front_2(j,i) = .true.
1037 if ( (maske(j,i)==2_i2b) &
1038 .and. (maske(j+1,i)==3_i2b) &
1040 flag_calving_front_2(j,i) = .true.
1043 if ( (maske(j,i)==2_i2b) &
1044 .and. (maske(j-1,i)==3_i2b) &
1046 flag_calving_front_2(j,i) = .true.
1053 if ( (maske(j,i)==2_i2b) &
1054 .and. (maske(j,i+1)==3_i2b) &
1056 flag_calving_front_2(j,i) = .true.
1059 if ( (maske(j,i)==2_i2b) &
1060 .and. (maske(j,i-1)==3_i2b) &
1062 flag_calving_front_2(j,i) = .true.
1067 stop
' sico_init: MARGIN must be either 1, 2 or 3!'
1072 #if (GRID==0 || GRID==1)
1076 dist_dxdy(jr,ir) = sqrt( (
real(ir,dp)*dxi)**2 + (
real(jr,dp)*deta)**2 )
1087 #if (DYNAMICS==1 || DYNAMICS==2)
1092 #if (MARGIN==3 || DYNAMICS==2)
1093 call
calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
1108 stop
' sico_init: DYNAMICS must be either 0, 1 or 2!'
1111 call
calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
1115 #if (CALCMOD==0 || CALCMOD==-1)
1125 enth_t(kt,j,i) = enth_c(0,j,i)
1140 if ( (maske(j,i) == 0_i2b).and.(n_cts(j,i) == 1_i2b) )
then
1146 enth_t(kt,j,i) = enth_c(0,j,i)
1153 #elif (CALCMOD==2 || CALCMOD==3)
1163 enth_t(kt,j,i) = enth_c(0,j,i)
1171 stop
' sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!'
1179 open(12, iostat=ios, &
1180 file=outpath//
'/'//trim(runname)//
'.ser', &
1182 if (ios /= 0) stop
' sico_init: Error when opening the ser file!'
1187 #if (!defined(OUTSER_V_A) || OUTSER_V_A==1)
1189 1102
format(
' t(a) D_Ts(deg C) z_sl(m)',/, &
1190 ' H_max(m) zs_max(m) V_g(m^3)', &
1191 ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1192 ' A_g(m^2) A_t(m^2) V_bm(m^3/a)', &
1193 ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1194 1103
format(
'----------------------------------------------------', &
1195 '---------------------------------------')
1197 #elif (OUTSER_V_A==2)
1199 1102
format(
' t(a) D_Ts(deg C) z_sl(m)',/, &
1200 ' V(m^3) V_g(m^3) V_f(m^3)', &
1201 ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1202 ' H_max(m) zs_max(m)', &
1203 ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1204 ' A_t(m^2) V_bm(m^3/a)', &
1205 ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1206 1103
format(
'----------------------------------------------------', &
1207 '---------------------------------------')
1210 stop
' sico_init: OUTSER_V_A must be either 1 or 2!'
1215 open(15, iostat=ios, &
1216 file=outpath//
'/'//trim(runname)//
'.sed', &
1218 if (ios /= 0) stop
' sico_init: Error when opening the sed file!'
1223 1108
format(
' t(a) D_Ts(deg C) z_sl(m)',/, &
1224 ' H_ave(m) Tbh_ave(C) Atb(m^2)')
1225 1109
format(
'----------------------------------------------------')
1231 open(13, iostat=ios, &
1232 file=outpath//
'/'//trim(runname)//
'.thk', &
1234 if (ios /= 0) stop
' sico_init: Error when opening the thk file!'
1239 1104
format(
' t(a) D_Ts(deg C) z_sl(m)',/, &
1240 ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
1241 ' IMAX+1 values [i = 0 (1) IMAX] in each record')
1242 1105
format(
'----------------------------------------------------')
1252 allocate(lambda_core(n_core), phi_core(n_core), &
1253 x_core(n_core), y_core(n_core))
1255 lambda_core(1) = 0.0_dp
1256 phi_core(1) = 0.0_dp
1257 x_core(1) = 3900.0_dp *1.0e+03_dp
1258 y_core(1) = 2000.0_dp *1.0e+03_dp
1260 lambda_core(2) = 0.0_dp
1261 phi_core(2) = 0.0_dp
1262 x_core(2) = 3800.0_dp *1.0e+03_dp
1263 y_core(2) = 2000.0_dp *1.0e+03_dp
1265 lambda_core(3) = 0.0_dp
1266 phi_core(3) = 0.0_dp
1267 x_core(3) = 3700.0_dp *1.0e+03_dp
1268 y_core(3) = 2000.0_dp *1.0e+03_dp
1270 lambda_core(4) = 0.0_dp
1271 phi_core(4) = 0.0_dp
1272 x_core(4) = 3500.0_dp *1.0e+03_dp
1273 y_core(4) = 2000.0_dp *1.0e+03_dp
1275 lambda_core(5) = 0.0_dp
1276 phi_core(5) = 0.0_dp
1277 x_core(5) = 3200.0_dp *1.0e+03_dp
1278 y_core(5) = 2000.0_dp *1.0e+03_dp
1280 lambda_core(6) = 0.0_dp
1281 phi_core(6) = 0.0_dp
1282 x_core(6) = 2900.0_dp *1.0e+03_dp
1283 y_core(6) = 2000.0_dp *1.0e+03_dp
1285 lambda_core(7) = 0.0_dp
1286 phi_core(7) = 0.0_dp
1287 x_core(7) = 2600.0_dp *1.0e+03_dp
1288 y_core(7) = 2000.0_dp *1.0e+03_dp
1290 open(14, iostat=ios, &
1291 file=outpath//
'/'//trim(runname)//
'.core', &
1293 if (ios /= 0) stop
' sico_init: Error when opening the core file!'
1295 if (forcing_flag == 1)
then
1300 1106
format(
' t(a) D_Ts(C) z_sl(m)',/, &
1301 ' H_P1(m) H_P2(m) H_P3(m) H_P4(m)', &
1302 ' H_P5(m) H_P6(m) H_P7(m)',/, &
1303 ' v_P1(m/a) v_P2(m/a) v_P3(m/a) v_P4(m/a)', &
1304 ' v_P5(m/a) v_P6(m/a) v_P7(m/a)',/, &
1305 ' T_P1(C) T_P2(C) T_P3(C) T_P4(C)', &
1306 ' T_P5(C) T_P6(C) T_P7(C)',/, &
1307 ' Rb_P1(W/m2) Rb_P2(W/m2) Rb_P3(W/m2) Rb_P4(W/m2)', &
1308 ' Rb_P5(W/m2) Rb_P6(W/m2) Rb_P7(W/m2)')
1309 1107
format(
'-----------------------------------------------------------------', &
1310 '---------------------------------------')
1321 flag_3d_output = .false.
1323 flag_3d_output = .true.
1325 stop
' sico_init: ERGDAT must be either 0 or 1!'
1329 call
output1(runname, time_init, delta_ts, glac_index, z_sl, &
1330 flag_3d_output, ndat2d, ndat3d)
1332 call
output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1333 flag_3d_output, ndat2d, ndat3d)
1335 stop
' sico_init: Parameter NETCDF must be either 1 or 2!'
1340 if (time_output(1) <= time_init+eps)
then
1343 flag_3d_output = .false.
1345 flag_3d_output = .true.
1347 stop
' sico_init: ERGDAT must be either 0 or 1!'
1351 call
output1(runname, time_init, delta_ts, glac_index, z_sl, &
1352 flag_3d_output, ndat2d, ndat3d)
1354 call
output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1355 flag_3d_output, ndat2d, ndat3d)
1357 stop
' sico_init: Parameter NETCDF must be either 1 or 2!'
1364 flag_3d_output = .false.
1367 call
output1(runname, time_init, delta_ts, glac_index, z_sl, &
1368 flag_3d_output, ndat2d, ndat3d)
1370 call
output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1371 flag_3d_output, ndat2d, ndat3d)
1373 stop
' sico_init: Parameter NETCDF must be either 1 or 2!'
1376 if (time_output(1) <= time_init+eps)
then
1378 flag_3d_output = .true.
1381 call
output1(runname, time_init, delta_ts, glac_index, z_sl, &
1382 flag_3d_output, ndat2d, ndat3d)
1384 call
output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1385 flag_3d_output, ndat2d, ndat3d)
1387 stop
' sico_init: Parameter NETCDF must be either 1 or 2!'
1393 stop
' sico_init: OUTPUT must be either 1, 2 or 3!'
1396 call
output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1398 call
output3(time_init, delta_ts, glac_index, z_sl)
1401 call
output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
subroutine calc_vz_ssa(z_sl, dxi, deta, dzeta_c)
Computation of the vertical velocity vz in the shallow shelf approximation.
subroutine calc_vxy_static()
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy etc. for static ice.
subroutine calc_vxy_sia(dzeta_c, dzeta_t)
Computation of the shear stresses txz, tyz, the effective shear stress sigma, the depth-averaged flui...
subroutine topography1(dxi, deta)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
subroutine, public ice_mat_eqs_pars(RF_table, R_T_val, KAPPA_table, C_table, n_temp_min, n_temp_max, RHO_I_val, RHO_C_val, KAPPA_C_val, C_C_val)
Setting of required physical parameters.
subroutine output2(time, dxi, deta, delta_ts, glac_index, z_sl)
Writing of time-series data on file in ASCII format.
subroutine phys_para()
Reading of physical parameters.
Declarations of kind types for SICOPOLIS.
subroutine topography3(dxi, deta, z_sl, anfdatname)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
subroutine output3(time, delta_ts, glac_index, z_sl)
Writing of time-series data of the ice thickness field on file in ASCII format.
subroutine topography3_nc(dxi, deta, z_sl, anfdatname)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
subroutine calc_temp_melt()
Computation of the melting temperatures.
subroutine output_nc(runname, time, delta_ts, glac_index, z_sl, flag_3d_output, ndat2d, ndat3d)
Writing of time-slice files in NetCDF format.
subroutine 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).
subroutine topography2(dxi, deta)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
subroutine calc_vz_static()
Computation of the vertical velocity vz for static ice.
subroutine calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
Computation of all components of the strain-rate tensor, the full effective strain rate and the shear...
subroutine calc_enhance(time)
Computation of the flow enhancement factor.
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.
subroutine calc_vxy_b_sia(time, z_sl)
Computation of the basal horizontal velocity vx_b, vy_b in the shallow ice approximation.
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)
Initialisations for SICOPOLIS.
subroutine boundary(time, dtime, dxi, deta, delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
Computation of the surface temperature (must be less than 0 deg C!) and of the accumulation-ablation ...
subroutine init_temp_water_age()
Initial temperature, water content and age.
Material quantities of ice: Rate factor, heat conductivity, specific heat (heat capacity).
subroutine 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.
Declarations of global variables for SICOPOLIS.
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_...
Conversion from temperature (temp) and water content (omega) to enthalpy (enth) and vice versa...
subroutine calc_vz_sia(dxi, deta, dzeta_c, dzeta_t)
Computation of the vertical velocity vz in the shallow ice approximation.
subroutine output1(runname, time, delta_ts, glac_index, z_sl, flag_3d_output, ndat2d, ndat3d)
Writing of time-slice files in binary format.
real(dp) function, public enth_fct_temp_omega(temp_val, omega_val)
Enthalpy as a function of temperature and water content.