*DK ADDHEEL                                                               ADHEEL.1
      SUBROUTINE ADHEEL                                                   ADHEEL.2
C     ****                                                                ADHEEL.3
C     ****     MODIFY STENCILS AND RHS SO AS TO INSERT HEELIS             ADHEEL.4
C     ****       POTENTIAL FOR GEOMAGNETIC LATITUDE .GE. 60.0 DEG.        ADHEEL.5
C     ****                                                                ADHEEL.6
      include "param.h"                                                   ADHEEL.7
      INCLUDE 'epara.h'                                                   ADHEEL.8
      include "const.h"                                                   ADHEEL.9
      include "eceee.h"                                                   ADHEEL.10
      include "eterp.h"                                                   ADHEEL.11
      include "edphi.h"                                                   ADHEEL.12
      include "econs.h"                                                   ADHEEL.13
      COMMON /IONCR/ ISTAR,IHEM,THETA0(2),OFFA(2),OFFC(2),DSKOFA(2),      ADHEEL.14
     1  DSKOFC(2),PHID(2),PHIN(2),PHIDP0(2),PHIDM0(2),PHINP0(2),          ADHEEL.15
     2  PHINM0(2),PSIM(2),PSIE(2),PCEN(2),RR1(2),COLATC(IMAXMP,JMX0),     ADHEEL.16
     3  P(IMX0,JMX0),CRIT(2)                                              ADHEEL.17
      CRIT(1)=0.261799387
      CRIT(2)=0.523598775                                                  ADHEEL.18
C     ****                                                                ADHEEL.19
C     ****     FILL ARRAY COLATC WITH NORTHERN AURORAL COLATITUDES        ADHEEL.20
C     ****       CORRESPONDING TO EACH NORTHERN GEOMAGNETIC GRID          ADHEEL.21
C     ****       POINT.  VALUES IN COLATC(IMX0,JMX0).                     ADHEEL.22
C     ****       ALSO CALCULATE FRACTIONAL PRESENCE OF DYNAMO             ADHEEL.23
C     ****       EQUATION GIVEN CRITICAL COLATITUDES CRIT(2).             ADHEEL.24
C     ****       VALUES IN P(IMX0,JMX0).                                  ADHEEL.25
C     ****                                                                ADHEEL.26
      CALL COLATH                                                         ADHEEL.27
C     ****                                                                ADHEEL.28
C     ****     CALCULATE PHIHM, THE HEELIS POTENTIAL IN                   ADHEEL.29
C     ****       GEOMAGNETIC COORDINATES.                                 ADHEEL.30
C     ****                                                                ADHEEL.31
      ISTAR=0                                                             ADHEEL.32
      CALL POTM                                                           ADHEEL.33
C     ****                                                                ADHEEL.34
C     ****     MODIFY STENCILS                                            ADHEEL.35
C     ****                                                                ADHEEL.36
      NCC=1                                                               ADHEEL.37
      IMX=IMX0                                                            ADHEEL.38
      JMX=JMX0                                                            ADHEEL.39
      DO N=1,5                                                            ADHEEL.40
       CALL STENMOD(IMX,JMX,CEE(NCC),PHIHM(1,JMX0),P,COLATC)              ADHEEL.41
       NCC=NCC+9*IMX*JMX                                                  ADHEEL.42
       IF(N.EQ.1)NCC=NCC+IMX*JMX                                          ADHEEL.43
       IMX=(IMX+1)/2                                                      ADHEEL.44
       JMX=(JMX+1)/2                                                      ADHEEL.45
      ENDDO                                                               ADHEEL.46
      RETURN                                                              ADHEEL.47
      END                                                                 ADHEEL.48
C                                                                         ADHEEL.49
*DK CLEAR                                                                  CLEAR.1
      SUBROUTINE CLEAR(CEE,IMX0,JMX0)                                      CLEAR.2
      IMPLICIT REAL*8 (A-H,O-Z)                                            CLEAR.3
C     ****                                                                 CLEAR.4
C     ****     CLEAR C ARRAYS FOR STENCIL COEFFICIENTS                     CLEAR.5
C     ****                                                                 CLEAR.6
C     ****     CEE WILL CONTAIN                                            CLEAR.7
C     ****       C0(IMX0,JMX0,10), C1(IMX1,JMX1,9), C2(IMX2,JMX2,9),       CLEAR.8
C     ****       C3(IMX3,JMX3,9), C4(IMX4,JMX4,9)                          CLEAR.9
C     ****                                                                 CLEAR.10
C     ****       WHERE:                                                    CLEAR.11
C     ****         IMX1 = (IMX0+1)/2                                       CLEAR.12
C     ****         JMX1 = (JMX0+1)/2                                       CLEAR.13
C     ****         AND SO ON                                               CLEAR.14
C     ****                                                                 CLEAR.15
      DIMENSION CEE(1)                                                     CLEAR.16
C     ****                                                                 CLEAR.17
C     ****     COMPUTE TOTAL SIZE OF CEE, N                                CLEAR.18
C     ****                                                                 CLEAR.19
      IMX=IMX0                                                             CLEAR.20
      JMX=JMX0                                                             CLEAR.21
      N=0                                                                  CLEAR.22
      DO M=1,5                                                             CLEAR.23
       N=N+IMX*JMX                                                         CLEAR.24
       IMX=(IMX+1)/2                                                       CLEAR.25
       JMX=(JMX+1)/2                                                       CLEAR.26
      ENDDO                                                                CLEAR.27
      N=9*N+IMX0*JMX0                                                      CLEAR.28
C     ****                                                                 CLEAR.29
C     ****     CLEAR ARRAY CEE                                             CLEAR.30
C     ****                                                                 CLEAR.31
      DO I=1,N                                                             CLEAR.32
       CEE(I)=0.                                                           CLEAR.33
      ENDDO                                                                CLEAR.34
      RETURN                                                               CLEAR.35
      END                                                                  CLEAR.36
C                                                                          CLEAR.37
*DK CNM                                                                      CNM.1
      SUBROUTINE CNM(ARRAY,CS,IMX0,JMX0,IMX,JMX,C,NCOEF,WK)                  CNM.2
      IMPLICIT REAL*8 (A-H,O-Z)                                              CNM.3
C     ****                                                                   CNM.4
C     ****     COMPUTE CONTRIBUTION TO STENCIL FROM ZIGM(NCOEF) ON GRID      CNM.5
C     ****       IMX BY JMX.  FINEST GRID IS IMX0 BY JMX0.                   CNM.6
C     ****                                                                   CNM.7
C     ****     INPUT:                                                        CNM.8
C     ****       ARRAY(-15:IMX0+16,JMX0)                                     CNM.9
C     ****         ARRAY OF ZIGM ON GRID IMX0 BY JMX0 EXTENDED USING         CNM.10
C     ****         SYMMETRY AND PERIODICITY                                  CNM.11
C     ****       CS(JMX0)                                                    CNM.12
C     ****         ARRAY OF COSINES OF LATITUDE VALUES                       CNM.13
C     ****       NCOEF = NUMBER IDENTIFYING COEFFICIENT                      CNM.14
C     ****         NCOEF = 1 FOR ZIGM11                                      CNM.15
C     ****         NCOEF = 2 FOR ZIGM12 (=ZIGMC+ZIGM2)                       CNM.16
C     ****         NCOEF = 3 FOR ZIGM21 (=ZIGMC-ZIGM2)                       CNM.17
C     ****         NCOEF = 4 FOR ZIGM22                                      CNM.18
C     ****                                                                   CNM.19
C     ****     OUTPUT:                                                       CNM.20
C     ****       C(IMX,JMX,9)                                                CNM.21
C     ****         ARRAY FOR GRID POINT STENCILS AT RESOLUTION               CNM.22
C     ****           IMX BY JMX                                              CNM.23
C     ****                                                                   CNM.24
      DIMENSION ARRAY(IMX0+32,1),CS(1),C(IMX,JMX,10),WK(IMX0,3)               CNM.25
C     ****                                                                   CNM.26
C     ****     COMPUTE SEPARATION OF GRID POINTS OF RESOLUTION               CNM.27
C     ****       IMX BY JMX WITHIN GRID OF RESOLUTION IMX0 BY JMX0.          CNM.28
C     ****     EVALUATE DLON AND DLAT, GRID SPACING OF IMX BY JMX.           CNM.29
C     ****                                                                   CNM.30
      NINT=(IMX0-1)/(IMX-1)                                                  CNM.31
C     ****                                                                   CNM.32
C     ****     SCAN ARRAY IMX BY JMX CALCULATING AND ADDING                  CNM.33
C     ****       CONTRIBUTIONS TO STENCIL FROM ZIGM(NCOEF)                   CNM.34
C     ****                                                                   CNM.35
      I0=1-NINT                                                              CNM.36
      J0=1-NINT                                                              CNM.37
C     ****                                                                   CNM.38
C     ****     INDENTIFY COEFFICIENT UNDER CONSIDERATION                     CNM.39
C     ****                                                                   CNM.40
      IF(NCOEF.EQ.1)THEN                                                     CNM.41
C      ****                                                                  CNM.42
C      ****     CALCULATIONS FOR ZIGM11                                      CNM.43
C      ****                                                                  CNM.44
       DO J=2,JMX-1                                                          CNM.45
        DO I=1,IMX                                                           CNM.46
         C(I,J,1)=C(I,J,1)+.5*(ARRAY(I0+I*NINT,J0+J*NINT)+                   CNM.47
     1            ARRAY(I0+(I+1)*NINT,J0+J*NINT))                            CNM.48
         C(I,J,5)=C(I,J,5)+.5*(ARRAY(I0+I*NINT,J0+J*NINT)+                   CNM.49
     1             ARRAY(I0+(I-1)*NINT,J0+J*NINT))                           CNM.51
         C(I,J,9)=C(I,J,9)-.5*(ARRAY(I0+(I+1)*NINT,J0+J*NINT)+               CNM.51
     1            2.*ARRAY(I0+I*NINT,J0+J*NINT)+                             CNM.52
     2            ARRAY(I0+(I-1)*NINT,J0+J*NINT))                            CNM.53
        ENDDO                                                                CNM.54
       ENDDO                                                                 CNM.55
      ELSE IF(NCOEF.EQ.2)THEN                                                CNM.56
C      ****                                                                  CNM.57
C      ****     CALCULATIONS FOR ZIGM12  (=ZIGMC+ZIGM2)                      CNM.58
C      ****                                                                  CNM.59
       DO J=2,JMX-1                                                          CNM.60
        DO I=1,IMX                                                           CNM.61
         C(I,J,2)=C(I,J,2)+.5*(ARRAY(I0+I*NINT,J0+J*NINT)+                   CNM.62
     1            ARRAY(I0+(I+1)*NINT,J0+J*NINT))                            CNM.63
         C(I,J,4)=C(I,J,4)-.5*(ARRAY(I0+I*NINT,J0+J*NINT)+                   CNM.64
     1            ARRAY(I0+(I-1)*NINT,J0+J*NINT))                            CNM.65
         C(I,J,6)=C(I,J,6)+.5*(ARRAY(I0+I*NINT,J0+J*NINT)+                   CNM.66
     1            ARRAY(I0+(I-1)*NINT,J0+J*NINT))                            CNM.67
         C(I,J,8)=C(I,J,8)-.5*(ARRAY(I0+I*NINT,J0+J*NINT)+                   CNM.68
     1            ARRAY(I0+(I+1)*NINT,J0+J*NINT))                            CNM.69
         WK(I,1)=.5*(ARRAY(I0+(I+1)*NINT,J0+J*NINT)-                         CNM.70
     1           ARRAY(I0+(I-1)*NINT,J0+J*NINT))                             CNM.71
         WK(I,2)=(C(I,J,3)+WK(I,1))*(C(I,J,7)-WK(I,1))                       CNM.72
         WK(I,3)=DSIGN(WK(I,1),C(I,J,3)+C(I,J,7))                            CNM.73
         IF(WK(I,2).GE.0.)WK(I,3)=0.                                         CNM.74
         C(I,J,3)=C(I,J,3)+WK(I,1)+WK(I,3)                                   CNM.75
         C(I,J,7)=C(I,J,7)-WK(I,1)+WK(I,3)                                   CNM.76
         C(I,J,9)=C(I,J,9)-2.*WK(I,3)                                        CNM.77
        ENDDO                                                                CNM.78
       ENDDO                                                                 CNM.79
      ELSE IF(NCOEF.EQ.3)THEN                                                CNM.80
C      ****                                                                  CNM.81
C      ****     CALCULATIONS FOR ZIGM21  (=ZIGMC-ZIGM2)                      CNM.82
C      ****                                                                  CNM.83
       DO J=2,JMX-1                                                          CNM.84
        DO I=1,IMX                                                           CNM.85
         C(I,J,2)=C(I,J,2)+.5*(ARRAY(I0+I*NINT,J0+J*NINT)+                   CNM.86
     1            ARRAY(I0+I*NINT,J0+(J+1)*NINT))                            CNM.87
         C(I,J,4)=C(I,J,4)-.5*(ARRAY(I0+I*NINT,J0+J*NINT)+                   CNM.88
     1            ARRAY(I0+I*NINT,J0+(J+1)*NINT))                            CNM.89
         C(I,J,6)=C(I,J,6)+.5*(ARRAY(I0+I*NINT,J0+J*NINT)+                   CNM.90
     1            ARRAY(I0+I*NINT,J0+(J-1)*NINT))                            CNM.91
         C(I,J,8)=C(I,J,8)-.5*(ARRAY(I0+I*NINT,J0+J*NINT)+                   CNM.92
     1            ARRAY(I0+I*NINT,J0+(J-1)*NINT))                            CNM.93
         WK(I,1)=.5*(ARRAY(I0+I*NINT,J0+(J+1)*NINT)-                         CNM.94
     1           ARRAY(I0+I*NINT,J0+(J-1)*NINT))                             CNM.95
         WK(I,2)=(C(I,J,1)+WK(I,1))*(C(I,J,5)-WK(I,1))                       CNM.96
         WK(I,3)=DSIGN(WK(I,1),C(I,J,1)+C(I,J,5))                            CNM.97
         IF(WK(I,2).GE.0.)WK(I,3)=0.                                         CNM.98
         C(I,J,1)=C(I,J,1)+WK(I,1)+WK(I,3)                                   CNM.99
         C(I,J,5)=C(I,J,5)-WK(I,1)+WK(I,3)                                   CNM.100
         C(I,J,9)=C(I,J,9)-2.*WK(I,3)                                        CNM.101
        ENDDO                                                                CNM.102
       ENDDO                                                                 CNM.103
C      ****                                                                  CNM.104
C      ****     LOW LATITUDE BOUNDARY CONDITION                              CNM.105
C      ****                                                                  CNM.106
       J=1                                                                   CNM.107
       DO I=1,IMX                                                            CNM.108
        C(I,J,2)=C(I,J,2)+.5*(ARRAY(I0+I*NINT,J0+J*NINT)+                    CNM.109
     1           ARRAY(I0+I*NINT,J0+(J+1)*NINT))                             CNM.110
        C(I,J,4)=C(I,J,4)-.5*(ARRAY(I0+I*NINT,J0+J*NINT)+                    CNM.111
     1            ARRAY(I0+I*NINT,J0+(J+1)*NINT))                            CNM.112
        WK(I,1)=.5*(ARRAY(I0+I*NINT,J0+J*NINT)+                              CNM.113
     1          ARRAY(I0+I*NINT,J0+(J+1)*NINT))                              CNM.114
        WK(I,2)=(C(I,J,1)+WK(I,1))*(C(I,J,5)-WK(I,1))                        CNM.115
        WK(I,3)=DSIGN(WK(I,1),C(I,J,1)+C(I,J,5))                             CNM.116
        IF(WK(I,2).GE.0.)WK(I,3)=0.                                          CNM.117
        C(I,J,1)=C(I,J,1)+WK(I,1)+WK(I,3)                                    CNM.118
        C(I,J,5)=C(I,J,5)-WK(I,1)+WK(I,3)                                    CNM.119
        C(I,J,9)=C(I,J,9)-2.*WK(I,3)                                         CNM.120
       ENDDO                                                                 CNM.121
      ELSE IF(NCOEF.EQ.4)THEN                                                CNM.122
C      ****                                                                  CNM.123
C      ****     CALCULATIONS FOR ZIGM22                                      CNM.124
C      ****                                                                  CNM.125
       DO J=2,JMX-1                                                          CNM.126
        DO I=1,IMX                                                           CNM.127
         C(I,J,3)=C(I,J,3)+.5*(ARRAY(I0+I*NINT,J0+J*NINT)                    CNM.128
     1            +ARRAY(I0+I*NINT,J0+(J+1)*NINT))                           CNM.129
         C(I,J,7)=C(I,J,7)+.5*(ARRAY(I0+I*NINT,J0+J*NINT)                    CNM.130
     1            +ARRAY(I0+I*NINT,J0+(J-1)*NINT))                           CNM.131
         C(I,J,9)=C(I,J,9)-.5*(ARRAY(I0+I*NINT,J0+(J-1)*NINT)                CNM.132
     1            +2*ARRAY(I0+I*NINT,J0+J*NINT)                              CNM.133
     2            +ARRAY(I0+I*NINT,J0+(J+1)*NINT))                           CNM.134
        ENDDO                                                                CNM.135
       ENDDO                                                                 CNM.136
C      ****                                                                  CNM.137
C      ****     LOW LATITUDE BOUNDARY CONDITION                              CNM.138
C      ****                                                                  CNM.139
       J=1                                                                   CNM.140
       DO I=1,IMX                                                            CNM.141
        C(I,J,3)=C(I,J,3)+.5*(ARRAY(I0+I*NINT,J0+J*NINT)                     CNM.142
     1           +ARRAY(I0+I*NINT,J0+(J+1)*NINT))                            CNM.143
        C(I,J,9)=C(I,J,9)-.5*(ARRAY(I0+I*NINT,J0+J*NINT)                     CNM.144
     2           +ARRAY(I0+I*NINT,J0+(J+1)*NINT))                            CNM.145
       ENDDO                                                                 CNM.146
      ENDIF                                                                  CNM.147
      RETURN                                                                 CNM.148
      END                                                                    CNM.149
C                                                                            CNM.150
*DK COLATH                                                                COLATH.1
      SUBROUTINE COLATH                                                   COLATH.2
C     ****                                                                COLATH.3
C     ****     FOR EACH NORTHERN HEMISPHERE GEOMAGNETIC GRID POINT,       COLATH.4
C     ****       DETERMINE COLATITUDE IN NORTHERN AURORAL CIRCLE          COLATH.5
C     ****       COORDINATES.  PLACE VALUES IN COLATC(IMAXMP,JMX0).       COLATH.6
C     ****       ALSO CALCULATE FRACTIONAL PRESENCE OF DYNAMO             COLATH.7
C     ****       EQUATION GIVEN CRITICAL COLATITUDES CRIT(2)              COLATH.8
C     ****                                                                COLATH.9
      include "param.h"                                                   COLATH.10
      INCLUDE 'epara.h'                                                   COLATH.11
      include "econs.h"                                                   COLATH.12
      COMMON/TRGM/RLATM(ZIMXP,ZJMX),RLONM(ZIMXP,ZJMX),                    COLATH.13
     1  DIPMAG(ZIMXP,ZJMX),DECMAG(ZIMXP,ZJMX),SNDEC(ZIMXP,ZJMX),          COLATH.14
     2  CSDEC(ZIMXP,ZJMX),SN2DEC(ZIMXP,ZJMX),SNCSDC(ZIMXP,ZJMX),          COLATH.15
     3  RLATMP(ZJMX),RLONMP(ZJMX),DLONS(ZJMX)                             COLATH.16
      COMMON /IONCR/ ISTAR,IHEM,THETA0(2),OFFA(2),OFFC(2),DSKOFA(2),      COLATH.17
     1  DSKOFC(2),PHID(2),PHIN(2),PHIDP0(2),PHIDM0(2),PHINP0(2),          COLATH.18
     2  PHINM0(2),PSIM(2),PSIE(2),PCEN(2),RR1(2),COLATC(IMAXMP,JMX0),     COLATH.19
     3  P(IMX0,JMX0),CRIT(2)                                              COLATH.20
C     ****                                                                COLATH.21
C     ****     PICK RELEVANT NORTHERN PARAMETERS FROM /IONCR/             COLATH.22
C     ****                                                                COLATH.23
      OFDC=DSQRT(OFFC(2)**2+DSKOFC(2)**2)                                 COLATH.24
      COSOFC=DCOS(OFDC)                                                   COLATH.25
      SINOFC=DSIN(OFDC)                                                   COLATH.26
      ASLONC=DASIN(DSKOFC(2)/OFDC)                                        COLATH.27
C     ****                                                                COLATH.28
C     ****     FILL ARRAY COLATC WITH REQUIRED COLATITUDE VALUES          COLATH.29
C     ****                                                                COLATH.30
      DO J=1,JMX0                                                         COLATH.31
       SINLAT=DSIN(DABS(YLATM(J+JMX0-1)))                                 COLATH.32
       COSLAT=DCOS(YLATM(J+JMX0-1))                                       COLATH.33
       DO I=1,IMAXMP                                                      COLATH.34
        COLATC(I,J)=DCOS(YLONM(I)-DLONS(1)+ASLONC)                        COLATH.35
        COLATC(I,J)=COSOFC*SINLAT-SINOFC*COSLAT*COLATC(I,J)               COLATH.36
        COLATC(I,J)=DACOS(COLATC(I,J))                                    COLATH.37
       ENDDO                                                              COLATH.38
      ENDDO                                                               COLATH.39
C     ****                                                                COLATH.40
C     ****     CALCULATE FRACTIONAL PRESENCE OF DYNAMO EQUATION           COLATH.41
C     ****       AT EACH NORTHERN HEMISPHERE GEOMAGNETIC GRID             COLATH.42
C     ****       POINT.  PLACE VALUES IN P(IMX0,JMX0)                     COLATH.43
C     ****                                                                COLATH.44
      DO J=1,JMX0                                                         COLATH.45
       DO I=1,IMX0                                                        COLATH.46
        P(I,J)=(COLATC(I,J)-CRIT(1))/(CRIT(2)-CRIT(1))                    COLATH.47
        IF(P(I,J).LT.0.)P(I,J)=0.                                         COLATH.48
        IF(P(I,J).GT.1.)P(I,J)=1.                                         COLATH.49
       ENDDO                                                              COLATH.50
      ENDDO                                                               COLATH.51
      RETURN                                                              COLATH.52
      END                                                                 COLATH.53
C                                                                         COLATH.54
*DK DIVIDE                                                                DIVIDE.1
      SUBROUTINE DIVIDE(C,IMX,JMX,IMX0,JMX0,CS)                           DIVIDE.2
      IMPLICIT REAL*8 (A-H,O-Z)                                           DIVIDE.3
C     ****                                                                DIVIDE.4
C     ****     DIVIDE STENCIL, C, BY COS(THETA(I,J))                      DIVIDE.5
C     ****                                                                DIVIDE.6
      DIMENSION C(IMX,JMX,10),CS(1)                                        DIVIDE.7
      NINT=(IMX0-1)/(IMX-1)                                               DIVIDE.8
      J0=1-NINT                                                           DIVIDE.9
      DO N=1,9                                                            DIVIDE.10
       DO J=1,JMX-1                                                       DIVIDE.11
        DO I=1,IMX                                                        DIVIDE.12
         C(I,J,N)=C(I,J,N)/(CS(J0+J*NINT)*NINT**2)                        DIVIDE.13
        ENDDO                                                             DIVIDE.14
       ENDDO                                                              DIVIDE.15
      ENDDO                                                               DIVIDE.16
      IF(NINT.EQ.1)THEN                                                   DIVIDE.17
       DO I=1,IMX                                                         DIVIDE.18
        C(I,1,10)=C(I,1,10)/CS(1)                                         DIVIDE.19
       ENDDO                                                              DIVIDE.20
      ENDIF                                                               DIVIDE.21
      RETURN                                                              DIVIDE.22
      END                                                                 DIVIDE.23
C                                                                         DIVIDE.24
*DK DYN2                                                                    DYN2.1
      SUBROUTINE DYN2                                                       DYN2.2
C     ****                                                                  DYN2.3
C     ****     SETS UP COEFFICIENTS FOR DYNAMO PDE                          DYN2.4
C     ****                                                                  DYN2.5
      include "param.h"                                                     DYN2.6
      INCLUDE 'epara.h'                                                     DYN2.7
      INCLUDE 'econs.h'                                                     DYN2.8
      include "ecoef.h"                                                     DYN2.9
      include "edyna.h"                                                     DYN2.10
      include "efild.h"                                                     DYN2.11
      include "eterp.h"                                                     DYN2.12
      include "eceee.h"                                                     DYN2.14
      include "edphi.h"                                                     DYN2.15
      include "start.h"                                                     DYN2.16
      DIMENSION ARRAY(-15:IMX0+16,JMX0),CS(JMX0)                            DYN2.17
      DIMENSION C0(IMX0,JMX0,10),C1(IMX1,JMX1,9),C2(IMX2,JMX2,9),           DYN2.18
     1          C3(IMX3,JMX3,9),C4(IMX4,JMX4,9)                             DYN2.19
      DIMENSION WK(IMX0,3)                                                  DYN2.20
      EQUIVALENCE (CEE,C0),(CEE(NC1),C1),(CEE(NC2),C2),(CEE(NC3),C3),       DYN2.21
     1            (CEE(NC4),C4)                                             DYN2.22
      common/eeeeee/zx(640,350,30,7),zy(640,350,30,32)
C     ****                                                                  DYN2.23
C     ****     TRANSFORM NEEDED FIELDS TO GEOMAGNETIC COORDINATES.          DYN2.24
C     ****       PERFORM FIELD-LINE INTEGRATIONS                            DYN2.25
C     ****       EVALUATE PDE COEFFICIENTS AND RHS                          DYN2.26
C     ****                                                                  DYN2.27
      CALL TRANSF2                                                          DYN2.28
C     ****                                                                  DYN2.29
C     ****     SET ARRAY NC AND MAGNETIC LATITUDE COSINE ARRAY              DYN2.30
C     ****                                                                  DYN2.31
      NC(1)=NC0                                                             DYN2.32
      NC(2)=NC1                                                             DYN2.33
      NC(3)=NC2                                                             DYN2.34
      NC(4)=NC3                                                             DYN2.35
      NC(5)=NC4                                                             DYN2.36
      NC(6)=NCEE                                                            DYN2.37
      DO J=1,JMX0                                                           DYN2.38
       CS(J)=DCOS(PI/2.-(JMX0-J)*DLATM)                                     DYN2.39
      ENDDO                                                                 DYN2.40
C     ****                                                                  DYN2.41
C     ****     SET UP DIFFERENCE COEFFICIENTS                               DYN2.42
C     ****                                                                  DYN2.43
C     ****       REPLACE ZIGM11 BY A,  ZIGM22 BY B,  ZIGMC BY C             DYN2.44
C     ****         AND ZIGM2 BY D                                           DYN2.45
C     ****                                                                  DYN2.46
      J0=JMX0-JMAXMH                                                        DYN2.47
      DO J=1,JMAXMH                                                         DYN2.48
       JJ=JMAXMH+J-1                                                        DYN2.49
       JJJ=JMAXMH-J+1                                                       DYN2.50
       DO I=1,IMAXMP                                                        DYN2.51
        ZIGMC(I,JJ)=(ZIGMC(I,JJ)+ZIGM2(I,JJ))/(4*DLATM*DLONM)               DYN2.52
        ZIGM2(I,JJ)=ZIGMC(I,JJ)-2.*ZIGM2(I,JJ)/(4*DLATM*DLONM)              DYN2.53
        ZIGM22(I,JJ)=ZIGM22(I,JJ)*CS(J0+J)/DLATM**2                         DYN2.54
        ZIGMC(I,JJJ)=-ZIGMC(I,JJ)                                           DYN2.55
        ZIGM2(I,JJJ)=-ZIGM2(I,JJ)                                           DYN2.56
        ZIGM22(I,JJJ)=ZIGM22(I,JJ)                                          DYN2.57
       ENDDO                                                                DYN2.58
       IF(J.NE.JMAXMH)THEN                                                  DYN2.59
        DO I=1,IMAXMP                                                       DYN2.60
         ZIGM11(I,JJ)=ZIGM11(I,JJ)/(CS(J0+J)*DLONM**2)                      DYN2.61
         ZIGM11(I,JJJ)=ZIGM11(I,JJ)                                         DYN2.62
        ENDDO                                                               DYN2.63
       ENDIF                                                                DYN2.64
      ENDDO                                                                 DYN2.65
C     ****                                                                  DYN2.66
C     ****     SET ZIGM11 TO ZERO AT MEGNETIC POLES TO AVOID FLOATING       DYN2.67
C     ****       EXCEPTION  ( VALUES AT POLES ARE NOT USED )                DYN2.68
C     ****                                                                  DYN2.69
      DO I=1,IMAXMP                                                         DYN2.70
        ZIGM11(I,1)=0.0                                                     DYN2.71
        ZIGM11(I,JMAXM)=0.0                                                 DYN2.72
      ENDDO                                                                 DYN2.73
C     ****                                                                  DYN2.74
C     ****     CLEAR ARRAY FOR DIFFERENCE STENCILS AT ALL LEVELS            DYN2.75
C     ****                                                                  DYN2.76
      CALL CLEAR(CEE,IMX0,JMX0)                                             DYN2.77
C     ****                                                                  DYN2.78
C     ****     CALCULATE CONTRIBUTION TO STENCILS FROM EACH PDE             DYN2.79
C     ****       COEFFICIENT                                                DYN2.80
C     ****                                                                  DYN2.81
      ARRAY(1,1)=1.                                                         DYN2.82
      SYM=1.                                                                DYN2.83
      CALL STENCIL(ZIGM11,CS(1),IMX0,JMX0,SYM,CEE,ARRAY(1,1),1,WK)          DYN2.84
      ARRAY(1,1)=4.                                                         DYN2.85
      SYM=1.                                                                DYN2.86
      CALL STENCIL(ZIGM22,CS(1),IMX0,JMX0,SYM,CEE,ARRAY(1,1),4,WK)          DYN2.87
      ARRAY(1,1)=2.                                                         DYN2.88
      SYM=-1.                                                               DYN2.89
      CALL STENCIL(ZIGMC,CS(1),IMX0,JMX0,SYM,CEE,ARRAY(1,1),2,WK)           DYN2.90
      ARRAY(1,1)=3.                                                         DYN2.91
      SYM=-1.                                                               DYN2.92
      CALL STENCIL(ZIGM2,CS(1),IMX0,JMX0,SYM,CEE,ARRAY(1,1),3,WK)           DYN2.93
C     ****                                                                  DYN2.94
C     ****     CALL RTHS TO INSERT RHS IN FINEST STENCIL                    DYN2.95
C     ****                                                                  DYN2.96
      CALL RTHS(RHS,IMX0,JMX0,JMAXMH,C0)                                    DYN2.97
C     ****                                                                  DYN2.98
C     ****      SET BOUNDARY CONDITION AT EQUATOR AND POLE                  DYN2.99
C     ****                                                                  DYN2.100
      CALL EDGES(C0,IMX0,JMX0)                                              DYN2.101
      CALL EDGES(C1,IMX1,JMX1)                                              DYN2.102
      CALL EDGES(C2,IMX2,JMX2)                                              DYN2.103
      CALL EDGES(C3,IMX3,JMX3)                                              DYN2.104
      CALL EDGES(C4,IMX4,JMX4)                                              DYN2.105
C     ****                                                                  DYN2.106
C     ****     DIVIDE STENCILS BY COS(THETA)                                DYN2.107
C     ****                                                                  DYN2.108
      CALL DIVIDE(C0,IMX0,JMX0,IMX0,JMX0,CS(1))                             DYN2.109
      CALL DIVIDE(C1,IMX1,JMX1,IMX0,JMX0,CS(1))                             DYN2.110
      CALL DIVIDE(C2,IMX2,JMX2,IMX0,JMX0,CS(1))                             DYN2.111
      CALL DIVIDE(C3,IMX3,JMX3,IMX0,JMX0,CS(1))                             DYN2.112
      CALL DIVIDE(C4,IMX4,JMX4,IMX0,JMX0,CS(1))                             DYN2.113
C     ****                                                                  DYN2.114
C     ****     SET VALUE OF SOLUTION TO 1. AT POLE                          DYN2.115
C     ****                                                                  DYN2.116
      DO I=1,IMX0                                                           DYN2.117
       C0(I,JMX0,10)=1.                                                     DYN2.118
      ENDDO                                                                 DYN2.119
C     ****                                                                  DYN2.120
C     ****     MODIFY STENCILS AND RHS SO THAT HEELIS POTENTIAL             DYN2.121
C     ****       ID INSERTED AT HIGH LATITUDE                               DYN2.122
C     ****                                                                  DYN2.123
      CALL ADHEEL                                                           DYN2.124
      JNTL=0                                                                DYN2.125
C     ****                                                                  DYN2.126
C     ****     CALL MUDPAC TO SOLVE PDE                                     DYN2.127
C     ****       SOLUTION RETURNED IN RIM(1)                                DYN2.128
C     ****                                                                  DYN2.129
      CALL MUD(PHIM,JNTL)                                                   DYN2.130
C     ****                                                                  DYN2.126
C     ****     MAP TWO-DIMENSIONAL POTENTAIL TO A 3-D ARRAY                 DYN2.127
C     ****                                                                  DYN2.129
      CALL THREED                                                           DYN2.139
      RETURN                                                                DYN2.140
      END                                                                   DYN2.141
C                                                                           DYN2.142
*DK EDGES                                                                  EDGES.1
      SUBROUTINE EDGES(C,IMX,JMX)                                          EDGES.2
      IMPLICIT REAL*8 (A-H,O-Z)                                            EDGES.3
C     ****                                                                 EDGES.4
C     ****     INSERT EQUATORIAL AND POLAR BOUNDARY CONDITIONS IN          EDGES.5
C     ****       STENCIL C(IMX,JMX,9)                                      EDGES.6
C     ****                                                                 EDGES.7
      DIMENSION C(IMX,JMX,10)                                               EDGES.8
C     ****                                                                 EDGES.9
C     ****     POLE                                                        EDGES.10
C     ****                                                                 EDGES.11
      DO N=1,8                                                             EDGES.12
       DO I=1,IMX                                                          EDGES.13
        C(I,JMX,N)=0.                                                      EDGES.14
       ENDDO                                                               EDGES.15
      ENDDO                                                                EDGES.16
      DO I=1,IMX                                                           EDGES.17
       C(I,JMX,9)=1.                                                       EDGES.18
      ENDDO                                                                EDGES.19
      RETURN                                                               EDGES.20
      END                                                                  EDGES.21
C                                                                          EDGES.22
*DK FLOW                                                                    FLOW.1
      SUBROUTINE FLOWION(DCUSP,DALFA,DALFA2,DFLUX,DFLUX2,DDRIZL)            FLOW.2
C     ****                                                                  FLOW.4
C     ****     COMPUTE ELECTRIC FIELD COMPONENTS BY SPACE CENTERED          FLOW.5
C     ****     TECHNIQUE(VOLTS/M)                                           FLOW.6
C     ****                                                                  FLOW.7
      INCLUDE 'param.h'                                                     FLOW.8
      INCLUDE 'const.h'                                                     FLOW.9
      INCLUDE 'phycs.h'                                                     FLOW.10
      COMMON/TRGM/RLATM(ZIMXP,ZJMX),RLONM(ZIMXP,ZJMX),                      FLOW.11
     1      DIPMAG(ZIMXP,ZJMX),DECMAG(ZIMXP,ZJMX),SNDEC(ZIMXP,ZJMX),        FLOW.12
     2      CSDEC(ZIMXP,ZJMX),SN2DEC(ZIMXP,ZJMX),SNCSDC(ZIMXP,ZJMX),        FLOW.13
     3      RLATMP(ZJMX),RLONMP(ZJMX),DLONS(ZJMX)                           FLOW.14
      DIMENSION DCUSP(ZIMX),DALFA(ZIMX),DFLUX(ZIMX),DDRIZL(ZIMX),           FLOW.15
     1          DALFA2(ZIMX),DFLUX2(ZIMX),DLAT(ZIMX),DLON(ZIMX)             FLOW.16
      COMMON/SUNANGLE/DLONSUN      
C     ****                                                                  FLOW.17
C     ****     SET UP DIPOLE COORDINATES                                    FLOW.18
C     ****                                                                  FLOW.19
      DO I=1,IMAX                                                           FLOW.20
       DLAT(I)=RLATM(I+2,J)                                                 FLOW.21
       DLON(I)=RLONM(I+2,J)-DLONS(J)                                        FLOW.22
      ENDDO                                                                 FLOW.23
      DLONSUN=DLONS(J)
C     ****                                                                  FLOW.24
C     ****     COMPUTE POTENTIAL                                            FLOW.25
C     ****                                                                  FLOW.26
      CALL FLOWV3(DLAT,DLON)                                                FLOW.27
      CALL POLHT(DCUSP,IMAX)                                                FLOW.28
      CALL AURHT(DLAT,DLON,DALFA,DFLUX,DDRIZL,DALFA2,DFLUX2,IMAX)           FLOW.29
      RETURN                                                                FLOW.30
      END                                                                   FLOW.31
C                                                                           FLOW.32
*DK FLOWV3                                                                FLOWV3.1
      SUBROUTINE FLOWV3(DLAT,DLON)                                        FLOWV3.2
C     ****                                                                FLOWV3.3
C     ****     COMPUTE ELECTRIC POTENTIAL                                 FLOWV3.4
C     ****                                                                FLOWV3.5
      INCLUDE 'param.h'                                                   FLOWV3.6
      INCLUDE 'const.h'                                                   FLOWV3.7
      DIMENSION DLAT(ZIMX),DLON(ZIMX)                                     FLOWV3.8
      COMMON /IONCR/ ISTAR,IHEM,THETA0(2),OFFA(2),OFFC(2),DSKOFA(2),      FLOWV3.9
     |DSKOFC(2),PHID(2),PHIN(2),PHIDP0(2),PHIDM0(2),PHINP0(2),PHINM0(2),  FLOWV3.10
     2 PSIM(2),PSIE(2),PCEN(2),R1(2),COLATC(IMAXMP,JMX0),
     3 P(IMX0,JMX0),CRIT(2)                                               FLOWV3.11
      COMMON /OVALPOS/ COSOFA(2),SINOFA(2),ASLONA(2),PI,PI2               FLOWV3.12
      COMMON/ANGLE/ALON(ZIMX),COLAT(ZIMX)                                 FLOWV3.13
      DIMENSION SINLAT(ZIMX),COSLAT(ZIMX),SINLON(ZIMX),COSLON(ZIMX)       FLOWV3.14
      DATA E/1.D-6/,PIH/1.5707963/                                        FLOWV3.16
C     ****                                                                FLOWV3.16
C     ****     CONVERT PARAMETERS TO RADIANS ETC.                         FLOWV3.17
C     ****                                                                FLOWV3.18
      IF(ISTAR.EQ.0)THEN                                                  FLOWV3.19
       ISTAR=1                                                            FLOWV3.20
       PI=4.*DATAN(1.D0)                                                  FLOWV3.21
       PI2=2.*PI                                                          FLOWV3.22
       PIH=.5*PI                                                          FLOWV3.23
       DO N=1,2                                                           FLOWV3.24
        IF(ABS(OFFA(N)+DSKOFA(N)).GT.0.01)THEN                            FLOWV3.25
         OFDA=DSQRT(OFFA(N)**2+DSKOFA(N)**2)                              FLOWV3.26
         COSOFA(N)=DCOS(OFDA)                                             FLOWV3.27
         SINOFA(N)=DSIN(OFDA)                                             FLOWV3.28
         ASLONA(N)=DASIN(DSKOFA(N)/OFDA)                                  FLOWV3.29
        ELSE                                                              FLOWV3.30
         OFDA=0.                                                          FLOWV3.31
         COSOFA(N)=1.                                                     FLOWV3.32
         SINOFA(N)=0.                                                     FLOWV3.33
         ASLONA(N)=0.                                                     FLOWV3.34
        ENDIF                                                             FLOWV3.35
       ENDDO                                                              FLOWV3.36
      ENDIF                                                               FLOWV3.37
C     ****                                                                FLOWV3.38
C     ****               SET IHEM = 1,2  FOR  S,N HEMISPHERE              FLOWV3.39
C     ****                                                                FLOWV3.40
      IMAXD2=MAX0(1,IMAX/2)                                               FLOWV3.41
      IHEM=IDINT(DLAT(IMAXD2)*2./3.1416+2.)                               FLOWV3.42
C     ****                                                                FLOWV3.43
C     ****     TRANSFORM TO AURORAL CIRCLE CORDINATES                     FLOWV3.44
C     ****     I. E. CALCULATE COLAT AND ALON FOR POLHT AND AURHT         FLOWV3.45
C     ****                                                                FLOWV3.46
      DO I=1,IMAX                                                         FLOWV3.47
       SINLAT(I)=DSIN(DABS(DLAT(I)))                                      FLOWV3.48
       COSLAT(I)=DCOS(DLAT(I))                                            FLOWV3.49
       SINLON(I)=DSIN(DLON(I)+ASLONA(IHEM))                               FLOWV3.50
       COSLON(I)=DCOS(DLON(I)+ASLONA(IHEM))                               FLOWV3.51
       COLAT(I)=COSOFA(IHEM)*SINLAT(I)-SINOFA(IHEM)*COSLAT(I)*            FLOWV3.52
     1            COSLON(I)                                               FLOWV3.53
       COLAT(I)=DACOS(COLAT(I))                                           FLOWV3.54
       ALON(I)=DMOD(DATAN2(+SINLON(I)*COSLAT(I),SINLAT(I)*SINOFA(IHEM)+   FLOWV3.55
     1         COSOFA(IHEM)*COSLAT(I)*COSLON(I))-ASLONA(IHEM)+            FLOWV3.56
     2         3.*PI,PI2)-PI                                              FLOWV3.57
      ENDDO                                                               FLOWV3.58
      RETURN                                                              FLOWV3.59
      END                                                                 FLOWV3.60
C                                                                        FLOWV32.61
*DK FLOWV32                                                              FLOWV32.1
      SUBROUTINE FLOWV32(POTEN)                                          FLOWV32.2
C     ****                                                               FLOWV32.3
C     ****     COMPUTE ELECTRIC POTENTIAL                                FLOWV32.4
C     ****                                                               FLOWV32.5
      INCLUDE 'param.h'                                                  FLOWV32.6
      INCLUDE 'const.h'                                                  FLOWV32.7
C     ****                                                               FLOWV32.8
C     ****     WORK SPACE                                                FLOWV32.9
C     ****                                                               FLOWV32.10
      COMMON/POTE/ DLAT(IMAXM),DLON(IMAXM),PSI(8),PHI(IMAXM,8),          FLOWV32.11
     1  IFN(IMAXM),PHIFUN(IMAXM),COLAT(IMAXM),ALON(IMAXM),               FLOWV32.12
     2  SINLAT(IMAXM),COSLAT(IMAXM),SINLON(IMAXM),COSLON(IMAXM),         FLOWV32.13
     3  WK1(IMAXM),WK2(IMAXM),WK3(IMAXM),PHIFN2(IMAXM)                   FLOWV32.14
      COMMON /IONCR/ ISTAR,IHEM,THETA0(2),OFFA(2),OFFC(2),DSKOFA(2),     FLOWV32.15
     2    DSKOFC(2),PHID(2),PHIN(2),PHIDP0(2),PHIDM0(2),PHINP0(2),       FLOWV32.16
     2    PHINM0(2),PSIM(2),PSIE(2),PCEN(2),R1(2),COLATC(IMAXMP,JMX0),
     3    P(IMX0,JMX0),CRIT(2)                                           FLOWV32.17
      COMMON /OVALPOS/ COSOFA(2),SINOFA(2),ASLONA(2),PI,PI2              FLOWV32.18
      COMMON/POTE1/PHDPMX(2),PHDMMX(2),PHNPMX(2),PHNMMX(2),              FLOWV32.19
     1          COSOFC(2),SINOFC(2),ASLONC(2)                            FLOWV32.20
      DIMENSION POTEN(IMAXMP)                                            FLOWV32.21
      DATA E/1.D-6/                                                      FLOWV32.22
C     ****                                                               FLOWV32.23
C     ****     CONVERT PARAMETERS TO RADIANS ETC.                        FLOWV32.24
C     ****                                                               FLOWV32.25
      IF(ISTAR.EQ.0)THEN                                                 FLOWV32.26
       ISTAR=1                                                           FLOWV32.27
       PI=4.*DATAN(1.D0)                                                 FLOWV32.28
       PI2=2.*PI                                                         FLOWV32.29
       PIH=.5*PI                                                         FLOWV32.30
       DO N=1,2                                                          FLOWV32.31
        IF(ABS(OFFA(N)+DSKOFA(N)).GT.0.01)THEN                           FLOWV32.32
         OFDA=DSQRT(OFFA(N)**2+DSKOFA(N)**2)                             FLOWV32.33
         COSOFA(N)=DCOS(OFDA)                                            FLOWV32.34
         SINOFA(N)=DSIN(OFDA)                                            FLOWV32.35
         ASLONA(N)=DASIN(DSKOFA(N)/OFDA)                                 FLOWV32.36
        ELSE                                                             FLOWV32.37
         OFDA=0.                                                         FLOWV32.38
         COSOFA(N)=1.                                                    FLOWV32.39
         SINOFA(N)=0.                                                    FLOWV32.40
         ASLONA(N)=0.                                                    FLOWV32.41
        ENDIF                                                            FLOWV32.42
        IF(ABS(OFFC(N)+DSKOFC(N)).GT.0.01)THEN                           FLOWV32.43
         OFDC=DSQRT(OFFC(N)**2+DSKOFC(N)**2)                             FLOWV32.44
         COSOFC(N)=DCOS(OFDC)                                            FLOWV32.45
         SINOFC(N)=DSIN(OFDC)                                            FLOWV32.46
         ASLONC(N)=DASIN(DSKOFC(N)/OFDC)                                 FLOWV32.47
        ELSE                                                             FLOWV32.48
         OFDC=0.                                                         FLOWV32.49
         COSOFC(N)=1.                                                    FLOWV32.50
         SINOFC(N)=0.                                                    FLOWV32.51
         ASLONC(N)=0.                                                    FLOWV32.52
        ENDIF                                                            FLOWV32.53
        IF(PHIN(N).LT.PHID(N))PHIN(N)=PHIN(N)+PI2                        FLOWV32.54
        PHDPMX(N)=.5*DMIN1(PI,(PHIN(N)-PHID(N)))                         FLOWV32.55
        PHNPMX(N)=.5*DMIN1(PI,(PHID(N)-PHIN(N)+PI2))                     FLOWV32.56
        PHNMMX(N)=PHDPMX(N)                                              FLOWV32.57
        PHDMMX(N)=PHNPMX(N)                                              FLOWV32.58
       ENDDO                                                             FLOWV32.59
      ENDIF                                                              FLOWV32.60
C     ****                                                               FLOWV32.61
C     ****               SET IHEM = 1,2  FOR  S,N HEMISPHERE             FLOWV32.62
C     ****                                                               FLOWV32.63
      IMAXD2=MAX0(1,IMAX/2)                                              FLOWV32.64
      IHEM=IDINT(DLAT(IMAXD2)*2./3.1416+2.)                              FLOWV32.65
      SINTH0=DSIN(THETA0(IHEM))                                          FLOWV32.66
C     ****                                                               FLOWV32.67
C     ****     AVERAGE AMIE RESULTS SHOW R1=-2.6 FOR 11.3 DEGREES        FLOWV32.68
C     ****       (0.1972 RAD) BEYOND THETA0.                             FLOWV32.69
C     ****                                                               FLOWV32.70
      SINTHR1=DSIN(THETA0(IHEM)+0.1972)                                  FLOWV32.71
      PSI(1)=PSIE(IHEM)                                                  FLOWV32.72
      PSI(3)=PSIM(IHEM)                                                  FLOWV32.73
      DO N=2,4,2                                                         FLOWV32.74
       PSI(N)=PSI(N-1)                                                   FLOWV32.75
      ENDDO                                                              FLOWV32.76
      DO N=1,4                                                           FLOWV32.77
       PSI(N+4)=PSI(N)                                                   FLOWV32.78
      ENDDO                                                              FLOWV32.79
C     ****                                                               FLOWV32.80
C     ****     TRANSFORM TO AURORAL CIRCLE CORDINATES                    FLOWV32.81
C     ****                                                               FLOWV32.82
      DO I=1,IMAXM                                                       FLOWV32.83
       SINLAT(I)=DSIN(DABS(DLAT(I)))                                     FLOWV32.84
       COSLAT(I)=DCOS(DLAT(I))                                           FLOWV32.85
       SINLON(I)=DSIN(DLON(I)+ASLONC(IHEM))                              FLOWV32.86
       COSLON(I)=DCOS(DLON(I)+ASLONC(IHEM))                              FLOWV32.87
       COLAT(I)=COSOFC(IHEM)*SINLAT(I) - SINOFC(IHEM)*COSLAT(I)*         FLOWV32.88
     1            COSLON(I)                                              FLOWV32.89
       ALON(I)=DMOD(DATAN2(+SINLON(I)*COSLAT(I),SINLAT(I)*SINOFC(IHEM)+  FLOWV32.90
     1         COSOFC(IHEM)*COSLAT(I)*COSLON(I))-ASLONC(IHEM)+           FLOWV32.91
     2         3.*PI,PI2)-PI                                             FLOWV32.92
       COLAT(I)=DACOS(COLAT(I))                                          FLOWV32.93
C      ****                                                              FLOWV32.94
C      ****     CALCULATE BOUNDARIES FOR LONGITUDINAL FUNCTION           FLOWV32.95
C      ****                                                              FLOWV32.96
       WK1(I)=((COLAT(I)-THETA0(IHEM))/THETA0(IHEM))**2                  FLOWV32.97
       PHI(I,4)=PHID(IHEM)+E-DMIN1(PHIDM0(IHEM)+WK1(I)*                  FLOWV32.98
     1         (PIH-PHIDM0(IHEM)),PHDMMX(IHEM))                          FLOWV32.99
       PHI(I,5)=PHID(IHEM)-E+DMIN1(PHIDP0(IHEM)+WK1(I)*                  FLOWV32.100
     1         (PIH-PHIDP0(IHEM)),PHDPMX(IHEM))                          FLOWV32.101
       PHI(I,6)=PHIN(IHEM)+E-DMIN1(PHINM0(IHEM)+WK1(I)*                  FLOWV32.102
     1         (PIH-PHINM0(IHEM)),PHNMMX(IHEM))                          FLOWV32.103
       PHI(I,7)=PHIN(IHEM)-E+DMIN1(PHINP0(IHEM)+WK1(I)*                  FLOWV32.104
     1         (PIH-PHINP0(IHEM)),PHNPMX(IHEM))                          FLOWV32.105
       PHI(I,1)=PHI(I,5)-PI2                                             FLOWV32.106
       PHI(I,2)=PHI(I,6)-PI2                                             FLOWV32.107
       PHI(I,3)=PHI(I,7)-PI2                                             FLOWV32.108
       PHI(I,8)=PHI(I,4)+PI2                                             FLOWV32.109
       PHIFUN(I)=0.                                                      FLOWV32.110
       PHIFN2(I)=0.                                                      FLOWV32.111
       IF((COLAT(I)-THETA0(IHEM)).GE.0.0)THEN                            FLOWV32.112
        IFN(I)=3                                                         FLOWV32.113
       ELSE                                                              FLOWV32.114
        IFN(I)=2                                                         FLOWV32.115
       ENDIF                                                             FLOWV32.116
C      ****                                                              FLOWV32.117
C      ****    ADD RING CURRENT ROTATION TO POTENTIAL (PHIRC)            FLOWV32.118
C      ****                                                              FLOWV32.119
       PHIRC=0.                                                          FLOWV32.120
       WK2(I)=DMOD(ALON(I)+PHIRC+2.*PI2+PI,PI2)-PI                       FLOWV32.121
       WK3(I)=DMOD(ALON(I)+PHIRC+3.*PI2,PI2)-PI                          FLOWV32.122
      ENDDO                                                              FLOWV32.123
C     ****                                                               FLOWV32.124
C     ****     CALCULATE LONGITUDINAL VARIATION                          FLOWV32.125
C     ****                                                               FLOWV32.126
      DO N=1,7                                                           FLOWV32.127
       DO I=1,IMAXM                                                      FLOWV32.128
        PHIFUN(I)=PHIFUN(I)+.25*(PSI(N)+PSI(N+1)+(PSI(N)                 FLOWV32.129
     1            -PSI(N+1))*DCOS(DMOD(PI*(WK2(I)-PHI(I,N))/(PHI(I,N+1)  FLOWV32.130
     2            -PHI(I,N)),PI2)))*(1.-DSIGN(1.D0,(WK2(I)-PHI(I,N))*    FLOWV32.131
     3            (WK2(I)-PHI(I,N+1))))                                  FLOWV32.132
        PHIFN2(I)=PHIFN2(I)+.25*(PSI(N)+PSI(N+1)+(PSI(N)                 FLOWV32.133
     1            -PSI(N+1))*DCOS(DMOD(PI*(WK3(I)-PHI(I,N))/(PHI(I,N+1)  FLOWV32.134
     2            -PHI(I,N)),PI2)))*(1.-DSIGN(1.D0,(WK3(I)-PHI(I,N))*    FLOWV32.135
     3            (WK3(I)-PHI(I,N+1))))                                  FLOWV32.136
       ENDDO                                                             FLOWV32.137
      ENDDO                                                              FLOWV32.138
C     ****                                                               FLOWV32.139
C     ****     EVALUATE TOTAL POTENTIAL                                  FLOWV32.140
C     ****                                                               FLOWV32.141
      DO I=1,IMAXM                                                       FLOWV32.142
C      ****                                                              FLOWV32.143
C      ****     NEW HEELIS MODEL (PERSONAL COMMUNICATION, 1987)          FLOWV32.144
C      ****                                                              FLOWV32.145
       SSS=(2.*(PCEN(IHEM)-PHIFUN(I))+(PHIFUN(I)-PHIFN2(I))              FLOWV32.146
     1     *0.75)*(COLAT(I)/THETA0(IHEM))**3 +                           FLOWV32.147
     2     (1.5*(PHIFUN(I)+PHIFN2(I))-3.*PCEN(IHEM))*(COLAT(I)/          FLOWV32.148
     3     THETA0(IHEM))**2+0.75*(PHIFUN(I)-PHIFN2(I))*(COLAT(I)/        FLOWV32.149
     4     THETA0(IHEM))+PCEN(IHEM)                                      FLOWV32.150
       TTT=PHIFUN(I)*(DMAX1(DSIN(COLAT(I)),SINTH0)/                      FLOWV32.151
     5     SINTH0)**R1(IHEM)*DEXP(7.*(1.-DMAX1(DSIN(COLAT(I)),SINTHR1)/  FLOWV32.152
     6     SINTHR1) )                                                    FLOWV32.153
       IF((DBLE(IFN(I))-2).EQ.0.)THEN                                    FLOWV32.154
        POTEN(I)=SSS                                                     FLOWV32.155
       ELSE                                                              FLOWV32.156
        POTEN(I)=TTT                                                     FLOWV32.157
       ENDIF                                                             FLOWV32.158
      ENDDO                                                              FLOWV32.159
      RETURN                                                             FLOWV32.160
      END                                                                FLOWV32.161
C                                                                        FLOWV32.162
*DK GRDTRP                                                                GRDTRP.1
      SUBROUTINE GRDTRP(AM,AG,IM,JM,DIM,DJM,LG,LM,IMAXG,JMAXG,IMAXM,      GRDTRP.2
     1JMAXM)                                                              GRDTRP.3
      IMPLICIT REAL*8 (A-H,O-Z)                                           GRDTRP.4
C     ****                                                                GRDTRP.5
C     ****     SUBROUTINE TRANSFORMS SCALAR FIELD GIVEN ON GEOMAGNETIC    GRDTRP.6
C     ****     GRID TO GEOGRAPHIC GRID USING COEFFICIENTS FOR BILINEAR    GRDTRP.7
C     ****     INTERPOLATION PREVIOUSLY GENERATED BY CALL TO GRDSET.      GRDTRP.8
C     ****                                                                GRDTRP.9
C     ****     GEOGRAPHIC GRID:                                           GRDTRP.10
C     ****       IMAXG EQUALLY SPACED LONGITUDE POINTS.  RANGE            GRDTRP.11
C     ****         -PI TO PI.                                             GRDTRP.12
C     ****       JMAXG LATITUDE POINTS, NOT NECESSARILY EQUALLY           GRDTRP.13
C     ****         SPACED.  RANGE -PI/2 TO P/2.                           GRDTRP.14
C     ****                                                                GRDTRP.15
C     ****     GEOMAGNETIC GRID:                                          GRDTRP.16
C     ****       IMAXM EQUALLY SPACED LONGITUDE POINTS.  RANGE            GRDTRP.17
C     ****         -PI TO PI.                                             GRDTRP.18
C     ****       JMAXM LATITUDE POINTS, NOT NECESSARILY EQUALLY           GRDTRP.19
C     ****         SPACED.  RANGE -PI/2 TO P/2.  MUST INCLUDE             GRDTRP.20
C     ****         POLES                                                  GRDTRP.21
C     ****                                                                GRDTRP.22
C     ****     INPUT:                                                     GRDTRP.23
C     ****                                                                GRDTRP.24
C     ****       AM(LM,1) IS 2-DIM SCALAR FIELD GIVEN ON GEOMAGNETIC      GRDTRP.25
C     ****         LATITUDE/LONGITUDE GRID.                               GRDTRP.26
C     ****         NB.  PERIODIC POINT IS REPEATED SO THAT                GRDTRP.27
C     ****         AM(IMAXM+1,JM) = AM(1,JM)                              GRDTRP.28
C     ****         (LM.GE.IMAXM+1)                                        GRDTRP.29
C     ****                                                                GRDTRP.30
C     ****       IMAXG, JMAXG  DIMENSIONS OF GEOGRAPHIC GRID              GRDTRP.31
C     ****                                                                GRDTRP.32
C     ****       IMAXM, JMAXM  DIMENSIONS OF GEOMAGNETIC GRID             GRDTRP.33
C     ****                                                                GRDTRP.34
C     ****       LM IS FIRST DIMENSION OF ARRAY AM INCALLING PROGRAM      GRDTRP.35
C     ****         (LM.GE.IMAXM+1)                                        GRDTRP.36
C     ****                                                                GRDTRP.37
C     ****       IM(LG,1) GIVES ROUNDED DOWN GEOMAGNETIC LONGITUDE        GRDTRP.38
C     ****         INDEX FOR EACH GEOGRAPHIC GRID POINT                   GRDTRP.39
C     ****                                                                GRDTRP.40
C     ****       JM(LG,1) GIVES ROUNDED DOWN GEOMAGNETIC LATITUDE         GRDTRP.41
C     ****         INDEX FOR EACH GEOGRAPHIC GRID POINT                   GRDTRP.42
C     ****                                                                GRDTRP.43
C     ****       DIM(LG,1) IS FRACTIONAL GEOMAGNETIC LONGITUDE INDEX      GRDTRP.44
C     ****         FOR EACH GEOGRAPHIC GRID POINT.  USED FOR BILINEAR     GRDTRP.45
C     ****         INTERPOLATION WITHIN GRID CELL.                        GRDTRP.46
C     ****                                                                GRDTRP.47
C     ****       DJM(LG,1) IS FRACTIONAL GEOMAGNETIC LATITUDE INDEX       GRDTRP.48
C     ****         FOR EACH GEOGRAPHIC GRID POINT.  USED FOR BILINEAR     GRDTRP.49
C     ****         INTERPOLATION WITHIN GRID CELL.                        GRDTRP.50
C     ****                                                                GRDTRP.51
C     ****       LG IS FIRST DIMENSION OF ARRAYS AM,IM,JM,DIM,DJM IN      GRDTRP.52
C     ****         CALLING PROGRAM.                                       GRDTRP.53
C     ****         (LG.GE.IMAXG)                                          GRDTRP.54
C     ****                                                                GRDTRP.55
C     ****       NOTE:  IM, JM, DIM, DJM ARE PRODUCED BY PREVIOUS CALL    GRDTRP.56
C     ****         TO GRDSET.                                             GRDTRP.57
C     ****                                                                GRDTRP.58
C     ****     OUTPUT:                                                    GRDTRP.59
C     ****                                                                GRDTRP.60
C     ****       AG(LG,1) IS 2-DIM SCALAR FIELD TRANSFORMED TO            GRDTRP.61
C     ****         GEOGRAPHIC LATITUDE/LONGITUDE GRID.                    GRDTRP.62
C     ****         (LG.GE.IMAXG)                                          GRDTRP.63
C     ****                                                                GRDTRP.64
C     ****     PERFORM BILINEAR INTERPOLATION                             GRDTRP.65
C     ****                                                                GRDTRP.66
      DIMENSION AM(LM,1),AG(LG,1),IM(LG,1),JM(LG,1),DIM(LG,1),DJM(LG,1)   GRDTRP.67
      DO JG=1,JMAXG                                                       GRDTRP.68
       DO IG=1,IMAXG                                                      GRDTRP.69
        AG(IG,JG)=                                                        GRDTRP.70
     1    AM(IM(IG,JG),JM(IG,JG))*(1.-DIM(IG,JG))*(1.-DJM(IG,JG))+        GRDTRP.71
     2    AM(IM(IG,JG)+1,JM(IG,JG))*DIM(IG,JG)*(1.-DJM(IG,JG))+           GRDTRP.72
     3    AM(IM(IG,JG),JM(IG,JG)+1)*(1.-DIM(IG,JG))*DJM(IG,JG)+           GRDTRP.73
     2    AM(IM(IG,JG)+1,JM(IG,JG)+1)*DIM(IG,JG)*DJM(IG,JG)               GRDTRP.74
       ENDDO                                                              GRDTRP.75
      ENDDO                                                               GRDTRP.76
      RETURN                                                              GRDTRP.77
      END                                                                 GRDTRP.78
C                                                                         GRDTRP.79
*DK HEELIS                                                                HEELIS.1
      SUBROUTINE HEELIS                                                   HEELIS.2
C     ****                                                                HEELIS.3
C     ****     GENERATE AURORAL CIRCLE FIELDS, UI, VI, WI                 HEELIS.4
C     ****                                                                HEELIS.5
      INCLUDE 'param.h'                                                   HEELIS.6
      INCLUDE 'vscrw.h'                                                   HEELIS.7
      INCLUDE 'const.h'                                                   HEELIS.8
      INCLUDE 'index.h'                                                   HEELIS.9
      INCLUDE 'start.h'                                                   HEELIS.10
      INCLUDE 'buffe.h'                                                   HEELIS.11
      INCLUDE 'phycs.h'                                                   HEELIS.12
      INCLUDE 'enest.h'
      COMMON/MAGFLD/BXM(ZIMXP,2),BX(ZIMXP,ZJMX),BXP(ZIMXP,4),             HEELIS.13
     1BY(ZIMXP,ZJMX),BYP(ZIMXP,4),BZ(ZIMXP,ZJMX),BZP(ZIMXP,4),            HEELIS.14
     2BMOD(ZIMXP,ZJMX),BMODP(ZIMXP,2)                                     HEELIS.15
      COMMON/IONDRIFTG/SUG(IMXG,ZKMXP,JMXG),SVG(IMXG,ZKMXP,JMXG),
     1                 SWG(IMXG,ZKMXP,JMXG)
C     ****                                                                HEELIS.16
C     ****     CALCULATE UI AND VI                                        HEELIS.17
C     ****                                                                HEELIS.18
C     ****     SCALE UI AND VI FROM INPUT                                 HEELIS.19
C     ****                                                                HEELIS.20
      NUIK=NUI                                                            HEELIS.21
      NVIK=NVI                                                            HEELIS.22
      NWIK=NWI                                                            HEELIS.23
      DO K=1,KMAXP1
       DO I=1,IMXG                                                        HEELIS.24
        F(I+2,NUIK)=100.*SUG(I,K,J)                                       HEELIS.25
        F(I+2,NVIK)=100.*SVG(I,K,J)                                       HEELIS.25
        F(I+2,NWIK)=100.*SWG(I,K,J)                                       HEELIS.25
       ENDDO
C      ****                                                               HEELIS.35
C      ****     INSERT PERIODIC POINTS                                    HEELIS.36
C      ****                                                               HEELIS.37
       DO I=1,2                                                           HEELIS.38
        F(I,NUIK)=F(I+IMAX,NUIK)                                          HEELIS.39
        F(I,NVIK)=F(I+IMAX,NVIK)                                          HEELIS.39
        F(I,NWIK)=F(I+IMAX,NWIK)                                          HEELIS.39
        F(I+IMAXP2,NUIK)=F(I+2,NUIK)                                      HEELIS.45
        F(I+IMAXP2,NVIK)=F(I+2,NVIK)                                      HEELIS.45
        F(I+IMAXP2,NWIK)=F(I+2,NWIK)                                      HEELIS.45
       ENDDO                                                              HEELIS.51
       NUIK=NUIK+1
       NVIK=NVIK+1
       NWIK=NWIK+1
      ENDDO                                                               HEELIS.28
      USERLA=(DBLE(J-JMAX)-.5)*C(2)+C(110)/2.                             HEELIS.29
      ISHUNK=1                                                            HEELIS.30
      IF(DABS(USERLA).LT.C(110)/6.)ISHUNK=0                               HEELIS.31
      IF(ISHUNK.EQ.0)GO TO 24                                             HEELIS.32
      JJ=((2*J-JMAX-1)/IABS(2*J-JMAX-1)+3)/2                              HEELIS.33
      CALL FLOWION(T5(3),T6(3),T7(3),T8(3),T9(3),T3(3))                   HEELIS.34
C     ****                                                                HEELIS.35
C     ****     INSERT PERIODIC POINTS                                     HEELIS.36
C     ****                                                                HEELIS.37
      DO I=1,2                                                            HEELIS.38
       T3(I)=T3(I+IMAX)                                                   HEELIS.39
       T5(I)=T5(I+IMAX)                                                   HEELIS.40
       T6(I)=T6(I+IMAX)                                                   HEELIS.41
       T7(I)=T7(I+IMAX)                                                   HEELIS.42
       T8(I)=T8(I+IMAX)                                                   HEELIS.43
       T9(I)=T9(I+IMAX)                                                   HEELIS.44
       T3(I+IMAXP2)=T3(I+2)                                               HEELIS.45
       T5(I+IMAXP2)=T5(I+2)                                               HEELIS.46
       T6(I+IMAXP2)=T6(I+2)                                               HEELIS.47
       T7(I+IMAXP2)=T7(I+2)                                               HEELIS.48
       T8(I+IMAXP2)=T8(I+2)                                               HEELIS.49
       T9(I+IMAXP2)=T9(I+2)                                               HEELIS.50
      ENDDO                                                               HEELIS.51
      GO TO 25                                                            HEELIS.52
   24 CONTINUE                                                            HEELIS.53
      DO I=1,LEN1                                                         HEELIS.54
       T3(I)=0.                                                           HEELIS.55
       T5(I)=0.                                                           HEELIS.56
       T6(I)=0.                                                           HEELIS.57
       T7(I)=0.                                                           HEELIS.58
       T8(I)=0.                                                           HEELIS.59
       T9(I)=0.                                                           HEELIS.60
      ENDDO                                                               HEELIS.61
   25 CONTINUE                                                            HEELIS.62
      RETURN                                                              HEELIS.63
      END                                                                 HEELIS.64
C                                                                         HEELIS.65
*DK HTRPEX                                                                HTRPEX.1
      SUBROUTINE HTRPEX(COEFF,ARRAY,IMX0,JMX0,SYM)                        HTRPEX.2
      IMPLICIT REAL*8 (A-H,O-Z)                                           HTRPEX.3
C     ****                                                                HTRPEX.4
C     ****     PERFORM HALF WAY INTERPOLATION ON ARRAY COEFF AND EXTEND   HTRPEX.5
C     ****       OVER 16 GRID POINTS .  RESULT RETURNED IN ARRAY          HTRPEX.6
C     ****                                                                HTRPEX.7
C     ****     INPUT:                                                     HTRPEX.8
C     ****       COEFF((IMX0+1)/2,(JMX0+1)/2)  CONTAINS ARRAY OF ONE OF   HTRPEX.9
C     ****         PDE COEFFICIENTS.                                      HTRPEX.10
C     ****       SYM =1.  IF COEFFICIENT SYMMETRIC W.R.T. EQUATOR         HTRPEX.11
C     ****           = -1.  FOR ANTISYMMETRY                              HTRPEX.12
C     ****                                                                HTRPEX.13
C     ****     OUTPUT:                                                    HTRPEX.14
C     ****       ARRAY(-15:IMX0+16,JMX0)  INTERPOLATED AND                HTRPEX.15
C     ****         EXTENDED COEFFICIENT RETURNED IN THIS ARRAY.           HTRPEX.16
C     ****                                                                HTRPEX.17
      DIMENSION COEFF(IMX0,JMX0),ARRAY(IMX0+32,1)                         HTRPEX.18
C     ****                                                                HTRPEX.19
C     ****     COPY COEFF INTO POSITIONS IN ARRAY                         HTRPEX.20
C     ****                                                                HTRPEX.21
      DO J=1,JMX0                                                         HTRPEX.22
       JJ=JMX0-J+1                                                        HTRPEX.23
       DO I=1,IMX0                                                        HTRPEX.24
        ARRAY(I,J)=SYM*COEFF(I,JJ)                                        HTRPEX.25
       ENDDO                                                              HTRPEX.26
      ENDDO                                                               HTRPEX.27
C     ****                                                                HTRPEX.28
C     ****     NOW EXTEND OVER 32 GRID SPACES TO ALLOW FOR A TOTAL OF     HTRPEX.29
C     ****       FIVE GRID LEVELS                                         HTRPEX.30
C     ****                                                                HTRPEX.31
      DO I=1,16                                                           HTRPEX.32
       DO J=1,JMX0                                                        HTRPEX.33
        ARRAY(1-I,J)=ARRAY(IMX0-I,J)                                      HTRPEX.34
        ARRAY(IMX0+I,J)=ARRAY(1+I,J)                                      HTRPEX.35
       ENDDO                                                              HTRPEX.36
      ENDDO                                                               HTRPEX.37
      RETURN                                                              HTRPEX.38
      END                                                                 HTRPEX.39
C                                                                         HTRPEX.40
*DK INTGRLS                                                              INTGRLS.1
      SUBROUTINE INTGRL2(JM)                                             INTGRLS.2
C     ****                                                               INTGRLS.3
C     ****     PERFORMS APPROXIMATE FIELD LINE QUADRATURES AT ALL        INTGRLS.4
C     ****       GEOMAGNETIC GRID POINTS                                 INTGRLS.5
C     ****                                                               INTGRLS.6
      INCLUDE "param.h"                                                  INTGRLS.7
      INCLUDE 'const.h'
      INCLUDE 'epara.h'                                                  INTGRLS.8
      INCLUDE 'econs.h'                                                  INTGRLS.9
      INCLUDE 'ecoef.h'                                                  INTGRLS.10
C     ****                                                               INTGRLS.11
C     ****     FIELDS TO BE TRANSFORMED TO GEOMAGNETIC SPACE             INTGRLS.12
C     ****       EXTRAPOLATED DOWN TO 90 KM                              INTGRLS.13
C     ****                                                               INTGRLS.14
       COMMON/TRNSFR/ZZ(IMAXGP,0:JMAXGP,2),                              INTGRLS.15
     1               SSIGMA1(IMAXGP,0:JMAXGP),                           INTGRLS.16
     2               SSIGMA2(IMAXGP,0:JMAXGP),                           INTGRLS.17
     3               ADOTV(IMAXGP,0:JMAXGP,2),                           INTGRLS.18
     4               ADOTA(IMAXGP,0:JMAXGP,2),                           INTGRLS.19
     5               A1DTA2(IMAXGP,0:JMAXGP),                            INTGRLS.20
     6               SINI(IMAXGP,0:JMAXGP),                              INTGRLS.21
     7               BE3(IMAXGP,0:JMAXGP),                               INTGRLS.22
C     ****                                                               INTGRLS.23
C     ****     THESE SAME QUANTITIES IN GEOMAGNETIC SPACE                INTGRLS.24
C     ****       (ONE LATITUDE AT A TIME)                                INTGRLS.25
C     ****                                                               INTGRLS.26
     1              ZM(IMAXMP,2),                                        INTGRLS.27
     2              SIGMA1M(IMAXMP),                                     INTGRLS.28
     3              SIGMA2M(IMAXMP),                                     INTGRLS.29
     4              ADOTVM(IMAXMP,2),                                    INTGRLS.30
     5              ADOTAM(IMAXMP,2),                                    INTGRLS.31
     6              A1A2M(IMAXMP),                                       INTGRLS.32
     7              SINIAM(IMAXMP),                                      INTGRLS.33
     8              BMODM(IMAXMP),                                       INTGRLS.34
C     ****                                                               INTGRLS.35
C     ****     99/2/23 Note: be3,BMODM are at 90 km, whereas BMOD        INTGRLS.36
C     ****     (in fieldz.h) is at 130 km.                               INTGRLS.37
C     ****                                                               INTGRLS.38
C     ****     Equatorial quantities                                     INTGRLS.39
C     ****                                                               INTGRLS.40
     1              SIGMA1E(IMAXMP),                                     INTGRLS.41
     2              SIGMA2E(IMAXMP),                                     INTGRLS.42
     3              ADOTVE(IMAXMP,2),                                    INTGRLS.43
C     ****                                                               INTGRLS.44
C     ****     WORK SPACE                                                INTGRLS.45
C     ****                                                               INTGRLS.46
     1              TINT1(JMAXM),TINT2(JMAXM),TINT3(JMAXM)               INTGRLS.47
C     ****                                                               INTGRLS.48
C     ****     QUANTITIES needed for field-line integrations             INTGRLS.49
C     ****                                                               INTGRLS.50
      DIMENSION AAM(IMAXMP),RRM(IMAXMP,2),                               INTGRLS.51
     1	        RTRAMRM(IMAXMP,2)                                        INTGRLS.52
      DIMENSION HTFUNC(IMAXMP),HTFUNC2(IMAXMP)                           INTGRLS.53
      DATA Z0/90.D5/                                                     INTGRLS.54
      SINLM=DSIN(YLATM(JM))                                              INTGRLS.55
      CLM2=1.-SINLM*SINLM                                                INTGRLS.56
      RA=R0/CLM2                                                         INTGRLS.57
      SQOMRRA=DSQRT(1.D0-R0/RA)                                          INTGRLS.58
      SQRRA=DSQRT(R0/RA)                                                 INTGRLS.59
      AFAC=2.*DSQRT(RA-R0)                                               INTGRLS.60
      HTFAC=DSQRT(RA-.75*R0)                                             INTGRLS.61
      DO I=1,IMAXM                                                       INTGRLS.62
       AAM(I)=AFAC/DABS(SINIAM(I))                                       INTGRLS.63
      ENDDO                                                              INTGRLS.64
C     ****                                                               INTGRLS.65
C     ****     3-D FIELDS                                                INTGRLS.66
C     ****                                                               INTGRLS.67
      DO K=1,2                                                           INTGRLS.68
       DO I=1,IMAXM                                                      INTGRLS.69
C       ****                                                             INTGRLS.70
C       ****     RR = R0+Z-Z0                                            INTGRLS.71
C       ****                                                             INTGRLS.72
        RRM(I,K)=R0+ZM(I,K)-Z0                                           INTGRLS.73
C       ****                                                             INTGRLS.74
C       ****     RTRAMR = RA-R IF +IVE, ZERO OTHERWISE                   INTGRLS.75
C       ****                                                             INTGRLS.76
        RTRAMRM(I,K)=DMAX1(0.D0,RA-RRM(I,K))                             INTGRLS.77
C       ****                                                             INTGRLS.78
C       ****     RTRAMR = SQRT(RA-R                                      INTGRLS.79
C       ****                                                             INTGRLS.80
        RTRAMRM(I,K)=DSQRT(RTRAMRM(I,K))                                 INTGRLS.81
       ENDDO                                                             INTGRLS.82
      ENDDO                                                              INTGRLS.83
C     ****                                                               INTGRLS.84
C     ****     INTERPOLATE TO MIDDLE OF HORIZONTAL SLICES                INTGRLS.85
C     ****                                                               INTGRLS.86
      DO I=1,IMAXM                                                       INTGRLS.88
       RRM(I,1)=.5*(RRM(I,1)+RRM(I,2))                                   INTGRLS.89
       RTRAMRM(I,1)=RTRAMRM(I,1)-RTRAMRM(I,2)                            INTGRLS.90
C      ****                                                              INTGRLS.91
C      **** htfunc = factor by which to multiply AAM(I)*d(sqrt(ra-r))    INTGRLS.92
C      ****                                                              INTGRLS.93
       HTFUNC(I)=DSQRT(RA-.75D0*RRM(I,K))/HTFAC                          INTGRLS.94
       HTFUNC2(I)=HTFUNC(I)**2                                           INTGRLS.95
      ENDDO                                                              INTGRLS.96
C     ****                                                               INTGRLS.98
C     ****     COMPUTE INTEGRALS                                         INTGRLS.99
C     ****                                                               INTGRLS.100
      DO I=1,IMAXM                                                       INTGRLS.101
       ZIGM11(I,JM)=0.                                                   INTGRLS.102
       ZIGM22(I,JM)=0.                                                   INTGRLS.103
       ZIGM2(I,JM)=0.                                                    INTGRLS.104
       ZIGMC(I,JM)=0.                                                    INTGRLS.105
       RIM(I,JM,1)=0.                                                    INTGRLS.106
       RIM(I,JM,2)=0.                                                    INTGRLS.107
      ENDDO                                                              INTGRLS.108
      DO I=1,IMAXM                                                       INTGRLS.110
       RORA=DMIN1(1.D0,RRM(I,1)/RA)                                      INTGRLS.111
       DEL=(SQOMRRA*DSQRT(RORA)-SQRRA*DSQRT(1.D0-RORA))/                 INTGRLS.112
     1     DABS(YLATM(JM))                                               INTGRLS.113
       OMDEL=1.-DEL                                                      INTGRLS.114
C      ****     Interpolate conductivities and winds in latitude along   INTGRLS.115
C      ****     field line, assuming linear variation between foot of    INTGRLS.116
C      ****     field line and magnetic equator. (For field lines other  INTGRLS.117
C      ****     than those near the magnetic equator, del is nearly      INTGRLS.118
C      ****     zero, so that the interpolated values are essentially    INTGRLS.119
C      ****     the values for the latitude of the foot of the field     INTGRLS.120
C      ****     line; inaccuracy of the assumption of linear variation   INTGRLS.121
C      ****     is thus unimportant for these field lines.)              INTGRLS.122
C      ****                                                              INTGRLS.123
       SIG1=OMDEL*SIGMA1M(I)+DEL*SIGMA1E(I)                              INTGRLS.124
       SIG2=OMDEL*SIGMA2M(I)+DEL*SIGMA2E(I)                              INTGRLS.125
       UE1=OMDEL*ADOTVM(I,1)+DEL*ADOTVE(I,1)                             INTGRLS.126
       UE2=OMDEL*ADOTVM(I,2)+DEL*ADOTVE(I,2)                             INTGRLS.127
       ZIGM11(I,JM)=ZIGM11(I,JM)+SIG1*RTRAMRM(I,1)                         INTGRLS.128
       ZIGM22(I,JM)=ZIGM22(I,JM)+SIG1*RTRAMRM(I,1)*HTFUNC2(I)              INTGRLS.129
C      ****                                                              INTGRLS.130
C      ****      Following uses Richmond convention for sign of ZIGMC    INTGRLS.131
C      ****     (= -ZIGMC of Ridley)                                     INTGRLS.132
C      ****                                                              INTGRLS.133
       ZIGMC(I,JM)=ZIGMC(I,JM)+SIG1*RTRAMRM(I,1)*HTFUNC(I)                 INTGRLS.134
       ZIGM2(I,JM)=ZIGM2(I,JM)+SIG2*RTRAMRM(I,1)*HTFUNC(I)                 INTGRLS.135
       RIM(I,JM,1)=RIM(I,JM,1)+(SIG1*ADOTAM(I,1)*UE2+                    INTGRLS.136
     2             (SIG2-SIG1*A1A2M(I))*HTFUNC(I)*UE1)*RTRAMRM(I,1)        INTGRLS.137
       RIM(I,JM,2)=RIM(I,JM,2)+                                          INTGRLS.138
     1             (SIG1*ADOTAM(I,2)*HTFUNC2(I)*UE1-                     INTGRLS.139
     2             (SIG2+SIG1*A1A2M(I))*HTFUNC(I)*UE2)*RTRAMRM(I,1)        INTGRLS.140
      ENDDO                                                              INTGRLS.141
      DO I=1,IMAXM                                                       INTGRLS.143
C      ****                                                              INTGRLS.144
C      ****     ZIGM's are in S, RIM's are in A/m                        INTGRLS.145
C      ****                                                              INTGRLS.146
C      ****     At this point, ZIGM11 is int[sigP*d1**2/D], i.e.,        INTGRLS.147
C      ****     Sigma_(phi phi)/abs(sin Im)                              INTGRLS.148
C      ****                                                              INTGRLS.149
       ZIGM11(I,JM)=1.D-2*ZIGM11(I,JM)*AAM(I)*ADOTAM(I,1)                INTGRLS.150
C      ****                                                              INTGRLS.151
C      ****     At this point, ZIGM22 is int[sigP*d2**2/D], i.e.,        INTGRLS.152
C      ****     Sigma_(lam lam)*abs(sin Im)                              INTGRLS.153
C      ****                                                              INTGRLS.154
       ZIGM22(I,JM)=1.D-2*ZIGM22(I,JM)*AAM(I)*ADOTAM(I,2)                INTGRLS.155
       ZIGMC(I,JM)=1.D-2*ZIGMC(I,JM)*AAM(I)*A1A2M(I)                     INTGRLS.156
       ZIGM2(I,JM)=1.D-2*ZIGM2(I,JM)*AAM(I)                              INTGRLS.157
C      ****                                                              INTGRLS.158
C      ****      At this point, RIM1 is K_(m phi)^D/abs(sin Im)          INTGRLS.159
C      ****                                                              INTGRLS.160
       RIM(I,JM,1)=1.D-2*RIM(I,JM,1)*AAM(I)*BMODM(I)                     INTGRLS.161
C      ****                                                              INTGRLS.162
C      ****     At this point, RIM2 is K_(m lam)^D                       INTGRLS.163
C      ****                                                              INTGRLS.164
       RIM(I,JM,2)=1.D-2*RIM(I,JM,2)*AAM(I)*BMODM(I)                     INTGRLS.165
      ENDDO                                                              INTGRLS.166
      RETURN                                                             INTGRLS.167
      END                                                                INTGRLS.168
*DK LAMDAS                                                                LAMDAS.1
      SUBROUTINE LAMDAS                                                   LAMDAS.2
C     ****                                                                LAMDAS.3
C     ****     COMPUTE ION DRAG COEFFICIENTS                              LAMDAS.4
C     ****                                                                LAMDAS.5
C     ****     MAJOR SPECIES AND IONS ARE NUMBERED THUS                   LAMDAS.6
C     ****       O2  -  1            O2+  -  1                            LAMDAS.7
C     ****       O   -  2            O+   -  2                            LAMDAS.8
C     ****       N2  -  3            N2+  -  3                            LAMDAS.9
C     ****                                                                LAMDAS.10
      INCLUDE 'param.h'                                                   LAMDAS.11
      INCLUDE 'const.h'                                                   LAMDAS.12
      INCLUDE 'index.h'                                                   LAMDAS.13
      INCLUDE 'buffe.h'                                                   LAMDAS.14
      INCLUDE 'phycs.h'                                                   LAMDAS.15
      INCLUDE 'edyna.h'                                                   LAMDAS.16
      COMMON/TRGM/RLATM(ZIMXP,ZJMX),RLONM(ZIMXP,ZJMX),                    LAMDAS.17
     1      DIPMAG(ZIMXP,ZJMX),DECMAG(ZIMXP,ZJMX),SNDEC(ZIMXP,ZJMX),      LAMDAS.18
     2      CSDEC(ZIMXP,ZJMX),SN2DEC(ZIMXP,ZJMX),SNCSDC(ZIMXP,ZJMX),      LAMDAS.19
     3      RLATMP(ZJMX),RLONMP(ZJMX),DLONS(ZJMX)                         LAMDAS.20
      COMMON/MAGFLD/BXM(ZIMXP,2),BX(ZIMXP,ZJMX),BXP(ZIMXP,4),             LAMDAS.21
     1              BY(ZIMXP,ZJMX),BYP(ZIMXP,4),BZ(ZIMXP,ZJMX),           LAMDAS.22
     2              BZP(ZIMXP,4),BMOD(ZIMXP,ZJMX),BMODP(ZIMXP,2)          LAMDAS.23
      DIMENSION RM(ZIMXP,ZKMXP,3,3),V1(ZIMXP,ZKMXP),V2(ZIMXP,ZKMXP),      LAMDAS.24
     1          V3(ZIMXP,ZKMXP),V4(ZIMXP,ZKMXP),V5(ZIMXP,ZKMXP),          LAMDAS.25
     2          V6(ZIMXP,ZKMXP),T1(ZIMXP),T2(ZIMXP),T3(ZIMXP),T4(ZIMXP)   LAMDAS.26
      DIMENSION RMS(3)                                                    LAMDAS.27
      DATA RMS/32.,16.,30./,QE/1.602D-19/                                 LAMDAS.28
C     ****                                                                LAMDAS.29
C     ****     T1 = BGAUSS                                                LAMDAS.30
C     ****     T2 = 9.6489D3 * BGAUS    ( FOR OMEGA(I) )                  LAMDAS.31
C     ****     T3 = 1.7588028D7 * BGAUSS  ( FOR OMEGA(E) )                LAMDAS.32
C     ****     T4 = SIN(DIP)                                              LAMDAS.33
C     ****                                                                LAMDAS.34
      DO I=1,LEN1                                                         LAMDAS.35
       T1(I)=BMOD(I,J)                                                    LAMDAS.36
       T2(I)=9.6489D3*T1(I)                                               LAMDAS.37
       T3(I)=1.7588028D7*T1(I)                                            LAMDAS.38
       T4(I)=DSIN(DIPMAG(I,J))                                            LAMDAS.39
      ENDDO                                                               LAMDAS.40
C     ****                                                                LAMDAS.41
C     ****     V1 = PSI(O2)/M(O2), V2 = PSI(O)/M(O), V3 = PSI2)/M(N2)     LAMDAS.42
C     ****      V6 = TR = 0.5*(TI+TN)                 AT (K+1/2)          LAMDAS.43
C     ****                                                                LAMDAS.44
      NPS1K=NJ+NPS                                                        LAMDAS.45
      NPS2K=NJ+NPS2                                                       LAMDAS.46
      NTIK=NJ+NTI                                                         LAMDAS.47
      NTK=NJ+NT                                                           LAMDAS.48
      DO I=1,LEN2                                                         LAMDAS.49
       V1(I,1)=F(I,NPS1K)/RMASS(1)                                        LAMDAS.50
       V2(I,1)=F(I,NPS2K)/RMASS(2)                                        LAMDAS.51
       V3(I,1)=(1.-F(I,NPS1K)-F(I,NPS2K))/RMASS(3)                        LAMDAS.52
       V6(I,1)=0.5*(F(I,NTIK)+F(I,NTK))                                   LAMDAS.53
      ENDDO                                                               LAMDAS.54
C     ****                                                                LAMDAS.55
C     ****     RM = NEW(L,M)                                              LAMDAS.56
C     ****                                                                LAMDAS.57
C     ****     GENERATE NUMERICAL FACTORS                                 LAMDAS.58
C     ****                                                                LAMDAS.59
      CALL NEW(V6,RM)                                                     LAMDAS.60
C     ****                                                                LAMDAS.61
C     ****     MULTIPLY BY MAJOR SPECIES NUMBER DENSITIES                 LAMDAS.62
C     ****                                                                LAMDAS.63
      DO L=1,3                                                            LAMDAS.64
       DO I=1,LEN2                                                        LAMDAS.65
        RM(I,1,L,1)=RM(I,1,L,1)*V1(I,1)                                   LAMDAS.66
        RM(I,1,L,2)=RM(I,1,L,2)*V2(I,1)                                   LAMDAS.67
        RM(I,1,L,3)=RM(I,1,L,3)*V3(I,1)                                   LAMDAS.68
       ENDDO                                                              LAMDAS.69
      ENDDO                                                               LAMDAS.70
C     ****                                                                LAMDAS.71
C     ****     V1(L) = SUM(M=1,3) RM1(L,M)    ( K+1/2 )                   LAMDAS.72
C     ****                                                                LAMDAS.73
      DO I=1,LEN2                                                         LAMDAS.74
       V1(I,1)=0.                                                         LAMDAS.75
       V2(I,1)=0.                                                         LAMDAS.76
       V3(I,1)=0.                                                         LAMDAS.77
      ENDDO                                                               LAMDAS.78
      DO M=1,3                                                            LAMDAS.79
       DO I=1,LEN2                                                        LAMDAS.80
        V1(I,1)=V1(I,1)+RM(I,1,1,M)                                       LAMDAS.81
        V2(I,1)=V2(I,1)+RM(I,1,2,M)                                       LAMDAS.82
        V3(I,1)=V3(I,1)+RM(I,1,3,M)                                       LAMDAS.83
       ENDDO                                                              LAMDAS.84
      ENDDO                                                               LAMDAS.85
C     ****                                                                LAMDAS.86
C     ****     RM(1) = N(O2+), RM(2) = N(O+), RM(3) = N(NO+)              LAMDAS.87
C     ****                      ( AT (K+1/2) )                            LAMDAS.88
C     ****                                                                LAMDAS.89
      NO2PK=NJ+NO2P                                                       LAMDAS.90
      NOPK=NJ+NOP                                                         LAMDAS.91
      NNOPK=NNOP                                                          LAMDAS.92
      DO I=1,LEN2                                                         LAMDAS.93
       RM(I,1,1,1)=F(I,NO2PK)                                             LAMDAS.94
       RM(I,1,2,1)=F(I,NOPK)                                              LAMDAS.95
       RM(I,1,3,1)=F(I,NNOPK)                                             LAMDAS.96
      ENDDO                                                               LAMDAS.97
C     ****                                                                LAMDAS.98
C     ****     V6 = N * MBAR            ( K+1/2 )                         LAMDAS.99
C     ****                                                                LAMDAS.100
      NTK=NJ+NT-1                                                         LAMDAS.101
      NMSK=NJ+NMS-1                                                       LAMDAS.102
      DO K=1,KMAX                                                         LAMDAS.103
       NTK=NTK+1                                                          LAMDAS.104
       NMSK=NMSK+1                                                        LAMDAS.105
       DO I=1,LEN1                                                        LAMDAS.106
        V6(I,K)=C(81)*EXPS(K)*.5*(F(I,NMSK)+F(I,NMSK+1))/(C(84)*          LAMDAS.107
     1            F(I,NTK))                                               LAMDAS.108
       ENDDO                                                              LAMDAS.109
      ENDDO                                                               LAMDAS.110
C     ****                                                                LAMDAS.111
C     ****     V1(L) = V1(L)*(N*MBAR)/OMEGAI(L)    ( L=1,3 )  ( K=1/2 )   LAMDAS.112
C     ****                                                                LAMDAS.113
      DO K=1,KMAX                                                         LAMDAS.114
       DO I=1,LEN1                                                        LAMDAS.115
        V1(I,K)=V1(I,K)*V6(I,K)*RMS(1)/T2(I)                              LAMDAS.116
        V2(I,K)=V2(I,K)*V6(I,K)*RMS(2)/T2(I)                              LAMDAS.117
        V3(I,K)=V3(I,K)*V6(I,K)*RMS(3)/T2(I)                              LAMDAS.118
       ENDDO                                                              LAMDAS.119
      ENDDO                                                               LAMDAS.120
C     ****                                                                LAMDAS.121
C     ****     V4 = NU/OMEGA (E)                                          LAMDAS.122
C     ****                                                                LAMDAS.123
      NPS1K=NJ+NPS                                                        LAMDAS.124
      NPS2K=NJ+NPS2                                                       LAMDAS.125
      NTEK=NJ+NTE                                                         LAMDAS.126
      DO I=1,LEN2                                                         LAMDAS.127
       V5(I,1)=DSQRT(F(I,NTEK))                                           LAMDAS.128
       V4(I,1)=V6(I,1)*(2.33D-11*(1.-F(I,NPS1K)-F(I,NPS2K))/RMASS(3)*     LAMDAS.129
     1         F(I,NTEK)*(1.-1.21D-4*F(I,NTEK))+1.82D-10*F(I,NPS1K)/      LAMDAS.130
     2         RMASS(1)*V5(I,1)*(1.+3.6D-2*V5(I,1))+8.9D-11*              LAMDAS.131
     3         F(I,NPS2K)/RMASS(2)*V5(I,1)*(1.+5.7D-4*F(I,NTEK)))         LAMDAS.132
      ENDDO                                                               LAMDAS.133
C     ****                                                                LAMDAS.134
C     ****     V4 = NU(EN)/OMEGAE                                         LAMDAS.135
C     ****                                                                LAMDAS.136
      DO K=1,KMAX                                                         LAMDAS.137
       DO I=1,LEN1                                                        LAMDAS.138
        V4(I,K)=V4(I,K)/T3(I)                                             LAMDAS.139
       ENDDO                                                              LAMDAS.140
      ENDDO                                                               LAMDAS.141
C     ****                                                                LAMDAS.142
C     ****     RM2 = PEDERSEN FACTORS                                     LAMDAS.143
C     ****     RM3 = HALL FACTORS                                         LAMDAS.145
C     ****                                                                LAMDAS.146
      DO I=1,LEN2                                                         LAMDAS.147
       RM(I,1,1,2)=V4(I,1)/(1.+V4(I,1)**2)+V1(I,1)/(1.+V1(I,1)**2)        LAMDAS.148
       RM(I,1,2,2)=V4(I,1)/(1.+V4(I,1)**2)+V2(I,1)/(1.+V2(I,1)**2)        LAMDAS.149
       RM(I,1,3,2)=V4(I,1)/(1.+V4(I,1)**2)+V3(I,1)/(1.+V3(I,1)**2)        LAMDAS.140
       RM(I,1,1,3)=1./(1.+V4(I,1)**2)-1./(1+V1(I,1)**2)                   LAMDAS.150
       RM(I,1,2,3)=1./(1.+V4(I,1)**2)-1./(1+V2(I,1)**2)                   LAMDAS.151
       RM(I,1,3,3)=1./(1.+V4(I,1)**2)-1./(1+V3(I,1)**2)                   LAMDAS.152
      ENDDO                                                               LAMDAS.153
C     ****                                                                LAMDAS.154
C     ****     V1 = SUM(L=1,3) PEDERSEN FACTORS * ION DENSITIES           LAMDAS.155
C     ****     V2 = SUM(L=1,3) HALL FACTORS * ION DENSITIES               LAMDAS.156
C     ****                                                                LAMDAS.157
      DO I=1,LEN2                                                         LAMDAS.158
       V1(I,1)=0.                                                         LAMDAS.159
       V2(I,1)=0.                                                         LAMDAS.160
      ENDDO                                                               LAMDAS.161
      DO L=1,3                                                            LAMDAS.162
       DO I=1,LEN2                                                        LAMDAS.163
        V1(I,1)=V1(I,1)+RM(I,1,L,2)*RM(I,1,L,1)                           LAMDAS.164
        V2(I,1)=V2(I,1)+RM(I,1,L,3)*RM(I,1,L,1)                           LAMDAS.165
       ENDDO                                                              LAMDAS.166
      ENDDO                                                               LAMDAS.167
C     ****                                                                LAMDAS.168
C     ****     SET LXX, LYY, LXY                                          LAMDAS.169
C     ****                                                                LAMDAS.170
      DO K=1,KMAX                                                         LAMDAS.171
       DO I=1,LEN1                                                        LAMDAS.172
        V3(I,K)=V1(I,K)*T1(I)*1.D-1*QE*C(85)/V6(I,K)                      LAMDAS.173
        V4(I,K)=V3(I,K)*T4(I)**2                                          LAMDAS.174
        V5(I,K)=V2(I,K)*T1(I)*1.D-1*QE*C(85)/V6(I,K)                      LAMDAS.175
       ENDDO                                                              LAMDAS.176
      ENDDO                                                               LAMDAS.177
      NLXXK=NLXX+1                                                        LAMDAS.178
      NLYYK=NLYY+1                                                        LAMDAS.179
      NLXYK=NLXY+1                                                        LAMDAS.180
      DO I=1,LEN2-LEN1                                                    LAMDAS.181
       F(I,NLXXK)=DSQRT(DABS(V3(I,1)*V3(I,2)))                            LAMDAS.182
       F(I,NLYYK)=DSQRT(DABS(V4(I,1)*V4(I,2)))                            LAMDAS.183
       F(I,NLXYK)=DSQRT(DABS(V5(I,1)*V5(I,2)))                            LAMDAS.184
      ENDDO                                                               LAMDAS.185
C     ****                                                                LAMDAS.186
C     ****     BOUNDARIES                                                 LAMDAS.187
C     ****                                                                LAMDAS.188
      NLXXK=NLXX                                                          LAMDAS.189
      NLYYK=NLYY                                                          LAMDAS.190
      NLXYK=NLXY                                                          LAMDAS.191
      DO I=1,LEN1                                                         LAMDAS.192
       F(I,NLXXK)=DSQRT(V3(I,1)**3/V3(I,2))                               LAMDAS.193
       F(I,NLYYK)=DSQRT(V4(I,1)**3/V4(I,2))                               LAMDAS.194
       F(I,NLXYK)=DSQRT(V5(I,1)**3/V5(I,2))                               LAMDAS.195
       F(I,NLXXK+KMAX)=DSQRT(DABS(V3(I,KMAX)**3/V3(I,KMAXM1)))            LAMDAS.196
       F(I,NLYYK+KMAX)=DSQRT(DABS(V4(I,KMAX)**3/V4(I,KMAXM1)))            LAMDAS.197
       F(I,NLXYK+KMAX)=DSQRT(DABS(V5(I,KMAX)**3/V5(I,KMAXM1)))            LAMDAS.198
      ENDDO                                                               LAMDAS.199
      NLXYK=NLXY-1                                                        LAMDAS.200
      DO K=1,KMAXP1                                                       LAMDAS.201
       NLXYK=NLXYK+1                                                      LAMDAS.202
       DO I=1,IMAXP4                                                      LAMDAS.203
        F(I,NLXYK)=F(I,NLXYK)*T4(I)                                       LAMDAS.204
       ENDDO                                                              LAMDAS.205
      ENDDO                                                               LAMDAS.206
C     ****                                                                LAMDAS.207
C     ****     COPY SIGMA1, SIGMA2, Z, H, U, V AND W TO COMMON/DYNAMO/    LAMDAS.208
C     ****                                                                LAMDAS.209
      NZK=NJ+NZ-1                                                         LAMDAS.210
      NTK=NJ+NT-1                                                         LAMDAS.211
      NMSK=NJ+NMS-1                                                       LAMDAS.212
      NUK=NJ+NU-1                                                         LAMDAS.213
      NVK=NJ+NV-1                                                         LAMDAS.214
      NWK=NJ+NW-1                                                         LAMDAS.215
      DO K=1,KMAX                                                         LAMDAS.216
       NZK=NZK+1                                                          LAMDAS.217
       NTK=NTK+1                                                          LAMDAS.218
       NMSK=NMSK+1                                                        LAMDAS.219
       NUK=NUK+1                                                          LAMDAS.220
       NVK=NVK+1                                                          LAMDAS.221
       NWK=NWK+1                                                          LAMDAS.222
       DO I=1,IMAX+1                                                      LAMDAS.223
        SIGMA1(I,J,K)=V1(I+2,K)*1.d10*QE/T1(I+2)                          LAMDAS.224
        SIGMA2(I,J,K)=V2(I+2,K)*1.d10*QE/T1(I+2)                          LAMDAS.225
        Z(I,J,K)=F(I+2,NZK)                                               LAMDAS.226
        H=C(57)*F(I+2,NTK)/(.5d0*(F(I+2,NMSK)+F(I+2,NMSK+1))*C(54))       LAMDAS.227
        U(I,J,K)=F(I+2,NUK)                                               LAMDAS.229
        V(I,J,K)=F(I+2,NVK)                                               LAMDAS.230
        W(I,J,K)=.5*(F(I+2,NWK)+F(I+2,NWK+1))*H                           LAMDAS.231
       ENDDO                                                              LAMDAS.232
      ENDDO                                                               LAMDAS.233
      RETURN                                                              LAMDAS.238
      END                                                                 LAMDAS.239
C                                                                         LAMDAS.240
*DK POTM                                                                    POTM.1
      SUBROUTINE POTM                                                       POTM.2
C     ****                                                                  POTM.3
C     ****     CALCULATE HEELIS POTENTIAL IN GEOMAGNETIC COORDINATES        POTM.4
C     ****                                                                  POTM.5
      include "param.h"                                                     POTM.6
      INCLUDE 'epara.h'                                                     POTM.7
      include "const.h"                                                     POTM.8
      include "econs.h"                                                     POTM.9
      include "edphi.h"                                                     POTM.10
C     ****                                                                  POTM.11
C     ****     WORK SPACE                                                   POTM.12
C     ****                                                                  POTM.13
      COMMON/POTE/ DLAT(IMAXM),DLON(IMAXM),PSI(8),PHI(IMAXM,8),             POTM.14
     1  IFN(IMAXM),PHIFUN(IMAXM),COLAT(IMAXM),ALON(IMAXM),                  POTM.15
     2  SINLAT(IMAXM),COSLAT(IMAXM),SINLON(IMAXM),COSLON(IMAXM),            POTM.16
     3  WK1(IMAXM),WK2(IMAXM),WK3(IMAXM),PHIFN2(IMAXM)                      POTM.17
C     ****                                                                  POTM.18
C     ****     /TRGM/ CONTAINS MAGNETIC INFO                                POTM.19
C     ****                                                                  POTM.20
      COMMON/TRGM/RLATM(ZIMXP,ZJMX),RLONM(ZIMXP,ZJMX),                      POTM.21
     1  DIPMAG(ZIMXP,ZJMX),DECMAG(ZIMXP,ZJMX),SNDEC(ZIMXP,ZJMX),            POTM.22
     2  CSDEC(ZIMXP,ZJMX),SN2DEC(ZIMXP,ZJMX),SNCSDC(ZIMXP,ZJMX),            POTM.23
     3  RLATMP(ZJMX),RLONMP(ZJMX),DLONS(ZJMX)                               POTM.24
C     ****                                                                  POTM.25
C     ****     LOOP OVER MAGNETIC LATITUDE                                  POTM.26
C     ****                                                                  POTM.27
      DO J=1,JMAXM                                                          POTM.28
C      ****                                                                 POTM.29
C      ****     FILL ARRAYS IFLAG, DLAT, DLON, RATIO                        POTM.30
C      ****                                                                 POTM.31
       DO I=1,IMAXM                                                         POTM.32
        DLAT(I)=YLATM(J)                                                    POTM.33
        DLON(I)=YLONM(I)-DLONS(1)                                           POTM.34
       ENDDO                                                                POTM.35
C      ****                                                                 POTM.36
C      ****     CALCULATE HEELIS POTENTIAL IN PHIHM(IMAXMP,JMAXM)           POTM.37
C      ****                                                                 POTM.38
       IF(DABS(YLATM(J)).GT.C(110)/6.)THEN                                  POTM.39
        CALL FLOWV32(PHIHM(1,J))                                            POTM.40
       ELSE                                                                 POTM.41
        DO I=1,IMAXM                                                        POTM.42
         PHIHM(I,J)=0.                                                      POTM.43
        ENDDO                                                               POTM.44
       ENDIF                                                                POTM.45
      ENDDO                                                                 POTM.46
C     ****                                                                  POTM.47
C     ****    PERIODIC POINTS                                               POTM.48
C     ****                                                                  POTM.49
      DO J=1,JMAXM                                                          POTM.50
       PHIHM(IMAXMP,J)=PHIHM(1,J)                                           POTM.51
      ENDDO                                                                 POTM.52
      RETURN                                                                POTM.53
      END                                                                   POTM.54
C                                                                           POTM.55
*DK RTHS                                                                    RTHS.1
      SUBROUTINE RTHS(RHS,IMX0,JMX0,JMAXMH,C0)                              RTHS.2
      IMPLICIT REAL*8 (A-H,O-Z)                                             RTHS.3
C     ****                                                                  RTHS.4
C     ****     INSERT RHS IN STENCIL C0(10)                                 RTHS.5
C     ****                                                                  RTHS.6
C     ****     INPUT:                                                       RTHS.7
C     ****       RHS((IMX0+1)/2,(JMX0+1)/2)   RIGHT HAND SIDE OF PDE        RTHS.8
C     ****                                                                  RTHS.9
C     ****       CS(JMX0) COSINE ARRAY OF JMX0 LATITUDE VALUES              RTHS.10
C     ****                                                                  RTHS.11
C     ****       IMX0, JMX0   DIMENSION OF FINEST GRID IS IMX0 BY JMX0      RTHS.12
C     ****                                                                  RTHS.13
C     ****     OUTPUT:                                                      RTHS.14
C     ****       C0(IMX0,JMX0,10)                                           RTHS.15
C     ****         IN C0(I,J,10) FOR I=1,IMX0 AND J=1,JMX0                  RTHS.16
C     ****                                                                  RTHS.17
      DIMENSION RHS(IMX0,1),C0(IMX0,JMX0,10)                                RTHS.18
C     ****                                                                  RTHS.19
C     ****     CALCULATE NUMERICAL QUADRATURES OF WEIGHTED RHS AND          RTHS.20
C     ****       WEIGHTS.                                                   RTHS.21
C     ****                                                                  RTHS.22
      DO J=1,JMX0                                                           RTHS.23
       JJ=JMAXMH-JMX0+J                                                     RTHS.24
       DO I=1,IMX0                                                          RTHS.25
        C0(I,J,10)=RHS(I,JJ)                                                RTHS.26
       ENDDO                                                                RTHS.27
      ENDDO                                                                 RTHS.28
      RETURN                                                                RTHS.29
      END                                                                   RTHS.30
C                                                                           RTHS.31
*DK STENCIL                                                              STENCIL.1
      SUBROUTINE STENCIL(ZIGM,CS,IMX0,JMX0,SYM,CEE,ARRAY,NCOEF,WK)       STENCIL.2
      IMPLICIT REAL*8 (A-H,O-Z)                                          STENCIL.3
C     ****                                                               STENCIL.4
C     ****     CALCULATES CONTRIBUTION TO 3 BY 3 STENCIL FROM            STENCIL.5
C     ****       COEFFICIENT ,ZIGM, AT EACH GRID POINT AND LEVEL         STENCIL.6
C     ****                                                               STENCIL.7
C     ****     INPUT:                                                    STENCIL.8
C     ****       ZIGM((IMX0+1)/2,(JMX0+1)/2)   ARRAY OF COEFFICIENT      STENCIL.9
C     ****         VALUES.                                               STENCIL.10
C     ****       CS(-15:JMX0) ARRAY OF COSINES OF LATITUDE VALUES        STENCIL.11
C     ****       IMX0, JMX0  LONGITUDE AND LATITUDE DIMENSIONS           STENCIL.12
C     ****         OF FINEST GRID LEVEL                                  STENCIL.13
C     ****       SYM = 1.  IF ZIGM SYMMETRIC W.R.T. EQUATOR              STENCIL.14
C     ****           = -1.  FOR ANTISYMMETRY                             STENCIL.15
C     ****       NCOEF = NUMBER IDENTIFYING COEFFICIENT                  STENCIL.16
C     ****         NCOEF = 1 FOR ZIGM11                                  STENCIL.17
C     ****         NCOEF = 2 FOR ZIGM12 (=ZIGMC+ZIGM2)                   STENCIL.18
C     ****         NCOEF = 3 FOR ZIGM21 (=ZIGMC-ZIGM2)                   STENCIL.19
C     ****         NCOEF = 4 FOR ZIGM22                                  STENCIL.20
C     ****                                                               STENCIL.21
C     ****     OUTPUT:                                                   STENCIL.22
C     ****       CEE    STENCIL ARRAY CONSISTING OF                      STENCIL.23
C     ****         C0(IMX0,JMX0,10), C1(IMX1,JMX1,9), C2(IMX2,JMX2,9),   STENCIL.24
C     ****         C3(IMX3,JMX3,9), C4(IMX4,JMX4,9)                      STENCIL.25
C     ****           WHERE:   IMX1 = (IMX0+1)/2, JMX1 = (JMX0+1)/2,      STENCIL.26
C     ****                    IMX2 = (IMX1+1)/2, JMX2 = (JMX1+1)/2,      STENCIL.27
C     ****              AND SO ON                                        STENCIL.28
C     ****                                                               STENCIL.29
C     ****     WORK SPACE:                                               STENCIL.30
C     ****       ARRAY(-15:IMX0+16,-15:JMX0+2)                           STENCIL.31
C     ****         USED FOR INTERPOLATED AND EXTENDED ARRAY OF           STENCIL.32
C     ****         COEFFICIENT ZIGM                                      STENCIL.33
C     ****                                                               STENCIL.34
      DIMENSION ZIGM(IMX0,JMX0),CS(1),CEE(1),ARRAY(IMX0+32,1)            STENCIL.35
C     ****                                                               STENCIL.36
C     ****     PERFORM HALF WAY INTERPOLATION AND EXTEND ZIGM IN ARRAY   STENCIL.37
C     ****                                                               STENCIL.38
      CALL HTRPEX(ZIGM,ARRAY,IMX0,JMX0,SYM)                              STENCIL.39
C     ****                                                               STENCIL.40
C     ****     CALCULATE CONTRIBUTION TO STENCIL FOR EACH GRID POINT     STENCIL.41
C     ****       AND LEVEL                                               STENCIL.42
C     ****                                                               STENCIL.43
      NC=1                                                               STENCIL.44
      IMX=IMX0                                                           STENCIL.45
      JMX=JMX0                                                           STENCIL.46
      DO N=1,5                                                           STENCIL.47
       CALL CNM(ARRAY,CS,IMX0,JMX0,IMX,JMX,CEE(NC),NCOEF,WK)             STENCIL.48
       NC=NC+9*IMX*JMX                                                   STENCIL.49
       IF(N.EQ.1)NC=NC+IMX*JMX                                           STENCIL.50
       IMX=(IMX+1)/2                                                     STENCIL.51
       JMX=(JMX+1)/2                                                     STENCIL.52
      ENDDO                                                              STENCIL.53
      RETURN                                                             STENCIL.54
      END                                                                STENCIL.55
C                                                                        STENCIL.56
*DK STENMOD                                                              STENMOD.1
      SUBROUTINE STENMOD(IMX,JMX,C,PHIHM,P,COLATC)                       STENMOD.2
C     ****                                                               STENMOD.3
C     ****     MODIFY STENCIL TO SET POTENTIAL TO HEELIS VALUE           STENMOD.4
C     ****       WITHIN AURORAL CIRCLE                                   STENMOD.5
C     ****                                                               STENMOD.6
      include "param.h"                                                  STENMOD.7
      INCLUDE 'epara.h'                                                  STENMOD.8
      INCLUDE "econs.h"                                                  STENMOD.9
      DIMENSION C(IMX,JMX,10),PHIHM(IMX0,JMX0),COLATC(IMX0,JMX0),         STENMOD.10
     1          P(IMX0,JMX0)                                             STENMOD.11
C     ****                                                               STENMOD.12
C     ****     COMPUTE SEPARATION OF GRID POINTS FOR THIS RESOLUTION     STENMOD.13
C     ****                                                               STENMOD.14
      NINT=(IMX0-1)/(IMX-1)                                              STENMOD.15
      I0=1-NINT                                                          STENMOD.16
      J0=1-NINT                                                          STENMOD.17
C     ****                                                               STENMOD.18
C     ****     IF (NINT.EQ.1), WE ARE AT HIGHEST RESOLUTION.             STENMOD.18
C     ****       CORRECT RHS, WHICH IS IN C(10)                          STENMOD.20
C     ****                                                               STENMOD.21
      IF(NINT.EQ.1)THEN                                                  STENMOD.22
       DO J=1,JMX                                                        STENMOD.23
        DO I=1,IMX                                                       STENMOD.24
         C(I,J,10)=P(I,J)*C(I,J,10)+(1.-P(I,J))*C(I,J,9)*                STENMOD.25
     1              (DLATM/(10.*DTR))**2*PHIHM(I,J)                      STENMOD.26
        ENDDO                                                            STENMOD.27
       ENDDO                                                             STENMOD.28
      ENDIF                                                              STENMOD.29
C     ****                                                               STENMOD.30
C     ****     MODIFY STENCIL, C(I,J,N),N=1,9)                           STENMOD.31
C     ****                                                               STENMOD.32
      DO J=1,JMX                                                         STENMOD.33
       JJ=J0+J*NINT                                                      STENMOD.34
       DO N=1,8                                                          STENMOD.35
        DO I=1,IMX                                                       STENMOD.36
         C(I,J,N)=C(I,J,N)*P(I0+I*NINT,JJ)                               STENMOD.37
        ENDDO                                                            STENMOD.38
       ENDDO                                                             STENMOD.39
       DO I=1,IMX                                                        STENMOD.40
        C(I,J,9)=C(I,J,9)*P(I0+I*NINT,JJ)+(1-P(I0+I*NINT,JJ))*           STENMOD.41
     1           C(I,J,9)*(DLATM*DBLE(NINT)/(10.*DTR))**2                STENMOD.42
       ENDDO                                                             STENMOD.43
      ENDDO                                                              STENMOD.44
      RETURN                                                             STENMOD.45
      END                                                                STENMOD.46
C                                                                        STENMOD.47
*DK THREED                                                                THREED.1
      SUBROUTINE THREED                                                   THREED.2
C     ****                                                                THREED.3
C     ****     CALCULATE 3-D DYNAMO POTENTIAL IN GEOMAGNETIC COORDINATES  THREED.4
C     ****                                                                THREED.5
C     ****     PHIM(IMAXMP,JMAXM) = 2-D POTENTIAL IN GEOMAGNETIC          THREED.6
C     ****      COORDINATES OUTPUT BY DYNAMO CODE.                        THREED.7
C     ****                                                                THREED.8
C     ****     PHIM3D(IMAXMP,JMAXM,KMAXP) = 3-D DYNAMO POTENTIAL IN       THREED.9
C     ****       GEOMAGNETIC COORDINATES GENERATED BY THIS SUBROUTINE.    THREED.10
C     ****                                                                THREED.11
      include "param.h"                                                   THREED.12
      include 'const.h'
      INCLUDE 'epara.h'                                                   THREED.13
      include "econs.h"                                                   THREED.14
      include "edphi.h"                                                   THREED.15
      include "eterp.h"                                                   THREED.16
      DIMENSION THETAM(IMAXMP),ISLOT(IMAXMP),JSLOT(IMAXMP),               THREED.17
     1  PSLOT(IMAXMP),QSLOT(IMAXMP),PHIM3D(IMAXMP,JMAXM),                  THREED.18
     2  ZZM(IMAXMP,JMAXM)
C     ****                                                                THREED.19
C     ****     FOR EACH GEOMAGNETIC GRID POINT OVER THE 3-D GRID,         THREED.20
C     ****       DETERMINE THE LOCATION OF THE FOOT OF THE MAGNETIC       THREED.21
C     ****       FIELD PASSING THROUGH THE GRID POINT.  ASSUME THAT       THREED.22
C     ****       THE FIELD LINES ARE DIPOLE IN SHAPE.  RECALLING THAT     THREED.23
C     ****       THE MODEL ASSUMES CONSTANT ELECTRIC POTENTIAL ALONG      THREED.24
C     ****       EACH FIELD LINE, INTERPOLATION IN THE 2-D ARRAY PHIM     THREED.25
C     ****       GIVES THE POTENTIAL AT THE ORIGINAL GRID POINT.          THREED.26
C     ****                                                                THREED.27
C     ****     LOOP OVER GEOMAGNETIC GRID.                                THREED.28
C     ****                                                                THREED.29     
      DO K=1,KMAXP1
       READ(14)((ZZM(I,J),I=1,IMAXMP),J=1,JMAXM)
       DO J=2,JMAXM-1                                                     THREED.31
C       ****                                                              THREED.32
C       ****     COSLTM = COS(THETAS)                                     THREED.33
C       ****                                                              THREED.34
        COSLTM=DCOS(YLATM(J))                                             THREED.35
        DO I=1,IMAXM                                                      THREED.36
C        ****                                                             THREED.37
C        ****     THETAM(I) = THETAS FOR EACH VALUE OF I                  THREED.38
C        ****                                                             THREED.39
         THETAM(I)=DACOS(DSQRT(R0/(RE+ZZM(I,J)))*COSLTM)                  THREED.40
C        ****                                                             THREED.41
C        ****     FIND POSITION OF THETAM IN TABLE OF                     THREED.42
C        ****       THETA0 VS. THETAS                                     THREED.43
C        ****                                                             THREED.44
         PSLOT(I)=THETAM(I)*180./PI+1.                                    THREED.45
         ISLOT(I)=PSLOT(I)                                                THREED.46
         PSLOT(I)=PSLOT(I)-DBLE(ISLOT(I))                                 THREED.47
C        ****                                                             THREED.48
C        ****     DETERMINE THATA0 BY INTERPOLATION IN TABLE              THREED.49
C        ****                                                             THREED.50
         THETAM(I)=(1.-PSLOT(I))*TABLE(ISLOT(I),2)+PSLOT(I)*              THREED.51
     1             TABLE(ISLOT(I)+1,2)                                    THREED.52
C        ****                                                             THREED.53
C        ****     NOW IDENTIFY GEOMAGNETIC GRID ELEMENT IN WHICH          THREED.54
C        ****       THE POINT LIES IN PHIM                                THREED.55
C        ****                                                             THREED.56
         ISLOT(I)=I                                                       THREED.57
         PSLOT(I)=0.                                                      THREED.58
         QSLOT(I)=(THETAM(I)+PI/2)/DLATM+1.                               THREED.59
         JSLOT(I)=QSLOT(I)                                                THREED.60
         QSLOT(I)=QSLOT(I)-DBLE(JSLOT(I))                                 THREED.61
C        ****                                                             THREED.62
C        ****     INTERPOLATE IN PHIM USING BILINEAR INTERPOLATION        THREED.63
C        ****                                                             THREED.64
         PHIM3D(I,J)=(1-QSLOT(I))*((1-PSLOT(I))*                          THREED.65
     1                 PHIM(ISLOT(I),JSLOT(I))+PSLOT(I)*                  THREED.66
     2                 PHIM(ISLOT(I)+1,JSLOT(I)))+QSLOT(I)*               THREED.67
     3                 ((1.-PSLOT(I))*PHIM(ISLOT(I),JSLOT(I)+1)+          THREED.68
     4                 PSLOT(I)*PHIM(ISLOT(I)+1,JSLOT(I)+1))              THREED.69
        ENDDO                                                             THREED.70
       ENDDO                                                              THREED.71
C      ****                                                               THREED.73
C      ****     CALCULATE VALUE OF PHIM3D AT GEOMAGNETIC POLES            THREED.74
C      ****                                                               THREED.75
       PHIMS=0.                                                           THREED.76
       PHIMN=0.                                                           THREED.76
       DO I=1,IMAXM
        PHIMS=PHIMS+PHIM(I,1)                                             THREED.77
        PHIMn=PHIMn+PHIM(I,JMAXM)                                         THREED.77
       ENDDO
       DO I=1,IMAXM                                                       THREED.81
        PHIM3D(I,1)=PHIMS/DBLE(IMAXM)                                     THREED.82
        PHIM3D(I,JMAXM)=PHIMN/DBLE(IMAXM)                                 THREED.83
       ENDDO                                                              THREED.84
C      ****                                                               THREED.86
C      ****     PERIODIC POINT                                            THREED.87
C      ****                                                               THREED.88
       DO J=1,JMAXM                                                       THREED.90
        PHIM3D(IMAXMP,J)=PHIM3D(1,J)                                      THREED.91
       ENDDO                                                              THREED.92
C      ****                                                               THREED.86
C      ****     WRITE POTENTIAL TO A TEMPORAL FILE                        THREED.87
C      ****                                                               THREED.88
       WRITE(11)((PHIM3D(I,J),I=1,IMAXMP),J=1,JMAXM)
      ENDDO                                                               THREED.93
      REWIND(11)
      REWIND(14)
      RETURN                                                              THREED.122
      END                                                                 THREED.123
C                                                                         THREED.124
*DK                                                                     TRANSFR2.1
      SUBROUTINE TRANSF2                                                TRANSFR2.3
C     ****                                                              TRANSFR2.4
C     ****     CALCULATE QUANTITIES NEEDED FOR FIELD-LINE INTEGRATIONS  TRANSFR2.5
C     ****       AND TRANSFORM TO GEOMAGNETIC COORDINATE SYSTEM.        TRANSFR2.6
C     ****                                                              TRANSFR2.7
      include "param.h"                                                 TRANSFR2.8
      INCLUDE 'epara.h'                                                 TRANSFR2.9
      include "ecoef.h"                                                 TRANSFR2.10
      include 'const.h'
      include "econs.h"                                                 TRANSFR2.11
      include "eterp.h"                                                 TRANSFR2.12
      include "edphi.h"                                                 TRANSFR2.14
      include "efild.h"                                                 TRANSFR2.15
      include "start.h"                                                 TRANSFR2.16
C     ****                                                              TRANSFR2.17
C     ****     FIELDS TO BE TRANSFORMED TO GEOMAGNETIC SPACE            TRANSFR2.18
C     ****       EXTRAPOLATED DOWN TO 90 KM                             TRANSFR2.19
C     ****                                                              TRANSFR2.20
       COMMON/TRNSFR/ZZ(IMAXGP,0:JMAXGP,2),                             TRANSFR2.21
     1               SSIGMA1(IMAXGP,0:JMAXGP),                          TRANSFR2.22
     2               SSIGMA2(IMAXGP,0:JMAXGP),                          TRANSFR2.23
     3               ADOTV(IMAXGP,0:JMAXGP,2),                          TRANSFR2.24
     4               ADOTA(IMAXGP,0:JMAXGP,2),                          TRANSFR2.25
     5               A1DTA2(IMAXGP,0:JMAXGP),                           TRANSFR2.26
     6               SINI(IMAXGP,0:JMAXGP),                             TRANSFR2.27
     7               BE3(IMAXGP,0:JMAXGP),                              TRANSFR2.28
C     ****                                                              TRANSFR2.29
C     ****     THESE SAME QUANTITIES IN GEOMAGNETIC SPACE               TRANSFR2.30
C     ****       (ONE LATITUDE AT A TIME)                               TRANSFR2.31
C     ****                                                              TRANSFR2.32
     1              ZM(IMAXMP,2),                                       TRANSFR2.33
     2              SIGMA1M(IMAXMP),                                    TRANSFR2.34
     3              SIGMA2M(IMAXMP),                                    TRANSFR2.35
     4              ADOTVM(IMAXMP,2),                                   TRANSFR2.36
     5              ADOTAM(IMAXMP,2),                                   TRANSFR2.37
     6              A1A2M(IMAXMP),                                      TRANSFR2.38
     7              SINIAM(IMAXMP),                                     TRANSFR2.39
     8              BMODM(IMAXMP),                                      TRANSFR2.40
C     ****                                                              TRANSFR2.41
C     ****     99/2/23 Note: be3,BMODM are at 90 km, whereas BMOD       TRANSFR2.42
C     ****     (in fieldz.h) is at 130 km.                              TRANSFR2.43
C     ****                                                              TRANSFR2.44
C     ****     Equatorial quantities                                    TRANSFR2.45
C     ****                                                              TRANSFR2.46
     1              SIGMA1E(IMAXMP),                                    TRANSFR2.47
     2              SIGMA2E(IMAXMP),                                    TRANSFR2.48
     3              ADOTVE(IMAXMP,2),                                   TRANSFR2.49
C     ****                                                              TRANSFR2.50
C     ****     WORK SPACE                                               TRANSFR2.51
C     ****                                                              TRANSFR2.52
     1              TINT1(JMAXM),TINT2(JMAXM),TINT3(JMAXM)              TRANSFR2.53
C     ****                                                              TRANSFR2.54
C     ****     SET CONSTANTS                                            TRANSFR2.55
C     ****       RL1,RL2 ARE RATES AT WHICH SIGMA1 AND SIGMA2 DECAY     TRANSFR2.56
C     ****         WITH HEIGHT BELOW BOTTOM OF MODEL                    TRANSFR2.57
C     ****       Z0 IS LOWEST Z LEVEL FOR START OF FIELD LINE           TRANSFR2.58
C     ****         INTEGRATIONS                                         TRANSFR2.59
C     ****                                                              TRANSFR2.60
      DATA RL1,RL2,Z0/5.D5,3.D5,90.D5/                                  TRANSFR2.61
      DIMENSION ZZM(IMAXMP,JMAXM),CLM2(IMAXGP,JMAXG),
     1          U(IMAXGP,0:JMAXGP),V(IMAXGP,0:JMAXGP),
     2          W(IMAXGP,0:JMAXGP)        
C     ****                                                              TRANSFR2.138
C     ****     LINEAR INTERPOLATION OF GLOBAL CONDUCTIVITIES AND        TRANSFR2.139
C     ****     PRESSURE LEVEL HEIGHT TO ELECTRODYNAMICS GRID            TRANSFR2.140
C     ****     SLOT NESTED GRID VALUES INTO THE RELATED FIELDS          TRANSFR2.140
C     ****                                                              TRANSFR2.141
      CALL  LAMDAIN
      DO J=1,JMAXG                                                      TRANSFR2.63
       DO I=1,IMAXGP                                                    TRANSFR2.64
        SINALAT=DSIN(ALATM(I,J))                                        TRANSFR2.65
        CLM2(I,J)=1.-SINALAT*SINALAT                                    TRANSFR2.66
        SINI(I,J)=ZZB(I,J)/BMOD(I,J)                                    TRANSFR2.67
C       ****                                                            TRANSFR2.68
C       ****     be3 is in T                                            TRANSFR2.69
C       ****                                                            TRANSFR2.70
        BE3(I,J)=1.D-9*BE3ARR(I,J)                                      TRANSFR2.71
       ENDDO                                                            TRANSFR2.72
      ENDDO                                                             TRANSFR2.73
C     ****                                                              TRANSFR2.138
C     ****     ADOTA(n) = d(n)**2/D                                     TRANSFR2.139
C     ****     A1DTA2 = (d(1) dot d(2)) /D                              TRANSFR2.140
C     ****                                                              TRANSFR2.141
      DO J=1,JMAXG                                                      TRANSFR2.142
       DO I=1,IMAXGP                                                    TRANSFR2.143
        DDD=DDDARR(I,J)                                                 TRANSFR2.144
        ADOTA(I,J,1)=(DVEC(I,J,1,1)**2+DVEC(I,J,2,1)**2                 TRANSFR2.145
     1               +DVEC(I,J,3,1)**2)/DDD                             TRANSFR2.146
        ADOTA(I,J,2)=(DVEC(I,J,1,2)**2+DVEC(I,J,2,2)**2                 TRANSFR2.147
     1               +DVEC(I,J,3,2)**2)/DDD                             TRANSFR2.148
        A1DTA2(I,J)=(DVEC(I,J,1,1)*DVEC(I,J,1,2)+                       TRANSFR2.149
     1              DVEC(I,J,2,1)*DVEC(I,J,2,2)+                        TRANSFR2.150
     2              DVEC(I,J,3,1)*DVEC(I,J,3,2))/DDD                    TRANSFR2.151
       ENDDO                                                            TRANSFR2.152
      ENDDO                                                             TRANSFR2.153
C     ****                                                              TRANSFR2.175
C     ****     READ IN LOWEST HEIGHT VALUES                             TRANSFR2.176
C     ****                                                              TRANSFR2.180
      READ(12)((ZZ(I,J,2),I=1,IMAXGP),J=0,JMAXGP)
      DO J=2,JMAXM-1                                                    TRANSFR2.190
       CALL GRDINT(ZM(1,2),ZZ(1,0,2),IG,JG,WT,IMAXGP,                   TRANSFR2.192
     1             IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J)                  TRANSFR2.193
C      ****                                                             TRANSFR2.195
C      ****     STORE ZM IN 3-D ARRAY ZZM(IMAXMP,JMAXM,-2:KMAXP)        TRANSFR2.196
C      ****       FOR USE IN THREED                                     TRANSFR2.197
C      ****                                                             TRANSFR2.198
       DO I=1,IMAXM                                                      TRANSFR2.200
        ZZM(I,J)=ZM(I,2)                                                 TRANSFR2.201
       ENDDO                                                             TRANSFR2.202
      ENDDO
C     ****                                                               TRANSFR2.261
C     ****     COMPUTE POLAR VALUES OF ZM                                TRANSFR2.262
C     ****                                                               TRANSFR2.263
      RATE=3.*DBLE(IMAXM)
      VSDOT1=0.
      VSDOT2=0.
      DO I=1,IMAXM
       VSDOT1=VSDOT1+ZZM(I,2)
       VSDOT2=VSDOT2+ZZM(I,3)
      ENDDO        
      ZZM(1,1)=(4.*VSDOT1-VSDOT2)/RATE                                   TRANSFR2.267
      VSDOT1=0.
      VSDOT2=0.
      DO I=1,IMAXM
       VSDOT1=VSDOT1+ZZM(I,JMAXM-1)
       VSDOT2=VSDOT2+ZZM(I,JMAXM-2)
      ENDDO  
      ZZM(1,JMAXM)=(4.*VSDOT1-VSDOT2)/RATE                              TRANSFR2.270
C     ****                                                              TRANSFR2.272
C     ****     EXTEND OVER LONGITUDE                                    TRANSFR2.273
C     ****                                                              TRANSFR2.274
      DO I=2,IMAXM                                                      TRANSFR2.276
       ZZM(I,1)=ZZM(1,1)                                                TRANSFR2.277
       ZZM(I,JMAXM)=ZZM(1,JMAXM)                                        TRANSFR2.278
      ENDDO                                                             TRANSFR2.279
C     ****                                                              TRANSFR2.281
C     ****     PERIODIC POINTS                                          TRANSFR2.282
C     ****                                                              TRANSFR2.283
      DO J=1,JMAXM                                                      TRANSFR2.285
       ZZM(IMAXMP,J)=ZZM(1,J)                                           TRANSFR2.286
      ENDDO                                                             TRANSFR2.287
      WRITE(14)((ZZM(I,J),I=1,IMAXMP),J=1,JMAXM)
C     ****                                                               TRANSFR2.175
C     ****     TRANSFORM NEEDED FIELDS TO GEOMAGNETIC COORDINATE SYSTEM  TRANSFR2.176
C     ****       ONE MAGNETIC  LATITUDE AT A TIME                        TRANSFR2.177
C     ****                                                               TRANSFR2.180
      DO K=-2,KMAXP1
       DO J=0,JMAXGP
        DO I=1,IMAXGP
         ZZ(I,J,1)=ZZ(I,J,2)
        ENDDO
       ENDDO
       READ(12)((ZZ(I,J,2),I=1,IMAXGP),J=0,JMAXGP)       
       DO J=0,JMAXGP
        DO I=1,IMACGP
         READ(13)SSIGMA1(I,J),SSIGMA2(I,J),U(I,J),V(I,J),W(I,J)
        ENDDO
       ENDDO
       write(6,*)k
C      ****                                                              TRANSFR2.112
C      ****     ADOTV = ue1,ue2 (m/s)                                    TRANSFR2.113
C      ****                                                              TRANSFR2.114
       if(K.NE.KMAXP1)THEN
        IF(KK.LT.1)KK=1                                                  TRANSFR2.117
        DO J=1,JMAXG                                                     TRANSFR2.118
         DO I=1,IMAXGP                                                   TRANSFR2.119
          ROR=R0/(R0+.5*(ZZ(I,J,1)+ZZ(I,J,2))-Z0)                        TRANSFR2.120
          RAT=1.D-2*R0R**1.5                                             TRANSFR2.121
          ADOTV(I,J,1)=RAT*(DVEC(I,J,1,1)*U(I,J)+                        TRANSFR2.122
     2                   DVEC(I,J,2,1)*V(I,J)+                           TRANSFR2.123
     3                   DVEC(I,J,3,1)*W(I,J))                           TRANSFR2.124
          RAT=RAT*DSQRT((4.-3.*CLM2(I,J))/(4.-3.*R0R*CLM2(I,J)))         TRANSFR2.125
C         ****                                                           TRANSFR2.126
C         ****     Note: clm2 is being used here to represent the        TRANSFR2.127
C         ****     squared cosine of thequasi-dipole latitude, not of    TRANSFR2.128
C         ****     the M(90) latitude, since the windvalues are aligned  TRANSFR2.129
C         ****     vertically, not along the field line.                 TRANSFR2.130
C         ****                                                           TRANSFR2.131
          ADOTV(I,J,2)=RAT*(DVEC(I,J,1,2)*U(I,J)+                        TRANSFR2.132
     2                   DVEC(I,J,2,2)*V(I,J)+                           TRANSFR2.133
     3                   DVEC(I,J,3,2)*W(I,J))                           TRANSFR2.134
         ENDDO                                                           TRANSFR2.135
        ENDDO                                  
       ENDIF
C      ****                                                              TRANSFR2.178
C      ****     EQUATORIAL VALUES                                        TRANSFR2.179
C      ****                                                              TRANSFR2.180
       J=JMAXM/2+1                                                       TRANSFR2.181
        CALL GRDINT(SIGMA1E(1),SSIGMA1(1,0),IG,JG,WT,IMAXGP,            TRANSFR2.192
     1              IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J)                 TRANSFR2.193
        CALL GRDINT(SIGMA2E(1),SSIGMA2(1,0),IG,JG,WT,IMAXGP,            TRANSFR2.192
     1              IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J)                 TRANSFR2.193
        CALL GRDINT(ADOTVE(1,1),ADOTV(1,0,1),IG,JG,WT,IMAXGP,            TRANSFR2.192
     1              IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J)                 TRANSFR2.193
        CALL GRDINT(ADOTVE(1,2),ADOTV(1,0,2),IG,JG,WT,IMAXGP,            TRANSFR2.192
     1               IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J)                 TRANSFR2.193
       SIGMA1E(IMAXMP)=SIGMA1E(1)                                       TRANSFR2.185
       SIGMA2E(IMAXMP)=SIGMA2E(1)                                       TRANSFR2.185
       ADOTVE(IMAXMP,1)=ADOTVE(1,1)                                       TRANSFR2.185
       ADOTVE(IMAXMP,2)=ADOTVE(1,2)                                       TRANSFR2.185
C      ****                                                              TRANSFR2.187
C      ****     OTHER LATITUDES                                          TRANSFR2.188
C      ****                                                              TRANSFR2.189
       DO J=2,JMAXM-1                                                    TRANSFR2.190
        IF(K.EQ.-2)THEN
         CALL GRDINT(ZM(1,1),ZZ(1,0,1),IG,JG,WT,IMAXGP,                  TRANSFR2.192
     1               IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J)                 TRANSFR2.193        
         CALL GRDINT(ZM(1,2),ZZ(1,0,2),IG,JG,WT,IMAXGP,                  TRANSFR2.192
     1               IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J)                 TRANSFR2.193
        ELSE
         DO I=1,IMAXMP
          ZM(I,1)=ZM(I,2)
         ENDDO
         CALL GRDINT(ZM(1,2),ZZ(1,0,2),IG,JG,WT,IMAXGP,                  TRANSFR2.192
     1               IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J)                 TRANSFR2.193
         CALL GRDINT(SIGMA1M(1),SSIGMA1(1,0),IG,JG,WT,IMAXGP,            TRANSFR2.192
     1               IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J)                 TRANSFR2.193
         CALL GRDINT(SIGMA2M(1),SSIGMA2(1,0),IG,JG,WT,IMAXGP,            TRANSFR2.192
     1               IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J)                 TRANSFR2.193
         CALL GRDINT(ADOTVM(1,1),ADOTV(1,0,1),IG,JG,WT,IMAXGP,            TRANSFR2.192
     1               IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J)                 TRANSFR2.193
         CALL GRDINT(ADOTVM(1,2),ADOTV(1,0,2),IG,JG,WT,IMAXGP,            TRANSFR2.192
     1               IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J)                 TRANSFR2.193
         CALL GRDINT(ADOTAM(1,1),ADOTA(1,0,1),IG,JG,WT,IMAXGP,            TRANSFR2.192
     1               IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J)                 TRANSFR2.193
         CALL GRDINT(ADOTAM(1,2),ADOTA(1,0,2),IG,JG,WT,IMAXGP,            TRANSFR2.192
     1               IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J)                 TRANSFR2.193
         CALL GRDINT(A1A2M(1),A1DTA2(1,0),IG,JG,WT,IMAXGP,               TRANSFR2.192
     1               IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J)                 TRANSFR2.193
         CALL GRDINT(SINIAM(1),SINI(1,0),IG,JG,WT,IMAXGP,                TRANSFR2.192
     1               IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J)                 TRANSFR2.193
         CALL GRDINT(BMODM(1),BE3(1,0),IG,JG,WT,IMAXGP,                  TRANSFR2.192
     1               IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J)                 TRANSFR2.193
        ENDIF
C       ****                                                             TRANSFR2.195
C       ****     STORE ZM IN 3-D ARRAY ZZM(IMAXMP,JMAXM,-2:KMAXP)        TRANSFR2.196
C       ****       FOR USE IN THREED                                     TRANSFR2.197
C       ****                                                             TRANSFR2.198
        DO I=1,IMAXM                                                    TRANSFR2.200
         ZZM(I,J)=ZM(I,2)                                                 TRANSFR2.201
        ENDDO                                                           TRANSFR2.202
        IF(J.NE.JMAXM/2+1)THEN                                           TRANSFR2.204
C        ****                                                            TRANSFR2.205
C        ****     CALL INTGRL2 TO PERFORM FIELD LINE INTEGRATIONS        TRANSFR2.207
C        ****                                                            TRANSFR2.208
         CALL INTGRL2(J)                                                 TRANSFR2.209
        ENDIF                                                            TRANSFR2.210
       ENDDO                                                             TRANSFR2.211
C      ****                                                              TRANSFR2.261
C      ****     COMPUTE POLAR VALUES OF ZM                               TRANSFR2.262
C      ****                                                              TRANSFR2.263
       RATE=3.*DBLE(IMAXM)
       VSDOT1=0.
       VSDOT2=0.
       DO I=1,IMAXM
        VSDOT1=VSDOT1+ZZM(I,2)
        VSDOT2=VSDOT2+ZZM(I,3)
       ENDDO        
       ZZM(1,1)=(4.*VSDOT1-VSDOT2)/RATE                                 TRANSFR2.267
       VSDOT1=0.
       VSDOT2=0.
       DO I=1,IMAXM
        VSDOT1=VSDOT1+ZZM(I,JMAXM-1)
        VSDOT2=VSDOT2+ZZM(I,JMAXM-2)
       ENDDO  
       ZZM(1,JMAXM)=(4.*VSDOT1-VSDOT2)/RATE                             TRANSFR2.270
C      ****                                                             TRANSFR2.272
C      ****     EXTEND OVER LONGITUDE                                   TRANSFR2.273
C      ****                                                             TRANSFR2.274
       DO I=2,IMAXM                                                     TRANSFR2.276
        ZZM(I,1)=ZZM(1,1)                                               TRANSFR2.277
        ZZM(I,JMAXM)=ZZM(1,JMAXM)                                       TRANSFR2.278
       ENDDO                                                            TRANSFR2.279
C      ****                                                             TRANSFR2.281
C      ****     PERIODIC POINTS                                         TRANSFR2.282
C      ****                                                             TRANSFR2.283
       DO J=1,JMAXM                                                     TRANSFR2.285
        ZZM(IMAXMP,J)=ZZM(1,J)                                          TRANSFR2.286
       ENDDO                                                            TRANSFR2.287
       WRITE(14)((ZZM(I,J),I=1,IMAXMP),J=1,JMAXM)
      ENDDO                                                             TRANSFR2.211
C     ****                                                              TRANSFR2.212
C     ****     Equatorial values:                                       TRANDFR2.213
C     ****     Assume that quantities primarily dependent on Pedersen   TRANDFR2.214
C     ****       conductivity have field-line integrals 1/4 as large as TRANDFR2.215
C     ****       the averages for next-higherfield lines; quantities    TRANDFR2.216
C     ****       primarily dependent on Hall conductivity have field-   TRANDFR2.217
C     ****       line integrals 0.12 as large as the averages for next- TRANDFR2.218
C     ****       higher field lines.  Exact values chosen should not be TRANDFR2.219
C     ****       important for potential calculation, as long as they   TRANDFR2.220
C     ****       are physically reasonable and not too different from   TRANDFR2.221
C     ****       adjacent values.                                       TRANDFR2.222
C     ****                                                              TRANSFR2.223
      J=JMAXM/2+1                                                       TRANSFR2.224
      DO I=1,IMAXM                                                      TRANSFR2.225
       ZIGM11(I,J)=.125*(ZIGM11(I,J-1)+ZIGM11(I,J+1))                   TRANSFR2.226
       ZIGM22(I,J)=.125*(ZIGM22(I,J-1)+ZIGM22(I,J+1))                   TRANSFR2.227
       ZIGMC(I,J)=.125*(ZIGMC(I,J-1)+ZIGMC(I,J+1))                      TRANSFR2.228
       ZIGM2(I,J)=.06*(ZIGM2(I,J-1)+ZIGM2(I,J+1))                       TRANSFR2.229
       RIM(I,J,1)=.06*(RIM(I,J-1,1)+RIM(I,J+1,1))                       TRANSFR2.230
       RIM(I,J,2)=.06*(RIM(I,J-1,2)+RIM(I,J+1,2))                       TRANSFR2.231
      ENDDO                                                             TRANSFR2.232
C     ****                                                              TRANSFR2.233
C     ****     Using notation of Richmond (1995) on right-hand side     TRANSFR2.234
C     ****     below:                                                   TRANSFR2.235
C     ****     Sigma_(phi phi) = ZIGM11*abs(sin I_m)                    TRANSFR2.236
C     ****     Sigma_(lam lam) = ZIGM22/abs(sin I_m)                    TRANSFR2.237
C     ****     Sigma_(phi lam) = +-(ZIGM2-ZIGMC)                        TRANSFR2.238
C     ****     Sigma_(lam phi) = -+(ZIGM2+ZIGMC)                        TRANSFR2.239
C     ****     K_(m phi)^D     = RIM(1)*abs(sin I_m)                    TRANSFR2.240
C     ****     K_(m lam)^D     = +-RIM(2)                               TRANSFR2.241
C     ****                                                              TRANSFR2.242
C     ****     SCALE quantities to modified (0) magnetic latitude       TRANSFR2.243
C     ****     system, multiplying or dividing by abs(sin I_m)          TRANSFR2.244
C     ****     [inverse contained in DT1DTS] as necessary.              TRANSFR2.245
C     ****     Sign of K_(m lam)^D in southern hemisphere remains       TRANSFR2.246
C     ****     reversed.                                                TRANSFR2.247
C     ****                                                              TRANSFR2.248
      DO J=2,JMAXM-1                                                    TRANSFR2.249
       CORFAC=RCOS0S(J)/DT1DTS(J)                                       TRANSFR2.250
       DO I=1,IMAXM                                                     TRANSFR2.251
        ZIGM11(I,J)=ZIGM11(I,J)*CORFAC                                  TRANSFR2.252
        ZIGM22(I,J)=ZIGM22(I,J)/CORFAC                                  TRANSFR2.253
        RIM(I,J,1)=RIM(I,J,1)/DT1DTS(J)                                 TRANSFR2.254
        RIM(I,J,2)=RIM(I,J,2)/RCOS0S(J)                                 TRANSFR2.255
       ENDDO                                                            TRANSFR2.256
      ENDDO                                                             TRANSFR2.257
      DO J=1,6*JMAXM                                                    TRANSFR2.258
       ZIGM11(IMAXMP,J)=ZIGM11(1,J)                                     TRANSFR2.259
      ENDDO                                                             TRANSFR2.260
C     ****                                                              TRANSFR2.289
C     ****     POLEs                                                    TRANSFR2.290
C     ****                                                              TRANSFR2.291
      RATE=3.*DBLE(IMAXM)
      VSDOT1=0.
      VSDOT2=0.
      DO I=1,IMAXM
       VSDOT1=VSDOT1+ZIGM11(I,2)
       VSDOT2=VSDOT2+ZIGM11(I,3)
      ENDDO 
      ZIGM11(1,1)=(4.*VSDOT1-VSDOT)/RATE                               TRANSFR2.294
      VSDOT1=0.
      VSDOT2=0.
      DO I=1,IMAXM
       VSDOT1=VSDOT1+ZIGM11(I,JMAXM-1)
       VSDOT2=VSDOT2+ZIGM11(I,JMAXM-2)
      ENDDO 
      ZIGM11(1,JMAXM)=(4.*VSDOT1-VSDOT2)/RATE                          TRANSFR2.297
      VSDOT1=0.
      VSDOT2=0.
      DO I=1,IMAXM
       VSDOT1=VSDOT1+ZIGMC(I,2)
       VSDOT2=VSDOT2+ZIGMC(I,3)
      ENDDO 
      ZIGMC(1,1)=(4.*VSDOT1-VSDOT2)/RATE                               TRANSFR2.300
      VSDOT1=0.
      VSDOT2=0.
      DO I=1,IMAXM
       VSDOT1=VSDOT1+ZIGMC(I,JMAXM-1)
       VSDOT2=VSDOT2+ZIGMC(I,JMAXM-2)
      ENDDO 
      ZIGMC(1,JMAXM)=(4.*VSDOT1-VSDOT2)/RATE                            TRANSFR2.303
      VSDOT1=0.
      VSDOT2=0.
      DO I=1,IMAXM
       VSDOT1=VSDOT1+ZIGM2(I,2)
       VSDOT2=VSDOT2+ZIGM2(I,3)
      ENDDO 
      ZIGM2(1,1)=(4.*VSDOT1-VSDOT2)/RATE                                TRANSFR2.306
      VSDOT1=0.
      VSDOT2=0.
      DO I=1,IMAXM
       VSDOT1=VSDOT1+ZIGM2(I,JMAXM-1)
       VSDOT2=VSDOT2+ZIGM2(I,JMAXM-2)
      ENDDO
      ZIGM2(1,JMAXM)=(4.*VSDOT1-VSDOT2)/RATE                            TRANSFR2.309
      VSDOT1=0.
      VSDOT2=0.
      DO I=1,IMAXM
       VSDOT1=VSDOT1+ZIGM22(I,2)
       VSDOT2=VSDOT2+ZIGM22(I,3)
      ENDDO 
      ZIGM22(1,1)=(4.*VSDOT1-VSDOT2)/RATE                               TRANSFR2.312
      VSDOT1=0.
      VSDOT2=0.
      DO I=1,IMAXM
       VSDOT1=VSDOT1+ZIGM22(I,JMAXM-1)
       VSDOT2=VSDOT2+ZIGM22(I,JMAXM-2)
      ENDDO
      ZIGM22(1,JMAXM)=(4.*VSDOT1-VSDOT2)/RATE                           TRANSFR2.315
C     ****                                                              TRANSFR2.316
C     ****     EXTEND OVER LONGITUDE                                    TRANSFR2.317
C     ****                                                              TRANSFR2.318
      DO I=2,IMAXM                                                      TRANSFR2.319
      ZIGM11(I,    1)=ZIGM11(1,    1)                                   TRANSFR2.320
      ZIGM11(I,JMAXM)=ZIGM11(1,JMAXM)                                   TRANSFR2.321
      ZIGMC(I,    1)=ZIGMC(1,    1)                                     TRANSFR2.322
      ZIGMC(I,JMAXM)=ZIGMC(1,JMAXM)                                     TRANSFR2.323
      ZIGM2(I,    1)=ZIGM2(1,    1)                                     TRANSFR2.324
      ZIGM2(I,JMAXM)=ZIGM2(1,JMAXM)                                     TRANSFR2.325
      ZIGM22(I,    1)=ZIGM22(1,    1)                                   TRANSFR2.326
      ZIGM22(I,JMAXM)=ZIGM22(1,JMAXM)                                   TRANSFR2.327
      ENDDO                                                             TRANSFR2.328
C     ****                                                              TRANSFR2.329
C     ****     VECTOR (I1,I2)                                           TRANSFR2.330
C     ****                                                              TRANSFR2.331
      DO I=1,IMAXM                                                      TRANSFR2.332
       RIM(I,1,1)=.5*(RIM(I,2,1)-RIM(1+MOD(I-1+IMAXM/2,IMAXM),2,1))     TRANSFR2.333
       RIM(I,JMAXM,1)=.5*(RIM(I,JMAXM-1,1)-RIM(1+MOD(I-1+IMAXM/2,       TRANSFR2.334
     1                IMAXM),JMAXM-1,1))                                TRANSFR2.335
       RIM(I,1,2)=.5*(RIM(I,2,2)-RIM(1+MOD(I-1+IMAXM/2,IMAXM),2,2))     TRANSFR2.336
       RIM(I,JMAXM,2)=.5*(RIM(I,JMAXM-1,2)-RIM(1+MOD(I-1+IMAXM/2,       TRANSFR2.337
     1                IMAXM),JMAXM-1,2))                                TRANSFR2.338
      ENDDO                                                             TRANSFR2.339
C     ****                                                              TRANSFR2.340
C     ****     PERIODIC POINT                                           TRANSFR2.341
C     ****                                                              TRANSFR2.342
      DO J=1,6*JMAXM                                                    TRANSFR2.343
       ZIGM11(IMAXMP,J)=ZIGM11(1,J)                                     TRANSFR2.344
      ENDDO                                                             TRANSFR2.345
C     ****                                                              TRANSFR2.346
C     ****     FOLD SOUTHERN HEMISPHERE OVER ON TO NORTHERN             TRANSFR2.347
C     ****                                                              TRANSFR2.348
      DO J=1,JMAXMH                                                     TRANSFR2.349
       DO I=1,IMAXMP                                                    TRANSFR2.350
        ZIGM11(I,JMAXM+1-J)=(ZIGM11(I,JMAXM+1-J)+ZIGM11(I,J))           TRANSFR2.351
        ZIGMC(I,JMAXM+1-J)=ZIGMC(I,JMAXM+1-J)+ZIGMC(I,J)                TRANSFR2.352
C       ****                                                            TRANSFR2.353
C       ****     Reverse sign of ZIGMC to be compatible with Cicely's   TRANSFR2.354
C       ****                                                            TRANSFR2.355
        ZIGMC(I,JMAXM+1-J)=-ZIGMC(I,JMAXM+1-J)                          TRANSFR2.356
        ZIGM2(I,JMAXM+1-J)=ZIGM2(I,JMAXM+1-J)+ZIGM2(I,J)                TRANSFR2.357
        ZIGM22(I,JMAXM+1-J)=(ZIGM22(I,JMAXM+1-J)+ZIGM22(I,J))           TRANSFR2.358
        RIM(I,JMAXM+1-J,1)=(RIM(I,JMAXM+1-J,1)+RIM(I,J,1))              TRANSFR2.359
        RIM(I,JMAXM+1-J,2)=(RIM(I,JMAXM+1-J,2)+RIM(I,J,2))              TRANSFR2.360
       ENDDO                                                            TRANSFR2.361
      ENDDO                                                             TRANSFR2.362
C     ****                                                              TRANSFR2.363
C     ****     CALCULATE RHS OF PDE FROM RIM(1) AND RIM(2)              TRANSFR2.364
C     ****                                                              TRANSFR2.365
C     ****     FILL ARRAY TINT1 WITH COS(THETA0)                        TRANSFR2.366
C     ****                                                              TRANSFR2.367
      DO J=1,JMAXM                                                      TRANSFR2.368
       TINT1(J)=DCOS(-PI/2.+(J-1)*DLATM)                                TRANSFR2.369
      ENDDO                                                             TRANSFR2.370
C     ****                                                              TRANSFR2.371
C     ****     TRANSFER RIM(1) TO TINT3 IN /WORK/                       TRANSFR2.372
C     ****                                                              TRANSFR2.373
      JE=JMAXMH                                                         TRANSFR2.374
      DO J=2,JE-1                                                       TRANSFR2.375
       JJ=J+JE-1                                                        TRANSFR2.376
       DO I=1,IMAXM                                                     TRANSFR2.377
        TINT3(I)=RIM(I,JJ,1)                                            TRANSFR2.378
       ENDDO                                                            TRANSFR2.379
C      ****                                                             TRANSFR2.380
C      ****     PERIODIC POINTS FOR TINT2                               TRANSFR2.381
C      ****                                                             TRANSFR2.382
       DO I=1,2                                                         TRANSFR2.383
        TINT3(I-2)=TINT3(I-2+IMAXM)                                     TRANSFR2.384
        TINT3(I+IMAXM)=TINT3(I)                                         TRANSFR2.385
       ENDDO                                                            TRANSFR2.386
C      ****                                                             TRANSFR2.387
C      ****     PERFORM DIFFERENTIATION OF RIM(1) W.R.T. LAMDA          TRANSFR2.388
C      ****                                                             TRANSFR2.389
       DO I=1,IMAXM                                                     TRANSFR2.390
        RHS(I,J)=1./(DLONM*TINT1(JE+J-1))*.5*(TINT3(I+1)-               TRANSFR2.391
     1           TINT3(I-1))                                            TRANSFR2.392
       ENDDO                                                            TRANSFR2.393
      ENDDO                                                             TRANSFR2.394
C     ****                                                              TRANSFR2.395
C     ****     PERFORM DIFFERENTIATION OF RIM(2) W.R.T. THETA0          TRANSFR2.396
C     ****                                                              TRANSFR2.397
      DO J=JE+1,JMAXM-1                                                 TRANSFR2.398
       JJ=J-JE+1                                                        TRANSFR2.399
       DO I=1,IMAXM                                                     TRANSFR2.400
        RHS(I,JJ)=RHS(I,JJ)+1./(DLATM*TINT1(J))*                        TRANSFR2.401
     1            .5*(RIM(I,J+1,2)*TINT1(J+1)-                          TRANSFR2.402
     2            RIM(I,J-1,2)*TINT1(J-1))                              TRANSFR2.403
       ENDDO                                                            TRANSFR2.404
      ENDDO                                                             TRANSFR2.405
C     ****                                                              TRANSFR2.406
C     ****     NOW DEAL WITH  JMAXM                                     TRANSFR2.407
C     ****                                                              TRANSFR2.408
      VSDOT1=0.
      DO I=1,IMAXM
       VSDOT1=VSDOT1+RIM(I,JMAXM-1,2)
      ENDDO                                                             TRANSFR2.409
      RHS(1,JMAXMH)=-2/DBLE(IMAXM)*VSDOT1/TINT1(JMAXM-1)                TRANSFR2.411
C     ****                                                              TRANSFR2.412
C     ****     AND EQUATOR                                              TRANSFR2.413
C     ****                                                              TRANSFR2.414
      DO I=1,IMAXM                                                      TRANSFR2.415
       RHS(I,1)=1./DLATM*.5*TINT1(JMAXMH+1)*RIM(I,JMAXMH+1,2)           TRANSFR2.416
      ENDDO                                                             TRANSFR2.417
C     ****                                                              TRANSFR2.418
C     ****     EXTEND OVER LONGITUDE                                    TRANSFR2.419
C     ****                                                              TRANSFR2.420
      DO I=2,IMAXM                                                      TRANSFR2.421
       RHS(I,JMAXMH)=RHS(1,JMAXMH)                                      TRANSFR2.422
      ENDDO                                                             TRANSFR2.423
C     ****                                                              TRANSFR2.424
C     ****     PERIODIC POINTS                                          TRANSFR2.425
C     ****                                                              TRANSFR2.426
      DO J=1,JMAXMH                                                     TRANSFR2.427
       RHS(IMAXMP,J)=RHS(1,J)                                           TRANSFR2.428
      ENDDO                                                             TRANSFR2.429
C     ****                                                              TRANSFR2.430
C     ****     SCALE  (MULTIPLY BY Earth radius in meters = R0*1.E-2)   TRANSFR2.431
C     ****                                                              TRANSFR2.432
      DO J=1,JMAXMH                                                     TRANSFR2.433
       DO I=1,IMAXMP                                                    TRANSFR2.434
        RHS(I,J)=RHS(I,J)*R0*1.D-2                                      TRANSFR2.435
       ENDDO                                                            TRANSFR2.436
      ENDDO                                                             TRANSFR2.437
      REWIND(12)
      REWIND(13)
      REWIND(14)
      RETURN                                                            TRANSFR2.438
      END                                                               TRANSFR2.439
C                                                                       TRANSFR2.440
*DK VDRIFT                                                                VDRIFT.1
      SUBROUTINE VDRIFT                                                   VDRIFT.2
C     ****                                                                VDRIFT.3
C     ****     Calculate 3-d electric field by numerical differentiation  VDRIFT.4
C     ****     of the 3-d dynamo potential in geomagnetic coordinates.    VDRIFT.5
C     ****                                                                VDRIFT.6
C     ****     Transform the three components of the field to geographic  VDRIFT.7
C     ****     space.                                                     VDRIFT.8
C     ****                                                                VDRIFT.9
C     ****     Assemble needed common blocks:                             VDRIFT.10
C     ****                                                                VDRIFT.11
      INCLUDE "param.h"                                                   VDRIFT.12
      INCLUDE 'epara.h'                                                   VDRIFT.13
      INCLUDE "const.h"                                                   VDRIFT.14
      INCLUDE "econs.h"                                                   VDRIFT.15
C     ****                                                                VDRIFT.18
C     ****     Work space in /  /                                         VDRIFT.19
C     ****                                                                VDRIFT.20
      common/phidd/PHI3D(0:IMAXMP2,0:JMAXMP,3)
      DIMENSION EX(IMAXGP,JMAXG,3),                                       VDRIFT.21
     1   EY(IMAXGP,JMAXG,3),EZ(IMAXGP,JMAXG,3),HT(IMAXGP,JMAXG,4)         VDRIFT.22
C     ****                                                                VDRIFT.24
C     ****     READ IN FIRST ALTITUDE (K=1) POTENTAIL                     VDRIFT.25
C     ****                                                                VDRIFT.27
      READ(11)((PHI3D(I,J,3),I=1,IMAXMP),J=1,JMAXM)
      DO K=1,4
       READ(12)((HT(I,J,4),I=1,IMAXGP),J=1,JMAXG)      
      ENDDO
C     ****                                                                VDRIFT.35
C     ****     Add points corresponding to J = 0 and JMAX+1 performing    VDRIFT.36
C     ****     extrapolation across the poles.                            VDRIFT.37
C     ****                                                                VDRIFT.38
      IMAXMH=IMAXM/2                                                      VDRIFT.39
      DO I=1,IMAXMP                                                       VDRIFT.41
       PHI3D(I,0,3)=PHI3D(1+MOD(I-1+IMAXMH,IMAXM),2,3)                    VDRIFT.42
       PHI3D(I,JMAXM+1,3)=PHI3D(1+MOD(I-1+IMAXMH,IMAXM),JMAXM-1,3)        VDRIFT.43
      ENDDO                                                               VDRIFT.45
C     ****                                                                VDRIFT.46
C     ****     Add periodic points to PHI3D                               VDRIFT.47
C     ****                                                                VDRIFT.48
      DO J=0,JMAXM+1                                                      VDRIFT.50
       PHI3D(0,J,3)=PHI3D(IMAXM,J,3)                                      VDRIFT.51
       PHI3D(IMAXMP+1,J,3)=PHI3D(2,J,3)                                   VDRIFT.52
      ENDDO                                                               VDRIFT.54
C     ****                                                                VDRIFT.24
C     ****     K CYCLE                                                    VDRIFT.25
C     ****                                                                VDRIFT.27
      K1=1
      K2=2
      K3=3
      DO K=1,KMAXP1
       IF(K.EQ.1)THEN
C       ****                                                                VDRIFT.24
C       ****     COPY POTENTIALS TO LOWER HEIGHT BUFFER                     VDRIFT.25
C       ****
        DO J=0,JMAXMP
         DO I=0,IMAXMP+1
          PHI3D(I,J,2)=PHI3D(I,J,3)
         ENDDO
        ENDDO
        DO J=1,JMAXG
         DO I=1,IMAXGP
          HT(I,J,3)=HT(I,J,2)
         ENDDO
        ENDDO
       ELSE IF(K.NE.KMAXP1)THEN
C       ****                                                                VDRIFT.24
C       ****     COPY POTENTIALS TO LOWER HEIGHT BUFFER                     VDRIFT.25
C       ****
        DO J=0,JMAXMP
         DO I=0,IMAXMP+1
          PHI3D(I,J,1)=PHI3D(I,J,2)
          PHI3D(I,J,2)=PHI3D(I,J,3)
         ENDDO
        ENDDO
        DO J=1,JMAXG
         DO I=1,IMAXGP
          HT(I,J,1)=HT(I,J,2)
          HT(I,J,2)=HT(I,J,3)
          HT(I,J,3)=HT(I,J,4)
         ENDDO
        ENDDO
       ENDIF
C      ****                                                                VDRIFT.24
C      ****     READ IN NEXT ALTITUDE (K) POTENTAIL                        VDRIFT.25
C      ****
       IF(K.NE.KMAXP1)READ(11)((PHI3D(I,J,3),I=1,IMAXMP),J=1,JMAXM)
       IF(K.NE.KMAXP1)READ(12)((HT(I,J,4),I=1,IMAXGP),J=1,JMAXG)
C      ****                                                                VDRIFT.35
C      ****     Add points corresponding to J = 0 and JMAX+1 performing    VDRIFT.36
C      ****     extrapolation across the poles.                            VDRIFT.37
C      ****                                                                VDRIFT.38
       IMAXMH=IMAXM/2                                                      VDRIFT.39
       DO I=1,IMAXMP                                                       VDRIFT.41
        PHI3D(I,0,3)=PHI3D(1+MOD(I-1+IMAXMH,IMAXM),2,3)                    VDRIFT.42
        PHI3D(I,JMAXM+1,3)=PHI3D(1+MOD(I-1+IMAXMH,IMAXM),JMAXM-1,3)        VDRIFT.43
       ENDDO                                                               VDRIFT.45
C      ****                                                                VDRIFT.46
C      ****     Add periodic points to PHI3D                               VDRIFT.47
C      ****                                                                VDRIFT.48
       DO J=0,JMAXM+1                                                      VDRIFT.50
        PHI3D(0,J,3)=PHI3D(IMAXM,J,3)                                      VDRIFT.51
        PHI3D(IMAXMP+1,J,3)=PHI3D(2,J,3)                                   VDRIFT.52
       ENDDO                                                               VDRIFT.54
C      ****
C      ****
C      ****
       IF(K.EQ.1)THEN
        CALL EMFIELD(EX(1,1,K1),EY(1,1,K1),EZ(1,1,K1),1)
       ELSE IF(K.EQ.2)THEN
        CALL EMFIELD(EX(1,1,K2),EY(1,1,K2),EZ(1,1,K2),0)
       ELSE IF(K.EQ.3)THEN
        CALL EMFIELD(EX(1,1,K3),EY(1,1,K3),EZ(1,1,K3),0)
        CALL IDRIFT(EX(1,1,K1),EY(1,1,K1),EZ(1,1,K3),EZ(1,1,K2),1,HT,1)
        CALL IDRIFT(EX(1,1,K2),EY(1,1,K2),EZ(1,1,K2),EZ(1,1,K2),2,HT,2)
        CALL IDRIFT(EX(1,1,K3),EY(1,1,K3),EZ(1,1,K3),EZ(1,1,K2),0,HT,3)
       ELSE IF(K.GE.4.AND.K.LE.KMAX-1)THEN
        CALL EMFIELD(EX(1,1,K1),EY(1,1,K1),EZ(1,1,K1),0)
        CALL IDRIFT(EX(1,1,K1),EY(1,1,K1),EZ(1,1,K1),EZ(1,1,K2),0,HT,K)
       ELSE IF(K.EQ.KMAX)THEN
        CALL EMFIELD(EX(1,1,K1),EY(1,1,K1),EZ(1,1,K2),0)
        CALL IDRIFT(EX(1,1,K1),EY(1,1,K1),EZ(1,1,K2),EZ(1,1,K2),0,
     1              HT,KMAX)
       ELSE IF(K.EQ.KMAXP1)THEN
        CALL EMFIELD(EX(1,1,K1),EY(1,1,K1),EZ(1,1,K3),1)       
        CALL IDRIFT(EX(1,1,K1),EY(1,1,K1),EZ(1,1,K2),EZ(1,1,K1),-1,
     1              HT,KMAXP1)
       ENDIF
      ENDDO
      REWIND(11)
      REWIND(12)
      RETURN
      END
*DK EMFIELD
      SUBROUTINE EMFIELD(EEX,EEY,EEZ,IN)
      INCLUDE "param.h"                                                   VDRIFT.12
      INCLUDE 'epara.h'                                                   VDRIFT.13
      INCLUDE "const.h"                                                   VDRIFT.14
      INCLUDE "econs.h"                                                   VDRIFT.15
      INCLUDE "efild.h"                                                   VDRIFT.17
C     ****                                                                VDRIFT.18
C     ****     Work space in /  /                                         VDRIFT.19
C     ****                                                                VDRIFT.20
      common/phidd/PHI3D(0:IMAXMP2,0:JMAXMP,3)
      DIMENSION EEX(IMAXGP,JMAXG) ,                                       VDRIFT.21
     1   EEY(IMAXGP,JMAXG),EEZ(IMAXGP,JMAXG),THETAM(IMAXG),               VDRIFT.22
     2   QSLOT(IMAXG),PSLOT(IMAXG),ISLOT(IMAXG),JSLOT(IMAXG)              VDRIFT.23
      DIMENSION EMX(IMAXMP,JMAXM),EMY(IMAXMP,JMAXM),EMZ(IMAXMP,JMAXM)
C     ****                                                                VDRIFT.55
C     ****     Calculate EMY, the meridional component of the electric    VDRIFT.56
C     ****     field, for all geomagnetic grid points.                    VDRIFT.57
C     ****                                                                VDRIFT.58
C     ****     Note: The factor DT0DTS converts from the latitudinally    VDRIFT.59
C     ****     distorted geomagnetic grid to standard apex latitudes.     VDRIFT.60
C     ****                                                                VDRIFT.61
      DO J=1,JMAXM                                                        VDRIFT.63
       DO I=1,IMAXMP                                                      VDRIFT.64
        EMY(I,J)=-(PHI3D(I,J+1,2)-PHI3D(I,J-1,2))/(2*DLATM)*DT0DTS(J)      
       ENDDO                                                              VDRIFT.67
      ENDDO                                                               VDRIFT.69
C     ****                                                                VDRIFT.70
C     ****     EMX = zonal electric field for all geomagnetic grid        VDRIFT.71
C     ****     points                                                     VDRIFT.72
C     ****                                                                VDRIFT.73
C     ****     Note: The factor RCOS0S converts from the latitudinally    VDRIFT.74
C     ****     distortded geomagnetic grid to standard apex latitudes.    VDRIFT.75
C     ****                                                                VDRIFT.76
      DO J=2,JMAXM-1                                                      VDRIFT.78
       CSTH0=DCOS(-C(110)/2.+(J-1)*DLATM)                                 VDRIFT.79
       DO I=1,IMAXMP                                                      VDRIFT.80
        EMX(I,J)=-(PHI3D(I+1,J,2)-PHI3D(I-1,J,2))/                        VDRIFT.81
     1           (2*DLONM*CSTH0)*RCOS0S(J)                                VDRIFT.82
       ENDDO                                                              VDRIFT.83
      ENDDO                                                               VDRIFT.84
C     ****                                                                VDRIFT.86
C     ****     Polar values of EMX                                        VDRIFT.87
C     ****                                                                VDRIFT.88
      IMXM3Q=(3*IMAXM)/4                                                  VDRIFT.89
      IMXM1Q=IMAXM/4                                                      VDRIFT.90
      DO I=1,IMAXMP                                                       VDRIFT.92
       EMX(I,1)=EMY(1+MOD(I-1+IMXM1Q,IMAXM),1)                            VDRIFT.93
       EMX(I,JMAXM)=EMY(1+MOD(I-1+IMXM3Q,IMAXM),JMAXM)                    VDRIFT.94
      ENDDO                                                               VDRIFT.96
C     ****                                                                VDRIFT.97
C     ****     Now calculate EMZ = d(PHI)/dz                              VDRIFT.98
C     ****                                                                VDRIFT.99
      IF(IN.NE.0)THEN
       DO J=1,JMAXM                                                       VDRIFT.101
        DO I=1,IMAXMP                                                     VDRIFT.102
         EMZ(I,J)=0.
        ENDDO                                                             VDRIFT.104
       ENDDO                                                              VDRIFT.105
      ELSE                                                                VDRIFT.100
       DO J=1,JMAXM                                                       VDRIFT.101
        DO I=1,IMAXMP                                                     VDRIFT.102
         EMZ(I,J)=-(PHI3D(I,J,3)-PHI3D(I,J,1))                            VDRIFT.103
        ENDDO                                                             VDRIFT.104
       ENDDO                                                              VDRIFT.105
      ENDIF                                                               VDRIFT.106
C     ****                                                                VDRIFT.107
C     ****     Transform EMX, EMY, EMZ to geographic space.               VDRIFT.108
C     ****                                                                VDRIFT.109
C     ****      For each geographic grid point, determine the             VDRIFT.110
C     ****      corresponding theta0 in the distorted geomagnetic         VDRIFT.111
C     ****      latitude grid.                                            VDRIFT.112
C     ****                                                                VDRIFT.113
C     ****      ALATM(I,J) is the thetas (standard apex latitude)         VDRIFT.114
C     ****      corresponding to geographic grid point (I,J).  Find the   VDRIFT.115
C     ****      slot frr this value in the table of theta0 vs.thetas,     VDRIFT.116
C     ****      which has 1 degree spacing from 0. to 90. degrees.        VDRIFT.117
C     ****                                                                VDRIFT.118
      DO J=1,JMAXG                                                        VDRIFT.119
       DO I=1,IMAXG                                                       VDRIFT.120
        PSLOT(I)=2.*DABS(ALATM(I,J))*180./C(110)+1.                           VDRIFT.121
        ISLOT(I)=PSLOT(I)                                                 VDRIFT.122
        PSLOT(I)=PSLOT(I)-DBLE(ISLOT(I))                                  VDRIFT.123
C       ****                                                              VDRIFT.124
C       ****     Now interpolate for theta0 in table using linear         VDRIFT.125
C       ****     interpolation.                                           VDRIFT.126
C       ****                                                              VDRIFT.127
        THETAM(I)=DSIGN((1.-PSLOT(I))*TABLE(ISLOT(I),2)+                  VDRIFT.128
     1            PSLOT(I)*TABLE(ISLOT(I)+1,2),ALATM(I,J))                VDRIFT.129
C       ****                                                              VDRIFT.130
C       ****     Locate magnetic grid element containing geographic       VDRIFT.131
C       ****     grid point (I,J)                                         VDRIFT.132
C       ****                                                              VDRIFT.133
        PSLOT(I)=(ALONM(I,J)+C(110))/DLONM+1.                             VDRIFT.134
        ISLOT(I)=PSLOT(I)                                                 VDRIFT.135
        PSLOT(I)=PSLOT(I)-DBLE(ISLOT(I))                                  VDRIFT.136
        QSLOT(I)=(THETAM(I)+C(110)/2.)/DLATM+1.                           VDRIFT.137
        JSLOT(I)=QSLOT(I)                                                 VDRIFT.138
        QSLOT(I)=QSLOT(I)-DBLE(JSLOT(I))                                  VDRIFT.139
       ENDDO                                                              VDRIFT.140
       DO I=1,IMAXG                                                       VDRIFT.142
C       ****                                                              VDRIFT.143
C       ****     Transform EMX, EMY, EMZ to geographic frame as           VDRIFT.144
C       ****     EX, EY and EZ                                            VDRIFT.145
C       ****                                                              VDRIFT.146
        EEX(I,J)=(1.-PSLOT(I))*(1.-QSLOT(I))*EMX(ISLOT(I),JSLOT(I))       VDRIFT.147
     1             +PSLOT(I)*(1.-QSLOT(I))*EMX(ISLOT(I)+1,JSLOT(I))       VDRIFT.148
     2             +PSLOT(I)*QSLOT(I)*EMX(ISLOT(I)+1,JSLOT(I)+1)          VDRIFT.149
     3             +(1.-PSLOT(I))*QSLOT(I)*EMX(ISLOT(I),JSLOT(I)+1)       VDRIFT.150
        EEY(I,J)=(1.-PSLOT(I))*(1.-QSLOT(I))*EMY(ISLOT(I),JSLOT(I))       VDRIFT.151
     1             +PSLOT(I)*(1.-QSLOT(I))*EMY(ISLOT(I)+1,JSLOT(I))       VDRIFT.152
     2             +PSLOT(I)*QSLOT(I)*EMY(ISLOT(I)+1,JSLOT(I)+1)          VDRIFT.153
     3             +(1.-PSLOT(I))*QSLOT(I)*EMY(ISLOT(I),JSLOT(I)+1)       VDRIFT.154
       ENDDO                                                              VDRIFT.155
       EEX(IMAXG+1,J)=EEX(1,J)                                            VDRIFT.156
       EEY(IMAXG+1,J)=EEY(1,J)                                            VDRIFT.157
       DO I=1,IMAXG                                                       VDRIFT.160
        EEZ(I,J)=(1.-PSLOT(I))*(1.-QSLOT(I))*EMZ(ISLOT(I),JSLOT(I))       VDRIFT.161
     2           +PSLOT(I)*(1.-QSLOT(I))*EMZ(ISLOT(I)+1,JSLOT(I))         VDRIFT.162
     3           +PSLOT(I)*QSLOT(I)*EMZ(ISLOT(I)+1,JSLOT(I)+1)            VDRIFT.163
     4           +(1.-PSLOT(I))*QSLOT(I)*EMZ(ISLOT(I),JSLOT(I)+1)         VDRIFT.164
       ENDDO                                                              VDRIFT.165
       EEZ(IMAXG+1,J)=EEZ(1,J)                                            VDRIFT.166
      ENDDO                                                               VDRIFT.167
      RETURN                                                              VDRIFT.169
      END                                                                 VDRIFT.170
C                                                                         VDRIFT.171
*DK IDRIFT                                                               VDRIFT2.1
      SUBROUTINE IDRIFT(EGX,EGY,EGZ1,EGZ2,IN,HHT,KM)                     VDRIFT2.2
C     ****                                                               VDRIFT2.3
C     ****     Copy ion drift velocities for latitude J from arrays UI,  VDRIFT2.4
C     ****     VI and WI in /DYNPHI/ to slots NUI, NVI and NWI in the    VDRIFT2.5
C     ****     F-array.                                                  VDRIFT2.6
C     ****                                                               VDRIFT2.7
      INCLUDE "param.h"                                                  VDRIFT2.8
      INCLUDE 'epara.h'                                                  VDRIFT2.9
      INCLUDE "buffe.h"                                                  VDRIFT2.10
      INCLUDE "const.h"                                                  VDRIFT2.11
      INCLUDE "econs.h"                                                  VDRIFT2.12
      INCLUDE "index.h"                                                  VDRIFT2.15
      INCLUDE 'efild.h'
      INCLUDE "phycs.h"                                                  VDRIFT2.16
      INCLUDE 'nestcnt.h'
      DIMENSION EGX(IMAXGP,JMAXG),EGY(IMAXGP,JMAXG),EGZ1(IMAXGP,JMAXG),
     1          EGZ2(IMAXGP,JMAXG),HHT(IMAXGP,JMAXG,4)                   VDRIFT2.17
      DIMENSION EX(IMAXGP,JMAXG),EY(IMAXGP,JMAXG),EZ(IMAXGP,JMAXG)       VDRIFT2.17
      COMMON/IOND/UION(IMAXGP,JMAXG),VION(IMAXGP,JMAXG),
     1          WION(IMAXGP,JMAXG)
      DO J=1,JMAXG
       DO I=1,IMAXGP                                                     VDRIFT2.21
C       ****                                                             VDRIFT2.22
C       ****     For latitude J, rotate EX and EY to geographic          VDRIFT2.23
C       ****     orientation using Jacobian.  Divide by distance from    VDRIFT2.24
C       ****     center of earth.                                        VDRIFT2.25
C       ****                                                             VDRIFT2.26
C       ****      EEX = rotated EX                                       VDRIFT2.27
C       ****      EEY = rotated EY                                       VDRIFT2.28
C       ****                                                             VDRIFT2.29
        EX(I,J)=(RJAC(I,J,1,1)*EGX(I,J)+                                 VDRIFT2.30
     1           RJAC(I,J,2,1)*EGY(I,J))/(R00+HHT(I,J,2))                VDRIFT2.31
        EY(I,J)=(RJAC(I,J,1,2)*EGX(I,J)+                                 VDRIFT2.32
     1           RJAC(I,J,2,2)*EGY(I,J))/(R00+HHT(I,J,2))                VDRIFT2.33
       ENDDO                                                             VDRIFT2.34
      ENDDO                                                              VDRIFT2.35
C     ****                                                               VDRIFT2.36
C     ****     For K = 2,KMAX divide EZ by (Z(K+1) - Z(K-1))             VDRIFT2.37
C     ****                                                               VDRIFT2.38
      IF(IN.EQ.1)THEN
       DO J=1,JMAXG
        DO I=1,IMAXGP
         EZ1=EGZ1(I,J)/(HHT(I,J,4)-HHT(I,J,2))
         EZ2=EGZ2(I,J)/(HHT(I,J,3)-HHT(I,J,1))
         EZ(I,J)=2.*EZ2-EZ1           
       ENDDO
       ENDDO
      ELSE IF(IN.EQ.-1)THEN
       DO J=1,JMAXG
        DO I=1,IMAXGP
         EZ1=EGZ1(I,J)/(HHT(I,J,4)-HHT(I,J,2))
         EZ2=EGZ2(I,J)/(HHT(I,J,3)-HHT(I,J,1))
         EZ(I,J)=2.*EZ1-EZ2   
        ENDDO
       ENDDO
      ELSE IF(IN.EQ.2)THEN
       DO J=1,JMAXG
        DO I=1,IMAXGP
         EZ(I,J)=EGZ1(I,J)/(HHT(I,J,3)-HHT(I,J,1))
        ENDDO
       ENDDO
      ELSE 
       DO J=1,JMAXG
        DO I=1,IMAXGP
         EZ(I,J)=EGZ1(I,J)/(HHT(I,J,4)-HHT(I,J,2))
        ENDDO
       ENDDO
      ENDIF
C     ****                                                               VDRIFT2.53
C     ****     VI = (E X B/B**2)                                         VDRIFT2.54
C     ****     Multiply by 1.E6 for results in m/sec                     VDRIFT2.55
C     ****                                                               VDRIFT2.56
      DO J=1,JMAXG                                                       VDRIFT2.60
       DO I=1,IMAXGP                                                     VDRIFT2.64
C       ****                                                             VDRIFT2.65
C       ****     ui = x-component of ion drift velocity                  VDRIFT2.66
C       ****                                                             VDRIFT2.67
        UION(I,J)=-(EY(I,J)*ZZB(I,J)+EZ(I,J)*XB(I,J))*1.D6/              VDRIFT2.68
     1              BMOD(I,J)**2                                         VDRIFT2.69
C       ****                                                             VDRIFT2.70
C       ****     vi = y-component of ion drift velocity                  VDRIFT2.71
C       ****                                                             VDRIFT2.72
        VION(I,J)=(EZ(I,J)*YB(I,J)+EX(I,J)*ZZB(I,J))*1.D6/               VDRIFT2.73
     1               BMOD(I,J)**2                                        VDRIFT2.74
C       ****                                                             VDRIFT2.75
C       ****     wi = Z-component of ion drift velocity                  VDRIFT2.76
C       ****                                                             VDRIFT2.77
        WION(I,J)=(EX(I,J)*XB(I,J)-EY(I,J)*YB(I,J))*1.D6/                VDRIFT2.78
     1             BMOD(I,J)**2                                          VDRIFT2.79
       ENDDO                                                             VDRIFT2.80
      ENDDO                                                              VDRIFT2.81
C     ****                                                               VDRIFT2.53
C     ****     SLOT ION DRIFT VELOCITIES INTO BOTH GLOBAL AND NESTED     VDRIFT2.54
C     ****     GRIDS                                                     VDRIFT2.55
C     ****                                                               VDRIFT2.56
      CALL VISLOT(KM)
      RETURN                                                             VDRIFT2.98
      END                                                                VDRIFT2.99
C                                                                        VDRIFT2.100
*DK VSLOT
      SUBROUTINE VISLOT(KK)
      INCLUDE 'param.h'
      INCLUDE 'epara.h'
      INCLUDE 'enest.h'
      INCLUDE 'const.h'
      INCLUDE 'nestcnt.h'
      COMMON/IOND/UI(IMAXGP,JMAXG),VI(IMAXGP,JMAXG),WI(IMAXGP,JMAXG)      
      COMMON/IONDRIFTG/SUG(IMXG,ZKMXP,JMXG),SVG(IMXG,ZKMXP,JMXG),
     1                 SWG(IMXG,ZKMXP,JMXG)
      COMMON/IONDRIFT1/SU1(IMXN1,ZKMXP,JMXN1),SV1(IMXN1,ZKMXP,JMXN1),
     1                 SW1(IMXN1,ZKMXP,JMXN1)
      COMMON/IONDRIFT2/SU2(IMXN2,ZKMXP,JMXN2),SV2(IMXN2,ZKMXP,JMXN2),
     1                 SW2(IMXN2,ZKMXP,JMXN2)
C     ****
C     ****     FIRST PUT ION DRIFT INTO GLOBAL GRID
C     ****
      DO J=1,JMAX
       JJ=5+(J-1)*9
       DO I=1,IMXG
        II=I+(I-1)*5
        SUG(I,KK,J)=UI(II,JJ)
        SVG(I,KK,J)=VI(II,JJ)
        SWG(I,KK,J)=WI(II,JJ)
       ENDDO
      ENDDO
C     ****
C     ****     THEN PUT ION DRIFT INTO FIRST LEVEL OF NESTED GRID
C     ****
      J1=5+(NESTJ1(1)-1)*9
      DO J=1,JMXN1
       JJ=J1+(J-1)*3
       DO I=1,IMXN1
        PHI=-180.+DBLE(NESTI1(1)-1)*5.+DBLE(I-1)*5./3.
        PHI=PHI+180.
        ISLOT=PHI
        QSLOT=PHI-ISLOT
        SU1(I,KK,J)=QSLOT*UI(ISLOT,JJ)+(1.-QSLOT)*UI(ISLOT+1,JJ)
        SV1(I,KK,J)=QSLOT*VI(ISLOT,JJ)+(1.-QSLOT)*VI(ISLOT+1,JJ)
        SW1(I,KK,J)=QSLOT*WI(ISLOT,JJ)+(1.-QSLOT)*WI(ISLOT+1,JJ)
       ENDDO
      ENDDO
C     ****
C     ****     FINALLY PUT ION DRIFT INTO SECOND LEVEL OF NESTED GRID
C     ****
      J2=5+(NESTJ1(1)-1)*9+(NESTJ1(2)-1)*3
      DO J=1,JMXN2
       JJ=J2+J-1
       DO I=1,IMXN2
        PHI=-180.+DBLE(NESTI1(1)-1)*5.+DBLE(NESTI1(2)-1)*5./3.+
     1      DBLE(I-1)*5./9.
        PHI=PHI+180.
        ISLOT=PHI
        QSLOT=PHI-ISLOT
        SU2(I,KK,J)=QSLOT*UI(ISLOT,JJ)+(1.-QSLOT)*UI(ISLOT+1,JJ)
        SV2(I,KK,J)=QSLOT*VI(ISLOT,JJ)+(1.-QSLOT)*VI(ISLOT+1,JJ)
        SW2(I,KK,J)=QSLOT*WI(ISLOT,JJ)+(1.-QSLOT)*WI(ISLOT+1,JJ)
       ENDDO
      ENDDO      
      RETURN
      END
C
*DK N1VIINTER
      SUBROUTINE N1VIINTER(LEVEL)
      INCLUDE 'n1param.h'
      INCLUDE 'n1const.h'
      INCLUDE 'n1array.h'
      INCLUDE 'index.h'
      INCLUDE 'n1phycs.h'
      INCLUDE 'nestcnt.h'
      COMMON/IONDRIFT/SU(IMX,KMXP,JMX),SV(IMX,KMXP,JMX),
     1                SW(IMX,KMXP,JMX)
      CALL N1UIVI(SU,SV,SW,LEVEL)     
      RETURN
      END
*DK N1UIVI
      SUBROUTINE N1UIVI(U,V,W,LEVEL)
      INCLUDE 'param.h'
      INCLUDE 'const.h'
      INCLUDE 'index.h'
      INCLUDE 'buffe.h'
      INCLUDE 'n1nestg.h'
      INCLUDE 'nestcnt.h'
      INCLUDE 'phycs.h'
      DIMENSION U(NNX,ZKMXP,NJX),V(NNX,ZKMXP,NJX),W(NNX,ZKMXP,NJX)
      NUIK=NUI                                                            HEELIS.21
      NVIK=NVI                                                            HEELIS.22
      NWIK=NWI                                                            HEELIS.23
      JJ=(J-NESTJ1(LEVEL))*NESTRI(LEVEL)+1
      RATIO=DBLE(NESTRI(LEVEL))
      DO K=1,KMAXP1
       DO I=1,INX
        II=NESTRI(LEVEL)*(I-1)+1
        IC=NESTI1(LEVEL)+I-1
        IF(I.LT.INX)THEN
         U(II,K,JJ)=F(IC,NUIK)
         V(II,K,JJ)=F(IC,NVIK)
         W(II,K,JJ)=F(IC,NWIK)
         DO L=1,NESTRI(LEVEL)-1
          IL=II+L
          U(IL,K,JJ)=F(IC,NUIK)+
     1                DBLE(L)*(F(IC+1,NUIK)-F(IC,NUIK))/RATIO
          V(IL,K,JJ)=F(IC,NVIK)+
     1                DBLE(L)*(F(IC+1,NVIK)-F(IC,NVIK))/RATIO
          W(IL,K,JJ)=F(IC,NWIK)+
     1                DBLE(L)*(F(IC+1,NWIK)-F(IC,NWIK))/RATIO
         ENDDO
        ELSE
         U(II,K,JJ)=F(IC,NUIK)
         V(II,K,JJ)=F(IC,NVIK)
         W(II,K,JJ)=F(IC,NWIK)
        ENDIF         
       ENDDO
       NUIK=NUIK+1                                                            HEELIS.21
       NVIK=NVIK+1                                                            HEELIS.22
       NWIK=NWIK+1
      ENDDO       
      RETURN
      END
*DK N1IONDRIFT
      SUBROUTINE N1IONDRIFT(LEVEL)
      INCLUDE 'n1param.h'
      INCLUDE 'n1array.h' 
      INCLUDE 'n1const.h'
      INCLUDE 'n1nestg.h'
      INCLUDE 'nestcnt.h'
      COMMON/IONDRIFT/SU(IMX,KMXP,JMX),SV(IMX,KMXP,JMX),
     1                SW(IMX,KMXP,JMX)
      RATIO=DBLE(NESTRI(LEVEL))
      DO J=1,JNX-1
       JJ=NESTRI(LEVEL)*(J-1)+1
       DO K=1,KMAXP1
        DO I=1,IMAX
         DO L=1,NESTRI(LEVEL)-1
          JL=JJ+L
          SU(I,K,JL)=SU(I,K,JJ)+
     1               DBLE(L)*(SU(I,K,JJ+1)-SU(I,K,JL))/RATIO
          SV(I,K,JL)=SV(I,K,JJ)+
     1               DBLE(L)*(SV(I,K,JJ+1)-SV(I,K,JL))/RATIO
          SW(I,K,JL)=SW(I,K,JJ)+
     1               DBLE(L)*(SW(I,K,JJ+1)-SW(I,K,JL))/RATIO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      RETURN
      END
*DK N1VIINTER
      SUBROUTINE N2VIINTER(LEVEL)
      INCLUDE 'n2param.h'
      INCLUDE 'n2const.h'
      INCLUDE 'index.h'
      COMMON/IONDRIFT2/SU(IMX,KMXP,JMX),SV(IMX,KMXP,JMX),
     1                SW(IMX,KMXP,JMX)
      write(6,*)'fff'
      CALL N2UIVI(SU,SV,SW,LEVEL)     
      RETURN
      END
*DK N1UIVI
      SUBROUTINE N2UIVI(U,V,W,LEVEL)
      INCLUDE 'n1param.h'
      INCLUDE 'n1array.h'
      INCLUDE 'n1const.h'
      INCLUDE 'index.h'
      INCLUDE 'buffe.h'
      INCLUDE 'n2nestg.h'
      INCLUDE 'nestcnt.h'
      INCLUDE 'n1phycs.h'
      DIMENSION U(NNX,KMXP,NJX),V(NNX,KMXP,NJX),W(NNX,KMXP,NJX)
      NUIK=NUI                                                            HEELIS.21
      NVIK=NVI                                                            HEELIS.22
      NWIK=NWI                                                            HEELIS.23
      JJ=(J-NESTJ1(LEVEL))*NESTRI(LEVEL)+1
      RATIO=DBLE(NESTRI(LEVEL))
      DO K=1,KMAXP1
       DO I=1,INX
        II=NESTRI(LEVEL)*(I-1)+1
        IC=NESTI1(LEVEL)+I-1
        IF(I.LT.INX)THEN
         U(II,K,JJ)=F(IC,NUIK)
         V(II,K,JJ)=F(IC,NVIK)
         W(II,K,JJ)=F(IC,NWIK)
         DO L=1,NESTRI(LEVEL)-1
          IL=II+L
          U(IL,K,JJ)=F(IC,NUIK)+
     1                DBLE(L)*(F(IC+1,NUIK)-F(IC,NUIK))/RATIO
          V(IL,K,JJ)=F(IC,NVIK)+
     1                DBLE(L)*(F(IC+1,NVIK)-F(IC,NVIK))/RATIO
          W(IL,K,JJ)=F(IC,NWIK)+
     1                DBLE(L)*(F(IC+1,NWIK)-F(IC,NWIK))/RATIO
         ENDDO
        ELSE
         U(II,K,JJ)=F(IC,NUIK)
         V(II,K,JJ)=F(IC,NVIK)
         W(II,K,JJ)=F(IC,NWIK)
        ENDIF         
       ENDDO
       NUIK=NUIK+1                                                            HEELIS.21
       NVIK=NVIK+1                                                            HEELIS.22
       NWIK=NWIK+1
      ENDDO       
      RETURN
      END
*DK N1IONDRIFT
      SUBROUTINE N2IONDRIFT(LEVEL)
      INCLUDE 'n2param.h'
      INCLUDE 'n2array.h' 
      INCLUDE 'n2const.h'
      INCLUDE 'n2nestg.h'
      INCLUDE 'nestcnt.h'
      COMMON/IONDRIFT2/SU(IMX,KMXP,JMX),SV(IMX,KMXP,JMX),
     1                SW(IMX,KMXP,JMX)
      write(6,*)ttt
      RATIO=DBLE(NESTRI(LEVEL))
      DO J=1,JNX-1
       JJ=NESTRI(LEVEL)*(J-1)+1
       DO K=1,KMAXP1
        DO I=1,IMAX
         DO L=1,NESTRI(LEVEL)-1
          JL=JJ+L
          SU(I,K,JL)=SU(I,K,JJ)+
     1               DBLE(L)*(SU(I,K,JJ+1)-SU(I,K,JL))/RATIO
          SV(I,K,JL)=SV(I,K,JJ)+
     1               DBLE(L)*(SV(I,K,JJ+1)-SV(I,K,JL))/RATIO
          SW(I,K,JL)=SW(I,K,JJ)+
     1               DBLE(L)*(SW(I,K,JJ+1)-SW(I,K,JL))/RATIO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      RETURN
      END
*DK LAMDAIN
      SUBROUTINE LAMDAIN
      INCLUDE 'param.h'
      INCLUDE 'const.h'
      INCLUDE 'edyna.h'
      INCLUDE 'enest.h'
      INCLUDE 'epara.h'
      COMMON/N1DYNAMO/SIGMA1N1(IMXN1,JMXN1,ZKMX),
     1                SIGMA2N1(IMXN1,JMXN1,ZKMX),                                      dynamo.h.4
     1                ZN1(IMXN1,JMXN1,ZKMX),UN1(IMXN1,JMXN1,ZKMX),                     dynamo.h.5
     2                VN1(IMXN1,JMXN1,ZKMX),WN1(IMXN1,JMXN1,ZKMX)                       dyna
      COMMON/N2DYNAMO/SIGMA1N2(IMXN2,JMXN2,ZKMX),
     1                SIGMA2N2(IMXN2,JMXN2,ZKMX),                                      dynamo.h.4
     1                ZN2(IMXN2,JMXN2,ZKMX),UN2(IMXN2,JMXN2,ZKMX),                     dynamo.h.5
     2                VN2(IMXN2,JMXN2,ZKMX),WN2(IMXN2,JMXN2,ZKMX)                       dynamo.h.6
      DIMENSION SIGMA1M(IMAXGP,0:JMAXGP),SIGMA2M(IMAXGP,0:JMAXGP),
     1          ZM(IMAXGP,0:JMAXGP),UM(IMAXGP,0:JMAXGP),
     2          VM(IMAXGP,0:JMAXGP),WM(IMAXGP,0:JMAXGP),
     3          ZZ(IMAXGP,0:JMAXGP,-2:2),WK1(IMXN1,JMXN1),
     4          WK2(IMXN2,JMXN2)
      DATA RL1,RL2,Z0/5.D5,3.D5,90.D5/                                  TRANSFR2.61
C     ****
C     ****     FIRST INTERPOLATION IN HEIGHT
C     ****
      DO K=1,KMAXP1
       DO J=1,JMAX
C       ****
C       ****     FIRST LINE IN LATITUDE
C       ****
        IF(J.EQ.1)THEN
         JJ=5+9*(J-1)
         DO I=1,IMAX+1
          II=5*(I-1)+1
          IF(I.LT.IMAX+1)THEN
           ZM(II,JJ)=Z(I,J,K)
           DO L=1,4
            IL=II+L
            ZM(IL,JJ)=Z(I,J,K)+DBLE(L)*(Z(I+1,J,K)-Z(I,J,K))/5.
           ENDDO
          ELSE
           ZM(II,JJ)=Z(I,J,K)
          ENDIF
         ENDDO
        ELSE
C        ****
C        ****     NEXT LINE IN LATITUDE
C        ****
         J2=5+9*(J-1)
         J1=5+9*(J-2)
         DO I=1,IMAX+1
          II=5*(I-1)+1
          IF(I.LT.IMAX+1)THEN
           ZM(II,J2)=Z(I,J,K)
           DO L=1,4
            IL=II+L
            ZM(IL,J2)=Z(I,J,K)+DBLE(L)*(Z(I+1,J,K)-Z(I,J,K))/5.
           ENDDO
          ELSE
           ZM(II,J2)=Z(I,J,K)
          ENDIF
         ENDDO
C        ****
C        ****     CALCULATE LATITUDINAL LINES IN BETWEEN
C        ****
         DO I=1,IMAXGP
          DO L=1,8
           ZM(I,J2-L)=ZM(I,J2)-DBLE(L)*(ZM(I,J2)-ZM(I,J1))/9.
          ENDDO
         ENDDO
        ENDIF
       ENDDO
C      ****
C      ****     POLES
C      ****
       J1=320
       DO J=1,4
        JJ=J1+J
        DO I=1,IMAXGP
         ZM(I,J)=ZM(I,5)-DBLE(5-J)*(ZM(I,5)-
     1           ZM(1+MOD(I-1+IMAXG/2,IMAXG),5))/DBLE(9)
         ZM(I,JJ)=ZM(I,J1)-DBLE(5-J)*(ZM(I,J1)-
     1            ZM(1+MOD(I-1+IMAXG/2,IMAXG),J1))/DBLE(9)
        ENDDO
       ENDDO
       CALL LAMDAPL(ZM,IMAXGP,JMAXGP)
C      ****
C      ****     PERIODIC POINTS 
C      ****
       DO J=0,JMAXGP
        ZM(IMAXGP,J)=ZM(1,J)
       ENDDO
C      ****
C      ****     SLOT IN FIRST LEVEL NESTED GRID 
C      ****
       DO J=1,JMXN1
        DO I=1,IMXN1
         WK1(I,J)=ZN1(I,J,K)
        ENDDO
       ENDDO
       CALL LAMDAIN1(ZM,WK1)
C      ****
C      ****     SLOT IN SECOND LEVEL NESTED GRID
C      ****
       DO J=1,JMXN2
        DO I=1,IMXN2
         WK2(I,J)=ZN2(I,J,K)
        ENDDO
       ENDDO
       CALL LAMDAIN2(ZM,WK2)
C      ****
C      ****     
C      ****
       IF(K.EQ.1)THEN
        DO KK=-2,0
         DO J=0,JMAXGP
          DO I=1,IMAXGP
           ZZ(I,J,KK)=Z0+DBLE(KK+2)*(ZM(I,J)-Z0)/3.
          ENDDO
         ENDDO
         WRITE(12)((ZZ(I,J,K),I=1,IMAXGP),J=0,JMAXGP)
        ENDDO
        DO J=0,JMAXGP
         DO I=1,IMAXGP
          ZZ(I,J,K)=ZM(I,J)
         ENDDO
        ENDDO
       ENDIF
       IF(K.EQ.2)THEN
        DO J=0,JMAXGP
         DO I=1,IMAXGP
          ZZ(I,J,K)=ZM(I,J)
         ENDDO
        ENDDO
       ENDIF
       WRITE(12)((ZM(I,J),I=1,IMAXGP),J=0,JMAXGP)
      ENDDO
C     ****
C     ****     EXTRA LEVEL FOR HEIGH (K=KMAXP1), ITS VALUES ARE THE SAME AS LEVEL
C     ****     KMAX
C     ****
      WRITE(12)((ZM(I,J),I=1,IMAXGP),J=0,JMAXGP)
C     ****
C     ****     OTHER FIELDS
C     ****
      DO K=1,KMAXP1
       DO J=1,JMAX
C       ****
C       ****     FIRST LINE IN LATITUDE
C       ****
        IF(J.EQ.1)THEN
         JJ=5+9*(J-1)
         DO I=1,IMAX+1
          II=5*(I-1)+1
          IF(I.LT.IMAX+1)THEN
           SIGMA1M(II,JJ)=SIGMA1(I,J,K)
           SIGMA2M(II,JJ)=SIGMA2(I,J,K)
           UM(II,JJ)=U(I,J,K)
           VM(II,JJ)=V(I,J,K)
           WM(II,JJ)=W(I,J,K)
           DO L=1,4
            IL=II+L
            SIGMA1M(IL,JJ)=SIGMA1(I,J,K)+DBLE(L)*(SIGMA1(I+1,J,K)-
     1                     SIGMA1(I,J,K))/5.
            SIGMA2M(IL,JJ)=SIGMA2(I,J,K)+DBLE(L)*(SIGMA2(I+1,J,K)-
     1                     SIGMA2(I,J,K))/5.
            UM(IL,JJ)=U(I,J,K)+DBLE(L)*(U(I+1,J,K)-U(I,J,K))/5.
            VM(IL,JJ)=V(I,J,K)+DBLE(L)*(V(I+1,J,K)-V(I,J,K))/5.
            WM(IL,JJ)=W(I,J,K)+DBLE(L)*(W(I+1,J,K)-W(I,J,K))/5.
           ENDDO
          ELSE
           SIGMA1M(II,JJ)=SIGMA1(I,J,K)
           SIGMA2M(II,JJ)=SIGMA2(I,J,K)
           UM(II,JJ)=U(I,J,K)
           VM(II,JJ)=V(I,J,K)
           WM(II,JJ)=W(I,J,K)
          ENDIF
         ENDDO
        ELSE
C        ****
C        ****     NEXT LINE IN LATITUDE
C        ****
         J2=5+9*(J-1)
         J1=5+9*(J-2)
         DO I=1,IMAX+1
          II=5*(I-1)+1
          IF(I.LT.IMAX+1)THEN
           SIGMA1M(II,J2)=SIGMA1(I,J,K)
           SIGMA2M(II,J2)=SIGMA2(I,J,K)
           UM(II,J2)=U(I,J,K)
           VM(II,J2)=V(I,J,K)
           WM(II,J2)=W(I,J,K)
           DO L=1,4
            IL=II+L
            SIGMA1M(IL,J2)=SIGMA1(I,J,K)+
     1                     DBLE(L)*(SIGMA1(I+1,J,K)-SIGMA1(I,J,K))/5.
            SIGMA2M(IL,J2)=SIGMA2(I,J,K)+
     1                     DBLE(L)*(SIGMA2(I+1,J,K)-SIGMA2(I,J,K))/5.
            UM(IL,J2)=U(I,J,K)+DBLE(L)*(U(I+1,J,K)-U(I,J,K))/5.
            VM(IL,J2)=V(I,J,K)+DBLE(L)*(V(I+1,J,K)-V(I,J,K))/5.
            WM(IL,J2)=W(I,J,K)+DBLE(L)*(W(I+1,J,K)-W(I,J,K))/5.
           ENDDO
          ELSE
           SIGMA1M(II,J2)=SIGMA1(I,J,K)
           SIGMA2M(II,J2)=SIGMA2(I,J,K)
           UM(II,J2)=U(I,J,K)
           VM(II,J2)=V(I,J,K)
           WM(II,J2)=W(I,J,K)
          ENDIF
         ENDDO
C        ****
C        ****     CALCULATE LATITUDINAL LINES IN BETWEEN
C        ****
         DO I=1,IMAXGP
          DO L=1,8
           SIGMA1M(I,J2-L)=SIGMA1M(I,J2)-
     1                     DBLE(L)*(SIGMA1M(I,J2)-SIGMA1M(I,J1))/9.
           SIGMA2M(I,J2-L)=SIGMA2M(I,J2)-
     1                     DBLE(L)*(SIGMA2M(I,J2)-SIGMA2M(I,J1))/9.
           UM(I,J2-L)=UM(I,J2)-DBLE(L)*(UM(I,J2)-UM(I,J1))/9.
           VM(I,J2-L)=VM(I,J2)-DBLE(L)*(VM(I,J2)-VM(I,J1))/9.
           WM(I,J2-L)=WM(I,J2)-DBLE(L)*(WM(I,J2)-WM(I,J1))/9.
          ENDDO
         ENDDO
        ENDIF
       ENDDO
C      ****
C      ****     POLES
C      ****
       J1=320
       DO J=1,4
        JJ=J1+J
        DO I=1,IMAXGP
         SIGMA1M(I,J)=SIGMA1M(I,5)-DBLE(5-J)*(SIGMA1M(I,5)-
     1                SIGMA1M(1+MOD(I-1+IMAXG/2,IMAXG),5))/DBLE(9)
         SIGMA2M(I,J)=SIGMA2M(I,5)-DBLE(5-J)*(SIGMA2M(I,5)-
     1                SIGMA2M(1+MOD(I-1+IMAXG/2,IMAXG),5))/DBLE(9)
         UM(I,J)=UM(I,5)-DBLE(5-J)*(UM(I,5)-
     1           UM(1+MOD(I-1+IMAXG/2,IMAXG),5))/DBLE(9)
         VM(I,J)=VM(I,5)-DBLE(5-J)*(VM(I,5)-
     1           VM(1+MOD(I-1+IMAXG/2,IMAXG),5))/DBLE(9)
         WM(I,J)=WM(I,5)-DBLE(5-J)*(WM(I,5)-
     1           WM(1+MOD(I-1+IMAXG/2,IMAXG),5))/DBLE(9)
         SIGMA1M(I,JJ)=SIGMA1M(I,J1)-DBLE(5-J)*(SIGMA1M(I,J1)-
     1                 SIGMA1M(1+MOD(I-1+IMAXG/2,IMAXG),J1))/DBLE(9)
         SIGMA2M(I,JJ)=SIGMA2M(I,J1)-DBLE(5-J)*(SIGMA2M(I,J1)-
     1                 SIGMA2M(1+MOD(I-1+IMAXG/2,IMAXG),J1))/DBLE(9)
         UM(I,JJ)=UM(I,J1)-DBLE(5-J)*(UM(I,J1)-
     1            UM(1+MOD(I-1+IMAXG/2,IMAXG),J1))/DBLE(9)
         VM(I,JJ)=VM(I,J1)-DBLE(5-J)*(VM(I,J1)-
     1            VM(1+MOD(I-1+IMAXG/2,IMAXG),J1))/DBLE(9)
         WM(I,JJ)=WM(I,J1)-DBLE(5-J)*(WM(I,J1)-
     1            WM(1+MOD(I-1+IMAXG/2,IMAXG),J1))/DBLE(9) 
        ENDDO
       ENDDO
       CALL LAMDAPL(SIGMA1M,IMAXGP,JMAXGP)
       CALL LAMDAPL(SIGMA2M,IMAXGP,JMAXGP)
       CALL LAMDAPL(UM,IMAXGP,JMAXGP)
       CALL LAMDAPL(VM,IMAXGP,JMAXGP)
       CALL LAMDAPL(WM,IMAXGP,JMAXGP)
C      ****
C      ****     PERIODIC POINTS 
C      ****
       DO J=0,JMAXGP
        SIGMA1M(IMAXGP,J)=SIGMA1M(1,J)
        SIGMA2M(IMAXGP,J)=SIGMA2M(1,J)
        UM(IMAXGP,J)=UM(1,J)
        VM(IMAXGP,J)=VM(1,J)
        WM(IMAXGP,J)=WM(1,J)        
       ENDDO
C      ****
C      ****     SLOT IN FIRST LEVEL NESTED GRID 
C      ****
       DO J=1,JMXN1
        DO I=1,IMXN1
         WK1(I,J)=SIGMA1N1(I,J,K)
        ENDDO
       ENDDO
       CALL LAMDAIN1(SIGMA1M,WK1)
       DO J=1,JMXN1
        DO I=1,IMXN1
         WK1(I,J)=SIGMA2N1(I,J,K)
        ENDDO
       ENDDO
       CALL LAMDAIN1(SIGMA2M,WK1)
       DO J=1,JMXN1
        DO I=1,IMXN1
         WK1(I,J)=UN1(I,J,K)
        ENDDO
       ENDDO
       CALL LAMDAIN1(UM,WK1)
       DO J=1,JMXN1
        DO I=1,IMXN1
         WK1(I,J)=VN1(I,J,K)
        ENDDO
       ENDDO
       CALL LAMDAIN1(VM,WK1)
       DO J=1,JMXN1
        DO I=1,IMXN1
         WK1(I,J)=WN1(I,J,K)
        ENDDO
       ENDDO
       CALL LAMDAIN1(WM,WK1)
C      ****
C      ****     SLOT IN SECOND LEVEL NESTED GRID
C      ****
       DO J=1,JMXN2
        DO I=1,IMXN2
         WK2(I,J)=SIGMA1N2(I,J,K)
        ENDDO
       ENDDO
       CALL LAMDAIN2(SIGMA2M,WK1)
       DO J=1,JMXN2
        DO I=1,IMXN2
         WK2(I,J)=SIGMA2N2(I,J,K)
        ENDDO
       ENDDO
       CALL LAMDAIN2(SIGMA2M,WK1)
       DO J=1,JMXN2
        DO I=1,IMXN2
         WK2(I,J)=UN2(I,J,K)
        ENDDO
       ENDDO
       CALL LAMDAIN2(UM,WK1)
       DO J=1,JMXN2
        DO I=1,IMXN2
         WK2(I,J)=VN2(I,J,K)
        ENDDO
       ENDDO
       CALL LAMDAIN2(VM,WK1)
       DO J=1,JMXN2
        DO I=1,IMXN2
         WK2(I,J)=WN2(I,J,K)
        ENDDO
       ENDDO
       CALL LAMDAIN2(WM,WK1)
C      ****
C      ****    ADD THREE EXTRAL LEVELS AT LOWER BOUNDARY
C      ****
       IF(K.EQ.1)THEN
        DO KK=-2,0
         DO J=0,JMAXGP
          DO I=1,IMAXGP       
           SIGMA1M(I,J)=SIGMA1M(I,J)*DEXP((ZZ(I,J,KK)+ZZ(I,J,KK+1)-
     1                  ZZ(I,J,1)-ZZ(I,J,2))/(2.*RL1))
           SIGMA2M(I,J)=SIGMA1M(I,J)*DEXP((ZZ(I,J,KK)+ZZ(I,J,KK+1)-
     1                  ZZ(I,J,1)-ZZ(I,J,2))/(2.*RL2))
          ENDDO
         ENDDO
C        ****
C        ****    WRITE DATA TO A TEMPORAL FILE FOR K=-2,-1 AND 0
C        ****
         DO J=0,JMAXGP
          DO I=1,IMACGP
           WRITE(13)SIGMA1M(I,J),SIGMA2M(I,J),UM(I,J),VM(I,J),WM(I,J)
          ENDDO
         ENDDO
        ENDDO
       ENDIF
C      ****
C      ****    WRITE DATA TO A TEMPORAL FILE FOR K=1 TO KMAX
C      ****
       DO J=0,JMAXGP
        DO I=1,IMACGP
         WRITE(13)SIGMA1M(I,J),SIGMA2M(I,J),UM(I,J),VM(I,J),WM(I,J)
        ENDDO
       ENDDO
      ENDDO
      REWIND(12)
      REWIND(13)
      RETURN
      END
C
*DK LAMDAPL
      SUBROUTINE LAMDAPL(W,IM,JM)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION W(IM,0:JM) 
      TEMPA=0.
      TEMPB=0.
      IMM1=IM-1
      JMM1=JM-1
      DO I=1,IMM1
       TEMPA=TEMPA+W(I,1)
       TEMPB=TEMPB+W(I,2)
      ENDDO
      DO I=1,IM
       W(I,0)=(9.*TEMPA-TEMPB)/(8.*IMM1)
      ENDDO
      TEMPA=0.
      TEMPB=0.
      DO I=1,IMM1
       TEMPA=TEMPA+W(I,JMM1)
       TEMPB=TEMPB+W(I,JMM1-1)
      ENDDO
      DO I=1,IM
       W(I,JM)=(9.*TEMPA-TEMPB)/(8.*IMM1)
      ENDDO       
      RETURN
      END
C
*DK LAMDAIN1
      SUBROUTINE LAMDAIN1(S,W)
      INCLUDE 'param.h'
      INCLUDE 'epara.h'
      INCLUDE 'enest.h'
      INCLUDE 'nestcnt.h'
      DIMENSION S(IMAXGP,0:JMAXGP),W(IMXN1,JMXN1)
C     ****
C     ****    FIND GRID INDEX IN LONGITUDE AND LATITUDE, ONE POINT INSIDE
C     ****    THE NESTED GRID, DISREGARD THE BOUNDARY POINTS
C     ****
      RATIO=5./DBLE(NESTRI(1))
      PHI1=DBLE(NESTI1(1)-3+1)*5.
      IPHI1=PHI1
      PHI2=DBLE(NESTIM(1)-3-1)*5
      IPHI2=PHI2
      J1=5+(NESTJ1(1)-1+1)*9
      JM=5+(NESTJM(1)-2)*9
      DO J=J1,JM
       THETA=-90.+(J-0.5)*5./9.
       PSLOT=(THETA-(-87.5+DBLE(NESTJ1(1)-1)*5.))/RATIO+1.
       JSLOT=PSLOT
       PSLOT=PSLOT-JSLOT
       DO I=IPHI1,IPHI2
        PHI=DBLE(I)
        QSLOT=(PHI-DBLE(NESTI1(1)-3)*5.)/RATIO+1.
        ISLOT=QSLOT
        QSLOT=QSLOT-ISLOT
        S(I,J)=(1.-PSLOT)*(1.-QSLOT)*W(ISLOT,JSLOT)                       VDRIFT.161
     2         +PSLOT*(1.-QSLOT)*W(ISLOT,JSLOT+1)                         VDRIFT.162
     3         +PSLOT*QSLOT*W(ISLOT+1,JSLOT+1)                            VDRIFT.163
     4         +(1.-PSLOT)*QSLOT*W(ISLOT+1,JSLOT)                         VDRIFT.164       
       ENDDO
      ENDDO 
      RETURN
      END
C
*DK LAMDAIN2
      SUBROUTINE LAMDAIN2(S,W)
      INCLUDE 'param.h'
      INCLUDE 'epara.h'
      INCLUDE 'enest.h'
      INCLUDE 'nestcnt.h'
      DIMENSION S(IMAXGP,0:JMAXGP),W(IMXN2,JMXN2)
C     ****
C     ****    FIND GRID INDEX IN LONGITUDE AND LATITUDE, ONE POINT INSIDE
C     ****    THE NESTED GRID, DISREGARD THE BOUNDARY POINTS
C     ****
      RATIO=5./DBLE(NESTRI(1))
      PHI1=DBLE(NESTI1(1)-3)*5.+DBLE(NESTI1(2)-1+1)*RATIO
      IPHI1=PHI1
      PHI2=DBLE(NESTI1(1)-3)*5.+DBLE(NESTIM(2)-2)*RATIO
      IPHI2=PHI2
      J1=5+(NESTJ1(1)-1)*9+(NESTJ1(2)-1+1)*3
      JM=5+(NESTJ1(1)-1)*9+(NESTJM(2)-2)*3
      DO J=J1,JM
       JJ=J+4-J1
       DO I=IPHI1,IPHI2
        PHI=DBLE(I)
        QSLOT=(PHI-PHI1)/RATIO*3+4.
        ISLOT=QSLOT
        QSLOT=QSLOT-ISLOT
       ENDDO
      ENDDO 
      RETURN
      END
C