7 !! Computation of the surface temperature (must be less than 0 deg C!)
8 !! and of the accumulation-ablation function.
12 !! Copyright 2009-2013 Ralf Greve
16 !! This file is part of SICOPOLIS.
18 !! SICOPOLIS is free software: you can redistribute it and/or modify
19 !! it under the terms of the GNU General Public License as published by
20 !! the Free Software Foundation, either version 3 of the License, or
21 !! (at your option) any later version.
23 !! SICOPOLIS is distributed in the hope that it will be useful,
24 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
25 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26 !! GNU General Public License for more details.
28 !! You should have received a copy of the GNU General Public License
29 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
31 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
34 !> Computation of the surface temperature (must be less than 0 deg C!)
35 !! and of the accumulation-ablation function.
36 !<------------------------------------------------------------------------------
38 delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
45 real(dp),
intent(in) :: time, dtime, dxi, deta
47 real(dp),
intent(out) :: delta_ts, glac_index, dzsl_dtau, z_mar
48 real(dp),
intent(inout) :: z_sl
57 integer(i4b) :: i_gr, i_kl
59 real(dp) :: z_sl_min, t1, t2, t3, t4, t5, t6
60 real(dp) :: time_gr, time_kl
61 real(dp),
dimension(0:JMAX,0:IMAX) :: dist
62 logical,
dimension(0:JMAX,0:IMAX) :: check_point
81 delta_ts = sine_amplit &
82 *cos(2.0_dp*pi*time/(sine_period*year_sec)) &
89 if (time/year_sec.lt.
real(grip_time_min,dp)) then
90 delta_ts = griptemp(0)
91 else if (time/year_sec.lt.
real(grip_time_max,dp)) then
93 i_kl = floor(((time/year_sec) &
94 -
real(grip_time_min,dp))/
real(grip_time_stp,dp))
97 i_gr = ceiling(((time/year_sec) &
98 -
real(grip_time_min,dp))/
real(grip_time_stp,dp))
99 i_gr = min(i_gr, ndata_grip)
101 if (i_kl.eq.i_gr)
then
103 delta_ts = griptemp(i_kl)
107 time_kl = (grip_time_min + i_kl*grip_time_stp) *year_sec
108 time_gr = (grip_time_min + i_gr*grip_time_stp) *year_sec
110 delta_ts = griptemp(i_kl) &
111 +(griptemp(i_gr)-griptemp(i_kl)) &
112 *(time-time_kl)/(time_gr-time_kl)
118 delta_ts = griptemp(ndata_grip)
121 delta_ts = delta_ts * grip_temp_fact
140 t1 = -250000.0_dp *year_sec
141 t2 = -140000.0_dp *year_sec
142 t3 = -125000.0_dp *year_sec
143 t4 = -21000.0_dp *year_sec
144 t5 = -8000.0_dp *year_sec
145 t6 = 0.0_dp *year_sec
149 else if (time.lt.t2)
then
150 z_sl = z_sl_min*(time-t1)/(t2-t1)
151 else if (time.lt.t3)
then
152 z_sl = -z_sl_min*(time-t3)/(t3-t2)
153 else if (time.lt.t4)
then
154 z_sl = z_sl_min*(time-t3)/(t4-t3)
155 else if (time.lt.t5)
then
156 z_sl = -z_sl_min*(time-t5)/(t5-t4)
157 else if (time.lt.t6)
then
167 if (time/year_sec.lt.
real(specmap_time_min,dp)) then
168 z_sl = specmap_zsl(0)
169 else if (time/year_sec.lt.
real(specmap_time_max,dp)) then
171 i_kl = floor(((time/year_sec) &
172 -
real(specmap_time_min,dp))/
real(specmap_time_stp,dp))
175 i_gr = ceiling(((time/year_sec) &
176 -
real(specmap_time_min,dp))/
real(specmap_time_stp,dp))
177 i_gr = min(i_gr, ndata_specmap)
179 if (i_kl.eq.i_gr)
then
181 z_sl = specmap_zsl(i_kl)
185 time_kl = (specmap_time_min + i_kl*specmap_time_stp) *year_sec
186 time_gr = (specmap_time_min + i_gr*specmap_time_stp) *year_sec
188 z_sl = specmap_zsl(i_kl) &
189 +(specmap_zsl(i_gr)-specmap_zsl(i_kl)) &
190 *(time-time_kl)/(time_gr-time_kl)
196 z_sl = specmap_zsl(ndata_specmap)
203 if ( z_sl_old > -999999.9_dp )
then
204 dzsl_dtau = (z_sl-z_sl_old)/dtime
213 #if ( MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3 )
215 #elif ( MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 )
216 z_mar = fact_z_mar*z_sl
217 #elif ( MARINE_ICE_CALVING==6 || MARINE_ICE_CALVING==7 )
218 if (z_sl >= -80.0_dp)
then
221 z_mar = 10.25_dp*(z_sl+80.0_dp)-200.0_dp
223 z_mar = fact_z_mar*z_mar
235 check_point(j,i) = .false.
241 if (maske(j,i).ge.2)
then
242 check_point(j ,i ) = .true.
243 check_point(j ,i+1) = .true.
244 check_point(j ,i-1) = .true.
245 check_point(j+1,i ) = .true.
246 check_point(j-1,i ) = .true.
253 if (check_point(j,i))
then
263 if (check_point(j,i))
then
264 maske(j,i) = maske_neu(j,i)
273 dist(j,i) = sqrt( (xi(i)-x_hat)**2 + (eta(j)-y_hat)**2 )
280 temp_s(j,i) = temp_min + s_t*dist(j,i)
281 temp_s(j,i) = temp_s(j,i) + delta_ts
283 if (temp_s(j,i).gt. -0.001_dp) temp_s(j,i) = -0.001_dp
286 accum_present(j,i) = b_max
287 accum(j,i) = accum_present(j,i)
289 as_perp(j,i) = s_b*(eld-dist(j,i))
290 if (as_perp(j,i).gt.accum(j,i)) as_perp(j,i) = accum(j,i)