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))
73 real(dp),
intent(in) :: time, dtime, dxi, deta
75 real(dp),
intent(out) :: delta_ts, glac_index, dzsl_dtau, z_mar
76 real(dp),
intent(inout) :: z_sl
84 integer(i4b) :: i, j, n
85 integer(i4b) :: i_gr, i_kl
86 integer(i4b) :: nrec_temp_precip
87 integer(i4b),
save :: nrec_temp_precip_save = -1
89 real(dp) :: z_sl_min, t1, t2, t3, t4, t5, t6
90 real(dp) :: year_sec_inv, time_in_years
91 real(dp) :: time_gr, time_kl
92 real(dp) :: z_sle_present, z_sle_help
93 real(dp),
dimension(0:JMAX,0:IMAX,0:12) :: precip
94 real(dp),
dimension(0:JMAX,0:IMAX) :: &
95 snowfall, rainfall, melt, melt_star
96 real(dp),
dimension(0:JMAX,0:IMAX,12) :: temp_mm
97 real(dp),
dimension(0:JMAX,0:IMAX) :: temp_ma, temp_ampl
98 real(dp),
dimension(0:JMAX,0:IMAX) :: temp_ma_anom, temp_mj_anom, precip_ma_anom
99 real(dp),
dimension(0:IMAX,0:JMAX),
save :: temp_ma_anom_tra, temp_mj_anom_tra, &
104 real(dp),
dimension(12) :: temp_mm_help
105 real(dp) :: temp_jja_help
106 real(dp),
dimension(0:JMAX,0:IMAX) :: ET
107 real(dp) :: theta_ma, c_ma, gamma_ma, &
108 theta_ma_1, c_ma_1, gamma_ma_1, &
109 theta_ma_2, c_ma_2, gamma_ma_2, &
110 theta_ma_3, c_ma_3, gamma_ma_3, &
111 zs_sep_1, zs_sep_2, &
112 theta_mj, c_mj, gamma_mj
113 real(dp) :: sine_factor
114 real(dp) :: gamma_p, zs_thresh, &
115 alpha_p, beta_p, temp_0, alpha_t, beta_t, &
116 temp_inv, temp_inv_present, &
117 temp_rain, temp_snow, &
118 inv_delta_temp_rain_snow, coeff(0:5), inv_sqrt2_s_stat, &
119 precip_fact, frac_solid
120 real(dp) :: s_stat, beta1, beta2, Pmax, mu, lambda_lti, temp_lti
121 logical,
dimension(0:JMAX,0:IMAX) :: check_point
122 logical,
save :: firstcall = .true.
127 integer(i4b) :: nc3cor(3)
129 integer(i4b) :: nc3cnt(3)
133 real(dp),
parameter :: &
134 inv_twelve = 1.0_dp/12.0_dp, one_third = 1.0_dp/3.0_dp
136 character(len=64),
parameter :: thisroutine =
'boundary' 138 year_sec_inv = 1.0_dp/year_sec
139 time_in_years = time*year_sec_inv
158 delta_ts = sine_amplit &
159 *cos(2.0_dp*
pi*time/(sine_period*year_sec)) &
170 i_kl = floor(((time/year_sec) &
174 i_gr = ceiling(((time/year_sec) &
178 if (i_kl == i_gr)
then 189 *(time-time_kl)/(time_gr-time_kl)
198 delta_ts = delta_ts * grip_temp_fact
207 else if (time/year_sec <
real(
gi_time_max,dp)) then
209 i_kl = floor(((time/year_sec) &
213 i_gr = ceiling(((time/year_sec) &
217 if (i_kl == i_gr)
then 223 time_kl = (
gi_time_min + i_kl*gi_time_stp) *year_sec
224 time_gr = (
gi_time_min + i_gr*gi_time_stp) *year_sec
228 *(time-time_kl)/(time_gr-time_kl)
240 #elif ( (TSURFACE==6) && (ACCSURFACE==6) ) 251 if ( nrec_temp_precip < 0 )
then 252 stop
' >>> boundary: nrec_temp_precip < 0, not allowed!' 254 stop
' >>> boundary: nrec_temp_precip > ndata_temp_precip, not allowed!' 257 if ( nrec_temp_precip /= nrec_temp_precip_save )
then 261 nc3cor(3) = nrec_temp_precip + 1
269 start=nc3cor, count=nc3cnt), thisroutine )
274 start=nc3cor, count=nc3cnt), thisroutine )
279 start=nc3cor, count=nc3cnt), thisroutine )
283 temp_ma_anom = transpose(temp_ma_anom_tra)
284 temp_mj_anom = transpose(temp_mj_anom_tra)
285 precip_ma_anom = transpose(precip_ma_anom_tra) *(1.0e-03_dp/year_sec)*(
rho_w/
rho)
288 nrec_temp_precip_save = nrec_temp_precip
303 t1 = -250000.0_dp *year_sec
304 t2 = -140000.0_dp *year_sec
305 t3 = -125000.0_dp *year_sec
306 t4 = -21000.0_dp *year_sec
307 t5 = -8000.0_dp *year_sec
308 t6 = 0.0_dp *year_sec
312 else if (time < t2)
then 313 z_sl = z_sl_min*(time-t1)/(t2-t1)
314 else if (time < t3)
then 315 z_sl = -z_sl_min*(time-t3)/(t3-t2)
316 else if (time < t4)
then 317 z_sl = z_sl_min*(time-t3)/(t4-t3)
318 else if (time < t5)
then 319 z_sl = -z_sl_min*(time-t5)/(t5-t4)
320 else if (time < t6)
then 334 i_kl = floor(((time/year_sec) &
338 i_gr = ceiling(((time/year_sec) &
342 if (i_kl == i_gr)
then 353 *(time-time_kl)/(time_gr-time_kl)
366 if ( z_sl_old > -999999.9_dp )
then 367 dzsl_dtau = (z_sl-z_sl_old)/dtime
376 #if ( MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3 ) 378 #elif ( MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 ) 379 z_mar = fact_z_mar*z_sl
380 #elif ( MARINE_ICE_CALVING==6 || MARINE_ICE_CALVING==7 ) 381 if (z_sl >= -80.0_dp)
then 384 z_mar = 10.25_dp*(z_sl+80.0_dp)-200.0_dp
386 z_mar = fact_z_mar*z_mar
398 check_point(j,i) = .false.
404 if (
maske(j,i) >= 2)
then 405 check_point(j ,i ) = .true.
406 check_point(j ,i+1) = .true.
407 check_point(j ,i-1) = .true.
408 check_point(j+1,i ) = .true.
409 check_point(j-1,i ) = .true.
416 if (check_point(j,i))
then 426 if (check_point(j,i))
then 434 #if TEMP_PRESENT_PARA == 1 /* Parameterisation by Fortuin and Oerlemans */ 438 gamma_ma = -9.14e-03_dp
442 gamma_mj = -6.92e-03_dp
445 #elif TEMP_PRESENT_PARA == 2 /* Parameterisation by Fortuin and Oerlemans */ 452 theta_ma_1 = 49.642_dp
456 theta_ma_2 = 36.689_dp
457 gamma_ma_2 = -5.102e-03_dp
460 theta_ma_3 = 7.405_dp
461 gamma_ma_3 = -14.285e-03_dp
465 gamma_mj = -6.92e-03_dp
470 stop
' >>> boundary: Parameter TEMP_PRESENT_PARA must be either 1 or 2!' 479 #if TEMP_PRESENT_PARA == 1 482 #elif TEMP_PRESENT_PARA == 2 483 if (
zs(j,i) <= zs_sep_1 )
then 486 else if (
zs(j,i) <= zs_sep_2 )
then 507 #elif (TSURFACE == 5) 514 #elif (TSURFACE == 6) 519 temp_ma(j,i) =
temp_ma_present(j,i) + temp_ma_anom_fact*temp_ma_anom(j,i)
520 temp_mm(j,i,7) =
temp_mj_present(j,i) + temp_mj_anom_fact*temp_mj_anom(j,i)
526 temp_ampl(j,i) = temp_mm(j,i,7) - temp_ma(j,i)
528 if (temp_ampl(j,i) <
eps)
then 539 sine_factor = sin((
real(n,dp)-4.0_dp)*
pi/6.0_dp)
543 temp_mm(j,i,n) = temp_ma(j,i) + sine_factor*temp_ampl(j,i)
551 #if (INITMIP_CONST_SMB==1) /* SMB held constant at its initial distribution */ 555 #if (ACCSURFACE <= 3) 557 #if (ELEV_DESERT == 1) 559 gamma_p = gamma_p*1.0e-03_dp
561 zs_thresh = zs_thresh
565 #elif (ACCSURFACE == 4) 575 #if (SOLID_PRECIP == 1) /* Marsiat (1994) */ 582 inv_delta_temp_rain_snow = 1.0_dp/(temp_rain-temp_snow)
584 #elif (SOLID_PRECIP == 2) /* Bales et al. (2009) */ 591 coeff(0) = 5.4714e-01_dp
592 coeff(1) = -9.1603e-02_dp
593 coeff(2) = -3.314e-03_dp
594 coeff(3) = 4.66e-04_dp
595 coeff(4) = 3.8e-05_dp
596 coeff(5) = 6.0e-07_dp
598 #elif (SOLID_PRECIP == 3) /* Huybrechts and de Wolde (1999) */ 602 temp_snow = temp_rain
608 inv_sqrt2_s_stat = 1.0_dp/(sqrt(2.0_dp)*s_stat)
612 #if (ABLSURFACE==1 || ABLSURFACE==2) 623 #elif (ABLSURFACE==3) 625 lambda_lti = lambda_lti *(0.001_dp/year_sec)*(
rho_w/
rho)
637 #if (ACCSURFACE <= 3) 641 #if (ELEV_DESERT == 0) 645 #elif (ELEV_DESERT == 1) 647 if (
zs_ref(j,i) < zs_thresh)
then 649 = exp(gamma_p*(max(
zs(j,i),zs_thresh)-zs_thresh))
652 = exp(gamma_p*(max(
zs(j,i),zs_thresh)-
zs_ref(j,i)))
656 stop
' >>> boundary: Parameter ELEV_DESERT must be either 0 or 1!' 668 precip_fact = accfact
670 precip_fact = 1.0_dp + gamma_s*delta_ts
672 precip_fact = exp(gamma_s*delta_ts)
675 #if (ACCSURFACE <= 3) 677 precip(j,i,0) = 0.0_dp
680 precip(j,i,n) = precip(j,i,n)*precip_fact
681 precip(j,i,0) = precip(j,i,0) + precip(j,i,n)*inv_twelve
685 #elif (ACCSURFACE == 4) 687 precip(j,i,0) = 0.0_dp
689 temp_inv = alpha_t * (temp_ma(j,i)+temp_0) + beta_t
692 precip_fact = exp(alpha_p*(temp_0/temp_inv_present-temp_0/temp_inv)) &
693 *(temp_inv_present/temp_inv)**2 &
694 *(1.0_dp+beta_p*(temp_inv-temp_inv_present))
698 precip(j,i,0) = precip(j,i,0) + precip(j,i,n)*inv_twelve
702 #elif (ACCSURFACE == 5) 704 precip(j,i,0) = 0.0_dp
708 #if (PRECIP_ANOM_INTERPOL==1) 711 #elif (PRECIP_ANOM_INTERPOL==2) 717 precip(j,i,0) = precip(j,i,0) + precip(j,i,n)*inv_twelve
721 #elif (ACCSURFACE == 6) 729 precip(j,i,n) = precip(j,i,0)
737 accum(j,i) = precip(j,i,0)
739 snowfall(j,i) = 0.0_dp
743 #if (SOLID_PRECIP == 1) /* Marsiat (1994) */ 745 if (temp_mm(j,i,n) >= temp_rain)
then 747 else if (temp_mm(j,i,n) <= temp_snow)
then 750 frac_solid = (temp_rain-temp_mm(j,i,n))*inv_delta_temp_rain_snow
753 #elif (SOLID_PRECIP == 2) /* Bales et al. (2009) */ 755 if (temp_mm(j,i,n) >= temp_rain)
then 757 else if (temp_mm(j,i,n) <= temp_snow)
then 760 frac_solid = coeff(0) + temp_mm(j,i,n) * ( coeff(1) &
761 + temp_mm(j,i,n) * ( coeff(2) &
762 + temp_mm(j,i,n) * ( coeff(3) &
763 + temp_mm(j,i,n) * ( coeff(4) &
764 + temp_mm(j,i,n) * coeff(5) ) ) ) )
768 #elif (SOLID_PRECIP == 3) /* Huybrechts and de Wolde (1999) */ 770 frac_solid = 1.0_dp &
771 - 0.5_dp*erfc((temp_rain-temp_mm(j,i,n))*inv_sqrt2_s_stat)
775 snowfall(j,i) = snowfall(j,i) + precip(j,i,n)*frac_solid*inv_twelve
779 rainfall(j,i) = precip(j,i,0) - snowfall(j,i)
781 if (snowfall(j,i) < 0.0_dp) snowfall(j,i) = 0.0_dp
782 if (rainfall(j,i) < 0.0_dp) rainfall(j,i) = 0.0_dp
788 #if (ABLSURFACE==1 || ABLSURFACE==2) 793 temp_mm_help(n) = temp_mm(j,i,n)
796 call pdd(temp_mm_help, s_stat, et(j,i))
803 if ((beta1*et(j,i)) <= (pmax*snowfall(j,i)))
then 804 melt_star(j,i) = beta1*et(j,i)
806 runoff(j,i) = melt(j,i)+rainfall(j,i)
808 melt_star(j,i) = pmax*snowfall(j,i)
809 melt(j,i) = beta2*(et(j,i)-melt_star(j,i)/beta1)
810 runoff(j,i) = melt(j,i)+rainfall(j,i)
813 #elif (ABLSURFACE==2) 815 if ( rainfall(j,i) <= (pmax*snowfall(j,i)) )
then 817 if ( (rainfall(j,i)+beta1*et(j,i)) <= (pmax*snowfall(j,i)) )
then 818 melt_star(j,i) = rainfall(j,i)+beta1*et(j,i)
822 melt_star(j,i) = pmax*snowfall(j,i)
824 *(et(j,i)-(melt_star(j,i)-rainfall(j,i))/beta1)
830 melt_star(j,i) = pmax*snowfall(j,i)
831 melt(j,i) = beta2*et(j,i)
832 runoff(j,i) = melt(j,i) + rainfall(j,i)-pmax*snowfall(j,i)
838 #elif (ABLSURFACE==3) 840 temp_jja_help = one_third*(temp_mm(j,i,6)+temp_mm(j,i,7)+temp_mm(j,i,8))
842 melt_star(j,i) = 0.0_dp
843 melt(j,i) = lambda_lti*max((temp_jja_help-temp_lti), 0.0_dp)
844 runoff(j,i) = melt(j,i) + rainfall(j,i)
855 #if (INITMIP_CONST_SMB==1) /* SMB held constant at its initial distribution */ 863 #if (defined(INITMIP_SMB_ANOM_FILE)) 867 if ((time_in_years > 0.0_dp).and.(time_in_years <= 40.0_dp))
then 868 as_perp =
as_perp + 0.025_dp*floor(time_in_years) * as_anom_initmip
869 else if (time_in_years > 40.0_dp)
then 879 where (melt_star >= melt)
880 temp_s = temp_ma + mu*(melt_star-melt)
893 && (marine_ice_formation==2) \
894 && (marine_ice_calving==9))
901 if (firstcall) firstcall = .false.
real(dp), dimension(0:jmax, 0:imax) phi
phi(j,i): Geographic latitude of grid point (i,j)
integer(i4b) ndata_specmap
ndata_specmap: Number of data values for the sea level
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) temp_mj_present
temp_mj_present(j,i): Present-day mean summer (northern hemisphere: July, southern hemisphere: Januar...
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), 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) beta1_0
BETA1_0: Degree-day factor for snow.
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!) ...
real(dp), dimension(0:jmax, 0:imax) precip_ma_present
precip_ma_present(j,i): Present-day mean annual precipitation rate at the ice surface ...
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
real(dp) beta2_0
BETA2_0: Degree-day factor for ice.
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
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(i2b), 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) 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...
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.
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 ...
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