'reinit'
'open rad_com_1.ctl'


* Function Arguments:
*    sndtemp - temperature data (Celsius) as a function of pressure
*    snddewp - dewpoint data (Celsius) as a function of pressure      
*
* NOTE:  Make sure to set the vertical range of the plot before running.
*        I.e., "SET LEV 1050 150", for example.   This does not have to 
*        be limited to the pressure range of your data.

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

*** VARIAVEIS A CONSIDERAR ***
'define sndtemp = TC'
'define snddewp = TD'
'define startp = 1000'
'define endp = 100'

*** SERÁ?????? ***
'define pmax =  1000'
'define pmin = 100'

*** DOMINIO ***
'set lon -7.75 -7.05'
'set lat 38.05 38.75'
'set lev 1000 100'



dis_t = 1
while (dis_t < 72)
    'set t 'dis_t
*************************************************************************

pw=precipw(sndtemp,snddewp,pmax,pmin)

*************************************************************************
dis_t = dis_t + 1

endwhile


function precipw(sndtemp,snddewp,startp,endp)

*-----------------------------------------------------------------------
* Function to calculate the precipitable water (cm) in a sounding
* starting at pressure level startp and ending at pressure level endp.
*-----------------------------------------------------------------------

ppold=-999
ttold=-999
ddold=-999
delp=10
Int=0
mix=0 
pres=startp
logpp=log(pres)
logppm=log(pres-delp)
while (pres >= endp)
   tt=interp(sndtemp,pres)
   dd=interp(snddewp,pres)
   if (tt>-900 & ttold>-900 & dd>-900 & ddold>-900) 
      e=satvap2(dd)
      mix=mixratio(e,pres)
      mixavg=(logpp*mix+logppm*mixold)/(logpp+logppm)
      Int=Int+1.020408*mixavg*delp
   endif
   ttold=tt
   ddold=dd
   ppold=pp
   mixold=mix
   pres=pres-delp
   logpp=logppm
   logppm=log(pres-delp)
endwhile

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

