boundary.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  b o u n d a r y
!
!> @file
!!
!! Computation of the surface temperature (must be less than 0 deg C!)
!! and of the accumulation-ablation function.
!!
!! @section Copyright
!!
!! Copyright 2009-2016 Ralf Greve
!!
!! @section License
!!
!! This file is part of SICOPOLIS.
!!
!! SICOPOLIS is free software: you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! SICOPOLIS is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with SICOPOLIS.  If not, see <http://www.gnu.org/licenses/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------------------------------------------------------------------------------
!> Computation of the surface temperature (must be less than 0 deg C!)
!! and of the accumulation-ablation function.
!<------------------------------------------------------------------------------
subroutine boundary(time, dtime, dxi, deta, &
                    delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)

use sico_types
use sico_variables
use sico_vars
#if (NETCDF > 1)
use netcdf
use nc_check
#endif

implicit none

real(dp), intent(in) :: time, dtime, dxi, deta

real(dp), intent(out)   :: delta_ts, glac_index, dzsl_dtau, z_mar
real(dp), intent(inout) :: z_sl

! Further return variables
! (defined as global variables in module sico_variables):
!
!    accum(j,i), evap(j,i), runoff(j,i), as_perp(j,i), temp_s(j,i)

integer(i4b) :: i, j, n
integer(i4b) :: i_gr, i_kl
integer(i4b) :: nrec_temp_precip
integer(i4b), save :: nrec_temp_precip_save = -1
real(dp) :: z_sl_old
real(dp) :: z_sl_min, t1, t2, t3, t4, t5, t6
real(dp) :: time_gr, time_kl
real(dp) :: z_sle_present, z_sle_help
real(dp), dimension(0:JMAX,0:IMAX,0:12) :: precip
real(dp), dimension(0:JMAX,0:IMAX) :: &
          snowfall, rainfall, melt, melt_star
real(dp), dimension(0:JMAX,0:IMAX,12) :: temp_mm
real(dp), dimension(0:JMAX,0:IMAX) :: temp_ma, temp_ampl
real(dp), dimension(0:JMAX,0:IMAX) :: temp_ma_anom, temp_mj_anom, precip_ma_anom
real(dp), dimension(0:IMAX,0:JMAX), save :: temp_ma_anom_tra, temp_mj_anom_tra, &
                                            precip_ma_anom_tra
                                    ! These arrays are transposed with respect
                                    ! to the usual SICOPOLIS convention; 
                                    ! they have i as first and j as second index
real(dp), dimension(12) :: temp_mm_help
real(dp) :: temp_jja_help
real(dp), dimension(0:JMAX,0:IMAX) :: et
real(dp) :: theta_ma, c_ma, kappa_ma, gamma_ma, &
            theta_mj, c_mj, kappa_mj, gamma_mj
real(dp) :: sine_factor
real(dp) :: gamma_p, zs_thresh, &
            temp_rain, temp_snow, &
            inv_delta_temp_rain_snow, coeff(0:5), inv_sqrt2_s_stat, &
            precip_fact, frac_solid
real(dp) :: s_stat, &
            phi_sep, temp_lt, temp_ht, inv_delta_temp_ht_lt, &
            beta1_lt, beta1_ht, beta2_lt, beta2_ht, &
            beta1, beta2, pmax, mu, lambda_lti, temp_lti
real(dp), dimension(256) :: pdd_mod_lat, delta_pdd_mod_lat_inv, &
                            pdd_mod_fac_w, pdd_mod_fac_e, pdd_mod_fac
real(dp) :: lon_w_e_sep, pdd_mod_fac_interpol
logical, dimension(0:JMAX,0:IMAX) :: check_point

#if (NETCDF > 1)
integer(i4b) :: ncv
!     ncv:       Variable ID
integer(i4b) :: nc3cor(3)
!     nc3cor(3): Corner of a 3-d array
integer(i4b) :: nc3cnt(3)
!     nc3cnt(3): Count of a 3-d array
#endif

real(dp), parameter :: &
          inv_twelve = 1.0_dp/12.0_dp, one_third = 1.0_dp/3.0_dp

character(len=64), parameter :: thisroutine = 'boundary'

!-------- Initialization of intent(out) variables --------

z_sl_old   = z_sl

delta_ts   = 0.0_dp
glac_index = 0.0_dp
z_sl       = 0.0_dp
dzsl_dtau  = 0.0_dp
z_mar      = 0.0_dp

!-------- Surface-temperature deviation from present values --------

#if TSURFACE==1
delta_ts = delta_ts0
!                           ! Steady state with prescribed constant
!                           ! air-temperature deviation
#elif TSURFACE==3
delta_ts = sine_amplit &
           *cos(2.0_dp*pi*time/(sine_period*year_sec)) &
           -sine_amplit
!                           ! Sinusoidal air-temperature forcing
#elif TSURFACE==4

!  ------ delta_ts from ice-core record

if (time/year_sec < real(grip_time_min,dp)) then
   delta_ts = griptemp(0)
else if (time/year_sec < real(grip_time_max,dp)) then

   i_kl = floor(((time/year_sec) &
          -real(grip_time_min,dp))/real(grip_time_stp,dp))
   i_kl = max(i_kl, 0)

   i_gr = ceiling(((time/year_sec) &
          -real(grip_time_min,dp))/real(grip_time_stp,dp))
   i_gr = min(i_gr, ndata_grip)

   if (i_kl == i_gr) then

      delta_ts = griptemp(i_kl)

   else

      time_kl = (grip_time_min + i_kl*grip_time_stp) *year_sec
      time_gr = (grip_time_min + i_gr*grip_time_stp) *year_sec

      delta_ts = griptemp(i_kl) &
                +(griptemp(i_gr)-griptemp(i_kl)) &
                *(time-time_kl)/(time_gr-time_kl)
                 ! linear interpolation of the ice-core data

   end if

else
   delta_ts  = griptemp(ndata_grip)
end if

delta_ts = delta_ts * grip_temp_fact
!                ! modification by constant factor

!-------- Glacial index --------

#elif TSURFACE==5

if (time/year_sec < real(gi_time_min,dp)) then
   glac_index = glacial_index(0)
else if (time/year_sec < real(gi_time_max,dp)) then

   i_kl = floor(((time/year_sec) &
          -real(gi_time_min,dp))/real(gi_time_stp,dp))
   i_kl = max(i_kl, 0)

   i_gr = ceiling(((time/year_sec) &
          -real(gi_time_min,dp))/real(gi_time_stp,dp))
   i_gr = min(i_gr, ndata_gi)

   if (i_kl == i_gr) then

      glac_index = glacial_index(i_kl)

   else

      time_kl = (gi_time_min + i_kl*gi_time_stp) *year_sec
      time_gr = (gi_time_min + i_gr*gi_time_stp) *year_sec

      glac_index = glacial_index(i_kl) &
                +(glacial_index(i_gr)-glacial_index(i_kl)) &
                *(time-time_kl)/(time_gr-time_kl)
                 ! linear interpolation of the glacial-index data

   end if

else
   glac_index  = glacial_index(ndata_gi)
end if

!-------- Reading of surface temperature and precipitation
!                                       directly from data --------

#elif ( (TSURFACE==6) && (ACCSURFACE==6) )

if (time/year_sec <= temp_precip_time_min) then
   nrec_temp_precip = 0
else if (time/year_sec < temp_precip_time_max) then
   nrec_temp_precip = nint( ((time/year_sec)-temp_precip_time_min) &
                            / temp_precip_time_stp )
else
   nrec_temp_precip = ndata_temp_precip
end if

if ( nrec_temp_precip < 0 ) then
   stop ' boundary: nrec_temp_precip < 0, not allowed!'
else if ( nrec_temp_precip > ndata_temp_precip ) then
   stop ' boundary: nrec_temp_precip > ndata_temp_precip, not allowed!'
end if

if ( nrec_temp_precip /= nrec_temp_precip_save ) then

   nc3cor(1) = 1
   nc3cor(2) = 1
   nc3cor(3) = nrec_temp_precip + 1
   nc3cnt(1) = imax + 1
   nc3cnt(2) = jmax + 1
   nc3cnt(3) = 1

   call check( nf90_inq_varid(ncid_temp_precip, 'annualtemp_anom', ncv), &
               thisroutine )
   call check( nf90_get_var(ncid_temp_precip, ncv, temp_ma_anom_tra, &
                            start=nc3cor, count=nc3cnt), thisroutine )

   call check( nf90_inq_varid(ncid_temp_precip, 'julytemp_anom', ncv), &
               thisroutine )
   call check( nf90_get_var(ncid_temp_precip, ncv, temp_mj_anom_tra, &
                            start=nc3cor, count=nc3cnt), thisroutine )

   call check( nf90_inq_varid(ncid_temp_precip, 'precipitation_anom', ncv), &
               thisroutine )
   call check( nf90_get_var(ncid_temp_precip, ncv, precip_ma_anom_tra, &
                            start=nc3cor, count=nc3cnt), thisroutine )

end if

temp_ma_anom   = transpose(temp_ma_anom_tra)
temp_mj_anom   = transpose(temp_mj_anom_tra)
precip_ma_anom = transpose(precip_ma_anom_tra) *(1.0e-03_dp/year_sec)*(rho_w/rho)
                                      ! mm/a water equiv. -> m/s ice equiv.

nrec_temp_precip_save = nrec_temp_precip

#endif

!-------- Sea level --------

#if SEA_LEVEL==1
!  ------ constant sea level
z_sl = z_sl0

#elif SEA_LEVEL==2
!  ------ saw-tooth-shaped palaeoclimatic sea-level forcing

z_sl_min = -130.0_dp

t1 = -250000.0_dp *year_sec
t2 = -140000.0_dp *year_sec
t3 = -125000.0_dp *year_sec
t4 =  -21000.0_dp *year_sec
t5 =   -8000.0_dp *year_sec
t6 =       0.0_dp *year_sec

if (time < t1) then
   z_sl = 0.0_dp
else if (time < t2) then
   z_sl = z_sl_min*(time-t1)/(t2-t1)
else if (time < t3) then
   z_sl = -z_sl_min*(time-t3)/(t3-t2)
else if (time < t4) then
   z_sl = z_sl_min*(time-t3)/(t4-t3)
else if (time < t5) then
   z_sl = -z_sl_min*(time-t5)/(t5-t4)
else if (time < t6) then
   z_sl = 0.0_dp
else
   z_sl = 0.0_dp
end if

#elif SEA_LEVEL==3

!  ------ z_sl from the SPECMAP record

if (time/year_sec < real(specmap_time_min,dp)) then
   z_sl = specmap_zsl(0)
else if (time/year_sec < real(specmap_time_max,dp)) then

   i_kl = floor(((time/year_sec) &
          -real(specmap_time_min,dp))/real(specmap_time_stp,dp))
   i_kl = max(i_kl, 0)

   i_gr = ceiling(((time/year_sec) &
          -real(specmap_time_min,dp))/real(specmap_time_stp,dp))
   i_gr = min(i_gr, ndata_specmap)

   if (i_kl == i_gr) then

      z_sl = specmap_zsl(i_kl)

   else

      time_kl = (specmap_time_min + i_kl*specmap_time_stp) *year_sec
      time_gr = (specmap_time_min + i_gr*specmap_time_stp) *year_sec

      z_sl = specmap_zsl(i_kl) &
                +(specmap_zsl(i_gr)-specmap_zsl(i_kl)) &
                *(time-time_kl)/(time_gr-time_kl)
                 ! linear interpolation of the sea-level data

   end if

else
   z_sl  = specmap_zsl(ndata_specmap)
end if

#endif

!  ------ Time derivative of the sea level

if ( z_sl_old > -999999.9_dp ) then
   dzsl_dtau = (z_sl-z_sl_old)/dtime
else   ! only dummy value for z_sl_old available
   dzsl_dtau = 0.0_dp
end if

!  ------ Minimum bedrock elevation for extent of marine ice

#if MARGIN==2

#if ( MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3 )
z_mar = z_mar
#elif ( MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 )
z_mar = fact_z_mar*z_sl
#elif ( MARINE_ICE_CALVING==6 || MARINE_ICE_CALVING==7 )
if (z_sl >= -80.0_dp) then
   z_mar = 2.5_dp*z_sl
else
   z_mar = 10.25_dp*(z_sl+80.0_dp)-200.0_dp
end if
z_mar = fact_z_mar*z_mar
#endif

#endif

!  ------ Update of the mask according to the sea level

!    ---- Check all sea and floating-ice points and their direct
!         neighbours

do i=0, imax
do j=0, jmax
   check_point(j,i) = .false.
end do
end do

do i=1, imax-1
do j=1, jmax-1
   if (maske(j,i) >= 2) then
      check_point(j  ,i  ) = .true.
      check_point(j  ,i+1) = .true.
      check_point(j  ,i-1) = .true.
      check_point(j+1,i  ) = .true.
      check_point(j-1,i  ) = .true.
   end if
end do
end do

do i=1, imax-1
do j=1, jmax-1
   if (check_point(j,i)) then
      maske_neu(j,i) = mask_update(z_sl, i, j)
   end if
end do
end do

!    ---- Assign new values of the mask

do i=1, imax-1
do j=1, jmax-1
   if (check_point(j,i)) then
      maske(j,i) = maske_neu(j,i)
   end if
end do
end do

!-------- Surface air temperatures --------

#if TEMP_PRESENT_PARA == 1   /* Parameterization by Ritz et al. (1997) */

theta_ma = 49.13_dp
gamma_ma = -7.992e-03_dp
c_ma     = -0.7576_dp
kappa_ma =  0.0_dp

theta_mj = 30.38_dp
gamma_mj = -6.277e-03_dp
c_mj     = -0.3262_dp
kappa_mj =  0.0_dp

#elif TEMP_PRESENT_PARA == 2   /* Parameterization by Fausto et al. (2009) */

theta_ma = 41.83_dp
gamma_ma = -6.309e-03_dp
c_ma     = -0.7189_dp
kappa_ma = -0.0672_dp

theta_mj = 14.70_dp
gamma_mj = -5.426e-03_dp
c_mj     = -0.1585_dp
kappa_mj = -0.0518_dp

#else

stop ' boundary: Parameter TEMP_PRESENT_PARA must be either 1 or 2!'

#endif

do i=0, imax
do j=0, jmax

!  ------ Present-day mean-annual air temperature

   temp_ma_present(j,i) = theta_ma &
                  + gamma_ma*zs(j,i) &
                  + c_ma*phi(j,i)*pi_180_inv &
                  + kappa_ma*(modulo(lambda(j,i)+pi,2.0_dp*pi)-pi)*pi_180_inv
                              ! west longitudes counted negatively

!  ------ Present-day mean-July (summer) air temperature

   temp_mj_present(j,i) = theta_mj &
                    + gamma_mj*zs(j,i) &
                    + c_mj*phi(j,i)*pi_180_inv &
                    + kappa_mj*(modulo(lambda(j,i)+pi,2.0_dp*pi)-pi)*pi_180_inv
                                ! west longitudes counted negatively

#if (TSURFACE <= 4)

!  ------ Correction with deviation delta_ts

   temp_ma(j,i)   = temp_ma_present(j,i) + delta_ts
   temp_mm(j,i,7) = temp_mj_present(j,i) + delta_ts

#elif (TSURFACE == 5)

!  ------ Correction with LGM anomaly and glacial index

   temp_ma(j,i)   = temp_ma_present(j,i) + glac_index*temp_ma_lgm_anom(j,i)
   temp_mm(j,i,7) = temp_mj_present(j,i) + glac_index*temp_mj_lgm_anom(j,i)

#elif (TSURFACE == 6)

!  ------ Mean-annual and mean-July (summer) surface temperatures
!                                            from data read above --------

   temp_ma(j,i)   = temp_ma_present(j,i) + temp_ma_anom_fact*temp_ma_anom(j,i)
   temp_mm(j,i,7) = temp_mj_present(j,i) + temp_mj_anom_fact*temp_mj_anom(j,i)

#endif

!  ------ Amplitude of the annual cycle

   temp_ampl(j,i) = temp_mm(j,i,7) - temp_ma(j,i)

   if (temp_ampl(j,i) < eps) then
      temp_ampl(j,i) = eps   ! Correction of amplitude, if required
   end if

end do
end do

!  ------ Monthly temperatures

do n=1, 12   ! month counter

   sine_factor = sin((real(n,dp)-4.0_dp)*pi/6.0_dp)

   do i=0, imax
   do j=0, jmax
      temp_mm(j,i,n) = temp_ma(j,i) + sine_factor*temp_ampl(j,i)
   end do
   end do

end do

!-------- Accumulation-ablation function as_perp --------

#if (ELEV_DESERT == 1)

gamma_p   = gamma_p*1.0e-03_dp   ! Precipitation lapse rate
                                 ! for elevation desertification, in m^(-1)
zs_thresh = zs_thresh            ! Elevation threshold, in m

#endif

#if (SOLID_PRECIP == 1)     /* Marsiat (1994) */

temp_rain =    7.0_dp   ! Threshold monthly mean temperature for
                        ! precipitation = 100% rain, in deg C
temp_snow =  -10.0_dp   ! Threshold monthly mean temperature for &
                        ! precipitation = 100% snow, in deg C

inv_delta_temp_rain_snow = 1.0_dp/(temp_rain-temp_snow)

#elif (SOLID_PRECIP == 2)   /* Bales et al. (2009) */

temp_rain =    7.2_dp   ! Threshold monthly mean temperature for &
                        ! precipitation = 100% rain, in deg C
temp_snow =  -11.6_dp   ! Threshold monthly mean temperature for &
                        ! precipitation = 100% snow, in deg C

coeff(0) =  5.4714e-01_dp   ! Coefficients
coeff(1) = -9.1603e-02_dp   ! of
coeff(2) = -3.314e-03_dp    ! the
coeff(3) =  4.66e-04_dp     ! fifth-order
coeff(4) =  3.8e-05_dp      ! polynomial
coeff(5) =  6.0e-07_dp      ! fit

#elif (SOLID_PRECIP == 3)   /* Huybrechts and de Wolde (1999) */

temp_rain = 2.0_dp      ! Threshold instantaneous temperature for &
                        ! precipitation = 100% rain, in deg C
temp_snow = temp_rain   ! Threshold instantaneous temperature for &
                        ! precipitation = 100% snow, in deg C

s_stat    = s_stat_0    ! Standard deviation of the air termperature
                        ! (same parameter as in the PDD model)

inv_sqrt2_s_stat = 1.0_dp/(sqrt(2.0_dp)*s_stat)

#endif

#if (ABLSURFACE==1 || ABLSURFACE==2)

s_stat   = s_stat_0

phi_sep  = phi_sep_0*pi_180   ! separates different domains for computation of
                              ! degree-day factors beta1 and beta2
                              ! deg N --> rad
temp_lt  = -1.0_dp
temp_ht  = 10.0_dp
inv_delta_temp_ht_lt = 1.0_dp/(temp_ht-temp_lt)

beta1_lt = beta1_lt_0  *(0.001_dp/86400.0_dp)*(rho_w/rho)
                           ! (mm WE)/(d*deg C) --> (m IE)/(s*deg C)
beta1_ht = beta1_ht_0  *(0.001_dp/86400.0_dp)*(rho_w/rho)
                           ! (mm WE)/(d*deg C) --> (m IE)/(s*deg C)
beta2_lt = beta2_lt_0  *(0.001_dp/86400.0_dp)*(rho_w/rho)
                           ! (mm WE)/(d*deg C) --> (m IE)/(s*deg C)
beta2_ht = beta2_ht_0  *(0.001_dp/86400.0_dp)*(rho_w/rho)
                           ! (mm WE)/(d*deg C) --> (m IE)/(s*deg C)
pmax     = pmax_0
mu       = mu_0        *(1000.0_dp*86400.0_dp)*(rho/rho_w)
                           ! (d*deg C)/(mm WE) --> (s*deg C)/(m IE)
#if (PDD_MODIFIER==2)

lon_w_e_sep = lon_w_e_sep *pi_180   ! deg N -> rad

if (lon_w_e_sep < 0.0_dp) then
   lon_w_e_sep = lon_w_e_sep + 2.0_dp*pi
else if (lon_w_e_sep >= (2.0_dp*pi)) then
   lon_w_e_sep = lon_w_e_sep - 2.0_dp*pi
end if

pdd_mod_lat   = 0.0_dp   ! initialisation
pdd_mod_fac_w = 0.0_dp   ! initialisation
pdd_mod_fac_e = 0.0_dp   ! initialisation

pdd_mod_lat(   1) = pdd_mod_lat_01 *pi_180   ! deg N -> rad
pdd_mod_fac_w( 1) = pdd_mod_fac_w_01
pdd_mod_fac_e( 1) = pdd_mod_fac_e_01
pdd_mod_lat(   2) = pdd_mod_lat_02 *pi_180   ! deg N -> rad
pdd_mod_fac_w( 2) = pdd_mod_fac_w_02
pdd_mod_fac_e( 2) = pdd_mod_fac_e_02
pdd_mod_lat(   3) = pdd_mod_lat_03 *pi_180   ! deg N -> rad
pdd_mod_fac_w( 3) = pdd_mod_fac_w_03
pdd_mod_fac_e( 3) = pdd_mod_fac_e_03
pdd_mod_lat(   4) = pdd_mod_lat_04 *pi_180   ! deg N -> rad
pdd_mod_fac_w( 4) = pdd_mod_fac_w_04
pdd_mod_fac_e( 4) = pdd_mod_fac_e_04
pdd_mod_lat(   5) = pdd_mod_lat_05 *pi_180   ! deg N -> rad
pdd_mod_fac_w( 5) = pdd_mod_fac_w_05
pdd_mod_fac_e( 5) = pdd_mod_fac_e_05
pdd_mod_lat(   6) = pdd_mod_lat_06 *pi_180   ! deg N -> rad
pdd_mod_fac_w( 6) = pdd_mod_fac_w_06
pdd_mod_fac_e( 6) = pdd_mod_fac_e_06
pdd_mod_lat(   7) = pdd_mod_lat_07 *pi_180   ! deg N -> rad
pdd_mod_fac_w( 7) = pdd_mod_fac_w_07
pdd_mod_fac_e( 7) = pdd_mod_fac_e_07
pdd_mod_lat(   8) = pdd_mod_lat_08 *pi_180   ! deg N -> rad
pdd_mod_fac_w( 8) = pdd_mod_fac_w_08
pdd_mod_fac_e( 8) = pdd_mod_fac_e_08
pdd_mod_lat(   9) = pdd_mod_lat_09 *pi_180   ! deg N -> rad
pdd_mod_fac_w( 9) = pdd_mod_fac_w_09
pdd_mod_fac_e( 9) = pdd_mod_fac_e_09
pdd_mod_lat(  10) = pdd_mod_lat_10 *pi_180   ! deg N -> rad
pdd_mod_fac_w(10) = pdd_mod_fac_w_10
pdd_mod_fac_e(10) = pdd_mod_fac_e_10

delta_pdd_mod_lat_inv = 0.0_dp   ! initialisation

do n=1, n_pdd_mod-1
   delta_pdd_mod_lat_inv(n) = 1.0_dp/(pdd_mod_lat(n+1)-pdd_mod_lat(n))
end do

#endif

#elif (ABLSURFACE==3)

lambda_lti = lambda_lti  *(0.001_dp/year_sec)*(rho_w/rho)
                         ! (mm WE)/(a*deg C) --> (m IE)/(s*deg C)
temp_lti   = temp_lti
mu         = 0.0_dp      ! no superimposed ice considered

#endif

do i=0, imax
do j=0, jmax

#if (ABLSURFACE==1 || ABLSURFACE==2)

   if (phi(j,i) <= phi_sep) then

         beta1 = beta1_ht
         beta2 = beta2_ht

   else

      if (temp_mm(j,i,7) >= temp_ht) then
         beta1 = beta1_ht
         beta2 = beta2_ht
      else if (temp_mm(j,i,7) <= temp_lt) then
         beta1 = beta1_lt
         beta2 = beta2_lt
      else
         beta1 = beta1_lt &
                 + (beta1_ht-beta1_lt) &
                   *inv_delta_temp_ht_lt*(temp_mm(j,i,7)-temp_lt)
         beta2 = beta2_ht &
                 + (beta2_lt-beta2_ht) &
                   *(inv_delta_temp_ht_lt*(temp_ht-temp_mm(j,i,7)))**3
      end if

   end if

#if (PDD_MODIFIER==2)

   if ( lambda(j,i) <= lon_w_e_sep ) then   ! West Greenland
      pdd_mod_fac = pdd_mod_fac_w
   else          ! lambda(j,i) > lon_w_e_sep, East Greenland
      pdd_mod_fac = pdd_mod_fac_e
   end if

   if (phi(j,i) <= pdd_mod_lat(1)) then

      beta1 = beta1 * pdd_mod_fac(1)
      beta2 = beta2 * pdd_mod_fac(1)

   else if (phi(j,i) >= pdd_mod_lat(n_pdd_mod)) then

      beta1 = beta1 * pdd_mod_fac(n_pdd_mod)
      beta2 = beta2 * pdd_mod_fac(n_pdd_mod)

   else

      do n=1, n_pdd_mod-1

         if ( (phi(j,i) >= pdd_mod_lat(n)) &
              .and. &
              (phi(j,i) <= pdd_mod_lat(n+1)) ) then

            pdd_mod_fac_interpol = pdd_mod_fac(n) &
                                   + (pdd_mod_fac(n+1)-pdd_mod_fac(n)) &
                                     *(phi(j,i)-pdd_mod_lat(n)) &
                                     *delta_pdd_mod_lat_inv(n)

            beta1 = beta1 * pdd_mod_fac_interpol
            beta2 = beta2 * pdd_mod_fac_interpol

            exit

         end if

      end do

   end if

#endif

#endif

!  ------ Accumulation

#if (ACCSURFACE <= 3)

!    ---- Elevation desertification of precipitation

#if (ELEV_DESERT == 0)

   precip_fact = 1.0_dp   ! no elevation desertification

#elif (ELEV_DESERT == 1)

   if (zs_ref(j,i) < zs_thresh) then
      precip_fact &
         = exp(gamma_p*(max(zs(j,i),zs_thresh)-zs_thresh))
   else
      precip_fact &
         = exp(gamma_p*(max(zs(j,i),zs_thresh)-zs_ref(j,i)))
   end if

#else
   stop ' boundary: Parameter ELEV_DESERT must be either 0 or 1!'
#endif

   do n=1, 12   ! month counter
      precip(j,i,n) = precip_present(j,i,n)*precip_fact
   end do

#endif

!    ---- Precipitation change related to changing climate

#if ACCSURFACE==1
   precip_fact = accfact
#elif ACCSURFACE==2
   precip_fact = 1.0_dp + gamma_s*delta_ts
#elif ACCSURFACE==3
   precip_fact = exp(gamma_s*delta_ts)
#endif

#if (ACCSURFACE <= 3)

   precip(j,i,0) = 0.0_dp   ! initialisation value for mean annual precip

   do n=1, 12   ! month counter
      precip(j,i,n) = precip(j,i,n)*precip_fact   ! monthly precip
      precip(j,i,0) = precip(j,i,0) + precip(j,i,n)*inv_twelve
                                                  ! mean annual precip
   end do

#elif (ACCSURFACE == 5)

   precip(j,i,0) = 0.0_dp   ! initialisation value for mean annual precip

   do n=1, 12   ! month counter

#if (PRECIP_ANOM_INTERPOL==1)
      precip_fact = 1.0_dp-glac_index+glac_index*precip_lgm_anom(j,i,n)
                    ! interpolation with a linear function
#elif (PRECIP_ANOM_INTERPOL==2)
      precip_fact = exp(-glac_index*gamma_precip_lgm_anom(j,i,n))
                    ! interpolation with an exponential function
#endif

      precip(j,i,n) = precip_present(j,i,n)*precip_fact   ! monthly precip
      precip(j,i,0) = precip(j,i,0) + precip(j,i,n)*inv_twelve
                                                      ! mean annual precip
   end do

#elif (ACCSURFACE == 6)

!      -- Mean-annual precipitation from data already read above

   precip(j,i,0) = precip_ma_present(j,i) + precip_anom_fact*precip_ma_anom(j,i)
                                      ! mean annual precip

   do n=1, 12   ! month counter
      precip(j,i,n) = precip(j,i,0)   ! monthly precip, assumed to be equal
                                      ! to the mean annual precip
   end do

#endif

!    ---- Annual accumulation, snowfall and rainfall rates

   accum(j,i) = precip(j,i,0)

   snowfall(j,i) = 0.0_dp   ! initialisation value

   do n=1, 12   ! month counter

#if (SOLID_PRECIP == 1)     /* Marsiat (1994) */

      if (temp_mm(j,i,n) >= temp_rain) then
         frac_solid = 0.0_dp
      else if (temp_mm(j,i,n) <= temp_snow) then
         frac_solid = 1.0_dp
      else
         frac_solid = (temp_rain-temp_mm(j,i,n))*inv_delta_temp_rain_snow
      end if

#elif (SOLID_PRECIP == 2)   /* Bales et al. (2009) */

      if (temp_mm(j,i,n) >= temp_rain) then
         frac_solid = 0.0_dp
      else if (temp_mm(j,i,n) <= temp_snow) then
         frac_solid = 1.0_dp
      else
         frac_solid = coeff(0) + temp_mm(j,i,n) * ( coeff(1) &
                               + temp_mm(j,i,n) * ( coeff(2) &
                               + temp_mm(j,i,n) * ( coeff(3) &
                               + temp_mm(j,i,n) * ( coeff(4) &
                               + temp_mm(j,i,n) *   coeff(5) ) ) ) )
               ! evaluation of 5th-order polynomial by Horner scheme
      end if

#elif (SOLID_PRECIP == 3)   /* Huybrechts and de Wolde (1999) */

      frac_solid = 1.0_dp &
                   - 0.5_dp*erfcc((temp_rain-temp_mm(j,i,n))*inv_sqrt2_s_stat)

#endif

      snowfall(j,i) = snowfall(j,i) + precip(j,i,n)*frac_solid*inv_twelve

   end do

   rainfall(j,i) = precip(j,i,0) - snowfall(j,i)

   if (snowfall(j,i) < 0.0_dp) snowfall(j,i) = 0.0_dp   ! correction of
   if (rainfall(j,i) < 0.0_dp) rainfall(j,i) = 0.0_dp   ! negative values

!  ------ Ablation

!    ---- Runoff

#if (ABLSURFACE==1 || ABLSURFACE==2)

!      -- Temperature excess ET

   do n=1, 12   ! month counter
      temp_mm_help(n) = temp_mm(j,i,n)
   end do

   call pdd(temp_mm_help, s_stat, et(j,i))

!      -- Formation rate of superimposed ice (melt_star), melt rate (melt)
!         and runoff rate (runoff)

#if (ABLSURFACE==1)

   if ((beta1*et(j,i)) <= (pmax*snowfall(j,i))) then
      melt_star(j,i) = beta1*et(j,i)
      melt(j,i)      = 0.0_dp
      runoff(j,i)    = melt(j,i)+rainfall(j,i)
   else
      melt_star(j,i) = pmax*snowfall(j,i)
      melt(j,i)      = beta2*(et(j,i)-melt_star(j,i)/beta1)
      runoff(j,i)    = melt(j,i)+rainfall(j,i)
   end if

#elif (ABLSURFACE==2)

   if ( rainfall(j,i) <= (pmax*snowfall(j,i)) ) then

      if ( (rainfall(j,i)+beta1*et(j,i)) <= (pmax*snowfall(j,i)) ) then
         melt_star(j,i) = rainfall(j,i)+beta1*et(j,i)
         melt(j,i)      = 0.0_dp
         runoff(j,i)    = melt(j,i)
      else
         melt_star(j,i) = pmax*snowfall(j,i)
         melt(j,i)      = beta2 &
                          *(et(j,i)-(melt_star(j,i)-rainfall(j,i))/beta1)
         runoff(j,i)    = melt(j,i)
      end if

   else

      melt_star(j,i) = pmax*snowfall(j,i)
      melt(j,i)      = beta2*et(j,i)
      runoff(j,i)    = melt(j,i) + rainfall(j,i)-pmax*snowfall(j,i)

   end if

#endif

#elif (ABLSURFACE==3)

   temp_jja_help  = one_third*(temp_mm(j,i,6)+temp_mm(j,i,7)+temp_mm(j,i,8))

   melt_star(j,i) = 0.0_dp   ! no superimposed ice considered
   melt(j,i)      = lambda_lti*max((temp_jja_help-temp_lti), 0.0_dp)
   runoff(j,i)    = melt(j,i) + rainfall(j,i)

#endif

!    ---- Evaporation

   evap(j,i) = 0.0_dp

!  ------ Accumulation minus ablation

   as_perp(j,i) = accum(j,i) - evap(j,i) - runoff(j,i)

!  ------ Ice-surface temperature (10-m firn temperature) temp_s,
!         including empirical firn-warming correction due to
!         refreezing meltwater when superimposed ice is formed

   if (melt_star(j,i) >= melt(j,i)) then
      temp_s(j,i) = temp_ma(j,i) &
                    +mu*(melt_star(j,i)-melt(j,i))
   else
      temp_s(j,i) = temp_ma(j,i)
   end if

   if (temp_s(j,i) > -0.001_dp) temp_s(j,i) = -0.001_dp
                               ! Cut-off of positive air temperatures

end do
end do

end subroutine boundary
!