anchor.diff:      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, local_heights
  +         layer_press, clima_layer_press
       REAL    (KIND=r8) :: thish2omxr, Mwet, Rwet, detlnp
   
       ! -----------------------
  
positive.code:           ! Ap(1:CmEp1) & Bp(1:CmEp1) are coeff at layer boundary for pressure
            ! calculation
            DO ilevel = 0, CmETA
               lpre(CmETA-ilevel) = (Ap(ilevel+1) + local_psurf(ixtrack,itimes) * Bp(ilevel+1))*1.E2 ! Pa
               clima_lpre(CmETA-ilevel) = (Ap(ilevel+1) + clima_psurf * Bp(ilevel+1))*1.E2 ! Pa
            END DO
  
  
            ! layer_press(1:CmETA) are pressure at layer centers in Pa
            ! local_heights are pressures in hPa at layer centers (I should get rid of this variable)
            DO ilevel = 1, CmETA ! for layer center
               layer_press(ilevel) = (lpre(ilevel-1) + lpre(ilevel))/2.0 ! Pa
               local_heights(ilevel)  = layer_press(ilevel)/1.E2 ! hPa
               clima_layer_press(ilevel) = (clima_lpre(ilevel-1) + clima_lpre(ilevel))/2.0 ! Pa
               clima_local_heights(ilevel)  = clima_layer_press(ilevel)/1.E2 ! hPa
               re_tmp(ilevel) = Temperature(idx_lon,idx_lat,CmETA+1-ilevel)
               re_gas(ilevel) = Gas_profiles(idx_lon,idx_lat,CmETA+1-ilevel)
               re_h2o(ilevel) = H2O_profiles(idx_lon,idx_lat,CmETA+1-ilevel)
  
negative.code:              detlnp = log(lpre(n)) - log(lpre(n1))
               lhgt = Rwet * ltmp(n) * detlnp / gplanet ! meter
               
               rho                = DEXP(layer_press(n)) / ltmp(n) / Rstar 
               aircolumn          = rho*lhgt*Navogadro*1.0E-4 ! # air/cm^2
               
               ! -------------------------------------------------------------
  