diff --git a/src/gsi/read_anowbufr.f90 b/src/gsi/read_anowbufr.f90 index e2b744eb6a..1873d0b877 100644 --- a/src/gsi/read_anowbufr.f90 +++ b/src/gsi/read_anowbufr.f90 @@ -307,6 +307,7 @@ subroutine read_anowbufr(nread,ndata,nodata,gstime,& ndata=ndata+1 nodata=nodata+1 + if(ndata>maxobs) exit cdata_all(iconc,ndata) = conc ! pm2_5 obs cdata_all(ierror,ndata) = obserror ! pm2_5 obs error diff --git a/src/gsi/read_dbz_nc.f90 b/src/gsi/read_dbz_nc.f90 index cddbd14de4..f6ac9aa112 100644 --- a/src/gsi/read_dbz_nc.f90 +++ b/src/gsi/read_dbz_nc.f90 @@ -417,6 +417,7 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no !#################### Data thinning ################### icntpnt=icntpnt+1 + if(icntpnt>maxobs) exit if(ithin > 0)then diff --git a/src/gsi/read_gmi.f90 b/src/gsi/read_gmi.f90 index 0f45aa7e28..f59529662a 100644 --- a/src/gsi/read_gmi.f90 +++ b/src/gsi/read_gmi.f90 @@ -528,6 +528,7 @@ subroutine read_gmi(mype,val_gmi,ithin,rmesh,jsatid,gstime,& flgch = 0 iobs=iobs+1 + if(iobs>maxobs) exit end do read_loop end do read_subset 690 continue diff --git a/src/gsi/read_goesglm.f90 b/src/gsi/read_goesglm.f90 index e0124abbf2..bf8639c72d 100644 --- a/src/gsi/read_goesglm.f90 +++ b/src/gsi/read_goesglm.f90 @@ -276,6 +276,7 @@ subroutine read_goesglm(nread,ndata,nodata,infile,obstype,lunout,twindin,sis) icntpnt=icntpnt+1 ndata=ndata+1 + if(ndata>maxobs) exit nodata=nodata+1 iout=ndata isort(icntpnt)=iout diff --git a/src/gsi/read_obs.F90 b/src/gsi/read_obs.F90 index 8e451c75be..cb4a7c4b8f 100644 --- a/src/gsi/read_obs.F90 +++ b/src/gsi/read_obs.F90 @@ -436,10 +436,10 @@ subroutine read_obs_check (lexist,filename,jsatid,dtype,minuse,nread) end do nread = nread + 1 end do airploop - else if(trim(filename) == 'satwndbufr')then + else if(index(filename,'satwnd') /=0 .or. index(filename,'satwhr') /=0) then lexist = .false. loop: do while(ireadmg(lnbufr,subset,idate2) >= 0) -! 5 GOES-R AMVs (NC005030, NC005031, NC005032, NC005034 and NC005039) +! 5 GOES-R AMVs (NC005030, NC005031, NC005032, NC005034, NC005039, NC005099) ! are added as the GOES-R bufr file provide do not contain other winds. ! May not be necessary with the operational satwnd BUFR if(trim(subset) == 'NC005010' .or. trim(subset) == 'NC005011' .or.& @@ -450,6 +450,7 @@ subroutine read_obs_check (lexist,filename,jsatid,dtype,minuse,nread) trim(subset) == 'NC005030' .or. trim(subset) == 'NC005031' .or.& trim(subset) == 'NC005032' .or. trim(subset) == 'NC005034' .or.& trim(subset) == 'NC005039' .or. & + trim(subset) == 'NC005099' .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.& @@ -1503,7 +1504,7 @@ subroutine read_obs(ndata,mype) else if(obstype == 'uv' .or. obstype == 'wspd10m' .or. & obstype == 'uwnd10m' .or. obstype == 'vwnd10m') then ! Process satellite winds which seperate from prepbufr - if ( index(infile,'satwnd') /=0 ) then + if ( index(infile,'satwnd') /=0 .or. index(infile,'satwhr') /=0 ) then call read_satwnd(nread,npuse,nouse,infile,obstype,lunout,gstime,twind,sis,& prsl_full,nobs_sub1(1,i)) string='READ_SATWND' diff --git a/src/gsi/read_radar.f90 b/src/gsi/read_radar.f90 index 40c77a7ee2..9ce156e736 100644 --- a/src/gsi/read_radar.f90 +++ b/src/gsi/read_radar.f90 @@ -2911,6 +2911,7 @@ subroutine read_radar(nread,ndata,nodata,infile,lunout,obstype,twind,sis,hgtl_fu !#################### Data thinning ################### icntpnt=icntpnt+1 + if(icntpnt>maxobs) exit if(ithin > 0)then if(zflag == 0)then @@ -4031,6 +4032,7 @@ subroutine read_radar_l2rw(ndata,nodata,lunout,obstype,sis,nobs,hgtl_full) end if !#################### Data thinning ################### icntpnt=icntpnt+1 + if(icntpnt>maxobs) exit ithin=1 !number of obs to keep per grid box if(radar_no_thinning) then ithin=-1 diff --git a/src/gsi/read_radar_wind_ascii.f90 b/src/gsi/read_radar_wind_ascii.f90 index 2e1b06a50c..604d9d0eca 100644 --- a/src/gsi/read_radar_wind_ascii.f90 +++ b/src/gsi/read_radar_wind_ascii.f90 @@ -509,6 +509,7 @@ subroutine read_radar_wind_ascii(nread,ndata,nodata,infile,lunout,obstype,sis,hg !#################### Data thinning ################### icntpnt=icntpnt+1 + if(icntpnt>maxobs) exit if(ithin > 0)then if(zflag == 0)then diff --git a/src/gsi/read_satwnd.f90 b/src/gsi/read_satwnd.f90 index 1679708787..943cf4d47b 100644 --- a/src/gsi/read_satwnd.f90 +++ b/src/gsi/read_satwnd.f90 @@ -18,6 +18,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis ! 253: EUMETSAT IR winds, 254: EUMETSAT WV deep layer winds ! 257,258,259: MODIS IR,WV cloud top, WV deep layer winds ! 260: VIIR IR winds +! 241: CIMSS enhanced AMV winds ! respectively ! For satellite subtype: 50-80 from EUMETSAT geostationary satellites(METEOSAT) ! 100-199 from JMA geostationary satellites(MTSAT) @@ -77,6 +78,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis ! 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 +! 2022-12-10 Bi - added code for CIMSS enhanced AMVs in new BUFR ! ! ! input argument list: @@ -155,7 +157,6 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis real(r_kind),parameter:: r799=799.0_r_kind real(r_kind),parameter:: r1200= 1200.0_r_kind real(r_kind),parameter:: r10000= 10000.0_r_kind - real(r_double),parameter:: rmiss=10d7 ! Declare local variables @@ -212,7 +213,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis real(r_double),dimension(13):: hdrdat real(r_double),dimension(4):: obsdat - real(r_double),dimension(2) :: hdrdat_test + real(r_double),dimension(2) :: hdrdat_test,hdrdat_005099 real(r_double),dimension(3,5) :: heightdat real(r_double),dimension(6,4) :: derdwdat real(r_double),dimension(3,12) :: qcdat @@ -509,7 +510,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis !GOES-R section of the 'if' statement over 'subsets' else if(trim(subset) == 'NC005030' .or. trim(subset) == 'NC005031' .or. trim(subset) == 'NC005032' .or. & - trim(subset) == 'NC005034' .or. trim(subset) == 'NC005039') then + trim(subset) == 'NC005034' .or. trim(subset) == 'NC005039' .or. trim(subset) == 'NC005099') then ! Commented out, because we need clarification for SWCM/hdrdat(9) from Yi Song ! NOTE: Once it is confirmed that SWCM values are sensible, apply this logic and replace lines 685-702 ! if(hdrdat(9) == one) then @@ -537,6 +538,8 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis itype=246 else if(trim(subset) == 'NC005031') then ! WV clear sky/deep layer itype=247 + else if(trim(subset) == 'NC005099') then + itype=241 endif else ! wind is not recognised and itype is not assigned cycle loop_report @@ -735,7 +738,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis do_qc = subset(1:7)=='NC00503'.and.nint(hdrdat(1))>=270 do_qc = do_qc.or.subset(1:7)=='NC00501' do_qc = do_qc.or.subset=='NC005081'.or.subset=='NC005091' - do_qc = do_qc.or.qcret>0 + do_qc = do_qc.or.qcret>0 ! assign types and get quality info: start @@ -1051,7 +1054,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis endif ! get quality information THIS SECTION NEEDS TO BE TESTED!!! call ufbint(lunin,rep_array,1,1,iret, '{AMVIVR}') - irep_array = int(rep_array) + irep_array = max(1,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 @@ -1175,9 +1178,12 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis endif ! Extra block for VIIRS NOAA20: End ! Extra block for GOES-R winds: Start - else if(trim(subset) == 'NC005030' .or. trim(subset) == 'NC005031' .or. trim(subset) == 'NC005032' .or. & !IR(LW) / CS WV / VIS GOES-R like winds - trim(subset) == 'NC005034' .or. trim(subset) == 'NC005039' ) then !CT WV / IR(SW) GOES-R like winds + else if(trim(subset) == 'NC005030' .or. trim(subset) == 'NC005031' .or. trim(subset) == 'NC005032' .or. & !IR(LW) / CS WV / VIS GOES-R like winds + trim(subset) == 'NC005034' .or. trim(subset) == 'NC005039' .or. trim(subset) == 'NC005099') then !CT WV / IR(SW) GOES-R like winds + if ( trim(subset) == 'NC005099' ) then + hdrdat(10)=61.23 ! set zenith angle for CIMSS AMVs to 67 to pass QC, no value in origional data + end if if(hdrdat(1) >=r250 .and. hdrdat(1) <=r299 ) then ! the range of NESDIS satellite IDs ! The sample newBUFR has SAID=259 (GOES-15) ! When GOES-R SAID is assigned, pls check @@ -1209,6 +1215,10 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis c_station_id='WV'//stationid c_sprvstg='WV' !write(6,*)'itype= ',itype + else if(trim(subset) == 'NC005099') then ! WV clear sky/deep layer + itype=241 + c_station_id='IR'//stationid + c_sprvstg='IR' endif ! call ufbint(lunin,rep_array,1,1,iret, '{AMVAHA}') @@ -1223,6 +1233,8 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis ! call ufbrep(lunin,amviii,12,irep_array,iret, 'LTDS SCLF SAID SIID CHNM SCCF ORBN SAZA BEARAZ EHAM PRLC TMDBST') ! deallocate( amviii ) + if (itype /= 241) then + call ufbint(lunin,rep_array,1,1,iret, '{AMVIVR}') irep_array = int(rep_array) allocate( amvivr(2,irep_array)) @@ -1253,7 +1265,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis if(wrf_nmm_regional) then ! type 251 has been determine not suitable to be subjected to pct1 range check - if(itype==240 .or. itype==245 .or. itype==246) then + if(itype==240 .or. itype==245 .or. itype==246 .or. itype==241) then if (pct1 < 0.04_r_kind) qm=15 if (pct1 > 0.50_r_kind) qm=15 elseif (itype==251) then @@ -1279,6 +1291,16 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis if (isflg == 1 .and. ppb > 850.0_r_kind) qm=15 ! low over land endif + else ! Assign values for the mnemonics/variables missing in original datafile for type 241 + + call ufbint(lunin,hdrdat_005099,2,1,iret, 'GNAPS PCCF'); + qifn=hdrdat_005099(2); + qm=2.0 ! do not reject the wind + pct1=0.4 ! do not reject the wind + ee=1.0 ! do not reject the wind + + endif + ! winds rejected by qc dont get used if (qm == 15) usage=r100 if (qm == 3 .or. qm ==7) woe=woe*r1_2 @@ -1288,9 +1310,12 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis if(itype==246 ) then; c_prvstg='GOESR' ; c_sprvstg='WVCT' ; endif if(itype==247 ) then; c_prvstg='GOESR' ; c_sprvstg='WVCS' ; endif if(itype==251 ) then; c_prvstg='GOESR' ; c_sprvstg='VIS' ; endif + if(itype==241 ) then; c_prvstg='GOESR' ; c_sprvstg='IR' ; endif !to be revisited I.Genkova + endif ! Extra block for GOES-R winds: End else ! wind is not recognised and itype is not assigned + write(6,*) 'read_satwnd: WIND IS NOT RECOGNIZEd and we are in hell' cycle loop_readsb endif @@ -1338,7 +1363,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis ! 3 snow ! 4 mixed if( .not. twodvar_regional) then - if(itype ==245 .or. itype ==252 .or. itype ==253 .or. itype ==240) then + if(itype ==245 .or. itype ==252 .or. itype ==253 .or. itype ==240 .or. itype ==241) then if(hdrdat(2) >20.0_r_kind) then call deter_sfc_type(dlat_earth,dlon_earth,t4dv,isflg,tsavg) if(isflg /= 0) cycle loop_readsb @@ -1465,7 +1490,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis ! GOES-R wind are identified/recognised here by subset, but it could be done by itype or SAID ! After completing the evaluation of GOES-R winds, REVISE this section!!! if(trim(subset) == 'NC005030' .or. trim(subset) == 'NC005031' .or. trim(subset) == 'NC005032' .or. & - trim(subset) == 'NC005034' .or. trim(subset) == 'NC005039' ) then + trim(subset) == 'NC005034' .or. trim(subset) == 'NC005039' .or. trim(subset) == 'NC005099') then obserr=obserr/two endif diff --git a/src/gsi/setupuwnd10m.f90 b/src/gsi/setupuwnd10m.f90 index 4552a7e81a..24a4e3d4f7 100644 --- a/src/gsi/setupuwnd10m.f90 +++ b/src/gsi/setupuwnd10m.f90 @@ -428,7 +428,7 @@ subroutine setupuwnd10m(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_d lowlevelsat=itype==242.or.itype==243.or.itype==245.or.itype==246.or. & itype==247.or.itype==250.or.itype==251.or.itype==252.or. & itype==253.or.itype==254.or.itype==257.or.itype==258.or. & - itype==259 + itype==259.or.itype==241 if (lowlevelsat .and. twodvar_regional) then call windfactor(presw,factw) data(iuob,i)=factw*data(iuob,i) diff --git a/src/gsi/setupvwnd10m.f90 b/src/gsi/setupvwnd10m.f90 index 0c601e716b..0f5b46900a 100644 --- a/src/gsi/setupvwnd10m.f90 +++ b/src/gsi/setupvwnd10m.f90 @@ -428,7 +428,7 @@ subroutine setupvwnd10m(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_d lowlevelsat=itype==242.or.itype==243.or.itype==245.or.itype==246.or. & itype==247.or.itype==250.or.itype==251.or.itype==252.or. & itype==253.or.itype==254.or.itype==257.or.itype==258.or. & - itype==259 + itype==259.or.itype==241 if (lowlevelsat .and. twodvar_regional) then call windfactor(presw,factw) data(iuob,i)=factw*data(iuob,i) diff --git a/src/gsi/setupw.f90 b/src/gsi/setupw.f90 index 087c3c34ab..784df1dfbe 100644 --- a/src/gsi/setupw.f90 +++ b/src/gsi/setupw.f90 @@ -877,7 +877,7 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav lowlevelsat=itype==242.or.itype==243.or.itype==245.or.itype==246.or. & itype==247.or.itype==250.or.itype==251.or.itype==252.or. & itype==253.or.itype==254.or.itype==257.or.itype==258.or. & - itype==259 + itype==259.or.itype==241 if (lowlevelsat .and. twodvar_regional) then call windfactor(presw,factw) data(iuob,i)=factw*data(iuob,i) @@ -1146,7 +1146,7 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav if(itype ==244) then ! AVHRR, use same as MODIS qcgross=r0_7*cgross(ikx) endif - if( itype == 245 .or. itype ==246) then + if( itype == 245 .or. itype ==246 .or. itype ==241) then if(presw <400.0_r_kind .and. presw >300.0_r_kind ) qcgross=r0_7*cgross(ikx) endif if(itype == 253 .or. itype ==254) then diff --git a/src/gsi/setupwspd10m.f90 b/src/gsi/setupwspd10m.f90 index 22618fbf9e..ad50c5b0c1 100644 --- a/src/gsi/setupwspd10m.f90 +++ b/src/gsi/setupwspd10m.f90 @@ -635,7 +635,7 @@ subroutine setupwspd10m(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_d lowlevelsat=itype==242.or.itype==243.or.itype==245.or.itype==246.or. & itype==247.or.itype==250.or.itype==251.or.itype==252.or. & itype==253.or.itype==254.or.itype==257.or.itype==258.or. & - itype==259 + itype==259.or.itype==241 if (lowlevelsat .and. twodvar_regional) then call windfactor(presw,factw) data(iuob,i)=factw*data(iuob,i)