C ******************************************************************* C ******************************************************************* subroutine traj(alons,alats,utstrt) C ******************************************************************* C ******************************************************************* C C set parameters to match size of input files parameter (izut=121,imx=73,kmx=25,jmx=36) parameter (nmx=5,imxint=73*nmx,jmxint=36*nmx) parameter (imx2=(imx-1)*nmx+1) parameter (jmx2=(jmx-1)*nmx+1,jmx3=jmx2+2) C components of the wind real u(imx2,jmxint),v(imx2,jmxint) C parcel local time real xlt(izut) real xltg(izut) C packet velocities real upkt(izut),vpkt(izut) C packet latitude and longitude real alat(izut),alon(izut) real alata(izut),alona(izut) C kuhlp? packet zonal coordinate integer kuhlp2,ipkt,ipku C packet meridional coordinate integer jpkt,jpku C packet vertical coordinate integer kpkt,kpku real utpkt(izut) real utpkt2(izut) real uvelin(imx2,jmx3,izut) real vvelin(imx2,jmx3,izut) common /velsin/uvelin,vvelin 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 common /trajout2/xltg C Also save the winds/ion drifts? along the trajectory path C common /trajout2/upkt,vpkt,wpkt C longitude grid real glon(361) C latitudfe grid real glat(36*5) real w(4) real pmlong(2) real dipole(2) real mag(4) DATA MAG/-74.5,127.,79.,-70./ PI=180./RTOD DIPOLE(1)=90.+MAG(1) DIPOLE(2)=90.-MAG(3) DTR=ATAN(1.)/45. CHMAG1 = MAG(3) CHMAG2 = MAG(4) IF (ALATS/ABS(ALATS) .LT. 0.)THEN CHMAG1 = MAG(1) CHMAG2 = MAG(2) ENDIF PMLONG(1)=CHMAG2-180 PMLONG(2)=CHMAG2 ibck=1 DTOR=ATAN(1.)/45. RTD=1./DTR pi = 180./rtd jeomgn = 1 astp = 1. timstp=60. re = 6371.E+5 do i = 1,361,1 glon(i) = -180. + (i-1) enddo do j = 1,180,1 glat(j) = -89.5 + (j-1) enddo write(*,*)'alons alats utstrt',alons,alats,utstrt C Time loop do 2034 jut=1,121,1 do 1033 i = 1,imx2,1 do 1033 j = 1,jmx3,1 u(i,j+2) = uvelin(i,j,jut) v(i,j+2) = vvelin(i,j,jut) 1033 continue C if(jut .eq. 1)write(35,*)u ut = utstrt-jut C C+++++++++++++++++++++++++++++++++++++++++++++ C CALCULATE TRAJECTORY LOCI C+++++++++++++++++++++++++++++++++++++++++++++ kuhlp2=0 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 FIND NEAREST NEIGHBOUR POINTS using the initial conditions for lat lon IPK1 = (ALONS+180.) + 1 IPK1 = MOD(IPK1+360,360) IF(IPK1 .EQ. 0)IPK1=IPK1 + 360 IPKHOLD = IPK1 IPK2 = MOD(IPK1+1,360) IF(IPK2 .EQ. 0)IPK2 = IPK2 + 360 JPK1 = (ALATS+90.5) JPKHOLD = JPK1 JPK2 = JPK1+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+181,360) IF(IPK3 .EQ. 0)IPK3=IPK3 + 360 IPK4 = MOD(IPK1+180,360) IF(IPK4 .EQ. 0)IPK4 = IPK4+360 ENDIF IF(JPK2 .GT. 180) THEN JPK2 = 180 MFLAG = 2 IPK3 = MOD(IPK1+181,360) IF(IPK3 .EQ. 0)IPK3=IPK3 + 360 IPK4 = MOD(IPK1+180,360) IF(IPK4 .EQ. 0)IPK4 = IPK4+360 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. 360)I=I-360 IF(I .LT. 1)I=I+360 J=JPKHOLD + JJJ - 3 AJPT = 1. IF(J .GT. 180)THEN I = MOD(IPK1+181,360) + 3 - III AJPT = -1. ENDIF IF(J .LT. 1)THEN I = MOD(IPK1+181,360) + 3 - III AJPT = -1. ENDIF IF(I .GT. 360)I=I-360 IF(I .LT. 1)I=I+360 IF(J .GT. 180)J=361-J IF(J .LT. 1)J=1-J u(i,j) = ajpt*u(i,j) v(i,j) = ajpt*v(i,j) 2027 CONTINUE c calculate weighting terms for nearest neighbour points: C write(*,*)ipk2,glon(ipk2),alons C write(*,*)jpk2,glat(jpk2),alats X1 = GLON(IPK2) - ALONS Y1 = GLAT(JPK2) - ALATS write(*,*)'glat alats',glat(jpk2),alats IF (MFLAG .GT. 0) THEN Y1 = 1. - GLAT(JPK1) + ALATS X1 = 1. - ALONS + GLON(IPK1) ENDIF WRITE(6,5040)X1,Y1 5040 FORMAT('X1 X2',2F6.2) S1 = SQRT( (1.0-X1)**2 + (1.0-Y1)**2 + 0.0000001) S2 = SQRT( (1.0-X1)**2 + Y1**2 + 0.000001) S3 = SQRT( X1**2 + (1.0-Y1)**2 +0.000001) S4 = SQRT( X1**2 + Y1**2+0.000001) C write(*,*)s1,s2,s3,s4,aaa AAA = 1.0 / ( (1./S1) + (1./S2) + (1./S3) + (1./S4)) C write(*,*)s1,s2,s3,s4,aaa C Get vertical information W1 = AAA/S1 W2 = AAA/S2 W3 = AAA/S3 W4 = AAA/S4 write(*,*)'wts',w1,w2,w3,w4,ipk1,ipk2,ipk3,ipk4, 1 jpk1,jpk2,jut,u(ipk1,jpk1),v(ipk1,jpk1) UPKT(JUT) = (W1*U(IPK1,JPK1)+W2* 1 U(IPK3,JPK2) + 2 W3*U(IPK2,JPK1)+W4* 3 U(IPK4,JPK2)) VPKT(JUT) = (W1*V(IPK1,JPK1)+W2* 1 V(IPK3,JPK2) + 2 W3*V(IPK2,JPK1)+W4* 3 V(IPK4,JPK2)) C write(*,*)'us',upkt(jut),vpkt(jut),ipk1,jpk1,u(ipk1,jpk1) IPKT=IPK1 IPKU=IPK1 JPKT=(ALATS +90.5) JPKU=JPK1 C 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 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/100. + ALON(JUT)/15. - DHOUR IF(XLT(JUT) .GT. 24.)XLT(JUT)=XLT(JUT)-24. IF(XLT(JUT) .LT. 0. )XLT(JUT)=XLT(JUT)+24. XLTG(JUT) = UT/100. + ALONA(JUT)/15. IF(XLTG(JUT) .GT. 24.)XLTG(JUT)=XLTG(JUT)-24. IF(XLTG(JUT) .LT. 0. )XLTG(JUT)=XLTG(JUT)+24. IHEMI=1 IF(ALAT(JUT) .LT. 0.)IHEMI=-1 C write(*,*)'xlt',xlt(jut) write(*,*)'alons alats utstrt',alons,alats,utstrt GO TO 2034 2024 CONTINUE C 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) UPNEW=UPOLD VPNEW=VPOLD IHEM=AGLAT/ABS(AGLAT) AHEM=FLOAT(IHEM) DO 6103 JNJK=1,2,1 ALNGX=AGLON ALATGX=AGLAT UPLEN = (UPOLD + UPNEW)/2. VPLEN = (VPOLD + VPNEW)/2. ATMSTP = TIMSTP C write(*,*)'timstp',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 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.) + 1 IPK1 = MOD(IPK1+360,360) IPKHOLD = IPK1 IF(IPK1 .EQ. 0)IPK1=IPK1 + 360 IPK2 = MOD(IPK1+1,360) IF(IPK2 .EQ. 0)IPK2 = IPK2 + 360 JPK1 = ((ALATGX)+90.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+181,360) IF(IPK3 .EQ. 0)IPK3=IPK3 + 360 IPK4 = MOD(IPK1+180,360) IF(IPK4 .EQ. 0)IPK4 = IPK4+360 ENDIF IF(JPK2 .GT. 180) THEN JPK2 = 180 MFLAG = 2 IPK3 = MOD(IPK1+181,360) IF(IPK3 .EQ. 0)IPK3=IPK3 + 360 IPK4 = MOD(IPK1+180,360) IF(IPK4 .EQ. 0)IPK4 = IPK4+360 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. 360)I=I-360 IF(I .LT. 1)I=I+360 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. 180)THEN I = MOD(IPK1+181,360) + 3 - III AJPT = -1. ENDIF IF(J .LT. 1)THEN I = MOD(IPK1+181,360) + 3 - III AJPT = -1. ENDIF IF(I .GT. 360)I=I-360 IF(I .LT. 1)I=I+360 IF(J .GT. 180)J=361-J IF(J .LT. 1)J=1-J u(i,j) = ajpt*u(i,j) v(i,j) = ajpt*v(i,j) 2039 CONTINUE c calculate weighting terms for nearest neighbour points: X1 = GLON(IPK2) - ALNGX Y1 = GLAT(JPK2) - ALATGX IF (MFLAG .GT. 0) THEN Y1 = 1. - GLAT(JPK1) + ALATGX X1 = 1. - ALNGX + GLON(IPK1) ENDIF C write(*,*)'x1 y1',x1,y1 S1 = SQRT( (1.0-X1)**2 + (1.0-Y1)**2 + 0.0000001) S2 = SQRT( (1.0-X1)**2 + Y1**2 + 0.000001) S3 = SQRT( X1**2 + (1.0-Y1)**2 +0.000001) S4 = SQRT( X1**2 + Y1**2+0.000001) AAA = 1.0 / ( (1./S1) + (1./S2) + (1./S3) + (1./S4)) 5065 FORMAT(1X,'SSSSS',9F7.3,I3) W1 = AAA/S1 W2 = AAA/S2 W3 = AAA/S3 W4 = AAA/S4 5070 FORMAT(1X,'WWS',4F8.3) C Find the velocities acting on the parcel UPLEN = (W1*U(IPK1,JPK1)+W2* 1 U(IPK3,JPK2) + 2 W3*U(IPK2,JPK1)+W4* 3 U(IPK4,JPK2)) VPLEN = (W1*V(IPK1,JPK1)+W2* 1 V(IPK3,JPK2) + 2 W3*V(IPK2,JPK1)+W4* 3 V(IPK4,JPK2)) 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) 6103 CONTINUE UPKT(JUT) = UPLEN VPKT(JUT) = VPLEN 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: C C FIND NEAREST NEIGHBOUR POINTS C UTPKT(JUT)=UT UTPKT2(JUT) = UTPKT2(JUT-1) + AUSTP*TIMSTP/3600. IPKT=IPK1 JPKT=((ALATGX) + 90.5) 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/100. + ALON(JUT)/15. - DHOUR IF(XLT(JUT) .GT. 24.)XLT(JUT)=XLT(JUT)-24. IF(XLT(JUT) .LT. 0. )XLT(JUT)=XLT(JUT)+24. XLTG(JUT) = UT/100. + ALONA(JUT)/15. IF(XLTG(JUT) .GT. 24.)XLTG(JUT)=XLTG(JUT)-24. IF(XLTG(JUT) .LT. 0. )XLTG(JUT)=XLTG(JUT)+24. C C C 2036 CONTINUE 2034 CONTINUE return end