Skip to content

Commit

Permalink
Error in processing VAD winds. (NOAA-EMC#617)
Browse files Browse the repository at this point in the history
  • Loading branch information
@jderber-NOAA
jderber-NOAA authored Sep 1, 2023
1 parent 9e5aa09 commit be4a3d9
Show file tree
Hide file tree
Showing 4 changed files with 91 additions and 91 deletions.
2 changes: 2 additions & 0 deletions src/gsi/constants.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,6 +74,7 @@ module constants
public :: psv_a, psv_b, psv_c, psv_d
public :: ef_alpha, ef_beta, ef_gamma
public :: max_varname_length
public :: max_filename_length
public :: z_w_max,tfrozen
public :: qmin,qcmin,tgmin
public :: i_missing, r_missing
Expand All @@ -91,6 +92,7 @@ module constants
! Declare derived constants
integer(i_kind):: huge_i_kind
integer(i_kind), parameter :: max_varname_length=20
integer(i_kind), parameter :: max_filename_length=80
real(r_single):: tiny_single, huge_single
real(r_kind):: xai, xa, xbi, xb, dldt, rozcon,ozcon,fv, tpwcon,eps, rd_over_g
real(r_kind):: el2orc, g_over_rd, rd_over_cp, cpr, omeps, epsm1, factor2
Expand Down
16 changes: 8 additions & 8 deletions src/gsi/gsi_rfv3io_mod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ module gsi_rfv3io_mod

use kinds, only: r_kind,i_kind
use gridmod, only: nlon_regional,nlat_regional
use constants, only:max_varname_length
use constants, only:max_varname_length,max_filename_length
use gsi_bundlemod, only : gsi_bundle
use general_sub2grid_mod, only: sub2grid_info
use gridmod, only: fv3_io_layout_y
Expand Down Expand Up @@ -2207,7 +2207,7 @@ subroutine gsi_fv3ncdf_read(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin)
real(r_kind),dimension(1,grd_ionouv%nlat,grd_ionouv%nlon,grd_ionouv%kbegin_loc:grd_ionouv%kend_alloc):: hwork
character(len=max_varname_length) :: varname,vgsiname
character(len=max_varname_length) :: name
character(len=max_varname_length) :: filenamein2
character(len=max_filename_length) :: filenamein2
real(r_kind),allocatable,dimension(:,:):: uu2d_tmp
integer(i_kind) :: countloc_tmp(3),startloc_tmp(3)

Expand Down Expand Up @@ -2381,7 +2381,7 @@ subroutine gsi_fv3ncdf_read_v1(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin)
type(gsi_bundle),intent(inout) :: cstate_nouv
real(r_kind),allocatable,dimension(:,:):: uu2d
real(r_kind),dimension(1,grd_ionouv%nlat,grd_ionouv%nlon,grd_ionouv%kbegin_loc:grd_ionouv%kend_alloc):: hwork
character(len=max_varname_length) :: filenamein2
character(len=max_filename_length) :: filenamein2
character(len=max_varname_length) :: varname,vgsiname


Expand Down Expand Up @@ -2482,7 +2482,7 @@ subroutine gsi_fv3ncdf_readuv(grd_uv,ges_u,ges_v,fv3filenamegin)
character(:), allocatable:: filenamein
real(r_kind),allocatable,dimension(:,:):: u2d,v2d
real(r_kind),allocatable,dimension(:,:):: uc2d,vc2d
character(len=max_varname_length) :: filenamein2
character(len=max_filename_length) :: filenamein2
character(len=max_varname_length) :: varname,vgsiname
real(r_kind),allocatable,dimension(:,:,:,:):: worksub
integer(i_kind) u_grd_VarId,v_grd_VarId
Expand Down Expand Up @@ -2658,7 +2658,7 @@ subroutine gsi_fv3ncdf_readuv_v1(grd_uv,ges_u,ges_v,fv3filenamegin)
real(r_kind),allocatable,dimension(:,:):: us2d,vw2d
real(r_kind),allocatable,dimension(:,:):: uorv2d
real(r_kind),allocatable,dimension(:,:,:,:):: worksub
character(len=max_varname_length) :: filenamein2
character(len=max_filename_length) :: filenamein2
character(len=max_varname_length) :: varname
integer(i_kind) nlatcase,nloncase
integer(i_kind) kbgn,kend
Expand Down Expand Up @@ -2778,7 +2778,7 @@ subroutine gsi_fv3ncdf_read_ens_parallel_over_ens(filenamein,fv3filenamegin, &
real(r_kind),dimension(nlat,nlon,nsig),intent(out),optional:: delp,tsen,w,q,oz,ql,qr,qs,qi,qg,dbz
character(len=max_varname_length) :: varname
character(len=max_varname_length) :: name
character(len=max_varname_length), allocatable,dimension(:) :: varname_files
character(len=max_filename_length), allocatable,dimension(:) :: varname_files

integer(i_kind) nlatcase,nloncase,nxcase,nycase,countloc(3),startloc(3),countloc_tmp(3),startloc_tmp(3)
integer(i_kind) ilev,ilevtot,inative,ivar
Expand Down Expand Up @@ -4097,7 +4097,7 @@ subroutine gsi_fv3ncdf_write(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3file
character(len=:), allocatable, intent(in) :: filenamein
type (type_fv3regfilenameg),intent (in) :: fv3filenamegin
real(r_kind),dimension(1,grd_ionouv%nlat,grd_ionouv%nlon,grd_ionouv%kbegin_loc:grd_ionouv%kend_alloc):: hwork
character(len=max_varname_length) :: filenamein2
character(len=max_filename_length) :: filenamein2
character(len=max_varname_length) :: varname,vgsiname,name

integer(i_kind) nlatcase,nloncase,nxcase,nycase,countloc(3),startloc(3)
Expand Down Expand Up @@ -4321,7 +4321,7 @@ subroutine gsi_fv3ncdf_write_v1(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3f
character(*),intent(in):: filenamein
type (type_fv3regfilenameg),intent (in) :: fv3filenamegin
real(r_kind),dimension(1,grd_ionouv%nlat,grd_ionouv%nlon,grd_ionouv%kbegin_loc:grd_ionouv%kend_alloc):: hwork
character(len=max_varname_length) :: filenamein2
character(len=max_filename_length) :: filenamein2

integer(i_kind) kbgn,kend
integer(i_kind) inative,ilev,ilevtot
Expand Down
158 changes: 81 additions & 77 deletions src/gsi/read_prepbufr.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2003,6 +2003,85 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
! don't use MESONET psfc obs if 8th character of station id is "x")
if( kx==188 .and. psob .and. sidchr(8)=='x' ) usage=r100

! Set inflate_error logical based on qm flag
inflate_error=.false.
if (qm==3 .or. qm==7) inflate_error=.true.

if(uvob) then
selev=stnelev
oelev=obsdat(4,k)
if(kx >= 280 .and. kx < 300 )then
if (twodvar_regional.and.(kx==288.or.kx==295)) then
oelev=windsensht+selev !windsensht: read in from prepbufr
else
oelev=r10+selev
endif
if (kx == 280 )then
it29=nint(hdr(8))
if(it29 == 522 .or. it29 == 523 .or. it29 == 531)then
! oelev=r20+selev
oelev=r20
end if
end if

if (kx == 282) oelev=r20+selev
if (kx == 285 .or. kx == 289 .or. kx == 290) then
oelev=selev
selev=zero
endif
else
if((kx >= 221 .and. kx <= 229) &
.and. selev >= oelev) oelev=r10+selev
end if

! Rotate winds to rotated coordinate
uob=obsdat(5,k)
vob=obsdat(6,k)
!* thin new VAD wind and generate VAD superob
if(kx==224.and.newvad)then
klev=k+5 !*average over 6 points
! klev=k !* no average
if(klev>levs) cycle loop_readsb
diffuu=obsdat(5,k)-fcstdat(1,k)
diffvv=obsdat(6,k)-fcstdat(2,k)
if(sqrt(diffuu**2+diffvv**2)>10.0_r_kind) cycle loop_k_levs
if(abs(diffvv)>8.0_r_kind) cycle loop_k_levs
!if(abs(diffvv)>5.0.and.oelev<5000.0.and.fcstdat(3,k)>276.3) cycle loop_k_levs
if(oelev>7000.0_r_kind) cycle loop_k_levs
if(abs(diffvv)>5.0_r_kind.and.oelev<5000.0_r_kind) cycle loop_k_levs
! write(6,*)'sliu diffuu,vv::',diffuu, diffvv
uob=0.0
vob=0.0
oelev=0.0
tkk=0
do ikkk=k,klev
diffhgt=obsdat(4,ikkk)-obsdat(4,k)
if(diffhgt<301.0_r_kind)then
uob=uob+obsdat(5,ikkk)
vob=vob+obsdat(6,ikkk)
oelev=oelev+obsdat(4,ikkk)
tkk=tkk+1
end if
end do
uob=uob/tkk
vob=vob/tkk
oelev=oelev/tkk

diffuu=5.0_r_kind;diffvv=5.0_r_kind
diffhgt=0.0_r_kind
do ikkk=k,klev
diffuu=abs(obsdat(5,ikkk)-uob)
if(diffhgt<diffuu)diffhgt=diffuu
diffvv=abs(obsdat(6,ikkk)-vob)
if(diffhgt<diffvv)diffhgt=diffvv
end do

if(diffhgt>5.0_r_kind)cycle LOOP_K_LEVS !* if u-u_avg>5.0, reject
if(tkk<3) cycle LOOP_K_LEVS !* obs numb<3, reject
!* unreasonable observation, will fix this in QC package
if(sqrt(uob**2+vob**2)>60.0_r_kind)cycle LOOP_readsb
end if
end if

! Get information from surface file necessary for conventional data here

Expand Down Expand Up @@ -2088,9 +2167,6 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
! Extract pressure level and quality marks
dlnpob=log(plevs(k)) ! ln(pressure in cb)

! Set inflate_error logical based on qm flag
inflate_error=.false.
if (qm==3 .or. qm==7) inflate_error=.true.

! Temperature
if(tob) then
Expand Down Expand Up @@ -2143,6 +2219,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&

! Winds
else if(uvob) then

if (aircraftobs .and. aircraft_t_bc .and. acft_profl_file) then
call errormod_aircraft(pqm,wqm,levs,plevs,errout,k,presl,dpres,nsig,lim_qm,hdr3)
else
Expand All @@ -2151,80 +2228,6 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
woe=obserr(5,k)*errout
if (inflate_error) woe=woe*r1_2
if(obsdat(1,k) < r50)woe=woe*r1_2
selev=stnelev
oelev=obsdat(4,k)
if(kx >= 280 .and. kx < 300 )then
if (twodvar_regional.and.(kx==288.or.kx==295)) then
oelev=windsensht+selev !windsensht: read in from prepbufr
else
oelev=r10+selev
endif
if (kx == 280 )then
it29=nint(hdr(8))
if(it29 == 522 .or. it29 == 523 .or. it29 == 531)then
! oelev=r20+selev
oelev=r20
end if
end if

if (kx == 282) oelev=r20+selev
if (kx == 285 .or. kx == 289 .or. kx == 290) then
oelev=selev
selev=zero
endif
else
if((kx >= 221 .and. kx <= 229) &
.and. selev >= oelev) oelev=r10+selev
end if

! Rotate winds to rotated coordinate
uob=obsdat(5,k)
vob=obsdat(6,k)
!* thin new VAD wind and generate VAD superob
if(kx==224.and.newvad)then
klev=k+5 !*average over 6 points
! klev=k !* no average
if(klev>levs) cycle loop_readsb
diffuu=obsdat(5,k)-fcstdat(1,k)
diffvv=obsdat(6,k)-fcstdat(2,k)
if(sqrt(diffuu**2+diffvv**2)>10.0_r_kind) cycle loop_k_levs
if(abs(diffvv)>8.0_r_kind) cycle loop_k_levs
!if(abs(diffvv)>5.0.and.oelev<5000.0.and.fcstdat(3,k)>276.3) cycle loop_k_levs
if(oelev>7000.0_r_kind) cycle loop_k_levs
if(abs(diffvv)>5.0_r_kind.and.oelev<5000.0_r_kind) cycle loop_k_levs
! write(6,*)'sliu diffuu,vv::',diffuu, diffvv
uob=0.0
vob=0.0
oelev=0.0
tkk=0
do ikkk=k,klev
diffhgt=obsdat(4,ikkk)-obsdat(4,k)
if(diffhgt<301.0_r_kind)then
uob=uob+obsdat(5,ikkk)
vob=vob+obsdat(6,ikkk)
oelev=oelev+obsdat(4,ikkk)
tkk=tkk+1
end if
end do
uob=uob/tkk
vob=vob/tkk
oelev=oelev/tkk

diffuu=5.0_r_kind;diffvv=5.0_r_kind
diffhgt=0.0_r_kind
do ikkk=k,klev
diffuu=abs(obsdat(5,ikkk)-uob)
if(diffhgt<diffuu)diffhgt=diffuu
diffvv=abs(obsdat(6,ikkk)-vob)
if(diffhgt<diffvv)diffhgt=diffvv
end do

if(diffhgt>5.0_r_kind)cycle LOOP_K_LEVS !* if u-u_avg>5.0, reject
if(tkk<3) cycle LOOP_K_LEVS !* obs numb<3, reject
!* unreasonable observation, will fix this in QC package
if(sqrt(uob**2+vob**2)>60.0_r_kind)cycle LOOP_readsb
end if

if(regional .and. .not. fv3_regional)then
u0=uob
v0=vob
Expand All @@ -2237,6 +2240,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
end if
endif


cdata_all(1,iout)=woe ! wind error
cdata_all(2,iout)=dlon ! grid relative longitude
cdata_all(3,iout)=dlat ! grid relative latitude
Expand Down
6 changes: 0 additions & 6 deletions src/gsi/setupw.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -419,14 +419,8 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
if (netcdf_diag) call init_netcdf_diag_
end if

num_bad_ikx=0
do i=1,nobs
muse(i)=nint(data(iuse,i)) <= jiter
ikx=nint(data(ikxx,i))
if(ikx < 1 .or. ikx > nconvtype) then
num_bad_ikx=num_bad_ikx+1
if(num_bad_ikx<=10) write(6,*)' in setupw ',ikx,i,nconvtype,mype
end if
end do
! If HD raobs available move prepbufr version to monitor
if(nhduv > 0)then
Expand Down

0 comments on commit be4a3d9

Please sign in to comment.