#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 "fieldz.h"
      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)
#include "trgm.h"
!
! Local:
      real :: sin10,cos10
      integer :: i,j
C     ****
C     ****     FILL /TGCM/ AND /MAGFLD/
C     ****
      SIN10=1.E-6 ! snoe
!     SIN10=0.17  ! kibo
      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(ZZB(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) = -ZZB(I,J)/BMOD(I,J)

!         write(6,"('magdyn: i=',i3,' j=',i3,' zb(i,j)=',e12.4,
!    |      ' bmod(i,j)=',e12.4,' bz(i,j)=',e12.4)") 
!    |      i,j,zzb(i,j),bmod(i,j),bz(i,j)  

          BMOD2(I,J) = BMOD(I,J)
C         ****     SET MINIMUM DIP TO 10 DEGREES
          if (ABS(BZ(I,J))-SIN10 < 0.) then
            BX(I,J)=BX(I,J)*(COS10/SQRT(1.-BZ(I,J)**2))
            BY(I,J)=BY(I,J)*(COS10/SQRT(1.-BZ(I,J)**2))
            BZ(I,J)=SIGN(SIN10,BZ(I,J))
          endif
    1 CONTINUE

!     do j=1,zjmx
!       write(6,"('j=',i3,' bz(:,j)=',/,(6e12.4))") j,bz(:,j)
!     enddo

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

!     do j=1,zjmx
!       write(6,"('j=',i3,' bz(:,j)=',/,(6e12.4))") j,bz(:,j)
!     enddo

      RETURN
      END