Vertintegral not working

I am using vertintegral, from version ioapi3.1, and it gives me the error that ZF and DENS are not found in the METCRO3D file:

 "METCRO3D" opened as OLD:READ-ONLY   
 File name "/pln10/STATE_CMAQ/CMAQv5.0.2/data/aqmp16/mcip/carb_met/METCRO3D_2012_07"
 File type GRDDED3 
 Execution ID "mcip"
 Grid name "METCRO_State_107"
 Dimensions: 97 rows, 107 cols, 18 lays, 15 vbles
 NetCDF ID:    131072  opened as READONLY            
 Starting date and time  2012183:060000 (6:00:00   July 1, 2012)
 Timestep                          010000 (1:00:00 hh:mm:ss)
 Maximum current record number       749
  Variable ZF not found in met file  METCRO3D
 Variable DENS not found in met file  METCRO3D
 
 *** ERROR ABORT in subroutine VERTINTEGRAL
 Fatal setup error(s)

But the variables are in the file (see below). Any clue of what could be happening?

netcdf METCRO3D_2012_07 {
dimensions:
TSTEP = UNLIMITED ; // (749 currently)
DATE-TIME = 2 ;
LAY = 18 ;
VAR = 15 ;
ROW = 97 ;
COL = 107 ;
variables:
int TFLAG(TSTEP, VAR, DATE-TIME) ;
TFLAG:units = “<YYYYDDD,HHMMSS>” ;
TFLAG:long_name = "TFLAG " ;
TFLAG:var_desc = "Timestep-valid flags: (1) YYYYDDD or (2) HHMMSS " ;
float JACOBF(TSTEP, LAY, ROW, COL) ;
JACOBF:long_name = "JACOBF " ;
JACOBF:units = "M " ;
JACOBF:var_desc = "total Jacobian at layer face " ;
float JACOBM(TSTEP, LAY, ROW, COL) ;
JACOBM:long_name = "JACOBM " ;
JACOBM:units = "M " ;
JACOBM:var_desc = "total Jacobian at layer middle " ;
float DENSA_J(TSTEP, LAY, ROW, COL) ;
DENSA_J:long_name = "DENSA_J " ;
DENSA_J:units = "KG/M2 " ;
DENSA_J:var_desc = "J-weighted air density: MM5-total dens, WRF-dry dens " ;
float WHAT_JD(TSTEP, LAY, ROW, COL) ;
WHAT_JD:long_name = "WHAT_JD " ;
WHAT_JD:units = "KG/(M*S) " ;
WHAT_JD:var_desc = "Jacobian and density weighted vert contravariant-W " ;
float TA(TSTEP, LAY, ROW, COL) ;
TA:long_name = "TA " ;
TA:units = "K " ;
TA:var_desc = "air temperature " ;
float QV(TSTEP, LAY, ROW, COL) ;
QV:long_name = "QV " ;
QV:units = "KG/KG " ;
QV:var_desc = "water vapor mixing ratio " ;
float PRES(TSTEP, LAY, ROW, COL) ;
PRES:long_name = "PRES " ;
PRES:units = "Pa " ;
PRES:var_desc = "pressure " ;
float DENS(TSTEP, LAY, ROW, COL) ;
DENS:long_name = "DENS " ;
DENS:units = "KG/M
3 " ;
DENS:var_desc = "density of air: MM5-total density, WRF-dry density " ;
float ZH(TSTEP, LAY, ROW, COL) ;
ZH:long_name = "ZH " ;
ZH:units = "M " ;
ZH:var_desc = "mid-layer height above ground " ;
float ZF(TSTEP, LAY, ROW, COL) ;
ZF:long_name = "ZF " ;
ZF:units = "M " ;
ZF:var_desc = "full-layer height above ground " ;

It looks to me that there is a mistake in the vertintegral.f routine. When checking if ZF and DENS are present in the Met file, the if statement should be IF(0 .GE. INDEX1…) instead of IF(0 .LE. INDEX1…) . Anyone had the same issue?

Here is the portion of the code that seems wrong:

!!...............   Get METCRO3D input meteorology file

IF ( .NOT. DESC3( METFILE ) ) THEN
    EFLAG = .TRUE.
    MESG = 'Could not get description for ' // METFILE
    CALL M3MESG( MESG )
 ELSE IF ( .NOT.FILCHK3( METFILE,  GRDDED3,                     &
                         NCOLS, NROWS, NLAYS, NTHIK ) ) THEN
    EFLAG = .TRUE.
    MESG = 'Inconsistent dimensions  for ' // METFILE
    CALL M3MESG( MESG )
ELSE IF ( .NOT.GRDCHK3( METFILE,                                &
                        P_ALP, P_BET, P_GAM, XCENT, YCENT,      &
                        XORIG, YORIG, XCELL, YCELL,             &
                        NLAYS, VGTYP, VGTOP, VGLEV ) ) THEN
    EFLAG = .TRUE.
    MESG = 'Inconsistent coord/grid  for ' // METFILE
    CALL M3MESG( MESG )
ELSE
    IF( **0 .LE. INDEX1**( 'ZF', NVARS3D, VNAME3D ) ) THEN
        EFLAG = .TRUE.
        MESG = 'Variable ZF not found in met file  '// METFILE
        CALL M3MESG( MESG )
    END IF
    IF( **0 .**LE**. INDEX1**( 'DENS', NVARS3D, VNAME3D ) ) THEN
        EFLAG = .TRUE.
        MESG = 'Variable DENS not found in met file  '// METFILE
        CALL M3MESG( MESG )
    END IF
END IF