module strong_slow_global_mod
!$$$   module documentation block
!                .      .    .                                       .
! module:    strong_slow_global_mod
!
! abstract: 
!          
! program history log:
!   2008-04-04  safford - add module doc block
!
! subroutines included;
!    init_strongvars_1                    --
!    initialize_strong_slow_global        --
!    strong_bal_correction_slow_global    -- compute balance adjustment
!    strong_bal_correction_slow_global_ad -- adjoint of strong_bal_correction
!    inmi_sub2grid                        -- special subdomain to horizontal grid
!    inmi_grid2sub                        -- reverse of inmi_sub2grid
!    get_mode_number                      -- define mode interleaving
!    gather_rmstends                      -- get BAL diagnostics
!
! variable definitions:
!     mode_number  - defines vertical mode layout for current processor
!     mode_number0 - defines how vertical modes are interleaved across processors
!                    in horizontal slab arrays
!
! attributes:
!   language:  f90
!   machine:   ibm RS/6000 SP
!$$$ end documentation block

  use kinds, only: i_kind
  use gridmod, only: nsig
  use mpimod, only: nuvlevs,nnnuvlevs
  implicit none

! set default to private
  private
! set subroutines to public
  public :: init_strongvars_1
  public :: initialize_strong_slow_global
  public :: strong_bal_correction_slow_global
  public :: strong_bal_correction_slow_global_ad
  public :: inmi_sub2grid
  public :: inmi_grid2sub
  public :: get_mode_number
  public :: gather_rmstends

  integer(i_kind),allocatable::mode_number(:),mode_number0(:)

contains

  subroutine init_strongvars_1(mype)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    init_strongvars_1
!
!   prgrmmr:
!
! abstract:
!
! 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

    call initialize_strong_slow_global(mype)

  end subroutine init_strongvars_1

subroutine initialize_strong_slow_global(mype)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    initialize_strong_slow_global
!
!   prgrmmr:
!
! abstract:
!
! 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

  allocate(mode_number(nuvlevs),mode_number0(nsig))
  call get_mode_number(mype)
!---------mode_number > 0 for 1st copy of u,v,m, reserved for correction delu,delv,delm
!---------mode_number < 0 for 2nd copy of u,v,m, reserved for grav part of tends, u_t_g,v_t_g,m_t_g

end subroutine initialize_strong_slow_global

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

!$$$  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-04-04  safford - rm unused vars
!   2008-10-08  derber  - modify to output streamfunction and vel. pot. and to not update time derivatives
!   2009-11-27  parrish - add uv_hyb_ens.  if present and true, then
!                          input/output variables psi=u, chi=v.
!   2010-03-31  treadon - replace specmod components with sp_a structure
!
!   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 psi on gaussian grid (subdomains)
!     chi      - output balanced perturbation chi 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 kinds, only: r_kind
  use mod_vtrans, only: depths,nvmodes_keep,vtrans,vtrans_inv
  use mod_inmi, only: m,gspeed,mmax,dinmi,gproj
  use gridmod, only: nlat,nlon,lat2,lon2,sp_a,grd_a
  use constants, only: izero,ione,zero
  use hybrid_ensemble_parameters, only: uv_hyb_ens
  implicit none

  integer(i_kind)                       ,intent(in   ) :: mype
  logical                               ,intent(in   ) :: bal_diagnostic,update,fullfield
  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(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,nsig)::delu,delv,delt
  real(r_kind),dimension(lat2,lon2)::delps
  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(lat2,lon2,nvmodes_keep)::utilde_t_g,vtilde_t_g,mtilde_t_g
  real(r_kind),dimension(nlat,nlon,nuvlevs)::uwork,vwork,mwork
  real(r_kind),dimension(nuvlevs)::rmstend_loc_uf,rmstend_g_loc_uf
  real(r_kind),dimension(nuvlevs)::rmstend_loc_f,rmstend_g_loc_f
  real(r_kind),dimension(sp_a%nc)::divhat,vorthat,mhat,deldivhat,delvorthat,delmhat
  real(r_kind) rmstend_all_uf,rmstend_all_g_uf,rmstend_all_f,rmstend_all_g_f

  integer(i_kind) i,j,k,kk,iad,mode
  logical filtered

  filtered=.true.

  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)

                      
  call inmi_sub2grid(utilde_t,utilde_t,uwork)
  call inmi_sub2grid(vtilde_t,vtilde_t,vwork)
  call inmi_sub2grid(mtilde_t,mtilde_t,mwork)

  rmstend_loc_uf=zero
  rmstend_g_loc_uf=zero
  rmstend_loc_f=zero
  rmstend_g_loc_f=zero
  do kk=1,nuvlevs
     mode=mode_number(kk)
     if(mode == izero) then
        do j=1,nlon
           do i=1,nlat
              uwork(i,j,kk)=zero
              vwork(i,j,kk)=zero
              mwork(i,j,kk)=zero
           end do
        end do
        cycle
     end if
      

!   3.  uwork,vwork,mwork    --> divhat,vorthat,mhat  (spectral transform)
!       (slabs)                  (slabs)

     vorthat=zero ; divhat=zero ; mhat=zero
     call general_uvg2zds(grd_a,sp_a,vorthat,divhat,uwork(1,1,kk),vwork(1,1,kk))
     call general_g2s0(grd_a,sp_a,mhat,mwork(1,1,kk))

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

     gspeed=sqrt(depths(abs(mode)))
     iad=ione
     do m=0,sp_a%jcap
        if(mode >  izero) then
!              here, delvorthat, etc contain field corrections necessary to zero out gravity component
!                                         of tendencies
           call dinmi(vorthat  (iad),divhat  (iad),mhat  (iad),&
                    delvorthat(iad),deldivhat(iad),delmhat(iad))
        else
!               here, delvorthat, etc contain gravity component of tendencies
           if(bal_diagnostic) &
              call gproj(vorthat(iad),divhat(iad),mhat(iad),delvorthat(iad),deldivhat(iad),delmhat(iad), &
                      rmstend_loc_uf(kk),rmstend_g_loc_uf(kk),.not.filtered)
           call gproj(vorthat(iad),divhat(iad),mhat(iad),delvorthat(iad),deldivhat(iad),delmhat(iad), &
                      rmstend_loc_f(kk),rmstend_g_loc_f(kk),filtered)
        end if
        iad=iad+2*(sp_a%jcap-m+ione)
     end do

!   5.  deldivhat,delvorthat,delmhat    -->  uwork,vwork,mwork   (spectral inverse transform)
!          (slabs)                             (slabs)

     do j=1,nlon
        do i=1,nlat
           uwork(i,j,kk)=zero
           vwork(i,j,kk)=zero
           mwork(i,j,kk)=zero
        end do
     end do
     if(uv_hyb_ens) then
        call general_zds2uvg(grd_a,sp_a,delvorthat,deldivhat,uwork(1,1,kk),vwork(1,1,kk))
     else
        call general_zds2pcg(grd_a,sp_a,delvorthat,deldivhat,uwork(1,1,kk),vwork(1,1,kk))
     end if
     call general_s2g0(grd_a,sp_a,delmhat,mwork(1,1,kk))

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

        if (fullfield) then
           write(6,*) 'STRONG_SLOW_GLOBAL:   FULL FIELD BALANCE DIAGNOSTICS --  '
        else
           write(6,*) 'STRONG_SLOW_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)
           write(6,'(" mode,rmstend_uf,rmstend_g_uf,rat = ",i5,2e14.4,2f10.4)') &
                              i,rmstend_uf(i),rmstend_g_uf(i),&
                              rmstend_g_uf(i)/(rmstend_uf(i)-rmstend_g_uf(i)), &
                              rmstend_g_uf(i)/(rmstend_uf(1)-rmstend_g_uf(1))
        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)
           write(6,'(" mode,rmstend_f,rmstend_g_f,rat = ",i5,2e14.4,2f10.4)') &
                              i,rmstend_f(i),rmstend_g_f(i),&
                              rmstend_g_f(i)/(rmstend_f(i)-rmstend_g_f(i)), &
                              rmstend_g_f(i)/(rmstend_f(1)-rmstend_g_f(1))
        end do
        write(6,'(" rmstend_all_uf,g_uf,rat = ",2e14.4,f10.4)') rmstend_all_uf,rmstend_all_g_uf, &
                                                 rmstend_all_g_uf/(rmstend_all_uf-rmstend_all_g_uf) 
        write(6,'(" rmstend_all_f,g_f,rat = ",2e14.4,f10.4)') rmstend_all_f,rmstend_all_g_f, &
                                                 rmstend_all_g_f/(rmstend_all_f-rmstend_all_g_f) 
     end if
  end if

  delutilde=zero
  delvtilde=zero
  delmtilde=zero
  call inmi_grid2sub(delutilde,utilde_t_g,uwork)
  call inmi_grid2sub(delvtilde,vtilde_t_g,vwork)
  call inmi_grid2sub(delmtilde,mtilde_t_g,mwork)
 
!   7.  delutilde,delvtilde,delmtilde  -->  delu,delv,delt,delps   (vertical mode inverse transform)
!       (subdomains)                      (subdomains)

  call vtrans_inv(delutilde,delvtilde,delmtilde,delu,delv,delt,delps)
!????????????????????????????????in here, can insert diagnostic based on utilde_t_g and utilde_t,etc.
! call vtrans_inv(utilde_t_g,vtilde_t_g,mtilde_t_g,u_t_g,v_t_g,t_t_g,ps_t_g)


!  update u,v,t,ps


  if(update) then
     do k=1,nsig
        do j=1,lon2
           do i=1,lat2
              psi(i,j,k)=psi(i,j,k)+delu(i,j,k)
              chi(i,j,k)=chi(i,j,k)+delv(i,j,k)
              t(i,j,k)=t(i,j,k)+delt(i,j,k)
           end do
        end do
     end do
     do j=1,lon2
        do i=1,lat2
           ps(i,j)=ps(i,j)+delps(i,j)
        end do
     end do
  end if

end subroutine strong_bal_correction_slow_global

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

!$$$  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
!   2009-11-27  parrish - add variable uv_hyb_ens.  if present and true, then
!                          input/output variables psi=u, chi=v.
!   2010-03-31  treadon - replace specmod components with sp_a structure
!
!   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
!     u        - input perturbation psi on gaussian grid (subdomains)
!     v        - 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 psi tendency on gaussian grid (subdomains)
!     v_t      - output perturbation chi 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 kinds, only: r_kind
  use mod_vtrans, only: depths,nvmodes_keep,vtrans_ad,vtrans_inv_ad
  use mod_inmi, only: m,gspeed,mmax,dinmi_ad,gproj_ad
  use gridmod, only: nlat,nlon,lat2,lon2,sp_a,grd_a
  use constants, only: izero,ione,zero
  use hybrid_ensemble_parameters, only: uv_hyb_ens
  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

  real(r_kind),dimension(lat2,lon2,nsig)::delu,delv,delt
  real(r_kind),dimension(lat2,lon2)::delps
  real(r_kind),dimension(lat2,lon2,nvmodes_keep)::utilde_t,vtilde_t,mtilde_t
  real(r_kind),dimension(lat2,lon2,nvmodes_keep)::utilde_t2,vtilde_t2,mtilde_t2
  real(r_kind),dimension(lat2,lon2,nvmodes_keep)::delutilde,delvtilde,delmtilde
  real(r_kind),dimension(lat2,lon2,nvmodes_keep)::utilde_t_g,vtilde_t_g,mtilde_t_g
  real(r_kind),dimension(nlat,nlon,nuvlevs)::uwork,vwork,mwork
  real(r_kind),dimension(sp_a%nc)::divhat,vorthat,mhat,deldivhat,delvorthat,delmhat

  integer(i_kind) mode,iad
  integer(i_kind) i,j,k,kk

  mmax=sp_a%jcap

!  adjoint of update u,v,t,ps

  do j=1,lon2
     do i=1,lat2
        delps(i,j)=ps(i,j)
     end do
  end do

  do k=1,nsig
     do j=1,lon2
        do i=1,lat2
           delu(i,j,k)=psi(i,j,k)
           delv(i,j,k)=chi(i,j,k)
           delt(i,j,k)=t(i,j,k)
        end do
     end do
  end do

!   7.  adjoint of delutilde,delvtilde,delmtilde  -->  delu,delv,delt,delps  (vert mode inverse transform)
!       (subdomains)                      (subdomains)

  delutilde=zero ; delvtilde=zero ; delmtilde=zero
  utilde_t_g=zero ; vtilde_t_g=zero ; mtilde_t_g=zero
  call vtrans_inv_ad(delutilde,delvtilde,delmtilde,delu,delv,delt,delps)

  call inmi_sub2grid(delutilde,utilde_t_g,uwork)
  call inmi_sub2grid(delvtilde,vtilde_t_g,vwork)
  call inmi_sub2grid(delmtilde,mtilde_t_g,mwork)

  do kk=1,nuvlevs
     mode=mode_number(kk)
     if(mode == izero) then
        do j=1,nlon
           do i=1,nlat
              uwork(i,j,kk)=zero
              vwork(i,j,kk)=zero
              mwork(i,j,kk)=zero
           end do
        end do
        cycle
     end if


!   5.  adjoint of deldivhat,delvorthat,delmhat    -->  uwork,vwork,mwork   (spectral inverse transform)
!                     (slabs)                             (slabs)

     call general_s2g0_ad(grd_a,sp_a,delmhat,mwork(1,1,kk))
     if(uv_hyb_ens) then
        call general_zds2uvg_ad(grd_a,sp_a,delvorthat,deldivhat,uwork(1,1,kk),vwork(1,1,kk))
     else
        call general_zds2pcg_ad(grd_a,sp_a,delvorthat,deldivhat,uwork(1,1,kk),vwork(1,1,kk))
     end if

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

     gspeed=sqrt(depths(abs(mode)))
     iad=ione
     vorthat=zero ; divhat=zero ; mhat=zero
     do m=0,sp_a%jcap
        if(mode >  izero) then
           call dinmi_ad(vorthat(iad),divhat(iad),mhat(iad),&
                       delvorthat(iad)   ,   deldivhat(iad),   delmhat(iad))
        else
           call gproj_ad(vorthat(iad),divhat(iad),mhat(iad),delvorthat(iad),deldivhat(iad),delmhat(iad))
        end if
        iad=iad+2*(sp_a%jcap-m+ione)
     end do

!   3.  adjoint of uwork,vwork,mwork    --> divhat,vorthat,mhat  (spectral transform)
!                       (slabs)                  (slabs)

     call general_g2s0_ad(grd_a,sp_a,mhat,mwork(1,1,kk))
     call general_uvg2zds_ad(grd_a,sp_a,vorthat,divhat,uwork(1,1,kk),vwork(1,1,kk))

  end do

!???????????????????????this part needs some thought--probably need to have two distinct arrays,
!  which are then added together after these calls
  call inmi_grid2sub(utilde_t,utilde_t2,uwork)
  call inmi_grid2sub(vtilde_t,vtilde_t2,vwork)
  call inmi_grid2sub(mtilde_t,mtilde_t2,mwork)
  utilde_t=utilde_t+utilde_t2
  vtilde_t=vtilde_t+vtilde_t2
  mtilde_t=mtilde_t+mtilde_t2
!
!    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)

end subroutine strong_bal_correction_slow_global_ad

subroutine inmi_sub2grid(utilde,utilde2,uwork)

!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    inmi_sub2grid  special subdomain to horizontal grid
!   prgmmr: parrish          org: np23                date: 2006-08-03
!
! abstract: convert two subdomain 3d fields utilde,utilde2 from subdomains
!           to horizontal slab array uwork.  note that must have 
!           2*nvmodes_keep <= nsig.
!
! program history log:
!   2006-08-03  parrish
!   2008-04-04  safford  - rm unused uses
!
!   input argument list:
!     utilde   - vertical tranformed variable on subdomains
!     utilde2  - vertical tranformed variable on subdomains
!
!   output argument list:
!     uwork    - output fields in horizontal slab mode
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

  use kinds, only: r_kind
  use constants, only: izero,zero
  use gridmod, only: lat2,iglobal,lon1,itotsub,lon2,lat1,ltosi,ltosj,nlon,nlat,&
       strip,reorder
  use mpimod, only: ierror,mpi_comm_world,&
       isduv_g,iscuv_g,irduv_g,ircuv_g,mpi_rtype
  use mod_vtrans, only: nvmodes_keep
  implicit none

! Declare passed variables
  real(r_kind),dimension(lat2,lon2,nvmodes_keep),intent(in   ) :: utilde,utilde2
  real(r_kind),dimension(nlat,nlon,nuvlevs)     ,intent(  out) :: uwork

  real(r_kind),dimension(lat2,lon2,nsig):: u

! Declare local variables
  integer(i_kind) i,j,k,isize,mode

  real(r_kind),dimension(lat1,lon1,nsig):: usm
  real(r_kind),dimension(itotsub,nuvlevs):: work3

! Initialize variables
  isize=max(iglobal,itotsub)

! zero out work array
  do k=1,nuvlevs
     do j=1,isize
        work3(j,k)=zero
     end do
  end do

!  copy input array into bigger internal array
  do k=1,nsig
     do j=1,lon2
        do i=1,lat2
           u(i,j,k)=zero
        end do
     end do
  end do
  do k=1,nsig
     mode=mode_number0(k)
     if(mode == izero) cycle
     if(mode >  izero) then
        do j=1,lon2
           do i=1,lat2
              u(i,j,k)=utilde(i,j,mode)
           end do
        end do
     else
        do j=1,lon2
           do i=1,lat2
              u(i,j,k)=utilde2(i,j,-mode)
           end do
        end do
     end if
    
  end do

!   begin strmfctn / vp to u/v section
!   strip off endpoints of input arrays on subdomains
!   note that right now, place in usm,vsm
  call strip(u,usm,nsig)

!   put on global slabs
  call mpi_alltoallv(usm(1,1,1),iscuv_g,isduv_g,&
       mpi_rtype,work3(1,1),ircuv_g,irduv_g,mpi_rtype,&
       mpi_comm_world,ierror)

!   reorder work arrays and transfer to output array
  call reorder(work3,nuvlevs,nnnuvlevs)
  do k=1,nnnuvlevs
     do i=1,iglobal
        uwork(ltosi(i),ltosj(i),k)=work3(i,k)
     end do
  end do

end subroutine inmi_sub2grid

subroutine inmi_grid2sub(utilde,utilde2,uwork)

!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    inmi_grid2sub  reverse of inmi_sub2grid
!   prgmmr: parrish          org: np23                date: 2006-08-03
!
! abstract: reverse of inmi_sub2grid
!
! program history log:
!   2006-08-03  parrish
!   2008-04-04  safford - rm unused vars and uses
!
!   input argument list:
!     uwork    - input fields in horizontal slab mode
!
!   output argument list:
!     utilde   - vertical tranformed variable on subdomains
!     utilde2  - vertical tranformed variable on subdomains
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

  use kinds, only: r_kind
  use gridmod, only: lat2,iglobal,itotsub,lon2,nlat,nlon,ltosi_s,ltosj_s,&
       reorder2
  use mpimod, only: iscuv_s,ierror,mpi_comm_world,irduv_s,ircuv_s,&
       isduv_s,mpi_rtype
  use mod_vtrans, only: nvmodes_keep
  use constants, only: izero
  implicit none

! Declare passed variables
  real(r_kind),dimension(lat2,lon2,nvmodes_keep),intent(  out) :: utilde,utilde2
  real(r_kind),dimension(nlat,nlon,nuvlevs)     ,intent(in   ) :: uwork

  real(r_kind),dimension(lat2,lon2,nsig):: u

! Declare local variables
  integer(i_kind) i,j,k,isize,mode

  real(r_kind),dimension(itotsub,nuvlevs):: work3

! Initialize variables
  isize=max(iglobal,itotsub)

  do k=1,nnnuvlevs
     do i=1,itotsub
        work3(i,k)=uwork(ltosi_s(i),ltosj_s(i),k)
     end do
  end do
!   reorder the work array for the mpi communication
  call reorder2(work3,nuvlevs,nnnuvlevs)

!   get u back on subdomains
  call mpi_alltoallv(work3(1,1),iscuv_s,isduv_s,&
       mpi_rtype,u(1,1,1),ircuv_s,irduv_s,mpi_rtype,&
       mpi_comm_world,ierror)
  do k=1,nsig
     mode=mode_number0(k)
     if(mode == izero) cycle
     if(mode >  izero) then
        do j=1,lon2
           do i=1,lat2
              utilde(i,j,mode)=u(i,j,k)
           end do
        end do
     else
        do j=1,lon2
           do i=1,lat2
              utilde2(i,j,-mode)=u(i,j,k)
           end do
        end do
     end if
  end do

end subroutine inmi_grid2sub

subroutine get_mode_number(mype)

!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    get_mode_number  define mode interleaving
!   prgmmr: parrish          org: np23                date: 2006-08-03
!
! abstract: define how modes are interleaved across processors
!           when converting vertical transformed variables between
!           subdomains and horizontal slab storage
!
! program history log:
!   2006-08-03  parrish
!
!   input argument list:
!     mype     - current processor number
!
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

  use kinds, only: r_kind
  use mpimod, only: ierror,mpi_comm_world,mpi_integer,npe
  use mod_vtrans, only: nvmodes_keep
  use constants, only: izero,ione
  implicit none

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

  integer(i_kind) i,ii
  integer(i_kind) nuvlevs0(npe),ndisp(npe+ione)
  integer(i_kind) nuvlev_use,kchk

  if(nvmodes_keep*2 >  nsig) then
     write(6,*)' error in get_mode_number, currently necessary for nvmodes_keep*2 <= nsig '
     call stop2(89)
  end if
  if (mod(nsig,npe)==izero) then
     kchk=npe
  else
     kchk=mod(nsig,npe)
  end if
  if(mype+ione <= kchk) then
     nuvlev_use=nuvlevs
  else
     nuvlev_use=nuvlevs-ione
  end if
  ii=izero
  mode_number=izero
  do i=1,2*nvmodes_keep
     if(mod(i-ione,npe) == mype.and.ii <  nuvlev_use) then
        ii=ii+ione
        if(i <= nvmodes_keep) then
           mode_number(ii)=i
        else
           mode_number(ii)=nvmodes_keep-i
        end if
     end if
  end do
  call mpi_allgather(nuvlev_use,ione,mpi_integer,nuvlevs0,ione,mpi_integer,mpi_comm_world,ierror)
  ndisp(1)=izero
  do i=2,npe+ione
     ndisp(i)=ndisp(i-ione)+nuvlevs0(i-ione)
  end do
  call mpi_allgatherv(mode_number,nuvlev_use,mpi_integer, &
                   mode_number0,nuvlevs0,ndisp,mpi_integer,mpi_comm_world,ierror)

end subroutine get_mode_number

subroutine gather_rmstends(rmstend_loc,rmstend,mype)

!$$$  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
!
!   input argument list:
!     rmstend_loc - previously computed energy norms of vertical modes
!                   as distributed on local processors
!     mype     - current processor number
!
!   output argument list:
!     rmstend  -  all vertical modes of rmstend assembled across all processors
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

  use kinds, only: r_kind
  use mpimod, only: ierror,mpi_comm_world,mpi_integer,mpi_rtype,npe
  use mod_vtrans, only: nvmodes_keep
  use constants, only: izero,ione
  implicit none

  real(r_kind)   ,intent(in   ) :: rmstend_loc(nuvlevs)
  real(r_kind)   ,intent(  out) :: rmstend(nvmodes_keep)
  integer(i_kind),intent(in   ) :: mype

  integer(i_kind) i
  integer(i_kind) nuvlevs0(npe),ndisp(npe+ione)
  integer(i_kind) nuvlev_use,kchk
  real(r_kind) work(nsig)

  if (mod(nsig,npe)==izero) then
     kchk=npe
  else
     kchk=mod(nsig,npe)
  end if
  if(mype+ione <= kchk) then
     nuvlev_use=nuvlevs
  else
     nuvlev_use=nuvlevs-ione
  end if
  call mpi_allgather(nuvlev_use,ione,mpi_integer,nuvlevs0,ione,mpi_integer,mpi_comm_world,ierror)
  ndisp(1)=izero
  do i=2,npe+ione
     ndisp(i)=ndisp(i-ione)+nuvlevs0(i-ione)
  end do
  call mpi_allgatherv(rmstend_loc,nuvlev_use,mpi_rtype, &
                     work,nuvlevs0,ndisp,mpi_rtype,mpi_comm_world,ierror)
  do i=1,nsig
     if(mode_number0(i) <  izero) rmstend(-mode_number0(i))=work(i)
  end do

end subroutine gather_rmstends

end module strong_slow_global_mod