sub_ice_shelf_melting_param.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  s u b _ i c e _ s h e l f _ m e l t i n g _ p a r a m
!
!> @file
!!
!! Sub-ice-shelf melting parameterisation for Antarctica.
!!
!! @section Copyright
!!
!! Copyright 2013-2016 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/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------------------------------------------------------------------------------
!> Sub-ice-shelf melting parameterisation for Antarctica.
!<------------------------------------------------------------------------------
subroutine sub_ice_shelf_melting_param(z_sl, i, j, Q_bm_floating)

use sico_types
use sico_variables
use sico_vars

implicit none

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

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

real(dp) :: lon_d, lat_d, tmb, phi_par
real(dp) :: toc, omega, alpha, draft0, draft
real(dp) :: h_w_now, q_bm_scaling_factor

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

real(dp), parameter :: euler = 2.718281828459045_dp

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

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 (defined(ANT))   /* Antarctic ice sheet */

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

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

!-------- 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 should not be called for maske(j,i) < 2!'
   stop
end if

tmb = -beta_sw*draft - delta_tm_sw

#if (defined(H_W_0))

if (h_w_0 > eps) then
   q_bm_scaling_factor = tanh(euler*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

q_bm_floating  = q_bm_scaling_factor &
                 *phi_par*omega*(toc-tmb)*(draft/draft0)**alpha

#else   /* not Antarctic ice sheet */

write(6, fmt='(a)') ' sub_ice_shelf_melting_param:'
write(6, fmt='(a)') '     Parameterisation only defined for Antarctica!'
stop

#endif

end subroutine sub_ice_shelf_melting_param
!