   
   module module_trajectory

   use module_driver_constants,  only : max_domains
   use module_state_description, only : num_chem

   implicit none

   private
   public :: trajectory_init
#ifdef NETCDF
   public :: trajectory_driver
   public :: trajectory_dchm_tstep_init
   public :: trajectory_dchm_tstep_set
#endif
   public :: traject
   public :: traj_cnt

   integer, parameter :: vals_max = 1000
   integer, parameter :: traj_max = 1000
   integer, parameter :: var_max  = 100
   integer, parameter :: pkg_max  = 6
   integer, parameter :: dyn_max  = 8
   integer, parameter :: chm_pkg  = 1
   integer, parameter :: hyd_pkg  = 2
   integer, parameter :: trc_pkg  = 3
   integer, parameter :: dyn_pkg  = 4
   integer, parameter :: msc_pkg  = 5
   integer, parameter :: dchm_pkg = 6
   real, parameter    :: missing_val = -9999.
   real, parameter    :: zero_val    = 0.
   
   integer :: n_dom        ! domain count
   integer, pointer :: n_vals                         ! wrking count of time points in buffer
   integer, target  :: n_vals_dm(max_domains)         ! count of time points in buffer per domain
   integer, target  :: n_dchm_dm(max_domains)         ! count of dchm variables in  buffer per domain
   integer, allocatable, target :: num_msc_dm(:)      ! number misc species
   integer, allocatable, target :: dchm_buf_ndx_dm(:,:)  ! dchm buffer indices
   integer, pointer :: dchm_buf_ndx(:)
   integer :: offset       ! into variable arrays

   type default
     character(len=19) :: start_time
     character(len=19) :: stop_time
     character(len=32) :: chm_name(var_max)
     character(len=32) :: hyd_name(var_max)
     character(len=32) :: trc_name(var_max)
     character(len=32) :: dyn_name(var_max)
     character(len=32) :: msc_name(var_max)
     character(len=32) :: dchm_name(var_max)
   end type default

   type base
     real    :: lon
     real    :: lat
     real    :: lev
     real    :: x, y
     real    :: traj_var(var_max)
     integer :: n_chm_var
     integer :: n_ct_var
     integer :: n_hyd_var
     integer :: n_trc_var
     integer :: n_dyn_var
     integer :: n_msc_var
     integer :: n_dchm_var
     integer :: chm_ndx(var_max)
     integer :: hyd_ndx(var_max)
     integer :: trc_ndx(var_max)
     integer :: dyn_ndx(var_max)
     integer :: msc_ndx(var_max)
     integer :: dchm_ndx(var_max)
     character(len=19) :: start_time
     character(len=19) :: stop_time
     character(len=32) :: chm_spc(var_max)
     character(len=32) :: hyd_spc(var_max)
     character(len=32) :: trc_spc(var_max)
     character(len=32) :: dyn_var(var_max)
     character(len=32) :: msc_var(var_max)
     character(len=32) :: dchm_spc(var_max)
     logical :: in_dom
     logical :: in_patch
     logical :: is_stationary
     logical :: chm_is_gas(var_max)
   end type base

   type buffer
     real    :: trj_i(vals_max)
     real    :: trj_j(vals_max)
     real    :: trj_k(vals_max)
     real    :: trj_lons(vals_max)
     real    :: trj_lats(vals_max)
     real, allocatable :: chm_vals(:,:)
     real, allocatable :: trc_vals(:,:)
     real, allocatable :: hyd_vals(:,:)
     real, allocatable :: dyn_vals(:,:)
     real, allocatable :: msc_vals(:,:)
     real, allocatable :: dchm_vals(:,:)
     character(len=19) :: times(vals_max)
   end type buffer

   type statevar
     character(len=80) :: Varname
     character(len=80) :: Description
     character(len=80) :: Units
     character(len=10) :: MemOrd
     character(len=10) :: Stagger
     integer           :: Ndim
     real, pointer     :: rfield_2d(:,:)
     real, pointer     :: rfield_3d(:,:,:)
   end type statevar

   integer,    allocatable           :: traj_cnt(:)
   type(base), allocatable, target   :: traject(:,:)
   type(base), pointer               :: trjects(:)
   type(buffer), allocatable, target :: trj_buff(:,:)
   type(buffer), pointer             :: trj_pbf(:)
   type(statevar), allocatable, target :: St_Vars_dm(:,:)
   type(statevar), pointer             :: St_Vars(:)

   real, allocatable  :: dchm_buff(:,:,:,:)

   character(len=256) :: dyn_var_lst(dyn_max)
   character(len=32)  :: dyn_var_desc_att(dyn_max)
   character(len=32)  :: dyn_var_unit_att(dyn_max)
   character(len=4)   :: pkg_tag(pkg_max) = (/ 'chm ', 'hyd ', 'trc ', 'dyn ', 'msc ', 'dchm' /)

   logical                      :: do_chemstep
   logical, allocatable, target :: pkg_has_vars_dm(:,:,:)
   logical, pointer             :: pkg_has_vars(:,:)
   logical, allocatable, target :: trj_msk_dm(:,:,:,:)
   logical, pointer             :: trj_msk(:,:)
   logical, allocatable, target :: St_Vars_msk_dm(:,:)
   logical, pointer             :: St_Vars_msk(:)
   logical, allocatable, target :: dchm_msk_dm(:,:)
   logical, pointer             :: dchm_msk(:)
   logical, allocatable         :: is_initialized(:)

   CONTAINS

   subroutine trajectory_init( grid, config_flags, &
                               ims,ime, jms,jme, kms,kme )

   use module_domain
   use module_llxy,              only : proj_info, latlon_to_ij
   use module_configure,         only : grid_config_rec_type
   use module_state_description, only : no_trajectory, param_first_scalar, num_chem, num_moist, num_tracer
   use module_scalar_tables,     only : chem_dname_table, moist_dname_table, tracer_dname_table
   use module_model_constants,   only : g
   use module_domain_type,       only : fieldlist
#ifdef DM_PARALLEL
   use module_dm,                only : wrf_dm_sum_integer, wrf_dm_max_real
#endif

!-----------------------------------------------------------------------------
!  dummy arguments
!-----------------------------------------------------------------------------
   integer, intent(in)      :: ims,ime, jms,jme, kms,kme
   type(domain), intent(inout)            :: grid
   type(grid_config_rec_type), intent(in) :: config_flags

!-----------------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------------
   integer :: astat
   integer :: dm
   integer, pointer :: num_msc      ! number misc species
   integer, pointer :: n_dchm       ! number dchm buffer species
   integer :: ierr, ios
   integer :: i, j, k, k1, m, m1, m2
   integer :: i_end, j_end, u_lim
   integer :: n, pkg, trj
   integer :: n_traj, n_traj_1
   integer :: p_size
   integer :: unitno
   integer :: ids,ide, jds,jde, kds,kde
   integer :: ips,ipe, jps,jpe, kps,kpe
   integer :: dbg_lvl
   integer :: target
   integer :: n_def_var(pkg_max)
   real    :: x, y
   real    :: z_dm_bot, z_dm_top
   real    :: z(kms:kme-1)
   real    :: z_at_w(ims:ime,kms:kme,jms:jme)
   character(len=256), allocatable :: msc_tbl(:)
   character(len=32)  :: var_name
   character(len=128) :: filename
   character(len=256) :: err_mes
   character(len=32)  :: wrk_var_name
   character(len=32)  :: wrk_chr(var_max)
   character(len=32)  :: wrk_def_name(var_max)
   logical :: exists
   logical :: is_root_proc
   logical :: rstrt
   logical :: mask(var_max)
   logical :: trj_mask(traj_max)

   type(grid_config_rec_type) :: config_temp
   type(proj_info)            :: proj

   type(default) :: traj_def
   type(base)    :: traj_type(traj_max)

   logical, external :: wrf_dm_on_monitor
   integer, external :: get_unused_unit

   namelist / traj_spec /    traj_type
   namelist / traj_default / traj_def

#ifndef NETCDF
   call wrf_error_fatal( 'trajectory_init: requires netcdf' )
#endif

   offset = param_first_scalar - 1
!-----------------------------------------------------------------------------
!  check for trajectory option
!-----------------------------------------------------------------------------
   if( grid%traj_opt == no_trajectory ) then
     write(err_mes,'(''trajectory_init('',i2.2,''): traj_opt = no_trajectory'')') grid%id
     call wrf_message( trim(err_mes) )
     return
   endif
!-----------------------------------------------------------------------------
!  set domain count, restarting?
!-----------------------------------------------------------------------------
   dm = grid%id
   call nl_get_restart( dm,rstrt )
   if( dm == 1 .and. .not. allocated( traject ) ) then
     call nl_get_max_dom( 1,n_dom )
!-----------------------------------------------------------------------------
!  allocate module variables
!-----------------------------------------------------------------------------
     allocate( traject(traj_max,n_dom),traj_cnt(n_dom), &
               pkg_has_vars_dm(traj_max,pkg_max,n_dom), &
               num_msc_dm(n_dom), is_initialized(n_dom),stat=astat )
     if( astat /= 0 ) then
       write(err_mes,'(''trajectory_init('',i2.2,''): failed to allocate traject...num_msc_dm; error = '',i6)') dm,astat
       call wrf_error_fatal( trim( err_mes  ) )
     endif
     is_initialized(:) = .false.
   endif

   if( is_initialized(dm) ) then
     return
   endif

   trjects => traject(:,dm)

   call get_ijk_from_grid( grid ,                   &
                           ids, ide, jds, jde, kds, kde, &
                           n,n,  n,n, n,n,          &
                           ips, ipe, jps, jpe, kps, kpe    )

   n_vals      => n_vals_dm(dm)
   n_vals      = 0
   traj_cnt(dm) = 0
   is_root_proc = wrf_dm_on_monitor()
   pkg_has_vars => pkg_has_vars_dm(:,:,dm)

   dyn_var_lst(1:dyn_max)      = (/ 'p       ', 'T', 'z', 'u', 'v', 'w', 'rainprod', 'evapprod' /)
   dyn_var_unit_att(1:dyn_max) = (/ 'hPa', 'K', 'm', 'm/s', 'm/s', 'm/s', 's-1', 's-1' /)
   dyn_var_desc_att(1:dyn_max) = (/ 'pressure            ', 'temperature', 'height', &
                                    'x wind component', 'y wind component', 'z wind component', &
                                    'rain production rate', 'rain evap rate' /)
!-----------------------------------------------------------------------------
!  master proc
!-----------------------------------------------------------------------------
!master_proc: &
!   if( is_root_proc ) then
      write(filename,'(''wrfinput_traj_d'',i2.2)',iostat=ios) dm
      if( ios /= 0 ) then
        write(err_mes,'(''trajectory_init('',i2.2,''): failed to set filename: error = '',i6)') dm,ios
        call wrf_error_fatal( trim( err_mes  ) )
      endif
      inquire( file=trim(filename),exist=exists )
input_file: &
      if( exists ) then
        unitno = get_unused_unit()
        if( unitno <= 0 ) then
          call wrf_error_fatal( 'trajectory_init: failed to get unit number' )
        endif
!-----------------------------------------------------------------------------
!  open file
!-----------------------------------------------------------------------------
        open( unit = unitno,file=trim(filename),iostat=ios )
        if( ios /= 0 ) then
          write(err_mes,'(''trajectory_init('',i2.2,''): failed to open '',a,''; error = '',i6)') dm,trim(filename),ios
          call wrf_error_fatal( trim( err_mes  ) )
        endif
!-----------------------------------------------------------------------------
!  initialize trajectories
!-----------------------------------------------------------------------------
        traj_def%start_time = ' '
        traj_def%stop_time  = ' '
        traj_def%chm_name(:) = ' '
        traj_def%hyd_name(:) = ' '
        traj_def%trc_name(:) = ' '
        traj_def%dyn_name(:) = ' '
        traj_def%msc_name(:) = ' '
        traj_def%dchm_name(:) = ' '

        do trj = 1,traj_max
          traj_type(trj)%start_time = ' '
          traj_type(trj)%stop_time  = ' '
          traj_type(trj)%chm_spc(:) = ' '
          traj_type(trj)%dchm_spc(:) = ' '
          traj_type(trj)%hyd_spc(:) = ' '
          traj_type(trj)%trc_spc(:) = ' '
          traj_type(trj)%dyn_var(:) = ' '
          traj_type(trj)%msc_var(:) = ' '
          traj_type(trj)%chm_ndx(:) = 0
          traj_type(trj)%hyd_ndx(:) = 0
          traj_type(trj)%trc_ndx(:) = 0
          traj_type(trj)%dyn_ndx(:) = 0
          traj_type(trj)%msc_ndx(:) = 0
          traj_type(trj)%dchm_ndx(:) = 0
          traj_type(trj)%n_chm_var = 0 ; traj_type(trj)%n_ct_var = 0
          traj_type(trj)%n_hyd_var = 0 ; traj_type(trj)%n_trc_var = 0 
          traj_type(trj)%n_dyn_var = 0 ; traj_type(trj)%n_msc_var = 0
          traj_type(trj)%n_dchm_var = 0
          traj_type(trj)%is_stationary = .false.
          traj_type(trj)%chm_is_gas(:) = .true.
          traj_type(trj)%lon           = missing_val
          traj_type(trj)%lat           = missing_val
          traj_type(trj)%lev           = missing_val
        end do
!-----------------------------------------------------------------------------
!  read file
!-----------------------------------------------------------------------------
        read(unit=unitno,nml=traj_default,iostat=ios)
        if( ios /= 0 ) then
          close( unit=unitno )
          write(err_mes,'(''trajectory_init('',i2.2,''): failed to read '',a,''; error = '',i6)') dm,trim(filename),ios
          call wrf_error_fatal( trim( err_mes  ) )
        endif
        read(unit=unitno,nml=traj_spec,iostat=ios)
        if( ios /= 0 ) then
          close( unit=unitno )
          write(err_mes,'(''trajectory_init('',i2.2,''): failed to read '',a,''; error = '',i6)') dm,trim(filename),ios
          call wrf_error_fatal( trim( err_mes  ) )
        endif
        close( unit=unitno )
      else input_file
        write(err_mes,'(''trajectory_init('',i2.2,''): failed to find '',a)') dm,trim(filename)
        call wrf_error_fatal( trim( err_mes ) )
      endif input_file

!-----------------------------------------------------------------------------
!  process the namelist input
!-----------------------------------------------------------------------------
      do pkg = 1,pkg_max
        select case( trim(pkg_tag(pkg)) )
          case( 'chm' )
            wrk_def_name(:) = traj_def%chm_name(:)
          case( 'hyd' )
            wrk_def_name(:) = traj_def%hyd_name(:)
          case( 'trc' )
            wrk_def_name(:) = traj_def%trc_name(:)
          case( 'dyn' )
            wrk_def_name(:) = traj_def%dyn_name(:)
          case( 'msc' )
            wrk_def_name(:) = traj_def%msc_name(:)
          case( 'dchm' )
            wrk_def_name(:) = traj_def%dchm_name(:)
        end select
        do m = 1,var_max
          if( wrk_def_name(m) == ' ' ) then
            exit
          endif
        end do
        n_def_var(pkg) = m - 1
      end do

      if( traj_def%start_time /= ' ' ) then
        write(err_mes,'(''trajectory_init('',i2.2,''): default start time = '',a)') dm,traj_def%start_time
        call wrf_message( trim(err_mes) )
      endif
      if( traj_def%stop_time /= ' ' ) then
        write(err_mes,'(''trajectory_init('',i2.2,''): default stop  time = '',a)') dm,traj_def%stop_time
        call wrf_message( trim(err_mes) )
      endif

      do pkg = 1,pkg_max
        if( n_def_var(pkg) > 0 ) then
          write(*,*) ' '
          write(*,'(''trajectory_init('',i2.2,''): default '',a,'' variables'')') dm,pkg_tag(pkg)
          select case( trim(pkg_tag(pkg)) )
            case( 'chm' )
              wrk_def_name(:) = traj_def%chm_name(:)
            case( 'hyd' )
              wrk_def_name(:) = traj_def%hyd_name(:)
            case( 'trc' )
              wrk_def_name(:) = traj_def%trc_name(:)
            case( 'dyn' )
              wrk_def_name(:) = traj_def%dyn_name(:)
            case( 'msc' )
              wrk_def_name(:) = traj_def%msc_name(:)
            case( 'dchm' )
              wrk_def_name(:) = traj_def%dchm_name(:)
          end select
          write(*,*) wrk_def_name(:n_def_var(pkg))
        endif
      end do

      do n_traj = 1,traj_max
        if( traj_type(n_traj)%lon == missing_val .or. &
            traj_type(n_traj)%lat == missing_val .or. &
            traj_type(n_traj)%lev == missing_val ) then
          exit
        endif
      end do
      n_traj = n_traj - 1

has_trajectories: &
      if( n_traj > 0 ) then
!-----------------------------------------------------------------------------
!  set individual trajectories to default if specified
!-----------------------------------------------------------------------------
        if( traj_def%start_time /= ' ' ) then
          do trj = 1,n_traj
            if( traj_type(trj)%start_time == ' ' ) then
              traj_type(trj)%start_time = traj_def%start_time
            endif
          end do
        endif
        if( traj_def%stop_time /= ' ' ) then
          do trj = 1,n_traj
            if( traj_type(trj)%stop_time == ' ' ) then
              traj_type(trj)%stop_time = traj_def%stop_time
            endif
          end do
        endif

        do pkg = 1,pkg_max
          select case( trim(pkg_tag(pkg)) )
            case( 'chm' )
              wrk_def_name(:) = traj_def%chm_name(:)
            case( 'hyd' )
              wrk_def_name(:) = traj_def%hyd_name(:)
            case( 'trc' )
              wrk_def_name(:) = traj_def%trc_name(:)
            case( 'dyn' )
              wrk_def_name(:) = traj_def%dyn_name(:)
            case( 'msc' )
              wrk_def_name(:) = traj_def%msc_name(:)
            case( 'dchm' )
              wrk_def_name(:) = traj_def%dchm_name(:)
          end select
          do trj = 1,n_traj
            select case( trim(pkg_tag(pkg)) )
              case( 'chm' )
                wrk_var_name = traj_type(trj)%chm_spc(1)
              case( 'hyd' )
                wrk_var_name = traj_type(trj)%hyd_spc(1)
              case( 'trc' )
                wrk_var_name = traj_type(trj)%trc_spc(1)
              case( 'dyn' )
                wrk_var_name = traj_type(trj)%dyn_var(1)
              case( 'msc' )
                wrk_var_name = traj_type(trj)%msc_var(1)
              case( 'dchm' )
                wrk_var_name = traj_type(trj)%dchm_spc(1)
            end select
            if( wrk_var_name == ' ' ) then
              m1 = n_def_var(pkg)
              select case( trim(pkg_tag(pkg)) )
                case( 'chm' )
                  traj_type(trj)%chm_spc(:m1) = traj_def%chm_name(:m1)
                case( 'hyd' )
                  traj_type(trj)%hyd_spc(:m1) = traj_def%hyd_name(:m1)
                case( 'trc' )
                  traj_type(trj)%trc_spc(:m1) = traj_def%trc_name(:m1)
                case( 'dyn' )
                  traj_type(trj)%dyn_var(:m1) = traj_def%dyn_name(:m1)
                case( 'msc' )
                  traj_type(trj)%msc_var(:m1) = traj_def%msc_name(:m1)
                case( 'dchm' )
                  traj_type(trj)%dchm_spc(:m1) = traj_def%dchm_name(:m1)
              end select
            endif
          end do
        end do
!-----------------------------------------------------------------------------
!  scan registry real, 2d, 3d variables
!-----------------------------------------------------------------------------
        call reg_scan( grid )
        num_msc => num_msc_dm(dm)

        call get_wrf_debug_level( dbg_lvl )
        if( dbg_lvl > 200 ) then
          write(*,*) ' '
          write(*,'(''trajectory_init('',i2.2,''): Registry 2d,3d variables'')') dm
          do n = 1,num_msc
            write(*,*) trim(St_Vars(n)%varname)
          end do
          n = count( St_Vars(:num_msc)%Stagger == 'X' )
          if( n > 0 ) then
            write(*,*) ' '
            write(*,*) 'Registry 2d,3d variables with staggered X'
            do n = 1,num_msc
              if( St_Vars(n)%Stagger == 'X' ) then
                write(*,*) trim(St_Vars(n)%varname)
              endif
            end do
          endif
          n = count( St_Vars(:num_msc)%Stagger == 'Y' )
          if( n > 0 ) then
            write(*,*) ' '
            write(*,*) 'trajectory_init: Registry 2d,3d variables with staggered Y'
            do n = 1,num_msc
              if( St_Vars(n)%Stagger == 'Y' ) then
                write(*,*) trim(St_Vars(n)%varname)
              endif
            end do
          endif
          n = count( St_Vars(:num_msc)%Stagger == 'Z' )
          if( n > 0 ) then
            write(*,*) ' '
            write(*,*) 'trajectory_init: Registry 2d,3d variables with staggered Z'
            do n = 1,num_msc
              if( St_Vars(n)%Stagger == 'Z' ) then
                write(*,*) trim(St_Vars(n)%varname)
              endif
            end do
          endif
        endif
!-----------------------------------------------------------------------------
!  get variable counts
!-----------------------------------------------------------------------------
        do trj = 1,n_traj
          if( num_chem > 1 ) then
            call get_var_cnt( traj_type(trj)%n_chm_var, traj_type(trj)%chm_spc )
            call get_var_cnt( traj_type(trj)%n_dchm_var, traj_type(trj)%dchm_spc )
          else
            traj_type(trj)%n_chm_var  = 0
            traj_type(trj)%n_dchm_var = 0
          endif
          if( num_moist > 1 ) then
            call get_var_cnt( traj_type(trj)%n_hyd_var, traj_type(trj)%hyd_spc )
          else
            traj_type(trj)%n_hyd_var = 0
          endif
          if( num_tracer > 1 ) then
            call get_var_cnt( traj_type(trj)%n_trc_var, traj_type(trj)%trc_spc )
          else
            traj_type(trj)%n_trc_var = 0
          endif
          if( num_msc > 1 ) then
            call get_var_cnt( traj_type(trj)%n_msc_var, traj_type(trj)%msc_var )
          else
            traj_type(trj)%n_msc_var = 0
          endif
          call get_var_cnt( traj_type(trj)%n_dyn_var, traj_type(trj)%dyn_var )
        end do

        if( any( traj_type(:n_traj)%n_msc_var > 0 ) ) then
          allocate( msc_tbl(num_msc),stat=astat )
          if( astat /= 0 ) then
            call wrf_error_fatal( 'trajectory_init: failed to find allocate msc_tbl' )
          endif
          do m = 1,num_msc
            msc_tbl(m) = trim( St_Vars(m)%Varname )
          end do
        endif
!-----------------------------------------------------------------------------
!  check for trajectory variables in simulation
!-----------------------------------------------------------------------------
        do trj = 1,n_traj
          if( num_chem > 1 ) then
            if( traj_type(trj)%n_chm_var > 0 ) then
              mask(:) = .false.
              call scan_vars( traj_type(trj)%n_chm_var, traj_type(trj)%chm_spc, traj_type(trj)%chm_ndx, &
                              num_chem, chem_dname_table(dm,:), &
                              traj_type(trj)%chm_is_gas, 'chm' )
            endif
            if( traj_type(trj)%n_dchm_var > 0 ) then
              mask(:) = .false.
              call scan_vars( traj_type(trj)%n_dchm_var, traj_type(trj)%dchm_spc, traj_type(trj)%dchm_ndx, &
                              num_chem, chem_dname_table(dm,:), &
                              traj_type(trj)%chm_is_gas, 'chm' )
            endif
          endif
          if( traj_type(trj)%n_hyd_var > 0 .and. num_moist > 1 ) then
            mask(:) = .false.
            call scan_vars( traj_type(trj)%n_hyd_var, traj_type(trj)%hyd_spc, traj_type(trj)%hyd_ndx, &
                            num_moist, moist_dname_table(dm,:), &
                            traj_type(trj)%chm_is_gas, 'hyd' )
          endif
          if( traj_type(trj)%n_trc_var > 0 .and. num_tracer > 1 ) then
            mask(:) = .false.
            call scan_vars( traj_type(trj)%n_trc_var, traj_type(trj)%trc_spc, traj_type(trj)%trc_ndx, &
                            num_tracer, tracer_dname_table(dm,:), &
                            traj_type(trj)%chm_is_gas, 'trc' )
          endif
          if( traj_type(trj)%n_dyn_var > 0 ) then
            mask(:) = .false.
            call scan_vars( traj_type(trj)%n_dyn_var, traj_type(trj)%dyn_var, traj_type(trj)%dyn_ndx, &
                            dyn_max, dyn_var_lst(:), &
                            traj_type(trj)%chm_is_gas, 'dyn' )
          endif
          if( traj_type(trj)%n_msc_var > 0 ) then
            mask(:) = .false.
            call scan_vars( traj_type(trj)%n_msc_var, traj_type(trj)%msc_var, traj_type(trj)%msc_ndx, &
                            num_msc, msc_tbl(:), &
                            traj_type(trj)%chm_is_gas, 'msc' )
          endif
        end do

        do trj = 1,n_traj
          if( traj_type(trj)%n_msc_var > 0 ) then
            do m = 1,traj_type(trj)%n_msc_var
              St_Vars_msk(traj_type(trj)%msc_ndx(m)-offset) = .true.
            end do
          endif
        end do
        m = count( St_Vars_msk(:num_msc) )

        if( allocated( msc_tbl ) ) then
          deallocate( msc_tbl )
        endif

!-----------------------------------------------------------------------------
!  remove trajectories with no variables
!-----------------------------------------------------------------------------
        n_traj_1 = count( (traj_type(:n_traj)%n_chm_var + traj_type(:n_traj)%n_hyd_var &
                           + traj_type(:n_traj)%n_trc_var + traj_type(:n_traj)%n_dyn_var &
                           + traj_type(:n_traj)%n_msc_var + traj_type(:n_traj)%n_dchm_var) > 0 )
        if( n_traj_1 > 0 ) then
          if( n_traj_1 /= n_traj ) then
            trj_mask(1:n_traj) = traj_type(1:n_traj)%in_dom
            m = 1
            do trj = 1,n_traj
              if( trj_mask(trj) ) then
                if( trj /= m ) then
                  traj_type(m) = traj_type(trj)
                  m = m + 1
                endif
              endif
            end do
            n_traj = n_traj_1
          endif
        else
          grid%traj_opt = no_trajectory
          return
        endif
!-----------------------------------------------------------------------------
!  allocate buffer type
!-----------------------------------------------------------------------------
        if( is_root_proc ) then
          if( dm == 1 ) then
            allocate( trj_buff(traj_max,n_dom),stat=astat )
            if( astat /= 0 ) then
              write(err_mes,'(''trajectory_init: failed to allocate traj_buff; error = '',i6)') astat
              call wrf_error_fatal( trim( err_mes  ) )
            endif
          endif
          trj_pbf => trj_buff(:,dm)
        endif
      endif has_trajectories
!   endif master_proc
!-----------------------------------------------------------------------------
!   for dchm package make sure dchm_ndx is ordered list
!-----------------------------------------------------------------------------
    do trj = 1,n_traj
      n = traj_type(trj)%n_dchm_var
      if( n > 0 .and.  any( traj_type(trj)%dchm_ndx(1:n-1) > traj_type(trj)%dchm_ndx(2:n) ) ) then
        trj_mask(:num_chem) = .false. 
        do m = 1,n
          trj_mask(traj_type(trj)%dchm_ndx(m)) = .true. 
        end do
        traj_type(trj)%dchm_ndx(:n) = pack( (/ (m,m=1,num_chem) /),trj_mask(:num_chem) )
        do m = 1,n
          m1 = traj_type(trj)%dchm_ndx(m)
          traj_type(trj)%dchm_spc(m) = trim(chem_dname_table(dm,m1))
        end do
      endif
    end do
!-----------------------------------------------------------------------------
!  if initial run, overwrite existing grid trajectory variables
!-----------------------------------------------------------------------------
#ifdef DM_PARALLEL
!  call wrf_dm_bcast_integer( n_traj,1 )
!  if( n_traj > 0 ) then
!    p_size = (5+var_max)*RWORDSIZE + (2+var_max)*LWORDSIZE &
!           + (5+4*var_max)*IWORDSIZE + 4*32*var_max + 38
!    do trj = 1,n_traj
!      call wrf_dm_bcast_bytes( traj_type(trj),p_size )
!    end do
!  endif
#endif
!  call wrf_abort( 'Debugging' )
is_cold_start: &
   if( .not. rstrt ) then
!-----------------------------------------------------------------------------
!  calc and check trajectory start point
!-----------------------------------------------------------------------------
has_trajectories_a: &
     if( n_traj > 0 ) then
       config_temp = config_flags
       call trajmapproj( grid, config_temp, proj )
       i_end = min(ipe,ide-1)
       j_end = min(jpe,jde-1)
       do j = jps,j_end
         do k = kps,kpe
           z_at_w(ips:i_end,k,j) = (grid%ph_2(ips:i_end,k,j) + grid%phb(ips:i_end,k,j))/g
         end do
       end do
!-----------------------------------------------------------------------------
!  first check x,y
!-----------------------------------------------------------------------------
traj_loop: &
       do trj = 1,n_traj
         if( traj_type(trj)%lat /= missing_val .and. &
             traj_type(trj)%lon /= missing_val ) then
           call latlon_to_ij( proj, traj_type(trj)%lat, traj_type(trj)%lon, x, y )
           traj_type(trj)%in_dom = &
                  (x >= real(ids) .and. x <= real(ide-1) .and. &
                   y >= real(jds) .and. y <= real(jde-1))
#ifdef SW_DEBUG
          write(err_mes,'(''trajectory_init('',i2.2,''): x,y = '',1p,2g15.7)') dm,x,y
          call wrf_debug( 0,trim(err_mes) )
#endif
         else
           traj_type(trj)%in_dom = .false.
         endif
#ifdef SW_DEBUG
         write(err_mes,'(''trajectory_init('',i2.2,''): traj '',i4,'' in dom  = '',l)') dm,trj,traj_type(trj)%in_dom
         call wrf_debug( 0,trim(err_mes) )
#endif
!-----------------------------------------------------------------------------
!  then check z
!-----------------------------------------------------------------------------
is_in_domain: &
         if( traj_type(trj)%in_dom ) then
           i = nint( x )
           j = nint( y )
           traj_type(trj)%in_patch = &
                  (i >= ips .and. i <= min(ipe,ide-1) .and. &
                   j >= jps .and. j <= min(jpe,jde-1))
is_in_patch: &
           if( traj_type(trj)%in_patch ) then
             k1 = kde - 1
!            z_dm_bot = .5*(z_at_w(i,kds,j) + z_at_w(i,kds+1,j))
!            z_dm_top = .5*(z_at_w(i,k1,j) + z_at_w(i,k1+1,j))
             z_dm_bot = z_at_w(i,kds,j)
             z_dm_top = z_at_w(i,k1,j)
             write(err_mes,'(''trajectory_init('',i2.2,''): traj '',i1,'' i,j,z_bot,z_top,lev = '',2i4,1p3g15.7)') &
                dm,trj,i,j,z_dm_bot,z_dm_top,traj_type(trj)%lev
             call wrf_debug( 0,trim(err_mes) )
             if( traj_type(trj)%lev >= z_dm_bot .and. &
                 traj_type(trj)%lev <= z_dm_top ) then
               traj_type(trj)%in_dom = .true.
               traj_type(trj)%x      = x
               traj_type(trj)%y      = y
             else
               traj_type(trj)%in_dom = .false.
             endif
           else is_in_patch
             traj_type(trj)%x = missing_val
             traj_type(trj)%y = missing_val
           endif is_in_patch
         endif is_in_domain
       end do traj_loop
!      n_traj = count( traj_type(:n_traj)%in_dom )
       n_traj_1 = count( traj_type(:n_traj)%in_dom )
     else has_trajectories_a
       n_traj_1 = 0
     endif has_trajectories_a

!-----------------------------------------------------------------------------
!  remove out of domain trajectories
!-----------------------------------------------------------------------------
     write(err_mes,'(''trajectory_init('',i2.2,''): traj cnt = '',2i6)') dm,n_traj_1,n_traj
     call wrf_debug( 0,trim(err_mes) )
     if( n_traj_1 /= n_traj ) then
       trj_mask(1:n_traj) = traj_type(1:n_traj)%in_dom
       m = 1
       do trj = 1,n_traj
         if( trj_mask(trj) ) then
           if( trj /= m ) then
             traj_type(m) = traj_type(trj)
             m = m + 1
           endif
         endif
       end do
       n_traj = n_traj_1
     endif

has_trajectories_b: &
     if( n_traj_1 > 0 ) then
       grid%traj_i(:) = missing_val
       grid%traj_j(:) = missing_val
       grid%traj_k(:) = missing_val
       grid%traj_long(:) = missing_val
       grid%traj_lat(:)  = missing_val
!-----------------------------------------------------------------------------
!  set initial trajectory spatial coordinates
!-----------------------------------------------------------------------------
       k1 = kde - 1
       do trj = 1,n_traj
         if( traj_type(trj)%in_patch ) then
           grid%traj_i(trj) = traj_type(trj)%x
           grid%traj_j(trj) = traj_type(trj)%y
           grid%traj_long(trj) = traj_type(trj)%lon
           grid%traj_lat(trj)  = traj_type(trj)%lat
           i = nint( traj_type(trj)%x )
           j = nint( traj_type(trj)%y )
!          z(kds:k1) = .5*(z_at_w(i,kds:k1,j) + z_at_w(i,kds+1:k1+1,j))
           z(kds:k1) = z_at_w(i,kds:k1,j)
           do k = kds+1,k1
             if( traj_type(trj)%lev <= z(k) ) then
               grid%traj_k(trj) = real(k - 1) &
                              + (traj_type(trj)%lev - z(k-1))/(z(k) - z(k-1))
               exit
             endif
           end do
           write(err_mes,'(''trajectory_init('',i2.2,''): trj,k,z(k-1:k),traj_k = '',2i3,1p3g15.7)') &
              dm,trj,k,z(k-1:k),grid%traj_k(trj)
           call wrf_debug( 0,trim(err_mes) )
         endif
       end do
     else has_trajectories_b
       grid%traj_opt = no_trajectory
       return
     endif has_trajectories_b
   else is_cold_start
     if( n_traj > 0 ) then
       call set_in_dom
     else
       traj_cnt(dm) = n_traj
       return
     endif
   endif is_cold_start

!-----------------------------------------------------------------------------
!  transfer from local variable to module variable
!-----------------------------------------------------------------------------
   traj_cnt(dm) = n_traj
   do trj = 1,n_traj
     trjects(trj) = traj_type(trj)
   end do
!-----------------------------------------------------------------------------
!  create netcdf trajectory file
!-----------------------------------------------------------------------------
   if( .not. rstrt ) then
     do trj = 1,n_traj
#ifdef DM_PARALLEL
       grid%traj_i(trj) = wrf_dm_max_real( grid%traj_i(trj) )
       grid%traj_j(trj) = wrf_dm_max_real( grid%traj_j(trj) )
       grid%traj_k(trj) = wrf_dm_max_real( grid%traj_k(trj) )
       grid%traj_long(trj) = wrf_dm_max_real( grid%traj_long(trj) )
       grid%traj_lat(trj)  = wrf_dm_max_real( grid%traj_lat(trj) )
#else
       grid%traj_i(trj) = grid%traj_i(trj)
       grid%traj_j(trj) = grid%traj_j(trj)
       grid%traj_k(trj) = grid%traj_k(trj)
       grid%traj_long(trj) = grid%traj_long(trj)
       grid%traj_lat(trj)  = grid%traj_lat(trj)
#endif
     end do
   endif

   do pkg = 1,pkg_max
     select case( trim(pkg_tag(pkg)) )
       case( 'chm' )
         pkg_has_vars(:n_traj,pkg) = trjects(:n_traj)%n_chm_var > 0
       case( 'hyd' )
         pkg_has_vars(:n_traj,pkg) = trjects(:n_traj)%n_hyd_var > 0
       case( 'trc' )
         pkg_has_vars(:n_traj,pkg) = trjects(:n_traj)%n_trc_var > 0
       case( 'dyn' )
         pkg_has_vars(:n_traj,pkg) = trjects(:n_traj)%n_dyn_var > 0
       case( 'msc' )
         pkg_has_vars(:n_traj,pkg) = trjects(:n_traj)%n_msc_var > 0
       case( 'dchm' )
         pkg_has_vars(:n_traj,pkg) = trjects(:n_traj)%n_dchm_var > 0
     end select
   end do

!-----------------------------------------------------------------------------
!  setup for  dchm buffer
!-----------------------------------------------------------------------------
   n_dchm_dm(dm) = 0
   if( any( pkg_has_vars(:n_traj,dchm_pkg) ) ) then
     if( .not. allocated( dchm_msk_dm ) ) then
       allocate( dchm_msk_dm(num_chem,n_dom),stat=astat )
       if( astat /= 0 ) then
         write(err_mes,'(''trajectory_init('',i2.2,''): failed to allocate dchm_msk_dm; error = '',i6)') dm,astat
         call wrf_error_fatal( trim( err_mes  ) )
       endif
     endif
     if( .not. allocated( dchm_buf_ndx_dm ) ) then
       allocate( dchm_buf_ndx_dm(num_chem,n_dom),stat=astat )
       if( astat /= 0 ) then
         write(err_mes,'(''trajectory_init('',i2.2,''): failed to allocate dchm_buf_ndx_dm; error = '',i6)') dm,astat
         call wrf_error_fatal( trim( err_mes  ) )
       endif
     endif
     n_dchm   => n_dchm_dm(dm)
     dchm_msk => dchm_msk_dm(:,dm)
     dchm_msk(:) = .false.
     dchm_buf_ndx => dchm_buf_ndx_dm(:,dm)
     dchm_buf_ndx(:) = 0
     do trj = 1,n_traj
       do m = 1,trjects(trj)%n_dchm_var
         dchm_msk(trjects(trj)%dchm_ndx(m)) = .true.
       end do
     end do
     n_dchm = count( dchm_msk )
     dchm_buf_ndx(:n_dchm) = pack( (/ (m,m=1,num_chem) /),dchm_msk(:num_chem) )
     do trj = 1,n_traj
       if( trjects(trj)%n_dchm_var > 0 ) then
         do m1 = 1,trjects(trj)%n_dchm_var
           target = trjects(trj)%dchm_ndx(m1)
           do m2 = 1,n_dchm
             if( target == dchm_buf_ndx(m2) ) then
               trjects(trj)%dchm_ndx(m1) = m2 + offset
               exit
             endif
           end do
         end do
       endif
     end do
   endif

master_proc_a: &
   if( is_root_proc ) then
     if( dm == 1 ) then
       n = max(num_chem,num_moist,num_tracer,num_msc+offset,dyn_max+offset)
       if( n  > offset .and. .not. allocated(trj_msk_dm) ) then
         allocate( trj_msk_dm(traj_max,n,pkg_max,n_dom),stat=astat )
         if( astat /= 0 ) then
           write(err_mes,'(''trajectory_init: failed to allocate trj_msk_dm; error = '',i6)') astat
           call wrf_error_fatal( trim( err_mes  ) )
         endif
         trj_msk_dm(:,:,:,:) = .false.
       endif
     endif
is_initial: &
     if( .not. is_initialized(dm) ) then
!-----------------------------------------------------------------------------
!  allocate buffer arrays
!-----------------------------------------------------------------------------
trj_loop: &
       do trj = 1,n_traj
         if( trjects(trj)%in_dom ) then
pkg_loop:  do pkg = 1,pkg_max
             select case( trim(pkg_tag(pkg)) )
               case( 'chm' )
                 trj_msk => trj_msk_dm(:,:,chm_pkg,dm)
                 m1 = trjects(trj)%n_chm_var
                 if( m1 > 0 ) then
                   allocate( trj_pbf(trj)%chm_vals(vals_max,m1),stat=astat)
                   do m = 1,m1
                     trj_msk(trj,trjects(trj)%chm_ndx(m)) = .true.
                   end do
                 endif
               case( 'hyd' )
                 trj_msk => trj_msk_dm(:,:,hyd_pkg,dm)
                 m1 = trjects(trj)%n_hyd_var
                 if( m1 > 0 ) then
                   allocate( trj_pbf(trj)%hyd_vals(vals_max,m1),stat=astat)
                   do m = 1,m1
                     trj_msk(trj,trjects(trj)%hyd_ndx(m)) = .true.
                   end do
                 endif
               case( 'trc' )
                 trj_msk => trj_msk_dm(:,:,trc_pkg,dm)
                 m1 = trjects(trj)%n_trc_var
                 if( m1 > 0 ) then
                   allocate( trj_pbf(trj)%trc_vals(vals_max,m1),stat=astat)
                   do m = 1,m1
                     trj_msk(trj,trjects(trj)%trc_ndx(m)) = .true.
                   end do
                 endif
               case( 'dyn' )
                 trj_msk => trj_msk_dm(:,:,dyn_pkg,dm)
                 m1 = trjects(trj)%n_dyn_var
                 if( m1 > 0 ) then
                   allocate( trj_pbf(trj)%dyn_vals(vals_max,m1),stat=astat)
                   do m = 1,m1
                     trj_msk(trj,trjects(trj)%dyn_ndx(m)) = .true.
                   end do
                 endif
               case( 'msc' )
                 trj_msk => trj_msk_dm(:,:,msc_pkg,dm)
                 m1 = trjects(trj)%n_msc_var
                 if( m1 > 0 ) then
                   allocate( trj_pbf(trj)%msc_vals(vals_max,m1),stat=astat)
                   do m = 1,m1
                     trj_msk(trj,trjects(trj)%msc_ndx(m)) = .true.
                   end do
                 endif
               case( 'dchm' )
                 trj_msk => trj_msk_dm(:,:,dchm_pkg,dm)
                 m1 = trjects(trj)%n_dchm_var
                 if( m1 > 0 ) then
                   allocate( trj_pbf(trj)%dchm_vals(vals_max,m1),stat=astat)
                   do m = 1,m1
                     trj_msk(trj,trjects(trj)%dchm_ndx(m)) = .true.
                   end do
                 endif
             end select
             if( astat /= 0 ) then
               write(err_mes,'(''trajectory_init: failed to allocate buffer%'',a,''; error = '',i6)') &
                   pkg_tag(pkg),astat
               call wrf_error_fatal( trim( err_mes  ) )
             endif
           end do pkg_loop
         endif
       end do trj_loop

       do pkg = 1,pkg_max
         select case( trim(pkg_tag(pkg)) )
           case( 'chm' )
             trj_msk => trj_msk_dm(:,:,chm_pkg,dm)
             u_lim = num_chem
           case( 'hyd' )
             trj_msk => trj_msk_dm(:,:,hyd_pkg,dm)
             u_lim = num_moist
           case( 'trc' )
             trj_msk => trj_msk_dm(:,:,trc_pkg,dm)
             u_lim = num_tracer
           case( 'dyn' )
             trj_msk => trj_msk_dm(:,:,dyn_pkg,dm)
             u_lim = dyn_max + offset
           case( 'msc' )
             trj_msk => trj_msk_dm(:,:,msc_pkg,dm)
             u_lim = num_msc + offset
           case( 'dchm' )
             trj_msk => trj_msk_dm(:,:,dchm_pkg,dm)
             u_lim = num_chem
         end select
         do trj = 1,n_traj
           trj_msk(trj,1) = any( trj_msk(trj,param_first_scalar:u_lim) )
         end do
       end do
       is_initialized(dm) = .true.
       if( .not. rstrt ) then
         call trajectory_create_file( grid, n_traj )
       endif
     endif is_initial
   endif master_proc_a
#ifdef DM_PARALLEL
   call wrf_dm_bcast_logical( is_initialized,n_dom )
#endif

   CONTAINS

   subroutine reg_scan( grid )

   use module_domain_type, only : domain, fieldlist

!-----------------------------------------------------------------------------
!  dummy arguments
!-----------------------------------------------------------------------------
   type(domain), intent(in) :: grid

!-----------------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------------
   integer         :: astat
   integer         :: cnt
   integer         :: dm, dm_ndx
   integer         :: m, n
   type(fieldlist), pointer :: p
   logical         :: valid
   type(statevar), allocatable :: St_Vars_wrk(:,:)
   logical, allocatable :: St_Vars_msk_wrk(:,:)

!-----------------------------------------------------------------------------
!  just get the count
!-----------------------------------------------------------------------------
   dm = grid%id
   p => grid%head_statevars%next ; cnt = 0
   do while( associated(p) )
     valid = (p%Type == 'R' .or. p%Type == 'r') .and. (p%Ndim == 2 .or. p%Ndim == 3)
     if( valid ) then
       valid = (p%Ndim == 2 .and. p%MemoryOrder(:2) == 'XY') .or. &
               (p%Ndim == 3 .and. p%MemoryOrder(:3) == 'XZY')
       if( valid .and. p%Ndim == 3 ) then
         valid = p%em2 == kme
       endif
     endif
     if( valid ) then
       cnt = cnt + 1
     endif
     p => p%next
   end do

   write(*,'(''reg_scan('',i2.2,''): found '',i4,'' state variables '')') dm,cnt

!-----------------------------------------------------------------------------
!  now allocate and set St_Vars
!-----------------------------------------------------------------------------
   num_msc_dm(dm) = cnt
   if( cnt > 0 ) then
     if( .not. allocated( St_Vars_dm ) ) then
       allocate( St_Vars_dm(cnt,n_dom),St_Vars_msk_dm(cnt,n_dom),stat=astat )
       if( astat /= 0 ) then
         call wrf_error_fatal( 'reg_scan: failed to allocate St_Vars,St_Vars_msk' )
       endif
     elseif( cnt > maxval(num_msc_dm(1:dm-1)) ) then
       n = size( St_vars_dm,dim=1 )
       allocate( St_Vars_wrk(n,dm-1),St_Vars_msk_wrk(n,dm-1),stat=astat )
       if( astat /= 0 ) then
         call wrf_error_fatal( 'reg_scan: failed to allocate St_Vars,St_Vars_msk wrk arrays' )
       endif
       do dm_ndx = 1,dm-1
         do m = 1,n
           St_Vars_wrk(m,dm_ndx) = St_Vars_dm(m,dm_ndx)
           St_Vars_msk_wrk(m,dm_ndx) = St_Vars_msk_dm(m,dm_ndx)
         end do
       end do
       deallocate( St_vars_dm,St_Vars_msk_dm )
       allocate( St_Vars_dm(cnt,n_dom),St_Vars_msk_dm(cnt,n_dom),stat=astat )
       if( astat /= 0 ) then
         call wrf_error_fatal( 'reg_scan: failed to allocate St_Vars,St_Vars_msk' )
       endif
       do dm_ndx = 1,dm-1
         do m = 1,n
           St_Vars_dm(m,dm_ndx) = St_Vars_wrk(m,dm_ndx)
           St_Vars_msk_dm(m,dm_ndx) = St_Vars_msk_wrk(m,dm_ndx)
         end do
       end do
       deallocate( St_vars_wrk,St_Vars_msk_wrk )
     endif

     St_Vars     => St_Vars_dm(:,dm)
     St_Vars_msk => St_Vars_msk_dm(:,dm)

     St_Vars_msk(:cnt) = .false.
     p => grid%head_statevars%next ; cnt = 0
     do while( associated(p) )
       valid = (p%Type == 'R' .or. p%Type == 'r') .and. (p%Ndim == 2 .or. p%Ndim == 3)
       if( valid ) then
         valid = (p%Ndim == 2 .and. p%MemoryOrder(:2) == 'XY') .or. &
                 (p%Ndim == 3 .and. p%MemoryOrder(:3) == 'XZY')
         if( valid .and. p%Ndim == 3 ) then
           valid = p%em2 == kme
         endif
         if( valid ) then
           cnt = cnt + 1
           St_Vars(cnt)%Varname = p%Varname 
           St_Vars(cnt)%Description = p%Description 
           St_Vars(cnt)%Units   = p%Units 
           St_Vars(cnt)%MemOrd  = p%MemoryOrder 
           St_Vars(cnt)%Stagger = p%Stagger
           St_Vars(cnt)%Ndim = p%Ndim 
           if( p%Ndim == 3 ) then
             St_Vars(cnt)%rfield_3d => p%rfield_3d
             nullify( St_Vars(cnt)%rfield_2d )
           else
             St_Vars(cnt)%rfield_2d => p%rfield_2d
             nullify( St_Vars(cnt)%rfield_3d )
           endif
         endif
       endif
       p => p%next
     end do
   endif

   end subroutine reg_scan

   subroutine get_var_cnt( n_vars, vars )

!-----------------------------------------------------------------------------
!  dummy arguments
!-----------------------------------------------------------------------------
   integer, intent(inout) :: n_vars
   character(len=32), intent(inout)  :: vars(:)

   mask(:) = vars(:) /= ' '
   wrk_chr(:) = ' '
   m1 = 0
   do n = 1,var_max
     if( mask(n) ) then
       m1 = m1 + 1
       wrk_chr(m1) = vars(n)
     endif
   end do

   n_vars  = count( wrk_chr(:) /= ' ' )
   vars(:) = wrk_chr(:)

   end subroutine get_var_cnt

   subroutine scan_vars( n_vars, vars, var_ndx, n_tbl, tbl, &
                         is_gas, var_type )

!-----------------------------------------------------------------------------
!  dummy arguments
!-----------------------------------------------------------------------------
   integer, intent(inout) :: n_vars
   integer, intent(in)    :: n_tbl
   integer, intent(inout) :: var_ndx(:)
   logical, intent(inout) :: is_gas(:)
   character(len=*), intent(in)      :: var_type
   character(len=32), intent(inout)  :: vars(:)
   character(len=256), intent(in)    :: tbl(:)

!-----------------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------------
   integer :: ml
   integer :: ndx(var_max)

var_loop: &
   do n = 1,n_vars
     var_name = vars(n)
     if( trim(var_type) == 'chm' ) then
       i = index( '(a)', var_name )
       if( i == 0 ) then
         i = index( '(A)', var_name )
       endif
       if( i > 0 ) then
         is_gas(n) = .false.
         var_name(i:) = ' '
         vars(n)(i:)  = ' '
       endif
     endif
     if( trim(var_type) == 'dyn' .or. trim(var_type) == 'msc' ) then
       ml = 1
     else
       ml = param_first_scalar
     endif
!    if( trim(var_type) /= 'dyn' ) then
!      ml = param_first_scalar
!    else
!      ml = 1
!    endif
     do m1 = ml,n_tbl
       if( trim(var_name) == trim(tbl(m1)) ) then
         mask(n) = .true.
         ndx(n)  = m1
         exit
       endif 
     end do
     if( .not. mask(n) ) then
       write(err_mes,'(''scan_vars('',i2.2,''): '',a,'' not in '',a,'' opt'')') dm,trim(var_name),var_type
       call wrf_message( trim(err_mes) )
     endif
   end do var_loop

   if( trim(var_type) == 'dyn' .or. trim(var_type) == 'msc' ) then
     ndx(1:n_vars) = ndx(1:n_vars) + offset
   endif

   wrk_chr(:) = ' '
   m1 = 0
   do n = 1,var_max
     if( mask(n) ) then
       m1 = m1 + 1
       wrk_chr(m1) = vars(n)
     endif
   end do

   var_ndx(:count(mask)) = pack( ndx(:),mask=mask)
   vars(:) = wrk_chr(:)
   n_vars  = count( mask )

   end subroutine scan_vars

   subroutine set_in_dom

!-----------------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------------
   integer :: ncid
   integer :: ios
   integer :: time_id
   integer :: varid
   integer :: time_ndx

   real    :: traj_i(n_traj)
   real    :: traj_j(n_traj)
   real    :: traj_k(n_traj)

   character(len=256) :: err_mes
   character(len=256) :: filename

include 'netcdf.inc'

!---------------------------------------------------------------------
!  open netcdf file
!---------------------------------------------------------------------
   write(filename,'(''wrfout_traj_d'',i2.2)',iostat=ios) dm
   if( ios /= 0 ) then
     write(err_mes,'(''set_in_dom: failed to set filename: error = '',i6)') ios
     call wrf_error_fatal( trim( err_mes  ) )
   endif
   ios = nf_open( trim(filename), nf_nowrite, ncid )
   if( ios /= 0 ) then
     write(err_mes,'(''set_in_dom: failed to open '',a,'': error = '',i6)') trim(filename),ios
     call wrf_error_fatal( trim( err_mes  ) )
   endif

!---------------------------------------------------------------------
!  get current time index
!---------------------------------------------------------------------
   err_mes = 'set_in_dom: failed to get time id'
   call handle_ncerr( nf_inq_dimid( ncid, 'time', time_id ),trim(err_mes) )
   err_mes = 'set_in_dom: failed to get time dimension'
   call handle_ncerr( nf_inq_dimlen( ncid, time_id, time_ndx ),trim(err_mes) )

!---------------------------------------------------------------------
!  read in last traj_{i,j,k} variables
!---------------------------------------------------------------------
   err_mes = 'set_in_dom: failed to get traj_i id'
   call handle_ncerr( nf_inq_varid( ncid, 'traj_i', varid ),trim(err_mes) )
   err_mes = 'set_in_dom: failed to get traj_i'
   call handle_ncerr( nf_get_vara_real( ncid, varid, (/ 1,time_ndx /), (/ n_traj,1 /), &
                                        traj_i ),trim(err_mes) )
   err_mes = 'set_in_dom: failed to get traj_j id'
   call handle_ncerr( nf_inq_varid( ncid, 'traj_j', varid ),trim(err_mes) )
   err_mes = 'set_in_dom: failed to get traj_j'
   call handle_ncerr( nf_get_vara_real( ncid, varid, (/ 1,time_ndx /), (/ n_traj,1 /), &
                                        traj_j ),trim(err_mes) )
   err_mes = 'set_in_dom: failed to get traj_k id'
   call handle_ncerr( nf_inq_varid( ncid, 'traj_k', varid ),trim(err_mes) )
   err_mes = 'set_in_dom: failed to get traj_k'
   call handle_ncerr( nf_get_vara_real( ncid, varid, (/ 1,time_ndx /), (/ n_traj,1 /), &
                                        traj_k ),trim(err_mes) )
   traj_type(1:n_traj)%in_dom = traj_i(1:n_traj) /= missing_val .and.  traj_j(1:n_traj) /= missing_val &
                                                                .and.  traj_k(1:n_traj) /= missing_val

   ios = nf_close( ncid )

   end subroutine set_in_dom

   end subroutine trajectory_init

#ifdef NETCDF
   subroutine trajectory_driver( grid )

#ifdef DM_PARALLEL
   use module_dm, only : &
                  local_communicator, mytask, ntasks, ntasks_x, ntasks_y                   &
                 ,local_communicator_periodic, wrf_dm_max_real, wrf_dm_max_int
   use module_comm_dm, only : halo_em_chem_e_3_sub, halo_em_moist_e_3_sub
   use module_comm_dm, only : halo_em_tracer_e_3_sub
#endif
   use module_domain
   use module_date_time
   use module_state_description, only : num_chem
   use module_state_description, only : num_moist, num_tracer, param_first_scalar

!-----------------------------------------------------------------------------
!  dummy arguments
!-----------------------------------------------------------------------------
   type(domain), intent(in) :: grid

!-----------------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------------
   integer :: ims,ime, jms,jme, kms,kme
   integer :: ids,ide, jds,jde, kds,kde
   integer :: ips,ipe, jps,jpe, kps,kpe
   integer :: dm
   integer :: j, k
   integer :: il, iu, ios, jl, ju, kl, ku
   integer :: m, mu, n, ndx, n_vars, n_traj
   integer :: pkg_var_cnt(traj_max)
   integer :: ncid
   integer :: pkg, trj
   integer :: n_msc_buf
   integer :: num_chem_sav
   integer, pointer :: num_msc      ! number misc species
   integer, pointer :: n_dchm       ! number dchm buffer species
   integer, pointer :: dchm_buf_ndx(:)
   integer :: St_Vars_ndx
   integer, allocatable :: St_Vars_buf_ndx(:)
#ifndef DM_PARALLEL
   integer :: mytask
#endif
   integer :: traj_proc(traj_max), glb_traj_proc(traj_max)
   real :: dchm_fill_val(traj_max)
   real :: x, y, zi
   real :: delsx, delsy, o_delsx, o_delsy
   real :: delsz, o_delsz
   real :: max_conc
   real :: horz_conc(2)
   real, pointer :: traj_conc(:)
   real, target  :: traj_val(var_max,traj_max)
   real, pointer :: wrk4d(:,:,:,:)
   real, allocatable, target :: chem(:,:,:,:)
   real, allocatable, target :: moist(:,:,:,:)
   real, allocatable, target :: tracer(:,:,:,:)
   character(len=19)  :: current_timestr, next_timestr
   character(len=32)  :: var_name(var_max)
   character(len=256) :: err_mes
   logical :: has_dchm
   logical :: is_root_proc
   logical :: is_in_patch_gap
   logical :: flsh_buff
   logical :: found
   logical :: traj_is_active(traj_max)
   logical :: pkg_is_active(traj_max,pkg_max)
   logical, pointer :: pkg_has_vars(:,:)

   logical, external :: wrf_dm_on_monitor

   type(WRFU_Time) :: current_time, next_time
   type(WRFU_Time) :: start_time, stop_time

include 'netcdf.inc'
#ifdef DM_PARALLEL
include 'mpif.h'
#endif

#ifndef DM_PARALLEL
   mytask = 0
#endif
   dm = grid%id
   n_traj = traj_cnt(dm)
has_trajectories: &
   if( n_traj > 0 ) then
     St_Vars => St_Vars_dm(:,dm)
     St_Vars_msk => St_Vars_msk_dm(:,dm)
     num_msc => num_msc_dm(dm)
     trjects => traject(:,dm)
     n_vals  => n_vals_dm(dm)
     is_root_proc = wrf_dm_on_monitor()
     if( is_root_proc ) then
       trj_pbf => trj_buff(:,dm)
     endif
     has_dchm = any( trjects(:n_traj)%n_dchm_var > 0 )
     if( has_dchm ) then
       n_dchm => n_dchm_dm(dm)
       dchm_buf_ndx => dchm_buf_ndx_dm(:,dm)
     endif

     call get_ijk_from_grid( grid ,                   &
                             ids, ide, jds, jde, kds, kde,    &
                             ims, ime, jms, jme, kms, kme,    &
                             ips, ipe, jps, jpe, kps, kpe    )

!-----------------------------------------------------------------------------
!  is trajectory in time interval?
!-----------------------------------------------------------------------------
     call domain_clock_get( grid, current_time=current_time, current_timestr=current_timestr)
     call geth_newdate( next_timestr, current_timestr, int(grid%dt) )
     call wrf_atotime( next_timestr(1:19), next_time )
     do trj = 1,n_traj
       call wrf_atotime( traject(trj,dm)%start_time(1:19), start_time )
       call wrf_atotime( traject(trj,dm)%stop_time(1:19), stop_time )
       traj_is_active(trj) = next_time .ge. start_time .and. next_time .le. stop_time
     end do
!-----------------------------------------------------------------------------
!  is trajectory in domain?
!-----------------------------------------------------------------------------
     pkg_has_vars => pkg_has_vars_dm(:,:,dm)
     do trj = 1,n_traj
       if( traj_is_active(trj) ) then
         trjects(trj)%in_dom = grid%traj_i(trj) >= real(ids) .and. grid%traj_i(trj) <= real(ide-1)
         if( trjects(trj)%in_dom ) then
           trjects(trj)%in_dom = grid%traj_j(trj) >= real(jds) .and. grid%traj_j(trj) <= real(jde-1)
         endif
         if( trjects(trj)%in_dom ) then
           trjects(trj)%in_dom = grid%traj_k(trj) >= real(kps) .and. grid%traj_k(trj) <= real( min( kpe,kde-1 ) )
         endif
         traj_is_active(trj) = trjects(trj)%in_dom
       endif
     end do
     do pkg = 1,pkg_max
       pkg_is_active(:n_traj,pkg) = traj_is_active(:n_traj) .and. pkg_has_vars(:n_traj,pkg)
     end do
!-----------------------------------------------------------------------------
!  check whether this is a chemistry time step
!-----------------------------------------------------------------------------
     dchm_fill_val(:n_traj) = missing_val
     if( .not. do_chemstep ) then
       do trj = 1,n_traj
         if( pkg_is_active(trj,dchm_pkg) ) then
           dchm_fill_val(trj) = zero_val
         endif
       end do
       pkg_is_active(:n_traj,dchm_pkg) = .false.
     endif
!-----------------------------------------------------------------------------
!  is trajectory in mpi process?
!-----------------------------------------------------------------------------
     traj_proc(:n_traj) = -1
     do trj = 1,n_traj
       if( traj_is_active(trj) ) then
         trjects(trj)%in_patch = grid%traj_i(trj) >= real(ips) .and. grid%traj_i(trj) <= real( min( ipe,ide-1 ) )
         if( trjects(trj)%in_patch ) then
           trjects(trj)%in_patch = grid%traj_j(trj) >= real(jps) .and. grid%traj_j(trj) <= real( min( jpe,jde-1 ) )
         endif
         if( trjects(trj)%in_patch ) then
           trjects(trj)%in_patch = grid%traj_k(trj) >= real(kps) .and. grid%traj_k(trj) <= real( min( kpe,kde-1 ) )
         endif
         if( trjects(trj)%in_patch ) then
           traj_proc(trj) = mytask + 1
         else
           traj_proc(trj) = 0
         endif
       endif
     end do
     do trj = 1,n_traj
       glb_traj_proc(trj) = wrf_dm_max_int( traj_proc(trj) )
     end do
!-----------------------------------------------------------------------------
!  any trajectories in "gap" between patches?
!-----------------------------------------------------------------------------
     do trj = 1,n_traj
       if( traj_is_active(trj) .and. glb_traj_proc(trj) == 0 ) then
         trjects(trj)%in_patch = grid%traj_i(trj) >= real(ips) .and. grid%traj_i(trj) <= real( min( ipe+1,ide-1 ) )
         if( trjects(trj)%in_patch ) then
           trjects(trj)%in_patch = grid%traj_j(trj) >= real(jps) .and. grid%traj_j(trj) <= real( min( jpe+1,jde-1 ) )
         endif
         if( trjects(trj)%in_patch ) then
           trjects(trj)%in_patch = grid%traj_k(trj) >= real(kps) .and. grid%traj_k(trj) <= real( min( kme,kde-1 ) )
         endif
         if( trjects(trj)%in_patch ) then
           traj_proc(trj) = mytask + 1
         else
           traj_proc(trj) = 0
         endif
         if( traj_proc(trj) /= 0 ) then
           call wrf_debug( 0,'^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^')
           write(err_mes,'(''Gapper '',i5,''; x,y,zi = '',1p,3g15.7)') trj,grid%traj_i(trj),grid%traj_j(trj),grid%traj_k(trj)
           call wrf_debug( 0,trim(err_mes) )
           write(err_mes,'(''Gapper ips,ipe,jps,jpe = '',4i5)') ips,ipe,jps,jpe
           call wrf_debug( 0,trim(err_mes) )
           call wrf_debug( 0,'^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^')
         endif
       endif
     end do

!-----------------------------------------------------------------------------
!  buffer traj time and position
!-----------------------------------------------------------------------------
     if( is_root_proc ) then
       n_vals = n_vals + 1
       if( grid%itimestep > 0 ) then
         trj_pbf(:n_traj)%times(n_vals) = next_timestr
       else
         trj_pbf(:n_traj)%times(n_vals) = current_timestr
       endif
       do trj = 1,n_traj
         trj_pbf(trj)%trj_i(n_vals) = grid%traj_i(trj)
         trj_pbf(trj)%trj_j(n_vals) = grid%traj_j(trj)
         trj_pbf(trj)%trj_k(n_vals) = grid%traj_k(trj)
         trj_pbf(trj)%trj_lons(n_vals) = grid%traj_long(trj)
         trj_pbf(trj)%trj_lats(n_vals) = grid%traj_lat(trj)
       end do
     endif

     do trj = 1,n_traj
       traj_val(:,trj) = missing_val
     end do

pkg_loop: &
     do pkg = 1,pkg_max 
pkg_has_active_traj: &
       if( any( pkg_is_active(:n_traj,pkg) ) ) then
!-----------------------------------------------------------------------------
!  allocate working data array
!-----------------------------------------------------------------------------
         select case( trim(pkg_tag(pkg)) )
           case( 'chm' )
             allocate( chem(ims:ime,kms:kme,jms:jme,num_chem),stat=ios )
           case( 'hyd' )
             allocate( moist(ims:ime,kms:kme,jms:jme,num_moist),stat=ios )
           case( 'trc' )
             allocate( tracer(ims:ime,kms:kme,jms:jme,num_tracer),stat=ios )
           case( 'dyn' )
             allocate( chem(ims:ime,kms:kme,jms:jme,dyn_max+offset+2),stat=ios )
           case( 'msc' )
             m = count( St_Vars_msk(:num_msc) )
             if( m > 0 ) then
               allocate( chem(ims:ime,kms:kme,jms:jme,m+offset), &
                         St_Vars_buf_ndx(m+offset),stat=ios )
             endif
           case( 'dchm' )
             if( n_dchm > 0 ) then
               allocate( chem(ims:ime,kms:kme,jms:jme,n_dchm+offset),stat=ios )
             endif
         end select
         if( ios /= 0 ) then
           write(err_mes,'(''trajectory_driver: failed to allocate wrk4d: error = '',i6)') ios
           call wrf_error_fatal( trim( err_mes  ) )
         endif
         select case( trim(pkg_tag(pkg)) )
           case( 'chm' )
             do m = 1,num_chem
               do j = jps,jpe
                 do k = kps,kpe
                   chem(ips:ipe,k,j,m) = grid%chem(ips:ipe,k,j,m)
                 end do
               end do
             end do
           case( 'hyd' )
             do m = 1,num_moist
               do j = jps,jpe
                 do k = kps,kpe
                   moist(ips:ipe,k,j,m) = grid%moist(ips:ipe,k,j,m)
                 end do
               end do
             end do
           case( 'trc' )
             do m = 1,num_tracer
               do j = jps,jpe
                 do k = kps,kpe
                   tracer(ips:ipe,k,j,m) = grid%tracer(ips:ipe,k,j,m)
                 end do
               end do
             end do
           case( 'dyn' )
             call pack_dyn_vals
           case( 'msc' )
             n_msc_buf = 1
             do m = 1,num_msc
               if( St_Vars_msk(m) ) then
                 n_msc_buf = n_msc_buf + 1
                 St_Vars_buf_ndx(n_msc_buf) = m
                 do j = jps,jpe
                   do k = kps,kpe
                     chem(ips:ipe,k,j,n_msc_buf) = St_Vars(m)%rfield_3d(ips:ipe,k,j)
                   end do
                 end do
               endif
             end do
           case( 'dchm' )
             do m = param_first_scalar,n_dchm+offset
               do j = jps,jpe
                 do k = kps,kpe
                   chem(ips:ipe,k,j,m) = dchm_buff(ips:ipe,k,j,m)
                 end do
               end do
             end do
         end select
#ifdef DM_PARALLEL
!-----------------------------------------------------------------------------
!  any trajectories in "gap" between patches?
!-----------------------------------------------------------------------------
             is_in_patch_gap = any( glb_traj_proc(:n_traj) == 0 .and. pkg_is_active(:n_traj,pkg) )
             if( is_in_patch_gap ) then
!              write(err_mes,'(''glb_traj_proc mask = '',10i5)') glb_traj_proc(:n_traj)
!              call wrf_debug( 0,trim(err_mes) )
               select case( trim(pkg_tag(pkg)) )
                 case( 'chm' )
#      include "HALO_EM_CHEM_E_3.inc"
                 case( 'hyd' )
#      include "HALO_EM_MOIST_E_3.inc"
                 case( 'trc' )
#      include "HALO_EM_TRACER_E_3.inc"
                 case( 'dyn' )
                   num_chem_sav = num_chem
                   num_chem     = dyn_max + offset + 2
#      include "HALO_EM_CHEM_E_3.inc"
                   num_chem = num_chem_sav
                 case( 'msc' )
                   num_chem_sav = num_chem
                   num_chem     = n_msc_buf
#      include "HALO_EM_CHEM_E_3.inc"
                   num_chem = num_chem_sav
                 case( 'dchm' )
                   num_chem_sav = num_chem
                   num_chem     = n_dchm
#      include "HALO_EM_CHEM_E_3.inc"
                   num_chem = num_chem_sav
               end select
             endif
#endif

traj_loop: &
         do trj = 1,n_traj
           select case( trim(pkg_tag(pkg)) )
             case( 'chm' )
               n_vars = traject(trj,dm)%n_chm_var
             case( 'hyd' )
               n_vars = traject(trj,dm)%n_hyd_var
             case( 'trc' )
               n_vars = traject(trj,dm)%n_trc_var
             case( 'dyn' )
               n_vars = traject(trj,dm)%n_dyn_var
             case( 'msc' )
               n_vars = traject(trj,dm)%n_msc_var
             case( 'dchm' )
               n_vars = traject(trj,dm)%n_dchm_var
           end select
pkg_is_active_in_traj: &
           if( pkg_is_active(trj,pkg) ) then
             select case( trim(pkg_tag(pkg)) )
               case( 'chm', 'dyn', 'msc', 'dchm' )
                 wrk4d => chem
               case( 'hyd' )
                 wrk4d => moist
               case( 'trc' )
                 wrk4d => tracer
             end select
in_patch:    if( traj_proc(trj) == mytask+1 ) then
               x = grid%traj_i(trj)
               y = grid%traj_j(trj)
               zi = grid%traj_k(trj)
               il = int( x ) ; iu = il + 1
               jl = int( y ) ; ju = jl + 1
               kl = int( zi ) ; ku = kl + 1
               delsx = x - floor(x) ; o_delsx = 1. - delsx
               delsy = y - floor(y) ; o_delsy = 1. - delsy
               delsz = zi - floor(zi) ; o_delsz = 1. - delsz
var_loop:      do n = 1,n_vars
                 found = .true.
                 select case( trim(pkg_tag(pkg)) )
                   case( 'chm' )
                     ndx = traject(trj,dm)%chm_ndx(n)
                   case( 'hyd' )
                     ndx = traject(trj,dm)%hyd_ndx(n)
                   case( 'trc' )
                     ndx = traject(trj,dm)%trc_ndx(n)
                   case( 'dyn' )
                     ndx = 1
                     call set_dyn_vals
                   case( 'msc' )
                     found = .false.
                     St_Vars_ndx = trjects(trj)%msc_ndx(n) - offset
                     do ndx = param_first_scalar,n_msc_buf
                       if( St_Vars_ndx == St_Vars_buf_ndx(ndx) ) then
                         found = .true.
                         exit
                       endif
                     end do
                     if( found ) then
                       call set_msc_vals
                     endif
                     ndx = 1
                   case( 'dchm' )
                     ndx = trjects(trj)%dchm_ndx(n)
                 end select
                 if( found ) then
                   horz_conc(1) = o_delsx*o_delsy*wrk4d(il,kl,jl,ndx) + o_delsy*delsx*wrk4d(iu,kl,jl,ndx) &
                                + delsy*o_delsx*wrk4d(il,kl,ju,ndx) + delsx*delsy*wrk4d(iu,kl,ju,ndx)
                   horz_conc(2) = o_delsx*o_delsy*wrk4d(il,ku,jl,ndx) + o_delsy*delsx*wrk4d(iu,ku,jl,ndx) &
                                + delsy*o_delsx*wrk4d(il,ku,ju,ndx) + delsx*delsy*wrk4d(iu,ku,ju,ndx)
                   traject(trj,dm)%traj_var(n) = delsz*horz_conc(2) + o_delsz*horz_conc(1)
                 else
                   traject(trj,dm)%traj_var(n) = missing_val
                 endif
               end do var_loop
             else in_patch
               traject(trj,dm)%traj_var(:n_vars) = missing_val
             endif in_patch
             traj_conc => traj_val(:,trj)
             do n = 1,n_vars
#ifdef DM_PARALLEL
               max_conc = wrf_dm_max_real( traject(trj,dm)%traj_var(n) )
#else
               max_conc = traject(trj,dm)%traj_var(n)
#endif
               if( is_root_proc ) then
                 traj_conc(n) = max_conc
               endif
             end do
           else pkg_is_active_in_traj
             if( is_root_proc .and. n_vars > 0 ) then
               traj_conc => traj_val(:,trj)
               traj_conc(:n_vars) = missing_val
             endif
           endif pkg_is_active_in_traj
!-----------------------------------------------------------------------------
!  buffer traj chm,trc,hyb,msc,dchm variables
!-----------------------------------------------------------------------------
           if( is_root_proc .and. n_vars > 0 ) then
             select case( pkg_tag(pkg) )
               case( 'chm' )
                 trj_pbf(trj)%chm_vals(n_vals,:n_vars) = traj_conc(:n_vars)
               case( 'hyd' )
                 trj_pbf(trj)%hyd_vals(n_vals,:n_vars) = traj_conc(:n_vars)
               case( 'trc' )
                 trj_pbf(trj)%trc_vals(n_vals,:n_vars) = traj_conc(:n_vars)
               case( 'dyn' )
                 trj_pbf(trj)%dyn_vals(n_vals,:n_vars) = traj_conc(:n_vars)
               case( 'msc' )
                 trj_pbf(trj)%msc_vals(n_vals,:n_vars) = traj_conc(:n_vars)
               case( 'dchm' )
                 trj_pbf(trj)%dchm_vals(n_vals,:n_vars) = traj_conc(:n_vars)
             end select
           endif
         end do traj_loop
         if( allocated( chem ) ) then
           deallocate( chem )
         endif
         if( allocated( moist ) ) then
           deallocate( moist )
         endif
         if( allocated( tracer ) ) then
           deallocate( tracer )
         endif
         if( allocated( St_Vars_buf_ndx ) ) then
           deallocate( St_Vars_buf_ndx )
         endif
         if( pkg == dchm_pkg .and. allocated( dchm_buff ) ) then
           deallocate( dchm_buff )
         endif
       else pkg_has_active_traj
         if( is_root_proc ) then
           do trj = 1,n_traj
             select case( trim(pkg_tag(pkg)) )
               case( 'chm' )
                 n_vars = traject(trj,dm)%n_chm_var
               case( 'hyd' )
                 n_vars = traject(trj,dm)%n_hyd_var
               case( 'trc' )
                 n_vars = traject(trj,dm)%n_trc_var
               case( 'dyn' )
                 n_vars = traject(trj,dm)%n_dyn_var
               case( 'msc' )
                 n_vars = traject(trj,dm)%n_msc_var
               case( 'dchm' )
                 n_vars = traject(trj,dm)%n_dchm_var
             end select
             if( n_vars > 0 ) then
               select case( pkg_tag(pkg) )
                 case( 'chm' )
                   trj_pbf(trj)%chm_vals(n_vals,:n_vars) = missing_val
                 case( 'hyd' )
                   trj_pbf(trj)%hyd_vals(n_vals,:n_vars) = missing_val
                 case( 'trc' )
                   trj_pbf(trj)%trc_vals(n_vals,:n_vars) = missing_val
                 case( 'dyn' )
                   trj_pbf(trj)%dyn_vals(n_vals,:n_vars) = missing_val
                 case( 'msc' )
                   trj_pbf(trj)%msc_vals(n_vals,:n_vars) = missing_val
                 case( 'dchm' )
                   trj_pbf(trj)%dchm_vals(n_vals,:n_vars) = dchm_fill_val(trj)
               end select
             endif
           end do
         endif
       endif pkg_has_active_traj
     end do pkg_loop
!-----------------------------------------------------------------------------
!  output trajectory buffer
!-----------------------------------------------------------------------------
     if( is_root_proc ) then
       flsh_buff = n_vals == vals_max .or. domain_last_time_step( grid )
       if( flsh_buff ) then
         call trajectory_write_file( n_traj, n_vals, dm )
       endif
     endif
   endif has_trajectories

   CONTAINS

   subroutine pack_dyn_vals

   integer :: mp1

   do m = 1,dyn_max+2
     mp1 = m + 1
     select case( m )
       case( 1 )
         do j = jps,jpe
           do k = kps,kpe
             chem(ips:ipe,k,j,mp1) = grid%p(ips:ipe,k,j)
           end do
         end do
       case( 2 )
         do j = jps,jpe
           do k = kps,kpe
             chem(ips:ipe,k,j,mp1) = grid%t_2(ips:ipe,k,j)
           end do
         end do
       case( 3 )
         do j = jps,jpe
           do k = kps,kpe
             chem(ips:ipe,k,j,mp1) = grid%pb(ips:ipe,k,j)
           end do
         end do
       case( 4 )
         do j = jps,jpe
           do k = kps,kpe
             chem(ips:ipe,k,j,mp1) = grid%u_2(ips:ipe,k,j)
           end do
         end do
       case( 5 )
         do j = jps,jpe
           do k = kps,kpe
             chem(ips:ipe,k,j,mp1) = grid%v_2(ips:ipe,k,j)
           end do
         end do
       case( 6 )
         do j = jps,jpe
           do k = kps,kpe
             chem(ips:ipe,k,j,mp1) = grid%w_2(ips:ipe,k,j)
           end do
         end do
       case( 7 )
         do j = jps,jpe
           do k = kps,kpe
             chem(ips:ipe,k,j,mp1) = grid%ph_2(ips:ipe,k,j)
           end do
         end do
       case( 8 )
         do j = jps,jpe
           do k = kps,kpe
             chem(ips:ipe,k,j,mp1) = grid%phb(ips:ipe,k,j)
           end do
         end do
       case( 9 )
         do j = jps,jpe
           do k = kps,kpe
             chem(ips:ipe,k,j,mp1) = grid%rainprod(ips:ipe,k,j)
           end do
         end do
       case( 10 )
         do j = jps,jpe
           do k = kps,kpe
             chem(ips:ipe,k,j,mp1) = grid%evapprod(ips:ipe,k,j)
           end do
         end do
     end select
   end do

   end subroutine pack_dyn_vals

   subroutine set_dyn_vals

   use module_model_constants, only : g, t0, p1000mb, rcp

   integer :: ilp1, iup1, jlp1, jup1, klp1, kup1
   real    :: ginv, pinv

   select case( traject(trj,dm)%dyn_var(n) )
     case( 'p' )
       wrk4d(il,kl,jl,1) = chem(il,kl,jl,2) + chem(il,kl,jl,4)
       wrk4d(iu,kl,jl,1) = chem(iu,kl,jl,2) + chem(iu,kl,jl,4)
       wrk4d(il,ku,jl,1) = chem(il,ku,jl,2) + chem(il,ku,jl,4)
       wrk4d(iu,ku,jl,1) = chem(iu,ku,jl,2) + chem(iu,ku,jl,4)
       wrk4d(il,kl,ju,1) = chem(il,kl,ju,2) + chem(il,kl,ju,4)
       wrk4d(iu,kl,ju,1) = chem(iu,kl,ju,2) + chem(iu,kl,ju,4)
       wrk4d(il,ku,ju,1) = chem(il,ku,ju,2) + chem(il,ku,ju,4)
       wrk4d(iu,ku,ju,1) = chem(iu,ku,ju,2) + chem(iu,ku,ju,4)
     case( 'T' )
       wrk4d(il,kl,jl,1) = chem(il,kl,jl,2) + chem(il,kl,jl,4)
       wrk4d(iu,kl,jl,1) = chem(iu,kl,jl,2) + chem(iu,kl,jl,4)
       wrk4d(il,ku,jl,1) = chem(il,ku,jl,2) + chem(il,ku,jl,4)
       wrk4d(iu,ku,jl,1) = chem(iu,ku,jl,2) + chem(iu,ku,jl,4)
       wrk4d(il,kl,ju,1) = chem(il,kl,ju,2) + chem(il,kl,ju,4)
       wrk4d(iu,kl,ju,1) = chem(iu,kl,ju,2) + chem(iu,kl,ju,4)
       wrk4d(il,ku,ju,1) = chem(il,ku,ju,2) + chem(il,ku,ju,4)
       wrk4d(iu,ku,ju,1) = chem(iu,ku,ju,2) + chem(iu,ku,ju,4)

       pinv = 1./p1000mb
       wrk4d(il,kl,jl,1) = (chem(il,kl,jl,3) + t0)*(wrk4d(il,kl,jl,1)*pinv)**rcp
       wrk4d(iu,kl,jl,1) = (chem(iu,kl,jl,3) + t0)*(wrk4d(iu,kl,jl,1)*pinv)**rcp
       wrk4d(il,ku,jl,1) = (chem(il,ku,jl,3) + t0)*(wrk4d(il,ku,jl,1)*pinv)**rcp
       wrk4d(iu,ku,jl,1) = (chem(iu,ku,jl,3) + t0)*(wrk4d(iu,ku,jl,1)*pinv)**rcp
       wrk4d(il,kl,ju,1) = (chem(il,kl,ju,3) + t0)*(wrk4d(il,kl,ju,1)*pinv)**rcp
       wrk4d(iu,kl,ju,1) = (chem(iu,kl,ju,3) + t0)*(wrk4d(iu,kl,ju,1)*pinv)**rcp
       wrk4d(il,ku,ju,1) = (chem(il,ku,ju,3) + t0)*(wrk4d(il,ku,ju,1)*pinv)**rcp
       wrk4d(iu,ku,ju,1) = (chem(iu,ku,ju,3) + t0)*(wrk4d(iu,ku,ju,1)*pinv)**rcp
     case( 'z' )
       klp1 = kl + 1 ; kup1 = ku + 1
       ginv = 1./g
       wrk4d(il,kl,jl,1) = .5*(chem(il,kl,jl,8) + chem(il,klp1,jl,8) &
                               + chem(il,kl,jl,9) + chem(il,klp1,jl,9))*ginv
       wrk4d(iu,kl,jl,1) = .5*(chem(iu,kl,jl,8) + chem(iu,klp1,jl,8) &
                               + chem(iu,kl,jl,9) + chem(iu,klp1,jl,9))*ginv
       wrk4d(il,ku,jl,1) = .5*(chem(il,ku,jl,8) + chem(il,kup1,jl,8) &
                               + chem(il,ku,jl,9) + chem(il,kup1,jl,9))*ginv
       wrk4d(iu,ku,jl,1) = .5*(chem(iu,ku,jl,8) + chem(iu,kup1,jl,8) &
                               + chem(iu,ku,jl,9) + chem(iu,kup1,jl,9))*ginv
       wrk4d(il,kl,ju,1) = .5*(chem(il,kl,ju,8) + chem(il,klp1,ju,8) &
                               + chem(il,kl,ju,9) + chem(il,klp1,ju,9))*ginv
       wrk4d(iu,kl,ju,1) = .5*(chem(iu,kl,ju,8) + chem(iu,klp1,ju,8) &
                               + chem(iu,kl,ju,9) + chem(iu,klp1,ju,9))*ginv
       wrk4d(il,ku,ju,1) = .5*(chem(il,ku,ju,8) + chem(il,kup1,ju,8) &
                               + chem(il,ku,ju,9) + chem(il,kup1,ju,9))*ginv
       wrk4d(iu,ku,ju,1) = .5*(chem(iu,ku,ju,8) + chem(iu,kup1,ju,8) &
                               + chem(iu,ku,ju,9) + chem(iu,kup1,ju,9))*ginv
     case( 'u' )
       ilp1 = il + 1 ; iup1 = iu + 1
       wrk4d(il,kl,jl,1) = .5*(chem(il,kl,jl,5) + chem(ilp1,kl,jl,5))
       wrk4d(iu,kl,jl,1) = .5*(chem(iu,kl,jl,5) + chem(iup1,kl,jl,5))
       wrk4d(il,ku,jl,1) = .5*(chem(il,ku,jl,5) + chem(ilp1,ku,jl,5))
       wrk4d(iu,ku,jl,1) = .5*(chem(iu,ku,jl,5) + chem(iup1,ku,jl,5))
       wrk4d(il,kl,ju,1) = .5*(chem(il,kl,ju,5) + chem(ilp1,kl,ju,5))
       wrk4d(iu,kl,ju,1) = .5*(chem(iu,kl,ju,5) + chem(iup1,kl,ju,5))
       wrk4d(il,ku,ju,1) = .5*(chem(il,ku,ju,5) + chem(ilp1,ku,ju,5))
       wrk4d(iu,ku,ju,1) = .5*(chem(iu,ku,ju,5) + chem(iup1,ku,ju,5))
     case( 'v' )
       jlp1 = jl + 1 ; jup1 = ju + 1
       wrk4d(il,kl,jl,1) = .5*(chem(il,kl,jl,6) + chem(il,kl,jlp1,6))
       wrk4d(iu,kl,jl,1) = .5*(chem(iu,kl,jl,6) + chem(iu,kl,jlp1,6))
       wrk4d(il,ku,jl,1) = .5*(chem(il,ku,jl,6) + chem(il,ku,jlp1,6))
       wrk4d(iu,ku,jl,1) = .5*(chem(iu,ku,jl,6) + chem(iu,ku,jlp1,6))
       wrk4d(il,kl,ju,1) = .5*(chem(il,kl,ju,6) + chem(il,kl,jup1,6))
       wrk4d(iu,kl,ju,1) = .5*(chem(iu,kl,ju,6) + chem(iu,kl,jup1,6))
       wrk4d(il,ku,ju,1) = .5*(chem(il,ku,ju,6) + chem(il,ku,jup1,6))
       wrk4d(iu,ku,ju,1) = .5*(chem(iu,ku,ju,6) + chem(iu,ku,jup1,6))
     case( 'w' )
       klp1 = kl + 1 ; kup1 = ku + 1
       wrk4d(il,kl,jl,1) = .5*(chem(il,kl,jl,7) + chem(il,klp1,jl,7))
       wrk4d(iu,kl,jl,1) = .5*(chem(iu,kl,jl,7) + chem(iu,klp1,jl,7))
       wrk4d(il,ku,jl,1) = .5*(chem(il,ku,jl,7) + chem(il,kup1,jl,7))
       wrk4d(iu,ku,jl,1) = .5*(chem(iu,ku,jl,7) + chem(iu,kup1,jl,7))
       wrk4d(il,kl,ju,1) = .5*(chem(il,kl,ju,7) + chem(il,klp1,ju,7))
       wrk4d(iu,kl,ju,1) = .5*(chem(iu,kl,ju,7) + chem(iu,klp1,ju,7))
       wrk4d(il,ku,ju,1) = .5*(chem(il,ku,ju,7) + chem(il,kup1,ju,7))
       wrk4d(iu,ku,ju,1) = .5*(chem(iu,ku,ju,7) + chem(iu,kup1,ju,7))
     case( 'rainprod' )
       wrk4d(il,kl,jl,1) = chem(il,kl,jl,10)
       wrk4d(iu,kl,jl,1) = chem(iu,kl,jl,10)
       wrk4d(il,ku,jl,1) = chem(il,ku,jl,10)
       wrk4d(iu,ku,jl,1) = chem(iu,ku,jl,10)
       wrk4d(il,kl,ju,1) = chem(il,kl,ju,10)
       wrk4d(iu,kl,ju,1) = chem(iu,kl,ju,10)
       wrk4d(il,ku,ju,1) = chem(il,ku,ju,10)
       wrk4d(iu,ku,ju,1) = chem(iu,ku,ju,10)
     case( 'evapprod' )
       wrk4d(il,kl,jl,1) = chem(il,kl,jl,11)
       wrk4d(iu,kl,jl,1) = chem(iu,kl,jl,11)
       wrk4d(il,ku,jl,1) = chem(il,ku,jl,11)
       wrk4d(iu,ku,jl,1) = chem(iu,ku,jl,11)
       wrk4d(il,kl,ju,1) = chem(il,kl,ju,11)
       wrk4d(iu,kl,ju,1) = chem(iu,kl,ju,11)
       wrk4d(il,ku,ju,1) = chem(il,ku,ju,11)
       wrk4d(iu,ku,ju,1) = chem(iu,ku,ju,11)
   end select

   end subroutine set_dyn_vals

   subroutine set_msc_vals

   integer :: ilp1, iup1, jlp1, jup1, klp1, kup1

   select case( trim(St_Vars(St_Vars_ndx)%Stagger) )
     case( 'X' )
       ilp1 = il + 1 ; iup1 = iu + 1
       wrk4d(il,kl,jl,1) = .5*(chem(il,kl,jl,ndx) + chem(ilp1,kl,jl,ndx))
       wrk4d(iu,kl,jl,1) = .5*(chem(iu,kl,jl,ndx) + chem(iup1,kl,jl,ndx))
       wrk4d(il,ku,jl,1) = .5*(chem(il,ku,jl,ndx) + chem(ilp1,ku,jl,ndx))
       wrk4d(iu,ku,jl,1) = .5*(chem(iu,ku,jl,ndx) + chem(iup1,ku,jl,ndx))
       wrk4d(il,kl,ju,1) = .5*(chem(il,kl,ju,ndx) + chem(ilp1,kl,ju,ndx))
       wrk4d(iu,kl,ju,1) = .5*(chem(iu,kl,ju,ndx) + chem(iup1,kl,ju,ndx))
       wrk4d(il,ku,ju,1) = .5*(chem(il,ku,ju,ndx) + chem(ilp1,ku,ju,ndx))
       wrk4d(iu,ku,ju,1) = .5*(chem(iu,ku,ju,ndx) + chem(iup1,ku,ju,ndx))
     case( 'Y' )
       jlp1 = jl + 1 ; jup1 = ju + 1
       wrk4d(il,kl,jl,1) = .5*(chem(il,kl,jl,ndx) + chem(il,kl,jlp1,ndx))
       wrk4d(iu,kl,jl,1) = .5*(chem(iu,kl,jl,ndx) + chem(iu,kl,jlp1,ndx))
       wrk4d(il,ku,jl,1) = .5*(chem(il,ku,jl,ndx) + chem(il,ku,jlp1,ndx))
       wrk4d(iu,ku,jl,1) = .5*(chem(iu,ku,jl,ndx) + chem(iu,ku,jlp1,ndx))
       wrk4d(il,kl,ju,1) = .5*(chem(il,kl,ju,ndx) + chem(il,kl,jup1,ndx))
       wrk4d(iu,kl,ju,1) = .5*(chem(iu,kl,ju,ndx) + chem(iu,kl,jup1,ndx))
       wrk4d(il,ku,ju,1) = .5*(chem(il,ku,ju,ndx) + chem(il,ku,jup1,ndx))
       wrk4d(iu,ku,ju,1) = .5*(chem(iu,ku,ju,ndx) + chem(iu,ku,jup1,ndx))
     case( 'Z' )
       klp1 = kl + 1 ; kup1 = ku + 1
       wrk4d(il,kl,jl,1) = .5*(chem(il,kl,jl,ndx) + chem(il,klp1,jl,ndx))
       wrk4d(iu,kl,jl,1) = .5*(chem(iu,kl,jl,ndx) + chem(iu,klp1,jl,ndx))
       wrk4d(il,ku,jl,1) = .5*(chem(il,ku,jl,ndx) + chem(il,kup1,jl,ndx))
       wrk4d(iu,ku,jl,1) = .5*(chem(iu,ku,jl,ndx) + chem(iu,kup1,jl,ndx))
       wrk4d(il,kl,ju,1) = .5*(chem(il,kl,ju,ndx) + chem(il,klp1,ju,ndx))
       wrk4d(iu,kl,ju,1) = .5*(chem(iu,kl,ju,ndx) + chem(iu,klp1,ju,ndx))
       wrk4d(il,ku,ju,1) = .5*(chem(il,ku,ju,ndx) + chem(il,kup1,ju,ndx))
       wrk4d(iu,ku,ju,1) = .5*(chem(iu,ku,ju,ndx) + chem(iu,kup1,ju,ndx))
     case default
       wrk4d(il,kl,jl,1) = chem(il,kl,jl,ndx)
       wrk4d(iu,kl,jl,1) = chem(iu,kl,jl,ndx)
       wrk4d(il,ku,jl,1) = chem(il,ku,jl,ndx)
       wrk4d(iu,ku,jl,1) = chem(iu,ku,jl,ndx)
       wrk4d(il,kl,ju,1) = chem(il,kl,ju,ndx)
       wrk4d(iu,kl,ju,1) = chem(iu,kl,ju,ndx)
       wrk4d(il,ku,ju,1) = chem(il,ku,ju,ndx)
       wrk4d(iu,ku,ju,1) = chem(iu,ku,ju,ndx)
   end select

   end subroutine set_msc_vals

   end subroutine trajectory_driver

   subroutine trajectory_dchm_tstep_init( grid, is_chemstep )
!-----------------------------------------------------------------------------
!  initialize dchm buffer
!-----------------------------------------------------------------------------
   use module_domain, only : domain, get_ijk_from_grid
   use module_state_description, only : num_chem

!-----------------------------------------------------------------------------
!  dummy arguments
!-----------------------------------------------------------------------------
   logical, intent(in)      :: is_chemstep
   type(domain), intent(in) :: grid

!-----------------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------------
   integer :: astat
   integer :: j, k, m
   integer :: chm_ndx
   integer :: dm
   integer :: ims,ime, jms,jme, kms,kme
   integer :: ids,ide, jds,jde, kds,kde
   integer :: ips,ipe, jps,jpe, kps,kpe
   integer, pointer :: n_dchm       ! number dchm buffer species
   integer, pointer :: dchm_buf_ndx(:)
   character(len=256) :: err_mes

   do_chemstep = is_chemstep

   if( is_chemstep ) then
     dm = grid%id
     n_dchm => n_dchm_dm(dm)
     if( n_dchm > 0 ) then
       call get_ijk_from_grid( grid ,                   &
                               ids, ide, jds, jde, kds, kde,    &
                               ims, ime, jms, jme, kms, kme,    &
                               ips, ipe, jps, jpe, kps, kpe    )
       if( allocated( dchm_buff ) ) then
         deallocate( dchm_buff )
       endif
       allocate( dchm_buff(ims:ime,kms:kme,jms:jme,n_dchm+offset),stat=astat )
       if( astat /= 0 ) then
         write(err_mes,'(''trajectory_dchm_tstep_init('',i2.2,''): failed to allocate wrk4d: error = '',i6)') dm,astat
         call wrf_error_fatal( trim( err_mes  ) )
       endif
       dchm_buf_ndx => dchm_buf_ndx_dm(:,dm)
       do m = 1,n_dchm
         chm_ndx = dchm_buf_ndx(m)
         do j = jps,jpe
           do k = kps,kpe
             dchm_buff(ips:ipe,k,j,m+offset) = grid%chem(ips:ipe,k,j,chm_ndx)
           end do
         end do
       end do
     endif
   endif

   end subroutine trajectory_dchm_tstep_init

   subroutine trajectory_dchm_tstep_set( grid )
!-----------------------------------------------------------------------------
!  set dchm buffer
!-----------------------------------------------------------------------------
   use module_domain, only : domain, get_ijk_from_grid
   use module_state_description, only : num_chem

!-----------------------------------------------------------------------------
!  dummy arguments
!-----------------------------------------------------------------------------
   type(domain), intent(in) :: grid

!-----------------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------------
   integer :: j, k, m, mp1
   integer :: chm_ndx
   integer :: dm
   integer :: ims,ime, jms,jme, kms,kme
   integer :: ids,ide, jds,jde, kds,kde
   integer :: ips,ipe, jps,jpe, kps,kpe
   integer, pointer :: n_dchm       ! number dchm buffer species
   integer, pointer :: dchm_buf_ndx(:)

   dm = grid%id
   n_dchm => n_dchm_dm(dm)
   if( n_dchm > 0 ) then
     call get_ijk_from_grid( grid ,                   &
                             ids, ide, jds, jde, kds, kde,    &
                             ims, ime, jms, jme, kms, kme,    &
                             ips, ipe, jps, jpe, kps, kpe    )
     dchm_buf_ndx => dchm_buf_ndx_dm(:,dm)
     do m = 1,n_dchm
       mp1 = m + offset
       chm_ndx = dchm_buf_ndx(m)
       do j = jps,jpe
         do k = kps,kpe
           dchm_buff(ips:ipe,k,j,mp1) = grid%chem(ips:ipe,k,j,chm_ndx) - dchm_buff(ips:ipe,k,j,mp1)
         end do
       end do
     end do
   endif

   end subroutine trajectory_dchm_tstep_set

   subroutine trajectory_create_file( grid, n_traj )
!-----------------------------------------------------------------------------
!  create trajectory netcdf file
!-----------------------------------------------------------------------------
   use module_domain
   use module_state_description, only : param_first_scalar, num_chem, num_moist, num_tracer
   use module_scalar_tables,     only : chem_dname_table, moist_dname_table, tracer_dname_table

!-----------------------------------------------------------------------------
!  dummy arguments
!-----------------------------------------------------------------------------
   integer, intent(in)       :: n_traj
   type(domain), intent(in)  :: grid

!-----------------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------------
   integer :: dm
   integer :: ncid, ios
   integer :: traj_dim, time_dim, Times_dim
   integer :: varid
   integer, pointer :: num_msc
   integer :: var_dims(2)
   integer :: m, n, trj, pkg
   character(len=10)  :: coord_name(5) = (/ 'traj_i    ', 'traj_j    ', 'traj_k    ', &
                                            'traj_long ', 'traj_lat  ' /)
   character(len=256) :: filename
   character(len=256) :: var_name
   character(len=256) :: err_mes
   character(len=256) :: description
   character(len=256) :: units

   logical, external :: wrf_dm_on_monitor

include 'netcdf.inc'

master_proc: &
   if( wrf_dm_on_monitor() ) then
     dm = grid%id
     write(filename,'(''wrfout_traj_d'',i2.2)',iostat=ios) dm
     if( ios /= 0 ) then
       write(err_mes,'(''trajectory_create_file: failed to set filename: error = '',i6)') ios
       call wrf_error_fatal( trim( err_mes  ) )
     endif
!    ios = nf_create( trim(filename), or(nf_clobber,nf_netcdf4), ncid )
     ios = nf_create( trim(filename), nf_clobber, ncid )
     if( ios /= nf_noerr ) then
       write(err_mes,'(''trajectory_create_file: failed to create '',a,'': error = '',i6)') trim(filename),ios
       call wrf_error_fatal( trim( err_mes  ) )
     endif
!-----------------------------------------------------------------------------
!  define dimensions
!-----------------------------------------------------------------------------
     err_mes = 'trajectory_create_file: failed to create traj dimension'
     call handle_ncerr( nf_def_dim( ncid, 'traj', n_traj, traj_dim ), trim(err_mes) )
     err_mes = 'trajectory_create_file: failed to create time dimension'
     call handle_ncerr( nf_def_dim( ncid, 'time', nf_unlimited, time_dim ), trim(err_mes) )
     err_mes = 'trajectory_create_file: failed to create Times dimension'
     call handle_ncerr( nf_def_dim( ncid, 'DateStrLen', 19, Times_dim ), trim(err_mes) )
!-----------------------------------------------------------------------------
!  define variables
!-----------------------------------------------------------------------------
     var_dims(:) = (/ Times_dim,time_dim /)
     err_mes = 'trajectory_create_file: failed to create Times variable'
     call handle_ncerr( nf_def_var( ncid, 'Times', nf_char, 2, var_dims, varid ), trim(err_mes) )

!-----------------------------------------------------------------------------
!  first the coordinate variables
!-----------------------------------------------------------------------------
     var_dims(:) = (/ traj_dim,time_dim /)
     do m = 1,5
       var_name = coord_name(m)
       err_mes = 'trajectory_create_file: failed to create ' // trim(var_name) // ' variable'
       call handle_ncerr( nf_def_var( ncid, trim(var_name), nf_real, 2, var_dims, varid ), trim(err_mes) )
     end do
!-----------------------------------------------------------------------------
!  then the species variables
!-----------------------------------------------------------------------------
     num_msc => num_msc_dm(dm)
pgk_loop: &
     do pkg = 1,pkg_max
       select case( pkg_tag(pkg) )
         case( 'chm' )
           trj_msk => trj_msk_dm(:,:,chm_pkg,dm)
           if( any( trj_msk(:n_traj,1) ) ) then
             call def_vars( 'chm' )
           endif
         case( 'hyd' )
           trj_msk => trj_msk_dm(:,:,hyd_pkg,dm)
           if( any( trj_msk(:n_traj,1) ) ) then
             call def_vars( 'hyd' )
           endif
         case( 'trc' )
           trj_msk => trj_msk_dm(:,:,trc_pkg,dm)
           if( any( trj_msk(:n_traj,1) ) ) then
             call def_vars( 'trc' )
           endif
         case( 'dyn' )
           trj_msk => trj_msk_dm(:,:,dyn_pkg,dm)
           if( any( trj_msk(:n_traj,1) ) ) then
             call def_vars( 'dyn' )
           endif
         case( 'msc' )
           trj_msk => trj_msk_dm(:,:,msc_pkg,dm)
           if( any( trj_msk(:n_traj,1) ) ) then
             call def_vars( 'msc' )
           endif
         case( 'dchm' )
           trj_msk => trj_msk_dm(:,:,dchm_pkg,dm)
           if( any( trj_msk(:n_traj,1) ) ) then
             call def_vars( 'dchm' )
           endif
       end select
     end do pgk_loop

     err_mes = 'trajectory_create_file: failed to end definition for file ' // trim(filename)
     call handle_ncerr( nf_enddef( ncid ), trim(err_mes) )
     err_mes = 'trajectory_create_file: failed to close file ' // trim(filename)
     call handle_ncerr( nf_close( ncid ), trim(err_mes) )
   endif master_proc

   CONTAINS

   subroutine def_vars( var_type )

   character(len=*), intent(in)  :: var_type

   integer :: m, ndx, trj
   integer, pointer  :: n_dchm
   character(len=32) :: spc_name

   select case( var_type )
     case( 'chm' )
       trj_msk => trj_msk_dm(:,:,chm_pkg,dm)
       do n = param_first_scalar,num_chem
         if( any( trj_msk(:n_traj,n) ) ) then
           spc_name = chem_dname_table(dm,n)
           write(var_name,'(a,''_traj'')') trim(spc_name)
           err_mes = 'def_vars: failed to create ' // trim(var_name) // ' variable'
           call handle_ncerr( nf_def_var( ncid, trim(var_name), nf_real, 2, var_dims, varid ), trim(err_mes) )
           description = trim(var_name) // ' mixing ratio'
           call handle_ncerr( nf_put_att_text( ncid, varid, 'description', len_trim(description), description ), trim(err_mes) )
           units = 'ppmv'
trj_loop:  do trj = 1,n_traj
             if( trj_msk(trj,n) ) then
               do m = 1,trjects(trj)%n_chm_var
                 if( trim(trjects(trj)%chm_spc(m)) == trim(spc_name) ) then
                   if( .not. trjects(trj)%chm_is_gas(m) ) then
                     units = 'ug/kg-dryair'
                   endif
                   exit trj_loop
                 endif
               end do
               exit trj_loop
             endif
           end do trj_loop
           call handle_ncerr( nf_put_att_text( ncid, varid, 'units', len_trim(units), units ), trim(err_mes) )
         endif
       end do
     case( 'hyd' )
       trj_msk => trj_msk_dm(:,:,hyd_pkg,dm)
       do n = param_first_scalar,num_moist
         if( any( trj_msk(:n_traj,n) ) ) then
           write(var_name,'(a,''_traj'')') trim(moist_dname_table(dm,n))
           err_mes = 'def_vars: failed to create ' // trim(var_name) // ' variable'
           call handle_ncerr( nf_def_var( ncid, trim(var_name), nf_real, 2, var_dims, varid ), trim(err_mes) )
           description = trim(var_name) // ' mixing ratio'
           call handle_ncerr( nf_put_att_text( ncid, varid, 'description', len_trim(description), description ), trim(err_mes) )
           units = 'ug/kg-dryair'
           call handle_ncerr( nf_put_att_text( ncid, varid, 'units', len_trim(units), units ), trim(err_mes) )
         endif
       end do
     case( 'trc' )
       trj_msk => trj_msk_dm(:,:,trc_pkg,dm)
       do n = param_first_scalar,num_tracer
         if( any( trj_msk(:n_traj,n) ) ) then
           write(var_name,'(a,''_traj'')') trim(tracer_dname_table(dm,n))
           err_mes = 'def_vars: failed to create ' // trim(var_name) // ' variable'
           call handle_ncerr( nf_def_var( ncid, trim(var_name), nf_real, 2, var_dims, varid ), trim(err_mes) )
           description = trim(var_name) // ' mixing ratio'
           call handle_ncerr( nf_put_att_text( ncid, varid, 'description', len_trim(description), description ), trim(err_mes) )
           units = ' '
           call handle_ncerr( nf_put_att_text( ncid, varid, 'units', len_trim(units), units ), trim(err_mes) )
         endif
       end do
     case( 'dyn' )
       trj_msk => trj_msk_dm(:,:,dyn_pkg,dm)
       do n = param_first_scalar,dyn_max + offset
         if( any( trj_msk(:n_traj,n) ) ) then
           write(var_name,'(a,''_traj'')') trim(dyn_var_lst(n-1))
           err_mes = 'def_vars: failed to create ' // trim(var_name) // ' variable'
           call handle_ncerr( nf_def_var( ncid, trim(var_name), nf_real, 2, var_dims, varid ), trim(err_mes) )
           description = trim(dyn_var_desc_att(n-1))
           call handle_ncerr( nf_put_att_text( ncid, varid, 'description', len_trim(description), description ), trim(err_mes) )
           units = trim(dyn_var_unit_att(n-1))
           call handle_ncerr( nf_put_att_text( ncid, varid, 'units', len_trim(units), units ), trim(err_mes) )
         endif
       end do
     case( 'msc' )
       trj_msk => trj_msk_dm(:,:,msc_pkg,dm)
       do n = param_first_scalar,num_msc + offset
         if( any( trj_msk(:n_traj,n) ) ) then
           write(var_name,'(a,''_traj'')') trim(St_Vars(n-1)%Varname)
           err_mes = 'def_vars: failed to create ' // trim(var_name) // ' variable'
           call handle_ncerr( nf_def_var( ncid, trim(var_name), nf_real, 2, var_dims, varid ), trim(err_mes) )
           description = trim(St_Vars(n-1)%Description)
           call handle_ncerr( nf_put_att_text( ncid, varid, 'description', len_trim(description), description ), trim(err_mes) )
           units = trim(St_Vars(n-1)%Units)
           call handle_ncerr( nf_put_att_text( ncid, varid, 'units', len_trim(units), units ), trim(err_mes) )
         endif
       end do
     case( 'dchm' )
       trj_msk => trj_msk_dm(:,:,dchm_pkg,dm)
       n_dchm  => n_dchm_dm(dm)
       dchm_buf_ndx => dchm_buf_ndx_dm(:,dm)
       do n = 1,n_dchm
         if( any( trj_msk(:n_traj,n+offset) ) ) then
           ndx = dchm_buf_ndx(n)
           spc_name = 'dchm_' // trim(chem_dname_table(dm,ndx))
           write(var_name,'(a,''_traj'')') trim(spc_name)
           err_mes = 'def_vars: failed to create ' // trim(var_name) // ' variable'
           call handle_ncerr( nf_def_var( ncid, trim(var_name), nf_real, 2, var_dims, varid ), trim(err_mes) )
           description = trim(var_name) // ' mixing ratio'
           call handle_ncerr( nf_put_att_text( ncid, varid, 'description', len_trim(description), description ), trim(err_mes) )
           units = 'ppmv'
           call handle_ncerr( nf_put_att_text( ncid, varid, 'units', len_trim(units), units ), trim(err_mes) )
         endif
       end do
   end select

   end subroutine def_vars

   end subroutine trajectory_create_file

   subroutine trajectory_write_file( n_traj, n_vals, dm )
!-----------------------------------------------------------------------------
!  create trajectory netcdf file
!-----------------------------------------------------------------------------
   use module_domain
   use module_state_description, only : param_first_scalar, num_chem, num_moist, num_tracer
   use module_scalar_tables,     only : chem_dname_table, moist_dname_table, tracer_dname_table

!-----------------------------------------------------------------------------
!  dummy arguments
!-----------------------------------------------------------------------------
   integer, intent(in)        :: n_traj
   integer, intent(inout)     :: n_vals
   integer, intent(in)        :: dm

!-----------------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------------
   integer :: ncid
   integer :: astat, ios
   integer :: time_id
   integer :: varid
   integer :: l, m, n, trj, pkg, spc, spcp1
   integer :: time_ndx
   integer :: buf_ndx
   integer :: ndx
   integer, pointer :: num_msc      ! number misc species
   integer, pointer :: n_dchm       ! total number of dchm species in domain
   real, allocatable :: holder(:,:)
   character(len=10)  :: coord_name(5) = (/ 'traj_i    ', 'traj_j    ', 'traj_k    ', &
                                            'traj_long ', 'traj_lat  ' /)
   character(len=256) :: var_name
   character(len=256) :: err_mes
   character(len=256) :: filename
   character(len=256) :: spcname

   logical :: found

include 'netcdf.inc'

!---------------------------------------------------------------------
!  open netcdf file
!---------------------------------------------------------------------
   write(filename,'(''wrfout_traj_d'',i2.2)',iostat=ios) dm
   if( ios /= 0 ) then
     write(err_mes,'(''trajectory_write_file: failed to set filename: error = '',i6)') ios
     call wrf_error_fatal( trim( err_mes  ) )
   endif
   ios = nf_open( trim(filename), nf_write, ncid )
   if( ios /= 0 ) then
     write(err_mes,'(''trajectory_write_file: failed to open '',a,'': error = '',i6)') trim(filename),ios
     call wrf_error_fatal( trim( err_mes  ) )
   endif

!---------------------------------------------------------------------
!  allocate wrking array
!---------------------------------------------------------------------
   allocate( holder(n_traj,n_vals),stat=astat )
   if( astat /= 0 ) then
     write(err_mes,'(''trajectory_write_file: failed to allocate holder; error = '',i6)') astat
     call wrf_error_fatal( trim( err_mes  ) )
   endif

!---------------------------------------------------------------------
!  get current time index
!---------------------------------------------------------------------
   err_mes = 'trajectory_write_file: failed to get time id'
   call handle_ncerr( nf_inq_dimid( ncid, 'time', time_id ),trim(err_mes) )
   err_mes = 'trajectory_write_file: failed to get time dimension'
   call handle_ncerr( nf_inq_dimlen( ncid, time_id, time_ndx ),trim(err_mes) )
   time_ndx = time_ndx + 1

!---------------------------------------------------------------------
!  write out trajectory times
!---------------------------------------------------------------------
   err_mes = 'trajectory_write_file: failed to get Times id'
   call handle_ncerr( nf_inq_varid( ncid, 'Times', varid ),trim(err_mes) )
   err_mes = 'trajectory_write_file: failed to write Times'
   call handle_ncerr( nf_put_vara_text( ncid, varid, (/ 1,time_ndx /), (/ 19,n_vals /), trj_pbf(1)%times(:n_vals) ), trim(err_mes) )

!---------------------------------------------------------------------
!  write out trajectory coordinates
!---------------------------------------------------------------------
coord_loop: &
   do l = 1,5
     var_name = coord_name(l)
     err_mes = 'trajectory_write_file: failed to get '// trim(var_name) // ' id'
     call handle_ncerr( nf_inq_varid( ncid, trim(var_name), varid ),trim(err_mes) )
     select case( l )
       case( 1 )
         do n = 1,n_traj
           holder(n,:n_vals) = trj_pbf(n)%trj_i(:n_vals)
         end do
       case( 2 )
         do n = 1,n_traj
           holder(n,:n_vals) = trj_pbf(n)%trj_j(:n_vals)
         end do
       case( 3 )
         do n = 1,n_traj
           holder(n,:n_vals) = trj_pbf(n)%trj_k(:n_vals)
         end do
       case( 4 )
         do n = 1,n_traj
           holder(n,:n_vals) = trj_pbf(n)%trj_lons(:n_vals)
         end do
       case( 5 )
         do n = 1,n_traj
           holder(n,:n_vals) = trj_pbf(n)%trj_lats(:n_vals)
         end do
       end select
       err_mes = 'trajectory_write_file: failed to write ' // trim(var_name)
       call handle_ncerr( nf_put_vara_real( ncid, varid, (/ 1,time_ndx /), (/ n_traj,n_vals /), &
                                            holder ), trim(err_mes) )
   end do coord_loop

   St_Vars => St_Vars_dm(:,dm)
   St_Vars_msk => St_Vars_msk_dm(:,dm)
   num_msc => num_msc_dm(dm)
!---------------------------------------------------------------------
!  write out trajectory variables
!---------------------------------------------------------------------
pkg_loop: &
   do pkg = 1,pkg_max
     select case( pkg_tag(pkg) )
       case( 'chm' )
         trj_msk => trj_msk_dm(:,:,chm_pkg,dm)
         do spc = param_first_scalar,num_chem
           if( any( trj_msk(:n_traj,spc) ) ) then
             holder(:,:) = missing_val
             do trj = 1,n_traj
               if( trj_msk(trj,spc) ) then
                 buf_ndx = get_spc_buf_ndx( trjects(trj)%n_chm_var, trjects(trj)%chm_spc, chem_dname_table(dm,spc) )
                 if( buf_ndx > 0 ) then
                   holder(trj,:n_vals) = trj_pbf(trj)%chm_vals(:n_vals,buf_ndx)
                 endif
               endif
             end do
             write(var_name,'(a,''_traj'')') trim(chem_dname_table(dm,spc))
             err_mes = 'trajectory_write_file: failed to get '// trim(var_name) // ' id'
             call handle_ncerr( nf_inq_varid( ncid, trim(var_name), varid ),trim(err_mes) )
             err_mes = 'trajectory_write_file: failed to write ' // trim(var_name)
             call handle_ncerr( nf_put_vara_real( ncid, varid, (/ 1,time_ndx /), (/ n_traj,n_vals /), &
                                                  holder ),trim(err_mes) )
           endif
         end do
       case( 'hyd' )
         trj_msk => trj_msk_dm(:,:,hyd_pkg,dm)
         do spc = param_first_scalar,num_moist
           if( any( trj_msk(:n_traj,spc) ) ) then
             holder(:,:) = missing_val
             do trj = 1,n_traj
               if( trj_msk(trj,spc) ) then
                 buf_ndx = get_spc_buf_ndx( trjects(trj)%n_hyd_var, trjects(trj)%hyd_spc, moist_dname_table(dm,spc) )
                 if( buf_ndx > 0 ) then
                   holder(trj,:n_vals) = trj_pbf(trj)%hyd_vals(:n_vals,buf_ndx)
                 endif
               endif
             end do
             write(var_name,'(a,''_traj'')') trim(moist_dname_table(dm,spc))
             err_mes = 'trajectory_write_file: failed to get '// trim(var_name) // ' id'
             call handle_ncerr( nf_inq_varid( ncid, trim(var_name), varid ),trim(err_mes) )
             err_mes = 'trajectory_write_file: failed to write ' // trim(var_name)
             call handle_ncerr( nf_put_vara_real( ncid, varid, (/ 1,time_ndx /), (/ n_traj,n_vals /), &
                                                  holder ),trim(err_mes) )
           endif
         end do
       case( 'trc' )
         trj_msk => trj_msk_dm(:,:,trc_pkg,dm)
         do spc = param_first_scalar,num_tracer
           if( any( trj_msk(:n_traj,spc) ) ) then
             holder(:,:) = missing_val
             do trj = 1,n_traj
               if( trj_msk(trj,spc) ) then
                 buf_ndx = get_spc_buf_ndx( trjects(trj)%n_trc_var, trjects(trj)%trc_spc, tracer_dname_table(dm,spc) )
                 if( buf_ndx > 0 ) then
                   holder(trj,:n_vals) = trj_pbf(trj)%trc_vals(:n_vals,buf_ndx)
                 endif
               endif
             end do
             write(var_name,'(a,''_traj'')') trim(tracer_dname_table(dm,spc))
             err_mes = 'trajectory_write_file: failed to get '// trim(var_name) // ' id'
             call handle_ncerr( nf_inq_varid( ncid, trim(var_name), varid ),trim(err_mes) )
             err_mes = 'trajectory_write_file: failed to write ' // trim(var_name)
             call handle_ncerr( nf_put_vara_real( ncid, varid, (/ 1,time_ndx /), (/ n_traj,n_vals /), &
                                                  holder ),trim(err_mes) )
           endif
         end do
       case( 'dyn' )
         trj_msk => trj_msk_dm(:,:,dyn_pkg,dm)
         do spc = param_first_scalar,dyn_max+offset
           if( any( trj_msk(:n_traj,spc) ) ) then
             holder(:,:) = missing_val
             do trj = 1,n_traj
               if( trj_msk(trj,spc) ) then
                 buf_ndx = get_spc_buf_ndx( trjects(trj)%n_dyn_var, trjects(trj)%dyn_var, dyn_var_lst(spc-offset) )
                 if( buf_ndx > 0 ) then
                   holder(trj,:n_vals) = trj_pbf(trj)%dyn_vals(:n_vals,buf_ndx)
                 endif
               endif
             end do
             write(var_name,'(a,''_traj'')') trim(dyn_var_lst(spc-offset))
             err_mes = 'trajectory_write_file: failed to get '// trim(var_name) // ' id'
             call handle_ncerr( nf_inq_varid( ncid, trim(var_name), varid ),trim(err_mes) )
             err_mes = 'trajectory_write_file: failed to write ' // trim(var_name)
             call handle_ncerr( nf_put_vara_real( ncid, varid, (/ 1,time_ndx /), (/ n_traj,n_vals /), &
                                                  holder ),trim(err_mes) )
           endif
         end do
       case( 'msc' )
         trj_msk => trj_msk_dm(:,:,msc_pkg,dm)
         do spc = param_first_scalar,num_msc+offset
           if( any( trj_msk(:n_traj,spc) ) ) then
             holder(:,:) = missing_val
             do trj = 1,n_traj
               if( trj_msk(trj,spc) ) then
                 found = .false.
                 do buf_ndx = 1,traject(trj,dm)%n_msc_var
                   if( traject(trj,dm)%msc_ndx(buf_ndx) == spc ) then
                     found = .true.
                     exit
                   endif
                 end do
                 if( found ) then
                   holder(trj,:n_vals) = trj_pbf(trj)%msc_vals(:n_vals,buf_ndx)
                 endif
               endif
             end do
             write(var_name,'(a,''_traj'')') trim(St_Vars(spc-offset)%Varname)
             err_mes = 'trajectory_write_file: failed to get '// trim(var_name) // ' id'
             call handle_ncerr( nf_inq_varid( ncid, trim(var_name), varid ),trim(err_mes) )
             err_mes = 'trajectory_write_file: failed to write ' // trim(var_name)
             call handle_ncerr( nf_put_vara_real( ncid, varid, (/ 1,time_ndx /), (/ n_traj,n_vals /), &
                                                  holder ),trim(err_mes) )
           endif
         end do
       case( 'dchm' )
         dchm_buf_ndx => dchm_buf_ndx_dm(:,dm)
         trj_msk => trj_msk_dm(:,:,dchm_pkg,dm)
         n_dchm  => n_dchm_dm(dm)
         do spc = 1,n_dchm
           spcp1 = spc + 1
           if( any( trj_msk(:n_traj,spcp1) ) ) then
             holder(:,:) = missing_val
             do trj = 1,n_traj
               if( trj_msk(trj,spcp1) ) then
                 buf_ndx = get_dchm_buf_ndx( trjects(trj)%n_dchm_var, trjects(trj)%dchm_ndx, spcp1 )
                 if( buf_ndx > 0 ) then
                   holder(trj,:n_vals) = trj_pbf(trj)%dchm_vals(:n_vals,buf_ndx)
                 endif
               endif
             end do
             ndx = dchm_buf_ndx(spc)
             var_name = 'dchm_' // trim(chem_dname_table(dm,ndx)) // '_traj'
!            write(var_name,'(a,''_traj'')') trim(chem_dname_table(dm,ndx))
             err_mes = 'trajectory_write_file: failed to get '// trim(var_name) // ' id'
             call handle_ncerr( nf_inq_varid( ncid, trim(var_name), varid ),trim(err_mes) )
             err_mes = 'trajectory_write_file: failed to write ' // trim(var_name)
             call handle_ncerr( nf_put_vara_real( ncid, varid, (/ 1,time_ndx /), (/ n_traj,n_vals /), &
                                                  holder ),trim(err_mes) )
           endif
         end do
     end select
   end do pkg_loop

   n_vals = 0

   ios = nf_close( ncid )

   if( allocated( holder ) ) then
     deallocate( holder )
   endif

   end subroutine trajectory_write_file

   integer function get_spc_buf_ndx( ncnt, list, match_name )

   integer, intent(in)          :: ncnt
   character(len=*), intent(in) :: match_name
   character(len=*), intent(in) :: list(:)

   integer :: spc

   get_spc_buf_ndx = -1
   do spc = 1,ncnt
     if( trim(match_name) == trim(list(spc)) ) then
       get_spc_buf_ndx = spc
       exit
     endif
   end do

   end function get_spc_buf_ndx

   integer function get_dchm_buf_ndx( ncnt, list, match_ndx )

   integer, intent(in)          :: ncnt
   integer, intent(in)          :: match_ndx
   integer, intent(in)          :: list(:)

   integer :: spc

   get_dchm_buf_ndx = -1
   do spc = 1,ncnt
     if( match_ndx == list(spc) ) then
       get_dchm_buf_ndx = spc
       exit
     endif
   end do

   end function get_dchm_buf_ndx

   subroutine handle_ncerr( ret, mes )
!---------------------------------------------------------------------
!	... netcdf error handling routine
!---------------------------------------------------------------------
!---------------------------------------------------------------------
!	... dummy arguments
!---------------------------------------------------------------------
   integer, intent(in) :: ret
   character(len=*), intent(in) :: mes

include 'netcdf.inc'

   if( ret /= nf_noerr ) then
      call wrf_message( trim(mes) )
      call wrf_message( trim(nf_strerror(ret)) )
      call wrf_abort
   endif

   end subroutine handle_ncerr
#endif

   subroutine trajmapproj (grid,config_flags,ts_proj)

   use module_domain
   use module_llxy
   use module_configure, only : grid_config_rec_type, model_config_rec
   use module_dm, only : wrf_dm_min_real

   IMPLICIT NONE


!------------------------------------------------------------------------
! Arguments
!------------------------------------------------------------------------
   TYPE(domain), INTENT(IN) :: grid
   TYPE(grid_config_rec_type) , INTENT(IN)  :: config_flags
   TYPE(PROJ_INFO), INTENT(out) :: ts_proj

!------------------------------------------------------------------------
! Local variables
!------------------------------------------------------------------------
   REAL :: ts_rx, ts_ry, ts_xlat, ts_xlong, ts_hgt
   REAL :: known_lat, known_lon

   INTEGER :: ids, ide, jds, jde, kds, kde,        &
              ims, ime, jms, jme, kms, kme,        &
              ips, ipe, jps, jpe, kps, kpe

   TYPE (grid_config_rec_type)               :: config_flags_temp

   config_flags_temp = config_flags

   call get_ijk_from_grid ( grid ,                               &
                            ids, ide, jds, jde, kds, kde,        &
                            ims, ime, jms, jme, kms, kme,        &
                            ips, ipe, jps, jpe, kps, kpe )

   call model_to_grid_config_rec ( grid%id , model_config_rec , config_flags_temp )


!------------------------------------------------------------------------
! Set up map transformation structure
!------------------------------------------------------------------------
   call map_init( ts_proj )

   IF (ips <= 1 .AND. 1 <= ipe .AND. jps <= 1 .AND. 1 <= jpe) THEN
      known_lat = grid%xlat(1,1)
      known_lon = grid%xlong(1,1)
   ELSE
      known_lat = 9999.
      known_lon = 9999.
   END IF
   known_lat = wrf_dm_min_real(known_lat)
   known_lon = wrf_dm_min_real(known_lon)


   select case( config_flags%map_proj )
!------------------------------------------------------------------------
! Mercator
!------------------------------------------------------------------------
     case( PROJ_MERC )
       call map_set(PROJ_MERC, ts_proj,               &
                    truelat1 = config_flags%truelat1, &
                    lat1     = known_lat,             &
                    lon1     = known_lon,             &
                    knowni   = 1.,                    &
                    knownj   = 1.,                    &
                    dx       = config_flags%dx)
!------------------------------------------------------------------------
! Lambert conformal
!------------------------------------------------------------------------
     case( PROJ_LC )
       call map_set(PROJ_LC, ts_proj,                  &
                    truelat1 = config_flags%truelat1,  &
                    truelat2 = config_flags%truelat2,  &
                    stdlon   = config_flags%stand_lon, &
                    lat1     = known_lat,              &
                    lon1     = known_lon,              &
                    knowni   = 1.,                     &
                    knownj   = 1.,                     &
                    dx       = config_flags%dx)
!------------------------------------------------------------------------
! Polar stereographic
!------------------------------------------------------------------------
     case( PROJ_PS )
       call map_set(PROJ_PS, ts_proj,                  &
                    truelat1 = config_flags%truelat1,  &
                    stdlon   = config_flags%stand_lon, &
                    lat1     = known_lat,              &
                    lon1     = known_lon,              &
                    knowni   = 1.,                     &
                    knownj   = 1.,                     &
                    dx       = config_flags%dx)
!------------------------------------------------------------------------
! Cassini (global ARW)
!------------------------------------------------------------------------
     case( PROJ_CASSINI )
       call map_set(PROJ_CASSINI, ts_proj,                            &
                    latinc   = grid%dy*360.0/(2.0*EARTH_RADIUS_M*PI), &
                    loninc   = grid%dx*360.0/(2.0*EARTH_RADIUS_M*PI), &
                    lat1     = known_lat,                             &
                    lon1     = known_lon,                             &
! We still need to get POLE_LAT and POLE_LON metadata variables before
!   this will work for rotated poles.
                    lat0     = 90.0,                                  &
                    lon0     = 0.0,                                   &
                    knowni   = 1.,                                    &
                    knownj   = 1.,                                    &
                    stdlon   = config_flags%stand_lon)
!------------------------------------------------------------------------
! Rotated latitude-longitude
!------------------------------------------------------------------------
     case( PROJ_ROTLL )
       call map_set(PROJ_ROTLL, ts_proj,                      &
                    ixdim    = grid%e_we-1,                   &
                    jydim    = grid%e_sn-1,                   &
                    phi      = real(grid%e_sn-2)*grid%dy/2.0, &
                    lambda   = real(grid%e_we-2)*grid%dx,     &
                    lat1     = config_flags%cen_lat,          &
                    lon1     = config_flags%cen_lon,          &
                    latinc   = grid%dy,                       &
                    loninc   = grid%dx,                       &
                    stagger  = HH)
   end select

   end subroutine trajmapproj

   end module module_trajectory
