anchor.diff:            ! ----------------------------------------------------
             ! Just selecting the closest location to lat lon pixel
             ! ----------------------------------------------------
  -          idx_lat = MINVAL(MINLOC(ABS(latvals(1:Cmlat) - lat(ixtrack,itimes) )))
  -          idx_lon = MINVAL(MINLOC(ABS(lonvals(1:Cmlon) - lon(ixtrack,itimes) )))
  +          idx_lat = MINVAL(MINLOC(ABS(latvals(1:Cmlat) - REAL(lat(ixtrack,itimes),KIND=8) )))
  +          idx_lon = MINVAL(MINLOC(ABS(lonvals(1:Cmlon) - REAL(lon(ixtrack,itimes),KIND=8) )))
   
             ! ----------------------------------------------
             ! Convert pixel terrain height to pressure using
             ! USSA (good below 50 km)
             ! ----------------------------------------------
  -          CALL ussa_press(REAL(terrain_height(ixtrack,itimes),KIND=r8) / 1000.0_r8, &
  -               local_psurf(ixtrack,itimes) )
  -          clima_psurf = REAL(Psurface(idx_lon, idx_lat), KIND=r8)
  +!!$          CALL ussa_press(REAL(terrain_height(ixtrack,itimes),KIND=r8) / 1000.0_r8, &
  +!!$               local_psurf(ixtrack,itimes) )
  +          
  +          ! ---------------------------------------------------------
  +          ! For OMH2O I should use terrain heigh from the climatology
  +          ! ---------------------------------------------------------
  +          local_psurf(ixtrack,itimes) = linInterpol(Cmlon, Cmlat, &
  +               lonvals(1:Cmlat), latvals(1:Cmlon), &
  +               Psurface(1:Cmlon,1:Cmlat), &
  +               REAL(lon(ixtrack,itimes),KIND=8), REAL(lat(ixtrack,itimes),KIND=8), &
  +               status=status)
  +          clima_psurf = Psurface(idx_lon, idx_lat)
             
             ! lpre(0:CmETA) is pressure at layer boundaries, same unit as
             ! local_surf, convert from hPa to Pa, it is a CmETA+1 = CmEp1 array
  
positive.code:        RETURN
      END IF
  
      ! -------------------------------
      ! Read dimension-defining arrays
      ! -------------------------------
      he5_start_1d = zerocl ; he5_stride_1d = onecl ; he5_edge_1d = Cmlatcl
      he5stat = HE5_SWrdfld ( swath_id, cli_lat_field, &
           he5_start_1d, he5_stride_1d, he5_edge_1d, latvals(1:Cmlat) )
      he5_start_1d = zerocl ; he5_stride_1d = onecl ; he5_edge_1d = Cmloncl
      he5stat = HE5_SWrdfld ( swath_id, cli_lon_field, &
           he5_start_1d, he5_stride_1d, he5_edge_1d, lonvals(1:Cmlon) )
      he5_start_1d = zerocl ; he5_stride_1d = onecl ; he5_edge_1d = CmEp1cl
      IF ( he5stat /= pge_errstat_ok ) &
           CALL error_check ( he5stat, OMI_S_SUCCESS, pge_errstat_error, OMSAO_E_PREFITCOL, &
  
negative.code:     he5stat = HE5_SWrdfld (                                &
           swath_id, TRIM(ADJUSTL(gasdatafieldname)),        &
           he5_start_3d, he5_stride_3d, he5_edge_3d,         &
           dummy_h2opr(1:Cmlon,1:Cmlat,1:CmETA) )
  
      ! -----------------------------------------
      ! Read gas datafield scale factor attribute
  