c @(#)trackav.F 1.1 11/1/94 program trackav c weighted averages of data along track c c < link trackav,'mlslib','epiclib','epslib','cdflib',sys$library:vaxcrtl/lib > c c written 1-JUL-1994 by mls c modified 14-JUL-1994 by mls to compute distance with great circle routine c from E.D.Cokelet c modified 16-AUG-1994 by mls; fixed variable distance intervals where the c target ending time is before the end of the file c #ifdef unix #include #else implicit none include 'ep_epslib:epsystem.fh' #endif c max number of points in depth integer mzdim parameter (mzdim=5000) c max number of points in time integer mscan parameter (mscan=50000) c max number of points total (z-t array) integer mpt parameter (mpt=100000) c max number of variables integer mv parameter (mv=20) c max number of time pairs integer mtarget parameter (mtarget = 500) c limit number of data files to read integer mfiles parameter (mfiles=1000) c max length of strings integer smax parameter (smax=255) character*(smax) epvx(mv) character*(smax) ccsname(mv),cclname(mv),ccgname(mv) character*(smax) ccvfrmt(mv),ccvunit(mv) character*(smax) czname(mv),czfrmt(mv),czunit(mv),cztype(mv) integer lepvx(mv),nv,pcnt(mv),ivar(mv),nvar,idvx(mv),nz(mv) integer bloc(mscan),nloc(mscan),nprofile(mscan),idz(mv) integer nds(mzdim,mv),midept(2,mscan) integer btrkept(2,mscan),etrkept(2,mscan) integer ntarget,btarget(2,mtarget),etarget(2,mtarget) real enssec(mscan),beglon(mscan),beglat(mscan) real midlon(mscan),midlat(mscan),endlon(mscan),endlat(mscan) real middist(mscan),ensdist(mscan),gapdist(mscan) real beglen(mscan),endlen(mscan),trklen(mscan),trkbeg(mscan) real btrklon(mscan),btrklat(mscan) real mtrklon(mscan),mtrklat(mscan) real etrklon(mscan),etrklat(mscan) real ztmp(mzdim),vtmpkp(mzdim,mv),zval(mzdim,mv) real uavg(mzdim,mv) real vtmp(mpt) real*8 uds(mzdim,mv),ds(mzdim,mv) logical vlabget(mv),zlabget(mv),havevx(mv) logical beginp(mscan) character*(smax) ptrinfile,datinfile character*(smax) ptroutfile,datoutfile,ptrname,outname character*(smax) csname,clname,cgname,cvfrmt,cvunit character*(smax) cxname,cxfrmt,cxunit character*(smax) cyname,cyfrmt,cyunit character*(smax) string,cfld(2),c1,c2 integer lfld(2),nfld,ldat,lname,l1,l2,llen integer icnt,nfiles,in,out,iscr,jscr,itab,istat integer ivtype,ixtype,iytype,ivfst,ivlst,iv integer i,k,l,n,nn,nk,idt,idx,idy,idxb,idyb,idxe,idye,nv_use integer epfin,epfout,lci(4),uci(4),origin(4),ifdim(4),imatch(4) integer ichoice,itmp,ierr,nptz,nptt,navg,lstavg,navg_in_file integer year,month,day,hour,minute,iavg,nprf,iprf,ilst,ntim integer meptkp(2),ept1(2),ept2(2),ept3(2),tdel(2),tval(2) real dtarget,rsec,lonmin,lonmax,ddbm,ddme,ddeb real blonkp,mlonkp,elonkp,blatkp,mlatkp,elatkp,eseckp real trkdist,trkaccm,avgend,rzspacing,trkmid real d1,d2,dx,t1,t2,tx,x1,x2,y1,y2 logical uniform_distance,ok,default,all_vars,lon180 logical lastfile,need_bloc,pinterp,zget,zsame logical west_positive,missing_data c functions real mytsec,maklone,makwest logical mytcmp,chkwest,chkept in = 5 out= 6 all_vars = .true. c explain write(out,10) 10 format(/,2x,'TRACKAV - averages interval data along track',//, 1 2x,'The averages are weighted by the distance', 1 ' represented by the profile',/,2x,'interval', 1 ' along the track. Input data files must' 1 ' contain start/mid/end',/,2x,'positions and', 1 ' times describing each profile interval.',/) 11 format(a) c log program usage call logger('trackav') c see if symbol was set by calling procedure llen = 0 call rtl_lib_get_symbol('EPIC_FILE',ptrinfile,llen,itab,istat) c prompt for pointer file name if (llen.le.0) call ask_ptrfilespec(in,out,ptrinfile) c set error options call eperropts(EPPRINTERR+EPCONTW) c open pointer file call epopen(epfin,ptrinfile,EPREAD) c prompt for type of average write(out,20) 20 format(2x,'Enter 1 for variable distance averages', 1 ' (select times to average between)',/, 1 8x,'[2] for uniform distance averages (select', 1 ' distance to average over)',/, 1 2x,'Enter choice > ',$) ichoice = 2 string = ' ' read(in,11) string if (string.ne.' ') then read(string,*,iostat=ierr) itmp if (ierr.eq.0) ichoice = itmp endif uniform_distance = (ichoice.ne.1) c prompt for average interval length as distance or time if (uniform_distance) then do 36 nn = 1, 10 write(out,32) 32 format(2x,'Enter distance to average (km) > ',$) string = ' ' read(in,11) string if (string.ne.' ') then read(string,*,iostat=ierr) dtarget if (ierr.eq.0) goto 38 endif 36 continue call stop_msg(out,'error reading input','error') 38 continue else write(out,42) 42 format(2x,'Enter 1 to read times from file',/, 1 2x,' [2] to type times from keyboard',/, 1 2x,'Enter choice > ',$) ichoice = 2 string = ' ' read(in,11) string if (string.ne.' ') then read(string,*,iostat=ierr) itmp if (ierr.eq.0) ichoice = itmp endif if (ichoice.eq.1) then c read times from file write(out,44) 44 format(2x,'Use an ASCII file with entries', 1 ' of the form DAY MONTH YEAR HOUR MINUTE SECOND',/, 1 2x,'example times: 1-JAN-84 10:00, 16 12 1992', 1 ' 23:43:01.00, 6 JUL 82 0600',/, 1 2x,'Each line should contain two times which', 1 ' data are to be averaged between.',/,2x,'The', 1 ' two times should be separated by a comma' 1 ' (e.g. time1, time2).',/,2x,'The lines in the file', 1 ' should be arranged chronologically.',/,2x, 1 'Enter filename > ',$) ok = .false. string = ' ' read(in,11) string call mystr_trim(smax,string,llen) if (llen.lt.1) goto 48 call open_oldlst(itmp,string(1:llen),ierr) if (ierr.ne.0) goto 48 ntarget = 0 do 46 nn = 1, mtarget string = ' ' 45 read(itmp,11,end=47,err=48) string if (string.eq.' ') goto 45 ok = .false. call mystr_parse(smax,smax,2,',',.false.,.true., 1 string,cfld,lfld,nfld) if (nfld.lt.2) goto 48 call mystr_timdmy2numeric(lfld(1),cfld(1),ok, 1 year,month,day,hour,minute,rsec) if (.not. ok) goto 48 call mdyhmstoeptime(month,day,year,hour,minute,rsec, 1 btarget(1,ntarget+1)) call mystr_timdmy2numeric(lfld(2),cfld(2),ok, 1 year,month,day,hour,minute,rsec) if (.not. ok) goto 48 call mdyhmstoeptime(month,day,year,hour,minute,rsec, 1 etarget(1,ntarget+1)) ntarget = ntarget + 1 46 continue 47 continue 48 if (.not. ok) call stop_msg(out, 1 'error reading time file','ERROR') close(itmp) else c prompt for times write(out,50) 50 format(2x,'Times are DAY MON YEAR [HHMM]',/, 1 2x,'Example: 3 Jul 86, 24 Dec 89 0600',/, 1 2x,'Example: 3 7 86, 24 12 89 0600',/, 1 2x,'Enter time1, time2 [or to terminate]',/) ntarget = 0 do 56 nn = 1, mtarget ok = .false. write(out,52) 52 format(2x,'Enter time1, time2 > ',$) string = ' ' read(in,11) string if (string.eq.' ') goto 58 call mystr_parse(smax,smax,2,',',.false.,.true., 1 string,cfld,lfld,nfld) if (nfld.lt.2) goto 54 call mystr_timdmy2numeric(lfld(1),cfld(1),ok, 1 year,month,day,hour,minute,rsec) if (.not. ok) goto 54 call mdyhmstoeptime(month,day,year,hour,minute,rsec, 1 btarget(1,ntarget+1)) call mystr_timdmy2numeric(lfld(2),cfld(2),ok, 1 year,month,day,hour,minute,rsec) if (.not. ok) goto 54 call mdyhmstoeptime(month,day,year,hour,minute,rsec, 1 etarget(1,ntarget+1)) ntarget = ntarget + 1 54 if (.not. ok) write(out,*) ' *** error reading input' 56 continue 58 continue endif write(out,60) ntarget 60 format(2x,'Number of time pairs read = ',i4) if (ntarget.lt.1) call stop_msg(out, 1 'error - no time pairs read','error') c show times read and check chronology do 64 nn = 1, ntarget call eptimetostr(btarget(1,nn), 1 'DD-MMM-YYYY hh:mm:ss.ff ',cfld(1),lfld(1)) call eptimetostr(etarget(1,nn), 1 'DD-MMM-YYYY hh:mm:ss.ff ',cfld(2),lfld(2)) write(out,62) nn,cfld(1),cfld(2) 62 format(2x,'Time pair ',i4,' > ',2a25) ok = .true. if (nn.gt.1) ok = 1 mytcmp(etarget(1,nn-1),'le',btarget(1,nn)) if (ok) ok = (mytcmp(btarget(1,nn),'lt',etarget(1,nn))) if (.not. ok) call stop_msg(out, 1 'error - times not chronological','error') 64 continue endif c prompt for variables to plot string = 'Process all variables in the first file?' llen = max(1,index(string,'?')) default = all_vars call ask_yn(in,out,string,llen,default,all_vars) if (.not. all_vars) then write(out,80) 80 format(2x,'Select the variable(s) for averaging') call ask_classicvars(in,out,smax,mv,' ',epvx,lepvx,nv) endif c c done prompting user c do 130 iv = 1, mv pcnt(iv) = 0 vlabget(iv) = .true. zlabget(iv) = .true. call szero(mzdim,ds(1,iv),uds(1,iv),nds(1,iv)) 130 continue do 132 n = 1, mscan nloc(n) = 0 nprofile(n) = 0 132 continue nfiles = 0 iprf = 0 nprf = 0 c open scratch files call open_scrunf(iscr,'scratchi.tmp') call open_scrunf(jscr,'scratchj.tmp') do 300 icnt = 1, mfiles if (eplastdatafile(epfin)) goto 302 c set to next data file call epsetnext(epfin,datinfile) call mystr_trim(smax,datinfile,llen) write(out,*)' reading file ',datinfile(:max(1,llen)) c if processing all vars then for first file get the variable codes of the data if (all_vars .and. nfiles.eq.0) then c get the id codes for the reserved variables if (.not. epfindvar(epfin,'ens_length',idt)) idt = 0 if (.not. epfindvar(epfin,'start_lon',idxb)) idxb = 0 if (.not. epfindvar(epfin,'start_lat',idyb)) idyb = 0 if (.not. epfindvar(epfin,'lon',idx)) idx = 0 if (.not. epfindvar(epfin,'lat',idy)) idy = 0 if (.not. epfindvar(epfin,'end_lon',idxe)) idxe = 0 if (.not. epfindvar(epfin,'end_lat',idye)) idye = 0 c get the variable list call epgetvarlist(epfin,ivar,mv,nvar) if (nvar.lt.1) call stop_msg(out, 1 'no variables in file','error') c select all the variables which aren't reserved variables nv = 0 do 180 i = 1, nvar itmp = ivar(i) if ( itmp.ne.idt .and. 1 itmp.ne.idxb .and. itmp.ne.idxe .and. 1 itmp.ne.idyb .and. itmp.ne.idye .and. 1 itmp.ne.idx .and. itmp.ne.idy ) then nv = nv + 1 write(epvx(nv),178) 'eps',itmp call mystr_collapse(10,epvx(nv),lepvx(nv)) 178 format(a,i6) endif 180 continue endif c check if ensemble length is in the file if (.not. epfindvar(epfin,'ens_length',idt)) call stop_msg 1 (out,'ens_length not found in file','error') c check if start and end positions are in the file ok = ( epfindvar(epfin,'start_lon',idxb) .and. 1 epfindvar(epfin,'start_lat',idyb) .and. 1 epfindvar(epfin,'end_lon',idxe) .and. 1 epfindvar(epfin,'end_lat',idye) ) if (.not. ok) call stop_msg(out, 1 'start/end positions not found in file','error') c check if mid positions are in the file as variables ok = ( epfindvar(epfin,'lon',idx) .and. 1 epfindvar(epfin,'lat',idy) ) if (.not. ok) call stop_msg(out, 1 'mid positions not found in file','error') c check variable(s) in file nv_use = 0 do 200 iv = 1, nv if (.not. epfindvar(epfin,epvx(iv),idvx(iv)) ) then write(out,192) epvx(iv)(:lepvx(iv)) 192 format(2x,'file is missing variable ',a, 1 ' -- skipping this variable') havevx(iv) = .false. else call epinqvar(epfin,idvx(iv),csname,clname,cvfrmt, 1 cvunit,ivtype) if (ivtype.ne.EPREAL) then write(out,196)epvx(iv)(:lepvx(iv)),idvx(iv),ivtype 196 format(2x,'variable ',a,' ( code = ',i6,' )', 1 ' is type ',i6,' -- skipping this variable') havevx(iv) = .false. else havevx(iv) = .true. nv_use = nv_use + 1 c get labels for output file if (vlabget(iv)) then cclname(iv) = clname ccsname(iv) = csname ccvfrmt(iv) = cvfrmt ccvunit(iv) = cvunit vlabget(iv) = .false. c get the generic name csname = ' ' cgname = ' ' call epinqcode(idvx(iv),csname,clname,cgname, 1 cvfrmt,cvunit) if (cgname.eq.' ') then ccgname(iv) = csname else ccgname(iv) = cgname endif endif endif endif 200 continue c skip file if no selected variables if (nv_use.lt.1) then write(out,214) 214 format(2x,'> no variables selected', 1 ' -- skipping this file') goto 300 endif c c here, we are using this file (unless not uniform_distance averages, c and the file is not in the time range specified) c nfiles = nfiles + 1 c get ivfst, index of first variable we want that is in the file do 216 iv = 1, nv if (havevx(iv)) then ivfst = iv goto 217 endif 216 continue 217 continue c check that variables are all on same time grid imatch(1) = 0 imatch(2) = 0 imatch(3) = 0 imatch(4) = 1 nn = ivfst+1 do 218 iv = nn, nv if (havevx(iv)) then if (.not. epcheckgrid(epfin,idvx(ivfst), 1 idvx(iv),imatch) ) call stop_msg(out, 1 'variables not on same time grid','error') endif 218 continue c check that variables on same time grid as position, ens length ok = (epcheckgrid(epfin,idvx(ivfst),idt,imatch) .and. 1 epcheckgrid(epfin,idvx(ivfst),idx,imatch) .and. 1 epcheckgrid(epfin,idvx(ivfst),idy,imatch) .and. 1 epcheckgrid(epfin,idvx(ivfst),idxb,imatch).and. 1 epcheckgrid(epfin,idvx(ivfst),idyb,imatch).and. 1 epcheckgrid(epfin,idvx(ivfst),idxe,imatch).and. 1 epcheckgrid(epfin,idvx(ivfst),idye,imatch) ) if (.not. ok) call stop_msg(out, 1 'vars not on same time grid as position or ens-len','error') c get the shape of the 1-d time variables; set to 1-d in x=y=z c time grid has been checked; all 1-d variables are on same time axis call epgetvarshap(epfin,idt,lci,uci) do 228 n = 1, 3 lci(n) = 1 uci(n) = 1 ifdim(n) = 1 228 continue ifdim(4) = (uci(4) - lci(4)) + 1 nptt = ifdim(4) if (nptt+2 .gt. mscan) call stop_msg(out, 1 'mscan size too small','ERROR') c read the time axis into midept call epgettaxis(epfin,idt,lci,uci,midept,nptt) c read ens_length into enssec call epinqvar(epfin,idt,csname,clname,cvfrmt,cvunit,ivtype) if (ivtype.ne.EPREAL) call stop_msg(out, 1 'error with ens_length data type','ERROR') call epgetvar(epfin,idt,lci,uci,enssec,ifdim) c read start positions into beglon,beglat call epinqvar(epfin,idxb,csname,cxname,cxfrmt,cxunit,ixtype) call epinqvar(epfin,idyb,csname,cyname,cyfrmt,cyunit,iytype) if (ixtype.ne.EPREAL.or.iytype.ne.EPREAL) call stop_msg 1 (out,'error with start lat or lon data type','ERROR') west_positive = chkwest(idxb,cxname,cxunit) call epgetvar(epfin,idxb,lci,uci,beglon,ifdim) call epgetvar(epfin,idyb,lci,uci,beglat,ifdim) do 232 n = 1, nptt beglon(n) = maklone(beglon(n),west_positive) 232 continue c read end positions into endlon,endlat call epinqvar(epfin,idxe,csname,cxname,cxfrmt,cxunit,ixtype) call epinqvar(epfin,idye,csname,cyname,cyfrmt,cyunit,iytype) if (ixtype.ne.EPREAL.or.iytype.ne.EPREAL) call stop_msg 1 (out,'error with end lat or lon data type','ERROR') west_positive = chkwest(idxe,cxname,cxunit) call epgetvar(epfin,idxe,lci,uci,endlon,ifdim) call epgetvar(epfin,idye,lci,uci,endlat,ifdim) do 236 n = 1, nptt endlon(n) = maklone(endlon(n),west_positive) 236 continue c read mid positions into midlon,midlat call epinqvar(epfin,idx,csname,cxname,cxfrmt,cxunit,ixtype) call epinqvar(epfin,idy,csname,cyname,cyfrmt,cyunit,iytype) if (ixtype.ne.EPREAL.or.iytype.ne.EPREAL) call stop_msg 1 (out,'error with lat or lon data type','ERROR') west_positive = chkwest(idx,cxname,cxunit) call epgetvar(epfin,idx,lci,uci,midlon,ifdim) call epgetvar(epfin,idy,lci,uci,midlat,ifdim) do 244 n = 1, nptt midlon(n) = maklone(midlon(n),west_positive) 244 continue c first time, check if longitudes should be fixed for crossing greenwich if (nfiles.eq.1) then call fix_lon(nptt,midlon,lonmin,lonmax) lon180 = (lonmin.lt.0) endif c change to +/- 180 if first file was changed if (lon180) then do 246 n = 1, nptt if (beglon(n).gt.180) beglon(n) = beglon(n) - 360.0 if (midlon(n).gt.180) midlon(n) = midlon(n) - 360.0 if (endlon(n).gt.180) endlon(n) = endlon(n) - 360.0 246 continue endif if (nfiles.eq.1) then c this is the first file, get the carry-over to use with next file blonkp = beglon(nptt) mlonkp = midlon(nptt) elonkp = endlon(nptt) blatkp = beglat(nptt) mlatkp = midlat(nptt) elatkp = endlat(nptt) eseckp = enssec(nptt) meptkp(1) = midept(1,nptt) meptkp(2) = midept(2,nptt) nptt = nptt - 1 else c add the last point from the previous file to the start of this file c and keep the last point in this file to use with next file call rshft(nptt,beglon,blonkp) call rshft(nptt,midlon,mlonkp) call rshft(nptt,endlon,elonkp) call rshft(nptt,beglat,blatkp) call rshft(nptt,midlat,mlatkp) call rshft(nptt,endlat,elatkp) call rshft(nptt,enssec,eseckp) call tshft(nptt,midept,meptkp) endif c check if this is the last data file lastfile = eplastdatafile(epfin) if (lastfile) then c last file, use the carry-over and pad with the end point of the last profile nptt = nptt + 1 beglon(nptt+1) = endlon(nptt) midlon(nptt+1) = endlon(nptt) endlon(nptt+1) = endlon(nptt) beglat(nptt+1) = endlat(nptt) midlat(nptt+1) = endlat(nptt) endlat(nptt+1) = endlat(nptt) enssec(nptt+1) = 0.0 call mytget(2,enssec(nptt),midept(1,nptt), 1 midept(1,nptt+1)) endif c check that at least one point to process if (nptt.lt.1) goto 251 c check chronology if (.not. chkept(midept,nptt) ) call stop_msg(out, 1 'error, times not chronological','error') c compute distances c mid dist = beg to mid c ens dist = beg to mid + mid to end c gap dist = end to start of next do 250 n = 1, nptt call grt_circle(beglon(n),beglat(n), 1 midlon(n),midlat(n),ddbm) call grt_circle(midlon(n),midlat(n), 1 endlon(n),endlat(n),ddme) call grt_circle(endlon(n),endlat(n), 1 beglon(n+1),beglat(n+1),ddeb) middist(n) = ddbm/1000.0 ensdist(n) = (ddbm + ddme)/1000.0 gapdist(n) = ddeb/1000.0 250 continue 251 continue c ======================================================================= c ======================================================================= c ======================================================================= c c find avg interval boundaries for this file c if (nfiles .eq. 1) then ntim = 1 navg = 1 lstavg = 0 trkdist = 0.0 trkaccm = 0.0 bloc(1) = 1 beglen(1) = 0.0 trkbeg(1) = 0.0 beginp(1) = .false. btrklon(1) = beglon(1) btrklat(1) = beglat(1) call mytget(1,enssec(1),midept(1,1),btrkept(1,1)) need_bloc = .true. else bloc(navg) = 0 endif navg_in_file = 0 c check that at least one point to process if (nptt.lt.1) goto 272 c ======================================================================= if (uniform_distance) then c ======================================================================= c uniform-distance averages do 258 n = 1, nptt nprf = nprf + 1 trkdist = trkdist + ensdist(n) + gapdist(n) nn = 0 c --------------------- c return to here at the end of the avg interval 252 continue c --------------------- nprofile(navg) = nprofile(navg) + 1 write(jscr) navg,nprofile(navg),nprf, 1 beglon(n),beglat(n), 1 midlon(n),midlat(n),endlon(n),endlat(n), 1 middist(n),ensdist(n),gapdist(n) if (trkdist.ge.dtarget) then navg_in_file = navg_in_file + 1 write(jscr) navg,nprofile(navg),nprf, 1 beglon(n+1),beglat(n+1), 1 midlon(n+1),midlat(n+1),endlon(n+1),endlat(n+1), 1 middist(n+1),ensdist(n+1),gapdist(n+1) c nloc is the current profile, where the avg interval boundary lies nloc(navg) = nprf bloc(navg+1) = nprf c avgend is the distance from the beginning of the current profile c to the avg interval boundary avgend = dtarget - trkaccm if (ensdist(n) .lt. avgend) then c avg interval ends within a gap, not within a profile beginp(navg+1) = .false. beglen(navg+1) = 0.0 endlen(navg) = ensdist(n) if (gapdist(n).gt.dtarget) then c its a big gap, don't interpolate pinterp = .false. avgend = ensdist(n) trkbeg(navg+1) = 0.0 else c its a small gap, setup to interpolate across the gap pinterp = .true. trkbeg(navg+1) = ensdist(n) + gapdist(n) - avgend d1 = 0 d2 = gapdist(n) dx = avgend - ensdist(n) x1 = endlon(n) x2 = beglon(n+1) y1 = endlat(n) y2 = beglat(n+1) call mytget(2,enssec(n),midept(1,n),ept1) call mytget(1,enssec(n+1),midept(1,n+1),ept2) call eptimesub(ept1,ept2,tval) t1 = 0 t2 = mytsec(tval) endif else c avg interval ends within a profile beginp(navg+1) = .true. beglen(navg+1) = ensdist(n) - avgend trkbeg(navg+1) = ensdist(n) + gapdist(n) - avgend endlen(navg) = avgend pinterp = .true. if (middist(n).gt.avgend) then c avg interval ends within first half of profile d1 = 0 d2 = middist(n) dx = avgend x1 = beglon(n) x2 = midlon(n) y1 = beglat(n) y2 = midlat(n) call mytget(1,enssec(n),midept(1,n),ept1) t1 = 0 t2 = enssec(n)/2. else c avg interval ends within last half of profile d1 = 0 d2 = ensdist(n) - middist(n) dx = avgend - middist(n) x1 = midlon(n) x2 = endlon(n) y1 = midlat(n) y2 = endlat(n) ept1(1) = midept(1,n) ept1(2) = midept(2,n) t1 = 0. t2 = enssec(n)/2. endif endif c adjust endlen if avg interval is entirely within the profile if (trkaccm.lt.0.0) endlen(navg)=endlen(navg)+trkaccm if (pinterp) then trklen(navg) = dtarget c interpolate to get position/time call lintrp(etrklon(navg),dx,x1,d1,x2,d2) call lintrp(etrklat(navg),dx,y1,d1,y2,d2) call lintrp(tx,dx,t1,d1,t2,d2) tdel(1) = tx/86400.0 tdel(2) = (tx - tdel(1)*86400 )*1000 call mytadd(ept1,tdel,etrkept(1,navg)) c next avg interval begins where this avg interval ends btrklon(navg+1) = etrklon(navg) btrklat(navg+1) = etrklat(navg) btrkept(1,navg+1) = etrkept(1,navg) btrkept(2,navg+1) = etrkept(2,navg) else trklen(navg) = trkaccm + ensdist(n) c assign position/time etrklon(navg) = endlon(n) etrklat(navg) = endlat(n) call mytget(2,enssec(n),midept(1,n),etrkept(1,navg)) c next avg interval begins with the next profile btrklon(navg+1) = beglon(n+1) btrklat(navg+1) = beglat(n+1) call mytget(1,enssec(n+1),midept(1,n+1), 1 btrkept(1,navg+1)) endif c move to next avg interval and reset the track distance accumulator navg = navg + 1 if (navg+2.ge.mscan) call stop_msg(out, 1 'error, mscan too small for navg','error') trkdist = trkbeg(navg) c reset beginp/beglen if the avg interval ends within the current profile/gap if (dtarget .le. trkdist) then nn = nn + 1 if (nn.gt.1000) call stop_msg(out, 1 'too many avgs within profile interval','error') beginp(navg) = .false. beglen(navg) = 0.0 trkbeg(navg) = 0.0 trkaccm = -avgend endif c loop back goto 252 endif trkaccm = trkdist 258 continue c ======================================================================= else c ======================================================================= c variable-distance averages (avg between input times btarget and etarget) do 268 n = 1, nptt nprf = nprf + 1 nn = 0 call mytget(1,enssec(n),midept(1,n),ept1) call mytget(2,enssec(n),midept(1,n),ept2) call mytget(1,enssec(n+1),midept(1,n+1),ept3) c --------------------- c return to here at the end of an avg interval to check c if the new avg interval is within this profile 262 continue c --------------------- missing_data = .false. if ( mytcmp(etarget(1,ntim),'lt',ept1)) then c the avg interval ends before this profile begins c move to the next avg interval ntim = ntim + 1 c pop out if all avg intervals have been defined if (ntim.gt.ntarget) goto 271 c make sure we are looking for the start of an average interval if (.not. need_bloc) call stop_msg(out, 1 'avg interval ends before start of profile','error') c loop back to check this new avg interval goto 262 elseif ( mytcmp(btarget(1,ntim),'gt',ept2) .and. 1 mytcmp(etarget(1,ntim),'lt',ept3)) then c the avg interval is in the gap between this profile and the next profile missing_data = .true. nprofile(navg) = nprofile(navg) + 1 write(jscr) navg,nprofile(navg),nprf, 1 beglon(n),beglat(n), 1 midlon(n),midlat(n),endlon(n),endlat(n), 1 middist(n),ensdist(n),gapdist(n) elseif (mytcmp(btarget(1,ntim),'le',ept2)) then c the avg interval begins before this profile ends, so, c this profile either begins the avg interval or is within the avg interval nprofile(navg) = nprofile(navg) + 1 write(jscr) navg,nprofile(navg),nprf, 1 beglon(n),beglat(n), 1 midlon(n),midlat(n),endlon(n),endlat(n), 1 middist(n),ensdist(n),gapdist(n) if (need_bloc) then c this profile begins the avg interval need_bloc = .false. bloc(navg) = nprf beginp(navg) = .true. if (mytcmp(btarget(1,ntim),'ge',midept(1,n))) then c avg interval begins within last half of profile pinterp = .true. t1 = 0 t2 = enssec(n)/2. call eptimesub(midept(1,n),btarget(1,ntim),tval) tx = mytsec(tval) x1 = midlon(n) x2 = endlon(n) y1 = midlat(n) y2 = endlat(n) d1 = middist(n) d2 = ensdist(n) elseif (mytcmp(btarget(1,ntim),'ge',ept1)) then c avg interval begins within first half of profile pinterp = .true. t1 = 0 t2 = enssec(n)/2. call eptimesub(ept1,btarget(1,ntim),tval) tx = mytsec(tval) x1 = beglon(n) x2 = midlon(n) y1 = beglat(n) y2 = midlat(n) d1 = 0.0 d2 = middist(n) else c avg interval begins before start of profile c use the start time of this profile as the start of the avg interval pinterp = .false. btrklon(navg) = beglon(n) btrklat(navg) = beglat(n) beglen(navg) = ensdist(n) btrkept(1,navg) = ept1(1) btrkept(2,navg) = ept1(2) endif if (pinterp) then c interpolate to get position/distance call lintrp(btrklon(navg),tx,x1,t1,x2,t2) call lintrp(btrklat(navg),tx,y1,t1,y2,t2) call lintrp(dx,tx,d1,t1,d2,t2) beglen(navg) = ensdist(n) - dx c target beginning time is start of avg interval btrkept(1,navg) = btarget(1,ntim) btrkept(2,navg) = btarget(2,ntim) endif trkbeg(navg) = beglen(navg) + gapdist(n) c reset the track distance accumulator trkaccm = 0.0 trkdist = beglen(navg) + gapdist(n) else c this profile is within the avg interval c accumulate the track distance trkdist = trkdist + ensdist(n) + gapdist(n) endif endif if (mytcmp(etarget(1,ntim),'lt',ept3)) then c the avg interval ends before the start of the next profile, so, c this profile is the last profile in the avg interval c finish this avg interval and move to the next avg interval nloc(navg) = nprf navg_in_file = navg_in_file + 1 write(jscr) navg,nprofile(navg),nprf, 1 beglon(n+1),beglat(n+1), 1 midlon(n+1),midlat(n+1),endlon(n+1),endlat(n+1), 1 middist(n+1),ensdist(n+1),gapdist(n+1) if (missing_data) then c there's no data for this avg interval pinterp = .false. beginp(navg) = .false. beglen(navg) = 0.0 endlen(navg) = 0.0 c make the avg interval the start/end of gap btrklon(navg) = endlon(n) btrklat(navg) = endlat(n) btrkept(1,navg) = ept2(1) btrkept(2,navg) = ept2(2) etrklon(navg) = beglon(n+1) etrklat(navg) = beglat(n+1) etrkept(1,navg) = ept3(1) etrkept(2,navg) = ept3(2) trklen(navg) = gapdist(n) trkbeg(navg) = gapdist(n) elseif (mytcmp(etarget(1,ntim),'gt',ept2)) then c avg interval ends after end of profile c use the end time of this profile as the end of the avg interval pinterp = .false. etrklon(navg) = endlon(n) etrklat(navg) = endlat(n) etrkept(1,navg) = ept2(1) etrkept(2,navg) = ept2(2) endlen(navg) = ensdist(n) trklen(navg) = trkaccm + endlen(navg) elseif (mytcmp(etarget(1,ntim),'ge',midept(1,n))) then c avg interval ends within last half of profile pinterp = .true. t1 = 0 t2 = enssec(n)/2. call eptimesub(midept(1,n),etarget(1,ntim),tval) tx = mytsec(tval) x1 = midlon(n) x2 = endlon(n) y1 = midlat(n) y2 = endlat(n) d1 = middist(n) d2 = ensdist(n) elseif (mytcmp(etarget(1,ntim),'ge',ept1)) then c avg interval ends within first half of profile pinterp = .true. t1 = 0 t2 = enssec(n)/2. call eptimesub(ept1,etarget(1,ntim),tval) tx = mytsec(tval) x1 = beglon(n) x2 = midlon(n) y1 = beglat(n) y2 = midlat(n) d1 = 0.0 d2 = middist(n) endif if (pinterp) then c interpolate to get position/distance call lintrp(etrklon(navg),tx,x1,t1,x2,t2) call lintrp(etrklat(navg),tx,y1,t1,y2,t2) call lintrp(endlen(navg),tx,d1,t1,d2,t2) trklen(navg) = trkaccm + endlen(navg) c target ending time is end of avg interval etrkept(1,navg) = etarget(1,ntim) etrkept(2,navg) = etarget(2,ntim) endif c adjust lengths if avg interval is entirely in this profile if ( (bloc(navg).eq.nloc(navg)) .and. 1 (.not. missing_data) ) then trklen(navg) = beglen(navg)+endlen(navg) - ensdist(n) endlen(navg) = trklen(navg) beglen(navg) = 0.0 beginp(navg) = .false. endif ntim = ntim + 1 c pop out if all avg intervals have been defined if (ntim.gt.ntarget) goto 271 c move to next avg interval and reset the track distance accumulator navg = navg + 1 if (navg+2.ge.mscan) call stop_msg(out, 1 'error, mscan too small for navg','error') trkdist = 0.0 need_bloc = .true. c loop back to check if new avg interval is within the current profile nn = nn + 1 if (nn.gt.1000) call stop_msg(out, 1 'too many avgs within profile interval','error') goto 262 endif trkaccm = trkdist 268 continue 271 continue c ======================================================================= endif c ======================================================================= 272 continue c this is the "last" file if we have all the time intervals defined if ( (.not. uniform_distance) .and. 1 (ntim.gt.ntarget) ) lastfile = .true. c terminate final interval on last data file if (lastfile) then if ( (.not.uniform_distance) .and. (need_bloc) .and. 1 (ntim.le.ntarget) ) then c variable distance, but never found beginning for current interval; discard it navg = navg - 1 elseif( (.not. uniform_distance) .and. 1 (ntim.gt.ntarget) ) then c variable distance, but we have all the intervals defined already continue else navg_in_file = navg_in_file + 1 nloc(navg) = nprf etrklon(navg) = endlon(nptt) etrklat(navg) = endlat(nptt) call mytget(2,enssec(nptt),midept(1,nptt),etrkept(1,navg)) trklen(navg) = trkaccm endlen(navg) = ensdist(nptt) c adjust endlen if avg interval is entirely within the profile if (uniform_distance .and. endlen(navg).gt.dtarget) then endlen(navg) = trkaccm beginp(navg) = .false. beglen(navg) = 0.0 trkbeg(navg) = 0.0 endif write(jscr) navg,nprofile(navg),nprf, 1 beglon(nptt+1),beglat(nptt+1), 1 midlon(nptt+1),midlat(nptt+1), 1 endlon(nptt+1),endlat(nptt+1), 1 middist(nptt+1),ensdist(nptt+1),gapdist(nptt+1) endif ntarget = navg need_bloc = .false. endif c ======================================================================= c ======================================================================= c ======================================================================= c c by here, the avg intervals are defined for this data file c read in vars, accumulate sums, avg data c do 290 iv = 1, nv c skip if the variable is not in the file if (.not. havevx(iv)) goto 290 c get the shape call epgetvarshap(epfin,idvx(iv),lci,uci) uci(1) = lci(1) uci(2) = lci(2) nptz = uci(3) - lci(3) + 1 nptt = uci(4) - lci(4) + 1 c check array sizes if (nptz*nptt .gt. mpt) call stop_msg(out, 1 'vtmp array size too small','ERROR') if (nptz .gt. mzdim) call stop_msg(out, 1 'mzdim size too small','ERROR') c set dimension ifdim(1) = 1 ifdim(2) = 1 ifdim(3) = nptz ifdim(4) = nptt c get the depth axis if (iv.eq.ivfst) then zget = .true. else imatch(1) = 0 imatch(2) = 0 imatch(3) = 1 imatch(4) = 0 if(epcheckgrid(epfin,idvx(iv),idvx(ivlst),imatch))then zget = .false. else zget = .true. endif endif if (zget) call epgetaxis(epfin,idvx(iv),3, 1 lci,uci,ztmp,mzdim) if (zlabget(iv)) then zlabget(iv) = .false. call epinqaxis(epfin,idvx(iv),3, 1 idz(iv),czname(iv),czfrmt(iv), 1 czunit(iv),cztype(iv),rzspacing) nz(iv) = nptz do 276 k = 1, nptz zval(k,iv) = ztmp(k) 276 continue else if (nptz.ne.nz(iv)) call stop_msg(out, 1 'var z-grid does not match','error') do 278 k = 1, nptz if (ztmp(k).ne.zval(k,iv)) call stop_msg(out, 1 'var z-grid does not match','error') 278 continue endif c get the variable call epgetvar(epfin,idvx(iv),lci,uci,vtmp,ifdim) c get the carry-over between files if (nfiles.eq.1) then nn = 1 + (nptt-1)*nptz do 282 k = 1, nptz vtmpkp(k,iv) = vtmp(nn) nn = nn + 1 282 continue nptt = nptt - 1 else call vshft(nptz,nptt,vtmp,vtmpkp(1,iv)) endif c use the carry-over if this is the last file if (lastfile) nptt = nptt + 1 c get the averages for this variable for this file iavg = lstavg + 1 if (nptt.lt.1) goto 289 do 288 l = 1, nptt if ( (.not. uniform_distance) .and. 1 (iavg.eq.navg) .and. (need_bloc) ) goto 289 nn = 1 + (l-1)*nptz nk = nn + nptz - 1 if (l+iprf.lt.bloc(iavg)) goto 288 c check if at avg interval boundary if (l+iprf.eq.nloc(iavg)) then c avg interval ended within this profile, add final contribution call myaccm(nptz,vtmp(nn), 1 endlen(iavg),ds(1,iv),uds(1,iv),nds(1,iv)) c compute avg do 284 k = 1, nptz if (nds(k,iv).eq.0) then uavg(k,iv) = 1.e35 elseif (ds(k,iv).le.1.d-5) then uavg(k,iv) = 0.0 else uavg(k,iv) = uds(k,iv)/ds(k,iv) endif 284 continue c write output write(iscr) iavg,idvx(iv),nptz, 1 (uavg(k,iv),k=1,nptz) call mypcnt(nptz,uavg(k,iv),pcnt(iv)) c move to next avg interval 286 continue iavg = iavg + 1 if (iavg.gt.navg) then if (.not. lastfile) call stop_msg(out, 1 'internal error-1 tracking navg','error') if (uniform_distance) then if (l.ne.nptt) call stop_msg(out, 1 'internal error-2 tracking navg','error') else if (iavg.le.ntarget) call stop_msg(out, 1 'internal error-3 tracking navg','error') c here, we've read the last profile for the final average (not necessarily c the last profile in the file); skip to process next variable goto 289 endif else c check if next avg interval is entirely within the current profile/gap c if so, there's no averaging to be done, just use the current profile if (l+iprf.eq.nloc(iavg)) then c but make sure we're not entirely in the gap if (endlen(iavg).le.0.0) then write(iscr) iavg,idvx(iv),nptz, 1 (1.e35,k=nn,nk) else write(iscr) iavg,idvx(iv),nptz, 1 (vtmp(k),k=nn,nk) call mypcnt(nptz,vtmp(nn),pcnt(iv)) endif goto 286 endif c zero ds uds nds arrays call szero(mzdim,ds(1,iv),uds(1,iv),nds(1,iv)) c check if next avg interval begins in this profile c and if so, add starting contribution if (l+iprf.eq.bloc(iavg) .and. beginp(iavg)) 1 call myaccm(nptz,vtmp(nn), 1 beglen(iavg),ds(1,iv),uds(1,iv),nds(1,iv)) endif elseif (l+iprf.eq.bloc(iavg)) then c if avg interval begins within this profile, add starting contribution if (beginp(iavg)) call myaccm(nptz,vtmp(nn), 1 beglen(iavg),ds(1,iv),uds(1,iv),nds(1,iv)) else c within an avg interval; accumulate sums call myaccm(nptz,vtmp(nn),ensdist(l), 1 ds(1,iv),uds(1,iv),nds(1,iv)) endif c loop back to process next profile 288 continue 289 continue ivlst = iv c loop back to process next variable 290 continue lstavg = lstavg + navg_in_file iprf = iprf + nptt if ((.not.uniform_distance).and.(lstavg.ge.ntarget)) goto 302 c loop back to process next file 300 continue c c by here, done reading input files c 302 continue if (nfiles.lt.1) call stop_msg(out,'no data files read','error') write(out,*)' done reading input data' navg = lstavg write(out,*)' number of avg intervals = ',navg if (navg.lt.1) call stop_msg(out,'no avgs computed','error') c get ivfst, index of first averaged variable that has some valid data ivfst = 0 do 306 iv = 1, nv if (pcnt(iv).gt.0) then ivfst = iv goto 307 endif 306 continue 307 continue if (ivfst.eq.0) call stop_msg(out, 1 'no valid data averaged','error') c c compute midept (mid time) c and enssec (length of interval in seconds) c from the beginning,ending times c do 310 nn = 1, navg call eptimesub(btrkept(1,nn),etrkept(1,nn),tval) call mythalf(tval,tdel) call mytadd(btrkept(1,nn),tdel,midept(1,nn)) enssec(nn) = mytsec(tval) 310 continue c c now get the positions at the mid point of the avg intervals c rewind(jscr) trkaccm = 0.0 ilst = 0 do 340 nn = 1, navg c read in positions, distances for profiles bracketing the avg interval nprf = nprofile(nn) + 1 do 320 n = 1, nprf read(jscr) iavg,iprf,itmp,beglon(n),beglat(n), 1 midlon(n),midlat(n),endlon(n),endlat(n), 1 middist(n),ensdist(n),gapdist(n) if (iavg.ne.nn) call stop_msg(out, 1 'internal error-1 getting mid-positions','error') 320 continue nprf = nprofile(nn) if (iprf.ne.nprf) call stop_msg(out, 1 'internal error-2 getting mid-positions','error') c reset track accumulator if this is a new profile if (itmp.ne.ilst) trkaccm = 0.0 ilst = itmp c target mid distance is half the track length for the avg interval trkmid = trklen(nn)/2.0 c adjust trkmid relative to start of profile if (uniform_distance .and. nprf.eq.1) then trkmid = trkmid + trkaccm else trkmid = trkmid + ensdist(1) + gapdist(1) - trkbeg(nn) endif trkdist = 0 trkaccm = trkaccm + endlen(nn) pinterp = .false. do 330 n = 1, nprf if (trkdist + middist(n) .ge. trkmid) then pinterp = .true. d1 = trkdist d2 = trkdist + middist(n) dx = trkmid x1 = beglon(n) x2 = midlon(n) y1 = beglat(n) y2 = midlat(n) elseif (trkdist + ensdist(n) .ge. trkmid) then pinterp = .true. d1 = trkdist + middist(n) d2 = trkdist + ensdist(n) dx = trkmid x1 = midlon(n) x2 = endlon(n) y1 = midlat(n) y2 = endlat(n) elseif (trkdist + ensdist(n)+gapdist(n) .ge. trkmid) then pinterp = .true. d1 = trkdist + ensdist(n) d2 = trkdist + ensdist(n) + gapdist(n) dx = trkmid x1 = endlon(n) x2 = beglon(n+1) y1 = endlat(n) y2 = beglat(n+1) endif if (pinterp) then call lintrp(mtrklon(nn),dx,x1,d1,x2,d2) call lintrp(mtrklat(nn),dx,y1,d1,y2,d2) goto 338 endif trkdist = trkdist + ensdist(n) + gapdist(n) 330 continue call stop_msg(out,'error getting avg mid lat/lon','error') 338 continue 340 continue close(jscr) c c by here, we have the avg interval beginning, middle, ending positions/times c c change to +west for longitudes do 350 n = 1, navg btrklon(n) = makwest(btrklon(n),lon180) mtrklon(n) = makwest(mtrklon(n),lon180) etrklon(n) = makwest(etrklon(n),lon180) 350 continue ptrname = ptrinfile call make_name(ptrname,llen,outname,lname,c1,l1,c2,l2) c open the output pointer file ptroutfile = outname(:lname)//'_trackav.dat' call mystr_trim(smax,ptroutfile,llen) call epopen(epfout,ptroutfile,EPCREATE+EPCDF) datoutfile = outname(:lname)//'.cdf' call mystr_trim(smax,datoutfile,ldat) call epsetnext(epfout,datoutfile) write(out,*)' writing output file ',datoutfile(:ldat) c copy attributes from last input file to the output file call epcopyatts(epfin,epfout) c close input pointer file call epclose(epfin) c write global attribute if (uniform_distance) then write(string,362) dtarget 362 format('Target distance = ',g15.7,' km') else write(string,364) 364 format('Averaging along track between input times',a) endif call mystr_compress(smax,string,llen) call epputattvalc(epfout,'TRACKAV_COMMENT',llen,string) c write out the 1-d variables do 380 n = 1, 4 lci(n) = 1 uci(n) = 1 origin(n) = 1 ifdim(n) = 1 380 continue uci(4) = navg ifdim(4) = navg c ensemble length idt = 1200 call epsetvar(epfout,idt,' ',' ',' ',' ',' ',EPREAL) c dummy x axis idx = 501 cvfrmt = ' ' cvunit = ' ' call epinqcode(idx,csname,clname,cgname,cvfrmt,cvunit) call epsetaxis(epfout,idt,1,0,'x-start',cvfrmt,cvunit,EPEVEN) call epputaxis(epfout,idt,1,origin,lci,uci,mtrklon) c dummy y axis idy = 500 cvfrmt = ' ' cvunit = ' ' call epinqcode(idy,csname,clname,cgname,cvfrmt,cvunit) call epsetaxis(epfout,idt,2,0,'y-start',cvfrmt,cvunit,EPEVEN) call epputaxis(epfout,idt,2,origin,lci,uci,mtrklat) c dummy z axis call epsetaxis(epfout,idt,3,0,'z-start', 1 czfrmt(ivfst),czunit(ivfst),EPEVEN) call epputaxis(epfout,idt,3,origin,lci,uci,zval(1,ivfst)) c actual t axis call epsettaxis(epfout,idt,'t-axis',' ',' ',EPUNEVEN) call epputtaxis(epfout,idt,origin,lci,uci,midept) c put data call epputvar(epfout,idt,origin,lci,uci,enssec,ifdim) c start longitude idxb = 1501 call epsetvar(epfout,idxb,' ',' ',' ',' ',' ',EPREAL) c copy grid call epcopygrid(epfout,idt,epfout,idxb) c put data call epputvar(epfout,idxb,origin,lci,uci,btrklon,ifdim) c start latitude idyb = 1500 call epsetvar(epfout,idyb,' ',' ',' ',' ',' ',EPREAL) c copy grid call epcopygrid(epfout,idt,epfout,idyb) c put data call epputvar(epfout,idyb,origin,lci,uci,btrklat,ifdim) c mid longitude idx = 501 call epsetvar(epfout,idx,' ',' ',' ',' ',' ',EPREAL) c copy grid call epcopygrid(epfout,idt,epfout,idx) c put data call epputvar(epfout,idx,origin,lci,uci,mtrklon,ifdim) c mid latitude idy = 500 call epsetvar(epfout,idy,' ',' ',' ',' ',' ',EPREAL) c copy grid call epcopygrid(epfout,idt,epfout,idy) c put data call epputvar(epfout,idy,origin,lci,uci,mtrklat,ifdim) c end longitude idxe = 1511 call epsetvar(epfout,idxe,' ',' ',' ',' ',' ',EPREAL) c copy grid call epcopygrid(epfout,idt,epfout,idxe) c put data call epputvar(epfout,idxe,origin,lci,uci,etrklon,ifdim) c end latitude idye = 1510 call epsetvar(epfout,idye,' ',' ',' ',' ',' ',EPREAL) c copy grid call epcopygrid(epfout,idt,epfout,idye) c put data call epputvar(epfout,idye,origin,lci,uci,etrklat,ifdim) c check if variables are on the same z-grid zsame = .false. nptz = nz(ivfst) nn = ivfst + 1 do 394 iv = nn, nv if (pcnt(iv).gt.0) then if (nz(iv).ne.nz(ivfst)) goto 396 do 392 k = 1, nptz if (zval(k,iv).ne.zval(k,ivfst)) goto 396 392 continue endif 394 continue zsame = .true. 396 continue c setup array of bad-val do 398 k = 1, mzdim vtmpkp(k,1) = 1.e35 398 continue c loop to write each variable to the output file ivlst = 0 do 480 iv = 1, nv c print msg and skip if no good data read for this variable if (vlabget(iv) .or. zlabget(iv) .or. pcnt(iv).lt.1) then write(out,408) epvx(iv)(:lepvx(iv)) 408 format(2x,'> no valid data read for ',a,' ***') goto 480 endif rewind(iscr) c set the name call epsetvar(epfout,idvx(iv),ccsname(iv),cclname(iv), 1 ccgname(iv),ccvfrmt(iv),ccvunit(iv),EPREAL) c set the z dimension uci(3) = nz(iv) ifdim(3) = nz(iv) c copy grid if all vars were on the same z-grid if (zsame .and. ivlst.ne.0) then call epcopygrid(epfout,idvx(ivlst),epfout,idvx(iv)) else c copy the x-y-t grid from ensemble_length call epcopygrid(epfout,idt,epfout,idvx(iv)) c set the z grid call epsetaxis(epfout,idvx(iv),3,idz(iv),czname(iv), 1 czfrmt(iv),czunit(iv),cztype(iv)) call epputaxis(epfout,idvx(iv),3,origin,lci,uci, 1 zval(1,iv)) endif c read data off scratch unit into vtmp iavg = 0 do 420 n = 1, navg if (iavg.lt.n) then 418 read(iscr) iavg,itmp,nptz,(ztmp(k),k=1,nptz) if (itmp.ne.idvx(iv)) goto 418 if (nptz.ne.nz(iv)) call stop_msg(out, 1 'internal error-1 reading avgs','error') if (iavg.lt.n) call stop_msg(out, 1 'internal error-2 reading avgs','error') endif if (iavg.eq.n) then call myfill(nptz,n,ztmp,vtmp) else c here, data are missing (not in input file) for this avg interval for this var call myfill(nptz,n,vtmpkp,vtmp) endif 420 continue c put data call epputvar(epfout,idvx(iv),origin,lci,uci,vtmp,ifdim) ivlst = iv c loop back to process next variable 480 continue c close scratch unit close(iscr) c close output pointer file call epclose(epfout) call stop_msg(out,'done writing data','ok') end c ----------------------------------------------------------------------- subroutine mytadd(time1,delta,time2) c ----------------------------------------------------------------------- c like eptimeadd, but handles negative delta time implicit none integer time1(2),delta(2),time2(2) integer mydel(2),rem,msecperday parameter (msecperday=86400000) rem = delta(2)/msecperday mydel(1) = delta(1) + rem mydel(2) = delta(2) - rem*msecperday time2(1) = time1(1) + mydel(1) time2(2) = time1(2) + mydel(2) if (time2(2) .ge. msecperday) then time2(1) = time2(1) + 1 time2(2) = time2(2) - msecperday endif if (time2(2) .lt. 0) then time2(1) = time2(1) - 1 time2(2) = time2(2) + msecperday endif return end c ----------------------------------------------------------------------- subroutine mythalf(timdiff,timhalf) c ----------------------------------------------------------------------- c divides timdiff by 2 implicit none integer timdiff(2),timhalf(2) timhalf(1) = timdiff(1)/2 timhalf(2) = timdiff(2)/2 if (mod(timdiff(1),2).ne.0) timhalf(2) = timhalf(2)+43200000 if (timhalf(2).gt.86400000) then timhalf(2) = timhalf(2) - 86400000 timhalf(1) = timhalf(1) + 1 endif return end c ----------------------------------------------------------------------- subroutine mytget(iflag,rsec,tmid,tout) c ----------------------------------------------------------------------- c iflag = 1 get tout = start of interval c iflag = 2 get tout = end of interval c given mid time (tmid) and length of interval in seconds (rsec) implicit none integer iflag,tmid(2),tout(2) real rsec,rval integer tval(2) if (rsec.eq.0.0) then tout(1) = tmid(1) tout(2) = tmid(2) return endif rval = abs(rsec)/2.0 if (iflag.eq.1) rval = -rval tval(1) = rval / 86400.0 tval(2) = (rval - float(tval(1))*86400.0) * 1000.0 call mytadd(tmid,tval,tout) return end c ----------------------------------------------------------------------- real function mytsec(tval) c ----------------------------------------------------------------------- c changes eps delta time to seconds implicit none integer tval(2) mytsec = float(tval(1))*86400.0 + float(tval(2))/1000.0 return end c ----------------------------------------------------------------------- logical function mytcmp(t1,cmp,t2) c ----------------------------------------------------------------------- c compares eps times t1 and t2 c cmp is how to compare, e.g. 'gt' 'ge' 'eq' 'le' 'lt' implicit none character*(*) cmp integer t1(2),t2(2) logical less,equal,great character*4 str str = cmp call mystr_lowercase(str,4) equal = ( t1(1).eq.t2(1) .and. t1(2).eq.t2(2) ) less = ( (t1(1).lt.t2(1)) .or. 1 (t1(1).eq.t2(1) .and. t1(2).lt.t2(2)) ) great = (.not. less .and. .not. equal) if ( (equal .and. (index(str,'e').ne.0) ) .or. 1 (less .and. (index(str,'l').ne.0) ) .or. 1 (great .and. (index(str,'g').ne.0) ) ) then mytcmp = .true. else mytcmp = .false. endif return end c ----------------------------------------------------------------------- logical function chkept(ept,nt) c ----------------------------------------------------------------------- c checks chronology of eps time array c returns true if chronological implicit none integer ept(2,*),nt integer n,n1 n1 = 1 do 100 n = 2, nt if ( (ept(1,n).lt.ept(1,n1)) .or. 1 (ept(1,n).eq.ept(1,n1) .and. 1 ept(2,n).lt.ept(2,n1)) ) then chkept = .false. return endif n1 = n 100 continue chkept = .true. return end c ----------------------------------------------------------------------- real function maklone(xval,west_positive) c ----------------------------------------------------------------------- c changes to +e longitude 0-360 implicit none real xval,xtmp logical west_positive xtmp = xval if (xtmp.lt.0) xtmp = xtmp + 360.0 if (west_positive) xtmp = 360.0 - xtmp if (xtmp.ge.360.0) xtmp = xtmp - 360.0 maklone = xtmp return end c ----------------------------------------------------------------------- real function makwest(xval,do180) c ----------------------------------------------------------------------- c changes to +w longitude implicit none real xval,xtmp logical do180 xtmp = xval if (xtmp.lt.0) xtmp = xtmp + 360.0 xtmp = 360.0 - xtmp if (do180 .and. xtmp.gt.180) xtmp = xtmp - 360.0 if (xtmp.ge.360.0) xtmp = xtmp - 360.0 makwest = xtmp return end c ----------------------------------------------------------------------- subroutine lintrp(xx,dx,x1,d1,x2,d2) c ----------------------------------------------------------------------- c linearly interpolates to get xx at dx given x1,d1 x2,d2 implicit none real xx,dx,x1,d1,x2,d2 if (abs(d2-d1) .le. 1.e-5) then xx = x1 else xx = x1 + ( (dx-d1)/(d2-d1) )*(x2-x1) endif return end c ----------------------------------------------------------------------- subroutine szero(n,ds,uds,nds) c ----------------------------------------------------------------------- c zeros out arrays ds,uds,nds implicit none integer i,n,nds(*) real*8 ds(*),uds(*) do 100 i = 1, n ds(i) = 0.d0 uds(i) = 0.d0 nds(i) = 0 100 continue return end c ----------------------------------------------------------------------- subroutine myaccm(nz,v,d,ds,uds,nds) c ----------------------------------------------------------------------- c accumulates ds, uds, and nds given nz, v() and d implicit none integer k,nz,nds(*) real v(*),d real*8 ds(*),uds(*) do 100 k = 1, nz if (v(k).lt.1.e35) then ds(k) = ds(k) + d uds(k) = uds(k) + v(k)*d nds(k) = nds(k) + 1 endif 100 continue return end c ----------------------------------------------------------------------- subroutine rshft(n,x,xkp) c ----------------------------------------------------------------------- c inserts input xkp as first element of x (shifts elements in x() down by 1) c returns new xkp as last element of x implicit none integer i,n real x(*),xkp,xnxt,xlst xlst = xkp do 100 i = 1, n xnxt = x(i) x(i) = xlst xlst = xnxt 100 continue x(n+1) = xlst xkp = x(n+1) return end c ----------------------------------------------------------------------- subroutine vshft(nz,nt,v,vkp) c ----------------------------------------------------------------------- c inserts input profile vkp as first profile in v array c (moves profiles in v array over by 1 element) c returns new vkp as last profile of v array implicit none integer i,k,nz,nt real v(nz,*),vkp(*) real xnxt,xlst do 200 k = 1, nz xlst = vkp(k) do 100 i = 1, nt xnxt = v(k,i) v(k,i) = xlst xlst = xnxt 100 continue v(k,nt+1) = xlst vkp(k) = v(k,nt+1) 200 continue return end c ----------------------------------------------------------------------- subroutine tshft(n,t,tkp) c ----------------------------------------------------------------------- c inserts input tkp as first element of t (shifts elements in t() down by 1) c returns new tkp as last element of t implicit none integer i,n,t(2,*),tkp(2),tnxt(2),tlst(2) tlst(1) = tkp(1) tlst(2) = tkp(2) do 100 i = 1, n tnxt(1) = t(1,i) tnxt(2) = t(2,i) t(1,i) = tlst(1) t(2,i) = tlst(2) tlst(1) = tnxt(1) tlst(2) = tnxt(2) 100 continue t(1,n+1) = tlst(1) t(2,n+1) = tlst(2) tkp(1) = t(1,n+1) tkp(2) = t(2,n+1) return end c ----------------------------------------------------------------------- subroutine mypcnt(nz,v,pcnt) c ----------------------------------------------------------------------- c increments pcnt if any element of v is valid (< 1.e35) implicit none integer nz,pcnt,k real v(*) do 100 k = 1, nz if (v(k).lt.1.e35) then pcnt = pcnt + 1 return endif 100 continue return end c ----------------------------------------------------------------------- subroutine myfill(nz,n,z,v) c ----------------------------------------------------------------------- c puts z-profile into v array as the profile n implicit none integer nz,n,nn,k real z(*), v(*) nn = 1 + (n-1)*nz do 100 k = 1, nz v(nn) = z(k) nn = nn + 1 100 continue return end