!++mmb ! module module_trajectory module trajectory !----------------------------------------------------------------------------- ! NOTE! This code assumes only one domain at this time. The pkg_has_vars ! array will need to have a domain dimension. !----------------------------------------------------------------------------- implicit none private public :: trajectory_init #ifdef NETCDF public :: trajectory_driver #endif public :: traject integer, parameter :: vals_max = 1000 integer, parameter :: traj_max = 1000 integer, parameter :: var_max = 100 integer, parameter :: pkg_max = 3 real, parameter :: missing_val = -9999. integer :: n_dom integer :: n_vals type default character(len=19) :: start_time character(len=19) :: stop_time character(len=32) :: chm_name(var_max) character(len=32) :: hyd_name(var_max) character(len=32) :: trc_name(var_max) end type default type base real :: lon real :: lat real :: lev real :: x, y real :: traj_var(var_max) integer :: n_chm_var integer :: n_ct_var integer :: n_hyd_var integer :: n_trc_var integer :: n_dyn_var integer :: chm_ndx(var_max) integer :: hyd_ndx(var_max) integer :: trc_ndx(var_max) character(len=19) :: start_time character(len=19) :: stop_time character(len=32) :: chm_spc(var_max) character(len=32) :: hyd_spc(var_max) character(len=32) :: trc_spc(var_max) logical :: in_dom logical :: in_patch logical :: chm_is_gas(var_max) end type base type buffer real :: trj_i(vals_max) real :: trj_j(vals_max) real :: trj_k(vals_max) real :: trj_lons(vals_max) real :: trj_lats(vals_max) real, allocatable :: chm_vals(:,:) real, allocatable :: trc_vals(:,:) real, allocatable :: hyd_vals(:,:) character(len=19) :: times(vals_max) end type buffer integer, allocatable :: traj_cnt(:) type(base), allocatable, target :: traject(:,:) type(base), pointer :: trjects(:) type(buffer), allocatable, target :: trj_buff(:,:) type(buffer), pointer :: trj_pbf(:) character(len=3) :: pkg_tag(pkg_max) = (/ 'chm', 'hyd', 'trc' /) logical, allocatable, target :: pkg_has_vars_dm(:,:,:) CONTAINS subroutine trajectory_init( grid, config_flags, & ims,ime, jms,jme, kms,kme ) use module_domain use module_llxy, only : proj_info, latlon_to_ij use module_configure, only : grid_config_rec_type use module_state_description, only : no_trajectory, param_first_scalar, num_chem, num_moist, num_tracer use module_scalar_tables, only : chem_dname_table, moist_dname_table, tracer_dname_table use module_model_constants, only : g #ifdef DM_PARALLEL use module_dm, only : wrf_dm_sum_integer, wrf_dm_max_real #endif !----------------------------------------------------------------------------- ! dummy arguments !----------------------------------------------------------------------------- integer, intent(in) :: ims,ime, jms,jme, kms,kme type(domain), intent(inout) :: grid type(grid_config_rec_type), intent(in) :: config_flags !----------------------------------------------------------------------------- ! local variables !----------------------------------------------------------------------------- integer :: astat integer :: dm integer :: ierr, ios integer :: i, j, k, k1, m, m1 integer :: i_end, j_end integer :: n, pkg, trj integer :: n_traj, n_traj_1 integer :: p_size integer :: unitno integer :: ids,ide, jds,jde, kds,kde integer :: ips,ipe, jps,jpe, kps,kpe integer :: n_def_var(pkg_max) real :: x, y real :: z_dm_bot, z_dm_top real :: z(kms:kme-1) real :: z_at_w(ims:ime,kms:kme,jms:jme) character(len=32) :: var_name character(len=128) :: filename character(len=256) :: err_mes character(len=32) :: wrk_var_name character(len=32) :: wrk_chr(var_max) character(len=32) :: wrk_def_name(var_max) logical :: exists logical :: is_root_proc logical :: rstrt logical :: mask(var_max) logical, pointer :: pkg_has_vars(:,:) type(grid_config_rec_type) :: config_temp type(proj_info) :: proj type(default) :: traj_def type(base) :: traj_type(traj_max) logical, external :: wrf_dm_on_monitor integer, external :: get_unused_unit namelist / traj_spec / traj_type namelist / traj_default / traj_def #ifndef NETCDF call wrf_error_fatal( 'trajectory_init: requires netcdf' ) #endif !----------------------------------------------------------------------------- ! check for trajectory option !----------------------------------------------------------------------------- if( grid%traj_opt == no_trajectory ) then write(err_mes,'(''trajectory_init: traj_opt = no_trajectory for grid '',i2.2)') grid%id call wrf_message( trim(err_mes) ) return endif !----------------------------------------------------------------------------- ! set domain count, restarting? !----------------------------------------------------------------------------- dm = grid%id call nl_get_restart( dm,rstrt ) if( dm == 1 ) then call nl_get_max_dom( 1,n_dom ) !----------------------------------------------------------------------------- ! allocate module variables !----------------------------------------------------------------------------- allocate( traject(traj_max,n_dom),traj_cnt(n_dom), & pkg_has_vars_dm(traj_max,pkg_max,n_dom),stat=astat ) if( astat /= 0 ) then write(err_mes,'(''trajectory_init: failed to allocate traject,traj_cnt; error = '',i6)') astat call wrf_error_fatal( trim( err_mes ) ) endif endif trjects => traject(:,dm) call get_ijk_from_grid( grid , & ids, ide, jds, jde, kds, kde, & n,n, n,n, n,n, & ips, ipe, jps, jpe, kps, kpe ) n_vals = 0 traj_cnt(:) = 0 is_root_proc = wrf_dm_on_monitor() pkg_has_vars => pkg_has_vars_dm(:,:,dm) !----------------------------------------------------------------------------- ! master proc !----------------------------------------------------------------------------- master_proc: & if( is_root_proc ) then write(filename,'(''wrfinput_traj_d'',i2.2)',iostat=ios) grid%id if( ios /= 0 ) then write(err_mes,'(''trajectory_init: failed to set filename: error = '',i6)') ios call wrf_error_fatal( trim( err_mes ) ) endif inquire( file=trim(filename),exist=exists ) input_file: & if( exists ) then unitno = get_unused_unit() if( unitno <= 0 ) then call wrf_error_fatal( 'trajectory_init: failed to get unit number' ) endif !----------------------------------------------------------------------------- ! open file !----------------------------------------------------------------------------- open( unit = unitno,file=trim(filename),iostat=ios ) if( ios /= 0 ) then write(err_mes,'(''trajectory_init: failed to open '',a,''; error = '',i6)') trim(filename),ios call wrf_error_fatal( trim( err_mes ) ) endif !----------------------------------------------------------------------------- ! initialize trajectories !----------------------------------------------------------------------------- traj_def%start_time = ' ' traj_def%stop_time = ' ' traj_def%chm_name(:) = ' ' traj_def%hyd_name(:) = ' ' traj_def%trc_name(:) = ' ' do trj = 1,traj_max traj_type(trj)%start_time = ' ' traj_type(trj)%stop_time = ' ' traj_type(trj)%chm_spc(:) = ' ' traj_type(trj)%hyd_spc(:) = ' ' traj_type(trj)%trc_spc(:) = ' ' traj_type(trj)%chm_ndx(:) = 0 traj_type(trj)%hyd_ndx(:) = 0 traj_type(trj)%trc_ndx(:) = 0 traj_type(trj)%n_chm_var = 0 ; traj_type(trj)%n_ct_var = 0 traj_type(trj)%n_hyd_var = 0 ; traj_type(trj)%n_trc_var = 0 ; traj_type(trj)%n_dyn_var = 0 traj_type(trj)%chm_is_gas(:) = .true. traj_type(trj)%lon = missing_val traj_type(trj)%lat = missing_val traj_type(trj)%lev = missing_val end do !----------------------------------------------------------------------------- ! read file !----------------------------------------------------------------------------- read(unit=unitno,nml=traj_default,iostat=ios) if( ios /= 0 ) then close( unit=unitno ) write(err_mes,'(''trajectory_init: failed to read '',a,''; error = '',i6)') trim(filename),ios call wrf_error_fatal( trim( err_mes ) ) endif read(unit=unitno,nml=traj_spec,iostat=ios) if( ios /= 0 ) then close( unit=unitno ) write(err_mes,'(''trajectory_init: failed to read '',a,''; error = '',i6)') trim(filename),ios call wrf_error_fatal( trim( err_mes ) ) endif close( unit=unitno ) else input_file write(err_mes,'(''trajectory_init: failed to find '',a)') trim(filename) call wrf_error_fatal( trim( err_mes ) ) endif input_file !----------------------------------------------------------------------------- ! process the namelist input !----------------------------------------------------------------------------- do pkg = 1,pkg_max select case( trim(pkg_tag(pkg)) ) case( 'chm' ) wrk_def_name(:) = traj_def%chm_name(:) case( 'hyd' ) wrk_def_name(:) = traj_def%hyd_name(:) case( 'trc' ) wrk_def_name(:) = traj_def%trc_name(:) end select do m = 1,var_max if( wrk_def_name(m) == ' ' ) then exit endif end do n_def_var(pkg) = m - 1 end do if( traj_def%start_time /= ' ' ) then write(*,*) 'trajectory_init: default start time = ',traj_def%start_time endif if( traj_def%stop_time /= ' ' ) then write(*,*) 'trajectory_init: default stop time = ',traj_def%stop_time endif do pkg = 1,pkg_max if( n_def_var(pkg) > 0 ) then write(*,*) ' ' write(*,*) 'trajectory_init: default ',pkg_tag(pkg),' variables' select case( trim(pkg_tag(pkg)) ) case( 'chm' ) wrk_def_name(:) = traj_def%chm_name(:) case( 'hyd' ) wrk_def_name(:) = traj_def%hyd_name(:) case( 'trc' ) wrk_def_name(:) = traj_def%trc_name(:) end select write(*,*) wrk_def_name(:n_def_var(pkg)) endif end do do n_traj = 1,traj_max if( traj_type(n_traj)%lon == missing_val .or. & traj_type(n_traj)%lat == missing_val .or. & traj_type(n_traj)%lev == missing_val ) then exit endif end do n_traj = n_traj - 1 has_trajectories: & if( n_traj > 0 ) then !----------------------------------------------------------------------------- ! set individual trajectories to default if specified !----------------------------------------------------------------------------- if( traj_def%start_time /= ' ' ) then do trj = 1,n_traj if( traj_type(trj)%start_time == ' ' ) then traj_type(trj)%start_time = traj_def%start_time endif end do endif if( traj_def%stop_time /= ' ' ) then do trj = 1,n_traj if( traj_type(trj)%stop_time == ' ' ) then traj_type(trj)%stop_time = traj_def%stop_time endif end do endif do pkg = 1,pkg_max select case( trim(pkg_tag(pkg)) ) case( 'chm' ) wrk_def_name(:) = traj_def%chm_name(:) case( 'hyd' ) wrk_def_name(:) = traj_def%hyd_name(:) case( 'trc' ) wrk_def_name(:) = traj_def%trc_name(:) end select do trj = 1,n_traj select case( trim(pkg_tag(pkg)) ) case( 'chm' ) wrk_var_name = traj_type(trj)%chm_spc(1) case( 'hyd' ) wrk_var_name = traj_type(trj)%hyd_spc(1) case( 'trc' ) wrk_var_name = traj_type(trj)%trc_spc(1) end select if( wrk_var_name == ' ' ) then m1 = n_def_var(pkg) select case( trim(pkg_tag(pkg)) ) case( 'chm' ) traj_type(trj)%chm_spc(:m1) = traj_def%chm_name(:m1) case( 'hyd' ) traj_type(trj)%hyd_spc(:m1) = traj_def%hyd_name(:m1) case( 'trc' ) traj_type(trj)%trc_spc(:m1) = traj_def%trc_name(:m1) end select endif end do end do !----------------------------------------------------------------------------- ! get variable counts !----------------------------------------------------------------------------- do trj = 1,n_traj if( num_chem > 1 ) then call get_var_cnt( traj_type(trj)%n_chm_var, traj_type(trj)%chm_spc ) else traj_type(trj)%n_chm_var = 0 endif if( num_moist > 1 ) then call get_var_cnt( traj_type(trj)%n_hyd_var, traj_type(trj)%hyd_spc ) else traj_type(trj)%n_hyd_var = 0 endif if( num_tracer > 1 ) then call get_var_cnt( traj_type(trj)%n_trc_var, traj_type(trj)%trc_spc ) else traj_type(trj)%n_trc_var = 0 endif end do !----------------------------------------------------------------------------- ! check for trajectory variables in simulation !----------------------------------------------------------------------------- do trj = 1,n_traj if( traj_type(trj)%n_chm_var > 0 .and. num_chem > 1 ) then mask(:) = .false. call scan_vars( traj_type(trj)%n_chm_var, traj_type(trj)%chm_spc, traj_type(trj)%chm_ndx, & num_chem, chem_dname_table(dm,:), & traj_type(trj)%chm_is_gas, 'chm' ) endif if( traj_type(trj)%n_hyd_var > 0 .and. num_moist > 1 ) then mask(:) = .false. call scan_vars( traj_type(trj)%n_hyd_var, traj_type(trj)%hyd_spc, traj_type(trj)%hyd_ndx, & num_moist, moist_dname_table(dm,:), & traj_type(trj)%chm_is_gas, 'hyd' ) endif if( traj_type(trj)%n_trc_var > 0 .and. num_tracer > 1 ) then mask(:) = .false. call scan_vars( traj_type(trj)%n_trc_var, traj_type(trj)%trc_spc, traj_type(trj)%trc_ndx, & num_tracer, tracer_dname_table(dm,:), & traj_type(trj)%chm_is_gas, 'trc' ) endif end do n_traj = count( (traj_type(:n_traj)%n_chm_var+traj_type(:n_traj)%n_hyd_var+traj_type(:n_traj)%n_trc_var) > 0 ) !----------------------------------------------------------------------------- ! allocate buffer type !----------------------------------------------------------------------------- if( dm == 1 ) then allocate( trj_buff(traj_max,n_dom),stat=astat ) if( astat /= 0 ) then write(err_mes,'(''trajectory_init: failed to allocate traj_buff; error = '',i6)') astat call wrf_error_fatal( trim( err_mes ) ) endif endif trj_pbf => trj_buff(:,dm) endif has_trajectories endif master_proc !----------------------------------------------------------------------------- ! if initial run, overwrite existing grid trajectory variables !----------------------------------------------------------------------------- #ifdef DM_PARALLEL call wrf_dm_bcast_integer( n_traj,1 ) if( n_traj > 0 ) then p_size = (5+var_max)*RWORDSIZE + (2+var_max)*LWORDSIZE & + (5+3*var_max)*IWORDSIZE + 3*32*var_max + 38 do trj = 1,n_traj call wrf_dm_bcast_bytes( traj_type(trj),p_size ) end do endif #endif is_cold_start: & if( .not. rstrt ) then !----------------------------------------------------------------------------- ! calc and check trajectory start point !----------------------------------------------------------------------------- has_trajectories_a: & if( n_traj > 0 ) then config_temp = config_flags call trajmapproj( grid, config_temp, proj ) i_end = min(ipe,ide-1) j_end = min(jpe,jde-1) do j = jps,j_end do k = kps,kpe z_at_w(ips:i_end,k,j) = (grid%ph_2(ips:i_end,k,j) + grid%phb(ips:i_end,k,j))/g end do end do !----------------------------------------------------------------------------- ! first check x,y !----------------------------------------------------------------------------- traj_loop: & do trj = 1,n_traj if( traj_type(trj)%lat /= missing_val .and. & traj_type(trj)%lon /= missing_val ) then call latlon_to_ij( proj, traj_type(trj)%lat, traj_type(trj)%lon, x, y ) traj_type(trj)%in_dom = & (x >= real(ids) .and. x <= real(ide-1) .and. & y >= real(jds) .and. y <= real(jde-1)) write(err_mes,'(''trajectory_init: x,y = '',1p,2g15.7)') x,y call wrf_debug( 0,trim(err_mes) ) else traj_type(trj)%in_dom = .false. endif write(err_mes,'(''trajectory_init: traj '',i1,'' in dom = '',l)') trj,traj_type(trj)%in_dom call wrf_debug( 0,trim(err_mes) ) !----------------------------------------------------------------------------- ! then check z !----------------------------------------------------------------------------- is_in_domain: & if( traj_type(trj)%in_dom ) then i = nint( x ) j = nint( y ) traj_type(trj)%in_patch = & (i >= ips .and. i <= min(ipe,ide-1) .and. & j >= jps .and. j <= min(jpe,jde-1)) is_in_patch: & if( traj_type(trj)%in_patch ) then k1 = kde - 1 ! z_dm_bot = .5*(z_at_w(i,kds,j) + z_at_w(i,kds+1,j)) ! z_dm_top = .5*(z_at_w(i,k1,j) + z_at_w(i,k1+1,j)) z_dm_bot = z_at_w(i,kds,j) z_dm_top = z_at_w(i,k1,j) write(err_mes,'(''trajectory_init: traj '',i1,'' i,j,z_bot,z_top,lev = '',2i4,1p3g15.7)') & trj,i,j,z_dm_bot,z_dm_top,traj_type(trj)%lev call wrf_debug( 0,trim(err_mes) ) if( traj_type(trj)%lev >= z_dm_bot .and. & traj_type(trj)%lev <= z_dm_top ) then traj_type(trj)%in_dom = .true. traj_type(trj)%x = x traj_type(trj)%y = y else traj_type(trj)%in_dom = .false. endif else is_in_patch traj_type(trj)%x = missing_val traj_type(trj)%y = missing_val endif is_in_patch endif is_in_domain end do traj_loop ! n_traj = count( traj_type(:n_traj)%in_dom ) n_traj_1 = count( traj_type(:n_traj)%in_dom ) write(err_mes,'(''trajectory_init: traj cnt = '',2i6)') n_traj_1,n_traj call wrf_debug( 0,trim(err_mes) ) endif has_trajectories_a has_trajectories_b: & if( n_traj_1 > 0 ) then grid%traj_i(:) = missing_val grid%traj_j(:) = missing_val grid%traj_k(:) = missing_val grid%traj_long(:) = missing_val grid%traj_lat(:) = missing_val !----------------------------------------------------------------------------- ! set initial trajectory spatial coordinates !----------------------------------------------------------------------------- k1 = kde - 1 do trj = 1,n_traj if( traj_type(trj)%in_patch ) then grid%traj_i(trj) = traj_type(trj)%x grid%traj_j(trj) = traj_type(trj)%y grid%traj_long(trj) = traj_type(trj)%lon grid%traj_lat(trj) = traj_type(trj)%lat i = nint( traj_type(trj)%x ) j = nint( traj_type(trj)%y ) ! z(kds:k1) = .5*(z_at_w(i,kds:k1,j) + z_at_w(i,kds+1:k1+1,j)) z(kds:k1) = z_at_w(i,kds:k1,j) do k = kds+1,k1 if( traj_type(trj)%lev <= z(k) ) then grid%traj_k(trj) = real(k - 1) & + (traj_type(trj)%lev - z(k-1))/(z(k) - z(k-1)) exit endif end do write(err_mes,'(''trajectory_init: trj,k,z(k-1:k),traj_k = '',2i3,1p3g15.7)') & trj,k,z(k-1:k),grid%traj_k(trj) call wrf_debug( 0,trim(err_mes) ) endif end do else has_trajectories_b grid%traj_opt = no_trajectory return endif has_trajectories_b endif is_cold_start !----------------------------------------------------------------------------- ! transfer from local variable to module variable !----------------------------------------------------------------------------- traj_cnt(dm) = n_traj do trj = 1,n_traj traject(trj,dm) = traj_type(trj) end do !----------------------------------------------------------------------------- ! create netcdf trajectory file !----------------------------------------------------------------------------- if( .not. rstrt ) then do trj = 1,n_traj #ifdef DM_PARALLEL grid%traj_i(trj) = wrf_dm_max_real( grid%traj_i(trj) ) grid%traj_j(trj) = wrf_dm_max_real( grid%traj_j(trj) ) grid%traj_k(trj) = wrf_dm_max_real( grid%traj_k(trj) ) grid%traj_long(trj) = wrf_dm_max_real( grid%traj_long(trj) ) grid%traj_lat(trj) = wrf_dm_max_real( grid%traj_lat(trj) ) #else grid%traj_i(trj) = grid%traj_i(trj) grid%traj_j(trj) = grid%traj_j(trj) grid%traj_k(trj) = grid%traj_k(trj) grid%traj_long(trj) = grid%traj_long(trj) grid%traj_lat(trj) = grid%traj_lat(trj) #endif end do call trajectory_create_file( grid ) endif do pkg = 1,pkg_max select case( trim(pkg_tag(pkg)) ) case( 'chm' ) pkg_has_vars(:n_traj,pkg) = traject(:n_traj,dm)%n_chm_var > 0 case( 'hyd' ) pkg_has_vars(:n_traj,pkg) = traject(:n_traj,dm)%n_hyd_var > 0 case( 'trc' ) pkg_has_vars(:n_traj,pkg) = traject(:n_traj,dm)%n_trc_var > 0 end select end do if( is_root_proc ) then !----------------------------------------------------------------------------- ! allocate buffer arrays !----------------------------------------------------------------------------- do trj = 1,n_traj if( traject(trj,dm)%in_dom ) then do pkg = 1,pkg_max select case( trim(pkg_tag(pkg)) ) case( 'chm' ) m1 = traject(trj,dm)%n_chm_var if( m1 > 0 ) then allocate( trj_pbf(trj)%chm_vals(vals_max,m1),stat=astat) endif case( 'hyd' ) m1 = traject(trj,dm)%n_hyd_var if( m1 > 0 ) then allocate( trj_pbf(trj)%hyd_vals(vals_max,m1),stat=astat) endif case( 'trc' ) m1 = traject(trj,dm)%n_trc_var if( m1 > 0 ) then allocate( trj_pbf(trj)%trc_vals(vals_max,m1),stat=astat) endif end select if( astat /= 0 ) then write(err_mes,'(''trajectory_init: failed to allocate buffer%'',a,''; error = '',i6)') & pkg_tag(pkg),astat call wrf_error_fatal( trim( err_mes ) ) endif end do endif end do endif CONTAINS subroutine get_var_cnt( n_vars, vars ) !----------------------------------------------------------------------------- ! dummy arguments !----------------------------------------------------------------------------- integer, intent(inout) :: n_vars character(len=32), intent(inout) :: vars(:) mask(:) = vars(:) /= ' ' wrk_chr(:) = ' ' m1 = 0 do n = 1,var_max if( mask(n) ) then m1 = m1 + 1 wrk_chr(m1) = vars(n) endif end do n_vars = count( wrk_chr(:) /= ' ' ) vars(:) = wrk_chr(:) end subroutine get_var_cnt subroutine scan_vars( n_vars, vars, var_ndx, n_tbl, tbl, & is_gas, var_type ) !----------------------------------------------------------------------------- ! dummy arguments !----------------------------------------------------------------------------- integer, intent(inout) :: n_vars integer, intent(in) :: n_tbl integer, intent(inout) :: var_ndx(:) logical, intent(inout) :: is_gas(:) character(len=*), intent(in) :: var_type character(len=32), intent(inout) :: vars(:) character(len=256), intent(in) :: tbl(:) !----------------------------------------------------------------------------- ! local variables !----------------------------------------------------------------------------- integer :: ndx(var_max) var_loop: & do n = 1,n_vars var_name = vars(n) if( trim(var_type) == 'chm' ) then i = index( '(a)', var_name ) if( i == 0 ) then i = index( '(A)', var_name ) endif if( i > 0 ) then is_gas(n) = .false. var_name(i:) = ' ' vars(n)(i:) = ' ' endif endif do m1 = param_first_scalar,n_tbl if( trim(var_name) == trim(tbl(m1)) ) then mask(n) = .true. ndx(n) = m1 exit endif end do if( .not. mask(n) ) then write(err_mes,'(''trajectory_init: '',a,'' not in '',a,'' opt'')') trim(var_name),var_type call wrf_message( trim(err_mes) ) endif end do var_loop wrk_chr(:) = ' ' m1 = 0 do n = 1,var_max if( mask(n) ) then m1 = m1 + 1 wrk_chr(m1) = vars(n) endif end do var_ndx(:count(mask)) = pack( ndx(:),mask=mask) vars(:) = wrk_chr(:) n_vars = count( mask ) end subroutine scan_vars end subroutine trajectory_init #ifdef NETCDF subroutine trajectory_driver( grid ) #ifdef DM_PARALLEL use module_dm, only : & local_communicator, mytask, ntasks, ntasks_x, ntasks_y & ,local_communicator_periodic, wrf_dm_max_real, wrf_dm_max_int use module_comm_dm, only : halo_em_chem_e_3_sub, halo_em_moist_e_3_sub use module_comm_dm, only : halo_em_tracer_e_3_sub #endif use module_domain use module_date_time use module_state_description, only : num_chem, num_moist, num_tracer !----------------------------------------------------------------------------- ! dummy arguments !----------------------------------------------------------------------------- type(domain), intent(in) :: grid !----------------------------------------------------------------------------- ! local variables !----------------------------------------------------------------------------- integer :: ims,ime, jms,jme, kms,kme integer :: ids,ide, jds,jde, kds,kde integer :: ips,ipe, jps,jpe, kps,kpe integer :: dm integer :: j, k integer :: il, iu, ios, jl, ju, kl, ku integer :: m, mu, n, ndx, n_vars, n_traj integer :: pkg_var_cnt(traj_max) integer :: ncid integer :: pkg, trj #ifndef DM_PARALLEL integer :: mytask #endif integer :: traj_proc(traj_max), glb_traj_proc(traj_max) real :: x, y, zi real :: delsx, delsy, o_delsx, o_delsy real :: delsz, o_delsz real :: max_conc real :: horz_conc(2) real, pointer :: traj_conc(:) real, target :: traj_val(var_max,traj_max) real, pointer :: wrk4d(:,:,:,:) real, allocatable, target :: chem(:,:,:,:) real, allocatable, target :: moist(:,:,:,:) real, allocatable, target :: tracer(:,:,:,:) character(len=19) :: current_timestr, next_timestr character(len=32) :: var_name(var_max) character(len=256) :: err_mes logical :: is_root_proc logical :: is_in_patch_gap logical :: flsh_buff logical :: traj_is_active(traj_max) logical :: pkg_is_active(traj_max,pkg_max) logical, pointer :: pkg_has_vars(:,:) logical, external :: wrf_dm_on_monitor type(WRFU_Time) :: current_time, next_time type(WRFU_Time) :: start_time, stop_time include 'netcdf.inc' #ifdef DM_PARALLEL include 'mpif.h' #endif #ifndef DM_PARALLEL mytask = 0 #endif dm = grid%id n_traj = traj_cnt(dm) has_trajectories: & if( n_traj > 0 ) then trjects => traject(:,dm) is_root_proc = wrf_dm_on_monitor() if( is_root_proc ) then trj_pbf => trj_buff(:,dm) endif call get_ijk_from_grid( grid , & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe ) !----------------------------------------------------------------------------- ! is trajectory in time interval? !----------------------------------------------------------------------------- call domain_clock_get( grid, current_time=current_time, current_timestr=current_timestr) call geth_newdate( next_timestr, current_timestr, int(grid%dt) ) call wrf_atotime( next_timestr, next_time ) do trj = 1,n_traj call wrf_atotime( traject(trj,dm)%start_time, start_time ) call wrf_atotime( traject(trj,dm)%stop_time, stop_time ) traj_is_active(trj) = next_time .ge. start_time .and. next_time .le. stop_time end do !----------------------------------------------------------------------------- ! is trajectory in domain? !----------------------------------------------------------------------------- pkg_has_vars => pkg_has_vars_dm(:,:,dm) do trj = 1,n_traj if( traj_is_active(trj) ) then traject(trj,dm)%in_dom = grid%traj_i(trj) >= real(ids) .and. grid%traj_i(trj) <= real(ide-1) if( traject(trj,dm)%in_dom ) then traject(trj,dm)%in_dom = grid%traj_j(trj) >= real(jds) .and. grid%traj_j(trj) <= real(jde-1) endif if( traject(trj,dm)%in_dom ) then traject(trj,dm)%in_dom = grid%traj_k(trj) >= real(kps) .and. grid%traj_k(trj) <= real( min( kpe,kde-1 ) ) endif traj_is_active(trj) = traject(trj,dm)%in_dom endif end do do pkg = 1,pkg_max pkg_is_active(:n_traj,pkg) = traj_is_active(:n_traj) .and. pkg_has_vars(:n_traj,pkg) end do !----------------------------------------------------------------------------- ! is trajectory in mpi process? !----------------------------------------------------------------------------- traj_proc(:n_traj) = -1 do trj = 1,n_traj if( traj_is_active(trj) ) then traject(trj,dm)%in_patch = grid%traj_i(trj) >= real(ips) .and. grid%traj_i(trj) <= real( min( ipe,ide-1 ) ) if( traject(trj,dm)%in_patch ) then traject(trj,dm)%in_patch = grid%traj_j(trj) >= real(jps) .and. grid%traj_j(trj) <= real( min( jpe,jde-1 ) ) endif if( traject(trj,dm)%in_patch ) then traject(trj,dm)%in_patch = grid%traj_k(trj) >= real(kps) .and. grid%traj_k(trj) <= real( min( kpe,kde-1 ) ) endif if( traject(trj,dm)%in_patch ) then traj_proc(trj) = mytask + 1 else traj_proc(trj) = 0 endif endif end do do trj = 1,n_traj glb_traj_proc(trj) = wrf_dm_max_int( traj_proc(trj) ) end do !----------------------------------------------------------------------------- ! any trajectories in "gap" between patches? !----------------------------------------------------------------------------- do trj = 1,n_traj if( traj_is_active(trj) .and. glb_traj_proc(trj) == 0 ) then traject(trj,dm)%in_patch = grid%traj_i(trj) >= real(ips) .and. grid%traj_i(trj) <= real( min( ipe+1,ide-1 ) ) if( traject(trj,dm)%in_patch ) then traject(trj,dm)%in_patch = grid%traj_j(trj) >= real(jps) .and. grid%traj_j(trj) <= real( min( jpe+1,jde-1 ) ) endif if( traject(trj,dm)%in_patch ) then traject(trj,dm)%in_patch = grid%traj_k(trj) >= real(kps) .and. grid%traj_k(trj) <= real( min( kme,kde-1 ) ) endif if( traject(trj,dm)%in_patch ) then traj_proc(trj) = mytask + 1 else traj_proc(trj) = 0 endif if( traj_proc(trj) /= 0 ) then call wrf_debug( 0,'^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^') write(err_mes,'(''Gapper '',i5,''; x,y,zi = '',1p,3g15.7)') trj,grid%traj_i(trj),grid%traj_j(trj),grid%traj_k(trj) call wrf_debug( 0,trim(err_mes) ) write(err_mes,'(''Gapper ips,ipe,jps,jpe = '',4i5)') ips,ipe,jps,jpe call wrf_debug( 0,trim(err_mes) ) call wrf_debug( 0,'^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^') endif endif end do !----------------------------------------------------------------------------- ! buffer traj time and position !----------------------------------------------------------------------------- if( is_root_proc ) then n_vals = n_vals + 1 if( grid%itimestep > 0 ) then trj_pbf(:n_traj)%times(n_vals) = next_timestr else trj_pbf(:n_traj)%times(n_vals) = current_timestr endif do trj = 1,n_traj trj_pbf(trj)%trj_i(n_vals) = grid%traj_i(trj) trj_pbf(trj)%trj_j(n_vals) = grid%traj_j(trj) trj_pbf(trj)%trj_k(n_vals) = grid%traj_k(trj) trj_pbf(trj)%trj_lons(n_vals) = grid%traj_long(trj) trj_pbf(trj)%trj_lats(n_vals) = grid%traj_lat(trj) end do endif do trj = 1,n_traj traj_val(:,trj) = missing_val end do pkg_loop: & do pkg = 1,pkg_max pkg_has_active_traj: & if( any( pkg_is_active(:n_traj,pkg) ) ) then !----------------------------------------------------------------------------- ! allocate working data array !----------------------------------------------------------------------------- select case( trim(pkg_tag(pkg)) ) case( 'chm' ) allocate( chem(ims:ime,kms:kme,jms:jme,num_chem),stat=ios ) case( 'hyd' ) allocate( moist(ims:ime,kms:kme,jms:jme,num_moist),stat=ios ) case( 'trc' ) allocate( tracer(ims:ime,kms:kme,jms:jme,num_tracer),stat=ios ) end select if( ios /= 0 ) then write(err_mes,'(''trajectory_driver: failed to allocate wrk4d: error = '',i6)') ios call wrf_error_fatal( trim( err_mes ) ) endif select case( trim(pkg_tag(pkg)) ) case( 'chm' ) do m = 1,num_chem do j = jps,jpe do k = kps,kpe chem(ips:ipe,k,j,m) = grid%chem(ips:ipe,k,j,m) end do end do end do case( 'hyd' ) do m = 1,num_moist do j = jps,jpe do k = kps,kpe moist(ips:ipe,k,j,m) = grid%moist(ips:ipe,k,j,m) end do end do end do case( 'trc' ) do m = 1,num_tracer do j = jps,jpe do k = kps,kpe tracer(ips:ipe,k,j,m) = grid%tracer(ips:ipe,k,j,m) end do end do end do end select #ifdef DM_PARALLEL !----------------------------------------------------------------------------- ! any trajectories in "gap" between patches? !----------------------------------------------------------------------------- is_in_patch_gap = any( glb_traj_proc(:n_traj) == 0 .and. pkg_is_active(:n_traj,pkg) ) if( is_in_patch_gap ) then ! write(err_mes,'(''glb_traj_proc mask = '',10i5)') glb_traj_proc(:n_traj) ! call wrf_debug( 0,trim(err_mes) ) select case( trim(pkg_tag(pkg)) ) case( 'chm' ) # include "HALO_EM_CHEM_E_3.inc" case( 'hyd' ) # include "HALO_EM_MOIST_E_3.inc" case( 'trc' ) # include "HALO_EM_TRACER_E_3.inc" end select endif #endif traj_loop: & do trj = 1,n_traj select case( trim(pkg_tag(pkg)) ) case( 'chm' ) n_vars = traject(trj,dm)%n_chm_var case( 'hyd' ) n_vars = traject(trj,dm)%n_hyd_var case( 'trc' ) n_vars = traject(trj,dm)%n_trc_var end select pkg_is_active_in_traj: & if( pkg_is_active(trj,pkg) ) then select case( trim(pkg_tag(pkg)) ) case( 'chm' ) wrk4d => chem case( 'hyd' ) wrk4d => moist case( 'trc' ) wrk4d => tracer end select in_patch: if( traj_proc(trj) == mytask+1 ) then x = grid%traj_i(trj) y = grid%traj_j(trj) zi = grid%traj_k(trj) il = int( x ) ; iu = il + 1 jl = int( y ) ; ju = jl + 1 kl = int( zi ) ; ku = kl + 1 delsx = x - floor(x) ; o_delsx = 1. - delsx delsy = y - floor(y) ; o_delsy = 1. - delsy delsz = zi - floor(zi) ; o_delsz = 1. - delsz var_loop: do n = 1,n_vars select case( trim(pkg_tag(pkg)) ) case( 'chm' ) ndx = traject(trj,dm)%chm_ndx(n) case( 'hyd' ) ndx = traject(trj,dm)%hyd_ndx(n) case( 'trc' ) ndx = traject(trj,dm)%trc_ndx(n) end select horz_conc(1) = o_delsx*o_delsy*wrk4d(il,kl,jl,ndx) + o_delsy*delsx*wrk4d(iu,kl,jl,ndx) & + delsy*o_delsx*wrk4d(il,kl,ju,ndx) + delsx*delsy*wrk4d(iu,kl,ju,ndx) horz_conc(2) = o_delsx*o_delsy*wrk4d(il,ku,jl,ndx) + o_delsy*delsx*wrk4d(iu,ku,jl,ndx) & + delsy*o_delsx*wrk4d(il,ku,ju,ndx) + delsx*delsy*wrk4d(iu,ku,ju,ndx) traject(trj,dm)%traj_var(n) = delsz*horz_conc(2) + o_delsz*horz_conc(1) end do var_loop else in_patch traject(trj,dm)%traj_var(:n_vars) = missing_val endif in_patch traj_conc => traj_val(:,trj) do n = 1,n_vars #ifdef DM_PARALLEL max_conc = wrf_dm_max_real( traject(trj,dm)%traj_var(n) ) #else max_conc = traject(trj,dm)%traj_var(n) #endif if( is_root_proc ) then traj_conc(n) = max_conc endif end do else pkg_is_active_in_traj if( is_root_proc .and. n_vars > 0 ) then traj_conc => traj_val(:,trj) traj_conc(:n_vars) = missing_val endif endif pkg_is_active_in_traj !----------------------------------------------------------------------------- ! buffer traj chm,trc,hyb variables !----------------------------------------------------------------------------- if( is_root_proc .and. n_vars > 0 ) then select case( pkg_tag(pkg) ) case( 'chm' ) trj_pbf(trj)%chm_vals(n_vals,:n_vars) = traj_conc(:n_vars) case( 'hyd' ) trj_pbf(trj)%hyd_vals(n_vals,:n_vars) = traj_conc(:n_vars) case( 'trc' ) trj_pbf(trj)%trc_vals(n_vals,:n_vars) = traj_conc(:n_vars) end select endif end do traj_loop if( allocated( chem ) ) then deallocate( chem ) endif if( allocated( moist ) ) then deallocate( moist ) endif if( allocated( tracer ) ) then deallocate( tracer ) endif else pkg_has_active_traj if( is_root_proc ) then do trj = 1,n_traj select case( trim(pkg_tag(pkg)) ) case( 'chm' ) n_vars = traject(trj,dm)%n_chm_var case( 'hyd' ) n_vars = traject(trj,dm)%n_hyd_var case( 'trc' ) n_vars = traject(trj,dm)%n_trc_var end select if( n_vars > 0 ) then select case( pkg_tag(pkg) ) case( 'chm' ) trj_pbf(trj)%chm_vals(n_vals,:n_vars) = missing_val case( 'hyd' ) trj_pbf(trj)%hyd_vals(n_vals,:n_vars) = missing_val case( 'trc' ) trj_pbf(trj)%trc_vals(n_vals,:n_vars) = missing_val end select endif end do endif endif pkg_has_active_traj end do pkg_loop !----------------------------------------------------------------------------- ! output trajectory buffer !----------------------------------------------------------------------------- if( is_root_proc ) then flsh_buff = n_vals == vals_max .or. domain_last_time_step( grid ) if( flsh_buff ) then call trajectory_write_file( n_traj, n_vals, dm ) endif endif endif has_trajectories end subroutine trajectory_driver subroutine trajectory_create_file( grid ) !----------------------------------------------------------------------------- ! create trajectory netcdf file !----------------------------------------------------------------------------- use module_domain !----------------------------------------------------------------------------- ! dummy arguments !----------------------------------------------------------------------------- type(domain), intent(in) :: grid !----------------------------------------------------------------------------- ! local variables !----------------------------------------------------------------------------- integer :: dm integer :: ncid, ios integer :: time_dim, Times_dim integer :: varid integer :: var_dims(2) integer :: m, n, trj character(len=256) :: filename character(len=256) :: var_name character(len=256) :: err_mes character(len=256) :: description character(len=256) :: units logical, external :: wrf_dm_on_monitor include 'netcdf.inc' master_proc: & if( wrf_dm_on_monitor() ) then dm = grid%id write(filename,'(''wrfout_traj_d'',i2.2)',iostat=ios) dm if( ios /= 0 ) then write(err_mes,'(''trajectory_create_file: failed to set filename: error = '',i6)') ios call wrf_error_fatal( trim( err_mes ) ) endif ios = nf_create( trim(filename), or(nf_clobber,nf_netcdf4), ncid ) if( ios /= 0 ) then write(err_mes,'(''trajectory_create_file: failed to create '',a,'': error = '',i6)') trim(filename),ios call wrf_error_fatal( trim( err_mes ) ) endif !----------------------------------------------------------------------------- ! define dimensions !----------------------------------------------------------------------------- err_mes = 'trajectory_create_file: failed to create time dimension' call handle_ncerr( nf_def_dim( ncid, 'time', nf_unlimited, time_dim ), trim(err_mes) ) err_mes = 'trajectory_create_file: failed to create Times dimension' call handle_ncerr( nf_def_dim( ncid, 'DateStrLen', 19, Times_dim ), trim(err_mes) ) !----------------------------------------------------------------------------- ! define variables !----------------------------------------------------------------------------- var_dims(:) = (/ Times_dim,time_dim /) err_mes = 'trajectory_create_file: failed to create Times variable' call handle_ncerr( nf_def_var( ncid, 'Times', nf_char, 2, var_dims, varid ), trim(err_mes) ) var_dims(1) = time_dim traj_loop: & do trj = 1,traj_cnt(dm) write(var_name,'(''traj_i_'',i4.4)') trj err_mes = 'trajectory_create_file: failed to create ' // trim(var_name) // ' variable' call handle_ncerr( nf_def_var( ncid, trim(var_name), nf_real, 1, var_dims, varid ), trim(err_mes) ) write(var_name,'(''traj_j_'',i4.4)') trj err_mes = 'trajectory_create_file: failed to create ' // trim(var_name) // ' variable' call handle_ncerr( nf_def_var( ncid, trim(var_name), nf_real, 1, var_dims, varid ), trim(err_mes) ) write(var_name,'(''traj_k_'',i4.4)') trj err_mes = 'trajectory_create_file: failed to create ' // trim(var_name) // ' variable' call handle_ncerr( nf_def_var( ncid, trim(var_name), nf_real, 1, var_dims, varid ), trim(err_mes) ) write(var_name,'(''traj_long_'',i4.4)') trj err_mes = 'trajectory_create_file: failed to create ' // trim(var_name) // ' variable' call handle_ncerr( nf_def_var( ncid, trim(var_name), nf_real, 1, var_dims, varid ), trim(err_mes) ) write(var_name,'(''traj_lat_'',i4.4)') trj err_mes = 'trajectory_create_file: failed to create ' // trim(var_name) // ' variable' call handle_ncerr( nf_def_var( ncid, trim(var_name), nf_real, 1, var_dims, varid ), trim(err_mes) ) if( num_chem > 1 .and. traject(trj,dm)%n_chm_var > 0 ) then call def_vars( traject(trj,dm)%n_chm_var, traject(trj,dm)%chm_spc, & 'chm', is_gas=traject(trj,dm)%chm_is_gas ) endif if( num_moist > 1 .and. traject(trj,dm)%n_hyd_var > 0 ) then call def_vars( traject(trj,dm)%n_hyd_var, traject(trj,dm)%hyd_spc, 'hyd' ) endif if( num_tracer > 1 .and. traject(trj,dm)%n_trc_var > 0 ) then call def_vars( traject(trj,dm)%n_trc_var, traject(trj,dm)%trc_spc, 'trc' ) endif end do traj_loop err_mes = 'trajectory_create_file: failed to end definition for file ' // trim(filename) call handle_ncerr( nf_enddef( ncid ), trim(err_mes) ) err_mes = 'trajectory_create_file: failed to close file ' // trim(filename) call handle_ncerr( nf_close( ncid ), trim(err_mes) ) endif master_proc CONTAINS subroutine def_vars( n_vars, var_names, var_type, is_gas ) integer, intent(in) :: n_vars character(len=32), intent(in) :: var_names(:) character(len=*), intent(in) :: var_type logical, optional, intent(in) :: is_gas(:) do n = 1,n_vars write(var_name,'(a,''_traj_'',i4.4)') trim(var_names(n)),trj err_mes = 'def_vars: failed to create ' // trim(var_name) // ' variable' call handle_ncerr( nf_def_var( ncid, trim(var_name), nf_real, 1, var_dims, varid ), trim(err_mes) ) description = trim(var_name) // ' mixing ratio' call handle_ncerr( nf_put_att_text( ncid, varid, 'description', len_trim(description), description ), trim(err_mes) ) units = 'ppmv' if( present( is_gas ) ) then if( .not. is_gas(n) ) then units = 'ug/kg-dryair' endif elseif( trim(var_type) == 'hyd' ) then units = 'kg/kg' elseif( trim(var_type) == 'trc' ) then units = ' ' endif call handle_ncerr( nf_put_att_text( ncid, varid, 'units', len_trim(units), units ), trim(err_mes) ) end do end subroutine def_vars end subroutine trajectory_create_file subroutine trajectory_write_file( n_traj, n_vals, dm ) !----------------------------------------------------------------------------- ! create trajectory netcdf file !----------------------------------------------------------------------------- use module_domain !----------------------------------------------------------------------------- ! dummy arguments !----------------------------------------------------------------------------- integer, intent(in) :: n_traj integer, intent(inout) :: n_vals integer, intent(in) :: dm !----------------------------------------------------------------------------- ! local variables !----------------------------------------------------------------------------- integer :: ncid integer :: ios integer :: time_id integer :: varid integer :: l, m, n, trj, pkg integer :: time_ndx character(len=10) :: suffix(5) = (/ 'traj_i_ ', 'traj_j_ ', 'traj_k_ ', & 'traj_long_', 'traj_lat_ ' /) character(len=256) :: var_name character(len=256) :: err_mes character(len=256) :: filename include 'netcdf.inc' !--------------------------------------------------------------------- ! open netcdf file !--------------------------------------------------------------------- write(filename,'(''wrfout_traj_d'',i2.2)',iostat=ios) dm if( ios /= 0 ) then write(err_mes,'(''trajectory_write_file: failed to set filename: error = '',i6)') ios call wrf_error_fatal( trim( err_mes ) ) endif ios = nf_open( trim(filename), nf_write, ncid ) if( ios /= 0 ) then write(err_mes,'(''trajectory_write_file: failed to open '',a,'': error = '',i6)') trim(filename),ios call wrf_error_fatal( trim( err_mes ) ) endif !--------------------------------------------------------------------- ! get current time index !--------------------------------------------------------------------- err_mes = 'trajectory_write_file: failed to get time id' call handle_ncerr( nf_inq_dimid( ncid, 'time', time_id ),trim(err_mes) ) err_mes = 'trajectory_write_file: failed to get time dimension' call handle_ncerr( nf_inq_dimlen( ncid, time_id, time_ndx ),trim(err_mes) ) time_ndx = time_ndx + 1 !--------------------------------------------------------------------- ! write out trajectory times !--------------------------------------------------------------------- err_mes = 'trajectory_write_file: failed to get Times id' call handle_ncerr( nf_inq_varid( ncid, 'Times', varid ),trim(err_mes) ) err_mes = 'trajectory_write_file: failed to write Times' call handle_ncerr( nf_put_vara_text( ncid, varid, (/ 1,time_ndx /), (/ 19,n_vals /), trj_pbf(1)%times(:n_vals) ), trim(err_mes) ) traj_loop: & do m = 1,n_traj !--------------------------------------------------------------------- ! write out trajectory coordinates !--------------------------------------------------------------------- pos_loop: & do l = 1,5 write(var_name,'(a,i4.4)') trim(suffix(l)),m err_mes = 'trajectory_write_file: failed to get '// trim(var_name) // ' id' call handle_ncerr( nf_inq_varid( ncid, trim(var_name), varid ),trim(err_mes) ) err_mes = 'trajectory_write_file: failed to write ' // trim(var_name) select case( l ) case( 1 ) call handle_ncerr( nf_put_vara_real( ncid, varid, (/ time_ndx /), (/ n_vals /), & trj_pbf(m)%trj_i(:n_vals) ), trim(err_mes) ) case( 2 ) call handle_ncerr( nf_put_vara_real( ncid, varid, (/ time_ndx /), (/ n_vals /), & trj_pbf(m)%trj_j(:n_vals) ), trim(err_mes) ) case( 3 ) call handle_ncerr( nf_put_vara_real( ncid, varid, (/ time_ndx /), (/ n_vals /), & trj_pbf(m)%trj_k(:n_vals) ), trim(err_mes) ) case( 4 ) call handle_ncerr( nf_put_vara_real( ncid, varid, (/ time_ndx /), (/ n_vals /), & trj_pbf(m)%trj_lons(:n_vals) ), trim(err_mes) ) case( 5 ) call handle_ncerr( nf_put_vara_real( ncid, varid, (/ time_ndx /), (/ n_vals /), & trj_pbf(m)%trj_lats(:n_vals) ), trim(err_mes) ) end select end do pos_loop !--------------------------------------------------------------------- ! write out trajectory variables !--------------------------------------------------------------------- pkg_loop: & do pkg = 1,pkg_max select case( pkg_tag(pkg) ) case( 'chm' ) if( trjects(m)%n_chm_var > 0 ) then do n = 1,trjects(m)%n_chm_var write(var_name,'(a,''_traj_'',i4.4)') trim(trjects(m)%chm_spc(n)),m err_mes = 'trajectory_write_file: failed to get '// trim(var_name) // ' id' call handle_ncerr( nf_inq_varid( ncid, trim(var_name), varid ),trim(err_mes) ) err_mes = 'trajectory_write_file: failed to write ' // trim(var_name) call handle_ncerr( nf_put_vara_real( ncid, varid, (/ time_ndx /), (/ n_vals /), & trj_pbf(m)%chm_vals(:n_vals,n) ),trim(err_mes) ) end do endif case( 'hyd' ) if( trjects(m)%n_hyd_var > 0 ) then do n = 1,trjects(m)%n_hyd_var write(var_name,'(a,''_traj_'',i4.4)') trim(trjects(m)%hyd_spc(n)),m err_mes = 'trajectory_write_file: failed to get '// trim(var_name) // ' id' call handle_ncerr( nf_inq_varid( ncid, trim(var_name), varid ),trim(err_mes) ) err_mes = 'trajectory_write_file: failed to write ' // trim(var_name) call handle_ncerr( nf_put_vara_real( ncid, varid, (/ time_ndx /), (/ n_vals /), & trj_pbf(m)%hyd_vals(:n_vals,n) ),trim(err_mes) ) end do endif case( 'trc' ) if( trjects(m)%n_trc_var > 0 ) then do n = 1,trjects(m)%n_trc_var write(var_name,'(a,''_traj_'',i4.4)') trim(trjects(m)%trc_spc(n)),m err_mes = 'trajectory_write_file: failed to get '// trim(var_name) // ' id' call handle_ncerr( nf_inq_varid( ncid, trim(var_name), varid ),trim(err_mes) ) err_mes = 'trajectory_write_file: failed to write ' // trim(var_name) call handle_ncerr( nf_put_vara_real( ncid, varid, (/ time_ndx /), (/ n_vals /), & trj_pbf(m)%trc_vals(:n_vals,n) ),trim(err_mes) ) end do endif end select end do pkg_loop end do traj_loop n_vals = 0 ios = nf_close( ncid ) end subroutine trajectory_write_file subroutine handle_ncerr( ret, mes ) !--------------------------------------------------------------------- ! ... netcdf error handling routine !--------------------------------------------------------------------- !--------------------------------------------------------------------- ! ... dummy arguments !--------------------------------------------------------------------- integer, intent(in) :: ret character(len=*), intent(in) :: mes include 'netcdf.inc' if( ret /= nf_noerr ) then call wrf_message( trim(mes) ) call wrf_message( trim(nf_strerror(ret)) ) call wrf_abort endif end subroutine handle_ncerr #endif subroutine trajmapproj (grid,config_flags,ts_proj) use module_domain use module_llxy use module_configure, only : grid_config_rec_type, model_config_rec use module_dm, only : wrf_dm_min_real IMPLICIT NONE !------------------------------------------------------------------------ ! Arguments !------------------------------------------------------------------------ TYPE(domain), INTENT(IN) :: grid TYPE(grid_config_rec_type) , INTENT(IN) :: config_flags TYPE(PROJ_INFO), INTENT(out) :: ts_proj !------------------------------------------------------------------------ ! Local variables !------------------------------------------------------------------------ REAL :: ts_rx, ts_ry, ts_xlat, ts_xlong, ts_hgt REAL :: known_lat, known_lon INTEGER :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe TYPE (grid_config_rec_type) :: config_flags_temp config_flags_temp = config_flags call get_ijk_from_grid ( grid , & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe ) call model_to_grid_config_rec ( grid%id , model_config_rec , config_flags_temp ) !------------------------------------------------------------------------ ! Set up map transformation structure !------------------------------------------------------------------------ call map_init( ts_proj ) IF (ips <= 1 .AND. 1 <= ipe .AND. jps <= 1 .AND. 1 <= jpe) THEN known_lat = grid%xlat(1,1) known_lon = grid%xlong(1,1) ELSE known_lat = 9999. known_lon = 9999. END IF known_lat = wrf_dm_min_real(known_lat) known_lon = wrf_dm_min_real(known_lon) select case( config_flags%map_proj ) !------------------------------------------------------------------------ ! Mercator !------------------------------------------------------------------------ case( PROJ_MERC ) call map_set(PROJ_MERC, ts_proj, & truelat1 = config_flags%truelat1, & lat1 = known_lat, & lon1 = known_lon, & knowni = 1., & knownj = 1., & dx = config_flags%dx) !------------------------------------------------------------------------ ! Lambert conformal !------------------------------------------------------------------------ case( PROJ_LC ) call map_set(PROJ_LC, ts_proj, & truelat1 = config_flags%truelat1, & truelat2 = config_flags%truelat2, & stdlon = config_flags%stand_lon, & lat1 = known_lat, & lon1 = known_lon, & knowni = 1., & knownj = 1., & dx = config_flags%dx) !------------------------------------------------------------------------ ! Polar stereographic !------------------------------------------------------------------------ case( PROJ_PS ) call map_set(PROJ_PS, ts_proj, & truelat1 = config_flags%truelat1, & stdlon = config_flags%stand_lon, & lat1 = known_lat, & lon1 = known_lon, & knowni = 1., & knownj = 1., & dx = config_flags%dx) !------------------------------------------------------------------------ ! Cassini (global ARW) !------------------------------------------------------------------------ case( PROJ_CASSINI ) call map_set(PROJ_CASSINI, ts_proj, & latinc = grid%dy*360.0/(2.0*EARTH_RADIUS_M*PI), & loninc = grid%dx*360.0/(2.0*EARTH_RADIUS_M*PI), & lat1 = known_lat, & lon1 = known_lon, & ! We still need to get POLE_LAT and POLE_LON metadata variables before ! this will work for rotated poles. lat0 = 90.0, & lon0 = 0.0, & knowni = 1., & knownj = 1., & stdlon = config_flags%stand_lon) !------------------------------------------------------------------------ ! Rotated latitude-longitude !------------------------------------------------------------------------ case( PROJ_ROTLL ) call map_set(PROJ_ROTLL, ts_proj, & ixdim = grid%e_we-1, & jydim = grid%e_sn-1, & phi = real(grid%e_sn-2)*grid%dy/2.0, & lambda = real(grid%e_we-2)*grid%dx, & lat1 = config_flags%cen_lat, & lon1 = config_flags%cen_lon, & latinc = grid%dy, & loninc = grid%dx, & stagger = HH) end select end subroutine trajmapproj !++mmb ! end module module_trajectory end module trajectory !--mmb