share/module_trajectory.F

   1
   2    module module_trajectory
   3
   4    use module_driver_constants,  only : max_domains
   5 #if( EM_CORE == 1 )
   6    use module_state_description, only : num_chem
   7 #endif
   8
   9    implicit none
  10
  11    private
  12    public :: trajectory_init
  13 #ifdef NETCDF
  14    public :: trajectory_driver
  15    public :: trajectory_dchm_tstep_init
  16    public :: trajectory_dchm_tstep_set
  17 #endif
  18    public :: traject
  19    public :: traj_cnt
  20
  21    integer, parameter :: vals_max = 1000
  22    integer, parameter :: traj_max = 1000
  23    integer, parameter :: var_max  = 100
  24    integer, parameter :: pkg_max  = 6
  25    integer, parameter :: dyn_max  = 8
  26    integer, parameter :: chm_pkg  = 1
  27    integer, parameter :: hyd_pkg  = 2
  28    integer, parameter :: trc_pkg  = 3
  29    integer, parameter :: dyn_pkg  = 4
  30    integer, parameter :: msc_pkg  = 5
  31    integer, parameter :: dchm_pkg = 6
  32    real, parameter    :: missing_val = -9999.
  33    real, parameter    :: zero_val    = 0.
  34
  35    integer :: n_dom        ! domain count
  36    integer, pointer :: n_vals                         ! wrking count of time points in buffer
  37    integer, target  :: n_vals_dm(max_domains)         ! count of time points in buffer per domain
  38    integer, target  :: n_dchm_dm(max_domains)         ! count of dchm variables in  buffer per domain
  39    integer, allocatable, target :: num_msc_dm(:)      ! number misc species
  40    integer, allocatable, target :: dchm_buf_ndx_dm(:,:)  ! dchm buffer indices
  41    integer, pointer :: dchm_buf_ndx(:)
  42    integer :: offset       ! into variable arrays
  43
  44    type default
  45      character(len=19) :: start_time
  46      character(len=19) :: stop_time
  47      character(len=32) :: chm_name(var_max)
  48      character(len=32) :: hyd_name(var_max)
  49      character(len=32) :: trc_name(var_max)
  50      character(len=32) :: dyn_name(var_max)
  51      character(len=32) :: msc_name(var_max)
  52      character(len=32) :: dchm_name(var_max)
  53    end type default
  54
  55    type base
  56      real    :: lon
  57      real    :: lat
  58      real    :: lev
  59      real    :: x, y
  60      real    :: traj_var(var_max)
  61      integer :: n_chm_var
  62      integer :: n_ct_var
  63      integer :: n_hyd_var
  64      integer :: n_trc_var
  65      integer :: n_dyn_var
  66      integer :: n_msc_var
  67      integer :: n_dchm_var
  68      integer :: chm_ndx(var_max)
  69      integer :: hyd_ndx(var_max)
  70      integer :: trc_ndx(var_max)
  71      integer :: dyn_ndx(var_max)
  72      integer :: msc_ndx(var_max)
  73      integer :: dchm_ndx(var_max)
  74      character(len=19) :: start_time
  75      character(len=19) :: stop_time
  76      character(len=32) :: chm_spc(var_max)
  77      character(len=32) :: hyd_spc(var_max)
  78      character(len=32) :: trc_spc(var_max)
  79      character(len=32) :: dyn_var(var_max)
  80      character(len=32) :: msc_var(var_max)
  81      character(len=32) :: dchm_spc(var_max)
  82      logical :: in_dom
  83      logical :: in_patch
  84      logical :: is_stationary
  85      logical :: z_is_agl
  86      logical :: chm_is_gas(var_max)
  87    end type base
  88
  89    type buffer
  90      real    :: trj_i(vals_max)
  91      real    :: trj_j(vals_max)
  92      real    :: trj_k(vals_max)
  93      real    :: trj_lons(vals_max)
  94      real    :: trj_lats(vals_max)
  95      real, allocatable :: chm_vals(:,:)
  96      real, allocatable :: trc_vals(:,:)
  97      real, allocatable :: hyd_vals(:,:)
  98      real, allocatable :: dyn_vals(:,:)
  99      real, allocatable :: msc_vals(:,:)
 100      real, allocatable :: dchm_vals(:,:)
 101      character(len=19) :: times(vals_max)
 102    end type buffer
 103
 104    type statevar
 105      character(len=80) :: Varname
 106      character(len=80) :: Description
 107      character(len=80) :: Units
 108      character(len=10) :: MemOrd
 109      character(len=10) :: Stagger
 110      integer           :: Ndim
 111      real, pointer     :: rfield_2d(:,:)
 112      real, pointer     :: rfield_3d(:,:,:)
 113    end type statevar
 114
 115    integer,    allocatable           :: traj_cnt(:)
 116    type(base), allocatable, target   :: traject(:,:)
 117    type(base), pointer               :: trjects(:)
 118    type(buffer), allocatable, target :: trj_buff(:,:)
 119    type(buffer), pointer             :: trj_pbf(:)
 120    type(statevar), allocatable, target :: St_Vars_dm(:,:)
 121    type(statevar), pointer             :: St_Vars(:)
 122
 123    real, allocatable  :: dchm_buff(:,:,:,:)
 124
 125    character(len=256) :: dyn_var_lst(dyn_max)
 126    character(len=32)  :: dyn_var_desc_att(dyn_max)
 127    character(len=32)  :: dyn_var_unit_att(dyn_max)
 128    character(len=4)   :: pkg_tag(pkg_max) = (/ 'chm ', 'hyd ', 'trc ', 'dyn ', 'msc ', 'dchm' /)
 129
 130    logical                      :: do_chemstep
 131    logical, allocatable, target :: pkg_has_vars_dm(:,:,:)
 132    logical, pointer             :: pkg_has_vars(:,:)
 133    logical, allocatable, target :: trj_msk_dm(:,:,:,:)
 134    logical, pointer             :: trj_msk(:,:)
 135    logical, allocatable, target :: St_Vars_msk_dm(:,:)
 136    logical, pointer             :: St_Vars_msk(:)
 137    logical, allocatable, target :: dchm_msk_dm(:,:)
 138    logical, pointer             :: dchm_msk(:)
 139    logical, allocatable         :: is_initialized(:)
 140    logical, allocatable         :: dm_has_traj(:)
 141
 142    CONTAINS
 143
 144    subroutine trajectory_init( grid, config_flags, &
 145                                ims,ime, jms,jme, kms,kme )
 146
 147    use module_domain
 148    use module_llxy,              only : proj_info, latlon_to_ij
 149    use module_configure,         only : grid_config_rec_type
 150    use module_state_description, only : no_trajectory, param_first_scalar, num_chem, num_moist, num_tracer
 151    use module_scalar_tables,     only : chem_dname_table, moist_dname_table, tracer_dname_table
 152    use module_model_constants,   only : g
 153    use module_domain_type,       only : fieldlist
 154 #ifdef DM_PARALLEL
 155    use module_dm,                only : wrf_dm_sum_integer, wrf_dm_max_real
 156 #endif
 157
 158 !-----------------------------------------------------------------------------
 159 !  dummy arguments
 160 !-----------------------------------------------------------------------------
 161    integer, intent(in)      :: ims,ime, jms,jme, kms,kme
 162    type(domain), intent(inout)            :: grid
 163    type(grid_config_rec_type), intent(in) :: config_flags
 164
 165 !-----------------------------------------------------------------------------
 166 !  local variables
 167 !-----------------------------------------------------------------------------
 168    integer :: astat
 169    integer :: dm
 170    integer, pointer :: num_msc      ! number misc species
 171    integer, pointer :: n_dchm       ! number dchm buffer species
 172    integer :: ierr, ios
 173    integer :: i, j, k, k1, m, m1, m2
 174    integer :: i_end, j_end, u_lim
 175    integer :: n, pkg, trj
 176    integer :: n_traj, n_traj_1
 177    integer :: p_size
 178    integer :: unitno
 179    integer :: ids,ide, jds,jde, kds,kde
 180    integer :: ips,ipe, jps,jpe, kps,kpe
 181    integer :: dbg_lvl
 182    integer :: target
 183    integer :: n_def_var(pkg_max)
 184    real    :: x, y
 185    real    :: z_dm_bot, z_dm_top
 186    real    :: z(kms:kme-1)
 187    real    :: z_at_w(ims:ime,kms:kme,jms:jme)
 188    character(len=256), allocatable :: msc_tbl(:)
 189    character(len=32)  :: var_name
 190    character(len=128) :: filename
 191    character(len=256) :: err_mes
 192    character(len=32)  :: wrk_var_name
 193    character(len=32)  :: wrk_chr(var_max)
 194    character(len=32)  :: wrk_def_name(var_max)
 195    logical :: exists
 196    logical :: is_root_proc
 197    logical :: rstrt
 198    logical :: mask(var_max)
 199    logical :: trj_mask(traj_max)
 200
 201    type(grid_config_rec_type) :: config_temp
 202    type(proj_info)            :: proj
 203
 204    type(default) :: traj_def
 205    type(base)    :: traj_type(traj_max)
 206
 207    logical, external :: wrf_dm_on_monitor
 208    integer, external :: get_unused_unit
 209
 210    namelist / traj_spec /    traj_type
 211    namelist / traj_default / traj_def
 212
 213 #ifndef NETCDF
 214    call wrf_error_fatal( 'trajectory_init: requires netcdf' )
 215 #endif
 216
 217    offset = param_first_scalar - 1
 218
 219    if( .not. allocated( dm_has_traj ) ) then
 220      call nl_get_max_dom( 1,n_dom )
 221      allocate( dm_has_traj(n_dom),traj_cnt(n_dom),stat=astat )
 222      if( astat /= 0 ) then
 223        write(err_mes,'(''trajectory_init('',i2.2,''): failed to allocate dm_has_traj,traj_cnt; error = '',i6)') dm,astat
 224        call wrf_error_fatal( trim( err_mes  ) )
 225      endif
 226      traj_cnt(:) = 0
 227    endif
 228
 229 !-----------------------------------------------------------------------------
 230 !  check for trajectory option
 231 !-----------------------------------------------------------------------------
 232    if( grid%traj_opt == no_trajectory ) then
 233      write(err_mes,'(''trajectory_init('',i2.2,''): traj_opt = no_trajectory'')') grid%id
 234      call wrf_message( trim(err_mes) )
 235      dm_has_traj(:) = .false.
 236      return
 237    endif
 238
 239 !-----------------------------------------------------------------------------
 240 !  domain requests trajectories
 241 !-----------------------------------------------------------------------------
 242    dm = grid%id
 243    if( .not. config_flags%dm_has_traj ) then
 244      write(err_mes,'(''trajectory_init('',i2.2,''): no trajectory calculation for domain '',i2.2)') dm,dm
 245      call wrf_message( trim( err_mes ) )
 246      dm_has_traj(dm) = .false.
 247      return
 248    else
 249      dm_has_traj(dm) = .true.
 250    endif
 251
 252 !-----------------------------------------------------------------------------
 253 !  set domain count, restarting?
 254 !-----------------------------------------------------------------------------
 255    call nl_get_restart( dm,rstrt )
 256    if( .not. allocated( traject ) ) then
 257 !-----------------------------------------------------------------------------
 258 !  allocate module variables
 259 !-----------------------------------------------------------------------------
 260      allocate( traject(traj_max,n_dom), &
 261                pkg_has_vars_dm(traj_max,pkg_max,n_dom), &
 262                num_msc_dm(n_dom), is_initialized(n_dom),stat=astat )
 263      if( astat /= 0 ) then
 264        write(err_mes,'(''trajectory_init('',i2.2,''): failed to allocate traject...num_msc_dm; error = '',i6)') dm,astat
 265        call wrf_error_fatal( trim( err_mes  ) )
 266      endif
 267      is_initialized(:) = .false.
 268    endif
 269
 270    if( is_initialized(dm) ) then
 271      return
 272    endif
 273
 274    trjects => traject(:,dm)
 275
 276    call get_ijk_from_grid( grid ,                   &
 277                            ids, ide, jds, jde, kds, kde, &
 278                            n,n,  n,n, n,n,          &
 279                            ips, ipe, jps, jpe, kps, kpe    )
 280
 281    n_vals      => n_vals_dm(dm)
 282    n_vals      = 0
 283    is_root_proc = wrf_dm_on_monitor()
 284    pkg_has_vars => pkg_has_vars_dm(:,:,dm)
 285
 286    dyn_var_lst(1:dyn_max)      = (/ 'p       ', 'T       ', 'z       ', 'u       ', &
 287                                     'v       ', 'w       ', 'rainprod', 'evapprod' /)
 288 #if( WRF_CHEM == 1 )
 289    if( grid%wetscav_onoff < 1 ) then
 290      dyn_var_lst(dyn_max-1:dyn_max) = (/ 'is_blank', 'is_blank' /)
 291    endif
 292 #endif
 293    dyn_var_unit_att(1:dyn_max) = (/ 'hPa', 'K  ', 'm  ', 'm/s', 'm/s', 'm/s', 's-1', 's-1' /)
 294    dyn_var_desc_att(1:dyn_max) = (/ 'pressure            ', 'temperature         ', 'height              ', &
 295                                     'x wind component    ', 'y wind component    ', 'z wind component    ', &
 296                                     'rain production rate', 'rain evap rate      ' /)
 297 !-----------------------------------------------------------------------------
 298 !  master proc
 299 !-----------------------------------------------------------------------------
 300 !master_proc: &
 301 !   if( is_root_proc ) then
 302       write(filename,'(''wrfinput_traj_d'',i2.2)',iostat=ios) dm
 303       if( ios /= 0 ) then
 304         write(err_mes,'(''trajectory_init('',i2.2,''): failed to set filename: error = '',i6)') dm,ios
 305         call wrf_error_fatal( trim( err_mes  ) )
 306       endif
 307       inquire( file=trim(filename),exist=exists )
 308 input_file: &
 309       if( exists ) then
 310         unitno = get_unused_unit()
 311         if( unitno <= 0 ) then
 312           call wrf_error_fatal( 'trajectory_init: failed to get unit number' )
 313         endif
 314 !-----------------------------------------------------------------------------
 315 !  open file
 316 !-----------------------------------------------------------------------------
 317         open( unit = unitno,file=trim(filename),iostat=ios )
 318         if( ios /= 0 ) then
 319           write(err_mes,'(''trajectory_init('',i2.2,''): failed to open '',a,''; error = '',i6)') dm,trim(filename),ios
 320           call wrf_error_fatal( trim( err_mes  ) )
 321         endif
 322 !-----------------------------------------------------------------------------
 323 !  initialize trajectories
 324 !-----------------------------------------------------------------------------
 325         traj_def%start_time = ' '
 326         traj_def%stop_time  = ' '
 327         traj_def%chm_name(:) = ' '
 328         traj_def%hyd_name(:) = ' '
 329         traj_def%trc_name(:) = ' '
 330         traj_def%dyn_name(:) = ' '
 331         traj_def%msc_name(:) = ' '
 332         traj_def%dchm_name(:) = ' '
 333
 334         do trj = 1,traj_max
 335           traj_type(trj)%start_time = ' '
 336           traj_type(trj)%stop_time  = ' '
 337           traj_type(trj)%chm_spc(:) = ' '
 338           traj_type(trj)%dchm_spc(:) = ' '
 339           traj_type(trj)%hyd_spc(:) = ' '
 340           traj_type(trj)%trc_spc(:) = ' '
 341           traj_type(trj)%dyn_var(:) = ' '
 342           traj_type(trj)%msc_var(:) = ' '
 343           traj_type(trj)%chm_ndx(:) = 0
 344           traj_type(trj)%hyd_ndx(:) = 0
 345           traj_type(trj)%trc_ndx(:) = 0
 346           traj_type(trj)%dyn_ndx(:) = 0
 347           traj_type(trj)%msc_ndx(:) = 0
 348           traj_type(trj)%dchm_ndx(:) = 0
 349           traj_type(trj)%n_chm_var = 0 ; traj_type(trj)%n_ct_var = 0
 350           traj_type(trj)%n_hyd_var = 0 ; traj_type(trj)%n_trc_var = 0
 351           traj_type(trj)%n_dyn_var = 0 ; traj_type(trj)%n_msc_var = 0
 352           traj_type(trj)%n_dchm_var = 0
 353           traj_type(trj)%is_stationary = .false.
 354           traj_type(trj)%chm_is_gas(:) = .true.
 355           traj_type(trj)%z_is_agl      = .true.
 356           traj_type(trj)%lon           = missing_val
 357           traj_type(trj)%lat           = missing_val
 358           traj_type(trj)%lev           = missing_val
 359         end do
 360 !-----------------------------------------------------------------------------
 361 !  read file
 362 !-----------------------------------------------------------------------------
 363         read(unit=unitno,nml=traj_default,iostat=ios)
 364         if( ios /= 0 ) then
 365           close( unit=unitno )
 366           write(err_mes,'(''trajectory_init('',i2.2,''): failed to read '',a,''; error = '',i6)') dm,trim(filename),ios
 367           call wrf_error_fatal( trim( err_mes  ) )
 368         endif
 369         read(unit=unitno,nml=traj_spec,iostat=ios)
 370         if( ios /= 0 ) then
 371           close( unit=unitno )
 372           write(err_mes,'(''trajectory_init('',i2.2,''): failed to read '',a,''; error = '',i6)') dm,trim(filename),ios
 373           call wrf_error_fatal( trim( err_mes  ) )
 374         endif
 375         close( unit=unitno )
 376       else input_file
 377         write(err_mes,'(''trajectory_init('',i2.2,''): no '',a,'' file'')') dm,trim(filename)
 378         call wrf_message( trim( err_mes ) )
 379         traj_cnt(dm) = 0
 380         return
 381       endif input_file
 382
 383 !-----------------------------------------------------------------------------
 384 !  process the namelist input
 385 !-----------------------------------------------------------------------------
 386       do pkg = 1,pkg_max
 387         select case( trim(pkg_tag(pkg)) )
 388           case( 'chm' )
 389             wrk_def_name(:) = traj_def%chm_name(:)
 390           case( 'dchm' )
 391             wrk_def_name(:) = traj_def%dchm_name(:)
 392           case( 'hyd' )
 393             wrk_def_name(:) = traj_def%hyd_name(:)
 394           case( 'trc' )
 395             wrk_def_name(:) = traj_def%trc_name(:)
 396           case( 'dyn' )
 397             wrk_def_name(:) = traj_def%dyn_name(:)
 398           case( 'msc' )
 399             wrk_def_name(:) = traj_def%msc_name(:)
 400         end select
 401         do m = 1,var_max
 402           if( wrk_def_name(m) == ' ' ) then
 403             exit
 404           endif
 405         end do
 406         n_def_var(pkg) = m - 1
 407       end do
 408
 409       if( traj_def%start_time /= ' ' ) then
 410         write(err_mes,'(''trajectory_init('',i2.2,''): default start time = '',a)') dm,traj_def%start_time
 411         call wrf_message( trim(err_mes) )
 412       endif
 413       if( traj_def%stop_time /= ' ' ) then
 414         write(err_mes,'(''trajectory_init('',i2.2,''): default stop  time = '',a)') dm,traj_def%stop_time
 415         call wrf_message( trim(err_mes) )
 416       endif
 417
 418       do pkg = 1,pkg_max
 419         if( n_def_var(pkg) > 0 ) then
 420           write(*,*) ' '
 421           write(*,'(''trajectory_init('',i2.2,''): default '',a,'' variables'')') dm,pkg_tag(pkg)
 422           select case( trim(pkg_tag(pkg)) )
 423 #if( WRF_CHEM == 1 )
 424             case( 'chm' )
 425               wrk_def_name(:) = traj_def%chm_name(:)
 426             case( 'dchm' )
 427               wrk_def_name(:) = traj_def%dchm_name(:)
 428 #endif
 429             case( 'hyd' )
 430               wrk_def_name(:) = traj_def%hyd_name(:)
 431             case( 'trc' )
 432               wrk_def_name(:) = traj_def%trc_name(:)
 433             case( 'dyn' )
 434               wrk_def_name(:) = traj_def%dyn_name(:)
 435             case( 'msc' )
 436               wrk_def_name(:) = traj_def%msc_name(:)
 437           end select
 438           write(*,*) wrk_def_name(:n_def_var(pkg))
 439         endif
 440       end do
 441
 442       do n_traj = 1,traj_max
 443         if( traj_type(n_traj)%lon == missing_val .or. &
 444             traj_type(n_traj)%lat == missing_val .or. &
 445             traj_type(n_traj)%lev == missing_val ) then
 446           exit
 447         endif
 448       end do
 449       n_traj = n_traj - 1
 450
 451 has_trajectories: &
 452       if( n_traj > 0 ) then
 453 !-----------------------------------------------------------------------------
 454 !  set individual trajectories to default if specified
 455 !-----------------------------------------------------------------------------
 456         if( traj_def%start_time /= ' ' ) then
 457           do trj = 1,n_traj
 458             if( traj_type(trj)%start_time == ' ' ) then
 459               traj_type(trj)%start_time = traj_def%start_time
 460             endif
 461           end do
 462         endif
 463         if( traj_def%stop_time /= ' ' ) then
 464           do trj = 1,n_traj
 465             if( traj_type(trj)%stop_time == ' ' ) then
 466               traj_type(trj)%stop_time = traj_def%stop_time
 467             endif
 468           end do
 469         endif
 470
 471         do pkg = 1,pkg_max
 472           select case( trim(pkg_tag(pkg)) )
 473 #if( WRF_CHEM == 1 )
 474             case( 'chm' )
 475               wrk_def_name(:) = traj_def%chm_name(:)
 476             case( 'dchm' )
 477               wrk_def_name(:) = traj_def%dchm_name(:)
 478 #endif
 479             case( 'hyd' )
 480               wrk_def_name(:) = traj_def%hyd_name(:)
 481             case( 'trc' )
 482               wrk_def_name(:) = traj_def%trc_name(:)
 483             case( 'dyn' )
 484               wrk_def_name(:) = traj_def%dyn_name(:)
 485             case( 'msc' )
 486               wrk_def_name(:) = traj_def%msc_name(:)
 487           end select
 488           do trj = 1,n_traj
 489             select case( trim(pkg_tag(pkg)) )
 490 #if( WRF_CHEM == 1 )
 491               case( 'chm' )
 492                 wrk_var_name = traj_type(trj)%chm_spc(1)
 493               case( 'dchm' )
 494                 wrk_var_name = traj_type(trj)%dchm_spc(1)
 495 #endif
 496               case( 'hyd' )
 497                 wrk_var_name = traj_type(trj)%hyd_spc(1)
 498               case( 'trc' )
 499                 wrk_var_name = traj_type(trj)%trc_spc(1)
 500               case( 'dyn' )
 501                 wrk_var_name = traj_type(trj)%dyn_var(1)
 502               case( 'msc' )
 503                 wrk_var_name = traj_type(trj)%msc_var(1)
 504             end select
 505             if( wrk_var_name == ' ' ) then
 506               m1 = n_def_var(pkg)
 507               select case( trim(pkg_tag(pkg)) )
 508 #if( WRF_CHEM == 1 )
 509                 case( 'chm' )
 510                   traj_type(trj)%chm_spc(:m1) = traj_def%chm_name(:m1)
 511                 case( 'dchm' )
 512                   traj_type(trj)%dchm_spc(:m1) = traj_def%dchm_name(:m1)
 513 #endif
 514                 case( 'hyd' )
 515                   traj_type(trj)%hyd_spc(:m1) = traj_def%hyd_name(:m1)
 516                 case( 'trc' )
 517                   traj_type(trj)%trc_spc(:m1) = traj_def%trc_name(:m1)
 518                 case( 'dyn' )
 519                   traj_type(trj)%dyn_var(:m1) = traj_def%dyn_name(:m1)
 520                 case( 'msc' )
 521                   traj_type(trj)%msc_var(:m1) = traj_def%msc_name(:m1)
 522               end select
 523             endif
 524           end do
 525         end do
 526 !-----------------------------------------------------------------------------
 527 !  scan registry real, 2d, 3d variables
 528 !-----------------------------------------------------------------------------
 529         call reg_scan( grid )
 530         num_msc => num_msc_dm(dm)
 531
 532         call get_wrf_debug_level( dbg_lvl )
 533         if( dbg_lvl > 200 ) then
 534           write(*,*) ' '
 535           write(*,'(''trajectory_init('',i2.2,''): Registry 2d,3d variables'')') dm
 536           do n = 1,num_msc
 537             write(*,*) trim(St_Vars(n)%varname)
 538           end do
 539           n = count( St_Vars(:num_msc)%Stagger == 'X' )
 540           if( n > 0 ) then
 541             write(*,*) ' '
 542             write(*,*) 'Registry 2d,3d variables with staggered X'
 543             do n = 1,num_msc
 544               if( St_Vars(n)%Stagger == 'X' ) then
 545                 write(*,*) trim(St_Vars(n)%varname)
 546               endif
 547             end do
 548           endif
 549           n = count( St_Vars(:num_msc)%Stagger == 'Y' )
 550           if( n > 0 ) then
 551             write(*,*) ' '
 552             write(*,*) 'trajectory_init: Registry 2d,3d variables with staggered Y'
 553             do n = 1,num_msc
 554               if( St_Vars(n)%Stagger == 'Y' ) then
 555                 write(*,*) trim(St_Vars(n)%varname)
 556               endif
 557             end do
 558           endif
 559           n = count( St_Vars(:num_msc)%Stagger == 'Z' )
 560           if( n > 0 ) then
 561             write(*,*) ' '
 562             write(*,*) 'trajectory_init: Registry 2d,3d variables with staggered Z'
 563             do n = 1,num_msc
 564               if( St_Vars(n)%Stagger == 'Z' ) then
 565                 write(*,*) trim(St_Vars(n)%varname)
 566               endif
 567             end do
 568           endif
 569         endif
 570 !-----------------------------------------------------------------------------
 571 !  get variable counts
 572 !-----------------------------------------------------------------------------
 573         do trj = 1,n_traj
 574           if( num_chem > 1 ) then
 575 #if( WRF_CHEM == 1 )
 576             call get_var_cnt( traj_type(trj)%n_chm_var, traj_type(trj)%chm_spc )
 577             call get_var_cnt( traj_type(trj)%n_dchm_var, traj_type(trj)%dchm_spc )
 578 #endif
 579           else
 580             traj_type(trj)%n_chm_var  = 0
 581             traj_type(trj)%n_dchm_var = 0
 582           endif
 583           if( num_moist > 1 ) then
 584             call get_var_cnt( traj_type(trj)%n_hyd_var, traj_type(trj)%hyd_spc )
 585           else
 586             traj_type(trj)%n_hyd_var = 0
 587           endif
 588           if( num_tracer > 1 ) then
 589             call get_var_cnt( traj_type(trj)%n_trc_var, traj_type(trj)%trc_spc )
 590           else
 591             traj_type(trj)%n_trc_var = 0
 592           endif
 593           if( num_msc > 1 ) then
 594             call get_var_cnt( traj_type(trj)%n_msc_var, traj_type(trj)%msc_var )
 595           else
 596             traj_type(trj)%n_msc_var = 0
 597           endif
 598           call get_var_cnt( traj_type(trj)%n_dyn_var, traj_type(trj)%dyn_var )
 599         end do
 600
 601         if( any( traj_type(:n_traj)%n_msc_var > 0 ) ) then
 602           allocate( msc_tbl(num_msc),stat=astat )
 603           if( astat /= 0 ) then
 604             call wrf_error_fatal( 'trajectory_init: failed to find allocate msc_tbl' )
 605           endif
 606           do m = 1,num_msc
 607             msc_tbl(m) = trim( St_Vars(m)%Varname )
 608           end do
 609         endif
 610 !-----------------------------------------------------------------------------
 611 !  check for trajectory variables in simulation
 612 !-----------------------------------------------------------------------------
 613         do trj = 1,n_traj
 614 #if( WRF_CHEM == 1 )
 615           if( num_chem > 1 ) then
 616             if( traj_type(trj)%n_chm_var > 0 ) then
 617               mask(:) = .false.
 618               call scan_vars( traj_type(trj)%n_chm_var, traj_type(trj)%chm_spc, traj_type(trj)%chm_ndx, &
 619                               num_chem, chem_dname_table(dm,:), &
 620                               traj_type(trj)%chm_is_gas, 'chm' )
 621             endif
 622             if( traj_type(trj)%n_dchm_var > 0 ) then
 623               mask(:) = .false.
 624               call scan_vars( traj_type(trj)%n_dchm_var, traj_type(trj)%dchm_spc, traj_type(trj)%dchm_ndx, &
 625                               num_chem, chem_dname_table(dm,:), &
 626                               traj_type(trj)%chm_is_gas, 'chm' )
 627             endif
 628           endif
 629 #endif
 630           if( traj_type(trj)%n_hyd_var > 0 .and. num_moist > 1 ) then
 631             mask(:) = .false.
 632             call scan_vars( traj_type(trj)%n_hyd_var, traj_type(trj)%hyd_spc, traj_type(trj)%hyd_ndx, &
 633                             num_moist, moist_dname_table(dm,:), &
 634                             traj_type(trj)%chm_is_gas, 'hyd' )
 635           endif
 636           if( traj_type(trj)%n_trc_var > 0 .and. num_tracer > 1 ) then
 637             mask(:) = .false.
 638             call scan_vars( traj_type(trj)%n_trc_var, traj_type(trj)%trc_spc, traj_type(trj)%trc_ndx, &
 639                             num_tracer, tracer_dname_table(dm,:), &
 640                             traj_type(trj)%chm_is_gas, 'trc' )
 641           endif
 642           if( traj_type(trj)%n_dyn_var > 0 ) then
 643             mask(:) = .false.
 644             call scan_vars( traj_type(trj)%n_dyn_var, traj_type(trj)%dyn_var, traj_type(trj)%dyn_ndx, &
 645                             dyn_max, dyn_var_lst(:), &
 646                             traj_type(trj)%chm_is_gas, 'dyn' )
 647           endif
 648           if( traj_type(trj)%n_msc_var > 0 ) then
 649             mask(:) = .false.
 650             call scan_vars( traj_type(trj)%n_msc_var, traj_type(trj)%msc_var, traj_type(trj)%msc_ndx, &
 651                             num_msc, msc_tbl(:), &
 652                             traj_type(trj)%chm_is_gas, 'msc' )
 653           endif
 654         end do
 655
 656         do trj = 1,n_traj
 657           if( traj_type(trj)%n_msc_var > 0 ) then
 658             do m = 1,traj_type(trj)%n_msc_var
 659               St_Vars_msk(traj_type(trj)%msc_ndx(m)-offset) = .true.
 660             end do
 661           endif
 662         end do
 663         m = count( St_Vars_msk(:num_msc) )
 664
 665         if( allocated( msc_tbl ) ) then
 666           deallocate( msc_tbl )
 667         endif
 668
 669 !-----------------------------------------------------------------------------
 670 !  remove trajectories with no variables
 671 !-----------------------------------------------------------------------------
 672         n_traj_1 = count( (traj_type(:n_traj)%n_chm_var + traj_type(:n_traj)%n_hyd_var &
 673                            + traj_type(:n_traj)%n_trc_var + traj_type(:n_traj)%n_dyn_var &
 674                            + traj_type(:n_traj)%n_msc_var + traj_type(:n_traj)%n_dchm_var) > 0 )
 675         if( n_traj_1 > 0 ) then
 676           if( n_traj_1 /= n_traj ) then
 677             trj_mask(1:n_traj) = traj_type(1:n_traj)%in_dom
 678             m = 1
 679             do trj = 1,n_traj
 680               if( trj_mask(trj) ) then
 681                 if( trj /= m ) then
 682                   traj_type(m) = traj_type(trj)
 683                   m = m + 1
 684                 endif
 685               endif
 686             end do
 687             n_traj = n_traj_1
 688           endif
 689         else
 690           dm_has_traj(dm) = .false.
 691           return
 692         endif
 693 !-----------------------------------------------------------------------------
 694 !  allocate buffer type
 695 !-----------------------------------------------------------------------------
 696         if( is_root_proc ) then
 697           if( dm == 1 ) then
 698             allocate( trj_buff(traj_max,n_dom),stat=astat )
 699             if( astat /= 0 ) then
 700               write(err_mes,'(''trajectory_init: failed to allocate traj_buff; error = '',i6)') astat
 701               call wrf_error_fatal( trim( err_mes  ) )
 702             endif
 703           endif
 704           trj_pbf => trj_buff(:,dm)
 705         endif
 706       endif has_trajectories
 707 !   endif master_proc
 708 #if( WRF_CHEM == 1 )
 709 !-----------------------------------------------------------------------------
 710 !   for dchm package make sure dchm_ndx is ordered list
 711 !-----------------------------------------------------------------------------
 712     do trj = 1,n_traj
 713       n = traj_type(trj)%n_dchm_var
 714       if( n > 0 .and.  any( traj_type(trj)%dchm_ndx(1:n-1) > traj_type(trj)%dchm_ndx(2:n) ) ) then
 715         trj_mask(:num_chem) = .false.
 716         do m = 1,n
 717           trj_mask(traj_type(trj)%dchm_ndx(m)) = .true.
 718         end do
 719         traj_type(trj)%dchm_ndx(:n) = pack( (/ (m,m=1,num_chem) /),trj_mask(:num_chem) )
 720         do m = 1,n
 721           m1 = traj_type(trj)%dchm_ndx(m)
 722           traj_type(trj)%dchm_spc(m) = trim(chem_dname_table(dm,m1))
 723         end do
 724       endif
 725     end do
 726 #endif
 727 !-----------------------------------------------------------------------------
 728 !  if initial run, overwrite existing grid trajectory variables
 729 !-----------------------------------------------------------------------------
 730 #ifdef DM_PARALLEL
 731 !  call wrf_dm_bcast_integer( n_traj,1 )
 732 !  if( n_traj > 0 ) then
 733 !    p_size = (5+var_max)*RWORDSIZE + (2+var_max)*LWORDSIZE &
 734 !           + (5+4*var_max)*IWORDSIZE + 4*32*var_max + 38
 735 !    do trj = 1,n_traj
 736 !      call wrf_dm_bcast_bytes( traj_type(trj),p_size )
 737 !    end do
 738 !  endif
 739 #endif
 740 !  call wrf_abort( 'Debugging' )
 741 is_cold_start: &
 742    if( .not. rstrt ) then
 743 !-----------------------------------------------------------------------------
 744 !  calc and check trajectory start point
 745 !-----------------------------------------------------------------------------
 746 has_trajectories_a: &
 747      if( n_traj > 0 ) then
 748        config_temp = config_flags
 749        call trajmapproj( grid, config_temp, proj )
 750        i_end = min(ipe,ide-1)
 751        j_end = min(jpe,jde-1)
 752        do j = jps,j_end
 753          do k = kps,kpe
 754            z_at_w(ips:i_end,k,j) = (grid%ph_2(ips:i_end,k,j) + grid%phb(ips:i_end,k,j))/g
 755          end do
 756        end do
 757 !-----------------------------------------------------------------------------
 758 !  first check x,y
 759 !-----------------------------------------------------------------------------
 760 traj_loop: &
 761        do trj = 1,n_traj
 762          if( traj_type(trj)%lat /= missing_val .and. &
 763              traj_type(trj)%lon /= missing_val ) then
 764            call latlon_to_ij( proj, traj_type(trj)%lat, traj_type(trj)%lon, x, y )
 765            traj_type(trj)%in_dom = &
 766                   (x >= real(ids) .and. x <= real(ide-1) .and. &
 767                    y >= real(jds) .and. y <= real(jde-1))
 768 #ifdef SW_DEBUG
 769           write(err_mes,'(''trajectory_init('',i2.2,''): x,y = '',1p,2g15.7)') dm,x,y
 770           call wrf_debug( 0,trim(err_mes) )
 771 #endif
 772          else
 773            traj_type(trj)%in_dom = .false.
 774          endif
 775 #ifdef SW_DEBUG
 776          write(err_mes,'(''trajectory_init('',i2.2,''): traj '',i4,'' in dom  = '',l)') dm,trj,traj_type(trj)%in_dom
 777          call wrf_debug( 0,trim(err_mes) )
 778 #endif
 779 !-----------------------------------------------------------------------------
 780 !  then check z
 781 !-----------------------------------------------------------------------------
 782 is_in_domain: &
 783          if( traj_type(trj)%in_dom ) then
 784            i = nint( x )
 785            j = nint( y )
 786            traj_type(trj)%in_patch = &
 787                   (i >= ips .and. i <= min(ipe,ide-1) .and. &
 788                    j >= jps .and. j <= min(jpe,jde-1))
 789 is_in_patch: &
 790            if( traj_type(trj)%in_patch ) then
 791              k1 = kde - 1
 792              if( traj_type(trj)%z_is_agl ) then
 793                traj_type(trj)%lev = traj_type(trj)%lev + z_at_w(i,kds,j)
 794              endif
 795              z_dm_bot = z_at_w(i,kds,j)
 796              z_dm_top = z_at_w(i,k1,j)
 797              write(err_mes,'(''trajectory_init('',i2.2,''): traj '',i3.3,'' i,j,z_bot,z_top,lev = '',2i4,1p3g15.7)') &
 798                 dm,trj,i,j,z_dm_bot,z_dm_top,traj_type(trj)%lev
 799              call wrf_debug( 0,trim(err_mes) )
 800              if( traj_type(trj)%lev >= z_dm_bot .and. &
 801                  traj_type(trj)%lev <= z_dm_top ) then
 802                traj_type(trj)%in_dom = .true.
 803                traj_type(trj)%x      = x
 804                traj_type(trj)%y      = y
 805              else
 806                traj_type(trj)%in_dom = .false.
 807              endif
 808            else is_in_patch
 809              traj_type(trj)%x = missing_val
 810              traj_type(trj)%y = missing_val
 811            endif is_in_patch
 812          endif is_in_domain
 813        end do traj_loop
 814        n_traj_1 = count( traj_type(:n_traj)%in_dom )
 815      else has_trajectories_a
 816        n_traj_1 = 0
 817      endif has_trajectories_a
 818
 819 !-----------------------------------------------------------------------------
 820 !  remove out of domain trajectories
 821 !-----------------------------------------------------------------------------
 822      write(err_mes,'(''trajectory_init('',i2.2,''): traj cnt = '',2i6)') dm,n_traj_1,n_traj
 823      call wrf_debug( 0,trim(err_mes) )
 824      if( n_traj_1 /= n_traj ) then
 825        trj_mask(1:n_traj) = traj_type(1:n_traj)%in_dom
 826        m = 1
 827        do trj = 1,n_traj
 828          if( trj_mask(trj) ) then
 829            if( trj /= m ) then
 830              traj_type(m) = traj_type(trj)
 831              m = m + 1
 832            endif
 833          endif
 834        end do
 835        n_traj = n_traj_1
 836      endif
 837
 838 has_trajectories_b: &
 839      if( n_traj_1 > 0 ) then
 840        grid%traj_i(:) = missing_val
 841        grid%traj_j(:) = missing_val
 842        grid%traj_k(:) = missing_val
 843        grid%traj_long(:) = missing_val
 844        grid%traj_lat(:)  = missing_val
 845 !-----------------------------------------------------------------------------
 846 !  set initial trajectory spatial coordinates
 847 !-----------------------------------------------------------------------------
 848        k1 = kde - 1
 849        do trj = 1,n_traj
 850          if( traj_type(trj)%in_patch ) then
 851            grid%traj_i(trj) = traj_type(trj)%x
 852            grid%traj_j(trj) = traj_type(trj)%y
 853            grid%traj_long(trj) = traj_type(trj)%lon
 854            grid%traj_lat(trj)  = traj_type(trj)%lat
 855            i = nint( traj_type(trj)%x )
 856            j = nint( traj_type(trj)%y )
 857 !          z(kds:k1) = .5*(z_at_w(i,kds:k1,j) + z_at_w(i,kds+1:k1+1,j))
 858            z(kds:k1) = z_at_w(i,kds:k1,j)
 859            do k = kds+1,k1
 860              if( traj_type(trj)%lev <= z(k) ) then
 861                grid%traj_k(trj) = real(k - 1) &
 862                               + (traj_type(trj)%lev - z(k-1))/(z(k) - z(k-1))
 863                exit
 864              endif
 865            end do
 866            write(err_mes,'(''trajectory_init('',i2.2,''): trj,k,z(k-1:k),traj_k = '',2i3,1p3g15.7)') &
 867               dm,trj,k,z(k-1:k),grid%traj_k(trj)
 868            call wrf_debug( 0,trim(err_mes) )
 869          endif
 870        end do
 871      else has_trajectories_b
 872        dm_has_traj(dm) = .false.
 873        return
 874      endif has_trajectories_b
 875    else is_cold_start
 876      if( n_traj > 0 ) then
 877        call set_in_dom
 878      else
 879        traj_cnt(dm) = n_traj
 880        return
 881      endif
 882    endif is_cold_start
 883
 884 !-----------------------------------------------------------------------------
 885 !  transfer from local variable to module variable
 886 !-----------------------------------------------------------------------------
 887    traj_cnt(dm) = n_traj
 888    do trj = 1,n_traj
 889      trjects(trj) = traj_type(trj)
 890    end do
 891 !-----------------------------------------------------------------------------
 892 !  create netcdf trajectory file
 893 !-----------------------------------------------------------------------------
 894    if( .not. rstrt ) then
 895      do trj = 1,n_traj
 896 #ifdef DM_PARALLEL
 897        grid%traj_i(trj) = wrf_dm_max_real( grid%traj_i(trj) )
 898        grid%traj_j(trj) = wrf_dm_max_real( grid%traj_j(trj) )
 899        grid%traj_k(trj) = wrf_dm_max_real( grid%traj_k(trj) )
 900        grid%traj_long(trj) = wrf_dm_max_real( grid%traj_long(trj) )
 901        grid%traj_lat(trj)  = wrf_dm_max_real( grid%traj_lat(trj) )
 902 #else
 903        grid%traj_i(trj) = grid%traj_i(trj)
 904        grid%traj_j(trj) = grid%traj_j(trj)
 905        grid%traj_k(trj) = grid%traj_k(trj)
 906        grid%traj_long(trj) = grid%traj_long(trj)
 907        grid%traj_lat(trj)  = grid%traj_lat(trj)
 908 #endif
 909      end do
 910    endif
 911
 912    do pkg = 1,pkg_max
 913      select case( trim(pkg_tag(pkg)) )
 914 #if( WRF_CHEM == 1 )
 915        case( 'chm' )
 916          pkg_has_vars(:n_traj,pkg) = trjects(:n_traj)%n_chm_var > 0
 917        case( 'dchm' )
 918          pkg_has_vars(:n_traj,pkg) = trjects(:n_traj)%n_dchm_var > 0
 919 #endif
 920        case( 'hyd' )
 921          pkg_has_vars(:n_traj,pkg) = trjects(:n_traj)%n_hyd_var > 0
 922        case( 'trc' )
 923          pkg_has_vars(:n_traj,pkg) = trjects(:n_traj)%n_trc_var > 0
 924        case( 'dyn' )
 925          pkg_has_vars(:n_traj,pkg) = trjects(:n_traj)%n_dyn_var > 0
 926        case( 'msc' )
 927          pkg_has_vars(:n_traj,pkg) = trjects(:n_traj)%n_msc_var > 0
 928      end select
 929    end do
 930
 931 #if( WRF_CHEM == 1 )
 932 !-----------------------------------------------------------------------------
 933 !  setup for  dchm buffer
 934 !-----------------------------------------------------------------------------
 935    n_dchm_dm(dm) = 0
 936    if( any( pkg_has_vars(:n_traj,dchm_pkg) ) ) then
 937      if( .not. allocated( dchm_msk_dm ) ) then
 938        allocate( dchm_msk_dm(num_chem,n_dom),stat=astat )
 939        if( astat /= 0 ) then
 940          write(err_mes,'(''trajectory_init('',i2.2,''): failed to allocate dchm_msk_dm; error = '',i6)') dm,astat
 941          call wrf_error_fatal( trim( err_mes  ) )
 942        endif
 943      endif
 944      if( .not. allocated( dchm_buf_ndx_dm ) ) then
 945        allocate( dchm_buf_ndx_dm(num_chem,n_dom),stat=astat )
 946        if( astat /= 0 ) then
 947          write(err_mes,'(''trajectory_init('',i2.2,''): failed to allocate dchm_buf_ndx_dm; error = '',i6)') dm,astat
 948          call wrf_error_fatal( trim( err_mes  ) )
 949        endif
 950      endif
 951      n_dchm   => n_dchm_dm(dm)
 952      dchm_msk => dchm_msk_dm(:,dm)
 953      dchm_msk(:) = .false.
 954      dchm_buf_ndx => dchm_buf_ndx_dm(:,dm)
 955      dchm_buf_ndx(:) = 0
 956      do trj = 1,n_traj
 957        do m = 1,trjects(trj)%n_dchm_var
 958          dchm_msk(trjects(trj)%dchm_ndx(m)) = .true.
 959        end do
 960      end do
 961      n_dchm = count( dchm_msk )
 962      dchm_buf_ndx(:n_dchm) = pack( (/ (m,m=1,num_chem) /),dchm_msk(:num_chem) )
 963      do trj = 1,n_traj
 964        if( trjects(trj)%n_dchm_var > 0 ) then
 965          do m1 = 1,trjects(trj)%n_dchm_var
 966            target = trjects(trj)%dchm_ndx(m1)
 967            do m2 = 1,n_dchm
 968              if( target == dchm_buf_ndx(m2) ) then
 969                trjects(trj)%dchm_ndx(m1) = m2 + offset
 970                exit
 971              endif
 972            end do
 973          end do
 974        endif
 975      end do
 976    endif
 977 #endif
 978
 979 master_proc_a: &
 980    if( is_root_proc ) then
 981      if( dm == 1 ) then
 982        n = max(num_chem,num_moist,num_tracer,num_msc+offset,dyn_max+offset)
 983        if( n  > offset .and. .not. allocated(trj_msk_dm) ) then
 984          allocate( trj_msk_dm(traj_max,n,pkg_max,n_dom),stat=astat )
 985          if( astat /= 0 ) then
 986            write(err_mes,'(''trajectory_init: failed to allocate trj_msk_dm; error = '',i6)') astat
 987            call wrf_error_fatal( trim( err_mes  ) )
 988          endif
 989          trj_msk_dm(:,:,:,:) = .false.
 990        endif
 991      endif
 992 is_initial: &
 993      if( .not. is_initialized(dm) ) then
 994 !-----------------------------------------------------------------------------
 995 !  allocate buffer arrays
 996 !-----------------------------------------------------------------------------
 997 trj_loop: &
 998        do trj = 1,n_traj
 999 pkg_loop:  do pkg = 1,pkg_max
1000              astat = 0
1001              select case( trim(pkg_tag(pkg)) )
1002 #if( WRF_CHEM == 1 )
1003                case( 'chm' )
1004                  trj_msk => trj_msk_dm(:,:,chm_pkg,dm)
1005                  m1 = trjects(trj)%n_chm_var
1006                  if( m1 > 0 ) then
1007                    allocate( trj_pbf(trj)%chm_vals(vals_max,m1),stat=astat)
1008                    do m = 1,m1
1009                      trj_msk(trj,trjects(trj)%chm_ndx(m)) = trjects(trj)%in_dom
1010                    end do
1011                  endif
1012                case( 'dchm' )
1013                  trj_msk => trj_msk_dm(:,:,dchm_pkg,dm)
1014                  m1 = trjects(trj)%n_dchm_var
1015                  if( m1 > 0 ) then
1016                    allocate( trj_pbf(trj)%dchm_vals(vals_max,m1),stat=astat)
1017                    do m = 1,m1
1018                      trj_msk(trj,trjects(trj)%dchm_ndx(m)) = trjects(trj)%in_dom
1019                    end do
1020                  endif
1021 #endif
1022                case( 'hyd' )
1023                  trj_msk => trj_msk_dm(:,:,hyd_pkg,dm)
1024                  m1 = trjects(trj)%n_hyd_var
1025                  if( m1 > 0 ) then
1026                    allocate( trj_pbf(trj)%hyd_vals(vals_max,m1),stat=astat)
1027                    do m = 1,m1
1028                      trj_msk(trj,trjects(trj)%hyd_ndx(m)) = trjects(trj)%in_dom
1029                    end do
1030                  endif
1031                case( 'trc' )
1032                  trj_msk => trj_msk_dm(:,:,trc_pkg,dm)
1033                  m1 = trjects(trj)%n_trc_var
1034                  if( m1 > 0 ) then
1035                    allocate( trj_pbf(trj)%trc_vals(vals_max,m1),stat=astat)
1036                    do m = 1,m1
1037                      trj_msk(trj,trjects(trj)%trc_ndx(m)) = trjects(trj)%in_dom
1038                    end do
1039                  endif
1040                case( 'dyn' )
1041                  trj_msk => trj_msk_dm(:,:,dyn_pkg,dm)
1042                  m1 = trjects(trj)%n_dyn_var
1043                  if( m1 > 0 ) then
1044                    allocate( trj_pbf(trj)%dyn_vals(vals_max,m1),stat=astat)
1045                    do m = 1,m1
1046                      trj_msk(trj,trjects(trj)%dyn_ndx(m)) = trjects(trj)%in_dom
1047                    end do
1048                  endif
1049                case( 'msc' )
1050                  trj_msk => trj_msk_dm(:,:,msc_pkg,dm)
1051                  m1 = trjects(trj)%n_msc_var
1052                  if( m1 > 0 ) then
1053                    allocate( trj_pbf(trj)%msc_vals(vals_max,m1),stat=astat)
1054                    do m = 1,m1
1055                      trj_msk(trj,trjects(trj)%msc_ndx(m)) = trjects(trj)%in_dom
1056                    end do
1057                  endif
1058              end select
1059              if( astat /= 0 ) then
1060                write(err_mes,'(''trajectory_init: failed to allocate buffer%'',a,''; error = '',i6)') &
1061                    pkg_tag(pkg),astat
1062                call wrf_error_fatal( trim( err_mes  ) )
1063              endif
1064            end do pkg_loop
1065        end do trj_loop
1066
1067        do pkg = 1,pkg_max
1068          select case( trim(pkg_tag(pkg)) )
1069 #if( WRF_CHEM == 1 )
1070            case( 'chm' )
1071              trj_msk => trj_msk_dm(:,:,chm_pkg,dm)
1072              u_lim = num_chem
1073            case( 'dchm' )
1074              trj_msk => trj_msk_dm(:,:,dchm_pkg,dm)
1075              u_lim = num_chem
1076 #endif
1077            case( 'hyd' )
1078              trj_msk => trj_msk_dm(:,:,hyd_pkg,dm)
1079              u_lim = num_moist
1080            case( 'trc' )
1081              trj_msk => trj_msk_dm(:,:,trc_pkg,dm)
1082              u_lim = num_tracer
1083            case( 'dyn' )
1084              trj_msk => trj_msk_dm(:,:,dyn_pkg,dm)
1085              u_lim = dyn_max + offset
1086            case( 'msc' )
1087              trj_msk => trj_msk_dm(:,:,msc_pkg,dm)
1088              u_lim = num_msc + offset
1089          end select
1090          do trj = 1,n_traj
1091            trj_msk(trj,1) = any( trj_msk(trj,param_first_scalar:u_lim) )
1092          end do
1093        end do
1094        is_initialized(dm) = .true.
1095        if( .not. rstrt ) then
1096          call trajectory_create_file( grid, n_traj )
1097        endif
1098      endif is_initial
1099    endif master_proc_a
1100 #ifdef DM_PARALLEL
1101    call wrf_dm_bcast_logical( is_initialized,n_dom )
1102 #endif
1103
1104    CONTAINS
1105
1106    subroutine reg_scan( grid )
1107
1108    use module_domain_type, only : domain, fieldlist
1109
1110 !-----------------------------------------------------------------------------
1111 !  dummy arguments
1112 !-----------------------------------------------------------------------------
1113    type(domain), intent(in) :: grid
1114
1115    integer, parameter :: nVerbotten = 8
1116 !-----------------------------------------------------------------------------
1117 !  local variables
1118 !-----------------------------------------------------------------------------
1119    integer         :: astat
1120    integer         :: cnt
1121    integer         :: dm, dm_ndx
1122    integer         :: m, n
1123    type(fieldlist), pointer :: p
1124    logical         :: valid
1125    type(statevar), allocatable :: St_Vars_wrk(:,:)
1126    logical, allocatable :: St_Vars_msk_wrk(:,:)
1127    character(len=80) :: tstring
1128    character(len=9) :: Verbotten(nVerbotten) = (/ 'zx       ', 'zy       ', &
1129                                                   'RUNDGTEN ', 'RVNDGTEN ', &
1130                                                   'U_NDG_OLD', 'V_NDG_OLD', &
1131                                                   'U_NDG_NEW', 'V_NDG_NEW'/)
1132
1133 !-----------------------------------------------------------------------------
1134 !  just get the count
1135 !-----------------------------------------------------------------------------
1136    dm = grid%id
1137    p => grid%head_statevars%next ; cnt = 0
1138    do while( associated(p) )
1139      valid = (p%Type == 'R' .or. p%Type == 'r') .and. &
1140              (p%Ndim == 3 .or. (p%Ndim == 4 .and. p%scalar_array))
1141      if( valid ) then
1142        valid = p%MemoryOrder(:3) == 'XZY' .and. (p%em2 == kde-1 .or. p%em2 == kde)
1143        if( valid ) then
1144          if( p%Ndim == 3 ) then
1145            do m = 1,nVerbotten
1146              if( trim(Verbotten(m)) == trim(p%Varname) ) then
1147                valid = .false.
1148                exit
1149              endif
1150            enddo
1151            if( valid ) then
1152              cnt = cnt + 1
1153            endif
1154          elseif( p%Ndim == 4 ) then
1155            do n = param_first_scalar,p%num_table(dm)
1156              valid = .true.
1157              do m = 1,nVerbotten
1158                tstring = trim(Verbotten(m))
1159                call upcase( tstring )
1160                if( trim(tstring) == trim(p%dname_table(dm,n)) ) then
1161                  valid = .false.
1162                  exit
1163                endif
1164              enddo
1165              if( valid ) then
1166                cnt = cnt + 1
1167              endif
1168            enddo
1169          endif
1170        endif
1171      endif
1172      p => p%next
1173    end do
1174
1175    write(*,'(''reg_scan('',i2.2,''): found '',i4,'' state variables '')') dm,cnt
1176
1177 !-----------------------------------------------------------------------------
1178 !  now allocate and set St_Vars
1179 !-----------------------------------------------------------------------------
1180    num_msc_dm(dm) = cnt
1181    if( cnt > 0 ) then
1182      if( .not. allocated( St_Vars_dm ) ) then
1183        allocate( St_Vars_dm(cnt,n_dom),St_Vars_msk_dm(cnt,n_dom),stat=astat )
1184        if( astat /= 0 ) then
1185          call wrf_error_fatal( 'reg_scan: failed to allocate St_Vars,St_Vars_msk' )
1186        endif
1187      elseif( cnt > maxval(num_msc_dm(1:dm-1)) ) then
1188        n = size( St_vars_dm,dim=1 )
1189        allocate( St_Vars_wrk(n,dm-1),St_Vars_msk_wrk(n,dm-1),stat=astat )
1190        if( astat /= 0 ) then
1191          call wrf_error_fatal( 'reg_scan: failed to allocate St_Vars,St_Vars_msk wrk arrays' )
1192        endif
1193        do dm_ndx = 1,dm-1
1194          do m = 1,n
1195            St_Vars_wrk(m,dm_ndx) = St_Vars_dm(m,dm_ndx)
1196            St_Vars_msk_wrk(m,dm_ndx) = St_Vars_msk_dm(m,dm_ndx)
1197          end do
1198        end do
1199        deallocate( St_vars_dm,St_Vars_msk_dm )
1200        allocate( St_Vars_dm(cnt,n_dom),St_Vars_msk_dm(cnt,n_dom),stat=astat )
1201        if( astat /= 0 ) then
1202          call wrf_error_fatal( 'reg_scan: failed to allocate St_Vars,St_Vars_msk' )
1203        endif
1204        do dm_ndx = 1,dm-1
1205          do m = 1,n
1206            St_Vars_dm(m,dm_ndx) = St_Vars_wrk(m,dm_ndx)
1207            St_Vars_msk_dm(m,dm_ndx) = St_Vars_msk_wrk(m,dm_ndx)
1208          end do
1209        end do
1210        deallocate( St_vars_wrk,St_Vars_msk_wrk )
1211      endif
1212
1213      St_Vars     => St_Vars_dm(:,dm)
1214      St_Vars_msk => St_Vars_msk_dm(:,dm)
1215
1216      St_Vars_msk(:cnt) = .false.
1217      p => grid%head_statevars%next ; cnt = 0
1218
1219      do while( associated(p) )
1220        valid = (p%Type == 'R' .or. p%Type == 'r') .and. &
1221                (p%Ndim == 3 .or. (p%Ndim == 4 .and. p%scalar_array))
1222        if( valid ) then
1223          valid = p%MemoryOrder(:3) == 'XZY' .and. (p%em2 == kde-1 .or. p%em2 == kde)
1224          if( valid ) then
1225            if( p%Ndim == 3 ) then
1226              do m = 1,nVerbotten
1227                if( trim(Verbotten(m)) == trim(p%Varname) ) then
1228                  valid = .false.
1229                  exit
1230                endif
1231              enddo
1232              if( valid ) then
1233                cnt = cnt + 1
1234                St_Vars(cnt)%Varname = p%Varname
1235                St_Vars(cnt)%Description = p%Description
1236                St_Vars(cnt)%Units   = p%Units
1237                St_Vars(cnt)%MemOrd  = p%MemoryOrder
1238                St_Vars(cnt)%Stagger = p%Stagger
1239                St_Vars(cnt)%Ndim = p%Ndim
1240                St_Vars(cnt)%rfield_3d => p%rfield_3d
1241              endif
1242            elseif( p%Ndim == 4 ) then
1243              do n = param_first_scalar,p%num_table(dm)
1244                valid = .true.
1245                do m = 1,nVerbotten
1246                  tstring = trim(Verbotten(m))
1247                  call upcase( tstring )
1248                  if( trim(tstring) == trim(p%dname_table(dm,n)) ) then
1249                    valid = .false.
1250                    exit
1251                  endif
1252                enddo
1253                if( valid ) then
1254                  cnt = cnt + 1
1255                  St_Vars(cnt)%Varname = trim(p%dname_table(dm,n))
1256                  St_Vars(cnt)%Description = p%Description
1257                  St_Vars(cnt)%Units   = p%Units
1258                  St_Vars(cnt)%MemOrd  = p%MemoryOrder
1259                  St_Vars(cnt)%Stagger = p%Stagger
1260                  St_Vars(cnt)%Ndim = 3
1261                  St_Vars(cnt)%rfield_3d => p%rfield_4d(:,:,:,p%index_table(n,dm))
1262                endif
1263              enddo
1264            endif
1265          endif
1266        endif
1267        p => p%next
1268      end do
1269    endif
1270
1271    end subroutine reg_scan
1272
1273    subroutine get_var_cnt( n_vars, vars )
1274
1275 !-----------------------------------------------------------------------------
1276 !  dummy arguments
1277 !-----------------------------------------------------------------------------
1278    integer, intent(inout) :: n_vars
1279    character(len=32), intent(inout)  :: vars(:)
1280
1281    mask(:) = vars(:) /= ' '
1282    wrk_chr(:) = ' '
1283    m1 = 0
1284    do n = 1,var_max
1285      if( mask(n) ) then
1286        m1 = m1 + 1
1287        wrk_chr(m1) = vars(n)
1288      endif
1289    end do
1290
1291    n_vars  = count( wrk_chr(:) /= ' ' )
1292    vars(:) = wrk_chr(:)
1293
1294    end subroutine get_var_cnt
1295
1296    subroutine scan_vars( n_vars, vars, var_ndx, n_tbl, tbl, &
1297                          is_gas, var_type )
1298
1299 !-----------------------------------------------------------------------------
1300 !  dummy arguments
1301 !-----------------------------------------------------------------------------
1302    integer, intent(inout) :: n_vars
1303    integer, intent(in)    :: n_tbl
1304    integer, intent(inout) :: var_ndx(:)
1305    logical, intent(inout) :: is_gas(:)
1306    character(len=*), intent(in)      :: var_type
1307    character(len=32), intent(inout)  :: vars(:)
1308    character(len=256), intent(in)    :: tbl(:)
1309
1310 !-----------------------------------------------------------------------------
1311 !  local variables
1312 !-----------------------------------------------------------------------------
1313    integer :: ml
1314    integer :: ndx(var_max)
1315
1316 var_loop: &
1317    do n = 1,n_vars
1318      var_name = vars(n)
1319      if( trim(var_type) == 'chm' ) then
1320        i = index( '(a)', var_name )
1321        if( i == 0 ) then
1322          i = index( '(A)', var_name )
1323        endif
1324        if( i > 0 ) then
1325          is_gas(n) = .false.
1326          var_name(i:) = ' '
1327          vars(n)(i:)  = ' '
1328        endif
1329      endif
1330      if( trim(var_type) == 'dyn' .or. trim(var_type) == 'msc' ) then
1331        ml = 1
1332      else
1333        ml = param_first_scalar
1334      endif
1335 !    if( trim(var_type) /= 'dyn' ) then
1336 !      ml = param_first_scalar
1337 !    else
1338 !      ml = 1
1339 !    endif
1340      do m1 = ml,n_tbl
1341        if( trim(var_name) == trim(tbl(m1)) ) then
1342          mask(n) = .true.
1343          ndx(n)  = m1
1344          exit
1345        endif
1346      end do
1347      if( .not. mask(n) ) then
1348        write(err_mes,'(''scan_vars('',i2.2,''): '',a,'' not in '',a,'' opt'')') dm,trim(var_name),var_type
1349        call wrf_message( trim(err_mes) )
1350      endif
1351    end do var_loop
1352
1353    if( trim(var_type) == 'dyn' .or. trim(var_type) == 'msc' ) then
1354      ndx(1:n_vars) = ndx(1:n_vars) + offset
1355    endif
1356
1357    wrk_chr(:) = ' '
1358    m1 = 0
1359    do n = 1,var_max
1360      if( mask(n) ) then
1361        m1 = m1 + 1
1362        wrk_chr(m1) = vars(n)
1363      endif
1364    end do
1365
1366    var_ndx(:count(mask)) = pack( ndx(:),mask=mask)
1367    vars(:) = wrk_chr(:)
1368    n_vars  = count( mask )
1369
1370    end subroutine scan_vars
1371
1372    subroutine set_in_dom
1373
1374 !-----------------------------------------------------------------------------
1375 !  local variables
1376 !-----------------------------------------------------------------------------
1377    integer :: ncid
1378    integer :: ios
1379    integer :: time_id
1380    integer :: varid
1381    integer :: time_ndx
1382
1383    real    :: traj_i(n_traj)
1384    real    :: traj_j(n_traj)
1385    real    :: traj_k(n_traj)
1386
1387    character(len=256) :: err_mes
1388    character(len=256) :: filename
1389
1390 include 'netcdf.inc'
1391
1392 !---------------------------------------------------------------------
1393 !  open netcdf file
1394 !---------------------------------------------------------------------
1395    write(filename,'(''wrfout_traj_d'',i2.2)',iostat=ios) dm
1396    if( ios /= 0 ) then
1397      write(err_mes,'(''set_in_dom: failed to set filename: error = '',i6)') ios
1398      call wrf_error_fatal( trim( err_mes  ) )
1399    endif
1400    ios = nf_open( trim(filename), nf_nowrite, ncid )
1401    if( ios /= 0 ) then
1402      write(err_mes,'(''set_in_dom: failed to open '',a,'': error = '',i6)') trim(filename),ios
1403      call wrf_error_fatal( trim( err_mes  ) )
1404    endif
1405
1406 !---------------------------------------------------------------------
1407 !  get current time index
1408 !---------------------------------------------------------------------
1409    err_mes = 'set_in_dom: failed to get time id'
1410    call handle_ncerr( nf_inq_dimid( ncid, 'time', time_id ),trim(err_mes) )
1411    err_mes = 'set_in_dom: failed to get time dimension'
1412    call handle_ncerr( nf_inq_dimlen( ncid, time_id, time_ndx ),trim(err_mes) )
1413
1414 !---------------------------------------------------------------------
1415 !  read in last traj_{i,j,k} variables
1416 !---------------------------------------------------------------------
1417    err_mes = 'set_in_dom: failed to get traj_i id'
1418    call handle_ncerr( nf_inq_varid( ncid, 'traj_i', varid ),trim(err_mes) )
1419    err_mes = 'set_in_dom: failed to get traj_i'
1420    call handle_ncerr( nf_get_vara_real( ncid, varid, (/ 1,time_ndx /), (/ n_traj,1 /), &
1421                                         traj_i ),trim(err_mes) )
1422    err_mes = 'set_in_dom: failed to get traj_j id'
1423    call handle_ncerr( nf_inq_varid( ncid, 'traj_j', varid ),trim(err_mes) )
1424    err_mes = 'set_in_dom: failed to get traj_j'
1425    call handle_ncerr( nf_get_vara_real( ncid, varid, (/ 1,time_ndx /), (/ n_traj,1 /), &
1426                                         traj_j ),trim(err_mes) )
1427    err_mes = 'set_in_dom: failed to get traj_k id'
1428    call handle_ncerr( nf_inq_varid( ncid, 'traj_k', varid ),trim(err_mes) )
1429    err_mes = 'set_in_dom: failed to get traj_k'
1430    call handle_ncerr( nf_get_vara_real( ncid, varid, (/ 1,time_ndx /), (/ n_traj,1 /), &
1431                                         traj_k ),trim(err_mes) )
1432    traj_type(1:n_traj)%in_dom = traj_i(1:n_traj) /= missing_val .and.  traj_j(1:n_traj) /= missing_val &
1433                                                                 .and.  traj_k(1:n_traj) /= missing_val
1434
1435    ios = nf_close( ncid )
1436
1437    end subroutine set_in_dom
1438
1439    end subroutine trajectory_init
1440
1441 #ifdef NETCDF
1442    subroutine trajectory_driver( grid )
1443
1444 #ifdef DM_PARALLEL
1445    use module_dm, only : &
1446                   local_communicator, mytask, ntasks, ntasks_x, ntasks_y                   &
1447                  ,local_communicator_periodic, wrf_dm_max_real, wrf_dm_max_int
1448    use module_comm_dm, only : halo_em_chem_e_3_sub, halo_em_moist_e_3_sub
1449    use module_comm_dm, only : halo_em_tracer_e_3_sub
1450 #endif
1451    use module_domain
1452    use module_date_time
1453    use module_state_description, only : num_chem
1454    use module_state_description, only : num_moist, num_tracer, param_first_scalar
1455
1456 !-----------------------------------------------------------------------------
1457 !  dummy arguments
1458 !-----------------------------------------------------------------------------
1459    type(domain), intent(in) :: grid
1460
1461 !-----------------------------------------------------------------------------
1462 !  local variables
1463 !-----------------------------------------------------------------------------
1464    integer :: ims,ime, jms,jme, kms,kme
1465    integer :: ids,ide, jds,jde, kds,kde
1466    integer :: ips,ipe, jps,jpe, kps,kpe
1467    integer :: dm
1468    integer :: j, k
1469    integer :: il, iu, ios, jl, ju, kl, ku
1470    integer :: m, mu, n, ndx, n_vars, n_traj
1471    integer :: pkg_var_cnt(traj_max)
1472    integer :: ncid
1473    integer :: pkg, trj
1474    integer :: n_msc_buf
1475    integer :: num_chem_sav
1476    integer, pointer :: num_msc      ! number misc species
1477    integer, pointer :: n_dchm       ! number dchm buffer species
1478    integer, pointer :: dchm_buf_ndx(:)
1479    integer :: St_Vars_ndx
1480    integer, allocatable :: St_Vars_buf_ndx(:)
1481 #ifndef DM_PARALLEL
1482    integer :: mytask
1483 #endif
1484    integer :: traj_proc(traj_max), glb_traj_proc(traj_max)
1485    real :: dchm_fill_val(traj_max)
1486    real :: x, y, zi
1487    real :: delsx, delsy, o_delsx, o_delsy
1488    real :: delsz, o_delsz
1489    real :: max_conc
1490    real :: horz_conc(2)
1491    real, pointer :: traj_conc(:)
1492    real, target  :: traj_val(var_max,traj_max)
1493    real, pointer :: wrk4d(:,:,:,:)
1494    real, allocatable, target :: chem(:,:,:,:)
1495    real, allocatable, target :: moist(:,:,:,:)
1496    real, allocatable, target :: tracer(:,:,:,:)
1497    character(len=19)  :: current_timestr, next_timestr
1498    character(len=32)  :: var_name(var_max)
1499    character(len=256) :: err_mes
1500    logical :: has_dchm
1501    logical :: is_root_proc
1502    logical :: is_in_patch_gap
1503    logical :: flsh_buff
1504    logical :: found
1505    logical :: traj_is_active(traj_max)
1506    logical :: pkg_is_active(traj_max,pkg_max)
1507    logical, pointer :: pkg_has_vars(:,:)
1508
1509    logical, external :: wrf_dm_on_monitor
1510
1511    type(WRFU_Time) :: current_time, next_time
1512    type(WRFU_Time) :: start_time, stop_time
1513
1514 include 'netcdf.inc'
1515
1516 #ifndef DM_PARALLEL
1517    mytask = 0
1518 #endif
1519    dm = grid%id
1520    n_traj = traj_cnt(dm)
1521 has_trajectories: &
1522    if( dm_has_traj(dm) .and. n_traj > 0 ) then
1523      St_Vars => St_Vars_dm(:,dm)
1524      St_Vars_msk => St_Vars_msk_dm(:,dm)
1525      num_msc => num_msc_dm(dm)
1526      trjects => traject(:,dm)
1527      n_vals  => n_vals_dm(dm)
1528      is_root_proc = wrf_dm_on_monitor()
1529      if( is_root_proc ) then
1530        trj_pbf => trj_buff(:,dm)
1531      endif
1532      has_dchm = any( trjects(:n_traj)%n_dchm_var > 0 )
1533      if( has_dchm ) then
1534        n_dchm => n_dchm_dm(dm)
1535        dchm_buf_ndx => dchm_buf_ndx_dm(:,dm)
1536      endif
1537
1538      call get_ijk_from_grid( grid ,                   &
1539                              ids, ide, jds, jde, kds, kde,    &
1540                              ims, ime, jms, jme, kms, kme,    &
1541                              ips, ipe, jps, jpe, kps, kpe    )
1542
1543 !-----------------------------------------------------------------------------
1544 !  is trajectory in time interval?
1545 !-----------------------------------------------------------------------------
1546      call domain_clock_get( grid, current_time=current_time, current_timestr=current_timestr)
1547      call geth_newdate( next_timestr, current_timestr, int(grid%dt) )
1548      call wrf_atotime( next_timestr(1:19), next_time )
1549      do trj = 1,n_traj
1550        call wrf_atotime( traject(trj,dm)%start_time(1:19), start_time )
1551        call wrf_atotime( traject(trj,dm)%stop_time(1:19), stop_time )
1552        traj_is_active(trj) = next_time .ge. start_time .and. next_time .le. stop_time
1553      end do
1554 !-----------------------------------------------------------------------------
1555 !  is trajectory in domain?
1556 !-----------------------------------------------------------------------------
1557      pkg_has_vars => pkg_has_vars_dm(:,:,dm)
1558      do trj = 1,n_traj
1559        if( traj_is_active(trj) ) then
1560          trjects(trj)%in_dom = grid%traj_i(trj) >= real(ids) .and. grid%traj_i(trj) <= real(ide-1)
1561          if( trjects(trj)%in_dom ) then
1562            trjects(trj)%in_dom = grid%traj_j(trj) >= real(jds) .and. grid%traj_j(trj) <= real(jde-1)
1563          endif
1564          if( trjects(trj)%in_dom ) then
1565            trjects(trj)%in_dom = grid%traj_k(trj) >= real(kps) .and. grid%traj_k(trj) <= real( min( kpe,kde-1 ) )
1566          endif
1567          traj_is_active(trj) = trjects(trj)%in_dom
1568        endif
1569      end do
1570      do pkg = 1,pkg_max
1571        pkg_is_active(:n_traj,pkg) = traj_is_active(:n_traj) .and. pkg_has_vars(:n_traj,pkg)
1572      end do
1573 !-----------------------------------------------------------------------------
1574 !  check whether this is a chemistry time step
1575 !-----------------------------------------------------------------------------
1576      dchm_fill_val(:n_traj) = missing_val
1577      if( .not. do_chemstep ) then
1578        do trj = 1,n_traj
1579          if( pkg_is_active(trj,dchm_pkg) ) then
1580            dchm_fill_val(trj) = zero_val
1581          endif
1582        end do
1583        pkg_is_active(:n_traj,dchm_pkg) = .false.
1584      endif
1585 !-----------------------------------------------------------------------------
1586 !  is trajectory in mpi process?
1587 !-----------------------------------------------------------------------------
1588      traj_proc(:n_traj) = -1
1589      do trj = 1,n_traj
1590        if( traj_is_active(trj) ) then
1591          trjects(trj)%in_patch = grid%traj_i(trj) >= real(ips) .and. grid%traj_i(trj) <= real( min( ipe,ide-1 ) )
1592          if( trjects(trj)%in_patch ) then
1593            trjects(trj)%in_patch = grid%traj_j(trj) >= real(jps) .and. grid%traj_j(trj) <= real( min( jpe,jde-1 ) )
1594          endif
1595          if( trjects(trj)%in_patch ) then
1596            trjects(trj)%in_patch = grid%traj_k(trj) >= real(kps) .and. grid%traj_k(trj) <= real( min( kpe,kde-1 ) )
1597          endif
1598          if( trjects(trj)%in_patch ) then
1599            traj_proc(trj) = mytask + 1
1600          else
1601            traj_proc(trj) = 0
1602          endif
1603        endif
1604      end do
1605 #ifdef DM_PARALLEL
1606      do trj = 1,n_traj
1607        glb_traj_proc(trj) = wrf_dm_max_int( traj_proc(trj) )
1608      end do
1609 #else
1610      glb_traj_proc(1:n_traj) = traj_proc(1:n_traj)
1611 #endif
1612 !-----------------------------------------------------------------------------
1613 !  any trajectories in "gap" between patches?
1614 !-----------------------------------------------------------------------------
1615      do trj = 1,n_traj
1616        if( traj_is_active(trj) .and. glb_traj_proc(trj) == 0 ) then
1617          trjects(trj)%in_patch = grid%traj_i(trj) >= real(ips) .and. grid%traj_i(trj) <= real( min( ipe+1,ide-1 ) )
1618          if( trjects(trj)%in_patch ) then
1619            trjects(trj)%in_patch = grid%traj_j(trj) >= real(jps) .and. grid%traj_j(trj) <= real( min( jpe+1,jde-1 ) )
1620          endif
1621          if( trjects(trj)%in_patch ) then
1622            trjects(trj)%in_patch = grid%traj_k(trj) >= real(kps) .and. grid%traj_k(trj) <= real( min( kme,kde-1 ) )
1623          endif
1624          if( trjects(trj)%in_patch ) then
1625            traj_proc(trj) = mytask + 1
1626          else
1627            traj_proc(trj) = 0
1628          endif
1629          if( traj_proc(trj) /= 0 ) then
1630            call wrf_debug( 0,'^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^')
1631            write(err_mes,'(''Gapper '',i5,''; x,y,zi = '',1p,3g15.7)') trj,grid%traj_i(trj),grid%traj_j(trj),grid%traj_k(trj)
1632            call wrf_debug( 0,trim(err_mes) )
1633            write(err_mes,'(''Gapper ips,ipe,jps,jpe = '',4i5)') ips,ipe,jps,jpe
1634            call wrf_debug( 0,trim(err_mes) )
1635            call wrf_debug( 0,'^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^')
1636          endif
1637        endif
1638      end do
1639
1640 !-----------------------------------------------------------------------------
1641 !  buffer traj time and position
1642 !-----------------------------------------------------------------------------
1643      if( is_root_proc ) then
1644        n_vals = n_vals + 1
1645        if( grid%itimestep > 0 ) then
1646          trj_pbf(:n_traj)%times(n_vals) = next_timestr
1647        else
1648          trj_pbf(:n_traj)%times(n_vals) = current_timestr
1649        endif
1650        do trj = 1,n_traj
1651          trj_pbf(trj)%trj_i(n_vals) = grid%traj_i(trj)
1652          trj_pbf(trj)%trj_j(n_vals) = grid%traj_j(trj)
1653          trj_pbf(trj)%trj_k(n_vals) = grid%traj_k(trj)
1654          trj_pbf(trj)%trj_lons(n_vals) = grid%traj_long(trj)
1655          trj_pbf(trj)%trj_lats(n_vals) = grid%traj_lat(trj)
1656        end do
1657      endif
1658
1659      do trj = 1,n_traj
1660        traj_val(:,trj) = missing_val
1661      end do
1662
1663 pkg_loop: &
1664      do pkg = 1,pkg_max
1665 pkg_has_active_traj: &
1666        if( any( pkg_is_active(:n_traj,pkg) ) ) then
1667 !-----------------------------------------------------------------------------
1668 !  allocate working data array
1669 !-----------------------------------------------------------------------------
1670          select case( trim(pkg_tag(pkg)) )
1671 #if( WRF_CHEM == 1 )
1672            case( 'chm' )
1673              allocate( chem(ims:ime,kms:kme,jms:jme,num_chem),stat=ios )
1674            case( 'dchm' )
1675              if( n_dchm > 0 ) then
1676                allocate( chem(ims:ime,kms:kme,jms:jme,n_dchm+offset),stat=ios )
1677              endif
1678 #endif
1679            case( 'hyd' )
1680              allocate( moist(ims:ime,kms:kme,jms:jme,num_moist),stat=ios )
1681            case( 'trc' )
1682              allocate( tracer(ims:ime,kms:kme,jms:jme,num_tracer),stat=ios )
1683            case( 'dyn' )
1684              allocate( chem(ims:ime,kms:kme,jms:jme,dyn_max+offset+2),stat=ios )
1685            case( 'msc' )
1686              m = count( St_Vars_msk(:num_msc) )
1687              if( m > 0 ) then
1688                allocate( chem(ims:ime,kms:kme,jms:jme,m+offset), &
1689                          St_Vars_buf_ndx(m+offset),stat=ios )
1690              endif
1691            case default
1692              ios = 0
1693          end select
1694          if( ios /= 0 ) then
1695            write(err_mes,'(''trajectory_driver: failed to allocate wrk4d: error = '',i6)') ios
1696            call wrf_error_fatal( trim( err_mes  ) )
1697          endif
1698          select case( trim(pkg_tag(pkg)) )
1699 #if( WRF_CHEM == 1 )
1700            case( 'chm' )
1701              do m = 1,num_chem
1702                do j = jps,jpe
1703                  do k = kps,kpe
1704                    chem(ips:ipe,k,j,m) = grid%chem(ips:ipe,k,j,m)
1705                  end do
1706                end do
1707              end do
1708            case( 'dchm' )
1709              do m = param_first_scalar,n_dchm+offset
1710                do j = jps,jpe
1711                  do k = kps,kpe
1712                    chem(ips:ipe,k,j,m) = dchm_buff(ips:ipe,k,j,m)
1713                  end do
1714                end do
1715              end do
1716 #endif
1717            case( 'hyd' )
1718              do m = 1,num_moist
1719                do j = jps,jpe
1720                  do k = kps,kpe
1721                    moist(ips:ipe,k,j,m) = grid%moist(ips:ipe,k,j,m)
1722                  end do
1723                end do
1724              end do
1725            case( 'trc' )
1726              do m = 1,num_tracer
1727                do j = jps,jpe
1728                  do k = kps,kpe
1729                    tracer(ips:ipe,k,j,m) = grid%tracer(ips:ipe,k,j,m)
1730                  end do
1731                end do
1732              end do
1733            case( 'dyn' )
1734              call pack_dyn_vals
1735            case( 'msc' )
1736              n_msc_buf = 1
1737              do m = 1,num_msc
1738                if( St_Vars_msk(m) ) then
1739                  n_msc_buf = n_msc_buf + 1
1740                  St_Vars_buf_ndx(n_msc_buf) = m
1741                  do j = jps,jpe
1742                    do k = kps,kpe
1743                      chem(ips:ipe,k,j,n_msc_buf) = St_Vars(m)%rfield_3d(ips:ipe,k,j)
1744                    end do
1745                  end do
1746                endif
1747              end do
1748          end select
1749 #ifdef DM_PARALLEL
1750 !-----------------------------------------------------------------------------
1751 !  any trajectories in "gap" between patches?
1752 !-----------------------------------------------------------------------------
1753              is_in_patch_gap = any( glb_traj_proc(:n_traj) == 0 .and. pkg_is_active(:n_traj,pkg) )
1754              if( is_in_patch_gap ) then
1755 !              write(err_mes,'(''glb_traj_proc mask = '',10i5)') glb_traj_proc(:n_traj)
1756 !              call wrf_debug( 0,trim(err_mes) )
1757                select case( trim(pkg_tag(pkg)) )
1758 #if( WRF_CHEM == 1 )
1759                  case( 'chm' )
1760 #      include "HALO_EM_CHEM_E_3.inc"
1761                  case( 'dchm' )
1762                    num_chem_sav = num_chem
1763                    num_chem     = n_dchm
1764 #      include "HALO_EM_CHEM_E_3.inc"
1765                    num_chem = num_chem_sav
1766 #endif
1767                  case( 'hyd' )
1768 #      include "HALO_EM_MOIST_E_3.inc"
1769                  case( 'trc' )
1770 #      include "HALO_EM_TRACER_E_3.inc"
1771                  case( 'dyn' )
1772                    num_chem_sav = num_chem
1773                    num_chem     = dyn_max + offset + 2
1774 #      include "HALO_EM_CHEM_E_3.inc"
1775                    num_chem = num_chem_sav
1776                  case( 'msc' )
1777                    num_chem_sav = num_chem
1778                    num_chem     = n_msc_buf
1779 #      include "HALO_EM_CHEM_E_3.inc"
1780                    num_chem = num_chem_sav
1781                end select
1782              endif
1783 #endif
1784
1785 traj_loop: &
1786          do trj = 1,n_traj
1787            select case( trim(pkg_tag(pkg)) )
1788              case( 'chm' )
1789                n_vars = traject(trj,dm)%n_chm_var
1790              case( 'hyd' )
1791                n_vars = traject(trj,dm)%n_hyd_var
1792              case( 'trc' )
1793                n_vars = traject(trj,dm)%n_trc_var
1794              case( 'dyn' )
1795                n_vars = traject(trj,dm)%n_dyn_var
1796              case( 'msc' )
1797                n_vars = traject(trj,dm)%n_msc_var
1798              case( 'dchm' )
1799                n_vars = traject(trj,dm)%n_dchm_var
1800            end select
1801 pkg_is_active_in_traj: &
1802            if( pkg_is_active(trj,pkg) ) then
1803              select case( trim(pkg_tag(pkg)) )
1804 #if( WRF_CHEM == 1 )
1805                case( 'chm', 'dchm' )
1806                  wrk4d => chem
1807 #endif
1808                case( 'dyn', 'msc' )
1809                  wrk4d => chem
1810                case( 'hyd' )
1811                  wrk4d => moist
1812                case( 'trc' )
1813                  wrk4d => tracer
1814              end select
1815 in_patch:    if( traj_proc(trj) == mytask+1 ) then
1816                x = grid%traj_i(trj)
1817                y = grid%traj_j(trj)
1818                zi = grid%traj_k(trj)
1819                il = int( x ) ; iu = il + 1
1820                jl = int( y ) ; ju = jl + 1
1821                kl = int( zi ) ; ku = kl + 1
1822                delsx = x - floor(x) ; o_delsx = 1. - delsx
1823                delsy = y - floor(y) ; o_delsy = 1. - delsy
1824                delsz = zi - floor(zi) ; o_delsz = 1. - delsz
1825 var_loop:      do n = 1,n_vars
1826                  found = .true.
1827                  select case( trim(pkg_tag(pkg)) )
1828                    case( 'chm' )
1829                      ndx = traject(trj,dm)%chm_ndx(n)
1830                    case( 'hyd' )
1831                      ndx = traject(trj,dm)%hyd_ndx(n)
1832                    case( 'trc' )
1833                      ndx = traject(trj,dm)%trc_ndx(n)
1834                    case( 'dyn' )
1835                      ndx = 1
1836                      call set_dyn_vals
1837                    case( 'msc' )
1838                      found = .false.
1839                      St_Vars_ndx = trjects(trj)%msc_ndx(n) - offset
1840                      do ndx = param_first_scalar,n_msc_buf
1841                        if( St_Vars_ndx == St_Vars_buf_ndx(ndx) ) then
1842                          found = .true.
1843                          exit
1844                        endif
1845                      end do
1846                      if( found ) then
1847                        call set_msc_vals
1848                      endif
1849                      ndx = 1
1850                    case( 'dchm' )
1851                      ndx = trjects(trj)%dchm_ndx(n)
1852                  end select
1853                  if( found ) then
1854                    horz_conc(1) = o_delsx*o_delsy*wrk4d(il,kl,jl,ndx) + o_delsy*delsx*wrk4d(iu,kl,jl,ndx) &
1855                                 + delsy*o_delsx*wrk4d(il,kl,ju,ndx) + delsx*delsy*wrk4d(iu,kl,ju,ndx)
1856                    horz_conc(2) = o_delsx*o_delsy*wrk4d(il,ku,jl,ndx) + o_delsy*delsx*wrk4d(iu,ku,jl,ndx) &
1857                                 + delsy*o_delsx*wrk4d(il,ku,ju,ndx) + delsx*delsy*wrk4d(iu,ku,ju,ndx)
1858                    traject(trj,dm)%traj_var(n) = delsz*horz_conc(2) + o_delsz*horz_conc(1)
1859                  else
1860                    traject(trj,dm)%traj_var(n) = missing_val
1861                  endif
1862                end do var_loop
1863              else in_patch
1864                traject(trj,dm)%traj_var(:n_vars) = missing_val
1865              endif in_patch
1866              traj_conc => traj_val(:,trj)
1867              do n = 1,n_vars
1868 #ifdef DM_PARALLEL
1869                max_conc = wrf_dm_max_real( traject(trj,dm)%traj_var(n) )
1870 #else
1871                max_conc = traject(trj,dm)%traj_var(n)
1872 #endif
1873                if( is_root_proc ) then
1874                  traj_conc(n) = max_conc
1875                endif
1876              end do
1877            else pkg_is_active_in_traj
1878              if( is_root_proc .and. n_vars > 0 ) then
1879                traj_conc => traj_val(:,trj)
1880                traj_conc(:n_vars) = missing_val
1881              endif
1882            endif pkg_is_active_in_traj
1883 !-----------------------------------------------------------------------------
1884 !  buffer traj chm,trc,hyb,msc,dchm variables
1885 !-----------------------------------------------------------------------------
1886            if( is_root_proc .and. n_vars > 0 ) then
1887              select case( pkg_tag(pkg) )
1888 #if( WRF_CHEM == 1 )
1889                case( 'chm' )
1890                  trj_pbf(trj)%chm_vals(n_vals,:n_vars) = traj_conc(:n_vars)
1891                case( 'dchm' )
1892                  trj_pbf(trj)%dchm_vals(n_vals,:n_vars) = traj_conc(:n_vars)
1893 #endif
1894                case( 'hyd' )
1895                  trj_pbf(trj)%hyd_vals(n_vals,:n_vars) = traj_conc(:n_vars)
1896                case( 'trc' )
1897                  trj_pbf(trj)%trc_vals(n_vals,:n_vars) = traj_conc(:n_vars)
1898                case( 'dyn' )
1899                  trj_pbf(trj)%dyn_vals(n_vals,:n_vars) = traj_conc(:n_vars)
1900                case( 'msc' )
1901                  trj_pbf(trj)%msc_vals(n_vals,:n_vars) = traj_conc(:n_vars)
1902              end select
1903            endif
1904          end do traj_loop
1905          if( allocated( chem ) ) then
1906            deallocate( chem )
1907          endif
1908          if( allocated( moist ) ) then
1909            deallocate( moist )
1910          endif
1911          if( allocated( tracer ) ) then
1912            deallocate( tracer )
1913          endif
1914          if( allocated( St_Vars_buf_ndx ) ) then
1915            deallocate( St_Vars_buf_ndx )
1916          endif
1917          if( pkg == dchm_pkg .and. allocated( dchm_buff ) ) then
1918            deallocate( dchm_buff )
1919          endif
1920        else pkg_has_active_traj
1921          if( is_root_proc ) then
1922            do trj = 1,n_traj
1923              select case( trim(pkg_tag(pkg)) )
1924 #if( WRF_CHEM == 1 )
1925                case( 'chm' )
1926                  n_vars = traject(trj,dm)%n_chm_var
1927                case( 'dchm' )
1928                  n_vars = traject(trj,dm)%n_dchm_var
1929 #endif
1930                case( 'hyd' )
1931                  n_vars = traject(trj,dm)%n_hyd_var
1932                case( 'trc' )
1933                  n_vars = traject(trj,dm)%n_trc_var
1934                case( 'dyn' )
1935                  n_vars = traject(trj,dm)%n_dyn_var
1936                case( 'msc' )
1937                  n_vars = traject(trj,dm)%n_msc_var
1938              end select
1939              if( n_vars > 0 ) then
1940                select case( pkg_tag(pkg) )
1941 #if( WRF_CHEM == 1 )
1942                  case( 'chm' )
1943                    trj_pbf(trj)%chm_vals(n_vals,:n_vars) = missing_val
1944                  case( 'dchm' )
1945                    trj_pbf(trj)%dchm_vals(n_vals,:n_vars) = dchm_fill_val(trj)
1946 #endif
1947                  case( 'hyd' )
1948                    trj_pbf(trj)%hyd_vals(n_vals,:n_vars) = missing_val
1949                  case( 'trc' )
1950                    trj_pbf(trj)%trc_vals(n_vals,:n_vars) = missing_val
1951                  case( 'dyn' )
1952                    trj_pbf(trj)%dyn_vals(n_vals,:n_vars) = missing_val
1953                  case( 'msc' )
1954                    trj_pbf(trj)%msc_vals(n_vals,:n_vars) = missing_val
1955                end select
1956              endif
1957            end do
1958          endif
1959        endif pkg_has_active_traj
1960      end do pkg_loop
1961 !-----------------------------------------------------------------------------
1962 !  output trajectory buffer
1963 !-----------------------------------------------------------------------------
1964      if( is_root_proc ) then
1965        flsh_buff = n_vals == vals_max .or. domain_last_time_step( grid )
1966        if( flsh_buff ) then
1967          call trajectory_write_file( n_traj, n_vals, dm )
1968        endif
1969      endif
1970    endif has_trajectories
1971
1972    CONTAINS
1973
1974    subroutine pack_dyn_vals
1975
1976    integer :: mp1
1977
1978    do m = 1,dyn_max+2
1979      mp1 = m + 1
1980      select case( m )
1981        case( 1 )
1982          do j = jps,jpe
1983            do k = kps,kpe
1984              chem(ips:ipe,k,j,mp1) = grid%p(ips:ipe,k,j)
1985            end do
1986          end do
1987        case( 2 )
1988          do j = jps,jpe
1989            do k = kps,kpe
1990              chem(ips:ipe,k,j,mp1) = grid%t_2(ips:ipe,k,j)
1991            end do
1992          end do
1993        case( 3 )
1994          do j = jps,jpe
1995            do k = kps,kpe
1996              chem(ips:ipe,k,j,mp1) = grid%pb(ips:ipe,k,j)
1997            end do
1998          end do
1999        case( 4 )
2000          do j = jps,jpe
2001            do k = kps,kpe
2002              chem(ips:ipe,k,j,mp1) = grid%u_2(ips:ipe,k,j)
2003            end do
2004          end do
2005        case( 5 )
2006          do j = jps,jpe
2007            do k = kps,kpe
2008              chem(ips:ipe,k,j,mp1) = grid%v_2(ips:ipe,k,j)
2009            end do
2010          end do
2011        case( 6 )
2012          do j = jps,jpe
2013            do k = kps,kpe
2014              chem(ips:ipe,k,j,mp1) = grid%w_2(ips:ipe,k,j)
2015            end do
2016          end do
2017        case( 7 )
2018          do j = jps,jpe
2019            do k = kps,kpe
2020              chem(ips:ipe,k,j,mp1) = grid%ph_2(ips:ipe,k,j)
2021            end do
2022          end do
2023        case( 8 )
2024          do j = jps,jpe
2025            do k = kps,kpe
2026              chem(ips:ipe,k,j,mp1) = grid%phb(ips:ipe,k,j)
2027            end do
2028          end do
2029 #if( WRF_CHEM == 1 )
2030        case( 9 )
2031          do j = jps,jpe
2032            do k = kps,kpe
2033              chem(ips:ipe,k,j,mp1) = grid%wetscav_frcing(ips:ipe,k,j,p_rainprod)
2034            end do
2035          end do
2036        case( 10 )
2037          do j = jps,jpe
2038            do k = kps,kpe
2039              chem(ips:ipe,k,j,mp1) = grid%wetscav_frcing(ips:ipe,k,j,p_evapprod)
2040            end do
2041          end do
2042 #endif
2043      end select
2044    end do
2045
2046    end subroutine pack_dyn_vals
2047
2048    subroutine set_dyn_vals
2049
2050    use module_model_constants, only : g, t0, p1000mb, rcp
2051
2052    integer :: ilp1, iup1, jlp1, jup1, klp1, kup1
2053    real    :: ginv, pinv
2054
2055    select case( traject(trj,dm)%dyn_var(n) )
2056      case( 'p' )
2057        wrk4d(il,kl,jl,1) = chem(il,kl,jl,2) + chem(il,kl,jl,4)
2058        wrk4d(iu,kl,jl,1) = chem(iu,kl,jl,2) + chem(iu,kl,jl,4)
2059        wrk4d(il,ku,jl,1) = chem(il,ku,jl,2) + chem(il,ku,jl,4)
2060        wrk4d(iu,ku,jl,1) = chem(iu,ku,jl,2) + chem(iu,ku,jl,4)
2061        wrk4d(il,kl,ju,1) = chem(il,kl,ju,2) + chem(il,kl,ju,4)
2062        wrk4d(iu,kl,ju,1) = chem(iu,kl,ju,2) + chem(iu,kl,ju,4)
2063        wrk4d(il,ku,ju,1) = chem(il,ku,ju,2) + chem(il,ku,ju,4)
2064        wrk4d(iu,ku,ju,1) = chem(iu,ku,ju,2) + chem(iu,ku,ju,4)
2065      case( 'T' )
2066        wrk4d(il,kl,jl,1) = chem(il,kl,jl,2) + chem(il,kl,jl,4)
2067        wrk4d(iu,kl,jl,1) = chem(iu,kl,jl,2) + chem(iu,kl,jl,4)
2068        wrk4d(il,ku,jl,1) = chem(il,ku,jl,2) + chem(il,ku,jl,4)
2069        wrk4d(iu,ku,jl,1) = chem(iu,ku,jl,2) + chem(iu,ku,jl,4)
2070        wrk4d(il,kl,ju,1) = chem(il,kl,ju,2) + chem(il,kl,ju,4)
2071        wrk4d(iu,kl,ju,1) = chem(iu,kl,ju,2) + chem(iu,kl,ju,4)
2072        wrk4d(il,ku,ju,1) = chem(il,ku,ju,2) + chem(il,ku,ju,4)
2073        wrk4d(iu,ku,ju,1) = chem(iu,ku,ju,2) + chem(iu,ku,ju,4)
2074
2075        pinv = 1./p1000mb
2076        wrk4d(il,kl,jl,1) = (chem(il,kl,jl,3) + t0)*(wrk4d(il,kl,jl,1)*pinv)**rcp
2077        wrk4d(iu,kl,jl,1) = (chem(iu,kl,jl,3) + t0)*(wrk4d(iu,kl,jl,1)*pinv)**rcp
2078        wrk4d(il,ku,jl,1) = (chem(il,ku,jl,3) + t0)*(wrk4d(il,ku,jl,1)*pinv)**rcp
2079        wrk4d(iu,ku,jl,1) = (chem(iu,ku,jl,3) + t0)*(wrk4d(iu,ku,jl,1)*pinv)**rcp
2080        wrk4d(il,kl,ju,1) = (chem(il,kl,ju,3) + t0)*(wrk4d(il,kl,ju,1)*pinv)**rcp
2081        wrk4d(iu,kl,ju,1) = (chem(iu,kl,ju,3) + t0)*(wrk4d(iu,kl,ju,1)*pinv)**rcp
2082        wrk4d(il,ku,ju,1) = (chem(il,ku,ju,3) + t0)*(wrk4d(il,ku,ju,1)*pinv)**rcp
2083        wrk4d(iu,ku,ju,1) = (chem(iu,ku,ju,3) + t0)*(wrk4d(iu,ku,ju,1)*pinv)**rcp
2084      case( 'z' )
2085        klp1 = kl + 1 ; kup1 = ku + 1
2086        ginv = 1./g
2087        wrk4d(il,kl,jl,1) = .5*(chem(il,kl,jl,8) + chem(il,klp1,jl,8) &
2088                                + chem(il,kl,jl,9) + chem(il,klp1,jl,9))*ginv
2089        wrk4d(iu,kl,jl,1) = .5*(chem(iu,kl,jl,8) + chem(iu,klp1,jl,8) &
2090                                + chem(iu,kl,jl,9) + chem(iu,klp1,jl,9))*ginv
2091        wrk4d(il,ku,jl,1) = .5*(chem(il,ku,jl,8) + chem(il,kup1,jl,8) &
2092                                + chem(il,ku,jl,9) + chem(il,kup1,jl,9))*ginv
2093        wrk4d(iu,ku,jl,1) = .5*(chem(iu,ku,jl,8) + chem(iu,kup1,jl,8) &
2094                                + chem(iu,ku,jl,9) + chem(iu,kup1,jl,9))*ginv
2095        wrk4d(il,kl,ju,1) = .5*(chem(il,kl,ju,8) + chem(il,klp1,ju,8) &
2096                                + chem(il,kl,ju,9) + chem(il,klp1,ju,9))*ginv
2097        wrk4d(iu,kl,ju,1) = .5*(chem(iu,kl,ju,8) + chem(iu,klp1,ju,8) &
2098                                + chem(iu,kl,ju,9) + chem(iu,klp1,ju,9))*ginv
2099        wrk4d(il,ku,ju,1) = .5*(chem(il,ku,ju,8) + chem(il,kup1,ju,8) &
2100                                + chem(il,ku,ju,9) + chem(il,kup1,ju,9))*ginv
2101        wrk4d(iu,ku,ju,1) = .5*(chem(iu,ku,ju,8) + chem(iu,kup1,ju,8) &
2102                                + chem(iu,ku,ju,9) + chem(iu,kup1,ju,9))*ginv
2103      case( 'u' )
2104        ilp1 = il + 1 ; iup1 = iu + 1
2105        wrk4d(il,kl,jl,1) = .5*(chem(il,kl,jl,5) + chem(ilp1,kl,jl,5))
2106        wrk4d(iu,kl,jl,1) = .5*(chem(iu,kl,jl,5) + chem(iup1,kl,jl,5))
2107        wrk4d(il,ku,jl,1) = .5*(chem(il,ku,jl,5) + chem(ilp1,ku,jl,5))
2108        wrk4d(iu,ku,jl,1) = .5*(chem(iu,ku,jl,5) + chem(iup1,ku,jl,5))
2109        wrk4d(il,kl,ju,1) = .5*(chem(il,kl,ju,5) + chem(ilp1,kl,ju,5))
2110        wrk4d(iu,kl,ju,1) = .5*(chem(iu,kl,ju,5) + chem(iup1,kl,ju,5))
2111        wrk4d(il,ku,ju,1) = .5*(chem(il,ku,ju,5) + chem(ilp1,ku,ju,5))
2112        wrk4d(iu,ku,ju,1) = .5*(chem(iu,ku,ju,5) + chem(iup1,ku,ju,5))
2113      case( 'v' )
2114        jlp1 = jl + 1 ; jup1 = ju + 1
2115        wrk4d(il,kl,jl,1) = .5*(chem(il,kl,jl,6) + chem(il,kl,jlp1,6))
2116        wrk4d(iu,kl,jl,1) = .5*(chem(iu,kl,jl,6) + chem(iu,kl,jlp1,6))
2117        wrk4d(il,ku,jl,1) = .5*(chem(il,ku,jl,6) + chem(il,ku,jlp1,6))
2118        wrk4d(iu,ku,jl,1) = .5*(chem(iu,ku,jl,6) + chem(iu,ku,jlp1,6))
2119        wrk4d(il,kl,ju,1) = .5*(chem(il,kl,ju,6) + chem(il,kl,jup1,6))
2120        wrk4d(iu,kl,ju,1) = .5*(chem(iu,kl,ju,6) + chem(iu,kl,jup1,6))
2121        wrk4d(il,ku,ju,1) = .5*(chem(il,ku,ju,6) + chem(il,ku,jup1,6))
2122        wrk4d(iu,ku,ju,1) = .5*(chem(iu,ku,ju,6) + chem(iu,ku,jup1,6))
2123      case( 'w' )
2124        klp1 = kl + 1 ; kup1 = ku + 1
2125        wrk4d(il,kl,jl,1) = .5*(chem(il,kl,jl,7) + chem(il,klp1,jl,7))
2126        wrk4d(iu,kl,jl,1) = .5*(chem(iu,kl,jl,7) + chem(iu,klp1,jl,7))
2127        wrk4d(il,ku,jl,1) = .5*(chem(il,ku,jl,7) + chem(il,kup1,jl,7))
2128        wrk4d(iu,ku,jl,1) = .5*(chem(iu,ku,jl,7) + chem(iu,kup1,jl,7))
2129        wrk4d(il,kl,ju,1) = .5*(chem(il,kl,ju,7) + chem(il,klp1,ju,7))
2130        wrk4d(iu,kl,ju,1) = .5*(chem(iu,kl,ju,7) + chem(iu,klp1,ju,7))
2131        wrk4d(il,ku,ju,1) = .5*(chem(il,ku,ju,7) + chem(il,kup1,ju,7))
2132        wrk4d(iu,ku,ju,1) = .5*(chem(iu,ku,ju,7) + chem(iu,kup1,ju,7))
2133 #if( WRF_CHEM == 1 )
2134      case( 'rainprod' )
2135        wrk4d(il,kl,jl,1) = chem(il,kl,jl,10)
2136        wrk4d(iu,kl,jl,1) = chem(iu,kl,jl,10)
2137        wrk4d(il,ku,jl,1) = chem(il,ku,jl,10)
2138        wrk4d(iu,ku,jl,1) = chem(iu,ku,jl,10)
2139        wrk4d(il,kl,ju,1) = chem(il,kl,ju,10)
2140        wrk4d(iu,kl,ju,1) = chem(iu,kl,ju,10)
2141        wrk4d(il,ku,ju,1) = chem(il,ku,ju,10)
2142        wrk4d(iu,ku,ju,1) = chem(iu,ku,ju,10)
2143      case( 'evapprod' )
2144        wrk4d(il,kl,jl,1) = chem(il,kl,jl,11)
2145        wrk4d(iu,kl,jl,1) = chem(iu,kl,jl,11)
2146        wrk4d(il,ku,jl,1) = chem(il,ku,jl,11)
2147        wrk4d(iu,ku,jl,1) = chem(iu,ku,jl,11)
2148        wrk4d(il,kl,ju,1) = chem(il,kl,ju,11)
2149        wrk4d(iu,kl,ju,1) = chem(iu,kl,ju,11)
2150        wrk4d(il,ku,ju,1) = chem(il,ku,ju,11)
2151        wrk4d(iu,ku,ju,1) = chem(iu,ku,ju,11)
2152 #endif
2153    end select
2154
2155    end subroutine set_dyn_vals
2156
2157    subroutine set_msc_vals
2158
2159    integer :: ilp1, iup1, jlp1, jup1, klp1, kup1
2160
2161    select case( trim(St_Vars(St_Vars_ndx)%Stagger) )
2162      case( 'X' )
2163        ilp1 = il + 1 ; iup1 = iu + 1
2164        wrk4d(il,kl,jl,1) = .5*(chem(il,kl,jl,ndx) + chem(ilp1,kl,jl,ndx))
2165        wrk4d(iu,kl,jl,1) = .5*(chem(iu,kl,jl,ndx) + chem(iup1,kl,jl,ndx))
2166        wrk4d(il,ku,jl,1) = .5*(chem(il,ku,jl,ndx) + chem(ilp1,ku,jl,ndx))
2167        wrk4d(iu,ku,jl,1) = .5*(chem(iu,ku,jl,ndx) + chem(iup1,ku,jl,ndx))
2168        wrk4d(il,kl,ju,1) = .5*(chem(il,kl,ju,ndx) + chem(ilp1,kl,ju,ndx))
2169        wrk4d(iu,kl,ju,1) = .5*(chem(iu,kl,ju,ndx) + chem(iup1,kl,ju,ndx))
2170        wrk4d(il,ku,ju,1) = .5*(chem(il,ku,ju,ndx) + chem(ilp1,ku,ju,ndx))
2171        wrk4d(iu,ku,ju,1) = .5*(chem(iu,ku,ju,ndx) + chem(iup1,ku,ju,ndx))
2172      case( 'Y' )
2173        jlp1 = jl + 1 ; jup1 = ju + 1
2174        wrk4d(il,kl,jl,1) = .5*(chem(il,kl,jl,ndx) + chem(il,kl,jlp1,ndx))
2175        wrk4d(iu,kl,jl,1) = .5*(chem(iu,kl,jl,ndx) + chem(iu,kl,jlp1,ndx))
2176        wrk4d(il,ku,jl,1) = .5*(chem(il,ku,jl,ndx) + chem(il,ku,jlp1,ndx))
2177        wrk4d(iu,ku,jl,1) = .5*(chem(iu,ku,jl,ndx) + chem(iu,ku,jlp1,ndx))
2178        wrk4d(il,kl,ju,1) = .5*(chem(il,kl,ju,ndx) + chem(il,kl,jup1,ndx))
2179        wrk4d(iu,kl,ju,1) = .5*(chem(iu,kl,ju,ndx) + chem(iu,kl,jup1,ndx))
2180        wrk4d(il,ku,ju,1) = .5*(chem(il,ku,ju,ndx) + chem(il,ku,jup1,ndx))
2181        wrk4d(iu,ku,ju,1) = .5*(chem(iu,ku,ju,ndx) + chem(iu,ku,jup1,ndx))
2182      case( 'Z' )
2183        klp1 = kl + 1 ; kup1 = ku + 1
2184        wrk4d(il,kl,jl,1) = .5*(chem(il,kl,jl,ndx) + chem(il,klp1,jl,ndx))
2185        wrk4d(iu,kl,jl,1) = .5*(chem(iu,kl,jl,ndx) + chem(iu,klp1,jl,ndx))
2186        wrk4d(il,ku,jl,1) = .5*(chem(il,ku,jl,ndx) + chem(il,kup1,jl,ndx))
2187        wrk4d(iu,ku,jl,1) = .5*(chem(iu,ku,jl,ndx) + chem(iu,kup1,jl,ndx))
2188        wrk4d(il,kl,ju,1) = .5*(chem(il,kl,ju,ndx) + chem(il,klp1,ju,ndx))
2189        wrk4d(iu,kl,ju,1) = .5*(chem(iu,kl,ju,ndx) + chem(iu,klp1,ju,ndx))
2190        wrk4d(il,ku,ju,1) = .5*(chem(il,ku,ju,ndx) + chem(il,kup1,ju,ndx))
2191        wrk4d(iu,ku,ju,1) = .5*(chem(iu,ku,ju,ndx) + chem(iu,kup1,ju,ndx))
2192      case default
2193        wrk4d(il,kl,jl,1) = chem(il,kl,jl,ndx)
2194        wrk4d(iu,kl,jl,1) = chem(iu,kl,jl,ndx)
2195        wrk4d(il,ku,jl,1) = chem(il,ku,jl,ndx)
2196        wrk4d(iu,ku,jl,1) = chem(iu,ku,jl,ndx)
2197        wrk4d(il,kl,ju,1) = chem(il,kl,ju,ndx)
2198        wrk4d(iu,kl,ju,1) = chem(iu,kl,ju,ndx)
2199        wrk4d(il,ku,ju,1) = chem(il,ku,ju,ndx)
2200        wrk4d(iu,ku,ju,1) = chem(iu,ku,ju,ndx)
2201    end select
2202
2203    end subroutine set_msc_vals
2204
2205    end subroutine trajectory_driver
2206
2207    subroutine trajectory_dchm_tstep_init( grid, is_chemstep )
2208 !-----------------------------------------------------------------------------
2209 !  initialize dchm buffer
2210 !-----------------------------------------------------------------------------
2211    use module_domain, only : domain, get_ijk_from_grid
2212    use module_state_description, only : num_chem
2213
2214 !-----------------------------------------------------------------------------
2215 !  dummy arguments
2216 !-----------------------------------------------------------------------------
2217    logical, intent(in)      :: is_chemstep
2218    type(domain), intent(in) :: grid
2219
2220 !-----------------------------------------------------------------------------
2221 !  local variables
2222 !-----------------------------------------------------------------------------
2223    integer :: astat
2224    integer :: j, k, m
2225    integer :: chm_ndx
2226    integer :: dm
2227    integer :: ims,ime, jms,jme, kms,kme
2228    integer :: ids,ide, jds,jde, kds,kde
2229    integer :: ips,ipe, jps,jpe, kps,kpe
2230    integer, pointer :: n_dchm       ! number dchm buffer species
2231    integer, pointer :: dchm_buf_ndx(:)
2232    character(len=256) :: err_mes
2233
2234    do_chemstep = is_chemstep
2235
2236    if( is_chemstep ) then
2237      dm = grid%id
2238      n_dchm => n_dchm_dm(dm)
2239      if( n_dchm > 0 ) then
2240        call get_ijk_from_grid( grid ,                   &
2241                                ids, ide, jds, jde, kds, kde,    &
2242                                ims, ime, jms, jme, kms, kme,    &
2243                                ips, ipe, jps, jpe, kps, kpe    )
2244        if( allocated( dchm_buff ) ) then
2245          deallocate( dchm_buff )
2246        endif
2247        allocate( dchm_buff(ims:ime,kms:kme,jms:jme,n_dchm+offset),stat=astat )
2248        if( astat /= 0 ) then
2249          write(err_mes,'(''trajectory_dchm_tstep_init('',i2.2,''): failed to allocate wrk4d: error = '',i6)') dm,astat
2250          call wrf_error_fatal( trim( err_mes  ) )
2251        endif
2252        dchm_buf_ndx => dchm_buf_ndx_dm(:,dm)
2253        do m = 1,n_dchm
2254          chm_ndx = dchm_buf_ndx(m)
2255          do j = jps,jpe
2256            do k = kps,kpe
2257              dchm_buff(ips:ipe,k,j,m+offset) = grid%chem(ips:ipe,k,j,chm_ndx)
2258            end do
2259          end do
2260        end do
2261      endif
2262    endif
2263
2264    end subroutine trajectory_dchm_tstep_init
2265
2266    subroutine trajectory_dchm_tstep_set( grid )
2267 !-----------------------------------------------------------------------------
2268 !  set dchm buffer
2269 !-----------------------------------------------------------------------------
2270    use module_domain, only : domain, get_ijk_from_grid
2271    use module_state_description, only : num_chem
2272
2273 !-----------------------------------------------------------------------------
2274 !  dummy arguments
2275 !-----------------------------------------------------------------------------
2276    type(domain), intent(in) :: grid
2277
2278 !-----------------------------------------------------------------------------
2279 !  local variables
2280 !-----------------------------------------------------------------------------
2281    integer :: j, k, m, mp1
2282    integer :: chm_ndx
2283    integer :: dm
2284    integer :: ims,ime, jms,jme, kms,kme
2285    integer :: ids,ide, jds,jde, kds,kde
2286    integer :: ips,ipe, jps,jpe, kps,kpe
2287    integer, pointer :: n_dchm       ! number dchm buffer species
2288    integer, pointer :: dchm_buf_ndx(:)
2289
2290    dm = grid%id
2291    n_dchm => n_dchm_dm(dm)
2292    if( n_dchm > 0 ) then
2293      call get_ijk_from_grid( grid ,                   &
2294                              ids, ide, jds, jde, kds, kde,    &
2295                              ims, ime, jms, jme, kms, kme,    &
2296                              ips, ipe, jps, jpe, kps, kpe    )
2297      dchm_buf_ndx => dchm_buf_ndx_dm(:,dm)
2298      do m = 1,n_dchm
2299        mp1 = m + offset
2300        chm_ndx = dchm_buf_ndx(m)
2301        do j = jps,jpe
2302          do k = kps,kpe
2303            dchm_buff(ips:ipe,k,j,mp1) = grid%chem(ips:ipe,k,j,chm_ndx) - dchm_buff(ips:ipe,k,j,mp1)
2304          end do
2305        end do
2306      end do
2307    endif
2308
2309    end subroutine trajectory_dchm_tstep_set
2310
2311    subroutine trajectory_create_file( grid, n_traj )
2312 !-----------------------------------------------------------------------------
2313 !  create trajectory netcdf file
2314 !-----------------------------------------------------------------------------
2315    use module_domain
2316    use module_state_description, only : param_first_scalar, num_chem, num_moist, num_tracer
2317    use module_scalar_tables,     only : chem_dname_table, moist_dname_table, tracer_dname_table
2318
2319 !-----------------------------------------------------------------------------
2320 !  dummy arguments
2321 !-----------------------------------------------------------------------------
2322    integer, intent(in)       :: n_traj
2323    type(domain), intent(in)  :: grid
2324
2325 !-----------------------------------------------------------------------------
2326 !  local variables
2327 !-----------------------------------------------------------------------------
2328    integer :: dm
2329    integer :: ncid, ios
2330    integer :: traj_dim, time_dim, Times_dim
2331    integer :: varid
2332    integer, pointer :: num_msc
2333    integer :: var_dims(2)
2334    integer :: m, n, trj, pkg
2335    character(len=10)  :: coord_name(5) = (/ 'traj_i    ', 'traj_j    ', 'traj_k    ', &
2336                                             'traj_long ', 'traj_lat  ' /)
2337    character(len=256) :: filename
2338    character(len=256) :: var_name
2339    character(len=256) :: err_mes
2340    character(len=256) :: description
2341    character(len=256) :: units
2342
2343    logical, external :: wrf_dm_on_monitor
2344
2345 include 'netcdf.inc'
2346
2347 master_proc: &
2348    if( wrf_dm_on_monitor() ) then
2349      dm = grid%id
2350      write(filename,'(''wrfout_traj_d'',i2.2)',iostat=ios) dm
2351      if( ios /= 0 ) then
2352        write(err_mes,'(''trajectory_create_file: failed to set filename: error = '',i6)') ios
2353        call wrf_error_fatal( trim( err_mes  ) )
2354      endif
2355 !    ios = nf_create( trim(filename), or(nf_clobber,nf_netcdf4), ncid )
2356      ios = nf_create( trim(filename), nf_clobber, ncid )
2357      if( ios /= nf_noerr ) then
2358        write(err_mes,'(''trajectory_create_file: failed to create '',a,'': error = '',i6)') trim(filename),ios
2359        call wrf_error_fatal( trim( err_mes  ) )
2360      endif
2361 !-----------------------------------------------------------------------------
2362 !  define dimensions
2363 !-----------------------------------------------------------------------------
2364      err_mes = 'trajectory_create_file: failed to create traj dimension'
2365      call handle_ncerr( nf_def_dim( ncid, 'traj', n_traj, traj_dim ), trim(err_mes) )
2366      err_mes = 'trajectory_create_file: failed to create time dimension'
2367      call handle_ncerr( nf_def_dim( ncid, 'time', nf_unlimited, time_dim ), trim(err_mes) )
2368      err_mes = 'trajectory_create_file: failed to create Times dimension'
2369      call handle_ncerr( nf_def_dim( ncid, 'DateStrLen', 19, Times_dim ), trim(err_mes) )
2370 !-----------------------------------------------------------------------------
2371 !  define variables
2372 !-----------------------------------------------------------------------------
2373      var_dims(:) = (/ Times_dim,time_dim /)
2374      err_mes = 'trajectory_create_file: failed to create Times variable'
2375      call handle_ncerr( nf_def_var( ncid, 'Times', nf_char, 2, var_dims, varid ), trim(err_mes) )
2376
2377 !-----------------------------------------------------------------------------
2378 !  first the coordinate variables
2379 !-----------------------------------------------------------------------------
2380      var_dims(:) = (/ traj_dim,time_dim /)
2381      do m = 1,5
2382        var_name = coord_name(m)
2383        err_mes = 'trajectory_create_file: failed to create ' // trim(var_name) // ' variable'
2384        call handle_ncerr( nf_def_var( ncid, trim(var_name), nf_real, 2, var_dims, varid ), trim(err_mes) )
2385      end do
2386 !-----------------------------------------------------------------------------
2387 !  then the species variables
2388 !-----------------------------------------------------------------------------
2389      num_msc => num_msc_dm(dm)
2390 pgk_loop: &
2391      do pkg = 1,pkg_max
2392        select case( pkg_tag(pkg) )
2393          case( 'chm' )
2394            trj_msk => trj_msk_dm(:,:,chm_pkg,dm)
2395            if( any( trj_msk(:n_traj,1) ) ) then
2396              call def_vars( 'chm' )
2397            endif
2398          case( 'hyd' )
2399            trj_msk => trj_msk_dm(:,:,hyd_pkg,dm)
2400            if( any( trj_msk(:n_traj,1) ) ) then
2401              call def_vars( 'hyd' )
2402            endif
2403          case( 'trc' )
2404            trj_msk => trj_msk_dm(:,:,trc_pkg,dm)
2405            if( any( trj_msk(:n_traj,1) ) ) then
2406              call def_vars( 'trc' )
2407            endif
2408          case( 'dyn' )
2409            trj_msk => trj_msk_dm(:,:,dyn_pkg,dm)
2410            if( any( trj_msk(:n_traj,1) ) ) then
2411              call def_vars( 'dyn' )
2412            endif
2413          case( 'msc' )
2414            trj_msk => trj_msk_dm(:,:,msc_pkg,dm)
2415            if( any( trj_msk(:n_traj,1) ) ) then
2416              call def_vars( 'msc' )
2417            endif
2418          case( 'dchm' )
2419            trj_msk => trj_msk_dm(:,:,dchm_pkg,dm)
2420            if( any( trj_msk(:n_traj,1) ) ) then
2421              call def_vars( 'dchm' )
2422            endif
2423        end select
2424      end do pgk_loop
2425
2426      err_mes = 'trajectory_create_file: failed to end definition for file ' // trim(filename)
2427      call handle_ncerr( nf_enddef( ncid ), trim(err_mes) )
2428      err_mes = 'trajectory_create_file: failed to close file ' // trim(filename)
2429      call handle_ncerr( nf_close( ncid ), trim(err_mes) )
2430    endif master_proc
2431
2432    CONTAINS
2433
2434    subroutine def_vars( var_type )
2435
2436    character(len=*), intent(in)  :: var_type
2437
2438    integer :: m, ndx, trj
2439    integer, pointer  :: n_dchm
2440    character(len=32) :: spc_name
2441
2442    select case( var_type )
2443      case( 'chm' )
2444        trj_msk => trj_msk_dm(:,:,chm_pkg,dm)
2445        do n = param_first_scalar,num_chem
2446          if( any( trj_msk(:n_traj,n) ) ) then
2447            spc_name = chem_dname_table(dm,n)
2448            write(var_name,'(a,''_traj'')') trim(spc_name)
2449            err_mes = 'def_vars: failed to create ' // trim(var_name) // ' variable'
2450            call handle_ncerr( nf_def_var( ncid, trim(var_name), nf_real, 2, var_dims, varid ), trim(err_mes) )
2451            description = trim(var_name) // ' mixing ratio'
2452            call handle_ncerr( nf_put_att_text( ncid, varid, 'description', len_trim(description), description ), trim(err_mes) )
2453            units = 'ppmv'
2454 trj_loop:  do trj = 1,n_traj
2455              if( trj_msk(trj,n) ) then
2456                do m = 1,trjects(trj)%n_chm_var
2457                  if( trim(trjects(trj)%chm_spc(m)) == trim(spc_name) ) then
2458                    if( .not. trjects(trj)%chm_is_gas(m) ) then
2459                      units = 'ug/kg-dryair'
2460                    endif
2461                    exit trj_loop
2462                  endif
2463                end do
2464                exit trj_loop
2465              endif
2466            end do trj_loop
2467            call handle_ncerr( nf_put_att_text( ncid, varid, 'units', len_trim(units), units ), trim(err_mes) )
2468          endif
2469        end do
2470      case( 'hyd' )
2471        trj_msk => trj_msk_dm(:,:,hyd_pkg,dm)
2472        do n = param_first_scalar,num_moist
2473          if( any( trj_msk(:n_traj,n) ) ) then
2474            write(var_name,'(a,''_traj'')') trim(moist_dname_table(dm,n))
2475            err_mes = 'def_vars: failed to create ' // trim(var_name) // ' variable'
2476            call handle_ncerr( nf_def_var( ncid, trim(var_name), nf_real, 2, var_dims, varid ), trim(err_mes) )
2477            description = trim(var_name) // ' mixing ratio'
2478            call handle_ncerr( nf_put_att_text( ncid, varid, 'description', len_trim(description), description ), trim(err_mes) )
2479            units = 'ug/kg-dryair'
2480            call handle_ncerr( nf_put_att_text( ncid, varid, 'units', len_trim(units), units ), trim(err_mes) )
2481          endif
2482        end do
2483      case( 'trc' )
2484        trj_msk => trj_msk_dm(:,:,trc_pkg,dm)
2485        do n = param_first_scalar,num_tracer
2486          if( any( trj_msk(:n_traj,n) ) ) then
2487            write(var_name,'(a,''_traj'')') trim(tracer_dname_table(dm,n))
2488            err_mes = 'def_vars: failed to create ' // trim(var_name) // ' variable'
2489            call handle_ncerr( nf_def_var( ncid, trim(var_name), nf_real, 2, var_dims, varid ), trim(err_mes) )
2490            description = trim(var_name) // ' mixing ratio'
2491            call handle_ncerr( nf_put_att_text( ncid, varid, 'description', len_trim(description), description ), trim(err_mes) )
2492            units = ' '
2493            call handle_ncerr( nf_put_att_text( ncid, varid, 'units', len_trim(units), units ), trim(err_mes) )
2494          endif
2495        end do
2496      case( 'dyn' )
2497        trj_msk => trj_msk_dm(:,:,dyn_pkg,dm)
2498        do n = param_first_scalar,dyn_max + offset
2499          if( any( trj_msk(:n_traj,n) ) ) then
2500            write(var_name,'(a,''_traj'')') trim(dyn_var_lst(n-1))
2501            err_mes = 'def_vars: failed to create ' // trim(var_name) // ' variable'
2502            call handle_ncerr( nf_def_var( ncid, trim(var_name), nf_real, 2, var_dims, varid ), trim(err_mes) )
2503            description = trim(dyn_var_desc_att(n-1))
2504            call handle_ncerr( nf_put_att_text( ncid, varid, 'description', len_trim(description), description ), trim(err_mes) )
2505            units = trim(dyn_var_unit_att(n-1))
2506            call handle_ncerr( nf_put_att_text( ncid, varid, 'units', len_trim(units), units ), trim(err_mes) )
2507          endif
2508        end do
2509      case( 'msc' )
2510        trj_msk => trj_msk_dm(:,:,msc_pkg,dm)
2511        do n = param_first_scalar,num_msc + offset
2512          if( any( trj_msk(:n_traj,n) ) ) then
2513            write(var_name,'(a,''_traj'')') trim(St_Vars(n-1)%Varname)
2514            err_mes = 'def_vars: failed to create ' // trim(var_name) // ' variable'
2515            call handle_ncerr( nf_def_var( ncid, trim(var_name), nf_real, 2, var_dims, varid ), trim(err_mes) )
2516            description = trim(St_Vars(n-1)%Description)
2517            call handle_ncerr( nf_put_att_text( ncid, varid, 'description', len_trim(description), description ), trim(err_mes) )
2518            units = trim(St_Vars(n-1)%Units)
2519            call handle_ncerr( nf_put_att_text( ncid, varid, 'units', len_trim(units), units ), trim(err_mes) )
2520          endif
2521        end do
2522      case( 'dchm' )
2523        trj_msk => trj_msk_dm(:,:,dchm_pkg,dm)
2524        n_dchm  => n_dchm_dm(dm)
2525        dchm_buf_ndx => dchm_buf_ndx_dm(:,dm)
2526        do n = 1,n_dchm
2527          if( any( trj_msk(:n_traj,n+offset) ) ) then
2528            ndx = dchm_buf_ndx(n)
2529            spc_name = 'dchm_' // trim(chem_dname_table(dm,ndx))
2530            write(var_name,'(a,''_traj'')') trim(spc_name)
2531            err_mes = 'def_vars: failed to create ' // trim(var_name) // ' variable'
2532            call handle_ncerr( nf_def_var( ncid, trim(var_name), nf_real, 2, var_dims, varid ), trim(err_mes) )
2533            description = trim(var_name) // ' mixing ratio'
2534            call handle_ncerr( nf_put_att_text( ncid, varid, 'description', len_trim(description), description ), trim(err_mes) )
2535            units = 'ppmv'
2536            call handle_ncerr( nf_put_att_text( ncid, varid, 'units', len_trim(units), units ), trim(err_mes) )
2537          endif
2538        end do
2539    end select
2540
2541    end subroutine def_vars
2542
2543    end subroutine trajectory_create_file
2544
2545    subroutine trajectory_write_file( n_traj, n_vals, dm )
2546 !-----------------------------------------------------------------------------
2547 !  create trajectory netcdf file
2548 !-----------------------------------------------------------------------------
2549    use module_domain
2550    use module_state_description, only : param_first_scalar, num_chem, num_moist, num_tracer
2551    use module_scalar_tables,     only : chem_dname_table, moist_dname_table, tracer_dname_table
2552
2553 !-----------------------------------------------------------------------------
2554 !  dummy arguments
2555 !-----------------------------------------------------------------------------
2556    integer, intent(in)        :: n_traj
2557    integer, intent(inout)     :: n_vals
2558    integer, intent(in)        :: dm
2559
2560 !-----------------------------------------------------------------------------
2561 !  local variables
2562 !-----------------------------------------------------------------------------
2563    integer :: ncid
2564    integer :: astat, ios
2565    integer :: time_id
2566    integer :: varid
2567    integer :: l, m, n, trj, pkg, spc, spcp1
2568    integer :: time_ndx
2569    integer :: buf_ndx
2570    integer :: ndx
2571    integer, pointer :: num_msc      ! number misc species
2572    integer, pointer :: n_dchm       ! total number of dchm species in domain
2573    real, allocatable :: holder(:,:)
2574    character(len=10)  :: coord_name(5) = (/ 'traj_i    ', 'traj_j    ', 'traj_k    ', &
2575                                             'traj_long ', 'traj_lat  ' /)
2576    character(len=256) :: var_name
2577    character(len=256) :: err_mes
2578    character(len=256) :: filename
2579    character(len=256) :: spcname
2580
2581    logical :: found
2582
2583 include 'netcdf.inc'
2584
2585 !---------------------------------------------------------------------
2586 !  open netcdf file
2587 !---------------------------------------------------------------------
2588    write(filename,'(''wrfout_traj_d'',i2.2)',iostat=ios) dm
2589    if( ios /= 0 ) then
2590      write(err_mes,'(''trajectory_write_file: failed to set filename: error = '',i6)') ios
2591      call wrf_error_fatal( trim( err_mes  ) )
2592    endif
2593    ios = nf_open( trim(filename), nf_write, ncid )
2594    if( ios /= 0 ) then
2595      write(err_mes,'(''trajectory_write_file: failed to open '',a,'': error = '',i6)') trim(filename),ios
2596      call wrf_error_fatal( trim( err_mes  ) )
2597    endif
2598
2599 !---------------------------------------------------------------------
2600 !  allocate wrking array
2601 !---------------------------------------------------------------------
2602    allocate( holder(n_traj,n_vals),stat=astat )
2603    if( astat /= 0 ) then
2604      write(err_mes,'(''trajectory_write_file: failed to allocate holder; error = '',i6)') astat
2605      call wrf_error_fatal( trim( err_mes  ) )
2606    endif
2607
2608 !---------------------------------------------------------------------
2609 !  get current time index
2610 !---------------------------------------------------------------------
2611    err_mes = 'trajectory_write_file: failed to get time id'
2612    call handle_ncerr( nf_inq_dimid( ncid, 'time', time_id ),trim(err_mes) )
2613    err_mes = 'trajectory_write_file: failed to get time dimension'
2614    call handle_ncerr( nf_inq_dimlen( ncid, time_id, time_ndx ),trim(err_mes) )
2615    time_ndx = time_ndx + 1
2616
2617 !---------------------------------------------------------------------
2618 !  write out trajectory times
2619 !---------------------------------------------------------------------
2620    err_mes = 'trajectory_write_file: failed to get Times id'
2621    call handle_ncerr( nf_inq_varid( ncid, 'Times', varid ),trim(err_mes) )
2622    err_mes = 'trajectory_write_file: failed to write Times'
2623    call handle_ncerr( nf_put_vara_text( ncid, varid, (/ 1,time_ndx /), (/ 19,n_vals /), trj_pbf(1)%times(:n_vals) ), trim(err_mes) )
2624
2625 !---------------------------------------------------------------------
2626 !  write out trajectory coordinates
2627 !---------------------------------------------------------------------
2628 coord_loop: &
2629    do l = 1,5
2630      var_name = coord_name(l)
2631      err_mes = 'trajectory_write_file: failed to get '// trim(var_name) // ' id'
2632      call handle_ncerr( nf_inq_varid( ncid, trim(var_name), varid ),trim(err_mes) )
2633      select case( l )
2634        case( 1 )
2635          do n = 1,n_traj
2636            holder(n,:n_vals) = trj_pbf(n)%trj_i(:n_vals)
2637          end do
2638        case( 2 )
2639          do n = 1,n_traj
2640            holder(n,:n_vals) = trj_pbf(n)%trj_j(:n_vals)
2641          end do
2642        case( 3 )
2643          do n = 1,n_traj
2644            holder(n,:n_vals) = trj_pbf(n)%trj_k(:n_vals)
2645          end do
2646        case( 4 )
2647          do n = 1,n_traj
2648            holder(n,:n_vals) = trj_pbf(n)%trj_lons(:n_vals)
2649          end do
2650        case( 5 )
2651          do n = 1,n_traj
2652            holder(n,:n_vals) = trj_pbf(n)%trj_lats(:n_vals)
2653          end do
2654        end select
2655        err_mes = 'trajectory_write_file: failed to write ' // trim(var_name)
2656        call handle_ncerr( nf_put_vara_real( ncid, varid, (/ 1,time_ndx /), (/ n_traj,n_vals /), &
2657                                             holder ), trim(err_mes) )
2658    end do coord_loop
2659
2660    St_Vars => St_Vars_dm(:,dm)
2661    St_Vars_msk => St_Vars_msk_dm(:,dm)
2662    num_msc => num_msc_dm(dm)
2663 !---------------------------------------------------------------------
2664 !  write out trajectory variables
2665 !---------------------------------------------------------------------
2666 pkg_loop: &
2667    do pkg = 1,pkg_max
2668      select case( pkg_tag(pkg) )
2669        case( 'chm' )
2670          trj_msk => trj_msk_dm(:,:,chm_pkg,dm)
2671          do spc = param_first_scalar,num_chem
2672            if( any( trj_msk(:n_traj,spc) ) ) then
2673              holder(:,:) = missing_val
2674              do trj = 1,n_traj
2675                if( trj_msk(trj,spc) ) then
2676                  buf_ndx = get_spc_buf_ndx( trjects(trj)%n_chm_var, trjects(trj)%chm_spc, chem_dname_table(dm,spc) )
2677                  if( buf_ndx > 0 ) then
2678                    holder(trj,:n_vals) = trj_pbf(trj)%chm_vals(:n_vals,buf_ndx)
2679                  endif
2680                endif
2681              end do
2682              write(var_name,'(a,''_traj'')') trim(chem_dname_table(dm,spc))
2683              err_mes = 'trajectory_write_file: failed to get '// trim(var_name) // ' id'
2684              call handle_ncerr( nf_inq_varid( ncid, trim(var_name), varid ),trim(err_mes) )
2685              err_mes = 'trajectory_write_file: failed to write ' // trim(var_name)
2686              call handle_ncerr( nf_put_vara_real( ncid, varid, (/ 1,time_ndx /), (/ n_traj,n_vals /), &
2687                                                   holder ),trim(err_mes) )
2688            endif
2689          end do
2690        case( 'hyd' )
2691          trj_msk => trj_msk_dm(:,:,hyd_pkg,dm)
2692          do spc = param_first_scalar,num_moist
2693            if( any( trj_msk(:n_traj,spc) ) ) then
2694              holder(:,:) = missing_val
2695              do trj = 1,n_traj
2696                if( trj_msk(trj,spc) ) then
2697                  buf_ndx = get_spc_buf_ndx( trjects(trj)%n_hyd_var, trjects(trj)%hyd_spc, moist_dname_table(dm,spc) )
2698                  if( buf_ndx > 0 ) then
2699                    holder(trj,:n_vals) = trj_pbf(trj)%hyd_vals(:n_vals,buf_ndx)
2700                  endif
2701                endif
2702              end do
2703              write(var_name,'(a,''_traj'')') trim(moist_dname_table(dm,spc))
2704              err_mes = 'trajectory_write_file: failed to get '// trim(var_name) // ' id'
2705              call handle_ncerr( nf_inq_varid( ncid, trim(var_name), varid ),trim(err_mes) )
2706              err_mes = 'trajectory_write_file: failed to write ' // trim(var_name)
2707              call handle_ncerr( nf_put_vara_real( ncid, varid, (/ 1,time_ndx /), (/ n_traj,n_vals /), &
2708                                                   holder ),trim(err_mes) )
2709            endif
2710          end do
2711        case( 'trc' )
2712          trj_msk => trj_msk_dm(:,:,trc_pkg,dm)
2713          do spc = param_first_scalar,num_tracer
2714            if( any( trj_msk(:n_traj,spc) ) ) then
2715              holder(:,:) = missing_val
2716              do trj = 1,n_traj
2717                if( trj_msk(trj,spc) ) then
2718                  buf_ndx = get_spc_buf_ndx( trjects(trj)%n_trc_var, trjects(trj)%trc_spc, tracer_dname_table(dm,spc) )
2719                  if( buf_ndx > 0 ) then
2720                    holder(trj,:n_vals) = trj_pbf(trj)%trc_vals(:n_vals,buf_ndx)
2721                  endif
2722                endif
2723              end do
2724              write(var_name,'(a,''_traj'')') trim(tracer_dname_table(dm,spc))
2725              err_mes = 'trajectory_write_file: failed to get '// trim(var_name) // ' id'
2726              call handle_ncerr( nf_inq_varid( ncid, trim(var_name), varid ),trim(err_mes) )
2727              err_mes = 'trajectory_write_file: failed to write ' // trim(var_name)
2728              call handle_ncerr( nf_put_vara_real( ncid, varid, (/ 1,time_ndx /), (/ n_traj,n_vals /), &
2729                                                   holder ),trim(err_mes) )
2730            endif
2731          end do
2732        case( 'dyn' )
2733          trj_msk => trj_msk_dm(:,:,dyn_pkg,dm)
2734          do spc = param_first_scalar,dyn_max+offset
2735            if( any( trj_msk(:n_traj,spc) ) ) then
2736              holder(:,:) = missing_val
2737              do trj = 1,n_traj
2738                if( trj_msk(trj,spc) ) then
2739                  buf_ndx = get_spc_buf_ndx( trjects(trj)%n_dyn_var, trjects(trj)%dyn_var, dyn_var_lst(spc-offset) )
2740                  if( buf_ndx > 0 ) then
2741                    holder(trj,:n_vals) = trj_pbf(trj)%dyn_vals(:n_vals,buf_ndx)
2742                  endif
2743                endif
2744              end do
2745              write(var_name,'(a,''_traj'')') trim(dyn_var_lst(spc-offset))
2746              err_mes = 'trajectory_write_file: failed to get '// trim(var_name) // ' id'
2747              call handle_ncerr( nf_inq_varid( ncid, trim(var_name), varid ),trim(err_mes) )
2748              err_mes = 'trajectory_write_file: failed to write ' // trim(var_name)
2749              call handle_ncerr( nf_put_vara_real( ncid, varid, (/ 1,time_ndx /), (/ n_traj,n_vals /), &
2750                                                   holder ),trim(err_mes) )
2751            endif
2752          end do
2753        case( 'msc' )
2754          trj_msk => trj_msk_dm(:,:,msc_pkg,dm)
2755          do spc = param_first_scalar,num_msc+offset
2756            if( any( trj_msk(:n_traj,spc) ) ) then
2757              holder(:,:) = missing_val
2758              do trj = 1,n_traj
2759                if( trj_msk(trj,spc) ) then
2760                  found = .false.
2761                  do buf_ndx = 1,traject(trj,dm)%n_msc_var
2762                    if( traject(trj,dm)%msc_ndx(buf_ndx) == spc ) then
2763                      found = .true.
2764                      exit
2765                    endif
2766                  end do
2767                  if( found ) then
2768                    holder(trj,:n_vals) = trj_pbf(trj)%msc_vals(:n_vals,buf_ndx)
2769                  endif
2770                endif
2771              end do
2772              write(var_name,'(a,''_traj'')') trim(St_Vars(spc-offset)%Varname)
2773              err_mes = 'trajectory_write_file: failed to get '// trim(var_name) // ' id'
2774              call handle_ncerr( nf_inq_varid( ncid, trim(var_name), varid ),trim(err_mes) )
2775              err_mes = 'trajectory_write_file: failed to write ' // trim(var_name)
2776              call handle_ncerr( nf_put_vara_real( ncid, varid, (/ 1,time_ndx /), (/ n_traj,n_vals /), &
2777                                                   holder ),trim(err_mes) )
2778            endif
2779          end do
2780        case( 'dchm' )
2781          dchm_buf_ndx => dchm_buf_ndx_dm(:,dm)
2782          trj_msk => trj_msk_dm(:,:,dchm_pkg,dm)
2783          n_dchm  => n_dchm_dm(dm)
2784          do spc = 1,n_dchm
2785            spcp1 = spc + 1
2786            if( any( trj_msk(:n_traj,spcp1) ) ) then
2787              holder(:,:) = missing_val
2788              do trj = 1,n_traj
2789                if( trj_msk(trj,spcp1) ) then
2790                  buf_ndx = get_dchm_buf_ndx( trjects(trj)%n_dchm_var, trjects(trj)%dchm_ndx, spcp1 )
2791                  if( buf_ndx > 0 ) then
2792                    holder(trj,:n_vals) = trj_pbf(trj)%dchm_vals(:n_vals,buf_ndx)
2793                  endif
2794                endif
2795              end do
2796              ndx = dchm_buf_ndx(spc)
2797              var_name = 'dchm_' // trim(chem_dname_table(dm,ndx)) // '_traj'
2798 !            write(var_name,'(a,''_traj'')') trim(chem_dname_table(dm,ndx))
2799              err_mes = 'trajectory_write_file: failed to get '// trim(var_name) // ' id'
2800              call handle_ncerr( nf_inq_varid( ncid, trim(var_name), varid ),trim(err_mes) )
2801              err_mes = 'trajectory_write_file: failed to write ' // trim(var_name)
2802              call handle_ncerr( nf_put_vara_real( ncid, varid, (/ 1,time_ndx /), (/ n_traj,n_vals /), &
2803                                                   holder ),trim(err_mes) )
2804            endif
2805          end do
2806      end select
2807    end do pkg_loop
2808
2809    n_vals = 0
2810
2811    ios = nf_close( ncid )
2812
2813    if( allocated( holder ) ) then
2814      deallocate( holder )
2815    endif
2816
2817    end subroutine trajectory_write_file
2818
2819    integer function get_spc_buf_ndx( ncnt, list, match_name )
2820
2821    integer, intent(in)          :: ncnt
2822    character(len=*), intent(in) :: match_name
2823    character(len=*), intent(in) :: list(:)
2824
2825    integer :: spc
2826
2827    get_spc_buf_ndx = -1
2828    do spc = 1,ncnt
2829      if( trim(match_name) == trim(list(spc)) ) then
2830        get_spc_buf_ndx = spc
2831        exit
2832      endif
2833    end do
2834
2835    end function get_spc_buf_ndx
2836
2837    integer function get_dchm_buf_ndx( ncnt, list, match_ndx )
2838
2839    integer, intent(in)          :: ncnt
2840    integer, intent(in)          :: match_ndx
2841    integer, intent(in)          :: list(:)
2842
2843    integer :: spc
2844
2845    get_dchm_buf_ndx = -1
2846    do spc = 1,ncnt
2847      if( match_ndx == list(spc) ) then
2848        get_dchm_buf_ndx = spc
2849        exit
2850      endif
2851    end do
2852
2853    end function get_dchm_buf_ndx
2854
2855    subroutine handle_ncerr( ret, mes )
2856 !---------------------------------------------------------------------
2857 !       ... netcdf error handling routine
2858 !---------------------------------------------------------------------
2859 !---------------------------------------------------------------------
2860 !       ... dummy arguments
2861 !---------------------------------------------------------------------
2862    integer, intent(in) :: ret
2863    character(len=*), intent(in) :: mes
2864
2865 include 'netcdf.inc'
2866
2867    if( ret /= nf_noerr ) then
2868       call wrf_message( trim(mes) )
2869       call wrf_message( trim(nf_strerror(ret)) )
2870       call wrf_abort
2871    endif
2872
2873    end subroutine handle_ncerr
2874 #endif
2875
2876    subroutine trajmapproj (grid,config_flags,ts_proj)
2877
2878    use module_domain
2879    use module_llxy
2880    use module_configure, only : grid_config_rec_type, model_config_rec
2881    use module_dm, only : wrf_dm_min_real
2882
2883    IMPLICIT NONE
2884
2885
2886 !------------------------------------------------------------------------
2887 ! Arguments
2888 !------------------------------------------------------------------------
2889    TYPE(domain), INTENT(IN) :: grid
2890    TYPE(grid_config_rec_type) , INTENT(IN)  :: config_flags
2891    TYPE(PROJ_INFO), INTENT(out) :: ts_proj
2892
2893 !------------------------------------------------------------------------
2894 ! Local variables
2895 !------------------------------------------------------------------------
2896    REAL :: ts_rx, ts_ry, ts_xlat, ts_xlong, ts_hgt
2897    REAL :: known_lat, known_lon
2898
2899    INTEGER :: ids, ide, jds, jde, kds, kde,        &
2900               ims, ime, jms, jme, kms, kme,        &
2901               ips, ipe, jps, jpe, kps, kpe
2902
2903    TYPE (grid_config_rec_type)               :: config_flags_temp
2904
2905    config_flags_temp = config_flags
2906
2907    call get_ijk_from_grid ( grid ,                               &
2908                             ids, ide, jds, jde, kds, kde,        &
2909                             ims, ime, jms, jme, kms, kme,        &
2910                             ips, ipe, jps, jpe, kps, kpe )
2911
2912    call model_to_grid_config_rec ( grid%id , model_config_rec , config_flags_temp )
2913
2914
2915 !------------------------------------------------------------------------
2916 ! Set up map transformation structure
2917 !------------------------------------------------------------------------
2918    call map_init( ts_proj )
2919
2920    IF (ips <= 1 .AND. 1 <= ipe .AND. jps <= 1 .AND. 1 <= jpe) THEN
2921       known_lat = grid%xlat(1,1)
2922       known_lon = grid%xlong(1,1)
2923    ELSE
2924       known_lat = 9999.
2925       known_lon = 9999.
2926    END IF
2927    known_lat = wrf_dm_min_real(known_lat)
2928    known_lon = wrf_dm_min_real(known_lon)
2929
2930
2931    select case( config_flags%map_proj )
2932 !------------------------------------------------------------------------
2933 ! Mercator
2934 !------------------------------------------------------------------------
2935      case( PROJ_MERC )
2936        call map_set(PROJ_MERC, ts_proj,               &
2937                     truelat1 = config_flags%truelat1, &
2938                     lat1     = known_lat,             &
2939                     lon1     = known_lon,             &
2940                     knowni   = 1.,                    &
2941                     knownj   = 1.,                    &
2942                     dx       = config_flags%dx)
2943 !------------------------------------------------------------------------
2944 ! Lambert conformal
2945 !------------------------------------------------------------------------
2946      case( PROJ_LC )
2947        call map_set(PROJ_LC, ts_proj,                  &
2948                     truelat1 = config_flags%truelat1,  &
2949                     truelat2 = config_flags%truelat2,  &
2950                     stdlon   = config_flags%stand_lon, &
2951                     lat1     = known_lat,              &
2952                     lon1     = known_lon,              &
2953                     knowni   = 1.,                     &
2954                     knownj   = 1.,                     &
2955                     dx       = config_flags%dx)
2956 !------------------------------------------------------------------------
2957 ! Polar stereographic
2958 !------------------------------------------------------------------------
2959      case( PROJ_PS )
2960        call map_set(PROJ_PS, ts_proj,                  &
2961                     truelat1 = config_flags%truelat1,  &
2962                     stdlon   = config_flags%stand_lon, &
2963                     lat1     = known_lat,              &
2964                     lon1     = known_lon,              &
2965                     knowni   = 1.,                     &
2966                     knownj   = 1.,                     &
2967                     dx       = config_flags%dx)
2968 !------------------------------------------------------------------------
2969 ! Cassini (global ARW)
2970 !------------------------------------------------------------------------
2971      case( PROJ_CASSINI )
2972        call map_set(PROJ_CASSINI, ts_proj,                            &
2973                     latinc   = grid%dy*360.0/(2.0*EARTH_RADIUS_M*PI), &
2974                     loninc   = grid%dx*360.0/(2.0*EARTH_RADIUS_M*PI), &
2975                     lat1     = known_lat,                             &
2976                     lon1     = known_lon,                             &
2977 ! We still need to get POLE_LAT and POLE_LON metadata variables before
2978 !   this will work for rotated poles.
2979                     lat0     = 90.0,                                  &
2980                     lon0     = 0.0,                                   &
2981                     knowni   = 1.,                                    &
2982                     knownj   = 1.,                                    &
2983                     stdlon   = config_flags%stand_lon)
2984 !------------------------------------------------------------------------
2985 ! Rotated latitude-longitude
2986 !------------------------------------------------------------------------
2987      case( PROJ_ROTLL )
2988        call map_set(PROJ_ROTLL, ts_proj,                      &
2989                     ixdim    = grid%e_we-1,                   &
2990                     jydim    = grid%e_sn-1,                   &
2991                     phi      = real(grid%e_sn-2)*grid%dy/2.0, &
2992                     lambda   = real(grid%e_we-2)*grid%dx,     &
2993                     lat1     = config_flags%cen_lat,          &
2994                     lon1     = config_flags%cen_lon,          &
2995                     latinc   = grid%dy,                       &
2996                     loninc   = grid%dx,                       &
2997                     stagger  = HH)
2998    end select
2999
3000    end subroutine trajmapproj
3001
3002    subroutine UPCASE( lstring )
3003 !----------------------------------------------------------------------
3004 !       ... Convert character string lstring to upper case
3005 !----------------------------------------------------------------------
3006    implicit none
3007
3008 !-----------------------------------------------------------------------
3009 !       ... Dummy args
3010 !-----------------------------------------------------------------------
3011    character(len=*), intent(inout) ::  lstring
3012
3013 !-----------------------------------------------------------------------
3014 !       ... Local variables
3015 !-----------------------------------------------------------------------
3016    integer :: i
3017
3018    do i = 1,LEN_TRIM( lstring )
3019      if( ICHAR(lstring(i:i)) >= 97 .and.  ICHAR(lstring(i:i)) <= 122 ) then
3020        lstring(i:i) = CHAR(ICHAR(lstring(i:i)) - 32)
3021      end if
3022    end do
3023
3024    end subroutine UPCASE
3025
3026    end module module_trajectory