% C *******************************************************************
% C *******************************************************************
%                   subroutine traj(ut0,u,v)
% C *******************************************************************
% C *******************************************************************
% C
% C set parameters to match size of input files
%       parameter (ntraj=20,izut=20,imx=73,kmx=25,jmx=36)
% C Three components of the wind
%       real u(imx,kmx,jmx),v(imx,kmx,jmx),w(imx,kmx,jmx)
% C Mass mixing ratio
%       real poo2(imx,kmx,jmx),poo(imx,kmx,jmx)
% C Temperature
%       real T(imx,kmx,jmx)
% C ion velocities
%       real uvidyn(imx,jmx,2),velew(imx,jmx),velns(imx,jmx)
% C packet height and starting pressure surface
%       real zpkt(izut),azpkt(izut)
% C parcel local time
%       real xlt(izut)
% C packet velocities
%       real upkt(izut),vpkt(izut),wpkt(izut)
% C packet latitude and longitude
%       real alat(izut),alon(izut)
%       real alata(izut),alona(izut)
% C kuhlp? packet zonal coordinate
%       real kuhlp2,ipkt,ipku
% C packet meridional coordinate 
%       real jpkt,jpku
% C packet vertical coordinate
%       real kpkt,kpku
%       real utpkt(izut)
%       real utpkt2(izut)
%       common /trajin1/alons,alats,akstrt,zpkt
% C When Jeomgn=1 alon and alat are in geomagnetic coordinates
% C alata and alona are in geographic coordinates
%       common /trajout1/alon,alat,alona,alata,xlt
% C Also save the winds/ion drifts? along the trajectory path
% C      common /trajout2/upkt,vpkt,wpkt
% C longitude grid
%       real glon(73)
% C latitudfe grid
%       real glat(36)
      ibck=1;
      DTR=atan(1.0)/45;
      RTD=1./DTR;
      jeomgn = 1;
      astp = 1.
      timstp=240.
      re = 6371.E+5;
      for i = 1:144
       glon(i) = -180. + (i-1)*2.5
      end
      for j = 1:72
       glat(j) = -88.75 + (j-1)*2.5
      end
% C
% C+++++++++++++++++++++++++++++++++++++++++++++
% C  CALCULATE TRAJECTORY LOCI
% C+++++++++++++++++++++++++++++++++++++++++++++
      JUT=IU+1;
      UTOO=UTO;
      UT = UT0;
% C
% C   IF FIRST POINT INITIALISE
%       IF(JUT .EQ. 1)GO TO 2021
% C   IF UT=SPECIAL SET IN DATA STATEMENT, INITIALISE
% C   OTHERWISE SKIP INITIALISATION:
%       GO TO 2024
%  2021 CONTINUE
      KUHLP2=1;
% C Read in the initial conditions for pressure surface
      ZPKT(JUT)=AKSTRT;
      AZPKT(JUT)=ZPKT(JUT);
% C FIND NEAREST NEIGHBOUR POINTS using the initial conditions for lat lon
      IPK1 = (ALONS+180.)/2 + 1;
      IPK1 = MOD(IPK1+144,144)
      IF(IPK1 .EQ. 0)IPK1=IPK1 + 144
      IPKHOLD = IPK1
      IPK2 = MOD(IPK1+1,144)
      IF(IPK2 .EQ. 0)IPK2 = IPK2 + 144
      JPK1 = (ALATS+92.5)/5.
      JPKHOLD = JPK1
      JPK2 = JPK1+1
      KPK1 = AKSTRT
      KPK2 = KPK1 + 1
% C IPK3 and IPK4 are used to make transpolar calculations easy
      IPK3 = IPK1
      IPK4 = IPK2
      MFLAG = 0
% C The next 2 if blocks are invoked for points close to the pole
      IF(JPK1 .EQ. 0) THEN
         JPK1 = 1
         MFLAG = 1
         IPK3 = MOD(IPK1+37,144)
         IF(IPK3 .EQ. 0)IPK3=IPK3 + 144
         IPK4 = MOD(IPK1+72,144)
         IF(IPK4 .EQ. 0)IPK4 = IPK4+144
        ENDIF
      IF(JPK2 .GT. 72) THEN
         JPK2 = 72
         MFLAG = 2
         IPK3 = MOD(IPK1+37,144)
         IF(IPK3 .EQ. 0)IPK3=IPK3 + 144
         IPK4 = MOD(IPK1+72,144)
         IF(IPK4 .EQ. 0)IPK4 = IPK4+144
        ENDIF
% C Find U,V,W near the trajectory path and allow for transpolar points
      DO 2027 JJJ=1,6
      DO 2027 KKK=1,4
      DO 2027 III=1,6
      I=IPKHOLD + III - 3
      IF(I .GT. 144)I=I-144
      IF(I .LT. 1)I=I+144
      J=JPKHOLD + JJJ - 3
      AJPT = 1.
      IF(J .GT. 72)THEN
        I = MOD(IPK1+37,144) + 3 - III
        AJPT = -1.
      ENDIF
      IF(J .LT. 1)THEN
        I = MOD(IPK1+37,144) + 3 - III
        AJPT = -1.
      ENDIF
      IF(I .GT. 144)I=I-144
      IF(I .LT. 1)I=I+144
      IF(J .GT. 72)J=73-J
      IF(J .LT. 1)J=1-J
      K=KPK1+KKK-2
      WRITE(6,5035)U(I,K,J)
 5035 FORMAT(//1X,'U1 ',F8.2/)
% C ADJUST VERTICAL WINDS using the mean molecular mass and scale height
      AMBAR=1./((POO2(I,K,J)/32.)+(POO(I,K,J)/16.)+((1.-POO2(I,K,J)-
     1         POO(I,K,J))/28.))
      SHGT=8314.*T(I,K,J)/(AMBAR*8.7)
      W(I,K,J)=W(I,K,J)/SHGT
% C  Ion velocities must be used for ion trajectories
      IF(JEOMGN .EQ. 1)THEN
      ENDIF
 2027 CONTINUE
% c  calculate weighting terms for nearest neighbour points:
      X1 = GLON(IPK2) - ALONS
      Y1 = GLAT(JPK2) - ALATS
      IF (MFLAG .GT.  0) THEN
        Y1 = 5. - GLAT(JPK1) + ALATS
        X1 = 5. - ALONS + GLON(IPK1)
        ENDIF
      WRITE(6,5040)X1,Y1
%  5040 FORMAT('X1 X2',2F6.2)
      S1 = SQRT( (5.0-X1)**2 + (5.0-Y1)**2  + 0.0000001)
      S2 = SQRT( (5.0-X1)**2 + Y1**2  + 0.000001)
      S3 = SQRT( X1**2 + (5.0-Y1)**2 +0.000001)
      S4 = SQRT( X1**2 + Y1**2+0.000001)
      S5 = SQRT( (5.0-X1)**2 + (5.0-Y1)**2 +0.000001)
      S6 = SQRT( (5.0-X1)**2 + Y1**2 + 0.000001)
      S7 = SQRT( X1**2 + (5.0-Y1)**2 + 0.000001)
      S8 = SQRT( X1**2 + Y1**2 + 0.000001)
      AAA = 1.0 / ( (1./S1) + (1./S2) + (1./S3) + (1./S4))
      WA = AAA/S1
      WB = AAA/S2
      WC = AAA/S3
      WD = AAA/S4
% C Get vertical information
      Z1 = FLOAT(KPK2)-ZPKT(JUT)
      W1 = AAA/S1 *(Z1)
      W2 = AAA/S2 *(Z1)
      W3 = AAA/S3 *(Z1)
      W4 = AAA/S4 *(Z1)
      W5 = AAA/S5 *(1.-Z1)
      W6 = AAA/S6 *(1.-Z1)
      W7 = AAA/S7 *(1.-Z1)
      W8 = AAA/S8 *(1.-Z1)
      UPKT(JUT) = (W1*U(IPK1,KPK1,JPK1)+W2*
                     U(IPK3,KPK1,JPK2) +
                     W3*U(IPK2,KPK1,JPK1)+W4*
                     U(IPK4,KPK1,JPK2) +
                     W5*U(IPK1,KPK2,JPK1)+W6*
                     U(IPK3,KPK2,JPK2) +
                     W7*U(IPK2,KPK2,JPK1)+W8*
                     U(IPK4,KPK2,JPK2))
      VPKT(JUT) = (W1*V(IPK1,KPK1,JPK1)+W2*
                     V(IPK3,KPK1,JPK2) +
                     W3*V(IPK2,KPK1,JPK1)+W4*
                     V(IPK4,KPK1,JPK2) +
                     W5*V(IPK1,KPK2,JPK1)+W6*
                     V(IPK3,KPK2,JPK2) +
                     W7*V(IPK2,KPK2,JPK1)+W8*
                     V(IPK4,KPK2,JPK2))
      WPKT(JUT) = (W1*W(IPK1,KPK1,JPK1)+W2*W(IPK3,KPK1,JPK2)+
                     W3*W(IPK2,KPK1,JPK1)+W4*W(IPK4,KPK1,JPK2)+
                     W5*W(IPK1,KPK2,JPK1)+W6*W(IPK3,KPK2,JPK2)+
                     W7*W(IPK2,KPK2,JPK1)+W8*W(IPK4,KPK2,JPK2))
                     
      IF(JEOMGN .EQ. 1)THEN
      WPKT(JUT) = 0.
      ENDIF
      IPKT=IPK1
      IPKU=IPK1
      JPKT=(ALATS +92.5)/5.
      JPKU=JPK1
      KPKT=KPK1
      KPKU=KPK1

      UTPKT(JUT)=UT
      AUSTP = ASTP
      IF(IBCK .EQ. 1)AUSTP = -1.
      UTPKT2(JUT)=UT
      ALAT(JUT) = ALATS
      ALON(JUT) = ALONS
      AGLAT = ALATS
      AGLON = ALONS
      ALONA(JUT) = ALONS
      ALATA(JUT) = ALATS
% C For ion trajectories or neutral trajectories in geomagnetic coords
      IF(JEOMGN .EQ. 1 .OR. KEOMGN .EQ. 1)THEN
        CALL GTM (ALAT(JUT)*DTOR,ALON(JUT)*DTOR,CHMAG1*DTOR,
     1 CHMAG2*DTOR,A1,A2,DIP1,DEC1,W,1)
        ALON(JUT) = A2/DTOR
        ALAT(JUT) = A1/DTOR
      ENDIF
% C Dhour is set to allow for local time displacement between geomagnetic
% C and geographic poles
      ZKP1 = AKSTRT
      DHOUR=0.
      CHGHR = 4.667
      IF(JEOMGN .EQ. 1 .OR. KEOMGN .EQ. 1)THEN
      DHOUR = CHGHR
      IF (ALATS/ABS(ALATS) .LT. 0.) THEN
      CHGHR = -8.467
      DHOUR = CHGHR
      ENDIF
      ENDIF
      XLT(JUT)  = UT + ALON(JUT)/15. - DHOUR
      IF(XLT(JUT) .GT. 24.)XLT(JUT)=XLT(JUT)-24.
      IF(XLT(JUT) .LT. 0. )XLT(JUT)=XLT(JUT)+24.
      IHEMI=1
      IF(ALAT(JUT) .LT. 0.)IHEMI=-1
 2024 CONTINUE
C
      GO TO 2034
 2030 CONTINUE
C find next trajectory point,assume height=300km for now AGB
      REI = RE/1.E2 + 300.*1.E3
      IF(KUHLP2 .EQ. 1)GO TO 2036
C  SOLVE SPHERICAL TRIANGLE TO OBTAIN COORDS OF NEW TRAJ POINT
      ALNGX=AGLON
      ALATGX=AGLAT
      UPOLD=UPKT(JUT-1)
      VPOLD=VPKT(JUT-1)
      WPOLD=WPKT(JUT-1)
      UPNEW=UPOLD
      VPNEW=VPOLD
      WPNEW=WPOLD
      CZKP1=ZKP1
      IHEM=AGLAT/ABS(AGLAT)
      AHEM=FLOAT(IHEM)
      DO 6103 JNJK=1,2,1
      ALNGX=AGLON
      ALATGX=AGLAT
      ZKP1 = CZKP1
      UPLEN = (UPOLD + UPNEW)/2.
      VPLEN = (VPOLD + VPNEW)/2.
      WPLEN = (WPOLD + WPNEW)/2.
      ATMSTP = TIMSTP
% C  SOLVE SPHERICAL TRIANGLE TO OBTAIN COORDS OF NEW TRAJ POINT
      VMAG = SQRT(VPLEN**2 + UPLEN**2)
% C ANGULAR DISTANCE TRAVELLED:
      BETA = VMAG * ATMSTP / REI
      ALMONE = -(ALNGX - AGLON) * DTR
      UCLEN = UPLEN * COS(ALMONE) - AHEM*VPLEN * SIN(ALMONE)
      VCLEN = AHEM*UPLEN * SIN(ALMONE) + VPLEN * COS(ALMONE)
% C If backward trajectories are being calculated, wind velocities affect
% C  the parcel in the opposite sense
      ABCK=1.
      IF(IBCK .EQ. 1)ABCK=-1.
      ALPHA = ATAN2(ABCK*UCLEN,ABCK*VCLEN)
      THNEW = ASIN( SIN(ALATGX*DTR) * COS(BETA)
     1            + COS(ALATGX*DTR) *SIN(BETA) * COS(ALPHA) )
C
      XXX = SIN(BETA) * SIN(ALPHA)
      YYY = COS(ALATGX*DTR) * COS(BETA)
     1      - ( SIN(ALATGX*DTR) * SIN(BETA) * COS(ALPHA) )
      GLNEW = ALNGX + ATAN2(XXX,YYY)/DTR
      ALATGX = THNEW/DTR
      ALNGX = GLNEW
      ZKP1=ZKP1+ABCK*(TIMSTP/2.0)*(WPLEN)
C
C
      IF(ALATGX .GT. 90.)ALNGX=180.+ALNGX
      IF(ALATGX .GT. 90.)ALATGX=180.-ALATGX
      IF(ALATGX .LT. -90.)ALNGX=180.+ALNGX
      IF(ALATGX .LT.-90.)ALATGX=-180.-ALATGX
      IF(ALNGX .LT. -180.)ALNGX=ALNGX+360.
      IF(ALNGX .GT. 180.)ALNGX=ALNGX-360.
C Define the parcel location in i,j (longitude and latitude) grids
      IPK1 = (ALNGX+180.)/5. + 1
      IPK1 = MOD(IPK1+144,144)
      IPKHOLD = IPK1
      IF(IPK1 .EQ. 0)IPK1=IPK1 + 144
      IPK2 = MOD(IPK1+1,144)
      IF(IPK2 .EQ. 0)IPK2 = IPK2 + 144
      JPK1 = ((ALATGX)+92.5)/5.
      JPKHOLD = JPK1
      JPK2 = JPK1+1
      KPK1 = ZKP1
      KPK2 = KPK1+1
      IPK3 = IPK1
      IPK4 = IPK2
      MFLAG = 0
C  This section deals with transpolar points
      IF(JPK1 .EQ. 0) THEN
         JPK1 = 1
         MFLAG = 1
         IPK3 = MOD(IPK1+37,144)
         IF(IPK3 .EQ. 0)IPK3=IPK3 + 144
         IPK4 = MOD(IPK1+72,144)
         IF(IPK4 .EQ. 0)IPK4 = IPK4+144
        ENDIF
      IF(JPK2 .GT. 72) THEN
         JPK2 = 72
         MFLAG = 2
         IPK3 = MOD(IPK1+37,144)
         IF(IPK3 .EQ. 0)IPK3=IPK3 + 144
         IPK4 = MOD(IPK1+72,144)
         IF(IPK4 .EQ. 0)IPK4 = IPK4+144
        ENDIF
C Find U,V,W near the trajectory path
      DO 2039 JJJ=1,6
      DO 2039 KKK=1,4
      DO 2039 III=1,6
      I=IPKHOLD + III - 3
      IF(I .GT. 144)I=I-144
      IF(I .LT. 1)I=I+144
      J=JPKHOLD + JJJ - 3
C  AJPT inverts the sense of the velocity for transpolar points
      AJPT = 1.
C  These sections deal with transpolar points
      IF(J .GT. 72)THEN
        I = MOD(IPK1+37,144) + 3 - III
        AJPT = -1.
      ENDIF
      IF(J .LT. 1)THEN
        I = MOD(IPK1+37,144) + 3 - III
        AJPT = -1.
      ENDIF
      IF(I .GT. 144)I=I-144
      IF(I .LT. 1)I=I+144
      IF(J .GT. 72)J=73-J
      IF(J .LT. 1)J=1-J
      K=KPK1+KKK-2
C ADJUST VERTICAL WINDS to w/H, where H is the scale height
      AMBAR=1./((POO2(I,K,J)/32.)+(POO(I,K,J)/16.)+((1.-POO2(I,K,J)-
     1         POO(I,K,J))/28.))
      SHGT=8314.*T(I,K,J)/(AMBAR*8.7)
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C note
      W(I,K,J)=W(I,K,J)/SHGT
C  If ion trajectories are to be calculated use ion winds
 2039 CONTINUE
c  calculate weighting terms for nearest neighbour points:
      X1 = GLON(IPK2) - ALNGX
      Y1 = GLAT(JPK2) - ALATGX
      IF (MFLAG .GT.  0) THEN
        Y1 = 5. - GLAT(JPK1) + ALATGX
        X1 = 5. - ALNGX + GLON(IPK1)
        ENDIF
      IF(MFLAG .NE. 0)WRITE(6,5060)X1,Y1,ALNGX,ALATGX,MFLAG,IPK1,JPK1
 5060 FORMAT(1X,4F8.3,'MFLAG',I3,'IPK1 JPK1',2I3)
      S1 = SQRT( (5.0-X1)**2 + (5.0-Y1)**2  + 0.0000001)
      S2 = SQRT( (5.0-X1)**2 + Y1**2  + 0.000001)
      S3 = SQRT( X1**2 + (5.0-Y1)**2 +0.000001)
      S4 = SQRT( X1**2 + Y1**2+0.000001)
      S5 = SQRT( (5.0-X1)**2 + (5.0-Y1)**2 +0.000001)
      S6 = SQRT( (5.0-X1)**2 + Y1**2 + 0.000001)
      S7 = SQRT( X1**2 + (5.0-Y1)**2 + 0.000001)
      S8 = SQRT( X1**2 + Y1**2 + 0.000001)
      AAA = 1.0 / ( (1./S1) + (1./S2) + (1./S3) + (1./S4))
C  See earlier comment
      Z1 = FLOAT(KPK2)-ZKP1
 5065 FORMAT(1X,'SSSSS',9F7.3,I3)
      W1 = AAA/S1 *Z1
      W2 = AAA/S2 *Z1
      W3 = AAA/S3 *Z1
      W4 = AAA/S4 *Z1
      W5 = AAA/S5 *(1.-Z1)
      W6 = AAA/S6 *(1.-Z1)
      W7 = AAA/S7 *(1.-Z1)
      W8 = AAA/S8 *(1.-Z1)
 5070 FORMAT(1X,'WWS',4F8.3)
C Find the velocities acting on the parcel
      UPLEN = (W1*U(IPK1,KPK1,JPK1)+W2*
     1                U(IPK3,KPK1,JPK2) +
     2                W3*U(IPK2,KPK1,JPK1)+W4*
     3                U(IPK4,KPK1,JPK2) +
     4                W5*U(IPK1,KPK2,JPK1)+W6*
     5                U(IPK3,KPK2,JPK2) +
     6                W7*U(IPK2,KPK2,JPK1)+W8*
     7                U(IPK4,KPK2,JPK2))
      VPLEN = (W1*V(IPK1,KPK1,JPK1)+W2*
     1                V(IPK3,KPK1,JPK2) +
     2                W3*V(IPK2,KPK1,JPK1)+W4*
     3                V(IPK4,KPK1,JPK2) +
     4                W5*V(IPK1,KPK2,JPK1)+W6*
     5                V(IPK3,KPK2,JPK2) +
     6                W7*V(IPK2,KPK2,JPK1)+W8*
     7                V(IPK4,KPK2,JPK2))
      WPLEN = (W1*W(IPK1,KPK1,JPK1)+W2*W(IPK3,KPK1,JPK2)+
     1                W3*W(IPK2,KPK1,JPK1)+W4*W(IPK4,KPK1,JPK2)+
     2                W5*W(IPK1,KPK2,JPK1)+W6*W(IPK3,KPK2,JPK2)+
     3                W7*W(IPK2,KPK2,JPK1)+W8*W(IPK4,KPK2,JPK2))
C No vertical drifts for ion trajectories
      IF(JEOMGN .EQ. 1)THEN
      WPLEN = 0.
      ENDIF
      ALMONE = (ALNGX - AGLON) * DTR
C Adjust so that the newly calculated velocitiesare arranged in the
C zonal and meridional directions at the point of calculation.
      UPNEW = UPLEN * COS(ALMONE) - AHEM*VPLEN * SIN(ALMONE)
      VPNEW = AHEM*UPLEN * SIN(ALMONE) + VPLEN * COS(ALMONE)
      WPNEW = WPLEN
 6103 CONTINUE
      UPKT(JUT) = UPLEN
      VPKT(JUT) = VPLEN
      WPKT(JUT) = WPLEN
      ALAT(JUT) = ALATGX
      ALON(JUT) = ALNGX
      AGLAT = ALATGX
      AGLON = ALNGX
      ALONA(JUT) = AGLON
      ALATA(JUT) = AGLAT
      IF(JEOMGN .EQ. 1 .OR. KEOMGN .EQ. 1)THEN
        CALL GTM (ALAT(JUT)*DTOR,ALON(JUT)*DTOR,CHMAG1*DTOR,
     1 CHMAG2*DTOR,A1,A2,DIP1,DEC1,W,1)
        ALON(JUT) = A2/DTOR
        ALAT(JUT) = A1/DTOR
      ENDIF
C ANGULAR DISTANCE TRAVELLED:
      ZPKT(JUT)=ZKP1
      AZPKT(JUT)=ZPKT(JUT)
C
C FIND NEAREST NEIGHBOUR POINTS
C
      UTPKT(JUT)=UT
      UTPKT2(JUT) = UTPKT2(JUT-1) + AUSTP*TIMSTP/3600.
      IPKT=IPK1
      JPKT=((ALATGX) + 92.5)/5.
      KPKT=KPK1
C  DHOUR fixes local time differences between the geomagnetic and
C  geographic poles
      DHOUR=0.
      CHGHR = 4.667
      IF(JEOMGN .EQ. 1 .OR. KEOMGN .EQ. 1)THEN
      DHOUR = CHGHR
      IF (ALATS/ABS(ALATS) .LT. 0.) THEN
      CHGHR = -8.467
      DHOUR = CHGHR
      ENDIF
      ENDIF
      XLT(JUT) = UT + ALON(JUT)/15. - DHOUR
      IF(XLT(JUT) .GT. 24.)XLT(JUT)=XLT(JUT)-24.
      IF(XLT(JUT) .LT. 0. )XLT(JUT)=XLT(JUT)+24.
 5080 FORMAT(' DEBUG:',2F6.1,2F6.1,6F6.2)
C
C
C
 2036 CONTINUE
 2034 CONTINUE
      return
      end