38 delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
50 real(dp),
intent(in) :: time, dtime, dxi, deta
52 real(dp),
intent(out) :: delta_ts, glac_index, dzsl_dtau, z_mar
53 real(dp),
intent(inout) :: z_sl
60 integer(i4b) :: i, j, n
61 integer(i4b) :: i_gr, i_kl
62 integer(i4b) :: nrec_temp_precip
63 integer(i4b),
save :: nrec_temp_precip_save = -1
65 real(dp) :: z_sl_min, t1, t2, t3, t4, t5, t6
66 real(dp) :: time_gr, time_kl
67 real(dp) :: z_sle_present, z_sle_help
68 real(dp),
dimension(0:JMAX,0:IMAX,0:12) :: precip
69 real(dp),
dimension(0:JMAX,0:IMAX) :: &
70 snowfall, rainfall, melt, melt_star
71 real(dp),
dimension(0:JMAX,0:IMAX,12) :: temp_mm
72 real(dp),
dimension(0:JMAX,0:IMAX) :: temp_ma, temp_ampl
73 real(dp),
dimension(0:JMAX,0:IMAX) :: temp_ma_anom, temp_mj_anom, precip_ma_anom
74 real(dp),
dimension(0:IMAX,0:JMAX),
save :: temp_ma_anom_tra, temp_mj_anom_tra, &
79 real(dp),
dimension(12) :: temp_mm_help
80 real(dp) :: temp_jja_help
81 real(dp),
dimension(0:JMAX,0:IMAX) :: et
82 real(dp) :: theta_ma, c_ma, gamma_ma, &
83 theta_ma_1, c_ma_1, gamma_ma_1, &
84 theta_ma_2, c_ma_2, gamma_ma_2, &
85 theta_ma_3, c_ma_3, gamma_ma_3, &
87 theta_mj, c_mj, gamma_mj
88 real(dp) :: sine_factor
89 real(dp) :: gamma_p, zs_thresh, &
90 alpha_p, beta_p, temp_0, alpha_t, beta_t, &
91 temp_inv, temp_inv_present, &
92 temp_rain, temp_snow, &
93 inv_delta_temp_rain_snow, coeff(0:5), inv_sqrt2_s_stat, &
94 precip_fact, frac_solid
95 real(dp) :: s_stat, beta1, beta2, pmax, mu, lambda_lti, temp_lti
96 logical,
dimension(0:JMAX,0:IMAX) :: check_point
101 integer(i4b) :: nc3cor(3)
103 integer(i4b) :: nc3cnt(3)
107 real(dp),
parameter :: &
108 inv_twelve = 1.0_dp/12.0_dp, one_third = 1.0_dp/3.0_dp
110 character(len=64),
parameter :: thisroutine =
'boundary'
129 delta_ts = sine_amplit &
130 *cos(2.0_dp*pi*time/(sine_period*year_sec)) &
137 if (time/year_sec <
real(grip_time_min,dp)) then
138 delta_ts = griptemp(0)
139 else if (time/year_sec <
real(grip_time_max,dp)) then
141 i_kl = floor(((time/year_sec) &
142 -
real(grip_time_min,dp))/
real(grip_time_stp,dp))
145 i_gr = ceiling(((time/year_sec) &
146 -
real(grip_time_min,dp))/
real(grip_time_stp,dp))
147 i_gr = min(i_gr, ndata_grip)
149 if (i_kl == i_gr)
then
151 delta_ts = griptemp(i_kl)
155 time_kl = (grip_time_min + i_kl*grip_time_stp) *year_sec
156 time_gr = (grip_time_min + i_gr*grip_time_stp) *year_sec
158 delta_ts = griptemp(i_kl) &
159 +(griptemp(i_gr)-griptemp(i_kl)) &
160 *(time-time_kl)/(time_gr-time_kl)
166 delta_ts = griptemp(ndata_grip)
169 delta_ts = delta_ts * grip_temp_fact
176 if (time/year_sec <
real(gi_time_min,dp)) then
177 glac_index = glacial_index(0)
178 else if (time/year_sec <
real(gi_time_max,dp)) then
180 i_kl = floor(((time/year_sec) &
181 -
real(gi_time_min,dp))/
real(gi_time_stp,dp))
184 i_gr = ceiling(((time/year_sec) &
185 -
real(gi_time_min,dp))/
real(gi_time_stp,dp))
186 i_gr = min(i_gr, ndata_gi)
188 if (i_kl == i_gr)
then
190 glac_index = glacial_index(i_kl)
194 time_kl = (gi_time_min + i_kl*gi_time_stp) *year_sec
195 time_gr = (gi_time_min + i_gr*gi_time_stp) *year_sec
197 glac_index = glacial_index(i_kl) &
198 +(glacial_index(i_gr)-glacial_index(i_kl)) &
199 *(time-time_kl)/(time_gr-time_kl)
205 glac_index = glacial_index(ndata_gi)
211 #elif ( (TSURFACE==6) && (ACCSURFACE==6) )
213 if (time/year_sec <= temp_precip_time_min)
then
215 else if (time/year_sec < temp_precip_time_max)
then
216 nrec_temp_precip = nint( ((time/year_sec)-temp_precip_time_min) &
217 / temp_precip_time_stp )
219 nrec_temp_precip = ndata_temp_precip
222 if ( nrec_temp_precip < 0 )
then
223 stop
' boundary: nrec_temp_precip < 0, not allowed!'
224 else if ( nrec_temp_precip > ndata_temp_precip )
then
225 stop
' boundary: nrec_temp_precip > ndata_temp_precip, not allowed!'
228 if ( nrec_temp_precip /= nrec_temp_precip_save )
then
232 nc3cor(3) = nrec_temp_precip + 1
237 call
check( nf90_inq_varid(ncid_temp_precip,
'annualtemp_anom', ncv), &
239 call
check( nf90_get_var(ncid_temp_precip, ncv, temp_ma_anom_tra, &
240 start=nc3cor, count=nc3cnt), thisroutine )
242 call
check( nf90_inq_varid(ncid_temp_precip,
'januarytemp_anom', ncv), &
244 call
check( nf90_get_var(ncid_temp_precip, ncv, temp_mj_anom_tra, &
245 start=nc3cor, count=nc3cnt), thisroutine )
247 call
check( nf90_inq_varid(ncid_temp_precip,
'precipitation_anom', ncv), &
249 call
check( nf90_get_var(ncid_temp_precip, ncv, precip_ma_anom_tra, &
250 start=nc3cor, count=nc3cnt), thisroutine )
254 temp_ma_anom = transpose(temp_ma_anom_tra)
255 temp_mj_anom = transpose(temp_mj_anom_tra)
256 precip_ma_anom = transpose(precip_ma_anom_tra) *(1.0e-03_dp/year_sec)*(rho_w/rho)
259 nrec_temp_precip_save = nrec_temp_precip
274 t1 = -250000.0_dp *year_sec
275 t2 = -140000.0_dp *year_sec
276 t3 = -125000.0_dp *year_sec
277 t4 = -21000.0_dp *year_sec
278 t5 = -8000.0_dp *year_sec
279 t6 = 0.0_dp *year_sec
283 else if (time < t2)
then
284 z_sl = z_sl_min*(time-t1)/(t2-t1)
285 else if (time < t3)
then
286 z_sl = -z_sl_min*(time-t3)/(t3-t2)
287 else if (time < t4)
then
288 z_sl = z_sl_min*(time-t3)/(t4-t3)
289 else if (time < t5)
then
290 z_sl = -z_sl_min*(time-t5)/(t5-t4)
291 else if (time < t6)
then
301 if (time/year_sec <
real(specmap_time_min,dp)) then
302 z_sl = specmap_zsl(0)
303 else if (time/year_sec <
real(specmap_time_max,dp)) then
305 i_kl = floor(((time/year_sec) &
306 -
real(specmap_time_min,dp))/
real(specmap_time_stp,dp))
309 i_gr = ceiling(((time/year_sec) &
310 -
real(specmap_time_min,dp))/
real(specmap_time_stp,dp))
311 i_gr = min(i_gr, ndata_specmap)
313 if (i_kl == i_gr)
then
315 z_sl = specmap_zsl(i_kl)
319 time_kl = (specmap_time_min + i_kl*specmap_time_stp) *year_sec
320 time_gr = (specmap_time_min + i_gr*specmap_time_stp) *year_sec
322 z_sl = specmap_zsl(i_kl) &
323 +(specmap_zsl(i_gr)-specmap_zsl(i_kl)) &
324 *(time-time_kl)/(time_gr-time_kl)
330 z_sl = specmap_zsl(ndata_specmap)
337 if ( z_sl_old > -999999.9_dp )
then
338 dzsl_dtau = (z_sl-z_sl_old)/dtime
347 #if ( MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3 )
349 #elif ( MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 )
350 z_mar = fact_z_mar*z_sl
351 #elif ( MARINE_ICE_CALVING==6 || MARINE_ICE_CALVING==7 )
352 if (z_sl >= -80.0_dp)
then
355 z_mar = 10.25_dp*(z_sl+80.0_dp)-200.0_dp
357 z_mar = fact_z_mar*z_mar
369 check_point(j,i) = .false.
375 if (maske(j,i) >= 2)
then
376 check_point(j ,i ) = .true.
377 check_point(j ,i+1) = .true.
378 check_point(j ,i-1) = .true.
379 check_point(j+1,i ) = .true.
380 check_point(j-1,i ) = .true.
387 if (check_point(j,i))
then
397 if (check_point(j,i))
then
398 maske(j,i) = maske_neu(j,i)
405 #if TEMP_PRESENT_PARA == 1 /* Parameterisation by Fortuin and Oerlemans */
409 gamma_ma = -9.14e-03_dp
413 gamma_mj = -6.92e-03_dp
416 #elif TEMP_PRESENT_PARA == 2 /* Parameterisation by Fortuin and Oerlemans */
423 theta_ma_1 = 49.642_dp
427 theta_ma_2 = 36.689_dp
428 gamma_ma_2 = -5.102e-03_dp
431 theta_ma_3 = 7.405_dp
432 gamma_ma_3 = -14.285e-03_dp
436 gamma_mj = -6.92e-03_dp
441 stop
' boundary: Parameter TEMP_PRESENT_PARA must be either 1 or 2!'
450 #if TEMP_PRESENT_PARA == 1
451 temp_ma_present(j,i) = theta_ma + gamma_ma*zs(j,i) &
452 + c_ma*abs(phi(j,i))*pi_180_inv
453 #elif TEMP_PRESENT_PARA == 2
454 if ( zs(j,i) <= zs_sep_1 )
then
455 temp_ma_present(j,i) = theta_ma_1 + gamma_ma_1*zs(j,i) &
456 + c_ma_1*abs(phi(j,i))*pi_180_inv
457 else if ( zs(j,i) <= zs_sep_2 )
then
458 temp_ma_present(j,i) = theta_ma_2 + gamma_ma_2*zs(j,i) &
459 + c_ma_2*abs(phi(j,i))*pi_180_inv
461 temp_ma_present(j,i) = theta_ma_3 + gamma_ma_3*zs(j,i) &
462 + c_ma_3*abs(phi(j,i))*pi_180_inv
468 temp_mj_present(j,i) = theta_mj + gamma_mj*zs(j,i) &
469 + c_mj*abs(phi(j,i))*pi_180_inv
475 temp_ma(j,i) = temp_ma_present(j,i) + delta_ts
476 temp_mm(j,i,7) = temp_mj_present(j,i) + delta_ts
478 #elif (TSURFACE == 5)
482 temp_ma(j,i) = temp_ma_present(j,i) + glac_index*temp_ma_lgm_anom(j,i)
483 temp_mm(j,i,7) = temp_mj_present(j,i) + glac_index*temp_mj_lgm_anom(j,i)
485 #elif (TSURFACE == 6)
490 temp_ma(j,i) = temp_ma_present(j,i) + temp_ma_anom_fact*temp_ma_anom(j,i)
491 temp_mm(j,i,7) = temp_mj_present(j,i) + temp_mj_anom_fact*temp_mj_anom(j,i)
497 temp_ampl(j,i) = temp_mm(j,i,7) - temp_ma(j,i)
499 if (temp_ampl(j,i) < eps)
then
510 sine_factor = sin((
real(n,dp)-4.0_dp)*pi/6.0_dp)
514 temp_mm(j,i,n) = temp_ma(j,i) + sine_factor*temp_ampl(j,i)
522 #if (ACCSURFACE <= 3)
524 #if (ELEV_DESERT == 1)
526 gamma_p = gamma_p*1.0e-03_dp
528 zs_thresh = zs_thresh
532 #elif (ACCSURFACE == 4)
542 #if (SOLID_PRECIP == 1) /* Marsiat (1994) */
549 inv_delta_temp_rain_snow = 1.0_dp/(temp_rain-temp_snow)
551 #elif (SOLID_PRECIP == 2) /* Bales et al. (2009) */
558 coeff(0) = 5.4714e-01_dp
559 coeff(1) = -9.1603e-02_dp
560 coeff(2) = -3.314e-03_dp
561 coeff(3) = 4.66e-04_dp
562 coeff(4) = 3.8e-05_dp
563 coeff(5) = 6.0e-07_dp
565 #elif (SOLID_PRECIP == 3) /* Huybrechts and de Wolde (1999) */
569 temp_snow = temp_rain
575 inv_sqrt2_s_stat = 1.0_dp/(sqrt(2.0_dp)*s_stat)
579 #if (ABLSURFACE==1 || ABLSURFACE==2)
582 beta1 = beta1_0 *(0.001_dp/86400.0_dp)*(rho_w/rho)
584 beta2 = beta2_0 *(0.001_dp/86400.0_dp)*(rho_w/rho)
587 mu = mu_0 *(1000.0_dp*86400.0_dp)*(rho/rho_w)
590 #elif (ABLSURFACE==3)
592 lambda_lti = lambda_lti *(0.001_dp/year_sec)*(rho_w/rho)
604 #if (ACCSURFACE <= 3)
608 #if (ELEV_DESERT == 0)
612 #elif (ELEV_DESERT == 1)
614 if (zs_ref(j,i) < zs_thresh)
then
616 = exp(gamma_p*(max(zs(j,i),zs_thresh)-zs_thresh))
619 = exp(gamma_p*(max(zs(j,i),zs_thresh)-zs_ref(j,i)))
623 stop
' boundary: Parameter ELEV_DESERT must be either 0 or 1!'
627 precip(j,i,n) = precip_present(j,i,n)*precip_fact
635 precip_fact = accfact
637 precip_fact = 1.0_dp + gamma_s*delta_ts
639 precip_fact = exp(gamma_s*delta_ts)
642 #if (ACCSURFACE <= 3)
644 precip(j,i,0) = 0.0_dp
647 precip(j,i,n) = precip(j,i,n)*precip_fact
648 precip(j,i,0) = precip(j,i,0) + precip(j,i,n)*inv_twelve
652 #elif (ACCSURFACE == 4)
654 precip(j,i,0) = 0.0_dp
656 temp_inv = alpha_t * (temp_ma(j,i)+temp_0) + beta_t
657 temp_inv_present = alpha_t * (temp_ma_present(j,i)+temp_0) + beta_t
659 precip_fact = exp(alpha_p*(temp_0/temp_inv_present-temp_0/temp_inv)) &
660 *(temp_inv_present/temp_inv)**2 &
661 *(1.0_dp+beta_p*(temp_inv-temp_inv_present))
664 precip(j,i,n) = precip_present(j,i,n)*precip_fact
665 precip(j,i,0) = precip(j,i,0) + precip(j,i,n)*inv_twelve
669 #elif (ACCSURFACE == 5)
671 precip(j,i,0) = 0.0_dp
675 #if (PRECIP_ANOM_INTERPOL==1)
676 precip_fact = 1.0_dp-glac_index+glac_index*precip_lgm_anom(j,i,n)
678 #elif (PRECIP_ANOM_INTERPOL==2)
679 precip_fact = exp(-glac_index*gamma_precip_lgm_anom(j,i,n))
683 precip(j,i,n) = precip_present(j,i,n)*precip_fact
684 precip(j,i,0) = precip(j,i,0) + precip(j,i,n)*inv_twelve
688 #elif (ACCSURFACE == 6)
692 precip(j,i,0) = precip_ma_present(j,i) + precip_anom_fact*precip_ma_anom(j,i)
696 precip(j,i,n) = precip(j,i,0)
704 accum(j,i) = precip(j,i,0)
706 snowfall(j,i) = 0.0_dp
710 #if (SOLID_PRECIP == 1) /* Marsiat (1994) */
712 if (temp_mm(j,i,n) >= temp_rain)
then
714 else if (temp_mm(j,i,n) <= temp_snow)
then
717 frac_solid = (temp_rain-temp_mm(j,i,n))*inv_delta_temp_rain_snow
720 #elif (SOLID_PRECIP == 2) /* Bales et al. (2009) */
722 if (temp_mm(j,i,n) >= temp_rain)
then
724 else if (temp_mm(j,i,n) <= temp_snow)
then
727 frac_solid = coeff(0) + temp_mm(j,i,n) * ( coeff(1) &
728 + temp_mm(j,i,n) * ( coeff(2) &
729 + temp_mm(j,i,n) * ( coeff(3) &
730 + temp_mm(j,i,n) * ( coeff(4) &
731 + temp_mm(j,i,n) * coeff(5) ) ) ) )
735 #elif (SOLID_PRECIP == 3) /* Huybrechts and de Wolde (1999) */
737 frac_solid = 1.0_dp &
738 - 0.5_dp*
erfcc((temp_rain-temp_mm(j,i,n))*inv_sqrt2_s_stat)
742 snowfall(j,i) = snowfall(j,i) + precip(j,i,n)*frac_solid*inv_twelve
746 rainfall(j,i) = precip(j,i,0) - snowfall(j,i)
748 if (snowfall(j,i) < 0.0_dp) snowfall(j,i) = 0.0_dp
749 if (rainfall(j,i) < 0.0_dp) rainfall(j,i) = 0.0_dp
755 #if (ABLSURFACE==1 || ABLSURFACE==2)
760 temp_mm_help(n) = temp_mm(j,i,n)
763 call
pdd(temp_mm_help, s_stat, et(j,i))
770 if ((beta1*et(j,i)) <= (pmax*snowfall(j,i)))
then
771 melt_star(j,i) = beta1*et(j,i)
773 runoff(j,i) = melt(j,i)+rainfall(j,i)
775 melt_star(j,i) = pmax*snowfall(j,i)
776 melt(j,i) = beta2*(et(j,i)-melt_star(j,i)/beta1)
777 runoff(j,i) = melt(j,i)+rainfall(j,i)
780 #elif (ABLSURFACE==2)
782 if ( rainfall(j,i) <= (pmax*snowfall(j,i)) )
then
784 if ( (rainfall(j,i)+beta1*et(j,i)) <= (pmax*snowfall(j,i)) )
then
785 melt_star(j,i) = rainfall(j,i)+beta1*et(j,i)
787 runoff(j,i) = melt(j,i)
789 melt_star(j,i) = pmax*snowfall(j,i)
791 *(et(j,i)-(melt_star(j,i)-rainfall(j,i))/beta1)
792 runoff(j,i) = melt(j,i)
797 melt_star(j,i) = pmax*snowfall(j,i)
798 melt(j,i) = beta2*et(j,i)
799 runoff(j,i) = melt(j,i) + rainfall(j,i)-pmax*snowfall(j,i)
805 #elif (ABLSURFACE==3)
807 temp_jja_help = one_third*(temp_mm(j,i,6)+temp_mm(j,i,7)+temp_mm(j,i,8))
809 melt_star(j,i) = 0.0_dp
810 melt(j,i) = lambda_lti*max((temp_jja_help-temp_lti), 0.0_dp)
811 runoff(j,i) = melt(j,i) + rainfall(j,i)
821 as_perp(j,i) = accum(j,i) - evap(j,i) - runoff(j,i)
827 if (melt_star(j,i) >= melt(j,i))
then
828 temp_s(j,i) = temp_ma(j,i) &
829 +mu*(melt_star(j,i)-melt(j,i))
831 temp_s(j,i) = temp_ma(j,i)
834 if (temp_s(j,i) > -0.001_dp) temp_s(j,i) = -0.001_dp
integer(i2b) function mask_update(z_sl, i, j)
Update of the topography mask due to changes of the sea level.
Declarations of kind types for SICOPOLIS.
subroutine pdd(temp_mm, s_stat, ET)
Computation of the positive degree days (PDD) with statistical temperature fluctuations; based on sem...
subroutine check(status, ch_calling_routine)
NetCDF error capturing.
real(dp) function erfcc(x)
Computation of the complementary error function erfc(x) = 1-erf(x) with a fractional error everywhere...
Declarations of global variables for SICOPOLIS (for the ANT domain).
subroutine boundary(time, dtime, dxi, deta, delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
Computation of the surface temperature (must be less than 0 deg C!) and of the accumulation-ablation ...
Declarations of global variables for SICOPOLIS.