discharge_workers_m.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Module :  d i s c h a r g e _ w o r k e r s _ m
!
!> @file
!!
!! Ice discharge parameterization for the Greenland ice sheet.
!!
!! References:
!! @li Calov, R. and A. Robinson,  and M. Perrette, \n
!!     and A. Ganopolski 2015.\n
!!     Simulating the Greenland ice sheet under present-day and 
!!     palaeo constraints including a new discharge parameterization.\n
!!     The Cryosphere 9: 179-196.
!!
!! @section Copyright
!!
!! Copyright 2011-2017 Reinhard Calov, Andrey Ganopolski
!! Potsdam Institute for Climate Impact Research  
!!
!! @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/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------------------------------------------------------------------------------
!> Ice discharge parameterization for the Greenland ice sheet.
!<------------------------------------------------------------------------------
module discharge_workers_m

  use sico_types_m
  use sico_variables_m
  use sico_vars_m

  use compare_float_m

  implicit none

  integer(i4b), private :: disc
  integer(i4b), private :: n_discharge_call, iter_mar_coa
  real(dp),     private :: c_dis_0, s_dis, c_dis_fac
  real(dp),     private :: t_sub_pd, alpha_sub, alpha_o, m_h, m_d, r_mar_eff
  real(dp),     private :: t_sea_freeze
  real(dp),     public  :: dt_glann, dt_sub

  integer(i2b), dimension(0:JMAX,0:IMAX), public  :: mask_mar
  real(dp),     dimension(0:JMAX,0:IMAX), private :: c_dis
  real(dp),     dimension(0:JMAX,0:IMAX), public  :: cst_dist, cos_grad_tc
  real(dp),     dimension(0:JMAX,0:IMAX), public  :: dis_perp

  private
  public :: disc_param, disc_fields, calc_c_dis_0, discharge

contains

!-------------------------------------------------------------------------------
!> Ice discharge parameters (Greenland).
!! [Assign ice discharge parameters.]
!<------------------------------------------------------------------------------
  subroutine disc_param(dtime)

  ! Author: Reinhard Calov
  ! Institution: Potsdam Institute for Climate Impact Research  
  ! Date: 10.6.16    

  ! Purpose: Calculate discharge parameters.

  implicit none

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

  real(dp) :: dtime_mar_coa

  disc = disc 

  c_dis_0   = c_dis_0
  c_dis_fac = c_dis_fac

  m_h = m_h
  m_d = m_d

  r_mar_eff = r_mar_eff *1000.0_dp   ! km -> m

#if (defined(S_DIS))
  s_dis = s_dis
#else
  s_dis = 1.0_dp   ! default value
#endif

#if (defined(ALPHA_SUB))
  alpha_sub = alpha_sub
#else
  alpha_sub = 0.5_dp   ! default value
#endif

#if (defined(ALPHA_O))
  alpha_o = alpha_o
#else
  alpha_o = 1.0_dp   ! default value
#endif

  t_sea_freeze = 0.0_dp - delta_tm_sw   ! freezing temperature of sea water

  n_discharge_call = -1

#if (defined(DTIME_MAR_COA0))
  dtime_mar_coa = dtime_mar_coa0*year_sec   ! a -> s
#else
  stop ' >>> disc_param: DTIME_MAR_COA0 not defined in header file!'
#endif

  if (.not.approx_integer_multiple(dtime_mar_coa, dtime, eps_sp_dp)) &
     stop ' >>> disc_param: dtime_mar_coa must be a multiple of dtime!'

  iter_mar_coa = nint(dtime_mar_coa/dtime)

#if (GRID > 1)
  stop ' >>> disc_param: GRID==2 not allowed for this application!'
#endif

#if (ANF_DAT != 1  && defined(EXEC_MAKE_C_DIS_0))
  stop ' >>> disc_param: EXEC_MAKE_C_DIS_0 defined requires ANF_DAT==1!'
#endif

  end subroutine disc_param

!-------------------------------------------------------------------------------
!> Constant in ice discharge parameterization (Greenland).
!! [Determine (amount of magnitude of) constant in ice discharge
!! parameterization.]
!<------------------------------------------------------------------------------
  subroutine calc_c_dis_0(dxi, deta)

  ! Author: Reinhard Calov
  ! Institution: Potsdam Institute for Climate Impact Research
  ! Date: 8.6.16    
  ! Purpose: Compute c_dis_0 indirectly from present-day topography 
  !          from the condition that the parameterized total discharge 
  !          equals the discharge give by dis_target. Note: This c_dis_0 could
  !          differ from the optimal one yielded from the siumualted ice thickness.
  !
  !          This is very useful, because c_dis_0 can vary a lot for different powers 
  !          m_D and m_H. This way we easier yield the approximately right value
  !          for c_dis_0 for given powers m_D and m_H.

  ! This routine is not to be used regularly, and it is only executed if the
  ! parameter EXEC_MAKE_C_DIS_0 is defined in the header file.

  implicit none

  real(dp), intent(in) :: dxi, deta

  integer(i4b) :: i, j
  real(dp) :: disc_target, disc_tot, dT_sub

  ! targeted total ice discharge
  ! 350 Gt/yr Calov et al. (2015)   !!! 400 Gt/yr van den Broeke et al (2016)
  disc_target = 350.0_dp ! in Gt/yr  

  h_c=zs-zb; h_t=0.0_dp

  c_dis_0   =  1.0_dp
  c_dis_fac =  1.0_dp
  s_dis     =  1.0_dp
  call disc_fields()

  ! ensure that we get present-day discharge here (disc=1). 
  dt_glann = 0.0_dp
  n_discharge_call = -1
  call discharge(dxi, deta)

  !  disc_tot=0.0_dp
  !  do i=0, IMAX
  !  do j=0, JMAX
  !    if(maske(j,i).eq.0.or.maske(j,i).eq.3) then
  !      disc_tot=disc_tot+dis_perp(j,i)*area(j,i)
  !    end if
  !  end do
  !  end do
  !
  !  disc_tot = disc_tot*YEAR_SEC         ! m^3/s -> m^3/a

  disc_tot = 0.0_dp

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

     if (mask_mar(j,i) == 1) then
        disc_tot = disc_tot + dis_perp(j,i)*area(j,i)
     end if

  end do
  end do

  disc_tot = disc_tot *year_sec *(rho/rho_w)
                  ! m^3/s ice equiv. -> m^3/a water equiv.

  disc_target = disc_target*1.0e+12_dp/rho_w   ! Gt/a -> m^3/a water equiv.

  write(6, fmt='(a)') ' '
  write(6, fmt='(a)') ' * * * * * * * * * * * * * * * * * * * * * * * * * * * * '

  write(6, fmt='(a)') ' '
  write(6, fmt='(3x,a,es12.4)') 'c_dis_0_init = ', c_dis_0

  c_dis_0 = c_dis_0 * disc_target/disc_tot

  write(6, fmt='(a)') ' '
  write(6, fmt='(3x,a,es12.4)') 'disc_target  = ', disc_target
  write(6, fmt='(3x,a,es12.4)') 'disc_tot     = ', disc_tot

  write(6, fmt='(a)') ' '
  write(6, fmt='(3x,a,es12.4)') '--> c_dis_0  = ', c_dis_0

  end subroutine calc_c_dis_0

!-------------------------------------------------------------------------------
!> Dependence of ice discharge coefficient on latitude (Greenland).
!! [Determine dependence of ice discharge coefficient on latitude.
!! This can be improved. For now I recommend s_dis=1.]
!<------------------------------------------------------------------------------
  subroutine disc_fields()

  ! Author: Reinhard Calov
  ! Institution: Potsdam Institute for Climate Impact Research  
  ! Date: 8.6.16    

  !   Purpose: Assign fields relevant for ice discharge parameterization
  !            1. Discharge coefficient field using linear dependence on 
  !               latitude.
  !            2. Subsurface temperature dependence on longitude and latitude.

  implicit none

  integer(i4b) :: i, j
  real(dp)     :: delta_phi=25.0_dp ! approximately the phi
                                    ! spanning the middle of domain

  c_dis = c_dis_0*c_dis_fac &
                 *(1.0_dp-(1.0_dp-s_dis)*(phi*pi_180_inv-60.0_dp)/delta_phi)

  c_dis = max(c_dis, 0.0_dp)

  t_sub_pd = 3.0_dp ! Can depend on lambda, phi later on

  end subroutine disc_fields

!-------------------------------------------------------------------------------
!> Ice discharge parameterization main formula, controler (general).
!! [Compute ice discharge via a parameterization using distance of ice 
!! margin to coast and ice thickness as parameters.]
!<------------------------------------------------------------------------------
  subroutine discharge(dxi, deta)

  ! Authors: Reinhard Calov, Andrey Ganopolski
  ! Institution: Potsdam Institute for Climate Impact Research  
  ! Date: 11.6.16

  implicit none

  real(dp), intent(in)    :: dxi, deta

  real(dp), parameter :: alpha_tc=60.0_dp ! maximal allowed angle
  real(dp) :: cos_tc

  n_discharge_call = n_discharge_call + 1

  cos_tc=dcos(alpha_tc*pi_180)

  if ((mod(n_discharge_call, iter_mar_coa)==0).or.(n_discharge_call==1)) then
     write(6, fmt='(10x,a)') 'Computation of mask_mar, cst_dist, cos_grad_tc'
     call marginal_ring(dxi, deta)
     call coastal_distance(dxi, deta)
  end if

  if(disc.ge.1) then !------------------ disc >= 1
    where(mask_mar.eq.1.and.cos_grad_tc.ge.cos_tc) 
      dis_perp=c_dis*(h_c+h_t)**m_h/cst_dist**m_d
    elsewhere
      dis_perp=0.0_dp
    end where
    if(disc.eq.2) then !----------------- disc = 2
      call calc_sub_oc_dt(t_sub_pd, dt_glann, dt_sub)
      dis_perp=dis_perp*(1.0_dp+alpha_sub*dt_sub) ! actual discharge present-day one corrected 
                                                ! via sub-ocean temperature anomaly
    end if
    dis_perp=max(0.0_dp, dis_perp) ! ensure positive values
  else if(disc.eq.0) then !-------------- disc = 0
    dis_perp=0.0_dp
  end if

  end subroutine discharge

!-------------------------------------------------------------------------------
!> Distance to the coast (general).
!! [Compute distance to the coast for every land point.]
!<------------------------------------------------------------------------------
  subroutine coastal_distance(dxi, deta)

  ! Author: Reinhard Calov
  ! Institution: Potsdam Institute for Climate Impact Research  
  ! Date: 2.11.11
  ! Methode: Options are:
  !          i_search=1: Brute force.
  !          i_search=2: Defining a square with two grid distances side length around the 
  !                      actual grid point and enlanging the square successively by
  !                      one grid. Because this should have been a circle, the
  !                      square is enlarged again if the coast line was found; this time 
  !                      enlargement is by a rectangles a the respective four side of 
  !                      the square. Smaller side lenght equals the rectangle root of 2 of the
  !                      larger side. 

  implicit none

  real(dp), intent(in) :: dxi, deta

  ! --------------------------------------------------
  real(dp), parameter :: grid_dist=5000.0_dp
  real(dp), parameter :: c_smooth=0.2_dp
  integer(i4b), parameter :: i_search=2
  ! --------------------------------------------------
  integer(i4b) :: i, j, i_pos, j_pos, l, l_e, d_l
  real(dp) :: dxi_inv, deta_inv
  real(dp) :: cst_dist_tmp, qrt_1, qrt_2
  logical :: leave_loop

  real(dp), allocatable, dimension(:,:) :: zs_tmp, grad_zs_x, grad_zs_y, &
  grad_dist_x, grad_dist_y

  allocate(zs_tmp(0:jmax,0:imax), grad_zs_x(0:jmax,0:imax), grad_zs_y(0:jmax,0:imax), &
  grad_dist_x(0:jmax,0:imax), grad_dist_y(0:jmax,0:imax))

    select case (i_search)

    case(1)
    ! brute force computation of minimal distance to coast for every land point

      do i_pos=0, imax
      do j_pos=0, jmax
        if(maske(j_pos,i_pos).ne.2) then
          cst_dist(j_pos,i_pos)=1.d+20
          do i=0, imax
          do j=0, jmax
            if(maske(j,i).eq.2) then
              cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
              if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
                cst_dist(j_pos,i_pos)=cst_dist_tmp
              end if
            end if
          end do
          end do
        else
          cst_dist(j_pos,i_pos)=grid_dist
        end if
      end do
      end do

    case(2)

      do i_pos=0, imax
      do j_pos=0, jmax
        if(maske(j_pos,i_pos).ne.2) then
          leave_loop=.false.
          cst_dist(j_pos,i_pos)=1.d+20
          do l=1, max(imax,jmax)
            do i=max(i_pos-l,0), min(i_pos+l,imax)
              j=max(j_pos-l,0); j=min(j,jmax)
              if(maske(j,i).eq.2) then
                leave_loop=.true.
                cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
                if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
                  cst_dist(j_pos,i_pos)=cst_dist_tmp
                end if
              end if
              j=min(j_pos+l,jmax)
              if(maske(j,i).eq.2) then
                leave_loop=.true.
                cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
                if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
                  cst_dist(j_pos,i_pos)=cst_dist_tmp
                end if
              end if
            end do
            do j=max(j_pos-l+1,0), min(j_pos+l-1,jmax)
              i=max(i_pos-l,0); i=min(i,imax)
              if(maske(j,i).eq.2) then
                leave_loop=.true.
                cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
                if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
                  cst_dist(j_pos,i_pos)=cst_dist_tmp
                end if
              end if
              i=min(i_pos+l,imax)
              if(maske(j,i).eq.2) then
                leave_loop=.true.
                cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
                if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
                  cst_dist(j_pos,i_pos)=cst_dist_tmp
                end if
              end if
            end do
            if(leave_loop) then
              l_e=l
              d_l=ceiling(l_e*sqrt(2.0_dp)+0.001_dp)
              exit   ! leave loop in l
            end if
          end do
          ! left
          do i=max(i_pos-(l_e+d_l),0),min(i_pos-(l_e-1),imax)
          do j=max(j_pos-l_e,0),min(j_pos+l_e,jmax)
            if(maske(j,i).eq.2) then
              cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
              if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
                cst_dist(j_pos,i_pos)=cst_dist_tmp
              end if
            end if
          end do
          end do
          ! right
          do i=max(i_pos+l_e+1,0),min(i_pos+l_e+d_l,imax)
          do j=max(j_pos-l_e,0),min(j_pos+l_e,jmax)
            if(maske(j,i).eq.2) then
              cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
              if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
                cst_dist(j_pos,i_pos)=cst_dist_tmp
              end if
            end if
          end do
          end do
          ! lower
          do i=max(i_pos-l_e,0),min(i_pos+l_e,imax)
          do j=max(j_pos-(l_e+d_l),0),min(j_pos-(l_e-1),jmax)
            if(maske(j,i).eq.2) then
              cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
              if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
                cst_dist(j_pos,i_pos)=cst_dist_tmp
              end if
            end if
          end do
          end do
          ! upper
          do i=max(i_pos-l_e,0),min(i_pos+l_e,imax)
          do j=max(j_pos+l_e+1,0),min(j_pos+l_e+d_l,jmax)
            if(maske(j,i).eq.2) then
              cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
              if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
                cst_dist(j_pos,i_pos)=cst_dist_tmp
              end if
            end if
          end do
          end do
        else
          cst_dist(j_pos,i_pos)=grid_dist
        end if
      end do
      end do

    end select

  ! test GIS only

    if(.false.) then
      do i=0, imax
      do j=0, jmax
        if(xi(i).le.-350000.0_dp &
          .and.-2800000.0_dp.le.eta(j) &
          .and.eta(j).le.-2250000.0_dp &
          .or. &
          xi(i).ge.100000.0_dp &
          .and.-2280000.0_dp.le.eta(j) &
          .and.eta(j).le.-1880000.0_dp &
          .or. &
          eta(j).ge.-1250000.0_dp &
          .and.0.0_dp.le.xi(i) &
          .and.xi(i).le.230000.0_dp) then
             cst_dist(j,i)=0.3_dp*cst_dist(j,i)
        end if
      end do
      end do
    end if

  ! forbid zero, something like grid distance is reasonable here.
  ! Dependents of coordinate system, for spherical coordinate system: 
  ! take length on Earth's surface. For now, simply just 5 km are taken.

    cst_dist=max(grid_dist, cst_dist)

  ! forbit to large values
    cst_dist=min(1.0e10_dp, cst_dist)

    zs_tmp=zs

  ! smoothing the surface topography
    do i=1, imax-1
    do j=1, jmax-1
       zs_tmp(j,i) = c_smooth*(zs(j,i-1)+zs(j,i+1)+zs(j-1,i)+zs(j+1,i)) &
                     +(1.0_dp-4.0_dp*c_smooth)*zs(j,i)
    end do
    end do

  ! gradients (2. order for now ...)

    dxi_inv  = 1.0_dp/dxi
    deta_inv = 1.0_dp/deta

!  ------ x-derivatives

    do i=1, imax-1
    do j=0, jmax
      grad_zs_x(j,i)    = (zs_tmp(j,i+1)-zs_tmp(j,i-1))*0.5_dp*dxi_inv*insq_g11_g(j,i)
      grad_dist_x(j,i)  = (cst_dist(j,i+1)-cst_dist(j,i-1))*0.5_dp*dxi_inv*insq_g11_g(j,i)
    end do
    end do

    do j=0, jmax
      grad_zs_x(j,0)       = (zs_tmp(j,1)-zs_tmp(j,0))*dxi_inv*insq_g11_g(j,0)
      grad_zs_x(j,imax)    = (zs_tmp(j,imax)-zs_tmp(j,imax-1))*dxi_inv*insq_g11_g(j,imax)
      grad_dist_x(j,0)     = (cst_dist(j,1)-cst_dist(j,0))*dxi_inv*insq_g11_g(j,0)
      grad_dist_x(j,imax)  = (cst_dist(j,imax)-cst_dist(j,imax-1))*dxi_inv*insq_g11_g(j,imax)
    end do

!  ------ y-derivatives

    do i=0, imax
    do j=1, jmax-1
      grad_zs_y(j,i)    = (zs_tmp(j+1,i)-zs_tmp(j-1,i))*0.5_dp*deta_inv*insq_g22_g(j,i)
      grad_dist_y(j,i)  = (cst_dist(j+1,i)-cst_dist(j-1,i))*0.5_dp*deta_inv*insq_g22_g(j,i)
    end do
    end do

    do i=0, imax
      grad_zs_y(0,i)       = (zs_tmp(1,i)-zs_tmp(0,i))*deta_inv*insq_g22_g(0,i)
      grad_zs_y(jmax,i)    = (zs_tmp(jmax,i)-zs_tmp(jmax-1,i))*deta_inv*insq_g22_g(jmax,i)
      grad_dist_y(0,i)     = (cst_dist(1,i)-cst_dist(0,i))*deta_inv*insq_g22_g(0,i)
      grad_dist_y(jmax,i)  = (cst_dist(jmax,i)-cst_dist(jmax-1,i))*deta_inv*insq_g22_g(jmax,i)
    end do

  ! angle between topographical gradient and coastal distance gradient

    do i=0, imax
    do j=0, jmax
      qrt_1=grad_zs_x(j,i)*grad_zs_x(j,i)+grad_zs_y(j,i)*grad_zs_y(j,i)
      qrt_2=grad_dist_x(j,i)*grad_dist_x(j,i)+grad_dist_y(j,i)*grad_dist_y(j,i)
      if(qrt_1.ne.0.0_dp.and.qrt_2.ne.0.0_dp) then
  !!! top_cst_alpha(j,i)=dacos( (grad_zs_x(j,i)*grad_dist_x(j,i)+grad_zs_y(j,i)*grad_dist_y(j,i))/ &
  !!! (sqrt(qrt_1)*sqrt(qrt_2)) )

        cos_grad_tc(j,i)= (grad_zs_x(j,i)*grad_dist_x(j,i)+grad_zs_y(j,i)*grad_dist_y(j,i))/ &
        (sqrt(qrt_1)*sqrt(qrt_2))
      else
        cos_grad_tc(j,i)=-1.0_dp
      end if
    end do
    end do

    deallocate(zs_tmp, grad_zs_x, grad_zs_y, grad_dist_x, grad_dist_y)

  end subroutine coastal_distance

!-------------------------------------------------------------------------------
!> Ring along an ice sheet margin (general).
!! [Compute marginal ring.]
!<------------------------------------------------------------------------------
  subroutine marginal_ring(dxi, deta)

  ! Author: Reinhard Calov
  ! Institution: Potsdam Institute for Climate Impact Research  
  ! Date: 28.10.11

  ! Purpose: Determine an r_max_eff wide ring arong ice margins 
  ! towards the interior of the ice sheet. Small stripes of land are
  ! not considered as land. For the logics De Morgan is used often.

  ! Methode: Two staggered loops in i,j. The inner i,j loop acts inside
  ! a rectangle defined by r_mar_eff. r_mar_eff should be small compared
  ! to the domain size. 

    implicit none

    real(dp), intent(in) :: dxi, deta

    integer(i4b) :: i, j, i_pos, j_pos, di_eff, dj_eff
    integer(i4b) :: count_tmp
    real(dp) :: r_p

    if(r_mar_eff.le.1.0e6_dp) then
      mask_mar=0
      do i_pos=1, imax-1
      do j_pos=1, jmax-1
        if((maske(j_pos,i_pos).eq.1.or.maske(j_pos,i_pos).eq.2).and.     &
           .not.(maske(j_pos,i_pos).eq.1.and.maske(j_pos,i_pos-1).eq.0   &
                 .and.maske(j_pos,i_pos+1).eq.0.or. &
                  maske(j_pos,i_pos).eq.1.and.maske(j_pos-1,i_pos).eq.0  &
                  .and.maske(j_pos+1,i_pos).eq.0.or. &
                   maske(j_pos,i_pos).eq.1.and.maske(j_pos,i_pos-1).eq.3 &
                   .and.maske(j_pos,i_pos+1).eq.3.or. &
                   maske(j_pos,i_pos).eq.1.and.maske(j_pos-1,i_pos).eq.3 &
                   .and.maske(j_pos+1,i_pos).eq.3.or. &
                   maske(j_pos,i_pos).eq.1.and.maske(j_pos,i_pos-1).eq.3 &
                   .and.maske(j_pos,i_pos+1).eq.0.or. &
                   maske(j_pos,i_pos).eq.1.and.maske(j_pos-1,i_pos).eq.0 &
                   .and.maske(j_pos+1,i_pos).eq.3)) then ! outside ice sheet, exclude isolated land stripes
          di_eff=int(r_mar_eff/dxi)+1; dj_eff=int(r_mar_eff/deta)+1 ! only for grid=0, 1 yet!
          do i=max(i_pos-di_eff,0), min(i_pos+di_eff,imax)
          do j=max(j_pos-dj_eff,0), min(j_pos+dj_eff,jmax)
            r_p=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
            if(r_p.le.r_mar_eff.and.(maske(j,i).eq.0.or.maske(j,i).eq.3)) then
              mask_mar(j,i)=1
            end if
          end do
          end do
        end if
      end do
      end do
      where(maske.ge.1)
        mask_mar=1
      end where
    else ! the ring encompassed entire Greenland 
      mask_mar=1
    end if

  end subroutine marginal_ring

!-------------------------------------------------------------------------------
!> Anomaly of subsurface temperature (general).
!! [Compute anomaly of subsurface temperature with a simple parameterization.]
!<------------------------------------------------------------------------------
  subroutine calc_sub_oc_dt(T_sub_PD, dT_glann, dT_sub)

  ! Author: Reinhard Calov
  ! Institution: Potsdam Institute for Climate Impact Research  
  ! Purpose: Compute anomaly of subsurface temperature
  !          with a simple parameterization
  !          using global annual temperature
  !          anomaly (e.g. from CLIMBER)

  implicit none

  real(dp), intent(in)  :: T_sub_PD, dT_glann
  real(dp), intent(out) :: dT_sub

  real(dp) :: T_sub

  dt_sub=alpha_o*dt_glann                    ! Sub-ocean temperature anomaly
  t_sub=max(t_sea_freeze, t_sub_pd)+dt_sub   ! Actual sub-ocean temperature
  t_sub=max(t_sea_freeze, t_sub)             ! Ocean temperature should stay
                                             ! above freezing point
  dt_sub=t_sub-t_sub_pd

  end subroutine calc_sub_oc_dt

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

end module discharge_workers_m
!