ufo_gnssro_bndgsi_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 to handle gnssro bending angle observations following 
!> the NCEP/GSI (2018 Aug) implementation

module ufo_gnssro_bndgsi_mod
  use iso_c_binding
  use kinds
  use ufo_vars_mod
  use ufo_geovals_mod
  use ufo_geovals_mod_c,   only: ufo_geovals_registry
  use vert_interp_mod
  use ufo_basis_mod,      only: ufo_basis
  use lag_interp_mod,     only: lag_interp_const, lag_interp_smthweights
  use obsspace_mod

  implicit none
  public             :: ufo_gnssro_bndgsi
  private

  !> Fortran derived type for gnssro trajectory
  type, extends(ufo_basis) :: ufo_gnssro_bndgsi
  contains
   procedure :: simobs    => ufo_gnssro_bndgsi_simobs
  end type ufo_gnssro_bndgsi

  contains
! ------------------------------------------------------------------------------
  subroutine ufo_gnssro_bndgsi_simobs(self, geovals, hofx, obss)
      use gnssro_mod_constants
      use gnssro_mod_transform
      use gnssro_mod_grids, only: get_coordinate_value
      use obsspace_mod, only: obsspace_get_db
      implicit none
      class(ufo_gnssro_bndGSI), 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_bndgsi_simobs"
      logical, parameter              :: use_compress = .true.
      real(kind_real), parameter      :: r1em6 = 1.0e-6_kind_real
      real(kind_real), parameter      :: r1em3 = 1.0e-3_kind_real
      real(kind_real), parameter      :: six   = 6.0_kind_real
      real(kind_real), parameter      :: ds    = 10000.0_kind_real
      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
      integer, parameter              :: max_string    = 800
      real(kind_real), parameter      :: miss_values   = -99999.00
      real(kind_real), parameter      :: crit_gradRefr = 157.0_kind_real !criteria for the refractivity gradient

      character(max_string)           :: err_msg
      integer                         :: iobs, ilev, igrd, klev
      integer                         :: nlev, nlev1, nobs, nlevExt, nlevCheck
      real(kind_real)                 :: rnlevExt
      real(kind_real)                 :: w4(4), dw4(4)
      integer                         :: ierr
      type(ufo_geoval), pointer       :: t, gphi, prsi, q
      real(kind_real), allocatable    :: gesT(:,:), gesZi(:,:), gesPi(:,:), gesQ(:,:)
      real(kind_real), allocatable    :: refr(:,:), radius(:,:)
      real(kind_real), allocatable    :: refrIndex(:), refrXrad(:), geomzi(:), refrXrad_new(:)
      real(kind_real), allocatable    :: lagConst(:,:)
      real(kind_real), allocatable    :: obsLat(:), obsImpP(:), obsLocR(:), obsGeoid(:)
      real(kind_real)                 :: specHmean, tmean
      real(kind_real)                 :: d_refrXrad
      real(kind_real)                 :: derivRefr_s(ngrd),grids(ngrd),refrXrad_s(ngrd)
      real(kind_real)                 :: sIndx, sIndxExt
      integer                         :: indx
      real(kind_real)                 :: bndIntgd, bendingAngle, obsImpH, gradRefr
      logical                         :: obs_check, qc_layer_SR
      integer                         :: nobs_outIntgl
      integer                         :: count_SR, top_layer_SP, top_layer_SR, bot_layer_SR !for super refraction
      integer                         :: count_rejection
       real(kind_real)                 :: jacob    
      nobs = 0
      nlev = 0
      nlev1 = 0
      nlevext = 0
      sindxext = one
      hofx(:) = miss_values

      ! check if nobs is consistent in geovals & hofx
      if (geovals%nobs /= size(hofx)) then
        write(err_msg,*) myname_, ' error: nobs inconsistent!'
        call abor1_ftn(err_msg)
      endif
      ! check if prsi (pressure at model interface levels) variable is in geovals and get it
      call ufo_geovals_get_var(geovals, var_prsi, prsi,status=ierr)
      if (ierr/=0) then
         write(err_msg,*) myname_, trim(var_prsi), ' doesnt exist'
         call abor1_ftn(err_msg)
      endif
      ! check if t (temperature) variable is in geovals and get it
      call ufo_geovals_get_var(geovals, var_t, t,status=ierr)  
      if (ierr/=0) then
         write(err_msg,*) myname_, trim(var_t), ' doesnt exist'
         call abor1_ftn(err_msg)
      endif
      ! check if q(specific humidity) variable is in geovals and get it
      call ufo_geovals_get_var(geovals, var_q, q,status=ierr)
       if (ierr/=0) then
         write(err_msg,*) myname_, trim(var_q), ' doesnt exist'
         call abor1_ftn(err_msg)
      endif
      ! check if gphi (geopotential height at model interface levels) variable is in geovals and get it
      call ufo_geovals_get_var(geovals, var_zi, gphi,status=ierr) 
      if (ierr/=0) then
         write(err_msg,*) myname_, trim(var_zi), ' doesnt exist'
         call abor1_ftn(err_msg)
      endif

 
      nlev  = t%nval ! number of model levels
      nlev1 = prsi%nval ! number of model interface levels 
      if (nlev1 /= nlev+1) then
         write(err_msg,*) myname_, ' Numbers of vertical profiles, nlev1 and nlev, dont match'
         call abor1_ftn(err_msg)
      endif
      nlevext = nlev + nlevadd
      nlevcheck = min(23, nlev) !number of levels to check super refraction
      nobs  = geovals%nobs ! number of observations

      allocate(gespi(nlev1,nobs)) 
      allocate(geszi(nlev1,nobs)) 
      allocate(gest(nlev,nobs)) 
      allocate(gesq(nlev,nobs)) 

      !FV3 background is from top to bottom. Reserse the vertical order 
      do ilev=1, nlev
         gest(ilev,:) = t%vals(nlev-ilev+1,:)
         gesq(ilev,:) = q%vals(nlev-ilev+1,:)
      enddo

      do ilev=1, nlev1
         gespi(ilev,:) = prsi%vals(nlev1-ilev+1,:)
         geszi(ilev,:) = gphi%vals(nlev1-ilev+1,:)
      enddo

      do igrd = 0, ngrd-1
          grids(igrd+1) = igrd * ds
      end do 

      allocate(geomzi(nlev))      !geometric height at interface model levels
      allocate(radius(nlev,nobs)) !distance between earth center to model interface level

      allocate(refr(nlevext,nobs))    !refractivity N 
      allocate(refrindex(nlev))       !refractivity index n
      allocate(refrxrad(0:nlevext+1)) !x=nr, r: radius
      allocate(lagconst(3,nlevext))   !x=nr, r: radius
      allocate(refrxrad_new(nlevext+newadd))

      allocate(obslat(nobs))
      allocate(obsimpp(nobs))
      allocate(obslocr(nobs))
      allocate(obsgeoid(nobs))

      call obsspace_get_db(obss, "Metadata", "Latitude", obslat)
      call obsspace_get_db(obss, "Metadata", "IMPP", obsimpp)   !observed impact parameter; meter
      call obsspace_get_db(obss, "Metadata", "ELRC", obslocr)   !local radius of earth; meter
      call obsspace_get_db(obss, "Metadata", "GEODU", obsgeoid) !Geoid; meter

      nobs_outintgl = 0 !initialize count of observations out of integral grids  
      count_rejection = 0

      obs_loop: do iobs = 1, nobs 
         do ilev = 1,nlev 

         ! compute guess geometric height from geopotential height at model interface levels
            call geop2geometric( obslat(iobs), geszi(ilev,iobs), geomzi(ilev), jacob)

         ! compute guess radius 
            radius(ilev,iobs) = geomzi(ilev) + obsgeoid(iobs) + obslocr(iobs)   ! radius r

         ! compute guess refractivity and refractivity index at model interface levels
            if(ilev > 1) then
               spechmean = (gesq(ilev,iobs) + gesq(ilev-1,iobs))/two
               tmean = (gest(ilev,iobs) + gest(ilev-1,iobs) )/two
            else
               spechmean = gesq(1,iobs)
               tmean = gest(1,iobs)
            endif
            call compute_refractivity(tmean, spechmean, gespi(ilev,iobs), refr(ilev,iobs),use_compress) !refr N
            refrindex(ilev) = one + (r1em6*refr(ilev,iobs))         ! refr index n
            refrxrad(ilev)  = refrindex(ilev) * radius(ilev,iobs)   ! x=nr

         end do 

         ! data rejection based on model background !
         ! (1) skip data below the model levels
         call get_coordinate_value(obsimpp(iobs), sindx,refrxrad(1),nlev,"increasing")
         if (sindx < one .or. sindx > float(nlev)) then 
             count_rejection = count_rejection + 1
             cycle
         end if

         ! (2) super-refraction
         qc_layer_sr=.false.
         count_sr=0
         top_layer_sr=0
         bot_layer_sr=0
         obsimph = (obsimpp(iobs) - obslocr(iobs)) * r1em3 !impact heigt: a-r_earth
         if (obsimph <= six) then
            do klev=nlevcheck,1,-1

               ! check for model SR layer 
               gradrefr = 1000.0_kind_real * (refr(klev+1,iobs)-refr(klev,iobs)) / (radius(klev+1,iobs)-radius(klev,iobs))

               if (abs(gradrefr)>= half*crit_gradrefr) then  !Super refractivity - likely, to be used in obs SR qc
                  qc_layer_sr=.true.   !SR-likely layer detected
               endif

               !if (((ref_rad(klev+1)-ref_rad(klev))/(radius(klev+1,iobs)-radius(klev,iobs))) < zero) then
               if (abs(gradrefr) >= 0.75_kind_real*crit_gradrefr) then  !relax to close-to-SR conditions
                  count_sr=count_sr+1 ! layers of SR

                  if (count_sr > 1 ) then
                     !if(abs(bot_layer_SR-klev) > 1) write(6,*) 'WARNING GPSRO: non-consecutive SR layers'
                     bot_layer_sr=klev
                  else
                     top_layer_sr=klev
                     bot_layer_sr=top_layer_sr
                  endif

               endif
            end do
            if (top_layer_sr >= 1 .and. obsimpp(iobs) <= refrxrad(top_layer_sr+2)) then !obs inside model SR layer
               count_rejection = count_rejection + 1
               cycle
            end if
          endif


         ! Extend atmosphere above interface level nlev
         d_refrxrad = refrxrad(nlev) - refrxrad(nlev-1)
         do ilev=1,nlevadd
            refrxrad(nlev+ilev)=refrxrad(nlev)+ ilev*d_refrxrad ! extended x_i
            refr(nlev+ilev,iobs)=refr(nlev+ilev-1,iobs)**2/refr(nlev+ilev-2,iobs) ! exended N_i
         end do

         refrxrad(0)=refrxrad(3)
         refrxrad(nlevext+1)=refrxrad(nlevext-2)

         do ilev = 1,nlevext
            call lag_interp_const(lagconst(:,ilev),refrxrad(ilev-1:ilev+1),3)
         enddo

         ! Set up a new equally-spaced vertical grid for integral 

         derivrefr_s = zero
         grids_loop: do igrd =1,ngrd
         !use the new grids (s) for bending angle computation
           refrxrad_s(igrd)=sqrt(grids(igrd)**2 + obsimpp(iobs)**2) !x_s^2=s^2+a^2

           call get_coordinate_value(refrxrad_s(igrd), sindx,refrxrad(1:nlevext),nlevext,"increasing")

           rnlevext = float(nlevext)

           if (sindx < rnlevext) then  !obs inside the new grid
           !HS 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(refrxrad(indx-1:indx+2),refrxrad_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
              !need make sure dot_product is available or add the code
              derivrefr_s(igrd)=dot_product(dw4,refr(indx-1:indx+2,iobs)) !derivative dN/dx_s
              derivrefr_s(igrd)=max(zero,abs(derivrefr_s(igrd)))

           else
              obs_check=.true.
              nobs_outintgl=nobs_outintgl+1
              d_refrxrad=refrxrad(nlev)-refrxrad(nlev-1)
              do klev=1,newadd
                 refrxrad_new(nlevext+klev)=refrxrad(nlevext)+ klev*d_refrxrad ! extended x_i
              end do
              do klev=1,nlevext
                 refrxrad_new(klev)=refrxrad(klev)
              enddo
              call get_coordinate_value(refrxrad_s(igrd), sindx,refrxrad_new(1:nlevext+newadd),nlevext+newadd,"increasing")
              sindxext=max(sindx,sindxext)
           endif !obs in new grid
         end do grids_loop

!        bending angle (radians)
         bendingangle = ds*derivrefr_s(1)/refrxrad_s(1)
         do igrd = 2,ngrd
            bndintgd     = ds*derivrefr_s(igrd)/refrxrad_s(igrd)
            bendingangle = bendingangle + two*bndintgd
         end do
         bendingangle=r1em6 * obsimpp(iobs) * bendingangle
         hofx(iobs) = bendingangle

      end do obs_loop

      write(6,*) 'bndGSI: hofx ', &
                 'min = ', minval(hofx, mask=hofx > miss_values), 'min index = ', minloc(hofx), &
                 'max = ', maxval(hofx, mask=hofx > miss_values), 'max index = ', maxloc(hofx)
      write(6,*) 'bndGSI: ', count_rejection, ' out of ', nobs, ' rejected due to model vertical range and super refraction'

      !for tuning the nlevExt. New grids (s) should be in range  with nlevExt. If not, adjust the hardwired 
      if (nobs_outintgl>=1) then
         write(6,*)'bndGSI: Warning',nobs_outintgl,'obs outside integration grid. Increase nlevExt to',&
         int(sindxext)
      endif

      deallocate(obslat)
      deallocate(obsimpp)
      deallocate(obslocr)
      deallocate(obsgeoid)

      deallocate(gespi) 
      deallocate(geszi) 
      deallocate(gest) 
      deallocate(gesq) 
      deallocate(refr) 
      deallocate(refrindex) 
      deallocate(refrxrad) 
      deallocate(geomzi) 
      deallocate(radius) 
      deallocate(lagconst) 
      deallocate(refrxrad_new) 

   end subroutine ufo_gnssro_bndgsi_simobs
! ------------------------------------------------------------------------------
end module ufo_gnssro_bndgsi_mod