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
96 integer(i2b),
dimension(0:JMAX,0:IMAX) :: maske_ref
97 real(dp) :: dtime0, dtime_temp0, dtime_wss0, dtime_out0, dtime_ser0
98 real(dp) :: time_init0, time_end0, time_output0(100)
100 character(len=100) :: anfdatname
101 character(len=256) :: filename_with_path
102 character(len=256) :: shell_command
103 character :: ch_dummy
104 logical :: flag_init_output, flag_3d_output
106 character(len=64),
parameter :: fmt1 =
'(a)', &
111 character(len= 8) :: ch_imax
112 character(len=128) :: fmt4
114 write(ch_imax, fmt=
'(i8)') imax
115 write(fmt4, fmt=
'(a)')
'('//trim(adjustl(ch_imax))//
'(i1),i1)' 117 write(unit=6, fmt=
'(a)')
' ' 118 write(unit=6, fmt=
'(a)')
' -------- sico_init --------' 119 write(unit=6, fmt=
'(a)')
' ' 131 #elif (defined(EMTP2SGE)) 139 #elif (defined(NHEM)) 143 #elif (defined(SCAND)) 147 #elif (defined(TIBET)) 151 #elif (defined(NMARS)) 155 #elif (defined(SMARS)) 167 stop
' >>> sico_init: No valid domain specified!' 188 #if (CALCTHK==3 || CALCTHK==6 || MARGIN==3 || DYNAMICS==2) 189 call lis_initialize(ierr)
205 if ( trim(version) /= trim(sico_version) ) &
206 stop
' >>> sico_init: SICOPOLIS version not compatible with header file!' 210 #if (!defined(DYNAMICS)) 211 stop
' >>> sico_init: DYNAMICS not defined in the header file!' 214 #if (!defined(CALCMOD)) 215 stop
' >>> sico_init: CALCMOD not defined in the header file!' 218 #if (defined(ENTHMOD)) 219 write(6, fmt=
'(a)')
' >>> sico_init: ENTHMOD must not be defined any more.' 220 write(6, fmt=
'(a)')
' Please update your header file!' 227 #if (GRID==0 || GRID==1) 229 stop
' >>> sico_init: GRID==0 or 1 not allowed for the Tibet application!' 234 if ( (approx_equal(dlambda, 1.0_dp/3.0_dp,
eps_sp_dp)) &
235 .and.(approx_equal(dphi , 1.0_dp/3.0_dp,
eps_sp_dp)) )
then 237 if ((imax /= 135).or.(jmax /= 51))
then 238 stop
' >>> sico_init: IMAX and/or JMAX wrong!' 241 else if ( (approx_equal(dlambda, 1.0_dp/6.0_dp,
eps_sp_dp)) &
242 .and.(approx_equal(dphi , 1.0_dp/6.0_dp,
eps_sp_dp)) )
then 244 if ((imax /= 270).or.(jmax /= 102))
then 245 stop
' >>> sico_init: IMAX and/or JMAX wrong!' 248 else if ( (approx_equal(dlambda, 1.0_dp/12.0_dp,
eps_sp_dp)) &
249 .and.(approx_equal(dphi , 1.0_dp/12.0_dp,
eps_sp_dp)) )
then 251 if ((imax /= 540).or.(jmax /= 204))
then 252 stop
' >>> sico_init: IMAX and/or JMAX wrong!' 256 stop
' >>> sico_init: DLAMBDA / DPHI wrong!' 265 #if (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1) 268 write(6, fmt=
'(a)')
' >>> sico_init: For options CALCMOD==0, 2, 3 or -1,' 269 write(6, fmt=
'(a)')
' the separate kt domain is redundant.' 270 write(6, fmt=
'(a)')
' Therefore, consider setting KTMAX to 2.' 271 write(6, fmt=
'(a)')
' ' 280 #if (TSURFACE == 5 && ACCSURFACE != 5) 281 stop
' >>> sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!' 284 #if (TSURFACE != 5 && ACCSURFACE == 5) 285 stop
' >>> sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!' 292 stop
' >>> sico_init: Option ADV_HOR==1 (central differences) not defined!' 300 #if ((MARGIN==3 && DYNAMICS==1) || DYNAMICS==2) /* requires SSA or SStA */ 302 write(6, fmt=
'(a)')
' >>> sico_init: Distortion correction not implemented' 303 write(6, fmt=
'(a)')
' for the shallow shelf approximation (SSA)' 304 write(6, fmt=
'(a)')
' or the shelfy stream approximation (SStA)' 305 write(6, fmt=
'(a)')
' -> GRID==0 required!' 316 #elif (TSURFACE == 5) 325 dtime_temp0 = dtime_temp0
327 dtime_wss0 = dtime_wss0
332 dzeta_c = 1.0_dp/
real(kcmax,dp)
333 dzeta_t = 1.0_dp/
real(ktmax,dp)
334 dzeta_r = 1.0_dp/
real(krmax,dp)
343 if (deform >=
eps)
then 423 if ((i <= imax).and.(j <= jmax))
then 435 anfdatname = anfdatname
437 #if (defined(YEAR_ZERO)) 443 time_init0 = time_init0
444 time_end0 = time_end0
445 dtime_ser0 = dtime_ser0
447 #if (OUTPUT==1 || OUTPUT==3) 448 dtime_out0 = dtime_out0
451 #if (OUTPUT==2 || OUTPUT==3) 453 time_output0( 1) = time_out0_01
454 time_output0( 2) = time_out0_02
455 time_output0( 3) = time_out0_03
456 time_output0( 4) = time_out0_04
457 time_output0( 5) = time_out0_05
458 time_output0( 6) = time_out0_06
459 time_output0( 7) = time_out0_07
460 time_output0( 8) = time_out0_08
461 time_output0( 9) = time_out0_09
462 time_output0(10) = time_out0_10
463 time_output0(11) = time_out0_11
464 time_output0(12) = time_out0_12
465 time_output0(13) = time_out0_13
466 time_output0(14) = time_out0_14
467 time_output0(15) = time_out0_15
468 time_output0(16) = time_out0_16
469 time_output0(17) = time_out0_17
470 time_output0(18) = time_out0_18
471 time_output0(19) = time_out0_19
472 time_output0(20) = time_out0_20
477 shell_command =
'if [ ! -d' 478 shell_command = trim(shell_command)//
' '//outpath
479 shell_command = trim(shell_command)//
' '//
'] ; then mkdir' 480 shell_command = trim(shell_command)//
' '//outpath
481 shell_command = trim(shell_command)//
' '//
'; fi' 482 call system(trim(shell_command))
485 filename_with_path = trim(outpath)//
'/'//trim(runname)//
'.log' 487 open(10, iostat=ios, file=trim(filename_with_path), status=
'new')
489 if (ios /= 0) stop
' >>> sico_init: Error when opening the log file!' 491 write(10, fmt=trim(fmt1))
'Computational domain:' 493 write(10, fmt=trim(fmt1))
' ' 495 write(10, fmt=trim(fmt2a))
'GRID = ', grid
496 write(10, fmt=trim(fmt1))
' ' 498 write(10, fmt=trim(fmt2))
'imax =', imax
499 write(10, fmt=trim(fmt2))
'jmax =', jmax
500 write(10, fmt=trim(fmt2))
'kcmax =', kcmax
501 write(10, fmt=trim(fmt2))
'ktmax =', ktmax
502 write(10, fmt=trim(fmt2))
'krmax =', krmax
503 write(10, fmt=trim(fmt1))
' ' 505 write(10, fmt=trim(fmt3))
'a =',
aa 506 write(10, fmt=trim(fmt1))
' ' 508 #if (GRID==0 || GRID==1) 509 stop
' >>> sico_init: GRID==0 or 1 not allowed for this application!' 511 write(10, fmt=trim(fmt3))
'lambda0 =', lambda_0
512 write(10, fmt=trim(fmt3))
'phi0 =', phi_0
513 write(10, fmt=trim(fmt3))
'dlambda =', dlambda
514 write(10, fmt=trim(fmt3))
'dphi =', dphi
516 write(10, fmt=trim(fmt1))
' ' 518 write(10, fmt=trim(fmt3))
'year_zero =',
year_zero 519 write(10, fmt=trim(fmt3))
'time_init =', time_init0
520 write(10, fmt=trim(fmt3))
'time_end =', time_end0
521 write(10, fmt=trim(fmt3))
'dtime =', dtime0
522 write(10, fmt=trim(fmt3))
'dtime_temp =', dtime_temp0
524 write(10, fmt=trim(fmt3))
'dtime_wss =', dtime_wss0
526 write(10, fmt=trim(fmt1))
' ' 528 write(10, fmt=trim(fmt2))
'ANF_DAT =', anf_dat
529 write(10, fmt=trim(fmt1))
'zs_present file = '//zs_present_file
531 #if (defined(ZB_PRESENT_FILE)) 532 write(10, fmt=trim(fmt1))
'zb_present file = '//zb_present_file
534 write(10, fmt=trim(fmt1))
'zl_present file = '//zl_present_file
536 write(10, fmt=trim(fmt1))
'zl0 file = '//zl0_file
537 write(10, fmt=trim(fmt1))
'mask_present file = '//mask_present_file
538 #if (ANF_DAT==1 && defined(TEMP_INIT)) 539 write(10, fmt=trim(fmt2))
'TEMP_INIT =', temp_init
542 write(10, fmt=trim(fmt1))
'Initial-value file = '//anfdatname
543 write(10, fmt=trim(fmt1))
'Path to initial-value file = '//anfdatpath
545 write(10, fmt=trim(fmt1))
' ' 547 write(10, fmt=trim(fmt1))
'Physical-parameter file = '//phys_para_file
548 write(10, fmt=trim(fmt1))
' ' 550 #if (CALCTHK==2 || CALCTHK==3 || CALCTHK==5 || CALCTHK==6) 551 write(10, fmt=trim(fmt3))
'ovi_weight =', ovi_weight
552 #if (CALCTHK==2 || CALCTHK==5) 553 write(10, fmt=trim(fmt3))
'omega_sor =', omega_sor
555 write(10, fmt=trim(fmt2a))
'iter_max_sor = ', iter_max_sor
558 write(10, fmt=trim(fmt1))
' ' 561 write(10, fmt=trim(fmt1))
'temp_mm_present file = '//temp_mm_present_file
563 write(10, fmt=trim(fmt3))
'delta_ts0 =', delta_ts0
565 write(10, fmt=trim(fmt3))
'sine_amplit =', sine_amplit
566 write(10, fmt=trim(fmt3))
'sine_period =', sine_period
568 write(10, fmt=trim(fmt1))
'GRIP file = '//grip_temp_file
569 write(10, fmt=trim(fmt3))
'grip_temp_fact =', grip_temp_fact
571 write(10, fmt=trim(fmt1))
'Glacial-index file = '//glac_ind_file
572 write(10, fmt=trim(fmt1))
'temp_mm_anom file = '//temp_mm_anom_file
573 write(10, fmt=trim(fmt3))
'temp_mm_anom fact = ', temp_mm_anom_fact
576 write(10, fmt=trim(fmt1))
'precip_mm_present file = '//precip_mm_present_file
578 write(10, fmt=trim(fmt3))
'accfact =', accfact
579 #elif (ACCSURFACE==2 || ACCSURFACE==3) 580 write(10, fmt=trim(fmt3))
'gamma_s =', gamma_s
581 #elif (ACCSURFACE==5) 582 write(10, fmt=trim(fmt1))
'precip_mm_anom file = '//precip_mm_anom_file
583 write(10, fmt=trim(fmt3))
'precip_mm_anom fact = ', precip_mm_anom_fact
585 #if (ACCSURFACE <= 3) 586 write(10, fmt=trim(fmt2))
'ELEV_DESERT =', elev_desert
587 #if (ELEV_DESERT == 1) 588 write(10, fmt=trim(fmt3))
'gamma_p =', gamma_p
589 write(10, fmt=trim(fmt3))
'zs_thresh =', zs_thresh
594 write(10, fmt=trim(fmt3))
'lambda_lti =', lambda_lti
595 write(10, fmt=trim(fmt3))
'temp_lti =', temp_lti
598 write(10, fmt=trim(fmt2a))
'SEA_LEVEL = ', sea_level
600 write(10, fmt=trim(fmt3))
'z_sl0 =', z_sl0
602 write(10, fmt=trim(fmt1))
'sea-level file = '//sea_level_file
604 write(10, fmt=trim(fmt1))
' ' 607 #if (MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3) 608 write(10, fmt=trim(fmt3))
'z_mar =', z_mar
609 write(10, fmt=trim(fmt1))
' ' 610 #elif (MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \ 611 || marine_ice_calving==6 || marine_ice_calving==7)
612 write(10, fmt=trim(fmt3))
'fact_z_mar =', fact_z_mar
613 write(10, fmt=trim(fmt1))
' ' 614 #elif (MARINE_ICE_FORMATION==2 && MARINE_ICE_CALVING==9) 615 write(10, fmt=trim(fmt3))
'calv_uw_coeff =', calv_uw_coeff
616 write(10, fmt=trim(fmt3))
'r1_calv_uw =', r1_calv_uw
617 write(10, fmt=trim(fmt3))
'r2_calv_uw =', r2_calv_uw
618 write(10, fmt=trim(fmt1))
' ' 621 #if (ICE_SHELF_CALVING==2) 622 write(10, fmt=trim(fmt3))
'H_calv =', h_calv
623 write(10, fmt=trim(fmt1))
' ' 627 #if (defined(BASAL_HYDROLOGY)) 628 write(10, fmt=trim(fmt2a))
'BASAL_HYDROLOGY = ', basal_hydrology
630 write(10, fmt=trim(fmt2a))
'SLIDE_LAW = ', slide_law
631 write(10, fmt=trim(fmt3))
'c_slide =',
c_slide 632 write(10, fmt=trim(fmt3))
'gamma_slide =', gamma_slide
633 write(10, fmt=trim(fmt2a))
'p_weert = ',
p_weert 634 write(10, fmt=trim(fmt2a))
'q_weert = ',
q_weert 635 #if (defined(TIME_RAMP_UP_SLIDE)) 636 write(10, fmt=trim(fmt3))
'time_ramp_up_slide =', time_ramp_up_slide
639 write(10, fmt=trim(fmt3))
'red_pres_limit_fact =', red_pres_limit_fact
641 #if (BASAL_HYDROLOGY==1) 642 write(10, fmt=trim(fmt3))
'c_Hw_slide =', c_hw_slide
643 write(10, fmt=trim(fmt3))
'Hw0_slide =', hw0_slide
645 write(10, fmt=trim(fmt1))
' ' 647 write(10, fmt=trim(fmt2a))
'Q_GEO_MOD = ', q_geo_mod
649 write(10, fmt=trim(fmt3))
'q_geo =',
q_geo 651 write(10, fmt=trim(fmt1))
'q_geo file = '//q_geo_file
653 write(10, fmt=trim(fmt1))
' ' 655 #if (defined(MARINE_ICE_BASAL_MELTING)) 656 write(10, fmt=trim(fmt2))
'MARINE_ICE_BASAL_MELTING =', marine_ice_basal_melting
657 #if (MARINE_ICE_BASAL_MELTING==2 || MARINE_ICE_BASAL_MELTING==3) 658 write(10, fmt=trim(fmt3))
'qbm_marine =', qbm_marine
660 write(10, fmt=trim(fmt1))
' ' 663 write(10, fmt=trim(fmt2))
'REBOUND =', rebound
665 write(10, fmt=trim(fmt3))
'frac_llra =', frac_llra
667 #if (REBOUND==1 || REBOUND==2) 668 write(10, fmt=trim(fmt2))
'TIME_LAG_MOD =', time_lag_mod
669 #if (TIME_LAG_MOD==1) 670 write(10, fmt=trim(fmt3))
'time_lag =', time_lag
671 #elif (TIME_LAG_MOD==2) 672 write(10, fmt=trim(fmt1))
'time_lag_file = '//time_lag_file
674 stop
' >>> sico_init: TIME_LAG_MOD must be either 1 or 2!' 678 write(10, fmt=trim(fmt2))
'FLEX_RIG_MOD =', flex_rig_mod
679 #if (FLEX_RIG_MOD==1) 680 write(10, fmt=trim(fmt3))
'flex_rig =', flex_rig
681 #elif (FLEX_RIG_MOD==2) 682 write(10, fmt=trim(fmt1))
'flex_rig_file = '//flex_rig_file
684 stop
' >>> sico_init: FLEX_RIG_MOD must be either 1 or 2!' 687 write(10, fmt=trim(fmt2))
'Q_LITHO =', q_litho
688 write(10, fmt=trim(fmt1))
' ' 691 write(10, fmt=trim(fmt3))
'gr_size =', gr_size
692 write(10, fmt=trim(fmt1))
' ' 695 write(10, fmt=trim(fmt3))
'sigma_res =', sigma_res
696 write(10, fmt=trim(fmt1))
' ' 699 write(10, fmt=trim(fmt2))
'ENHMOD =', enhmod
700 #if (ENHMOD==1 || ENHMOD==2 || ENHMOD==3) 701 write(10, fmt=trim(fmt3))
'enh_fact =', enh_fact
703 #if (ENHMOD==2 || ENHMOD==3) 704 write(10, fmt=trim(fmt3))
'enh_intg =', enh_intg
707 write(10, fmt=trim(fmt3))
'age_trans =', age_trans_0
710 write(10, fmt=trim(fmt3))
'date_trans1 =', date_trans1_0
711 write(10, fmt=trim(fmt3))
'date_trans2 =', date_trans2_0
712 write(10, fmt=trim(fmt3))
'date_trans3 =', date_trans3_0
714 #if (ENHMOD==4 || ENHMOD==5) 715 write(10, fmt=trim(fmt3))
'enh_compr =', enh_compr
716 write(10, fmt=trim(fmt3))
'enh_shear =', enh_shear
718 #if ((ENHMOD==1 || ENHMOD==2 || ENHMOD==3 || ENHMOD==4) && MARGIN==3) 719 write(10, fmt=trim(fmt3))
'enh_shelf =', enh_shelf
721 write(10, fmt=trim(fmt1))
' ' 723 write(10, fmt=trim(fmt3))
'numdiff_H_t =', numdiff_h_t
724 write(10, fmt=trim(fmt3))
'tau_cts =', tau_cts
725 write(10, fmt=trim(fmt3))
'vh_max =', vh_max
726 write(10, fmt=trim(fmt3))
'hd_min =', hd_min
727 write(10, fmt=trim(fmt3))
'hd_max =', hd_max
728 #if (defined(QBM_MIN)) 729 write(10, fmt=trim(fmt3))
'qbm_min =', qbm_min
730 #elif (defined(QB_MIN)) /* obsolete */ 731 write(10, fmt=trim(fmt3))
'qb_min =', qb_min
733 #if (defined(QBM_MAX)) 734 write(10, fmt=trim(fmt3))
'qbm_max =', qbm_max
735 #elif (defined(QB_MAX)) /* obsolete */ 736 write(10, fmt=trim(fmt3))
'qb_max =', qb_max
738 write(10, fmt=trim(fmt3))
'age_min =', age_min
739 write(10, fmt=trim(fmt3))
'age_max =', age_max
740 write(10, fmt=trim(fmt3))
'mean_accum =', mean_accum
742 write(10, fmt=trim(fmt3))
'age_diff =', agediff
744 write(10, fmt=trim(fmt1))
' ' 746 write(10, fmt=trim(fmt2a))
'DYNAMICS = ', dynamics
747 #if ((DYNAMICS==1 && MARGIN==3) || DYNAMICS==2) 748 #if (defined(LIS_OPTS)) 749 write(10, fmt=trim(fmt1))
'lis_opts = '//lis_opts
751 #if (defined(N_ITER_SSA)) 752 write(10, fmt=trim(fmt2a))
'n_iter_ssa = ', n_iter_ssa
754 #if (defined(RELAX_FACT_SSA)) 755 write(10, fmt=trim(fmt3))
'relax_fact_ssa =', relax_fact_ssa
758 #if (DYNAMICS==2 && defined(RATIO_SL_THRESH)) 759 write(10, fmt=trim(fmt3))
'ratio_sl_thresh =', ratio_sl_thresh
761 #if (DYNAMICS==2 && defined(SSTA_SIA_WEIGH_FCT)) 762 write(10, fmt=trim(fmt2a))
'SSTA_SIA_WEIGH_FCT = ', ssta_sia_weigh_fct
764 write(10, fmt=trim(fmt1))
' ' 766 write(10, fmt=trim(fmt2))
'CALCMOD =', calcmod
767 #if (CALCMOD==-1 && defined(TEMP_CONST)) 768 write(10, fmt=trim(fmt3))
'TEMP_CONST =', temp_const
770 #if (CALCMOD==-1 && defined(AGE_CONST)) 771 write(10, fmt=trim(fmt3))
'AGE_CONST =', age_const
773 #if (CALCMOD==1 && defined(CTS_MELTING_FREEZING)) 774 write(10, fmt=trim(fmt2))
'CTS_MELTING_FREEZING =', cts_melting_freezing
776 write(10, fmt=trim(fmt2))
'FLOW_LAW =', flow_law
777 write(10, fmt=trim(fmt2))
'FIN_VISC =', fin_visc
778 write(10, fmt=trim(fmt2))
'MARGIN =', margin
780 write(10, fmt=trim(fmt2))
'MARINE_ICE_FORMATION =', marine_ice_formation
781 write(10, fmt=trim(fmt2))
'MARINE_ICE_CALVING =', marine_ice_calving
783 write(10, fmt=trim(fmt2))
'ICE_SHELF_CALVING =', ice_shelf_calving
785 #if (defined(THK_EVOL)) 786 write(10, fmt=trim(fmt2))
'THK_EVOL =', thk_evol
788 stop
' >>> sico_init: Define THK_EVOL in header file!' 790 #if (defined(CALCTHK)) 791 write(10, fmt=trim(fmt2))
'CALCTHK =', calcthk
793 stop
' >>> sico_init: Define CALCTHK in header file!' 795 write(10, fmt=trim(fmt2))
'ADV_HOR =', adv_hor
796 write(10, fmt=trim(fmt2))
'ADV_VERT =', adv_vert
797 write(10, fmt=trim(fmt2))
'TOPOGRAD =', topograd
798 write(10, fmt=trim(fmt2))
'TSURFACE =', tsurface
799 write(10, fmt=trim(fmt2))
'ACCSURFACE =', accsurface
801 write(10, fmt=trim(fmt2))
'PRECIP_ANOM_INTERPOL =', precip_anom_interpol
803 write(10, fmt=trim(fmt2))
'SOLID_PRECIP =', solid_precip
804 write(10, fmt=trim(fmt2))
'ABLSURFACE =', ablsurface
805 #if (defined(MB_ACCOUNT)) 806 write(10, fmt=trim(fmt2a))
'MB_ACCOUNT = ', mb_account
808 write(10, fmt=trim(fmt1))
' ' 810 #if (defined(OUT_TIMES)) 811 write(10, fmt=trim(fmt2a))
'OUT_TIMES = ', out_times
813 #if (defined(OUTPUT_INIT)) 814 write(10, fmt=trim(fmt2a))
'OUTPUT_INIT = ', output_init
816 write(10, fmt=trim(fmt2a))
'OUTPUT = ', output
817 #if (OUTPUT==1 || OUTPUT==3) 818 write(10, fmt=trim(fmt3))
'dtime_out =' , dtime_out0
820 write(10, fmt=trim(fmt3))
'dtime_ser =' , dtime_ser0
821 #if (OUTPUT==1 || OUTPUT==2) 822 write(10, fmt=trim(fmt2a))
'ERGDAT = ', ergdat
824 write(10, fmt=trim(fmt2a))
'NETCDF = ', netcdf
825 #if (OUTPUT==2 || OUTPUT==3) 826 write(10, fmt=trim(fmt2a))
'n_output = ', n_output
828 write(10, fmt=trim(fmt3))
'time_output =' , time_output0(n)
831 write(10, fmt=trim(fmt1))
' ' 833 write(10, fmt=trim(fmt1))
'Program version and date: '//version//
' / '//date
835 close(10, status=
'keep')
840 time_init = time_init0*year_sec
841 time_end = time_end0*year_sec
842 dtime = dtime0*year_sec
843 dtime_temp = dtime_temp0*year_sec
845 dtime_wss = dtime_wss0*year_sec
847 dtime_ser = dtime_ser0*year_sec
848 #if (OUTPUT==1 || OUTPUT==3) 849 dtime_out = dtime_out0*year_sec
851 #if (OUTPUT==2 || OUTPUT==3) 853 time_output(n) = time_output0(n)*year_sec
857 if (.not.approx_integer_multiple(dtime_temp, dtime,
eps_sp_dp)) &
858 stop
' >>> sico_init: dtime_temp must be a multiple of dtime!' 861 if (.not.approx_integer_multiple(dtime_wss, dtime,
eps_sp_dp)) &
862 stop
' >>> sico_init: dtime_wss must be a multiple of dtime!' 865 if (.not.approx_integer_multiple(dtime_ser, dtime,
eps_sp_dp)) &
866 stop
' >>> sico_init: dtime_ser must be a multiple of dtime!' 868 #if (OUTPUT==1 || OUTPUT==3) 869 if (.not.approx_integer_multiple(dtime_out, dtime,
eps_sp_dp)) &
870 stop
' >>> sico_init: dtime_out must be a multiple of dtime!' 877 #if (GRID==0 || GRID==1) 879 stop
' >>> sico_init: GRID==0 or 1 not allowed for the Tibet application!' 884 trim(precip_mm_present_file)
886 open(21, iostat=ios, file=trim(filename_with_path), recl=
rcl1, status=
'old')
890 if (ios /= 0) stop
' >>> sico_init: Error when opening the precipitation file!' 892 do m=1, 6;
read(21, fmt=
'(a)') ch_dummy;
end do 895 do m=1, 3;
read(21, fmt=
'(a)') ch_dummy;
end do 901 close(21, status=
'keep')
915 trim(precip_mm_anom_file)
917 open(21, iostat=ios, file=trim(filename_with_path), recl=
rcl1, status=
'old')
921 if (ios /= 0) stop
' >>> sico_init: Error when opening the precip anomaly file!' 923 do m=1, 6;
read(21, fmt=
'(a)') ch_dummy;
end do 926 do m=1, 3;
read(21, fmt=
'(a)') ch_dummy;
end do 932 close(21, status=
'keep')
939 #if (PRECIP_ANOM_INTERPOL==1) 943 #elif (PRECIP_ANOM_INTERPOL==2) 948 stop
' >>> sico_init: Wrong value of switch PRECIP_ANOM_INTERPOL!' 958 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(
rho_w/
rho)
963 #if (GRID==0 || GRID==1) 965 stop
' >>> sico_init: GRID==0 or 1 not allowed for the Tibet application!' 970 trim(mask_present_file)
972 open(24, iostat=ios, file=trim(filename_with_path), recl=
rcl2, status=
'old')
976 if (ios /= 0) stop
' >>> sico_init: Error when opening the mask file!' 978 do m=1, 6;
read(24, fmt=
'(a)') ch_dummy;
end do 981 read(24, fmt=trim(fmt4)) (maske_ref(j,i), i=0,imax)
984 close(24, status=
'keep')
991 trim(temp_mm_present_file)
993 open(21, iostat=ios, file=trim(filename_with_path), recl=
rcl1, status=
'old')
997 if (ios /= 0) stop
' >>> sico_init: Error when opening the temperature file!' 999 do m=1, 6;
read(21, fmt=
'(a)') ch_dummy;
end do 1002 do m=1, 3;
read(21, fmt=
'(a)') ch_dummy;
end do 1008 close(21, status=
'keep')
1016 filename_with_path = trim(inpath)//
'/'//trim(
ch_domain_short)//
'/'// &
1017 trim(temp_mm_anom_file)
1019 open(21, iostat=ios, file=trim(filename_with_path), recl=
rcl1, status=
'old')
1023 if (ios /= 0) stop
' >>> sico_init: Error when opening the temperature anomaly file!' 1025 do m=1, 6;
read(21, fmt=
'(a)') ch_dummy;
end do 1028 do m=1, 3;
read(21, fmt=
'(a)') ch_dummy;
end do 1034 close(21, status=
'keep')
1043 #if (GRID==0 || GRID==1) 1045 stop
' >>> sico_init: GRID==0 or 1 not allowed for the Tibet application!' 1049 filename_with_path = trim(inpath)//
'/'//trim(
ch_domain_short)//
'/'// &
1050 trim(zs_present_file)
1052 open(21, iostat=ios, file=trim(filename_with_path), recl=
rcl1, status=
'old')
1056 if (ios /= 0) stop
' >>> sico_init: Error when opening the zs file!' 1058 do m=1, 6;
read(21, fmt=
'(a)') ch_dummy;
end do 1061 read(21, fmt=*) (
zs_ref(j,i), i=0,imax)
1064 close(21, status=
'keep')
1071 if (maske_ref(j,i) >= 2)
zs_ref(j,i) = 0.0_dp
1079 filename_with_path = trim(inpath)//
'/general/'//trim(grip_temp_file)
1081 open(21, iostat=ios, file=trim(filename_with_path), status=
'old')
1083 if (ios /= 0) stop
' >>> sico_init: Error when opening the data file for delta_ts!' 1087 if (ch_dummy /=
'#')
then 1088 write(6, fmt=*)
' >>> sico_init: grip_time_min, grip_time_stp, grip_time_max' 1089 write(6, fmt=*)
' not defined in data file for delta_ts!' 1098 read(21, fmt=*) d_dummy,
griptemp(n)
1101 close(21, status=
'keep')
1109 filename_with_path = trim(inpath)//
'/general/'//trim(glac_ind_file)
1111 open(21, iostat=ios, file=trim(filename_with_path), status=
'old')
1113 if (ios /= 0) stop
' >>> sico_init: Error when opening the glacial-index file!' 1117 if (ch_dummy /=
'#')
then 1118 write(6, fmt=*)
' >>> sico_init: gi_time_min, gi_time_stp, gi_time_max' 1119 write(6, fmt=*)
' not defined in glacial-index file!' 1131 close(21, status=
'keep')
1139 filename_with_path = trim(inpath)//
'/general/'//trim(sea_level_file)
1141 open(21, iostat=ios, file=trim(filename_with_path), status=
'old')
1143 if (ios /= 0) stop
' >>> sico_init: Error when opening the data file for z_sl!' 1147 if (ch_dummy /=
'#')
then 1148 write(6, fmt=*)
' >>> sico_init: specmap_time_min, specmap_time_stp, specmap_time_max' 1149 write(6, fmt=*)
' not defined in data file for z_sl!' 1161 close(21, status=
'keep')
1177 #elif (Q_GEO_MOD==2) 1181 filename_with_path = trim(inpath)//
'/'//trim(
ch_domain_short)//
'/'// &
1184 open(21, iostat=ios, file=trim(filename_with_path), recl=
rcl1, status=
'old')
1186 if (ios /= 0) stop
' >>> sico_init: Error when opening the qgeo file!' 1188 do m=1, 6;
read(21, fmt=
'(a)') ch_dummy;
end do 1191 read(21, fmt=*) (
q_geo(j,i), i=0,imax)
1194 close(21, status=
'keep')
1206 #if (REBOUND==0 || REBOUND==1) 1219 #if (REBOUND==1 || REBOUND==2) 1221 #if (TIME_LAG_MOD==1) 1225 #elif (TIME_LAG_MOD==2) 1227 filename_with_path = trim(inpath)//
'/'//trim(
ch_domain_short)//
'/'// &
1230 open(29, iostat=ios, file=trim(filename_with_path), recl=
rcl1, status=
'old')
1232 if (ios /= 0) stop
' >>> sico_init: Error when opening the time-lag file!' 1234 do n=1, 6;
read(29, fmt=
'(a)') ch_dummy;
end do 1240 close(29, status=
'keep')
1256 #if (FLEX_RIG_MOD==1) 1260 #elif (FLEX_RIG_MOD==2) 1262 filename_with_path = trim(inpath)//
'/'//trim(
ch_domain_short)//
'/'// &
1265 open(29, iostat=ios, file=trim(filename_with_path), recl=
rcl1, status=
'old')
1267 if (ios /= 0) stop
' >>> sico_init: Error when opening the flex-rig file!' 1269 do n=1, 6;
read(29, fmt=
'(a)') ch_dummy;
end do 1275 close(29, status=
'keep')
1279 #elif (REBOUND==0 || REBOUND==1) 1291 stop
' >>> sico_init: ANF_DAT==1 not allowed for Tibet application!' 1301 call boundary(time_init, dtime, dxi, deta, &
1302 delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1326 stop
' >>> sico_init: Parameter ENHMOD must be between 1 and 5!' 1335 call boundary(time_init, dtime, dxi, deta, &
1336 delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1367 #if (MARGIN==1 || MARGIN==2) 1376 if ( (
maske(j,i)==0_i2b) &
1378 ( (
maske(j,i+1)==3_i2b) &
1379 .or.(
maske(j,i-1)==3_i2b) &
1380 .or.(
maske(j+1,i)==3_i2b) &
1381 .or.(
maske(j-1,i)==3_i2b) ) &
1385 if ( (
maske(j,i)==3_i2b) &
1387 ( (
maske(j,i+1)==0_i2b) &
1388 .or.(
maske(j,i-1)==0_i2b) &
1389 .or.(
maske(j+1,i)==0_i2b) &
1390 .or.(
maske(j-1,i)==0_i2b) ) &
1394 if ( (
maske(j,i)==3_i2b) &
1396 ( (
maske(j,i+1)==2_i2b) &
1397 .or.(
maske(j,i-1)==2_i2b) &
1398 .or.(
maske(j+1,i)==2_i2b) &
1399 .or.(
maske(j-1,i)==2_i2b) ) &
1403 if ( (
maske(j,i)==2_i2b) &
1405 ( (
maske(j,i+1)==3_i2b) &
1406 .or.(
maske(j,i-1)==3_i2b) &
1407 .or.(
maske(j+1,i)==3_i2b) &
1408 .or.(
maske(j-1,i)==3_i2b) ) &
1418 if ( (
maske(j,i)==2_i2b) &
1419 .and. (
maske(j+1,i)==3_i2b) &
1424 if ( (
maske(j,i)==2_i2b) &
1425 .and. (
maske(j-1,i)==3_i2b) &
1434 if ( (
maske(j,i)==2_i2b) &
1435 .and. (
maske(j,i+1)==3_i2b) &
1440 if ( (
maske(j,i)==2_i2b) &
1441 .and. (
maske(j,i-1)==3_i2b) &
1448 stop
' >>> sico_init: MARGIN must be either 1, 2 or 3!' 1459 + (
sq_g22_g(jmax/2,imax/2)*
real(jr,dp)*deta)**2 )
1473 #if (DYNAMICS==1 || DYNAMICS==2) 1478 #if (MARGIN==3 || DYNAMICS==2) 1479 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
1494 stop
' >>> sico_init: DYNAMICS must be either 0, 1 or 2!' 1497 call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
1501 #if (CALCMOD==0 || CALCMOD==-1) 1526 if ( (
maske(j,i) == 0_i2b).and.(
n_cts(j,i) == 1_i2b) )
then 1539 #elif (CALCMOD==2 || CALCMOD==3) 1557 stop
' >>> sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!' 1565 filename_with_path = trim(outpath)//
'/'//trim(runname)//
'.ser' 1567 open(12, iostat=ios, file=trim(filename_with_path), status=
'new')
1569 if (ios /= 0) stop
' >>> sico_init: Error when opening the ser file!' 1576 1102
format(
' t(a) D_Ts(deg C) z_sl(m)',/, &
1577 ' V(m^3) V_g(m^3) V_f(m^3)', &
1578 ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1579 ' V_sle(m) V_t(m^3)', &
1581 ' H_max(m) H_t_max(m)', &
1582 ' zs_max(m) vs_max(m/a) Tbh_max(C)')
1583 1103
format(
'----------------------------------------------------', &
1584 '---------------------------------------')
1591 1112
format(
' t(a) glac_ind(1) z_sl(m)',/, &
1592 ' V(m^3) V_g(m^3) V_f(m^3)', &
1593 ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1594 ' V_sle(m) V_t(m^3)', &
1596 ' H_max(m) H_t_max(m)', &
1597 ' zs_max(m) vs_max(m/a) Tbh_max(C)')
1598 1113
format(
'----------------------------------------------------', &
1599 '---------------------------------------')
1607 filename_with_path = trim(outpath)//
'/'//trim(runname)//
'.core' 1609 open(14, iostat=ios, file=trim(filename_with_path), status=
'new')
1611 write(14,
'(1x,a)')
'---------------------' 1612 write(14,
'(1x,a)')
'No boreholes defined.' 1613 write(14,
'(1x,a)')
'---------------------' 1617 #if (defined(OUTPUT_INIT)) 1619 #if (OUTPUT_INIT==0) 1620 flag_init_output = .false.
1621 #elif (OUTPUT_INIT==1) 1622 flag_init_output = .true.
1624 stop
' >>> sico_init: OUTPUT_INIT must be either 0 or 1!' 1628 flag_init_output = .true.
1634 flag_3d_output = .false.
1636 flag_3d_output = .true.
1638 stop
' >>> sico_init: ERGDAT must be either 0 or 1!' 1641 if (flag_init_output) &
1642 call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1643 flag_3d_output, ndat2d, ndat3d)
1647 if (time_output(1) <= time_init+
eps)
then 1650 flag_3d_output = .false.
1652 flag_3d_output = .true.
1654 stop
' >>> sico_init: ERGDAT must be either 0 or 1!' 1657 call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1658 flag_3d_output, ndat2d, ndat3d)
1664 flag_3d_output = .false.
1666 if (flag_init_output) &
1667 call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1668 flag_3d_output, ndat2d, ndat3d)
1670 if (time_output(1) <= time_init+
eps)
then 1672 flag_3d_output = .true.
1674 call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1675 flag_3d_output, ndat2d, ndat3d)
1680 stop
' >>> sico_init: OUTPUT must be either 1, 2 or 3!' 1683 if (flag_init_output)
then 1684 call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1685 call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1700 real(dp),
intent(out) :: dxi, deta
1706 dxi=0.0_dp; deta=0.0_dp
1708 write(6,
'(/1x,a)')
'>>> topography1: ANF_DAT==1 not defined for this domain!' 1720 #if (GRID==0 || GRID==1) 1729 real(dp),
intent(out) :: dxi, deta
1731 integer(i4b) :: i, j, n
1733 real(dp) :: xi0, eta0
1734 character :: ch_dummy
1736 character(len= 8) :: ch_imax
1737 character(len=128) :: fmt4
1738 character(len=256) :: filename_with_path
1740 write(ch_imax, fmt=
'(i8)') imax
1741 write(fmt4, fmt=
'(a)')
'('//trim(adjustl(ch_imax))//
'(i1),i1)' 1745 filename_with_path = trim(inpath)//
'/'//trim(
ch_domain_short)//
'/'// &
1748 open(23, iostat=ios, file=trim(filename_with_path), recl=
rcl1, status=
'old')
1750 if (ios /= 0) stop
' >>> topography2: Error when opening the zl0 file!' 1752 filename_with_path = trim(inpath)//
'/'//trim(
ch_domain_short)//
'/'// &
1753 trim(mask_present_file)
1755 open(24, iostat=ios, file=trim(filename_with_path), recl=
rcl2, status=
'old')
1757 if (ios /= 0) stop
' >>> topography2: Error when opening the mask file!' 1759 do n=1, 6;
read(23, fmt=
'(a)') ch_dummy;
end do 1762 read(23, fmt=*) (
zl0(j,i), i=0,imax)
1765 do n=1, 6;
read(24, fmt=
'(a)') ch_dummy;
end do 1768 read(24, fmt=trim(fmt4)) (
maske(j,i), i=0,imax)
1771 close(23, status=
'keep')
1772 close(24, status=
'keep')
1785 if (
maske(j,i) <= 1)
then 1792 #if (MARGIN==1 || MARGIN==2) 1802 xi(i) = xi0 +
real(i,
dp)*dxi
1803 eta(j) = eta0 +
real(j,
dp)*deta
1830 #if (GRID==0 || GRID==1) /* Stereographic projection */ 1835 #elif (GRID==2) /* Geographic coordinates */ 1868 subroutine topography3(dxi, deta, z_sl, anfdatname)
1872 #if (GRID==0 || GRID==1) 1881 character(len=100),
intent(in) :: anfdatname
1883 real(dp),
intent(out) :: dxi, deta, z_sl
1885 integer(i4b) :: i, j, n
1887 character(len=256) :: filename_with_path
1888 character :: ch_dummy
1896 filename_with_path = trim(inpath)//
'/'//trim(
ch_domain_short)//
'/'// &
1899 open(23, iostat=ios, file=trim(filename_with_path), recl=
rcl1, status=
'old')
1901 if (ios /= 0) stop
' >>> topography3: Error when opening the zl0 file!' 1903 do n=1, 6;
read(23, fmt=
'(a)') ch_dummy;
end do 1906 read(23, fmt=*) (
zl0(j,i), i=0,imax)
1909 close(23, status=
'keep')
1922 #if (GRID==0 || GRID==1) /* Stereographic projection */ 1927 #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.
subroutine, public read_kei()
Reading of the tabulated kei function.
logical, dimension(0:jmax, 0:imax) flag_inner_point
flag_inner_point(j,i): Inner-point flag. .true.: inner point, .false.: margin point ...
integer(i4b) grip_time_stp
grip_time_stp: Time step of the data values for the surface temperature anomaly
integer(i2b), dimension(0:jmax, 0:imax) n_cts
n_cts(j,i): Mask for thermal conditions. -1: cold ice base, 0: temperate ice base with cold ice above...
integer(i4b) ndata_specmap
ndata_specmap: Number of data values for the sea level
subroutine, public calc_c_int_inv_table()
Computation of the inverse of the temperature integral of the specific heat of ice as a table (c_int_...
real(dp), dimension(0:kcmax) eaz_c
eaz_c(kc): Abbreviation for exp(aa*zeta(kc))
real(dp), dimension(0:jmax, 0:imax) q_geo
q_geo(j,i): Geothermal heat flux
real(dp), dimension(0:jmax, 0:imax, 12) temp_mm_present
temp_mm_present(j,i,n): Present-day mean monthly surface temperature
subroutine sico_init(delta_ts, glac_index, mean_accum, dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, time, time_init, time_end, time_output, dxi, deta, dzeta_c, dzeta_t, dzeta_r, z_sl, dzsl_dtau, z_mar, ii, jj, nn, ndat2d, ndat3d, n_output, runname)
Main routine of sico_init_m: Initialisations for SICOPOLIS.
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)
real(dp), dimension(:), allocatable glacial_index
glacial_index(n): Data values for the glacial index
Computation of all components of the strain-rate tensor, the full effective strain rate and the shear...
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...
real(dp), dimension(0:jmax, 0:imax) zs_ref
zs_ref(j,i): Reference elevation for precip_present, temp_ma_present and temp_mj_present ...
subroutine, public calc_enhance_5()
Computation of the flow enhancement factor. Case ENHMOD==5: minimal anisotropic enhancement factor fo...
real(dp) b
B: Semi-minor axis of the planet.
real(dp), dimension(0:jmax, 0:imax) zl
zl(j,i): Coordinate z of the lithosphere surface
subroutine, public calc_vxy_static()
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy etc. for static ice.
subroutine, public output4(time, dxi, deta, delta_ts, glac_index, z_sl)
Writing of time-series data of the deep ice cores on file in ASCII format (and optionally in NetCDF f...
integer(i2b), dimension(0:jmax, 0:imax) kc_cts
kc_cts(j,i): Position kc of the CTS (for COLD, ENTC, ENTM)
integer(i4b) grip_time_max
grip_time_max: Maximum time of the data values for the surface temperature anomaly ...
real(dp), dimension(0:jmax, 0:imax) q_tld
Q_tld(j,i): Water drainage rate from the temperate layer.
subroutine, public calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy and the flux across ...
Declarations of global variables for SICOPOLIS (for the ANT domain).
integer(i4b), parameter rcl1
rcl1: Maximum length of record for input files (with factor 3 safety margin)
Initialisations for SICOPOLIS.
integer(i4b) ndata_gi
ndata_gi: Number of data values for the glacial index
Computation of the flow enhancement factor.
integer(i4b) ndata_grip
ndata_grip: Number of data values for the surface temperature anomaly
Computation of the horizontal velocity vx, vy.
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.
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
integer(i4b) gi_time_min
gi_time_min: Minimum time of the data values for the glacial index
real(dp), dimension(0:jmax, 0:imax) h_c
H_c(j,i): Thickness of ice in the upper (kc) domain (thickness of the cold-ice layer for POLY...
real(dp) kei_r_max
kei_r_max: Maximum value of the argument r of the tabulated kei function
real(dp) year_zero
year_zero: SICOPOLIS year zero in astronomical year numbering [ = signed year CE (AD) ] ...
real(dp), dimension(:), allocatable specmap_zsl
specmap_zsl(n): Data values for the sea level
Computation of the melting and basal temperatures.
subroutine topograd_2(dxi, deta, n_switch)
Calculation of topography gradients on the staggered grid and on the grid points (the latter by fourt...
real(dp), dimension(-190:10) kappa
KAPPA(n): Tabulated values for the heat conductivity of ice.
Computation of the vertical velocity vz.
integer(i4b) gi_time_stp
gi_time_stp: Time step of the data values for the glacial index
real(dp), dimension(0:jmax, 0:imax) dzm_dtau
dzm_dtau(j,i): Derivative of zm by tau (time)
real(dp), dimension(0:krmax) zeta_r
zeta_r(kr): Sigma coordinate zeta_r of grid point kr
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) omega_c
omega_c(kc,j,i): Water content in the upper (kc) ice domain
integer(i2b), dimension(0:jmax, 0:imax) maske_old
maske_old(j,i): Old value of maske (at the previous time step)
character(len=64) ch_domain_long
ch_domain_long: Long name of the computational domain
Definition of the components g11 and g22 of the metric tensor of the applied numerical coordinates...
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.
integer(i4b) gi_time_max
gi_time_max: Maximum time of the data values for the glacial index
real(dp) phi0
PHI0: Standard parallel of the stereographic projection.
subroutine, public ice_mat_eqs_pars(RF_table, R_T_val, KAPPA_table, C_table, n_tmp_min, n_tmp_max, RHO_I_val, RHO_C_val, KAPPA_C_val, C_C_val)
Setting of required physical parameters.
real(dp), dimension(0:jmax, 0:imax, 12) gamma_precip_lgm_anom
gamma_precip_lgm_anom(j,i,n): negative natural logarithm of precip_lgm_anom(j,i,n) ...
Comparison of single- or double-precision floating-point numbers.
subroutine, public calc_enhance_1()
Computation of the flow enhancement factor. Case ENHMOD==1: constant for grounded ice...
real(dp), dimension(0:jmax, 0:imax) zm
zm(j,i): Coordinate z of the bottom of the upper (kc) ice domain = top of the lower (kt) ice domain (...
real(dp), dimension(-190:10) rf
RF(n): Tabulated values for the rate factor of cold ice.
real(dp), dimension(0:jmax, 0:imax) q_b_tot
Q_b_tot(j,i): Sum of Q_bm and Q_tld.
real(dp), dimension(0:jmax, 0:imax) q_bm
Q_bm(j,i): Basal melting rate.
real(dp) function, public enth_fct_temp_omega(temp_val, omega_val)
Enthalpy as a function of temperature and water content.
subroutine, public calc_vxy_sia(dzeta_c, dzeta_t)
Computation of the shear stresses txz, tyz, the effective shear stress sigma, the depth-averaged flui...
real(dp), dimension(-190:10) c
C(n): Tabulated values for the specific heat of ice.
real(dp), dimension(0:kcmax) eaz_c_quotient
eaz_c_quotient(kc): Abbreviation for (eaz_c(kc)-1.0)/(ea-1.0)
real(dp), dimension(-10000:10000) kei
kei(n): Tabulated values of the kei function (Kelvin function of zero order)
subroutine boundary(time, dtime, dxi, deta, delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
Main routine of boundary_m: Computation of the surface temperature (must be less than 0 deg C!) and o...
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...
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(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) 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 ...
real(dp) l
L: Latent heat of ice.
real(dp), parameter eps_sp_dp
eps_sp_dp: Small number to single-precision accuracy in double precision
integer(i4b) n_data_kei
n_data_kei: Number of tabulated values of the kei function
real(dp) rho
RHO: Density of ice.
subroutine topography3(dxi, deta, z_sl, anfdatname)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
integer(i4b), dimension(0:jmax, 0:imax) p_weert
p_weert(j,i): Weertman exponent for the basal shear stress
integer(i4b) n_core
n_core: Number of positions to be considered in the time-series file for deep boreholes ...
integer, parameter dp
Double-precision reals.
subroutine topography2(dxi, deta)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
real(dp), dimension(0:jmax, 0:imax) zb
zb(j,i): Coordinate z of the ice base
logical, dimension(0:jmax, 0:imax) flag_grounding_line_2
flag_grounding_line_2(j,i): Grounding line flag. .true.: grounding line point (floating ice point wit...
real(dp), dimension(0:jmax, 0:imax, 12) precip_present
precip_present(j,i,n): Present-day mean monthly precipitation rate at the ice surface ...
Writing of output data on files.
real(dp), dimension(0:jmax, 0:imax) dzb_dtau
dzb_dtau(j,i): Derivative of zb by tau (time)
real(dp), dimension(0:kcmax) zeta_c
zeta_c(kc): Sigma coordinate zeta_c of grid point kc
real(dp) lambda0
LAMBDA0: Reference longitude (central meridian) of the stereographic projection.
subroutine, public calc_c_int_table(c_table, n_tmp_min, n_tmp_max, L_val)
Computation of the temperature integral of the specific heat of ice as a table (c_int_table). Further, definition of the latent heat of ice.
real(dp), dimension(0:jmax, 0:imax) h_t
H_t(j,i): Thickness of ice in the lower (kt) domain (thickness of the temperate layer for POLY...
real(dp), dimension(0:jmax, 0:imax, 12) precip_lgm_anom
precip_lgm_anom(j,i,n): LGM anomaly (ratio LGM/present) of the mean monthly precipitation rate at the...
subroutine topography1(dxi, deta)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
real(dp), dimension(0:jmax, 0:imax, 12) temp_mm_lgm_anom
temp_mm_lgm_anom(j,i,n): LGM anomaly (difference LGM - present) of the mean monthly surface temperatu...
real(dp), dimension(0:jmax, 0:imax) h_w
H_w(j,i): Thickness of the water column under the ice base.
subroutine, public calc_vxy_b_sia(time, z_sl)
Computation of the basal horizontal velocity vx_b, vy_b in the shallow ice approximation.
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...