      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 /sight/ sigp(73,0:37),sigh(73,0:37),qjn(73,0:37),
     | tec(73,0:37),qwind(73,0:37),qamie(73,0:37),wtot(73,0:37),
     | work(73,0:37),fwindu(73,0:37),fwindv(73,0:37),
     | famieu(73,0:37),famiev(73,0:37),ftot(73,0:37),
     | famie(73,0:37),fwind(73,0:37)
      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  tm1(zimxp),tm2(zimxp),tm3(zimxp),sndip(zimxp)
      real :: sigma1_sech(zimxp,zkmxp),sigma2_sech(zimxp,zkmxp)

      DIMENSION RMS(3)
      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
	if (f(i,ntek) .le. 0.) then
	 write (6,"(1x,'neg TE: j i Te=',2i4,e12.4)") j,i,f(i,ntek)
	 call clsgks
	 stop
	endif
      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
	if (v3(i,1,1)*v3(i,2,1) .le. 0. .or. v4(i,1,1)*v4(i,2,1)
     |   .le. 0. .or. v5(i,1,1)*v5(i,2,1) .le. 0.) then
	 write (6,"(1x,'neg lmd1; j i xx,yy,xy=',2i4,6e12.4)") j,i,
     |    v3(i,1,1),v3(i,2,1),v4(i,1,1),v4(i,2,1),v5(i,1,1),v5(i,2,1)
	 call clsgks
	 stop
	endif
      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
	if (v3(i,1,1)*v3(i,2,1) .le. 0. .or. v4(i,1,1)*v4(i,2,1)
     |   .le. 0. .or. v5(i,1,1)*v5(i,2,1) .le. 0.) then
	 write (6,"(1x,'neg lmd2; j i xx,yy,xy=',2i4,6e12.4)") j,i,
     |    v3(i,1,1),v3(i,2,1),v4(i,1,1),v4(i,2,1),v5(i,1,1),v5(i,2,1)
	 call clsgks
	 stop
	endif
      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  Get ht integral of Ped, Hall conductance
C      write (6,"('LAMDAS: ht integration of sigp and sigh for 
C    |        j = ',i3)") j
	      do i=1,imax+1
		sigp(i,j) = 0.
		sigh(i,j) = 0.
	      enddo
              nzk = nj+nz-1
              do k=1,kmax
		nzk = nzk+1
		do 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
                enddo
              enddo
!
! Add sigmas to secondary histories:
!     subroutine addfsech(fname,f,idim1,idim2,ndim2,ilat)
! Field passed to addfsech must have ZIMXP (76) in 1st dimension,
!   and ZKMXP in 2nd dimension.
! Since sigma1,2 have zimxp1 (73) in 1st dimension, so transfer
!   to local arrays sigma1_sech and sigma2_sech w/ periodic points.
! Since ndim2 < idim2 by 1 (ZKMX < ZKMXP) in this case, addfsech will 
!   pad the top level with spval.
! ZIMXP = 76, ZIMXP1 = 72+1 = 73, ZKMX = 28, ZKMXP = ZKMX+1 = 29
! From dynamo.h: SIGMA1(ZIMXP1,ZJMX,ZKMX),SIGMA2(ZIMXP1,ZJMX,ZKMX)
!
      sigma1_sech(1:zimxp1,1:zkmx) = sigma1(:,j,:)
      sigma1_sech(zimxp1+1,:) = sigma1_sech(2,:) ! 74 <- 2
      sigma1_sech(zimxp1+2,:) = sigma1_sech(3,:) ! 75 <- 3
      sigma1_sech(zimxp1+3,:) = sigma1_sech(4,:) ! 76 <- 4
      call addfsech('SIGMA1',sigma1_sech,zimxp,kmaxp1,kmax,j)
      sigma2_sech(1:zimxp1,1:zkmx) = sigma2(:,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
      call addfsech('SIGMA2',sigma2_sech,zimxp,kmaxp1,kmax,j)

C calculate the Joule heating due to wind and ion drift
C qwind = Joule Heting due to wind 
C qamie = Joule Heating due to ion drift
C work = the mechanical work
C wtot = total electromagnetic power
      do i=1,len1
C sndip = sin(dipmag)
                sndip(i) = sin(dipmag(i,j))
        if (abs(sndip(i)) .lt. 1.e-3) 
     |     print *, 'sndip is zero at i= ',i,sndip(i)
C term1=(bx**2 + bz**2)
                tm1(i) = sn2dec(i,j)+(1.-sn2dec(i,j))*sndip(i)**2
C term2=(by**2 + bz**2)
                tm2(i) = (1.-sn2dec(i,j))+sn2dec(i,j)*sndip(i)**2
C term3=bxby
                tm3(i) = sndec(i,j)*csdec(i,j) *(1.-sndip(i)**2)
      enddo

C      write (6,"('LAMDAS: ht integration of qwind and qamie')")

      do i=1,imax+1
	qwind(i,j) = 0.
	qamie(i,j) = 0.
	wtot(i,j) = 0.
	work(i,j) = 0.
	fwindu(i,j) = 0.
	fwindv(i,j) = 0.
	famieu(i,j) = 0.
	famiev(i,j) = 0.
      enddo
      nuik=ndj+nui-1
      nvik=ndj+nvi-1
      nuk=nj+nu-1
      nvk=nj+nv-1
      nzk = nj+nz-1
      do k=1,kmax
        nuik = nuik + 1
        nvik = nvik + 1
        nuk = nuk + 1
        nvk = nvk + 1
        nzk = nzk+1
        do i=1,imax+1
C     write (6,"('LAMDAS: integration of qwind and qamie at 
C    |          i= ',i3)") i
C qwind is the Joule heating associated with nuetral wind terms
C qamie=sigp{(Bx^2+Bz^2)*Wx^2 + (By^2+Bz^2)*Wy^2-2.*Bx*By*Wx*Wy
C        -2.*B^2*Vx*Wx - 2.*B^2*Vy*Wy}
               qwind0 = sigma1(i,j,k)*bmod(i,j)**2*
     +         (tm2(i)*f(i,nuk)**2+tm1(i)*
     +           f(i,nvk)**2-2.*tm3(i)*f(i,nuk)*f(i,nvk)
     +           -2.*(f(i,nuk)*f(i,nuik)+f(i,nvk)*f(i,nvik)))
C qamie is the Joule heating withou nuetral wind
C qamie=sigp{(Bx^2+Bz^2)*Vx^2 + (By^2+Bz^2)*Vy^2+2.*Bx*By*Vx*Vy}/sin(dip)^2
               qamie0 = sigma1(i,j,k)*(tm1(i)*f(i,nuik)**2+tm2(i)*
     +           f(i,nvik)**2+2.*tm3(i)*f(i,nuik)*f(i,nvik))
               qamie0 = qamie0*bmod(i,j)**2 / sndip(i)**2

C wtot = total electric power
C wtot = sigp (E + UxB).E = Qamie + sigp*{(UxB).E}
               wtot0 = - sigma1(i,j,k)*bmod(i,j)**2*(f(i,nuk)*f(i,nuik)
     +           +f(i,nvk)*f(i,nvik)) + sigma2(i,j,k)*bmod(i,j)**2*
     +           (tm3(i)*(f(i,nuk)*
     +           f(i,nuik)-f(i,nvk)*f(i,nvik))-tm1(i)*f(i,nvk)*f(i,nuik)
     +           +tm2(i)*f(i,nuk)*f(i,nvik))/sndip(i)
               wtot0 = qamie0 + wtot0
C Calculate the mechanical work done by wind
               work0 = wtot0 - qwind0 - qamie0

C To convert to mW/m^2 by mutipling by e-11 (since sigp (mho/m),height(cm)
C  B(in Gauss), and velocity (cm/s)
		  qwind(i,j) = qwind(i,j) +
     +              (f(i,nzk+1)-f(i,nzk)) * qwind0 * 1.e-11
		  qamie(i,j) = qamie(i,j) +
     +              (f(i,nzk+1)-f(i,nzk)) * qamie0 * 1.e-11
		  wtot(i,j) = wtot(i,j) +
     +              (f(i,nzk+1)-f(i,nzk)) * wtot0 * 1.e-11
		  work(i,j) = work(i,j) +
     +              (f(i,nzk+1)-f(i,nzk)) * work0 * 1.e-11
C Calculate the horizontal currents
               fwindu0 = -tm1(i)*f(i,nvk)*sigma1(i,j,k)
     +            + tm3(i)*f(i,nuk)*sigma1(i,j,k)
     +            + sndip(i)*f(i,nuk)*sigma2(i,j,k)
               fwindv0 = tm2(i)*f(i,nuk)*sigma1(i,j,k)
     +            + sndip(i)*f(i,nvk)*sigma2(i,j,k)
     +            - tm3(i)*f(i,nvk)*sigma1(i,j,k)
C Convert wind velocity from cm/s to m/s by multipling 1.e-2 and
C Convert B from Gauss to T multipling 1.e-4 and 1.e-2 from cm to m
               fwindu0 = fwindu0*bmod(i,j)*1.e-8/sndip(i)
               fwindv0 = fwindv0*bmod(i,j)*1.e-8/sndip(i)
               famieu0 = -tm1(i)*sigma2(i,j,k)*
     +            f(i,nuik) + sndip(i)*sigma1(i,j,k)*f(i,nvik)
     +            - tm3(i)*sigma2(i,j,k)*f(i,nvik)
               famiev0 = -tm2(i)*sigma2(i,j,k)*
     +            f(i,nvik) - sndip(i)*sigma1(i,j,k)*f(i,nuik)
     +            - tm3(i)*sigma2(i,j,k)*f(i,nuik)
C Convert wind velocity from cm/s to m/s by multipling 1.e-2 and
C Convert B from Gauss to T multipling 1.e-4 and 1.e-2 from cm to m
               famieu0 = famieu0 * bmod(i,j)*1.e-8 / sndip(i)**2
               famiev0 = famiev0 * bmod(i,j)*1.e-8 / sndip(i)**2
		  fwindu(i,j) = fwindu(i,j) +
     +              (f(i,nzk+1)-f(i,nzk)) * fwindu0
		  fwindv(i,j) = fwindv(i,j) +
     +              (f(i,nzk+1)-f(i,nzk)) * fwindv0
		  famieu(i,j) = famieu(i,j) +
     +              (f(i,nzk+1)-f(i,nzk)) * famieu0
		  famiev(i,j) = famiev(i,j) +
     +              (f(i,nzk+1)-f(i,nzk)) * famiev0

        enddo
      enddo
C Set the values at the poles
      sumfaus = 0.
      sumfavs = 0.
      sumfwus = 0.
      sumfwvs = 0.
      sumfaun= 0.
      sumfavn= 0.
      sumfwun= 0.
      sumfwvn= 0.
      do i=1,imax
        sumfwus = sumfwus + fwindu(i,1)
        sumfwun = sumfwun + fwindu(i,36)
        sumfwvs = sumfwvs + fwindv(i,1)
        sumfwvn = sumfwvn + fwindv(i,36)
        sumfaus = sumfaus + famieu(i,1)
        sumfaun = sumfaun + famieu(i,36)
        sumfavs = sumfavs + famiev(i,1)
        sumfavn = sumfavn + famiev(i,36)
      enddo
      do i=1,imax+1
        fwindu(i,0) = sumfwus/72.
        fwindu(i,37) = sumfwun/72.
        fwindv(i,0) = sumfwvs/72.
        fwindv(i,37) = sumfwvn/72.
        famieu(i,0) = sumfaus/72.
        famieu(i,37) = sumfaun/72.
        famiev(i,0) = sumfavs/72.
        famiev(i,37) = sumfavn/72.
      enddo

C     write (6,"('LAMDAS: end of ht integration of qwind and qamie')")

      RETURN
      END
C
