!> @file
!>
!> @brief rtsig() reads and transforms sigma file.
!>
!> This subprogram reads a sigma file and transforms
!> the fields to a designated global grid.
!> Add Iredells subroutine to read sigma files.
!>
!> @param[out] lusig integer(sigio_intkind) sigma file unit number.
!> @param[out] head type(sigio_head) sigma file header.
!> @param[out] k1 integer first model level to return.
!> @param[out] k2 integer last model level to return.
!> @param[out] kgds integer (200) GDS to which to transform.
!> @param[out] ijo integer dimension of output fields.
!> @param[out] levs integer number of total vertical levels.
!> @param[out] ntrac integer number of output tracers.
!> @param[out] jcap integer number of waves.
!> @param[out] lnt2 integer (jcap+1)*(jcap+2).
!> @param[out] h real (ijo) surface orography (m).
!> @param[out] p real (ijo) surface pressure (Pa).
!> @param[out] px real (ijo) log surface pressure x-gradient (1/m).
!> @param[out] py real (ijo) log surface pressure y-gradient (1/m).
!> @param[out] t real (ijo,k1:k2) temperature (K).
!> @param[out] tx real (ijo,k1:k2) virtual temperature x-gradient (K/m).
!> @param[out] ty real (ijo,k1:k2) virtual temperature y-gradient (K/m).
!> @param[out] u real (ijo,k1:k2) x-component wind (m/s).
!> @param[out] v real (ijo,k1:k2) y-component wind (m/s).
!> @param[out] d real (ijo,k1:k2) wind divergence (1/s).
!> @param[out] trc real (ijo,k1:k2,ntrac) tracers.
!><pre>
!>                                   1 = specific humidity (kg/kg)
!>                                   2 = Ozone mixing ratio (kg/kg)
!>                                   3 = cloud condensate mixing ratio (kg/kg)
!>                                       atomic oxyge, oxygen etc
!>
!></pre>
!> @param[out] iret Integer return code.
!>
!> ### Program History Log
!> Date | Programmer | Comments
!> -----|------------|---------
!> 1999-10-18 | Mark Iredell      | Initial
!> 2013-04-19 | Jun  Wang         | Add option to get tmp and ps(in pascal) from enthalpy and ps(cb) option
!> 2013-05-06 | Shrinivas Moorthi | Initialize midea to 0
!> 2013-05-07 | Shrinivas Moorthi | Remove mo3, mct, midea and define io3, ict etc correctly and get correct cloud condensate.
!> 2013-08-02 | Shrinivas Moorthi | Rewrote the whole routine to read the sigma file differently and to read all tracers. Added sptezj for two 2d fields
!> 2014-02-20 | Shrinivas Moorthi | Modified conversion from spectral to grid taking advantage of threding in SP library. This really speeds up the code. Also threaded loop for Temperature from Tv
!>
!> @author Mark Iredell np23 @date 1999-10-18

! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!> rtsig() reads and transforms sigma file.
!>
!> @param[out] lusig integer(sigio_intkind) sigma file unit number.
!> @param[out] head type(sigio_head) sigma file header.
!> @param[out] k1 integer first model level to return.
!> @param[out] k2 integer last model level to return.
!> @param[out] kgds integer (200) GDS to which to transform.
!> @param[out] ijo integer dimension of output fields.
!> @param[out] levs integer number of total vertical levels.
!> @param[out] ntrac integer number of output tracers.
!> @param[out] jcap integer number of waves.
!> @param[out] lnt2 integer (jcap+1)*(jcap+2).
!> @param[in] me integer
!> @param[out] h real (ijo) surface orography (m).
!> @param[out] p real (ijo) surface pressure (Pa).
!> @param[out] px real (ijo) log surface pressure x-gradient (1/m).
!> @param[out] py real (ijo) log surface pressure y-gradient (1/m).
!> @param[out] t real (ijo,k1:k2) temperature (K).
!> @param[out] u real (ijo,k1:k2) x-component wind (m/s).
!> @param[out] v real (ijo,k1:k2) y-component wind (m/s).
!> @param[out] d real (ijo,k1:k2) wind divergence (1/s).
!> @param[out] trc real (ijo,k1:k2,ntrac) tracers.
!><pre>
!>                                   1 = specific humidity (kg/kg)
!>                                   2 = Ozone mixing ratio (kg/kg)
!>                                   3 = cloud condensate mixing ratio (kg/kg)
!>                                       atomic oxyge, oxygen etc
!>
!></pre>
!> @param[out] iret Integer return code.
!>
!> ### Program History Log
!> Date | Programmer | Comments
!> -----|------------|---------
!> 1999-10-18 | Mark Iredell      | Initial
!> 2013-04-19 | Jun  Wang         | Add option to get tmp and ps(in pascal) from enthalpy and ps(cb) option
!> 2013-05-06 | Shrinivas Moorthi | Initialize midea to 0
!> 2013-05-07 | Shrinivas Moorthi | Remove mo3, mct, midea and define io3, ict etc correctly and get correct cloud condensate.
!> 2013-08-02 | Shrinivas Moorthi | Rewrote the whole routine to read the sigma file differently and to read all tracers. Added sptezj for two 2d fields
!> 2014-02-20 | Shrinivas Moorthi | Modified conversion from spectral to grid taking advantage of threding in SP library. This really speeds up the code. Also threaded loop for Temperature from Tv
!>
!> @author Mark Iredell np23 @date 1999-10-18
subroutine rtsig(lusig,head,k1,k2,kgds,ijo,levs,ntrac,jcap,lnt2,me,     &
                 h,p,px,py,t,u,v,d,trc,iret)

  use sigio_module,   only : sigio_intkind, sigio_head
  use sigio_r_module, only : sigio_dati, sigio_rrdati
  use physcons_post,       only : con_omega, con_fvirt
  use omp_lib
  implicit none
  integer(sigio_intkind),intent(in)    :: lusig
  type(sigio_head),      intent(in)    :: head
  integer,intent(in)                   :: k1,k2,kgds(200),ijo,levs,ntrac,jcap,lnt2,me
  real,dimension(ijo),   intent(out)   :: h,p,px,py
  real,dimension(ijo,k1:k2),intent(out):: t,u,v,d
  real,dimension(ijo,k1:k2,ntrac),intent(out),target :: trc
  integer,intent(out) :: iret
!
  integer idrt,io,jo,iidea
! integer idrt,io,jo,mo3,mct,iidea,midea
  integer(sigio_intkind):: irets
! type(sigio_datm):: datm
  type(sigio_dati) dati
! type griddata
! real,dimension(:,:),pointer :: datm
! endtype griddata
! type(griddata),dimension(:),pointer :: datatrc
  real, target ::  trisca(lnt2,k1:k2+1), triscb(lnt2,k1:k2)
  real,dimension(:),  allocatable :: cpi
  real,dimension(:,:),allocatable :: wrk
  integer io3,ict,jct,n,i,k,jc,nt
  integer idvm, klen
  real    pmean,sumq,xcp
! integer, parameter :: latch=20
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  Determine output grid
  idrt = kgds(1)
  if(kgds(1) == 0 .and. kgds(4) < 90000) idrt = 256
  io = kgds(2)
  jo = kgds(3)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  Read and transform surface fields
  iret = 1
  if (me == 0) then
    print*,'Debug rtsig= ',lusig,k1,k2,ijo,kgds(1:20)
  endif

  idvm  = head%idvm
! jc = omp_get_num_threads()
! write(0,*)' in RTSIG lnt2=',lnt2,' threads=',jc,' latch=',latch,   &
!          ' jcap=',jcap,' io=',io,' jo=',jo,' ijo=',ijo
!
  if (k2 < k1) return

  dati%i = 1                                           ! hs
  dati%f => trisca(:,k1)
  call sigio_rrdati(lusig,head,dati,irets)
  if(irets /= 0) return

! call sptez(0,jcap,idrt,io,jo,trisca(1,k1),h,1)
! call sptez(0,jcap,idrt,io,jo,dats%hs,h,1)
! call sptez(0,jcap,idrt,io,jo,dats%ps,p,1)
! call sptezj(jcap,lnt2,1,idrt,io,jo,jc,trisca,h,latch,1)
!
  dati%i = 2                               ! Surface pressure
  dati%f => trisca(:,k1+1)
  call sigio_rrdati(lusig,head,dati,irets)
  if(irets /= 0) return
!
! call sptez(0,jcap,idrt,io,jo,trisca(1,k1),p,1)
! call sptezj(jcap,lnt2,1,idrt,io,jo,jc,trisca,p,latch,1)
!--
  allocate(wrk(ijo,2))
  call sptezm(0,jcap,idrt,io,jo,2,trisca(1,k1),wrk,1)
  if( mod(idvm,10) < 2) then
!$omp parallel do private(i)
    do i=1,ijo
      h(i) = wrk(i,1)
      p(i) = 1.e3*exp(wrk(i,2))
!     p(i) = 1.e3*exp(p(i))
    enddo
  elseif(mod(idvm,10) == 2) then
!$omp parallel do private(i)
    do i=1,ijo
      h(i) = wrk(i,1)
      p(i) = 1000.*wrk(i,2)
!     p(i) = 1000.*p(i)
    enddo
  endif
  if (allocated(wrk)) deallocate(wrk)

  call sptezd(0,jcap,idrt,io,jo,trisca(1,k1+1),pmean,px,py,1)
  iret = 0

! if (k2 < k1) return

! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  Read and transform fields on levels k1 through k2
  iret = 2
  if (k2 >= k1) then
    klen = k2-k1+1
    do k=k1,k2
      write(0,*)' retriving T for k=',k,' k1=',k1,' k2=',k2
      dati%i = k + 2                        ! Virtual Temperature or CpT
      dati%f => trisca(:,k)

      call sigio_rrdati(lusig,head,dati,iret)
    enddo
    call sptezm(0,jcap,idrt,io,jo,klen,trisca(1,k1),t(1,k1),1)
!   call sptezm(0,jcap,idrt,io,jo,klen,trisca,t,1)
    do k=k1,k2
      dati%i = levs + 2 + (k-1) * 2 + 1     ! Divergence
      dati%f => trisca(:,k)
      call sigio_rrdati(lusig,head,dati,irets)
      if(irets /= 0) return
      dati%i = levs + 2 + (k-1) * 2 + 2     ! Vorticity
      dati%f => triscb(:,k)
      call sigio_rrdati(lusig,head,dati,irets)
      if(irets /= 0) return
    enddo
    call sptezmv(0,jcap,idrt,io,jo,klen,trisca(1,k1),triscb(1,k1),  &
                 u(1,k1),v(1,k1),1)
    call sptezm(0,jcap,idrt,io,jo,klen,trisca(1,k1),d(1,k1),1)

!   call sptezm(0,jcap,idrt,io,jo,1,triscb,z(1,k),1)
    write(0,*)' retriving d/z for k=',k,' k1=',k1,' k2=',k2
!   datm%z(3,:) = datm%z(3,:)+2*con_omega/sqrt(1.5)
!   call sptezm(0,jcap,idrt,io,jo,klen,datm%z,z,1)
   write(0,*)' start get tracer'
    do nt=1,ntrac
      do k=k1,k2
        dati%i = levs * (2+nt) + 2 + k      ! Tracers starting with q
        dati%f => trisca(:,k)
        call sigio_rrdati(lusig,head,dati,irets)
      enddo
      call sptezm(0,jcap,idrt,io,jo,klen,trisca(1,k1),trc(1,k1,nt),1)
      write(0,*)' retriving d/z for nt=',nt,'ntrac=',ntrac,'k=',k,' k1=',k1,' k2=',k2
    enddo
    !t=t/(1+con_fvirt*sh)
   write(0,*)' end get tracer,idvm=',idvm,'ijo=',ijo,'ntrac=',ntrac
!
!-- get temp 
    if (mod(idvm/10,10) == 3) then ! Enthalpy case
      allocate(cpi(0:ntrac))
!     write(0,*)'aft read sig, cpi=',head%cpi
      cpi(0:ntrac) = head%cpi(1:ntrac+1)
!     write(0,*)'cpi=',cpi(0:ntrac)
!$omp parallel do private(k,i,xcp,sumq,n)
      do k=k1,k2
        do i=1,ijo
          xcp  = 0.0
          sumq = 0.0
          do n=1,ntrac
            if( cpi(n) /= 0.0 ) then
              xcp  = xcp  + cpi(n)*trc(i,k,n)
              sumq = sumq + trc(i,k,n)
            endif
          enddo
          xcp    = (1.-sumq)*cpi(0) + xcp
          t(i,k) = t(i,k) / xcp   ! Now g1 contains T
        enddo
      enddo
      if (allocated(cpi)) deallocate(cpi)
    else
!$omp parallel do private(i,k)
      do k=k1,k2
        do i=1,ijo
          t(i,k) = t(i,k) / (1+con_fvirt*trc(i,k,1)) !get temp from virtual temp
        enddo
      enddo
    endif
  endif
! write(0,*)'end comput t'
  iret=0
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
end subroutine  

! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!> modstuff() computes model coordinate dependent functions.
!>
!> This subprogram computes fields which depend on the model coordinate
!> such as pressure thickness and vertical velocity.
!>
!> @param[in] km integer number of levels.
!> @param[in] idvc integer vertical coordinate id (1 for sigma and 2 for hybrid).
!> @param[in] idsl integer type of sigma structure (1 for phillips or 2 for mean).
!> @param[in] nvcoord integer number of vertical coordinates.
!> @param[in] vcoord real (km+1,nvcoord) vertical coordinates.
!> @param[in] ps real surface pressure (Pa).
!> @param[in] psx real log surface pressure x-gradient (1/m).
!> @param[in] psy real log surface pressure y-gradient (1/m).
!> @param[in] d real (km) wind divergence (1/s).
!> @param[in] u real (km) x-component wind (m/s).
!> @param[in] v real (km) y-component wind (m/s).
!> @param[out] pi real (km+1) interface pressure (Pa).
!> @param[out] pm real (km) mid-layer pressure (Pa).
!> @param[out] om real (km) vertical velocity (Pa/s).
!>
!> ### Program History Log
!> Date | Programmer | Comments
!> -----|------------|---------
!> 1999-10-18 | Mark Iredell | Initial
!> 2013-04-19 | Jun Wang     | Add option to get pi by using 8 byte real computation
!>
!> @author Mark Iredell np23 @date 1999-10-18
  subroutine modstuff(km,idvc,idsl,nvcoord,vcoord,ps,psx,psy,d,u,v,&
                      pi,pm,om)
    use sigio_module, only: sigio_modprd
    implicit none
    integer,intent(in):: km,idvc,idsl,nvcoord
    real,intent(in):: vcoord(km+1,nvcoord)
    real,intent(in):: ps,psx,psy
    real,intent(in):: u(km),v(km),d(km)
!   real,intent(out):: pi(km+1),pm(km)
    real*8, intent(out):: pi(km+1),pm(km)
    real,intent(out):: om(km)
    real*8 ps8,pm8(km),pd8(km),vcoord8(km+1,nvcoord)
    real*8 dpmdps(km),dpddps(km),dpidps(km+1),pi8(km+1)
    real vgradp,pd(km),px(km),py(km),os
    integer k,iret,logk
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    ps8=ps
    vcoord8=vcoord
    call sigio_modprd(1,1,km,nvcoord,idvc,idsl,vcoord8,iret,&
                     ps=(/ps8/),&
                     pm=pm8,pd=pd8,dpmdps=dpmdps,dpddps=dpddps)
!
!jw: has to be 8 real for wam     
    pi8(1)=ps
    pm=pm8
!   pd=pd8
    dpidps(1)=1.
    do k=1,km
      pi8(k+1)=pi8(k)-pd8(k)
      dpidps(k+1)=dpidps(k)-dpddps(k)
!     if(pi(8)<0.) then
!        print *,'in modstuff,pi8=',pi8(k)
!     endif
    enddo
    pi=pi8
!
    os=0
    do k=km,1,-1
      vgradp=u(k)*psx+v(k)*psy
      os=os-vgradp*ps*(dpmdps(k)-dpidps(k+1))-d(k)*(pm(k)-pi(k+1))
      om(k)=vgradp*ps*dpmdps(k)+os
      os=os-vgradp*ps*(dpidps(k)-dpmdps(k))-d(k)*(pi(k)-pm(k))
    enddo
    px=ps*dpmdps*psx
    py=ps*dpmdps*psy
  end subroutine

!-------------------------------------------------------------------------------
!> modstuff2() computes model coordinate dependent functions.
!>
!> This subprogram computes fields which depend on the model coordinate
!> such as pressure thickness and vertical velocity.
!>
!> @param[in] im integer inner computational domain.
!> @param[in] ix integer maximum inner dimension.
!> @param[in] km integer number of levels.
!> @param[in] idvc integer vertical coordinate id (1 for sigma and 2 for hybrid).
!> @param[in] idsl integer type of sigma structure (1 for phillips or 2 for mean).
!> @param[in] nvcoord integer number of vertical coordinates.
!> @param[in] vcoord real (km+1,nvcoord) vertical coordinates.
!> @param[in] ps real surface pressure (Pa).
!> @param[in] psx real log surface pressure x-gradient (1/m).
!> @param[in] psy real log surface pressure y-gradient (1/m).
!> @param[in] d real (km) wind divergence (1/s).
!> @param[in] u real (km) x-component wind (m/s).
!> @param[in] v real (km) y-component wind (m/s).
!> @param[out] pi real (km+1) interface pressure (Pa).
!> @param[out] pm real (km) mid-layer pressure (Pa).
!> @param[out] om real (km) vertical velocity (Pa/s).
!> @param[in] me integer
!>
!> ### Program History Log
!> Date | Programmer | Comments
!> -----|------------|---------
!> 1999-10-18 | Mark Iredell      | Initial
!> 2013-04-19 | Jun  Wang         | Add option to get pi by using 8 byte real computation
!> 2013-08-13 | Shrinivas Moorthi | Modified to include im points and thread
!>
!> @author Mark Iredell np23 @date 1999-10-18
  subroutine modstuff2(im,ix,km,idvc,idsl,nvcoord,vcoord,ps,psx,psy,d,u,v,&
                       pi,pm,om,me)
    use sigio_module, only : sigio_modprd
    implicit none
    integer,                    intent(in)  :: im,ix,km,idvc,idsl,nvcoord,me
    real,                       intent(in)  :: vcoord(km+1,nvcoord)
    real,   dimension(ix),      intent(in)  :: ps,psx,psy
    real,   dimension(ix,km),   intent(in)  :: u,v,d
    real*8, dimension(ix,km+1), intent(out) :: pi
    real*8, dimension(ix,km),   intent(out) :: pm
    real,   dimension(ix,km),   intent(out) :: om
!   real*8, allocatable :: ps8(:), pm8(:,:), pd8(:,:),dpmdps(:,:),dpddps(:,:), &
!                          dpidps(:,:),pi8(:,:),vcoord8(:,:)
!   real,   allocatable :: os(:)
!   real,   allocatable :: pd(:,:),px(:,:), py(:,:), os(:)

!   real vgradpps

    real*8 ps8(ix),pm8(ix,km),pd8(ix,km),vcoord8(km+1,nvcoord)
    real*8 dpmdps(ix,km),dpddps(ix,km),dpidps(ix,km+1),pi8(ix,km+1)
    real vgradpps,pd(im,km),os(im)
!   real vgradpps,pd(im,km),px(im,km),py(im,km),os(im),tem
    integer i,k,iret,logk
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    ps8     = ps
    vcoord8 = vcoord
    call sigio_modprd(im,ix,km,nvcoord,idvc,idsl,vcoord8,iret,                 &
                     ps=ps8,pd=pd8,dpddps=dpddps,pm=pm8,dpmdps=dpmdps)

!
!   if (me == 0) then
!     write(0,*)' pd8=',pd8(1,60:64)
!     write(0,*)' pm8=',pm8(1,60:64)
!    endif
!jw: has to be 8 real for wam     

!$omp parallel do private(i)
    do i=1,im
      pi8(i,1)    = ps(i)
      dpidps(i,1) = 1.
      os(i)       = 0
      pi(i,1)     = pi8(i,1)
    enddo
    do k=1,km
!$omp parallel do private(i)
      do i=1,im
        pi8(i,k+1)    = pi8(i,k)    - pd8(i,k)
        dpidps(i,k+1) = dpidps(i,k) - dpddps(i,k)
!       if(pi(i,8)<0.) then
!          print *,'in modstuff,pi8=',pi8(i,k),' i=',i,' k=',k,' me=',me
!       endif
        pi(i,k+1) = pi8(i,k+1)
        pm(i,k)   = pm8(i,k)
      enddo
    enddo
!
    do k=km,1,-1
!$omp parallel do private(i,vgradpps)
      do i=1,im
        vgradpps = (u(i,k)*psx(i) + v(i,k)*psy(i)) * ps(i)

        os(i)    = os(i) - vgradpps*(dpmdps(i,k)-dpidps(i,k+1))             &
                         - d(i,k)*(pm(i,k)-pi(i,k+1))

        om(i,k)  = os(i) + vgradpps*dpmdps(i,k)

        os(i)    = os(i) - vgradpps*(dpidps(i,k)-dpmdps(i,k))               &
                         - d(i,k)*(pi(i,k)-pm(i,k))
      enddo
    enddo
  end subroutine

!-----------------------------------------------------------------------
!> trssc() transforms sigma spectral fields to grid.
!>
!> Transforms sigma spectral fields to grid and converts
!> log surface pressure to surface pressure and virtual temperature
!> to temperature.
!>
!> @param[in] jcap integer spectral truncation.
!> @param[in] nc integer first dimension (nc>=(jcap+1)*(jcap+2)).
!> @param[in] km integer number of levels.
!> @param[in] ntrac integer number of tracers.
!> @param[in] idvc integer 
!> @param[in] idvm integer mass variable id.
!> @param[in] idsl integer
!> @param[in] nvcoord integer number of vertical coordinates.
!> @param[in] vcoord real (km+1,nvcoord) vertical coordinates.
!> @param[in] cpi real
!> @param[in] idrt integer data representation type.
!> @param[in] lonb integer number of longitudes.
!> @param[in] latb integer number of latitudes.
!> @param[in] ijl integer horizontal dimension.
!> @param[in] ijn integer
!> @param[in] j1 integer first latitude.
!> @param[in] j2 integer last latitude.
!> @param[in] jc integer number of cpus.
!> @param[in] chgq0 integer
!> @param[in] szs real (nc) orography.
!> @param[in] sps real (nc) log surface pressure.
!> @param[in] st real (nc,levs) virtual temperature.
!> @param[in] sd real (nc,levs) divergence.
!> @param[in] sz real (nc,levs) vorticity.
!> @param[in] sq real (nc,levs*ntrac) tracers.
!> @param[inout] gfszs real (ijn)
!> @param[inout] gfsps real (ijn)
!> @param[inout] gfsp real (ijn, km)
!> @param[inout] gfsdp real (ijn, km)
!> @param[inout] gfst real (ijn, km)
!> @param[inout] gfsu real (ijn, km)
!> @param[inout] gfsv real (ijn, km)
!> @param[inout] gfsq real (ijn, km*ntrac)
!> @param[inout] gfsw real (ijn, km)
!>
!> ### Program History Log
!> Date | Programmer | Comments
!> -----|------------|---------
!> 1999-10-18 | Mark Iredell | Initial
!>
!> @author Mark Iredell w/nmc23 @date 1992-10-31
      subroutine trssc(jcap,nc,km,ntrac,idvc,idvm,idsl,nvcoord,vcoord,     &
                       cpi,idrt,lonb,latb,ijl,ijn,j1,j2,jc,chgq0,          &
                       szs,sps,st,sd,sz,sq,gfszs,gfsps,gfsp,gfsdp,         &
                       gfst,gfsu,gfsv,gfsq,gfsw)
      implicit none
      integer,intent(in)::jcap,nc,km,ntrac,idvc,idvm,idsl,nvcoord,idrt,lonb,latb
      integer,intent(in)::ijl,ijn,j1,j2,jc,chgq0
      real,intent(in):: szs(nc),sps(nc),st(nc,km),sd(nc,km),sz(nc,km),sq(nc,km*ntrac)
      real,intent(in):: cpi(0:ntrac)
      real*8,intent(in):: vcoord(km+1,nvcoord)
      real,dimension(ijn),intent(inout):: gfszs,gfsps
      real,dimension(ijn,km),intent(inout):: gfsp,gfsdp,gfst,gfsu,gfsv,gfsw
      real,dimension(ijn,km*ntrac),intent(inout):: gfsq
      real zs(ijl),ps(ijl),t(ijl,km),u(ijl,km),v(ijl,km),q(ijl,km*ntrac)
      real wi(ijl,km),pi(ijl,km),dpo(ijl,km)
      real tvcon,xcp,sumq
      integer thermodyn_id,jn,js,is,in
      integer jj,jjm,k,n,j,i,ij,lonb2
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  spectral transforms
      if(j1==732)print*,'sample input to trssc= ',jcap,nc,km,ntrac, &
                          idvc,idvm,idsl,nvcoord,    &
                       idrt,lonb,latb,ijl,ijn,j1,j2,jc,chgq0
      lonb2=lonb*2
      ij=lonb2*(j2-j1+1)
      in=1
      is=1+lonb
      call sptran(0,jcap,idrt,lonb,latb,1,1,1,lonb2,lonb2,nc,ijl,           &
                  j1,j2,jc,szs,zs(in),zs(is),1)
      call sptran(0,jcap,idrt,lonb,latb,1,1,1,lonb2,lonb2,nc,ijl,           &
                  j1,j2,jc,sps,ps(in),ps(is),1)
      call sptran(0,jcap,idrt,lonb,latb,km,1,1,lonb2,lonb2,nc,ijl,          &
                  j1,j2,jc,st,t(in,1),t(is,1),1)
      call sptranv(0,jcap,idrt,lonb,latb,km,1,1,lonb2,lonb2,nc,ijl,         &
                   j1,j2,jc,sd,sz,u(in,1),u(is,1),v(in,1),v(is,1),1)
      call sptran(0,jcap,idrt,lonb,latb,km*ntrac,1,1,lonb2,lonb2,nc,ijl,    &
                  j1,j2,jc,sq,q(in,1),q(is,1),1)
      if(j1==732)then
       do k=1,km
        do i=1,ijl
	 if(t(i,k)>400. .or. t(i,k)<100.)print*,'bad T from sptran',i,k,t(i,k)
	 if(q(i,k)>1.)print*,'bad Q  from sptran',i,k,q(i,k)
	 if(q(i,2*k)>1.)print*,'bad Q  from sptran',i,k,q(i,2*k)
	 if(q(i,3*k)>1.)print*,'bad Q  from sptran',i,k,q(i,3*k)
	end do
       end do	 
      end if  
      select case(mod(idvm,10))
      case(0,1)
        do i=1,ij
          ps(i)=1.e3*exp(ps(i))
        enddo
      case(2)
        do i=1,ij
          ps(i)=1.e3*ps(i)
        enddo
      case default
        do i=1,ij
          ps(i)=1.e3*exp(ps(i))
        enddo
      end select
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      thermodyn_id=mod(idvm/10,10)
      if (thermodyn_id == 3) then
        do k=1,km
         do i=1,ij
            t(i,k) = t(i,k)/cpi(0)   ! enthalpy (cpt/cpd)
         end do
        end do
!
      endif

      call getomega(jcap,nc,km,idvc,idvm,idrt,idsl,                    &
        nvcoord,vcoord,lonb,latb,ijl,j1,j2,1,sd,sps,                   &
        ps,t,u,v,wi,pi,dpo)
      if(j1==732)then
       do k=1,km
        do i=1,ijl
	 if(t(i,k)>400. .or. t(i,k)<100.)print*,'bad T after getomega',i,k,t(i,k)
	 if(q(i,k)>1. )print*,'bad Q  after getomega',i,k,q(i,k)
	 if(q(i,2*k)>1. )print*,'bad Q  after getomega',i,2*k,q(i,2*k)
	end do
       end do	 
      end if  
      if(thermodyn_id /= 2)then
!  convert to surface pressure and temperature
         if (thermodyn_id == 3) then
           do k=1,km
            do i=1,ij
               xcp  = 0.0
               sumq = 0.0
               do n=1,ntrac
                 if( cpi(n) /= 0.0 ) then
                   xcp  = xcp  + cpi(n)*q(i,k+(n-1)*km)
                   sumq = sumq + q(i,k+(n-1)*km)
                 endif
               enddo
               t(i,k)  = t(i,k)/((1.-sumq)*cpi(0)+xcp)
            end do
           end do

          else 
           tvcon=(461.50/287.05-1.)
           t(:,:) = t(:,:)/(1.+tvcon*q(:,1:km))
          endif
      end if
      if(j1==732)then
       do k=1,km
        do i=1,ijl
	 if(t(i,k)>400. .or. t(i,k)<100.)print*,'bad T after Tv to T',i,k,t(i,k)
	 if(q(i,k)>1.)print*,'bad Q  after Tv to T',i,k,q(i,k)
	 if(q(i,2*k)>1. )print*,'bad Q  after Tv to T',i,k,q(i,2*k)
	end do
       end do	 
      end if  
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!----force tracers to be positive
      if (chgq0 == 1) q = max(q, 0.0)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  pass data to gfsdatao
      do j=1,j2-j1+1
        jn=j+j1-1
        js=latb+1-jn
        jn=(jn-1)*lonb
        js=(js-1)*lonb
        jj=j*lonb
        jjm=(j-1)*lonb
        do i=1,lonb
          gfszs(i+jn) = zs(i+jjm)
          gfsps(i+jn) = ps(i+jjm)
        enddo
        do i=1,lonb
          gfszs(i+js) = zs(i+jj)
          gfsps(i+js) = ps(i+jj)
        enddo
        do k=1,km
         do i=1,lonb
          gfsdp(i+jn,k) = dpo(i+jjm,k)
          gfsp(i+jn,k)  = pi(i+jjm,k)
          gfst(i+jn,k)  = t(i+jjm,k)
          gfsu(i+jn,k)  = u(i+jjm,k)
          gfsv(i+jn,k)  = v(i+jjm,k)
          gfsw(i+jn,k)  = wi(i+jjm,k)
         enddo
         do i=1,lonb
          gfsdp(i+js,k) = dpo(i+jj,k)
          gfsp(i+js,k)  = pi(i+jj,k)
          gfst(i+js,k)  = t(i+jj,k)
          gfsu(i+js,k)  = u(i+jj,k)
          gfsv(i+js,k)  = v(i+jj,k)
          gfsw(i+js,k)  = wi(i+jj,k)
         enddo
        enddo
        do k=1,km*ntrac
          do i=1,lonb
            gfsq(i+jn,k) = q(i+jjm,k)
          enddo
          do i=1,lonb
            gfsq(i+js,k) = q(i+jj,k)
          enddo
        enddo
      enddo
      return
      end
!-----------------------------------------------------------------------
!> getomega() computes omega.
!>
!> @param[in] jcap integer spectral truncation.
!> @param[in] nc integer first dimension (nc>=(jcap+1)*(jcap+2)).
!> @param[in] km integer number of levels.
!> @param[in] idvc integer
!> @param[in] idvm integer mass variable id.
!> @param[in] idrt integer data representation type.
!> @param[in] idsl integer
!> @param[in] nvcoord integer number of vertical coordinates.
!> @param[in] vcoord real (km+1,nvcoord) vertical coordinates.
!> @param[in] lonb integer number of longitudes.
!> @param[in] latb integer number of latitudes.
!> @param[in] ijn integer
!> @param[in] j1 integer first latitude.
!> @param[in] j2 integer last latitude.
!> @param[in] jc integer number of cpus.
!> @param[in] sd real (nc,km) divergence.
!> @param[in] sps real (nc) log surface pressure.
!> @param[in] psi real (ijn)
!> @param[in] ti real (ijn,km)
!> @param[in] ui real (ijn,km)
!> @param[in] vi real (ijn,km)
!> @param[out] wi real (ijn,km)
!> @param[out] pm real (ijn,km)
!> @param[out] pd real (ijn,km)
!>
      subroutine getomega(jcap,nc,km,idvc,idvm,idrt,idsl,nvcoord,vcoord,  &
            lonb,latb,ijn,j1,j2,jc,sd,sps,psi,ti,ui,vi,wi,pm,pd)
!!!!!
      use sigio_module, only : sigio_modprd
      implicit none

      integer,intent(in):: jcap,nc,km,idvc,idvm,idrt,idsl,nvcoord
      integer,intent(in):: lonb,latb,j1,j2,jc,ijn
      real*8,intent(in):: vcoord(km+1,nvcoord)
      real,intent(in):: sd(nc,km),sps(nc)
      real,intent(in):: psi(ijn),ti(ijn,km),ui(ijn,km),vi(ijn,km)
      real,intent(out):: wi(ijn,km),pm(ijn,km),pd(ijn,km)
      real :: pi(ijn,km+1)
      real :: os
      real*8 psi8(ijn),ti8(ijn,km),pm8(ijn,km),pd8(ijn,km)
      real*8 dpmdps(ijn,km),dpddps(ijn,km),dpidps(ijn,km+1),vgradp,psmean
      real di(ijn,km),psx(ijn),psy(ijn)
      integer k,i,ij,lonb2,iret,is,in
!----1. spectral transform
      lonb2=lonb*2
      ij=lonb2*(j2-j1+1)
      in=1
      is=1+lonb
      call sptrand(0,jcap,idrt,lonb,latb,1,1,1,lonb2,lonb2,nc,ijn,    &
           j1,j2,jc,sps,psmean,psx(in),psx(is),psy(in),psy(is),1)

      call sptran(0,jcap,idrt,lonb,latb,km,1,1,lonb2,lonb2,nc,ijn,    &
                  j1,j2,jc,sd,di(in,1),di(is,1),1)
      psi8=psi	
      ti8=ti  

      call sigio_modprd(ijn,ijn,km,nvcoord,idvc,idsl,vcoord,iret,     &
                   ps=psi8,t=ti8,pm=pm8,pd=pd8,dpmdps=dpmdps,dpddps=dpddps)
      pm=pm8
      pd=pd8		   

      select case(mod(idvm,10))
      case(0,1)
          continue
      case(2)
          do i=1,ijn
           psx(i)=psx(i)/(psi(i)*1.0e-3)
           psy(i)=psy(i)/(psi(i)*1.0e-3)
          enddo
      case default
          do i=1,ijn
           psx(i)=psx(i)/psi(i)
           psy(i)=psy(i)/psi(i)
          enddo
      end select

!----3.omeda from modstuff
      do i=1,ijn
       pi(i,1)=psi(i)
       dpidps(i,1)=1.
      enddo
      do k=1,km
       do i=1,ijn
         pi(i,k+1)=pi(i,k)-pd(i,k)
         dpidps(i,k+1)=dpidps(i,k)-dpddps(i,k)
       enddo
      enddo
      do i=1,ijn
       os=0.
       do k=km,1,-1
        vgradp=ui(i,k)*psx(i)+vi(i,k)*psy(i)
        os=os-vgradp*psi(i)*(dpmdps(i,k)-dpidps(i,k+1))-                 &
           di(i,k)*(pm(i,k)-pi(i,k+1))
        wi(i,k)=vgradp*psi(i)*dpmdps(i,k)+os
        os=os-vgradp*psi(i)*(dpidps(i,k)-dpmdps(i,k))-                   &
           di(i,k)*(pi(i,k)-pm(i,k))
       enddo
!
      enddo
!---
       return
       end subroutine