!
      subroutine lamdas(tn,barm,o2,ox,ti,te,o2p,op,nop, lxx,lyy,lxy,lyx,
     |  lamda1,ped_out,hall_out,lev0,lev1,lon0,lon1,lat)
!
! Compute ion drag coefficients 
!      in geographic direction lxx, lyy, lxy, and lyx [1/s] (full levels)
!      in magnetic direction lamda1 [1/s] (full levels)
! Pedersen/ Hall conductivities [S/m] (half levels)
!
! 5/11/06 btf: Merging tiegcm1.8/lamdas (Astrid mods) for timegcm.
!              (lamda1 is now output, passed later to qjoule)
!
      use params_module,only: nlonp4
      use magfield_module,only: bmod2,sn2dec,sncsdc,csdec,sndec
      use cons_module,only: dipmin,rmassinv_o2,rmassinv_o1,rmassinv_n2,
     |  p0,expz,rmass_o2,rmass_o1,boltz,avo
      use magfield_module,only: dipmag
      use input_module,only: colfac
      use addfld_module,only: addfld
      implicit none
!
! Input args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  tn,   ! neutral temperature (deg K)
     |  barm, ! mean molecular mass
     |  o2,   ! molecular oxygen (mmr)
     |  ox,   ! atomic oxygen    (mmr)
     |  ti,   ! ion temperature  (deg K)
     |  te,   ! electron temperature (deg K)
     |  o2p,  ! O2+ number density (cm3)
     |  op,   ! O+ number density (cm3)
     |  nop   ! NO+ number density (cm3)
!
! Output args:
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(out) ::
     |  lxx,     ! lamda XX term
     |  lyy,     ! lamda YY term
     |  lxy,     ! lamda XY term
     |  lyx,     ! lamda YX term
     |  lamda1,  ! sigped*b**2/rho
     |  ped_out, ! Pedersen conductivity (siemens/m)
     |  hall_out ! Hall conductivity (siemens/m)
!
! Local:
      integer :: k,i,l,lonbeg,lonend
      integer ::i0,i1
      real,parameter :: 
     |    qe = 1.602e-19,           ! electron charge (coulomb)
     |    qeomeo10 =  1.7588028E7,  ! qe/m_e/10 [C/g]
     |    qeoNao10 =  9.6489E3      ! qe/N_a/10 [C/mol]
      real,parameter :: 
     |  rmass_op=16.,   ! mass of O+
     |  rmass_o2p=32.,  ! mass of O2+
     |  rmass_nop=30.,  ! mass of NO+
     |  rmi_nop=1./rmass_nop
      real :: sqrt_te   ! sqrt(te)
!
! Local (lon):
      real,dimension(lon0:lon1) :: 
     |  bgauss, qe_fac, sindip, cosdip, cos2dip, sin2dip, cos2dec,
     |  omega_o2p    ,omega_op    ,omega_nop,
     |  omega_o2p_inv,omega_op_inv,omega_nop_inv,
     |  omega_e      ,omega_e_inv
!
! Local (lev,lon):
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  xnmbar,     ! for conversion to volume density
     |  xn2,        ! N2 (1-o2-o mmr)
     |  tnti,       ! average of tn and ti
     |  o2_cm3, ox_cm3, n2_cm3, ! major species number densities (cm3)
     |  sigma_ped,  ! pedersen conductivity (siemens/m)
     |  sigma_hall, ! hall conductivity (siemens/m)
     |  ne,         ! electron density (assume o2p+op+nop)
     |  lamda2,     ! sighal*b**2/rho
     |  lamda1tmp,  ! temporary lamda1
     |  lamda2tmp,  ! temporary lamda2
     |  lxxnorot,   ! XX before rotation
     |  lyynorot,   ! YY before rotation
     |  lxynorot,   ! XY before rotation
     |  lyxnorot,   ! YX before rotation
!
! Ion-neutral momentum transfer collision frequencies:
!
     |  rnu_o2p_o2, ! O2+ ~ O2 collision freq (resonant, temperature dependent)
     |  rnu_op_o2 , ! O+  ~ O2 collision freq (non-resonant)
     |  rnu_nop_o2, ! NO+ ~ O2 collision freq (non-resonant)
!
     |  rnu_o2p_o,  ! O2+ ~ O  collision freq (non-resonant)
     |  rnu_op_o ,  ! O+  ~ O  collision freq (resonant, temperature dependent)
     |  rnu_nop_o,  ! NO+ ~ O  collision freq (non-resonant)
!
     |  rnu_o2p_n2, ! O2+ ~ N2 collision freq (non-resonant)
     |  rnu_op_n2 , ! O+  ~ N2 collision freq (non-resonant)
     |  rnu_nop_n2, ! NO+ ~ N2 collision freq (non-resonant)
!
     |  rnu_o2p,    ! [[o2p~o2]n(o2)+[o2p~o]n(o)+[o2p~n2]n(n2)]/w(o2p)
     |  rnu_op,     ! [[op ~o2]n(o2)+[op ~o]n(o)+[op ~n2]n(n2)]/w(op )
     |  rnu_nop,    ! [[nop~o2]n(o2)+[nop~o]n(o)+[nop~n2]n(n2)]/w(nop)
! 
     |  rnu_ne      ! electron~neutral
!
! Begin execution.
!
      i0 = lon0 ; i1 = lon1
!
! Save input ions to secondary history:
!
!     call addfld('O2P_ION',' ',' ',o2p(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('OP_ION' ,' ',' ',op (lev0:lev1-1,i0:i1),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('NOP_ION',' ',' ',nop(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!
      do i=lon0,lon1
        bgauss(i) = bmod2(i,lat) ! magnetic field strength
!
! 1.e10 for SI units = 1.e6 (cm3->m3) * 1.e4 (gauss->tesla)
!
        qe_fac(i) = qe*1.e10/bgauss(i) 
! 
! gyrofrequencies: omega_i = eB/m_i  [1/s]
!                  omega_e = eB/m_e  [1/s]
! with qeoNao10 = e/Na [C/mol T/kg g/Gauss]
!      qeomeo10 = e/m_e [C/g T/kg g/Gauss]
! 1/10 in qeoNao10 and qeomeo10 for conversion from Gauss/g to T/kg
!
        omega_op (i)    = qeoNao10*bgauss(i)*rmassinv_o1
        omega_o2p(i)    = qeoNao10*bgauss(i)*rmassinv_o2
        omega_nop(i)    = qeoNao10*bgauss(i)*rmi_nop
        omega_op_inv (i)= 1./omega_op(i)
        omega_o2p_inv(i)= 1./omega_o2p(i)
        omega_nop_inv(i)= 1./omega_nop(i)
        omega_e(i)      = qeomeo10*bgauss(i)
        omega_e_inv(i)  = 1./omega_e(i)
!
! Sin and cos for rotation of lamdas:
!   dipmag (magnetic dip angle) is in magfield module (magfield.F), 
!   dipmin (minimum dip angle) is in cons module (cons.F) approx. 10 deg.
!
        if (abs(dipmag(i,lat)) >= dipmin) then
          sindip(i) = sin(dipmag(i,lat))
          cosdip(i) = cos(dipmag(i,lat))
        else
          if (dipmag(i,lat) >= 0.) then
            sindip(i) = sin(dipmin)
            cosdip(i) = cos(dipmin)
          else
            sindip(i) = sin(-dipmin)
            cosdip(i) = cos(-dipmin)
          endif
        endif
        cos2dip(i) = cosdip(i)**2
        sin2dip(i) = sindip(i)**2
        cos2dec(i) = csdec(i,lat)**2
      enddo ! i=lon0,lon1
!
! Ion-neutral momentum transfer collision frequencies [cm^3/s]:
! (Defines "NU" matrix, formerly sub new.F)
!
      do i=lon0,lon1
        do k=lev0,lev1-1
          tnti(k,i) = 0.5*(ti(k,i)+tn(k,i))           ! ave of tn & ti
!
! O2 collision frequencies:
          rnu_o2p_o2(k,i) = 2.59E-11*sqrt(tnti(k,i))* ! O2+ ~ O2 (resonant)
     |      (1.-0.073*alog10(tnti(k,i)))**2
          rnu_op_o2 (k,i) = 6.64E-10                  ! O+  ~ O2
          rnu_nop_o2(k,i) = 4.27E-10                  ! NO+ ~ O2
!
! O collision frequencies:
          rnu_o2p_o(k,i) = 2.31E-10                   ! O2+ ~ O
          rnu_op_o (k,i) = 3.67e-11*sqrt(tnti(k,i))*  ! O+  ~ O (resonant)
     |      (1.-0.064*alog10(tnti(k,i)))**2*colfac    
          rnu_nop_o(k,i) = 2.44E-10                   ! NO+ ~ O
!
! N2 collision frequencies:
          rnu_o2p_n2(k,i) = 4.13E-10                  ! O2+ ~ N2
          rnu_op_n2 (k,i) = 6.82E-10                  ! O+  ~ N2
          rnu_nop_n2(k,i) = 4.34E-10                  ! NO+ ~ N2
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1
!
! Major species number densities [1/cm^3]:
! mass density [g/cm^3/mol]: xnmbar half levels
! pressure at half level: p0*expz(k) [dyn/cm^2]
! mean molecular weight at full level: barm [g/mol]
! Boltzmann constant [erg/K]
! neutral temperature [K]
! mass mixing rations[-]: o2, o1, xn2
! inverse of mean molecular weight [mol/g]:rmassinv_o2, rmassinv_o1, rmassinv_n2
!
      do i=lon0,lon1
        do k=lev0,lev1-1
          xnmbar(k,i) = p0*expz(k)*.5*(barm(k,i)+barm(k+1,i))/
     |      (boltz*tn(k,i))
          o2_cm3(k,i) = o2(k,i)*xnmbar(k,i)*rmassinv_o2
          ox_cm3(k,i) = ox(k,i)*xnmbar(k,i)*rmassinv_o1
          xn2(k,i) = (1.-o2(k,i)-ox(k,i))
          if (xn2(k,i) < 1.e-20) xn2(k,i) = 1.e-20
          n2_cm3(k,i) = xn2(k,i)*xnmbar(k,i)*rmassinv_n2
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1

!     call addfld('O2_CM3',' ',' ',o2_cm3(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('O_CM3' ,' ',' ',ox_cm3(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('N2_CM3',' ',' ',n2_cm3(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!
! Multiply collision freq by neutral number density and sum for each ion:
! nu_in is divided by gyrofrequency omega_i
! nu_in/omega_i [-]:
! rnu_o2p = [[o2p~o2]n(o2)+[o2p~o]n(o)+[o2p~n2]n(n2)]/w(o2p)
! rnu_op  = [[op ~o2]n(o2)+[op ~o]n(o)+[op ~n2]n(n2)]/w(op )
! rnu_nop = [[nop~o2]n(o2)+[nop~o]n(o)+[nop~n2]n(n2)]/w(nop)
!
      do i=lon0,lon1
        do k=lev0,lev1-1
          rnu_o2p(k,i) = (rnu_o2p_o2(k,i)*o2_cm3(k,i) +
     |                    rnu_o2p_o (k,i)*ox_cm3(k,i) +
     |                    rnu_o2p_n2(k,i)*n2_cm3(k,i))*omega_o2p_inv(i)
          rnu_op (k,i) = (rnu_op_o2 (k,i)*o2_cm3(k,i) +
     |                    rnu_op_o  (k,i)*ox_cm3(k,i) +
     |                    rnu_op_n2 (k,i)*n2_cm3(k,i))*omega_op_inv(i)
          rnu_nop(k,i) = (rnu_nop_o2(k,i)*o2_cm3(k,i) +
     |                    rnu_nop_o (k,i)*ox_cm3(k,i) +
     |                    rnu_nop_n2(k,i)*n2_cm3(k,i))*omega_nop_inv(i)
!
! neutral~electron collision frequency (from Banks & Kockards) nu_en
! divided by gyrofrequency omega_2:
! nu_en/omega_e [-]
!
          sqrt_te = sqrt(te(k,i))
          rnu_ne(k,i) = 
     |      (2.33e-11*n2_cm3(k,i)*te(k,i)*(1.-1.21e-4*te(k,i))+
     |       1.82e-10*o2_cm3(k,i)*sqrt_te*(1.+3.60e-2*sqrt_te)+
     |       8.90e-11*ox_cm3(k,i)*sqrt_te*(1.+5.70e-4*te(k,i)))*
     |       omega_e_inv(i)
!
!
! 6/2/06 btf: Multiply rnu_ne by 4, as per Richmond:
! The effective electron-neutral collision frequency is increased in 
! an an hoc manner by a factor of 4 in order for the model to produce
! electric fields and currents below 105 km that agree better with
! observations, as recommended by Gagnepain et al. (J. Atmos. Terr. 
! Phys., 39, 1119-1124, 1977).
!
          rnu_ne(k,i) = rnu_ne(k,i)*4.
!
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1

!     call addfld('RNU_O2P',' ',' ',rnu_o2p(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('RNU_OP' ,' ',' ',rnu_op(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('RNU_NOP',' ',' ',rnu_nop(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('RNU_NE' ,' ',' ',rnu_ne(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)

!
! Pedersen and Hall conductivities (siemens/m):
! Qe_fac includes conversion from CGS to SI units,
!  -> e/B [C/T 10^6 m^3/cm^3], see above.
! number densities [1/cm^3]
!
      do i=lon0,lon1
        do k=lev0,lev1-1
!
! ne = electron density assuming charge equilibrium (cm3):
          ne(k,i) = op(k,i)+o2p(k,i)+nop(k,i) 
!
! Pedersen conductivity (siemens/m):
          sigma_ped(k,i) = qe_fac(i)*
     |      ((op (k,i)*rnu_op (k,i)/(1.+rnu_op (k,i)**2))+
     |       (o2p(k,i)*rnu_o2p(k,i)/(1.+rnu_o2p(k,i)**2))+
     |       (nop(k,i)*rnu_nop(k,i)/(1.+rnu_nop(k,i)**2))+
     |       (ne (k,i)*rnu_ne (k,i)/(1.+rnu_ne (k,i)**2)))
!
! Hall conductivity (siemens/m):
          sigma_hall(k,i) = qe_fac(i)*
     |      (ne (k,i)/(1.+rnu_ne (k,i)**2)-
     |       op (k,i)/(1.+rnu_op (k,i)**2)-
     |       o2p(k,i)/(1.+rnu_o2p(k,i)**2)-
     |       nop(k,i)/(1.+rnu_nop(k,i)**2))

        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1

!     call addfld('ELECDEN' ,' ',' ',ne(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('SIGMAPED',' ',' ',sigma_ped(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('SIGMAHAL',' ',' ',sigma_hall(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!
! ion drag coefficients lamda1,2 [1/s] (full level)
! Pedersen/ Hall conductivity [S/m]: sigma_ped, sigma_hall (half level)
! mass density [g/cm^3/mol]: xnmbar (half level)
! magnetic field [Gauss]: bgauss
! Na Avoagdro number [1/mol]
! conversion: tesla = 1.e-4*bgauss
!             kg/m3 = 1.e-3*g/(1.e-6*cm^3) = 1.e3*g/cm3
!
      do i=lon0,lon1
        do k=lev0,lev1-1
          lamda1tmp(k,i) = (sigma_ped(k,i)*(1.e-4*bgauss(i))**2)*avo/
     |                     (1.e3*xnmbar(k,i))
          lamda2tmp(k,i) = (sigma_hall(k,i)*(1.e-4*bgauss(i))**2)*avo/
     |                     (1.e3*xnmbar(k,i))
        enddo ! k=lev0,lev1-1
!
! Lamda1,2 to full (interface) levels:
        do k=lev0,lev1-2
!
! 5/17/04 btf: FP invalid here on SGI dataproc due to negative lamda1tmp
!              at k==lev1 (reference to lamda1tmp at k+1):
!         write(6,"('lamdas: lat=',i3,' i=',i3,' k=',i3,' lamda1tmp=',
!    |      e12.4,' lamda1tmp(k+1,i)=',e12.4)") lat,i,k,
!    |      lamda1tmp(k,i),lamda1tmp(k+1,i)         

          lamda1(k+1,i) = sqrt(lamda1tmp(k,i)*lamda1tmp(k+1,i))
          lamda2(k+1,i) = sqrt(lamda2tmp(k,i)*lamda2tmp(k+1,i))
        enddo ! k=lev0,lev1-1
! Bottom boundary:
        lamda1(lev0,i) = sqrt(lamda1tmp(lev0,i)**3/lamda1tmp(lev0+1,i))
        lamda2(lev0,i) = sqrt(lamda2tmp(lev0,i)**3/lamda2tmp(lev0+1,i))
! Top boundary:
        lamda1(lev1,i)= sqrt(lamda1tmp(lev1-1,i)**3/lamda1tmp(lev1-2,i))
        lamda2(lev1,i)= sqrt(lamda2tmp(lev1-1,i)**3/lamda2tmp(lev1-2,i))
!
! Non-rotated lamdas:
        do k=lev0,lev1
          lxxnorot(k,i) = lamda1(k,i)
          lyynorot(k,i) = lamda1(k,i)*sin2dip(i)
          lxynorot(k,i) = lamda2(k,i)*sindip(i)
          lyxnorot(k,i) = lxynorot(k,i)
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1

!     call addfld('LAMDA1',' ',' ',lamda1, ! appears to need periodic pts
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('LAMDA2',' ',' ',lamda2,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('LXXNOROT',' ',' ',lxxnorot,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('LYYNOROT',' ',' ',lyynorot,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('LXYNOROT',' ',' ',lxynorot,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('LYXNOROT',' ',' ',lyxnorot,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!
! Rotate lamdas for displacement of geomagnetic and geographic poles:
!
      do i=lon0,lon1
        do k=lev0,lev1
          lxx(k,i)= lxxnorot(k,i)*cos2dec(i)+lyynorot(k,i)*sn2dec(i,lat)
          lyy(k,i)= lyynorot(k,i)*cos2dec(i)+lxxnorot(k,i)*sn2dec(i,lat)
          lyx(k,i)= lxynorot(k,i)-(lyynorot(k,i)-lxxnorot(k,i))*
     |      sndec(i,lat)*csdec(i,lat)
          lxy(k,i)= lxynorot(k,i)+(lyynorot(k,i)-lxxnorot(k,i))*
     |      sndec(i,lat)*csdec(i,lat)
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1
 
!     call addfld('LXX','LXX','Hz',lxx(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('LYY','LYY','Hz',lyy(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('LXY','LXY','Hz',lxy(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('LYX','LYX','Hz',lyx(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('LAMDA1','LAMDA1 ','Hz',lamda1,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)

      lonbeg = lon0
      if (lon0==1) lonbeg = 3
      lonend = lon1
      if (lon1==nlonp4) lonend = nlonp4-2
!
! Output conductivities [S/m] (half levels):
      do i=lonbeg,lonend
        do k=lev0,lev1
          ped_out(k,i) = sigma_ped(k,i)
          hall_out(k,i) = sigma_hall(k,i)
        enddo
      enddo

!     call addfld('SIGMAPED' ,'SIGMAPED','S/m',
!    |  ped_out(lev0:lev1-1,i0:i1),'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('SIGMAHAL','SIGMAHAL','S/m',
!    |  hall_out(lev0:lev1-1,i0:i1),'lev',lev0,lev1-1,'lon',i0,i1,lat)

      end subroutine lamdas
