Bring in HAFSv1 Related Changes (NOAA-EMC#608)

**Description**
Bring in HAFSv1 related maxobs changes and the capability of
assimilating GOES-R high-resolution AMVs.
Resolved NOAA-EMC#599 

**Type of change**
Use "maxobs" as a condition to check whether the number of observations
exceeds the limit, 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
This update also added the capability of assimilating the CIMSS enhanced
GOES-R AMVs in a new "satwhr" bufr file.
Please delete options that are not relevant.

- [x] Bug fix (non-breaking change which fixes an issue)
- [x] New feature (non-breaking change which adds functionality)

**How Has This Been Tested?**
This updates passed the HAFS related regression tests. All tests are
performed on Orion.
  
**Checklist**

- [x] My code follows the style guidelines of this project
- [x] I have performed a self-review of my own code
- [x] I have commented my code, particularly in hard-to-understand areas
- [x] New and existing tests pass with my changes
- [x] Any dependent changes have been merged and published

**DUE DATE for this PR is 9/21/2023.** If this PR is not merged into
develop by this date, the PR will be closed and returned to the
developer.

---------

Co-authored-by: jswhit <[email protected]>
Co-authored-by: [email protected] <[email protected]>
Co-authored-by: Jili Dong <[email protected]>
Co-authored-by: Bin Liu <[email protected]>
Co-authored-by: daprediction <[email protected]>
Co-authored-by: TingLei-NOAA <[email protected]>
Co-authored-by: TingLei-NOAA <[email protected]>
Co-authored-by: xulu <[email protected]>
Co-authored-by: Li Bi <[email protected]>
Co-authored-by: edward.safford <[email protected]>
Co-authored-by: MichaelLueken-NOAA <[email protected]>
Co-authored-by: RussTreadon-NOAA <[email protected]>
Co-authored-by: Michael Lueken <[email protected]>
Co-authored-by: AndrewEichmann-NOAA <[email protected]>
Co-authored-by: Rahul Mahajan <[email protected]>
Co-authored-by: Emily Liu <[email protected]>
Co-authored-by: emilyhcliu <[email protected]>
Co-authored-by: Andrew Collard <[email protected]>
Co-authored-by: RussTreadon-NOAA <[email protected]>

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

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
 

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

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

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'

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

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

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
 

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)

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)

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

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)