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