ufo_gnssro_bndropp2d_mod.F90

Go to the documentation of this file.

! (C) Copyright 2017-2018 UCAR
! 
! 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 for gnssro bending angle ropp2d forward operator
!> following the ROPP (2018 Aug) implementation
 
module ufo_gnssro_bndropp2d_mod
 
use fckit_configuration_module, only: fckit_configuration 
!use iso_c_binding
use kinds
use ufo_vars_mod
use ufo_geovals_mod
use ufo_geovals_mod_c, only: ufo_geovals_registry
use ufo_basis_mod,     only: ufo_basis
use vert_interp_mod
use lag_interp_mod,    only: lag_interp_const, lag_interp_smthweights
use obsspace_mod
use missing_values_mod
use ufo_gnssro_ropp2d_utils_mod
use ufo_gnssro_ropp1d_utils_mod
 
use gnssro_mod_conf
use fckit_log_module,  only : fckit_log
 
implicit none
public             :: ufo_gnssro_bndropp2d
private
  !> Fortran derived type for gnssro trajectory
type, extends(ufo_basis) :: ufo_gnssro_bndropp2d
  type(gnssro_conf)  :: roconf
  real(kind_real), allocatable  :: obslon2d(:), obslat2d(:)  !2d location
  contains
 
    procedure :: setup     => ufo_gnssro_bndropp2d_setup
    procedure :: simobs    => ufo_gnssro_bndropp2d_simobs
end type ufo_gnssro_bndropp2d
 
contains
 
! ------------------------------------------------------------------------------
subroutine ufo_gnssro_bndropp2d_setup(self, f_conf, c_size)
  implicit none
  class(ufo_gnssro_bndropp2d), intent(inout) :: self
  type(fckit_configuration),   intent(in)    :: f_conf
  integer,                     intent(in)    :: c_size ! 1d obsspace vector length
 
  call gnssro_conf_setup(self%roconf,f_conf)
 
  allocate(self%obsLon2d(c_size*self%roconf%n_horiz))
  allocate(self%obsLat2d(c_size*self%roconf%n_horiz))
 
  self%obsLon2d = 0.0
  self%obsLat2d = 0.0
 
end subroutine ufo_gnssro_bndropp2d_setup
 
! ------------------------------------------------------------------------------
subroutine ufo_gnssro_bndropp2d_simobs(self, geovals, hofx, obss)
  use ropp_fm_types, only: state2dfm, state1dfm
  use ropp_fm_types, only: obs1dbangle
  use typesizes,     only: wp => eightbytereal
  use datetimetypes, only: dp
 
  implicit none
  class(ufo_gnssro_bndropp2d), intent(in)    :: self
  type(ufo_geovals),           intent(in)    :: geovals
  real(kind_real),             intent(inout) :: hofx(:)
  type(c_ptr), value,          intent(in)    :: obss
  real(c_double)                             :: missing
 
  type(state2dfm)                    :: x
  type(state1dfm)                    :: x1d
  type(obs1dbangle)                  :: y
 
  character(len=*), parameter   :: myname_="ufo_gnssro_bndropp2d_simobs"
  integer, parameter            :: max_string = 800
  character(max_string)         :: err_msg
  integer                       :: nlev, nlocs, iobs, nvprof
  integer                       :: iflip
  type(ufo_geoval), pointer     :: t, q, prs, gph, gph_sfc
  real(kind_real), allocatable  :: obsImpP(:),obsLocR(:),obsGeoid(:),obsAzim(:) !nlocs
  real(kind_real), allocatable  :: obsLat(:),obsLon(:)                          !nlocs
  real(kind_real), allocatable  :: obsLonnh(:),obsLatnh(:)                      !n_horiz
  integer                       :: n_horiz
  real(kind_real)               :: dtheta
  real(kind_real)                    :: ob_time
  integer, allocatable, dimension(:) :: ichk
 
  n_horiz = self%roconf%n_horiz
  dtheta  = self%roconf%dtheta
 
  write(err_msg,*) "TRACE: ufo_gnssro_bndropp2d_simobs: begin"
  call fckit_log%info(err_msg)
! check if nlocs is consistent in geovals & hofx
  if (geovals%nlocs /= size(hofx)*n_horiz) then
      write(err_msg,*) myname_, ' error: 2d nlocs inconsistent! geovals%nlocs, size(hofx), &
                                  and n_horiz are', geovals%nlocs, size(hofx), n_horiz
      call abor1_ftn(err_msg)
  endif
 
! get variables from geovals
  call ufo_geovals_get_var(geovals, var_ts,    t)         ! temperature
  call ufo_geovals_get_var(geovals, var_q,     q)         ! specific humidity
  call ufo_geovals_get_var(geovals, var_prs,   prs)       ! pressure
  call ufo_geovals_get_var(geovals, var_z,     gph)       ! geopotential height
  call ufo_geovals_get_var(geovals, var_sfc_geomz, gph_sfc)   ! surface geopotential height
 
  missing = missing_value(missing)
  nlev    = t%nval ! number of model levels
  nlocs   = obsspace_get_nlocs(obss)
 
  iflip = 0
  if (prs%vals(1,1) .lt. prs%vals(prs%nval,1) ) then
    iflip = 1
    write(err_msg,'(a)') '  ufo_gnssro_bndropp2d_simobs:'//new_line('a')//                         &
                         '  Model vertical height profile is in descending order,'//new_line('a')// &
                         '  but ROPP requires it to be ascending order, need flip'
    call fckit_log%info(err_msg)
  end if
 
! set obs space struture
  allocate(obslon(nlocs))
  allocate(obslat(nlocs))
  allocate(obsimpp(nlocs))
  allocate(obslocr(nlocs))
  allocate(obsgeoid(nlocs))
  allocate(obsazim(nlocs))
  allocate(obslatnh(n_horiz))
  allocate(obslonnh(n_horiz))
 
  call obsspace_get_db(obss, "MetaData", "longitude",        obslon)
  call obsspace_get_db(obss, "MetaData", "latitude",         obslat)
  call obsspace_get_db(obss, "MetaData", "impact_parameter", obsimpp)
  call obsspace_get_db(obss, "MetaData", "earth_radius_of_curvature", obslocr)
  call obsspace_get_db(obss, "MetaData", "geoid_height_above_reference_ellipsoid", obsgeoid)
  call obsspace_get_db(obss, "MetaData", "sensor_azimuth_angle", obsazim)
 
  nvprof  = 1  ! no. of bending angles in profile
  ob_time = 0.0
  allocate(ichk(nvprof))
  ichk(:) = 0
  write(err_msg,*) "TRACE: ufo_gnssro_bndropp2d_simobs: begin observation loop, nlocs =  ", nlocs
  call fckit_log%info(err_msg)
 
! loop through the obs
  obs_loop: do iobs = 1, nlocs  
 
    if ( ( obsimpp(iobs)-obslocr(iobs)-obsgeoid(iobs) ) <= self%roconf%top_2d .and. &
           obsazim(iobs) /= missing ) then
 
      obslatnh = self%obsLat2d( (iobs-1)*n_horiz+1:iobs*n_horiz )
      obslonnh = self%obsLon2d( (iobs-1)*n_horiz+1:iobs*n_horiz )
      call init_ropp_2d_statevec(obslonnh, obslatnh,                  &
                      t%vals(:,(iobs-1)*n_horiz+1:iobs*n_horiz),      &
                      q%vals(:,(iobs-1)*n_horiz+1:iobs*n_horiz),      &
                    prs%vals(:,(iobs-1)*n_horiz+1:iobs*n_horiz),      &
                    gph%vals(:,(iobs-1)*n_horiz+1:iobs*n_horiz),      &
                    nlev,x,n_horiz,dtheta,iflip)
     
      call init_ropp_2d_obvec(nvprof,      &
                    obsimpp(iobs),         &
                     obslat(iobs),         &
                     obslon(iobs),         &
                    obslocr(iobs),         &
                   obsgeoid(iobs),         &
                             y)
 
      call ropp_fm_bangle_2d(x,y)
 
    else ! apply ropp1d above top_2d
 
      call init_ropp_1d_statevec(ob_time,            &
                               obslon(iobs),         &
                               obslat(iobs),         &
                               t%vals(:,(iobs-1)*n_horiz+1+(n_horiz-1)/2),       &
                               q%vals(:,(iobs-1)*n_horiz+1+(n_horiz-1)/2),       &
                               prs%vals(:,(iobs-1)*n_horiz+1+(n_horiz-1)/2),     &
                               gph%vals(:,(iobs-1)*n_horiz+1+(n_horiz-1)/2),     &
                               nlev,                                             &
                               gph_sfc%vals(1,(iobs-1)*n_horiz+1+(n_horiz-1)/2), &
                               x1d, iflip)
 
      call init_ropp_1d_obvec(nvprof,          &
                              obsimpp(iobs),   &
                              ichk, ob_time,   &
                              obslat(iobs),    &
                              obslon(iobs),    &
                              obslocr(iobs),   &
                              obsgeoid(iobs),  &
                              y)
 
      call ropp_fm_bangle_1d(x1d,y)
 
    end if
 
!   hack -- handling ropp missing value 
    if (y%bangle(nvprof) .lt. -900.0_wp ) then
       hofx(iobs) = missing
       y%bangle(nvprof) = missing
    else
       hofx(iobs) = y%bangle(nvprof)  ! nvprof is just one point
    end if
!   hack -- handling ropp missing value 
 
!   deallocate ropp structures  
    if ( ( obsimpp(iobs)-obslocr(iobs)-obsgeoid(iobs) ) <= self%roconf%top_2d ) then
      call ropp_tidy_up_2d(x,y)
    else
      call ropp_tidy_up_1d(x1d,y)
    end if
 
  end do obs_loop
 
  deallocate(obslat)
  deallocate(obslon)
  deallocate(obsimpp)
  deallocate(obslocr)
  deallocate(obsgeoid)
  deallocate(obsazim)
  deallocate(obslatnh)
  deallocate(obslonnh)
  deallocate(ichk)
 
  write(err_msg,*) "TRACE: ufo_gnssro_bndropp2d_simobs: completed"
  call fckit_log%info(err_msg)
     
  return
 
end subroutine ufo_gnssro_bndropp2d_simobs
! ------------------------------------------------------------------------------
 
end module ufo_gnssro_bndropp2d_mod