      SUBROUTINE DYN
C     ****
C     ****     SETS UP COEFFICIENTS FOR DYNAMO PDE
C     ****
      include "params.h"
      PARAMETER (IMAXG=ZIMX, JMAXG=ZJMX)
      PARAMETER (IMAXGP=IMAXG+1, JMAXGP=JMAXG+1)
      include "coefm.h"
      include "consts.h"
      include "dynamo.h"
      include "field.h"
      include "cterp.h"
      include "trig.h"
      include "blnk.h"
      include "ceee.h"
      include "dynphi.h"
      DIMENSION ARRAY(-15:IMX0+16,JMX0),CS(JMX0)
      EQUIVALENCE (F(1),ARRAY(-15,1))
      DIMENSION 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)
      DIMENSION WK(IMX0,3)
      dimension c0sv(imx0,jmx0,10)
C  FACE = FAC from electric fields only
C  FACVN = FAC from neutral winds only
C  FACT = total FAC
      common /facvne/face(imaxmp,jmaxmh),facvn(imaxmp,jmaxmh),
     |  fact(imaxmp,jmaxmh),fac3(imaxmp,jmaxmh,3)
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
	do 1100 i=1,imx0
	do 1100 k=1,jmx0
	facvn(i,k) = -50.*dt0dts(jmx0-1+k)*rhs(i,k)/(r0*r0)
	do 1100 m=1,10
	c0sv(i,k,m) = 50.*dt0dts(jmx0-1+k)*c0(i,k,m)/(r0*r0)
1100    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     ****
      CALL MUD(RIM,JNTL)
	do 1300 i=1,imx0
	face(i,1) = 0.
	fact(i,1) = facvn(i,1)
	im1 = i - 1
	if (i .eq. 1) im1 = imx0 - 1
	ip1 = i + 1
	if (i .eq. imx0) ip1 = 2
	sum90 = 0.
	do 1200 k=2,jmaxmh-1
	k0 = jmx0 + k - 1
	kp1 = jmx0 + k
	km1 = jmx0 + k - 2
	face(i,k) = c0sv(i,k,1)*rim(ip1,k0)+c0sv(i,k,2)*rim(ip1,kp1)
     |   + c0sv(i,k,3)*rim(i,kp1) + c0sv(i,k,4)*rim(im1,kp1)
     |   + c0sv(i,k,5)*rim(im1,k0) + c0sv(i,k,6)*rim(im1,km1)
     |   + c0sv(i,k,7)*rim(i,km1) + c0sv(i,k,8)*rim(ip1,km1)
     |   + c0sv(i,k,9)*rim(i,k0)
	fact(i,k) = face(i,k) + facvn(i,k)
1200    continue
	sum90 = sum90 + face(i,jmaxmh-1)
1300    continue
	do 1400 i=1,imx0
	fact(i,jmaxmh) = sum90 + facvn(i,jmaxmh)
1400    face(i,jmaxmh) = sum90
c     CALL EZCNTR (FACVN,IMX0,JMAXMH)
c     CALL EZCNTR (FACE,IMX0,JMAXMH)
c     CALL EZCNTR (FACT,IMX0,JMAXMH)
	do 1500 i=1,imaxmp
	do 1500 k=1,jmaxmh
	fac3(i,k,1) = facvn(i,jmaxmh+1-k)
	fac3(i,k,2) = face(i,jmaxmh+1-k)
	fac3(i,k,3) = fact(i,jmaxmh+1-k)
1500    continue
c     CALL EZCNTR (FAC3(1,1,1),IMX0,13)
c     CALL EZCNTR (FAC3(1,1,2),IMX0,13)
c     CALL EZCNTR (FAC3(1,1,3),IMX0,20)
c        write (6,"(1x,'facEvn(80,49) (A/m2) in mlon and mlat from ',
c    |    '90 to 45 for Peymirat:')")
c        write (6,*) ((fac3(i,k,3),i=1,imaxm),k=1,20)

C     ****
C     ****     COPY OUTPUT POTENTIAL FROM RIM TO PHIM(IMAXMP,JMAXM)
C     ****
      DO 33 J = 1,JMAXM
        DO 33 I = 1,IMAXMP
	  PHIM(I,J) = RIM(I,J) - PHIANT(I,J)
   33 CONTINUE
C  Plot of dynamo + hi-lat potential
C     CALL EZCNTR (PHIM,IMX0,JMAXM)
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
