*****************************************************************
*         HORIZONTAL FORWARD AND BACKWARD TRAJECTORIES          *
*           Department of Astronomy and Meteorology             *
*                    University of Barcelona                    *
*                                                               *
*               (C) Bernat Codina - spring 2001                 * 
*                      bcodina@am.ub.es                         *
*****************************************************************
*  Modified by Ted Pollock (University of Alberta, email:
*  edwardp@ualberta.ca)  to incorporate vertical motions into trajectory
*  calculations. Summer 2007. 
*  The program interpolates u, v, and w between sigma levels.

  function main(arg)
    _At=6371229
    _PI=3.141592654
    _D2R=_PI/180
    _R2D=180/_PI

* Wind components in the ctl file. Must be changed if other
* notation is employed.
    _U='U'
    _V='V'
    _W='W'
    _H='H'
    _TER='TER'

    if(inici())
      return
    endif
    if (arg='')
      say 'Click on the map to start a trajectory.'
      say 'Alternatively, you can also issue:'
      say '        traj lon_ini lat_ini z_ini'
      say ' '
      'q pos'
      x=subwrd(result,3)
      y=subwrd(result,4)
      'q xy2w 'x' 'y
      LON0=subwrd(result,3)
      LAT0=subwrd(result,6)
    else
      LON0=subwrd(arg,1)
      LAT0=subwrd(arg,2)
      j=subwrd(arg,3)
      say j
    endif
*   (j) is the starting level, (k) is a marker for recording the number of iterations the loop has gone through, (SUM) keeps track of the sum of vertical distances traveled via d=E(w*t), (hgt) is the height of the current trajectory calculation. Done by adding the SUM to the height of the current sigma level. fbswitch=2 for backwards only.
     fbswitch=0
* (oj) needed only in the case the script switches from forward to backward trajectories, (up,dn) set other than zero only in the case the script just switched sigma levels
     oj=j
     up=0
     dn=0
* Model sigma levels, and sigma values are needed
     level.1=0.99500
     level.2=0.98500
     level.3=0.97000
     level.4=0.94500
     level.5=0.91000
     level.6=0.87000
     level.7=0.82500
     level.8=0.77500
     level.9=0.72500
     level.10=0.67500
     level.11=0.62500
     level.12=0.57500
     level.13=0.52500
     level.14=0.47500
     level.15=0.42500
     level.16=0.37500
     level.17=0.32500
     level.18=0.27500
     level.19=0.22500
     level.20=0.17500
     level.21=0.12500
     level.22=0.07500
     level.23=0.02500
     'set lev ' level.j
* Defining the colours used for the sigma levels along the trajectory. (1-sigma)*255
'set rgb 51 11 11 11'
'set rgb 52 22 22 22'
'set rgb 53 33 33 33'
'set rgb 54 44 44 44'
'set rgb 55 55 55 55'
'set rgb 56 66 66 66'
'set rgb 57 77 77 77'
'set rgb 58 88 88 88'
'set rgb 59 99 99 99'
'set rgb 60 110 110 110'
'set rgb 61 121 121 121'
'set rgb 62 132 132 132'
'set rgb 63 143 143 143'
'set rgb 64 154 154 154'
'set rgb 65 165 165 165'
'set rgb 66 176 176 176'
'set rgb 67 187 187 187'
'set rgb 68 198 198 198'
'set rgb 69 209 209 209'
'set rgb 70 220 220 220'
'set rgb 71 231 231 231'
'set rgb 72 242 242 242'
'set rgb 73 254 254 254'

* volta=1 for fwd+bck volta=2 for bck only
    volta=2
    while (volta<=2)
      if (volta=1)
        say 'Forward trajectory (in red):'
        signe=1
        _color=2
        Tm=-999999
        TM=_Tmax
      else
        say 'Backward trajectory (in blue):'
        signe=-1
        _color=4
        TM=999999
        Tm=1
      endif
* added
        'set lev ' level.oj
        j=oj
*
      LONt=LON0
      LATt=LAT0
      aflag=1
      SUM=0
      hgt=0
      Zop=0
      Zom=0
      BFACT=1
      TFACT=0
      t=_Tini
      say t' 'LONt' 'LATt
      'set T ' t
**************************************************
******DISTANCE CALCULATOR
      Ut=interp(_U,LONt,LATt)  
      Vt=interp(_V,LONt,LATt)
      Wt=interp(_W,LONt,LATt)
***********************************************
      Utm1=Ut
      Vtm1=Vt
      Wtm1=Wt
* Determines the heights of the initial sigma levels 
        if(j=1)
          'set lev ' level.j
          Zom=interp(_H,LONt,LATt)
          j=j+1
          'set lev ' level.j
          Zup1m=interp(_H,LONt,LATt)
          j=j+1
          'set lev ' level.j
          Zup2m=interp(_H,LONt,LATt)
          j=j-2
          Zdn1m=-9999
          Zdn2m=-99999
        endif
*
        if(j=2)
          'set lev ' level.j
          Zom=interp(_H,LONt,LATt)
          j=j+1
          'set lev ' level.j
          Zup1m=interp(_H,LONt,LATt)
          j=j+1
          'set lev ' level.j
          Zup2m=interp(_H,LONt,LATt)
          j=j-3
          'set lev ' level.j
          Zdn1m=interp(_H,LONt,LATt)
          j=j+1
          Zdn2m=-99999
        endif
*
        if(j>2&j<22)
          'set lev ' level.j
          Zom=interp(_H,LONt,LATt)
          j=j+1
          'set lev ' level.j
          Zup1m=interp(_H,LONt,LATt)
          j=j+1
          'set lev ' level.j
          Zup2m=interp(_H,LONt,LATt)
          j=j-4
          'set lev ' level.j
          Zdn2m=interp(_H,LONt,LATt)
          j=j+1
          'set lev ' level.j
          Zdn1m=interp(_H,LONt,LATt)
          j=j+1
        endif
*
        if(j=22)
          'set lev ' level.j
          Zom=interp(_H,LONt,LATt)
          j=j+1
          'set lev ' level.j
          Zup1m=interp(_H,LONt,LATt)
          j=j-3
          'set lev ' level.j
          Zdn2m=interp(_H,LONt,LATt)
          j=j+1
          'set lev ' level.j
          Zdn1m=interp(_H,LONt,LATt)
          j=j+1
          Zup2m=999999
        endif
*
        if(j=23)
          'set lev ' level.j
          Zom=interp(_H,LONt,LATt)
          j=j-1
          'set lev ' level.j
          Zdn1m=interp(_H,LONt,LATt)
          j=j-1
          'set lev ' level.j
          Zdn2m=interp(_H,LONt,LATt)
          j=j+2
          Zup1m=999999
          Zup2m=9999999
        endif
********************************
      while (t<TM & t>Tm)
        Ut=Utm1
        Vt=Vtm1
        Wt=Wtm1
        'set T ' t+signe
        LON0tm1=LONt
        LAT0tm1=LATt
        DR=1000000
        iter=0
        if(t=3)
          fi(1)
        endif
*        say startwhile************************************************
        while(DR>0.1 & iter<15)
          Umig=signe*(Utm1+Ut)/2
          Vmig=signe*(Vtm1+Vt)/2
          Wmig=signe*(Wtm1+Wt)/2
          wdist=Wmig*_DT
          'd sqrt('Umig'*'Umig'+'Vmig'*'Vmig')'
          V=subwrd(result,4)
          D=V*_DT
          'd atan2('Vmig','Umig')'
          alpha=subwrd(result,4)*_R2D
          rec=mou(LONt,LATt,D,alpha)
          LON1tm1=subwrd(rec,1)
          LAT1tm1=subwrd(rec,2)
          DR=dist(LON1tm1,LAT1tm1,LON0tm1,LAT0tm1)
************
*********
        if(j=1&SUM<0)
         SUM=0
        endif
        if(j=23&SUM>=0)
         SUM=-0.0000001
        endif
************
***********************************************
* Determines the heights of the sigma levels 

        if(j=1)
          'set lev ' level.j
          Zop=interp(_H,LON1tm1,LAT1tm1)
          j=j+1
          'set lev ' level.j
          Zup1p=interp(_H,LON1tm1,LAT1tm1)
          j=j+1
          'set lev ' level.j
          Zup2p=interp(_H,LON1tm1,LAT1tm1)
          j=j-2
          Zdn1p=-9999
          Zdn2p=-99999
*          say t' 'Siglevs' 'Zdn2p' 'Zdn1p' 'Zop' 'Zup1p' 'Zup2p' 'j
        endif
*
        if(j=2)
          'set lev ' level.j
          Zop=interp(_H,LON1tm1,LAT1tm1)
          j=j+1
          'set lev ' level.j
          Zup1p=interp(_H,LON1tm1,LAT1tm1)
          j=j+1
          'set lev ' level.j
          Zup2p=interp(_H,LON1tm1,LAT1tm1)
          j=j-3
          'set lev ' level.j
          Zdn1p=interp(_H,LON1tm1,LAT1tm1)
          j=j+1
          Zdn2p=-99999
*          say t' 'Siglevs' 'Zdn2p' 'Zdn1p' 'Zop' 'Zup1p' 'Zup2p' 'j
        endif
*
        if(j>2&j<22)
          'set lev ' level.j
          Zop=interp(_H,LON1tm1,LAT1tm1)
          j=j+1
          'set lev ' level.j
          Zup1p=interp(_H,LON1tm1,LAT1tm1)
          j=j+1
          'set lev ' level.j
          Zup2p=interp(_H,LON1tm1,LAT1tm1)
          j=j-4
          'set lev ' level.j
          Zdn2p=interp(_H,LON1tm1,LAT1tm1)
          j=j+1
          'set lev ' level.j
          Zdn1p=interp(_H,LON1tm1,LAT1tm1)
          j=j+1
*          say t' 'Siglevs' 'Zdn2p' 'Zdn1p' 'Zop' 'Zup1p' 'Zup2p' 'j
        endif
*
        if(j=22)
          'set lev ' level.j
          Zop=interp(_H,LON1tm1,LAT1tm1)
          j=j+1
          'set lev ' level.j
          Zup1p=interp(_H,LON1tm1,LAT1tm1)
          j=j-3
          'set lev ' level.j
          Zdn2p=interp(_H,LON1tm1,LAT1tm1)
          j=j+1
          'set lev ' level.j
          Zdn1p=interp(_H,LON1tm1,LAT1tm1)
          j=j+1
          Zup2p=999999
*          say t' 'Siglevs' 'Zdn2p' 'Zdn1p' 'Zop' 'Zup1p' 'Zup2p' 'j
        endif
*
        if(j=23)
          'set lev ' level.j
          Zop=interp(_H,LON1tm1,LAT1tm1)
          j=j-1
          'set lev ' level.j
          Zdn1p=interp(_H,LON1tm1,LAT1tm1)
          j=j-1
          'set lev ' level.j
          Zdn2p=interp(_H,LON1tm1,LAT1tm1)
          j=j+2
          Zup1p=999999
          Zup2p=9999999
*          say t' 'Siglevs' 'Zdn2p' 'Zdn1p' 'Zop' 'Zup1p' 'Zup2p' 'j
        endif

        if(Zop<0)
         fi(1)
        endif


        zdist=Zop-Zom

***************************************
* below: level switch, if the sum is great enough to switch up or down one
* or two levels. In this example there are 21 levels. The 1st, 2nd, 20th,
* and 21st (i.e. top 2 and bottom 2) need to be treated differently than
* levels in the middle.
       x=0
       TEMP=SUM-zdist+wdist
       if(j=1&x=0)
          if((TEMP+Zop)>=Zup2p)
            TEMP=TEMP+Zop-Zup2p
            x=2
           else
            if((TEMP+Zop)>=Zup1p)
              TEMP=TEMP+Zop-Zup2p
              x=1
            endif
          endif
        endif
*
        if(j=2&x=0)
          if((TEMP+Zop)>=Zup2p)
            TEMP=TEMP+Zop-Zup2p
            x=2
           else
            if((TEMP+Zop)>=Zup1p)
              TEMP=TEMP+Zop-Zup1p
              x=1
            endif
          endif
          if((TEMP+Zop)<=Zdn1p)
            TEMP=TEMP+Zop-Zdn1p
            x=-1
          endif
        endif
*
        if(j>2&j<22&x=0)
          if((TEMP+Zop)>=Zup2p)
            TEMP=TEMP+Zop-Zup2p
            x=2
           else
            if((TEMP+Zop)>=Zup1p)
              TEMP=TEMP+Zop-Zup1p
              x=1
            endif
          endif
          if((TEMP+Zop)<=Zdn2p)
            TEMP=TEMP+Zop-Zdn2p
            x=-2
           else
            if((TEMP+Zop)<=Zdn1p)
              TEMP=TEMP+Zop-Zdn1p
              x=-1
            endif
          endif
        endif
*
        if(j=22&x=0)
          if((TEMP+Zop)<=Zdn2p)
            TEMP=TEMP+Zop-Zdn2p
            x=-2
           else
            if((TEMP+Zop)<=Zdn1p)
              TEMP=TEMP+Zop-Zdn1p
              x=-1
            endif
          endif
          if((TEMP+Zop)>=Zup1p)
            TEMP=TEMP+Zop-Zup1p
            x=1
          endif
        endif
*
        if(j=23&x=0)
          if((TEMP+Zop)<=Zdn2p)
            TEMP=TEMP+Zop-Zdn2p
            x=-2
            if((TEMP+Zop)<=Zdn1p)
              TEMP=TEMP+Zop-Zdn1p
              x=-1
            endif
          endif
        endif

        p=j+x
        'set lev ' level.p
        
************
************
* Determines the heights of the sigma levels 
       if(x!=0)
        if(p=1)
          'set lev ' level.p
          Zop=interp(_H,LON1tm1,LAT1tm1)
          p=p+1
          'set lev ' level.p
          Zup1p=interp(_H,LON1tm1,LAT1tm1)
          p=p+1
          'set lev ' level.p
          Zup2p=interp(_H,LON1tm1,LAT1tm1)
          p=p-2
          Zdn1p=-9999
          Zdn2p=-99999
*          say t' 'Siglevs' 'Zdn2p' 'Zdn1p' 'Zop' 'Zup1p' 'Zup2p' 'j
        endif
*
        if(p=2)
          'set lev ' level.p
          Zop=interp(_H,LON1tm1,LAT1tm1)
          p=p+1
          'set lev ' level.p
          Zup1p=interp(_H,LON1tm1,LAT1tm1)
          p=p+1
          'set lev ' level.p
          Zup2p=interp(_H,LON1tm1,LAT1tm1)
          p=p-3
          'set lev ' level.p
          Zdn1p=interp(_H,LON1tm1,LAT1tm1)
          p=p+1
          Zdn2p=-99999
*          say t' 'Siglevs' 'Zdn2p' 'Zdn1p' 'Zop' 'Zup1p' 'Zup2p' 'j
        endif
*
        if(p>2&p<22)
          'set lev ' level.p
          Zop=interp(_H,LON1tm1,LAT1tm1)
          p=p+1
          'set lev ' level.p
          Zup1p=interp(_H,LON1tm1,LAT1tm1)
          p=p+1
          'set lev ' level.p
          Zup2p=interp(_H,LON1tm1,LAT1tm1)
          p=p-4
          'set lev ' level.p
          Zdn2p=interp(_H,LON1tm1,LAT1tm1)
          p=p+1
          'set lev ' level.p
          Zdn1p=interp(_H,LON1tm1,LAT1tm1)
          p=p+1
*          say t' 'Siglevs' 'Zdn2p' 'Zdn1p' 'Zop' 'Zup1p' 'Zup2p' 'j
        endif
*
        if(p=22)
          'set lev ' level.p
          Zop=interp(_H,LON1tm1,LAT1tm1)
          p=p+1
          'set lev ' level.p
          Zup1p=interp(_H,LON1tm1,LAT1tm1)
          p=p-3
          'set lev ' level.p
          Zdn2p=interp(_H,LON1tm1,LAT1tm1)
          p=p+1
          'set lev ' level.p
          Zdn1p=interp(_H,LON1tm1,LAT1tm1)
          p=p+1
          Zup2p=999999
*          say t' 'Siglevs' 'Zdn2p' 'Zdn1p' 'Zop' 'Zup1p' 'Zup2p' 'j
        endif
*
        if(p=23)
          'set lev ' level.p
          Zop=interp(_H,LON1tm1,LAT1tm1)
          p=p-1
          'set lev ' level.p
          Zdn1p=interp(_H,LON1tm1,LAT1tm1)
          p=p-1
          'set lev ' level.p
          Zdn2p=interp(_H,LON1tm1,LAT1tm1)
          p=p+2
          Zup1p=999999
          Zup2p=9999999
*          say t' 'Siglevs' 'Zdn2p' 'Zdn1p' 'Zop' 'Zup1p' 'Zup2p' 'j
        endif

       endif

***************************************
********************************
* Require factors for vertical interpolation
       if(x=0)
        if((SUM-zdist+wdist)>=0)
         TFACT=(SUM-zdist+wdist)/(Zup1p-Zop)
         BFACT=1-TFACT
        else
         BFACT=(SUM-zdist+wdist)/(Zdn1p-Zop)
         TFACT=1-BFACT
        endif 
       else
        if(TEMP>=0)
         TFACT=TEMP/(Zup1p-Zop)
         BFACT=1-TFACT
        else
         BFACT=TEMP/(Zdn1p-Zop)
         TFACT=1-BFACT
        endif 
       endif
******Calculate vertical height change due to vertically interpolated
* vertical wind
        if(TEMP>=0)
         if(p<23)
          WTMPx=interp(_W,LON1tm1,LAT1tm1)
          p=p+1
          'set lev ' level.p
          WTMPy=interp(_W,LON1tm1,LAT1tm1)
          p=p-1
          'set lev ' level.p
          Wtm1=((WTMPx*BFACT)+(WTMPy*TFACT))
         else
          Wtm1=interp(_W,LON1tm1,LAT1tm1)
         endif
        endif
        if(TEMP<0)
         if(p>1)
          WTMPx=interp(_W,LON1tm1,LAT1tm1)
          p=p-1
          'set lev ' level.p
          WTMPy=interp(_W,LON1tm1,LAT1tm1)
          p=p+1
          'set lev ' level.p
          Wtm1=((WTMPy*BFACT)+(WTMPx*TFACT))
         else
          Wtm1=interp(_W,LON1tm1,LAT1tm1)
         endif
        endif
***********************
************************
********Modified to interpolate in the vertical
          if(TEMP>=0)
           if(p<23)
            Utm1a=interp(_U,LON1tm1,LAT1tm1)
            Vtm1a=interp(_V,LON1tm1,LAT1tm1)
            p=p+1
            'set lev ' level.p
            Utm1b=interp(_U,LON1tm1,LAT1tm1)
            Vtm1b=interp(_V,LON1tm1,LAT1tm1)
            p=p-1
            'set lev ' level.p
            Utm1=(BFACT*Utm1a)+(TFACT*Utm1b)
            Vtm1=(BFACT*Vtm1a)+(TFACT*Vtm1b)
           else
            Utm1=interp(_U,LON1tm1,LAT1tm1)
            Vtm1=interp(_V,LON1tm1,LAT1tm1)
           endif 
          else
           if(p>1)
            Utm1a=interp(_U,LON1tm1,LAT1tm1)
            Vtm1a=interp(_V,LON1tm1,LAT1tm1)
            p=p-1
            'set lev ' level.p
            Utm1b=interp(_U,LON1tm1,LAT1tm1)
            Vtm1b=interp(_V,LON1tm1,LAT1tm1)
            p=p+1
            'set lev ' level.p
            Utm1=(TFACT*Utm1a)+(BFACT*Utm1b)
            Vtm1=(TFACT*Vtm1a)+(BFACT*Vtm1b)
           else
            Utm1=interp(_U,LON1tm1,LAT1tm1)
            Vtm1=interp(_V,LON1tm1,LAT1tm1)
           endif 
          endif
*********************
*         say iter' 'p' 'j' 'zdist' 'wdist' 'Utm1' 'Vtm1' 'Wtm1
*        say TEMP' 'SUM' 'TFACT' 'BFACT' 'j
***********
          LON0tm1=LON1tm1
          LAT0tm1=LAT1tm1
          iter=iter+1
        endwhile
*******************************************************************
*          say t' 'Sigs' 'Zdn2m' 'Zdn1m' 'Zom' 'Zup1m' 'Zup2m' 'j
*          say t' 'Sigs' 'Zdn2p' 'Zdn1p' 'Zop' 'Zup1p' 'Zup2p' 'j

        j=j+x
        'set lev ' level.j

        if(x!=0)
         SUM=TEMP
         else
         SUM=SUM-zdist+wdist
        endif

        hgt=Zop+SUM
        terhgt=interp(_TER,LON1tm1,LAT1tm1)
        say t' 'LON1tm1' 'LAT1tm1' 'j' 'hgt' 'terhgt
*        say t' 'Zdn1m' 'Zom' 'Zup1m' 'SUM
*        say done loop done loop done loop done loop done loop
********Switch to determine whether back trajectory originates from
********Caspian sea
        if(LON1tm1<=60&LAT1tm1>=39&t>289&t<361&j>=12&aflag=1)
         aflag=0
         say aflag
        endif

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

        if (iter=15)
*          say '--> Lack of convergence: DR='DR' m.'
        endif

        tx=t
        linia(LONt,LATt,LON1tm1,LAT1tm1,j,tx,signe)

*****  Draw line and
****** END DISTANCE CALCULATOR

*****************************
*********
        if(j=1&SUM<0)
         SUM=0
        endif
        if(j=23&SUM>=0)
         SUM=-0.0000001
        endif
************

        LONt=LON1tm1
        LATt=LAT1tm1
        Zup2m=Zup2p
        Zup1m=Zup1p
        Zom=Zop
        Zdn1m=Zdn1p
        Zdn2m=Zdn2p
        
        t=t+signe
      endwhile
      volta=volta+1
    endwhile
    fi(0)
*    say hgt
  return

  function mou(lon0,lat0,dist,alpha)
* Coordinates of a point located a distance dist (in meters) and angle
* alpha (in degrees and mathematical convention) from the point (lon0,lat0).
    if (dist<1);return(lon0' 'lat0);endif
    lon0=lon0*_D2R
    lat0=lat0*_D2R
    alpha=90-alpha
    if (alpha<0)
      alpha=360+alpha
    endif
    if (alpha>=360)
      alpha=alpha-360
    endif
    A=alpha*_D2R
    b=dist/_At
    c=_PI/2-lat0
    'd acos(cos('b')*cos('c')+sin('b')*sin('c')*cos('A'))'
    a=subwrd(result,4)
    lat1=(_PI/2-a)*_R2D
    'd asin(sin('b')*sin('A')/sin('a'))'
    B=subwrd(result,4)
    lon1=(lon0+B)*_R2D
  return(lon1' 'lat1)

  function interp(camp,lon,lat)
* Interpolates variable camp in (lon,lat).
    'q w2gr 'lon' 'lat
    x=subwrd(result,3);y=subwrd(result,6)
    x.2=int(x);y.3=int(y)
    x.3=x.2;y.4=y.3
    x.1=x.2+1;y.1=y.3+1
    x.4=x.1;y.2=y.1
    num=0;den=0
    i=1
    while(i<=4)
      if(x.i<_Xmin|x.i>_Xmax|y.i<_Ymin|y.i>_Ymax)
        fi(1)
      endif
      'set x 'x.i
      'set y 'y.i
      'q gr2w 'x.i' 'y.i
      rec=sublin(result,1)
      lon.i=subwrd(rec,3)
      lat.i=subwrd(rec,6)
      'd 'camp
      rec=sublin(result,1)
      if(subwrd(rec,1)!='Result')
        rec=sublin(result,2)
      endif
      z.i=subwrd(rec,4)
      d.i=dist(lon,lat,lon.i,lat.i)
      if (d.i=0)
        return (z.i)
        break
      endif
      num=num+1/d.i*z.i
      den=den+1/d.i
      i=i+1
    endwhile
  return (num/den)

  function dist(lon1,lat1,lon2,lat2)
* Distance between two points on the Earth surface
    phi=lon1*_D2R;theta=(90-lat1)*_D2R
    'd sin('theta')*cos('phi')'
    x1=subwrd(result,4)
    'd sin('theta')*sin('phi')'
    y1=subwrd(result,4)
    'd cos('theta')'
    z1=subwrd(result,4)
    phi=lon2*_D2R;theta=(90-lat2)*_D2R
    'd sin('theta')*cos('phi')'
    x2=subwrd(result,4)
    'd sin('theta')*sin('phi')'
    y2=subwrd(result,4)
    'd cos('theta')'
    z2=subwrd(result,4)
    x=y1*z2-y2*z1
    y=x2*z1-x1*z2
    z=x1*y2-x2*y1
    d2=x*x+y*y+z*z
    'd asin(sqrt('d2'))'
  return (subwrd(result,4)*_At)

  function int(x)
* Integer part of x
    i=1
    while(i<17)
      if(substr(x,i,1)=".");break;endif
      i=i+1
    endwhile
    if (i=17); return(x);endif
    if (i=1); return('0');endif
  return (substr(x,1,i-1))

  function linia(lon0,lat0,lon1,lat1,phi,TT,sing)
* Draws a line limited by two circles in _color
  if ((lon1-lon0)>2.5|(lon0-lon1)>2.5)
    fi(1)
  endif;
    pi=phi+50
    _color=pi
    'set line '_color
* modified colour
    'q w2xy 'lon0' 'lat0
    x0=subwrd(result,3)
    y0=subwrd(result,6)
    'q w2xy 'lon1' 'lat1
    x1=subwrd(result,3)
    y1=subwrd(result,6)
    'draw line 'x0' 'y0' 'x1' 'y1
    'set line '2
    if(TT=433|TT=505|TT=577|TT=649)
      'draw mark 3 'x0' 'y0' 0.075'
      'set line '_color
    endif
    'set line '9
    if(TT=289|TT=361|TT=721|TT=793)
      'draw mark 3 'x0' 'y0' 0.075'
      'set line '_color
    endif
    'set line '0
    if(TT=73|TT=145|TT=217)
      'draw mark 3 'x0' 'y0' 0.075'
      'set line '_color
    endif
    return
  
  function inici
* Some verifications and variable initialization
  'q w2gr 1 1'
  if(subwrd(result,3)='environment')
    say 'Please open first a file and draw a chart.'
    return(1)
  endif
  'q dims'
  rec1=sublin(result,2)
  rec2=sublin(result,3)
  rec3=sublin(result,4)
  rec4=sublin(result,5)
  if(subwrd(rec1,3)='fixed' | subwrd(rec2,3)='fixed')
    say 'X and Y dimensions must be opened:'
    say rec1
    say rec2
    return(1)
  endif
  if(subwrd(rec3,3)!='fixed')
    say 'Vertical dimension must be fixed:'
    say rec3
    return(1)
  endif
  if(subwrd(rec4,3)!='fixed')
    say 'Time dimension must be fixed'
    say rec4
    return(1)
  endif
  _Xmin=subwrd(rec1,11)+0.1
  _Xmax=subwrd(rec1,13)-0.1
  _Ymin=subwrd(rec2,11)+0.1
  _Ymax=subwrd(rec2,13)-0.1
  _Tini=subwrd(rec4,9)
  'q file'
  rec=sublin(result,5)
  _Tmax=subwrd(rec,12)
  rec=sublin(result,2)
  fitxer=subwrd(rec,2)
  rc=0
  while(rc=0)
    rec=read(fitxer)
    rc=sublin(rec,1)
    lin=sublin(rec,2)
    if (subwrd(lin,1)='tdef' | subwrd(lin,1)='TDEF')
      interval=subwrd(lin,5)
      break
    endif
  endwhile
  if (rc!=0)
    say 'TDEF not found in 'fitxer
    return(1)
  endif
  rc=close(fitxer)
  i=1
  factor=0
  while(i<10)
    a=substr(interval,i,2)
    if(a='mn'|a='MN'|a='Mn')
      factor=60
      break
    endif
    if(a='hr'|a='HR'|a='Hr')
      factor=3600
      break
    endif
    i=i+1
  endwhile
  if(factor>0)
    _DT=factor*substr(interval,1,i-1)
  else
    say 'HR or MN missing in line TDEF of file 'fitxer
    return(1)
  endif
  return (0)

* Function to determine the height of the ground along the trajectory
  function HGTERR(LONA,LATA)
     'set lat ' LATA
     'set lon ' LONA
     'd ter'
     HterRRR=subwrd(result,10)
   return(HterRRR)

*Functions for determining the heights of sigma surfaces
   function hsig(LONA,LATA,level.p,ZZ)
     'set lev ' level.p
     'set lat ' LATA
     'set lon ' LONA
*     'd h'
*     Hsig=subwrd(result,10)
     Fsig=interp(ZZ,LONA,LATA)
   return(Fsig)
   
  function fi(i)
* Things are left as they were before execution of traj
    'set x ' _Xmin' '_Xmax
    'set y ' _Ymin' '_Ymax
    'set t ' _Tini
* Does anybody know how to nicely quit a GrADS script?
    if (i=1)
      say 'The trajectory can't be completed.'
      say 'Press Ctrl+C and Enter to quit.'
      pull nothing
      fi(1)
    endif
  return