shift_cts_downward.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  s h i f t _ c t s _ d o w n w a r d
!
!> @file
!!
!! Downward shifting of the CTS.
!!
!! @section Copyright
!!
!! Copyright 2009-2013 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/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------------------------------------------------------------------------------
!> Downward shifting of the CTS.
!<------------------------------------------------------------------------------
subroutine shift_cts_downward(at1, at2_1, at2_2, at3_1, at3_2, &
              at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
              aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
              ai1, ai2, ai3, ai4, mean_accum_inv, dzeta_t, &
              dtime_temp, dtt_2dxi, dtt_2deta, dtime_temp_inv, &
              i, j)

use sico_types
use sico_variables

implicit none
integer(i4b) :: i, j
real(dp) :: zm_shift
real(dp) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
            at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
            at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
            ai1(0:kcmax), ai2(0:kcmax), ai3, ai4, &
            atr1, am1, am2, alb1
real(dp) :: aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld
real(dp) :: mean_accum_inv
real(dp) :: dtt_2dxi, dtt_2deta, dtime_temp, dtime_temp_inv, dzeta_t
real(dp) :: difftemp_a, difftemp_b, interpol

zm_shift = 1.0_dp   ! CTS shift in intervals of 1 m

!-------- Temperature discrepancy from the computation of the main
!         program --------

difftemp_a = temp_c_neu(0,j,i)-(-beta*h_c_neu(j,i))
if (difftemp_a.ge.0.0_dp) return

!-------- Shift CTS downward until it is too low --------

   10 continue

   zm_neu(j,i)  = zm_neu(j,i) - zm_shift

!  ------ Special case: CTS too close to the base

   if (zm_neu(j,i).le.zb(j,i)) then

      zm_shift = (zm_neu(j,i)+zm_shift)-(zb(j,i)+0.001_dp)
      zm_neu(j,i)  = zb(j,i)+0.001_dp
      h_c_neu(j,i) = h_c_neu(j,i)+h_t_neu(j,i)-0.001_dp
      h_t_neu(j,i) = 0.001_dp
!                   ! CTS positioned 1 mm above ice base --------

      dh_t_dtau(j,i) = (zm_neu(j,i)-zm(j,i))*dtime_temp_inv
      dzm_dtau(j,i)  = dzb_dtau(j,i)+dh_t_dtau(j,i)
      dh_c_dtau(j,i) = dzs_dtau(j,i)-dzm_dtau(j,i)

      am_perp(j,i) = am_perp_st(j,i) + dzm_dtau(j,i)

      call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
                 at4_1, at4_2, at5, at6, at7, atr1, &
                 am1, am2, alb1, &
                 aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
                 ai1, ai2, ai3, mean_accum_inv, dzeta_t, &
                 dtime_temp, dtt_2dxi, dtt_2deta, i, j)

      difftemp_b = difftemp_a
      difftemp_a = temp_c_neu(0,j,i)-(-beta*h_c_neu(j,i))

      if (difftemp_a.ge.0.0_dp) then ! CTS remains above the base

!    ---- Interpolate the CTS position from the last (_a) and the
!         last but one (_b) value, weighed with the temperature
!         discrepancies at the CTS --------

         interpol = difftemp_a/(difftemp_a-difftemp_b)*zm_shift

         zm_neu(j,i)  = zm_neu(j,i)  + interpol
         h_c_neu(j,i) = h_c_neu(j,i) - interpol
         h_t_neu(j,i) = h_t_neu(j,i) + interpol

         dh_t_dtau(j,i) = (zm_neu(j,i)-zm(j,i))*dtime_temp_inv
         dzm_dtau(j,i)  = dzb_dtau(j,i)+dh_t_dtau(j,i)
         dh_c_dtau(j,i) = dzs_dtau(j,i)-dzm_dtau(j,i)

         am_perp(j,i) = am_perp_st(j,i) + dzm_dtau(j,i)

         call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
                 at4_1, at4_2, at5, at6, at7, atr1, &
                 am1, am2, alb1, &
                 aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
                 ai1, ai2, ai3, mean_accum_inv, dzeta_t, &
                 dtime_temp, dtt_2dxi, dtt_2deta, i, j)

      else   ! CTS disappears

         n_cts_neu(j,i) = 0
         zm_neu(j,i) = zb(j,i)
         h_c_neu(j,i) = h_c_neu(j,i)+h_t_neu(j,i)
         h_t_neu(j,i) = 0.0_dp

         dh_t_dtau(j,i) = (zm_neu(j,i)-zm(j,i))*dtime_temp_inv
         dzm_dtau(j,i)  = dzb_dtau(j,i)+dh_t_dtau(j,i)
         dh_c_dtau(j,i) = dzs_dtau(j,i)-dzm_dtau(j,i)

         am_perp(j,i) = am_perp_st(j,i) + dzm_dtau(j,i)

         call calc_temp2(at1, at2_1, at2_2, at3_1, at3_2, &
                 at4_1, at4_2, at5, at6, at7, atr1, alb1, &
                 ai1, ai2, ai4, mean_accum_inv, &
                 dtime_temp, dtt_2dxi, dtt_2deta, i, j)

      end if

      return
   end if

!  ------ End of treatment of special case

   h_c_neu(j,i) = h_c_neu(j,i) + zm_shift
   h_t_neu(j,i) = h_t_neu(j,i) - zm_shift

   dh_t_dtau(j,i) = (zm_neu(j,i)-zm(j,i))*dtime_temp_inv
   dzm_dtau(j,i)  = dzb_dtau(j,i)+dh_t_dtau(j,i)
   dh_c_dtau(j,i) = dzs_dtau(j,i)-dzm_dtau(j,i)

   am_perp(j,i) = am_perp_st(j,i) + dzm_dtau(j,i)

   call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
                 at4_1, at4_2, at5, at6, at7, atr1, &
                 am1, am2, alb1, &
                 aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
                 ai1, ai2, ai3, mean_accum_inv, dzeta_t, &
                 dtime_temp, dtt_2dxi, dtt_2deta, i, j)

   difftemp_b = difftemp_a
   difftemp_a = temp_c_neu(0,j,i)-(-beta*h_c_neu(j,i))

if (difftemp_a.lt.0.0_dp) go to 10

!-------- Interpolate the CTS position from the last (_a) and the
!         last but one (_b) value, weighed with the temperature
!         discrepancies at the CTS --------

interpol = difftemp_a/(difftemp_a-difftemp_b)*zm_shift

zm_neu(j,i)  = zm_neu(j,i)  + interpol
h_c_neu(j,i) = h_c_neu(j,i) - interpol
h_t_neu(j,i) = h_t_neu(j,i) + interpol

dh_t_dtau(j,i) = (zm_neu(j,i)-zm(j,i))*dtime_temp_inv
dzm_dtau(j,i)  = dzb_dtau(j,i)+dh_t_dtau(j,i)
dh_c_dtau(j,i) = dzs_dtau(j,i)-dzm_dtau(j,i)

am_perp(j,i) = am_perp_st(j,i) + dzm_dtau(j,i)

call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
                at4_1, at4_2, at5, at6, at7, atr1, &
                am1, am2, alb1, &
                aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
                ai1, ai2, ai3, mean_accum_inv, dzeta_t, &
                dtime_temp, dtt_2dxi, dtt_2deta, i, j)

end subroutine shift_cts_downward
!