#include "dims.h"
      SUBROUTINE MAGDYN
      implicit none
C     ****
C     ****     SET UP VERSIONS OF /TRGM/ AND /MAGFLD/ CALCULATED
C     ****       FROM TRUE FIELD AVAILABLE IN /FIELD/
C     ****
#include "params.h"
      integer,PARAMETER :: ZIMXP2=ZIMX+2,ZJMXP2=ZJMX+2
#include "field.h"
!
! Note this /MAGFLD/ is the same memory as /MAGFLD/ in magfld.h,
! therefore order is important:
      real :: bx,by,bz,bmod2
      COMMON/MAGFLD/BX(-1:ZIMXP2,-1:ZJMXP2),BY(-1:ZIMXP2,-1:ZJMXP2),
     1  BZ(-1:ZIMXP2,-1:ZJMXP2),BMOD2(-1:ZIMXP2,-1:ZJMXP2)
      real :: rlatm,rlonm,dipmag,decmag,sndec,csdec,sn2dec,sncsdc,
     |  dumdum
      COMMON/TRGM/RLATM(ZIMXP,ZJMX),RLONM(ZIMXP,ZJMX),
     1  DIPMAG(ZIMXP,ZJMX),DECMAG(ZIMXP,ZJMX),SNDEC(ZIMXP,ZJMX),
     2  CSDEC(ZIMXP,ZJMX),SN2DEC(ZIMXP,ZJMX),SNCSDC(ZIMXP,ZJMX),
     3  DUMDUM(ZJMX,3)
C     real,parameter :: SIN10=1.E-6
C     real,parameter :: SIN10=0.01
      real,parameter :: SIN10=0.17
      real :: cos10,temp
      integer :: i,j
C     ****
C     ****     FILL /TRGM/ AND /MAGFLD/
C     ****
      COS10=SQRT(1.-SIN10**2)
      DO 1 J = 1,ZJMX
        DO 1 I = 1,ZIMX
          RLATM(I+2,J) = ALATM(I,J)
          RLONM(I+2,J) = ALONM(I,J)
          DIPMAG(I+2,J) = ATAN(ZB(I,J)/SQRT(XB(I,J)**2+YB(I,J)**2))
          DECMAG(I+2,J) = -ATAN2(YB(I,J),XB(I,J))
          SNDEC(I+2,J) = SIN(DECMAG(I+2,J))
          CSDEC(I+2,J) = COS(DECMAG(I+2,J))
          SN2DEC(I+2,J) = SNDEC(I+2,J)**2
          SNCSDC(I+2,J) = SNDEC(I+2,J)*CSDEC(I+2,J)
          BX(I,J) = YB(I,J)/BMOD(I,J)
          BY(I,J) = XB(I,J)/BMOD(I,J)
          BZ(I,J) = -ZB(I,J)/BMOD(I,J)
          BMOD2(I,J) = BMOD(I,J)
C         ****     SET MINIMUM DIP TO 10 DEGREES
!         BX(I,J)=BX(I,J)*merge(1.,COS10/SQRT(1.-BZ(I,J)**2),
!    1      ABS(BZ(I,J))-SIN10>=0.)
!         BY(I,J)=BY(I,J)*merge(1.,COS10/SQRT(1.-BZ(I,J)**2),
!    1      ABS(BZ(I,J)) -SIN10>=0.)
!         BZ(I,J)=merge(BZ(I,J),SIGN(SIN10,BZ(I,J)),
!    1      ABS(BZ(I,J))-SIN10>=0.)
          if (abs(bz(i,j)) < sin10) then
            temp = cos10/sqrt(1.-bz(i,j)**2)
            bx(i,j) = bx(i,j)*temp
            by(i,j) = by(i,j)*temp
            bz(i,j) = sign(sin10,bz(i,j))
          endif
    1 CONTINUE
C     ****
C     ****     VALUES AT J = -1, 0, ZJMAXP1, ZJNMP2
C     ****
      DO 2 J = 1,2
CDIR$ IVDEP
        DO 2 I = 1,ZIMX
          BX(I,J-2) = -BX(1+MOD(I-1+ZIMX/2,ZIMX),3-J)
          BY(I,J-2) = -BY(1+MOD(I-1+ZIMX/2,ZIMX),3-J)
          BZ(I,J-2) = BZ(1+MOD(I-1+ZIMX/2,ZIMX),3-J)
          BMOD2(I,J-2) = BMOD2(1+MOD(I-1+ZIMX/2,ZIMX),3-J)
          BX(I,ZJMX+J) = -BX(1+MOD(I-1+ZIMX/2,ZIMX),ZJMX+1-J)
          BY(I,ZJMX+J) = -BY(1+MOD(I-1+ZIMX/2,ZIMX),ZJMX+1-J)
          BZ(I,ZJMX+J) = BZ(1+MOD(I-1+ZIMX/2,ZIMX),ZJMX+1-J)
          BMOD2(I,ZJMX+J) = BMOD2(1+MOD(I-1+ZIMX/2,ZIMX),ZJMX+1-J)
    2 CONTINUE
C     ****
C     ****     PERIODIC POINTS
C     ****
      DO  3 I = 1,2
        DO 3  J = 1,8*ZJMX
          RLATM(I,J) = RLATM(I+ZIMX,J)
          RLATM(I+ZIMXP2,J) = RLATM(I+2,J)
    3 CONTINUE
      DO 4 I = 1,2
        DO 4 J = -1,4*(ZJMX+4)-2
          BX(I-2,J) = BX(I-2+ZIMX,J)
          BX(I+ZIMX,J) = BX(I,J)
    4 CONTINUE
      RETURN
      END