* ========================================
* 2a. Calculate the temperature advection from gridded dataset
* ========================================
'set rgb  16    0    0  255'
'set rgb  17   55   55  255'
'set rgb  18  110  110  255'
'set rgb  19  165  165  255'
'set rgb  20  220  220  255'

'set rgb  21  255  220  220'
'set rgb  22  255  165  165'
'set rgb  23  255  110  110'
'set rgb  24  255   55   55'
'set rgb  25  255    0    0'

'reinit'
'open narr-a_221_19930314_0000_000.ctl'
'q dims'


lons=subwrd(xline,6)' 'subwrd(xline,8)
lats=subwrd(yline,6)' 'subwrd(yline,8)



'set lev 850'
'tmp=tmpprs'
'vwnd=vgrdprs'
'uwnd=ugrdprs'

'set gxout shaded'
'set mpdset hires'

'pi=3.14159'
'dtr=pi/180'
'a=6.37122e6'


* get dx and dy
'dy=cdiff(lat,y)*dtr*a'
'dx=cdiff(lon,x)*dtr*a*cos(lat*dtr)'
* calculating temp gradient
'dtdx=cdiff(tmp,x)/dx'
'dtdy=cdiff(tmp,y)/dy'

* calculating temp advection
'define tadv=(-uwnd*dtdx-vwnd*dtdy)'
'set grads off'
'set lon -77 -65'
'set lat 40 50'
'set gxout shaded'
'set clevs -8 -6 -4 -2 -1 1 2 4 6 8'
'set ccols 16 17 18 19 20 0 21 22 23 24 25'

*'colormaps -map b2r -custom -6 -4 -2 0 2 4 6 8 10'
'd tadv*3600'
'q w2xy -70.3097 43.6465'
* get PWM coordinates on the plot
x=subwrd(result,3)
y=subwrd(result,6)
'set line 14'
'draw mark 1 'x' 'y' 0.15'
'draw title 850 hPa Temp. Adv. 00Z Mar 14, 1993'
cbarn.gs
* get T advection value first
'set lon -70.3097'
'set lat 43.6465'
'd tadv*3600'
adv=subwrd(result,4)
say '===================='
say 'T advection = 'adv ' K/hr'
say '===================='


* ========================================
* Start calculation of 2b
* ========================================

* Use 800 and 900 hPa to calculate lapse rate
* Get temperature at 800
'set lev 800'
'tmp=tmpprs'
'set lon -70.3097'
'set lat 43.6465'
'd tmp'
T800=subwrd(result,4)

* Get temperature at 900
'set lev 900'
'tmp=tmpprs'
'set lon -70.3097'
'set lat 43.6465'
'd tmp'
T900=subwrd(result,4)

* Get temperature at 1000
'set lev 1000'
'tmp=tmpprs'
'set lon -70.3097'
'set lat 43.6465'
'd tmp'
T1000=subwrd(result,4)
say 'T800 = ' T800
say 'T900 = ' T900
say 'T1000= ' T1000

* Use the skew-T to find out the moist adiabatic rate
TM800=-14.8+273
TM900=-8.0+273

* get omega850 in hPa/s from the previous calculation
omega850=0.61*(-0.026)
say 'omega850 = ' omega850 ' hPa/s'

* Calculate the lapse rate for the environment and moist adiabat
DTDP   =( T800- T900)/(800-900)
DTDP_ma=(TM800-TM900)/(800-900)

say 'DTDP   = 'DTDP
say 'DTDP_ma= 'DTDP_ma
* Calculate the final local change of temperature in K/hr
omegaTerm= omega850*(DTDP-DTDP_ma)*3600
say 'omegaTerm= ' omegaTerm

dTdt= adv - omegaTerm
say 'local tendency of T= 'dTdt ' K/hr'