C
      SUBROUTINE THREED
      implicit none
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"
      integer,parameter :: kmax=zkmx, kmaxp=kmax+1
      include "consts.h"
      include "dynphi.h"
!     real :: thetam,pslot,qslot
!     integer :: islot,jslot
!     COMMON THETAM(IMAXMP),ISLOT(IMAXMP),JSLOT(IMAXMP),
!    1  PSLOT(IMAXMP),QSLOT(IMAXMP)
!
! 4/00: make these local (s.a., vdrift.f):
      real :: thetam(imaxmp),pslot(imaxmp),qslot(imaxmp)
      integer :: islot(imaxmp),jslot(imaxmp)
!
      include "cterp.h"
!
! Local:
      integer :: j,k,i
      real :: cosltm,phims,phimn
      character(len=80) :: title

! External:
      real,external :: sddot ! in util.F
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           ****
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 = sddot(IMAXM,UNIT,PHIM(1,1))/FLOAT(IMAXM)
      PHIMN = sddot(IMAXM,UNIT,PHIM(1,JMAXM))/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

!     do j=1,jmaxm
!       do k=1,zkmxp,zkmxp-1
!         write(6,"('threed: j=',i3,' k=',i3,' phim3d(:,j,k)=',/,
!    |      (6e12.4))") j,k,phim3d(:,j,k)
!       enddo
!     enddo
!
! Save magnetic potential on secondary histories:
      do j=1,jmaxm
        call addfsech('PHIM3D','PHIM3D',' ',
     |    phim3d(:,j,:),imaxmp,kmaxp+3,kmaxp+3,j)
      enddo
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
      RETURN
      END