55 subroutine boundary(time, dtime, dxi, deta, &
56 delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
59 && (marine_ice_formation==2) \
60 && (marine_ice_calving==9))
69 real(dp),
intent(in) :: time, dtime, dxi, deta
71 real(dp),
intent(out) :: delta_ts, glac_index, dzsl_dtau, z_mar
72 real(dp),
intent(inout) :: z_sl
80 integer(i4b) :: i, j, n
81 integer(i4b) :: i_gr, i_kl
83 real(dp) :: z_sl_min, t1, t2, t3, t4, t5, t6
84 real(dp) :: time_gr, time_kl
85 real(dp) :: z_sle_present, z_sle_help
86 real(dp),
dimension(0:JMAX,0:IMAX,0:12) :: precip
87 real(dp),
dimension(0:JMAX,0:IMAX,12) :: temp_mm
88 real(dp),
dimension(0:JMAX,0:IMAX) :: temp_ma
89 real(dp),
dimension(12) :: temp_mm_help
90 real(dp) :: temp_jja_help
91 real(dp) :: gamma_t, temp_diff
92 real(dp) :: gamma_p, zs_thresh, &
93 temp_rain, temp_snow, &
94 inv_delta_temp_rain_snow, coeff(0:5), inv_sqrt2_s_stat, &
95 precip_fact, frac_solid
96 real(dp) :: s_stat, beta1, beta2, pmax, mu, lambda_lti, temp_lti
97 logical,
dimension(0:JMAX,0:IMAX) :: check_point
98 logical,
save :: firstcall = .true.
100 real(dp),
parameter :: &
101 inv_twelve = 1.0_dp/12.0_dp, one_third = 1.0_dp/3.0_dp
120 delta_ts = sine_amplit &
121 *cos(2.0_dp*
pi*time/(sine_period*year_sec)) &
132 i_kl = floor(((time/year_sec) &
136 i_gr = ceiling(((time/year_sec) &
140 if (i_kl.eq.i_gr)
then 151 *(time-time_kl)/(time_gr-time_kl)
160 delta_ts = delta_ts * grip_temp_fact
171 i_kl = floor(((time/year_sec) &
175 i_gr = ceiling(((time/year_sec) &
179 if (i_kl == i_gr)
then 190 *(time-time_kl)/(time_gr-time_kl)
212 t1 = -250000.0_dp *year_sec
213 t2 = -140000.0_dp *year_sec
214 t3 = -125000.0_dp *year_sec
215 t4 = -21000.0_dp *year_sec
216 t5 = -8000.0_dp *year_sec
217 t6 = 0.0_dp *year_sec
221 else if (time.lt.t2)
then 222 z_sl = z_sl_min*(time-t1)/(t2-t1)
223 else if (time.lt.t3)
then 224 z_sl = -z_sl_min*(time-t3)/(t3-t2)
225 else if (time.lt.t4)
then 226 z_sl = z_sl_min*(time-t3)/(t4-t3)
227 else if (time.lt.t5)
then 228 z_sl = -z_sl_min*(time-t5)/(t5-t4)
229 else if (time.lt.t6)
then 243 i_kl = floor(((time/year_sec) &
247 i_gr = ceiling(((time/year_sec) &
251 if (i_kl.eq.i_gr)
then 262 *(time-time_kl)/(time_gr-time_kl)
275 if ( z_sl_old > -999999.9_dp )
then 276 dzsl_dtau = (z_sl-z_sl_old)/dtime
285 #if ( MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3 ) 287 #elif ( MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 ) 288 z_mar = fact_z_mar*z_sl
289 #elif ( MARINE_ICE_CALVING==6 || MARINE_ICE_CALVING==7 ) 290 if (z_sl >= -80.0_dp)
then 293 z_mar = 10.25_dp*(z_sl+80.0_dp)-200.0_dp
295 z_mar = fact_z_mar*z_mar
307 check_point(j,i) = .false.
313 if (
maske(j,i) >= 2_i1b)
then 314 check_point(j ,i ) = .true.
315 check_point(j ,i+1) = .true.
316 check_point(j ,i-1) = .true.
317 check_point(j+1,i ) = .true.
318 check_point(j-1,i ) = .true.
325 if (check_point(j,i))
then 335 if (check_point(j,i))
then 343 gamma_t = -6.5e-03_dp
353 temp_diff = gamma_t*(
zs(j,i)-
zs_ref(j,i)) + delta_ts
359 #elif (TSURFACE == 5) 364 temp_diff = gamma_t*(
zs(j,i)-
zs_ref(j,i))
376 temp_ma(j,i) = 0.0_dp
379 temp_ma(j,i) = temp_ma(j,i) + temp_mm(j,i,n)*inv_twelve
387 #if (ELEV_DESERT == 1) 389 gamma_p = gamma_p*1.0e-03_dp
391 zs_thresh = zs_thresh
395 #if (SOLID_PRECIP == 1) /* Marsiat (1994) */ 402 inv_delta_temp_rain_snow = 1.0_dp/(temp_rain-temp_snow)
404 #elif (SOLID_PRECIP == 2) /* Bales et al. (2009) */ 411 coeff(0) = 5.4714e-01_dp
412 coeff(1) = -9.1603e-02_dp
413 coeff(2) = -3.314e-03_dp
414 coeff(3) = 4.66e-04_dp
415 coeff(4) = 3.8e-05_dp
416 coeff(5) = 6.0e-07_dp
418 #elif (SOLID_PRECIP == 3) /* Huybrechts and de Wolde (1999) */ 422 temp_snow = temp_rain
428 inv_sqrt2_s_stat = 1.0_dp/(sqrt(2.0_dp)*s_stat)
432 #if (ABLSURFACE==1 || ABLSURFACE==2) 443 #elif (ABLSURFACE==3) 445 lambda_lti = lambda_lti *(0.001_dp/year_sec)*(
rho_w/
rho)
457 #if (ACCSURFACE <= 3) 461 #if (ELEV_DESERT == 0) 465 #elif (ELEV_DESERT == 1) 467 if (
zs_ref(j,i) < zs_thresh)
then 469 = exp(gamma_p*(max(
zs(j,i),zs_thresh)-zs_thresh))
472 = exp(gamma_p*(max(
zs(j,i),zs_thresh)-
zs_ref(j,i)))
476 errormsg =
' >>> boundary: Parameter ELEV_DESERT must be either 0 or 1!' 489 precip_fact = accfact
491 precip_fact = 1.0_dp + gamma_s*delta_ts
493 precip_fact = exp(gamma_s*delta_ts)
496 #if (ACCSURFACE <= 3) 498 precip(j,i,0) = 0.0_dp
501 precip(j,i,n) = precip(j,i,n)*precip_fact
502 precip(j,i,0) = precip(j,i,0) + precip(j,i,n)*inv_twelve
506 #elif (ACCSURFACE == 5) 508 precip(j,i,0) = 0.0_dp
512 #if (PRECIP_ANOM_INTERPOL==1) 515 #elif (PRECIP_ANOM_INTERPOL==2) 521 precip(j,i,0) = precip(j,i,0) + precip(j,i,n)*inv_twelve
529 accum(j,i) = precip(j,i,0)
535 #if (SOLID_PRECIP == 1) /* Marsiat (1994) */ 537 if (temp_mm(j,i,n) >= temp_rain)
then 539 else if (temp_mm(j,i,n) <= temp_snow)
then 542 frac_solid = (temp_rain-temp_mm(j,i,n))*inv_delta_temp_rain_snow
545 #elif (SOLID_PRECIP == 2) /* Bales et al. (2009) */ 547 if (temp_mm(j,i,n) >= temp_rain)
then 549 else if (temp_mm(j,i,n) <= temp_snow)
then 552 frac_solid = coeff(0) + temp_mm(j,i,n) * ( coeff(1) &
553 + temp_mm(j,i,n) * ( coeff(2) &
554 + temp_mm(j,i,n) * ( coeff(3) &
555 + temp_mm(j,i,n) * ( coeff(4) &
556 + temp_mm(j,i,n) * coeff(5) ) ) ) )
560 #elif (SOLID_PRECIP == 3) /* Huybrechts and de Wolde (1999) */ 562 frac_solid = 1.0_dp &
563 - 0.5_dp*erfc((temp_rain-temp_mm(j,i,n))*inv_sqrt2_s_stat)
580 #if (ABLSURFACE==1 || ABLSURFACE==2) 585 temp_mm_help(n) = temp_mm(j,i,n)
588 call pdd(temp_mm_help, s_stat,
et(j,i))
595 if ((beta1*
et(j,i)) <= (pmax*
snowfall(j,i)))
then 605 #elif (ABLSURFACE==2) 630 #elif (ABLSURFACE==3) 632 temp_jja_help = one_third*(temp_mm(j,i,6)+temp_mm(j,i,7)+temp_mm(j,i,8))
635 melt(j,i) = lambda_lti*max((temp_jja_help-temp_lti), 0.0_dp)
669 && (marine_ice_formation==2) \
670 && (marine_ice_calving==9))
677 if (firstcall) firstcall = .false.
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), dimension(0:jmax, 0:imax, 12) temp_mm_present
temp_mm_present(j,i,n): Present-day mean monthly surface temperature
real(dp) rho_w
RHO_W: Density of pure water.
subroutine error(error_message)
Main routine of error_m: Writing of error messages and stopping execution.
integer(i4b) specmap_time_max
specmap_time_max: Maximum time of the data values for the sea level
real(dp) mu_0
MU_0: Firn-warming correction.
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) rainfall
rainfall(j,i): Rainfall rate at the ice surface
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...
integer(i1b), dimension(0:jmax, 0:imax) maske_neu
maske_neu(j,i): New value of maske computed during an integration step
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), dimension(0:jmax, 0:imax) runoff
runoff(j,i): Runoff rate at the ice surface
real(dp), dimension(0:jmax, 0:imax) melt_star
melt_star(j,i): Superimposed ice formation 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 ...
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_s
temp_s(j,i): Ice surface temperature
real(dp), dimension(0:jmax, 0:imax) melt
melt(j,i): Melting rate at the ice surface
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) beta1_0
BETA1_0: Degree-day factor for snow.
real(dp), dimension(0:jmax, 0:imax) et
ET(j,i): Temperature excess at the ice surface (positive degree days divided by time) ...
real(dp), dimension(0:jmax, 0:imax) calv_uw_ice
Calving of "underwater ice".
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(i1b) 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...
integer(i4b) gi_time_min
gi_time_min: Minimum time of the data values for the glacial index
real(dp), dimension(:), allocatable specmap_zsl
specmap_zsl(n): Data values for the sea level
real(dp) beta2_0
BETA2_0: Degree-day factor for ice.
integer(i4b) gi_time_stp
gi_time_stp: Time step of the data values for the glacial index
Writing of error messages and stopping execution.
integer(i4b) gi_time_max
gi_time_max: Maximum time of the data values for the glacial index
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(i1b), 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 ...
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), parameter pi
pi: Constant pi
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 ...
character(len=256) errormsg
errormsg: Error-message string
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...
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...
Declarations of global variables for SICOPOLIS.
real(dp), dimension(:), allocatable griptemp
griptemp(n): Data values for the surface temperature anomaly