42 delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
53 real(dp),
intent(in) :: time, dtime, dxi, deta
55 real(dp),
intent(out) :: delta_ts, glac_index, dzsl_dtau, z_mar
56 real(dp),
intent(inout) :: z_sl
64 integer(i4b) :: i_gr, i_kl
65 integer(i4b) :: ndata_insol
67 real(dp) :: z_sl_min, t1, t2, t3, t4, t5, t6
68 real(dp) :: time_gr, time_kl
69 real(dp) :: z_sle_present, z_sle_help
70 real(dp),
dimension(0:JMAX,0:IMAX) :: temp_ma
71 real(dp) :: eld, g_mb, accum_factor
72 real(dp) :: erosion_chasm
73 real(dp),
dimension(0:JMAX,0:IMAX) :: et
74 real(dp) :: t_obliq_main, t_obliq_mod, &
75 obliq_ampl_max, obliq_ampl_min, obliq_ampl, obliq0, obliq, &
77 real(dp) :: insol_ma_90_now, insol_ma_90_present, &
78 obl_now, obl_present, ecc_now, ecc_present, &
79 ave_now, ave_present, cp_now, cp_present, &
80 temp_ma_90, temp_ma_90_present, zs_90
81 real(dp),
dimension(0:JMAX,0:IMAX) :: dist
82 real(dp) :: ave_data_i_kl, ave_data_i_gr
83 real(dp) :: q_geo_chasm
84 logical,
dimension(0:JMAX,0:IMAX) :: check_point
87 type (ins) :: temp_now, temp_present
90 real(dp),
parameter :: &
91 time_present = 0.0_dp, &
92 zs_90_present = -2.0e+03_dp
94 real(dp),
parameter :: &
99 real(dp),
parameter :: &
100 lambda_h2o = 2.86e+06_dp, &
104 real(dp),
parameter :: &
105 temp_s_min = -125.0_dp
125 delta_ts = -sine_amplit &
126 *cos(2.0_dp*pi*time/(sine_period*year_sec)) &
133 t_obliq_main = 1.25e+05_dp *year_sec
134 t_obliq_mod = 1.3e+06_dp *year_sec
136 obliq0 = 25.2_dp *pi_180
137 obliq_ampl_max = 10.0_dp *pi_180
138 obliq_ampl_min = 2.5_dp *pi_180
140 obliq_ampl = 0.5_dp*(obliq_ampl_max+obliq_ampl_min) &
141 -0.5_dp*(obliq_ampl_max-obliq_ampl_min) &
142 *cos(2.0_dp*pi*time/t_obliq_mod)
144 obliq = obliq0 + obliq_ampl*sin(2.0_dp*pi*time/t_obliq_main)
151 insol_ma_90_now = (sol/pi)*sin(obliq)/sqrt(1.0_dp-ecc**2)
153 insol_ma_90_present = (sol/pi)*sin(obliq0)/sqrt(1.0_dp-ecc0**2)
156 #elif ( TSURFACE==5 || TSURFACE==6 )
160 ndata_insol = (insol_time_max-insol_time_min)/insol_time_stp
162 if (time/year_sec.lt.
real(insol_time_min,dp)) then
164 insol_ma_90_now = insol_ma_90(0)
165 obl_now = obl_data(0)
166 ecc_now = ecc_data(0)
167 ave_now = ave_data(0)
170 else if (time/year_sec.lt.
real(insol_time_max,dp)) then
172 i_kl = floor(((time/year_sec) &
173 -
real(insol_time_min,dp))/
real(insol_time_stp,dp))
176 i_gr = ceiling(((time/year_sec) &
177 -
real(insol_time_min,dp))/
real(insol_time_stp,dp))
178 i_gr = min(i_gr, ndata_insol)
180 if (i_kl.eq.i_gr)
then
182 insol_ma_90_now = insol_ma_90(i_kl)
183 obl_now = obl_data(i_kl)
184 ecc_now = ecc_data(i_kl)
185 ave_now = ave_data(i_kl)
186 cp_now = cp_data(i_kl)
190 time_kl = (insol_time_min + i_kl*insol_time_stp) *year_sec
191 time_gr = (insol_time_min + i_gr*insol_time_stp) *year_sec
193 insol_ma_90_now = insol_ma_90(i_kl) &
194 +(insol_ma_90(i_gr)-insol_ma_90(i_kl)) &
195 *(time-time_kl)/(time_gr-time_kl)
196 obl_now = obl_data(i_kl) &
197 +(obl_data(i_gr)-obl_data(i_kl)) &
198 *(time-time_kl)/(time_gr-time_kl)
199 ecc_now = ecc_data(i_kl) &
200 +(ecc_data(i_gr)-ecc_data(i_kl)) &
201 *(time-time_kl)/(time_gr-time_kl)
203 if ( abs(ave_data(i_gr)-ave_data(i_kl)) < pi )
then
204 ave_data_i_kl = ave_data(i_kl)
205 ave_data_i_gr = ave_data(i_gr)
207 if ( ave_data(i_gr) > ave_data(i_kl) )
then
208 ave_data_i_kl = ave_data(i_kl) + 2.0_dp*pi
209 ave_data_i_gr = ave_data(i_gr)
211 ave_data_i_kl = ave_data(i_kl)
212 ave_data_i_gr = ave_data(i_gr) + 2.0_dp*pi
216 ave_now = ave_data_i_kl &
217 +(ave_data_i_gr-ave_data_i_kl) &
218 *(time-time_kl)/(time_gr-time_kl)
220 cp_now = cp_data(i_kl) &
221 +(cp_data(i_gr)-cp_data(i_kl)) &
222 *(time-time_kl)/(time_gr-time_kl)
228 insol_ma_90_now = insol_ma_90(ndata_insol)
229 obl_now = obl_data(ndata_insol)
230 ecc_now = ecc_data(ndata_insol)
231 ave_now = ave_data(ndata_insol)
232 cp_now = cp_data(ndata_insol)
238 if (time_present/year_sec.lt.
real(insol_time_min,dp)) then
240 insol_ma_90_present = insol_ma_90(0)
241 obl_present = obl_data(0)
242 ecc_present = ecc_data(0)
243 ave_present = ave_data(0)
244 cp_present = cp_data(0)
246 else if (time_present/year_sec.lt.
real(insol_time_max,dp)) then
248 i_kl = floor(((time_present/year_sec) &
249 -
real(insol_time_min,dp))/
real(insol_time_stp,dp))
252 i_gr = ceiling(((time_present/year_sec) &
253 -
real(insol_time_min,dp))/
real(insol_time_stp,dp))
254 i_gr = min(i_gr, ndata_insol)
256 if (i_kl.eq.i_gr)
then
258 insol_ma_90_present = insol_ma_90(i_kl)
259 obl_present = obl_data(i_kl)
260 ecc_present = ecc_data(i_kl)
261 ave_present = ave_data(i_kl)
262 cp_present = cp_data(i_kl)
266 time_kl = (insol_time_min + i_kl*insol_time_stp) *year_sec
267 time_gr = (insol_time_min + i_gr*insol_time_stp) *year_sec
269 insol_ma_90_present = insol_ma_90(i_kl) &
270 +(insol_ma_90(i_gr)-insol_ma_90(i_kl)) &
271 *(time_present-time_kl)/(time_gr-time_kl)
272 obl_present = obl_data(i_kl) &
273 +(obl_data(i_gr)-obl_data(i_kl)) &
274 *(time_present-time_kl)/(time_gr-time_kl)
275 ecc_present = ecc_data(i_kl) &
276 +(ecc_data(i_gr)-ecc_data(i_kl)) &
277 *(time_present-time_kl)/(time_gr-time_kl)
279 if ( abs(ave_data(i_gr)-ave_data(i_kl)) < pi )
then
280 ave_data_i_kl = ave_data(i_kl)
281 ave_data_i_gr = ave_data(i_gr)
283 if ( ave_data(i_gr) > ave_data(i_kl) )
then
284 ave_data_i_kl = ave_data(i_kl) + 2.0_dp*pi
285 ave_data_i_gr = ave_data(i_gr)
287 ave_data_i_kl = ave_data(i_kl)
288 ave_data_i_gr = ave_data(i_gr) + 2.0_dp*pi
292 ave_present = ave_data_i_kl &
293 +(ave_data_i_gr-ave_data_i_kl) &
294 *(time_present-time_kl)/(time_gr-time_kl)
296 cp_present = cp_data(i_kl) &
297 +(cp_data(i_gr)-cp_data(i_kl)) &
298 *(time_present-time_kl)/(time_gr-time_kl)
304 insol_ma_90_present = insol_ma_90(ndata_insol)
305 obl_present = obl_data(ndata_insol)
306 ecc_present = ecc_data(ndata_insol)
307 ave_present = ave_data(ndata_insol)
308 cp_present = cp_data(ndata_insol)
314 #if (TSURFACE==4 || TSURFACE==5)
318 temp_ma_90 = sqrt(sqrt(insol_ma_90_now*(1.0_dp-albedo)/(epsilon*sigma))) &
324 = sqrt(sqrt(insol_ma_90_present*(1.0_dp-albedo)/(epsilon*sigma))) &
329 delta_ts = temp_ma_90 - temp_ma_90_present
336 ecc = ecc_now, ave = ave_now*pi_180_inv, &
337 obl = obl_now*pi_180_inv, sa = albedo, ct = 148.7_dp)
339 temp_ma_90 =
instam(temp_now, 90.0_dp) + delta_ts0 &
345 ecc = ecc_present, ave = ave_present*pi_180_inv, &
346 obl = obl_present*pi_180_inv, sa = albedo, ct = 148.7_dp)
348 temp_ma_90_present =
instam(temp_present, 90.0_dp) + delta_ts0 &
353 delta_ts = temp_ma_90 - temp_ma_90_present
365 stop
' boundary: SEA_LEVEL==2 not allowed for nmars application!'
369 stop
' boundary: SEA_LEVEL==3 not allowed for nmars application!'
375 if ( z_sl_old > -999999.9_dp )
then
376 dzsl_dtau = (z_sl-z_sl_old)/dtime
385 #if ( MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3 )
387 #elif ( MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 )
388 z_mar = fact_z_mar*z_sl
389 #elif ( MARINE_ICE_CALVING==6 || MARINE_ICE_CALVING==7 )
390 if (z_sl >= -80.0_dp)
then
393 z_mar = 10.25_dp*(z_sl+80.0_dp)-200.0_dp
395 z_mar = fact_z_mar*z_mar
407 check_point(j,i) = .false.
413 if (maske(j,i).ge.2)
then
414 check_point(j ,i ) = .true.
415 check_point(j ,i+1) = .true.
416 check_point(j ,i-1) = .true.
417 check_point(j+1,i ) = .true.
418 check_point(j-1,i ) = .true.
425 if (check_point(j,i))
then
435 if (check_point(j,i))
then
436 maske(j,i) = maske_neu(j,i)
443 call
borehole(zs, 0.0_dp, 0.0_dp, dxi, deta,
'grid', zs_90)
451 #if (TSURFACE==1 || TSURFACE==2 || TSURFACE==3)
453 temp_ma_90_present = temp0_ma_90n
455 temp_ma(j,i) = temp_ma_90_present &
456 + (gamma_ma)*(zs(j,i)-zs_90_present) &
457 + (c_ma*pi_180_inv)*((90.0_dp*pi_180)-phi(j,i))
459 #elif (TSURFACE==4 || TSURFACE==5)
461 temp_ma(j,i) = temp_ma_90_present &
462 + (gamma_ma)*(zs(j,i)-zs_90_present) &
463 + (c_ma*pi_180_inv)*((90.0_dp*pi_180)-phi(j,i))
469 #if (TSURFACE==1 || TSURFACE==2 || TSURFACE==3)
471 temp_ma(j,i) = temp_ma(j,i) + delta_ts
473 #elif (TSURFACE==4 || TSURFACE==5)
475 temp_ma(j,i) = temp_ma(j,i) + delta_ts &
476 - (gamma_ma)*(zs_90-zs_90_present)
484 temp_ma(j,i) =
instam(temp_now, phi(j,i)*pi_180_inv) + delta_ts0 &
497 accum_factor = 1.0_dp+gamma_s*delta_ts
499 accum_factor = exp(gamma_s*delta_ts)
501 accum_factor = exp(lambda_h2o/(r_h2o*temp0_ref) &
502 -lambda_h2o/(r_h2o*(temp0_ref+delta_ts)))
507 accum(j,i) = accum_present(j,i)*accum_factor
508 if (accum(j,i).lt.0.0_dp) stop
' boundary: Negative accumulation rate!'
514 #if ( ABLSURFACE==1 || ABLSURFACE==2)
516 eld = eld_0 *1.0e+03_dp
519 g_mb = g_0 *(1.0e-06_dp/year_sec)*(rho_w/rho_i) &
520 *(1.0_dp/(1.0_dp-frac_dust))
523 g_mb = g_0 *(1.0e-06_dp/year_sec)
528 eld = eld*(1.0_dp+gamma_eld*delta_ts)
529 g_mb = g_mb*(gamma_g*accum_factor)
531 eld = eld*exp(gamma_eld*delta_ts)
532 g_mb = g_mb*(gamma_g*accum_factor)
535 if (eld.lt.eps) stop
' boundary: Negative equilibrium line distance eld!'
539 dist(j,i) = sqrt( xi(i)**2 + eta(j)**2 )
548 erosion_chasm = erosion_chasm *(1.0e-03_dp/year_sec)*(rho_w/rho_i) &
549 *(1.0_dp/(1.0_dp-frac_dust))
552 erosion_chasm = erosion_chasm *(1.0e-03_dp/year_sec)
561 #if ( ABLSURFACE==1 || ABLSURFACE==2 )
563 as_perp(j,i) = g_mb*(eld-dist(j,i))
564 if (as_perp(j,i).gt.accum(j,i)) as_perp(j,i) = accum(j,i)
568 evap(j,i) = accum(j,i) - as_perp(j,i)
575 if ( (maske_chasm(j,i) == 7) &
576 .and.(time >= time_chasm_init) &
577 .and.(time <= time_chasm_end) )
then
578 evap(j,i) = evap(j,i) + erosion_chasm
579 as_perp(j,i) = accum(j,i) - evap(j,i)
589 temp_s = min(temp_ma, -eps)
591 temp_s = max(temp_s, temp_s_min)
602 q_geo_chasm = q_geo_chasm *1.0e-03_dp
607 if ( (maske_chasm(j,i) == 7) &
608 .and.(time >= time_chasm_init) &
609 .and.(time <= time_chasm_end) )
then
610 q_geo(j,i) = q_geo_chasm
612 q_geo(j,i) = q_geo_normal(j,i)