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 #if (NETCDF > 1)
43 use netcdf
44 #endif
45 
46 implicit none
47 
48 real(dp), intent(in) :: time, dtime, dxi, deta
49 
50 real(dp), intent(out) :: delta_ts, glac_index, dzsl_dtau, z_mar
51 real(dp), intent(inout) :: z_sl
52 
53 ! Further return variables
54 ! (defined as global variables in module sico_variables):
55 !
56 ! accum(j,i), evap(j,i), runoff(j,i), as_perp(j,i), temp_s(j,i)
57 
58 integer(i2b) :: mask_update
59 integer(i4b) :: i, j, n
60 integer(i4b) :: i_gr, i_kl
61 integer(i4b) :: nrec_temp_precip
62 integer(i4b), save :: nrec_temp_precip_save = -1
63 real(dp) :: z_sl_old
64 real(dp) :: z_sl_min, t1, t2, t3, t4, t5, t6
65 real(dp) :: time_gr, time_kl
66 real(dp) :: z_sle_present, z_sle_help
67 real(dp), dimension(0:JMAX,0:IMAX,0:12) :: precip
68 real(dp), dimension(0:JMAX,0:IMAX) :: &
69  snowfall, rainfall, melt, melt_star
70 real(dp), dimension(0:JMAX,0:IMAX,12) :: temp_mm
71 real(dp), dimension(0:JMAX,0:IMAX) :: temp_ma, temp_ampl
72 real(dp), dimension(0:JMAX,0:IMAX) :: temp_ma_anom, temp_mj_anom, precip_ma_anom
73 real(dp), dimension(0:IMAX,0:JMAX), save :: temp_ma_anom_tra, temp_mj_anom_tra, &
74  precip_ma_anom_tra
75  ! These arrays are transposed with respect
76  ! to the usual SICOPOLIS convention;
77  ! they have i as first and j as second index
78 real(dp), dimension(12) :: temp_mm_help
79 real(dp) :: temp_jja_help
80 real(dp), dimension(0:JMAX,0:IMAX) :: et
81 real(dp) :: theta_ma, c_ma, gamma_ma, &
82  theta_ma_1, c_ma_1, gamma_ma_1, &
83  theta_ma_2, c_ma_2, gamma_ma_2, &
84  theta_ma_3, c_ma_3, gamma_ma_3, &
85  zs_sep_1, zs_sep_2, &
86  theta_mj, c_mj, gamma_mj
87 real(dp) :: sine_factor
88 real(dp) :: gamma_p, zs_thresh, &
89  alpha_p, beta_p, temp_0, alpha_t, beta_t, &
90  temp_inv, temp_inv_present, &
91  temp_rain, temp_snow, &
92  inv_delta_temp_rain_snow, coeff(0:5), inv_sqrt2_s_stat, &
93  precip_fact, frac_solid
94 real(dp) :: s_stat, beta1, beta2, pmax, mu, lambda_lti, temp_lti
95 real(dp) :: erfcc
96 logical, dimension(0:JMAX,0:IMAX) :: check_point
97 
98 #if (NETCDF > 1)
99 integer(i4b) :: ncv
100 ! ncv: Variable ID
101 integer(i4b) :: nc3cor(3)
102 ! nc3cor(3): Corner of a 3-d array
103 integer(i4b) :: nc3cnt(3)
104 ! nc3cnt(3): Count of a 3-d array
105 #endif
106 
107 real(dp), parameter :: &
108  inv_twelve = 1.0_dp/12.0_dp, one_third = 1.0_dp/3.0_dp
109 
110 !-------- Initialization of intent(out) variables --------
111 
112 z_sl_old = z_sl
113 
114 delta_ts = 0.0_dp
115 glac_index = 0.0_dp
116 z_sl = 0.0_dp
117 dzsl_dtau = 0.0_dp
118 z_mar = 0.0_dp
119 
120 !-------- Surface-temperature deviation from present values --------
121 
122 #if TSURFACE==1
123 delta_ts = delta_ts0
124 ! ! Steady state with prescribed constant
125 ! ! air-temperature deviation
126 #elif TSURFACE==3
127 delta_ts = sine_amplit &
128  *cos(2.0_dp*pi*time/(sine_period*year_sec)) &
129  -sine_amplit
130 ! ! Sinusoidal air-temperature forcing
131 #elif TSURFACE==4
132 
133 ! ------ delta_ts from ice-core record
134 
135 if (time/year_sec < real(grip_time_min,dp)) then
136  delta_ts = griptemp(0)
137 else if (time/year_sec < real(grip_time_max,dp)) then
138 
139  i_kl = floor(((time/year_sec) &
140  -real(grip_time_min,dp))/real(grip_time_stp,dp))
141  i_kl = max(i_kl, 0)
142 
143  i_gr = ceiling(((time/year_sec) &
144  -real(grip_time_min,dp))/real(grip_time_stp,dp))
145  i_gr = min(i_gr, ndata_grip)
146 
147  if (i_kl == i_gr) then
148 
149  delta_ts = griptemp(i_kl)
150 
151  else
152 
153  time_kl = (grip_time_min + i_kl*grip_time_stp) *year_sec
154  time_gr = (grip_time_min + i_gr*grip_time_stp) *year_sec
155 
156  delta_ts = griptemp(i_kl) &
157  +(griptemp(i_gr)-griptemp(i_kl)) &
158  *(time-time_kl)/(time_gr-time_kl)
159  ! linear interpolation of the ice-core data
160 
161  end if
162 
163 else
164  delta_ts = griptemp(ndata_grip)
165 end if
166 
167 delta_ts = delta_ts * grip_temp_fact
168 ! ! modification by constant factor
169 
170 !-------- Glacial index --------
171 
172 #elif TSURFACE==5
173 
174 if (time/year_sec < real(gi_time_min,dp)) then
175  glac_index = glacial_index(0)
176 else if (time/year_sec < real(gi_time_max,dp)) then
177 
178  i_kl = floor(((time/year_sec) &
179  -real(gi_time_min,dp))/real(gi_time_stp,dp))
180  i_kl = max(i_kl, 0)
181 
182  i_gr = ceiling(((time/year_sec) &
183  -real(gi_time_min,dp))/real(gi_time_stp,dp))
184  i_gr = min(i_gr, ndata_gi)
185 
186  if (i_kl == i_gr) then
187 
188  glac_index = glacial_index(i_kl)
189 
190  else
191 
192  time_kl = (gi_time_min + i_kl*gi_time_stp) *year_sec
193  time_gr = (gi_time_min + i_gr*gi_time_stp) *year_sec
194 
195  glac_index = glacial_index(i_kl) &
196  +(glacial_index(i_gr)-glacial_index(i_kl)) &
197  *(time-time_kl)/(time_gr-time_kl)
198  ! linear interpolation of the glacial-index data
199 
200  end if
201 
202 else
203  glac_index = glacial_index(ndata_gi)
204 end if
205 
206 !-------- Reading of surface temperature and precipitation
207 ! directly from data --------
208 
209 #elif ( (TSURFACE==6) && (ACCSURFACE==6) )
210 
211 if (time/year_sec <= temp_precip_time_min) then
212  nrec_temp_precip = 0
213 else if (time/year_sec < temp_precip_time_max) then
214  nrec_temp_precip = nint( ((time/year_sec)-temp_precip_time_min) &
215  / temp_precip_time_stp )
216 else
217  nrec_temp_precip = ndata_temp_precip
218 end if
219 
220 if ( nrec_temp_precip < 0 ) then
221  stop ' boundary: nrec_temp_precip < 0, not allowed!'
222 else if ( nrec_temp_precip > ndata_temp_precip ) then
223  stop ' boundary: nrec_temp_precip > ndata_temp_precip, not allowed!'
224 end if
225 
226 if ( nrec_temp_precip /= nrec_temp_precip_save ) then
227 
228  nc3cor(1) = 1
229  nc3cor(2) = 1
230  nc3cor(3) = nrec_temp_precip + 1
231  nc3cnt(1) = imax + 1
232  nc3cnt(2) = jmax + 1
233  nc3cnt(3) = 1
234 
235  call check( nf90_inq_varid(ncid_temp_precip, 'annualtemp_anom', ncv) )
236  call check( nf90_get_var(ncid_temp_precip, ncv, temp_ma_anom_tra, &
237  start=nc3cor, count=nc3cnt) )
238 
239  call check( nf90_inq_varid(ncid_temp_precip, 'januarytemp_anom', ncv) )
240  call check( nf90_get_var(ncid_temp_precip, ncv, temp_mj_anom_tra, &
241  start=nc3cor, count=nc3cnt) )
242 
243  call check( nf90_inq_varid(ncid_temp_precip, 'precipitation_anom', ncv) )
244  call check( nf90_get_var(ncid_temp_precip, ncv, precip_ma_anom_tra, &
245  start=nc3cor, count=nc3cnt) )
246 
247 end if
248 
249 temp_ma_anom = transpose(temp_ma_anom_tra)
250 temp_mj_anom = transpose(temp_mj_anom_tra)
251 precip_ma_anom = transpose(precip_ma_anom_tra) *(1.0e-03_dp/year_sec)*(rho_w/rho)
252  ! mm/a water equiv. -> m/s ice equiv.
253 
254 nrec_temp_precip_save = nrec_temp_precip
255 
256 #endif
257 
258 !-------- Sea level --------
259 
260 #if SEA_LEVEL==1
261 ! ------ constant sea level
262 z_sl = z_sl0
263 
264 #elif SEA_LEVEL==2
265 ! ------ saw-tooth-shaped palaeoclimatic sea-level forcing
266 
267 z_sl_min = -130.0_dp
268 
269 t1 = -250000.0_dp *year_sec
270 t2 = -140000.0_dp *year_sec
271 t3 = -125000.0_dp *year_sec
272 t4 = -21000.0_dp *year_sec
273 t5 = -8000.0_dp *year_sec
274 t6 = 0.0_dp *year_sec
275 
276 if (time < t1) then
277  z_sl = 0.0_dp
278 else if (time < t2) then
279  z_sl = z_sl_min*(time-t1)/(t2-t1)
280 else if (time < t3) then
281  z_sl = -z_sl_min*(time-t3)/(t3-t2)
282 else if (time < t4) then
283  z_sl = z_sl_min*(time-t3)/(t4-t3)
284 else if (time < t5) then
285  z_sl = -z_sl_min*(time-t5)/(t5-t4)
286 else if (time < t6) then
287  z_sl = 0.0_dp
288 else
289  z_sl = 0.0_dp
290 end if
291 
292 #elif SEA_LEVEL==3
293 
294 ! ------ z_sl from the SPECMAP record
295 
296 if (time/year_sec < real(specmap_time_min,dp)) then
297  z_sl = specmap_zsl(0)
298 else if (time/year_sec < real(specmap_time_max,dp)) then
299 
300  i_kl = floor(((time/year_sec) &
301  -real(specmap_time_min,dp))/real(specmap_time_stp,dp))
302  i_kl = max(i_kl, 0)
303 
304  i_gr = ceiling(((time/year_sec) &
305  -real(specmap_time_min,dp))/real(specmap_time_stp,dp))
306  i_gr = min(i_gr, ndata_specmap)
307 
308  if (i_kl == i_gr) then
309 
310  z_sl = specmap_zsl(i_kl)
311 
312  else
313 
314  time_kl = (specmap_time_min + i_kl*specmap_time_stp) *year_sec
315  time_gr = (specmap_time_min + i_gr*specmap_time_stp) *year_sec
316 
317  z_sl = specmap_zsl(i_kl) &
318  +(specmap_zsl(i_gr)-specmap_zsl(i_kl)) &
319  *(time-time_kl)/(time_gr-time_kl)
320  ! linear interpolation of the sea-level data
321 
322  end if
323 
324 else
325  z_sl = specmap_zsl(ndata_specmap)
326 end if
327 
328 #endif
329 
330 ! ------ Time derivative of the sea level
331 
332 if ( z_sl_old > -999999.9_dp ) then
333  dzsl_dtau = (z_sl-z_sl_old)/dtime
334 else ! only dummy value for z_sl_old available
335  dzsl_dtau = 0.0_dp
336 end if
337 
338 ! ------ Minimum bedrock elevation for extent of marine ice
339 
340 #if MARGIN==2
341 
342 #if ( MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3 )
343 z_mar = z_mar
344 #elif ( MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 )
345 z_mar = fact_z_mar*z_sl
346 #elif ( MARINE_ICE_CALVING==6 || MARINE_ICE_CALVING==7 )
347 if (z_sl >= -80.0_dp) then
348  z_mar = 2.5_dp*z_sl
349 else
350  z_mar = 10.25_dp*(z_sl+80.0_dp)-200.0_dp
351 end if
352 z_mar = fact_z_mar*z_mar
353 #endif
354 
355 #endif
356 
357 ! ------ Update of the mask according to the sea level
358 
359 ! ---- Check all sea and floating-ice points and their direct
360 ! neighbours
361 
362 do i=0, imax
363 do j=0, jmax
364  check_point(j,i) = .false.
365 end do
366 end do
367 
368 do i=1, imax-1
369 do j=1, jmax-1
370  if (maske(j,i) >= 2) then
371  check_point(j ,i ) = .true.
372  check_point(j ,i+1) = .true.
373  check_point(j ,i-1) = .true.
374  check_point(j+1,i ) = .true.
375  check_point(j-1,i ) = .true.
376  end if
377 end do
378 end do
379 
380 do i=1, imax-1
381 do j=1, jmax-1
382  if (check_point(j,i)) then
383  maske_neu(j,i) = mask_update(z_sl, i, j)
384  end if
385 end do
386 end do
387 
388 ! ---- Assign new values of the mask
389 
390 do i=1, imax-1
391 do j=1, jmax-1
392  if (check_point(j,i)) then
393  maske(j,i) = maske_neu(j,i)
394  end if
395 end do
396 end do
397 
398 !-------- Surface air temperatures --------
399 
400 #if TEMP_PRESENT_PARA == 1 /* Parameterisation by Fortuin and Oerlemans */
401  ! (1990) for the whole ice sheet
402 
403 theta_ma = 34.461_dp
404 gamma_ma = -9.14e-03_dp
405 c_ma = -0.688_dp
406 
407 theta_mj = 16.81_dp
408 gamma_mj = -6.92e-03_dp
409 c_mj = -0.27973_dp
410 
411 #elif TEMP_PRESENT_PARA == 2 /* Parameterisation by Fortuin and Oerlemans */
412  ! (1990), separately for three different
413  ! elevation ranges
414 
415 zs_sep_1 = 200.0_dp
416 zs_sep_2 = 1500.0_dp
417 
418 theta_ma_1 = 49.642_dp
419 gamma_ma_1 = 0.0_dp
420 c_ma_1 = -0.943_dp
421 
422 theta_ma_2 = 36.689_dp
423 gamma_ma_2 = -5.102e-03_dp
424 c_ma_2 = -0.725_dp
425 
426 theta_ma_3 = 7.405_dp
427 gamma_ma_3 = -14.285e-03_dp
428 c_ma_3 = -0.180_dp
429 
430 theta_mj = 16.81_dp
431 gamma_mj = -6.92e-03_dp
432 c_mj = -0.27937_dp
433 
434 #else
435 
436 stop ' boundary: Parameter TEMP_PRESENT_PARA must be either 1 or 2!'
437 
438 #endif
439 
440 do i=0, imax
441 do j=0, jmax
442 
443 ! ------ Present-day mean-annual air temperature
444 
445 #if TEMP_PRESENT_PARA == 1
446  temp_ma_present(j,i) = theta_ma + gamma_ma*zs(j,i) &
447  + c_ma*abs(phi(j,i))*pi_180_inv
448 #elif TEMP_PRESENT_PARA == 2
449  if ( zs(j,i) <= zs_sep_1 ) then
450  temp_ma_present(j,i) = theta_ma_1 + gamma_ma_1*zs(j,i) &
451  + c_ma_1*abs(phi(j,i))*pi_180_inv
452  else if ( zs(j,i) <= zs_sep_2 ) then
453  temp_ma_present(j,i) = theta_ma_2 + gamma_ma_2*zs(j,i) &
454  + c_ma_2*abs(phi(j,i))*pi_180_inv
455  else
456  temp_ma_present(j,i) = theta_ma_3 + gamma_ma_3*zs(j,i) &
457  + c_ma_3*abs(phi(j,i))*pi_180_inv
458  end if
459 #endif
460 
461 ! ------ Present-day mean-January (summer) air temperature
462 
463  temp_mj_present(j,i) = theta_mj + gamma_mj*zs(j,i) &
464  + c_mj*abs(phi(j,i))*pi_180_inv
465 
466 #if (TSURFACE <= 4)
467 
468 ! ------ Correction with deviation delta_ts
469 
470  temp_ma(j,i) = temp_ma_present(j,i) + delta_ts
471  temp_mm(j,i,7) = temp_mj_present(j,i) + delta_ts
472 
473 #elif (TSURFACE == 5)
474 
475 ! ------ Correction with LGM anomaly and glacial index
476 
477  temp_ma(j,i) = temp_ma_present(j,i) + glac_index*temp_ma_lgm_anom(j,i)
478  temp_mm(j,i,7) = temp_mj_present(j,i) + glac_index*temp_mj_lgm_anom(j,i)
479 
480 #elif (TSURFACE == 6)
481 
482 ! ------ Mean-annual and mean-January (summer) surface temperatures
483 ! from data read above --------
484 
485  temp_ma(j,i) = temp_ma_present(j,i) + temp_ma_anom_fact*temp_ma_anom(j,i)
486  temp_mm(j,i,7) = temp_mj_present(j,i) + temp_mj_anom_fact*temp_mj_anom(j,i)
487 
488 #endif
489 
490 ! ------ Amplitude of the annual cycle
491 
492  temp_ampl(j,i) = temp_mm(j,i,7) - temp_ma(j,i)
493 
494  if (temp_ampl(j,i) < eps) then
495  temp_ampl(j,i) = eps ! Correction of amplitude, if required
496  end if
497 
498 end do
499 end do
500 
501 ! ------ Monthly temperatures
502 
503 do n=1, 12 ! month counter
504 
505  sine_factor = sin((real(n,dp)-4.0_dp)*pi/6.0_dp)
506 
507  do i=0, imax
508  do j=0, jmax
509  temp_mm(j,i,n) = temp_ma(j,i) + sine_factor*temp_ampl(j,i)
510  end do
511  end do
512 
513 end do
514 
515 !-------- Accumulation-ablation function as_perp --------
516 
517 #if (ACCSURFACE <= 3)
518 
519 #if (ELEV_DESERT == 1)
520 
521 gamma_p = gamma_p*1.0e-03_dp ! Precipitation lapse rate
522  ! for elevation desertification, in m^(-1)
523 zs_thresh = zs_thresh ! Elevation threshold, in m
524 
525 #endif
526 
527 #elif (ACCSURFACE == 4)
528 
529 alpha_p = 22.47_dp
530 beta_p = 0.046_dp
531 temp_0 = 273.15_dp
532 alpha_t = 0.67_dp
533 beta_t = 88.9_dp
534 
535 #endif
536 
537 #if (SOLID_PRECIP == 1) /* Marsiat (1994) */
538 
539 temp_rain = 7.0_dp ! Threshold monthly mean temperature for
540  ! precipitation = 100% rain, in deg C
541 temp_snow = -10.0_dp ! Threshold monthly mean temperature for &
542  ! precipitation = 100% snow, in deg C
543 
544 inv_delta_temp_rain_snow = 1.0_dp/(temp_rain-temp_snow)
545 
546 #elif (SOLID_PRECIP == 2) /* Bales et al. (2009) */
547 
548 temp_rain = 7.2_dp ! Threshold monthly mean temperature for &
549  ! precipitation = 100% rain, in deg C
550 temp_snow = -11.6_dp ! Threshold monthly mean temperature for &
551  ! precipitation = 100% snow, in deg C
552 
553 coeff(0) = 5.4714e-01_dp ! Coefficients
554 coeff(1) = -9.1603e-02_dp ! of
555 coeff(2) = -3.314e-03_dp ! the
556 coeff(3) = 4.66e-04_dp ! fifth-order
557 coeff(4) = 3.8e-05_dp ! polynomial
558 coeff(5) = 6.0e-07_dp ! fit
559 
560 #elif (SOLID_PRECIP == 3) /* Huybrechts and de Wolde (1999) */
561 
562 temp_rain = 2.0_dp ! Threshold instantaneous temperature for &
563  ! precipitation = 100% rain, in deg C
564 temp_snow = temp_rain ! Threshold instantaneous temperature for &
565  ! precipitation = 100% snow, in deg C
566 
567 s_stat = s_stat_0 ! Standard deviation of the air termperature
568  ! (same parameter as in the PDD model)
569 
570 inv_sqrt2_s_stat = 1.0_dp/(sqrt(2.0_dp)*s_stat)
571 
572 #endif
573 
574 #if (ABLSURFACE==1 || ABLSURFACE==2)
575 
576 s_stat = s_stat_0
577 beta1 = beta1_0 *(0.001_dp/86400.0_dp)*(rho_w/rho)
578  ! (mm WE)/(d*deg C) --> (m IE)/(s*deg C)
579 beta2 = beta2_0 *(0.001_dp/86400.0_dp)*(rho_w/rho)
580  ! (mm WE)/(d*deg C) --> (m IE)/(s*deg C)
581 pmax = pmax_0
582 mu = mu_0 *(1000.0_dp*86400.0_dp)*(rho/rho_w)
583  ! (d*deg C)/(mm WE) --> (s*deg C)/(m IE)
584 
585 #elif (ABLSURFACE==3)
586 
587 lambda_lti = lambda_lti *(0.001_dp/year_sec)*(rho_w/rho)
588  ! (mm WE)/(a*deg C) --> (m IE)/(s*deg C)
589 temp_lti = temp_lti
590 mu = 0.0_dp ! no superimposed ice considered
591 
592 #endif
593 
594 do i=0, imax
595 do j=0, jmax
596 
597 ! ------ Accumulation
598 
599 #if (ACCSURFACE <= 3)
600 
601 ! ---- Elevation desertification of precipitation
602 
603 #if (ELEV_DESERT == 0)
604 
605  precip_fact = 1.0_dp ! no elevation desertification
606 
607 #elif (ELEV_DESERT == 1)
608 
609  if (zs_ref(j,i) < zs_thresh) then
610  precip_fact &
611  = exp(gamma_p*(max(zs(j,i),zs_thresh)-zs_thresh))
612  else
613  precip_fact &
614  = exp(gamma_p*(max(zs(j,i),zs_thresh)-zs_ref(j,i)))
615  end if
616 
617 #else
618  stop ' boundary: Parameter ELEV_DESERT must be either 0 or 1!'
619 #endif
620 
621  do n=1, 12 ! month counter
622  precip(j,i,n) = precip_present(j,i,n)*precip_fact
623  end do
624 
625 #endif
626 
627 ! ---- Precipitation change related to changing climate
628 
629 #if ACCSURFACE==1
630  precip_fact = accfact
631 #elif ACCSURFACE==2
632  precip_fact = 1.0_dp + gamma_s*delta_ts
633 #elif ACCSURFACE==3
634  precip_fact = exp(gamma_s*delta_ts)
635 #endif
636 
637 #if (ACCSURFACE <= 3)
638 
639  precip(j,i,0) = 0.0_dp ! initialisation value for mean annual precip
640 
641  do n=1, 12 ! month counter
642  precip(j,i,n) = precip(j,i,n)*precip_fact ! monthly precip
643  precip(j,i,0) = precip(j,i,0) + precip(j,i,n)*inv_twelve
644  ! mean annual precip
645  end do
646 
647 #elif (ACCSURFACE == 4)
648 
649  precip(j,i,0) = 0.0_dp ! initialisation value for mean annual precip
650 
651  temp_inv = alpha_t * (temp_ma(j,i)+temp_0) + beta_t ! in K
652  temp_inv_present = alpha_t * (temp_ma_present(j,i)+temp_0) + beta_t ! in K
653 
654  precip_fact = exp(alpha_p*(temp_0/temp_inv_present-temp_0/temp_inv)) &
655  *(temp_inv_present/temp_inv)**2 &
656  *(1.0_dp+beta_p*(temp_inv-temp_inv_present))
657 
658  do n=1, 12 ! month counter
659  precip(j,i,n) = precip_present(j,i,n)*precip_fact ! monthly precip
660  precip(j,i,0) = precip(j,i,0) + precip(j,i,n)*inv_twelve
661  ! mean annual precip
662  end do
663 
664 #elif (ACCSURFACE == 5)
665 
666  precip(j,i,0) = 0.0_dp ! initialisation value for mean annual precip
667 
668  do n=1, 12 ! month counter
669 
670 #if (PRECIP_ANOM_INTERPOL==1)
671  precip_fact = 1.0_dp-glac_index+glac_index*precip_lgm_anom(j,i,n)
672  ! interpolation with a linear function
673 #elif (PRECIP_ANOM_INTERPOL==2)
674  precip_fact = exp(-glac_index*gamma_precip_lgm_anom(j,i,n))
675  ! interpolation with an exponential function
676 #endif
677 
678  precip(j,i,n) = precip_present(j,i,n)*precip_fact ! monthly precip
679  precip(j,i,0) = precip(j,i,0) + precip(j,i,n)*inv_twelve
680  ! mean annual precip
681  end do
682 
683 #elif (ACCSURFACE == 6)
684 
685 ! -- Mean-annual precipitation from data already read above
686 
687  precip(j,i,0) = precip_ma_present(j,i) + precip_anom_fact*precip_ma_anom(j,i)
688  ! mean annual precip
689 
690  do n=1, 12 ! month counter
691  precip(j,i,n) = precip(j,i,0) ! monthly precip, assumed to be equal
692  ! to the mean annual precip
693  end do
694 
695 #endif
696 
697 ! ---- Annual accumulation, snowfall and rainfall rates
698 
699  accum(j,i) = precip(j,i,0)
700 
701  snowfall(j,i) = 0.0_dp ! initialisation value
702 
703  do n=1, 12 ! month counter
704 
705 #if (SOLID_PRECIP == 1) /* Marsiat (1994) */
706 
707  if (temp_mm(j,i,n) >= temp_rain) then
708  frac_solid = 0.0_dp
709  else if (temp_mm(j,i,n) <= temp_snow) then
710  frac_solid = 1.0_dp
711  else
712  frac_solid = (temp_rain-temp_mm(j,i,n))*inv_delta_temp_rain_snow
713  end if
714 
715 #elif (SOLID_PRECIP == 2) /* Bales et al. (2009) */
716 
717  if (temp_mm(j,i,n) >= temp_rain) then
718  frac_solid = 0.0_dp
719  else if (temp_mm(j,i,n) <= temp_snow) then
720  frac_solid = 1.0_dp
721  else
722  frac_solid = coeff(0) + temp_mm(j,i,n) * ( coeff(1) &
723  + temp_mm(j,i,n) * ( coeff(2) &
724  + temp_mm(j,i,n) * ( coeff(3) &
725  + temp_mm(j,i,n) * ( coeff(4) &
726  + temp_mm(j,i,n) * coeff(5) ) ) ) )
727  ! evaluation of 5th-order polynomial by Horner scheme
728  end if
729 
730 #elif (SOLID_PRECIP == 3) /* Huybrechts and de Wolde (1999) */
731 
732  frac_solid = 1.0_dp &
733  - 0.5_dp*erfcc((temp_rain-temp_mm(j,i,n))*inv_sqrt2_s_stat)
734 
735 #endif
736 
737  snowfall(j,i) = snowfall(j,i) + precip(j,i,n)*frac_solid*inv_twelve
738 
739  end do
740 
741  rainfall(j,i) = precip(j,i,0) - snowfall(j,i)
742 
743  if (snowfall(j,i) < 0.0_dp) snowfall(j,i) = 0.0_dp ! correction of
744  if (rainfall(j,i) < 0.0_dp) rainfall(j,i) = 0.0_dp ! negative values
745 
746 ! ------ Ablation
747 
748 ! ---- Runoff
749 
750 #if (ABLSURFACE==1 || ABLSURFACE==2)
751 
752 ! -- Temperature excess ET
753 
754  do n=1, 12 ! month counter
755  temp_mm_help(n) = temp_mm(j,i,n)
756  end do
757 
758  call pdd(temp_mm_help, s_stat, et(j,i))
759 
760 ! -- Formation rate of superimposed ice (melt_star), melt rate (melt)
761 ! and runoff rate (runoff)
762 
763 #if (ABLSURFACE==1)
764 
765  if ((beta1*et(j,i)) <= (pmax*snowfall(j,i))) then
766  melt_star(j,i) = beta1*et(j,i)
767  melt(j,i) = 0.0_dp
768  runoff(j,i) = melt(j,i)+rainfall(j,i)
769  else
770  melt_star(j,i) = pmax*snowfall(j,i)
771  melt(j,i) = beta2*(et(j,i)-melt_star(j,i)/beta1)
772  runoff(j,i) = melt(j,i)+rainfall(j,i)
773  end if
774 
775 #elif (ABLSURFACE==2)
776 
777  if ( rainfall(j,i) <= (pmax*snowfall(j,i)) ) then
778 
779  if ( (rainfall(j,i)+beta1*et(j,i)) <= (pmax*snowfall(j,i)) ) then
780  melt_star(j,i) = rainfall(j,i)+beta1*et(j,i)
781  melt(j,i) = 0.0_dp
782  runoff(j,i) = melt(j,i)
783  else
784  melt_star(j,i) = pmax*snowfall(j,i)
785  melt(j,i) = beta2 &
786  *(et(j,i)-(melt_star(j,i)-rainfall(j,i))/beta1)
787  runoff(j,i) = melt(j,i)
788  end if
789 
790  else
791 
792  melt_star(j,i) = pmax*snowfall(j,i)
793  melt(j,i) = beta2*et(j,i)
794  runoff(j,i) = melt(j,i) + rainfall(j,i)-pmax*snowfall(j,i)
795 
796  end if
797 
798 #endif
799 
800 #elif (ABLSURFACE==3)
801 
802  temp_jja_help = one_third*(temp_mm(j,i,6)+temp_mm(j,i,7)+temp_mm(j,i,8))
803 
804  melt_star(j,i) = 0.0_dp ! no superimposed ice considered
805  melt(j,i) = lambda_lti*max((temp_jja_help-temp_lti), 0.0_dp)
806  runoff(j,i) = melt(j,i) + rainfall(j,i)
807 
808 #endif
809 
810 ! ---- Evaporation
811 
812  evap(j,i) = 0.0_dp
813 
814 ! ------ Accumulation minus ablation
815 
816  as_perp(j,i) = accum(j,i) - evap(j,i) - runoff(j,i)
817 
818 ! ------ Ice-surface temperature (10-m firn temperature) temp_s,
819 ! including empirical firn-warming correction due to
820 ! refreezing meltwater when superimposed ice is formed
821 
822  if (melt_star(j,i) >= melt(j,i)) then
823  temp_s(j,i) = temp_ma(j,i) &
824  +mu*(melt_star(j,i)-melt(j,i))
825  else
826  temp_s(j,i) = temp_ma(j,i)
827  end if
828 
829  if (temp_s(j,i) > -0.001_dp) temp_s(j,i) = -0.001_dp
830  ! Cut-off of positive air temperatures
831 
832 end do
833 end do
834 
835 end subroutine boundary
836 !
Generated on Mon May 20 2013 10:18:34 for SICOPOLIS V3.0 by   doxygen 1.8.1