SICOPOLIS V5-dev  Revision 1173
boundary_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : b o u n d a r y _ m
4 !
5 !> @file
6 !!
7 !! Computation of the surface temperature (must be less than 0 deg C!)
8 !! and of the accumulation-ablation function.
9 !!
10 !! @section Copyright
11 !!
12 !! Copyright 2009-2017 Ralf Greve
13 !!
14 !! @section License
15 !!
16 !! This file is part of SICOPOLIS.
17 !!
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.
22 !!
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.
27 !!
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/>.
30 !<
31 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
32 
33 !-------------------------------------------------------------------------------
34 !> Computation of the surface temperature (must be less than 0 deg C!)
35 !! and of the accumulation-ablation function.
36 !<------------------------------------------------------------------------------
37 module boundary_m
38 
39  use sico_types_m
41  use sico_vars_m
42 
43  implicit none
44 
45  public
46 
47 contains
48 
49 !-------------------------------------------------------------------------------
50 !> Main routine of boundary_m:
51 !! Computation of the surface temperature (must be less than 0 deg C!)
52 !! and of the accumulation-ablation function.
53 !<------------------------------------------------------------------------------
54 subroutine boundary(time, dtime, dxi, deta, &
55  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
56 
57 #if ((MARGIN==2) \
58  && (marine_ice_formation==2) \
59  && (marine_ice_calving==9))
61 #endif
62 
64 
65 implicit none
66 
67 real(dp), intent(in) :: time, dtime, dxi, deta
68 
69 real(dp), intent(out) :: delta_ts, glac_index, dzsl_dtau, z_mar
70 real(dp), intent(inout) :: z_sl
71 
72 ! Further return variables
73 ! (defined as global variables in module sico_variables_m):
74 !
75 ! accum(j,i), evap(j,i), runoff(j,i), as_perp(j,i),
76 ! calv_grounded(j,i), temp_s(j,i)
77 
78 integer(i4b) :: i, j
79 integer(i4b) :: i_gr, i_kl
80 real(dp) :: z_sl_old
81 real(dp) :: z_sl_min, t1, t2, t3, t4, t5, t6
82 real(dp) :: time_gr, time_kl
83 real(dp), dimension(0:JMAX,0:IMAX) :: dist
84 logical, dimension(0:JMAX,0:IMAX) :: check_point
85 logical, save :: firstcall = .true.
86 
87 !-------- Initialization of intent(out) variables --------
88 
89 z_sl_old = z_sl
90 
91 delta_ts = 0.0_dp
92 glac_index = 0.0_dp
93 z_sl = 0.0_dp
94 dzsl_dtau = 0.0_dp
95 z_mar = 0.0_dp
96 
97 !-------- Surface-temperature deviation from present values --------
98 
99 #if (TSURFACE==1)
100 delta_ts = delta_ts0
101 ! ! Steady state with prescribed constant
102 ! ! air-temperature deviation
103 #elif (TSURFACE==3)
104 delta_ts = sine_amplit &
105  *cos(2.0_dp*pi*time/(sine_period*year_sec)) &
106  -sine_amplit
107 ! ! Sinusoidal air-temperature forcing
108 #elif (TSURFACE==4)
109 
110 ! ------ delta_ts from the GRIP record
111 
112 if (time/year_sec.lt.real(grip_time_min,dp)) then
113  delta_ts = griptemp(0)
114 else if (time/year_sec.lt.real(grip_time_max,dp)) then
115 
116  i_kl = floor(((time/year_sec) &
117  -real(grip_time_min,dp))/real(grip_time_stp,dp))
118  i_kl = max(i_kl, 0)
119 
120  i_gr = ceiling(((time/year_sec) &
121  -real(grip_time_min,dp))/real(grip_time_stp,dp))
122  i_gr = min(i_gr, ndata_grip)
123 
124  if (i_kl.eq.i_gr) then
125 
126  delta_ts = griptemp(i_kl)
127 
128  else
129 
130  time_kl = (grip_time_min + i_kl*grip_time_stp) *year_sec
131  time_gr = (grip_time_min + i_gr*grip_time_stp) *year_sec
132 
133  delta_ts = griptemp(i_kl) &
134  +(griptemp(i_gr)-griptemp(i_kl)) &
135  *(time-time_kl)/(time_gr-time_kl)
136  ! linear interpolation of the ice-core data
137 
138  end if
139 
140 else
141  delta_ts = griptemp(ndata_grip)
142 end if
143 
144 delta_ts = delta_ts * grip_temp_fact
145 ! ! modification by constant factor
146 
147 #elif (TSURFACE==5)
148 ! Enter here delta_ts scenario:
149 
150 #endif
151 
152 !-------- Sea level --------
153 
154 #if (SEA_LEVEL==1)
155 ! ------ constant sea level
156 z_sl = z_sl0
157 
158 #elif (SEA_LEVEL==2)
159 ! ------ saw-tooth-shaped palaeoclimatic sea-level forcing
160 
161 z_sl_min = -130.0_dp
162 
163 t1 = -250000.0_dp *year_sec
164 t2 = -140000.0_dp *year_sec
165 t3 = -125000.0_dp *year_sec
166 t4 = -21000.0_dp *year_sec
167 t5 = -8000.0_dp *year_sec
168 t6 = 0.0_dp *year_sec
169 
170 if (time.lt.t1) then
171  z_sl = 0.0_dp
172 else if (time.lt.t2) then
173  z_sl = z_sl_min*(time-t1)/(t2-t1)
174 else if (time.lt.t3) then
175  z_sl = -z_sl_min*(time-t3)/(t3-t2)
176 else if (time.lt.t4) then
177  z_sl = z_sl_min*(time-t3)/(t4-t3)
178 else if (time.lt.t5) then
179  z_sl = -z_sl_min*(time-t5)/(t5-t4)
180 else if (time.lt.t6) then
181  z_sl = 0.0_dp
182 else
183  z_sl = 0.0_dp
184 end if
185 
186 #elif (SEA_LEVEL==3)
187 
188 ! ------ z_sl from the SPECMAP record
189 
190 if (time/year_sec.lt.real(specmap_time_min,dp)) then
191  z_sl = specmap_zsl(0)
192 else if (time/year_sec.lt.real(specmap_time_max,dp)) then
193 
194  i_kl = floor(((time/year_sec) &
195  -real(specmap_time_min,dp))/real(specmap_time_stp,dp))
196  i_kl = max(i_kl, 0)
197 
198  i_gr = ceiling(((time/year_sec) &
199  -real(specmap_time_min,dp))/real(specmap_time_stp,dp))
200  i_gr = min(i_gr, ndata_specmap)
201 
202  if (i_kl.eq.i_gr) then
203 
204  z_sl = specmap_zsl(i_kl)
205 
206  else
207 
208  time_kl = (specmap_time_min + i_kl*specmap_time_stp) *year_sec
209  time_gr = (specmap_time_min + i_gr*specmap_time_stp) *year_sec
210 
211  z_sl = specmap_zsl(i_kl) &
212  +(specmap_zsl(i_gr)-specmap_zsl(i_kl)) &
213  *(time-time_kl)/(time_gr-time_kl)
214  ! linear interpolation of the sea-level data
215 
216  end if
217 
218 else
219  z_sl = specmap_zsl(ndata_specmap)
220 end if
221 
222 #endif
223 
224 ! ------ Time derivative of the sea level
225 
226 if ( z_sl_old > -999999.9_dp ) then
227  dzsl_dtau = (z_sl-z_sl_old)/dtime
228 else ! only dummy value for z_sl_old available
229  dzsl_dtau = 0.0_dp
230 end if
231 
232 ! ------ Minimum bedrock elevation for extent of marine ice
233 
234 #if (MARGIN==2)
235 
236 #if (MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3)
237 z_mar = z_mar
238 #elif (MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5)
239 z_mar = fact_z_mar*z_sl
240 #elif (MARINE_ICE_CALVING==6 || MARINE_ICE_CALVING==7)
241 if (z_sl >= -80.0_dp) then
242  z_mar = 2.5_dp*z_sl
243 else
244  z_mar = 10.25_dp*(z_sl+80.0_dp)-200.0_dp
245 end if
246 z_mar = fact_z_mar*z_mar
247 #endif
248 
249 #endif
250 
251 ! ------ Update of the mask according to the sea level
252 
253 ! ---- Check all sea and floating-ice points and their direct
254 ! neighbours
255 
256 do i=0, imax
257 do j=0, jmax
258  check_point(j,i) = .false.
259 end do
260 end do
261 
262 do i=1, imax-1
263 do j=1, jmax-1
264  if (maske(j,i).ge.2) then
265  check_point(j ,i ) = .true.
266  check_point(j ,i+1) = .true.
267  check_point(j ,i-1) = .true.
268  check_point(j+1,i ) = .true.
269  check_point(j-1,i ) = .true.
270  end if
271 end do
272 end do
273 
274 do i=1, imax-1
275 do j=1, jmax-1
276  if (check_point(j,i)) then
277  maske_neu(j,i) = mask_update_sea_level(z_sl, i, j)
278  end if
279 end do
280 end do
281 
282 ! ---- Assign new values of the mask
283 
284 do i=1, imax-1
285 do j=1, jmax-1
286  if (check_point(j,i)) then
287  maske(j,i) = maske_neu(j,i)
288  end if
289 end do
290 end do
291 
292 !-------- Surface air temperature and accumulation-ablation function --------
293 
294 #if (!defined(SURFACE_FORCING) || SURFACE_FORCING==1)
295 
296 do i=0, imax
297 do j=0, jmax
298  dist(j,i) = sqrt( (xi(i)-x_hat)**2 + (eta(j)-y_hat)**2 )
299 end do
300 end do
301 
302 do i=0, imax
303 do j=0, jmax
304 
305  temp_s(j,i) = temp_min + s_t*dist(j,i)
306  temp_s(j,i) = temp_s(j,i) + delta_ts
307  ! Correction with temperature deviation delta_ts
308  if (temp_s(j,i) > -0.001_dp) temp_s(j,i) = -0.001_dp
309  ! Cut-off of positive air temperatures
310 
311  accum_present(j,i) = b_max
312  accum(j,i) = accum_present(j,i)
313 
314  as_perp(j,i) = s_b*(eld-dist(j,i))
315  if (as_perp(j,i) > accum(j,i)) as_perp(j,i) = accum(j,i)
316 
317 end do
318 end do
319 
320 #elif (SURFACE_FORCING==2)
321 
322 do i=0, imax
323 do j=0, jmax
324 
325  temp_s(j,i) = temp_0 - gamma_t*zs(j,i)
326  temp_s(j,i) = temp_s(j,i) + delta_ts
327  ! Correction with temperature deviation delta_ts
328  if (temp_s(j,i) > -0.001_dp) temp_s(j,i) = -0.001_dp
329  ! Cut-off of positive air temperatures
330 
331  accum_present(j,i) = s_0
332  accum(j,i) = accum_present(j,i)
333 
334  as_perp(j,i) = m_0*(zs(j,i)-ela)
335  if (as_perp(j,i) > accum(j,i)) as_perp(j,i) = accum(j,i)
336 
337 end do
338 end do
339 
340 #endif
341 
342 runoff = accum - as_perp
343 evap = 0.0_dp
344 
345 !-------- Calving rate of grounded ice --------
346 
347 calv_grounded = 0.0_dp
348 
349 #if ((MARGIN==2) \
350  && (marine_ice_formation==2) \
351  && (marine_ice_calving==9))
352 
353 call calving_underwater_ice(z_sl)
355 
356 #endif
357 
358 if (firstcall) firstcall = .false.
359 
360 end subroutine boundary
361 
362 !-------------------------------------------------------------------------------
363 
364 end module boundary_m
365 !
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
integer(i4b) specmap_time_max
specmap_time_max: Maximum time of the data values for the sea level
real(dp) s_b
s_b: Gradient of accumulation rate change with horizontal distance
Definition: sico_vars_m.F90:54
real(dp), dimension(0:jmax, 0:imax) calv_grounded
calv_grounded(j,i): Calving rate of grounded ice
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
real(dp), dimension(0:jmax, 0:imax) as_perp
as_perp(j,i): Accumulation-ablation function at the ice surface (SMB)
real(dp), dimension(0:jmax, 0:imax) runoff
runoff(j,i): Runoff rate at the ice surface
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), 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 ...
real(dp) b_max
b_max: Maximum accumulation rate
Definition: sico_vars_m.F90:52
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars_m.F90:35
integer(i4b) ndata_grip
ndata_grip: Number of data values for the surface temperature anomaly
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), 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(:), allocatable specmap_zsl
specmap_zsl(n): Data values for the sea level
real(dp) eld
eld: Equilibrium line distance
Definition: sico_vars_m.F90:56
real(dp), dimension(0:jmax, 0:imax) accum_present
accum_present(j,i): Present-day accumulation rate at the ice surface (for EISMINT, ISMIP HEINO and the north and south polar caps of Mars)
integer(i2b), dimension(0:jmax, 0:imax) maske_neu
maske_neu(j,i): New value of maske computed during an integration step
real(dp) s_t
s_t: Gradient of surface temperature change with horizontal distance
Definition: sico_vars_m.F90:46
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...
Definition: boundary_m.F90:56
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 ...
real(dp), dimension(0:jmax) eta
eta(j): Coordinate eta (= y) of grid point j
real(dp) temp_min
temp_min: Minimum surface temperature
Definition: sico_vars_m.F90:44
real(dp), dimension(0:imax) xi
xi(i): Coordinate xi (= x) of grid point i
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 ...
Definition: boundary_m.F90:37
real(dp) y_hat
y_hat: Coordinate eta (= y) of the centre of the model domain
Definition: sico_vars_m.F90:50
real(dp) x_hat
x_hat: Coordinate xi (= x) of the centre of the model domain
Definition: sico_vars_m.F90:48
Update of the ice-land-ocean mask due to changes of the sea level.
integer, parameter dp
Double-precision reals.
Declarations of global variables for SICOPOLIS.
real(dp), dimension(:), allocatable griptemp
griptemp(n): Data values for the surface temperature anomaly