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  melt, melt_star
65 !real(dp), dimension(0:JMAX,0:IMAX) :: &
66 ! snowfall, rainfall, melt, melt_star
67 !snowfall, rainfall declared in sico_variables because of output5.F90
68 real(dp), dimension(0:JMAX,0:IMAX,12) :: temp_mm
69 real(dp), dimension(0:JMAX,0:IMAX) :: temp_ma
70 real(dp), dimension(12) :: temp_mm_help
71 real(dp) :: temp_jja_help
72 real(dp), dimension(0:JMAX,0:IMAX) :: et
73 real(dp) :: gamma_t, temp_diff
74 real(dp) :: gamma_p, zs_thresh, &
75  temp_rain, temp_snow, &
76  inv_delta_temp_rain_snow, coeff(0:5), inv_sqrt2_s_stat, &
77  precip_fact, frac_solid
78 real(dp) :: s_stat, beta1, beta2, pmax, mu, lambda_lti, temp_lti
79 real(dp) :: erfcc
80 logical, dimension(0:JMAX,0:IMAX) :: check_point
81 
82 real(dp), parameter :: &
83  inv_twelve = 1.0_dp/12.0_dp, one_third = 1.0_dp/3.0_dp
84 
85 
86 !-------- Initialization of intent(out) variables --------
87 
88 z_sl_old = z_sl
89 
90 delta_ts = 0.0_dp
91 glac_index = 0.0_dp
92 z_sl = 0.0_dp
93 dzsl_dtau = 0.0_dp
94 z_mar = 0.0_dp
95 
96 !-------- Surface-temperature deviation from present values --------
97 
98 #if TSURFACE==1
99 delta_ts = delta_ts0
100 
101 ! ! Steady state with prescribed constant
102 ! ! air-temperature deviation
103 #elif TSURFACE==3
104 delta_ts = sine_amplit &
105  *cos(2.0_dp*pi*time/(sine_period*year_sec)) &
106  -sine_amplit
107 ! ! Sinusoidal air-temperature forcing
108 #elif TSURFACE==4
109 
110 ! ------ delta_ts from ice-core record
111 
112 if (time/year_sec.lt.real(grip_time_min,dp)) then
113  delta_ts = griptemp(0)
114 else if (time/year_sec.lt.real(grip_time_max,dp)) then
115 
116  i_kl = floor(((time/year_sec)-real(grip_time_min,dp))/real(grip_time_stp,dp))
117  i_kl = max(i_kl, 0)
118 
119  i_gr = ceiling(((time/year_sec)-real(grip_time_min,dp))/real(grip_time_stp,dp))
120  i_gr = min(i_gr, ndata_grip)
121 
122  if (i_kl.eq.i_gr) then
123 
124  delta_ts = griptemp(i_kl)
125 
126  else
127 
128  time_kl = (grip_time_min + i_kl*grip_time_stp) *year_sec
129  time_gr = (grip_time_min + i_gr*grip_time_stp) *year_sec
130 
131  delta_ts = griptemp(i_kl) &
132  +(griptemp(i_gr)-griptemp(i_kl)) &
133  *(time-time_kl)/(time_gr-time_kl)
134  ! linear interpolation of the ice-core data
135 
136  end if
137 
138 else
139  delta_ts = griptemp(ndata_grip)
140 end if
141 
142 delta_ts = delta_ts * grip_temp_fact
143 ! ! modification by constant factor
144 
145 !-------- Glacial index --------
146 
147 #elif TSURFACE==5
148 
149 if (time/year_sec < real(gi_time_min,dp)) then
150  glac_index = glacial_index(0)
151 else if (time/year_sec < real(gi_time_max,dp)) then
152 
153  i_kl = floor(((time/year_sec)-real(gi_time_min,dp))/real(gi_time_stp,dp))
154  i_kl = max(i_kl, 0)
155 
156  i_gr = ceiling(((time/year_sec)-real(gi_time_min,dp))/real(gi_time_stp,dp))
157  i_gr = min(i_gr, ndata_gi)
158 
159  if (i_kl == i_gr) then
160 
161  glac_index = glacial_index(i_kl)
162 
163  else
164 
165  time_kl = (gi_time_min + i_kl*gi_time_stp) *year_sec
166  time_gr = (gi_time_min + i_gr*gi_time_stp) *year_sec
167 
168  glac_index = glacial_index(i_kl) &
169  +(glacial_index(i_gr)-glacial_index(i_kl)) &
170  *(time-time_kl)/(time_gr-time_kl)
171  ! linear interpolation of the glacial-index data
172 
173  end if
174 
175 else
176  glac_index = glacial_index(ndata_gi)
177 end if
178 
179 #endif
180 
181 !-------- Sea level --------
182 
183 #if SEA_LEVEL==1
184 ! ------ constant sea level
185 z_sl = z_sl0
186 
187 #elif SEA_LEVEL==2
188 ! ------ saw-tooth-shaped palaeoclimatic sea-level forcing
189 
190 z_sl_min = -130.0_dp
191 
192 t1 = -250000.0_dp *year_sec
193 t2 = -140000.0_dp *year_sec
194 t3 = -125000.0_dp *year_sec
195 t4 = -21000.0_dp *year_sec
196 t5 = -8000.0_dp *year_sec
197 t6 = 0.0_dp *year_sec
198 
199 if (time.lt.t1) then
200  z_sl = 0.0_dp
201 else if (time.lt.t2) then
202  z_sl = z_sl_min*(time-t1)/(t2-t1)
203 else if (time.lt.t3) then
204  z_sl = -z_sl_min*(time-t3)/(t3-t2)
205 else if (time.lt.t4) then
206  z_sl = z_sl_min*(time-t3)/(t4-t3)
207 else if (time.lt.t5) then
208  z_sl = -z_sl_min*(time-t5)/(t5-t4)
209 else if (time.lt.t6) then
210  z_sl = 0.0_dp
211 else
212  z_sl = 0.0_dp
213 end if
214 
215 #elif SEA_LEVEL==3
216 
217 ! ------ z_sl from the SPECMAP record
218 
219 if (time/year_sec.lt.real(specmap_time_min,dp)) then
220  z_sl = specmap_zsl(0)
221 else if (time/year_sec.lt.real(specmap_time_max,dp)) then
222 
223  i_kl = floor(((time/year_sec)-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)-real(specmap_time_min,dp))/real(specmap_time_stp,dp))
227  i_gr = min(i_gr, ndata_specmap)
228 
229  if (i_kl.eq.i_gr) then
230 
231  z_sl = specmap_zsl(i_kl)
232 
233  else
234 
235  time_kl = (specmap_time_min + i_kl*specmap_time_stp) *year_sec
236  time_gr = (specmap_time_min + i_gr*specmap_time_stp) *year_sec
237 
238  z_sl = specmap_zsl(i_kl) &
239  +(specmap_zsl(i_gr)-specmap_zsl(i_kl)) &
240  *(time-time_kl)/(time_gr-time_kl)
241  ! linear interpolation of the sea-level data
242 
243  end if
244 
245 else
246  z_sl = specmap_zsl(ndata_specmap)
247 end if
248 
249 #endif
250 
251 ! ------ Time derivative of the sea level
252 
253 if ( z_sl_old > -999999.9_dp ) then
254  dzsl_dtau = (z_sl-z_sl_old)/dtime
255 else ! only dummy value for z_sl_old available
256  dzsl_dtau = 0.0_dp
257 end if
258 
259 ! ------ Minimum bedrock elevation for extent of marine ice
260 
261 #if MARGIN==2
262 
263 #if ( MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3 )
264 z_mar = z_mar
265 #elif ( MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 )
266 z_mar = fact_z_mar*z_sl
267 #elif ( MARINE_ICE_CALVING==6 || MARINE_ICE_CALVING==7 )
268 if (z_sl >= -80.0_dp) then
269  z_mar = 2.5_dp*z_sl
270 else
271  z_mar = 10.25_dp*(z_sl+80.0_dp)-200.0_dp
272 end if
273 z_mar = fact_z_mar*z_mar
274 #endif
275 
276 #endif
277 
278 ! ------ Update of the mask according to the sea level
279 
280 ! ---- Check all sea and floating-ice points and their direct
281 ! neighbours
282 
283 do i=0, imax
284 do j=0, jmax
285  check_point(j,i) = .false.
286 end do
287 end do
288 
289 do i=1, imax-1
290 do j=1, jmax-1
291  if (maske(j,i).ge.2) then
292  check_point(j ,i ) = .true.
293  check_point(j ,i+1) = .true.
294  check_point(j ,i-1) = .true.
295  check_point(j+1,i ) = .true.
296  check_point(j-1,i ) = .true.
297  end if
298 end do
299 end do
300 
301 do i=1, imax-1
302 do j=1, jmax-1
303  if (check_point(j,i)) then
304  maske_neu(j,i) = mask_update(z_sl, i, j)
305  end if
306 end do
307 end do
308 
309 ! ---- Assign new values of the mask
310 
311 do i=1, imax-1
312 do j=1, jmax-1
313  if (check_point(j,i)) then
314  maske(j,i) = maske_neu(j,i)
315  end if
316 end do
317 end do
318 
319 !-------- Surface air temperatures --------
320 
321 gamma_t = -4.5e-03_dp ! atmospheric lapse rate
322 
323 do i=0, imax
324 do j=0, jmax
325 
326 #if (TSURFACE <= 4)
327 
328 ! ------ Correction of present monthly temperatures with elevation changes
329 ! and temperature deviation delta_ts
330 
331  temp_diff = gamma_t*(zs(j,i)-zs_ref(j,i)) + delta_ts
332 
333  do n=1, 12 ! month counter
334  temp_mm(j,i,n) = temp_mm_present(j,i,n) + temp_diff
335  end do
336 
337 #elif (TSURFACE == 5)
338 
339 ! ------ Correction of present monthly temperatures with LGM anomaly and
340 ! glacial index as well as elevation changes
341 
342  temp_diff = gamma_t*(zs(j,i)-zs_ref(j,i))
343 
344  do n=1, 12 ! month counter
345  temp_mm(j,i,n) = temp_mm_present(j,i,n) &
346  + glac_index*temp_mm_lgm_anom(j,i,n) &
347  + temp_diff
348  end do
349 
350 #endif
351 
352 ! ------ Mean annual air temperature
353 
354  temp_ma(j,i) = 0.0_dp ! initialisation value
355 
356  do n=1, 12 ! month counter
357  temp_ma(j,i) = temp_ma(j,i) + temp_mm(j,i,n)*inv_twelve
358  end do
359 
360 end do
361 end do
362 
363 !-------- Accumulation-ablation function as_perp --------
364 
365 #if (ELEV_DESERT == 1)
366 
367 gamma_p = gamma_p*1.0e-03_dp ! Precipitation lapse rate
368  ! for elevation desertification, in m^(-1)
369 zs_thresh = zs_thresh ! Elevation threshold, in m
370 
371 #endif
372 
373 #if (SOLID_PRECIP == 1) /* Marsiat (1994) */
374 
375 temp_rain = 7.0_dp ! Threshold monthly mean temperature for
376  ! precipitation = 100% rain, in deg C
377 temp_snow = -10.0_dp ! Threshold monthly mean temperature for &
378  ! precipitation = 100% snow, in deg C
379 
380 inv_delta_temp_rain_snow = 1.0_dp/(temp_rain-temp_snow)
381 
382 #elif (SOLID_PRECIP == 2) /* Bales et al. (2009) */
383 
384 temp_rain = 7.2_dp ! Threshold monthly mean temperature for &
385  ! precipitation = 100% rain, in deg C
386 temp_snow = -11.6_dp ! Threshold monthly mean temperature for &
387  ! precipitation = 100% snow, in deg C
388 
389 coeff(0) = 5.4714e-01_dp ! Coefficients
390 coeff(1) = -9.1603e-02_dp ! of
391 coeff(2) = -3.314e-03_dp ! the
392 coeff(3) = 4.66e-04_dp ! fifth-order
393 coeff(4) = 3.8e-05_dp ! polynomial
394 coeff(5) = 6.0e-07_dp ! fit
395 
396 #elif (SOLID_PRECIP == 3) /* Huybrechts and de Wolde (1999) */
397 
398 temp_rain = 2.0_dp ! Threshold instantaneous temperature for &
399  ! precipitation = 100% rain, in deg C
400 temp_snow = temp_rain ! Threshold instantaneous temperature for &
401  ! precipitation = 100% snow, in deg C
402 
403 s_stat = s_stat_0 ! Standard deviation of the air termperature
404  ! (same parameter as in the PDD model)
405 
406 inv_sqrt2_s_stat = 1.0_dp/(sqrt(2.0_dp)*s_stat)
407 
408 #endif
409 
410 #if (ABLSURFACE==1 || ABLSURFACE==2)
411 
412 s_stat = s_stat_0
413 beta1 = beta1_0 *(0.001_dp/86400.0_dp)*(rho_w/rho)
414  ! (mm WE)/(d*deg C) --> (m IE)/(s*deg C)
415 beta2 = beta2_0 *(0.001_dp/86400.0_dp)*(rho_w/rho)
416  ! (mm WE)/(d*deg C) --> (m IE)/(s*deg C)
417 
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 
569 ! ---- Formation rate of superimposed ice (melt_star), melt rate (melt)
570 ! and runoff rate (runoff)
571 
572 #if (ABLSURFACE==1)
573 
574  if ((beta1*et(j,i)) <= (pmax*snowfall(j,i))) then
575  melt_star(j,i) = beta1*et(j,i)
576  melt(j,i) = 0.0_dp
577  runoff(j,i) = melt(j,i)+rainfall(j,i)
578  else
579  melt_star(j,i) = pmax*snowfall(j,i)
580  melt(j,i) = beta2*(et(j,i)-melt_star(j,i)/beta1)
581  runoff(j,i) = melt(j,i)+rainfall(j,i)
582  end if
583 
584 #elif (ABLSURFACE==2)
585 
586  if ( rainfall(j,i) <= (pmax*snowfall(j,i)) ) then
587 
588  if ( (rainfall(j,i)+beta1*et(j,i)) <= (pmax*snowfall(j,i)) ) then
589  melt_star(j,i) = rainfall(j,i)+beta1*et(j,i)
590  melt(j,i) = 0.0_dp
591  runoff(j,i) = melt(j,i)
592  else
593  melt_star(j,i) = pmax*snowfall(j,i)
594  melt(j,i) = beta2 &
595  *(et(j,i)-(melt_star(j,i)-rainfall(j,i))/beta1)
596  runoff(j,i) = melt(j,i)
597  end if
598 
599  else
600 
601  melt_star(j,i) = pmax*snowfall(j,i)
602  melt(j,i) = beta2*et(j,i)
603  runoff(j,i) = melt(j,i) + rainfall(j,i)-pmax*snowfall(j,i)
604 
605  end if
606 
607 #endif
608 
609 #elif (ABLSURFACE==3)
610 
611  temp_jja_help = one_third*(temp_mm(j,i,6)+temp_mm(j,i,7)+temp_mm(j,i,8))
612 
613  melt_star(j,i) = 0.0_dp ! no superimposed ice considered
614  melt(j,i) = lambda_lti*max((temp_jja_help-temp_lti), 0.0_dp)
615  runoff(j,i) = melt(j,i) + rainfall(j,i)
616 
617 #endif
618 
619 ! ---- Evaporation
620 
621  evap(j,i) = 0.0_dp
622 
623 
624 ! ------ Accumulation minus runoff
625 
626  as_perp(j,i) = accum(j,i) - evap(j,i) - runoff(j,i)
627 
628  !hack for initial present-day conditions:
629  !set as_perp to zero, hold bedrock fixed and let zs evolve freely (initial zs smoothing)
630  !as_perp(j,i) = 0.0_dp
631 
632 ! ------ Ice-surface temperature (10-m firn temperature) temp_s,
633 ! including empirical firn-warming correction due to
634 ! refreezing meltwater when superimposed ice is formed
635 
636  if (melt_star(j,i).ge.melt(j,i)) then
637  temp_s(j,i) = temp_ma(j,i) &
638  +mu*(melt_star(j,i)-melt(j,i))
639  else
640  temp_s(j,i) = temp_ma(j,i)
641  end if
642 
643  if (temp_s(j,i) > -0.001_dp) temp_s(j,i) = -0.001_dp
644  ! Cut-off of positive air temperatures
645 
646 end do
647 end do
648 
649 end subroutine boundary
650 !
Generated on Mon May 20 2013 10:18:34 for SICOPOLIS V3.0 by   doxygen 1.8.1