* This program computs the residue of the thermodynamic energy equation
* The residue (or diabatic heating which is in K/s) is necessary to the comput of 
* genertion terms in the Lorenz energy cycle

'reinit'
'open uvtw.ctl'

***********************************************************
'q file'
rc=sublin(result,5)
nt=subwrd(rc,12)
nl=subwrd(rc,9)
nx=subwrd(rc,3)
ny=subwrd(rc,6)
rc=sublin(result,6)
'q time'
td=subwrd(result,3)
ano=substr (td, 9, 4)

***********************************************************

'set fwrite thermo.bin'
'set gxout fwrite'
'set x 1 360 '
'set y 1 181 '

* Note: Repeat 1000 and 100 to permit program diferentiate z levels.

string=' 1000 1000 975 950 925 900 850 800 750 700 650 600 550 500 450 400 350 300 250 200 150 100 100 '

*aa      --> m
*Rgas    --> 287.05 J/Kg K
*Cp      --> 1004   J/Kg K
*dt      --> s
*T       --> K
*u       --> m/s
*v       --> m/s
*P       --> Pa
*W       --> Pa/s

'define aa=6.37e6'
'define pi=2*asin(1)'
'define rd=pi/180'
'define clat=cos(lat*rd)'
'define dx=cdiff(lon,x)*rd'
'define dy=cdiff(lat,y)*rd'
'define kp=287.05/1004'
'define dt=3600*6*2'

k = 1
l = 2

while (k <= 36)

while (l <= 21+1)

'set t 'k

l1=l-1
l2=l+1

nivel  = subwrd(string,l)
nivel1 = subwrd(string,l1)
nivel2 = subwrd(string,l2)

* As vertical velocity is in units of Pa/sec and the levels (in CTL and nc data) 
* are in hPa, we must multiply the levels (hPa) by 100 to convert it in Pa units 
* (the same unit of vertical velocity)

'define p1='nivel1'*100'
'define p2='nivel2'*100'

say nivel' 'nivel1' 'nivel2' 'k' 'l

'set lev 'nivel

* Computing Temperature tendency term
*_____________ dT/dt _______________*
* Units must be in K/s

k1=k-1
k2=k+1


if(k!=1 | k!=nt)
'define ttend = (T(t='k2')-T(t='k1'))/dt'
endif

if(k=1)
'define ttend = (T(t='k2')-T(t='k'))/(dt/2)'
say 'k=1'
endif


if(k=nt)
'define ttend = (T(t='k')-T(t='k1'))/(dt/2)'
say 'k=nt'
endif

* Computing Temperature horizontal advection
* in hespherical coordinates
*_____________ V Grad(T) _______________*
* Units must be in K/s

'define advh = ((u*cdiff(T,x))/(clat*dx) + (v*cdiff(T,y))/dy)/aa'

* Computing Temperature vertical advection
* in hespherical coordinates
*_________ w Grad(T) or wdT/dp __________*
* Units must be in K/s

'define advv = w*(t(lev='nivel2')-t(lev='nivel1'))/(p2-p1)'

* Computing alfa term
* in hespherical coordinates
* _____RTW/CpP _____________*
* Units must be in K/s

'define qd   = -kp *  (w*(t))/(lev*100)'

* Computing ther esidual term or the diabatic term
* in hespherical coordinates
* _____ Q _____________*
* Units must be in K/s

'define res   = ttend + advh + advv + qd'

* Recording the diabatic term in units of K/day

'd ttend*86400'
'd advh*86400'
'd advv*86400'
'd qd*86400'
'd res*86400'



k = k
l = l + 1

endwhile

l = 2
k = k + 1

endwhile

'disable fwrite'

**********************************************************
write (thermo.ctl,'dset thermo.bin')
write (thermo.ctl,'undef  -9.96921e+36')
write (thermo.ctl,'title NCEP 1x1 grau')
write (thermo.ctl,"xdef "360" linear    0.0 "1.0"")
write (thermo.ctl,"ydef "181" linear  -90.0 "1.0"")
write (thermo.ctl,"zdef "21" levels  "1000" "975" "950" "925" "900" "850" "800" "750" "700" "650" "600" "550" "500" "450" "400" "350" "300" "250" "200" "150" "100" ")
write (thermo.ctl,"tdef "36" linear "00Z21MAR2004" "06hr" "   )
write (thermo.ctl,'vars  5')
write (thermo.ctl,'ttend 21  99  temperature tendency  [K/dia] ')
write (thermo.ctl,'advh  21  99  horiz. temperature advection  [K/dia] ')
write (thermo.ctl,'advv  21  99  vertical temperature advection  [K/dia] ')
write (thermo.ctl,'qd    21  99                                  [K/dia] ')
write (thermo.ctl,"res "21"  99  residuo of diabatic heating  (Q)  [K/dia] ")
write (thermo.ctl,'endvars')
**********************************************************

quit