ufo_gnssro_bndnbam_util_mod.F90

Go to the documentation of this file.

!  
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 
 
!> Fortran module of Gnssro NBAM (NCEP's Bending Angle Method) operator
 
module ufo_gnssro_bndnbam_util_mod
  use fckit_log_module, only: fckit_log
  use kinds
  use missing_values_mod
  use gnssro_mod_constants
  use vert_interp_mod
  use lag_interp_mod,     only: lag_interp_const, lag_interp_smthweights
  use gnssro_mod_transform
  use gnssro_mod_grids, only: get_coordinate_value
 
  implicit none
  public             :: ufo_gnssro_bndnbam_simobs_single
  private
 
  contains
! ------------------------------------------------------------------------------
subroutine ufo_gnssro_bndnbam_simobs_single( &
           obsLat, obsGeoid, obsLocR, obsImpP, &
           grids, ngrd, &
           nlev, nlevExt, nlevAdd, nlevCheck, &
           radius,ref,refIndex,refXrad, &
           bendingAngle)
! -------------------------------------------------------------------------------
  character(len=*), parameter    :: myname  = "ufo_gnssro_bndnbam_simobs_single"
 
  real(kind_real), intent(out)   :: bendingangle
 
  integer, intent(in)            :: nlev
  integer, intent(in)            :: nlevext
  integer, intent(in)            :: nlevadd
  integer, intent(in)            :: nlevcheck
  integer, intent(in)            :: ngrd
  real(kind_real), intent(in)    :: obslat, obsgeoid, obslocr, obsimpp ! obsspace
  real(kind_real), intent(in)    :: grids(ngrd)
 
  real(kind_real), intent(in)    :: radius(nlev)
  real(kind_real), intent(in)    :: refindex(nlev)
  real(kind_real), intent(inout) :: ref(nlevext)
  real(kind_real), intent(inout) :: refxrad(0:nlevext+1)
  real(kind_real)                :: lagconst(3,nlevext)
 
  real(kind_real)                 :: sindx
  real(kind_real)                 :: obsimph
  integer                         :: k, igrd, indx
  real(kind_real)                 :: d_refxrad
  real(kind_real)                 :: w4(4), dw4(4)
  real(kind_real)                 :: bndintgd
  real(kind_real)                 :: rnlevext
  real(kind_real)                 :: derivref_s(ngrd)
  real(kind_real)                 :: refxrad_s(ngrd)
 
!------------------------------------------------------------
! Extend atmosphere above interface level nlev
  d_refxrad = refxrad(nlev) - refxrad(nlev-1)
  do k = 1, nlevadd
     refxrad(nlev+k)=refxrad(nlev)+ k*d_refxrad    ! extended x_i
     ref(nlev+k)=ref(nlev+k-1)**2/ref(nlev+k-2) ! exended N_i
  end do
 
  refxrad(0)=refxrad(3)
  refxrad(nlevext+1)=refxrad(nlevext-2)
  do k = 1,nlevext
     call lag_interp_const(lagconst(:,k),refxrad(k-1:k+1),3)
  enddo
 
! integrate on a new set of equally-spaced vertical grid 
  derivref_s = zero
  grids_loop: do igrd =1,ngrd
    refxrad_s(igrd)=sqrt(grids(igrd)**2 + obsimpp**2) !x_s^2=s^2+a^2
    call get_coordinate_value(refxrad_s(igrd), sindx,refxrad(1:nlevext),nlevext,"increasing")
    rnlevext = float(nlevext)
 
    if (sindx > zero .and. sindx < rnlevext) then  !obs inside the new grid
       indx=sindx
!      Compute derivative at new grids (dN/ds) using Lagrange interpolators
       call lag_interp_smthweights(refxrad(indx-1:indx+2),refxrad_s(igrd),&
                                   lagconst(:,indx),lagconst(:,indx+1),   &
                                   w4,dw4,4)
       if (indx==1) then
          w4(4)=w4(4)+w4(1); w4(1:3)=w4(2:4);w4(4)=zero
          dw4(4)=dw4(4)+dw4(1);dw4(1:3)=dw4(2:4);dw4(4)=zero
          indx=indx+1
       endif
       if (indx==nlevext-1) then
          w4(1)=w4(1)+w4(4); w4(2:4)=w4(1:3);w4(1)=zero
          dw4(1)=dw4(1)+dw4(4); dw4(2:4)=dw4(1:3);dw4(1)=zero
          indx=indx-1
       endif
 
       derivref_s(igrd)=dot_product(dw4,ref(indx-1:indx+2)) !derivative dN/dx_s
       derivref_s(igrd)=max(zero,abs(derivref_s(igrd)))
 
    else
      return
    endif !obs in new grid
  end do grids_loop
 
! bending angle (radians)
  bendingangle = ds*derivref_s(1)/refxrad_s(1)
  do igrd = 2,ngrd
     bndintgd     = ds*derivref_s(igrd)/refxrad_s(igrd)
     bendingangle = bendingangle + two*bndintgd
  end do
  bendingangle=r1em6 * obsimpp * bendingangle
 
end subroutine ufo_gnssro_bndnbam_simobs_single
! ------------------------------------------------------------------------------
end module ufo_gnssro_bndnbam_util_mod