Skip to content

Commit

Permalink
ARW Trajectory Diagnostics (wrf-model#84)
Browse files Browse the repository at this point in the history 
KEYWORDS: trajectory

SOURCE: internal

DESCRIPTION OF CHANGES:

The original "trajectory" code in routine module_em.F only computed the trajectory
spatial coordinates per trajectory, per time step.  The new code, in share/module_trajectory.F,
takes trajectory coordinate values and linearly interpolates in x,y, and z to get the
value of up 100 variables per trajectory per time step.  The code will detect and handle
trajectory coordinate values that are between the grid parcels.  The user may have
up to 1000 trajectories per domain and as mentioned up to 100 variables per trajectory.

LIST OF MODIFIED FILES:

M       Registry/Registry.EM_COMMON
M       Registry/Registry.NMM
M       chem/chem_driver.F
M       chem/depend.chem
M       dyn_em/module_em.F
M       dyn_em/solve_em.F
M       dyn_em/start_em.F
M       external/RSL_LITE/module_dm.F
M       main/depend.common
M       share/Makefile
M       share/module_check_a_mundo.F
A       share/module_trajectory.F
M       share/solve_interface.F

TESTS CONDUCTED:

1. WTF reg test, version 3.04, has been successfully completed. 
2. The code has been tested in ACOM by Mary Barth,  Megan Bella, and Gabrielle Pfister.
  • Loading branch information
@stacywalters @davegill
stacywalters authored and davegill committed Jan 12, 2017
1 parent 0f00166 commit 8cb2192
Show file tree
Hide file tree
Showing 13 changed files with 3,165 additions and 49 deletions.
3 changes: 2 additions & 1 deletion Registry/Registry.EM_COMMON
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2094,7 +2094,7 @@ rconfig integer ocean_levels namelist,domains 1 30 i
rconfig real ocean_z namelist,domains max_ocean -1 - "vertical profile of layer depths for ocean" "m"
rconfig real ocean_t namelist,domains max_ocean -1 - "vertical profile of ocean temps" "K"
rconfig real ocean_s namelist,domains max_ocean -1 - "vertical profile of salinity"
rconfig integer num_traj namelist,domains 1 25 irh "num_traj" "#of trajectory" ""
rconfig integer num_traj namelist,domains 1 1000 irh "num_traj" "#of trajectory" ""

# variable for time series of vertical profile of U, V, Theta, GHT,a nd QVAPOR
rconfig integer max_ts_level namelist,domains 1 15 - "max_ts_level" "Highest model level for time series output"
Expand Down Expand Up @@ -2263,6 +2263,7 @@ rconfig integer pxlsm_smois_init namelist,physics max_domains 1
rconfig integer omlcall namelist,physics 1 0 h "omlcall" "temporary holder to allow checking for new name: oml_opt"
rconfig integer sf_ocean_physics namelist,physics 1 0 h "sf_ocean_physics" "activate ocean model 0=no, 1=1d mixed layer, 2=3D PWP" ""
rconfig integer traj_opt namelist,physics 1 0 h "traj_opt" "activate trajectory calculation 0=no, 1=on" ""
rconfig logical dm_has_traj namelist,physics max_domains .false. rh "has_traj" "activate trajectory calculation per domain" ""
rconfig integer tracercall namelist,physics 1 0 h "tracercall" "activate tracer calculation 0=no, 1=on" ""
rconfig real OMDT namelist,physics 1 1 h "OMDT" "Timestep of ocean model" "s"
rconfig real oml_hml0 namelist,physics 1 50 h "oml_hml0" "oml initial mixed layer depth value" "m"
Expand Down
5 changes: 5 additions & 0 deletions Registry/Registry.NMM
View file
Original file line number Diff line number Diff line change
Expand Up @@ -162,6 +162,11 @@ ifdef HWRF=1
include registry.tracker
endif

# For compilation only; any setting other than zero will cause the simulation to halt
rconfig integer traj_opt namelist,physics 1 0 h "traj_opt" "activate trajectory calculation 0=no, 1=on" ""
rconfig logical dm_has_traj namelist,physics max_domains .false. rh "has_traj" "activate trajectory calculation per domain" ""
rconfig integer num_traj namelist,domains 1 1000 irh "num_traj" "#of trajectory" ""

# Nest motion safeguard: don't let nest get close to parent boundary.
# Default values are lowest possible - anything lower would read
# outside of memory in intermediate domain.
Expand Down
9 changes: 8 additions & 1 deletion chem/chem_driver.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,7 @@ subroutine chem_driver ( grid , config_flags &
pcnst =>pcnst_runtime, numgas_mam, cam_mam_aerosols
USE module_cu_camzm_driver, only: zm_conv_tend_2
USE module_cam_mam_gas_wetdep_driver, only: cam_mam_gas_wetdep_driver
USE module_trajectory, only: trajectory_dchm_tstep_init, trajectory_dchm_tstep_set

IMPLICIT NONE

Expand Down Expand Up @@ -1135,7 +1136,9 @@ end SUBROUTINE sum_pm_driver
! Calculate rate of n2o5 hydrolysis
call wrf_debug(15,'calling calc_het_n2o5')


write(msg,'(''chem_driver('',i2.2,''): Calling dchm_tstep_init'')') grid%id
call wrf_debug( 200,trim(msg) )
call trajectory_dchm_tstep_init( grid, do_chemstep )

!
! For the chemistry tracer mode, only emissions and vertical mixing are done.
Expand Down Expand Up @@ -1226,6 +1229,10 @@ end SUBROUTINE sum_pm_driver
! chem(its:ite,kts:kte,jts:jte,p_h2o2)=vch2o2_old(its:ite,kts:kte,jts:jte)
endif

write(msg,'(''chem_driver('',i2.2,''): Calling dchm_tstep_set'')') grid%id
call wrf_debug( 200,trim(msg) )
call trajectory_dchm_tstep_set( grid )

IF(config_flags%conv_tr_aqchem == 0 ) THEN
so2so4_selecta: SELECT CASE(config_flags%chem_opt)
CASE (RADM2SORG,RADM2SORG_KPP,RACMSORG_KPP,RACM_SOA_VBS_KPP)
Expand Down
2 changes: 1 addition & 1 deletion chem/depend.chem
View file
Original file line number Diff line number Diff line change
Expand Up @@ -280,7 +280,7 @@ module_tropopause.o: module_interpolate.o

module_upper_bc_driver.o: module_tropopause.o

chem_driver.o: module_radm.o ../dyn_em/module_convtrans_prep.o module_chem_utilities.o module_data_radm2.o module_dep_simple.o module_bioemi_simple.o module_vertmx_wrf.o module_phot_mad.o module_aerosols_sorgam.o module_aerosols_soa_vbs.o module_aerosols_sorgam_vbs.o module_data_cbmz.o module_cbmz.o module_wetscav_driver.o dry_dep_driver.o emissions_driver.o module_input_tracer.o module_input_tracer_data.o module_tropopause.o module_upper_bc_driver.o module_ctrans_grell.o module_data_soa_vbs.o module_aer_opt_out.o module_data_sorgam.o module_gocart_so2so4.o ../phys/module_cu_camzm_driver.o module_cam_mam_gas_wetdep_driver.o module_dust_load.o module_chem_cup.o
chem_driver.o: module_radm.o ../dyn_em/module_convtrans_prep.o module_chem_utilities.o module_data_radm2.o module_dep_simple.o module_bioemi_simple.o module_vertmx_wrf.o module_phot_mad.o module_aerosols_sorgam.o module_aerosols_soa_vbs.o module_aerosols_sorgam_vbs.o module_data_cbmz.o module_cbmz.o module_wetscav_driver.o dry_dep_driver.o emissions_driver.o module_input_tracer.o module_input_tracer_data.o module_tropopause.o module_upper_bc_driver.o module_ctrans_grell.o module_data_soa_vbs.o module_aer_opt_out.o module_data_sorgam.o module_gocart_so2so4.o ../phys/module_cu_camzm_driver.o module_cam_mam_gas_wetdep_driver.o module_dust_load.o module_chem_cup.o ../share/module_trajectory.o

aerosol_driver.o: module_data_sorgam.o module_aerosols_sorgam.o module_data_soa_vbs.o module_aerosols_soa_vbs.o module_aerosols_sorgam_vbs.o module_mosaic_driver.o

Expand Down
134 changes: 90 additions & 44 deletions dyn_em/module_em.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2237,6 +2237,12 @@ subroutine trajectory ( grid,config_flags, &
its, ite, jts, jte, kts, kte )

!---------------------------------------------------------------------------------------------------

use module_date_time
use module_utility
use module_domain, only : domain_clock_get
use module_trajectory, only : traject, traj_cnt

implicit none
!---------------------------------------------------------------------------------------------------
! Subroutine trajectory calculates forward Lagrangian trajectory
Expand Down Expand Up @@ -2278,7 +2284,10 @@ subroutine trajectory ( grid,config_flags, &
real, dimension(ims:ime,kms:kme,jms:jme)::u,v,w
real, dimension(ims:ime,jms:jme),intent(in)::msft,msfu,msfv
real, dimension(ims:ime,jms:jme),intent(in)::muu,muv,mut
integer :: i_traj,j_traj,k_traj,tjk,k
integer :: j,k
integer :: i_beg,i_end
integer :: i_traj,j_traj,k_traj,tjk
integer :: dm
real :: traj_u,traj_v,traj_w
real :: rdx_grid,rdy_grid,rdz_grid
real :: deltx, delty, deltz,ax
Expand All @@ -2292,46 +2301,90 @@ subroutine trajectory ( grid,config_flags, &
real :: w_temp_upper,w_temp_lower
real :: eta_old, eta_new
integer :: keta, keta_temp
character(len=19) :: current_timestr, next_timestr, wrk_timestr
character(len=256) :: dbg_mes
logical :: has_proj_map
! varalbe for map projectory
real:: known_lat, known_lon
TYPE (grid_config_rec_type) :: config_flags_temp
TYPE (grid_config_rec_type) :: config_flags_temp
TYPE(WRFU_Time) :: current_time, next_time
TYPE(WRFU_Time) :: start_time, stop_time

config_flags_temp = config_flags
call trajmapproj(grid, config_flags_temp,proj)


! convert ru_m, rv_m and ww_m in u,v,w
const1=1.0/2.0/sqrt(2.0)
do k=kms,kme
u(:,k,:)=ru_m(:,k,:)/muu(:,:)*msfu(:,:)
v(:,k,:)=rv_m(:,k,:)/muv(:,:)*msfv(:,:)
w(:,k,:)=ww_m(:,k,:)/mut(:,:)*msft(:,:)
enddo
do tjk = 1,config_flags%num_traj
eta_old = 0.0
eta_new = 0.0
keta=0
keta_temp=0
if (traj_i(tjk) .ne. -9999.0) then
if ( ( proj%code .EQ. PROJ_LC ) .OR. &
( proj%code .EQ. PROJ_PS ) .OR. &

dm = grid%id

write(dbg_mes,'(''trajectory('',i2.2,''): Entering trajectory'')') dm
call wrf_debug(200,trim(dbg_mes) )

i_beg = max( 1,its-1 )
i_end = min( ite+2,ide-1 )
do j=max(1,jts-1),min(jte+2,jde-1)
do k=kms,kme-1
u(its:i_end,k,j)=ru_m(its:i_end,k,j)/MUU(its:i_end,j)*msfu(its:i_end,j)
v(its:i_end,k,j)=rv_m(its:i_end,k,j)/MUV(its:i_end,j)*msfv(its:i_end,j)
enddo
do k=kms,kme
w(its:i_end,k,j)=ww_m(its:i_end,k,j)/MUT(its:i_end,j)*msft(its:i_end,j)
enddo
enddo
has_proj_map = &
( proj%code .EQ. PROJ_LC ) .OR. &
( proj%code .EQ. PROJ_PS_WGS84 ) .OR. &
( proj%code .EQ. PROJ_ALBERS_NAD83 ) .OR. &
( proj%code .EQ. PROJ_MERC ) .OR. &
( proj%code .EQ. PROJ_LATLON ) .OR. &
( proj%code .EQ. PROJ_CYL ) .OR. &
( proj%code .EQ. PROJ_CASSINI ) .OR. &
( proj%code .EQ. PROJ_GAUSS ) .OR. &
( proj%code .EQ. PROJ_ROTLL ) ) THEN
call latlon_to_ij (proj, traj_lat(tjk),traj_long(tjk),traj_i(tjk),traj_j(tjk))
( proj%code .EQ. PROJ_ROTLL )

traj_loop: &
do tjk = 1,traj_cnt(dm)
!do tjk = 1,config_flags%num_traj
traj_is_active: &
if (traj_i(tjk) .ne. -9999.0) then
eta_old = 0.0
eta_new = 0.0
keta = 0
keta_temp = 0
if( has_proj_map ) then
call latlon_to_ij (proj, &
traj_lat(tjk),traj_long(tjk),traj_i(tjk),traj_j(tjk))
end if
i_traj=floor(traj_i(tjk)) ! find the lower_left_bottom corner for trajectory
j_traj=floor(traj_j(tjk)) !
k_traj=floor(traj_k(tjk)) !
if ((i_traj .ge. its .and. i_traj .le. ite .and. i_traj .lt. ide) .and. &
(j_traj .ge. jts .and. j_traj .le. jte .and. j_traj .lt. jde) .and. &
(k_traj .le. kte .and. k_traj .lt. kde)) then

i_traj = floor(traj_i(tjk)) ! find the lower_left_bottom corner for
!trajectory
j_traj = floor(traj_j(tjk)) !
k_traj = floor(traj_k(tjk)) !
traj_in_domain: &
if ((i_traj .ge. its .and. i_traj .le. ite .and. i_traj .lt. ide) .and. &
(j_traj .ge. jts .and. j_traj .le. jte .and. j_traj .lt. jde) .and. &
(k_traj .le. kte .and. k_traj .lt. kde)) then
!-----------------------------------------------------------------------------
! 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 )
wrk_timestr(1:19) = traject(tjk,dm)%start_time(1:19)
! write(*,*) ' '
! write(*,*) traject(tjk,dm)
! write(dbg_mes,'(''trajectory('',i2,2,''): tjk,start_time = '',i3,1x,a)') &
! dm,tjk,wrk_timestr
! call wrf_debug( 200,trim(dbg_mes) )
call wrf_atotime( wrk_timestr(1:19), start_time )
wrk_timestr(1:19) = traject(tjk,dm)%stop_time(1:19)
call wrf_atotime( wrk_timestr(1:19), stop_time )
is_in_time_interval: &
if( next_time .ge. start_time .and. next_time .le. stop_time &
.and. .not. traject(tjk,dm)%is_stationary ) then
! for u : check x stagger
if (traj_i(tjk)-real(floor(traj_i(tjk))) .ge. 0.5 ) then
i_u=floor(traj_i(tjk)) + 1
Expand Down Expand Up @@ -2472,32 +2525,25 @@ subroutine trajectory ( grid,config_flags, &
traj_k(tjk) = traj_k(tjk)-0.5
endif
!! convert i,j,k into lon, lat
if ( ( proj%code .EQ. PROJ_LC ) .OR. &
( proj%code .EQ. PROJ_PS ) .OR. &
( proj%code .EQ. PROJ_PS_WGS84 ) .OR. &
( proj%code .EQ. PROJ_ALBERS_NAD83 ) .OR. &
( proj%code .EQ. PROJ_MERC ) .OR. &
( proj%code .EQ. PROJ_LATLON ) .OR. &
( proj%code .EQ. PROJ_CYL ) .OR. &
( proj%code .EQ. PROJ_CASSINI ) .OR. &
( proj%code .EQ. PROJ_GAUSS ) .OR. &
( proj%code .EQ. PROJ_ROTLL ) ) THEN
call ij_to_latlon (proj, traj_i(tjk), traj_j(tjk),traj_lat(tjk),traj_long(tjk))
end if
else
if( has_proj_map ) then
call ij_to_latlon (proj, traj_i(tjk), &
traj_j(tjk),traj_lat(tjk),traj_long(tjk))
endif
endif is_in_time_interval
else traj_in_domain
traj_i(tjk) = -9999.0
traj_j(tjk) = -9999.0
traj_k(tjk) = -9999.0
traj_long(tjk) = -9999.0
traj_lat(tjk) = -9999.0
endif
endif
traj_i(tjk) = wrf_dm_max_real(traj_i(tjk))! save the information from highest fomain
traj_j(tjk) = wrf_dm_max_real(traj_j(tjk))
traj_k(tjk) = wrf_dm_max_real(traj_k(tjk))
traj_long(tjk) = wrf_dm_max_real(traj_long(tjk))
traj_lat(tjk) = wrf_dm_max_real(traj_lat(tjk))
enddo
endif traj_in_domain
endif traj_is_active
traj_i(tjk) = wrf_dm_max_real(traj_i(tjk))
traj_j(tjk) = wrf_dm_max_real(traj_j(tjk))
traj_k(tjk) = wrf_dm_max_real(traj_k(tjk))
traj_long(tjk) = wrf_dm_max_real(traj_long(tjk))
traj_lat(tjk) = wrf_dm_max_real(traj_lat(tjk))
enddo traj_loop
!end trajectory

END SUBROUTINE trajectory
Expand Down
2 changes: 1 addition & 1 deletion dyn_em/solve_em.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3418,7 +3418,7 @@ SUBROUTINE solve_em ( grid , config_flags &
call trajectory (grid,config_flags, &
grid%dt,grid%itimestep,grid%ru_m, grid%rv_m, grid%ww_m,&
grid%mut,grid%muu,grid%muv, &
grid%muts,muus,muvs, &
grid%rdx, grid%rdy, grid%rdn, grid%rdnw,grid%rdzw, &
grid%traj_i,grid%traj_j,grid%traj_k, &
grid%traj_long,grid%traj_lat, &
Expand Down
7 changes: 7 additions & 0 deletions dyn_em/start_em.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@ SUBROUTINE start_domain_em ( grid, allowed_to_read &
USE module_lightning_driver, ONLY : lightning_init
USE module_fr_fire_driver_wrf, ONLY : fire_driver_em_init
USE module_stoch, ONLY : setup_rand_perturb, rand_seed, update_stoch, initialize_stoch
USE module_trajectory, ONLY : trajectory_init
#if (WRF_CHEM == 1)
USE module_aerosols_sorgam, ONLY: sum_pm_sorgam
USE module_gocart_aerosols, ONLY: sum_pm_gocart
Expand Down Expand Up @@ -1909,6 +1910,12 @@ SUBROUTINE start_domain_em ( grid, allowed_to_read &
CALL wrf_debug ( 100 , 'start_domain_em: After call to fire_driver_em_init' )
endif
if( grid%traj_opt /= no_trajectory ) then
call trajectory_init( grid, config_flags, &
ims,ime, jms,jme, kms,kme )
CALL wrf_debug ( 100 , 'start_domain_em: After call to trajectory_init' )
endif
CALL wrf_debug ( 100 , 'start_domain_em: Returning' )
Expand Down
14 changes: 14 additions & 0 deletions external/RSL_LITE/module_dm.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1295,6 +1295,20 @@ INTEGER function getrealmpitype()
RETURN
END FUNCTION getrealmpitype
REAL FUNCTION wrf_dm_max_int ( inval )
IMPLICIT NONE
#ifndef STUBMPI
INCLUDE 'mpif.h'
INTEGER, intent(in) :: inval
INTEGER :: ierr, retval
CALL mpi_allreduce ( inval, retval , 1, MPI_INT, MPI_MAX, local_communicator, ierr )
wrf_dm_max_int = retval
#else
INTEGER, intent(in) :: inval
wrf_dm_max_int = inval
#endif
END FUNCTION wrf_dm_max_int
REAL FUNCTION wrf_dm_max_real ( inval )
IMPLICIT NONE
#ifndef STUBMPI
Expand Down
11 changes: 10 additions & 1 deletion main/depend.common
View file
Original file line number Diff line number Diff line change
Expand Up @@ -873,9 +873,18 @@ module_ra_aerosol.o :\

# DEPENDENCIES for share

module_trajectory.o: ../frame/module_domain.o \
../frame/module_configure.o \
../frame/module_dm.o \
../frame/module_comm_dm.o \
../frame/module_state_description.o \
module_model_constants.o \
module_date_time.o \
module_llxy.o

solve_interface.o: solve_em.int ../frame/module_domain.o ../frame/module_configure.o \
../frame/module_timing.o ../frame/module_driver_constants.o \
../frame/module_wrf_error.o
../frame/module_wrf_error.o module_trajectory.o

start_domain.o: start_domain_em.int wrf_timeseries.o track_driver.o ../frame/module_domain.o ../frame/module_configure.o ../share/module_llxy.o

Expand Down
1 change: 1 addition & 0 deletions share/Makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,7 @@ OBJS3 = \
landread.o \
track_driver.o \
track_input.o \
module_trajectory.o \
bobrand.o \
wrf_timeseries.o \
track_driver.o \
Expand Down
7 changes: 7 additions & 0 deletions share/module_check_a_mundo.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -584,6 +584,13 @@ SUBROUTINE check_nml_consistency
model_config_rec%num_traj = 0
END IF

#elif( NMM_CORE == 1 )
!----------------------------------------------------------------------------
! If NMM core and trajectories are on then halt.
!----------------------------------------------------------------------------
IF ( model_config_rec%traj_opt /= 0 ) THEN
call wrf_error_fatal( 'Trajectories not supported in NMM core' )
END IF
#endif

#if (EM_CORE == 1)
Expand Down
Loading

0 comments on commit 8cb2192

Please sign in to comment.