Feature/hafsv1 maxobs goesr meso amvs (#17)

* Drop off observations when the number of obs exceeds maxobs to avoid the out of bound/dimension issue in read_anowbufr.f90 read_dbz_nc.f90 read_gmi.f90 read_goesglm.f90 read_radar.f90 read_radar_wind_ascii.f90. (From @yonghuiweng) * Add the capability of assimilating the CIMSS enhanced GOES-R AMVs in a new satwhr bufr file. (From @lbi2018 and @yonghuiweng)
NOAA-EMC · Sep 15, 2023 · 2cabebd · 2cabebd
1 parent 008c63c
commit 2cabebd
Show file tree

Hide file tree

Showing 12 changed files with 54 additions and 18 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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,7 @@ 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 +214,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 +511,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 +539,9 @@ 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
+              !write(6,*) 'NC005099 readin'
            endif
         else ! wind is not recognised and itype is not assigned
            cycle loop_report
@@ -735,7 +740,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 +1056,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 +1180,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 +1217,11 @@ 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'
+                    !write(6,*)'itype= ',itype
                  endif
 
 !                 call ufbint(lunin,rep_array,1,1,iret, '{AMVAHA}')
@@ -1223,6 +1236,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 +1268,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 +1294,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 +1313,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 +1366,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 +1493,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
@@ -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
@@ -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
@@ -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
@@ -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)