Error in processing VAD winds.

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

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

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(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
 
@@ -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(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
@@ -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

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