Skip to content

Commit

Permalink
GitHub Issue #243. AMVs changes:add AMVQ mitigation flag, Metop-BC an…
Browse files Browse the repository at this point in the history 
…d Him-8 new BUFR, station_id bug fix
  • Loading branch information
@ilianagenkova
ilianagenkova committed Mar 7, 2022
1 parent 0b141d1 commit 04ecee7
Show file tree
Hide file tree
Showing 3 changed files with 82 additions and 26 deletions.
1 change: 1 addition & 0 deletions src/gsi/read_obs.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -448,6 +448,7 @@ subroutine read_obs_check (lexist,filename,jsatid,dtype,minuse,nread)
trim(subset) == 'NC005039' .or. &
trim(subset) == 'NC005090' .or. trim(subset) == 'NC005091' .or.&
trim(subset) == 'NC005067' .or. trim(subset) == 'NC005068' .or. trim(subset) == 'NC005069' .or.&
trim(subset) == 'NC005047' .or. trim(subset) == 'NC005048' .or. trim(subset) == 'NC005049' .or.&
trim(subset) == 'NC005081' .or. &
trim(subset) == 'NC005072' ) then
lexist = .true.
Expand Down
94 changes: 73 additions & 21 deletions src/gsi/read_satwnd.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis
! 2013-08-26 McCarty -modified to remove automatic rejection of AVHRR winds
! 2013-09-20 Su - set satellite ID as satellite wind subtype
! 2014-07-16 Su - read VIIRS winds

! 2014-10-16 Su -add optione for 4d thinning and option to keep thinned data
! 2015-02-23 Rancic/Thomas - add thin4d to time window logical
! 2015-02-26 su - add njqc as an option to choose new non linear qc
Expand All @@ -73,7 +74,11 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis
! or hilber curve downweighting
!
! 2020-05-04 wu - no rotate_wind for fv3_regional
!
! 2021-07-25 Genkova - read GOES-17 AMVQ flag:8-mitigated height
! 16-mit.target, 24-mit.target & height; write in diag
! 2021-07-25 Genkova - added code for Metop-B/C winds in new BUFR,NC005081 !
! 2022-01-20 Genkova - added missing station_id for polar winds
! 2022-01-20 Genkova - added code for Meteosat and Himawari AMVs in new BUFR
!
!
! input argument list:
Expand Down Expand Up @@ -207,7 +212,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis
real(r_kind),dimension(nsig):: presl

real(r_double),dimension(13):: hdrdat
real(r_double),dimension(4):: obsdat
real(r_double),dimension(5):: obsdat
real(r_double),dimension(2) :: hdrdat_test
real(r_double),dimension(3,5) :: heightdat
real(r_double),dimension(6,4) :: derdwdat
Expand Down Expand Up @@ -237,8 +242,8 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis
data hdrtr_v2 /'SAID CLATH CLONH YEAR MNTH DAYS HOUR MINU SWCM SAZA OGCE SCCF SWQM'/ ! OGCE replaces GCLONG, OGCE exists in old and new BUFR
! SWQM doesn't exist in the new BUFR, so qm is initialized to '2' manually

data obstr_v1 /'HAMD PRLC WDIR WSPD'/
data obstr_v2 /'EHAM PRLC WDIR WSPD'/
data obstr_v1 /'HAMD PRLC WDIR WSPD AMVQ'/
data obstr_v2 /'EHAM PRLC WDIR WSPD AMVQ'/
! data heightr/'MDPT '/
! data derdwtr/'TWIND'/
data qcstr /' OGCE GNAP PCCF'/
Expand Down Expand Up @@ -266,7 +271,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis
! Set lower limits for observation errors
werrmin=one
nsattype=0
nreal=26
nreal=27
if(perturb_obs ) nreal=nreal+2
ntread=1
ntmatch=0
Expand Down Expand Up @@ -374,6 +379,19 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis
itype=250
endif
endif
else if(trim(subset) == 'NC005047' .or. trim(subset) == 'NC005048' .or.&
trim(subset) == 'NC005049') then ! read new Him-8 BURF
if( hdrdat(1) >=r100 .and. hdrdat(1) <=r199 ) then ! the range of JMA satellite IDS
if(hdrdat(9) == one) then ! IR winds
itype=252
else if(hdrdat(9) == two) then ! visible winds
itype=242
else if(hdrdat(9) == three) then ! WV cloud top
itype=250
else if(hdrdat(9) >= four) then ! WV deep layer, monitored
itype=250
endif
endif
else if(trim(subset) == 'NC005010' .or. trim(subset) == 'NC005011' .or. &
trim(subset) == 'NC005012' ) then
if( hdrdat(1) >=r250 .and. hdrdat(1) <=r299 ) then ! the range of NESDIS satellite IDS
Expand Down Expand Up @@ -411,8 +429,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis
endif
endif
else if( trim(subset) == 'NC005081') then
if( hdrdat(1) <10.0_r_kind .or. (hdrdat(1) >= 200.0_r_kind .and. &
hdrdat(1) <=223.0_r_kind) ) then ! the range of EUMETSAT and NOAA polar orbit satellite IDs
if( hdrdat(1) <10.0_r_kind ) then ! the range of EUMETSAT polar orbit satellite IDs new BUFR
if(hdrdat(9) == one) then ! IR winds
itype=244
else
Expand Down Expand Up @@ -627,10 +644,10 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis
call ufbint(lunin,hdrdat_test,2,1,iret, 'CLAT CLON')
if ( hdrdat_test(1) > 100000000.0_r_kind .and. hdrdat_test(2) > 100000000.0_r_kind ) then
call ufbint(lunin,hdrdat,13,1,iret,hdrtr_v2)
call ufbint(lunin,obsdat,4,1,iret,obstr_v2)
call ufbint(lunin,obsdat,5,1,iret,obstr_v2)
else
call ufbint(lunin,hdrdat,13,1,iret,hdrtr_v1)
call ufbint(lunin,obsdat,4,1,iret,obstr_v1)
call ufbint(lunin,obsdat,5,1,iret,obstr_v1)
endif

ppb=obsdat(2)
Expand Down Expand Up @@ -951,10 +968,49 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis
endif
enddo
endif
! Extra block for new EUMETSAT BUFR: Start
if(qifn <85.0_r_kind ) then ! qifn, QI without forecast
qm=15
endif
! Extra block for new JMA BUFR: Start
else if(trim(subset) == 'NC005047' .or. trim(subset) == 'NC005048' .or. &
trim(subset) == 'NC005049') then ! read new JMA BURF
if( hdrdat(1) >=r100 .and. hdrdat(1) <=r199 ) then ! The range of satellite IDs
c_prvstg='JMA'
if(hdrdat(10) >68.0_r_kind) cycle loop_readsb ! reject data zenith angle >68.0 degree
if(hdrdat(9) == one) then ! IR winds
itype=252
c_station_id='IR'//stationid
c_sprvstg='IR'
else if(hdrdat(9) == two) then ! visible winds
itype=242
c_station_id='VI'//stationid
c_sprvstg='VI'
else if(hdrdat(9) == three) then ! WV cloud top
itype=250
c_station_id='WV'//stationid
c_sprvstg='WV'
else if(hdrdat(9) >= four) then ! WV deep layer,monitoring
itype=250
qm=9 ! quality mark as 9, means the observation error needed to be set
c_station_id='WV'//stationid
c_sprvstg='WV'
endif
! get quality information THIS SECTION NEEDS TO BE TESTED!!!
call ufbint(lunin,rep_array,1,1,iret, '{AMVIVR}')
irep_array = int(rep_array)
allocate( amvivr(2,irep_array))
call ufbrep(lunin,amvivr,2,irep_array,iret, 'TCOV CVWD')
pct1 = amvivr(2,1) ! use of pct1 (a new variable in the BUFR) is introduced by Nebuda/Genkova
deallocate( amvivr )
call ufbseq(lunin,amvqic,2,4,iret, 'AMVQIC') ! AMVQIC:: GNAPS PCCF
qifn = amvqic(2,2) ! QI w/ fcst does not exist in this BUFR
ee = amvqic(2,4) ! NOTE: GOES-R's ee is in [m/s]
if(qifn <85.0_r_kind ) then ! qifn, QI without forecast
qm=15
endif
endif
! Extra block for new JMA BUFR: End
! Extra block for new EUMETSAT BUFR: Start
else if(trim(subset) == 'NC005067' .or. trim(subset) == 'NC005068' .or. &
trim(subset) == 'NC005069') then ! read new EUM BURF
if( hdrdat(1) <r80 .and. hdrdat(1) >= r50 ) then ! The range of satellite IDs
Expand Down Expand Up @@ -985,7 +1041,6 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis
call ufbrep(lunin,amvivr,2,irep_array,iret, 'TCOV CVWD')
pct1 = amvivr(2,1) ! use of pct1 (a new variable in the BUFR) is introduced by Nebuda/Genkova
deallocate( amvivr )

call ufbseq(lunin,amvqic,2,4,iret, 'AMVQIC') ! AMVQIC:: GNAPS PCCF
qifn = amvqic(2,2) ! QI w/ fcst does not exist in this BUFR
ee = amvqic(2,4) ! NOTE: GOES-R's ee is in [m/s]
Expand All @@ -995,31 +1050,27 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis
endif
! Extra block for new EUMETSAT BUFR: End
! Extra block for new Metop/AVHRR BUFR: Start
else if(trim(subset) == 'NC005081') then ! Metop/AVHRR from NESDIS
if( hdrdat(1) <10.0_r_kind .or. (hdrdat(1) >= 200.0_r_kind .and. &
hdrdat(1) <=223.0_r_kind) ) then ! The range of satellite IDs
c_prvstg='AVHRR'
else if(trim(subset) == 'NC005081') then ! Metop-B/C from NESDIS
if( hdrdat(1) <10.0_r_kind ) then ! The range of satellite IDs
c_prvstg='METOP'
if(hdrdat(9) == one) then ! IRwinds
itype=244
c_station_id='IR'//stationid
c_sprvstg='IR'
else
write(6,*) 'READ_SATWND: wrong derived method value'
endif

call ufbint(lunin,rep_array,1,1,iret, '{AMVIVR}')
irep_array = int(rep_array)
allocate( amvivr(2,irep_array))
call ufbrep(lunin,amvivr,2,irep_array,iret, 'TCOV CVWD')
pct1 = amvivr(2,1) ! use of pct1 is limited to GOES-16/17) as introduced by Nebuda/Genkova
deallocate( amvivr )

call ufbseq(lunin,amvqic,2,4,iret, 'AMVQIC') ! AMVQIC:: GNAPS PCCF
qifn = amvqic(2,2) ! QI w/ fcst does not exist in this BUFR
ee = amvqic(2,4) ! NOTE: GOES-R's ee is in [m/s]
endif
! Extra block for new Metop/AVHRR BUFR: End

! Extra block for VIIRS NOAA-20: Start
else if(trim(subset) == 'NC005091') then
if( hdrdat(1) >=r200 .and. hdrdat(1) <=r250 ) then ! Use this range in v16.*
Expand Down Expand Up @@ -1534,11 +1585,12 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis
cdata_all(22,iout)=r_prvstg(1,1) ! provider name
cdata_all(23,iout)=r_sprvstg(1,1) ! subprovider name
cdata_all(25,iout)=var_jb ! non linear qc parameter
cdata_all(26,iout)=one ! hilbert curve weight
cdata_all(26,iout)=one ! hilbert curve weight
cdata_all(27,iout)=obsdat(5) ! AMVQ for GOES-17 mitig.AMVs

if(perturb_obs)then
cdata_all(27,iout)=ran01dom()*perturb_fact ! u perturbation
cdata_all(28,iout)=ran01dom()*perturb_fact ! v perturbation
cdata_all(28,iout)=ran01dom()*perturb_fact ! u perturbation
cdata_all(29,iout)=ran01dom()*perturb_fact ! v perturbation
endif

enddo loop_readsb
Expand Down
13 changes: 8 additions & 5 deletions src/gsi/setupw.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -218,6 +218,7 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
! level; they are now loaded by
! aircraftinfo.
! 2020-05-04 wu - no rotate_wind for fv3_regional
! 2021-07-25 Genkova - write AMVQ in diagnostic files
!
! REMARKS:
! language: f90
Expand Down Expand Up @@ -287,7 +288,7 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
integer(i_kind) ihgt,ier2,iuse,ilate,ilone
integer(i_kind) izz,iprvd,isprvd
integer(i_kind) idomsfc,isfcr,iskint,iff10
integer(i_kind) ibb,ikk,ihil,idddd
integer(i_kind) ibb,ikk,ihil,idddd,iamvq

integer(i_kind) num_bad_ikx,iprev_station

Expand Down Expand Up @@ -378,8 +379,9 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
icat=24 ! index of data level category
ijb=25 ! index of non linear qc parameter
ihil=26 ! index of hilbert curve weight
iptrbu=27 ! index of u perturbation
iptrbv=28 ! index of v perturbation
iamvq=27 ! index of AMVQ
iptrbu=28 ! index of u perturbation
iptrbv=29 ! index of v perturbation

mm1=mype+1
scale=one
Expand All @@ -395,7 +397,7 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
if(conv_diagsave)then
ii=0
nchar=1
ioff0=25
ioff0=26
nreal=ioff0
if (lobsdiagsave) nreal=nreal+7*miter+2
if (twodvar_regional) then; nreal=nreal+2; allocate(cprvstg(nobs),csprvstg(nobs)); endif
Expand Down Expand Up @@ -1802,7 +1804,8 @@ subroutine contents_netcdf_diag_(udiag,vdiag)
call nc_diag_metadata("Errinv_Input", sngl(errinv_input) )
call nc_diag_metadata("Errinv_Adjust", sngl(errinv_adjst) )
call nc_diag_metadata("Errinv_Final", sngl(errinv_final) )

! AMVQ Mitigated winds
call nc_diag_metadata("Mitigation_flag_AMVQ", sngl(data(iamvq,i)) )
call nc_diag_metadata("Wind_Reduction_Factor_at_10m", sngl(factw) )

if (.not. regional .or. fv3_regional) then
Expand Down

0 comments on commit 04ecee7

Please sign in to comment.