#include "dims.h"
C
      SUBROUTINE THREED
      use cons_module,only: kmax,kmaxp1,pi
      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,kmaxp1) = 3-D DYNAMO POTENTIAL IN
C     ****       GEOMAGNETIC COORDINATES GENERATED BY THIS SUBROUTINE.
C     ****
#include "params.h"
#include "consdyn.h"
#include "dynphi.h"
      real :: thetam,pslot,qslot
      integer :: islot,jslot
      COMMON THETAM(IMAXMP),ISLOT(IMAXMP),JSLOT(IMAXMP),
     1  PSLOT(IMAXMP),QSLOT(IMAXMP)
#include "cterp.h"
!
! Local:
      real :: EEE=1.E-10
      integer :: i,j,k
      real :: cosltm,phims,phimn
      real,external :: sddot
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 = ZKBM,kmaxp1
        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) = (RE+ZZM(I,J,ZKBM))/(RE+ZZM(I,J,K))
            THETAM(I) = ACOS(SQRT(THETAM(I))*COSLTM*(1.-EEE))
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 = ZKBM,kmaxp1
        DO 4 I = 1,IMAXM
          PHIM3D(I,1,K) = PHIMS
          PHIM3D(I,JMAXM,K) = PHIMN
    4 CONTINUE
C     ****
C     ****     PERIODIC POINT
C     ****
      DO K = ZKBM,kmaxp1
        DO J = 1,JMAXM
          PHIM3D(IMAXMP,J,K) = PHIM3D(1,J,K)
        ENDDO
C       CALL EZCNTR(PHIM3D(1,1,K),IMAXMP,JMAXM)
      ENDDO
C     ****
C     ****     EXTEND PHIM3D FROM LEVEL ZKBM TO LEVEL 1 KEEPING
C     ****     POTENTIAL CONSTANT WITH HEIGHT.
C     ****
      DO K = 1,ZKBM-1
        DO I = 1,IMAXMP
          DO J = 1,JMAXM
            PHIM3D(I,J,K) = PHIM3D(I,J,ZKBM)
          ENDDO
        ENDDO
      ENDDO
C     DO I = 1,IMAXMP
C       DO K = 1,kmaxp1
C         DO J = 1,JMAXM
C           PPLLOOTT(J,K) = PHIM3D(I,J,K)
C         ENDDO
C       ENDDO
C       CALL EZCNTR(PPLLOOTT,JMAXM,kmaxp1)
C     ENDDO
C     ****
C     ****     TRANSFORM PHIM3D TO GEOGRAPHIC COORDINATES IN
C     ****       DYNPOT(IMAXGP,0:JMAXGP,kmaxp1)
C     ****
      DO 3 K = 1,kmaxp1
        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,kmaxp1
        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