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.)/45. RTD=1./DTR pi = 180./rtd jeomgn = 1 astp = 1. timstp=240. re = 6371.E+5 do i = 1,73,1 glon(i) = -180. + (i-1)*5. enddo do j = 1,36,1 glat(j) = -87.5 + (j-1)*5. enddo 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.)/5. + 1 IPK1 = MOD(IPK1+72,72) IF(IPK1 .EQ. 0)IPK1=IPK1 + 72 IPKHOLD = IPK1 IPK2 = MOD(IPK1+1,72) IF(IPK2 .EQ. 0)IPK2 = IPK2 + 72 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,72) IF(IPK3 .EQ. 0)IPK3=IPK3 + 72 IPK4 = MOD(IPK1+36,72) IF(IPK4 .EQ. 0)IPK4 = IPK4+72 ENDIF IF(JPK2 .GT. 36) THEN JPK2 = 36 MFLAG = 2 IPK3 = MOD(IPK1+37,72) IF(IPK3 .EQ. 0)IPK3=IPK3 + 72 IPK4 = MOD(IPK1+36,72) IF(IPK4 .EQ. 0)IPK4 = IPK4+72 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. 72)I=I-72 IF(I .LT. 1)I=I+72 J=JPKHOLD + JJJ - 3 AJPT = 1. IF(J .GT. 36)THEN I = MOD(IPK1+37,72) + 3 - III AJPT = -1. ENDIF IF(J .LT. 1)THEN I = MOD(IPK1+37,72) + 3 - III AJPT = -1. ENDIF IF(I .GT. 72)I=I-72 IF(I .LT. 1)I=I+72 IF(J .GT. 36)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* 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)) VPKT(JUT) = (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)) WPKT(JUT) = (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)) IF(JEOMGN .EQ. 1)THEN WPKT(JUT) = 0. ENDIF IPKT=IPK1 IPKU=IPK1 JPKT=(ALATS +92.5)/5. JPKU=JPK1 KPKT=KPK1 KPKU=KPK1 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 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+72,72) IPKHOLD = IPK1 IF(IPK1 .EQ. 0)IPK1=IPK1 + 72 IPK2 = MOD(IPK1+1,72) IF(IPK2 .EQ. 0)IPK2 = IPK2 + 72 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,72) IF(IPK3 .EQ. 0)IPK3=IPK3 + 72 IPK4 = MOD(IPK1+36,72) IF(IPK4 .EQ. 0)IPK4 = IPK4+72 ENDIF IF(JPK2 .GT. 36) THEN JPK2 = 36 MFLAG = 2 IPK3 = MOD(IPK1+37,72) IF(IPK3 .EQ. 0)IPK3=IPK3 + 72 IPK4 = MOD(IPK1+36,72) IF(IPK4 .EQ. 0)IPK4 = IPK4+72 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. 72)I=I-72 IF(I .LT. 1)I=I+72 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. 36)THEN I = MOD(IPK1+37,72) + 3 - III AJPT = -1. ENDIF IF(J .LT. 1)THEN I = MOD(IPK1+37,72) + 3 - III AJPT = -1. ENDIF IF(I .GT. 72)I=I-72 IF(I .LT. 1)I=I+72 IF(J .GT. 36)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