calc_bas_melt_m.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Module :  c a l c _ b a s _ m e l t _ m
!
!> @file
!!
!! Computation of the basal melting rate.
!!
!! @section Copyright
!!
!! Copyright 2009-2017 Ralf Greve, Ben Galton-Fenzi, Tatsuru Sato
!!
!! @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 basal melting rate.
!<------------------------------------------------------------------------------

module calc_bas_melt_m

  use sico_types_m
  use sico_variables_m
  use sico_vars_m

  implicit none

  private
  public :: calc_qbm

contains

!-------------------------------------------------------------------------------
!> Computation of the basal melting rate Q_bm.
!! Summation of Q_bm and Q_tld (water drainage rate from the temperate layer).
!<------------------------------------------------------------------------------
subroutine calc_qbm(time, z_sl, dzeta_c, dzeta_r)

use ice_material_properties_m, only : kappa_val

implicit none

real(dp), intent(in) :: time, z_sl
real(dp), intent(in) :: dzeta_c, dzeta_r

integer(i4b) :: i, j
integer(i4b) :: n_ocean, n_float
real(dp) :: aqbm1, aqbm2, aqbm3a, aqbm3b, aqbm4
real(dp) :: frictional_heating
real(dp) :: year_sec_inv, rhow_rho_ratio
real(dp) :: Q_bm_grounded, Q_bm_marine, Q_bm_floating
real(dp) :: qbm_min, qbm_max

!-------- Term abbreviations --------

year_sec_inv   = 1.0_dp/year_sec
rhow_rho_ratio = rho_w/rho

if (flag_aa_nonzero) then
   aqbm1 = (ea-1.0_dp)/(rho*l*aa*dzeta_c)
else
   aqbm1 = 1.0_dp/(rho*l*dzeta_c)
end if

aqbm2  = kappa_r/(rho*l*h_r*dzeta_r)
aqbm3a = g/l
aqbm3b = 1.0_dp/(rho*l)
aqbm4  = beta/(rho*l)

!-------- Computation of Q_bm --------

q_bm = 0.0_dp   ! initialisation

do i=1, imax-1
do j=1, jmax-1

   if (maske(j,i)==0) then   ! grounded ice

      if (n_cts(j,i)==-1) then

         frictional_heating = 0.0_dp
         q_bm(j,i)          = 0.0_dp

      else if (n_cts(j,i)==0) then

#if (DYNAMICS==2)
         if (.not.flag_shelfy_stream(j,i)) then
#endif
            frictional_heating &
               = -aqbm3a*h_c(j,i) &
                        *0.5_dp*(vx_t(0,j,i)+vx_t(0,j,i-1)) &
                        *dzs_dxi_g(j,i) &
                 -aqbm3a*h_c(j,i) &
                        *0.5_dp*(vy_t(0,j,i)+vy_t(0,j-1,i)) &
                        *dzs_deta_g(j,i)
#if (DYNAMICS==2)
         else   ! flag_shelfy_stream(j,i) == .true.

            frictional_heating &
               = aqbm3b &
                    * c_drag(j,i) &
                    * sqrt(vx_b_g(j,i)**2  &
                          +vy_b_g(j,i)**2) &
                                    **(1.0_dp+p_weert_inv(j,i))
         end if
#endif
         q_bm(j,i) =   aqbm1*(temp_c(1,j,i)-temp_c(0,j,i))/h_c(j,i) &
                            *kappa_val(temp_c(0,j,i)) &
                     - aqbm2*(temp_r(krmax,j,i)-temp_r(krmax-1,j,i)) &
                     + frictional_heating

      else   ! n_cts(j,i)==1

#if (DYNAMICS==2)
         if (.not.flag_shelfy_stream(j,i)) then
#endif
            frictional_heating &
               = -aqbm3a*(h_c(j,i)+h_t(j,i)) &
                        *0.5_dp*(vx_t(0,j,i)+vx_t(0,j,i-1)) &
                        *dzs_dxi_g(j,i) &
                 -aqbm3a*(h_c(j,i)+h_t(j,i)) &
                        *0.5_dp*(vy_t(0,j,i)+vy_t(0,j-1,i)) &
                        *dzs_deta_g(j,i)
#if (DYNAMICS==2)
         else   ! flag_shelfy_stream(j,i) == .true.

            frictional_heating &
               = aqbm3b &
                    * c_drag(j,i) &
                    * sqrt(vx_b_g(j,i)**2  &
                          +vy_b_g(j,i)**2) &
                                    **(1.0_dp+p_weert_inv(j,i))
         end if
#endif
         q_bm(j,i) =   aqbm4*kappa_val(temp_t_m(0,j,i)) &
                     - aqbm2*(temp_r(krmax,j,i)-temp_r(krmax-1,j,i)) &
                     + frictional_heating
           
      end if

#if (MARGIN==1 || MARGIN==2)

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

      continue   ! do nothing

#elif (MARINE_ICE_BASAL_MELTING==2)

      if ( (zb(j,i) < z_sl) &          ! marine ice
           .and. &
           (     (maske(j,i+1)==2) &   ! at least one
             .or.(maske(j,i-1)==2) &   ! nearest neighbour
             .or.(maske(j+1,i)==2) &   ! is
             .or.(maske(j-1,i)==2) &   ! ocean
           ) &
         ) then

         q_bm(j,i) = qbm_marine *year_sec_inv*rhow_rho_ratio
                                ! m/a water equiv. -> m/s ice equiv.

      end if

#elif (MARINE_ICE_BASAL_MELTING==3)

      if ( (zb(j,i) < z_sl) &          ! marine ice
           .and. &
           (     (maske(j,i+1)==2) &   ! at least one
             .or.(maske(j,i-1)==2) &   ! nearest neighbour
             .or.(maske(j+1,i)==2) &   ! is
             .or.(maske(j-1,i)==2) &   ! ocean
           ) &
         ) then

         n_ocean = 0
         if (maske(j,i+1)==2) n_ocean = n_ocean+1
         if (maske(j,i-1)==2) n_ocean = n_ocean+1
         if (maske(j+1,i)==2) n_ocean = n_ocean+1
         if (maske(j-1,i)==2) n_ocean = n_ocean+1

         if ( n_ocean > 0 ) then

            q_bm_grounded = q_bm(j,i)
            q_bm_marine   = qbm_marine *year_sec_inv*rhow_rho_ratio
                                       ! m/a water equiv. -> m/s ice equiv.

            q_bm(j,i) = (1.0_dp-0.25_dp*real(n_ocean,dp)) * q_bm_grounded &
                               +0.25_dp*real(n_ocean,dp)  * Q_bm_marine
                      ! weighed average of grounded ice melting (computed)
                      ! and marine ice melting (prescribed)
         else
            write(6,'(a)') ' >>> calc_qbm: Marine ice margin point does not'
            write(6,'(a)') '               have an ocean neighbour!'
            stop
         end if

      end if

#else
      stop ' >>> calc_qbm: MARINE_ICE_BASAL_MELTING must be 1, 2 or 3!'
#endif

#elif (MARGIN==3)

      if (flag_grounding_line_1(j,i)) then
                                ! grounding line (grounded-ice side)

#if (FLOATING_ICE_BASAL_MELTING==1 || FLOATING_ICE_BASAL_MELTING==4 \
                                   || floating_ice_basal_melting==5)

         continue   ! do nothing

#elif (FLOATING_ICE_BASAL_MELTING==2)

         q_bm(j,i) = qbm_float_1 *year_sec_inv*rhow_rho_ratio
                                ! m/a water equiv. -> m/s ice equiv.

#elif (FLOATING_ICE_BASAL_MELTING==3)

         n_float = 0
         if (maske(j,i+1)>=2) n_float = n_float+1
         if (maske(j,i-1)>=2) n_float = n_float+1
         if (maske(j+1,i)>=2) n_float = n_float+1
         if (maske(j-1,i)>=2) n_float = n_float+1

         if ( n_float > 0 ) then

            q_bm_grounded = q_bm(j,i)
            q_bm_floating = qbm_float_1 *year_sec_inv*rhow_rho_ratio
                                        ! m/a water equiv. -> m/s ice equiv.

            q_bm(j,i) = (1.0_dp-0.25_dp*real(n_float,dp)) * q_bm_grounded &
                               +0.25_dp*real(n_float,dp)  * Q_bm_floating
                      ! weighed average of melting for grounded and floating ice
         else
            write(6,'(a)') ' >>> calc_qbm: Grounding line point does not have'
            write(6,'(a)') '               a floating ice or ocean neighbour!'
            stop
         end if

#else
         write(6, fmt='(a)') ' >>> calc_qbm: FLOATING_ICE_BASAL_MELTING'
         write(6, fmt='(a)') '               must be 1, 2, 3, 4 or 5!'
         stop
#endif

      else if ( (zb(j,i) < z_sl) &          ! marine ice margin
                .and. &
                (     (maske(j,i+1)>=2) &   !  (at least one
                  .or.(maske(j,i-1)>=2) &   !   nearest neighbour
                  .or.(maske(j+1,i)>=2) &   !   is
                  .or.(maske(j-1,i)>=2) &   !   ocean)
                ) &
              ) then

#if (MARINE_ICE_BASAL_MELTING==1)

         continue   ! do nothing

#elif (MARINE_ICE_BASAL_MELTING==2)

         q_bm(j,i) = qbm_marine *year_sec_inv*rhow_rho_ratio
                                ! m/a water equiv. -> m/s ice equiv.

#elif (MARINE_ICE_BASAL_MELTING==3)

         n_ocean = 0
         if (maske(j,i+1)>=2) n_ocean = n_ocean+1
         if (maske(j,i-1)>=2) n_ocean = n_ocean+1
         if (maske(j+1,i)>=2) n_ocean = n_ocean+1
         if (maske(j-1,i)>=2) n_ocean = n_ocean+1

         if ( n_ocean > 0 ) then

            q_bm_grounded = q_bm(j,i)
            q_bm_marine   = qbm_marine *year_sec_inv*rhow_rho_ratio
                                       ! m/a water equiv. -> m/s ice equiv.

            q_bm(j,i) = (1.0_dp-0.25_dp*real(n_ocean,dp)) * q_bm_grounded &
                               +0.25_dp*real(n_ocean,dp)  * Q_bm_marine
                      ! weighed average of grounded ice melting (computed)
                      ! and marine ice melting (prescribed)
         else
            write(6,'(a)') ' >>> calc_qbm: Marine ice margin point does not'
            write(6,'(a)') '          have a floating ice or ocean neighbour!'
            stop
         end if

#else
         stop ' >>> calc_qbm: MARINE_ICE_BASAL_MELTING must be 1, 2 or 3!'
#endif

      end if

   else if ( (maske(j,i)==2).or.(maske(j,i)==3) ) then   ! floating ice or ocean

#if (FLOATING_ICE_BASAL_MELTING<=3)

      if (flag_grounding_line_2(j,i)) then
                                ! grounding line (floating-ice side)

         q_bm(j,i) = qbm_float_1 *year_sec_inv*rhow_rho_ratio
                                 ! m/a water equiv. -> m/s ice equiv.

      else if ( zl(j,i) > z_abyss ) then   ! floating ice over continental shelf

         q_bm(j,i) = qbm_float_2 *year_sec_inv*rhow_rho_ratio
                                 ! m/a water equiv. -> m/s ice equiv.

      else   ! ( zl(j,i) <= Z_ABYSS ), floating ice over abyssal ocean

         q_bm(j,i) = qbm_float_3 *year_sec_inv*rhow_rho_ratio
                                 ! m/a water equiv. -> m/s ice equiv.

      end if

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

      if ( zl(j,i) > z_abyss ) then   ! floating ice over continental shelf
         call sub_ice_shelf_melting_param(time, z_sl, &
                                          year_sec_inv, rhow_rho_ratio, &
                                          i, j, q_bm_floating)
         q_bm(j,i) = q_bm_floating
      else   ! ( zl(j,i) <= Z_ABYSS ), floating ice over abyssal ocean
         q_bm(j,i) = qbm_float_3 *year_sec_inv*rhow_rho_ratio
                                 ! m/a water equiv. -> m/s ice equiv.
      end if

#else
      write(6, fmt='(a)') ' >>> calc_qbm: FLOATING_ICE_BASAL_MELTING'
      write(6, fmt='(a)') '               must be 1, 2, 3, 4 or 5!'
      stop
#endif

#else
      stop ' >>> calc_qbm: MARGIN must be 1, 2 or 3!'
#endif

   end if

end do
end do

!-------- Sum of Q_bm and Q_tld --------

q_b_tot = q_bm + q_tld

!-------- Limitation of Q_bm, Q_tld and Q_b_tot
!                       (only for grounded ice) --------

#if (defined(QBM_MIN))
   qbm_min = qbm_min *year_sec_inv*rhow_rho_ratio
                      ! m/a water equiv. -> m/s ice equiv.
#elif (defined(QB_MIN)) /* obsolete */
   qbm_min = qb_min   ! m/s ice equiv.
#else
   stop ' >>> calc_qbm: Limiter QBM_MIN is undefined!'
#endif

#if (defined(QBM_MAX))
   qbm_max = qbm_max *year_sec_inv*rhow_rho_ratio
                      ! m/a water equiv. -> m/s ice equiv.
#elif (defined(QB_MAX)) /* obsolete */
   qbm_max = qb_max   ! m/s ice equiv.
#else
   stop ' >>> calc_qbm: Limiter QBM_MAX is undefined!'
#endif

#if ( defined(MARINE_ICE_BASAL_MELTING) \
      && ( marine_ice_basal_melting==2 || marine_ice_basal_melting==3 ) )

if (qbm_marine*year_sec_inv*rhow_rho_ratio > qbm_max) then
   write (6,'(a)') ' >>> calc_qbm: QBM_MARINE'
   write (6,'(a)') '               (basal melting rate at the ice front)'
   write (6,'(a)') '               is larger than the limiter qbm_max!'
   stop
end if

#endif

#if ( MARGIN==3 \
      && defined(floating_ice_basal_melting) \
      && floating_ice_basal_melting<=3 )

if (qbm_float_1*year_sec_inv*rhow_rho_ratio > qbm_max) then
   write (6,'(a)') ' >>> calc_qbm: QBM_FLOAT_1'
   write (6,'(a)') '          (basal melting rate in the grounding line zone)'
   write (6,'(a)') '          is larger than the limiter qbm_max!'
   stop
end if

#endif

do i=0, imax
do j=0, jmax
   if (maske(j,i)==0) then   ! grounded ice
      if (q_bm(j,i)    < qbm_min) q_bm(j,i)    = 0.0_dp
      if (q_bm(j,i)    > qbm_max) q_bm(j,i)    = qbm_max
      if (q_tld(j,i)   < qbm_min) q_tld(j,i)   = 0.0_dp
      if (q_tld(j,i)   > qbm_max) q_tld(j,i)   = qbm_max
      if (q_b_tot(j,i) < qbm_min) q_b_tot(j,i) = 0.0_dp
      if (q_b_tot(j,i) > qbm_max) q_b_tot(j,i) = qbm_max
   end if
end do
end do

end subroutine calc_qbm

!-------------------------------------------------------------------------------
!> Sub-ice-shelf melting parameterisation.
!<------------------------------------------------------------------------------
subroutine sub_ice_shelf_melting_param(time, z_sl, &
                                       year_sec_inv, rhow_rho_ratio, &
                                       i, j, q_bm_floating)

implicit none

integer(i4b), intent(in)  :: i, j
real(dp),     intent(in)  :: time, z_sl
real(dp),     intent(in)  :: year_sec_inv, rhow_rho_ratio

real(dp),     intent(out) :: Q_bm_floating

real(dp) :: time_in_years
real(dp) :: Toc, Tmb, T_forcing, Omega, alpha, draft0, draft
real(dp) :: lon_d, lat_d, Phi_par
real(dp) :: H_w_now, Q_bm_scaling_factor

real(dp), parameter :: beta_sw = 7.61e-04_dp

#if (FLOATING_ICE_BASAL_MELTING==5)

real(dp), parameter :: c_sw = 3974.0_dp, &
                       g_t  =    5.0e-05_dp

#endif

time_in_years = time*year_sec_inv

#if (FLOATING_ICE_BASAL_MELTING==4)

toc   = temp_ocean
omega = omega_qbm *year_sec_inv*rhow_rho_ratio
                  ! m/[a*degC^alpha] water equiv. -> m/[s*degC^alpha] ice equiv.
alpha = alpha_qbm

#elif (FLOATING_ICE_BASAL_MELTING==5)

#if (defined(ANT))   /* Antarctic ice sheet */

!-------- Definition of the sectors --------

lon_d          = lambda(j,i) *pi_180_inv   ! rad -> deg
lat_d          = phi(j,i)    *pi_180_inv   ! rad -> deg

lon_d          = modulo(lon_d+180.0_dp, 360.0_dp)-180.0_dp
                                           ! mapping to interval
                                           ! [-180 deg, +180 deg)

if ((lon_d>=-10.0_dp).and.(lon_d<60.0_dp)) then
                                      ! Western East Antarctica
   n_sector(j,i) = 1_i2b

   toc    =   0.23_dp
   omega  =   0.0076_dp
   alpha  =   1.33_dp
   draft0 = 200.0_dp

else if ((lon_d>=60.0_dp).and.(lon_d<80.0_dp)) then
                                      ! Amery/Prydz Bay
   n_sector(j,i) = 2_i2b

   toc    =  -1.61_dp
   omega  =   0.0128_dp
   alpha  =   0.94_dp
   draft0 = 200.0_dp

else if ((lon_d>=80.0_dp).and.(lon_d<130.0_dp)) then
                                      ! Sabrina Coast/
                                      ! Aurora subglacial basin
   n_sector(j,i) = 3_i2b

   toc    =  -0.16_dp
   omega  =   0.0497_dp
   alpha  =   0.85_dp
   draft0 = 200.0_dp

else if ( ((lon_d>=130.0_dp).and.(lon_d<159.0_dp)) &
          .or. &
          ((lon_d>=159.0_dp).and.(lon_d<170.0_dp) &
                            .and.(lat_d>=-72.0_dp)) ) then
                                      ! George V Coast/
                                      ! Wilkes subglacial basin
   n_sector(j,i) = 4_i2b

   toc    =  -0.22_dp
   omega  =   0.0423_dp
   alpha  =   0.11_dp
   draft0 = 200.0_dp

else if ( (lon_d>=159.0_dp) &
          .or. &
          (lon_d<-140.0_dp) &
          .or. &
          ((lon_d>=-140.0_dp).and.(lon_d<-120.0_dp) &
                             .and.(lat_d<-77.0_dp)) ) then
                                      ! Ross Sea
   n_sector(j,i) = 5_i2b

   toc    =  -1.8_dp
   omega  =   0.0181_dp
   alpha  =   0.73_dp
   draft0 = 200.0_dp

else if ( ((lon_d>=-140.0_dp).and.(lon_d<-120.0_dp) &
                             .and.(lat_d>=-77.0_dp)) &
          .or. &
          ((lon_d>=-120.0_dp).and.(lon_d<-90.0_dp)) ) then
                                      ! Amundsen Sea
   n_sector(j,i) = 6_i2b

   toc    =   0.42_dp
   omega  =   0.1023_dp
   alpha  =   0.26_dp
   draft0 = 200.0_dp

else if ( ((lon_d>=-90.0_dp).and.(lon_d<-66.0_dp) &
                            .and.(lat_d>=-74.0_dp)) ) then
                                      ! Bellingshausen Sea
   n_sector(j,i) = 7_i2b

   toc    =   0.62_dp
   omega  =   0.0644_dp
   alpha  =   0.16_dp
   draft0 = 200.0_dp

else
                                      ! Weddell Sea
   n_sector(j,i) = 8_i2b

   toc    =  -1.55_dp
   omega  =   0.0197_dp
   alpha  =   0.26_dp
   draft0 = 200.0_dp

endif

#else   /* not Antarctic ice sheet */

q_bm_floating = 0.0_dp   ! dummy return value

write(6, fmt='(a)') ' >>> sub_ice_shelf_melting_param:'
write(6, fmt='(a)') '         Case FLOATING_ICE_BASAL_MELTING==5'
write(6, fmt='(a)') '         only defined for Antarctica!'
stop

#endif

#else

continue   ! do nothing

#endif

!-------- Computation of the sub-ice-shelf melting rate --------

if (maske(j,i)==2_i2b) then   ! ocean
   draft   = 0.0_dp
   h_w_now = max((z_sl-zl(j,i)), 0.0_dp)
else if (maske(j,i)==3_i2b) then   ! floating ice
   draft   = max((z_sl   -zb(j,i)), 0.0_dp)
   h_w_now = max((zb(j,i)-zl(j,i)), 0.0_dp)
else
   write(6, fmt='(a)') ' >>> sub_ice_shelf_melting_param:'
   write(6, fmt='(a)') '         Routine must not be called for maske(j,i) < 2!'
   stop
end if

tmb = -beta_sw*draft - delta_tm_sw   ! temperature at the ice shelf base

t_forcing = max((toc-tmb), 0.0_dp)   ! thermal forcing

#if (defined(H_W_0))

if (h_w_0 > eps) then
   q_bm_scaling_factor = tanh(h_w_now/h_w_0)
else
   q_bm_scaling_factor = 1.0_dp
end if

#else
q_bm_scaling_factor = 1.0_dp
#endif

#if (FLOATING_ICE_BASAL_MELTING==4)

q_bm_floating  = q_bm_scaling_factor*omega*t_forcing**alpha

#elif (FLOATING_ICE_BASAL_MELTING==5)

phi_par = rho_sw*c_sw*g_t/(rho*l)

q_bm_floating  = q_bm_scaling_factor &
                 *phi_par*omega*t_forcing*(draft/draft0)**alpha

#else

q_bm_floating = 0.0_dp   ! dummy return value

write(6, fmt='(a)') ' >>> sub_ice_shelf_melting_param:'
write(6, fmt='(a)') '         FLOATING_ICE_BASAL_MELTING must be 4 or 5!'
stop

#endif

#if (defined(INITMIP_BMB_ANOM_FILE))

!  ------ Correction for ISMIP InitMIP

if ((time_in_years > 0.0_dp).and.(time_in_years <= 40.0_dp)) then
   q_bm_floating = q_bm_floating &
                   + 0.025_dp*floor(time_in_years) * ab_anom_initmip(j,i)
else if (time_in_years > 40.0_dp) then
   q_bm_floating = q_bm_floating + ab_anom_initmip(j,i)
end if

#endif

end subroutine sub_ice_shelf_melting_param

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

end module calc_bas_melt_m
!