       SUBROUTINE LAMDAS
C     ****
C     ****     COMPUTE ION DRAG COEFFICIENTS
C     ****
C     ****     MAJOR SPECIES AND IONS ARE NUMBERED THUS
C     ****       O2  -  1            O2+  -  1
C     ****       O   -  2            O+   -  2
C     ****       N2  -  3            N2+  -  3
C     ****
      include "params.h"
      include "blnk.h"
      include "cons.h"
      include "index.h"
      include "buff.h"
      include "phys.h"
      include "dynamo.h"
      COMMON/TRGM/RLATM(ZIMXP,ZJMX),RLONM(ZIMXP,ZJMX),
     1  DIPMAG(ZIMXP,ZJMX),DECMAG(ZIMXP,ZJMX),SNDEC(ZIMXP,ZJMX),
     2  CSDEC(ZIMXP,ZJMX),SN2DEC(ZIMXP,ZJMX),SNCSDC(ZIMXP,ZJMX),
     3  DUMDUM(ZJMX,3)
      COMMON/MAGFLD/BXM(ZIMXP,2),BX(ZIMXP,ZJMX),BXP(ZIMXP,4),
     1  BY(ZIMXP,ZJMX),BYP(ZIMXP,4),BZ(ZIMXP,ZJMX),BZP(ZIMXP,4),
     2  BMOD(ZIMXP,ZJMX),BMODP(ZIMXP,2)
      COMMON/VSCR/RM1(ZIMXP,ZKMXP,3,1),RM2(ZIMXP,ZKMXP,3,1),
     1RM3(ZIMXP,ZKMXP,3,1),V1(ZIMXP,ZKMXP,1),V2(ZIMXP,ZKMXP,1),
     2V3(ZIMXP,ZKMXP,1),V4(ZIMXP,ZKMXP,1),V5(ZIMXP,ZKMXP,1),
     3V6(ZIMXP,ZKMXP,1),T1(ZIMXP),T2(ZIMXP),T3(ZIMXP),T4(ZIMXP),
     4T5(ZIMXP),T6(ZIMXP),T7(ZIMXP)
      DIMENSION RMS(3)
      DIMENSION SIGP(73,0:37),SIGH(73,0:37)
      REAL SIGMA1_SECH(ZIMXP,ZKMXP),SIGMA2_SECH(ZIMXP,ZKMXP)
      DATA RMS/32.,16.,30./,QE/1.602E-19/
      DATA DIPMIN/.20/
C     ****
C     ****     T1 = BGAUSS
C     ****     T2 = 9.6489E3 * BGAUS    ( FOR OMEGA(I) )
C     ****     T3 = 1.7588028E7 * BGAUSS  ( FOR OMEGA(E) )
C     ****     T4 = SIN(DIP)
C     ****
      DO 1 I=1,LEN1
      T1(I)=BMOD(I,J)
      T2(I)=9.6489E3*T1(I)
      T3(I)=1.7588028E7*T1(I)
      T4(I)=SIN(DIPMAG(I,J))
      T5(I) = T4(I)**2
    1 CONTINUE
C     ****
C     ****     V1 = PSI(O2)/M(O2), V2 = PSI(O)/M(O), V1 = PSI2)/M(N2)
C     ****                       AT (K+1/2)
C     ****
      NPS1K=NJ+NPS
      NPS2K=NJ+NPS2
      NTIK=NJ+NTI
      NTK=NJ+NT
      DO 2 I=1,LEN2
      V1(I,1,1)=F(I,NPS1K)/RMASS(1)
      V2(I,1,1)=F(I,NPS2K)/RMASS(2)
      V3(I,1,1)=(1.-F(I,NPS1K)-F(I,NPS2K))/RMASS(3)
C     ****
C     ****     V6 = TR = 0.5*(TI+TN)           ( K+1/2 )
C     ****
      V6(I,1,1)=0.5*(F(I,NTIK)+F(I,NTK))
    2 CONTINUE
C     ****
C     ****     RM1 = NEW(L,M)
C     ****
C     ****     GENERATE NUMERICAL FACTORS
      CALL NEW(V6,RM1)
C     ****     MULTIPLY BY MAJOR SPECIES NUMBER DENSITIES
      DO 3 L=1,3
      DO 3 M=1,3
      DO 3 I=1,LEN2
      RM1(I,1,L,M)=RM1(I,1,L,M)*V1(I,1,M)
    3 CONTINUE
C     ****
C     ****     V1(L) = SUM(M=1,3) RM1(L,M)    ( K+1/2 )
C     ****
      DO 4 L=1,3
      DO 5 I=1,LEN2
      V1(I,1,L)=0.
    5 CONTINUE
      DO 4 M=1,3
      DO 4 I=1,LEN2
      V1(I,1,L)=V1(I,1,L)+RM1(I,1,L,M)
    4 CONTINUE
C     ****
C     ****     RM1(1) = N(O2+), RM1(2) = N(O+), RM1(1) = N(NO+)
C     ****                      ( AT (K+1/2) )
C     ****
      NO2PK=NJ+NO2P
      NOPK=NJ+NOP
      NNOPK=NNOP
      DO 6 I=1,LEN2
      RM1(I,1,1,1)=F(I,NO2PK)
      RM1(I,1,2,1)=F(I,NOPK)
      RM1(I,1,3,1)=F(I,NNOPK)
    6 CONTINUE
C     ****
C     ****     V6 = N * MBAR            ( K+1/2 )
C     ****
      NTK=NJ+NT-1
      NMSK=NJ+NMS-1
      DO 7 K=1,KMAX
      NTK=NTK+1
      NMSK=NMSK+1
      DO 7 I=1,LEN1
      V6(I,K,1)=C(81)*EXPS(K)*.5*(F(I,NMSK)+F(I,NMSK+1))/(C(84)*
     1F(I,NTK))
    7 CONTINUE
C     ****
C     ****     V1(L) = V1(L)*(N*MBAR)/OMEGAI(L)    ( L=1,3 )  ( K=1/2 )
C     ****
      DO 8 L=1,3
      DO 8 K=1,KMAX
      DO 8 I=1,LEN1
      V1(I,K,L)=V1(I,K,L)*V6(I,K,1)*RMS(L)/T2(I)
    8 CONTINUE
C     ****
C     ****     V4 = NU/OMEGA (E)
C     ****
      NPS1K=NJ+NPS
      NPS2K=NJ+NPS2
      NTEK=NJ+NTE
      DO 9 I=1,LEN2
      V5(I,1,1)=SQRT(F(I,NTEK))
      V4(I,1,1)=V6(I,1,1)*(2.33E-11*(1.-F(I,NPS1K)-F(I,NPS2K))/RMASS(3)*
     1F(I,NTEK)*(1.-1.21E-4*F(I,NTEK))+1.82E-10*F(I,NPS1K)/RMASS(1)*
     2V5(I,1,1)*(1.+3.6E-2*V5(I,1,1))+8.9E-11*F(I,NPS2K)/RMASS(2)*
     3V5(I,1,1)*(1.+5.7E-4*F(I,NTEK)))
    9 CONTINUE
C     ****
C     ****     V4 = NU(EN)/OMEGAE
C     ****
      DO 10 K=1,KMAX
      DO 10 I=1,LEN1
      V4(I,K,1)=V4(I,K,1)/T3(I)
   10 CONTINUE
C     ****
C     ****     RM2 = PEDERSEN FACTORS
C     ****     RM3 = HALL FACTORS
C     ****
      DO 11 L=1,3
      DO 11 I=1,LEN2
      RM2(I,1,L,1)=V4(I,1,1)/(1.+V4(I,1,1)**2)+V1(I,1,L)/(1.+V1(I,1,L)
     1**2)
      RM3(I,1,L,1)=1./(1.+V4(I,1,1)**2)-1./(1.+V1(I,1,L)**2)
   11 CONTINUE
C     ****
C     ****     V1 = SUM(L=1,3) PEDERSEN FACTORS * ION DENSITIES
C     ****     V2 = SUM(L=1,3) HALL FACTORS * ION DENSITIES
C     ****
      DO 12 I=1,LEN2
      V1(I,1,1)=0.
      V2(I,1,1)=0.
   12 CONTINUE
      DO 13 L=1,3
      DO 13 I=1,LEN2
      V1(I,1,1)=V1(I,1,1)+RM2(I,1,L,1)*RM1(I,1,L,1)
      V2(I,1,1)=V2(I,1,1)+RM3(I,1,L,1)*RM1(I,1,L,1)
   13 CONTINUE
C     ****
C     ****     SET LXX, LYY, LXY
C     ****
      DO 14 K=1,KMAX
      DO 14 I=1,LEN1
        V3(I,K,1)=V1(I,K,1)*T1(I)*1.E-1*QE*C(85)/V6(I,K,1)
        V4(I,K,1)=V3(I,K,1)*T5(I)
        V5(I,K,1)=V2(I,K,1)*T1(I)*1.E-1*QE*C(85)/V6(I,K,1)
   14 CONTINUE
      NLXXK=NLXX+1
      NLYYK=NLYY+1
      NLXYK=NLXY+1
      DO 15 I=1,LEN2-LEN1
      F(I,NLXXK)=SQRT(V3(I,1,1)*V3(I,2,1))
      F(I,NLYYK)=SQRT(V4(I,1,1)*V4(I,2,1))
      F(I,NLXYK)=SQRT(V5(I,1,1)*V5(I,2,1))
   15 CONTINUE
C     ****
C     ****     BOUNDARIES
C     ****
      NLXXK=NLXX
      NLYYK=NLYY
      NLXYK=NLXY
      DO 16 I=1,LEN1
      F(I,NLXXK)=SQRT(V3(I,1,1)**3/V3(I,2,1))
      F(I,NLYYK)=SQRT(V4(I,1,1)**3/V4(I,2,1))
      F(I,NLXYK)=SQRT(V5(I,1,1)**3/V5(I,2,1))
      F(I,NLXXK+KMAX)=SQRT(V3(I,KMAX,1)**3/V3(I,KMAXM1,1))
      F(I,NLYYK+KMAX)=SQRT(V4(I,KMAX,1)**3/V4(I,KMAXM1,1))
      F(I,NLXYK+KMAX)=SQRT(V5(I,KMAX,1)**3/V5(I,KMAXM1,1))
   16 CONTINUE
      NLXYK = NLXY-1
      DO 19 K =1,KMAXP1
        NLXYK = NLXYK+1
        DO 19 I = 1,IMAXP4
          F(I,NLXYK) = F(I,NLXYK)*T4(I)
   19 CONTINUE
C     ****
C     ****     COPY SIGMA1, SIGMA2, Z, H, U, V AND W TO COMMON/DYNAMO/
C     ****
      NZK = NJ+NZ-1
      NTK = NJ+NT-1
      NMSK = NJ+NMS-1
      NUK = NJ+NU-1
      NVK = NJ+NV-1
      NWK = NJ+NW-1
      DO 17 K = 1,KMAX
        NZK = NZK+1
        NTK = NTK+1
        NMSK = NMSK+1
        NUK = NUK+1
        NVK = NVK+1
        NWK = NWK+1
        DO 17 I = 1,IMAX+1
          SIGMA1(I,J,K) = V1(I+2,K,1)*1.E10*QE/T1(I+2)
          SIGMA2(I,J,K) = V2(I+2,K,1)*1.E10*QE/T1(I+2)
          Z(I,J,K) = F(I+2,NZK)
          H(I,J,K) = C(57)*F(I+2,NTK)/(.5*(F(I+2,NMSK)+F(I+2,NMSK+1))
     1      *C(54))
          U(I,J,K) = F(I+2,NUK)
          V(I,J,K) = F(I+2,NVK)
          W(I,J,K) = .5*(F(I+2,NWK)+F(I+2,NWK+1))*H(I,J,K)
   17 CONTINUE
      NZK = NZK+1
      DO 18 I = 1,IMAX+1
        Z(I,J,KMAXP1) = F(I+2,NZK)
   18 CONTINUE
C
C
C  Get ht integral of Ped, Hall conductance
C      WRITE(6,"('LAMDAS: ht integration of sigp and sigh for
C    |        j = ',i3)") j
      DO 29 I=1,IMAX+1
        SIGP(I,J) = 0.
        SIGH(I,J) = 0.
 29   CONTINUE 
      NZK= NJ+NZ-1
      DO 21 K=1,KMAX
         NZK= NZK+1
         DO 20 I=1,IMAX+1
C  Times ht integ by 1.e-2 because appears off for siemens (1.e-2cm/m?)
           SIGP(I,J) = SIGP(I,J) +
     +              (F(I,NZK+1)-F(I,NZK)) * SIGMA1(I,J,K) * 1.E-2
           SIGH(I,J) = SIGH(I,J) +
     +              (F(I,NZK+1)-F(I,NZK)) * SIGMA2(I,J,K) * 1.E-2
 20   CONTINUE 
 21   CONTINUE 
CCCCCCCCCCCCCCCCCCC save secondary histories CCCCCCCCCCCCCCCCCCCCC
C
C Add SIGMAs to secondary histories:
C     subroutine addfSECH(fname,f,idim1,idim2,ndim2,ilat)
C Field passed to addfSECH must have ZIMXP (76) in 1st dimension,
C   and ZKMXP in 2nd dimension.
C Since SIGMA1,2 have ZIMXP1 (73) in 1st dimension, so transfer
C   to local arrays SIGMA1_SECH and SIGMA2_SECH w/ periodic points.
C Since ndim2 < idim2 by 1 (ZKMX < ZKMXP) in this case, addfSECH will
C   pad the top level with spval.
C ZIMXP = 76, ZIMXP1 = 72+1 = 73, ZKMX = 28, ZKMXP = ZKMX+1 = 29
C From dynamo.h: SIGMA1(ZIMXP1,ZJMX,ZKMX),SIGMA2(ZIMXP1,ZJMX,ZKMX)
C
!     SIGMA1_SECH(1:ZIMXP1,1:ZKMX) = SIGMA1(:,j,:)
!     SIGMA1_SECH(1:ZIMXP1,ZKMXP) = SIGP(:,j)
!     SIGMA1_SECH(ZIMXP1+1,:) = SIGMA1_SECH(2,:) ! 74 <- 2
C
!     SIGMA1_SECH(ZIMXP1+2,:) = SIGMA1_SECH(3,:) ! 75 <- 3
!     SIGMA1_SECH(ZIMXP1+3,:) = SIGMA1_SECH(4,:) ! 76 <- 4
C     CALL ADDFSECH('PEDERSEN',SIGMA1_SECH,ZIMXP,KMAXP1,KMAXP1,J)
!     SIGMA2_SECH(1:ZIMXP1,1:ZKMX) = SIGMA2(:,J,:)
!     SIGMA2_SECH(1:ZIMXP1,ZKMXP) = SIGH(:,J)
!     SIGMA2_SECH(ZIMXP1+1,:) = SIGMA2_SECH(2,:) ! 74 <- 2
!     SIGMA2_SECH(ZIMXP1+2,:) = SIGMA2_SECH(3,:) ! 75 <- 3
!     SIGMA2_SECH(ZIMXP1+3,:) = SIGMA2_SECH(4,:) ! 76 <- 4
C     CALL ADDFSECH('HALL',SIGMA2_SECH,ZIMXP,KMAXP1,KMAXP1,J)
C
      RETURN
      END
C
