calc_enhance.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  c a l c _ e n h a n c e
!
!> @file
!!
!! Computation of the flow enhancement factor.
!!
!! @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 flow enhancement factor.
!<------------------------------------------------------------------------------
subroutine calc_enhance(time)

use sico_types
use sico_variables
use sico_vars

implicit none

real(dp), intent(in) :: time

integer(i4b) :: i, j, kc, kt
real(dp) :: age_trans, date_trans1, date_trans2, date_trans3
real(dp) :: enh_shear_compr_diff, enh_mult

#if (ENHMOD==1)

enh_t = enh_fact
enh_c = enh_fact

#elif (ENHMOD==2)

age_trans = age_trans_0*year_sec

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

   do kt=0, ktmax
      if (age_t(kt,j,i) < age_trans) then
         enh_t(kt,j,i) = enh_intg   ! Holocene ice
      else
         enh_t(kt,j,i) = enh_fact   ! Pleistocene ice
      end if
   end do

   do kc=0, kcmax
      if (age_c(kc,j,i) < age_trans) then
         enh_c(kc,j,i) = enh_intg   ! Holocene ice
      else
         enh_c(kc,j,i) = enh_fact   ! Pleistocene ice
      end if
   end do

end do
end do

#elif (ENHMOD==3)

date_trans1 = date_trans1_0*year_sec
date_trans2 = date_trans2_0*year_sec
date_trans3 = date_trans3_0*year_sec

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

   do kt=0, ktmax
      if ( (time-age_t(kt,j,i)) < date_trans1 ) then
         enh_t(kt,j,i) = enh_fact   ! pre-Eemian ice
      else if ( ((time-age_t(kt,j,i)) >= date_trans1).and. &
                ((time-age_t(kt,j,i)) <  date_trans2) ) then
         enh_t(kt,j,i) = enh_intg   ! Eemian ice
      else if ( ((time-age_t(kt,j,i)) >= date_trans2).and. &
                ((time-age_t(kt,j,i)) <  date_trans3) ) then
         enh_t(kt,j,i) = enh_fact   ! Weichselian ice
      else
         enh_t(kt,j,i) = enh_intg   ! Holocene ice
      end if
   end do

   do kc=0, kcmax
      if ( (time-age_c(kc,j,i)) < date_trans1 ) then
         enh_c(kc,j,i) = enh_fact   ! pre-Eemian ice
      else if ( ((time-age_c(kc,j,i)) >= date_trans1).and. &
                ((time-age_c(kc,j,i)) <  date_trans2) ) then
         enh_c(kc,j,i) = enh_intg   ! Eemian ice
      else if ( ((time-age_c(kc,j,i)) >= date_trans2).and. &
                ((time-age_c(kc,j,i)) <  date_trans3) ) then
         enh_c(kc,j,i) = enh_fact   ! Weichselian ice
      else
         enh_c(kc,j,i) = enh_intg   ! Holocene ice
      end if
   end do

end do
end do

#elif (ENHMOD==4 || ENHMOD==5)

enh_shear_compr_diff = enh_shear-enh_compr

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

   do kt=0, ktmax
      enh_t(kt,j,i) = enh_compr &
                      + enh_shear_compr_diff*lambda_shear_t(kt,j,i)**2
   end do

   do kc=0, kcmax
      enh_c(kc,j,i) = enh_compr &
                      + enh_shear_compr_diff*lambda_shear_c(kc,j,i)**2
   end do

end do
end do

#else
   stop ' calc_enhance: Parameter ENHMOD must be either 1, 2, 3, 4 or 5!'
#endif

!-------- Modification for floating ice --------

#if (MARGIN==3)

#if (ENHMOD==1 || ENHMOD==2 || ENHMOD==3 || ENHMOD==4)

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

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

      do kt=0, ktmax
         enh_t(kt,j,i) = enh_shelf
      end do

      do kc=0, kcmax
         enh_c(kc,j,i) = enh_shelf
      end do

   end if

end do
end do

#elif (ENHMOD==5)

continue   ! no modification

#else
   stop ' calc_enhance: Parameter ENHMOD must be either 1, 2, 3, 4 or 5!'
#endif

#endif

!-------- Modification due to dust content (Martian polar caps only) --------

#if ( defined(NMARS) \
      || defined(smars) )   /* martian ice sheet */

#if (FLOW_LAW==1)
   enh_mult = exp(-3.0_dp*2.0_dp*frac_dust)
#elif (FLOW_LAW==2)
   enh_mult = exp(-1.8_dp*2.0_dp*frac_dust)
#elif (FLOW_LAW==3)
   enh_mult = exp(-4.0_dp*2.0_dp*frac_dust)
#elif (FLOW_LAW==4)
   stop ' calc_enhance: FLOW_LAW==4 has not been implemented yet!'
#endif

enh_t = enh_t * enh_mult
enh_c = enh_c * enh_mult

#endif

end subroutine calc_enhance
!