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 real(dp) :: rad_inv
85 logical, dimension(0:JMAX,0:IMAX) :: check_point
86 logical, save :: firstcall = .true.
87 
88 !-------- Initialization of intent(out) variables --------
89 
90 z_sl_old = z_sl
91 
92 delta_ts = 0.0_dp
93 glac_index = 0.0_dp
94 z_sl = 0.0_dp
95 dzsl_dtau = 0.0_dp
96 z_mar = 0.0_dp
97 
98 !-------- Surface-temperature deviation from present values --------
99 
100 #if (TSURFACE==1)
101 delta_ts = delta_ts0
102 ! ! Steady state with prescribed constant
103 ! ! air-temperature deviation
104 #elif (TSURFACE==3)
105 delta_ts = sine_amplit &
106  *cos(2.0_dp*pi*time/(sine_period*year_sec)) &
107  -sine_amplit
108 ! ! Sinusoidal air-temperature forcing
109 #elif (TSURFACE==4)
110 
111 ! ------ delta_ts from the GRIP record
112 
113 if (time/year_sec.lt.real(grip_time_min,dp)) then
114  delta_ts = griptemp(0)
115 else if (time/year_sec.lt.real(grip_time_max,dp)) then
116 
117  i_kl = floor(((time/year_sec) &
118  -real(grip_time_min,dp))/real(grip_time_stp,dp))
119  i_kl = max(i_kl, 0)
120 
121  i_gr = ceiling(((time/year_sec) &
122  -real(grip_time_min,dp))/real(grip_time_stp,dp))
123  i_gr = min(i_gr, ndata_grip)
124 
125  if (i_kl.eq.i_gr) then
126 
127  delta_ts = griptemp(i_kl)
128 
129  else
130 
131  time_kl = (grip_time_min + i_kl*grip_time_stp) *year_sec
132  time_gr = (grip_time_min + i_gr*grip_time_stp) *year_sec
133 
134  delta_ts = griptemp(i_kl) &
135  +(griptemp(i_gr)-griptemp(i_kl)) &
136  *(time-time_kl)/(time_gr-time_kl)
137  ! linear interpolation of the ice-core data
138 
139  end if
140 
141 else
142  delta_ts = griptemp(ndata_grip)
143 end if
144 
145 delta_ts = delta_ts * grip_temp_fact
146 ! ! modification by constant factor
147 
148 #elif (TSURFACE==5)
149 ! Enter here delta_ts scenario:
150 
151 #endif
152 
153 !-------- Sea level --------
154 
155 #if (SEA_LEVEL==1)
156 ! ------ constant sea level
157 z_sl = z_sl0
158 
159 #elif (SEA_LEVEL==2)
160 ! ------ saw-tooth-shaped palaeoclimatic sea-level forcing
161 
162 z_sl_min = -130.0_dp
163 
164 t1 = -250000.0_dp *year_sec
165 t2 = -140000.0_dp *year_sec
166 t3 = -125000.0_dp *year_sec
167 t4 = -21000.0_dp *year_sec
168 t5 = -8000.0_dp *year_sec
169 t6 = 0.0_dp *year_sec
170 
171 if (time.lt.t1) then
172  z_sl = 0.0_dp
173 else if (time.lt.t2) then
174  z_sl = z_sl_min*(time-t1)/(t2-t1)
175 else if (time.lt.t3) then
176  z_sl = -z_sl_min*(time-t3)/(t3-t2)
177 else if (time.lt.t4) then
178  z_sl = z_sl_min*(time-t3)/(t4-t3)
179 else if (time.lt.t5) then
180  z_sl = -z_sl_min*(time-t5)/(t5-t4)
181 else if (time.lt.t6) then
182  z_sl = 0.0_dp
183 else
184  z_sl = 0.0_dp
185 end if
186 
187 #elif (SEA_LEVEL==3)
188 
189 ! ------ z_sl from the SPECMAP record
190 
191 if (time/year_sec.lt.real(specmap_time_min,dp)) then
192  z_sl = specmap_zsl(0)
193 else if (time/year_sec.lt.real(specmap_time_max,dp)) then
194 
195  i_kl = floor(((time/year_sec) &
196  -real(specmap_time_min,dp))/real(specmap_time_stp,dp))
197  i_kl = max(i_kl, 0)
198 
199  i_gr = ceiling(((time/year_sec) &
200  -real(specmap_time_min,dp))/real(specmap_time_stp,dp))
201  i_gr = min(i_gr, ndata_specmap)
202 
203  if (i_kl.eq.i_gr) then
204 
205  z_sl = specmap_zsl(i_kl)
206 
207  else
208 
209  time_kl = (specmap_time_min + i_kl*specmap_time_stp) *year_sec
210  time_gr = (specmap_time_min + i_gr*specmap_time_stp) *year_sec
211 
212  z_sl = specmap_zsl(i_kl) &
213  +(specmap_zsl(i_gr)-specmap_zsl(i_kl)) &
214  *(time-time_kl)/(time_gr-time_kl)
215  ! linear interpolation of the sea-level data
216 
217  end if
218 
219 else
220  z_sl = specmap_zsl(ndata_specmap)
221 end if
222 
223 #endif
224 
225 ! ------ Time derivative of the sea level
226 
227 if ( z_sl_old > -999999.9_dp ) then
228  dzsl_dtau = (z_sl-z_sl_old)/dtime
229 else ! only dummy value for z_sl_old available
230  dzsl_dtau = 0.0_dp
231 end if
232 
233 ! ------ Minimum bedrock elevation for extent of marine ice
234 
235 #if (MARGIN==2)
236 
237 #if (MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3)
238 z_mar = z_mar
239 #elif (MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5)
240 z_mar = fact_z_mar*z_sl
241 #elif (MARINE_ICE_CALVING==6 || MARINE_ICE_CALVING==7)
242 if (z_sl >= -80.0_dp) then
243  z_mar = 2.5_dp*z_sl
244 else
245  z_mar = 10.25_dp*(z_sl+80.0_dp)-200.0_dp
246 end if
247 z_mar = fact_z_mar*z_mar
248 #endif
249 
250 #endif
251 
252 ! ------ Update of the mask according to the sea level
253 
254 ! ---- Check all sea and floating-ice points and their direct
255 ! neighbours
256 
257 do i=0, imax
258 do j=0, jmax
259  check_point(j,i) = .false.
260 end do
261 end do
262 
263 do i=1, imax-1
264 do j=1, jmax-1
265  if (maske(j,i).ge.2) then
266  check_point(j ,i ) = .true.
267  check_point(j ,i+1) = .true.
268  check_point(j ,i-1) = .true.
269  check_point(j+1,i ) = .true.
270  check_point(j-1,i ) = .true.
271  end if
272 end do
273 end do
274 
275 do i=1, imax-1
276 do j=1, jmax-1
277  if (check_point(j,i)) then
278  maske_neu(j,i) = mask_update_sea_level(z_sl, i, j)
279  end if
280 end do
281 end do
282 
283 ! ---- Assign new values of the mask
284 
285 do i=1, imax-1
286 do j=1, jmax-1
287  if (check_point(j,i)) then
288  maske(j,i) = maske_neu(j,i)
289  end if
290 end do
291 end do
292 
293 !-------- Surface air temperature and accumulation-ablation function --------
294 
295 do i=0, imax
296 do j=0, jmax
297  dist(j,i) = sqrt( (xi(i)-x_hat)**2 + (eta(j)-y_hat)**2 )
298 end do
299 end do
300 
301 rad_inv = 1.0_dp/rad
302 
303 do i=0, imax
304 do j=0, jmax
305 
306  temp_s(j,i) = temp_min + s_t*dist(j,i)**3
307  temp_s(j,i) = temp_s(j,i) + delta_ts
308  ! Correction with temperature deviation delta_ts
309  if (temp_s(j,i) > -0.001_dp) temp_s(j,i) = -0.001_dp
310  ! Cut-off of positive air temperatures
311 
312  accum_present(j,i) = b_min + (b_max-b_min)*rad_inv*dist(j,i)
313  accum(j,i) = accum_present(j,i)
314 
315  as_perp(j,i) = accum(j,i)
316 
317 end do
318 end do
319 
320 runoff = accum - as_perp
321 evap = 0.0_dp
322 
323 !-------- Calving rate of grounded ice --------
324 
325 calv_grounded = 0.0_dp
326 
327 #if ((MARGIN==2) \
328  && (marine_ice_formation==2) \
329  && (marine_ice_calving==9))
330 
331 call calving_underwater_ice(z_sl)
333 
334 #endif
335 
336 if (firstcall) firstcall = .false.
337 
338 end subroutine boundary
339 
340 !-------------------------------------------------------------------------------
341 
342 end module boundary_m
343 !
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), 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".
real(dp) b_min
b_min: Minimum accumulation rate
Definition: sico_vars_m.F90:58
Declarations of kind types for SICOPOLIS.
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), 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) rad
rad: Radius of the model domain
Definition: sico_vars_m.F90:56
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