kernels calculated based on 5 yrs ERA-Interim reanalysis and RRTMG

the units of the atmospgeric temperature and water vapor kernels are w/m2/k/100hpa
dR(Ta)=sum(Kt(i)*dTa(i)*dP(i)/100), i denotes each layer, dP(i) is the thickness of each layer
dR(q) =sum(Kq(i)*dlnq(i)*Ta(i)^2*Rv/Lv*dP(i)/100)

in the file dp.nc
player denotes the middle pressure of each layer we perturbed
plevel denotes the boundary pressure of each layer
dp denotes the thickness (in hPa) of each layer

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
Note that the atmospheric temperature kernel at surface has 25 layers but not 24 as at TOA, because the change in surface temperature contributes to the effect of the lowest layer on surface downward radiation. 

to use the kernels:
1. read kernels of atmospheric temperature (Kt), water vapor (Kq),surface temperature (Kts) and albedo (Kalb)
   read dp, player, plevel from dp.nc

2. read anomaly of Ta (dTa), q (dq), Ts (dTs), and albedo (dalb), and Ta, surface pressure (sp) 
   then interpolate to the resolution of the kernels

3. calculate the thickness of each layer over everywhere 

nlayer=24
dps=new((/nlayer,nlat,nlon/),"float")
dps=dps@_FillValue
dp1=dps(0,:,:) ; the thickness of the lowest layer 

do i=0,nlat-1
do j=0,nlon-1
  indyp=min(ind(sp(i,j)-player.ge.0))
  dps(indyp,i,j)=sp(i,j)-plevel(indyp+1)
  dps(indyp+1:,i,j)=dp(indyp+1:)
  dp1(i,j)=sp(i,j)-plevel(indyp+1)
end do
end do

4. calculate the feedbacks
dr_q =  dim_sum_n(Kq*dqn*dps/100.,0),0)

dr_Ta=  dim_sum_n(Kt(1:,:,:)*dTa*dps/100.,0),0)+Kt(0,:,:)*dTs*dp1/100.

dr_Ts=  Kts*dTs
dr_alb= Kalb*dalb


