calc_temp_r.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  c a l c _ t e m p _ r
!
!> @file
!!
!! Computation of temperature for an ice-free column.
!!
!! @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 temperature for an ice-free column.
!<------------------------------------------------------------------------------
subroutine calc_temp_r(atr1, alb1, i, j)

use sico_types
use sico_variables
use sico_vars
use sico_sle_solvers, only : tri_sle

#if (CALCMOD==2 || CALCMOD==3)
use enth_temp_omega, only : enth_fct_temp_omega
#endif

implicit none

integer(i4b), intent(in) :: i, j
real(dp), intent(in) :: atr1, alb1

integer(i4b) :: kc, kt, kr
real(dp) :: ctr1, clb1
real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
            lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
            lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
            lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
            lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
real(dp) :: enth_val

!-------- Abbreviations --------

ctr1 = atr1
clb1 = alb1*q_geo(j,i)

!-------- Set up the equations for the bedrock temperature --------

kr=0
lgs_a1(kr) = 1.0_dp
lgs_a2(kr) = -1.0_dp
lgs_b(kr)    = clb1

#if Q_LITHO==1
!   (coupled heat-conducting bedrock)

do kr=1, krmax-1
   lgs_a0(kr) = - ctr1
   lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
   lgs_a2(kr) = - ctr1
   lgs_b(kr)    = temp_r(kr,j,i)
end do

#elif Q_LITHO==0
!   (no coupled heat-conducting bedrock)

do kr=1, krmax-1
   lgs_a0(kr) = 1.0_dp
   lgs_a1(kr) = 0.0_dp
   lgs_a2(kr) = -1.0_dp
   lgs_b(kr)  = 2.0_dp*clb1
end do

#endif

kr=krmax
lgs_a0(kr) = 0.0_dp
lgs_a1(kr) = 1.0_dp
lgs_b(kr)   = temp_s(j,i)

!-------- Solve system of linear equations --------

call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)

!-------- Assign the result --------

do kr=0, krmax
   temp_r_neu(kr,j,i) = lgs_x(kr)
end do

!-------- Water content and age (enthalpy method: also enthalpy)
!                       in the non-existing lower (kt) ice layer --------

#if (CALCMOD==2 || CALCMOD==3)
enth_val = enth_fct_temp_omega(temp_s(j,i), 0.0_dp)
#endif

do kt=0, ktmax
   omega_t_neu(kt,j,i) = 0.0_dp
   age_t_neu(kt,j,i)   = 0.0_dp
#if (CALCMOD==2 || CALCMOD==3)
   enth_t_neu(kt,j,i)  = enth_val
#endif
end do

!-------- Temperature and age
!         (enthalpy method: also water content and enthalpy)
!                   in the non-existing upper (kc) ice layer --------

do kc=0, kcmax
   temp_c_neu(kc,j,i)  = temp_s(j,i)
   age_c_neu(kc,j,i)   = 0.0_dp
#if (CALCMOD==2 || CALCMOD==3)
   omega_c_neu(kc,j,i) = 0.0_dp
   enth_c_neu(kc,j,i)  = enth_val
#endif
end do

end subroutine calc_temp_r
!