Skip to content

Commit

Permalink
Added changes of GSI dev branch in read_obs.F90 and read_prepbufr.f90
Browse files Browse the repository at this point in the history
  • Loading branch information
@JingCheng-NOAA
JingCheng-NOAA committed Sep 11, 2023
1 parent 6be52b0 commit 2e20205
Show file tree
Hide file tree
Showing 2 changed files with 97 additions and 3 deletions.
19 changes: 16 additions & 3 deletions src/gsi/read_obs.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -192,6 +192,7 @@ subroutine read_obs_check (lexist,filename,jsatid,dtype,minuse,nread)
if ( .not. l_use_dbz_directDA) then
if(trim(dtype) == 'dbz' )return
end if
if(trim(dtype) == 'fed' )return

! Use routine as usual

Expand Down Expand Up @@ -912,7 +913,8 @@ subroutine read_obs(ndata,mype)
obstype == 'mitm' .or. obstype=='pmsl' .or. &
obstype == 'howv' .or. obstype=='tcamt' .or. &
obstype=='lcbas' .or. obstype=='cldch' .or. obstype == 'larcglb' .or. &
obstype=='uwnd10m' .or. obstype=='vwnd10m' .or. obstype=='dbz' ) then
obstype=='uwnd10m' .or. obstype=='vwnd10m' .or. obstype=='dbz' .or. &
obstype=='fed') then
ditype(i) = 'conv'
else if (obstype == 'swcp' .or. obstype == 'lwcp') then
ditype(i) = 'wcp'
Expand Down Expand Up @@ -1083,12 +1085,13 @@ subroutine read_obs(ndata,mype)
if (ii>npem1) ii=0
if(mype==ii)then
call gsi_inquire(lenbytes,lexist,trim(dfile(i)),mype)
!call read_obs_check (lexist,trim(dfile(i)),dplat(i),dtype(i),minuse,read_rec1(i))

if (is_extOzone(dfile(i),obstype,dplat(i))) then
print*,'reading ',trim(dfile(i)),' ',obstype,' ',trim(dplat(i)),lexist,lenbytes
else
call read_obs_check(lexist,trim(dfile(i)),dplat(i),dtype(i),minuse,read_rec1(i))
call read_obs_check (lexist,trim(dfile(i)),dplat(i),dtype(i),minuse,read_rec1(i))
endif

! If no data set starting record to be 999999. Note if this is not large
! enough code should still work - just does a bit more work.

Expand Down Expand Up @@ -1302,6 +1305,10 @@ subroutine read_obs(ndata,mype)
use_hgtl_full=.true.
if(belong(i))use_hgtl_full_proc=.true.
end if
if(obstype == 'fed')then
use_hgtl_full=.true.
if(belong(i))use_hgtl_full_proc=.true.
end if
if(obstype == 'sst')then
if(belong(i))use_sfc=.true.
endif
Expand Down Expand Up @@ -1639,6 +1646,12 @@ subroutine read_obs(ndata,mype)
endif
end if

! Process flash extent density
else if (obstype == 'fed' ) then
print *, "calling read_fed"
call read_fed(nread,npuse,nouse,infile,obstype,lunout,twind,sis,nobs_sub1(1,i))
string='READ_FED'

! Process lagrangian data
else if (obstype == 'lag') then
call read_lag(nread,npuse,nouse,infile,lunout,obstype,&
Expand Down
81 changes: 81 additions & 0 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 @@ -2143,6 +2222,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 Down Expand Up @@ -2237,6 +2317,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

0 comments on commit 2e20205

Please sign in to comment.