borehole.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  b o r e h o l e
!
!> @file
!!
!! Computation of an arbitrary field quantity for a given borehole
!! position x_pos, y_pos by weighed averaging of the corresponding
!! gridded 2-d field.
!!
!! @section Copyright
!!
!! Copyright 2009-2013 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 an arbitrary field quantity for a given borehole
!! position x_pos, y_pos by weighed averaging of the corresponding
!! gridded 2-d field.
!<------------------------------------------------------------------------------
subroutine borehole(field, x_pos, y_pos, dxi, deta, ch_grid, field_val)

use sico_types
use sico_variables

implicit none

real(dp), dimension(0:JMAX,0:IMAX), intent(in) :: field
real(dp), intent(in) :: x_pos, y_pos, dxi, deta
character (len=*), intent(in) :: ch_grid

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

integer(i4b) :: i_kl, i_gr, j_kl, j_gr
real(dp) :: real_i, real_j, dist1, dist2, dist3, dist4, &
            weight1, weight2, weight3, weight4

!-------- Neighbour points, weighing factors --------

real_i = (x_pos-xi(0)) /dxi
real_j = (y_pos-eta(0))/deta

i_kl = floor(real_i)
j_kl = floor(real_j)
i_gr = ceiling(real_i)
j_gr = ceiling(real_j)

if (i_kl <= 1) stop ' borehole: i_kl <= 1 not allowed!'
if (j_kl <= 1) stop ' borehole: j_kl <= 1 not allowed!'
if (i_gr >= imax-1) stop ' borehole: i_gr >= IMAX-1 not allowed!'
if (j_gr >= jmax-1) stop ' borehole: j_gr >= JMAX-1 not allowed!'

dist1 = sqrt( (x_pos-xi(i_kl))**2 + (y_pos-eta(j_kl))**2 )
dist2 = sqrt( (x_pos-xi(i_kl))**2 + (y_pos-eta(j_gr))**2 )
dist3 = sqrt( (x_pos-xi(i_gr))**2 + (y_pos-eta(j_kl))**2 )
dist4 = sqrt( (x_pos-xi(i_gr))**2 + (y_pos-eta(j_gr))**2 )

weight1 = 1.0_dp/(dist1+eps)
weight2 = 1.0_dp/(dist2+eps)
weight3 = 1.0_dp/(dist3+eps)
weight4 = 1.0_dp/(dist4+eps)

!-------- Inverse-distance weighing of the four adjacent grid values --------

if (ch_grid=='grid') then   ! field(j,i) defined on grid points

   field_val = ( weight1*field(j_kl,i_kl) &
                +weight2*field(j_gr,i_kl) &
                +weight3*field(j_kl,i_gr) &
                +weight4*field(j_gr,i_gr) ) &
               /(weight1 + weight2 + weight3 + weight4)

else if (ch_grid=='sg_x') then
                          ! field(j,i) defined on staggered grid in x direction

   field_val = ( weight1*0.5_dp*(field(j_kl,i_kl)+field(j_kl,i_kl-1)) &
                +weight2*0.5_dp*(field(j_gr,i_kl)+field(j_gr,i_kl-1)) &
                +weight3*0.5_dp*(field(j_kl,i_gr)+field(j_kl,i_gr-1)) &
                +weight4*0.5_dp*(field(j_gr,i_gr)+field(j_gr,i_gr-1)) ) &
               /(weight1 + weight2 + weight3 + weight4)

else if (ch_grid=='sg_y') then
                          ! field(j,i) defined on staggered grid in y direction

   field_val = ( weight1*0.5_dp*(field(j_kl,i_kl)+field(j_kl-1,i_kl)) &
                +weight2*0.5_dp*(field(j_gr,i_kl)+field(j_gr-1,i_kl)) &
                +weight3*0.5_dp*(field(j_kl,i_gr)+field(j_kl-1,i_gr)) &
                +weight4*0.5_dp*(field(j_gr,i_gr)+field(j_gr-1,i_gr)) ) &
               /(weight1 + weight2 + weight3 + weight4)

else

   stop ' borehole: Parameter ch_grid has undefined value!'

end if

end subroutine borehole
!