ufo_gnssro_bndnbam_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)
!> nonlinear operator
 
module ufo_gnssro_bndnbam_mod
  use, intrinsic::  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 obsspace_mod
  use missing_values_mod
  use gnssro_mod_conf
  use gnssro_mod_constants
  use gnssro_mod_transform
  use gnssro_mod_grids,  only : get_coordinate_value
  use fckit_log_module,  only : fckit_log
  use ufo_gnssro_bndnbam_util_mod
  use ufo_utils_mod, only: cmp_strings 
 
  implicit none
  public             :: ufo_gnssro_bndnbam
  private
 
  !> Fortran derived type for gnssro trajectory
  type, extends(ufo_basis) :: ufo_gnssro_bndnbam
   type(gnssro_conf) :: roconf
  contains
   procedure :: setup     => ufo_gnssro_bndnbam_setup
   procedure :: simobs    => ufo_gnssro_bndnbam_simobs
  end type ufo_gnssro_bndnbam
 
  contains
! ------------------------------------------------------------------------------
subroutine ufo_gnssro_bndnbam_setup(self, f_conf)
 
use fckit_configuration_module, only: fckit_configuration 
implicit none
class(ufo_gnssro_bndnbam), intent(inout) :: self
type(fckit_configuration), intent(in)   :: f_conf
 
call gnssro_conf_setup(self%roconf,f_conf)
 
end subroutine ufo_gnssro_bndnbam_setup
 
subroutine ufo_gnssro_bndnbam_simobs(self, geovals, hofx, obss)
  implicit none
  class(ufo_gnssro_bndnbam), 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"
  character(max_string)                   :: err_msg
  integer                                 :: nrecs, nlocs
  integer, parameter                      :: nlevAdd = 13 !num of additional levels on top of exsiting model levels
  integer, parameter                      :: ngrd    = 80 !num of new veritcal grids for bending angle computation
  integer                                 :: iobs, k, igrd, irec, icount, kk
  integer                                 :: nlev, nlev1, nlevExt, nlevCheck
  type(ufo_geoval), pointer               :: t, q, gph, prs, zs
  real(kind_real), allocatable            :: gesT(:,:), gesZ(:,:), gesP(:,:), gesQ(:,:), gesTv(:,:), gesZs(:)
  real(kind_real), allocatable            :: obsLat(:), obsImpP(:),obsLocR(:), obsGeoid(:), obsValue(:)
  integer(c_size_t), allocatable          :: obsRecnum(:)
  real(kind_real), allocatable            :: temperature(:)
  real(kind_real), allocatable            :: humidity(:),refractivity(:),pressure(:)
  real(kind_real)                         :: temp, geop
  real(kind_real)                         :: wf
  integer                                 :: wi, wi2
  real(kind_real)                         :: grids(ngrd)
  real(kind_real), allocatable            :: refIndex(:), refXrad(:), geomz(:)
  real(kind_real), allocatable            :: ref(:), radius(:)
  real(kind_real)                         :: sIndx
  integer                                 :: indx
  integer                                 :: iflip
  integer,allocatable                     :: nlocs_begin(:)
  integer,allocatable                     :: nlocs_end(:)
  real(c_double)                          :: missing
  integer,          allocatable           :: super_refraction_flag(:), super(:), obs_max(:)
  real(kind_real),  allocatable           :: toss_max(:)
  integer                                 :: sr_hgt_idx
  real(kind_real)                         :: gradRef, obsImpH
  integer,          allocatable           :: LayerIdx(:)
 
  write(err_msg,*) myname, ": begin"
  call fckit_log%info(err_msg)
 
  nlocs   = obsspace_get_nlocs(obss) ! number of observations
  nrecs   = obsspace_get_nrecs(obss) ! number of records/profiles
  write(err_msg,*) myname, ': nlocs from gelvals and hofx, nrecs', geovals%nlocs, nlocs, nrecs
  call fckit_log%info(err_msg)
  missing = missing_value(missing)
 
  allocate(temperature(nlocs))
  temperature = missing
  if (trim(self%roconf%output_diags) .eq. "true") then
     allocate(humidity(nlocs))      ! at obs location
     allocate(pressure(nlocs))      ! at obs location
     allocate(refractivity(nlocs))  ! at obs location
     humidity = missing
     pressure = missing
     refractivity = missing
  endif
  allocate(super_refraction_flag(nlocs))
  super_refraction_flag = 0
  allocate(layeridx(nlocs))
  layeridx = 0
 
  if (nlocs > 0) then ! check if ZERO OBS
 
! check if nobs is consistent in geovals & hofx
  if (geovals%nlocs /= size(hofx)) then
    write(err_msg,*) myname, ': nlocs inconsistent!', geovals%nlocs, size(hofx)
    call abor1_ftn(err_msg)
  endif
 
! 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 if (self%roconf%vertlayer .eq. "full") then
    call ufo_geovals_get_var(geovals, var_prsi,  prs)       ! pressure
    call ufo_geovals_get_var(geovals, var_zi,    gph)       ! geopotential height
  else
    call ufo_geovals_get_var(geovals, var_prsi,  prs)       ! pressure
    call ufo_geovals_get_var(geovals, var_zi,    gph)      
    write(err_msg,*) myname,': vertlayer has to be mass of full, '//new_line('a')// &
                            '  will use full layer anyway'
    call fckit_log%info(err_msg)
  end if
 
  nlev  = t%nval   ! number of model mass levels
  nlev1 = prs%nval ! number of model pressure/height levels 
 
  allocate(gesp(nlev1,nlocs)) 
  allocate(gesz(nlev1,nlocs)) 
  allocate(gest(nlev,nlocs)) 
  allocate(gestv(nlev,nlocs))
  allocate(gesq(nlev,nlocs)) 
  allocate(geszs(nlocs))
 
! copy geovals to local background arrays
  iflip = 0
  if (prs%vals(1,1) .lt. prs%vals(prs%nval,1) ) then
     iflip = 1
     write(err_msg,'(a)')'  ufo_gnssro_bndnbam_simobs:'//new_line('a')//                         &
                         '  Model vertical height profile is in descending order,'//new_line('a')// &
                         '  but bndNBAM 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,:)
       gestv(k,:)= gest(k,:)*(1+gesq(k,:)*(rv_over_rd-1))
    enddo
    do k=1, nlev1
       gesp(k,:) = prs%vals(nlev1-k+1,:)
       gesz(k,:) = gph%vals(nlev1-k+1,:)
    enddo
  else  ! not flipping
    do k=1, nlev
       gest(k,:)  = t%vals(k,:)
       gesq(k,:)  = q%vals(k,:)
       gestv(k,:) = gest(k,:)*(1+gesq(k,:)*(rv_over_rd-1))
    enddo
 
    do k=1, nlev1
       gesp(k,:) = prs%vals(k,:)
       gesz(k,:) = gph%vals(k,:)
    enddo
  end if
       geszs(:) = zs%vals(1,:)
 
! if background t and q are on mass layers, 
!    while p and z are on interface layers, take the mean of t and q
!       -- NBAM manner
  if ( nlev1 /= nlev ) then  
     do k = nlev, 2, -1
        gesq(k,:) = half* (gesq(k,:) + gesq(k-1,:))
        gestv(k,:) = half* (gestv(k,:) + gestv(k-1,:))
!       PLEASE KEEP this COMMENT:
!       to exactly reproduce nbam, t is converted to tv, tv mean is calcualted,
!       then tv mean is converted to t mean
        gest(k,:) = gestv(k,:)/(1+ gesq(k,:)*(rv_over_rd-1))
     enddo
  end if
 
! set obs space struture
  allocate(obslat(nlocs))
  allocate(obsimpp(nlocs))
  allocate(obslocr(nlocs))
  allocate(obsgeoid(nlocs))
  allocate(obsvalue(nlocs))
  allocate(obsrecnum(nlocs))
  allocate(nlocs_begin(nrecs))
  allocate(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_db(obss, "ObsValue", "bending_angle", obsvalue)
  call obsspace_get_recnum(obss, obsrecnum)
 
  nlocs_begin=1
  nlocs_end=1
  icount = 1
  do iobs = 1, nlocs-1
    if (obsrecnum(iobs+1) /= obsrecnum(iobs)) then
      icount = icount +1  !counting number of records
      nlocs_end(icount-1)= iobs
      nlocs_begin(icount) = iobs+1
    end if
  end do
  nlocs_end(nrecs)= nlocs
  if (icount /= nrecs) then
    write(err_msg,*) myname, ': record number is not consistent :', icount, nrecs
    call abor1_ftn(err_msg)
  end if
 
  nlevext   = nlev + nlevadd
  nlevcheck = min(23, nlev)   !number of levels to check super refaction
 
! define new integration grids
  do igrd = 0, ngrd-1
     grids(igrd+1) = igrd * ds
  end do 
 
! bending angle forward model starts
  allocate(geomz(nlev))    ! geometric height
  allocate(radius(nlev))   ! tangent point radisu to earth center
  allocate(ref(nlevext))   ! refractivity
  allocate(refindex(nlev))              !refactivity index n
  allocate(refxrad(0:nlevext+1))        !x=nr, model conuterpart impact parameter
 
  allocate(super(nlocs))
  allocate(toss_max(nrecs))
  allocate(obs_max(nrecs))
 
  iobs = 0
  hofx =  missing
  super = 0
  obs_max  = 0
  toss_max = 0
 
  rec_loop: do irec = 1, nrecs
 
    obs_loop: do icount = nlocs_begin(irec), nlocs_end(irec)
 
      iobs = icount
      do k = 1, nlev
!        compute guess geometric height from geopotential height
         call geop2geometric(obslat(iobs), gesz(k,iobs)-geszs(iobs), geomz(k))
         radius(k) = geomz(k) + 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))  cycle obs_loop
 
!     save the obs vertical location index (unit: model layer)
      layeridx(iobs) = min(max(1, int(sindx)), nlev)
 
!     calculating virtual temperature at obs location to obs space for BackgroundCheck RONBAM
      indx=sindx
      wi=min(max(1,indx),nlev)
      wi2=max(1,min(indx+1,nlev))
      wf=sindx-float(wi)
      wf=max(zero,min(wf,one))
      temperature(iobs)=gestv(wi,iobs)*(one-wf)+gestv(wi2,iobs)*wf
      if (trim(self%roconf%output_diags) .eq. "true") then
         humidity(iobs)= gesq(wi,iobs)*(one-wf)+gesq(wi2,iobs)*wf 
         temp          = gest(wi,iobs)*(one-wf)+gest(wi2,iobs)*wf 
         geop          = gesz(wi,iobs)*(one-wf)+gesz(wi2,iobs)*wf
         pressure(iobs)= gesp(wi,iobs)/exp(two*grav*(geop-gesz(wi,iobs))/(rd*(temperature(iobs)+gestv(wi,iobs))))
         call compute_refractivity(temp, humidity(iobs), pressure(iobs),   &
                                 refractivity(iobs), self%roconf%use_compress)
      end if
 
!     (2) super-refaction
!     (2.1) GSI style super refraction check
      if(cmp_strings(self%roconf%super_ref_qc, "NBAM")) then
 
        obsimph = (obsimpp(iobs) - obslocr(iobs)) * r1em3 !impact heigt: a-r_earth
 
        if (obsimph <= six) then
           do k = nlevcheck, 1, -1
 
!             N gradient
              gradref = 1000.0 * (ref(k+1)-ref(k))/(radius(k+1)-radius(k))
!             check for model SR layer
              if (abs(gradref) >= 0.75*crit_gradrefr .and. obsimpp(iobs) <= refxrad(k+2)) then
                  super_refraction_flag(iobs) = 1
                  cycle obs_loop
              endif
           end do
 
!          relax to close-to-SR conditions, and check if obs is inside model SR layer
           
           if (self%roconf%sr_steps > 1                 &
              .and. obsvalue(iobs) >= 0.03 ) then
 
               do k = nlevcheck, 1, -1
                  gradref = 1000.0 * (ref(k+1)-ref(k))/(radius(k+1)-radius(k))
                  if (abs(gradref) >= half*crit_gradrefr & 
                     .and. super(iobs) == 0                   &
                     .and. toss_max(irec) <= obsvalue(iobs)) then
                      obs_max(irec) = iobs
                      toss_max(irec)= max(toss_max(irec), obsvalue(iobs))
                      super(iobs) = 1
                  end if
              end do ! k
 
           end if   ! end if(self%roconf%sr_steps > 1 
        end if ! obsImpH <= six
 
!    ROPP style super refraction check
     else if(cmp_strings(self%roconf%super_ref_qc, "ECMWF")) then
 
       sr_hgt_idx = 1
       do k = nlev, 2, -1
          if (refxrad(k) - refxrad(k-1) < 10.0) THEN
             sr_hgt_idx = k
             exit
          end if
       end do
 
       if (obsimpp(iobs) < refxrad(sr_hgt_idx)) then
          super_refraction_flag(iobs) = 1
          cycle obs_loop
       end if
 
     else
       write(err_msg,*) myname, ': super refraction method has to be NBAM or ECMWF!'
       call abor1_ftn(err_msg)
     end if
 
     if (super_refraction_flag(iobs) .eq. 0) then
     call ufo_gnssro_bndnbam_simobs_single(   &
               obslat(iobs), obsgeoid(iobs), obslocr(iobs), obsimpp(iobs), &
               grids, ngrd, &
               nlev, nlevext, nlevadd, nlevcheck, &
               radius(1:nlev),ref(1:nlevext),refindex(1:nlev),refxrad(0:nlevext),  &
               hofx(iobs))
 
     end if
    end do obs_loop
  end do rec_loop
 
  if (cmp_strings(self%roconf%super_ref_qc, "NBAM") .and. self%roconf%sr_steps > 1 ) then
     rec_loop2: do irec = 1, nrecs
 
       if (obs_max(irec) > 0 ) then
 
          obs_loop2: do k = nlocs_begin(irec), nlocs_end(irec)
            obsimph = (obsimpp(k) - obslocr(k)) * r1em3
            if (obsimph<=six .and. obsimpp(k)<=obsimpp(obs_max(irec)) .and.  &
                hofx(k)/=missing .and. super_refraction_flag(k)==0) then
                super_refraction_flag(k)=2
            end if
          end do obs_loop2
 
       end if
 
     end do rec_loop2
  end if
 
  deallocate(obslat)
  deallocate(obsimpp)
  deallocate(obslocr)
  deallocate(obsgeoid)
  deallocate(gesp) 
  deallocate(gesz) 
  deallocate(gest) 
  deallocate(gestv) 
  deallocate(gesq)
  deallocate(geszs) 
  deallocate(ref)
  deallocate(refindex)
  deallocate(refxrad)
  deallocate(geomz)
  deallocate(radius)
  deallocate(obsrecnum)
  deallocate(nlocs_begin)
  deallocate(nlocs_end)
  deallocate(super)
  deallocate(toss_max)
  deallocate(obs_max)
 
  write(err_msg,*) myname, ": complete"
  call fckit_log%info(err_msg)
  end if ! end check if ZERO OBS
 
! putting virtual temeprature at obs location to obs space for BackgroundCheck RONBAM
  call obsspace_put_db(obss, "MetaData", "virtual_temperature", temperature)
! putting super refraction flag to obs space 
  call obsspace_put_db(obss, "SRflag",   "bending_angle", super_refraction_flag)
! saving obs vertical model layer postion for later
  call obsspace_put_db(obss, "LayerIdx",   "bending_angle", layeridx)
  if (trim(self%roconf%output_diags) .eq. "true") then
      call obsspace_put_db(obss, "ObsDiag", "specific_humidity", humidity)
      call obsspace_put_db(obss, "ObsDiag", "refractivity", refractivity)
      call obsspace_put_db(obss, "ObsDiag", "pressure", pressure)
      deallocate(humidity)
      deallocate(pressure)
      deallocate(refractivity)
  end if
  deallocate(super_refraction_flag)
  deallocate(temperature)
  deallocate(layeridx)
 
end subroutine ufo_gnssro_bndnbam_simobs
! ------------------------------------------------------------------------------
end module ufo_gnssro_bndnbam_mod