C SUBROUTINE THREED C **** C **** CALCULATE 3-D DYNAMO POTENTIAL IN GEOMAGNETIC COORDINATES C **** C **** PHIM(IMAXMP,JMAXM) = 2-D POTENTIAL IN GEOMAGNETIC C **** COORDINATES OUTPUT BY DYNAMO CODE. C **** C **** PHIM3D(IMAXMP,JMAXM,KMAXP) = 3-D DYNAMO POTENTIAL IN C **** GEOMAGNETIC COORDINATES GENERATED BY THIS SUBROUTINE. C **** include "params.h" PARAMETER (IMAXG = ZIMX, JMAXG = ZJMX, KMAX = ZKMX) PARAMETER (IMAXGP = IMAXG+1, JMAXGP = JMAXG+1, 1 JMXGP2 = JMAXG+2, KMAXP = KMAX+1) include "consts.h" include "dynphi.h" COMMON THETAM(IMAXMP),ISLOT(IMAXMP),JSLOT(IMAXMP), 1 PSLOT(IMAXMP),QSLOT(IMAXMP) include "cterp.h" C **** C **** FOR EACH GEOMAGNETIC GRID POINT OVER THE 3-D GRID, C **** DETERMINE THE LOCATION OF THE FOOT OF THE MAGNETIC C **** FIELD PASSING THROUGH THE GRID POINT. ASSUME THAT C **** THE FIELD LINES ARE DIPOLE IN SHAPE. RECALLING THAT C **** THE MODEL ASSUMES CONSTANT ELECTRIC POTENTIAL ALONG C **** EACH FIELD LINE, INTERPOLATION IN THE 2-D ARRAY PHIM C **** GIVES THE POTENTIAL AT THE ORIGINAL GRID POINT. C **** C **** LOOP OVER GEOMAGNETIC GRID. C **** DO 1 K = -2,KMAXP DO 1 J = 2,JMAXM-1 C **** COSLTM = COS(THETAS) COSLTM = COS(YLATM(J)) DO 1 I = 1,IMAXM C **** THETAM(I) = THETAS FOR EACH VALUE OF I THETAM(I) = ACOS(SQRT(R0/(R00+ZZM(I,J,K)))*COSLTM) C **** C **** FIND POSITION OF THETAM IN TABLE OF C **** THETA0 VS. THETAS C **** PSLOT(I) = THETAM(I)*180./PI+1. ISLOT(I) = PSLOT(I) PSLOT(I) = PSLOT(I)-FLOAT(ISLOT(I)) C **** C **** DETERMINE THATA0 BY INTERPOLATION IN TABLE C **** THETAM(I) = (1.-PSLOT(I))*TABLE(ISLOT(I),2)+ PSLOT(I)* 1 TABLE(ISLOT(I)+1,2) C ****put the right sign for THETAM (G. Lu Oct. 1997)*** C the0 = THETAM(I) THETAM(I) = SIGN(THETAM(I),YLATM(J)) C if (K .eq. KMAXP .AND. I .eq. 1) C 1 WRITE(6,"('THREED: THETAM,YLATM,the0 -- ',3f7.1)") C 2 THETAM(I)*180./PI,YLATM(J)*180./PI,the0*180./PI C **** C **** NOW IDENTIFY GEOMAGNETIC GRID ELEMENT IN WHICH C **** THE POINT LIES IN PHIM C **** ISLOT(I) = I PSLOT(I) = 0. QSLOT(I) = (THETAM(I)+PI/2.)/DLATM+1. JSLOT(I) = QSLOT(I) QSLOT(I) = QSLOT(I)-FLOAT(JSLOT(I)) C **** C **** INTERPOLATE IN PHIM USING BILINEAR INTERPOLATION C **** PHIM3D(I,J,K) = (1.-QSLOT(I))*((1.-PSLOT(I))* 1 PHIM(ISLOT(I),JSLOT(I))+PSLOT(I)* 2 PHIM(ISLOT(I)+1,JSLOT(I)))+QSLOT(I)*((1.-PSLOT(I))* 3 PHIM(ISLOT(I),JSLOT(I)+1)+PSLOT(I)* 4 PHIM(ISLOT(I)+1,JSLOT(I)+1)) 1 CONTINUE C **** C **** CALCULATE VALUE OF PHIM3D AT GEOMAGNETIC POLES C **** PHIMS = SDOT(IMAXM,UNIT,1,PHIM(1,1),1)/FLOAT(IMAXM) PHIMN = SDOT(IMAXM,UNIT,1,PHIM(1,JMAXM),1)/FLOAT(IMAXM) DO 4 K = -2,KMAXP DO 4 I = 1,IMAXM PHIM3D(I,1,K) = PHIMS PHIM3D(I,JMAXM,K) = PHIMN 4 CONTINUE C **** C **** PERIODIC POINT C **** DO 2 K = 1,KMAXP DO 2 J = 1,JMAXM PHIM3D(IMAXMP,J,K) = PHIM3D(1,J,K) 2 CONTINUE C **** C **** TRANSFORM PHIM3D TO GEOGRAPHIC COORDINATES IN C **** DYNPOT(IMAXGP,0:JMAXGP,KMAXP) C **** DO 3 K = 1,KMAXP CALL GRDTRP(PHIM3D(1,1,K),DYNPOT(1,0,K),IM(1,0),JM(1,0), 1 DIM(1,0),DJM(1,0),IMAXGP,IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM) 3 CONTINUE C **** C **** PERIODIC POINTS C **** DO 5 K = 1,KMAXP DO 5 J = 0,JMAXG+1 DYNPOT(IMAXGP,J,K) = DYNPOT(1,J,K) 5 CONTINUE C **** C **** TRANSFORM PHIHM.TO GEOGRAPHIC COORDINATES IN C **** PHIH(IMAXGP,0:JMAXGP) C **** CALL GRDTRP(PHIHM(1,1),PHIH(1,0),IM(1,0),JM(1,0), 1 DIM(1,0),DJM(1,0),IMAXGP,IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM) C **** C **** PERIODIC POINTS C **** DO 6 J = 0,JMAXG+1 PHIH(IMAXGP,J) = PHIH(1,J) 6 CONTINUE C CALL EZCNTR(PHIM,IMAXMP,JMAXM) C CALL EZCNTR(PHIM3D(1,1,KMAXP),IMAXMP,JMAXM) C CALL EZCNTR(DYNPOT(1,0,KMAXP),IMAXGP,JMAXGP/2+1) C CALL EZCNTR(PHIH,IMAXM,JMAXM) C CALL MKCON(PHIM3D(1:imaxmp,1:jmaxm,KMAXP),IMAXMP,JMAXM,0.,0, C 1 'MLONG','MLAT','THREED: PHIM3D') C CALL MKCON(DYNPOT(1:imaxgp,0:jmaxgp,KMAXP),IMAXGP,JMAXGP+1, C 1 0.,0,'GLONG','GLAT','THREED: DYNPOT') RETURN END