54 subroutine boundary(time, dtime, dxi, deta, &
55 delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
63 && (marine_ice_formation==2) \
64 && (marine_ice_calving==9))
77 real(dp),
intent(in) :: time, dtime, dxi, deta
79 real(dp),
intent(out) :: delta_ts, glac_index, dzsl_dtau, z_mar
80 real(dp),
intent(inout) :: z_sl
88 integer(i4b) :: i, j, n
89 integer(i4b) :: i_gr, i_kl
90 integer(i4b) :: nrec_temp_precip
91 integer(i4b),
save :: nrec_temp_precip_save = -1
93 real(dp) :: z_sl_min, t1, t2, t3, t4, t5, t6
94 real(dp) :: year_sec_inv, time_in_years
95 real(dp) :: time_gr, time_kl
96 real(dp) :: z_sle_present, z_sle_help
98 real(dp),
dimension(0:JMAX,0:IMAX,0:12) :: precip
99 real(dp),
dimension(0:JMAX,0:IMAX) :: &
101 real(dp),
dimension(0:JMAX,0:IMAX,12) :: temp_mm
102 real(dp),
dimension(0:JMAX,0:IMAX) :: temp_ma, temp_ampl
103 real(dp),
dimension(0:JMAX,0:IMAX) :: temp_ma_anom, temp_mj_anom, precip_ma_anom
104 real(dp),
dimension(0:IMAX,0:JMAX),
save :: temp_ma_anom_tra, temp_mj_anom_tra, &
109 real(dp),
dimension(12) :: temp_mm_help
110 real(dp) :: temp_jja_help
111 real(dp),
dimension(0:JMAX,0:IMAX) :: et
112 real(dp) :: theta_ma, c_ma, kappa_ma, gamma_ma, &
113 theta_mj, c_mj, kappa_mj, gamma_mj
114 real(dp) :: sine_factor
115 real(dp) :: gamma_p, zs_thresh, &
116 temp_rain, temp_snow, &
117 inv_delta_temp_rain_snow, coeff(0:5), inv_sqrt2_s_stat, &
118 precip_fact, frac_solid
119 real(dp) :: s_stat, &
120 phi_sep, temp_lt, temp_ht, inv_delta_temp_ht_lt, &
121 beta1_lt, beta1_ht, beta2_lt, beta2_ht, &
122 beta1, beta2, pmax, mu, lambda_lti, temp_lti
123 logical,
dimension(0:JMAX,0:IMAX) :: check_point
124 logical,
save :: firstcall = .true.
129 integer(i4b) :: nc3cor(3)
131 integer(i4b) :: nc3cnt(3)
135 real(dp),
parameter :: &
136 inv_twelve = 1.0_dp/12.0_dp, one_third = 1.0_dp/3.0_dp
138 character(len=64),
parameter :: thisroutine =
'boundary' 140 year_sec_inv = 1.0_dp/year_sec
141 time_in_years = time*year_sec_inv
160 delta_ts = sine_amplit &
161 *cos(2.0_dp*
pi*time_in_years/sine_period) &
172 i_kl = floor((time_in_years &
176 i_gr = ceiling((time_in_years &
180 if (i_kl == i_gr)
then 191 *(time-time_kl)/(time_gr-time_kl)
200 delta_ts = delta_ts * grip_temp_fact
211 i_kl = floor((time_in_years &
215 i_gr = ceiling((time_in_years &
219 if (i_kl == i_gr)
then 230 *(time-time_kl)/(time_gr-time_kl)
242 #elif (TSURFACE==6 && ACCSURFACE==6) 253 if ( nrec_temp_precip < 0 )
then 254 stop
' >>> boundary: nrec_temp_precip < 0, not allowed!' 256 stop
' >>> boundary: nrec_temp_precip > ndata_temp_precip, not allowed!' 259 if ( nrec_temp_precip /= nrec_temp_precip_save )
then 263 nc3cor(3) = nrec_temp_precip + 1
271 start=nc3cor, count=nc3cnt), thisroutine )
276 start=nc3cor, count=nc3cnt), thisroutine )
281 start=nc3cor, count=nc3cnt), thisroutine )
285 temp_ma_anom = transpose(temp_ma_anom_tra)
286 temp_mj_anom = transpose(temp_mj_anom_tra)
287 precip_ma_anom = transpose(precip_ma_anom_tra) &
288 *(1.0e-03_dp*year_sec_inv)*(
rho_w/
rho)
291 nrec_temp_precip_save = nrec_temp_precip
307 if (time_in_years <
real(glann_time_min,
dp)) then
309 else if (time_in_years <
real(glann_time_max,
dp)) then
311 i_kl = floor((time_in_years &
312 -
real(glann_time_min,
dp))/
real(glann_time_stp,
dp))
315 i_gr = ceiling((time_in_years &
316 -
real(glann_time_min,
dp))/
real(glann_time_stp,
dp))
317 i_gr = min(i_gr, ndata_glann)
319 if (i_kl == i_gr)
then 325 time_kl = (glann_time_min + i_kl*glann_time_stp) *year_sec
326 time_gr = (glann_time_min + i_gr*glann_time_stp) *year_sec
329 +(dt_glann_climber(i_gr)-dt_glann_climber(i_kl)) &
330 *(time-time_kl)/(time_gr-time_kl)
336 dt_glann = dt_glann_climber(ndata_glann)
352 t1 = -250000.0_dp *year_sec
353 t2 = -140000.0_dp *year_sec
354 t3 = -125000.0_dp *year_sec
355 t4 = -21000.0_dp *year_sec
356 t5 = -8000.0_dp *year_sec
357 t6 = 0.0_dp *year_sec
361 else if (time < t2)
then 362 z_sl = z_sl_min*(time-t1)/(t2-t1)
363 else if (time < t3)
then 364 z_sl = -z_sl_min*(time-t3)/(t3-t2)
365 else if (time < t4)
then 366 z_sl = z_sl_min*(time-t3)/(t4-t3)
367 else if (time < t5)
then 368 z_sl = -z_sl_min*(time-t5)/(t5-t4)
369 else if (time < t6)
then 383 i_kl = floor((time_in_years &
387 i_gr = ceiling((time_in_years &
391 if (i_kl == i_gr)
then 402 *(time-time_kl)/(time_gr-time_kl)
415 if ( z_sl_old > -999999.9_dp )
then 416 dzsl_dtau = (z_sl-z_sl_old)/dtime
425 #if (MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3) 427 #elif (MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5) 428 z_mar = fact_z_mar*z_sl
429 #elif (MARINE_ICE_CALVING==6 || MARINE_ICE_CALVING==7) 430 if (z_sl >= -80.0_dp)
then 433 z_mar = 10.25_dp*(z_sl+80.0_dp)-200.0_dp
435 z_mar = fact_z_mar*z_mar
447 check_point(j,i) = .false.
453 if (
maske(j,i) >= 2)
then 454 check_point(j ,i ) = .true.
455 check_point(j ,i+1) = .true.
456 check_point(j ,i-1) = .true.
457 check_point(j+1,i ) = .true.
458 check_point(j-1,i ) = .true.
465 if (check_point(j,i))
then 475 if (check_point(j,i))
then 483 #if (TEMP_PRESENT_PARA==1) /* Parameterization by Ritz et al. (1997) */ 486 gamma_ma = -7.992e-03_dp
491 gamma_mj = -6.277e-03_dp
495 #elif (TEMP_PRESENT_PARA==2) /* Parameterization by Fausto et al. (2009) */ 498 gamma_ma = -6.309e-03_dp
500 kappa_ma = -0.0672_dp
503 gamma_mj = -5.426e-03_dp
505 kappa_mj = -0.0518_dp
509 stop
' >>> boundary: Parameter TEMP_PRESENT_PARA must be either 1 or 2!' 551 temp_ma(j,i) =
temp_ma_present(j,i) + temp_ma_anom_fact*temp_ma_anom(j,i)
552 temp_mm(j,i,7) =
temp_mj_present(j,i) + temp_mj_anom_fact*temp_mj_anom(j,i)
558 temp_ampl(j,i) = temp_mm(j,i,7) - temp_ma(j,i)
560 if (temp_ampl(j,i) <
eps)
then 571 sine_factor = sin((
real(n,
dp)-4.0_dp)*
pi/6.0_dp)
575 temp_mm(j,i,n) = temp_ma(j,i) + sine_factor*temp_ampl(j,i)
583 #if (INITMIP_CONST_SMB==1) /* SMB held constant at its initial distribution */ 589 gamma_p = gamma_p*1.0e-03_dp
591 zs_thresh = zs_thresh
595 #if (SOLID_PRECIP==1) /* Marsiat (1994) */ 602 inv_delta_temp_rain_snow = 1.0_dp/(temp_rain-temp_snow)
604 #elif (SOLID_PRECIP==2) /* Bales et al. (2009) */ 611 coeff(0) = 5.4714e-01_dp
612 coeff(1) = -9.1603e-02_dp
613 coeff(2) = -3.314e-03_dp
614 coeff(3) = 4.66e-04_dp
615 coeff(4) = 3.8e-05_dp
616 coeff(5) = 6.0e-07_dp
618 #elif (SOLID_PRECIP==3) /* Huybrechts and de Wolde (1999) */ 622 temp_snow = temp_rain
628 inv_sqrt2_s_stat = 1.0_dp/(sqrt(2.0_dp)*s_stat)
632 #if (ABLSURFACE==1 || ABLSURFACE==2) 641 inv_delta_temp_ht_lt = 1.0_dp/(temp_ht-temp_lt)
655 #elif (ABLSURFACE==3) 657 lambda_lti = lambda_lti *(0.001_dp*year_sec_inv)*(
rho_w/
rho)
667 #if (ABLSURFACE==1 || ABLSURFACE==2) 669 if (
phi(j,i) <= phi_sep)
then 676 if (temp_mm(j,i,7) >= temp_ht)
then 679 else if (temp_mm(j,i,7) <= temp_lt)
then 684 + (beta1_ht-beta1_lt) &
685 *inv_delta_temp_ht_lt*(temp_mm(j,i,7)-temp_lt)
687 + (beta2_lt-beta2_ht) &
688 *(inv_delta_temp_ht_lt*(temp_ht-temp_mm(j,i,7)))**3
705 #elif (ELEV_DESERT==1) 707 if (
zs_ref(j,i) < zs_thresh)
then 709 = exp(gamma_p*(max(
zs(j,i),zs_thresh)-zs_thresh))
712 = exp(gamma_p*(max(
zs(j,i),zs_thresh)-
zs_ref(j,i)))
716 stop
' >>> boundary: Parameter ELEV_DESERT must be either 0 or 1!' 728 precip_fact = accfact
729 #elif (ACCSURFACE==2) 730 precip_fact = 1.0_dp + gamma_s*delta_ts
731 #elif (ACCSURFACE==3) 732 precip_fact = exp(gamma_s*delta_ts)
737 precip(j,i,0) = 0.0_dp
740 precip(j,i,n) = precip(j,i,n)*precip_fact
741 precip(j,i,0) = precip(j,i,0) + precip(j,i,n)*inv_twelve
745 #elif (ACCSURFACE==5) 747 precip(j,i,0) = 0.0_dp
751 #if (PRECIP_ANOM_INTERPOL==1) 754 #elif (PRECIP_ANOM_INTERPOL==2) 760 precip(j,i,0) = precip(j,i,0) + precip(j,i,n)*inv_twelve
764 #elif (ACCSURFACE==6) 768 precip(j,i,0) = 0.0_dp
772 + precip_anom_fact*precip_ma_anom(j,i)
775 precip(j,i,0) = precip(j,i,0) + precip(j,i,n)*inv_twelve
783 accum(j,i) = precip(j,i,0)
789 #if (SOLID_PRECIP==1) /* Marsiat (1994) */ 791 if (temp_mm(j,i,n) >= temp_rain)
then 793 else if (temp_mm(j,i,n) <= temp_snow)
then 796 frac_solid = (temp_rain-temp_mm(j,i,n))*inv_delta_temp_rain_snow
799 #elif (SOLID_PRECIP==2) /* Bales et al. (2009) */ 801 if (temp_mm(j,i,n) >= temp_rain)
then 803 else if (temp_mm(j,i,n) <= temp_snow)
then 806 frac_solid = coeff(0) + temp_mm(j,i,n) * ( coeff(1) &
807 + temp_mm(j,i,n) * ( coeff(2) &
808 + temp_mm(j,i,n) * ( coeff(3) &
809 + temp_mm(j,i,n) * ( coeff(4) &
810 + temp_mm(j,i,n) * coeff(5) ) ) ) )
814 #elif (SOLID_PRECIP==3) /* Huybrechts and de Wolde (1999) */ 816 frac_solid = 1.0_dp &
817 - 0.5_dp*erfc((temp_rain-temp_mm(j,i,n))*inv_sqrt2_s_stat)
834 #if (ABLSURFACE==1 || ABLSURFACE==2) 839 temp_mm_help(n) = temp_mm(j,i,n)
842 call pdd(temp_mm_help, s_stat, et(j,i))
849 if ((beta1*et(j,i)) <= (pmax*
snowfall(j,i)))
then 850 melt_star(j,i) = beta1*et(j,i)
855 melt(j,i) = beta2*(et(j,i)-melt_star(j,i)/beta1)
859 #elif (ABLSURFACE==2) 864 melt_star(j,i) =
rainfall(j,i)+beta1*et(j,i)
870 *(et(j,i)-(melt_star(j,i)-
rainfall(j,i))/beta1)
877 melt(j,i) = beta2*et(j,i)
884 #elif (ABLSURFACE==3) 886 temp_jja_help = one_third*(temp_mm(j,i,6)+temp_mm(j,i,7)+temp_mm(j,i,8))
888 melt_star(j,i) = 0.0_dp
889 melt(j,i) = lambda_lti*max((temp_jja_help-temp_lti), 0.0_dp)
901 #if (INITMIP_CONST_SMB==1) /* SMB held constant at its initial distribution */ 909 #if (defined(INITMIP_SMB_ANOM_FILE)) 913 if ((time_in_years > 0.0_dp).and.(time_in_years <= 40.0_dp))
then 914 as_perp =
as_perp + 0.025_dp*floor(time_in_years) * as_anom_initmip
915 else if (time_in_years > 40.0_dp)
then 925 where (melt_star >= melt)
926 temp_s = temp_ma + mu*(melt_star-melt)
939 && (marine_ice_formation==2) \
940 && (marine_ice_calving==9))
947 #if (DISC>0) /* Ice discharge parameterization */ 954 if (firstcall) firstcall = .false.
real(dp), dimension(0:jmax, 0:imax) phi
phi(j,i): Geographic latitude of grid point (i,j)
real(dp) phi_sep_0
PHI_SEP_0: Separation latitude for the computation of the degree-day factors beta1 and beta2: Equator...
integer(i4b) grip_time_stp
grip_time_stp: Time step of the data values for the surface temperature anomaly
integer(i4b) ndata_specmap
ndata_specmap: Number of data values for the sea level
real(dp) beta2_ht_0
BETA2_HT_0: Degree-day factor for ice at high summer temperatures.
real(dp) beta1_ht_0
BETA1_HT_0: Degree-day factor for snow at high summer temperatures.
real(dp) rho_w
RHO_W: Density of pure water.
integer(i4b) specmap_time_max
specmap_time_max: Maximum time of the data values for the sea level
real(dp) temp_precip_time_stp
temp_precip_time_stp: Time step of the surface-temperature and precipitation data ...
real(dp) mu_0
MU_0: Firn-warming correction.
real(dp), parameter eps
eps: Small number
real(dp) temp_precip_time_max
temp_precip_time_max: Maximum time of the surface-temperature and precipitation data ...
real(dp), dimension(0:jmax, 0:imax) calv_grounded
calv_grounded(j,i): Calving rate of grounded ice
real(dp), dimension(0:jmax, 0:imax) snowfall
snowfall(j,i): Snowfall rate at the ice surface
real(dp), dimension(0:jmax, 0:imax) temp_mj_present
temp_mj_present(j,i): Present-day mean summer (northern hemisphere: July, southern hemisphere: Januar...
real(dp), dimension(0:jmax, 0:imax) lambda
lambda(j,i): Geographic longitude of grid point (i,j)
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 ...
integer(i4b) specmap_time_min
specmap_time_min: Minimum time of the data values for the sea level
Computation of the positive degree days (PDD) with statistical temperature fluctuations; based on sem...
real(dp), dimension(0:jmax, 0:imax) as_perp
as_perp(j,i): Accumulation-ablation function at the ice surface (SMB)
real(dp), dimension(:), allocatable glacial_index
glacial_index(n): Data values for the glacial index
real(dp) s_stat_0
S_STAT_0: Standard deviation of the air termperature for the degree-day model.
real(dp) beta2_lt_0
BETA2_LT_0: Degree-day factor for ice at low summer temperatures.
real(dp), dimension(0:jmax, 0:imax) runoff
runoff(j,i): Runoff rate at the ice surface
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 ...
integer(i2b) function, public mask_update_sea_level(z_sl, i, j)
Main function of mask_update_m: Update of the ice-land-ocean mask due to changes of the sea level...
real(dp), parameter pi_180_inv
pi_180_inv: 180 divided by pi (-> rad to deg)
real(dp), dimension(0:jmax, 0:imax) temp_s
temp_s(j,i): Ice surface temperature
integer(i4b) grip_time_max
grip_time_max: Maximum time of the data values for the surface temperature anomaly ...
Declarations of global variables for SICOPOLIS (for the ANT domain).
integer(i4b) ndata_gi
ndata_gi: Number of data values for the glacial index
integer(i4b) ndata_grip
ndata_grip: Number of data values for the surface temperature anomaly
real(dp), dimension(0:jmax, 0:imax) calv_uw_ice
real(dp), parameter pi_180
pi_180: pi divided by 180 (-> deg to rad)
Calving of "underwater ice".
real(dp), dimension(0:jmax, 0:imax), public dis_perp
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) pmax_0
PMAX_0: Saturation factor for the formation of superimposed ice.
real(dp), dimension(0:jmax, 0:imax) evap
evap(j,i): Evaporation rate at the ice surface
real(dp), dimension(0:jmax, 0:imax) accum
accum(j,i): Accumulation rate at the ice surface (includes liquid precipitation = rainfall!) ...
integer(i4b) gi_time_min
gi_time_min: Minimum time of the data values for the glacial index
real(dp) temp_precip_time_min
temp_precip_time_min: Minimum time of the surface-temperature and precipitation data ...
real(dp), dimension(:), allocatable specmap_zsl
specmap_zsl(n): Data values for the sea level
integer(i4b) gi_time_stp
gi_time_stp: Time step of the data values for the glacial index
integer(i4b) ndata_temp_precip
ndata_temp_precip: Number of surface-temperature and precipitation data
real(dp), dimension(0:jmax, 0:imax) temp_ma_present
temp_ma_present(j,i): Present-day mean annual surface temperature
real(dp), dimension(0:jmax, 0:imax) rainfall
rainfall(j,i): Rainfall rate at the ice surface
integer(i4b) gi_time_max
gi_time_max: Maximum time of the data values for the glacial index
real(dp) beta1_lt_0
BETA1_LT_0: Degree-day factor for snow at low summer temperatures.
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) ...
integer(i2b), dimension(0:jmax, 0:imax) maske_neu
maske_neu(j,i): New value of maske computed during an integration step
real(dp), public dt_glann
real(dp), dimension(0:jmax, 0:imax) temp_mj_lgm_anom
temp_mj_lgm_anom(j,i): LGM anomaly (difference LGM - present) of the mean summer (northern hemisphere...
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
real(dp), dimension(0:jmax, 0:imax) temp_ma_lgm_anom
temp_ma_lgm_anom(j,i): LGM anomaly (difference LGM - present) of the mean annual surface temperature ...
subroutine check(status, ch_calling_routine)
NetCDF error capturing.
subroutine, public discharge(dxi, deta)
Ice discharge parameterization main formula, controler (general). [Compute ice discharge via a parame...
real(dp), parameter pi
pi: Constant pi
integer(i4b) ncid_temp_precip
ncid_temp_precip: ID of the NetCDF file containing the surface-temperature and precipitation data as ...
integer(i4b) grip_time_min
grip_time_min: Minimum time of the data values for the surface temperature anomaly ...
subroutine calving_underwater_ice(z_sl)
Main routine: Calving of "underwater ice".
Computation of the surface temperature (must be less than 0 deg C!) and of the accumulation-ablation ...
subroutine, public pdd(temp_mm, s_stat, ET)
Main subroutine of pdd_m: Computation of the positive degree days (PDD) with statistical temperature ...
Ice discharge parameterization for the Greenland ice sheet.
real(dp) rho
RHO: Density of ice.
Update of the ice-land-ocean mask due to changes of the sea level.
integer, parameter dp
Double-precision reals.
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 ...
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...
Declarations of global variables for SICOPOLIS.
real(dp), dimension(:), allocatable griptemp
griptemp(n): Data values for the surface temperature anomaly