ufo_gnssro_bndnbam_tlad_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 of Gnssro NBAM (NCEP's Bending Angle Method)
!> tlad operator
 
module ufo_gnssro_bndnbam_tlad_mod
use fckit_configuration_module, only: fckit_configuration 
use kinds
use ufo_vars_mod
use ufo_geovals_mod
use ufo_geovals_mod_c, only: ufo_geovals_registry
use obsspace_mod
use missing_values_mod
use lag_interp_mod
use ufo_basis_tlad_mod,only: ufo_basis_tlad
use gnssro_mod_conf
use gnssro_mod_constants
use fckit_log_module,  only : fckit_log
 
implicit none
real(c_double)                             :: missing
 
!>  Fortran derived type for gnssro trajectory
type, extends(ufo_basis_tlad) :: ufo_gnssro_bndnbam_tlad
  private
  integer                       :: nlev, nlev1, nlocs, iflip, nrecs
  real(kind_real), allocatable  :: jac_t(:,:), jac_prs(:,:), jac_q(:,:)
  integer, allocatable          :: nlocs_begin(:), nlocs_end(:)
  type(gnssro_conf)             :: roconf
 
  contains
    procedure :: setup      => ufo_gnssro_bndnbam_tlad_setup
    procedure :: delete     => ufo_gnssro_bndnbam_tlad_delete
    procedure :: settraj    => ufo_gnssro_bndnbam_tlad_settraj
    procedure :: simobs_tl  => ufo_gnssro_bndnbam_simobs_tl
    procedure :: simobs_ad  => ufo_gnssro_bndnbam_simobs_ad
end type ufo_gnssro_bndnbam_tlad
 
contains
! ------------------------------------------------------------------------------
subroutine ufo_gnssro_bndnbam_tlad_setup(self, f_conf)
  implicit none
  class(ufo_gnssro_bndnbam_tlad), intent(inout) :: self
  type(fckit_configuration), intent(in)        :: f_conf
 
  call gnssro_conf_setup(self%roconf,f_conf)
 
end subroutine ufo_gnssro_bndnbam_tlad_setup
 
! ------------------------------------------------------------------------------
subroutine ufo_gnssro_bndnbam_tlad_settraj(self, geovals, obss)
  use gnssro_mod_transform
  use gnssro_mod_grids, only: get_coordinate_value
  use, intrinsic::  iso_c_binding
  implicit none
  class(ufo_gnssro_bndnbam_tlad), intent(inout) :: self
  type(ufo_geovals),             intent(in)    :: geovals
  type(c_ptr), value,            intent(in)    :: obss
  character(len=*), parameter     :: myname  = "ufo_gnssro_bndnbam_tlad_settraj"
  character(max_string)           :: err_msg
  integer, parameter              :: nlevAdd = 13 !num of additional levels on top of exsiting model levels
  integer, parameter              :: newAdd  = 20 !num of additional levels on top of extended levels
  integer, parameter              :: ngrd    = 80 !num of new veritcal grids for bending angle computation
  type(ufo_geoval), pointer       :: t, q, gph, prs, zs
  integer                         :: iobs,k,j, klev, irec, icount
  integer                         :: nrecs
  integer                         :: nlev, nlev1, nlocs, nlevExt, nlevCheck
  real(kind_real)                 :: dw4(4), dw4_tl(4)
  real(kind_real)                 :: geomzi
  real(kind_real), allocatable    :: obsLat(:), obsImpP(:), obsLocR(:), obsGeoid(:)
  integer(c_size_t), allocatable  :: obsRecnum(:)
  real(kind_real)                 :: d_refXrad, gradRef
  real(kind_real)                 :: d_refXrad_tl
  real(kind_real)                 :: grids(ngrd)
  real(kind_real)                 :: sIndx 
  integer                         :: indx
  real(kind_real)                 :: p_coef, t_coef, q_coef
  real(kind_real)                 :: fv, pw
  real(kind_real)                 :: dbetaxi, dbetan
  real(kind_real), allocatable    :: lagConst(:,:), lagConst_tl(:,:)
  real(kind_real), allocatable    :: gesT(:,:), gesQ(:,:), gesP(:,:), gesH(:,:), gesZs(:)                     
 
  real(kind_real),allocatable     :: radius(:), dzdh(:), refIndex(:)
  real(kind_real), allocatable    :: dhdp(:), dhdt(:)
  real(kind_real), allocatable    :: ref(:)
  real(kind_real), allocatable    :: refXrad(:)
  real(kind_real), allocatable    :: refXrad_s(:)
  real(kind_real), allocatable    :: refXrad_tl(:), ref_tl(:)
  real(kind_real), allocatable    :: dndp(:,:), dndt(:,:), dndq(:,:)
  real(kind_real), allocatable    :: dxidp(:,:), dxidt(:,:), dxidq(:,:)
  real(kind_real), allocatable    :: dbenddxi(:), dbenddn(:)
  integer,         allocatable    :: super_refraction_flag(:)
  integer,         allocatable    :: obsSRflag(:)
  integer                         :: hasSRflag
 
  write(err_msg,*) myname, ": begin"
  call fckit_log%info(err_msg)
 
! Make sure nothing already allocated   
  call self%delete()
 
  missing = missing_value(missing)
  nlocs   = obsspace_get_nlocs(obss) ! number of observations
  nrecs   = obsspace_get_nrecs(obss)
  write(err_msg,*) myname, ': nlocs from gelvals and hofx, nrecs', nlocs, nrecs
  call fckit_log%info(err_msg)
 
if (nlocs > 0 ) then
! get variables from geovals
  call ufo_geovals_get_var(geovals, var_ts,  t)         ! air temperature
  call ufo_geovals_get_var(geovals, var_q,   q)         ! specific humidity
  call ufo_geovals_get_var(geovals, var_sfc_geomz, zs)  ! surface geopotential height/surface altitude
 
  if (self%roconf%vertlayer .eq. "mass") then
    call ufo_geovals_get_var(geovals, var_prs,   prs)       ! pressure
    call ufo_geovals_get_var(geovals, var_z,     gph)       ! geopotential height
  else
    call ufo_geovals_get_var(geovals, var_prsi,  prs)       ! pressure
    call ufo_geovals_get_var(geovals, var_zi,    gph)
  end if
 
  nlev  = t%nval   ! number of model mass levels
  nlev1 = prs%nval
  self%nlev  = t%nval   ! number of model mass levels
  self%nlev1 = prs%nval ! number of model pressure/height levels 
  self%nlocs = nlocs
  self%nrecs = nrecs
 
  nlevext = nlev + nlevadd
  nlevcheck = min(23, nlev) !number of levels to check super refraction
 
  allocate(gest(nlev,nlocs))
  allocate(gesq(nlev,nlocs))
  allocate(gesp(nlev1,nlocs))
  allocate(gesh(nlev1,nlocs))
  allocate(geszs(nlocs))
 
! copy geovals to local background arrays
  self%iflip = 0
  if (prs%vals(1,1) .lt. prs%vals(prs%nval,1) ) then
     self%iflip = 1
     write(err_msg,'(a)')'  ufo_gnssro_bndnbam_tlad_settraj:'//new_line('a')//                         &
                         '  Model vertical height profile is in descending order,'//new_line('a')// &
                         '  but NBAM requires it to be ascending order, need flip'
    call fckit_log%info(err_msg)
    do k = 1, nlev
       gest(k,:) = t%vals(nlev-k+1,:)
       gesq(k,:) = q%vals(nlev-k+1,:)
    enddo
    do k = 1, nlev1
       gesp(k,:) = prs%vals(nlev1-k+1,:)
       gesh(k,:) = gph%vals(nlev1-k+1,:)
    enddo
  else  ! not flipping
    call fckit_log%info(err_msg)
 
    do k = 1, nlev
       gest(k,:) = t%vals(k,:)
       gesq(k,:) = q%vals(k,:)
    enddo
    do k = 1, nlev1
       gesp(k,:) = prs%vals(k,:)
       gesh(k,:) = gph%vals(k,:)
    enddo
  end if
 
! if all fields t and q are on mass layers,
! while p and z are on interface layers -- NBAM manner
  if ( nlev1 /= nlev ) then
     do k = nlev, 2, -1
        gest(k,:) = half* (gest(k,:) + gest(k-1,:))
        gesq(k,:) = half* (gesq(k,:) + gesq(k-1,:))
     enddo
  end if
        geszs(:) = zs%vals(1,:)
 
! set obs space struture
  allocate(obslat(nlocs))
  allocate(obsimpp(nlocs))
  allocate(obslocr(nlocs))
  allocate(obsgeoid(nlocs))
  allocate(obsrecnum(nlocs))
  allocate(self%nlocs_begin(nrecs))
  allocate(self%nlocs_end(nrecs))
 
  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_recnum(obss, obsrecnum)
 
  if (  obsspace_has(obss,  "SR_flag",  "bending_angle") ) then
    allocate(obssrflag(nlocs))
    call obsspace_get_db(obss,  "SR_flag", "bending_angle", obssrflag)
    hassrflag = 1 ! SR_flag generated in hofx for real case run
  else
    hassrflag = 0 ! SR_flag does not exist in ctest
  end if
 
  self%nlocs_begin=1
  self%nlocs_end=1
 
  icount = 1
  do iobs = 1, nlocs-1
    if (obsrecnum(iobs+1) /= obsrecnum(iobs)) then
      icount = icount +1  !counting number of records
      self%nlocs_end(icount-1)= iobs
      self%nlocs_begin(icount) = iobs+1
    end if
  end do
  self%nlocs_end(nrecs)= nlocs
 
  if (icount /= nrecs) then
    write(err_msg,*) "record number is not consistent :", icount, nrecs
    call fckit_log%info(err_msg)
  end if
 
  allocate(dhdp(nlev))
  allocate(dhdt(nlev))
  allocate(dzdh(nlev))
  allocate(radius(nlev))
  allocate(refindex(nlev))
  allocate(ref(nlevext)) 
  allocate(ref_tl(nlevext))
  allocate(refxrad(0:nlevext+1)) 
  allocate(refxrad_tl(0:nlevext+1))
  allocate(refxrad_s(ngrd))
  allocate(dndp(nlev,nlev))
  allocate(dndq(nlev,nlev))
  allocate(dndt(nlev,nlev))
  allocate(dxidp(nlev,nlev))
  allocate(dxidt(nlev,nlev))
  allocate(dxidq(nlev,nlev))
  allocate(dbenddxi(nlev))
  allocate(dbenddn(nlev))
  allocate(lagconst_tl(3,nlevext))
  allocate(lagconst(3,nlevext))
  allocate(self%jac_t(nlev,nlocs))
  allocate(self%jac_q(nlev,nlocs))
  allocate(self%jac_prs(nlev1,nlocs))
 
! Inizialize some variables
  dxidt=zero; dxidp=zero; dxidq=zero
  dndt=zero;  dndq=zero;  dndp=zero
 
! tempprary manner to handle the missing hofx 
  self%jac_t = missing
 
  do j = 1, ngrd
     grids(j) = (j-1) * ds
  end do
 
! calculate jacobian
  call gnssro_ref_constants(self%roconf%use_compress)
 
  iobs = 0
  rec_loop: do irec = 1, nrecs
    obs_loop: do icount = self%nlocs_begin(irec), self%nlocs_end(irec)
 
      iobs = iobs + 1
 
      if (hassrflag == 1) then
         if (obssrflag(iobs) > 0)  cycle obs_loop
      end if
 
      dxidt=zero; dxidp=zero; dxidq=zero
      dndt=zero;  dndq=zero;  dndp=zero
 
      do k = 1, nlev
!        geometric height nad dzdh jacobian
         call geop2geometric( obslat(iobs), gesh(k,iobs)-geszs(iobs), geomzi, dzdh(k))
!        guess radius 
         radius(k) = geomzi + geszs(iobs) + obsgeoid(iobs) + obslocr(iobs)   ! radius r
!        guess refactivity, refactivity index,  and impact parameter
         call compute_refractivity(gest(k,iobs), gesq(k,iobs), gesp(k,iobs), &
                                 ref(k),self%roconf%use_compress) 
         refindex(k)= one + (r1em6*ref(k)) 
         refxrad(k) = refindex(k) * radius(k) 
      end do
 
!     data rejection based on model background !
!     (1) skip data beyond model levels
      call get_coordinate_value(obsimpp(iobs),sindx,refxrad(1),nlev,"increasing")
      if (sindx < one .or. sindx > float(nlev)) then
         cycle obs_loop
      end if
 
      do k = 1, nlev
 
!       jacobian for refractivity(N)
        fv    = rv_over_rd-one
        pw    = rd_over_rv+gesq(k,iobs)*(one-rd_over_rv)
        q_coef =  n_b *gesp(k,iobs)/(gest(k,iobs)**2*pw**2)*rd_over_rv + &
                  n_c *gesp(k,iobs)/(gest(k,iobs)*  pw**2)*rd_over_rv
        p_coef =  n_a/gest(k,iobs)   + &
                  n_b*gesq(k,iobs)/(gest(k,iobs)**2*pw) + &
                  n_c*gesq(k,iobs)/(gest(k,iobs)*pw)
        t_coef = -n_a*gesp(k,iobs)/gest(k,iobs)**2 -  &
                  n_b*two*gesq(k,iobs)*gesp(k,iobs)/(gest(k,iobs)**3*pw) - &
                  n_c*gesq(k,iobs)*gesp(k,iobs)/(gest(k,iobs)**2*pw)
 
        dhdp=zero; dhdt=zero 
        if(k > 1) then
           do j = 2, k
              dhdt(j-1)= rd_over_g*(log(gesp(j-1,iobs))-log(gesp(j,iobs)))
              dhdp(j)  = dhdp(j)-rd_over_g*(gest(j-1,iobs)/gesp(j,iobs))
              dhdp(j-1)= dhdp(j-1)+rd_over_g*(gest(j-1,iobs)/gesp(j-1,iobs))
           end do
        end if
        if(k == 1)then
           dndt(k,k)=dndt(k,k)+t_coef
           dndq(k,k)=dndq(k,k)+q_coef
           dndp(k,k)=dndp(k,k)+p_coef
        else
           dndt(k,k)=dndt(k,k)+half*t_coef
           dndt(k,k-1)=dndt(k,k-1)+half*t_coef
           dndq(k,k)=dndq(k,k)+half*q_coef
           dndq(k,k-1)=dndq(k,k-1)+half*q_coef
           dndp(k,k)=p_coef
        end if
        do j = 1, nlev
           dxidt(k,j)=r1em6*radius(k)*dndt(k,j) + refindex(k)*dzdh(k)*dhdt(j)
           dxidq(k,j)=r1em6*radius(k)*dndq(k,j)
           dxidp(k,j)=r1em6*radius(k)*dndp(k,j) + refindex(k)*dzdh(k)*dhdp(j)
        end do
      end do !nlev loop
 
      d_refxrad=refxrad(nlev)-refxrad(nlev-1)
    
      do k = 1, nlevadd
         refxrad(nlev+k) = refxrad(nlev) + k*d_refxrad
         ref(nlev+k)     = ref(nlev+k-1)**2/ref(nlev+k-2) 
      end do
 
      refxrad(0)=refxrad(3)
      refxrad(nlevext+1)=refxrad(nlevext-2)
 
      do klev = 1, nlev
         refxrad_tl      = zero
         refxrad_tl(klev)= one
         ref_tl          = zero
         ref_tl(klev)    = one
         lagconst        = zero
         lagconst_tl     = zero
         d_refxrad_tl    = refxrad_tl(nlev)-refxrad_tl(nlev-1)
 
         do k = 1, nlevadd
            refxrad_tl(nlev+k) = refxrad_tl(nlev)+ k*d_refxrad_tl
            ref_tl(nlev+k)     = (two*ref(nlev+k-1)*ref_tl(nlev+k-1)/ref(nlev+k-2))-&
                                 (ref(nlev+k-1)**2/ref(nlev+k-2)**2)*ref_tl(nlev+k-2)
         end do
         refxrad_tl(0)=refxrad_tl(3)
         refxrad_tl(nlevext+1)=refxrad_tl(nlevext-2)
 
         do k=1,nlevext
            call lag_interp_const_tl(lagconst(:,k),lagconst_tl(:,k),refxrad(k-1:k+1),refxrad_tl(k-1:k+1),3)
         end do
       
         intloop2: do j = 1, ngrd
            refxrad_s(j) = sqrt(grids(j)**2 + obsimpp(iobs)**2) !x_s^2=s^2+a^2
            call get_coordinate_value(refxrad_s(j),sindx,refxrad(1:nlevext),nlevext,"increasing")
            indx=sindx
            call lag_interp_smthweights_tl(refxrad(indx-1:indx+2),refxrad_tl(indx-1:indx+2), &
                                           refxrad_s(j), lagconst(:,indx),lagconst_tl(:,indx),&
                                          lagconst(:,indx+1),lagconst_tl(:,indx+1),dw4,dw4_tl,4)
            if ( indx < nlevext) then
              if(indx==1) then
                dw4(4)=dw4(4)+dw4(1);dw4(1:3)=dw4(2:4);dw4(4)=zero
                dw4_tl(4)=dw4_tl(4)+dw4_tl(1);dw4_tl(1:3)=dw4_tl(2:4);dw4_tl(4)=zero
                indx=indx+1
              endif
              if(indx==nlevext-1) then
                dw4(1)=dw4(1)+dw4(4); dw4(2:4)=dw4(1:3);dw4(1)=zero
                dw4_tl(1)=dw4_tl(1)+dw4_tl(4); dw4_tl(2:4)=dw4_tl(1:3);dw4_tl(1)=zero
                indx=indx-1
              endif
           
              dbetaxi=(r1em6/refxrad_s(j))*dot_product(dw4_tl,ref(indx-1:indx+2))
              dbetan =(r1em6/refxrad_s(j))*dot_product(dw4,ref_tl(indx-1:indx+2))
 
              if(j == 1)then
                dbenddxi(klev)=dbetaxi
                dbenddn(klev)=dbetan
              else
                dbenddxi(klev)=dbenddxi(klev)+two*dbetaxi
                dbenddn(klev)=dbenddn(klev)+two*dbetan
              end if
           else
              cycle  obs_loop
           end if ! obs inside the new "s" grids
         end do intloop2
         dbenddxi(klev)=-dbenddxi(klev)*ds*obsimpp(iobs)
         dbenddn(klev)=-dbenddn(klev)*ds*obsimpp(iobs)
 
       end do
     
       do k = 1, nlev
          self%jac_t(k,iobs)=zero
          self%jac_q(k,iobs)=zero
          self%jac_prs(k,iobs)=zero
          do j = 1, nlev
             self%jac_t(k,iobs)  = self%jac_t(k,iobs)+dbenddxi(j)*dxidt(j,k)+ &
                                                   dbenddn(j) * dndt(j,k)
             self%jac_q(k,iobs)  = self%jac_q(k,iobs)+dbenddxi(j)*dxidq(j,k)+ &
                                                     dbenddn(j) * dndq(j,k)
             self%jac_prs(k,iobs)= self%jac_prs(k,iobs)+dbenddxi(j)*dxidp(j,k)+ &
                                                    dbenddn(j) * dndp(j,k)
          end do
        end do
        if ( nlev /= nlev1)   self%jac_prs(nlev1,iobs)=  0.
    end do obs_loop
  end do rec_loop
 
 
  deallocate(obslat)
  deallocate(obsimpp)
  deallocate(obslocr)
  deallocate(obsgeoid)
  deallocate(gest)
  deallocate(gesq)
  deallocate(gesp)
  deallocate(gesh)
  deallocate(geszs)
  deallocate(dhdp)
  deallocate(dhdt)
  deallocate(radius)
  deallocate(dzdh)
  deallocate(refindex)
  deallocate(ref) 
  deallocate(ref_tl)
  deallocate(refxrad)
  deallocate(refxrad_tl) 
  deallocate(refxrad_s)
  deallocate(dndp)
  deallocate(dndq)
  deallocate(dndt)
  deallocate(dxidp)
  deallocate(dxidt)
  deallocate(dxidq)
  deallocate(dbenddxi)
  deallocate(dbenddn)
  deallocate(lagconst)
  deallocate(lagconst_tl)
  deallocate(obsrecnum)
  if (allocated(obssrflag)) deallocate(obssrflag)
 
end if
  self%ltraj = .true.
  write(err_msg,*) myname, ": complete"
  call fckit_log%info(err_msg)
 
end subroutine ufo_gnssro_bndnbam_tlad_settraj
!----------------------------------------------------------------
subroutine ufo_gnssro_bndnbam_simobs_tl(self, geovals, hofx, obss)
  use gnssro_mod_transform
  implicit none
  class(ufo_gnssro_bndnbam_tlad), intent(in)    :: self
  type(ufo_geovals),             intent(in)    :: geovals
  real(kind_real),               intent(inout) :: hofx(:)
  type(c_ptr), value,            intent(in)    :: obss
  character(len=*), parameter     :: myname = "ufo_gnssro_bndnbam_simobs_tl"
  character(max_string)           :: err_msg 
  integer                         :: nlev, nlev1, nlocs
  integer                         :: iobs, k, irec, icount
  type(ufo_geoval), pointer       :: t_tl, prs_tl, q_tl
  real(kind_real), allocatable    :: gesT_tl(:,:), gesP_tl(:,:), gesQ_tl(:,:)
  real(kind_real)                 :: sumIntgl
 
  write(err_msg,*) myname, ": begin"
  call fckit_log%info(err_msg)
 
! check if trajectory was set
  if (.not. self%ltraj) then
      write(err_msg,*) myname, ' trajectory was not set!'
      call abor1_ftn(err_msg)
  endif
 
 
if (geovals%nlocs > 0 ) then
  hofx = missing_value(missing)
 
! check if nlocs is consistent in geovals & hofx
  if (self%nlocs /= size(hofx)) then
      write(err_msg,*) myname, ' error: nlocs inconsistent!'
      call abor1_ftn(err_msg)
  endif
  
  call ufo_geovals_get_var(geovals, var_ts,  t_tl)         ! air temperature
  call ufo_geovals_get_var(geovals, var_q,   q_tl)         ! specific humidity
  if (self%roconf%vertlayer .eq. "mass") then
    call ufo_geovals_get_var(geovals, var_prs,   prs_tl)       ! pressure
  else
    call ufo_geovals_get_var(geovals, var_prsi,  prs_tl)       ! pressure
  end if
 
  nlocs = self%nlocs
  nlev  = self%nlev
  nlev1 = self%nlev1
  allocate(gest_tl(nlev, nlocs)) 
  allocate(gesq_tl(nlev, nlocs))
  allocate(gesp_tl(nlev1,nlocs)) 
 
  if( self%iflip == 1 ) then
    do k = 1, nlev
       gest_tl(k,:) = t_tl%vals(nlev-k+1,:)
       gesq_tl(k,:) = q_tl%vals(nlev-k+1,:)
    enddo
    do k = 1, nlev1
       gesp_tl(k,:) = prs_tl%vals(nlev1-k+1,:)
    enddo
  else  ! not flipping
    do k =1, nlev
       gest_tl(k,:) = t_tl%vals(k,:)
       gesq_tl(k,:) = q_tl%vals(k,:)
    enddo
    do k =1, nlev1
       gesp_tl(k,:) = prs_tl%vals(k,:)
    enddo
  end if
 
! if all fields t and q are on mass layers,
! while p and z are on interface layers -- NBAM manner
  if ( nlev1 /= nlev ) then
     do k = nlev, 2, -1
        gest_tl(k,:) = half* (gest_tl(k,:) + gest_tl(k-1,:))
        gesq_tl(k,:) = half* (gesq_tl(k,:) + gesq_tl(k-1,:))
     enddo
  end if
 
  iobs = 0
  rec_loop: do irec = 1, self%nrecs
     obs_loop: do icount = self%nlocs_begin(irec), self%nlocs_end(irec)
        iobs = iobs + 1
        if (self%jac_t(1,iobs) /= missing ) then
        sumintgl = 0.0
        do k = 1, nlev
            sumintgl = sumintgl + self%jac_t(k,iobs)*gest_tl(k,iobs) &
                     + self%jac_q(k,iobs)*gesq_tl(k,iobs) &
                     + self%jac_prs(k,iobs)*gesp_tl(k,iobs)
 
        end do
        hofx(iobs)=sumintgl
        end if
     end do obs_loop
  end do rec_loop
 
  deallocate(gest_tl)
  deallocate(gesp_tl)
  deallocate(gesq_tl)
 
end if
 
  write(err_msg,*) "TRACE: ufo_gnssro_bndnbam_simobs_tl: begin"
  call fckit_log%info(err_msg)
 
end subroutine ufo_gnssro_bndnbam_simobs_tl
 
!----------------------------------------------------------------
subroutine ufo_gnssro_bndnbam_simobs_ad(self, geovals, hofx, obss)
  implicit none
  class(ufo_gnssro_bndnbam_tlad), intent(in)    :: self
  type(ufo_geovals),             intent(inout) :: geovals
  real(kind_real),               intent(in)    :: hofx(:)
  type(c_ptr), value,            intent(in)    :: obss
  type(ufo_geoval), pointer       :: t_ad, q_ad, prs_ad
  real(kind_real), allocatable    :: gesT_ad(:,:), gesP_ad(:,:), gesQ_ad(:,:)
  character(len=*), parameter   :: myname = "ufo_gnssro_bndnbam_simobs_ad"
  character(max_string)         :: err_msg
  integer                       :: nlocs, iobs, k, nlev,nlev1, icount, irec
 
  write(err_msg,*) myname,": begin"
  call fckit_log%info(err_msg)
   
! check if trajectory was set
  if (.not. self%ltraj) then
     write(err_msg,*) myname, ' trajectory wasnt set!'
     call abor1_ftn(err_msg)
  endif
 
if (self%nlocs > 0 ) then
 
! check if nlocs is consistent in geovals & hofx
  if (geovals%nlocs /= size(hofx)) then
      write(err_msg,*) myname, ' error: nlocs inconsistent!'
      call abor1_ftn(err_msg)
  endif
 
  call ufo_geovals_get_var(geovals, var_ts,  t_ad)         ! air temperature
  call ufo_geovals_get_var(geovals, var_q,   q_ad)         ! specific humidity
 
  if (self%roconf%vertlayer .eq. "mass") then
    call ufo_geovals_get_var(geovals, var_prs,   prs_ad)       ! pressure
  else
    call ufo_geovals_get_var(geovals, var_prsi,  prs_ad)       ! pressure
  end if
 
  nlocs = self%nlocs
  nlev  = self%nlev
  nlev1 = self%nlev1
  allocate(gest_ad(nlev,nlocs))
  allocate(gesq_ad(nlev,nlocs))
  allocate(gesp_ad(nlev1,nlocs))
 
  gest_ad = 0.0_kind_real
  gesq_ad = 0.0_kind_real
  gesp_ad = 0.0_kind_real
 
  iobs = 0
  rec_loop: do irec = 1, self%nrecs
    obs_loop: do icount = self%nlocs_begin(irec), self%nlocs_end(irec)
      iobs = iobs + 1
      if (self%jac_t(1,iobs) /= missing .and. hofx(iobs) /= missing) then
 
          do k = 1,nlev1
             if (k == nlev + 1) then
                gesp_ad(k,iobs) = gesp_ad(k,iobs) + hofx(iobs)*self%jac_prs(k,iobs)
             else
                gest_ad(k,iobs) = gest_ad(k,iobs) + hofx(iobs)*self%jac_t(k,iobs)
                gesq_ad(k,iobs) = gesq_ad(k,iobs) + hofx(iobs)*self%jac_q(k,iobs)
                gesp_ad(k,iobs) = gesp_ad(k,iobs) + hofx(iobs)*self%jac_prs(k,iobs)
             end if
          end do 
 
          if ( nlev1 /= nlev ) then
            do k = 2, nlev
               gest_ad(k-1,iobs) = half*gest_ad(k,iobs) +  gest_ad(k-1,iobs)
               gest_ad(k,iobs)   = half*gest_ad(k,iobs)
               gesq_ad(k-1,iobs) = half*gesq_ad(k,iobs) +  gesq_ad(k-1,iobs)
               gesq_ad(k,iobs)   = half*gesq_ad(k,iobs)
            enddo
          end if
 
      end if
    end do  obs_loop
  end do   rec_loop
 
  if( self%iflip == 1 ) then
    do k = 1, nlev
       t_ad%vals(nlev-k+1,:) = gest_ad(k,:)
       q_ad%vals(nlev-k+1,:) = gesq_ad(k,:)
    enddo
    do k = 1, nlev1
       prs_ad%vals(nlev1-k+1,:) = gesp_ad(k,:)
    enddo
  else  ! not flipping
    do k =1, nlev
       t_ad%vals(k,:) = gest_ad(k,:)
       q_ad%vals(k,:) = gesq_ad(k,:)
    enddo
    do k =1, nlev1
       prs_ad%vals(k,:) = gesp_ad(k,:) 
    enddo
  end if
 
  deallocate(gest_ad)
  deallocate(gesp_ad)
  deallocate(gesq_ad)
end if
  write(err_msg,*) "TRACE: ufo_gnssro_bndnbam_simobs_ad: begin"
  call fckit_log%info(err_msg)
 
end subroutine ufo_gnssro_bndnbam_simobs_ad
! ------------------------------------------------------------------------------
 
subroutine ufo_gnssro_bndnbam_tlad_delete(self)
  implicit none
  class(ufo_gnssro_bndnbam_tlad), intent(inout) :: self
  character(len=*), parameter :: myname="ufo_gnssro_bndnbam_tlad_delete"
      
  self%nlocs = 0
  self%nrecs = 0
  self%nlev  = 0
  self%nlev1 = 0
  if (allocated(self%nlocs_begin)) deallocate(self%nlocs_begin)
  if (allocated(self%nlocs_end)) deallocate(self%nlocs_end)
  if (allocated(self%jac_t)) deallocate(self%jac_t)
  if (allocated(self%jac_prs)) deallocate(self%jac_prs)
  if (allocated(self%jac_q)) deallocate(self%jac_q)
 
  self%ltraj = .false.
 
end subroutine ufo_gnssro_bndnbam_tlad_delete
 
! ------------------------------------------------------------------------------
end module ufo_gnssro_bndnbam_tlad_mod