module strong_fast_global_mod
!$$$   module documentation block
!                .      .    .                                       .
! module:    strong_fast_global_mod
!
! abstract: contains all routines for fast strong balance constraint
!

! program history log:
!   2008-04-04  safford - add moddule doc block and missing subroutine doc blocks
!   2012-02-08  kleist  - add uvflag in place of uv_hyb_ens
!   2012-11-23  parrish - Replace calls to spanaly_ns and spsynth_ns with inline code.
!                          Remove subroutines spanaly_ns and spsynth_ns.
!   2013-07-02  parrish - change name of init_strongvars_2 to init_strongvars.
!
! subroutines included:
!    init_strongvars          --
!    strong_bal_correction_fast_global
!    strong_bal_correction_fast_global_ad -- adjoint of strong_bal_correction
!    gather_rmstends0         -- get bal diagnostics
!    gather_rmstends          -- get bal diagnostics
!    inmi_coupler_sd2ew0      --
!    inmi_coupler_sd2ew1      --
!    inmi_coupler_sd2ew       --
!    inmi_coupler_ew2sd1      --
!    inmi_coupler_ew2sd       --
!    inmi_ew_trans            --
!    inmi_ew_invtrans_ad      --
!    inmi_ew_invtrans         --
!    inmi_ew_trans_ad         --
!    inmi_coupler_ew2ns0      --
!    inmi_coupler_ew2ns1      --
!    inmi_coupler_ew2ns       --
!    inmi_coupler_ns2ew       --
!    inmi_nsuvm2zdm           --
!    spdz2uv_ns               -- compute winds from div and vort for 1 zonal wave number
!    spuv2dz_ns               -- compute div,vort from winds for one zonal wave number
!
! variable definitions:
!
! attributes:
!   language:  f90
!   machine:   ibm rs/6000 sp
!$$$ end documentation block

  use kinds, only: i_kind,r_kind
  use gridmod, only: nlat,nlon,lat2,lon2,nsig,sp_a,jstart,istart
  use gridmod, only: ilat1,jlon1
  use constants, only: zero,one,two,rearth
  use mpimod, only: ierror,mpi_comm_world,mpi_integer4,mpi_rtype,mpi_sum,npe
  use mod_vtrans, only: speeds,nvmodes_keep,vtrans,vtrans_inv,vtrans_ad,vtrans_inv_ad
  implicit none

! set default to private
  private
! set subroutines to public
  public :: init_strongvars
  public :: strong_bal_correction_fast_global
  public :: strong_bal_correction_fast_global_ad
  public :: gather_rmstends0
  public :: gather_rmstends
  public :: inmi_coupler_sd2ew0
  public :: inmi_coupler_sd2ew1
  public :: inmi_coupler_sd2ew
  public :: inmi_coupler_ew2sd1
  public :: inmi_coupler_ew2sd
  public :: inmi_ew_trans
  public :: inmi_ew_invtrans_ad
  public :: inmi_ew_invtrans
  public :: inmi_ew_trans_ad
  public :: inmi_coupler_ew2ns0
  public :: inmi_coupler_ew2ns1
  public :: inmi_coupler_ew2ns
  public :: inmi_coupler_ns2ew
  public :: inmi_nsuvm2zdm
  public :: spdz2uv_ns
  public :: spuv2dz_ns

  integer(i_kind),allocatable:: mode_list(:,:)
                                                   !  mode_list(1,j) = lat index for ew strip j
                                                   !  mode_list(2,j) = vert mode number for ew strip j
                                                   !  mode_list(3,j) = pe of this lat/vert mode strip 
  integer(i_kind),allocatable:: mmode_list(:,:)
                                                   !  mmode_list(1,j) = m1 (zonal wave number 1) for ns strip
                                                   !  mmode_list(2,j) = m2 (zonal wave number 2) for ns strip
                                                   !  mmode_list(3,j) = vert mode number 1 for ns strip j
                                                   !  mmode_list(4,j) = vert mode number 2 for ns strip j
                                                   !  mmode_list(5,j) = pe for ns strip j
  integer(i_kind) nlatm_0,nlatm_1,m_0,m_1
  integer(i_kind) mthis
  integer(i_kind) nallsend_ew2sd,nallrecv_ew2sd
  integer(i_kind) nallsend,nallrecv
  integer(i_kind) nallsend_sd2ew,nallrecv_sd2ew
  integer(i_kind),allocatable,dimension(:)::mthis0,ndisp,indexglob
  integer(i_kind),allocatable,dimension(:)::nsend_sd2ew,nrecv_sd2ew
  integer(i_kind),allocatable,dimension(:)::ndsend_sd2ew,ndrecv_sd2ew
  integer(i_kind),allocatable,dimension(:)::nsend_ew2sd,nrecv_ew2sd
  integer(i_kind),allocatable,dimension(:)::ndsend_ew2sd,ndrecv_ew2sd
  integer(i_kind),allocatable,dimension(:)::nsend,nrecv,ndsend,ndrecv
  integer(i_kind),allocatable,dimension(:,:)::info_send_sd2ew,info_recv_sd2ew
  integer(i_kind),allocatable,dimension(:,:)::info_send_ew2sd,info_recv_ew2sd
  integer(i_kind),allocatable,dimension(:,:)::info_send,info_recv

contains

subroutine init_strongvars(mype)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    init_strongvars
!
!   prgrmmr: parrish
!
! abstract: initialize strong fast global initialization
!
! program history log:
!   2008-04-04  safford -- add subprogram doc block
!   2010-03-31  treadon - replace specmod jcap with sp_a structure
!
!   input argument list:
!     mype     - mpi task id
!
!   output argument list:
!
! attributes:
!   language:  f90
!   machine:   ibm rs/6000 sp
!
!$$$

  implicit none

  integer(i_kind),intent(in   ) :: mype

  allocate(mode_list(3,nlat*nvmodes_keep),mmode_list(5,(sp_a%jcap+1)*nvmodes_keep))
  call inmi_coupler_sd2ew0(mype)
  call inmi_coupler_ew2ns0(mype)
  call gather_rmstends0
  call inmi_coupler_sd2ew1(mype)
  call inmi_coupler_ew2ns1(mype)
  call inmi_coupler_ew2sd1(mype)

  return
end subroutine init_strongvars


subroutine strong_bal_correction_fast_global(u_t,v_t,t_t,ps_t,mype,psi,chi,t,ps, &
                                              bal_diagnostic,fullfield,update,uvflag)

!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    strong_bal_correction  strong balance correction
!   prgmmr: parrish          org: np23                date: 2006-07-15
!
! abstract: given input perturbation tendencies of u,v,t,ps from tlm on gaussian grid,
!           and input perturbation u,v,t,ps, compute balance adjustment to u,v,t,ps
!           which zeroes out input gravity component of perturbation tendencies.
!           also output, for later use, input tendencies projected onto gravity modes.
!
! program history log:
!   2006-07-15  parrish
!   2007-04-16  kleist  - modified for full field or incremental diagnostics
!   2008-10-08  parrish/derber - modify to output streamfunction and vel. pot. and not update time derivatives
!   2009-11-27  parrish - add uvflag.  if uvflag=true, then
!                          input/output variables psi=u, chi=v.
!   2010-03-31  treadon - replace specmod jcap with sp_a structure
!   2012-02-08  kleist  - add uvflag in place of uv_hyb_ens
!   2013-10-26  todling - prevent division by zero when rms's are zero
!
!   input argument list:
!     u_t      - input perturbation u tendency on gaussian grid (subdomains)
!     v_t      - input perturbation v tendency on gaussian grid (subdomains)
!     t_t      - input perturbation t tendency on gaussian grid (subdomains)
!     ps_t     - input perturbation surface pressure tendency on gaussian grid (subdomains)
!     mype     - current processor
!     psi      - input perturbation psi on gaussian grid (subdomains)
!     chi      - input perturbation chi on gaussian grid (subdomains)
!     t        - input perturbation t on gaussian grid (subdomains)
!     ps       - input perturbation surface pressure on gaussian grid (subdomains)
!     bal_diagnostic - if true, then compute bal diagnostic, a measure of amplitude
!                      of balanced gravity mode tendencies
!     fullfield - if true, full field diagnostics
!                 if false, incremental 
!     update   - if false, then do not update u,v,t,ps with balance increment
!
!   output argument list:
!     psi      - output balanced perturbation u on gaussian grid (subdomains)
!     chi      - output balanced perturbation v on gaussian grid (subdomains)
!     t        - output balanced perturbation t on gaussian grid (subdomains)
!     ps       - output balanced perturbation surface pressure on gaussian grid (subdomains)
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!$$$

  use mod_strong, only: dinmi,gproj,gproj_diag,gproj_diag_update
  use constants, only: tiny_r_kind
  implicit none

  integer(i_kind)                       ,intent(in   ) :: mype
  logical                               ,intent(in   ) :: bal_diagnostic,update,fullfield,uvflag
  real(r_kind),dimension(lat2,lon2,nsig),intent(in   ) :: u_t,v_t,t_t
  real(r_kind),dimension(lat2,lon2)     ,intent(in   ) :: ps_t
  real(r_kind),dimension(lat2,lon2,nsig),intent(inout) :: psi,chi,t
  real(r_kind),dimension(lat2,lon2)     ,intent(inout) :: ps

  real(r_kind),dimension(nvmodes_keep)::rmstend_uf,rmstend_g_uf
  real(r_kind),dimension(nvmodes_keep)::rmstend_f,rmstend_g_f

  real(r_kind),dimension(lat2,lon2,nvmodes_keep)::utilde_t,vtilde_t,mtilde_t
  real(r_kind),dimension(lat2,lon2,nvmodes_keep)::delutilde,delvtilde,delmtilde
  real(r_kind),dimension(2,0:sp_a%jcap)::divhat,vorthat,mhat,deldivhat,delvorthat,delmhat
  real(r_kind),allocatable,dimension(:,:)::rmstend_loc_uf,rmstend_g_loc_uf
  real(r_kind),allocatable,dimension(:,:)::rmstend_loc_f,rmstend_g_loc_f
  real(r_kind),allocatable,dimension(:,:,:,:)::uvm_ew
  real(r_kind),allocatable,dimension(:,:,:,:,:)::uvm_ewtrans,uvm_ns,zdm_hat
  real(r_kind) rmstend_all_uf,rmstend_all_g_uf,rmstend_all_f,rmstend_all_g_f
  real(r_kind) del2inv,rn,gspeed
  real(r_kind) diff1,diffi

  integer(i_kind) i,ipair,kk,mode,n,mmax,m

  if(.not. (update .or. bal_diagnostic))return

  mmax=sp_a%jcap

!   1.  u,v,t,ps   -->    utilde,vtilde,mtilde  (vertical mode transform)
!       (subdomains)         (subdomains)

  call vtrans(u_t,v_t,t_t,ps_t,utilde_t,vtilde_t,mtilde_t)

  allocate(uvm_ew(nlon,2,3,nlatm_0:nlatm_1),uvm_ewtrans(2,0:sp_a%jcap,2,3,nlatm_0:nlatm_1))
  allocate(uvm_ns(3,2,nlat,2,m_0:m_1),zdm_hat(3,2,nlat,2,m_0:m_1))
  call inmi_coupler_sd2ew(utilde_t,vtilde_t,mtilde_t,utilde_t,vtilde_t,mtilde_t, &
                          uvm_ew,mype)
  call inmi_ew_trans(uvm_ew,uvm_ewtrans)
  call inmi_coupler_ew2ns(uvm_ewtrans,uvm_ns)
  call inmi_nsuvm2zdm(uvm_ns,zdm_hat)
  if(bal_diagnostic)then
    allocate(rmstend_loc_f(2,m_0:m_1))
    allocate(rmstend_g_loc_f(2,m_0:m_1))
    rmstend_loc_f=zero
    rmstend_g_loc_f=zero
    allocate(rmstend_loc_uf(2,m_0:m_1))
    allocate(rmstend_g_loc_uf(2,m_0:m_1))
    rmstend_g_loc_uf=zero
    rmstend_loc_uf=zero
  end if
  if(update)then
!$omp parallel do  schedule(dynamic,1) private(kk,ipair,m,mode,gspeed,i,n) &
!$omp private(vorthat,divhat,delvorthat,deldivhat,delmhat,mhat,rn,del2inv)
    do kk=m_0,m_1
       do ipair=1,2
          m=mmode_list(ipair,kk)
          mode=mmode_list(ipair+2,kk)
          gspeed=speeds(abs(mode))
          i=0
          do n=m,sp_a%jcap
             i=i+1
             vorthat(1,n)=zdm_hat(1,1,i,ipair,kk)
             vorthat(2,n)=zdm_hat(1,2,i,ipair,kk)
             divhat( 1,n)=zdm_hat(2,1,i,ipair,kk)
             divhat( 2,n)=zdm_hat(2,2,i,ipair,kk)
             mhat(   1,n)=zdm_hat(3,1,i,ipair,kk)
             mhat(   2,n)=zdm_hat(3,2,i,ipair,kk)
          end do

!   4.  divhat,vorthat,mhat --> deldivhat,delvorthat,delmhat   (inmi correction)
!          (slabs)                        (slabs)

          if(mode >  0) then
!               here, delvorthat, etc contain field corrections necessary to zero out gravity component
!                                         of tendencies
             call dinmi(vorthat(1,m),divhat(1,m),mhat(1,m),delvorthat(1,m),deldivhat(1,m),delmhat(1,m),&
                      m,mmax,gspeed)
          else
!               here, delvorthat, etc contain gravity component of tendencies
             if(bal_diagnostic) then   ! bal_diagnostic and update
                call gproj_diag_update(vorthat(1,m),divhat(1,m),mhat(1,m),delvorthat(1,m),deldivhat(1,m),delmhat(1,m), &
                          rmstend_loc_uf(ipair,kk),rmstend_g_loc_uf(ipair,kk),rmstend_loc_f(ipair,kk), &
                          rmstend_g_loc_f(ipair,kk),m,mmax,gspeed)
             else              !  update only
                call gproj(vorthat(1,m),divhat(1,m),mhat(1,m),delvorthat(1,m),deldivhat(1,m),delmhat(1,m), &
                          m,mmax,gspeed)
             end if
          end if

          if(.not. uvflag)then
             do n=m,sp_a%jcap
                if(n >  0) then
                   rn=real(n,r_kind)
                   del2inv=-rearth**2/(rn*(rn+one))
                else
                   del2inv=zero
                end if
                delvorthat(1,n)=delvorthat(1,n)*del2inv
                delvorthat(2,n)=delvorthat(2,n)*del2inv
                deldivhat(1,n)=deldivhat(1,n)*del2inv
                deldivhat(2,n)=deldivhat(2,n)*del2inv
             end do
          end if
             
          i=0
          do n=m,sp_a%jcap
             i=i+1
             zdm_hat(1,1,i,ipair,kk)=delvorthat(1,n)
             zdm_hat(1,2,i,ipair,kk)=delvorthat(2,n)
             zdm_hat(2,1,i,ipair,kk)=deldivhat(1,n)
             zdm_hat(2,2,i,ipair,kk)=deldivhat(2,n)
             zdm_hat(3,1,i,ipair,kk)=delmhat(1,n)
             zdm_hat(3,2,i,ipair,kk)=delmhat(2,n)
          end do
        end do
     end do
!   7.  delutilde,delvtilde,delmtilde  -->  psi,chi,t,ps   (vertical mode inverse transform)
!       (subdomains)                      (subdomains)

!  update u,v,t,ps

     if(uvflag) then
        call inmi_nszdm2uvm(uvm_ns,zdm_hat)
     else
        call inmi_nspcm_hat2pcm(uvm_ns,zdm_hat)
     end if
     call inmi_coupler_ns2ew(uvm_ewtrans,uvm_ns)
     call inmi_ew_invtrans(uvm_ew,uvm_ewtrans)
     call inmi_coupler_ew2sd(delutilde,delvtilde,delmtilde,utilde_t,vtilde_t,mtilde_t,uvm_ew,mype)
     call vtrans_inv(delutilde,delvtilde,delmtilde,psi,chi,t,ps)
!    u_t_g=zero;v_t_g=zero;t_t_g=zero;ps_t_g=zero
!    call vtrans_inv(utilde_t,vtilde_t,mtilde_t,u_t_g,v_t_g,t_t_g,ps_t_g)
  else            ! bal_diagnostic only no update
!$omp parallel do  schedule(dynamic,1) private(kk,ipair,m,mode,gspeed,i,n) &
!$omp private(vorthat,divhat,mhat)
     do kk=m_0,m_1
        do ipair=1,2
           mode=mmode_list(ipair+2,kk)
           if(mode > 0)cycle
           m=mmode_list(ipair,kk)
           gspeed=speeds(abs(mode))
           i=0
           do n=m,sp_a%jcap
              i=i+1
              vorthat(1,n)=zdm_hat(1,1,i,ipair,kk)
              vorthat(2,n)=zdm_hat(1,2,i,ipair,kk)
              divhat( 1,n)=zdm_hat(2,1,i,ipair,kk)
              divhat( 2,n)=zdm_hat(2,2,i,ipair,kk)
              mhat(   1,n)=zdm_hat(3,1,i,ipair,kk)
              mhat(   2,n)=zdm_hat(3,2,i,ipair,kk)
           end do

           call gproj_diag(vorthat(1,m),divhat(1,m),mhat(1,m), &
                     rmstend_loc_uf(ipair,kk),rmstend_g_loc_uf(ipair,kk),rmstend_loc_f(ipair,kk), &
                     rmstend_g_loc_f(ipair,kk),m,mmax,gspeed)
        end do
     end do
  end if

  deallocate(uvm_ew,uvm_ewtrans,uvm_ns,zdm_hat)


  if(bal_diagnostic) then
     call gather_rmstends(rmstend_loc_uf,  rmstend_uf)
     call gather_rmstends(rmstend_g_loc_uf,rmstend_g_uf)
     call gather_rmstends(rmstend_loc_f,   rmstend_f)
     call gather_rmstends(rmstend_g_loc_f, rmstend_g_f)
     if(mype == 0) then
        rmstend_all_uf=zero
        rmstend_all_g_uf=zero

        if (fullfield) then
           write(6,*) 'strong_fast_global:   full field balance diagnostics --  '
        else
           write(6,*) 'strong_fast_global:   incremental balance diagnostics --  '
        end if

        do i=1,nvmodes_keep
           rmstend_all_uf=rmstend_all_uf+rmstend_uf(i)
           rmstend_all_g_uf=rmstend_all_g_uf+rmstend_g_uf(i)
           diffi = rmstend_uf(i)-rmstend_g_uf(i) 
           diff1 = rmstend_uf(1)-rmstend_g_uf(1)
           if(abs(diffi)<tiny_r_kind.or.abs(diff1)<tiny_r_kind) then
             write(6,'(" mode,rmstend_uf,rmstend_g_uf,rat = ",i5,2e28.18)') &
                              i,rmstend_uf(i),rmstend_g_uf(i)
           else
             write(6,'(" mode,rmstend_uf,rmstend_g_uf,rat = ",i5,2e28.18,2f24.18)') &
                              i,rmstend_uf(i),rmstend_g_uf(i),&
                              rmstend_g_uf(i)/diffi, &
                              rmstend_g_uf(i)/diff1
           endif
        end do
        rmstend_all_f=zero
        rmstend_all_g_f=zero
        do i=1,nvmodes_keep
           rmstend_all_f=rmstend_all_f+rmstend_f(i)
           rmstend_all_g_f=rmstend_all_g_f+rmstend_g_f(i)
           diffi = rmstend_f(i)-rmstend_g_f(i)
           diff1 = rmstend_f(1)-rmstend_g_f(1)
           if(abs(diffi)<tiny_r_kind.or.abs(diff1)<tiny_r_kind) then
             write(6,'(" mode,rmstend_f,rmstend_g_f = ",i5,2e28.18)') &
                              i,rmstend_f(i),rmstend_g_f(i)
           else
             write(6,'(" mode,rmstend_f,rmstend_g_f,rat = ",i5,2e28.18,2f24.18)') &
                              i,rmstend_f(i),rmstend_g_f(i),&
                              rmstend_g_f(i)/diffi, &
                              rmstend_g_f(i)/diff1
           endif
        end do
        diffi = rmstend_all_uf-rmstend_all_g_uf
        diff1 = rmstend_all_f-rmstend_all_g_f
        if(abs(diffi)<tiny_r_kind.or.abs(diff1)<tiny_r_kind) then
          write(6,'(" rmstend_all_uf,g_uf,rat = ",2e28.18)') rmstend_all_uf,rmstend_all_g_uf
          write(6,'(" rmstend_all_f,g_f,rat = ",2e28.18)') rmstend_all_f,rmstend_all_g_f
        else
          write(6,'(" rmstend_all_uf,g_uf,rat = ",2e28.18,f24.18)') rmstend_all_uf,rmstend_all_g_uf, &
                                                   diffi
          write(6,'(" rmstend_all_f,g_f,rat = ",2e28.18,f24.18)') rmstend_all_f,rmstend_all_g_f, &
                                                   diff1
        endif
     end if
     deallocate(rmstend_loc_uf,rmstend_g_loc_uf)
     deallocate(rmstend_loc_f,rmstend_g_loc_f)
  end if
  return

end subroutine strong_bal_correction_fast_global

subroutine strong_bal_correction_fast_global_ad(u_t,v_t,t_t,ps_t,mype,psi,chi,t,ps,uvflag)

!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    strong_bal_correction_ad adjoint of strong_bal_correction
!   prgmmr: parrish          org: np23                date: 2006-07-15
!
! abstract: adjoint of strong_bal_correction
!
! program history log:
!   2006-07-15  parrish
!   2008-10-08  parrish/derber - modify to output streamfunction and vel. pot. and not update time derivatives
!   2009-11-27  parrish - add uvflag.  if present and true, then
!                          input/output variables psi=u, chi=v.
!   2010-03-31  treadon - replace specmod jcap with sp_a structure
!   2012-02-08  kleist  - add uvflag in place of uv_hyb_ens
!
!   input argument list:
!     u_t      - input perturbation u tendency on gaussian grid (subdomains)
!     v_t      - input perturbation v tendency on gaussian grid (subdomains)
!     t_t      - input perturbation t tendency on gaussian grid (subdomains)
!     ps_t     - input perturbation surface pressure tendency on gaussian grid (subdomains)
!     mype     - current processor
!     psi      - input perturbation psi on gaussian grid (subdomains)
!     chi      - input perturbation chi on gaussian grid (subdomains)
!     t        - input perturbation t on gaussian grid (subdomains)
!     ps       - input perturbation surface pressure on gaussian grid (subdomains)
!
!   output argument list:
!     u_t      - output perturbation u tendency on gaussian grid (subdomains)
!     v_t      - output perturbation v tendency on gaussian grid (subdomains)
!     t_t      - output perturbation t tendency on gaussian grid (subdomains)
!     ps_t     - output perturbation surface pressure tendency on gaussian grid (subdomains)
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!$$$

  use mod_strong, only: dinmi_ad,gproj_ad
  implicit none

  integer(i_kind)                       ,intent(in   ) :: mype
  real(r_kind),dimension(lat2,lon2,nsig),intent(inout) :: u_t,v_t,t_t
  real(r_kind),dimension(lat2,lon2)     ,intent(inout) :: ps_t
  real(r_kind),dimension(lat2,lon2,nsig),intent(in   ) :: psi,chi,t
  real(r_kind),dimension(lat2,lon2)     ,intent(in   ) :: ps
  logical, intent(in)                                  :: uvflag

  real(r_kind),dimension(lat2,lon2,nvmodes_keep)::utilde_t,vtilde_t,mtilde_t
  real(r_kind),dimension(lat2,lon2,nvmodes_keep)::delutilde,delvtilde,delmtilde
  real(r_kind),dimension(2,0:sp_a%jcap)::divhat,vorthat,mhat,deldivhat,delvorthat,delmhat
  real(r_kind),allocatable,dimension(:,:,:,:)::uvm_ew
  real(r_kind),allocatable,dimension(:,:,:,:,:)::uvm_ewtrans,uvm_ns,zdm_hat
  real(r_kind) del2inv,rn,gspeed

  integer(i_kind) i,ipair,kk,mode,n,m,mmax

  mmax=sp_a%jcap

!  adjoint of update u,v,t,ps

!   7.  adjoint of delutilde,delvtilde,delmtilde  -->  psi,chi,t,ps  (vert mode inverse transform)
!       (subdomains)                      (subdomains)

  delutilde=zero ; delvtilde=zero ; delmtilde=zero
  call vtrans_inv_ad(delutilde,delvtilde,delmtilde,psi,chi,t,ps)

  allocate(uvm_ew(nlon,2,3,nlatm_0:nlatm_1),uvm_ewtrans(2,0:sp_a%jcap,2,3,nlatm_0:nlatm_1))
  allocate(uvm_ns(3,2,nlat,2,m_0:m_1),zdm_hat(3,2,nlat,2,m_0:m_1))
  utilde_t=zero ; vtilde_t=zero ; mtilde_t=zero
  call inmi_coupler_sd2ew(delutilde,delvtilde,delmtilde,utilde_t,vtilde_t,mtilde_t,uvm_ew,mype)
  call inmi_ew_invtrans_ad(uvm_ew,uvm_ewtrans)
  call inmi_coupler_ew2ns(uvm_ewtrans,uvm_ns)
  if(uvflag) then
     call inmi_nszdm2uvm_ad(uvm_ns,zdm_hat)
  else
     call inmi_nspcm_hat2pcm_ad(uvm_ns,zdm_hat)
  end if
!$omp parallel do  schedule(dynamic,1) private(kk,ipair,m,mode,gspeed,i,n) &
!$omp private(vorthat,divhat,delvorthat,deldivhat,delmhat,mhat,rn,del2inv)
  do kk=m_0,m_1
     do ipair=1,2
        m=mmode_list(ipair,kk)
        mode=mmode_list(ipair+2,kk)
        gspeed=speeds(abs(mode))
        do n=m,sp_a%jcap
           vorthat(1,n)=zero
           vorthat(2,n)=zero
           divhat( 1,n)=zero
           divhat( 2,n)=zero
           mhat(   1,n)=zero
           mhat(   2,n)=zero
        end do
        i=0
        do n=m,sp_a%jcap
           i=i+1
           delvorthat(1,n)=zdm_hat(1,1,i,ipair,kk)
           delvorthat(2,n)=zdm_hat(1,2,i,ipair,kk)
           deldivhat(1,n)=zdm_hat(2,1,i,ipair,kk)
           deldivhat(2,n)=zdm_hat(2,2,i,ipair,kk)
           delmhat(1,n)=zdm_hat(3,1,i,ipair,kk)
           delmhat(2,n)=zdm_hat(3,2,i,ipair,kk)
        end do
        if(.not. uvflag) then
           do n=m,sp_a%jcap
              if(n >  0) then
                 rn=real(n,r_kind) 
                 del2inv=-rearth**2/(rn*(rn+one))
              else
                 del2inv=zero
              end if
              delvorthat(1,n)=delvorthat(1,n)*del2inv
              delvorthat(2,n)=delvorthat(2,n)*del2inv
              deldivhat(1,n)=deldivhat(1,n)*del2inv
              deldivhat(2,n)=deldivhat(2,n)*del2inv
           end do
        end if
        if(mode >  0) then
           call dinmi_ad(vorthat(1,m),divhat(1,m),mhat(1,m),&
                       delvorthat(1,m)   ,   deldivhat(1,m),   delmhat(1,m),m,mmax,gspeed)
        else
           call gproj_ad(vorthat(1,m),divhat(1,m),mhat(1,m),delvorthat(1,m),deldivhat(1,m), &
                 delmhat(1,m),m,mmax,gspeed)
        end if

        i=0
        do n=m,sp_a%jcap
           i=i+1
           zdm_hat(1,1,i,ipair,kk)=vorthat(1,n)
           zdm_hat(1,2,i,ipair,kk)=vorthat(2,n)
           zdm_hat(2,1,i,ipair,kk)=divhat(1,n)
           zdm_hat(2,2,i,ipair,kk)=divhat(2,n)
           zdm_hat(3,1,i,ipair,kk)=mhat(1,n)
           zdm_hat(3,2,i,ipair,kk)=mhat(2,n)
        end do

     end do
  end do

  call inmi_nsuvm2zdm_ad(uvm_ns,zdm_hat)
  call inmi_coupler_ns2ew(uvm_ewtrans,uvm_ns)
  call inmi_ew_trans_ad(uvm_ew,uvm_ewtrans)
  call inmi_coupler_ew2sd(utilde_t,vtilde_t,mtilde_t,delutilde,delvtilde,delmtilde,uvm_ew,mype)

  utilde_t=utilde_t+delutilde
  vtilde_t=vtilde_t+delvtilde
  mtilde_t=mtilde_t+delmtilde
!
!    1.  adjoint of u,v,t,ps   -->    utilde,vtilde,mtilde  (vertical mode transform)
!                      (subdomains)         (subdomains)

  call vtrans_ad(u_t,v_t,t_t,ps_t,utilde_t,vtilde_t,mtilde_t)
  deallocate(uvm_ew,uvm_ewtrans,uvm_ns,zdm_hat)

end subroutine strong_bal_correction_fast_global_ad

subroutine gather_rmstends0

!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    gather_rmstends0  get bal diagnostics
!   prgmmr: parrish          org: np23                date: 2006-08-03
!
! abstract: compute bal diagnostics which give amplitude of 
!           gravity projection of energy norm of tendencies
!
! program history log:
!   2006-08-03  parrish
!   2008-04-04  safford  - rm unused uses
!   2010-03-31  treadon - replace specmod jcap with sp_a structure
!
!   input argument list:
!
!   output argument list:
!     rmstend  -  all vertical modes of rmstend assembled across all processors
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!$$$

  implicit none

  integer(i_kind),dimension(m_0:m_1):: indexloc
  integer(i_kind) i
  
  allocate(mthis0(npe),ndisp(npe+1),indexglob((sp_a%jcap+1)*nvmodes_keep))
  do i=m_0,m_1
     indexloc(i)=i
  end do
  mthis=m_1-m_0+1
  call mpi_allgather(mthis,1,mpi_integer4,mthis0,1,mpi_integer4,mpi_comm_world,ierror)
  ndisp(1)=0
  do i=2,npe+1
     ndisp(i)=ndisp(i-1)+mthis0(i-1)
  end do
  call mpi_allgatherv(indexloc,mthis,mpi_integer4, &
                      indexglob,mthis0,ndisp,mpi_integer4,mpi_comm_world,ierror)

end subroutine gather_rmstends0

subroutine gather_rmstends(rmstend_loc,rmstend)

!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    gather_rmstends  get bal diagnostics
!   prgmmr: parrish          org: np23                date: 2006-08-03
!
! abstract: compute bal diagnostics which give amplitude of 
!           gravity projection of energy norm of tendencies
!
! program history log:
!   2006-08-03  parrish
!   2008-04-04  safford  - rm unused uses
!
!   input argument list:
!     rmstend_loc - previously computed energy norms of vertical modes
!                   as distributed on local processors
!
!   output argument list:
!     rmstend  -  all vertical modes of rmstend assembled across all processors
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!$$$

  implicit none

  real(r_kind),intent(in   ) :: rmstend_loc(2,m_0:m_1)
  real(r_kind),intent(  out) :: rmstend(nvmodes_keep)

  real(r_kind),dimension(2,(sp_a%jcap+1)*nvmodes_keep)::work
  integer(i_kind) i,ii,mode,mpi_string1
  
  call mpi_type_contiguous(2,mpi_rtype,mpi_string1,ierror)
  call mpi_type_commit(mpi_string1,ierror)
  call mpi_allgatherv(rmstend_loc,mthis,mpi_string1, &
                     work,mthis0,ndisp,mpi_string1,mpi_comm_world,ierror)
  call mpi_type_free(mpi_string1,ierror)
  rmstend=zero
  do i=1,ndisp(npe+1)
     ii=indexglob(i)
     mode=mmode_list(3,ii)
     if(mode <  0) rmstend(-mode)=work(1,ii)+rmstend(-mode)
     mode=mmode_list(4,ii)
     if(mode <  0) rmstend(-mode)=work(2,ii)+rmstend(-mode)
  end do

end subroutine gather_rmstends

subroutine inmi_coupler_sd2ew0(mype)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    inmi_coupler_sd2ew0
!   prgmmr:
!
! abstract:  create ew (lat strips) subdivision for use with inmi
!
! program history log:
!   2008-04-04  safford  - add doc block
!
!   input argument list:
!     mype        - current processor number
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!$$$

  implicit none

  integer(i_kind),intent(in   ) :: mype

  integer(i_kind) i,k,kchk,kk,n,nn,nlatm_this

  nlatm_this=nlat*nvmodes_keep/npe
  if(mod(nlat*nvmodes_keep,npe)==0) then
     kchk=npe
  else
     nlatm_this=nlatm_this+1
     kchk=mod(nlat*nvmodes_keep,npe)
  end if

  nn=0
  do k=1,nvmodes_keep
     do i=1,nlat
        nn=nn+1
        mode_list(1,nn)=i
        mode_list(2,nn)=k
        mode_list(3,nn)=-1
     end do
  end do
  
  nlatm_0=-1
  nlatm_1=-2
  nn=0
  do n=1,npe
     if(n <= kchk) then
        kk=nlatm_this
     else
        kk=nlatm_this-1
     end if
     if(kk >  0) then
        if(mype+1 == n) then
           nlatm_0=nn+1
           nlatm_1=nn+kk
        end if
        do k=1,kk
           nn=nn+1
           mode_list(3,nn)=n
        end do
     end if
  end do
end subroutine inmi_coupler_sd2ew0

subroutine inmi_coupler_sd2ew1(mype)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    inmi_coupler_sd2ew1
!
!   prgrmmr:
!
! abstract: use mpi_alltoallv to move u_sd,v_sd,m_sd (subdomains) to uvm_ew (lat strips)
!
! program history log:
!   2008-04-04  safford -- add subprogram doc block
!
!   input argument list:
!     mype     - mpi task id
!
!   output argument list:
!
! attributes:
!   language:  f90
!   machine:   ibm rs/6000 sp
!
!$$$

  implicit none


  integer(i_kind),intent(in   ) :: mype

  integer(i_kind) i,ii,ii0,ilat,imode,j,mm1,nlonloc,ipe,ilatm,ilon,mpi_string1
  real(r_kind),dimension(nlat,nvmodes_keep):: mode2_list

  allocate(nsend_sd2ew(npe),nrecv_sd2ew(npe))
  allocate(ndsend_sd2ew(npe+1),ndrecv_sd2ew(npe+1))
  mm1=mype+1

  mode2_list=0
  do i=1,nlat*nvmodes_keep
     ilat=mode_list(1,i)
     imode=mode_list(2,i)
     if(mode2_list(ilat,imode) /= 0) then
        if(mype == 0) write(6,*)' problem in inmi_coupler_sd2ew0'
        call mpi_finalize(i)
        stop
     end if
     mode2_list(ilat,imode)=i
  end do
  do imode=1,nvmodes_keep
     do ilat=1,nlat
        if(mode2_list(ilat,imode) == 0) then
           if(mype == 0) write(6,*)' problem in inmi_coupler_sd2ew0'
           call mpi_finalize(i)
           stop
        end if
     end do
  end do



!  obtain counts of points to send to each pe from this pe

  nsend_sd2ew=0
  nlonloc=lon2-2
  do imode=1,nvmodes_keep
     if(imode == 0) cycle
     do i=2,lat2-1
        ilat=i+istart(mm1)-2
        j=mode2_list(ilat,imode)
        ipe=mode_list(3,j)
        nsend_sd2ew(ipe)=nsend_sd2ew(ipe)+nlonloc
     end do
  end do

  ndsend_sd2ew(1)=0
  do i=2,npe+1
     ndsend_sd2ew(i)=ndsend_sd2ew(i-1)+nsend_sd2ew(i-1)
  end do
  nallsend_sd2ew=ndsend_sd2ew(npe+1)
  allocate(info_send_sd2ew(2,nallsend_sd2ew))
  nsend_sd2ew=0
  do imode=1,nvmodes_keep
     if(imode == 0) cycle
     do i=2,lat2-1
        ilat=i+istart(mm1)-2
        ilatm=mode2_list(ilat,imode)
        ipe=mode_list(3,ilatm)
        do ii=2,lon2-1
           ilon=ii+jstart(mm1)-2
           nsend_sd2ew(ipe)=nsend_sd2ew(ipe)+1
           ii0=ndsend_sd2ew(ipe)+nsend_sd2ew(ipe)
           info_send_sd2ew(1,ii0)=ilon
           info_send_sd2ew(2,ii0)=ilatm
        end do
     end do
  end do

  call mpi_alltoall(nsend_sd2ew,1,mpi_integer4,nrecv_sd2ew,1,mpi_integer4,mpi_comm_world,ierror)
  ndrecv_sd2ew(1)=0
  do i=2,npe+1
     ndrecv_sd2ew(i)=ndrecv_sd2ew(i-1)+nrecv_sd2ew(i-1)
  end do
  nallrecv_sd2ew=ndrecv_sd2ew(npe+1)
  allocate(info_recv_sd2ew(2,nallrecv_sd2ew))
  call mpi_type_contiguous(2,mpi_integer4,mpi_string1,ierror)
  call mpi_type_commit(mpi_string1,ierror)
  call mpi_alltoallv(info_send_sd2ew,nsend_sd2ew,ndsend_sd2ew,mpi_string1, &
                     info_recv_sd2ew,nrecv_sd2ew,ndrecv_sd2ew,mpi_string1,mpi_comm_world,ierror)
  call mpi_type_free(mpi_string1,ierror)
    
end subroutine inmi_coupler_sd2ew1

subroutine inmi_coupler_sd2ew(u_sd1,v_sd1,m_sd1,u_sd2,v_sd2,m_sd2,uvm_ew,mype)

!$$$  subprogram documentation block
!                .      .    .
! subprogram:    init_coupler_sd2ew
!
!   prgrmmr:
!
! abstract: use mpi_alltoallv to move u_sd,v_sd,m_sd (subdomains) to uvm_ew (lat strips)
!
! program history log:
!   2008-04-04  safford -- add subprogram doc block, rm unused vars
!
!   input argument list:
!     mype     - mpi task id
!     u_sd1    -
!     v_sd1    -
!     m_sd1    -
!     u_sd2    -
!     v_sd2    -
!     m_sd2    -
!
!   output argument list:
!     uvm_ew   -
!
! attributes:
!   language:  f90
!   machine:   ibm rs/6000 sp
!
!$$$

  implicit none


  integer(i_kind)                                 ,intent(in   ) :: mype
  real(r_kind),dimension(lat2,lon2,nvmodes_keep)  ,intent(in   ) :: u_sd1,v_sd1,m_sd1
  real(r_kind),dimension(lat2,lon2,nvmodes_keep)  ,intent(in   ) :: u_sd2,v_sd2,m_sd2
  real(r_kind),dimension(nlon,2,3,nlatm_0:nlatm_1),intent(  out) :: uvm_ew

  real(r_kind),allocatable,dimension(:,:,:)::sendbuf,recvbuf
  integer(i_kind) ilat,imode,j,mm1,ilatm,ilon,mpi_string1

  mm1=mype+1

  allocate(sendbuf(2,3,nallsend_sd2ew))
  do j=1,nallsend_sd2ew
     ilon=info_send_sd2ew(1,j)-jstart(mm1)+2
     ilatm=info_send_sd2ew(2,j)
     ilat=mode_list(1,ilatm)-istart(mm1)+2
     imode=mode_list(2,ilatm)
     sendbuf(1,1,j)=u_sd1(ilat,ilon,imode)
     sendbuf(1,2,j)=v_sd1(ilat,ilon,imode)
     sendbuf(1,3,j)=m_sd1(ilat,ilon,imode)
     sendbuf(2,1,j)=u_sd2(ilat,ilon,imode)
     sendbuf(2,2,j)=v_sd2(ilat,ilon,imode)
     sendbuf(2,3,j)=m_sd2(ilat,ilon,imode)
  end do
  allocate(recvbuf(2,3,nallrecv_sd2ew))
  call mpi_type_contiguous(6,mpi_rtype,mpi_string1,ierror)
  call mpi_type_commit(mpi_string1,ierror)
  call mpi_alltoallv(sendbuf,nsend_sd2ew,ndsend_sd2ew,mpi_string1, &
                     recvbuf,nrecv_sd2ew,ndrecv_sd2ew,mpi_string1,mpi_comm_world,ierror)
  call mpi_type_free(mpi_string1,ierror)
  deallocate(sendbuf)
    
  do j=1,nallrecv_sd2ew
     ilon=info_recv_sd2ew(1,j)
     ilatm=info_recv_sd2ew(2,j)
     uvm_ew(ilon,1,1,ilatm)=recvbuf(1,1,j)
     uvm_ew(ilon,1,2,ilatm)=recvbuf(1,2,j)
     uvm_ew(ilon,1,3,ilatm)=recvbuf(1,3,j)
     uvm_ew(ilon,2,1,ilatm)=recvbuf(2,1,j)
     uvm_ew(ilon,2,2,ilatm)=recvbuf(2,2,j)
     uvm_ew(ilon,2,3,ilatm)=recvbuf(2,3,j)
  end do
  deallocate(recvbuf)
    
end subroutine inmi_coupler_sd2ew

subroutine inmi_coupler_ew2sd1(mype)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    inmi_coupler_ew2sd1
!
!   prgrmmr:
!
! abstract: use mpi_alltoallv to move uvm_ew (lat strips) to u_sd,v_sd,m_sd (subdomains)
!
! program history log:
!   2008-04-04  safford -- add subprogram doc block, rm unused vars and uses
!
!   input argument list:
!     mype     - mpi task id
!
!   output argument list:
!
! attributes:
!   language:  f90
!   machine:   ibm rs/6000 sp
!
!$$$

  implicit none


  integer(i_kind),intent(in   ) :: mype

  integer(i_kind) i,ilat,imode,j,k,mm1,ipe,ilon,mpi_string1,nn

  allocate(nsend_ew2sd(npe),nrecv_ew2sd(npe))
  allocate(ndsend_ew2sd(npe+1),ndrecv_ew2sd(npe+1))
  mm1=mype+1

!      1.  for each pe, gather up list of points from this set of lat strips destined
!             for subdomain of pe
  do ipe=1,npe
     nn=0
     do k=nlatm_0,nlatm_1
        ilat=mode_list(1,k)
        imode=mode_list(2,k)
        i=ilat-istart(ipe)+2
        if(i <  1.or.i >  ilat1(ipe)+2) cycle
        do j=1,jlon1(ipe)+2
           ilon=j+jstart(ipe)-2
           if(ilon <  1) ilon=ilon+nlon
           if(ilon >  nlon) ilon=ilon-nlon
           nn=nn+1
        end do
     end do
     nsend_ew2sd(ipe)=nn
  end do

  ndsend_ew2sd(1)=0
  do i=2,npe+1
     ndsend_ew2sd(i)=ndsend_ew2sd(i-1)+nsend_ew2sd(i-1)
  end do
  nallsend_ew2sd=ndsend_ew2sd(npe+1)
  allocate(info_send_ew2sd(3,nallsend_ew2sd))
  nn=0
  do ipe=1,npe
     do k=nlatm_0,nlatm_1
        ilat=mode_list(1,k)
        imode=mode_list(2,k)
        i=ilat-istart(ipe)+2
        if(i <  1.or.i >  ilat1(ipe)+2) cycle
        do j=1,jlon1(ipe)+2
           ilon=j+jstart(ipe)-2
           if(ilon <  1) ilon=ilon+nlon
           if(ilon >  nlon) ilon=ilon-nlon
           nn=nn+1
           info_send_ew2sd(1,nn)=ilon
           info_send_ew2sd(2,nn)=j
           info_send_ew2sd(3,nn)=k
        end do
     end do
  end do

  call mpi_alltoall(nsend_ew2sd,1,mpi_integer4,nrecv_ew2sd,1,mpi_integer4,mpi_comm_world,ierror)
  ndrecv_ew2sd(1)=0
  do i=2,npe+1
     ndrecv_ew2sd(i)=ndrecv_ew2sd(i-1)+nrecv_ew2sd(i-1)
  end do
  nallrecv_ew2sd=ndrecv_ew2sd(npe+1)
  allocate(info_recv_ew2sd(3,nallrecv_ew2sd))
  call mpi_type_contiguous(3,mpi_integer4,mpi_string1,ierror)
  call mpi_type_commit(mpi_string1,ierror)
  call mpi_alltoallv(info_send_ew2sd,nsend_ew2sd,ndsend_ew2sd,mpi_string1, &
                     info_recv_ew2sd,nrecv_ew2sd,ndrecv_ew2sd,mpi_string1,mpi_comm_world,ierror)
  call mpi_type_free(mpi_string1,ierror)

end subroutine inmi_coupler_ew2sd1

subroutine inmi_coupler_ew2sd(u_sd1,v_sd1,m_sd1,u_sd2,v_sd2,m_sd2,uvm_ew,mype)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    inmi_coupler_ew2sd
!
!   prgrmmr:
!
! abstract: use mpi_alltoallv to move uvm_ew (lat strips) to u_sd,v_sd,m_sd (subdomains)
!
! program history log:
!   2008-04-04  safford -- add subprogram doc block, rm unused vars
!
!   input argument list:
!     mype     - mpi task id
!     uvm_ew   -
!
!   output argument list:
!     u_sd1    -
!     v_sd1    -
!     m_sd1    -
!     u_sd2    -
!     v_sd2    -
!     m_sd2    -
!
! attributes:
!   language:  f90
!   machine:   ibm rs/6000 sp
!
!$$$

  implicit none


  integer(i_kind)                                 ,intent(in   ) :: mype
  real(r_kind),dimension(lat2,lon2,nvmodes_keep)  ,intent(  out) :: u_sd1,v_sd1,m_sd1
  real(r_kind),dimension(lat2,lon2,nvmodes_keep)  ,intent(  out) :: u_sd2,v_sd2,m_sd2
  real(r_kind),dimension(nlon,2,3,nlatm_0:nlatm_1),intent(in   ) :: uvm_ew

  real(r_kind),allocatable,dimension(:,:,:)::sendbuf,recvbuf
  integer(i_kind) ilat,imode,j,mm1,ilatm,ilon,mpi_string1,ilonloc

  mm1=mype+1

  allocate(sendbuf(2,3,nallsend_ew2sd))
  do j=1,nallsend_ew2sd
     ilon=info_send_ew2sd(1,j)
     ilatm=info_send_ew2sd(3,j)
     sendbuf(1,1,j)=uvm_ew(ilon,1,1,ilatm)
     sendbuf(1,2,j)=uvm_ew(ilon,1,2,ilatm)
     sendbuf(1,3,j)=uvm_ew(ilon,1,3,ilatm)
     sendbuf(2,1,j)=uvm_ew(ilon,2,1,ilatm)
     sendbuf(2,2,j)=uvm_ew(ilon,2,2,ilatm)
     sendbuf(2,3,j)=uvm_ew(ilon,2,3,ilatm)
  end do
  allocate(recvbuf(2,3,nallrecv_ew2sd))
  call mpi_type_contiguous(6,mpi_rtype,mpi_string1,ierror)
  call mpi_type_commit(mpi_string1,ierror)
  call mpi_alltoallv(sendbuf,nsend_ew2sd,ndsend_ew2sd,mpi_string1, &
                     recvbuf,nrecv_ew2sd,ndrecv_ew2sd,mpi_string1,mpi_comm_world,ierror)
  call mpi_type_free(mpi_string1,ierror)
  deallocate(sendbuf)
  do j=1,nallrecv_ew2sd
     ilonloc=info_recv_ew2sd(2,j)
     ilatm=info_recv_ew2sd(3,j)
     ilat=mode_list(1,ilatm)-istart(mm1)+2
     imode=mode_list(2,ilatm)
     u_sd1(ilat,ilonloc,imode)=recvbuf(1,1,j)
     v_sd1(ilat,ilonloc,imode)=recvbuf(1,2,j)
     m_sd1(ilat,ilonloc,imode)=recvbuf(1,3,j)
     u_sd2(ilat,ilonloc,imode)=recvbuf(2,1,j)
     v_sd2(ilat,ilonloc,imode)=recvbuf(2,2,j)
     m_sd2(ilat,ilonloc,imode)=recvbuf(2,3,j)
!--------------check for north or south pole
     ilat=-1
     if(mode_list(1,ilatm) == nlat) ilat=nlat-istart(mm1)+3
     if(mode_list(1,ilatm) == 1) ilat=2-istart(mm1)
     if(ilat == -1) cycle
!-----------------do repeat rows for north/south pole
     u_sd1(ilat,ilonloc,imode)=recvbuf(1,1,j)
     v_sd1(ilat,ilonloc,imode)=recvbuf(1,2,j)
     m_sd1(ilat,ilonloc,imode)=recvbuf(1,3,j)
     u_sd2(ilat,ilonloc,imode)=recvbuf(2,1,j)
     v_sd2(ilat,ilonloc,imode)=recvbuf(2,2,j)
     m_sd2(ilat,ilonloc,imode)=recvbuf(2,3,j)
  end do
  deallocate(recvbuf)

end subroutine inmi_coupler_ew2sd

subroutine inmi_ew_trans(uvm_ew,uvm_ewtrans)

!$$$  subprogram documentation block
!                .      .    .
! subprogram:    inmi_ew_trans
!
!   prgrmmr:
!
! abstract:
!
! program history log:
!   2008-04-04  safford -- add subprogram doc block, rm unused uses
!   2010-03-31  treadon - replace specmod components with sp_a structure
!
!   input argument list:
!     uvm_ew   -
!
!   output argument list:
!     uvm_ewtrans -
!
! attributes:
!   language:  f90
!   machine:   ibm rs/6000 sp
!
!$$$

  implicit none

  real(r_kind),dimension(nlon,2,3,nlatm_0:nlatm_1)    ,intent(in   ) :: uvm_ew
  real(r_kind),dimension(2,0:sp_a%jcap,2,3,nlatm_0:nlatm_1),intent(  out) :: uvm_ewtrans

  integer(i_kind) i,j,k
  real(r_kind),dimension(2,0:nlon/2,2)::halfwave
  real(r_kind),dimension(50000+4*sp_a%imax)::tmpafft

!$omp parallel do  schedule(dynamic,1) private(k,j,i,tmpafft,halfwave)
  do k=nlatm_0,nlatm_1
     tmpafft=sp_a%afft
     do j=1,3
        call spffte(nlon,1+nlon/2,nlon,2,halfwave,uvm_ew(1,1,j,k),-1,tmpafft)
        do i=0,sp_a%jcap
           uvm_ewtrans(1,i,1,j,k)=halfwave(1,i,1)
           uvm_ewtrans(2,i,1,j,k)=halfwave(2,i,1)
           uvm_ewtrans(1,i,2,j,k)=halfwave(1,i,2)
           uvm_ewtrans(2,i,2,j,k)=halfwave(2,i,2)
        end do
     end do
  end do

end subroutine inmi_ew_trans

subroutine inmi_ew_invtrans_ad(uvm_ew,uvm_ewtrans)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    inmi_ew_invtrans_ad
!
!   prgrmmr:
!
! abstract:
!
! program history log:
!   2008-04-04  safford -- add subprogram doc block, rm unused uses
!   2010-03-31  treadon - replace specmod components with sp_a structure
!
!   input argument list:
!     uvm_ew   -
!
!   output argument list:
!     uvm_ewtrans -
!
! attributes:
!   language:  f90
!   machine:   ibm rs/6000 sp
!
!$$$

  implicit none

  real(r_kind),dimension(nlon,2,3,nlatm_0:nlatm_1)    ,intent(in   ) :: uvm_ew
  real(r_kind),dimension(2,0:sp_a%jcap,2,3,nlatm_0:nlatm_1),intent(  out) :: uvm_ewtrans

  integer(i_kind) i,j,k
  real(r_kind) fnlon,fnlon2
  real(r_kind),dimension(2,0:nlon/2,2)::halfwave
  real(r_kind),dimension(50000+4*sp_a%imax)::tmpafft

  fnlon=real(nlon,r_kind)
  fnlon2=two*fnlon

!$omp parallel do  schedule(dynamic,1) private(k,j,i,tmpafft,halfwave)
  do k=nlatm_0,nlatm_1
     tmpafft=sp_a%afft
     do j=1,3
        call spffte(nlon,1+nlon/2,nlon,2,halfwave,uvm_ew(1,1,j,k),-1,tmpafft)
        uvm_ewtrans(1,0,1,j,k)=halfwave(1,0,1)*fnlon
        uvm_ewtrans(2,0,1,j,k)=halfwave(2,0,1)*fnlon
        uvm_ewtrans(1,0,2,j,k)=halfwave(1,0,2)*fnlon
        uvm_ewtrans(2,0,2,j,k)=halfwave(2,0,2)*fnlon
        do i=1,sp_a%jcap
           uvm_ewtrans(1,i,1,j,k)=halfwave(1,i,1)*fnlon2
           uvm_ewtrans(2,i,1,j,k)=halfwave(2,i,1)*fnlon2
           uvm_ewtrans(1,i,2,j,k)=halfwave(1,i,2)*fnlon2
           uvm_ewtrans(2,i,2,j,k)=halfwave(2,i,2)*fnlon2
        end do
     end do
  end do

end subroutine inmi_ew_invtrans_ad

subroutine inmi_ew_invtrans(uvm_ew,uvm_ewtrans)

!$$$  subprogram documentation block
!                .      .    .
! subprogram:    inmi_ew_invtrans
!
!   prgrmmr:
!
! abstract:
!
! program history log:
!   2008-04-04  safford -- add subprogram doc block, rm unused vars
!   2010-03-31  treadon - replace specmod components with sp_a structure
!
!   input argument list:
!     uvm_ewtrans -
!
!   output argument list:
!     uvm_ew      -
!
! attributes:
!   language:  f90
!   machine:   ibm rs/6000 sp
!
!$$$

  implicit none

  real(r_kind),dimension(nlon,2,3,nlatm_0:nlatm_1)    ,intent(  out) :: uvm_ew
  real(r_kind),dimension(2,0:sp_a%jcap,2,3,nlatm_0:nlatm_1),intent(in   ) :: uvm_ewtrans

  integer(i_kind) i,j,k
  real(r_kind),dimension(2,0:nlon/2,2)::halfwave
  real(r_kind),dimension(50000+4*sp_a%imax)::tmpafft

!$omp parallel do  schedule(dynamic,1) private(k,j,i,tmpafft,halfwave)
  do k=nlatm_0,nlatm_1
     tmpafft=sp_a%afft
     do j=1,3
        do i=0,sp_a%jcap
           halfwave(1,i,1)=uvm_ewtrans(1,i,1,j,k)
           halfwave(2,i,1)=uvm_ewtrans(2,i,1,j,k)
           halfwave(1,i,2)=uvm_ewtrans(1,i,2,j,k)
           halfwave(2,i,2)=uvm_ewtrans(2,i,2,j,k)
        end do
        do i=sp_a%jcap+1,nlon/2
           halfwave(1,i,1)=zero
           halfwave(2,i,1)=zero
           halfwave(1,i,2)=zero
           halfwave(2,i,2)=zero
        end do
        call spffte(nlon,1+nlon/2,nlon,2,halfwave,uvm_ew(1,1,j,k),1,tmpafft)
     end do
  end do

end subroutine inmi_ew_invtrans

subroutine inmi_ew_trans_ad(uvm_ew,uvm_ewtrans)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    inmi_ew_trans_ad
!
!   prgrmmr:
!
! abstract:
!
! program history log:
!   2008-04-04  safford -- add subprogram doc block, rm unused vars
!   2010-03-31  treadon - replace specmod components with sp_a structure
!
!   input argument list:
!     uvm_ewtrans -
!
!   output argument list:
!   output argument list:
!     uvm_ew      -
!
! attributes:
!   language:  f90
!   machine:   ibm rs/6000 sp
!
!$$$


  implicit none

  real(r_kind),dimension(nlon,2,3,nlatm_0:nlatm_1)    ,intent(  out) :: uvm_ew
  real(r_kind),dimension(2,0:sp_a%jcap,2,3,nlatm_0:nlatm_1),intent(in   ) :: uvm_ewtrans

  integer(i_kind) i,j,k
  real(r_kind) invnlon,invnlon2
  real(r_kind),dimension(2,0:nlon/2,2):: halfwave
  real(r_kind),dimension(50000+4*sp_a%imax)::tmpafft

  invnlon=one/real(nlon,r_kind)
  invnlon2=one/(two*real(nlon,r_kind))
!$omp parallel do  schedule(dynamic,1) private(k,j,i,tmpafft,halfwave)
  do k=nlatm_0,nlatm_1
     tmpafft=sp_a%afft
     do j=1,3
        halfwave(1,0,1)=uvm_ewtrans(1,0,1,j,k)*invnlon
        halfwave(2,0,1)=zero
        halfwave(1,0,2)=uvm_ewtrans(1,0,2,j,k)*invnlon
        halfwave(2,0,2)=zero
        do i=1,sp_a%jcap
           halfwave(1,i,1)=uvm_ewtrans(1,i,1,j,k)*invnlon2
           halfwave(2,i,1)=uvm_ewtrans(2,i,1,j,k)*invnlon2
           halfwave(1,i,2)=uvm_ewtrans(1,i,2,j,k)*invnlon2
           halfwave(2,i,2)=uvm_ewtrans(2,i,2,j,k)*invnlon2
        end do
        do i=sp_a%jcap+1,nlon/2
           halfwave(1,i,1)=zero
           halfwave(2,i,1)=zero
           halfwave(1,i,2)=zero
           halfwave(2,i,2)=zero
        end do
        call spffte(nlon,1+nlon/2,nlon,2,halfwave,uvm_ew(1,1,j,k),1,tmpafft)
     end do
  end do

end subroutine inmi_ew_trans_ad

subroutine inmi_coupler_ew2ns0(mype)

!$$$  subprogram documentation block
!                .      .    .
! subprogram:    inmi_coupler_ew2ns0
!
!   prgrmmr:
!
! abstract:
!
! program history log:
!   2008-04-04  safford -- add subprogram doc block, rm unused uses
!   2010-03-31  treadon - replace specmod jcap with sp_a structure
!
!   input argument list:
!     mype     - mpi task id
!
!   output argument list:
!
! attributes:
!   language:  f90
!   machine:   ibm rs/6000 sp
!
!$$$

  implicit none

  integer(i_kind),intent(in   ) :: mype

  integer(i_kind) k,kk,m,n,num_per_pe,total_groups,nn,kchk

!   in laying out by zonal wave number/vertical mode, have two types of groupings:

!   1.  jcap odd: 

!     group zonal wave numbers in pairs as follows:

!      (0,(jcap+1)/2)    then    ( m,jcap+1-m)  for m=1,(jcap-1)/2

!   2.  jcap even:

!     0  (+,- mode pair together),  then   (m,jcap+1-m ) for m=1,jcap/2  single modes

  total_groups=(sp_a%jcap+1)*nvmodes_keep
  num_per_pe=total_groups/npe
  if(mod(total_groups,npe)/=0) num_per_pe=num_per_pe+1
  if(mod(total_groups,npe)==0) then
     kchk=npe
  else
     kchk=mod(total_groups,npe)
  end if

  if(mod(sp_a%jcap,2) /= 0) then

!    case  jcap odd:
   
     nn=0
     do k=1,nvmodes_keep
        nn=nn+1
        mmode_list(1,nn)=0
        mmode_list(2,nn)=(sp_a%jcap+1)/2
        mmode_list(3,nn)=k
        mmode_list(4,nn)=k
        mmode_list(5,nn)=-1
        do m=1,(sp_a%jcap-1)/2
           nn=nn+1
           mmode_list(1,nn)=m
           mmode_list(2,nn)=sp_a%jcap+1-m
           mmode_list(3,nn)=k
           mmode_list(4,nn)=k
           mmode_list(5,nn)=-1
        end do
     end do
     do k=1,nvmodes_keep
        nn=nn+1
        mmode_list(1,nn)=0
        mmode_list(2,nn)=(sp_a%jcap+1)/2
        mmode_list(3,nn)=-k
        mmode_list(4,nn)=-k
        mmode_list(5,nn)=-1
        do m=1,(sp_a%jcap-1)/2
           nn=nn+1
           mmode_list(1,nn)=m
           mmode_list(2,nn)=sp_a%jcap+1-m
           mmode_list(3,nn)=-k
           mmode_list(4,nn)=-k
           mmode_list(5,nn)=-1
        end do
     end do

  else

!    case  jcap even:
   
     nn=0
     do k=1,nvmodes_keep
        nn=nn+1
        mmode_list(1,nn)=0
        mmode_list(2,nn)=0
        mmode_list(3,nn)=k
        mmode_list(4,nn)=-k
        mmode_list(5,nn)=-1
        do m=1,sp_a%jcap/2
           nn=nn+1
           mmode_list(1,nn)=m
           mmode_list(2,nn)=sp_a%jcap+1-m
           mmode_list(3,nn)=k
           mmode_list(4,nn)=k
           mmode_list(5,nn)=-1
        end do
     end do
     do k=1,nvmodes_keep
        do m=1,sp_a%jcap/2
           nn=nn+1
           mmode_list(1,nn)=m
           mmode_list(2,nn)=sp_a%jcap+1-m
           mmode_list(3,nn)=-k
           mmode_list(4,nn)=-k
           mmode_list(5,nn)=-1
        end do
     end do

  end if

  m_0=-1
  m_1=-2
  nn=0
  do n=1,npe
     if(n <= kchk) then
        kk=num_per_pe
     else
        kk=num_per_pe-1
     end if
     if(kk >  0) then
        if(mype+1 == n) then
           m_0=nn+1
           m_1=nn+kk
        end if
        do k=1,kk
           nn=nn+1
           mmode_list(5,nn)=n
        end do
     end if
  end do

end subroutine inmi_coupler_ew2ns0

subroutine inmi_coupler_ew2ns1(mype)

!$$$  subprogram documentation block
!                .      .    .
! subprogram:    inmi_coupler_ew2ns1
!
!   prgrmmr:
!
! abstract:
!
! program history log:
!   2008-04-04  safford -- add subprogram doc block
!   2010-03-31  treadon - replace specmod jcap with sp_a structure
!
!   input argument list:
!     mype     - mpi task id
!
!   output argument list:
!
! attributes:
!   language:  f90
!   machine:   ibm rs/6000 sp
!
!$$$

  implicit none

  integer(i_kind),intent(in   )::mype

  integer(i_kind) i,ip12,ipe,j,k,m,nn,m1,m2,ilat,imode,imode1,imode2
  integer(i_kind) mpi_string1,ibad,ibad0,loop
  real(r_kind),dimension(0:sp_a%jcap,-nvmodes_keep:nvmodes_keep)::mmode2_list

  allocate(nsend(npe),nrecv(npe),ndsend(npe+1),ndrecv(npe+1))
  nn=0
                
  mmode2_list=0
  do j=1,(sp_a%jcap+1)*nvmodes_keep
     m1=mmode_list(1,j)
     m2=mmode_list(2,j)
     imode1=mmode_list(3,j)
     imode2=mmode_list(4,j)
     if(imode1 == imode2) then
        if(mmode2_list(m1,imode1) /= 0.or.mmode2_list(m2,imode1) /= 0) then
           if(mype == 0) write(6,*)' problem in inmi_coupler_ew2ns'
           call mpi_finalize(i)
           stop
        end if
        mmode2_list(m1,imode1)=j
        mmode2_list(m2,imode1)=j
     end if
     if(m1 == m2) then
        if(mmode2_list(m1,imode1) /= 0.or.mmode2_list(m1,imode2) /= 0) then
           if(mype == 0) write(6,*)' problem in inmi_coupler_ew2ns'
           call mpi_finalize(i)
           stop
        end if
        mmode2_list(m1,imode1)=j
        mmode2_list(m1,imode2)=j
     end if
  end do
  do imode=-nvmodes_keep,nvmodes_keep
     if(imode == 0) cycle
     do m=0,sp_a%jcap
        if(mmode2_list(m,imode) == 0) then
           if(mype == 0) write(6,*)' problem in inmi_coupler_ew2ns'
           call mpi_finalize(i)
           stop
        end if
     end do
  end do

!  obtain counts of points to send to each pe from this pe

  nsend=0
  do k=nlatm_0,nlatm_1
     imode=mode_list(2,k)
     do m=0,sp_a%jcap
        j=mmode2_list(m,imode)
        ipe=mmode_list(5,j)
        nsend(ipe)=nsend(ipe)+1
        j=mmode2_list(m,-imode)
        ipe=mmode_list(5,j)
        nsend(ipe)=nsend(ipe)+1
     end do
  end do

  ndsend(1)=0
  do i=2,npe+1
     ndsend(i)=ndsend(i-1)+nsend(i-1)
  end do
  nallsend=ndsend(npe+1)
  allocate(info_send(6,nallsend))
  nsend=0
  ibad =0
  do k=nlatm_0,nlatm_1
     ilat=mode_list(1,k)
     do loop=1,2
        imode=mode_list(2,k)
        if(loop == 2) imode=-mode_list(2,k)
        do m=0,sp_a%jcap
           j=mmode2_list(m,imode)
           m1=mmode_list(1,j)
           m2=mmode_list(2,j)
           imode1=mmode_list(3,j)
           imode2=mmode_list(4,j)
           ipe=mmode_list(5,j)
           ip12=0
           if(m1 == m2.and.imode == imode1) ip12=1
           if(m1 == m2.and.imode == imode2) ip12=2
           if(imode1 == imode2.and.m == m1) ip12=1
           if(imode1 == imode2.and.m == m2) ip12=2
           if(ip12 == 0) ibad=ibad+1
           ipe=mmode_list(5,j)
           nsend(ipe)=nsend(ipe)+1
           info_send(1,ndsend(ipe)+nsend(ipe))=k
           info_send(2,ndsend(ipe)+nsend(ipe))=ilat
           info_send(3,ndsend(ipe)+nsend(ipe))=imode
           info_send(4,ndsend(ipe)+nsend(ipe))=m
           info_send(5,ndsend(ipe)+nsend(ipe))=j
           info_send(6,ndsend(ipe)+nsend(ipe))=ip12
        end do
     end do
  end do

  call mpi_allreduce(ibad,ibad0,1,mpi_integer4,mpi_sum,mpi_comm_world,ierror)
  if(ibad0 >  0) then
     if(mype == 0) write(0,*)' ibad = ',ibad0,'  inconsistency in inmi_coupler_ew2ns1'
     call mpi_finalize(ierror)
     stop
  end if

  call mpi_alltoall(nsend,1,mpi_integer4,nrecv,1,mpi_integer4,mpi_comm_world,ierror)
  ndrecv(1)=0
  do i=2,npe+1
     ndrecv(i)=ndrecv(i-1)+nrecv(i-1)
  end do
  nallrecv=ndrecv(npe+1)
  allocate(info_recv(6,nallrecv))
  call mpi_type_contiguous(6,mpi_integer4,mpi_string1,ierror)
  call mpi_type_commit(mpi_string1,ierror)
  call mpi_alltoallv(info_send,nsend,ndsend,mpi_string1, &
                     info_recv,nrecv,ndrecv,mpi_string1,mpi_comm_world,ierror)
  call mpi_type_free(mpi_string1,ierror)

end subroutine inmi_coupler_ew2ns1

subroutine inmi_coupler_ew2ns(uvm_ewtrans,uvm_ns)

!$$$  subprogram documentation block
!                .      .    .
! subprogram:    inmi_coupler_ew2ns
!
!   prgrmmr:
!
! abstract:
!
! program history log:
!   2008-04-04  safford -- add subprogram doc block, rm unused vars and uses
!   2010-03-31  treadon - replace specmod components with sp_a structure
!
!   input argument list:
!     uvm_ewtrans -
!
!   output argument list:
!     uvm_ns      -
!
! attributes:
!   language:  f90
!   machine:   ibm rs/6000 sp
!
!$$$

  implicit none

  real(r_kind),dimension(2,0:sp_a%jcap,2,3,nlatm_0:nlatm_1),intent(in   ) :: uvm_ewtrans
  real(r_kind),dimension(3,2,nlat,2,m_0:m_1)          ,intent(  out) :: uvm_ns

  integer(i_kind) ip12,j,m,mm,ilat,ilatm,imode,mpi_string1,loop
  real(r_kind),allocatable,dimension(:,:,:)::sendbuf,recvbuf

  allocate(sendbuf(3,2,nallsend))
  do j=1,nallsend
     ilatm=info_send(1,j)
     imode=info_send(3,j)
     m=info_send(4,j)
     loop=1
     if(imode <  0) loop=2
     sendbuf(1,1,j)=uvm_ewtrans(1,m,loop,1,ilatm)
     sendbuf(2,1,j)=uvm_ewtrans(1,m,loop,2,ilatm)
     sendbuf(3,1,j)=uvm_ewtrans(1,m,loop,3,ilatm)
     sendbuf(1,2,j)=uvm_ewtrans(2,m,loop,1,ilatm)
     sendbuf(2,2,j)=uvm_ewtrans(2,m,loop,2,ilatm)
     sendbuf(3,2,j)=uvm_ewtrans(2,m,loop,3,ilatm)
  end do
  allocate(recvbuf(3,2,nallrecv))
  call mpi_type_contiguous(6,mpi_rtype,mpi_string1,ierror)
  call mpi_type_commit(mpi_string1,ierror)
  call mpi_alltoallv(sendbuf,nsend,ndsend,mpi_string1, &
                     recvbuf,nrecv,ndrecv,mpi_string1,mpi_comm_world,ierror)
  call mpi_type_free(mpi_string1,ierror)
  deallocate(sendbuf)

  do j=1,nallrecv
     ilat=info_recv(2,j)
     mm=info_recv(5,j)
     ip12=info_recv(6,j)
     uvm_ns(1,1,ilat,ip12,mm)=recvbuf(1,1,j)
     uvm_ns(2,1,ilat,ip12,mm)=recvbuf(2,1,j)
     uvm_ns(3,1,ilat,ip12,mm)=recvbuf(3,1,j)
     uvm_ns(1,2,ilat,ip12,mm)=recvbuf(1,2,j)
     uvm_ns(2,2,ilat,ip12,mm)=recvbuf(2,2,j)
     uvm_ns(3,2,ilat,ip12,mm)=recvbuf(3,2,j)
  end do
  deallocate(recvbuf)

end subroutine inmi_coupler_ew2ns

subroutine inmi_coupler_ns2ew(uvm_ewtrans,uvm_ns)

!$$$  subprogram documentation block
!                .      .    .
! subprogram:    inmi_coupler_ns2ew
!
!   prgrmmr:
!
! abstract:
!
! program history log:
!   2008-04-04  safford -- add subprogram doc block, rm unused vars and uses
!   2010-03-31  treadon - replace specmod jcap with sp_a structure
!
!   input argument list:
!     uvm_ns   -
!
!   output argument list:
!     uvm_ewtrans -
!
! attributes:
!   language:  f90
!   machine:   ibm rs/6000 sp
!
!$$$

  implicit none

  real(r_kind),dimension(2,0:sp_a%jcap,2,3,nlatm_0:nlatm_1),intent(  out) :: uvm_ewtrans
  real(r_kind),dimension(3,2,nlat,2,m_0:m_1)          ,intent(in   ) :: uvm_ns

  integer(i_kind) ip12,j,m,mm,ilat,ilatm,imode,mpi_string1,loop
  real(r_kind),allocatable,dimension(:,:,:)::sendbuf,recvbuf


  allocate(recvbuf(3,2,nallrecv))
  do j=1,nallrecv
     ilat=info_recv(2,j)
     mm=info_recv(5,j)
     ip12=info_recv(6,j)
     recvbuf(1,1,j)=uvm_ns(1,1,ilat,ip12,mm)
     recvbuf(2,1,j)=uvm_ns(2,1,ilat,ip12,mm)
     recvbuf(3,1,j)=uvm_ns(3,1,ilat,ip12,mm)
     recvbuf(1,2,j)=uvm_ns(1,2,ilat,ip12,mm)
     recvbuf(2,2,j)=uvm_ns(2,2,ilat,ip12,mm)
     recvbuf(3,2,j)=uvm_ns(3,2,ilat,ip12,mm)
  end do
  allocate(sendbuf(3,2,nallsend))
  call mpi_type_contiguous(6,mpi_rtype,mpi_string1,ierror)
  call mpi_type_commit(mpi_string1,ierror)
  call mpi_alltoallv(recvbuf,nrecv,ndrecv,mpi_string1, &
                     sendbuf,nsend,ndsend,mpi_string1,mpi_comm_world,ierror)
  call mpi_type_free(mpi_string1,ierror)
  deallocate(recvbuf)
  do j=1,nallsend
     ilatm=info_send(1,j)
     imode=info_send(3,j)
     m=info_send(4,j)
     loop=1
     if(imode <  0) loop=2
     uvm_ewtrans(1,m,loop,1,ilatm)=sendbuf(1,1,j)
     uvm_ewtrans(1,m,loop,2,ilatm)=sendbuf(2,1,j)
     uvm_ewtrans(1,m,loop,3,ilatm)=sendbuf(3,1,j)
     uvm_ewtrans(2,m,loop,1,ilatm)=sendbuf(1,2,j)
     uvm_ewtrans(2,m,loop,2,ilatm)=sendbuf(2,2,j)
     uvm_ewtrans(2,m,loop,3,ilatm)=sendbuf(3,2,j)
  end do
  deallocate(sendbuf)

end subroutine inmi_coupler_ns2ew

subroutine inmi_nsuvm2zdm(uvm_ns,zdm_hat)

!$$$  subprogram documentation block
!                .      .    .
! subprogram:    inmi_nsuv2zdm
!
!   prgrmmr:
!
! abstract:
!
! program history log:
!   2008-04-04  safford -- add subprogram doc block, rm unused uses
!   2010-03-31  treadon - replace specmod components with sp_a structure
!   2012-11-23  parrish - replace calls to spanaly_ns with inline code
!
!   input argument list:
!     uvm_ns   -
!
!   output argument list:
!     zdm_had  -
!
! attributes:
!   language:  f90
!   machine:   ibm rs/6000 sp
!   machine:   ibm rs/6000 sp
!
!$$$

  implicit none

  real(r_kind),dimension(3,2,nlat,2,m_0:m_1),intent(in   ) :: uvm_ns
  real(r_kind),dimension(3,2,nlat,2,m_0:m_1),intent(  out) :: zdm_hat

  integer(i_kind) i,ics,j,jnorth,jsouth,m,mm,n,ipair
  real(r_kind),dimension(2,0:sp_a%jcap):: spcz,spcd,spcp
  real(r_kind),dimension(2,0:sp_a%jcap+1):: spcu,spcv
  real(r_kind),dimension(0:sp_a%jcap+1):: plnloc
  real(r_kind),dimension(2,2):: fu,fv,fp
  real(r_kind) f11u,f21u,f12u,f22u
  real(r_kind) f11v,f21v,f12v,f22v
  real(r_kind) f11p,f21p,f12p,f22p
  real(r_kind):: c1,c2

!$omp parallel do  schedule(dynamic,1) private(mm,ipair,m,ics,i,n) &
!$omp private(spcz,spcd,spcp,spcu,spcv,j,jnorth,jsouth,plnloc,fu,fv,fp,c1,c2) &
!$omp private(f11u,f21u,f12u,f22u,f11v,f21v,f12v,f22v,f11p,f21p,f12p,f22p)
  do mm=m_0,m_1
     do ipair=1,2
        m=mmode_list(ipair,mm)
        ics=1+m*(2*sp_a%jcap+3-m)/2-m

        do n=m,sp_a%jcap
           spcp(1,n)=zero
           spcp(2,n)=zero
           spcu(1,n)=zero
           spcu(2,n)=zero
           spcv(1,n)=zero
           spcv(2,n)=zero
        end do
        spcu(1,sp_a%jcap+1)=zero
        spcu(2,sp_a%jcap+1)=zero
        spcv(1,sp_a%jcap+1)=zero
        spcv(2,sp_a%jcap+1)=zero

        do j=sp_a%jb,sp_a%je
           jsouth=1+j
           jnorth=nlat-j

           c1=sp_a%wlat(j)
           c2=c1/sp_a%clat(j)
           fu(1,1)=uvm_ns(1,1,jnorth,ipair,mm)*c2
           fu(2,1)=uvm_ns(1,2,jnorth,ipair,mm)*c2
           fu(1,2)=uvm_ns(1,1,jsouth,ipair,mm)*c2
           fu(2,2)=uvm_ns(1,2,jsouth,ipair,mm)*c2
           fv(1,1)=uvm_ns(2,1,jnorth,ipair,mm)*c2
           fv(2,1)=uvm_ns(2,2,jnorth,ipair,mm)*c2
           fv(1,2)=uvm_ns(2,1,jsouth,ipair,mm)*c2
           fv(2,2)=uvm_ns(2,2,jsouth,ipair,mm)*c2
           fp(1,1)=uvm_ns(3,1,jnorth,ipair,mm)*c1
           fp(2,1)=uvm_ns(3,2,jnorth,ipair,mm)*c1
           fp(1,2)=uvm_ns(3,1,jsouth,ipair,mm)*c1
           fp(2,2)=uvm_ns(3,2,jsouth,ipair,mm)*c1
!           create plnloc

   
           do n=m,sp_a%jcap
              plnloc(n)=sp_a%pln(ics+n,j)
           end do
           plnloc(sp_a%jcap+1)=sp_a%plntop(m+1,j)
 
           f11u=fu(1,1)+fu(1,2)
           f21u=fu(2,1)+fu(2,2)
           f12u=fu(1,1)-fu(1,2)
           f22u=fu(2,1)-fu(2,2)
           f11v=fv(1,1)+fv(1,2)
           f21v=fv(2,1)+fv(2,2)
           f12v=fv(1,1)-fv(1,2)
           f22v=fv(2,1)-fv(2,2)
           f11p=fp(1,1)+fp(1,2)
           f21p=fp(2,1)+fp(2,2)
           f12p=fp(1,1)-fp(1,2)
           f22p=fp(2,1)-fp(2,2)
           do n=m,sp_a%jcap+1,2
              spcu(1,n)=spcu(1,n)+plnloc(n)*f11u
              spcu(2,n)=spcu(2,n)+plnloc(n)*f21u
              spcv(1,n)=spcv(1,n)+plnloc(n)*f11v
              spcv(2,n)=spcv(2,n)+plnloc(n)*f21v
           end do
           do n=m+1,sp_a%jcap+1,2
              spcu(1,n)=spcu(1,n)+plnloc(n)*f12u
              spcu(2,n)=spcu(2,n)+plnloc(n)*f22u
              spcv(1,n)=spcv(1,n)+plnloc(n)*f12v
              spcv(2,n)=spcv(2,n)+plnloc(n)*f22v
           end do
           do n=m,sp_a%jcap,2
              spcp(1,n)=spcp(1,n)+plnloc(n)*f11p
              spcp(2,n)=spcp(2,n)+plnloc(n)*f21p
           end do
           do n=m+1,sp_a%jcap,2
              spcp(1,n)=spcp(1,n)+plnloc(n)*f12p
              spcp(2,n)=spcp(2,n)+plnloc(n)*f22p
           end do

        end do

        call spuv2dz_ns(sp_a%jcap,m,ics, &
                spcu(1,m),spcv(1,m),spcu(1,sp_a%jcap+1),spcv(1,sp_a%jcap+1),spcd(1,m),spcz(1,m))

        i=0
        do n=m,sp_a%jcap
           i=i+1
           zdm_hat(1,1,i,ipair,mm)=spcz(1,n)*sp_a%enn1(ics+n)
           zdm_hat(1,2,i,ipair,mm)=spcz(2,n)*sp_a%enn1(ics+n)
           zdm_hat(2,1,i,ipair,mm)=spcd(1,n)*sp_a%enn1(ics+n)
           zdm_hat(2,2,i,ipair,mm)=spcd(2,n)*sp_a%enn1(ics+n)
           zdm_hat(3,1,i,ipair,mm)=spcp(1,n)
           zdm_hat(3,2,i,ipair,mm)=spcp(2,n)
        end do

     end do

  end do

end subroutine inmi_nsuvm2zdm

subroutine inmi_nszdm2uvm_ad(uvm_ns,zdm_hat)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    inmi_nszdm2uvm_ad
!
!   prgrmmr:
!
! abstract:
!
! program history log:
!   2008-04-04  safford -- add subprogram doc block, rm unused uses
!   2010-03-31  treadon - replace specmod components with sp_a structure
!   2012-11-23  parrish - replace calls to spanaly_ns with inline code
!
!   input argument list:
!     uvm_ns   -
!
!   output argument list:
!     zdm_had  -
!
! attributes:
!   language:  f90
!   machine:   ibm rs/6000 sp
!
!$$$

  implicit none

  real(r_kind),dimension(3,2,nlat,2,m_0:m_1),intent(in   ) :: uvm_ns
  real(r_kind),dimension(3,2,nlat,2,m_0:m_1),intent(  out) :: zdm_hat

  integer(i_kind) i,ics,j,jnorth,jsouth,m,mm,n,ipair
  real(r_kind),dimension(3,2,nlat)::uvm_ns_temp
  real(r_kind),dimension(2,0:sp_a%jcap):: spcz,spcd,spcp
  real(r_kind),dimension(2,0:sp_a%jcap+1):: spcu,spcv
  real(r_kind),dimension(0:sp_a%jcap+1):: plnloc
  real(r_kind),dimension(2,2):: fu,fv,fp
  real(r_kind) f11u,f21u,f12u,f22u
  real(r_kind) f11v,f21v,f12v,f22v
  real(r_kind) f11p,f21p,f12p,f22p
  real(r_kind):: c1

!$omp parallel do  schedule(dynamic,1) private(mm,ipair,m,ics,i,n) &
!$omp private(spcz,spcd,spcp,spcu,spcv,j,jnorth,jsouth,plnloc,fu,fv,fp,c1,uvm_ns_temp) &
!$omp private(f11u,f21u,f12u,f22u,f11v,f21v,f12v,f22v,f11p,f21p,f12p,f22p)
  do mm=m_0,m_1
     do ipair=1,2
        m=mmode_list(ipair,mm)
        ics=1+m*(2*sp_a%jcap+3-m)/2-m

        do n=m,sp_a%jcap
           spcp(1,n)=zero
           spcp(2,n)=zero
           spcu(1,n)=zero
           spcu(2,n)=zero
           spcv(1,n)=zero
           spcv(2,n)=zero
        end do
        spcu(1,sp_a%jcap+1)=zero
        spcu(2,sp_a%jcap+1)=zero
        spcv(1,sp_a%jcap+1)=zero
        spcv(2,sp_a%jcap+1)=zero

!  adjoint of set pole values
        uvm_ns_temp(:,:,:)=uvm_ns(:,:,:,ipair,mm)
        if(m == 0) then
           uvm_ns_temp(3,1,     2)=uvm_ns_temp(3,1,   1)+uvm_ns_temp(3,1,     2)
           uvm_ns_temp(3,1,nlat-1)=uvm_ns_temp(3,1,nlat)+uvm_ns_temp(3,1,nlat-1)
        else if(m == 1) then
           uvm_ns_temp(1,1,     2)=uvm_ns_temp(1,1,   1)+uvm_ns_temp(1,1,     2)
           uvm_ns_temp(1,2,     2)=uvm_ns_temp(1,2,   1)+uvm_ns_temp(1,2,     2)
           uvm_ns_temp(2,1,     2)=uvm_ns_temp(2,1,   1)+uvm_ns_temp(2,1,     2)
           uvm_ns_temp(2,2,     2)=uvm_ns_temp(2,2,   1)+uvm_ns_temp(2,2,     2)
           uvm_ns_temp(1,1,nlat-1)=uvm_ns_temp(1,1,nlat)+uvm_ns_temp(1,1,nlat-1)
           uvm_ns_temp(1,2,nlat-1)=uvm_ns_temp(1,2,nlat)+uvm_ns_temp(1,2,nlat-1)
           uvm_ns_temp(2,1,nlat-1)=uvm_ns_temp(2,1,nlat)+uvm_ns_temp(2,1,nlat-1)
           uvm_ns_temp(2,2,nlat-1)=uvm_ns_temp(2,2,nlat)+uvm_ns_temp(2,2,nlat-1)
        end if

        do j=sp_a%jb,sp_a%je
           jsouth=1+j
           jnorth=nlat-j

           c1=one/sp_a%clat(j)
           fu(1,1)=uvm_ns_temp(1,1,jnorth)*c1
           fu(2,1)=uvm_ns_temp(1,2,jnorth)*c1
           fu(1,2)=uvm_ns_temp(1,1,jsouth)*c1
           fu(2,2)=uvm_ns_temp(1,2,jsouth)*c1
           fv(1,1)=uvm_ns_temp(2,1,jnorth)*c1
           fv(2,1)=uvm_ns_temp(2,2,jnorth)*c1
           fv(1,2)=uvm_ns_temp(2,1,jsouth)*c1
           fv(2,2)=uvm_ns_temp(2,2,jsouth)*c1
           fp(1,1)=uvm_ns_temp(3,1,jnorth)
           fp(2,1)=uvm_ns_temp(3,2,jnorth)
           fp(1,2)=uvm_ns_temp(3,1,jsouth)
           fp(2,2)=uvm_ns_temp(3,2,jsouth)

!           create plnloc

   
           do n=m,sp_a%jcap
              plnloc(n)=sp_a%pln(ics+n,j)
           end do
           plnloc(sp_a%jcap+1)=sp_a%plntop(m+1,j)

           f11u=fu(1,1)+fu(1,2)
           f21u=fu(2,1)+fu(2,2)
           f12u=fu(1,1)-fu(1,2)
           f22u=fu(2,1)-fu(2,2)
           f11v=fv(1,1)+fv(1,2)
           f21v=fv(2,1)+fv(2,2)
           f12v=fv(1,1)-fv(1,2)
           f22v=fv(2,1)-fv(2,2)
           f11p=fp(1,1)+fp(1,2)
           f21p=fp(2,1)+fp(2,2)
           f12p=fp(1,1)-fp(1,2)
           f22p=fp(2,1)-fp(2,2)
           do n=m,sp_a%jcap+1,2
              spcu(1,n)=spcu(1,n)+plnloc(n)*f11u
              spcu(2,n)=spcu(2,n)+plnloc(n)*f21u
              spcv(1,n)=spcv(1,n)+plnloc(n)*f11v
              spcv(2,n)=spcv(2,n)+plnloc(n)*f21v
           end do
           do n=m+1,sp_a%jcap+1,2
              spcu(1,n)=spcu(1,n)+plnloc(n)*f12u
              spcu(2,n)=spcu(2,n)+plnloc(n)*f22u
              spcv(1,n)=spcv(1,n)+plnloc(n)*f12v
              spcv(2,n)=spcv(2,n)+plnloc(n)*f22v
           end do
           do n=m,sp_a%jcap,2
              spcp(1,n)=spcp(1,n)+plnloc(n)*f11p
              spcp(2,n)=spcp(2,n)+plnloc(n)*f21p
           end do
           do n=m+1,sp_a%jcap,2
              spcp(1,n)=spcp(1,n)+plnloc(n)*f12p
              spcp(2,n)=spcp(2,n)+plnloc(n)*f22p
           end do

        end do

        call spuv2dz_ns(sp_a%jcap,m,ics, &
                spcu(1,m),spcv(1,m),spcu(1,sp_a%jcap+1),spcv(1,sp_a%jcap+1),spcd(1,m),spcz(1,m))

        i=0
        if(m == 0) then
           i=i+1
           zdm_hat(1,1,i,ipair,mm)=zero
           zdm_hat(1,2,i,ipair,mm)=zero
           zdm_hat(2,1,i,ipair,mm)=zero
           zdm_hat(2,2,i,ipair,mm)=zero
           zdm_hat(3,1,i,ipair,mm)=spcp(1,0)
           zdm_hat(3,2,i,ipair,mm)=spcp(2,0)
        end if
        do n=max(1,m),sp_a%jcap
           i=i+1
           zdm_hat(1,1,i,ipair,mm)=spcz(1,n)
           zdm_hat(1,2,i,ipair,mm)=spcz(2,n)
           zdm_hat(2,1,i,ipair,mm)=spcd(1,n)
           zdm_hat(2,2,i,ipair,mm)=spcd(2,n)
           zdm_hat(3,1,i,ipair,mm)=spcp(1,n)
           zdm_hat(3,2,i,ipair,mm)=spcp(2,n)
        end do

     end do

  end do

end subroutine inmi_nszdm2uvm_ad

subroutine inmi_nszdm2uvm(uvm_ns,zdm_hat)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    inmi_szdm2uvm
!
!   prgrmmr:
!
! abstract:
!
! program history log:
!   2008-04-04  safford -- add subprogram doc block, rm unused uses
!   2010-03-31  treadon - replace specmod components with sp_a structure
!   2012-11-23  parrish - replace calls to spsynth_ns with inline code
!
!   input argument list:
!     uvm_ns   -
!
!   output argument list:
!     zdm_had  -
!
! attributes:
!   language:  f90
!   machine:   ibm rs/6000 sp
!
!$$$

  implicit none

  real(r_kind),dimension(3,2,nlat,2,m_0:m_1),intent(  out) :: uvm_ns
  real(r_kind),dimension(3,2,nlat,2,m_0:m_1),intent(in   ) :: zdm_hat

  integer(i_kind) i,ics,j,jnorth,jsouth,m,mm,n,ipair
  real(r_kind),dimension(2,0:sp_a%jcap):: spcz,spcd,spcp
  real(r_kind),dimension(2,0:sp_a%jcap+1):: spcu,spcv
  real(r_kind),dimension(0:sp_a%jcap+1):: plnloc
  real(r_kind),dimension(2,2):: fu,fv,fp
  real(r_kind) f1ur,f1ui,f2ur,f2ui
  real(r_kind) f1vr,f1vi,f2vr,f2vi
  real(r_kind) f1pr,f1pi,f2pr,f2pi

  real(r_kind):: c1

!$omp parallel do  schedule(dynamic,1) private(mm,ipair,m,ics,i,n) &
!$omp private(spcz,spcd,spcp,spcu,spcv,j,jnorth,jsouth,plnloc,fu,fv,fp,c1) &
!$omp private(f1ur,f1ui,f2ur,f2ui,f1vr,f1vi,f2vr,f2vi,f1pr,f1pi,f2pr,f2pi)
  do mm=m_0,m_1
     do ipair=1,2
        m=mmode_list(ipair,mm)
        ics=1+m*(2*sp_a%jcap+3-m)/2-m

!           gather up spcz, spcd, spcp

        i=0
        do n=m,sp_a%jcap
           i=i+1
           spcz(1,n)=zdm_hat(1,1,i,ipair,mm)
           spcz(2,n)=zdm_hat(1,2,i,ipair,mm)
           spcd(1,n)=zdm_hat(2,1,i,ipair,mm)
           spcd(2,n)=zdm_hat(2,2,i,ipair,mm)
           spcp(1,n)=zdm_hat(3,1,i,ipair,mm)
           spcp(2,n)=zdm_hat(3,2,i,ipair,mm)
        end do

!           convert to spcu, spcv

        call spdz2uv_ns(sp_a%jcap,m,ics, &
                spcd(1,m),spcz(1,m),spcu(1,m),spcv(1,m),spcu(1,sp_a%jcap+1),spcv(1,sp_a%jcap+1))

        do j=sp_a%jb,sp_a%je
           jsouth=1+j
           jnorth=nlat-j

!           create plnloc

           do n=m,sp_a%jcap
              plnloc(n)=sp_a%pln(ics+n,j)
           end do
           plnloc(sp_a%jcap+1)=sp_a%plntop(m+1,j)
 
!          obtain f

! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  zero out fourier coefficients.
           f1ur=zero
           f1ui=zero
           f2ur=zero
           f2ui=zero
           f1vr=zero
           f1vi=zero
           f2vr=zero
           f2vi=zero
           f1pr=zero
           f1pi=zero
           f2pr=zero
           f2pi=zero
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  synthesis over finite latitude.
!  for each zonal wavenumber, synthesize terms over total wavenumber.
!  synthesize even and odd polynomials separately.
           do n=m,sp_a%jcap+1,2
              f1ur=f1ur+plnloc(n)*spcu(1,n)
              f1ui=f1ui+plnloc(n)*spcu(2,n)
              f1vr=f1vr+plnloc(n)*spcv(1,n)
              f1vi=f1vi+plnloc(n)*spcv(2,n)
           enddo
           do n=m+1,sp_a%jcap+1,2
              f2ur=f2ur+plnloc(n)*spcu(1,n)
              f2ui=f2ui+plnloc(n)*spcu(2,n)
              f2vr=f2vr+plnloc(n)*spcv(1,n)
              f2vi=f2vi+plnloc(n)*spcv(2,n)
           enddo
           do n=m,sp_a%jcap,2
              f1pr=f1pr+plnloc(n)*spcp(1,n)
              f1pi=f1pi+plnloc(n)*spcp(2,n)
           enddo
           do n=m+1,sp_a%jcap,2
              f2pr=f2pr+plnloc(n)*spcp(1,n)
              f2pi=f2pi+plnloc(n)*spcp(2,n)
           enddo
!  separate fourier coefficients from each hemisphere.
!  odd polynomials contribute negatively to the southern hemisphere.
           fu(1,1)=f1ur+f2ur
           fu(2,1)=f1ui+f2ui
           fu(1,2)=f1ur-f2ur
           fu(2,2)=f1ui-f2ui
           fv(1,1)=f1vr+f2vr
           fv(2,1)=f1vi+f2vi
           fv(1,2)=f1vr-f2vr
           fv(2,2)=f1vi-f2vi
           fp(1,1)=f1pr+f2pr
           fp(2,1)=f1pi+f2pi
           fp(1,2)=f1pr-f2pr
           fp(2,2)=f1pi-f2pi
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

!          scatter back to output pairs of lats
 
           c1=one/sp_a%clat(j)
           uvm_ns(1,1,jnorth,ipair,mm)=fu(1,1)*c1
           uvm_ns(1,2,jnorth,ipair,mm)=fu(2,1)*c1
           uvm_ns(1,1,jsouth,ipair,mm)=fu(1,2)*c1
           uvm_ns(1,2,jsouth,ipair,mm)=fu(2,2)*c1
           uvm_ns(2,1,jnorth,ipair,mm)=fv(1,1)*c1
           uvm_ns(2,2,jnorth,ipair,mm)=fv(2,1)*c1
           uvm_ns(2,1,jsouth,ipair,mm)=fv(1,2)*c1
           uvm_ns(2,2,jsouth,ipair,mm)=fv(2,2)*c1
           uvm_ns(3,1,jnorth,ipair,mm)=fp(1,1)
           uvm_ns(3,2,jnorth,ipair,mm)=fp(2,1)
           uvm_ns(3,1,jsouth,ipair,mm)=fp(1,2)
           uvm_ns(3,2,jsouth,ipair,mm)=fp(2,2)

        end do

!  set pole values
        if(m == 0) then
           uvm_ns(1,1,1,ipair,mm)=zero
           uvm_ns(1,2,1,ipair,mm)=zero
           uvm_ns(2,1,1,ipair,mm)=zero
           uvm_ns(2,2,1,ipair,mm)=zero
           uvm_ns(3,1,1,ipair,mm)=uvm_ns(3,1,2,ipair,mm)
           uvm_ns(3,2,1,ipair,mm)=zero
           uvm_ns(1,1,nlat,ipair,mm)=zero
           uvm_ns(1,2,nlat,ipair,mm)=zero
           uvm_ns(2,1,nlat,ipair,mm)=zero
           uvm_ns(2,2,nlat,ipair,mm)=zero
           uvm_ns(3,1,nlat,ipair,mm)=uvm_ns(3,1,nlat-1,ipair,mm)
           uvm_ns(3,2,nlat,ipair,mm)=zero
        else if(m == 1) then
           uvm_ns(1,1,1,ipair,mm)=uvm_ns(1,1,2,ipair,mm)
           uvm_ns(1,2,1,ipair,mm)=uvm_ns(1,2,2,ipair,mm)
           uvm_ns(2,1,1,ipair,mm)=uvm_ns(2,1,2,ipair,mm)
           uvm_ns(2,2,1,ipair,mm)=uvm_ns(2,2,2,ipair,mm)
           uvm_ns(3,1,1,ipair,mm)=zero
           uvm_ns(3,2,1,ipair,mm)=zero
           uvm_ns(1,1,nlat,ipair,mm)=uvm_ns(1,1,nlat-1,ipair,mm)
           uvm_ns(1,2,nlat,ipair,mm)=uvm_ns(1,2,nlat-1,ipair,mm)
           uvm_ns(2,1,nlat,ipair,mm)=uvm_ns(2,1,nlat-1,ipair,mm)
           uvm_ns(2,2,nlat,ipair,mm)=uvm_ns(2,2,nlat-1,ipair,mm)
           uvm_ns(3,1,nlat,ipair,mm)=zero
           uvm_ns(3,2,nlat,ipair,mm)=zero
        else
           uvm_ns(1,1,1,ipair,mm)=zero
           uvm_ns(1,2,1,ipair,mm)=zero
           uvm_ns(2,1,1,ipair,mm)=zero
           uvm_ns(2,2,1,ipair,mm)=zero
           uvm_ns(3,1,1,ipair,mm)=zero
           uvm_ns(3,2,1,ipair,mm)=zero
           uvm_ns(1,1,nlat,ipair,mm)=zero
           uvm_ns(1,2,nlat,ipair,mm)=zero
           uvm_ns(2,1,nlat,ipair,mm)=zero
           uvm_ns(2,2,nlat,ipair,mm)=zero
           uvm_ns(3,1,nlat,ipair,mm)=zero
           uvm_ns(3,2,nlat,ipair,mm)=zero
        end if

     end do

  end do

end subroutine inmi_nszdm2uvm

subroutine inmi_nspcm_hat2pcm(pcm_ns,pcm_hat)

!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    inmi_nspcm_hat2pcm
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-08-13  lueken - added subprogram doc block
!   2010-03-31  treadon - replace specmod components with sp_a structure
!   2012-11-23  parrish - replace calls to spsynth_ns with inline code
!
!   input argument list:
!    pcm_hat
!
!   output argument list:
!    pcm_ns
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

  implicit none

  real(r_kind),dimension(3,2,nlat,2,m_0:m_1),intent(  out) :: pcm_ns
  real(r_kind),dimension(3,2,nlat,2,m_0:m_1),intent(in   ) :: pcm_hat

  integer(i_kind) i,ics,j,jnorth,jsouth,m,mm,n,ipair
  real(r_kind),dimension(2,0:sp_a%jcap):: spcp,spcc,spcm
  real(r_kind),dimension(0:sp_a%jcap+1):: plnloc
  real(r_kind),dimension(2,2):: fp,fc,fm
  real(r_kind) f1pr,f1pi,f2pr,f2pi
  real(r_kind) f1cr,f1ci,f2cr,f2ci
  real(r_kind) f1mr,f1mi,f2mr,f2mi

!$omp parallel do  schedule(dynamic,1) private(mm,ipair,m,ics,i,n) &
!$omp private(spcc,spcm,spcp,j,jnorth,jsouth,plnloc,fc,fm,fp) &
!$omp private(f1pr,f1pi,f2pr,f2pi,f1cr,f1ci,f2cr,f2ci,f1mr,f1mi,f2mr,f2mi)
  do mm=m_0,m_1
     do ipair=1,2
        m=mmode_list(ipair,mm)
        ics=1+m*(2*sp_a%jcap+3-m)/2-m

!           gather up spcp, spcc, spcm

        i=0
        do n=m,sp_a%jcap
           i=i+1
           spcp(1,n)=pcm_hat(1,1,i,ipair,mm)
           spcp(2,n)=pcm_hat(1,2,i,ipair,mm)
           spcc(1,n)=pcm_hat(2,1,i,ipair,mm)
           spcc(2,n)=pcm_hat(2,2,i,ipair,mm)
           spcm(1,n)=pcm_hat(3,1,i,ipair,mm)
           spcm(2,n)=pcm_hat(3,2,i,ipair,mm)
        end do

        do j=sp_a%jb,sp_a%je
           jsouth=1+j
           jnorth=nlat-j

!           create plnloc

           do n=m,sp_a%jcap
              plnloc(n)=sp_a%pln(ics+n,j)
           end do

!          obtain f

! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  zero out fourier coefficients.
           f1pr=zero
           f1pi=zero
           f2pr=zero
           f2pi=zero
           f1cr=zero
           f1ci=zero
           f2cr=zero
           f2ci=zero
           f1mr=zero
           f1mi=zero
           f2mr=zero
           f2mi=zero
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  synthesis over finite latitude.
!  for each zonal wavenumber, synthesize terms over total wavenumber.
!  synthesize even and odd polynomials separately.
           do n=m,sp_a%jcap,2
              f1pr=f1pr+plnloc(n)*spcp(1,n)
              f1pi=f1pi+plnloc(n)*spcp(2,n)
              f1cr=f1cr+plnloc(n)*spcc(1,n)
              f1ci=f1ci+plnloc(n)*spcc(2,n)
              f1mr=f1mr+plnloc(n)*spcm(1,n)
              f1mi=f1mi+plnloc(n)*spcm(2,n)
           enddo
           do n=m+1,sp_a%jcap,2
              f2pr=f2pr+plnloc(n)*spcp(1,n)
              f2pi=f2pi+plnloc(n)*spcp(2,n)
              f2cr=f2cr+plnloc(n)*spcc(1,n)
              f2ci=f2ci+plnloc(n)*spcc(2,n)
              f2mr=f2mr+plnloc(n)*spcm(1,n)
              f2mi=f2mi+plnloc(n)*spcm(2,n)
           enddo
!  separate fourier coefficients from each hemisphere.
!  odd polynomials contribute negatively to the southern hemisphere.
           fp(1,1)=f1pr+f2pr
           fp(2,1)=f1pi+f2pi
           fp(1,2)=f1pr-f2pr
           fp(2,2)=f1pi-f2pi
           fc(1,1)=f1cr+f2cr
           fc(2,1)=f1ci+f2ci
           fc(1,2)=f1cr-f2cr
           fc(2,2)=f1ci-f2ci
           fm(1,1)=f1mr+f2mr
           fm(2,1)=f1mi+f2mi
           fm(1,2)=f1mr-f2mr
           fm(2,2)=f1mi-f2mi

!          scatter back to output pairs of lats

           pcm_ns(1,1,jnorth,ipair,mm)=fp(1,1)
           pcm_ns(1,2,jnorth,ipair,mm)=fp(2,1)
           pcm_ns(1,1,jsouth,ipair,mm)=fp(1,2)
           pcm_ns(1,2,jsouth,ipair,mm)=fp(2,2)
           pcm_ns(2,1,jnorth,ipair,mm)=fc(1,1)
           pcm_ns(2,2,jnorth,ipair,mm)=fc(2,1)
           pcm_ns(2,1,jsouth,ipair,mm)=fc(1,2)
           pcm_ns(2,2,jsouth,ipair,mm)=fc(2,2)
           pcm_ns(3,1,jnorth,ipair,mm)=fm(1,1)
           pcm_ns(3,2,jnorth,ipair,mm)=fm(2,1)
           pcm_ns(3,1,jsouth,ipair,mm)=fm(1,2)
           pcm_ns(3,2,jsouth,ipair,mm)=fm(2,2)
 
        end do

!  set pole values
        if(m == 0) then
           pcm_ns(1,1,1,ipair,mm)=pcm_ns(1,1,2,ipair,mm)
           pcm_ns(1,2,1,ipair,mm)=zero
           pcm_ns(1,1,nlat,ipair,mm)=pcm_ns(1,1,nlat-1,ipair,mm)
           pcm_ns(1,2,nlat,ipair,mm)=zero
           pcm_ns(2,1,1,ipair,mm)=pcm_ns(2,1,2,ipair,mm)
           pcm_ns(2,2,1,ipair,mm)=zero
           pcm_ns(2,1,nlat,ipair,mm)=pcm_ns(2,1,nlat-1,ipair,mm)
           pcm_ns(2,2,nlat,ipair,mm)=zero
           pcm_ns(3,1,1,ipair,mm)=pcm_ns(3,1,2,ipair,mm)
           pcm_ns(3,2,1,ipair,mm)=zero
           pcm_ns(3,1,nlat,ipair,mm)=pcm_ns(3,1,nlat-1,ipair,mm)
           pcm_ns(3,2,nlat,ipair,mm)=zero
        else
           pcm_ns(1,1,1,ipair,mm)=zero
           pcm_ns(1,2,1,ipair,mm)=zero
           pcm_ns(2,1,1,ipair,mm)=zero
           pcm_ns(2,2,1,ipair,mm)=zero
           pcm_ns(3,1,1,ipair,mm)=zero
           pcm_ns(3,2,1,ipair,mm)=zero
           pcm_ns(1,1,nlat,ipair,mm)=zero
           pcm_ns(1,2,nlat,ipair,mm)=zero
           pcm_ns(2,1,nlat,ipair,mm)=zero
           pcm_ns(2,2,nlat,ipair,mm)=zero
           pcm_ns(3,1,nlat,ipair,mm)=zero
           pcm_ns(3,2,nlat,ipair,mm)=zero
        end if

     end do

  end do

end subroutine inmi_nspcm_hat2pcm

subroutine inmi_nspcm_hat2pcm_ad(pcm_ns,pcm_hat)

!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    inmi_nspcm_hat2pcm_ad
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-08-13  lueken - added subprogram doc block
!   2010-03-31  treadon - replace specmod components with sp_a structure
!   2012-11-23  parrish - replace calls to spanaly_ns with inline code
!
!   input argument list:
!    pcm_ns
!
!   output argument list:
!    pcm_hat
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block
  implicit none

  real(r_kind),dimension(3,2,nlat,2,m_0:m_1),intent(in   ) :: pcm_ns
  real(r_kind),dimension(3,2,nlat,2,m_0:m_1),intent(  out) :: pcm_hat

  real(r_kind),dimension(3,2,nlat)::pcm_ns_temp
  integer(i_kind) i,ics,j,jnorth,jsouth,m,mm,n,ipair
  real(r_kind),dimension(2,0:sp_a%jcap):: spcp,spcc,spcm
  real(r_kind),dimension(0:sp_a%jcap+1):: plnloc
  real(r_kind),dimension(2,2):: fp,fc,fm
  real(r_kind) f11p,f21p,f12p,f22p
  real(r_kind) f11c,f21c,f12c,f22c
  real(r_kind) f11m,f21m,f12m,f22m

  pcm_hat=zero
!$omp parallel do  schedule(dynamic,1) private(mm,ipair,m,ics,i,n) &
!$omp private(spcc,spcm,spcp,j,jnorth,jsouth,plnloc,fc,fm,fp,pcm_ns_temp) &
!$omp private(f11p,f21p,f12p,f22p,f11c,f21c,f12c,f22c,f11m,f21m,f12m,f22m)
  do mm=m_0,m_1
     do ipair=1,2
        m=mmode_list(ipair,mm)
        ics=1+m*(2*sp_a%jcap+3-m)/2-m
 
        do n=m,sp_a%jcap
           spcp(1,n)=zero
           spcp(2,n)=zero
           spcc(1,n)=zero
           spcc(2,n)=zero
           spcm(1,n)=zero
           spcm(2,n)=zero
        end do

!  adjoint of set pole values

        pcm_ns_temp(:,:,:)=pcm_ns(:,:,:,ipair,mm)
        if(m == 0) then
           pcm_ns_temp(1,1,     2)=pcm_ns_temp(1,1,   1)+pcm_ns_temp(1,1,     2)
           pcm_ns_temp(1,1,nlat-1)=pcm_ns_temp(1,1,nlat)+pcm_ns_temp(1,1,nlat-1)
           pcm_ns_temp(2,1,     2)=pcm_ns_temp(2,1,   1)+pcm_ns_temp(2,1,     2)
           pcm_ns_temp(2,1,nlat-1)=pcm_ns_temp(2,1,nlat)+pcm_ns_temp(2,1,nlat-1)
           pcm_ns_temp(3,1,     2)=pcm_ns_temp(3,1,   1)+pcm_ns_temp(3,1,     2)
           pcm_ns_temp(3,1,nlat-1)=pcm_ns_temp(3,1,nlat)+pcm_ns_temp(3,1,nlat-1)
        end if

        do j=sp_a%jb,sp_a%je
           jsouth=1+j
           jnorth=nlat-j
 
!          adjoint of scatter back to output pairs of lats

           fp(1,1)=pcm_ns_temp(1,1,jnorth)
           fp(2,1)=pcm_ns_temp(1,2,jnorth)
           fp(1,2)=pcm_ns_temp(1,1,jsouth)
           fp(2,2)=pcm_ns_temp(1,2,jsouth)
           fc(1,1)=pcm_ns_temp(2,1,jnorth)
           fc(2,1)=pcm_ns_temp(2,2,jnorth)
           fc(1,2)=pcm_ns_temp(2,1,jsouth)
           fc(2,2)=pcm_ns_temp(2,2,jsouth)
           fm(1,1)=pcm_ns_temp(3,1,jnorth)
           fm(2,1)=pcm_ns_temp(3,2,jnorth)
           fm(1,2)=pcm_ns_temp(3,1,jsouth)
           fm(2,2)=pcm_ns_temp(3,2,jsouth)
 
!           create plnloc

           do n=m,sp_a%jcap
              plnloc(n)=sp_a%pln(ics+n,j)
           end do

!          adjoint of obtain f

           f11p=fp(1,1)+fp(1,2)
           f21p=fp(2,1)+fp(2,2)
           f12p=fp(1,1)-fp(1,2)
           f22p=fp(2,1)-fp(2,2)
           f11c=fc(1,1)+fc(1,2)
           f21c=fc(2,1)+fc(2,2)
           f12c=fc(1,1)-fc(1,2)
           f22c=fc(2,1)-fc(2,2)
           f11m=fm(1,1)+fm(1,2)
           f21m=fm(2,1)+fm(2,2)
           f12m=fm(1,1)-fm(1,2)
           f22m=fm(2,1)-fm(2,2)
           do n=m,sp_a%jcap,2
              spcp(1,n)=spcp(1,n)+plnloc(n)*f11p
              spcp(2,n)=spcp(2,n)+plnloc(n)*f21p
              spcc(1,n)=spcc(1,n)+plnloc(n)*f11c
              spcc(2,n)=spcc(2,n)+plnloc(n)*f21c
              spcm(1,n)=spcm(1,n)+plnloc(n)*f11m
              spcm(2,n)=spcm(2,n)+plnloc(n)*f21m
           end do
           do n=m+1,sp_a%jcap,2
              spcp(1,n)=spcp(1,n)+plnloc(n)*f12p
              spcp(2,n)=spcp(2,n)+plnloc(n)*f22p
              spcc(1,n)=spcc(1,n)+plnloc(n)*f12c
              spcc(2,n)=spcc(2,n)+plnloc(n)*f22c
              spcm(1,n)=spcm(1,n)+plnloc(n)*f12m
              spcm(2,n)=spcm(2,n)+plnloc(n)*f22m
           end do
 
        end do

!       adjoint of gather up spcp, spcc, spcm

        i=0
        do n=m,sp_a%jcap
           i=i+1
           pcm_hat(1,1,i,ipair,mm)=spcp(1,n)
           pcm_hat(1,2,i,ipair,mm)=spcp(2,n)
           pcm_hat(2,1,i,ipair,mm)=spcc(1,n)
           pcm_hat(2,2,i,ipair,mm)=spcc(2,n)
           pcm_hat(3,1,i,ipair,mm)=spcm(1,n)
           pcm_hat(3,2,i,ipair,mm)=spcm(2,n)
        end do

     end do

  end do

end subroutine inmi_nspcm_hat2pcm_ad

subroutine inmi_nsuvm2zdm_ad(uvm_ns,zdm_hat)

!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    inmi_nsuvm2zdm_ad
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-08-13  lueken - added subprogram doc block
!   2010-03-31  treadon - replace specmod components with sp_a structure
!   2012-11-23  parrish - replace calls to spsynth_ns with inline code
!
!   input argument list:
!    uvm_ns
!    zdm_hat
!
!   output argument list:
!    uvm_ns
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

  implicit none

  real(r_kind),dimension(3,2,nlat,2,m_0:m_1),intent(inout) :: uvm_ns
  real(r_kind),dimension(3,2,nlat,2,m_0:m_1),intent(in   ) :: zdm_hat

  integer(i_kind) i,ics,j,jnorth,jsouth,m,mm,n,ipair
  real(r_kind),dimension(2,0:sp_a%jcap):: spcz,spcd,spcp
  real(r_kind),dimension(2,0:sp_a%jcap+1):: spcu,spcv
  real(r_kind),dimension(0:sp_a%jcap+1):: plnloc
  real(r_kind),dimension(2,2):: fu,fv,fp
  real(r_kind) f1ur,f1ui,f2ur,f2ui
  real(r_kind) f1vr,f1vi,f2vr,f2vi
  real(r_kind) f1pr,f1pi,f2pr,f2pi

  real(r_kind):: c1,c2

!$omp parallel do  schedule(dynamic,1) private(mm,ipair,m,ics,i,n) &
!$omp private(spcz,spcd,spcp,spcu,spcv,j,jnorth,jsouth,plnloc,fu,fv,fp,c1,c2) &
!$omp private(f1ur,f1ui,f2ur,f2ui,f1vr,f1vi,f2vr,f2vi,f1pr,f1pi,f2pr,f2pi)
  do mm=m_0,m_1
     do ipair=1,2
        m=mmode_list(ipair,mm)
        ics=1+m*(2*sp_a%jcap+3-m)/2-m

!           gather up spcz, spcd, spcp

        i=0
        do n=m,sp_a%jcap
           i=i+1
           spcz(1,n)=zdm_hat(1,1,i,ipair,mm)*sp_a%enn1(ics+n)
           spcz(2,n)=zdm_hat(1,2,i,ipair,mm)*sp_a%enn1(ics+n)
           spcd(1,n)=zdm_hat(2,1,i,ipair,mm)*sp_a%enn1(ics+n)
           spcd(2,n)=zdm_hat(2,2,i,ipair,mm)*sp_a%enn1(ics+n)
           spcp(1,n)=zdm_hat(3,1,i,ipair,mm)
           spcp(2,n)=zdm_hat(3,2,i,ipair,mm)
        end do

!           convert to spcu, spcv

        call spdz2uv_ns(sp_a%jcap,m,ics, &
                spcd(1,m),spcz(1,m),spcu(1,m),spcv(1,m),spcu(1,sp_a%jcap+1),spcv(1,sp_a%jcap+1))

        do j=sp_a%jb,sp_a%je
           jsouth=1+j
           jnorth=nlat-j
 
!           create plnloc

           do n=m,sp_a%jcap
              plnloc(n)=sp_a%pln(ics+n,j)
           end do
           plnloc(sp_a%jcap+1)=sp_a%plntop(m+1,j)

!          obtain f

! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  zero out fourier coefficients.
           f1ur=zero
           f1ui=zero
           f2ur=zero
           f2ui=zero
           f1vr=zero
           f1vi=zero
           f2vr=zero
           f2vi=zero
           f1pr=zero
           f1pi=zero
           f2pr=zero
           f2pi=zero
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  synthesis over finite latitude.
!  for each zonal wavenumber, synthesize terms over total wavenumber.
!  synthesize even and odd polynomials separately.
           do n=m,sp_a%jcap+1,2
              f1ur=f1ur+plnloc(n)*spcu(1,n)
              f1ui=f1ui+plnloc(n)*spcu(2,n)
              f1vr=f1vr+plnloc(n)*spcv(1,n)
              f1vi=f1vi+plnloc(n)*spcv(2,n)
           enddo
           do n=m+1,sp_a%jcap+1,2
              f2ur=f2ur+plnloc(n)*spcu(1,n)
              f2ui=f2ui+plnloc(n)*spcu(2,n)
              f2vr=f2vr+plnloc(n)*spcv(1,n)
              f2vi=f2vi+plnloc(n)*spcv(2,n)
           enddo
           do n=m,sp_a%jcap,2
              f1pr=f1pr+plnloc(n)*spcp(1,n)
              f1pi=f1pi+plnloc(n)*spcp(2,n)
           enddo
           do n=m+1,sp_a%jcap,2
              f2pr=f2pr+plnloc(n)*spcp(1,n)
              f2pi=f2pi+plnloc(n)*spcp(2,n)
           enddo
!  separate fourier coefficients from each hemisphere.
!  odd polynomials contribute negatively to the southern hemisphere.
           fu(1,1)=f1ur+f2ur
           fu(2,1)=f1ui+f2ui
           fu(1,2)=f1ur-f2ur
           fu(2,2)=f1ui-f2ui
           fv(1,1)=f1vr+f2vr
           fv(2,1)=f1vi+f2vi
           fv(1,2)=f1vr-f2vr
           fv(2,2)=f1vi-f2vi
           fp(1,1)=f1pr+f2pr
           fp(2,1)=f1pi+f2pi
           fp(1,2)=f1pr-f2pr
           fp(2,2)=f1pi-f2pi

!          scatter back to output pairs of lats

           c1=sp_a%wlat(j)
           c2=c1/sp_a%clat(j)
           uvm_ns(1,1,jnorth,ipair,mm)=fu(1,1)*c2
           uvm_ns(1,2,jnorth,ipair,mm)=fu(2,1)*c2
           uvm_ns(1,1,jsouth,ipair,mm)=fu(1,2)*c2
           uvm_ns(1,2,jsouth,ipair,mm)=fu(2,2)*c2
           uvm_ns(2,1,jnorth,ipair,mm)=fv(1,1)*c2
           uvm_ns(2,2,jnorth,ipair,mm)=fv(2,1)*c2
           uvm_ns(2,1,jsouth,ipair,mm)=fv(1,2)*c2
           uvm_ns(2,2,jsouth,ipair,mm)=fv(2,2)*c2
           uvm_ns(3,1,jnorth,ipair,mm)=fp(1,1)*c1
           uvm_ns(3,2,jnorth,ipair,mm)=fp(2,1)*c1
           uvm_ns(3,1,jsouth,ipair,mm)=fp(1,2)*c1
           uvm_ns(3,2,jsouth,ipair,mm)=fp(2,2)*c1

        end do

     end do

  end do

end subroutine inmi_nsuvm2zdm_ad
      subroutine spdz2uv_ns(m,l,ics,d,z,u,v,utop,vtop)
!$$$  subprogram documentation block
!
! subprogram:    spdz2uv_ns  compute winds from div and vort for one zonal wave number
!   prgmmr: iredell          org: w/nmc23     date: 92-10-31
!
! abstract: computes the wind components from divergence and vorticity
!           in spectral space.
!           subprogram speps should be called already.
!           if l is the zonal wavenumber, n is the total wavenumber,
!           eps(l,n)=sqrt((n**2-l**2)/(4*n**2-1)) and a is earth radius,
!           then the zonal wind component u is computed as
!             u(l,n)=-i*l/(n*(n+1))*a*d(l,n)
!                    +eps(l,n+1)/(n+1)*a*z(l,n+1)-eps(l,n)/n*a*z(l,n-1)
!           and the meridional wind component v is computed as
!             v(l,n)=-i*l/(n*(n+1))*a*z(l,n)
!                    -eps(l,n+1)/(n+1)*a*d(l,n+1)+eps(l,n)/n*a*d(l,n-1)
!           where d is divergence and z is vorticity.
!           u and v are weighted by the cosine of latitude.
!           extra terms are computed over top of the spectral domain.
!           advantage is taken of the fact that eps(l,l)=0
!           in order to vectorize over the entire spectral domain.
!           triangular truncation only
!
! program history log:
!   91-10-31  mark iredell
!   2006-09-05 parrish -- modify to do one zonal wave number only for parallel
!                         computation across processors by zonal wave number.
!   2010-05-14 derber  - make triangular truncation only
!
! usage:    call spdz2uv_ns(m,l,elonn1,eon,eontop,d,z,u,v,utop,vtop)
!
!   input argument list:
!     m        - integer spectral truncation
!     l        - zonal wave number
!     ics      - starting point is full spectral array
!     d        - real (2,l:m) divergence for zonal wave number l
!     z        - real (2,l:m) vorticity for zonal wave number l
!
!   output argument list:
!     u        - real (2,l:m) zonal wind (times coslat) for zonal wave number l
!     v        - real (2,l:m) merid wind (times coslat) for zonal wave number l
!     utop     - real (2) zonal wind (times coslat) over top for zonal wave number l
!     vtop     - real (2) merid wind (times coslat) over top for zonal wave number l
!
! attributes:
!   language: cray fortran
!
!$$$
      implicit none

      integer(i_kind),intent(in   ) :: m,l,ics
      real(r_kind)   ,intent(in   ) :: d(2,l:m),z(2,l:m)
      real(r_kind)   ,intent(  out) :: u(2,l:m),v(2,l:m)
      real(r_kind)   ,intent(  out) :: utop(2),vtop(2)

      integer(i_kind) n
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  compute winds in the spectral domain

      do n=l,m
         u(1,n)= sp_a%elonn1(ics+n)*d(2,n)
         u(2,n)=-sp_a%elonn1(ics+n)*d(1,n)
         v(1,n)= sp_a%elonn1(ics+n)*z(2,n)
         v(2,n)=-sp_a%elonn1(ics+n)*z(1,n)
      end do
      do n=l,m-1
         u(1,n)=u(1,n)+sp_a%eon(ics+n+1)*z(1,n+1)
         u(2,n)=u(2,n)+sp_a%eon(ics+n+1)*z(2,n+1)
         v(1,n)=v(1,n)-sp_a%eon(ics+n+1)*d(1,n+1)
         v(2,n)=v(2,n)-sp_a%eon(ics+n+1)*d(2,n+1)
      end do
      do n=l+1,m
         u(1,n)=u(1,n)-sp_a%eon(ics+n)*z(1,n-1)
         u(2,n)=u(2,n)-sp_a%eon(ics+n)*z(2,n-1)
         v(1,n)=v(1,n)+sp_a%eon(ics+n)*d(1,n-1)
         v(2,n)=v(2,n)+sp_a%eon(ics+n)*d(2,n-1)
      end do
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  compute winds over top of the spectral domain
      utop(1)=-sp_a%eontop(l+1)*z(1,m)
      utop(2)=-sp_a%eontop(l+1)*z(2,m)
      vtop(1)= sp_a%eontop(l+1)*d(1,m)
      vtop(2)= sp_a%eontop(l+1)*d(2,m)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      return
      end subroutine spdz2uv_ns
!-----------------------------------------------------------------------
      subroutine spuv2dz_ns(m,l,ics,u,v,utop,vtop,d,z)
!$$$  subprogram documentation block
!
! subprogram:    spuv2dz_ns  compute div,vort from winds for one zonal wave number
!   prgmmr: iredell          org: w/nmc23     date: 92-10-31
!
! abstract: computes the divergence and vorticity from wind components
!           in spectral space.
!           subprogram speps should be called already.
!           if l is the zonal wavenumber, n is the total wavenumber,
!           eps(l,n)=sqrt((n**2-l**2)/(4*n**2-1)) and a is earth radius,
!           then the divergence d is computed as
!             d(l,n)=i*l*a*u(l,n)
!                    +eps(l,n+1)*n*a*v(l,n+1)-eps(l,n)*(n+1)*a*v(l,n-1)
!           and the vorticity z is computed as
!             z(l,n)=i*l*a*v(l,n)
!                    -eps(l,n+1)*n*a*u(l,n+1)+eps(l,n)*(n+1)*a*u(l,n-1)
!           where u is the zonal wind and v is the meridional wind.
!           u and v are weighted by the secant of latitude.
!           extra terms are used over top of the spectral domain.
!           advantage is taken of the fact that eps(l,l)=0
!           in order to vectorize over the entire spectral domain.
!           triangular truncation only
!
! program history log:
!   91-10-31  mark iredell
!   2006-09-05 parrish -- modify to do one zonal wave number only for parallel
!                         computation across processors by zonal wave number.
!   2010-05-14 derber  - triangular truncation only
!
! usage:    call spuv2dz_ns(m,l,ics,u,v,utop,vtop,d,z)
!
!   input argument list:
!     l        - zonal wave number
!     m        - integer spectral truncation
!     ics      - starting point is full spectral array
!     u        - real (2,l:m) zonal wind (over coslat)
!     v        - real (2,l:m) merid wind (over coslat)
!     utop     - real (2) zonal wind (over coslat) over top
!     vtop     - real (2) merid wind (over coslat) over top
!
!   output argument list:
!     d        - real (2,l:m) divergence
!     z        - real (2,l:m) vorticity
!
! attributes:
!   language: cray fortran
!
!$$$
      implicit none

      integer(i_kind),intent(in   ) :: m,l,ics
      real(r_kind)   ,intent(in   ) :: u(2,l:m),v(2,l:m)
      real(r_kind)   ,intent(in   ) :: utop(2),vtop(2)
      real(r_kind)   ,intent(  out) :: d(2,l:m),z(2,l:m)

      integer(i_kind) n
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  compute terms from the spectral domain

      do n=l,m
         d(1,n)=-sp_a%elonn1(ics+n)*u(2,n)
         d(2,n)= sp_a%elonn1(ics+n)*u(1,n)
         z(1,n)=-sp_a%elonn1(ics+n)*v(2,n)
         z(2,n)= sp_a%elonn1(ics+n)*v(1,n)
      end do
      do n=l,m-1
         d(1,n)=d(1,n)+sp_a%eon(ics+n+1)*v(1,n+1)
         d(2,n)=d(2,n)+sp_a%eon(ics+n+1)*v(2,n+1)
         z(1,n)=z(1,n)-sp_a%eon(ics+n+1)*u(1,n+1)
         z(2,n)=z(2,n)-sp_a%eon(ics+n+1)*u(2,n+1)
      end do
      do n=l+1,m
         d(1,n)=d(1,n)-sp_a%eon(ics+n)*v(1,n-1)
         d(2,n)=d(2,n)-sp_a%eon(ics+n)*v(2,n-1)
         z(1,n)=z(1,n)+sp_a%eon(ics+n)*u(1,n-1)
         z(2,n)=z(2,n)+sp_a%eon(ics+n)*u(2,n-1)
      end do

!  compute terms from over top of the spectral domain
      n=m
      d(1,n)=d(1,n)+sp_a%eontop(l+1)*vtop(1)
      d(2,n)=d(2,n)+sp_a%eontop(l+1)*vtop(2)
      z(1,n)=z(1,n)-sp_a%eontop(l+1)*utop(1)
      z(2,n)=z(2,n)-sp_a%eontop(l+1)*utop(2)

      return
      end subroutine spuv2dz_ns

end module strong_fast_global_mod