anchor.diff:      ! ---------------
       ! Local variables
       ! ---------------
  -    INTEGER (KIND=i4) :: itimes, ixtrack, spix, epix, idx_lat, idx_lon, ilevel, n, n1
  -    REAL    (KIND=r8)                      :: rho, lhgt, aircolumn
  -    REAL    (KIND=r8), DIMENSION (0:CmETA) :: lpre, level_press
  -    REAL    (KIND=r8), DIMENSION (1:CmETA) :: ltmp, layer_press
  +    INTEGER (KIND=i4) :: itimes, ixtrack, spix, epix, idx_lat, idx_lon, ilevel, n, n1, one, status
  +    REAL    (KIND=r8)                      :: rho, lhgt, aircolumn, clima_psurf
  +    REAL    (KIND=r8), DIMENSION (0:CmETA) :: lpre, level_press, clima_lpre
  +    REAL    (KIND=r8), DIMENSION (1:CmETA) :: ltmp, lh2o, lgas, re_tmp, re_h2o, re_gas, &
  +         layer_press, clima_layer_press, clima_local_heights
       REAL    (KIND=r8) :: thish2omxr, Mwet, Rwet, detlnp
   
       ! -----------------------
  
positive.code:           idx_lat = MINVAL(MINLOC(ABS(latvals(1:Cmlat) - lat(ixtrack,itimes) )))
            idx_lon = MINVAL(MINLOC(ABS(lonvals(1:Cmlon) - lon(ixtrack,itimes) )))
  
            local_temperature(ixtrack,itimes,1:CmETA) = Temperature(idx_lon,idx_lat,1:CmETA)
            ltmp(1:CmETA)                             = Temperature(idx_lon,idx_lat,1:CmETA)
  
            ! ----------------------------------------------
            ! Convert pixel terrain height to pressure using
            ! Xiong suggested to use pressure altitude:
            !  Z = -16 alog10 (P / Po) Z in km and P in hPa.
            ! ----------------------------------------------
            local_psurf(ixtrack,itimes) = 1013.0_r8 * &
                 (10.0_r8 ** (REAL(terrain_height(ixtrack,itimes),KIND=r8) / 1000.0_r8 / (-16.0_r8)))
  
            ! 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
            ! Ap(1:CmEp1) & Bp(1:CmEp1) are coeff at layer boundary for pressure
            ! calculation
            DO ilevel = 0, CmETA
               lpre(ilevel) = (Ap(ilevel+1) + local_psurf(ixtrack,itimes) * Bp(ilevel+1))*1.E2 ! Pa
            END DO
  
  
  
negative.code:     ! ----------------------
      locerrstat = pge_errstat_ok
      
      one = 1_i4
      ! ------------------------------------------------------------
      ! Find the Climatology corresponding to each lat and lon pixel
      ! ------------------------------------------------------------
  