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