SICOPOLIS V3.0
 All Classes Files Functions Variables Macros
boundary.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : b o u n d a r y
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-2013 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 subroutine boundary(time, dtime, dxi, deta, &
38  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
39 
40 use sico_types
41 use sico_variables
42 
43 implicit none
44 
45 real(dp), intent(in) :: time, dtime, dxi, deta
46 
47 real(dp), intent(out) :: delta_ts, glac_index, dzsl_dtau, z_mar
48 real(dp), intent(inout) :: z_sl
49 
50 ! Further return variables
51 ! (defined as global variables in module sico_variables):
52 !
53 ! accum(j,i), evap(j,i), runoff(j,i), as_perp(j,i), temp_s(j,i)
54 
55 integer(i2b) :: mask_update
56 integer(i4b) :: i, j, n
57 integer(i4b) :: i_gr, i_kl
58 real(dp) :: z_sl_old
59 real(dp) :: z_sl_min, t1, t2, t3, t4, t5, t6
60 real(dp) :: time_gr, time_kl
61 real(dp) :: z_sle_present, z_sle_help
62 real(dp), dimension(0:JMAX,0:IMAX,0:12) :: precip
63 real(dp), dimension(0:JMAX,0:IMAX) :: &
64  snowfall, rainfall, melt, melt_star
65 real(dp), dimension(0:JMAX,0:IMAX,12) :: temp_mm
66 real(dp), dimension(0:JMAX,0:IMAX) :: temp_ma
67 real(dp), dimension(12) :: temp_mm_help
68 real(dp) :: temp_jja_help
69 real(dp), dimension(0:JMAX,0:IMAX) :: et
70 real(dp) :: gamma_t, temp_diff
71 real(dp) :: gamma_p, zs_thresh, &
72  temp_rain, temp_snow, &
73  inv_delta_temp_rain_snow, coeff(0:5), inv_sqrt2_s_stat, &
74  precip_fact, frac_solid
75 real(dp) :: s_stat, beta1, beta2, pmax, mu, lambda_lti, temp_lti
76 real(dp) :: erfcc
77 logical, dimension(0:JMAX,0:IMAX) :: check_point
78 
79 real(dp), parameter :: &
80  inv_twelve = 1.0_dp/12.0_dp, one_third = 1.0_dp/3.0_dp
81 
82 !-------- Initialization of intent(out) variables --------
83 
84 z_sl_old = z_sl
85 
86 delta_ts = 0.0_dp
87 glac_index = 0.0_dp
88 z_sl = 0.0_dp
89 dzsl_dtau = 0.0_dp
90 z_mar = 0.0_dp
91 
92 !-------- Surface-temperature deviation from present values --------
93 
94 #if TSURFACE==1
95 delta_ts = delta_ts0
96 ! ! Steady state with prescribed constant
97 ! ! air-temperature deviation
98 #elif TSURFACE==3
99 delta_ts = sine_amplit &
100  *cos(2.0_dp*pi*time/(sine_period*year_sec)) &
101  -sine_amplit
102 ! ! Sinusoidal air-temperature forcing
103 #elif TSURFACE==4
104 
105 ! ------ delta_ts from ice-core record
106 
107 if (time/year_sec.lt.real(grip_time_min,dp)) then
108  delta_ts = griptemp(0)
109 else if (time/year_sec.lt.real(grip_time_max,dp)) then
110 
111  i_kl = floor(((time/year_sec) &
112  -real(grip_time_min,dp))/real(grip_time_stp,dp))
113  i_kl = max(i_kl, 0)
114 
115  i_gr = ceiling(((time/year_sec) &
116  -real(grip_time_min,dp))/real(grip_time_stp,dp))
117  i_gr = min(i_gr, ndata_grip)
118 
119  if (i_kl.eq.i_gr) then
120 
121  delta_ts = griptemp(i_kl)
122 
123  else
124 
125  time_kl = (grip_time_min + i_kl*grip_time_stp) *year_sec
126  time_gr = (grip_time_min + i_gr*grip_time_stp) *year_sec
127 
128  delta_ts = griptemp(i_kl) &
129  +(griptemp(i_gr)-griptemp(i_kl)) &
130  *(time-time_kl)/(time_gr-time_kl)
131  ! linear interpolation of the ice-core data
132 
133  end if
134 
135 else
136  delta_ts = griptemp(ndata_grip)
137 end if
138 
139 delta_ts = delta_ts * grip_temp_fact
140 ! ! modification by constant factor
141 
142 !-------- Glacial index --------
143 
144 #elif TSURFACE==5
145 
146 if (time/year_sec < real(gi_time_min,dp)) then
147  glac_index = glacial_index(0)
148 else if (time/year_sec < real(gi_time_max,dp)) then
149 
150  i_kl = floor(((time/year_sec) &
151  -real(gi_time_min,dp))/real(gi_time_stp,dp))
152  i_kl = max(i_kl, 0)
153 
154  i_gr = ceiling(((time/year_sec) &
155  -real(gi_time_min,dp))/real(gi_time_stp,dp))
156  i_gr = min(i_gr, ndata_gi)
157 
158  if (i_kl == i_gr) then
159 
160  glac_index = glacial_index(i_kl)
161 
162  else
163 
164  time_kl = (gi_time_min + i_kl*gi_time_stp) *year_sec
165  time_gr = (gi_time_min + i_gr*gi_time_stp) *year_sec
166 
167  glac_index = glacial_index(i_kl) &
168  +(glacial_index(i_gr)-glacial_index(i_kl)) &
169  *(time-time_kl)/(time_gr-time_kl)
170  ! linear interpolation of the glacial-index data
171 
172  end if
173 
174 else
175  glac_index = glacial_index(ndata_gi)
176 end if
177 
178 #endif
179 
180 !-------- Sea level --------
181 
182 #if SEA_LEVEL==1
183 ! ------ constant sea level
184 z_sl = z_sl0
185 
186 #elif SEA_LEVEL==2
187 ! ------ saw-tooth-shaped palaeoclimatic sea-level forcing
188 
189 z_sl_min = -130.0_dp
190 
191 t1 = -250000.0_dp *year_sec
192 t2 = -140000.0_dp *year_sec
193 t3 = -125000.0_dp *year_sec
194 t4 = -21000.0_dp *year_sec
195 t5 = -8000.0_dp *year_sec
196 t6 = 0.0_dp *year_sec
197 
198 if (time.lt.t1) then
199  z_sl = 0.0_dp
200 else if (time.lt.t2) then
201  z_sl = z_sl_min*(time-t1)/(t2-t1)
202 else if (time.lt.t3) then
203  z_sl = -z_sl_min*(time-t3)/(t3-t2)
204 else if (time.lt.t4) then
205  z_sl = z_sl_min*(time-t3)/(t4-t3)
206 else if (time.lt.t5) then
207  z_sl = -z_sl_min*(time-t5)/(t5-t4)
208 else if (time.lt.t6) then
209  z_sl = 0.0_dp
210 else
211  z_sl = 0.0_dp
212 end if
213 
214 #elif SEA_LEVEL==3
215 
216 ! ------ z_sl from the SPECMAP record
217 
218 if (time/year_sec.lt.real(specmap_time_min,dp)) then
219  z_sl = specmap_zsl(0)
220 else if (time/year_sec.lt.real(specmap_time_max,dp)) then
221 
222  i_kl = floor(((time/year_sec) &
223  -real(specmap_time_min,dp))/real(specmap_time_stp,dp))
224  i_kl = max(i_kl, 0)
225 
226  i_gr = ceiling(((time/year_sec) &
227  -real(specmap_time_min,dp))/real(specmap_time_stp,dp))
228  i_gr = min(i_gr, ndata_specmap)
229 
230  if (i_kl.eq.i_gr) then
231 
232  z_sl = specmap_zsl(i_kl)
233 
234  else
235 
236  time_kl = (specmap_time_min + i_kl*specmap_time_stp) *year_sec
237  time_gr = (specmap_time_min + i_gr*specmap_time_stp) *year_sec
238 
239  z_sl = specmap_zsl(i_kl) &
240  +(specmap_zsl(i_gr)-specmap_zsl(i_kl)) &
241  *(time-time_kl)/(time_gr-time_kl)
242  ! linear interpolation of the sea-level data
243 
244  end if
245 
246 else
247  z_sl = specmap_zsl(ndata_specmap)
248 end if
249 
250 #endif
251 
252 ! ------ Time derivative of the sea level
253 
254 if ( z_sl_old > -999999.9_dp ) then
255  dzsl_dtau = (z_sl-z_sl_old)/dtime
256 else ! only dummy value for z_sl_old available
257  dzsl_dtau = 0.0_dp
258 end if
259 
260 ! ------ Minimum bedrock elevation for extent of marine ice
261 
262 #if MARGIN==2
263 
264 #if ( MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3 )
265 z_mar = z_mar
266 #elif ( MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 )
267 z_mar = fact_z_mar*z_sl
268 #elif ( MARINE_ICE_CALVING==6 || MARINE_ICE_CALVING==7 )
269 if (z_sl >= -80.0_dp) then
270  z_mar = 2.5_dp*z_sl
271 else
272  z_mar = 10.25_dp*(z_sl+80.0_dp)-200.0_dp
273 end if
274 z_mar = fact_z_mar*z_mar
275 #endif
276 
277 #endif
278 
279 ! ------ Update of the mask according to the sea level
280 
281 ! ---- Check all sea and floating-ice points and their direct
282 ! neighbours
283 
284 do i=0, imax
285 do j=0, jmax
286  check_point(j,i) = .false.
287 end do
288 end do
289 
290 do i=1, imax-1
291 do j=1, jmax-1
292  if (maske(j,i).ge.2) then
293  check_point(j ,i ) = .true.
294  check_point(j ,i+1) = .true.
295  check_point(j ,i-1) = .true.
296  check_point(j+1,i ) = .true.
297  check_point(j-1,i ) = .true.
298  end if
299 end do
300 end do
301 
302 do i=1, imax-1
303 do j=1, jmax-1
304  if (check_point(j,i)) then
305  maske_neu(j,i) = mask_update(z_sl, i, j)
306  end if
307 end do
308 end do
309 
310 ! ---- Assign new values of the mask
311 
312 do i=1, imax-1
313 do j=1, jmax-1
314  if (check_point(j,i)) then
315  maske(j,i) = maske_neu(j,i)
316  end if
317 end do
318 end do
319 
320 !-------- Surface air temperatures --------
321 
322 gamma_t = -6.5e-03_dp ! atmospheric lapse rate
323 
324 do i=0, imax
325 do j=0, jmax
326 
327 #if (TSURFACE <= 4)
328 
329 ! ------ Correction of present monthly temperatures with elevation changes
330 ! and temperature deviation delta_ts
331 
332  temp_diff = gamma_t*(zs(j,i)-zs_ref(j,i)) + delta_ts
333 
334  do n=1, 12 ! month counter
335  temp_mm(j,i,n) = temp_mm_present(j,i,n) + temp_diff
336  end do
337 
338 #elif (TSURFACE == 5)
339 
340 ! ------ Correction of present monthly temperatures with LGM anomaly and
341 ! glacial index as well as elevation changes
342 
343  temp_diff = gamma_t*(zs(j,i)-zs_ref(j,i))
344 
345  do n=1, 12 ! month counter
346  temp_mm(j,i,n) = temp_mm_present(j,i,n) &
347  + glac_index*temp_mm_lgm_anom(j,i,n) &
348  + temp_diff
349  end do
350 
351 #endif
352 
353 ! ------ Mean annual air temperature
354 
355  temp_ma(j,i) = 0.0_dp ! initialisation value
356 
357  do n=1, 12 ! month counter
358  temp_ma(j,i) = temp_ma(j,i) + temp_mm(j,i,n)*inv_twelve
359  end do
360 
361 end do
362 end do
363 
364 !-------- Accumulation-ablation function as_perp --------
365 
366 #if (ELEV_DESERT == 1)
367 
368 gamma_p = gamma_p*1.0e-03_dp ! Precipitation lapse rate
369  ! for elevation desertification, in m^(-1)
370 zs_thresh = zs_thresh ! Elevation threshold, in m
371 
372 #endif
373 
374 #if (SOLID_PRECIP == 1) /* Marsiat (1994) */
375 
376 temp_rain = 7.0_dp ! Threshold monthly mean temperature for
377  ! precipitation = 100% rain, in deg C
378 temp_snow = -10.0_dp ! Threshold monthly mean temperature for &
379  ! precipitation = 100% snow, in deg C
380 
381 inv_delta_temp_rain_snow = 1.0_dp/(temp_rain-temp_snow)
382 
383 #elif (SOLID_PRECIP == 2) /* Bales et al. (2009) */
384 
385 temp_rain = 7.2_dp ! Threshold monthly mean temperature for &
386  ! precipitation = 100% rain, in deg C
387 temp_snow = -11.6_dp ! Threshold monthly mean temperature for &
388  ! precipitation = 100% snow, in deg C
389 
390 coeff(0) = 5.4714e-01_dp ! Coefficients
391 coeff(1) = -9.1603e-02_dp ! of
392 coeff(2) = -3.314e-03_dp ! the
393 coeff(3) = 4.66e-04_dp ! fifth-order
394 coeff(4) = 3.8e-05_dp ! polynomial
395 coeff(5) = 6.0e-07_dp ! fit
396 
397 #elif (SOLID_PRECIP == 3) /* Huybrechts and de Wolde (1999) */
398 
399 temp_rain = 2.0_dp ! Threshold instantaneous temperature for &
400  ! precipitation = 100% rain, in deg C
401 temp_snow = temp_rain ! Threshold instantaneous temperature for &
402  ! precipitation = 100% snow, in deg C
403 
404 s_stat = s_stat_0 ! Standard deviation of the air termperature
405  ! (same parameter as in the PDD model)
406 
407 inv_sqrt2_s_stat = 1.0_dp/(sqrt(2.0_dp)*s_stat)
408 
409 #endif
410 
411 #if (ABLSURFACE==1 || ABLSURFACE==2)
412 
413 s_stat = s_stat_0
414 beta1 = beta1_0 *(0.001_dp/86400.0_dp)*(rho_w/rho)
415  ! (mm WE)/(d*deg C) --> (m IE)/(s*deg C)
416 beta2 = beta2_0 *(0.001_dp/86400.0_dp)*(rho_w/rho)
417  ! (mm WE)/(d*deg C) --> (m IE)/(s*deg C)
418 pmax = pmax_0
419 mu = mu_0 *(1000.0_dp*86400.0_dp)*(rho/rho_w)
420  ! (d*deg C)/(mm WE) --> (s*deg C)/(m IE)
421 
422 #elif (ABLSURFACE==3)
423 
424 lambda_lti = lambda_lti *(0.001_dp/year_sec)*(rho_w/rho)
425  ! (mm WE)/(a*deg C) --> (m IE)/(s*deg C)
426 temp_lti = temp_lti
427 mu = 0.0_dp ! no superimposed ice considered
428 
429 #endif
430 
431 do i=0, imax
432 do j=0, jmax
433 
434 ! ------ Accumulation
435 
436 #if (ACCSURFACE <= 3)
437 
438 ! ---- Elevation desertification of precipitation
439 
440 #if (ELEV_DESERT == 0)
441 
442  precip_fact = 1.0_dp ! no elevation desertification
443 
444 #elif (ELEV_DESERT == 1)
445 
446  if (zs_ref(j,i) < zs_thresh) then
447  precip_fact &
448  = exp(gamma_p*(max(zs(j,i),zs_thresh)-zs_thresh))
449  else
450  precip_fact &
451  = exp(gamma_p*(max(zs(j,i),zs_thresh)-zs_ref(j,i)))
452  end if
453 
454 #else
455  stop ' boundary: Parameter ELEV_DESERT must be either 0 or 1!'
456 #endif
457 
458  do n=1, 12 ! month counter
459  precip(j,i,n) = precip_present(j,i,n)*precip_fact
460  end do
461 
462 #endif
463 
464 ! ---- Precipitation change related to changing climate
465 
466 #if ACCSURFACE==1
467  precip_fact = accfact
468 #elif ACCSURFACE==2
469  precip_fact = 1.0_dp + gamma_s*delta_ts
470 #elif ACCSURFACE==3
471  precip_fact = exp(gamma_s*delta_ts)
472 #endif
473 
474 #if (ACCSURFACE <= 3)
475 
476  precip(j,i,0) = 0.0_dp ! initialisation value for mean annual precip
477 
478  do n=1, 12 ! month counter
479  precip(j,i,n) = precip(j,i,n)*precip_fact ! monthly precip
480  precip(j,i,0) = precip(j,i,0) + precip(j,i,n)*inv_twelve
481  ! mean annual precip
482  end do
483 
484 #elif (ACCSURFACE == 5)
485 
486  precip(j,i,0) = 0.0_dp ! initialisation value for mean annual precip
487 
488  do n=1, 12 ! month counter
489 
490 #if (PRECIP_ANOM_INTERPOL==1)
491  precip_fact = 1.0_dp-glac_index+glac_index*precip_lgm_anom(j,i,n)
492  ! interpolation with a linear function
493 #elif (PRECIP_ANOM_INTERPOL==2)
494  precip_fact = exp(-glac_index*gamma_precip_lgm_anom(j,i,n))
495  ! interpolation with an exponential function
496 #endif
497 
498  precip(j,i,n) = precip_present(j,i,n)*precip_fact ! monthly precip
499  precip(j,i,0) = precip(j,i,0) + precip(j,i,n)*inv_twelve
500  ! mean annual precip
501  end do
502 
503 #endif
504 
505 ! ---- Annual accumulation, snowfall and rainfall rates
506 
507  accum(j,i) = precip(j,i,0)
508 
509  snowfall(j,i) = 0.0_dp ! initialisation value
510 
511  do n=1, 12 ! month counter
512 
513 #if (SOLID_PRECIP == 1) /* Marsiat (1994) */
514 
515  if (temp_mm(j,i,n) >= temp_rain) then
516  frac_solid = 0.0_dp
517  else if (temp_mm(j,i,n) <= temp_snow) then
518  frac_solid = 1.0_dp
519  else
520  frac_solid = (temp_rain-temp_mm(j,i,n))*inv_delta_temp_rain_snow
521  end if
522 
523 #elif (SOLID_PRECIP == 2) /* Bales et al. (2009) */
524 
525  if (temp_mm(j,i,n) >= temp_rain) then
526  frac_solid = 0.0_dp
527  else if (temp_mm(j,i,n) <= temp_snow) then
528  frac_solid = 1.0_dp
529  else
530  frac_solid = coeff(0) + temp_mm(j,i,n) * ( coeff(1) &
531  + temp_mm(j,i,n) * ( coeff(2) &
532  + temp_mm(j,i,n) * ( coeff(3) &
533  + temp_mm(j,i,n) * ( coeff(4) &
534  + temp_mm(j,i,n) * coeff(5) ) ) ) )
535  ! evaluation of 5th-order polynomial by Horner scheme
536  end if
537 
538 #elif (SOLID_PRECIP == 3) /* Huybrechts and de Wolde (1999) */
539 
540  frac_solid = 1.0_dp &
541  - 0.5_dp*erfcc((temp_rain-temp_mm(j,i,n))*inv_sqrt2_s_stat)
542 
543 #endif
544 
545  snowfall(j,i) = snowfall(j,i) + precip(j,i,n)*frac_solid*inv_twelve
546 
547  end do
548 
549  rainfall(j,i) = precip(j,i,0) - snowfall(j,i)
550 
551  if (snowfall(j,i) < 0.0_dp) snowfall(j,i) = 0.0_dp ! correction of
552  if (rainfall(j,i) < 0.0_dp) rainfall(j,i) = 0.0_dp ! negative values
553 
554 ! ------ Ablation
555 
556 ! ---- Runoff
557 
558 #if (ABLSURFACE==1 || ABLSURFACE==2)
559 
560 ! -- Temperature excess ET
561 
562  do n=1, 12 ! month counter
563  temp_mm_help(n) = temp_mm(j,i,n)
564  end do
565 
566  call pdd(temp_mm_help, s_stat, et(j,i))
567 
568 ! -- Formation rate of superimposed ice (melt_star), melt rate (melt)
569 ! and runoff rate (runoff)
570 
571 #if (ABLSURFACE==1)
572 
573  if ((beta1*et(j,i)) <= (pmax*snowfall(j,i))) then
574  melt_star(j,i) = beta1*et(j,i)
575  melt(j,i) = 0.0_dp
576  runoff(j,i) = melt(j,i)+rainfall(j,i)
577  else
578  melt_star(j,i) = pmax*snowfall(j,i)
579  melt(j,i) = beta2*(et(j,i)-melt_star(j,i)/beta1)
580  runoff(j,i) = melt(j,i)+rainfall(j,i)
581  end if
582 
583 #elif (ABLSURFACE==2)
584 
585  if ( rainfall(j,i) <= (pmax*snowfall(j,i)) ) then
586 
587  if ( (rainfall(j,i)+beta1*et(j,i)) <= (pmax*snowfall(j,i)) ) then
588  melt_star(j,i) = rainfall(j,i)+beta1*et(j,i)
589  melt(j,i) = 0.0_dp
590  runoff(j,i) = melt(j,i)
591  else
592  melt_star(j,i) = pmax*snowfall(j,i)
593  melt(j,i) = beta2 &
594  *(et(j,i)-(melt_star(j,i)-rainfall(j,i))/beta1)
595  runoff(j,i) = melt(j,i)
596  end if
597 
598  else
599 
600  melt_star(j,i) = pmax*snowfall(j,i)
601  melt(j,i) = beta2*et(j,i)
602  runoff(j,i) = melt(j,i) + rainfall(j,i)-pmax*snowfall(j,i)
603 
604  end if
605 
606 #endif
607 
608 #elif (ABLSURFACE==3)
609 
610  temp_jja_help = one_third*(temp_mm(j,i,6)+temp_mm(j,i,7)+temp_mm(j,i,8))
611 
612  melt_star(j,i) = 0.0_dp ! no superimposed ice considered
613  melt(j,i) = lambda_lti*max((temp_jja_help-temp_lti), 0.0_dp)
614  runoff(j,i) = melt(j,i) + rainfall(j,i)
615 
616 #endif
617 
618 ! ---- Evaporation
619 
620  evap(j,i) = 0.0_dp
621 
622 ! ------ Accumulation minus ablation
623 
624  as_perp(j,i) = accum(j,i) - evap(j,i) - runoff(j,i)
625 
626 ! ------ Ice-surface temperature (10-m firn temperature) temp_s,
627 ! including empirical firn-warming correction due to
628 ! refreezing meltwater when superimposed ice is formed
629 
630  if (melt_star(j,i).ge.melt(j,i)) then
631  temp_s(j,i) = temp_ma(j,i) &
632  +mu*(melt_star(j,i)-melt(j,i))
633  else
634  temp_s(j,i) = temp_ma(j,i)
635  end if
636 
637  if (temp_s(j,i) > -0.001_dp) temp_s(j,i) = -0.001_dp
638  ! Cut-off of positive air temperatures
639 
640 end do
641 end do
642 
643 end subroutine boundary
644 !
Generated on Mon May 20 2013 10:18:34 for SICOPOLIS V3.0 by   doxygen 1.8.1