program TRAJECT c modified january 3,1996 (glenn daughenbaugh) to adjust to pc fortran c problem of opening too many files. the file are opened when they are c needed and then closed. any changes have the 'glenn ' comment nearby. c c This program computes forward or backward trajectories c from/to any point for any u,v,w data set. c c Definitions : c l : number of points along east-west direction c m : number of points north-south direction c n : number of levels in the vertical direction c c gl : grid spacing c wl : western most longitude c sl : southern most latitude c gp : interval between vertical levels c (uniform thickness in this case) c nump : total number of trajectories c nseg : number of segments in each trajectory c kd : number of days for trajectories c ttm : total time in minutes c ttn : interval in minutes between two observational c data sets c nos : number of sequential data sets c dt : trajectory time step c is : unit no of wind data of state of trajectories. c itd : interval of trajectory data written. c td : switch for forward or backward trajectories c use td = 1 or -1 respectively c parameter (l = 17 , m = 13 , n = 10 ,gl = 2.5,gp = 100.) parameter (wl = -25., sl = -30.,bp = 100.,td = -1.,is = 41 ) parameter (kd = 6 ,kob = 4 ,itd= 6 ,dt = 10.,nump = 1 ) parameter (nseg = 721,ttm= 120.*60.,ttn = 6.*60.,nos = ttm/ttn+1) dimension u (l,m,n),v (l,m,n),w (l,m,n),dx(m),nq(nump) dimension up(l,m,n),vp(l,m,n),wp(l,m,n) dimension u1(l,m,n),v1(l,m,n),w1(l,m,n) dimension u2(l,m,n),v2(l,m,n),w2(l,m,n) dimension xs(nump,nseg),ys(nump,nseg),ps(nump,nseg) c c Open input files .file trj.d contains trajectory output c open (11,file='stp.d',status='old ',form='formatted') open (13,file='trj.d',status='unknown',form='formatted') c glenn open (14,file='traj.14',status='unknown') c c The following files contain the 3 components of the wind c field for each 6 hourly interval and for the entire considered c domain. c open(41,file='uvw.120',status='old',form='formatted') c c Compute the eastern,northern and upper limits of the domain c el = wl+gl*(float(l-1)) ul = sl+gl*(float(m-1)) ep = bp+gp*(float(n-1)) c c Compute the x-grid spacing for all latitudes . c call EPR (m,sl,gl,dx) dy = gl*111.1*1000. do 13200 np = 1,nump nq(np) = 1 13200 continue nsteps = ttn/dt c c Read data for terminal/starting point of trajectories c in degrees of lat/long and pressure in mb. c iu2 = is if (td.eq.1.) iu2 = 21 tot = 0. istp = 0 c c Read the coordinates of the initial points. c do 13202 i = 1, nump read (11,97) xs (i,1),ys (i,1),ps (i,1) write(6,97) xs (i,1),ys (i,1),ps (i,1) 13202 continue c c Read u,v and w data for all levels for initial position c do 13204 k = 1, n read (iu2,92) ((u1(i,j,k),i=1,l),j=1,m) read (iu2,92) ((v1(i,j,k),i=1,l),j=1,m) read (iu2,92) ((w1(i,j,k),i=1,l),j=1,m) 13204 continue c c Read u,v and w data for all levels for subsequent times c nosr = nos-1 do 13216 knos = 1, nosr iu2 = iu2+1*td c glenn c open and close current/previous file call openf(iu2,td) write(14,*) iu2 do 13206 k = 1, n read (iu2,92) ((u2(i,j,k),i=1,l),j=1,m) read (iu2,92) ((v2(i,j,k),i=1,l),j=1,m) read (iu2,92) ((w2(i,j,k),i=1,l),j=1,m) 13206 continue c c Computations for each trajectory step for the period c between two wind data sets. c 5 do 13208 nst = 1, nsteps tot = tot +dt istp = istp+1 c c Computation of u,v,w increments in one trajectory time step c do 13210 k = 1, 10 do 13210 j = 1, m do 13210 i = 1, l dudt = ((u2(i,j,k)-u1(i,j,k))*(dt/ttn)) dvdt = ((v2(i,j,k)-v1(i,j,k))*(dt/ttn)) dwdt = ((w2(i,j,k)-w1(i,j,k))*(dt/ttn)) c c Computation of u,v,w at next trajectory time step c u (i,j,k) = u1(i,j,k)+dudt*(nst-1) v (i,j,k) = v1(i,j,k)+dvdt*(nst-1) w (i,j,k) = w1(i,j,k)+dwdt*(nst-1) up(i,j,k) = u1(i,j,k)+dudt*(nst) vp(i,j,k) = v1(i,j,k)+dvdt*(nst) wp(i,j,k) = w1(i,j,k)+dwdt*(nst) 13210 continue c c Computation for each individual trajectory c do 13212 np = 1, nump sx1 = xs(np,istp) sy1 = ys(np,istp) sp1 = ps(np,istp) write(14,*) xs(np,istp),sx1,np,istp write(14,*) ys(np,istp),sy1,np,istp write(14,*) ps(np,istp),sp1,np,istp if(nq(np).eq.0) go to 8 c c For each trajectory for a step compute the terminal/ c and lower pressure corner of grid cube and its distances c xd , yd and pd from this point and interpolate u,v and w c data at this point. c call LOCATE (l,m,n,sx1,sy1,sp1,wl,sl,bp,gl,gp, & ir,jr,kr,xd,yd,pd,nump,np,nq) if (nq(np).eq.0) go to 8 call INTRPL (u,l,m,n,ir,jr,kr,xd,yd,pd,up1) call INTRPL (v,l,m,n,ir,jr,kr,xd,yd,pd,vp1) call INTRPL (w,l,m,n,ir,jr,kr,xd,yd,pd,wp1) c c Construct path to get next point at one time step difference c and to interpolate u,v,w data at that point and obtain mean c up,vp,wp for this segment c call DISPL (up1,vp1,wp1,sx1,sy1,sp1,dx,dy & ,gl,m,jr,dt,td,sx2,sy2,sp2) call LOCATE (l,m,n,sx2,sy2,sp2,wl,sl,bp,gl,gp, & ir,jr,kr,xd,yd,pd,nump,np,nq) if(nq(np).eq.0) go to 8 call INTRPL (up,l,m,n,ir,jr,kr,xd,yd,pd,up2) call INTRPL (vp,l,m,n,ir,jr,kr,xd,yd,pd,vp2) call INTRPL (wp,l,m,n,ir,jr,kr,xd,yd,pd,wp2) c um = (up1+up2)/2.0 vm = (vp1+vp2)/2.0 wm = (wp1+wp2)/2.0 c c Compute the path for one time step with mean u,v & w. c to get the position of this point and shift to it c call DISPL (um,vm,wm,sx1,sy1,sp1,dx,dy, & gl,m,jr,dt,td,sx2,sy2,sp2) c c house keeping work c if(sx2.lt.wl.or.sx2.gt.el) nq(np) = 0 if(sy2.lt.sl.or.sy2.gt.ul) nq(np) = 0 if(sp2.lt.bp.or.sp2.gt.ep) nq(np) = 0 8 if(nq(np).eq.1) go to 16 sx2 = xs(np,istp) sy2 = ys(np,istp) sp2 = ps(np,istp) 16 xs(np,istp+1)= sx2 ys(np,istp+1)= sy2 ps(np,istp+1)= sp2 9 continue 13212 continue 13208 continue do 13214 k = 1, n do 13214 j = 1, m do 13214 i = 1, l u1(i,j,k) = u2(i,j,k) v1(i,j,k) = v2(i,j,k) w1(i,j,k) = w2(i,j,k) 13214 continue 13216 continue istp = istp+1 c c Writting the computed data for trajectories c do 13218 i = 1, nump write (13,96) istp do 13218 k = 1, istp,itd kt = (k-1)/6 write(13,95) i,kt,ys(i,k),xs(i,k),ps(i,k) 13218 continue 90 format(3f9.2) 92 format(6e13.6) 93 format(1x,6i7) 94 format(i5,3f10.3) 95 format(2i5,3e13.6) 96 format(3i7) 97 format(3f9.2) stop end subroutine openf(iu2,td) integer iu2 real td if(td.eq.-1.) close (iu2+1) if(td.eq.1.) close (iu2-1) if(iu2.eq.21) 1 open(21,file='uvw.000',status='old',form='formatted') if(iu2.eq.22) 1 open(22,file='uvw.006',status='old',form='formatted') if(iu2.eq.23) 1 open(23,file='uvw.012',status='old',form='formatted') if(iu2.eq.24) 1 open(24,file='uvw.018',status='old',form='formatted') if(iu2.eq.25) 1 open(25,file='uvw.024',status='old',form='formatted') if(iu2.eq.26) 1 open(26,file='uvw.030',status='old',form='formatted') if(iu2.eq.27) 1 open(27,file='uvw.036',status='old',form='formatted') if(iu2.eq.28) 1 open(28,file='uvw.042',status='old',form='formatted') if(iu2.eq.29) 1 open(29,file='uvw.048',status='old',form='formatted') if(iu2.eq.30) 1 open(30,file='uvw.054',status='old',form='formatted') if(iu2.eq.31) 1 open(31,file='uvw.060',status='old',form='formatted') if(iu2.eq.32) 1 open(32,file='uvw.066',status='old',form='formatted') if(iu2.eq.33) 1 open(33,file='uvw.072',status='old',form='formatted') if(iu2.eq.34) 1 open(34,file='uvw.078',status='old',form='formatted') if(iu2.eq.35) 1 open(35,file='uvw.084',status='old',form='formatted') if(iu2.eq.36) 1 open(36,file='uvw.090',status='old',form='formatted') if(iu2.eq.37) 1 open(37,file='uvw.096',status='old',form='formatted') if(iu2.eq.38) 1 open(38,file='uvw.102',status='old',form='formatted') if(iu2.eq.39) 1 open(39,file='uvw.108',status='old',form='formatted') if(iu2.eq.40) 1 open(40,file='uvw.114',status='old',form='formatted') return end subroutine EPR (m,sl,gl,dx) c c This subroutines computes the grid spacing along c the x-axis for all latitudes. c dimension dx(m) rad = (4.*atan(1.0))/180.0 do 13200 j = 1, m dx(j) = gl*(cos(((j-1)*gl+sl)*rad))*111.1*1000.0 13200 continue return end subroutine INTRPL (u,l,m,n,ir,jr,kr,xd,yd,pd,up) c c This subroutine performs a trilinear interpolation c of the wind field components. c c three dimensional interpolation grid c 4.__._________.3 c |\ \ |\ c | \ \ | \ c | \ *d1 | \ c | 1._|\__.__|___.2 c | | | | | c | | *p | | c 8.__.|_|______.7 | c \ \ | \ | c \ |\| \ | c \| *d2 \| c .__\_________. c 5 6 c dimension u(l,m,n) u1 = u(ir ,jr ,kr ) u2 = u(ir+1,jr ,kr ) u3 = u(ir+1,jr+1,kr ) u4 = u(ir ,jr+1,kr ) u5 = u(ir ,jr ,kr+1) u6 = u(ir+1,jr ,kr+1) u7 = u(ir+1,jr+1,kr+1) u8 = u(ir ,jr+1,kr+1) ud2 = u5+xd*(u6-u5)+yd*(u8+xd*(u7-u8)-u5+xd*(u6-u5)) ud1 = u1+xd*(u2-u1)+yd*(u4+xd*(u3-u4)-u1+xd*(u2-u1)) up = ud1+pd*(ud2-ud1) return end subroutine DISPL (up,vp,wp,sx,sy,sp,dx,dy & ,gl,m,j,dt,td,sx2,sy2,sp2) c c This subroutine computes the displacement along c the trajectory in one time step c dimension dx(m) dsx2 = up*60.0*dt*gl/dx(j) dsy2 = vp*60.0*dt*gl/dy dsp2 = wp*60.0*dt sx2 = sx+dsx2*td sy2 = sy+dsy2*td sp2 = sp+dsp2*td return end subroutine LOCATE (l,m,n,sx,sy,sp,wl,sl,bp,gl,gp, & ir,jr,kr,xd,yd,pd,nump,np,nq) c c This subroutine computes the i,j,k indices of the soutwest and c lower pressure corner of the grid box and the distances of the c position from this corner along axes as a fraction of grid c length .It also ckecks if the point moves outside the grid box c and then sets nq = 0 and terminates the trajectory . c dimension nq(nump) l1 = l-1 m1 = m-1 n1 = n-1 ir = int((sx-wl)/gl)+1 xd = (sx-wl-(float(ir-1))*gl)/gl jr = int((sy-sl)/gl)+1 yd = (sy-sl-(float(jr-1))*gl)/gl kr = int((sp-bp)/gp)+1 pd = (sp-bp-(float(kr-1))*gp)/gp if (ir.lt.1.or.ir.gt.l1) nq(np) = 0 if (jr.lt.1.or.jr.gt.m1) nq(np) = 0 if (kr.lt.1.or.kr.gt.n1) nq(np) = 0 return end