From be4a3d91c92731d3611c46737ca695ce05cb9527 Mon Sep 17 00:00:00 2001 From: jderber-NOAA <75998838+jderber-NOAA@users.noreply.github.com> Date: Fri, 1 Sep 2023 12:34:05 -0400 Subject: [PATCH] Error in processing VAD winds. (#617) --- src/gsi/constants.f90 | 2 + src/gsi/gsi_rfv3io_mod.f90 | 16 ++-- src/gsi/read_prepbufr.f90 | 158 +++++++++++++++++++------------------ src/gsi/setupw.f90 | 6 -- 4 files changed, 91 insertions(+), 91 deletions(-) diff --git a/src/gsi/constants.f90 b/src/gsi/constants.f90 index b4cf775068..291a18ec97 100644 --- a/src/gsi/constants.f90 +++ b/src/gsi/constants.f90 @@ -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 @@ -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 diff --git a/src/gsi/gsi_rfv3io_mod.f90 b/src/gsi/gsi_rfv3io_mod.f90 index 4fcb2aba1d..05f679cb60 100644 --- a/src/gsi/gsi_rfv3io_mod.f90 +++ b/src/gsi/gsi_rfv3io_mod.f90 @@ -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 @@ -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) @@ -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 @@ -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 @@ -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 @@ -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 @@ -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) @@ -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 diff --git a/src/gsi/read_prepbufr.f90 b/src/gsi/read_prepbufr.f90 index 355441e209..2bf3a7d05d 100644 --- a/src/gsi/read_prepbufr.f90 +++ b/src/gsi/read_prepbufr.f90 @@ -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(diffhgt5.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 @@ -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 @@ -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 @@ -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(diffhgt5.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 @@ -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 diff --git a/src/gsi/setupw.f90 b/src/gsi/setupw.f90 index 62b58a0485..97ed1f8883 100644 --- a/src/gsi/setupw.f90 +++ b/src/gsi/setupw.f90 @@ -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