#include "defs.h"
      module dyndiag_module
      use params_module,only: nlonp1,nlonp4,nlat,nlevp1,nlatp1
      implicit none
!
      real, dimension(nlonp4,nlat) ::
     |  sigp,    ! height integrated pedersen conductivity
     |  sigh,    ! height integrated hall conductivity
     |  qwind,   ! joule heating associated with neutral winds
     |  qamie,   ! joule heating from amie
     |  wtot,    ! total electromagnetic power
     |  work,    ! mechanical work
     |  fwindu,fwindv, ! horizontal current
     |  famieu,famiev  ! horizontal current from amie
!
      contains
!-----------------------------------------------------------------------
      subroutine dyndiag(sigma1,sigma2,z,un,vn,ui,vi,
     |  lev0,lev1,lon0,lon1,lat)
      use magfield_module,only: dipmag,sn2dec,csdec,sndec,bmod
!
! Calculate dynamo related diagnostics for secondary histories.
! (this sub called from inside lat loop in dynamics, after lamdas)
!
! Args:
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  sigma1,  ! pederson conductivity (from lamdas)
     |  sigma2,  ! hall conductivity (from lamdas)
     |  z,       ! geopotential from addiag
     |  un,vn,   ! neutral velocities
     |  ui,vi    ! ion velocities
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
!
! Local:
      integer :: i,k
      real,dimension(nlonp4) :: tm1,tm2,tm3,sndip 
!
! Get height integral of Ped, Hall conductance
      do i=lon0,lon1
        sigp(i,lat) = 0.
        sigh(i,lat) = 0.
      enddo
      do i=lon0,lon1
        do k=lev0,lev1-1
!
!  Times ht integ by 1.e-2 because appears off for siemens (1.e-2cm/m?)
          sigp(i,lat) = sigp(i,lat) +
     |      (z(k+1,i)-z(k,i)) * sigma1(k,i) * 1.e-2
          sigh(i,lat) = sigh(i,lat) +
     |      (z(k+1,i)-z(k,i)) * sigma2(k,i) * 1.e-2
        enddo
      enddo
!
! Calculate the Joule heating due to wind and ion drift
!  qwind = Joule Heting due to wind
!  qamie = Joule Heating due to ion drift
!  work = the mechanical work
!  wtot = total electromagnetic power
!
      do i=lon0,lon1
        sndip(i) = sin(dipmag(i,lat))
!       if (abs(sndip(i)) .lt. 1.e-3)
!    |     print *, 'sndip is zero at i= ',i,sndip(i)
!
! term1=(bx**2 + bz**2)
! term2=(by**2 + bz**2)
! term3=bxby
!
        tm1(i) = sn2dec(i,lat)+(1.-sn2dec(i,lat))*sndip(i)**2
        tm2(i) = (1.-sn2dec(i,lat))+sn2dec(i,lat)*sndip(i)**2
        tm3(i) = sndec(i,lat)*csdec(i,lat) *(1.-sndip(i)**2)
      enddo
!
! Init:
      qwind(:,lat) = 0.
      qamie(:,lat) = 0.
      wtot(:,lat) = 0.
      work(:,lat) = 0.
      fwindu(:,lat) = 0.
      fwindv(:,lat) = 0.
      famieu(:,lat) = 0.
      famiev(:,lat) = 0.
!
! Longitude, level loops:
      do i=lon0,lon1
        do k=lev0,lev1-1
!
! qwind is the Joule heating associated with nuetral wind terms
! qamie=sigp{(Bx^2+Bz^2)*Wx^2 + (By^2+Bz^2)*Wy^2-2.*Bx*By*Wx*Wy
!        -2.*B^2*Vx*Wx - 2.*B^2*Vy*Wy}
!
          qwind(i,lat) = sigma1(k,i)*bmod(i,lat)**2*(tm2(i)*
     |      un(k,i)**2+tm1(i)*vn(k,i)**2-2.*tm3(i)*un(k,i)*
     |      vn(k,i)-2.*(un(k,i)*ui(k,i)+vn(k,i)*vi(k,i)))
!
! qamie is Joule heating without neutral wind
! qamie=sigp{(Bx^2+Bz^2)*Vx^2+(By^2+Bz^2)*Vy^2+2.*Bx*By*Vx*Vy}/
!       sin(dip)^2
!
          qamie(i,lat) = sigma1(k,i)*(tm1(i)*ui(k,i)**2+tm2(i)*
     |      vi(k,i)**2+2.*tm3(i)*ui(k,i)*vi(k,i))
          qamie(i,lat) = qamie(i,lat)*bmod(i,lat)**2 / sndip(i)**2
!
! wtot = total electric power
! wtot = sigp (E + UxB).E = Qamie + sigp*{(UxB).E}
!
          wtot(i,lat) = - sigma1(k,i)*bmod(i,lat)**2*(un(k,i)*ui(k,i)+
     |      vn(k,i)*vi(k,i)) + sigma2(k,i)*bmod(i,lat)**2*
     |      (tm3(i)*(un(k,i)*ui(k,i)-vn(k,i)*vi(k,i))-tm1(i)*
     |      vn(k,i)*ui(k,i)+tm2(i)*un(k,i)*vi(k,i))/sndip(i)
          wtot(i,lat) = qamie(i,lat) + wtot(i,lat)
!
! To convert to mW/m^2 by mutipling by e-11 (since sigp (mho/m),height(cm)
!  B(in Gauss), and velocity (cm/s)
!
          qwind(i,lat) = qwind(i,lat)+(z(k+1,i)-z(k,i))*qwind(i,lat)*
     |      1.e-11
          qamie(i,lat) = qamie(i,lat)+(z(k+1,i)-z(k,i))*qamie(i,lat)*
     |      1.e-11
          wtot(i,lat) = wtot(i,lat)  +(z(k+1,i)-z(k,i))*wtot(i,lat)*
     |      1.e-11
          work(i,lat) = work(i,lat)  +(z(k+1,i)-z(k,i))*work(i,lat)*
     |      1.e-11
!
! Calculate horizontal currents:
          fwindu(i,lat) = -tm1(i)*vn(k,i)*sigma1(k,i)+tm3(i)*un(k,i)*
     |      sigma1(k,i) + sndip(i)*un(k,i)*sigma2(k,i)
          fwindv(i,lat) = tm2(i)*un(k,i)*sigma1(k,i)+sndip(i)*vn(k,i)*
     |      sigma2(k,i) - tm3(i)*vn(k,i)*sigma1(k,i)
!
! Convert wind velocity from cm/s to m/s by multipling 1.e-2 and
! Convert B from Gauss to T multipling 1.e-4 and 1.e-2 from cm to m
!
          fwindu(i,lat) = fwindu(i,lat)*bmod(i,lat)*1.e-8/sndip(i)
          fwindv(i,lat) = fwindv(i,lat)*bmod(i,lat)*1.e-8/sndip(i)
          famieu(i,lat) = -tm1(i)*sigma2(k,i)*
     |       ui(k,i) + sndip(i)*sigma1(k,i)*vi(k,i)
     |       - tm3(i)*sigma2(k,i)*vi(k,i)
          famiev(i,lat) = -tm2(i)*sigma2(k,i)*
     |       vi(k,i) - sndip(i)*sigma1(k,i)*ui(k,i)
     |       - tm3(i)*sigma2(k,i)*ui(k,i)
!
! Convert wind velocity from cm/s to m/s by multipling 1.e-2 and
! Convert B from Gauss to T multipling 1.e-4 and 1.e-2 from cm to m
!
          famieu(i,lat) = famieu(i,lat)*bmod(i,lat)*1.e-8/sndip(i)**2
          famiev(i,lat) = famiev(i,lat)*bmod(i,lat)*1.e-8/sndip(i)**2
          fwindu(i,lat) = fwindu(i,lat)+(z(k+1,i)-z(k,i))*fwindu(i,lat)
          fwindv(i,lat) = fwindv(i,lat)+(z(k+1,i)-z(k,i))*fwindv(i,lat)
          famieu(i,lat) = famieu(i,lat)+(z(k+1,i)-z(k,i))*famieu(i,lat)
          famiev(i,lat) = famiev(i,lat)+(z(k+1,i)-z(k,i))*famiev(i,lat)
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1
      end subroutine dyndiag
!-----------------------------------------------------------------------
      subroutine dyndiag_sech(lon0,lon1,lat0,lat1)
!
! Save 2d (lon,lat) diagnostics to secondary histories:
! This is called from dynamics after lamdas lat loop.
!
! Note total fwind and famie are not calculated here (as they were in
!   tgcm15) because of the need for lat-1,lat+1,lon-1,lon+1. These can 
!   be calculated by a post-processor after reading the u and v components 
!   from the secondary histories.
!
      integer,intent(in) :: lon0,lon1,lat0,lat1
!
      call addfsech_ij('PEDERSEN',' ',' ',sigp(lon0:lon1,lat0:lat1),
     |  lon0,lon1,lat0,lat1)
      call addfsech_ij('HALL'    ,' ',' ',sigh(lon0:lon1,lat0:lat1),
     |  lon0,lon1,lat0,lat1)
      call addfsech_ij('QWIND'   ,' ',' ',
     |  qwind(lon0:lon1,lat0:lat1),lon0,lon1,lat0,lat1)
      call addfsech_ij('QAMIE'   ,' ',' ',
     |  qamie(lon0:lon1,lat0:lat1),lon0,lon1,lat0,lat1)
      call addfsech_ij('FWINDU'   ,' ',' ',
     |  fwindu(lon0:lon1,lat0:lat1),lon0,lon1,lat0,lat1)
      call addfsech_ij('FWINDV'   ,' ',' ',
     |  fwindv(lon0:lon1,lat0:lat1),lon0,lon1,lat0,lat1)
      call addfsech_ij('FAMIEU'   ,' ',' ',
     |  famieu(lon0:lon1,lat0:lat1),lon0,lon1,lat0,lat1)
      call addfsech_ij('FAMIEV'   ,' ',' ',
     |  famiev(lon0:lon1,lat0:lat1),lon0,lon1,lat0,lat1)
      end subroutine dyndiag_sech
!-----------------------------------------------------------------------
      end module dyndiag_module