Fix for reading GFS v16 fields

src/gridcheck_gfs.f90

@@ -72,7 +72,7 @@ subroutine gridcheck_gfs
   integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
   !HSO  end
   integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip
-  real :: sizesouth,sizenorth,xauxa,pint
+  real :: sizesouth,sizenorth,xauxa,pint,xsec18
   real :: akm_usort(nwzmax)
   real,parameter :: eps=spacing(2.0_4*360.0_4) 
 
@@ -148,6 +148,8 @@ subroutine gridcheck_gfs
   call grib_get_real4_array(igrib,'values',zsec4,iret)
   call grib_check(iret,gribFunction,gribErrorMsg)
 
+  xsec18 = real(isec1(8))
+
   else ! GRIB Edition 2
 
   !read the grib2 identifiers
@@ -167,19 +169,20 @@ subroutine gridcheck_gfs
   isec1(6)=-1
   isec1(7)=-1
   isec1(8)=-1
+  xsec18 = -1.0
   if ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U
     isec1(6)=33          ! indicatorOfParameter
     isec1(7)=100         ! indicatorOfTypeOfLevel
-    isec1(8)=valSurf/100 ! level, convert to hPa
+    xsec18=valSurf/100.0 ! level, convert to hPa
   elseif ((parCat.eq.3).and.(parNum.eq.5).and.(typSurf.eq.1)) then ! TOPO
     isec1(6)=7           ! indicatorOfParameter
     isec1(7)=1           ! indicatorOfTypeOfLevel
-    isec1(8)=0
+    xsec18=real(0)
   elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.1) &
        .and.(discipl.eq.2)) then ! LSM
     isec1(6)=81          ! indicatorOfParameter
     isec1(7)=1           ! indicatorOfTypeOfLevel
-    isec1(8)=0
+    xsec18=real(0)
   endif
 
   if (isec1(6).ne.-1) then
@@ -300,7 +303,7 @@ subroutine gridcheck_gfs
 
   if((isec1(6).eq.33).and.(isec1(7).eq.100)) then ! check for U wind
     iumax=iumax+1
-    pres(iumax)=real(isec1(8))*100.0
+    pres(iumax)=xsec18*100.0
   endif
 
 

src/readwind_gfs.f90

@@ -216,7 +216,8 @@ subroutine readwind_gfs(indj,n,uuh,vvh,wwh)
       isec1(6)=2           ! indicatorOfParameter
       isec1(7)=102         ! indicatorOfTypeOfLevel
       xsec18=real(0)
-    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then ! SP
+   elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1).and. &
+        (discipl.eq.0)) then ! SP
       isec1(6)=1           ! indicatorOfParameter
       isec1(7)=1           ! indicatorOfTypeOfLevel
       xsec18=real(0)
@@ -601,7 +602,8 @@ subroutine readwind_gfs(indj,n,uuh,vvh,wwh)
         temp=tth(i,j,k,n)
         if (temp .eq. 0.0) then 
           write (*, *) i, j, k, n
-          temp = 273.0
+!          temp = 273.0
+          stop
         endif
         plev1=akm(k)+bkm(k)*ps(i,j,1,n)
         elev=ew(temp)*help/100.0