50 subroutine sico_init(delta_ts, glac_index, &
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, &
57 ndat2d, ndat3d, n_output, &
79 integer(i4b),
intent(out) :: ii((imax+1)*(jmax+1)), &
80 jj((imax+1)*(jmax+1)), &
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, &
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
93 integer(i4b) :: i, j, kc, kt, kr, m, n, ir, jr
94 integer(i4b) :: ios, ios1, ios2, ios3, ios4
96 real(dp) :: dtime0, dtime_temp0, dtime_wss0, dtime_out0, dtime_ser0
97 real(dp) :: time_init0, time_end0, time_output0(100)
99 character(len=100) :: anfdatname
100 character(len=256) :: filename_with_path
101 character(len=256) :: shell_command
102 character :: ch_dummy
103 logical :: flag_init_output, flag_3d_output
105 character(len=64),
parameter :: fmt1 =
'(a)', &
110 character(len= 8) :: ch_imax
111 character(len=128) :: fmt4
113 write(ch_imax, fmt=
'(i8)') imax
114 write(fmt4, fmt=
'(a)')
'('//trim(adjustl(ch_imax))//
'(i1),i1)' 116 write(unit=6, fmt=
'(a)')
' ' 117 write(unit=6, fmt=
'(a)')
' -------- sico_init --------' 118 write(unit=6, fmt=
'(a)')
' ' 130 #elif (defined(EMTP2SGE)) 138 #elif (defined(NHEM)) 142 #elif (defined(SCAND)) 146 #elif (defined(TIBET)) 150 #elif (defined(NMARS)) 154 #elif (defined(SMARS)) 166 stop
' >>> sico_init: No valid domain specified!' 187 #if (CALCTHK==3 || CALCTHK==6 || MARGIN==3 || DYNAMICS==2) 188 call lis_initialize(ierr)
212 if ( trim(version) /= trim(sico_version) ) &
213 stop
' >>> sico_init: SICOPOLIS version not compatible with header file!' 217 #if (!defined(DYNAMICS)) 218 stop
' >>> sico_init: DYNAMICS not defined in the header file!' 221 #if (!defined(CALCMOD)) 222 stop
' >>> sico_init: CALCMOD not defined in the header file!' 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!' 236 if (approx_equal(dx, 50.0_dp,
eps_sp_dp))
then 238 if ((imax /= 80).or.(jmax /= 80))
then 239 stop
' >>> sico_init: IMAX and/or JMAX wrong!' 242 else if (approx_equal(dx, 25.0_dp,
eps_sp_dp))
then 244 if ((imax /= 160).or.(jmax /= 160))
then 245 stop
' >>> sico_init: IMAX and/or JMAX wrong!' 249 stop
' >>> sico_init: DX wrong!' 254 stop
' >>> sico_init: GRID==1 not allowed for this application!' 258 stop
' >>> sico_init: GRID==2 not allowed for this application!' 266 #if (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1) 269 write(6, fmt=
'(a)')
' >>> sico_init: For options CALCMOD==0, 2, 3 or -1,' 270 write(6, fmt=
'(a)')
' the separate kt domain is redundant.' 271 write(6, fmt=
'(a)')
' Therefore, consider setting KTMAX to 2.' 272 write(6, fmt=
'(a)')
' ' 281 stop
' >>> sico_init: Option ADV_HOR==1 (central differences) not defined!' 289 #if ((MARGIN==3 && DYNAMICS==1) || DYNAMICS==2) /* requires SSA or SStA */ 291 write(6, fmt=
'(a)')
' >>> sico_init: Distortion correction not implemented' 292 write(6, fmt=
'(a)')
' for the shallow shelf approximation (SSA)' 293 write(6, fmt=
'(a)')
' or the shelfy stream approximation (SStA)' 294 write(6, fmt=
'(a)')
' -> GRID==0 required!' 306 dtime_temp0 = dtime_temp0
310 dzeta_c = 1.0_dp/
real(kcmax,dp)
311 dzeta_t = 1.0_dp/
real(ktmax,dp)
312 dzeta_r = 1.0_dp/
real(krmax,dp)
321 if (deform >=
eps)
then 401 if ((i <= imax).and.(j <= jmax))
then 413 anfdatname = anfdatname
415 #if (defined(YEAR_ZERO)) 421 time_init0 = time_init0
422 time_end0 = time_end0
423 dtime_ser0 = dtime_ser0
425 #if (OUTPUT==1 || OUTPUT==3) 426 dtime_out0 = dtime_out0
429 #if (OUTPUT==2 || OUTPUT==3) 431 time_output0( 1) = time_out0_01
432 time_output0( 2) = time_out0_02
433 time_output0( 3) = time_out0_03
434 time_output0( 4) = time_out0_04
435 time_output0( 5) = time_out0_05
436 time_output0( 6) = time_out0_06
437 time_output0( 7) = time_out0_07
438 time_output0( 8) = time_out0_08
439 time_output0( 9) = time_out0_09
440 time_output0(10) = time_out0_10
441 time_output0(11) = time_out0_11
442 time_output0(12) = time_out0_12
443 time_output0(13) = time_out0_13
444 time_output0(14) = time_out0_14
445 time_output0(15) = time_out0_15
446 time_output0(16) = time_out0_16
447 time_output0(17) = time_out0_17
448 time_output0(18) = time_out0_18
449 time_output0(19) = time_out0_19
450 time_output0(20) = time_out0_20
455 shell_command =
'if [ ! -d' 456 shell_command = trim(shell_command)//
' '//outpath
457 shell_command = trim(shell_command)//
' '//
'] ; then mkdir' 458 shell_command = trim(shell_command)//
' '//outpath
459 shell_command = trim(shell_command)//
' '//
'; fi' 460 call system(trim(shell_command))
463 filename_with_path = trim(outpath)//
'/'//trim(runname)//
'.log' 465 open(10, iostat=ios, file=trim(filename_with_path), status=
'new')
467 if (ios /= 0) stop
' >>> sico_init: Error when opening the log file!' 469 write(10, fmt=trim(fmt1))
'Computational domain:' 471 write(10, fmt=trim(fmt1))
' ' 473 write(10, fmt=trim(fmt2a))
'GRID = ', grid
474 write(10, fmt=trim(fmt1))
' ' 476 write(10, fmt=trim(fmt2))
'imax =', imax
477 write(10, fmt=trim(fmt2))
'jmax =', jmax
478 write(10, fmt=trim(fmt2))
'kcmax =', kcmax
479 write(10, fmt=trim(fmt2))
'ktmax =', ktmax
480 write(10, fmt=trim(fmt2))
'krmax =', krmax
481 write(10, fmt=trim(fmt1))
' ' 483 write(10, fmt=trim(fmt3))
'a =',
aa 484 write(10, fmt=trim(fmt1))
' ' 486 #if (GRID==0 || GRID==1) 487 write(10, fmt=trim(fmt3))
'x0 =', x0
488 write(10, fmt=trim(fmt3))
'y0 =', y0
489 write(10, fmt=trim(fmt3))
'dx =', dx
491 write(10, fmt=trim(fmt3))
'dlambda =', dlambda
492 write(10, fmt=trim(fmt3))
'dphi =', dphi
494 write(10, fmt=trim(fmt1))
' ' 496 write(10, fmt=trim(fmt3))
'year_zero =',
year_zero 497 write(10, fmt=trim(fmt3))
'time_init =', time_init0
498 write(10, fmt=trim(fmt3))
'time_end =', time_end0
499 write(10, fmt=trim(fmt3))
'dtime =', dtime0
500 write(10, fmt=trim(fmt3))
'dtime_temp =', dtime_temp0
501 write(10, fmt=trim(fmt1))
' ' 503 write(10, fmt=trim(fmt2))
'ANF_DAT =', anf_dat
504 write(10, fmt=trim(fmt1))
'mask_present file = '//mask_present_file
505 #if (ANF_DAT==1 && defined(TEMP_INIT)) 506 write(10, fmt=trim(fmt2))
'TEMP_INIT =', temp_init
509 write(10, fmt=trim(fmt1))
'Initial-value file = '//anfdatname
510 write(10, fmt=trim(fmt1))
'Path to initial-value file = '//anfdatpath
512 write(10, fmt=trim(fmt1))
' ' 514 write(10, fmt=trim(fmt1))
'Physical-parameter file = '//phys_para_file
515 write(10, fmt=trim(fmt1))
' ' 517 #if (CALCTHK==2 || CALCTHK==3 || CALCTHK==5 || CALCTHK==6) 518 write(10, fmt=trim(fmt3))
'ovi_weight =', ovi_weight
519 #if (CALCTHK==2 || CALCTHK==5) 520 write(10, fmt=trim(fmt3))
'omega_sor =', omega_sor
522 write(10, fmt=trim(fmt2a))
'iter_max_sor = ', iter_max_sor
525 write(10, fmt=trim(fmt1))
' ' 528 write(10, fmt=trim(fmt3))
'temp_min =',
temp_min 529 write(10, fmt=trim(fmt3))
's_t =',
s_t 530 write(10, fmt=trim(fmt3))
'x_hat =',
x_hat 531 write(10, fmt=trim(fmt3))
'y_hat =',
y_hat 532 write(10, fmt=trim(fmt3))
'rad =',
rad 533 write(10, fmt=trim(fmt3))
'b_min =',
b_min 534 write(10, fmt=trim(fmt3))
'b_max =',
b_max 535 write(10, fmt=trim(fmt1))
' ' 538 write(10, fmt=trim(fmt3))
'delta_ts0 =', delta_ts0
540 write(10, fmt=trim(fmt3))
'sine_amplit =', sine_amplit
541 write(10, fmt=trim(fmt3))
'sine_period =', sine_period
543 write(10, fmt=trim(fmt1))
'GRIP file = '//grip_temp_file
544 write(10, fmt=trim(fmt3))
'grip_temp_fact =', grip_temp_fact
547 write(10, fmt=trim(fmt2a))
'SEA_LEVEL = ', sea_level
549 write(10, fmt=trim(fmt3))
'z_sl0 =', z_sl0
551 write(10, fmt=trim(fmt1))
'sea-level file = '//sea_level_file
553 write(10, fmt=trim(fmt1))
' ' 556 #if (MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3) 557 write(10, fmt=trim(fmt3))
'z_mar =', z_mar
558 write(10, fmt=trim(fmt1))
' ' 559 #elif (MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \ 560 || marine_ice_calving==6 || marine_ice_calving==7)
561 write(10, fmt=trim(fmt3))
'fact_z_mar =', fact_z_mar
562 write(10, fmt=trim(fmt1))
' ' 563 #elif (MARINE_ICE_FORMATION==2 && MARINE_ICE_CALVING==9) 564 write(10, fmt=trim(fmt3))
'calv_uw_coeff =', calv_uw_coeff
565 write(10, fmt=trim(fmt3))
'r1_calv_uw =', r1_calv_uw
566 write(10, fmt=trim(fmt3))
'r2_calv_uw =', r2_calv_uw
567 write(10, fmt=trim(fmt1))
' ' 570 #if (ICE_SHELF_CALVING==2) 571 write(10, fmt=trim(fmt3))
'H_calv =', h_calv
572 write(10, fmt=trim(fmt1))
' ' 576 write(10, fmt=trim(fmt1))
'Sediment-mask file = '//mask_sedi_file
577 write(10, fmt=trim(fmt1))
' ' 579 #if (defined(BASAL_HYDROLOGY)) 580 write(10, fmt=trim(fmt2a))
'BASAL_HYDROLOGY = ', basal_hydrology
582 write(10, fmt=trim(fmt2a))
'SLIDE_LAW = ', slide_law
583 write(10, fmt=trim(fmt3))
'c_slide =',
c_slide 584 write(10, fmt=trim(fmt3))
'gamma_slide =', gamma_slide
585 write(10, fmt=trim(fmt2a))
'p_weert = ',
p_weert 586 write(10, fmt=trim(fmt2a))
'q_weert = ',
q_weert 587 #if (defined(TIME_RAMP_UP_SLIDE)) 588 write(10, fmt=trim(fmt3))
'time_ramp_up_slide =', time_ramp_up_slide
591 write(10, fmt=trim(fmt3))
'red_pres_limit_fact =', red_pres_limit_fact
593 #if (BASAL_HYDROLOGY==1) 594 write(10, fmt=trim(fmt3))
'c_Hw_slide =', c_hw_slide
595 write(10, fmt=trim(fmt3))
'Hw0_slide =', hw0_slide
597 write(10, fmt=trim(fmt1))
' ' 599 write(10, fmt=trim(fmt3))
'c_slide_sedi =', c_slide_sedi
600 write(10, fmt=trim(fmt3))
'gamma_slide_sedi =', gamma_slide_sedi
601 write(10, fmt=trim(fmt2))
'p_weert_sedi =', p_weert_sedi
602 write(10, fmt=trim(fmt2))
'q_weert_sedi =', q_weert_sedi
603 write(10, fmt=trim(fmt1))
' ' 605 write(10, fmt=trim(fmt2a))
'Q_GEO_MOD = ', q_geo_mod
607 write(10, fmt=trim(fmt3))
'q_geo =',
q_geo 609 write(10, fmt=trim(fmt1))
' ' 611 write(10, fmt=trim(fmt2))
'REBOUND =', rebound
613 write(10, fmt=trim(fmt3))
'frac_llra =', frac_llra
615 #if (REBOUND==1 || REBOUND==2) 616 write(10, fmt=trim(fmt2))
'TIME_LAG_MOD =', time_lag_mod
617 #if (TIME_LAG_MOD==1) 618 write(10, fmt=trim(fmt3))
'time_lag =', time_lag
619 #elif (TIME_LAG_MOD==2) 620 write(10, fmt=trim(fmt1))
'time_lag_file = '//time_lag_file
622 stop
' >>> sico_init: TIME_LAG_MOD must be either 1 or 2!' 626 write(10, fmt=trim(fmt2))
'FLEX_RIG_MOD =', flex_rig_mod
627 #if (FLEX_RIG_MOD==1) 628 write(10, fmt=trim(fmt3))
'flex_rig =', flex_rig
629 #elif (FLEX_RIG_MOD==2) 630 write(10, fmt=trim(fmt1))
'flex_rig_file = '//flex_rig_file
632 stop
' >>> sico_init: FLEX_RIG_MOD must be either 1 or 2!' 635 write(10, fmt=trim(fmt2))
'Q_LITHO =', q_litho
636 write(10, fmt=trim(fmt1))
' ' 639 write(10, fmt=trim(fmt3))
'gr_size =', gr_size
640 write(10, fmt=trim(fmt1))
' ' 643 write(10, fmt=trim(fmt3))
'sigma_res =', sigma_res
644 write(10, fmt=trim(fmt1))
' ' 647 write(10, fmt=trim(fmt2))
'ENHMOD =', enhmod
648 #if (ENHMOD==1 || ENHMOD==2 || ENHMOD==3) 649 write(10, fmt=trim(fmt3))
'enh_fact =', enh_fact
651 #if (ENHMOD==2 || ENHMOD==3) 652 write(10, fmt=trim(fmt3))
'enh_intg =', enh_intg
655 write(10, fmt=trim(fmt3))
'age_trans =', age_trans_0
658 write(10, fmt=trim(fmt3))
'date_trans1 =', date_trans1_0
659 write(10, fmt=trim(fmt3))
'date_trans2 =', date_trans2_0
660 write(10, fmt=trim(fmt3))
'date_trans3 =', date_trans3_0
662 #if (ENHMOD==4 || ENHMOD==5) 663 write(10, fmt=trim(fmt3))
'enh_compr =', enh_compr
664 write(10, fmt=trim(fmt3))
'enh_shear =', enh_shear
666 #if ((ENHMOD==1 || ENHMOD==2 || ENHMOD==3 || ENHMOD==4) && MARGIN==3) 667 write(10, fmt=trim(fmt3))
'enh_shelf =', enh_shelf
669 write(10, fmt=trim(fmt1))
' ' 671 write(10, fmt=trim(fmt3))
'numdiff_H_t =', numdiff_h_t
672 write(10, fmt=trim(fmt3))
'tau_cts =', tau_cts
673 write(10, fmt=trim(fmt3))
'vh_max =', vh_max
674 write(10, fmt=trim(fmt3))
'hd_min =', hd_min
675 write(10, fmt=trim(fmt3))
'hd_max =', hd_max
676 #if (defined(QBM_MIN)) 677 write(10, fmt=trim(fmt3))
'qbm_min =', qbm_min
678 #elif (defined(QB_MIN)) /* obsolete */ 679 write(10, fmt=trim(fmt3))
'qb_min =', qb_min
681 #if (defined(QBM_MAX)) 682 write(10, fmt=trim(fmt3))
'qbm_max =', qbm_max
683 #elif (defined(QB_MAX)) /* obsolete */ 684 write(10, fmt=trim(fmt3))
'qb_max =', qb_max
686 write(10, fmt=trim(fmt3))
'age_min =', age_min
687 write(10, fmt=trim(fmt3))
'age_max =', age_max
688 write(10, fmt=trim(fmt3))
'mean_accum =', mean_accum
690 write(10, fmt=trim(fmt3))
'age_diff =', agediff
692 write(10, fmt=trim(fmt1))
' ' 694 write(10, fmt=trim(fmt2a))
'DYNAMICS = ', dynamics
695 #if ((DYNAMICS==1 && MARGIN==3) || DYNAMICS==2) 696 #if (defined(LIS_OPTS)) 697 write(10, fmt=trim(fmt1))
'lis_opts = '//lis_opts
699 #if (defined(N_ITER_SSA)) 700 write(10, fmt=trim(fmt2a))
'n_iter_ssa = ', n_iter_ssa
702 #if (defined(RELAX_FACT_SSA)) 703 write(10, fmt=trim(fmt3))
'relax_fact_ssa =', relax_fact_ssa
706 #if (DYNAMICS==2 && defined(RATIO_SL_THRESH)) 707 write(10, fmt=trim(fmt3))
'ratio_sl_thresh =', ratio_sl_thresh
709 #if (DYNAMICS==2 && defined(SSTA_SIA_WEIGH_FCT)) 710 write(10, fmt=trim(fmt2a))
'SSTA_SIA_WEIGH_FCT = ', ssta_sia_weigh_fct
712 write(10, fmt=trim(fmt1))
' ' 714 write(10, fmt=trim(fmt2))
'CALCMOD =', calcmod
715 #if (CALCMOD==-1 && defined(TEMP_CONST)) 716 write(10, fmt=trim(fmt3))
'TEMP_CONST =', temp_const
718 #if (CALCMOD==-1 && defined(AGE_CONST)) 719 write(10, fmt=trim(fmt3))
'AGE_CONST =', age_const
721 #if (CALCMOD==1 && defined(CTS_MELTING_FREEZING)) 722 write(10, fmt=trim(fmt2))
'CTS_MELTING_FREEZING =', cts_melting_freezing
724 write(10, fmt=trim(fmt2))
'FLOW_LAW =', flow_law
725 write(10, fmt=trim(fmt2))
'FIN_VISC =', fin_visc
726 write(10, fmt=trim(fmt2))
'MARGIN =', margin
728 write(10, fmt=trim(fmt2))
'MARINE_ICE_FORMATION =', marine_ice_formation
729 write(10, fmt=trim(fmt2))
'MARINE_ICE_CALVING =', marine_ice_calving
731 write(10, fmt=trim(fmt2))
'ICE_SHELF_CALVING =', ice_shelf_calving
733 #if (defined(THK_EVOL)) 734 write(10, fmt=trim(fmt2))
'THK_EVOL =', thk_evol
736 stop
' >>> sico_init: Define THK_EVOL in header file!' 738 #if (defined(CALCTHK)) 739 write(10, fmt=trim(fmt2))
'CALCTHK =', calcthk
741 stop
' >>> sico_init: Define CALCTHK in header file!' 743 write(10, fmt=trim(fmt2))
'ADV_HOR =', adv_hor
744 write(10, fmt=trim(fmt2))
'ADV_VERT =', adv_vert
745 write(10, fmt=trim(fmt2))
'TOPOGRAD =', topograd
746 write(10, fmt=trim(fmt2))
'TSURFACE =', tsurface
747 #if (defined(MB_ACCOUNT)) 748 write(10, fmt=trim(fmt2a))
'MB_ACCOUNT = ', mb_account
750 write(10, fmt=trim(fmt1))
' ' 752 #if (defined(OUT_TIMES)) 753 write(10, fmt=trim(fmt2a))
'OUT_TIMES = ', out_times
755 #if (defined(OUTPUT_INIT)) 756 write(10, fmt=trim(fmt2a))
'OUTPUT_INIT = ', output_init
758 write(10, fmt=trim(fmt2a))
'OUTPUT = ', output
759 #if (OUTPUT==1 || OUTPUT==3) 760 write(10, fmt=trim(fmt3))
'dtime_out =' , dtime_out0
762 write(10, fmt=trim(fmt3))
'dtime_ser =' , dtime_ser0
763 #if (OUTPUT==1 || OUTPUT==2) 764 write(10, fmt=trim(fmt2a))
'ERGDAT = ', ergdat
766 write(10, fmt=trim(fmt2a))
'NETCDF = ', netcdf
767 #if (OUTPUT==2 || OUTPUT==3) 768 write(10, fmt=trim(fmt2a))
'n_output = ', n_output
770 write(10, fmt=trim(fmt3))
'time_output =' , time_output0(n)
773 write(10, fmt=trim(fmt1))
' ' 775 write(10, fmt=trim(fmt1))
'Program version and date: '//version//
' / '//date
777 close(10, status=
'keep')
782 time_init = time_init0*year_sec
783 time_end = time_end0*year_sec
784 dtime = dtime0*year_sec
785 dtime_temp = dtime_temp0*year_sec
786 dtime_ser = dtime_ser0*year_sec
787 #if (OUTPUT==1 || OUTPUT==3) 788 dtime_out = dtime_out0*year_sec
790 #if (OUTPUT==2 || OUTPUT==3) 792 time_output(n) = time_output0(n)*year_sec
796 if (.not.approx_integer_multiple(dtime_temp, dtime,
eps_sp_dp)) &
797 stop
' >>> sico_init: dtime_temp must be a multiple of dtime!' 799 if (.not.approx_integer_multiple(dtime_ser, dtime,
eps_sp_dp)) &
800 stop
' >>> sico_init: dtime_ser must be a multiple of dtime!' 802 #if (OUTPUT==1 || OUTPUT==3) 803 if (.not.approx_integer_multiple(dtime_out, dtime,
eps_sp_dp)) &
804 stop
' >>> sico_init: dtime_out must be a multiple of dtime!' 811 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(
rho_w/
rho)
819 open(24, iostat=ios, file=trim(filename_with_path), recl=
rcl2, status=
'old')
821 if (ios /= 0) stop
' >>> sico_init: Error when opening the rock/sediment mask file!' 823 do n=1, 6;
read(24, fmt=
'(a)') ch_dummy;
end do 826 read(24, fmt=trim(fmt4)) (
maske_sedi(j,i), i=0,imax)
829 close(24, status=
'keep')
835 filename_with_path = trim(inpath)//
'/general/'//trim(grip_temp_file)
837 open(21, iostat=ios, file=trim(filename_with_path), status=
'old')
839 if (ios /= 0) stop
' >>> sico_init: Error when opening the data file for delta_ts!' 843 if (ch_dummy /=
'#')
then 844 write(6, fmt=*)
' >>> sico_init: grip_time_min, grip_time_stp, grip_time_max' 845 write(6, fmt=*)
' not defined in data file for delta_ts!' 854 read(21, fmt=*) d_dummy,
griptemp(n)
857 close(21, status=
'keep')
865 filename_with_path = trim(inpath)//
'/general/'//trim(sea_level_file)
867 open(21, iostat=ios, file=trim(filename_with_path), status=
'old')
869 if (ios /= 0) stop
' >>> sico_init: Error when opening the data file for z_sl!' 873 if (ch_dummy /=
'#')
then 874 write(6, fmt=*)
' >>> sico_init: specmap_time_min, specmap_time_stp, specmap_time_max' 875 write(6, fmt=*)
' not defined in data file for z_sl!' 887 close(21, status=
'keep')
905 stop
' >>> sico_init: Option Q_GEO_MOD==2 not available for this application!' 912 #if (REBOUND==1 || REBOUND==2) 914 #if (TIME_LAG_MOD==1) 918 #elif (TIME_LAG_MOD==2) 923 open(29, iostat=ios, file=trim(filename_with_path), recl=
rcl1, status=
'old')
925 if (ios /= 0) stop
' >>> sico_init: Error when opening the time-lag file!' 927 do n=1, 6;
read(29, fmt=
'(a)') ch_dummy;
end do 933 close(29, status=
'keep')
949 #if (FLEX_RIG_MOD==1) 953 #elif (FLEX_RIG_MOD==2) 958 open(29, iostat=ios, file=trim(filename_with_path), recl=
rcl1, status=
'old')
960 if (ios /= 0) stop
' >>> sico_init: Error when opening the flex-rig file!' 962 do n=1, 6;
read(29, fmt=
'(a)') ch_dummy;
end do 968 close(29, status=
'keep')
972 #elif (REBOUND==0 || REBOUND==1) 988 call boundary(time_init, dtime, dxi, deta, &
989 delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1005 #if (!defined(TEMP_INIT) || TEMP_INIT==1) 1007 #elif (TEMP_INIT==2) 1009 #elif (TEMP_INIT==3) 1011 #elif (TEMP_INIT==4) 1014 write(6, fmt=
'(a)')
' >>> sico_init:' 1015 write(6, fmt=
'(a)')
' TEMP_INIT must be set to either 1, 2, 3 or 4!' 1030 stop
' >>> sico_init: Parameter ENHMOD must be between 1 and 5!' 1041 call boundary(time_init, dtime, dxi, deta, &
1042 delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1066 stop
' >>> sico_init: Parameter ENHMOD must be between 1 and 5!' 1075 call boundary(time_init, dtime, dxi, deta, &
1076 delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1107 #if (MARGIN==1 || MARGIN==2) 1116 if ( (
maske(j,i)==0_i2b) &
1118 ( (
maske(j,i+1)==3_i2b) &
1119 .or.(
maske(j,i-1)==3_i2b) &
1120 .or.(
maske(j+1,i)==3_i2b) &
1121 .or.(
maske(j-1,i)==3_i2b) ) &
1125 if ( (
maske(j,i)==3_i2b) &
1127 ( (
maske(j,i+1)==0_i2b) &
1128 .or.(
maske(j,i-1)==0_i2b) &
1129 .or.(
maske(j+1,i)==0_i2b) &
1130 .or.(
maske(j-1,i)==0_i2b) ) &
1134 if ( (
maske(j,i)==3_i2b) &
1136 ( (
maske(j,i+1)==2_i2b) &
1137 .or.(
maske(j,i-1)==2_i2b) &
1138 .or.(
maske(j+1,i)==2_i2b) &
1139 .or.(
maske(j-1,i)==2_i2b) ) &
1143 if ( (
maske(j,i)==2_i2b) &
1145 ( (
maske(j,i+1)==3_i2b) &
1146 .or.(
maske(j,i-1)==3_i2b) &
1147 .or.(
maske(j+1,i)==3_i2b) &
1148 .or.(
maske(j-1,i)==3_i2b) ) &
1158 if ( (
maske(j,i)==2_i2b) &
1159 .and. (
maske(j+1,i)==3_i2b) &
1164 if ( (
maske(j,i)==2_i2b) &
1165 .and. (
maske(j-1,i)==3_i2b) &
1174 if ( (
maske(j,i)==2_i2b) &
1175 .and. (
maske(j,i+1)==3_i2b) &
1180 if ( (
maske(j,i)==2_i2b) &
1181 .and. (
maske(j,i-1)==3_i2b) &
1188 stop
' >>> sico_init: MARGIN must be either 1, 2 or 3!' 1193 #if (GRID==0 || GRID==1) 1197 dist_dxdy(jr,ir) = sqrt( (
real(ir,dp)*dxi)**2 + (
real(jr,dp)*deta)**2 )
1208 #if (DYNAMICS==1 || DYNAMICS==2) 1213 #if (MARGIN==3 || DYNAMICS==2) 1214 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
1229 stop
' >>> sico_init: DYNAMICS must be either 0, 1 or 2!' 1232 call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
1236 #if (CALCMOD==0 || CALCMOD==-1) 1261 if ( (
maske(j,i) == 0_i2b).and.(
n_cts(j,i) == 1_i2b) )
then 1274 #elif (CALCMOD==2 || CALCMOD==3) 1292 stop
' >>> sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!' 1300 filename_with_path = trim(outpath)//
'/'//trim(runname)//
'.ser' 1302 open(12, iostat=ios, file=trim(filename_with_path), status=
'new')
1304 if (ios /= 0) stop
' >>> sico_init: Error when opening the ser file!' 1309 1102
format(
' t(a) D_Ts(deg C) z_sl(m)',/, &
1310 ' V(m^3) V_g(m^3) V_f(m^3)', &
1311 ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1312 ' V_sle(m) V_t(m^3)', &
1314 ' H_max(m) H_t_max(m)', &
1315 ' zs_max(m) vs_max(m/a) Tbh_max(C)')
1316 1103
format(
'----------------------------------------------------', &
1317 '---------------------------------------')
1321 filename_with_path = trim(outpath)//
'/'//trim(runname)//
'.sed' 1323 open(15, iostat=ios, file=trim(filename_with_path), status=
'new')
1325 if (ios /= 0) stop
' >>> sico_init: Error when opening the sed file!' 1330 1108
format(
' t(a) D_Ts(deg C) z_sl(m)',/, &
1331 ' H_ave(m) Tbh_ave(C) Atb(m^2)')
1332 1109
format(
'----------------------------------------------------')
1344 x_core(1) = 3900.0_dp *1.0e+03_dp
1345 y_core(1) = 2000.0_dp *1.0e+03_dp
1350 x_core(2) = 3800.0_dp *1.0e+03_dp
1351 y_core(2) = 2000.0_dp *1.0e+03_dp
1356 x_core(3) = 3700.0_dp *1.0e+03_dp
1357 y_core(3) = 2000.0_dp *1.0e+03_dp
1362 x_core(4) = 3500.0_dp *1.0e+03_dp
1363 y_core(4) = 2000.0_dp *1.0e+03_dp
1368 x_core(5) = 3200.0_dp *1.0e+03_dp
1369 y_core(5) = 2000.0_dp *1.0e+03_dp
1374 x_core(6) = 2900.0_dp *1.0e+03_dp
1375 y_core(6) = 2000.0_dp *1.0e+03_dp
1380 x_core(7) = 2600.0_dp *1.0e+03_dp
1381 y_core(7) = 2000.0_dp *1.0e+03_dp
1383 filename_with_path = trim(outpath)//
'/'//trim(runname)//
'.core' 1385 open(14, iostat=ios, file=trim(filename_with_path), status=
'new')
1387 if (ios /= 0) stop
' >>> sico_init: Error when opening the core file!' 1394 1106
format(
' t(a) D_Ts(C) z_sl(m)',/, &
1395 ' H_P1(m) H_P2(m) H_P3(m) H_P4(m)', &
1396 ' H_P5(m) H_P6(m) H_P7(m)',/, &
1397 ' v_P1(m/a) v_P2(m/a) v_P3(m/a) v_P4(m/a)', &
1398 ' v_P5(m/a) v_P6(m/a) v_P7(m/a)',/, &
1399 ' T_P1(C) T_P2(C) T_P3(C) T_P4(C)', &
1400 ' T_P5(C) T_P6(C) T_P7(C)')
1401 1107
format(
'-----------------------------------------------------------------', &
1402 '---------------------------------------')
1408 #if (defined(OUTPUT_INIT)) 1410 #if (OUTPUT_INIT==0) 1411 flag_init_output = .false.
1412 #elif (OUTPUT_INIT==1) 1413 flag_init_output = .true.
1415 stop
' >>> sico_init: OUTPUT_INIT must be either 0 or 1!' 1419 flag_init_output = .true.
1425 flag_3d_output = .false.
1427 flag_3d_output = .true.
1429 stop
' >>> sico_init: ERGDAT must be either 0 or 1!' 1432 if (flag_init_output) &
1433 call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1434 flag_3d_output, ndat2d, ndat3d)
1438 if (time_output(1) <= time_init+
eps)
then 1441 flag_3d_output = .false.
1443 flag_3d_output = .true.
1445 stop
' >>> sico_init: ERGDAT must be either 0 or 1!' 1448 call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1449 flag_3d_output, ndat2d, ndat3d)
1455 flag_3d_output = .false.
1457 if (flag_init_output) &
1458 call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1459 flag_3d_output, ndat2d, ndat3d)
1461 if (time_output(1) <= time_init+
eps)
then 1463 flag_3d_output = .true.
1465 call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1466 flag_3d_output, ndat2d, ndat3d)
1471 stop
' >>> sico_init: OUTPUT must be either 1, 2 or 3!' 1474 if (flag_init_output)
then 1475 call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1476 call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1489 #if (GRID==0 || GRID==1) 1498 real(dp),
intent(out) :: dxi, deta
1500 integer(i4b) :: i, j, n
1501 integer(i4b) :: ios, n_dummy
1503 real(dp) :: xi0, eta0
1504 character :: ch_dummy
1506 real(dp),
parameter :: h_init = 1.0e-03_dp
1508 character(len= 8) :: ch_imax
1509 character(len=128) :: fmt4
1510 character(len=256) :: filename_with_path
1512 write(ch_imax, fmt=
'(i8)') imax
1513 write(fmt4, fmt=
'(a)')
'('//trim(adjustl(ch_imax))//
'(i1),i1)' 1519 filename_with_path = trim(inpath)//
'/'//trim(
ch_domain_short)//
'/'// &
1520 trim(mask_present_file)
1522 open(24, iostat=ios, file=trim(filename_with_path), recl=
rcl2, status=
'old')
1524 if (ios /= 0) stop
' >>> topography1: Error when opening the mask file!' 1526 do n=1, 6;
read(24, fmt=
'(a)') ch_dummy;
end do 1529 read(24, fmt=trim(fmt4)) (
maske(j,i), i=0,imax)
1532 close(24, status=
'keep')
1536 #if (GRID==0 || GRID==1) 1539 deta = dx *1000.0_dp
1542 eta0 = y0 *1000.0_dp
1546 stop
' >>> topography1: GRID==2 not allowed for this application!' 1557 if (
maske(j,i) <= 1)
then 1559 zs(j,i) =
zs(j,i) + h_init
1562 xi(i) = xi0 +
real(i,
dp)*dxi
1563 eta(j) = eta0 +
real(j,
dp)*deta
1590 #if (GRID==0 || GRID==1) /* Stereographic projection */ 1595 #elif (GRID==2) /* Geographic coordinates */ 1630 #if (GRID==0 || GRID==1) 1639 real(dp),
intent(out) :: dxi, deta
1641 integer(i4b) :: i, j, n
1643 real(dp) :: xi0, eta0
1644 character :: ch_dummy
1646 character(len= 8) :: ch_imax
1647 character(len=128) :: fmt4
1648 character(len=256) :: filename_with_path
1650 write(ch_imax, fmt=
'(i8)') imax
1651 write(fmt4, fmt=
'(a)')
'('//trim(adjustl(ch_imax))//
'(i1),i1)' 1657 filename_with_path = trim(inpath)//
'/'//trim(
ch_domain_short)//
'/'// &
1658 trim(mask_present_file)
1660 open(24, iostat=ios, file=trim(filename_with_path), recl=
rcl2, status=
'old')
1662 if (ios /= 0) stop
' >>> topography2: Error when opening the mask file!' 1664 do n=1, 6;
read(24, fmt=
'(a)') ch_dummy;
end do 1667 read(24, fmt=trim(fmt4)) (
maske(j,i), i=0,imax)
1670 close(24, status=
'keep')
1674 #if (GRID==0 || GRID==1) 1677 deta = dx *1000.0_dp
1680 eta0 = y0 *1000.0_dp
1684 stop
' >>> topography2: GRID==2 not allowed for this application!' 1695 xi(i) = xi0 +
real(i,
dp)*dxi
1696 eta(j) = eta0 +
real(j,
dp)*deta
1723 #if (GRID==0 || GRID==1) /* Stereographic projection */ 1728 #elif (GRID==2) /* Geographic coordinates */ 1761 subroutine topography3(dxi, deta, z_sl, anfdatname)
1765 #if (GRID==0 || GRID==1) 1774 character(len=100),
intent(in) :: anfdatname
1776 real(dp),
intent(out) :: dxi, deta, z_sl
1778 integer(i4b) :: i, j
1790 #if (GRID==0 || GRID==1) 1793 deta = dx *1000.0_dp
1797 stop
' >>> topography3: GRID==2 not allowed for this application!' 1807 #if (GRID==0 || GRID==1) /* Stereographic projection */ 1812 #elif (GRID==2) /* Geographic coordinates */ 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.
real(dp), dimension(:), allocatable phi_core
phi_core(n): Geographical latitude of the prescribed borehole positions
logical, dimension(0:jmax, 0:imax) flag_inner_point
flag_inner_point(j,i): Inner-point flag. .true.: inner point, .false.: margin point ...
integer(i4b) grip_time_stp
grip_time_stp: Time step of the data values for the surface temperature anomaly
integer(i2b), dimension(0:jmax, 0:imax) n_cts
n_cts(j,i): Mask for thermal conditions. -1: cold ice base, 0: temperate ice base with cold ice above...
integer(i4b) ndata_specmap
ndata_specmap: Number of data values for the sea level
subroutine, public calc_c_int_inv_table()
Computation of the inverse of the temperature integral of the specific heat of ice as a table (c_int_...
real(dp), dimension(0:kcmax) eaz_c
eaz_c(kc): Abbreviation for exp(aa*zeta(kc))
real(dp), dimension(0:jmax, 0:imax) q_geo
q_geo(j,i): Geothermal heat flux
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.
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.
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.
integer(i4b) specmap_time_max
specmap_time_max: Maximum time of the data values for the sea level
Material properties of ice: Rate factor, heat conductivity, specific heat (heat capacity), creep function, viscosity.
real(dp), parameter eps
eps: Small number
subroutine, public read_phys_para()
Reading of physical parameters.
integer(i2b) forcing_flag
forcing_flag: Flag for the forcing type. 1: forcing by a spatially constant surface temperature anoma...
logical, dimension(0:jmax, 0:imax) flag_grounding_line_1
flag_grounding_line_1(j,i): Grounding line flag. .true.: grounding line point (grounded ice point wit...
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) enth_t
enth_t(kt,j,i): Enthalpy in the lower (kt) ice domain
Conversion from temperature (temp) and water content (omega) to enthalpy (enth) and vice versa...
subroutine, public init_temp_water_age_2()
Initial temperature, water content and age (case ANF_DAT==2: ice-free conditions with relaxed bedrock...
real(dp), dimension(0:jmax, 0:imax) dh_t_dtau
dH_t_dtau(j,i): Derivative of H_t by tau (time)
real(dp), dimension(0:jmax, 0:imax) lambda
lambda(j,i): Geographic longitude of grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) dh_c_dtau
dH_c_dtau(j,i): Derivative of H_c by tau (time)
integer(i2b), dimension(0:jmax, 0:imax) maske
maske(j,i): Ice-land-ocean mask. 0: grounded ice, 1: ice-free land, 2: ocean, 3: floating ice ...
real(dp), dimension(0:jmax, 0:imax) p_b_w
p_b_w(j,i): Basal water pressure
integer(i4b) specmap_time_min
specmap_time_min: Minimum time of the data values for the sea level
real(dp), dimension(0:jmax, 0:imax) as_perp
as_perp(j,i): Accumulation-ablation function at the ice surface (SMB)
real(dp), dimension(0:jmax, 0:imax) area
area(j,i): Area of grid cell associated with grid point (i,j)
Computation of all components of the strain-rate tensor, the full effective strain rate and the shear...
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
integer(i2b), dimension(0:jmax, 0:imax) maske_sedi
maske_sedi(j,i): Sediment mask. 1: hard rock, 7: soft sediment, 2: ocean.
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...
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...
integer(i2b), dimension(0:jmax, 0:imax) kc_cts
kc_cts(j,i): Position kc of the CTS (for COLD, ENTC, ENTM)
integer(i4b) grip_time_max
grip_time_max: Maximum time of the data values for the surface temperature anomaly ...
real(dp), dimension(0:jmax, 0:imax) q_tld
Q_tld(j,i): Water drainage rate from the temperate layer.
subroutine, public calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy and the flux across ...
real(dp) b_max
b_max: Maximum accumulation rate
Declarations of global variables for SICOPOLIS (for the ANT domain).
integer(i4b), parameter rcl1
rcl1: Maximum length of record for input files (with factor 3 safety margin)
Initialisations for SICOPOLIS.
Computation of the flow enhancement factor.
integer(i4b) ndata_grip
ndata_grip: Number of data values for the surface temperature anomaly
Computation of the horizontal velocity vx, vy.
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), 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.
real(dp) b_min
b_min: Minimum accumulation rate
real(dp) aa
aa: Exponential stretch parameter of the non-equidistant vertical grid in the upper (kc) ice domain ...
Reading of several input data.
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), 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) year_zero
year_zero: SICOPOLIS year zero in astronomical year numbering [ = signed year CE (AD) ] ...
real(dp), dimension(:), allocatable specmap_zsl
specmap_zsl(n): Data values for the sea level
Computation of the melting and basal temperatures.
subroutine topograd_2(dxi, deta, n_switch)
Calculation of topography gradients on the staggered grid and on the grid points (the latter by fourt...
real(dp), dimension(-190:10) kappa
KAPPA(n): Tabulated values for the heat conductivity of ice.
Computation of the vertical velocity vz.
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
Definition of the components g11 and g22 of the metric tensor of the applied numerical coordinates...
Computation of the forward or inverse stereographic projection, alternatively for a spherical or an e...
real(dp), dimension(0:jmax, 0:imax) zl0
zl0(j,i): zl for isostatically relaxed ice-free conditions
subroutine, public calc_vz_static()
Computation of the vertical velocity vz for static ice.
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.
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) s_t
s_t: Gradient of surface temperature change with horizontal distance
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...
real(dp), dimension(-190:10) c
C(n): Tabulated values for the specific heat of ice.
real(dp), dimension(0:kcmax) eaz_c_quotient
eaz_c_quotient(kc): Abbreviation for (eaz_c(kc)-1.0)/(ea-1.0)
subroutine, public init_temp_water_age_1_3()
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==3: present-day initial topogr...
subroutine boundary(time, dtime, dxi, deta, delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
Main routine of boundary_m: Computation of the surface temperature (must be less than 0 deg C!) and o...
integer(i4b) specmap_time_stp
specmap_time_stp: Time step of the data values for the sea level
integer(i4b), parameter rcl2
rcl2: Maximum length of record for input mask files (with factor 3 safety margin) ...
subroutine, public output2(time, dxi, deta, delta_ts, glac_index, z_sl)
Writing of time-series data on file in ASCII format (and optionally in NetCDF format).
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...
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...
subroutine, public calc_enhance_3(time)
Computation of the flow enhancement factor. Case ENHMOD==3: two different values depending on time of...
integer(i4b) grip_time_min
grip_time_min: Minimum time of the data values for the surface temperature anomaly ...
real(dp), dimension(0:jmax, 0:imax) time_lag_asth
time_lag_asth(j,i): Time lag of the relaxing asthenosphere
real(dp), dimension(0:jmax) eta
eta(j): Coordinate eta (= y) of grid point j
real(dp), dimension(:), allocatable x_core
x_core(n): Coordinate xi (= x) of the prescribed borehole positions
real(dp) temp_min
temp_min: Minimum surface temperature
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.
subroutine, public calc_enhance_4()
Computation of the flow enhancement factor. Case ENHMOD==4: minimal anisotropic enhancement factor fo...
real(dp), dimension(0:ktmax) zeta_t
zeta_t(kt): Sigma coordinate zeta_t of grid point kt
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 ...
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) rad
rad: Radius of the model domain
real(dp) y_hat
y_hat: Coordinate eta (= y) of the centre of the model domain
real(dp) rho
RHO: Density of ice.
real(dp) x_hat
x_hat: Coordinate xi (= x) of the centre of the model domain
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...
Writing of output data on files.
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.
Declarations of global variables for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) dzs_dtau
dzs_dtau(j,i): Derivative of zs by tau (time)
real(dp), dimension(0:jmax, 0:imax) flex_rig_lith
flex_rig_lith(j,i): Flexural rigidity of the lithosphere
real(dp), dimension(:), allocatable griptemp
griptemp(n): Data values for the surface temperature anomaly
real(dp), dimension(0:jmax, 0:imax) dzl_dtau
dzl_dtau(j,i): Derivative of zl by tau (time)
real(dp), dimension(0:jmax, 0:imax) sq_g22_g
sq_g22_g(j,i): Square root of the coefficient g22 of the metric tensor on grid point (i...