WebSVN – lagranto.ecmwf – Blame – /trunk/select/select.f

Rev	Author	Line No.	Line
3	michaesp	1	`PROGRAM select`
		2
		3	`c **************************************************************`
		4	`c * Select trajectories from LSL file *`
		5	`c * Michael Sprenger / January, February 2008 *`
		6	`c **************************************************************`
		7
		8	`implicit none`
		9
		10	`c --------------------------------------------------------------`
		11	`c Declaration of parameters`
		12	`c --------------------------------------------------------------`
		13
		14	`c Maximum number of columns per trajectory`
		15	`integer maxcol`
		16	`parameter (maxcol=100)`
		17
		18	`c Maximum number of commands`
		19	`integer maxcmd`
		20	`parameter (maxcmd=10000)`
		21
19	michaesp	22	`c -------------------------------------------------------------`
3	michaesp	23	`c Declaration of variables`
		24	`c --------------------------------------------------------------`
		25
		26	`c Input and output files`
		27	`character*120 inp_lslfile ! Input lsl file`
		28	`character*120 out_lslfile ! Output lsl file`
		29	`character*120 inp_criteria ! Input file with criteria`
		30	`integer inpmode ! Format of input file`
		31	`integer outmode ! Format of output file`
		32	`character*80 outformat ! Trajectory/Boolean/Index of output`
		33	`character*80 regionf ! Name of the regionfile`
		34
		35	`c Trajectories`
		36	`integer ntim ! Number of times`
		37	`integer ncol ! Number of columns (including time...)`
		38	`integer ntrainp,ntraout ! Number of trajectories`
		39	`character*80 vars(maxcol) ! Names of trajectory columns`
		40	`integer basetime(6) ! Base time of trajectory (first line in lsl)`
		41	`real,allocatable, dimension (:,:,:) :: trainp ! Input trajectories`
		42	`real,allocatable, dimension (:,:,:) :: traout ! Output trajectories`
		43	`integer,allocatable,dimension (:) :: selflag ! Flag for selection`
		44	`real,allocatable, dimension (:) :: time ! Times of the trajectory`
		45	`integer,allocatable,dimension (:,:) :: trigger ! Trigger column`
		46	`integer itrigger ! Column index for trigger`
		47	`character*80 addtrigger ! Flag whether to add trigger column`
		48
		49	`c Command stack`
		50	`real cmd(maxcmd) ! Decoded slection criterion`
		51	`integer ncmd ! Number of commands`
		52
		53	`c Common block for initialisation of polygon check`
		54	`real tlonv(2000),vlat_c(2000),vlon_c(2000)`
		55	`real xlat_c,xlon_c`
		56	`integer ibndry,nv_c`
		57	`data ibndry/0/`
		58	`common /spolybndry/vlat_c,vlon_c,nv_c,xlat_c,xlon_c,tlonv,ibndry`
		59
		60	`c Auxiliary variables`
		61	`integer stat ! Logical (state) variable`
		62	`integer fid,fod,fcr ! File identifier for input and output`
		63	`integer i,j,k ! Index counter`
		64	`integer isok ! Flag for selected trajectory`
		65	`real param(1000) ! List of parameters`
		66	`integer nparam ! Number of parameters`
		67	`character*80 specialstr ! Name of special command`
		68	`integer len`
		69	`integer,allocatable,dimension (:) :: trigger1 ! Trigger column`
		70	`real,allocatable,dimension (:,:) :: trainp1 ! A single trajectory`
		71	`character ch`
		72
		73	`c --------------------------------------------------------------`
		74	`c Preparations`
		75	`c --------------------------------------------------------------`
		76
		77	`c Write start message`
		78	`print*,'========================================================='`
		79	`print,' START OF PROGRAM SELECT *'`
		80	`print*`
		81
		82	`c Read parameter file`
		83	`open(10,file='select.param')`
		84	`read(10,*) inp_lslfile`
		85	`read(10,*) out_lslfile`
		86	`read(10,*) outformat`
		87	`read(10,*) inp_criteria`
		88	`read(10,*) ntrainp`
		89	`read(10,*) ntim`
		90	`read(10,*) ncol`
		91	`read(10,*) regionf`
		92	`read(10,*) addtrigger`
		93	`close(10)`
		94
		95	`c Set the formats of the input and output files`
		96	`call mode_tra(inpmode,inp_lslfile)`
		97	`if (inpmode.eq.-1) inpmode=1`
		98	`call mode_tra(outmode,out_lslfile)`
		99	`if ( (outmode.eq.-1).and.(outformat.ne.'startf') ) then`
		100	`outmode=1`
		101	`endif`
		102
		103	`c Allocate memory for a single trajectory`
		104	`allocate(trainp(ntrainp,ntim,ncol),stat=stat)`
		105	`if (stat.ne.0) stop '* error allocating array trainp *'`
		106	`allocate(time(ntim),stat=stat)`
		107	`if (stat.ne.0) stop '* error allocating array time *'`
		108	`allocate(selflag(ntrainp),stat=stat)`
		109	`if (stat.ne.0) stop '* error allocating array selflag *'`
		110	`allocate(trigger(ntrainp,ntim),stat=stat)`
		111	`if (stat.ne.0) stop '* error allocating array trigger *'`
		112	`allocate(trigger1(ntim),stat=stat)`
		113	`if (stat.ne.0) stop '* error allocating array trigger1 *'`
		114	`allocate(trainp1(ntim,ncol),stat=stat)`
		115	`if (stat.ne.0) stop '* error allocating array trainp1 *'`
		116
		117	`c Read the input trajectory file`
		118	`call ropen_tra(fid,inp_lslfile,ntrainp,ntim,ncol,`
		119	`> basetime,vars,inpmode)`
		120	`call read_tra (fid,trainp,ntrainp,ntim,ncol,inpmode)`
		121	`call close_tra(fid,inpmode)`
		122
		123	`c Check that first four columns correspond to time,lon,lat,p`
		124	`if ( (vars(1).ne.'time' ).or.`
		125	`> (vars(2).ne.'xpos' ).and.(vars(2).ne.'lon' ).or.`
		126	`> (vars(3).ne.'ypos' ).and.(vars(3).ne.'lat' ).or.`
		127	`> (vars(4).ne.'ppos' ).and.(vars(4).ne.'p' ) )`
		128	`>then`
		129	`print*,' ERROR: problem with input trajectories ...'`
		130	`stop`
		131	`endif`
		132	`vars(1) = 'time'`
		133	`vars(2) = 'lon'`
		134	`vars(3) = 'lat'`
		135	`vars(4) = 'p'`
		136
		137	`c Get the trajectory times from first trajectory`
		138	`do i=1,ntim`
		139	`time(i)=trainp(1,i,1)`
		140	`enddo`
		141
		142	`c Init the trigger field - first check whether it is already available`
		143	`itrigger = 0`
		144	`do i=1,ncol`
		145	`if ( vars(i).eq.'TRIGGER' ) itrigger = i`
		146	`enddo`
		147
		148	`if ( itrigger.ne.0 ) then`
		149	`do i=1,ntrainp`
		150	`do j=1,ntim`
		151	`trigger(i,j) = nint( trainp(i,j,itrigger) )`
		152	`enddo`
		153	`enddo`
		154	`else`
		155	`do i=1,ntrainp`
		156	`do j=1,ntim`
		157	`trigger(i,j) = 0`
		158	`enddo`
		159	`enddo`
		160	`endif`
		161
		162	`c Write some info about the trajectory`
		163	`print*,'---- INPUT PARAMETERS -----------------------------------'`
		164	`write(,)`
		165	`write(,) 'Input file : ',trim(inp_lslfile)`
		166	`write(,) 'Output file : ',trim(out_lslfile)`
		167	`write(,) 'Output format : ',trim(outformat)`
		168	`write(,) 'Criteria file : ',trim(inp_criteria)`
		169	`write(,) '# tra : ',ntrainp`
		170	`write(,) '# time : ',ntim`
		171	`write(,) '# col : ',ncol`
		172	`write(,) 'Region file : ',trim(regionf)`
		173	`write(,) 'Add trigger : ',trim(addtrigger)`
		174
		175	`print*`
		176	`print*,'---- INPUT TRAJECTORY FILE ------------------------------'`
		177	`print*`
		178	`write(*,'(1x,a12,i4,a10)') 'Vars : ',1,trim(vars(1))`
		179	`do i=2,ncol`
		180	`write(*,'(1x,a12,i4,a10)') ' ',i,trim(vars(i))`
		181	`enddo`
		182	`print*`
		183	`write(*,'(1x,a12,i4,f10.2)') 'Time : ',1,time(1)`
		184	`do i=2,ntim`
		185	`write(*,'(1x,a12,i4,f10.2)') ' ',i,time(i)`
		186	`enddo`
		187	`print*`
		188	`write(*,'(1x,a12,i4,i10)') 'Base date : ',1,basetime(1)`
		189	`write(*,'(1x,a12,i4,i10)') ' ',2,basetime(2)`
		190	`write(*,'(1x,a12,i4,i10)') ' ',3,basetime(3)`
		191	`write(*,'(1x,a12,i4,i10)') ' ',4,basetime(4)`
		192	`write(*,'(1x,a12,i4,i10)') ' ',5,basetime(5)`
		193	`print*`
		194	`write(*,'(1x,a12,i4,i10)') 'Time range : ',6,basetime(6)`
		195	`print*`
		196	`if ( itrigger.ne.0 ) then`
		197	`print*,'TRIGGER FIELD FOUND IN COLUMN ',itrigger`
		198	`print*`
		199	`endif`
		200
		201	`c Read and decode the selection criterion`
		202	`fcr = 10`
		203	`open(fcr,file=inp_criteria)`
		204	`ncmd=maxcmd`
		205	`call decode(fcr,cmd,ncmd,vars,ncol,time,ntim,regionf)`
		206	`close(fcr)`
		207
		208	`print*`
		209	`print*,'---- PSEUDO CODE FOR SELECTION --------------------------'`
		210	`print*`
		211	`call dumpcode(cmd,ncmd,vars,ncol,time,ntim)`
		212
		213	`c --------------------------------------------------------------`
		214	`c Loop over all trajectories - selection`
		215	`c --------------------------------------------------------------`
		216
		217	`c Prepare string and parameters for SPECIAL commands`
		218	`if ( cmd(1).eq.0 ) then`
		219
		220	`c Get command string`
		221	`j = 2`
		222	`len = nint(cmd(j))`
		223	`specialstr = ''`
		224	`do k=1,len`
		225	`j = j + 1`
		226	`specialstr = trim(specialstr)//char(nint(cmd(j)))`
		227	`enddo`
		228
		229	`c Get paramters`
		230	`j = j + 1`
		231	`nparam = nint(cmd(j))`
		232	`do k=1,nparam`
		233	`j = j + 1`
		234	`param(k) = cmd(j)`
		235	`enddo`
		236
		237	`endif`
		238
		239	`c Init the counter for selected trajectories`
		240	`ntraout = 0`
		241
		242	`c Loop over all trajectories`
		243	`do i=1,ntrainp`
		244
		245	`c Copy a single trajectory to <trainp1> and <trigger1>`
		246	`do j=1,ntim`
		247	`do k=1,ncol`
		248	`trainp1(j,k) = trainp(i,j,k)`
		249	`enddo`
		250	`trigger1(j) = trigger(i,j)`
		251	`enddo`
		252
		253	`C Skip the trajectory if missing data are found for positions`
		254	`isok = 1`
		255	`do j=1,ntim`
		256	`if ( trainp1(j,4).lt.0. ) isok = 0`
		257	`enddo`
		258
		259	`c Decide whether the trajectory is selected (handle SPECIAL commands)`
		260	`if ( isok.eq.1 ) then`
		261	`if (cmd(1).ne.0 ) then`
		262	`call select_tra (isok,cmd,ncmd,trainp1,trigger1,ntim,ncol)`
		263
		264	`else`
		265	`call special (isok,specialstr,trainp1,ntim,ncol,`
		266	`> vars,time,param,nparam)`
		267	`endif`
		268	`endif`
		269
		270	`c The trigger might be changed in the selection - copy it`
		271	`do j=1,ntim`
		272	`trigger(i,j) = trigger1(j)`
		273	`enddo`
		274
		275	`c Set flag for selected trajectories`
		276	`if (isok.eq.1) then`
		277	`selflag(i) = 1`
		278	`ntraout = ntraout + 1`
		279	`else`
		280	`selflag(i) = 0`
		281	`endif`
		282
		283	`enddo`
		284
		285	`c --------------------------------------------------------------`
		286	`c Write output trajectories`
		287	`c --------------------------------------------------------------`
		288
		289	`c ------ Write output trajectories -----------------------------`
		290	`if ( outformat.eq.'trajectory' ) then`
		291
		292
		293	`c Define the trigger field if it is not yet defined`
		294	`if ( ( addtrigger.eq.'-trigger' ).and.(itrigger.eq.0) ) then`
		295	`ncol = ncol + 1`
		296	`vars(ncol) = 'TRIGGER'`
		297	`itrigger = ncol`
		298	`endif`
		299
		300	`c Allocate memory for output trajectory`
		301	`allocate(traout(ntraout,ntim,ncol),stat=stat)`
		302	`if (stat.ne.0) stop '* error allocating array apply *'`
		303
		304	`c Set output trajectories`
		305	`j = 0`
		306	`do i=1,ntrainp`
		307	`if (selflag(i).eq.1) then`
		308	`j = j + 1`
		309	`traout(j,1:ntim,1:ncol) = trainp(i,1:ntim,1:ncol)`
		310	`if ( itrigger.ne.0 ) then`
		311	`traout(j,1:ntim,itrigger) = real(trigger(i,1:ntim))`
		312	`endif`
		313	`endif`
		314	`enddo`
		315
		316	`c Write trajectories`
		317	`call wopen_tra(fod,out_lslfile,ntraout,ntim,ncol,`
		318	`> basetime,vars,outmode)`
		319	`call write_tra(fod,traout,ntraout,ntim,ncol,outmode)`
		320	`call close_tra(fod,outmode)`
		321
		322	`c ------ Write index list -------------------------------------`
		323	`elseif ( outformat.eq.'index' ) then`
		324
		325	`open(10,file=out_lslfile)`
		326	`do i=1,ntrainp`
		327	`if ( selflag(i).eq.1) write(10,*) i`
		328	`enddo`
		329	`close(10)`
		330
		331	`c ------ Write boolean list -----------------------------------`
		332	`elseif ( outformat.eq.'boolean' ) then`
		333
		334	`open(10,file=out_lslfile)`
		335	`do i=1,ntrainp`
		336	`write(10,'(i1)') selflag(i)`
		337	`enddo`
		338	`close(10)`
		339
		340	`c ------ Write count -------------------------------------------`
		341	`elseif ( outformat.eq.'count' ) then`
		342
		343	`open(10,file=out_lslfile)`
		344	`write(10,'(i7)') ntraout`
		345	`close(10)`
		346
		347	`c ------ Write starting positions -----------------------------`
		348	`elseif ( outformat.eq.'startf' ) then`
		349
		350	`c Allocate memory for output trajectory`
		351	`allocate(traout(ntraout,1,ncol),stat=stat)`
		352	`if (stat.ne.0) stop '* error allocating array apply *'`
		353
		354	`c Set output trajectories`
		355	`j = 0`
		356	`do i=1,ntrainp`
		357	`if (selflag(i).eq.1) then`
		358	`j = j + 1`
		359	`traout(j,1,:) = trainp(i,1,:)`
		360	`endif`
		361	`enddo`
		362
		363	`c Write trajectories`
		364	`if (outmode.ne.-1) then`
		365
		366	`call wopen_tra(fod,out_lslfile,ntraout,1,ncol,`
		367	`> basetime,vars,outmode)`
		368	`call write_tra(fod,traout,ntraout,1,ncol,outmode)`
		369	`call close_tra(fod,outmode)`
		370
		371	`c Output as a triple list (corresponding to <startf> file)`
		372	`else`
		373
		374	`fid = 10`
		375	`open(fid,file=out_lslfile)`
		376	`do i=1,ntraout`
		377	`write(fid,'(3f10.3)') traout(i,1,2), ! longitude`
		378	`> traout(i,1,3), ! latitude`
		379	`> traout(i,1,4) ! pressure`
		380	`enddo`
		381	`close(fid)`
		382
		383	`endif`
		384
		385	`endif`
		386
		387	`c Write some status information, and end of program message`
		388	`print*`
		389	`print*,'---- STATUS INFORMATION --------------------------------'`
		390	`print*`
		391	`print*,' # input trajectories : ',ntrainp`
		392	`print*,' # output trajectories : ',ntraout`
		393	`print*`
		394	`print,' END OF PROGRAM SELECT *'`
		395	`print*,'========================================================='`
		396
		397	`stop`
		398
		399	`c Exception handling`
		400	`100 stop 'select: First column in input trajectory must be <time>'`
		401	`101 stop 'select: Input trajectory file is empty'`
		402
		403	`end`
		404
		405	`c --------------------------------------------------------------`
		406	`c Dump the command list`
		407	`c --------------------------------------------------------------`
		408
		409	`subroutine dumpcode(out,n,vars,nvars,times,ntimes)`
		410
		411	`c Write the command list to screen. The command list is decoded`
		412	`c by call to <decode>`
		413
		414	`implicit none`
		415
		416	`c Declaration of subroutine parameters`
		417	`integer n`
		418	`real out(n)`
		419	`integer nvars`
		420	`character*80 vars(nvars)`
		421	`integer ntimes`
		422	`real times(ntimes)`
		423
		424	`c A single command`
		425	`character*80 cmd`
		426	`character*80 var,mode,strtim`
		427	`integer nval`
		428	`integer ntim`
		429	`integer ivar,imode,icmd,itime`
		430
		431	`c Auxiliary variables`
		432	`integer i,j`
		433
		434	`c Loop through the complete list`
		435	`i=0`
		436	`100 if (i.lt.n) then`
		437
		438	`write(,) '---------------------------------------'`
		439
		440	`c Get command`
		441	`i=i+1`
		442	`icmd=nint(out(i))`
		443
		444	`c Special handling of SPECIAL commands`
		445	`if ( icmd.eq.0 ) then`
		446
		447	`c Write 'header' for SPECIAL command`
		448	`write(*,'(i5,f15.4,10x,a10)') i,out(i),'SPECIAL'`
		449
		450	`c Write command string`
		451	`i = i + 1`
		452	`icmd = nint(out(i))`
		453	`write(*,'(i5,f15.4,10x,a10)') i,out(i),'LEN(CMD)'`
		454	`do j=1,icmd`
		455	`i = i + 1`
		456	`ivar = nint(out(i))`
		457	`write(*,'(i5,f15.4,10x,a10)') i,out(i),char(ivar)`
		458	`enddo`
		459
		460	`c Write parameters`
		461	`i = i + 1`
		462	`nval = nint(out(i))`
		463	`write(*,'(i5,f15.4,10x,a10)') i,out(i),'#PARAMETER'`
		464	`do j=1,nval`
		465	`i=i+1`
		466	`if ( var.ne.'INPOLYGON' ) then`
		467	`write(*,'(i5,f15.4)') i,out(i)`
		468	`else`
		469	`write(*,'(i5,f15.4,a2)') i,out(i),char(nint(out(i)))`
		470	`endif`
		471	`enddo`
		472
		473	`c Nothing else to do - exit`
		474	`goto 200`
		475
		476	`endif`
		477
		478	`c Set the command`
		479	`if (icmd.eq. 1) cmd='GT'`
		480	`if (icmd.eq. 2) cmd='LT'`
		481	`if (icmd.eq. 3) cmd='IN'`
		482	`if (icmd.eq. 4) cmd='OUT'`
		483	`if (icmd.eq. 5) cmd='EQ'`
		484	`if (icmd.eq. 6) cmd='TRUE'`
		485	`if (icmd.eq. 7) cmd='FALSE'`
		486	`if (icmd.eq. 8) cmd='ALL'`
		487	`if (icmd.eq. 9) cmd='ANY'`
		488	`if (icmd.eq. 10) cmd='NONE'`
		489	`if (icmd.eq. -1) cmd='BEGIN'`
		490	`if (icmd.eq. -2) cmd='END'`
		491	`if (icmd.eq. -3) cmd='AND'`
		492	`if (icmd.eq. -4) cmd='OR'`
		493	`write(*,'(i5,f15.4,10x,a10)') i,out(i),trim(cmd)`
		494	`if (icmd.lt.0) goto 100`
		495
		496	`c Get variable`
		497	`i=i+1`
		498	`ivar=nint(out(i))`
		499	`if ( ivar.eq. -1 ) then`
		500	`var = 'DIST'`
		501	`elseif ( ivar.eq. -2 ) then`
		502	`var = 'DIST0'`
		503	`elseif ( ivar.eq. -3 ) then`
		504	`var = 'INPOLYGON'`
		505	`elseif ( ivar.eq. -4 ) then`
		506	`var = 'INBOX'`
		507	`elseif ( ivar.eq. -5 ) then`
		508	`var = 'INCIRCLE'`
		509	`elseif ( ivar.eq. -6 ) then`
		510	`var = 'INREGION'`
		511	`elseif ( ivar.eq. -7 ) then`
		512	`var = 'TRIGGER'`
9	michaesp	513	`elseif ( ivar.eq. -8 ) then`
		514	`var = 'VERT0'`
3	michaesp	515	`else`
		516	`var=vars(ivar)`
		517	`endif`
		518	`write(*,'(i5,f15.4,10x,a10)') i,out(i),trim(var)`
		519
		520	`c Get variable mode`
		521	`i=i+1`
		522	`imode=nint(out(i))`
15	michaesp	523	`if (imode.eq.1 ) mode='VALUE'`
		524	`if (imode.eq.2 ) mode='MEAN'`
		525	`if (imode.eq.3 ) mode='MAX'`
		526	`if (imode.eq.4 ) mode='MIN'`
		527	`if (imode.eq.5 ) mode='VAR'`
		528	`if (imode.eq.6 ) mode='SUM'`
		529	`if (imode.eq.7 ) mode='CHANGE'`
		530	`if (imode.eq.8 ) mode='DIFF'`
		531	`if (imode.eq.9 ) mode='RANGE'`
		532	`if (imode.eq.10) mode='ABS'`
3	michaesp	533	`write(*,'(i5,f15.4,10x,a10)') i,out(i),trim(mode)`
		534
		535	`c Get values`
		536	`i=i+1`
		537	`nval=nint(out(i))`
		538	`write(*,'(i5,f15.4,10x,a10)') i,out(i),'#PARAMETER'`
		539	`do j=1,nval`
		540	`i=i+1`
		541
		542	`if ( var.ne.'INPOLYGON' ) then`
		543	`write(*,'(i5,f15.4)') i,out(i)`
		544	`else`
		545	`write(*,'(i5,f15.4,a2)') i,out(i),char(nint(out(i)))`
		546	`endif`
		547	`enddo`
		548
9	michaesp	549	`c Get the number of times`
3	michaesp	550	`i=i+1`
		551	`ntim=nint(out(i))`
		552
9	michaesp	553	`c the number of times is variable - depending on TRIGGER`
		554	`if ( ntim .eq. -993 ) then`
		555	`write(*,'(i5,f15.4,7x,7x,a15)') i,out(i),'TIMES @ TRIGGER'`
		556
		557	`c Get the detailed list of times`
		558	`else`
		559	`write(*,'(i5,f15.4,7x,7x,a6)') i,out(i),'#TIMES'`
		560	`do j=1,ntim`
		561	`i=i+1`
		562	`write(*,'(i5,f15.4,f7.0)') i,out(i),times(nint(out(i)))`
		563	`enddo`
		564	`endif`
		565
3	michaesp	566	`c Get time mode`
		567	`i=i+1`
		568	`itime=nint(out(i))`
		569	`if (itime.eq.1) strtim='ALL'`
		570	`if (itime.eq.2) strtim='ANY'`
		571	`if (itime.eq.3) strtim='NONE'`
		572	`if (itime.lt.0) strtim='TRIGGER'`
		573
		574	`if ( strtim.ne.'TRIGGER' ) then`
		575	`write(*,'(i5,f15.4,10x,a10)') i,out(i),trim(strtim)`
		576	`else`
		577	`write(*,'(i5,f15.4,10x,a10,a3,i3)') i,out(i),trim(strtim),`
		578	`> ' ->',abs(itime)`
		579	`endif`
		580
		581	`goto 100`
		582
		583	`endif`
		584
		585	`c Exit point`
		586	`200 continue`
		587
		588	`write(,) '---------------------------------------'`
		589
		590	`end`
		591
		592
		593	`c --------------------------------------------------------------`
		594	`c Read and decode a selection set`
		595	`c --------------------------------------------------------------`
		596
		597	`subroutine decode(fid,out,n,vars,nvars,times,ntimes,regionf)`
		598
		599	`c A selection file is opened with file id <fid> and transformed`
		600	`c into a set of commands applied to the trajectories. On input`
		601	`c <n> sets the maximum dimension of <out>, on output it gives the`
		602	`c total length of the command string. The output is a list of`
		603	`c commands with the following format:`
		604	`c`
		605	`c out(i) = Command`
		606	`c out(i+1) = Column index of variable`
		607	`c out(i+2) = Mode for variable`
		608	`c out(i+3) = Number of parameters (n)`
		609	`c out(i+4:i+4+n) = Parameters`
		610	`c out(i+5+n) = Number of times`
		611	`c out(i+6+n:i+6+n+m) = List of time indices (m)`
		612	`c out(i+7+n+m) = Time mode`
		613	`c`
		614	`c For SPECIAL commands (to be coded in <special.f>) the format is:`
		615	`c`
		616	`c out(i) = Length of command string (n)`
		617	`c out(i+1:i+n) = Command string`
		618	`c out(i+n+1) = Number of parameters (m)`
		619	`c out(i+n+2:i+n+1+m) = List of parameters`
		620	`c`
		621	`c The following coding is used`
		622	`c`
		623	`c Command Variable mode Time mode`
		624	`c ---------- ------------- ---------`
		625	`c GT 1 VALUE 1 ALL 1`
		626	`c LT 2 MEAN 2 ANY 2`
		627	`c IN 3 MAX 3 NONE 3`
		628	`c OUT 4 MIN 4 TRIGGER -i (i the trigger index)`
		629	`c EQ 5 VAR 5`
		630	`c BEGIN -1 SUM 6`
		631	`c END -2 CHANGE 7`
		632	`c AND -3 DIFF 8`
9	michaesp	633	`c OR -4 RANGE 9`
15	michaesp	634	`c TRUE 6 ABS 10`
3	michaesp	635	`c FALSE 7`
		636	`c SPECIAL 0`
		637	`c ALL 8 (TRIGGER)`
		638	`c ANY 9 (TRIGGER)`
		639	`c NONE 10 (TRIGGER)`
		640
		641
		642	`c Several "implicit variables" are supported - out(i+1):`
		643	`c`
		644	`c DIST -1 : Path length of the trajectory`
		645	`c DIST0 -2 : Distance from starting position`
		646	`c INPOLYGON -3 : Specified polygon region`
		647	`c INBOX -4 : Longitude/latitude rectangle`
		648	`c INCIRCLE -5 : Within a specified radius`
		649	`c INREGION -6 : Within a specified rehion on the region file`
		650	`c TRIGGER -7 : Trigger field`
9	michaesp	651	`c VERT0 -8 : Vertical distance from starting position`
3	michaesp	652	`c`
		653	`c For the special commands BEGIN, END, AND and OR, only one field`
		654	`c in <out> is used.`
		655
		656	`implicit none`
		657
		658	`c Declaration of subroutine parameters`
		659	`integer fid`
		660	`integer n`
		661	`real out(n)`
		662	`integer nvars`
		663	`character*80 vars(nvars)`
		664	`integer ntimes`
		665	`real times(ntimes)`
		666	`character*80 regionf`
		667
		668	`c Numerical epsilon`
		669	`real eps`
		670	`parameter (eps=0.001)`
		671
		672	`c A single command term`
		673	`character*80 cmd`
		674	`character*80 var,mode`
		675	`integer nval`
		676	`real val(1000)`
		677	`integer ntim`
		678	`real tim(1000)`
		679	`character*80 tmode`
		680
		681	`c Specification of a polygon`
		682	`character*80 filename`
		683	`integer pn ! Number of entries in lat/lon poly`
		684	`real latpoly(1000) ! List of polygon latitudes`
		685	`real lonpoly(1000) ! List of polygon longitudes`
		686	`real loninpoly,latinpoly ! Lon/lat inside polygon`
		687
		688	`c Specification of a region`
		689	`real xcorner(4)`
		690	`real ycorner(4)`
		691	`integer iregion`
		692	`character*80 string`
		693
		694	`c Transformation to UPN (handling of logical operators)`
		695	`integer nlogical`
		696	`integer ilogical(n)`
		697	`real tmp(n)`
		698	`integer mlogical`
		699	`integer isor`
		700
		701	`c Auxiliary variables`
		702	`integer i,j`
		703	`integer j1,j2`
		704	`integer flag(ntimes)`
		705	`integer count`
		706	`integer ok`
9	michaesp	707	`integer itrigger,ttrigger`
3	michaesp	708
		709	`c Common block for initialisation of polygon check`
		710	`real tlonv(1000),vlat_c(1000),vlon_c(1000)`
		711	`real xlat_c,xlon_c`
		712	`integer ibndry,nv_c`
		713	`common /spolybndry/vlat_c,vlon_c,nv_c,xlat_c,xlon_c,tlonv,ibndry`
		714
		715	`c ------ Decode single commands ---------------------`
		716
		717	`c Reset the filename for polygons`
		718	`filename='nil'`
		719
		720	`c Reset the counter for logical commands`
		721	`nlogical=0`
		722
		723	`c Loop through all commands`
		724	`100 continue`
		725
		726	`c Read next command (handle special cases)`
		727	`read(fid,*) cmd`
		728
		729	`c Special handling of SPECIAL commands`
		730	`if ( cmd.eq.'SPECIAL' ) then`
		731
		732	`c Set the flag for SPECIAL command`
		733	`n = 1`
		734	`out(n) = 0.`
		735
		736	`c Read the command string`
		737	`read(fid,*) var`
		738
		739	`c Add the command string`
		740	`n = n + 1`
		741	`out(n) = len_trim(var)`
		742	`do j=1,len_trim(var)`
		743	`n = n + 1`
		744	`out(n) = ichar(var(j:j))`
		745	`enddo`
		746
		747	`c Read the parameters`
		748	`read(fid,*) nval`
		749	`read(fid,*) (val(i),i=1,nval)`
		750
		751	`c Add the parameters`
		752	`n = n + 1`
		753	`out(n)=real(nval)`
		754	`do i=1,nval`
		755	`n=n+1`
		756	`out(n)=val(i)`
		757	`enddo`
		758
		759	`c Goto exit point - nothing more top do`
		760	`goto 350`
		761
		762	`endif`
		763
		764	`c Handle structure commands`
		765	`if ( cmd.eq.'BEGIN') then`
		766	`out(1)=-1.`
		767	`n=1`
		768	`nlogical=1`
		769	`ilogical(1)=n`
		770	`goto 200`
		771	`elseif ( cmd.eq.'AND' ) then`
		772	`n=n+1`
		773	`out(n)=-3.`
		774	`nlogical=nlogical+1`
		775	`ilogical(nlogical)=n`
		776	`goto 200`
		777	`elseif ( cmd.eq.'OR' ) then`
		778	`n=n+1`
		779	`out(n)=-4.`
		780	`nlogical=nlogical+1`
		781	`ilogical(nlogical)=n`
		782	`goto 200`
		783	`elseif ( cmd.eq.'END' ) then`
		784	`n=n+1`
		785	`out(n)=-2.`
		786	`nlogical=nlogical+1`
		787	`ilogical(nlogical)=n`
		788	`goto 300`
		789	`endif`
		790
		791	`c Read other fields associated with the command`
		792	`read(fid,*) var,mode`
		793
		794	`c Read parameter`
		795	`if ( var.eq.'INPOLYGON' ) then`
		796	`read(fid,*) nval`
		797	`read(fid,*) filename`
		798	`filename = trim(filename)`
		799	`else`
		800	`read(fid,*) nval`
		801	`read(fid,*) (val(i),i=1,nval)`
		802	`endif`
		803
9	michaesp	804	`c Read times (on request, change to special trigger times)`
3	michaesp	805	`read(fid,*) ntim`
9	michaesp	806	`if ( ntim.eq.-993 ) then`
		807	`ttrigger = 1`
		808	`ntim = 1`
		809	`tim(1) = -999.`
		810	`else`
		811	`ttrigger = 0`
		812	`read(fid,*) (tim(i),i=1,ntim)`
		813	`endif`
3	michaesp	814	`read(fid,*) tmode`
		815
		816	`c Bring CAPITAL "TIME,LAT,LON,P" into "time,lat,lon,p"`
		817	`if (var.eq.'TIME') var='time'`
		818	`if (var.eq.'LAT' ) var='lat'`
		819	`if (var.eq.'LON' ) var='lon'`
		820	`if (var.eq.'P' ) var='p'`
		821
		822	`c If the time mode is 'TRIGGER', all times of a trajectory`
		823	`c must be considered`
		824	`if ( tmode.eq.'TRIGGER' ) then`
		825	`itrigger = nint(tim(1))`
		826	`ntim = 1`
		827	`tim(1) = -999.`
		828	`endif`
		829
		830	`c Special times: transform into real time`
		831	`do i=1,ntim`
		832	`if ( abs(tim(i)+996.).lt.eps ) tim(i)=times(1)`
		833	`if ( abs(tim(i)+995.).lt.eps ) tim(i)=times(ntimes)`
		834	`enddo`
		835
		836	`c Check whether times are valid`
		837	`ok=0`
		838	`do i=1,ntim`
		839	`if ( (abs(tim(i)+994.).gt.eps).and.`
		840	`> (abs(tim(i)+999.).gt.eps) )`
		841	`> then`
		842	`do j=1,ntimes`
		843	`if ( abs(tim(i)-times(j)).lt.eps ) then`
		844	`ok=ok+1`
		845	`endif`
		846	`enddo`
		847	`else`
		848	`ok=ok+1`
		849	`endif`
		850	`enddo`
		851	`if (ok.ne.ntim) goto 400`
		852
		853	`c Select all times which are included in the criterion`
		854	`do i=1,ntimes`
		855	`flag(i)=0`
		856	`enddo`
		857	`i=1`
		858	`150 if (i.le.ntim) then`
		859
		860	`c A list of times`
		861	`if ( (abs(tim(i)+994.).lt.eps) ) then`
		862	`j1=0`
		863	`do j=1,ntimes`
		864	`if ( abs(tim(i-1)-times(j)).lt.eps ) then`
		865	`j1=j`
		866	`endif`
		867	`enddo`
		868	`j2=0`
		869	`do j=1,ntimes`
		870	`if ( abs(tim(i+1)-times(j)).lt.eps ) then`
		871	`j2=j`
		872	`endif`
		873	`enddo`
		874	`if ( (j1.eq.0).or.(j2.eq.0) ) goto 400`
		875	`do j=j1,j2`
19	michaesp	876	`flag(j)=i`
3	michaesp	877	`enddo`
		878	`i=i+1`
		879
		880	`c Explicitly given time value`
		881	`else`
		882	`do j=1,ntimes`
		883	`if ( abs(tim(i)-times(j)).lt.eps ) then`
19	michaesp	884	`flag(j)=i`
3	michaesp	885	`endif`
		886	`enddo`
		887
		888	`endif`
		889
		890	`i=i+1`
		891	`goto 150`
		892
		893	`endif`
		894
		895	`c Write command identifier`
		896	`n=n+1`
		897	`if (cmd.eq.'GT' ) out(n)= 1.`
		898	`if (cmd.eq.'LT' ) out(n)= 2.`
		899	`if (cmd.eq.'IN' ) out(n)= 3.`
		900	`if (cmd.eq.'OUT' ) out(n)= 4.`
		901	`if (cmd.eq.'EQ' ) out(n)= 5.`
		902	`if (cmd.eq.'TRUE' ) out(n)= 6.`
		903	`if (cmd.eq.'FALSE' ) out(n)= 7.`
		904	`if (cmd.eq.'ALL ' ) out(n)= 8.`
		905	`if (cmd.eq.'ANY ' ) out(n)= 9.`
		906	`if (cmd.eq.'NONE ' ) out(n)=10.`
		907
		908	`c Write index for variable - force implicit trigger`
		909	`ok=0`
		910	`do j=1,nvars`
		911	`if (vars(j).eq.var) ok=j`
		912	`enddo`
		913
		914	`if (var.eq.'TRIGGER') ok = 0`
		915
		916	`if (ok.eq.0) then`
		917	`if (var.eq.'DIST') then`
		918	`ok = -1`
		919	`elseif (var.eq.'DIST0') then`
		920	`ok = -2`
		921	`elseif (var.eq.'INPOLYGON') then`
		922	`ok = -3`
		923	`elseif (var.eq.'INBOX') then`
		924	`ok = -4`
		925	`elseif (var.eq.'INCIRCLE') then`
		926	`ok = -5`
		927	`elseif (var.eq.'INREGION') then`
		928	`ok = -6`
		929	`elseif (var.eq.'TRIGGER') then`
		930	`ok = -7`
9	michaesp	931	`elseif (var.eq.'VERT0') then`
		932	`ok = -8`
3	michaesp	933	`else`
		934	`goto 400`
		935	`endif`
		936	`endif`
		937	`n=n+1`
		938	`out(n)=real(ok)`
		939
		940	`c Write mode for variable`
		941	`ok=0`
		942	`if (mode.eq.'VALUE' ) ok=1`
		943	`if (mode.eq.'MEAN' ) ok=2`
		944	`if (mode.eq.'MAX' ) ok=3`
		945	`if (mode.eq.'MIN' ) ok=4`
		946	`if (mode.eq.'VAR' ) ok=5`
		947	`if (mode.eq.'SUM' ) ok=6`
		948	`if (mode.eq.'CHANGE' ) ok=7`
		949	`if (mode.eq.'DIFF' ) ok=8`
9	michaesp	950	`if (mode.eq.'RANGE' ) ok=9`
15	michaesp	951	`if (mode.eq.'ABS' ) ok=10`
3	michaesp	952	`if (ok.eq.0) goto 400`
		953	`n=n+1`
		954	`out(n)=real(ok)`
		955
		956	`c Write the parameter values: INPOLYGON`
		957	`if ( var.eq.'INPOLYGON' ) then`
		958
		959	`n = n+1`
		960	`out(n) = len_trim(filename)`
		961	`do j=1,len_trim(filename)`
		962	`n = n + 1`
		963	`out(n) = ichar(filename(j:j))`
		964	`enddo`
		965
		966	`c Write parameter value: INREGION`
		967	`elseif ( var.eq.'INREGION' ) then`
		968
		969	`iregion = nint(val(1))`
		970
		971	`open(fid+1,file=regionf)`
		972
		973	`50 read(fid+1,*,end=51) string`
		974
		975	`if ( string(1:1).ne.'#' ) then`
		976	`call regionsplit(string,i,xcorner,ycorner)`
		977	`if ( i.eq.iregion ) goto 52`
		978	`endif`
		979
		980	`goto 50`
		981
		982	`51 close(fid+1)`
		983
		984	`print*,' ERROR: region ',iregion,' not found on ',`
		985	`> trim(regionf)`
		986	`stop`
		987
		988	`52 continue`
		989
		990	`n = n + 1`
		991	`out(n) = 8 ! Number of parameters`
		992	`do i=1,4`
		993	`n = n + 1`
		994	`out(n) = xcorner(i)`
		995	`enddo`
		996	`do i=1,4`
		997	`n = n + 1`
		998	`out(n) = ycorner(i)`
		999	`enddo`
		1000
		1001	`c Write parameter values: all other cases`
		1002	`else`
		1003	`n=n+1`
		1004	`out(n)=real(nval)`
		1005	`do i=1,nval`
		1006	`n=n+1`
		1007	`out(n)=val(i)`
		1008	`enddo`
		1009	`endif`
		1010
9	michaesp	1011	`c Special time handling: only trigger times are cosidered`
		1012	`if ( ttrigger.eq.1 ) then`
		1013	`n = n+1`
		1014	`out(n)=-993.`
		1015
3	michaesp	1016	`c All times are selected`
9	michaesp	1017	`elseif ( abs(tim(1)+999.).lt.eps ) then`
3	michaesp	1018	`n=n+1`
		1019	`out(n)=real(ntimes)`
		1020	`do i=1,ntimes`
		1021	`n=n+1`
		1022	`out(n)=real(i)`
		1023	`enddo`
		1024
		1025	`c A selection of times is given`
		1026	`else`
		1027	`count=0`
		1028	`do i=1,ntimes`
19	michaesp	1029	`if ( flag(i).ne.0 ) then`
		1030	`count=count+1`
		1031	`endif`
3	michaesp	1032	`enddo`
		1033	`n=n+1`
		1034	`out(n)=real(count)`
19	michaesp	1035	`do i=1,count`
		1036	`do j=1,ntimes`
		1037	`if (flag(j).eq.i) then`
		1038	`n = n + 1`
		1039	`out(n) = real(j)`
		1040	`endif`
		1041	`enddo`
3	michaesp	1042	`enddo`
19	michaesp	1043
3	michaesp	1044	`endif`
		1045
		1046	`c Write the time mode`
		1047	`if ( tmode.eq.'ALL') then`
		1048	`n=n+1`
		1049	`out(n)=1.`
		1050	`elseif ( tmode.eq.'ANY') then`
		1051	`n=n+1`
		1052	`out(n)=2.`
		1053	`elseif ( tmode.eq.'NONE') then`
		1054	`n=n+1`
		1055	`out(n)=3.`
		1056	`elseif ( tmode.eq.'TRIGGER') then`
		1057	`n=n+1`
		1058	`out(n)=-real(itrigger)`
		1059	`endif`
		1060
		1061	`c End loop: handle single command`
		1062	`200 continue`
		1063	`goto 100`
		1064
		1065	`c End loop: loop over all commands`
		1066	`300 continue`
		1067
		1068	`c ------ Read polygon file, if requested -----------`
		1069	`if ( filename.ne.'nil' ) then`
		1070
		1071	`print*`
		1072	`print*,`
		1073	`> '---- POLYGON --------------------------------------------'`
		1074
		1075	`print*`
		1076	`print*,'Filename = ',trim(filename)`
		1077	`print*`
		1078
		1079	`c Read list of polygon coordinates from file`
		1080	`pn = 0`
		1081	`open(fid+1,file=filename)`
		1082	`read(fid+1,*) loninpoly,latinpoly`
		1083	`print*,'Inside (lon/lat) =',loninpoly,latinpoly`
		1084	`print*`
		1085	`510 continue`
		1086	`pn = pn + 1`
		1087	`read(fid+1,*,end=511) lonpoly(pn),`
		1088	`> latpoly(pn)`
		1089
		1090	`print*,pn,lonpoly(pn),latpoly(pn)`
		1091
		1092	`goto 510`
		1093	`511 continue`
		1094	`pn = pn - 1`
		1095	`close(fid+1)`
		1096
		1097	`c Define the polygon boundaries`
		1098	`call DefSPolyBndry(latpoly,lonpoly,pn,latinpoly,loninpoly)`
		1099
		1100	`endif`
		1101
		1102	`c ------ Transform to UPN --------------------------`
		1103
		1104	`c Check whether logical commands are ok`
		1105	`mlogical=nint(out(ilogical(1)))`
		1106	`if ( mlogical.ne.-1) goto 400`
		1107	`mlogical=nint(out(ilogical(nlogical)))`
		1108	`if ( mlogical.ne.-2) goto 400`
		1109
		1110	`c No transformation necessary if only one command`
		1111	`if (nlogical.eq.2) goto 350`
		1112
		1113	`c Copy the output to temporary list`
		1114	`do i=1,n`
		1115	`tmp(i)=out(i)`
		1116	`enddo`
		1117
		1118	`c Set BEGIN statement`
		1119	`n=1`
		1120	`out(n)=-1.`
		1121
		1122	`c Reorder commands and transform to UPN`
		1123	`isor=0`
		1124	`do i=1,nlogical-1`
		1125
		1126	`c Get the logical command`
		1127	`mlogical=nint(out(ilogical(i)))`
		1128
		1129	`c Connecting OR`
		1130	`if (mlogical.eq.-4) then`
		1131	`if (isor.eq.1) then`
		1132	`n=n+1`
		1133	`out(n)=-4.`
		1134	`else`
		1135	`isor=1`
		1136	`endif`
		1137	`endif`
		1138
		1139	`c Put the command onto the stack`
		1140	`do j=ilogical(i)+1,ilogical(i+1)-1`
		1141	`n=n+1`
		1142	`out(n)=tmp(j)`
		1143	`enddo`
		1144
		1145	`c Connecting AND operator`
		1146	`if ( mlogical.eq.-3 ) then`
		1147	`n=n+1`
		1148	`out(n)=-3.`
		1149	`endif`
		1150
		1151	`enddo`
		1152
		1153	`c Set final connecting OR`
		1154	`if (isor.eq.1) then`
		1155	`n=n+1`
		1156	`out(n)=-4.`
		1157	`endif`
		1158
		1159	`c Set END statement`
		1160	`n=n+1`
		1161	`out(n)=-2.`
		1162
		1163	`c ------ Exit point ---------------------------------`
		1164
		1165	`350 continue`
		1166	`return`
		1167
		1168	`c ----- Exception handling --------------------------`
		1169	`400 print*,'Invalid selection criterion... Stop'`
		1170	`stop`
		1171
		1172	`end`
		1173
		1174
		1175	`c --------------------------------------------------------------`
		1176	`c Decide whether a trajectory is selected or not`
		1177	`c --------------------------------------------------------------`
		1178
		1179	`subroutine select_tra (select,cmd,n,tra,trigger,ntim,ncol)`
		1180
		1181	`c Decide whether a single trajectory is selected (<select=1>) or`
		1182	`c is not selected <select=0> according to the selection criterion`
		1183	`c given in <cmd(ncmd)>. The selection criterion <cmd(ncmd)> is`
		1184	`c returned from the call to the subroutine <decode>. The trajectory`
		1185	`c is given in <tra(ntim,ncol)> where <ntim> is the number of times`
		1186	`c and <ncol> is the number of columns.`
		1187	`c`
		1188	`c Important note: the structure of <tra(ntim,ncol)> must match to the`
		1189	`c call parameter <vars,nvars,times,ntimes> in subroutine <decode>.`
		1190
		1191	`implicit none`
		1192
		1193	`c Declaration of subroutine parameters`
		1194	`integer select`
		1195	`integer n`
		1196	`real cmd(n)`
		1197	`integer ntim,ncol`
		1198	`real tra(ntim,ncol)`
		1199	`integer trigger(ntim)`
		1200
		1201	`c Numerical epsilon (for test of equality)`
		1202	`real eps`
		1203	`parameter (eps=0.000001)`
		1204
		1205	`c A single command and the associated field`
		1206	`integer icmd,ivar,imode,itime,nsel,nval`
		1207	`integer time(ntim)`
		1208	`real param(100)`
		1209	`real var (ntim)`
		1210	`integer intvar(ntim)`
		1211
		1212	`c Boolean values for a single time, a single command and build-up`
		1213	`integer stack(100)`
		1214	`integer nstack`
		1215	`integer istrue(ntim)`
		1216	`integer decision`
		1217
		1218	`c Auxiliary variables`
		1219	`integer i,j,k`
		1220	`real tmp,mea`
		1221	`integer istack1,istack2`
		1222	`real lat0,lon0,lat1,lon1`
		1223	`real length(ntim)`
		1224	`integer flag`
		1225	`real dist`
		1226	`real xcorner(4),ycorner(4)`
		1227	`integer iparam`
		1228	`character ch`
9	michaesp	1229	`real varmin,varmax`
		1230	`real lev0,lev1`
3	michaesp	1231
		1232	`c Common block for initialisation of polygon check`
		1233	`real tlonv(1000),vlat_c(1000),vlon_c(1000)`
		1234	`real xlat_c,xlon_c`
		1235	`integer ibndry,nv_c`
		1236	`common /spolybndry/vlat_c,vlon_c,nv_c,xlat_c,xlon_c,tlonv,ibndry`
		1237
		1238	`c Externals`
		1239	`real sdis ! Spherical distance`
		1240	`external sdis`
		1241	`integer inregion ! Boolean flag for regions`
		1242	`external inregion`
		1243
		1244	`c Reset the decision stack (with locical values)`
		1245	`nstack=0`
		1246
		1247	`c Loop through the complete command list`
		1248	`i=0`
		1249	`100 if (i.lt.n) then`
		1250
		1251	`c --- Get the command -------------------------------`
		1252	`i=i+1`
		1253	`icmd=nint(cmd(i))`
		1254
		1255	`c --- Handle structural commands (BEGIN, END, AND, OR)`
		1256
		1257	`c Handle BEGIN statement`
		1258	`if ( icmd.eq.-1) then`
		1259	`nstack=0`
		1260	`goto 200`
		1261	`endif`
		1262
		1263	`c Handle END statement`
		1264	`if (icmd.eq.-2) then`
		1265	`goto 300`
		1266	`endif`
		1267
		1268	`c Handle AND statement`
		1269	`if (icmd.eq.-3) then`
		1270	`istack1=stack(nstack)`
		1271	`nstack=nstack-1`
		1272	`istack2=stack(nstack)`
		1273	`if ((istack1.eq.1).and.(istack2.eq.1)) then`
		1274	`stack(nstack)=1`
		1275	`else`
		1276	`stack(nstack)=0`
		1277	`endif`
		1278	`goto 200`
		1279	`endif`
		1280
		1281	`c Handle OR statement`
		1282	`if (icmd.eq.-4) then`
		1283	`istack1=stack(nstack)`
		1284	`nstack=nstack-1`
		1285	`istack2=stack(nstack)`
		1286	`if ((istack1.eq.1).or.(istack2.eq.1)) then`
		1287	`stack(nstack)=1`
		1288	`else`
		1289	`stack(nstack)=0`
		1290	`endif`
		1291	`goto 200`
		1292	`endif`
		1293
		1294	`c --- Get all command details (parameters, modes, times)`
		1295
		1296	`c Get variable (<ivar> gets the column index in <tra>)`
		1297	`i=i+1`
		1298	`ivar=nint(cmd(i))`
		1299
		1300	`c Get variable mode`
		1301	`i=i+1`
		1302	`imode=nint(cmd(i))`
		1303
		1304	`c Get parameter values`
		1305	`i=i+1`
		1306	`nval=nint(cmd(i))`
		1307	`do j=1,nval`
		1308	`i=i+1`
		1309	`param(j)=cmd(i)`
		1310	`enddo`
		1311
		1312	`c Get times (<time(j)> gets the row indices of <tra>)`
		1313	`i=i+1`
		1314	`nsel=nint(cmd(i))`
9	michaesp	1315	`if ( nsel .eq. -993 ) then`
		1316	`nsel = 0`
		1317	`do k=1,ntim`
		1318	`if ( trigger(k).ne.0 ) then`
		1319	`nsel = nsel + 1`
		1320	`time(nsel) = k`
		1321	`endif`
		1322	`enddo`
		1323	`else`
		1324	`do j=1,nsel`
		1325	`i=i+1`
		1326	`time(j)=nint(cmd(i))`
		1327	`enddo`
		1328	`endif`
3	michaesp	1329
9	michaesp	1330	`c If no times are selected, exit with non-select status`
		1331	`if ( nsel.eq.0 ) then`
		1332	`stack(1) = 0`
		1333	`goto 300`
		1334	`endif`
		1335
3	michaesp	1336	`c Get time mode`
		1337	`i=i+1`
		1338	`itime=nint(cmd(i))`
		1339
		1340	`c --- Prepare field values for analysis -----------`
		1341
		1342	`c Implicit variable: DIST`
		1343	`if ( ivar.eq. -1 ) then`
		1344	`length(1) = 0.`
		1345	`do j=2,ntim`
		1346	`lon0 = tra(j-1,2)`
		1347	`lat0 = tra(j-1,3)`
		1348	`lon1 = tra(j ,2)`
		1349	`lat1 = tra(j ,3)`
		1350	`length(j) = length(j-1) + sdis(lon0,lat0,lon1,lat1)`
		1351	`enddo`
		1352	`do j=1,nsel`
		1353	`var(j) = length( time(j) )`
		1354	`enddo`
		1355
		1356
		1357	`c Implict variable: DIST0`
		1358	`elseif ( ivar.eq. -2 ) then`
		1359	`do j=1,nsel`
		1360	`lon0 = tra(1 ,2)`
		1361	`lat0 = tra(1 ,3)`
		1362	`lon1 = tra(time(j),2)`
		1363	`lat1 = tra(time(j),3)`
		1364	`var(j) = sdis(lon0,lat0,lon1,lat1)`
		1365	`enddo`
		1366
		1367	`c Implict variable: INPOLYGON`
		1368	`elseif ( ivar.eq. -3 ) then`
		1369	`do j=1,nsel`
		1370	`lon1 = tra(time(j),2)`
		1371	`lat1 = tra(time(j),3)`
		1372	`call LctPtRelBndry(lat1,lon1,flag)`
		1373	`if ( (flag.eq.1).or.(flag.eq.2) ) then`
		1374	`var(j) = 1.`
		1375	`else`
		1376	`var(j) = 0.`
		1377	`endif`
		1378	`enddo`
		1379
		1380	`c Implict variable: INBOX`
		1381	`elseif ( ivar.eq. -4 ) then`
		1382	`do j=1,nsel`
		1383	`lon1 = tra(time(j),2)`
		1384	`lat1 = tra(time(j),3)`
		1385	`if ( ( lon1.ge.param(1) ).and. ! lonmin`
		1386	`> ( lon1.le.param(2) ).and. ! lonmax`
		1387	`> ( lat1.ge.param(3) ).and. ! latmin`
		1388	`> ( lat1.le.param(4) ) ) ! latmax`
		1389	`> then`
		1390	`var(j) = 1`
		1391	`else`
		1392	`var(j) = 0`
		1393	`endif`
		1394	`enddo`
		1395
		1396	`c Implict variable: INCIRCLE (lonc=param(1),latc=param(2),radius=param(3))`
		1397	`elseif ( ivar.eq. -5 ) then`
		1398	`do j=1,nsel`
		1399	`lon1 = tra(time(j),2)`
		1400	`lat1 = tra(time(j),3)`
		1401
		1402	`dist = sdis( lon1,lat1,param(1),param(2) )`
		1403
		1404	`if ( dist.le.param(3) ) then`
		1405	`var(j) = 1`
		1406	`else`
		1407	`var(j) = 0`
		1408	`endif`
		1409	`enddo`
		1410
		1411	`c Implict variable: INREGION (xcorner=param(1..4),ycorner=param(5..8) )`
		1412	`elseif ( ivar.eq.-6 ) then`
		1413
		1414	`do j=1,4`
		1415	`xcorner(j) = param(j )`
		1416	`ycorner(j) = param(j+4)`
		1417	`enddo`
		1418
		1419	`do j=1,nsel`
		1420	`lon1 = tra(time(j),2)`
		1421	`lat1 = tra(time(j),3)`
		1422	`var(j) = inregion (lon1,lat1,xcorner,ycorner)`
		1423	`enddo`
		1424
		1425	`c Implict variable: TRIGGER`
		1426	`elseif ( ivar.eq. -7 ) then`
		1427	`do j=1,nsel`
		1428	`intvar(j) = trigger( time(j) )`
		1429	`enddo`
		1430
9	michaesp	1431	`c Implicit variable: VERT0`
		1432	`elseif ( ivar.eq. -8 ) then`
		1433	`do j=1,nsel`
		1434	`lev0 = tra(1 ,4)`
		1435	`lev1 = tra(time(j),4)`
		1436	`var(j) = lev0 - lev1`
		1437	`enddo`
		1438
3	michaesp	1439	`c Explicit variable (column index <ivar>)`
		1440	`else`
		1441	`do j=1,nsel`
		1442	`var(j) = tra(time(j),ivar)`
		1443	`enddo`
		1444
		1445	`endif`
		1446
		1447	`c Take MEAN of the variable (mean of selected times)`
		1448	`if (imode.eq.2) then`
		1449	`tmp=0.`
		1450	`do j=1,nsel`
		1451	`tmp=tmp+var(j)`
		1452	`enddo`
		1453	`var(1)=tmp/real(nsel)`
		1454	`nsel=1`
		1455
		1456	`c Take MAX of the variable (maximum of selected times)`
		1457	`elseif (imode.eq.3) then`
		1458	`tmp=var(1)`
		1459	`do j=2,nsel`
		1460	`if (var(j).gt.tmp) tmp=var(j)`
		1461	`enddo`
		1462	`var(1)=tmp`
		1463	`nsel=1`
		1464
		1465	`c Take MIN of the variable (minimum of selected times)`
		1466	`elseif (imode.eq.4) then`
		1467	`tmp=var(1)`
		1468	`do j=2,nsel`
		1469	`if (var(j).lt.tmp) tmp=var(j)`
		1470	`enddo`
		1471	`var(1)=tmp`
		1472	`nsel=1`
		1473
		1474	`c Take VAR of the variable (variance over all selected times)`
		1475	`elseif (imode.eq.5) then`
		1476	`tmp=0.`
		1477	`do j=1,nsel`
		1478	`tmp=tmp+var(j)`
		1479	`enddo`
		1480	`mea=tmp/real(nsel)`
		1481	`do j=1,nsel`
		1482	`tmp=tmp+(var(j)-mea)**2`
		1483	`enddo`
		1484	`var(1)=1./real(nsel-1)*tmp`
		1485	`nsel=1`
		1486
		1487	`c Take SUM of the variable (sum over all selected times)`
		1488	`elseif (imode.eq.6) then`
		1489	`tmp=0.`
		1490	`do j=1,nsel`
		1491	`tmp=tmp+var(j)`
		1492	`enddo`
		1493	`var(1)=tmp`
		1494	`nsel=1`
		1495
9	michaesp	1496	`c Take CHANGE of the variable (absolute difference between first and last time)`
3	michaesp	1497	`elseif (imode.eq.7) then`
		1498	`var(1)=abs(var(1)-var(nsel))`
		1499	`nsel=1`
		1500
		1501	`c Take DIFF of the variable (first minus last time)`
		1502	`elseif (imode.eq.8) then`
		1503	`var(1)=var(1)-var(nsel)`
		1504	`nsel=1`
		1505
9	michaesp	1506	`c Take RANGE of the variable`
		1507	`elseif (imode.eq.9) then`
		1508	`varmax=var(1)`
		1509	`varmin=var(1)`
		1510	`do j=2,nsel`
		1511	`if (var(j).gt.varmax) varmax=var(j)`
		1512	`if (var(j).lt.varmin) varmin=var(j)`
		1513	`enddo`
		1514	`var(1) = varmax - varmin`
		1515	`nsel=1`
15	michaesp	1516
		1517	`c Take the absolute value of the variable`
		1518	`elseif (imode.eq.10) then`
		1519	`do j=1,nsel`
		1520	`var(j)=abs( var(j) )`
		1521	`enddo`
3	michaesp	1522	`endif`
		1523
		1524	`c --- Apply the operators to the single values ---`
		1525
		1526	`do j=1,nsel`
		1527
		1528	`c GT`
		1529	`if (icmd.eq.1) then`
		1530	`if (var(j).gt.param(1)) then`
		1531	`istrue(j)=1`
		1532	`else`
		1533	`istrue(j)=0`
		1534	`endif`
		1535
		1536	`c LT`
		1537	`elseif (icmd.eq.2) then`
		1538	`if (var(j).lt.param(1)) then`
		1539	`istrue(j)=1`
		1540	`else`
		1541	`istrue(j)=0`
		1542	`endif`
		1543
		1544	`c IN`
		1545	`elseif (icmd.eq.3) then`
		1546	`if ( (var(j).gt.param(1)).and.`
		1547	`> (var(j).lt.param(2)) )`
		1548	`> then`
		1549	`istrue(j)=1`
		1550	`else`
		1551	`istrue(j)=0`
		1552	`endif`
		1553
		1554	`c OUT`
		1555	`elseif (icmd.eq.4) then`
		1556	`if ( (var(j).lt.param(1)).or.`
		1557	`> (var(j).gt.param(2)) )`
		1558	`> then`
		1559	`istrue(j)=1`
		1560	`else`
		1561	`istrue(j)=0`
		1562	`endif`
		1563
		1564	`c EQ`
		1565	`elseif (icmd.eq.5) then`
		1566	`if (abs(var(j)-param(1)).lt.eps) then`
		1567	`istrue(j)=1`
		1568	`else`
		1569	`istrue(j)=0`
		1570	`endif`
		1571
		1572	`c TRUE`
		1573	`elseif (icmd.eq.6) then`
		1574	`if (abs(var(j)).lt.eps) then`
		1575	`istrue(j)=0`
		1576	`else`
		1577	`istrue(j)=1`
		1578	`endif`
		1579
		1580	`c FALSE`
		1581	`elseif (icmd.eq.7) then`
		1582	`if (abs(var(j)).lt.eps) then`
		1583	`istrue(j)=1`
		1584	`else`
		1585	`istrue(j)=0`
		1586	`endif`
		1587
		1588	`c ALL`
		1589	`elseif (icmd.eq.8) then`
		1590	`istrue(j) = 1`
		1591	`do k=1,nval`
		1592	`iparam = nint(param(k))-1`
		1593	`if (btest(intvar(j),iparam).eqv..false.) then`
		1594	`istrue(j) = 0`
		1595	`endif`
		1596	`enddo`
		1597
		1598	`c ANY`
		1599	`elseif (icmd.eq.9) then`
		1600	`istrue(j) = 0`
		1601	`do k=1,nval`
		1602	`iparam = nint(param(k))-1`
		1603	`if (btest(intvar(j),iparam).eqv..true.) then`
		1604	`istrue(j) = 1`
		1605	`endif`
		1606	`enddo`
		1607
		1608	`c NONE`
		1609	`elseif (icmd.eq.10) then`
		1610	`istrue(j) = 1`
		1611	`do k=1,nval`
		1612	`iparam = nint(param(k))-1`
		1613	`if (btest(intvar(j),iparam).eqv..true.) then`
		1614	`istrue(j) = 0`
		1615	`endif`
		1616	`enddo`
		1617
		1618	`endif`
		1619
		1620	`enddo`
		1621
		1622	`c --- Determine the overall boolean value ----------`
		1623
		1624	`c ALL`
		1625	`if (itime.eq.1) then`
		1626	`decision=1`
		1627	`do j=1,nsel`
		1628	`if (istrue(j).eq.0) then`
		1629	`decision=0`
		1630	`goto 110`
		1631	`endif`
		1632	`enddo`
		1633	`110 continue`
		1634
		1635	`c ANY`
		1636	`elseif (itime.eq.2) then`
		1637	`decision=0`
		1638	`do j=1,nsel`
		1639	`if (istrue(j).eq.1) then`
		1640	`decision=1`
		1641	`goto 120`
		1642	`endif`
		1643	`enddo`
		1644	`120 continue`
		1645
		1646	`c NONE`
		1647	`elseif (itime.eq.3) then`
		1648	`decision=1`
		1649	`do j=1,nsel`
		1650	`if (istrue(j).eq.1) then`
		1651	`decision=0`
		1652	`goto 130`
		1653	`endif`
		1654	`enddo`
		1655	`130 continue`
		1656
		1657	`c TRIGGER`
		1658	`elseif (itime.lt.0) then`
		1659	`decision=1`
		1660	`do j=1,nsel`
		1661	`if (istrue(j).eq.1) then`
		1662	`trigger(j) = ior( trigger(j), 2**(abs(itime)-1) )`
		1663	`endif`
		1664	`enddo`
		1665
		1666	`endif`
		1667
		1668	`c --- Put the new boolean value onto the stack`
		1669
		1670	`nstack=nstack+1`
		1671	`stack(nstack)=decision`
		1672
		1673	`c Exit point for loop`
		1674	`200 continue`
		1675	`goto 100`
		1676
		1677	`endif`
		1678
		1679	`c Return the decision (selected or non-selected)`
		1680	`300 continue`
		1681
		1682	`select=stack(1)`
		1683
		1684	`end`
		1685
		1686
		1687	`c --------------------------------------------------------------------------`
		1688	`c Split a region string and get corners of the domain`
		1689	`c --------------------------------------------------------------------------`
		1690
		1691	`subroutine regionsplit(string,iregion,xcorner,ycorner)`
		1692
		1693	`c The region string comes either as <lonw,lone,lats,latn> or as <lon1,lat1,`
		1694	`c lon2,lat2,lon3,lat3,lon4,lat4>: split it into ints components and get the`
		1695	`c four coordinates for the region`
		1696
		1697	`implicit none`
		1698
		1699	`c Declaration of subroutine parameters`
		1700	`character*80 string`
		1701	`real xcorner(4),ycorner(4)`
		1702	`integer iregion`
		1703
		1704	`c Local variables`
		1705	`integer i,n`
		1706	`integer il,ir`
		1707	`real subfloat (80)`
		1708	`integer stat`
		1709	`integer len`
		1710
		1711	`c ------- Split the string`
		1712	`i = 1`
		1713	`n = 0`
		1714	`stat = 0`
		1715	`il = 1`
		1716	`len = len_trim(string)`
		1717
		1718	`100 continue`
		1719
		1720	`c Find start of a substring`
		1721	`do while ( stat.eq.0 )`
		1722	`if ( string(i:i).ne.' ' ) then`
		1723	`stat = 1`
		1724	`il = i`
		1725	`else`
		1726	`i = i + 1`
		1727	`endif`
		1728	`enddo`
		1729
		1730	`c Find end of substring`
		1731	`do while ( stat.eq.1 )`
		1732	`if ( ( string(i:i).eq.' ' ) .or. ( i.eq.len ) ) then`
		1733	`stat = 2`
		1734	`ir = i`
		1735	`else`
		1736	`i = i + 1`
		1737	`endif`
		1738	`enddo`
		1739
		1740	`c Convert the substring into a number`
		1741	`if ( stat.eq.2 ) then`
		1742	`n = n + 1`
		1743	`read(string(il:ir),*) subfloat(n)`
		1744	`stat = 0`
		1745	`endif`
		1746
		1747	`if ( i.lt.len ) goto 100`
		1748
		1749
		1750	`c -------- Get the region number`
		1751
		1752	`iregion = nint(subfloat(1))`
		1753
		1754	`c -------- Get the corners of the region`
		1755
		1756	`if ( n.eq.5 ) then ! lonw(2),lone(3),lats(4),latn(5)`
		1757
		1758	`xcorner(1) = subfloat(2)`
		1759	`ycorner(1) = subfloat(4)`
		1760
		1761	`xcorner(2) = subfloat(3)`
		1762	`ycorner(2) = subfloat(4)`
		1763
		1764	`xcorner(3) = subfloat(3)`
		1765	`ycorner(3) = subfloat(5)`
		1766
		1767	`xcorner(4) = subfloat(2)`
		1768	`ycorner(4) = subfloat(5)`
		1769
		1770	`elseif ( n.eq.9 ) then ! lon1,lat1,lon2,lat2,lon3,lon4,lat4`
		1771
		1772	`xcorner(1) = subfloat(2)`
		1773	`ycorner(1) = subfloat(3)`
		1774
		1775	`xcorner(2) = subfloat(4)`
		1776	`ycorner(2) = subfloat(5)`
		1777
		1778	`xcorner(3) = subfloat(6)`
		1779	`ycorner(3) = subfloat(7)`
		1780
		1781	`xcorner(4) = subfloat(8)`
		1782	`ycorner(4) = subfloat(9)`
		1783
		1784	`else`
		1785
		1786	`print*,' ERROR: invalid region specification '`
		1787	`print*,' ',trim(string)`
		1788	`stop`
		1789
		1790	`endif`
		1791
		1792
		1793	`end`
		1794
		1795	`c --------------------------------------------------------------------------`
		1796	`c Decide whether lat/lon point is in or out of region`
		1797	`c --------------------------------------------------------------------------`
		1798
		1799	`integer function inregion (lon,lat,xcorner,ycorner)`
		1800
		1801	`c Decide whether point (lon/lat) is in the region specified by <xcorner(1..4),`
		1802	`c ycorner(1..4).`
		1803
		1804	`implicit none`
		1805
		1806	`c Declaration of subroutine parameters`
		1807	`real lon,lat`
		1808	`real xcorner(4),ycorner(4)`
		1809
		1810	`c Local variables`
		1811	`integer flag`
		1812	`real xmin,xmax,ymin,ymax`
		1813	`integer i`
		1814
		1815	`c Reset the flag`
		1816	`flag = 0`
		1817
		1818	`c Set some boundaries`
		1819	`xmax = xcorner(1)`
		1820	`xmin = xcorner(1)`
		1821	`ymax = ycorner(1)`
		1822	`ymin = ycorner(1)`
		1823	`do i=2,4`
		1824	`if (xcorner(i).lt.xmin) xmin = xcorner(i)`
		1825	`if (xcorner(i).gt.xmax) xmax = xcorner(i)`
		1826	`if (ycorner(i).lt.ymin) ymin = ycorner(i)`
		1827	`if (ycorner(i).gt.ymax) ymax = ycorner(i)`
		1828	`enddo`
		1829
		1830	`c Do the tests - set flag=1 if all tests pased`
		1831	`if (lon.lt.xmin) goto 970`
		1832	`if (lon.gt.xmax) goto 970`
		1833	`if (lat.lt.ymin) goto 970`
		1834	`if (lat.gt.ymax) goto 970`
		1835
		1836	`if ((lon-xcorner(1))*(ycorner(2)-ycorner(1))-`
		1837	`> (lat-ycorner(1))*(xcorner(2)-xcorner(1)).gt.0.) goto 970`
		1838	`if ((lon-xcorner(2))*(ycorner(3)-ycorner(2))-`
		1839	`> (lat-ycorner(2))*(xcorner(3)-xcorner(2)).gt.0.) goto 970`
		1840	`if ((lon-xcorner(3))*(ycorner(4)-ycorner(3))-`
		1841	`> (lat-ycorner(3))*(xcorner(4)-xcorner(3)).gt.0.) goto 970`
		1842	`if ((lon-xcorner(4))*(ycorner(1)-ycorner(4))-`
		1843	`> (lat-ycorner(4))*(xcorner(1)-xcorner(4)).gt.0.) goto 970`
		1844
		1845	`flag = 1`
		1846
		1847	`c Return the value`
		1848	`970 continue`
		1849
		1850	`inregion = flag`
		1851
		1852	`return`
		1853
		1854	`end`
		1855
		1856
		1857	`c --------------------------------------------------------------------------`
		1858	`c Spherical distance between lat/lon points`
		1859	`c --------------------------------------------------------------------------`
		1860
		1861	`real function sdis(xp,yp,xq,yq)`
		1862	`c`
		1863	`c calculates spherical distance (in km) between two points given`
		1864	`c by their spherical coordinates (xp,yp) and (xq,yq), respectively.`
		1865	`c`
		1866	`real re`
		1867	`parameter (re=6370.)`
		1868	`real pi180`
		1869	`parameter (pi180=3.14159/180.)`
		1870	`real xp,yp,xq,yq,arg`
		1871
		1872	`arg=sin(pi180yp)sin(pi180*yq)+`
		1873	`> cos(pi180yp)cos(pi180yq)cos(pi180*(xp-xq))`
		1874	`if (arg.lt.-1.) arg=-1.`
		1875	`if (arg.gt.1.) arg=1.`
		1876
		1877	`sdis=re*acos(arg)`
		1878
		1879	`end`
		1880
		1881
		1882	`c ****************************************************************`
		1883	`c * Given some spherical polygon S and some point X known to be *`
		1884	`c * located inside S, these routines will determine if an arbit- *`
		1885	`c * -rary point P lies inside S, outside S, or on its boundary. *`
		1886	`c * The calling program must first call DefSPolyBndry to define *`
		1887	`c * the boundary of S and the point X. Any subsequent call to *`
		1888	`c * subroutine LctPtRelBndry will determine if some point P lies *`
		1889	`c * inside or outside S, or on its boundary. (Usually *`
		1890	`c * DefSPolyBndry is called once, then LctPrRelBndry is called *`
		1891	`c * many times). *`
		1892	`c * *`
		1893	`c * REFERENCE: Bevis, M. and Chatelain, J.-L. (1989) *`
		1894	`c * Maflaematical Geology, vol 21. *`
		1895	`c * VERSION 1.0 *`
		1896	`c ****************************************************************`
		1897
		1898	`Subroutine DefSPolyBndry(vlat,vlon,nv,xlat, xlon)`
		1899
		1900	`c ****************************************************************`
		1901	`c * This mmn entry point is used m define ~e spheric~ polygon S *`
		1902	`c * and the point X. *`
		1903	`c * ARGUMENTS: *`
		1904	`c * vlat,vlon (sent) ... vectors containing the latitude and *`
		1905	`c * longitude of each vertex of the *`
		1906	`c * spherical polygon S. The ith.vertex is *`
		1907	`c * located at [vlat(i),vlon(i)]. *`
		1908	`c * nv (sent) ... the number of vertices and sides in the *`
		1909	`c * spherical polygon S *`
		1910	`c * xlat,xlon (sent) ... latitude and longitude of some point X *`
		1911	`c * located inside S. X must not be located *`
		1912	`c * on any great circle that includes two *`
		1913	`c * vertices of S. *`
		1914	`c * *`
		1915	`c * UNITS AND SIGN CONVENTION: *`
		1916	`c * Latitudes and longitudes are specified in degrees. *`
		1917	`c * Latitudes are positive to the north and negative to the *`
		1918	`c * south. *`
		1919	`c * Longitudes are positive to the east and negative to the *`
		1920	`c * west. *`
		1921	`c * *`
		1922	`c * VERTEX ENUMERATION: *`
		1923	`c * The vertices of S should be numbered sequentially around the *`
		1924	`c * border of the spherical polygon. Vertex 1 lies between vertex*`
		1925	`c * nv and vertex 2. Neighbouring vertices must be seperated by *`
		1926	`c * less than 180 degrees. (In order to generate a polygon side *`
		1927	`c * whose arc length equals or exceeds 180 degrees simply *`
		1928	`c * introduce an additional (pseudo)vertex). Having chosen *`
		1929	`c * vertex 1, the user may number the remaining vertices in *`
		1930	`c * either direction. However if the user wishes to use the *`
		1931	`c * subroutine SPA to determine the area of the polygon S (Bevis *`
		1932	`c * & Cambareri, 1987, Math. Geol., v.19, p. 335-346) then he or *`
		1933	`c * she must follow the convention whereby in moving around the *`
		1934	`c * polygon border in the direction of increasing vertex number *`
		1935	`c * clockwise bends occur at salient vertices. A vertex is *`
		1936	`c * salient if the interior angle is less than 180 degrees. *`
		1937	`c * (In the case of a convex polygon this convention implies *`
		1938	`c * that vertices are numbered in clockwise sequence). *`
		1939	`c ****************************************************************`
		1940
		1941	`implicit none`
		1942
		1943	`integer mxnv,nv`
		1944
		1945	`c ----------------------------------------------------------------`
		1946	`c Edit next statement to increase maximum number of vertices that`
		1947	`c may be used to define the spherical polygon S`
		1948	`c The value of parameter mxnv in subroutine LctPtRelBndry must match`
		1949	`c that of parameter mxnv in this subroutine, as assigned above.`
		1950	`c ----------------------------------------------------------------`
		1951	`parameter (mxnv=2000)`
		1952
		1953	`real vlat(nv),vlon(nv),xlat,xlon,dellon`
		1954	`real tlonv(mxnv),vlat_c(mxnv),vlon_c(mxnv),xlat_c,xlon_c`
		1955	`integer i,ibndry,nv_c,ip`
		1956
		1957	`data ibndry/0/`
		1958
		1959	`common/spolybndry/vlat_c,vlon_c,nv_c,xlat_c,xlon_c,tlonv,ibndry`
		1960
		1961	`if (nv.gt.mxnv) then`
		1962	`print *,'nv exceeds maximum allowed value'`
		1963	`print *,'adjust parameter mxnv in subroutine DefSPolyBndry'`
		1964	`stop`
		1965	`endif`
		1966
		1967	`ibndry=1 ! boundary defined at least once (flag)`
		1968	`nv_c=nv ! copy for named common`
		1969	`xlat_c=xlat ! . . . .`
		1970	`xlon_c=xlon !`
		1971
		1972	`do i=1,nv`
		1973	`vlat_c(i)=vlat(i) ! "`
		1974	`vlon_c(i)=vlon(i) !`
		1975
		1976	`call TrnsfmLon(xlat,xlon,vlat(i),vlon(i),tlonv(i))`
		1977
		1978	`if (i.gt.1) then`
		1979	`ip=i-1`
		1980	`else`
		1981	`ip=nv`
		1982	`endif`
		1983
		1984	`if ((vlat(i).eq.vlat(ip)).and.(vlon(i).eq.vlon(ip))) then`
		1985	`print *,'DefSPolyBndry detects user error:'`
		1986	`print *,'vertices ',i,' and ',ip,' are not distinct'`
		1987	`print*,'lat ',i,ip,vlat(i),vlat(ip)`
		1988	`print*,'lon ',i,ip,vlon(i),vlon(ip)`
		1989	`stop`
		1990	`endif`
		1991
		1992	`if (tlonv(i).eq.tlonv(ip)) then`
		1993	`print *,'DefSPolyBndry detects user error:'`
		1994	`print *,'vertices ',i,' & ',ip,' on same gt. circle as X'`
		1995	`stop`
		1996	`endif`
		1997
		1998	`if (vlat(i).eq.(-vlat(ip))) then`
		1999	`dellon=vlon(i)-vlon(ip)`
		2000	`if (dellon.gt.+180.) dellon=dellon-360.`
		2001	`if (dellon.lt.-180.) dellon=dellon-360.`
		2002	`if ((dellon.eq.+180.0).or.(dellon.eq.-180.0)) then`
		2003	`print *,'DefSPolyBndry detects user error:'`
		2004	`print *,'vertices ',i,' and ',ip,' are antipodal'`
		2005	`stop`
		2006	`endif`
		2007	`endif`
		2008	`enddo`
		2009
		2010	`return`
		2011
		2012	`end`
		2013
		2014
		2015	`c ****************************************************************`
		2016
		2017	`Subroutine LctPtRelBndry(plat,plon,location)`
		2018
		2019	`c ****************************************************************`
		2020
		2021	`c ****************************************************************`
		2022	`c * This routine is used to see if some point P is located *`
		2023	`c * inside, outside or on the boundary of the spherical polygon *`
		2024	`c * S previously defined by a call to subroutine DefSPolyBndry. *`
		2025	`c * There is a single restriction on point P: it must not be *`
		2026	`c * antipodal to the point X defined in the call to DefSPolyBndry*`
		2027	`c * (ie.P and X cannot be seperated by exactly 180 degrees). *`
		2028	`c * ARGUMENTS: *`
		2029	`c * plat,plon (sent)... the latitude and longitude of point P *`
		2030	`c * location (returned)... specifies the location of P: *`
		2031	`c * location=0 implies P is outside of S *`
		2032	`c * location=1 implies P is inside of S *`
		2033	`c * location=2 implies P on boundary of S *`
		2034	`c * location=3 implies user error (P is *`
		2035	`c * antipodal to X) *`
		2036	`c * UNFfS AND SIGN CONVENTION: *`
		2037	`c * Latitudes and longitudes are specified in degrees. *`
		2038	`c * Latitudes are positive to the north and negative to the *`
		2039	`c * south. *`
		2040	`c * Longitudes are positive to the east and negative to the *`
		2041	`c * west. *`
		2042	`c ****************************************************************`
		2043
		2044	`implicit none`
		2045
		2046	`integer mxnv`
		2047
		2048	`c ----------------------------------------------------------------`
		2049	`c The statement below must match that in subroutine DefSPolyBndry`
		2050	`c ----------------------------------------------------------------`
		2051
		2052	`parameter (mxnv=2000)`
		2053
		2054	`real tlonv(mxnv),vlat_c(mxnv),vlon_c(mxnv),xlat_c,xlon_c`
		2055	`real plat,plon,vAlat,vAlon,vBlat,vBlon,tlonA,tlonB,tlonP`
		2056	`real tlon_X,tlon_P,tlon_B,dellon`
		2057	`integer i,ibndry,nv_c,location,icross,ibrngAB,ibrngAP,ibrngPB`
		2058	`integer ibrng_BX,ibrng_BP,istrike`
		2059
		2060	`common/spolybndry/vlat_c,vlon_c,nv_c,xlat_c,xlon_c,tlonv,ibndry`
		2061
		2062	`if (ibndry.eq.0) then ! user has never defined the bndry`
		2063	`print*,'Subroutine LctPtRelBndry detects user error:'`
		2064	`print*,'Subroutine DefSPolyBndry must be called before'`
		2065	`print*,'subroutine LctPtRelBndry can be called'`
		2066	`stop`
		2067	`endif`
		2068
		2069	`if (plat.eq.(-xlat_c)) then`
		2070	`dellon=plon-xlon_c`
		2071	`if (dellon.lt.(-180.)) dellon=dellon+360.`
		2072	`if (dellon.gt.+180.) dellon=dellon-360.`
		2073	`if ((dellon.eq.+180.0).or.(dellon.eq.-180.)) then`
		2074	`print*,'Warning: LctPtRelBndry detects case P antipodal`
		2075	`> to X'`
		2076	`print*,'location of P relative to S is undetermined'`
		2077	`location=3`
		2078	`return`
		2079	`endif`
		2080	`endif`
		2081
		2082	`location=0 ! default ( P is outside S)`
		2083	`icross=0 ! initialize counter`
		2084
		2085	`if ((plat.eq.xlat_c).and.(plon.eq.xlon_c)) then`
		2086	`location=1`
		2087	`return`
		2088	`endif`
		2089
		2090
		2091	`call TrnsfmLon (xlat_c,xlon_c,plat,plon,tlonP)`
		2092
		2093	`do i=1,nv_c ! start of loop over sides of S`
		2094
		2095	`vAlat=vlat_c(i)`
		2096	`vAlon=vlon_c(i)`
		2097	`tlonA=tlonv(i)`
		2098
		2099	`if (i.lt.nv_c) then`
		2100	`vBlat=vlat_c(i+1)`
		2101	`vBlon=vlon_c(i+1)`
		2102	`tlonB=tlonv(i+1)`
		2103	`else`
		2104	`vBlat=vlat_c(1)`
		2105	`vBlon=vlon_c(1)`
		2106	`tlonB=tlonv(1)`
		2107	`endif`
		2108
		2109	`istrike=0`
		2110
		2111	`if (tlonP.eq.tlonA) then`
		2112	`istrike=1`
		2113	`else`
		2114	`call EastOrWest(tlonA,tlonB,ibrngAB)`
		2115	`call EastOrWest(tlonA,tlonP,ibrngAP)`
		2116	`call EastOrWest(tlonP,tlonB,ibrngPB)`
		2117
		2118
		2119	`if((ibrngAP.eq.ibrngAB).and.(ibrngPB.eq.ibrngAB)) istrike=1`
		2120	`endif`
		2121
		2122
		2123	`if (istrike.eq.1) then`
		2124
		2125	`if ((plat.eq.vAlat).and.(plon.eq.vAlon)) then`
		2126	`location=2 ! P lies on a vertex of S`
		2127	`return`
		2128	`endif`
		2129	`call TrnsfmLon(vAlat,vAlon,xlat_c,xlon_c,tlon_X)`
		2130	`call TrnsfmLon(vAlat,vAlon,vBlat,vBlon,tlon_B)`
		2131	`call TrnsfmLon(vAlat,vAlon,plat,plon,tlon_P)`
		2132
		2133	`if (tlon_P.eq.tlon_B) then`
		2134	`location=2 ! P lies on side of S`
		2135	`return`
		2136	`else`
		2137	`call EastOrWest(tlon_B,tlon_X,ibrng_BX)`
		2138	`call EastOrWest(tlon_B,tlon_P,ibrng_BP)`
		2139	`if(ibrng_BX.eq.(-ibrng_BP)) icross=icross+1`
		2140	`endif`
		2141
		2142	`endif`
		2143	`enddo ! end of loop over the sides of S`
		2144
		2145
		2146	`c if the arc XP crosses the boundary S an even number of times then P`
		2147	`c is in S`
		2148
		2149	`if (mod(icross,2).eq.0) location=1`
		2150
		2151	`return`
		2152
		2153	`end`
		2154
		2155
		2156	`c ****************************************************************`
		2157
		2158	`subroutine TrnsfmLon(plat,plon,qlat,qlon,tranlon)`
		2159
		2160	`c ****************************************************************`
		2161	`c * This subroutine is required by subroutines DefSPolyBndry & *`
		2162	`c * LctPtRelBndry. It finds the 'longitude' of point Q in a *`
		2163	`c * geographic coordinate system for which point P acts as a *`
		2164	`c * 'north pole'. SENT: plat,plon,qlat,qlon, in degrees. *`
		2165	`c * RETURNED: tranlon, in degrees. *`
		2166	`c ****************************************************************`
		2167
		2168	`implicit none`
		2169
		2170	`real pi,dtr,plat,plon,qlat,qlon,tranlon,t,b`
		2171	`parameter (pi=3.141592654,dtr=pi/180.0)`
		2172
		2173	`if (plat.eq.90.) then`
		2174	`tranlon=qlon`
		2175	`else`
		2176	`t=sin((qlon-plon)dtr)cos(qlat*dtr)`
		2177	`b=sin(dtrqlat)cos(platdtr)-cos(qlatdtr)sin(platdtr)`
		2178	`> cos((qlon-plon)dtr)`
		2179	`tranlon=atan2(t,b)/dtr`
		2180	`endif`
		2181
		2182	`return`
		2183	`end`
		2184
		2185	`c ****************************************************************`
		2186
		2187	`subroutine EastOrWest(clon,dlon,ibrng)`
		2188
		2189	`c ****************************************************************`
		2190	`c * This subroutine is required by subroutine LctPtRelBndry. *`
		2191	`c * This routine determines if in travelling the shortest path *`
		2192	`c * from point C (at longitude clon) to point D (at longitude *`
		2193	`c * dlon) one is heading east, west or neither. *`
		2194	`c * SENT: clon,dlon; in degrees. RETURNED: ibrng *`
		2195	`c * (1=east,-1=west, 0=neither). *`
		2196	`c ****************************************************************`
		2197
		2198	`implicit none`
		2199	`real clon,dlon,del`
		2200	`integer ibrng`
		2201	`del=dlon-clon`
		2202	`if (del.gt.180.) del=del-360.`
		2203	`if (del.lt.-180.) del=del+360.`
		2204	`if ((del.gt.0.0).and.(del.ne.180.)) then`
		2205	`ibrng=-1 ! (D is west of C)`
		2206	`elseif ((del.lt.0.0).and.(del.ne.-180.)) then`
		2207	`ibrng=+1 ! (D is east of C)`
		2208	`else`
		2209	`ibrng=0 ! (D north or south of C)`
		2210	`endif`
		2211	`return`
		2212	`end`

Subversion Repositories lagranto.ecmwf

(root)/trunk/select/select.f @ 5 – Rev 19