pro mtrajmain,startarr,yrarr,ndays,latin,lonin,dt,forward, $ latsave,lonsave,qsave,temsave,tsave,yrsave,szsave,psave,err, $ source=source, noprint=noprint,theta=theta,$ save_memory=save_memory,usave=usave,vsave=vsave,plot=plot,dir=dir,scale=scale0,$ latpole=latpole,lonpole=lonpole,segment=segment,noplot=noplot ;+ ; NAME: ; mtrajmain ; PURPOSE: ; Trajectory Code Main Program ; CATEGORY: ; trajectory ;- common metfiles,dayfracs,lats,lons noprint=1 if forward eq 1 then start=min(startarr) else start=max(startarr) dayni=start nparcels=n_elements(latin) if n_elements(save_memory) eq 0 then save_memory=1 ; set up initial constants asteps=long(0) yrn=yrarr dtd=dt*86400. c1=!pi/180. c2=1./(6.37e6*c1) nstep=fix(ndays/dt) ; number of time steps ; if save_memory option then adjust the number of steps nstep2=nstep if save_memory ne 1 then $ nstep2 = save_memory*(fix(float(nstep)/save_memory)) if nstep2 ne nstep then begin print,'*** Warning **** Integration length extended due to conserve memory request' nstep = nstep2+1 endif if nstep eq 0 then begin print,'Integration length smaller than a time step' err=1 return endif mtraj_selectdata,start,ndays,dayfracs,forward=forward ; compute maximum integration length if err ne 0 then return ; set up initial vectors sm=nstep/save_memory+1 tsave=fltarr(sm) yrsave=intarr(sm) latsave=fltarr(n_elements(latin),sm) lonsave=latsave temsave=latsave qsave=latsave psave=latsave if n_elements(vsave) ne 0 then vsave=latsave if n_elements(usave) ne 0 then usave=latsave ; load initial fields latsave(0,0)=latin lonsave(0,0)=lonin curlat=latin curlon=lonin sth=size(theta) yrsave(0)=fix(yrn(0)) tsave(0)=dayni-day1900(string(yrn(0),'(i2)')+'0101')+1 mtraj_getdata,dayni,curlat,curlon,uout,vout,tout,qout,err,$ forward=forward,$ dayfracs1=dayfracs1,u1=u1,v1=v1,q1=q1,t1=t1,$ dayfracs2=dayfracs2,u2=u2,v2=v2,q2=q2,t2=t2,dir=dir temsave(0,0)=tout pout=1000.*(tout/theta)^(3.5) psave(0,0)=pout qsave(0,0)=qout ; begin print out of initial data if nparcels gt 100 then $ if not keyword_set(noprint) then print,'** Truncating the parcel print out list due to length' if not keyword_set(noprint) then print,' Initial Parcel Data' if n_elements(dif_start) le 1 then begin ; all the same start times if not keyword_set(noprint) then print,' No. Theta (K) Lat. Lon. Height(km) T PV' jxmax=nparcels if nparcels gt 100 then begin delta = nparcels/100. nfx=long(indgen(100))*long(delta) jxmax=n_elements(nfx) endif else nfx=indgen(nparcels) for jx=0L,jxmax-1 do begin th5 = theta if not keyword_set(noprint) then print,nfx(jx),th5,curlat(nfx(jx)),curlon(nfx(jx)),7.0*alog(1000./pout(nfx(jx))),$ tout(nfx(jx)),qout(nfx(jx)),format='(I7,F10.1,4(F10.1),E12.3)' endfor endif else begin ; different start times if not keyword_set(noprint) then print,' # YY:DDD Theta (K) Lat. Lon. Height(km) T PV' jxmax=nparcels if n_elements(yrn) lt 1024 then begin txx=float(startarr-day1900(string(yrn,'(i2)')+'0101')+1) endif else begin txx=startarr*0. for jtx=0L,n_elements(startarr)-1 do $ txx(jtx)=float(startarr(jtx))-day1900(string(yrn(jtx),'(i2)')+'0101')+1 endelse yxx=fix(yrn) if nparcels gt 100 then begin delta = nparcels/100. nfx=long(indgen(100))*long(delta) jxmax=n_elements(nfx) endif else nfx=indgen(nparcels) for jx=0L,jxmax-1 do begin th5 = theta ttim=' '+str(yxx(nfx(jx)))+':'+strtrim(string(txx(nfx(jx)),format='(f6.2)'),2)+' ' if not keyword_set(noprint) then print,nfx(jx),ttim,th5,curlat(nfx(jx)),curlon(nfx(jx)),7.0*alog(1000./pout(nfx(jx))),$ tout(nfx(jx)),qout(nfx(jx)),format='(I7,A7,F10.1,4(F10.1),E12.3)' endfor endelse if not keyword_set(noprint) then print,' ' ; begin time stepping il=1 for tsteps=1,nstep do begin dayni=dayni+forward*dt xcurlat=curlat xcurlon=curlon mtraj_runge,newlat,newlon,curlat,curlon,dayni,err,$ forward=forward,$ timestep=dt,uout=uout,vout=vout,tout=tout,$ dayfracs1=dayfracs1,u1=u1,v1=v1,q1=q1,t1=t1,$ dayfracs2=dayfracs2,u2=u2,v2=v2,q2=q2,t2=t2,dir=dir if err ne 0 then goto,error1 curlat=newlat curlon=newlon ; for parcels that haven't started copy their old positions into the new positions ; so that they looked fixed to the user if forward gt 0 then begin ns=where((startarr gt dayni)) endif else begin ns=where((startarr lt dayni)) endelse if ns(0) ge 0 then begin curlat(ns)=xcurlat(ns) curlon(ns)=xcurlon(ns) endif ; ----- now get u,v,t,p and pv at new positions mtraj_getdata,dayni,curlat,curlon,uout1,vout1,$ tout1,qout1,err,$ forward=forward,$ dayfracs1=dayfracs1,u1=u1,v1=v1,q1=q1,t1=t1,$ dayfracs2=dayfracs2,u2=u2,v2=v2,q2=q2,t2=t2,dir=dir if err ne 0 then goto,error1 vout=vout1 uout=uout1 tout=tout1 qout=qout1 pout=1000.*(tout/theta)^(3.5) ; okay we are basically done with the calculation - now plot values timer=tsteps*dt if not keyword_set(noplot) then mtraj_plot,curlat,curlon,theta,dayni,pout,tout,qout,timer=timer,scale=scale0,$ latpole=latpole,lonpole=lonpole,segment=segment ; save variables if (tsteps mod save_memory eq 0) then begin ; save memory by "mod-ing" output asteps=asteps+1 ; actual counter if n_elements(dif_start) le 1 then begin tsave(asteps)=float(dayni-day1900(string(yrn(0),'(i2)')+'0101')+1) yrsave(asteps)=fix(yrn(0)) endif else begin if forward gt 0 then begin ns=where(startarr le dayni) endif else begin ns=where(startarr ge dayni) endelse for nj=0L,n_elements(ns)-1 do begin tsave(ns(nj),asteps)=float(dayni-day1900(string(yrn(ns(nj)),'(i2)')+'0101')+1) yrsave(ns(nj),asteps)=fix(yrn(ns(nj))) endfor endelse latsave(0,asteps)=curlat lonsave(0,asteps)=curlon temsave(0,asteps)=tout psave(0,asteps)=pout qsave(0,asteps)=qout endif finished: endfor ; tsteps ;stop s=size(latsave) astepsx=s(2)-1 error1: ; compute solar zenith angle if n_elements(szsave) ge 0 then begin szsave=latsave*0. if n_elements(dif_start) le 1 then begin for jj=0,asteps do $ szsave(0,jj)=solar_zenith(latsave(*,jj),lonsave(*,jj),dnum=tsave(jj),year=yrsave(jj)) endif else begin for jj=0,asteps do begin t1=tsave(*,jj) lat1=latsave(*,jj) lon1=lonsave(*,jj) yr1=yrsave(*,jj) s1=lat1*0. n=where(t1 ge 0.) if n(0) ge 0 then $ s1(n)=solar_zenith(lat1(n),lon1(n),dnum=t1(n(0)),year=yr1(n(0))) szsave(0,jj)=s1 endfor endelse endif return end