#include "dims.h"
      SUBROUTINE DYN
      implicit none
C     ****
C     ****     SETS UP COEFFICIENTS FOR DYNAMO PDE
C     ****
#include "params.h"
#include "coefm.h"
#include "dynamo.h"
#include "fieldz.h"
#include "cterp.h"
#include "trig.h"
#include "fcom.h"
#include "ceee.h"
#include "dynphi.h"

! Local:
      real :: ARRAY(-15:IMX0+16,JMX0),CS(JMX0)
      EQUIVALENCE (F(1,1),ARRAY(-15,1))
      real :: C0(IMX0,JMX0,10),C1(IMX1,JMX1,9),C2(IMX2,JMX2,9),
     1  C3(IMX3,JMX3,9),C4(IMX4,JMX4,9)
      EQUIVALENCE (CEE,C0),(CEE(NC1),C1),(CEE(NC2),C2),(CEE(NC3),C3),
     1  (CEE(NC4),C4)
      real :: WK(IMX0,3),pi,dlonm,dlatm,sym
      integer :: j,j0,jj,jjj,i,jntl
C     ****
C     ****     TRANSFORM NEEDED FIELDS TO GEOMAGNETIC COORDINATES.
C     ****       PERFORM FIELD-LINE INTEGRATIONS
C     ****       EVALUATE PDE COEFFICIENTS AND RHS
C     ****
      CALL TRANSF
C     ****
C     ****     SET ARRAY NC AND MAGNETIC LATITUDE COSINE ARRAY
C     ****
      NC(1) = NC0
      NC(2) = NC1
      NC(3) = NC2
      NC(4) = NC3
      NC(5) = NC4
      NC(6) = NCEE
      PI = 4.*ATAN(1.)
      DLONM = 2.*PI/IMAXM
      DLATM = PI/(JMAXM-1)
      DO 10 J = 1,JMX0
        CS(J) = COS(PI/2.-(JMX0-J)*DLATM)
   10 CONTINUE
C     ****
C     ****     SET UP DIFFERENCE COEFFICIENTS
C     ****
C     ****       REPLACE ZIGM11 BY A,  ZIGM22 BY B,  ZIGMC BY C
C     ****         AND ZIGM2 BY D
C     ****
      J0 = JMX0-JMAXMH
      DO 30 J = 1,JMAXMH
        JJ = JMAXMH+J-1
        JJJ = JMAXMH -J+1
        DO 31 I = 1,IMAXMP
          ZIGMC(I,JJ) = (ZIGMC(I,JJ)+ZIGM2(I,JJ))/(4.*DLATM*DLONM)
          ZIGM2(I,JJ) = ZIGMC(I,JJ)-2.*ZIGM2(I,JJ)/(4.*DLATM*DLONM)
          ZIGM22(I,JJ) = ZIGM22(I,JJ)*CS(J0+J)/DLATM**2
          ZIGMC(I,JJJ) = -ZIGMC(I,JJ)
          ZIGM2(I,JJJ) = -ZIGM2(I,JJ)
          ZIGM22(I,JJJ) = ZIGM22(I,JJ)
   31   CONTINUE
        IF(J.NE.JMAXMH)THEN
          DO 32 I = 1,IMAXMP
            ZIGM11(I,JJ) = ZIGM11(I,JJ)/(CS(J0+J)*DLONM**2)
            ZIGM11(I,JJJ) = ZIGM11(I,JJ)
   32     CONTINUE
        ENDIF
   30 CONTINUE
C     ****
C     ****     SET ZIGM11 TO ZERO AT MEGNETIC POLES TO AVOID FLOATING
C     ****       EXCEPTION  ( VALUES AT POLES ARE NOT USED )
C     ****
      DO 34 I = 1,IMAXMP
        ZIGM11(I,1) = 0.0
        ZIGM11(I,JMAXM) = 0.0
   34 CONTINUE
C     ****
C     ****     CLEAR ARRAY FOR DIFFERENCE STENCILS AT ALL LEVELS
C     ****
      CALL CLEAR(CEE,IMX0,JMX0)
C     ****
C     ****     CALCULATE CONTRIBUTION TO STENCILS FROM EACH PDE
C     ****       COEFFICIENT
C     ****
      ARRAY(1,1) = 1.
      SYM = 1.
      CALL STENCIL(ZIGM11,CS(1),IMX0,JMX0,SYM,CEE,ARRAY(1,1),1,WK)
      ARRAY(1,1) = 4.
      SYM = 1.
      CALL STENCIL(ZIGM22,CS(1),IMX0,JMX0,SYM,CEE,ARRAY(1,1),4,WK)
      ARRAY(1,1) = 2.
      SYM = -1.
      CALL STENCIL(ZIGMC,CS(1),IMX0,JMX0,SYM,CEE,ARRAY(1,1),2,WK)
      ARRAY(1,1) = 3.
      SYM = -1.
      CALL STENCIL(ZIGM2,CS(1),IMX0,JMX0,SYM,CEE,ARRAY(1,1),3,WK)
C     ****
C     ****     CALL RTHS TO INSERT RHS IN FINEST STENCIL
C     ****
      CALL RTHS(RHS,IMX0,JMX0,JMAXMH,C0)
C     ****
C     ****      SET BOUNDARY CONDITION AT EQUATOR AND POLE
C     ****
      CALL EDGES(C0,IMX0,JMX0)
      CALL EDGES(C1,IMX1,JMX1)
      CALL EDGES(C2,IMX2,JMX2)
      CALL EDGES(C3,IMX3,JMX3)
      CALL EDGES(C4,IMX4,JMX4)
C     ****
C     ****     DIVIDE STENCILS BY COS(THETA)
C     ****
C     ****
      CALL DIVIDE(C0,IMX0,JMX0,IMX0,JMX0,CS(1))
      CALL DIVIDE(C1,IMX1,JMX1,IMX0,JMX0,CS(1))
      CALL DIVIDE(C2,IMX2,JMX2,IMX0,JMX0,CS(1))
      CALL DIVIDE(C3,IMX3,JMX3,IMX0,JMX0,CS(1))
      CALL DIVIDE(C4,IMX4,JMX4,IMX0,JMX0,CS(1))
C     ****
C     ****     SET VALUE OF SOLUTION TO 1. AT POLE
C     ****
      DO 20 I = 1,IMX0
        C0(I,JMX0,10) = 1.
   20 CONTINUE
C     ****
C     ****     MODIFY STENCILS AND RHS SO THAT HEELIS POTENTIAL
C     ****       ID INSERTED AT HIGH LATITUDE
C     ****
      CALL ADHEEL
      JNTL = 0
C     ****
C     ****     CALL MUDPAC TO SOLVE PDE
C     ****       SOLUTION RETURNED IN RIM(1)
C     ****
!     write(6,"('dyn before mud')")
      CALL MUD(RIM,JNTL)
!     write(6,"('dyn after mud')")
C     ****
C     ****     COPY OUTPUT POTENTIAL FROM RIM TO PHIM(IMAXMP,JMAXM)
C     ****
      DO J = 1,JMAXM
        DO I = 1,IMAXMP
!         PHIM(I,J) = RIM(I,J,1)
          PHIM(I,J) = RIM(I,J,1) - PHIANT(I,J)
        enddo
      enddo
C     ****
C     ****     CALL THREED TO GENERATE 3-D POTENTIAL ARRAY IN
C     ****       GEOMAGNETIC COORDINATES.  RESULT RETURNED IN
C     ****       PHIM3D(IMAXMP,JMAXM,-2:KMAXP).  ALSO TRANSFORM THIS
C     ****       FIELD TO GEOGRAPHIC SPACE IN ARRAY
C     ****       DYNPOT(IMAXGP,0:JMAXGP,KMAXP)
C     ****
      CALL THREED
      RETURN
      END
C