boundary_m.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Module :  b o u n d a r y _ m
!
!> @file
!!
!! Computation of the surface temperature (must be less than 0 deg C!)
!! and of the accumulation-ablation function.
!!
!! @section Copyright
!!
!! Copyright 2009-2019 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.
!<------------------------------------------------------------------------------
module boundary_m

  use sico_types_m
  use sico_variables_m
  use sico_vars_m
  use error_m

  implicit none

#if (INITMIP_CONST_SMB==1)  /* SMB held constant at its initial distribution */
  real(dp), dimension(0:JMAX,0:IMAX), save :: accum_save, evap_save, runoff_save
#endif

  public

contains

!-------------------------------------------------------------------------------
!> Main routine of boundary_m:
!! 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)

#if defined(ALLOW_OPENAD) /* OpenAD */
  use ctrl_m, only: myceiling, myfloor
  use sico_maths_m, only: my_erfc
#endif /* OpenAD */

#if (NETCDF > 1)
  use netcdf
  use nc_check_m
#endif

#if ((MARGIN==2) \
      && (marine_ice_formation==2) \
      && (marine_ice_calving==9))
  use calving_underwater_ice_m
#endif

#if (DISC > 0)
  use discharge_workers_m, only: discharge, dis_perp, dt_glann
#endif

  use mask_update_sea_level_m
  use pdd_m

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_m):
!
!    accum(j,i), evap(j,i), runoff(j,i), as_perp(j,i),
!    calv_grounded(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) :: year_sec_inv, time_in_years
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,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) :: theta_ma, theta_ma_offset, c_ma, kappa_ma, gamma_ma, &
            theta_mj, theta_mj_offset, 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
logical, dimension(0:JMAX,0:IMAX) :: check_point
logical, save                     :: firstcall = .true.

#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'

#if defined(ALLOW_OPENAD) /* OpenAD */
integer(i4b) :: i_time_in_years
real(dp)     :: temp_val
#endif /* OpenAD */

year_sec_inv  = 1.0_dp/year_sec
time_in_years = time*year_sec_inv

!-------- 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_in_years/sine_period) &
           -sine_amplit
!                           ! Sinusoidal air-temperature forcing
#elif (TSURFACE==4)

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

if (time_in_years < real(grip_time_min,dp)) then
   delta_ts = griptemp(0)
else if (time_in_years < real(grip_time_max,dp)) then

#if !defined(ALLOW_OPENAD) /* Normal */
   i_kl = floor((time_in_years &
          -real(grip_time_min,dp))/real(grip_time_stp,dp))
#else /* OpenAD */
   call myfloor((time_in_years &
                -real(grip_time_min,dp))/real(grip_time_stp,dp),i_kl)
#endif /* Normal vs. OpenAD */
   i_kl = max(i_kl, 0)

#if !defined(ALLOW_OPENAD) /* Normal */
   i_gr = ceiling((time_in_years &
          -real(grip_time_min,dp))/real(grip_time_stp,dp))
#else /* OpenAD */
   call myceiling(((time_in_years &
          -real(grip_time_min,dp))/real(grip_time_stp,dp)),i_gr)
#endif /* Normal vs. OpenAD */
   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_in_years < real(gi_time_min,dp)) then
   glac_index = glacial_index(0)
else if (time_in_years < real(gi_time_max,dp)) then

#if !defined(ALLOW_OPENAD)/* Normal */
   i_kl = floor((time_in_years &
          -real(gi_time_min,dp))/real(gi_time_stp,dp))
#else /* OpenAD */
   call myfloor(((time_in_years &
          -real(gi_time_min,dp))/real(gi_time_stp,dp)), i_kl)
#endif /* Normal vs. OpenAD */
   i_kl = max(i_kl, 0)

#if !defined(ALLOW_OPENAD)/* Normal */
   i_gr = ceiling((time_in_years &
          -real(gi_time_min,dp))/real(gi_time_stp,dp))
#else /* OpenAD */
   call myceiling(((time_in_years &
          -real(gi_time_min,dp))/real(gi_time_stp,dp)), i_gr)
#endif /* Normal vs. OpenAD */
   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_in_years <= temp_precip_time_min) then
   nrec_temp_precip = 0
else if (time_in_years < temp_precip_time_max) then
   nrec_temp_precip = nint( (time_in_years-temp_precip_time_min) &
                            / temp_precip_time_stp )
else
   nrec_temp_precip = ndata_temp_precip
end if

if ( nrec_temp_precip < 0 ) then
   errormsg = ' >>> boundary: nrec_temp_precip < 0, not allowed!'
   call error(errormsg)
else if ( nrec_temp_precip > ndata_temp_precip ) then
   errormsg = ' >>> boundary: ' &
            //'nrec_temp_precip > ndata_temp_precip, not allowed!'
   call error(errormsg)
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_inv)*(rho_w/rho)
                       ! mm/a water equiv. -> m/s ice equiv.

nrec_temp_precip_save = nrec_temp_precip

#endif

!-------- Global annual temperature anomaly
!         (for parameterising the sub-ocean temperature anomaly
!                             for the ice discharge parameterisation) --------

#if (DISC==1)

dt_glann = 0.0_dp   ! anomaly set to zero

#elif (DISC==2)

!  ------ dT_glann from CLIMBER simulations

if (time_in_years < real(glann_time_min,dp)) then
   dt_glann = dt_glann_climber(0)
else if (time_in_years < real(glann_time_max,dp)) then

#if !defined(ALLOW_OPENAD) /* Normal */
   i_kl = floor((time_in_years &
          -real(glann_time_min,dp))/real(glann_time_stp,dp))
   i_kl = max(i_kl, 0)

   i_gr = ceiling((time_in_years &
          -real(glann_time_min,dp))/real(glann_time_stp,dp))
   i_gr = min(i_gr, ndata_glann)
#else /* OpenAD */
   call myfloor(((time_in_years &
          -real(glann_time_min,dp))/real(glann_time_stp,dp)), i_kl)
   i_kl = max(i_kl, 0)

   call myceiling(((time_in_years &
          -real(glann_time_min,dp))/real(glann_time_stp,dp)), i_gr)
   i_gr = min(i_gr, ndata_glann)
#endif /* Normal vs. OpenAD */

   if (i_kl == i_gr) then

      dt_glann = dt_glann_climber(i_kl)

   else

      time_kl = (glann_time_min + i_kl*glann_time_stp) *year_sec
      time_gr = (glann_time_min + i_gr*glann_time_stp) *year_sec

      dt_glann = dt_glann_climber(i_kl) &
                +(dt_glann_climber(i_gr)-dt_glann_climber(i_kl)) &
                *(time-time_kl)/(time_gr-time_kl)
                 ! linear interpolation of the data

   end if

else
   dt_glann  = dt_glann_climber(ndata_glann)
end if

#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_in_years < real(specmap_time_min,dp)) then
   z_sl = specmap_zsl(0)
else if (time_in_years < real(specmap_time_max,dp)) then

#if !defined(ALLOW_OPENAD) /* Normal */
   i_kl = floor((time_in_years &
          -real(specmap_time_min,dp))/real(specmap_time_stp,dp))
   i_kl = max(i_kl, 0)

   i_gr = ceiling((time_in_years &
          -real(specmap_time_min,dp))/real(specmap_time_stp,dp))
   i_gr = min(i_gr, ndata_specmap)
#else /* OpenAD */
   call myfloor(((time_in_years &
          -real(specmap_time_min,dp))/real(specmap_time_stp,dp)), i_kl)
   i_kl = max(i_kl, 0)

   call myceiling(((time_in_years &
          -real(specmap_time_min,dp))/real(specmap_time_stp,dp)), i_gr)
   i_gr = min(i_gr, ndata_specmap)
#endif /* Normal vs. OpenAD */

   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_i1b) 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_sea_level(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) */

#if (defined(TEMP_PRESENT_OFFSET))
theta_ma_offset = temp_present_offset
theta_mj_offset = temp_present_offset
#else
theta_ma_offset = 0.0_dp
theta_mj_offset = 0.0_dp
#endif

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

theta_mj = 30.38_dp + theta_mj_offset
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) */

#if (defined(TEMP_PRESENT_OFFSET))
theta_ma_offset = temp_present_offset
theta_mj_offset = temp_present_offset
#else
theta_ma_offset = 0.0_dp
theta_mj_offset = 0.0_dp
#endif

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

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

#else

errormsg = ' >>> boundary: Parameter TEMP_PRESENT_PARA must be either 1 or 2!'
call error(errormsg)

#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 &
#if !defined(ALLOW_OPENAD) /* Normal */
                  + kappa_ma*(modulo(lambda(j,i)+pi,2.0_dp*pi)-pi)*pi_180_inv
                              ! west longitudes counted negatively
#else /* OpenAD */
                  + kappa_ma*((lambda(j,i)+pi - &
                              (2.0_dp*pi * int((lambda(j,i)+pi)/(2.0_dp*pi)))) &
                              - pi) &
                            *pi_180_inv
#endif /* Normal vs. OpenAD */

!  ------ 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 &
#if !defined(ALLOW_OPENAD) /* Normal */
                    + kappa_mj*(modulo(lambda(j,i)+pi,2.0_dp*pi)-pi)*pi_180_inv
                                ! west longitudes counted negatively
#else /* OpenAD */
                    + kappa_mj*((lambda(j,i)+pi - &
                                (2.0_dp*pi * int((lambda(j,i)+pi)/(2.0_dp*pi)))) &
                                - pi) &
                              *pi_180_inv
#endif /* Normal vs. OpenAD */

#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 (INITMIP_CONST_SMB==1)  /* SMB held constant at its initial distribution */
if (firstcall) then
#endif

#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)

#elif (ABLSURFACE==3)

lambda_lti = lambda_lti  *(0.001_dp*year_sec_inv)*(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

#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
   errormsg = ' >>> boundary: Parameter ELEV_DESERT must be either 0 or 1!'
   call error(errormsg)
#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)

!      -- Precipitation from data already read above

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

   do n=1, 12   ! month counter
      precip(j,i,n) = precip_present(j,i,n) &
                      + precip_anom_fact*precip_ma_anom(j,i)
                        ! monthly precip (anomaly assumed to be equal
                        !                to the mean annual precip anomaly)
      precip(j,i,0) = precip(j,i,0) + precip(j,i,n)*inv_twelve
                        ! 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) */

#if !defined(ALLOW_OPENAD) /* Normal */
      frac_solid = 1.0_dp &
                   - 0.5_dp*erfc((temp_rain-temp_mm(j,i,n))*inv_sqrt2_s_stat)
#else /* OpenAD */
      call my_erfc((temp_rain-temp_mm(j,i,n))*inv_sqrt2_s_stat, temp_val)
      frac_solid = 1.0_dp - 0.5_dp*temp_val
#endif /* Normal vs. OpenAD */

#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

end do
end do

#if (INITMIP_CONST_SMB==1)  /* SMB held constant at its initial distribution */

accum_save  = accum
evap_save   = evap
runoff_save = runoff

else   ! firstcall == .false.

accum  = accum_save
evap   = evap_save
runoff = runoff_save

end if   ! firstcall == .true./.false.

#endif

!  ------ SMB = accumulation minus ablation

as_perp = accum - evap - runoff

!    ---- Prescribed SMB correction

smb_corr_prescribed = smb_corr_in

#if (defined(INITMIP_SMB_ANOM_FILE))   /* Correction for ISMIP InitMIP */

#if defined(ALLOW_OPENAD) /* OpenAD */
call myfloor(time_in_years, i_time_in_years)
#endif /* OpenAD */

if ((time_in_years > 0.0_dp).and.(time_in_years <= 40.0_dp)) then
#if !defined(ALLOW_OPENAD) /* Normal */
   smb_corr_prescribed = smb_corr_prescribed &
                              + 0.025_dp*floor(time_in_years) * smb_anom_initmip
#else /* OpenAD */
   smb_corr_prescribed = smb_corr_prescribed &
                              + 0.025_dp*(i_time_in_years) * smb_anom_initmip
#endif /* Normal vs. OpenAD */
else if (time_in_years > 40.0_dp) then
   smb_corr_prescribed = smb_corr_prescribed &
                              + smb_anom_initmip
end if

#endif

as_perp = as_perp + smb_corr_prescribed

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

#if !defined(ALLOW_OPENAD) /* Normal */

where (melt_star >= melt)
   temp_s = temp_ma + mu*(melt_star-melt)
elsewhere
   temp_s = temp_ma
end where

where (temp_s > -0.001_dp) temp_s = -0.001_dp
                            ! Cut-off of positive air temperatures

#else /* OpenAD */

do i=0, imax
do j=0, jmax
   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
end do
end do

do i=0, imax
do j=0, jmax
   if (temp_s(j,i) > -0.001_dp) then
      temp_s(j,i) = -0.001_dp   ! Cut-off of positive air temperatures
   end if
end do
end do

#endif /* Normal vs. OpenAD */

!-------- Calving rate of grounded ice --------

calv_grounded = 0.0_dp

#if ((MARGIN==2) \
      && (marine_ice_formation==2) \
      && (marine_ice_calving==9))

call calving_underwater_ice(z_sl)
calv_grounded = calv_grounded + calv_uw_ice

#endif

#if (DISC>0)   /* Ice discharge parameterization */

call discharge(dxi, deta)
calv_grounded = calv_grounded + dis_perp

#endif

if (firstcall) firstcall = .false.

end subroutine boundary

!-------------------------------------------------------------------------------

end module boundary_m
!