init_temp_water_age.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  i n i t _ t e m p _ w a t e r _ a g e
!
!> @file
!!
!! Initial temperature, water content and age.
!!
!! @section Copyright
!!
!! Copyright 2009-2016 Ralf Greve, Thorben Dunse
!!
!! @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/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------------------------------------------------------------------------------
!> Initial temperature, water content and age.
!<------------------------------------------------------------------------------
subroutine init_temp_water_age()

use sico_types
use sico_variables
use sico_vars
use ice_material_quantities, only : kappa_val, c_val

implicit none

integer(i4b) :: i, j, kc, kt, kr
real(dp)     :: kappa_const_val, c_const_val
real(dp)     :: as_val, h_val, qgeo_val, k, z_above_base, erf_val_1, erf_val_2
real(dp)     :: temp_ice_base, temp_scale_factor

!-------- Initial ice temperature --------

!  ------ Present topography

#if (ANF_DAT==1)

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

#if (!defined(TEMP_INIT) || TEMP_INIT==1)

#if (defined(NMARS) || defined(SMARS))   /* Polar caps of Mars */

   do kc=0, kcmax
      temp_c(kc,j,i) = -100.0_dp
   end do

#else   /* all other domains */

   do kc=0, kcmax
      temp_c(kc,j,i) =  -10.0_dp
   end do

#endif

#elif (TEMP_INIT==2)

   do kc=0, kcmax
      temp_c(kc,j,i) = temp_s(j,i)
   end do

#elif (TEMP_INIT==3)

#if (defined(NMARS) || defined(SMARS))   /* Polar caps of Mars */
   kappa_const_val = kappa_val(-100.0_dp)
#else   /* all other domains */
   kappa_const_val = kappa_val(-10.0_dp)
#endif

   if (maske(j,i)<=2_i2b) then

      do kc=0, kcmax

      !%%% TDUNSE 05.12.2008: vertical temp profile depends on surface temp 
      !%%% and increases lineray with depth according to the geothermal heat flux;
      !%%% upper temperature limit of ice: pressure melting point
      !%%% KAPPA_ice = 2.072_dp and q_geo = 50 mW/m^2 => 0.024 Km^-1
      !%%% current description:
      !
      !    temp_c(kc,j,i) = temp_ma_present(j,i) &
      !                   + q_geo(j,i)/2.072_dp &
      !                     *H_c(j,i)*(1.0_dp-eaz_c_quotient(kc))
      ! 
      !    if (temp_c(kc,j,i) > -BETA*H_c(j,i)) then
      !        temp_c(kc,j,i) = -BETA*H_c(j,i)
      !    end if

         temp_c(kc,j,i) = temp_s(j,i) &
                          + (q_geo(j,i)/kappa_const_val) &
                            *h_c(j,i)*(1.0_dp-eaz_c_quotient(kc))
                          ! linear temperature distribution according to the
                          ! geothermal heat flux
      end do

      if (temp_c(0,j,i) >= -beta*h_c(j,i)) then

         temp_ice_base = -beta*h_c(j,i)

         do kc=0, kcmax
            temp_c(kc,j,i) = temp_s(j,i) &
                             + (temp_ice_base-temp_s(j,i)) &
                               *(1.0_dp-eaz_c_quotient(kc))
         end do

      end if

   else   ! maske(j,i)==3_i2b, floating ice

      temp_ice_base = -beta*h_c(j,i) - delta_tm_sw

      do kc=0, kcmax
         temp_c(kc,j,i) = temp_s(j,i) &
                          + (temp_ice_base-temp_s(j,i)) &
                            *(1.0_dp-eaz_c_quotient(kc))
      end do

   end if

#elif (TEMP_INIT==4)

#if (defined(NMARS) || defined(SMARS))   /* Polar caps of Mars */
   kappa_const_val = kappa_val(-100.0_dp)
       c_const_val =     c_val(-100.0_dp)
#else   /* all other domains */
   kappa_const_val = kappa_val(-10.0_dp)
       c_const_val =     c_val(-10.0_dp)
#endif

   if (maske(j,i)<=2_i2b) then
      as_val = max(as_perp(j,i), epsi)
   else   ! maske(j,i)==3_i2b, floating ice
      as_val = epsi   ! this will produce an almost linear temperature profile
   end if

   h_val    = max(h_c(j,i)  , eps)
   qgeo_val = max(q_geo(j,i), eps)

   k = sqrt( (kappa_const_val/(rho*c_const_val)) * (h_val/as_val) )

   erf_val_1 = erf(h_c(j,i)/(sqrt(2.0_dp)*k))

   do kc=0, kcmax
      z_above_base   = h_c(j,i)*eaz_c_quotient(kc)
      erf_val_2      = erf(z_above_base/(sqrt(2.0_dp)*k))
      temp_c(kc,j,i) = temp_s(j,i) &
                        + (qgeo_val/kappa_const_val) &
                          * sqrt(0.5_dp*pi)*k*(erf_val_1-erf_val_2)
   end do

   if ( (maske(j,i) <= 2_i2b).and.(temp_c(0,j,i) >= -beta*h_c(j,i)) ) then
      temp_ice_base     = -beta*h_c(j,i)
      temp_scale_factor = (temp_ice_base-temp_s(j,i)) &
                                /(temp_c(0,j,i)-temp_s(j,i))
   else if (maske(j,i) == 3_i2b) then
      temp_ice_base     = -beta*h_c(j,i)-delta_tm_sw
      temp_scale_factor = (temp_ice_base-temp_s(j,i)) &
                                /(temp_c(0,j,i)-temp_s(j,i))
   else
      temp_scale_factor = 1.0_dp
   end if

   do kc=0, kcmax
      temp_c(kc,j,i) = temp_s(j,i) &
                       + temp_scale_factor*(temp_c(kc,j,i)-temp_s(j,i))
   end do

#else

   write(6, fmt='(a)') ' ice_temp:'
   write(6, fmt='(a)') '     TEMP_INIT must be set to either 1, 2, 3 or 4!'
   stop

#endif

end do
end do

!  ------ Ice-free, relaxed bedrock

#elif (ANF_DAT==2)

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

   do kc=0, kcmax
      temp_c(kc,j,i) = temp_s(j,i)
   end do

end do
end do

#else

write(6, fmt='(a)') ' ice_temp:'
write(6, fmt='(a)') '     Only to be called for ANF_DAT==1 or ANF_DAT==2!'
stop

#endif

!-------- Initial lithosphere temperature --------

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

   do kr=0, krmax
      temp_r(kr,j,i) = temp_c(0,j,i) &
                       + (q_geo(j,i)/kappa_r) &
                         *h_r*(1.0_dp-zeta_r(kr))
            ! linear temperature distribution according to the
            ! geothermal heat flux
   end do

end do
end do

!-------- Initial water content --------

omega_c = 0.0_dp   ! only required for the enthalpy method
omega_t = 0.0_dp

!-------- Initial age --------

#if (defined(ASF))   /* Austfonna */

age_c =  3500.0_dp*year_sec
age_t =  3500.0_dp*year_sec

#else   /* all other domains */

age_c = 15000.0_dp*year_sec
age_t = 15000.0_dp*year_sec

#endif

end subroutine init_temp_water_age
!