subroutine trajectory ( grid,config_flags, & dt,itimestep,ru_m, rv_m, ww_m, & mut,muu,muv,c1h,c2h,c1f,c2f, & rdx, rdy, rdn, rdnw,rdzw, & traj_i,traj_j,traj_k, & traj_long,traj_lat, & xlong,xlat, & msft,msfu,msfv, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & 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 TYPE(proj_info) :: proj TYPE(domain), INTENT(IN) :: grid TYPE (grid_config_rec_type) , INTENT(IN) :: config_flags INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte INTEGER , INTENT(IN ) :: itimestep REAL , INTENT(IN ) :: rdx, & rdy, & dt REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, & rdnw REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: c1h, c2h, c1f, c2f REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , & INTENT(IN ) :: ru_m, & rv_m, & ww_m REAL , DIMENSION( ims:ime , jms:jme ) , & INTENT(IN ) ::xlong, xlat real, dimension(1:config_flags%num_traj), intent(inout) :: traj_i,traj_j,traj_k real, dimension(1:config_flags%num_traj), intent(inout) :: traj_long,traj_lat real, dimension(ims:ime,kms:kme,jms:jme),intent(in) :: rdzw 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 :: 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 real :: const1 real :: temp_i,temp_j integer :: i_u,j_v,k_w real :: u_a,u_b,u_c,u_d,v_a,v_b,v_c,v_d,w_a,w_b,w_c,w_d real :: d_a,d_b,d_c,d_d real :: u_temp_upper,u_temp_lower real :: v_temp_upper,v_temp_lower 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 real:: known_lat, known_lon 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) const1=1.0/2.0/sqrt(2.0) 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 ) traj_loop: & do tjk = 1,traj_cnt(dm) 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)) 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 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) 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 if (traj_i(tjk)-real(floor(traj_i(tjk))) .ge. 0.5 ) then i_u=floor(traj_i(tjk)) + 1 else i_u=floor(traj_i(tjk)) endif if (k_traj .ge. 1 ) then u_a=u(i_u ,k_traj,j_traj+1) u_b=u(i_u ,k_traj,j_traj ) u_c=u(i_u+1,k_traj,j_traj ) u_d=u(i_u+1,k_traj,j_traj+1) else u_a=0.0 u_b=0.0 u_c=0.0 u_d=0.0 endif d_a=abs((real(i_u+1)-(traj_i(tjk)+0.5))*(real(j_traj )-traj_j(tjk))) d_b=abs((real(i_u+1)-(traj_i(tjk)+0.5))*(real(j_traj+1)-traj_j(tjk))) d_c=abs((real(i_u )-(traj_i(tjk)+0.5))*(real(j_traj+1)-traj_j(tjk))) d_d=abs((real(i_u )-(traj_i(tjk)+0.5))*(real(j_traj )-traj_j(tjk))) u_temp_lower=(u_a*d_a+u_b*d_b+u_c*d_c+u_d*d_d)/(d_a+d_b+d_c+d_d) u_a=u(i_u ,k_traj+1,j_traj+1) u_b=u(i_u ,k_traj+1,j_traj ) u_c=u(i_u+1,k_traj+1,j_traj ) u_d=u(i_u+1,k_traj+1,j_traj+1) d_a=abs((real(i_u+1)-(traj_i(tjk)+0.5))*(real(j_traj )-traj_j(tjk))) d_b=abs((real(i_u+1)-(traj_i(tjk)+0.5))*(real(j_traj+1)-traj_j(tjk))) d_c=abs((real(i_u )-(traj_i(tjk)+0.5))*(real(j_traj+1)-traj_j(tjk))) d_d=abs((real(i_u )-(traj_i(tjk)+0.5))*(real(j_traj )-traj_j(tjk))) u_temp_upper=(u_a*d_a+u_b*d_b+u_c*d_c+u_d*d_d)/(d_a+d_b+d_c+d_d) traj_u=u_temp_upper*abs(real(k_traj)-traj_k(tjk))+u_temp_lower*abs(real(k_traj+1)-traj_k(tjk)) if (traj_j(tjk)-real(floor(traj_j(tjk))) .ge. 0.5 ) then j_v=floor(traj_j(tjk)) + 1 else j_v=floor(traj_j(tjk)) endif if (k_traj .ge. 1 ) then v_a=v(i_traj ,k_traj,j_v+1) v_b=v(i_traj ,k_traj,j_v ) v_c=v(i_traj+1,k_traj,j_v ) v_d=v(i_traj+1,k_traj,j_v+1) else v_a=0.0 v_b=0.0 v_c=0.0 v_d=0.0 endif d_a=abs((real(i_traj+1)-traj_i(tjk))*(real(j_v )-(traj_j(tjk)+0.5))) d_b=abs((real(i_traj+1)-traj_i(tjk))*(real(j_v+1)-(traj_j(tjk)+0.5))) d_c=abs((real(i_traj )-traj_i(tjk))*(real(j_v+1)-(traj_j(tjk)+0.5))) d_d=abs((real(i_traj )-traj_i(tjk))*(real(j_v )-(traj_j(tjk)+0.5))) v_temp_lower=(v_a*d_a+v_b*d_b+v_c*d_c+v_d*d_d)/(d_a+d_b+d_c+d_d) v_a=v(i_traj ,k_traj+1,j_v+1) v_b=v(i_traj ,k_traj+1,j_v ) v_c=v(i_traj+1,k_traj+1,j_v ) v_d=v(i_traj+1,k_traj+1,j_v+1) d_a=abs((real(i_traj+1)-traj_i(tjk))*(real(j_v )-(traj_j(tjk)+0.5))) d_b=abs((real(i_traj+1)-traj_i(tjk))*(real(j_v+1)-(traj_j(tjk)+0.5))) d_c=abs((real(i_traj )-traj_i(tjk))*(real(j_v+1)-(traj_j(tjk)+0.5))) d_d=abs((real(i_traj )-traj_i(tjk))*(real(j_v )-(traj_j(tjk)+0.5))) v_temp_upper=(v_a*d_a+v_b*d_b+v_c*d_c+v_d*d_d)/(d_a+d_b+d_c+d_d) traj_v=v_temp_upper*abs(real(k_traj)-traj_k(tjk))+v_temp_lower*abs(real(k_traj+1)-traj_k(tjk)) if (traj_k(tjk)-real(floor(traj_k(tjk))) .ge. 0.5 ) then k_w=floor(traj_k(tjk)) + 1 else k_w=floor(traj_k(tjk)) endif if (k_w .ge. 1) then w_b=w(i_traj ,k_w ,j_traj) w_c=w(i_traj+1,k_w ,j_traj) else w_b=0.0 w_c=0.0 endif w_a=w(i_traj ,k_w+1,j_traj) w_d=w(i_traj+1,k_w+1,j_traj) d_a=abs((real(i_traj+1)-traj_i(tjk))*(real(k_w )-(traj_k(tjk)+0.5))) d_b=abs((real(i_traj+1)-traj_i(tjk))*(real(k_w+1)-(traj_k(tjk)+0.5))) d_c=abs((real(i_traj )-traj_i(tjk))*(real(k_w+1)-(traj_k(tjk)+0.5))) d_d=abs((real(i_traj )-traj_i(tjk))*(real(k_w )-(traj_k(tjk)+0.5))) w_temp_lower=(w_a*d_a+w_b*d_b+w_c*d_c+w_d*d_d)/(d_a+d_b+d_c+d_d) if (k_w .ge. 1) then w_b=w(i_traj ,k_w ,j_traj+1) w_c=w(i_traj+1,k_w ,j_traj+1) else w_b=0.0 w_c=0.0 endif w_a=w(i_traj ,k_w+1,j_traj+1) w_d=w(i_traj+1,k_w+1,j_traj+1) d_a=abs((real(i_traj+1)-traj_i(tjk))*(real(k_w )-(traj_k(tjk)+0.5))) d_b=abs((real(i_traj+1)-traj_i(tjk))*(real(k_w+1)-(traj_k(tjk)+0.5))) d_c=abs((real(i_traj )-traj_i(tjk))*(real(k_w+1)-(traj_k(tjk)+0.5))) d_d=abs((real(i_traj )-traj_i(tjk))*(real(k_w )-(traj_k(tjk)+0.5))) w_temp_upper=(w_a*d_a+w_b*d_b+w_c*d_c+w_d*d_d)/(d_a+d_b+d_c+d_d) traj_w=w_temp_upper*abs(real(j_traj)-traj_j(tjk))+w_temp_lower*abs(real(j_traj+1)-traj_j(tjk)) eta_old=grid%znw(k_w+1)*abs(traj_k(tjk)+0.5-real(k_w))+grid%znw(k_w)*abs(traj_k(tjk)+0.5-real(k_w+1)) rdx_grid=rdx*msft(i_traj,j_traj) rdy_grid=rdy*msft(i_traj,j_traj) rdz_grid=rdnw(k_traj) deltx=traj_u*DT delty=traj_v*DT deltz=traj_w*DT traj_i(tjk)=traj_i(tjk)+deltx*rdx_grid traj_j(tjk)=traj_j(tjk)+delty*rdy_grid eta_new=eta_old+deltz keta_temp = 0 do keta=1, kme-1 if (eta_new .le. grid%znw(keta) .and. eta_new .gt. grid%znw(keta+1)) then keta_temp=keta endif enddo if (keta_temp .eq. 0) then traj_k(tjk) = traj_k(tjk) else traj_k(tjk) = (real(keta_temp)*abs(eta_new-grid%znw(keta_temp+1))+ & real(keta_temp+1)*abs(eta_new-grid%znw(keta_temp))) & /(grid%znw(keta_temp)-grid%znw(keta_temp+1)) traj_k(tjk) = traj_k(tjk)-0.5 endif 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 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 SUBROUTINE trajectory