pro mtraj,startday=startday,starttime=starttime,ndays=ndays,latin=latin,$ lonin=lonin,forward=forward,szsave=szsave,theta=theta, $ latsave=latsave,lonsave=lonsave,psave=psave,temsave=temsave,tsave=tsave,$ yrsave=yrsave,qsave=qsave,save_memory=save_memory,source=source,dir=dir,scale=scale0,$ latpole=latpole,lonpole=lonpole,segment=segment,noplot=noplot ;+ ; NAME: ; mtraj ; PURPOSE: ; Mini-trajectory Driver Program ; CATEGORY: ; trajectory ; CALLING SEQUENCE: ; ; INPUTS: only as keywords ; ; OPTIONAL INPUT PARAMETERS: ; none ; ; KEYWORD PARAMETERS: ; Mandatory keyword inputs ; startday (date for beginning of integration) string YYMMDD ; may be an array of dates, but dates shouls be close ; to each other ; starttime (real number) HH.mm where HH=hours (0-23) mm=minutes (0-59) ; may be an array of starting times [or] fraction of a day if the ; keyword /dayfrac is set ; ndays = number of days of integration, can be a fraction ; if multiple starttimes are used then ndays is the longest integration lengthr ; ; latin = vector of latitudes of points (negative for SH) ; lonin = vector of longitudes (same length as latin) ; theta = either a scalar or vector of the isentropic values for the points ; If a vector, then theta has same dimension as latin. ; ; Optional keyword inputs ; forward (+1 or -1), if lt 0 then back trajectory (default) ; delta_t the time step in days (default, 0.01) ; dayfrac = set this keyword if you set starttime in fractions of a day ; diabatic = switch to use diabatic cooling ; deltatheta = the theta "width" used to determine if new theta ; interpolation is needed in the diabatic case. The use ; of deltatheta greatly speeds up the code as new theta levels ; do not need to be generated for each small change in parcel theta ; save_memory = this variable indicates how often the data will be stored ; for example if save_memory=10 then every 10 time steps will be stored. ; Save once a day, = 50; default = 10 ; dtype = 'FX' or 'FIN' or 'FINR' default is 'FX' ; FX uses a 2km spaced theta grid for interpolation ; FIN uses a 0.5 km grid near the tropoause and a 1 km ; grid in the stratosphere and troposphere for finer resolution ; FINR used only for the NMC reanalysis which goes up to 10 mb ; other wise same as the FIN dataset ; dir_out = string which contains the path for the data interpolated onto ; potential surfaces. This is used if you wish to store test data sets ; overwrite= 0 use existing data files if found, =1 force an overwrite ; this is used if new heating rate data becomes available ; nodata_save = 0 use existing data files, =1 don't save files ; if nodata_save and overwrite are set then then new data files are computed but ; not saved. ; noheat_average =0 heating rate is globally averaged - else it is not ; asaved = number of days This is the auto-save feature. The history file is ; written out automatically every asaved days. Default is 10. Th history file is ; written as an idl save file command ; save,file='autosave.sav',tsave,latsave,lonsave,qsave,temsave,psave,yrsave ; ; ; Model control keywords ; source= data source character code, default 'NMC' ; create = call the parcel creation routine ftraj_newparts. ftraj_newparts calls a user ; defined routine: ; ftraj_np,latinx,loninx,theta,dayni,tsteps,index,lattemp,lontemp,theta_temp ; Input ; latinx=current parcel latitudes ; loninx=current parcel longitudes ; dayni = current day number ; tsteps = time step index ; Output ; index = long array which lists numerically the values of old parcels with a -1 ; where there are new parcels. ; lattemp = output array of latitudes (new +old parcels) ; lontemp = output array of longitudes (new +old parcels) ; theta_temp = output array of thetas (new +old parcels) ; Example: latinx = [20.,30.,40.], lontemp=[0.,0.,0.] ; The user wants to add parcels in between these points. The the output should be ; lattemp=[20.,25.,30.,35.,40.], lontemp=[0.,0.,0.,0.,0.], theta_temp= ... ; index= long([0,-1,1,-1,2]) ; Of course, if the user doesn't care about the order, then the following would also work ; lattemp=[20.,30.,40.,25.,35.], lontemp=[0.,0.,0.,0.,0.], theta_temp= ... ; index= long([0,1,2,-1,-1]) ; ** note that the addition of parcels will create a whole new history line for each parcel ; with values of zero until the parcel creation point. ; delete = delete parcels by calling the user defined routine ftraj_killparts with the call list: ; ftraj_killparts,latinx,loninx,theta,dayni,press,temp,pv,dead=dead ; Input ; the input values are as above (see create) plus parcel ; press=pressure ; temp=temperature ; pv=potential vorticity ; Output ; dead=parcel index array which has a value of -1 for parcels which are no longer ; to be used ; *** Note that the parcel is not actually eliminated - the history is saved, but the ; code no longer computes the parcel position. ; plot = number. Will call ftraj_plot when every 'plot' time steps. Ftraj_plot can also$ ; be a user defined ; routine. The call sequence is: ; ftraj_plot,curlat,curlon,curtheta,dayni,pout,tout,qout,dead=dead,u1=u1,$ ; t1=t1,dtype=dtype,scalar_val=scalar_val ; (see above for definition of the variables) ; ; File control keywords ; GRIDSPEC = (String) The code string of the gridding used. ; Default is no gridding. if ASM is the data default then ; Gridspec= 'GG2%5X2' ; special = string for certain met data sets - default = null ; for the NMC reanalysis set this to 'SREAN' ; seq = string for certain met data sets - default = null ; for the NMC reanalysis set this to 'E01' ; avespec = time averaging character code. ; Defaults to any character ; /xdr = create an xdr file that trajread can read. ; outfile = name of the output file ; ; Commonly used data source files ; UKMO data = source='UKM' mid 1991- present (beware of fall 1998) ; ASM data (DAO assimilation data) at 2x2.5 degrees resolution ; NMC data = source='NMC' mid 1979 - present (tropical winds are bogus) ; NMC Reanalysis = source='NMC,sequence='E01',special='SREAN' (no heating rates) 1954-1998 ; ; OUTPUTS: ; ; variables can be output through ; the call list as keywords: ; ; n = number of time steps, i = number of parcels = n_elements(latin) ; tsave(n) - time in days and fraction thereof in days ; if individual start times are used ; tsave(i,n) - for non relevant times tsave is set to -1 ; yrsave(n) - year values and fraction thereof in days ; if individual start times are used ; yrsave(i,n) - for non relevant times yrsave is set to 0 ; latsave(i,n) - latitudes in degrees ; lonsave(i,n) - longitude in degrees ; temsave(i,n) - temperature in K ; psave(i,n) - pressure in mb ; qsave(i,n) - PV in MKS units ; startarr(n) - starting daynumber for parcels ; if parcels are created, this variable has the creation date ; to convert the number into a date string use v=day1900(startarr,/string) ; scalar_val={a:r(n),b:l(n), ..} - this is a user assigned anonymous structure ; one common use is for scalar_val to be defined within ftraj_np ; ; Special meanings in these arrays ; if psave values are 1000. and temsave values are 300. then the parcel has been killed ; ; These variables are written out only if defined before the call ; szsave(i,n) - solar zenith angle in degrees ; The following variables are not written into the xdr files ; but can be output as a diagnostic through the call list ; hsave(i,n) - heating rate in K/day units ; usave(i,n) = zonal wind ; vsave(i,n) = meridional wind ; To trigger these options you must set szsave, usave, hsave and vsave to dummy variables in ; the call list e.g. ; usave=1 ; vsave=1 ; szsave=1 ; ftraj,....,usave=usave,vsave=vsave,szsave=szsave ; ; Multiple starting times ; In the case of multiple starting times or if parcels are injected the code keeps of ; when the parcels are started. This data is found is startarr(i) - which stores the ; day number (and day fraction) since Jan 1, 1900 (use date1900 to convert). ; The date/time that the code starts ; is stored in startdno. Thus the difference d=startday(i)-startdno is the time the ; parcel was injected after the code started to keep time. You will need this information ; since tsave stores the daynumber of the code and yrsave stores the year number. ; For parcels that haven't started running, the parcel history will either be a fixed ; latitude, longitude (equal to the starting value) until the parcel starts) or zer0 ; in the case of injected parcels. ; ; OPTIONAL OUTPUT PARAMETERS: ; COMMON BLOCKS: ; SIDE EFFECTS: ; ; Some Comments About Using the Fast Trajectory Model ; ; The fast trajectory model is built around interpolating data from ; fixed isentropic surfaces. The model uses precomputed arrays of the data ; on these fixed surfaces or ; will create and store them if they are not available. Heating rate ; arrays are automatically incorporated into the data arrays, but ; if no heating rates are available, these arrays are set to zero. ; ftraj will warn you if you use zero heating rate arrays when diabatic ; keyword is set. ; ; The model is designed to run fastest if a single isentropic ; surface is specified, but is also speedy (compared to trajdriver) ; if multiple surfaces or diabatic calcualtions are done. Timing ; estimates suggest that the model is 10-50 x faster than trajdriver ; in diabatic or multiple level isentropic mode. ; ; - Mark R. Schoeberl (schom@zephyr.gsfc.nasa.gov) ; ; RESTRICTIONS: ; PROCEDURE: ; This routine sets up the system for calls to ; ftrajmain ; ; REQUIRED ROUTINES: ; dayjul, day1900, date1900 ; MODIFICATION HISTORY: ; ; Trajdriver written by M. R. Schoeberl 7/3/89 ; Modified extensively by M. R. Schoeberl 7/3/91 ; ftraj routines created from trajdriver routines by M. Schoeberl 10/27/97 ; modified ftraj routines to correct for surface pressure MRS 8/31/99 ; Updated parcel injection/deletion routine, added plot routine MRS 10/4/99 ; 1999-11-29:nash:added avespec ; mtraj 1/6/00 ;- common ftraj_plot,firstcall,window1 firstcall=0 if n_elements(save_memory) eq 0 and n_elements(starttime) eq 1 then save_memory=10 if n_elements(latpole) eq 0 then latpole=90. if n_elements(lonpole) eq 0. then lonpole=0. if n_elements(starttime) gt 1 and n_elements(save_memory) eq 0 then save_memory=1 if n_elements(dir) eq 0 then dir='/app/ftp/pub/solve/met_data/mini_plan/' ; this is the location of the data files ; error messages if n_elements(theta) eq 0 then theta=480. if not (keyword_set(startday) and keyword_set(ndays) $ and keyword_set(latin) and keyword_set(lonin) ) then $ begin print,'Missing mandatory input variable(s):' if not keyword_set(startday) then print,' startday' if not keyword_set(starttime) then print,' starttime' if not keyword_set(ndays) then print,' ndays' if not keyword_set(latin) then print,' latin' if not keyword_set(lonin) then print,' lonin' print,'returning from mtraj' return endif ; check data consistency nlats=n_elements(latin) nlons=n_elements(lonin) if nlons ne nlats then begin print,'Number of latitude and longitude pts must be the same' print,' number of latitudes ',nlats print,' number of longitudes ',nlons err=-1 return endif ; float lat and lon arrays latin=float(latin) lonin=float(lonin) ; check theta s=size(theta) if n_elements(dt) eq 0 then dt=0.02 if n_elements(forward) eq 0 then forward=-1 ; set up time arrays sss=size(starttime) if sss(0) eq 0 then begin ; all the data start on the same day and time startdayx=replicate(startday,n_elements(latin)) if not keyword_set(dayfrac) then $ dayfract=(fix(starttime)+100.*(starttime-fix(starttime))/60.)/24. $ else dayfract=starttime start=double(day1900(startday))+double(dayfract) yr=fix(strmid(startday,0,2)) dayfractarr=fltarr(n_elements(latin))+dayfract startarr=double(fltarr(n_elements(latin)))+start yrarr=intarr(n_elements(latin))+yr endif else begin if sss(1) ne n_elements(latin) then begin stop,'array for starttime ne to latin size' endif dayfractarr=fltarr(n_elements(latin)) startarr=double(fltarr(n_elements(latin))) yrarr=intarr(n_elements(latin)) s2=size(startday) if s2(1) ne n_elements(latin) then begin startdayx=replicate(startday,n_elements(latin)) endif else startdayx=startday for jjs=long(0),sss(1)-1 do begin if not keyword_set(dayfrac) then $ dayfractarr(jjs)=(fix(starttime(jjs))+100.*(starttime(jjs)-fix(starttime(jjs)))/60.)/24. $ else dayfractarr=starttime startarr(jjs)=double(day1900(startdayx(jjs)))+double(dayfractarr(jjs)) yrarr(jjs)=fix(strmid(startdayx(jjs),0,2)) endfor endelse if not keyword_set(noprint) then print,' ***** GSFC Mini TRAJECTORY MODEL *****' if not keyword_set(noprint) then if forward lt 0 then print,'Back integration' else print,'Forward integration' if not keyword_set(noprint) then print,'Number of integration days requested = ',ndays if sss(0) eq 0 then begin if not keyword_set(noprint) then print,'Starting day, '+nicedate(startday) endif else begin if not keyword_set(noprint) then print,'Multiple starting times ranging from:' n=where( min(startarr) eq startarr) if not keyword_set(noprint) then if not keyword_set(noprint) then $ print,' '+nicedate(startdayx(n(0))),starttime(n(0)) if not keyword_set(noprint) then if not keyword_set(noprint) then print,'to:' n=where(max(startarr) eq startarr) if not keyword_set(noprint) then $ print,' '+nicedate(startdayx(n(0))),starttime(n(0)) endelse print,'Number of parcels = ',n_elements(latin) if save_memory ne 1 then if not keyword_set(noprint) then $ print,'Save Memory option ON, core reduction factor =',save_memory ; dt is set time step in fraction of a day print,'********' ;--------------------------------------------------------------- n=where(lonin lt 0.) if n(0) ge 0 then lonin(n)=lonin(n)+360. err=0 ; main trajectory computation program mtrajmain,startarr,yrarr,ndays,latin,lonin,dt,forward, $ latsave,lonsave,qsave,temsave,tsave,yrsave,szsave,psave,err, $ source=source,noprint=noprint,$ save_memory=save_memory,usave=usave,vsave=vsave,$ plot=plot,dir=dir,theta=theta,scale=scale0,latpole=latpole,$ segment=segment,lonpole=lonpole,noplot=noplot end