WebSVN – lagranto.ecmwf – Blame – /trunk/lib/iotra.f

Rev	Author	Line No.	Line
3	michaesp	1	`c ****************************************************************`
		2	`c * This package provides IO routines for trajectories. A file *`
		3	`c * is characterised by the filename <filename> and the file *`
		4	`c * identifier <fid>. Different modes <mode> are supported: *`
		5	`c * mode=1: ascii, sorted by trajectory; *`
		6	`c * mode=2: ascii, sorted by time; *`
		7	`c * mode=3: fortran (unformatted) *`
		8	`c * mode=4: IVE netcdf (for compatibiltzy reasons) *`
		9	`c * A trajectory set is given as 3d array <tra(ntra,ntim,ncol)> *`
		10	`c * where <ntra> is the number of trajectories, <ntim> the *`
		11	`c * number of times of each trajectory and <ncol> the number of *`
		12	`c * columns of the trajectory. The first 4 columns are: time, *`
		13	`c * longitude, latitude, pressure. The other columns are traced *`
		14	`c * fields. The complete list of all columns is given in the *`
		15	`c * array <vars(ncol)>. Finally, the reference date is given in *`
		16	`c * the array <time(6)=year,month,day,hour,time length of the *`
		17	`c * trajectory (hour,min)>. *`
		18	`c * *`
		19	`c * Author: Michael Sprenger, September 2008 *`
		20	`c ****************************************************************`
		21
		22	`c ----------------------------------------------------------------`
		23	`c Open a trajectory file for reading`
		24	`c ----------------------------------------------------------------`
		25
		26	`subroutine ropen_tra(fid,filename,ntra,ntim,ncol,time,vars,mode)`
		27
		28	`implicit none`
		29
		30	`c Declaration of subroutine parameters`
		31	`integer fid`
		32	`character*80 filename`
		33	`integer mode`
		34	`integer ntra,ntim,ncol`
		35	`integer time(6)`
		36	`character*80 vars(ncol)`
		37
		38	`c Auxiliary variables`
		39	`integer vardim(4)`
		40	`real varmin(4),varmax(4),stag(4)`
		41	`real mdv`
		42	`character*80 cfn`
		43	`integer ierr`
		44	`integer i`
		45	`integer nvars`
		46
		47	`c Open file`
		48	`if (mode.eq.1) then`
		49	`fid = 10`
		50	`open(fid,file=filename)`
		51	`elseif (mode.eq.2) then`
		52	`fid = 10`
		53	`open(fid,file=filename)`
		54	`elseif (mode.eq.3) then`
		55	`open(fid,file=filename,form='unformatted')`
		56	`elseif (mode.eq.4) then`
		57	`call cdfopn(filename,fid,ierr)`
5	michaesp	58	`elseif (mode.eq.5) then`
		59	`print*,' ERROR: Reading KML not supported'`
		60	`stop`
3	michaesp	61	`endif`
		62
		63	`c Read header information`
		64	`call read_hea(fid,time,vars,ntra,ntim,ncol,mode)`
		65
		66	`end`
		67
		68
		69	`c ----------------------------------------------------------------`
		70	`c Open a trajectory file for wrinting`
		71	`c ----------------------------------------------------------------`
		72
		73	`subroutine wopen_tra(fid,filename,ntra,ntim,ncol,time,vars,mode)`
		74
		75	`implicit none`
		76
		77	`c Declaration of subroutine parameters`
		78	`integer fid`
		79	`character*80 filename`
		80	`integer mode`
		81	`integer ntra,ntim,ncol`
		82	`integer time(6)`
		83	`character*80 vars(ncol)`
		84
		85	`c Auxiliary variables`
		86	`integer vardim(4)`
		87	`real varmin(4),varmax(4),stag(4)`
		88	`real mdv`
		89	`character*80 cfn`
		90	`integer ierr`
		91	`integer i`
		92	`character*80 varname`
		93	`real rtime(6)`
		94
		95	`c Open file`
		96	`if (mode.eq.1) then`
		97	`fid = 10`
		98	`open(fid,file=filename)`
		99	`elseif (mode.eq.2) then`
		100	`fid = 10`
		101	`open(fid,file=filename)`
		102	`elseif (mode.eq.3) then`
		103	`open(fid,file=filename,form='unformatted')`
		104	`elseif (mode.eq.4) then`
		105	`vardim(1)=ntra`
		106	`vardim(2)=1`
		107	`vardim(3)=1`
		108	`vardim(4)=1`
		109	`cfn =trim(filename)//'_cst'`
		110	`mdv =-999.98999`
		111	`call crecdf(filename,fid,varmin,varmax,3,cfn,ierr)`
5	michaesp	112	`elseif (mode.eq.5) then`
		113	`fid = 10`
		114	`open(fid,file=filename)`
3	michaesp	115	`endif`
		116
		117	`c Write header information`
		118	`call write_hea(fid,time,vars,ntra,ntim,ncol,mode)`
		119
		120	`end`
		121
		122
		123	`c ----------------------------------------------------------------`
		124	`c Read a trajectory`
		125	`c ----------------------------------------------------------------`
		126
		127	`subroutine read_tra(fid,tra,ntra,ntim,ncol,mode)`
		128
		129	`implicit none`
		130
		131	`c Declaration of subroutine parameters`
		132	`integer fid`
		133	`integer ntim`
		134	`integer ncol`
		135	`integer ntra`
		136	`real tra(ntra,ntim,ncol)`
		137	`integer mode`
		138
		139	`c Auxiliary variables`
		140	`integer i,j,n`
		141	`real arr(ntra)`
		142	`integer ntimes`
		143	`real times(1000)`
		144	`integer ierr`
		145	`character*80 vars(ncol+2)`
		146	`integer nvars`
		147
		148	`c Read ascii mode, sorted by trajectory (mode=1)`
		149	`if (mode.eq.1) then`
		150	`read(fid,*,end=100)`
		151	`do n=1,ntra`
		152	`do i=1,ntim`
		153	`read(fid,*,end=110) (tra(n,i,j),j=1,ncol)`
		154	`enddo`
		155	`enddo`
		156
		157	`c Read ascii mode, sorted by time (mode=2)`
		158	`elseif (mode.eq.2) then`
		159	`read(fid,*,end=100)`
		160	`do i=1,ntim`
		161	`do n=1,ntra`
		162	`read(fid,*,end=100) (tra(n,i,j),j=1,ncol)`
		163	`enddo`
		164	`enddo`
		165
		166	`c Read fortran mode (mode=3)`
		167	`elseif (mode.eq.3) then`
		168	`read(fid) tra`
		169
		170	`c Read IVE netcdf mode (mode=4)`
		171	`elseif (mode.eq.4) then`
		172	`call gettimes(fid,times,ntimes,ierr)`
		173	`call getvars(fid,nvars,vars,ierr)`
		174	`do i=1,ntim`
		175	`do j=1,ncol`
		176	`if (j.eq.1) then`
		177	`do n=1,ntra`
		178	`tra(n,i,1)=times(i)`
		179	`enddo`
		180	`else`
		181	`call getdat(fid,vars(j),times(i),0,arr,ierr)`
		182	`do n=1,ntra`
		183	`tra(n,i,j)=arr(n)`
		184	`enddo`
		185	`endif`
		186	`enddo`
		187	`enddo`
		188
		189	`endif`
		190
		191	`return`
		192
		193	`c End of file has been reached: set negative <fid>`
		194	`100 fid=-fid`
		195	`return`
		196
		197	`c Error: incomplete trajectory`
		198	`110 print*,'<read_tra>: Incomplete trajectory... Stop'`
		199	`stop`
		200
		201	`end`
		202
		203
		204	`c ----------------------------------------------------------------`
		205	`c Write a trajectory`
		206	`c ----------------------------------------------------------------`
		207
		208	`subroutine write_tra(fid,tra,ntra,ntim,ncol,mode)`
		209
		210	`implicit none`
		211
		212	`c Declaration of subroutine parameters`
		213	`integer fid`
		214	`integer ntim`
		215	`integer ncol`
		216	`integer ntra`
		217	`real tra(ntra,ntim,ncol)`
		218	`integer mode`
		219
		220	`c Auxiliary variables`
		221	`integer i,j,n`
		222	`real arr(ntra)`
		223	`integer ierr`
		224	`real time`
		225	`character*80 vars(ncol+2)`
		226	`integer nvars`
5	michaesp	227	`character*20 lonstr,latstr,levstr`
		228	`character*80 outstr`
		229	`real ref_z(3000),ref_p(3000),ref_t(3000)`
		230	`real lev`
		231	`character*80 path`
3	michaesp	232
		233	`c Write ascii mode, sorted by trajectory (mode=1)`
		234	`if (mode.eq.1) then`
		235	`do n=1,ntra`
		236	`write(fid,*)`
		237	`do i=1,ntim`
		238
		239	`c Avoid ugly *s or missing space in output`
		240	`do j=5,ncol`
		241	`if ( abs(tra(n,i,j)).gt.9999.) then`
		242	`print*,'Format problem : ',tra(n,i,j),' -> -999.99'`
		243	`tra(n,i,j) = -999.99`
		244	`endif`
		245	`enddo`
		246
		247	`write(fid,'(1f7.2,f9.2,f8.2,i6,100f10.3)')`
		248	`> (tra(n,i,j),j=1,3), ! time, lon, lat`
		249	`> nint(tra(n,i,4)), ! p`
		250	`> (tra(n,i,j),j=5,ncol) ! fields`
		251	`enddo`
		252	`enddo`
		253
		254	`c Write ascii mode, sorted by time (mode=2)`
		255	`elseif (mode.eq.2) then`
		256	`do i=1,ntim`
		257	`write(fid,*)`
		258	`do n=1,ntra`
		259
		260	`c Avoid ugly *s or missing space in output`
		261	`do j=5,ncol`
		262	`if ( abs(tra(n,i,j)).gt.9999.) then`
		263	`print*,'Format problem : ',tra(n,i,j),' -> -999.99'`
		264	`tra(n,i,j) = -999.99`
		265	`endif`
		266	`enddo`
		267
		268	`write(fid,'(1f7.2,f9.2,f8.2,i6,100f10.3)')`
		269	`> (tra(n,i,j),j=1,3), ! time, lon, lat`
		270	`> nint(tra(n,i,4)), ! p`
		271	`> (tra(n,i,j),j=5,ncol) ! fields`
		272	`enddo`
		273	`enddo`
		274
		275	`c Write fortran mode (mode=3)`
		276	`elseif (mode.eq.3) then`
		277	`write(fid) tra`
		278
		279	`c Write netcdf mode (mode=4)`
		280	`elseif (mode.eq.4) then`
		281	`call getvars(fid,nvars,vars,ierr)`
		282	`do i=1,ntim`
		283	`time=tra(1,i,1)`
		284	`do j=2,ncol`
		285	`do n=1,ntra`
		286	`arr(n)=tra(n,i,j)`
		287	`enddo`
		288	`call putdat(fid,vars(j),time,0,arr,ierr)`
		289	`enddo`
		290	`enddo`
5	michaesp	291
		292	`c Write KML mode (mode=5)`
		293	`elseif (mode.eq.5) then`
		294
		295	`call getenv('DYN_TOOLS',path)`
		296	`path = trim(path)//'/lagranto.ecmwf/goodies/'`
		297
		298	`open(fid+1,file=trim(path)//'reformat.refprof')`
		299
		300	`do n=1,6`
		301	`read(fid+1,*)`
		302	`enddo`
		303	`do n=1,3000`
		304	`read(fid+1,*) ref_z(n),ref_t(n),ref_p(n)`
		305	`ref_p(n) = 0.01 * ref_p(n)`
		306	`enddo`
		307
		308	`close(fid+1)`
		309
		310	`do n=1,ntra`
		311	`write(fid,"(A)") '<Placemark>'`
		312	`write(fid,"(A)") '<name>Absolute Extruded</name>'`
		313	`write(fid,"(A)") '<styleUrl>#yellowkLineGreenPoly</styleUrl>'`
		314	`write(fid,"(A)") '<LineString>'`
		315	`write(fid,"(A)") '<extrude>1</extrude>'`
		316	`write(fid,"(A)") '<tessellate>1</tessellate>'`
		317	`write(fid,"(A)") '<altitudeMode>absolute</altitudeMode>'`
		318	`write(fid,"(A)") '<coordinates>'`
		319
		320	`do i=1,ntim`
		321	`write(lonstr,*) tra(n,i,2)`
		322	`write(latstr,*) tra(n,i,3)`
		323
		324	`call binary(lev,tra(n,i,4),ref_z,ref_p)`
		325	`write(levstr,*) lev`
		326
		327	`outstr = trim(adjustl(lonstr))//','//`
		328	`> trim(adjustl(latstr))//','//`
		329	`> trim(adjustl(levstr))`
		330
		331	`write(fid,"(A)") outstr`
		332
		333	`enddo`
		334
		335	`write(fid,*) '</coordinates>'`
		336	`write(fid,*) '</LineString>'`
		337	`write(fid,*) '</Placemark>'`
		338	`enddo`
		339
3	michaesp	340	`endif`
		341
		342	`end`
		343
		344
		345	`c ----------------------------------------------------------------`
		346	`c Read header from trajectory file`
		347	`c ----------------------------------------------------------------`
		348
		349	`subroutine read_hea(fid,time,vars,ntra,ntim,ncol,mode)`
		350
		351	`implicit none`
		352
		353	`c Declaration of subroutine parameters`
		354	`integer fid`
		355	`integer time(6)`
		356	`integer ntra,ntim,ncol`
		357	`character*80 vars(ncol)`
		358	`integer mode`
		359
		360	`c Auxiliary variables`
		361	`integer i`
		362	`character ch(500)`
		363	`character*500 str`
		364	`integer ich(500)`
		365	`integer isstr,ileft,iright`
		366	`character*80 varname`
		367	`real rtime(6)`
		368	`integer ierr`
		369	`integer nvars`
		370	`character*15 str1`
		371	`character str2`
		372	`character*13 str3`
		373	`character*4 str4`
		374	`character*80 linestr`
		375	`integer itmp1,itmp2`
		376
		377	`c Read ascii format (mode=1,2)`
		378	`if ( (mode.eq.1).or.(mode.eq.2) ) then`
		379
		380	`c Read the time specification (old and new format)`
		381	`read(fid,'(a80)') linestr`
		382
		383	`if ( linestr(1:14).eq.'Reference date' ) then`
		384	`read(linestr,'(a15,i4,i2,i2,a1,i2,i2,a13,i8,a4)')`
		385	`> str1,`
		386	`> time(1),time(2),time(3),str2,time(4),time(5),`
		387	`> str3,time(6),str4`
		388
		389	`elseif ( linestr(1:11).eq.'time period' ) then`
		390	`read(linestr,'(a12,i4,i2,i2,a1,i2,a4,i6,a1,i2,a2)')`
		391	`> str1,`
		392	`> time(1),time(2),time(3),str2,time(4),`
		393	`> str3,itmp1,str3,itmp2,str4`
		394	`time(5) = 0`
		395	`time(6) = itmp1 * 60 + itmp2`
		396
		397	`endif`
		398
		399	`c Skip the empty line and read field names`
		400	`read(fid,*)`
		401	`read(fid,'(a500)',end=100) str`
		402	`do i=1,500`
		403	`ch(i)=str(i:i)`
		404	`enddo`
		405
		406	`c Split the input string`
		407	`isstr=0`
		408	`nvars=0`
		409	`do i=1,500`
		410	`if ( (isstr.eq.0).and.(ch(i).ne.' ') ) then`
		411	`isstr=1`
		412	`ileft=i`
		413	`elseif ( (isstr.eq.1).and.(ch(i).eq.' ') ) then`
		414	`isstr=0`
		415	`iright=i-1`
		416	`nvars=nvars+1`
		417	`vars(nvars)=str(ileft:iright)`
		418	`endif`
		419	`enddo`
		420
		421	`c Skip the empty line`
		422	`read(fid,*,end=100)`
		423
		424	`c Read fortran mode (mode=3)`
		425	`elseif (mode.eq.3) then`
		426	`read(fid) ntra,ntim,ncol`
		427	`read(fid) time`
		428	`read(fid) vars`
		429
		430	`c Read IVE netcdf mode (mode=4)`
		431	`elseif (mode.eq.4) then`
		432	`call getvars(fid,nvars,vars,ierr)`
		433	`varname='BASEDATE'`
		434	`call getdat(fid,varname,0.,0,rtime,ierr)`
		435	`do i=1,6`
		436	`time(i)=nint(rtime(i))`
		437	`enddo`
		438	`endif`
		439
		440	`return`
		441
		442	`c End of file has been reached`
		443	`100 fid=-fid`
		444	`return`
		445
		446	`c Excetion handling`
		447	`110 print*,'<read_hea>: Unexspected time format.... Stop'`
		448	`stop`
		449
		450	`end`
		451
		452
		453	`c ----------------------------------------------------------------`
		454	`c Write header to trajectory file (in ascii mode)`
		455	`c ----------------------------------------------------------------`
		456
		457	`subroutine write_hea(fid,time,vars,ntra,ntim,ncol,mode)`
		458
		459	`implicit none`
		460
		461	`c Declaration of subroutine parameters`
		462	`integer fid`
		463	`integer time(6)`
		464	`integer ntra,ntim,ncol`
		465	`character*80 vars(ncol)`
		466	`integer mode`
		467
		468	`c Auxiliary variables`
		469	`integer i`
		470	`character*500 str`
		471	`character*4 str1`
		472	`character*2 str2,str3,str4,str5,str6`
		473	`integer vardim(4)`
		474	`real varmin(4),varmax(4),stag(4)`
		475	`real mdv`
		476	`integer ierr`
		477	`character*80 varname`
		478	`real rtime(6)`
		479	`integer nvars`
		480
		481	`c Write ascii format (mode=1,2)`
		482	`if ( (mode.eq.1).or.(mode.eq.2) ) then`
		483
		484	`c Get the strings for output`
		485	`write(str1,'(i4)') time(1)`
		486	`write(str2,'(i2)') time(2)`
		487	`write(str3,'(i2)') time(3)`
		488	`write(str4,'(i2)') time(4)`
		489	`write(str5,'(i2)') time(5)`
		490	`if (time(2).eq. 0) str2(1:1)='0'`
		491	`if (time(3).eq. 0) str3(1:1)='0'`
		492	`if (time(4).eq. 0) str4(1:1)='0'`
		493	`if (time(5).eq. 0) str5(1:1)='0'`
		494	`if (time(2).lt.10) str2(1:1)='0'`
		495	`if (time(3).lt.10) str3(1:1)='0'`
		496	`if (time(4).lt.10) str4(1:1)='0'`
		497	`if (time(5).lt.10) str5(1:1)='0'`
		498
		499	`c Write the time specification`
		500	`write(fid,'(a15,a4,a2,a2,a1,a2,a2,a13,i8,a4)')`
		501	`> 'Reference date ',`
		502	`> str1,str2,str3,'_',str4,str5,`
		503	`> ' / Time range',time(6), ' min'`
		504	`write(fid,*)`
		505
		506	`c Write variable names`
		507	`str=''`
		508	`do i=1,ncol`
		509	`str=trim(str)//trim(vars(i))`
		510	`enddo`
		511	`write(fid,'(a6,a9,a8,a6,100a10)') (trim(vars(i)),i=1,ncol)`
		512	`write(fid,'(a6,a9,a8,a6,100a10)')`
		513	`> '------','---------','--------','------',`
		514	`> ('----------',i=5,ncol)`
		515
		516	`c Write fortran mode (mode=3)`
		517	`elseif (mode.eq.3) then`
		518	`write(fid) ntra,ntim,ncol`
		519	`write(fid) time`
		520	`write(fid) vars`
		521
		522	`c Write IVE netcdf format (mode=4)`
		523	`elseif (mode.eq.4) then`
		524	`vardim(1)=ntra`
		525	`vardim(2)=1`
		526	`vardim(3)=1`
		527	`vardim(4)=1`
		528	`mdv =-999.98999`
		529	`do i=2,ncol`
		530	`call putdef(fid,vars(i),4,mdv,vardim,`
		531	`> varmin,varmax,stag,ierr)`
		532	`enddo`
		533	`varname='BASEDATE'`
		534	`vardim(1)=6`
		535	`call putdef(fid,varname,4,mdv,vardim,`
		536	`> varmin,varmax,stag,ierr)`
		537	`do i=1,6`
		538	`rtime(i)=real(time(i))`
		539	`enddo`
		540	`call putdat(fid,varname,0.,0,rtime,ierr)`
		541
5	michaesp	542	`c Write KML format (mode=5)`
		543	`elseif (mode.eq.5) then`
		544
		545	`write(fid,"(A)") '<?xml version="1.0" encoding="UTF-8"?>'`
		546	`write(fid,"(A)") '<kml xmlns="http://www.opengis.net/kml/2.2">'`
		547	`write(fid,"(A)") '<Document>'`
		548	`write(fid,"(A)") '<name>Paths</name>'`
		549	`write(fid,"(A)") '<Style id="yellowLineGreenPoly">'`
		550	`write(fid,"(A)") '<LineStyle>'`
		551	`c write(fid,*) '<color>7f00ffff</color>' ! Yellow`
		552	`write(fid,"(A)") '<color>500A0A0A</color>' ! Black`
		553	`write(fid,"(A)") '<width>4</width>'`
		554	`write(fid,"(A)") '</LineStyle>'`
		555	`write(fid,"(A)") '<PolyStyle>'`
		556	`write(fid,"(A)") '<color>7f00ff00</color>'`
		557	`write(fid,"(A)") '</PolyStyle>'`
		558	`write(fid,"(A)") '</Style>'`
		559
3	michaesp	560	`endif`
		561
		562	`end`
		563
		564
		565	`c ----------------------------------------------------------------`
		566	`c Close a trajectory file`
		567	`c ----------------------------------------------------------------`
		568
		569	`subroutine close_tra(fid,mode)`
		570
		571	`implicit none`
		572
		573	`c Declaration of subroutine parameters`
		574	`integer fid`
		575	`integer mode`
		576
		577	`c Auxiliary variables`
		578	`integer ierr`
		579
		580	`c Close file`
		581	`if (mode.eq.1) then`
		582	`close(abs(fid))`
5	michaesp	583
3	michaesp	584	`elseif (mode.eq.2) then`
		585	`close(abs(fid))`
5	michaesp	586
3	michaesp	587	`elseif (mode.eq.3) then`
		588	`close(fid)`
5	michaesp	589
3	michaesp	590	`elseif (mode.eq.4) then`
		591	`call clscdf(fid,ierr)`
5	michaesp	592
		593	`elseif (mode.eq.5) then`
		594	`write(fid,"(A)") '</Document>'`
		595	`write(fid,"(A)") '</kml>'`
		596	`close(abs(fid))`
3	michaesp	597	`endif`
		598
		599	`end`
		600
		601
		602	`c ----------------------------------------------------------------`
		603	`c Determine the mode of a trajectory file`
		604	`c ----------------------------------------------------------------`
		605
		606	`subroutine mode_tra(mode,filename)`
		607
		608	`implicit none`
		609
		610	`c Declaration of subroutine parameters`
		611	`integer mode`
		612	`character*80 filename`
		613
		614	`c Auxiliary variables`
		615	`integer len`
5	michaesp	616	`character char0,char1,char2,char3,char4`
3	michaesp	617
		618	`c Get mode`
		619	`mode=-1`
		620
		621	`len = len_trim(filename)`
		622
5	michaesp	623	`c Mode specified by number`
3	michaesp	624	`char0 = filename((len-1):(len-1))`
		625	`char1 = filename(len:len)`
		626
		627	`if ( (char0.eq.'.').and.(char1.eq.'1') ) mode=1`
		628	`if ( (char0.eq.'.').and.(char1.eq.'2') ) mode=2`
		629	`if ( (char0.eq.'.').and.(char1.eq.'3') ) mode=3`
		630	`if ( (char0.eq.'.').and.(char1.eq.'4') ) mode=4`
5	michaesp	631	`if ( (char0.eq.'.').and.(char1.eq.'5') ) mode=5`
3	michaesp	632
5	michaesp	633	`if ( mode.gt.0 ) return`
		634
		635	`c Mode specified by appendix`
		636	`char0 = filename((len-3):(len-3))`
		637	`char1 = filename((len-2):(len-2))`
		638	`char2 = filename((len-1):(len-1))`
		639	`char3 = filename(len:len)`
		640	`if ( (char1.eq.'.').and.(char2.eq.'l').and.(char3.eq.'s') ) mode=1`
		641	`if ( (char1.eq.'.').and.(char2.eq.'t').and.(char3.eq.'i') ) mode=2`
		642	`if ( (char1.eq.'.').and.(char2.eq.'d').and.(char3.eq.'u') ) mode=3`
		643	`if ( (char1.eq.'.').and.(char2.eq.'n').and.(char3.eq.'c') ) mode=4`
		644
		645	`if ( (char0.eq.'.').and.(char1.eq.'k').and.`
		646	`> (char2.eq.'m').and.`
		647	`> (char3.eq.'l') ) mode = 5`
		648
3	michaesp	649	`end`
		650
		651
		652	`c ----------------------------------------------------------------`
		653	`c Get dimension of a trajectory file`
		654	`c ----------------------------------------------------------------`
		655
		656	`subroutine info_tra(filename,ntra,ntim,ncol,mode)`
		657
		658	`implicit none`
		659
		660	`c Declaration of subroutine parameters`
		661	`integer fid`
		662	`character*80 filename`
		663	`integer mode`
		664	`integer ntra,ntim,ncol`
		665
		666	`c Auxiliary variables`
		667	`integer vardim(4)`
		668	`real varmin(4),varmax(4),stag(4)`
		669	`real mdv`
		670	`character*80 cfn`
		671	`integer ierr`
		672	`integer i,ndim`
		673	`character*80 vars(100)`
		674	`integer nvars`
		675	`integer ntimes`
		676	`real times(100)`
		677	`character*500 str`
		678	`integer nline0,nline1,nline2`
		679	`integer isstr,isok`
		680	`character ch`
		681
		682	`c Open file`
		683	`if (mode.eq.1) then`
		684	`fid=10`
		685	`open(fid,file=filename)`
		686	`elseif (mode.eq.2) then`
		687	`fid=10`
		688	`open(fid,file=filename)`
		689	`elseif (mode.eq.3) then`
		690	`fid=10`
		691	`open(fid,file=filename,form='unformatted')`
		692	`elseif (mode.eq.4) then`
		693	`call cdfopn(filename,fid,ierr)`
		694	`endif`
		695
		696	`c Get dimension information`
		697	`if ( (mode.eq.1).or.(mode.eq.2) ) then`
		698	`read(fid,*)`
		699	`read(fid,*)`
		700	`read(fid,'(a500)') str`
		701	`read(fid,*)`
		702
		703	`c Get the number of columns`
		704	`isstr=0`
		705	`ncol =0`
		706	`do i=1,500`
		707	`ch = str(i:i)`
		708
		709	`if ( (isstr.eq.0).and.(ch.ne.' ') ) then`
		710	`isstr=1`
		711	`elseif ( (isstr.eq.1).and.(ch.eq.' ') ) then`
		712	`isstr=0`
		713	`ncol=ncol+1`
		714	`endif`
		715	`enddo`
		716
		717	`c Get the first data block`
		718	`nline0 = 5`
		719	`nline1 = 5`
		720	`read(fid,*)`
		721	`100 read(fid,'(a500)',end=110) str`
		722	`if (str.ne.'') then`
		723	`nline1 = nline1 + 1`
		724	`goto 100`
		725	`endif`
		726	`110 continue`
		727
		728	`c Get the total numbers of lines in the data block`
		729	`nline2 = nline1`
		730	`120 read(fid,*,end=130)`
		731	`nline2 = nline2 + 1`
		732	`goto 120`
		733	`130 nline2 = nline2 + 1`
		734
		735	`c Set the dimensions`
		736	`if (mode.eq.1) then`
		737	`ntim = nline1 - nline0`
		738	`ntra = (nline2-nline0+1)/(ntim+1)`
		739	`else`
		740	`ntra = nline1 - nline0`
		741	`ntim = (nline2-nline0+1)/(ntra+1)`
		742	`endif`
		743
		744	`elseif (mode.eq.3) then`
		745	`read(fid) ntra,ntim,ncol`
		746
		747	`elseif (mode.eq.4) then`
		748	`call gettimes(fid,times,ntimes,ierr)`
		749	`call getvars(fid,nvars,vars,ierr)`
		750	`call getdef(fid,trim(vars(2)),ndim,mdv,vardim,`
		751	`> varmin,varmax,stag,ierr)`
		752	`ntra = vardim(1)`
		753	`ntim = ntimes`
		754	`ncol = nvars-1`
		755
		756	`endif`
		757
		758	`c Close file`
		759	`if (mode.eq.1) then`
		760	`close(fid)`
		761	`elseif (mode.eq.2) then`
		762	`close(fid)`
		763	`elseif (mode.eq.3) then`
		764	`close(fid)`
		765	`elseif (mode.eq.4) then`
		766	`call clscdf(fid,ierr)`
		767	`endif`
		768
5	michaesp	769	`end`
		770
3	michaesp	771
5	michaesp	772	`c ----------------------------------------------------------------`
		773	`c Binary search algorithm`
		774	`c ----------------------------------------------------------------`
		775
		776	`subroutine binary (z,p,ref_z,ref_p)`
		777
		778	`implicit none`
		779
		780	`c Declaration of subroutine parameters`
		781	`real z`
		782	`real p`
		783	`real ref_z(3000)`
		784	`real ref_p(3000)`
		785
		786	`c Auxiliary variables`
		787	`integer i0,i1,im`
		788
		789	`c Binary search`
		790	`i0 = 1`
		791	`i1 = 3000`
		792	`100 continue`
		793	`im = (i0 + i1) / 2`
		794	`if ( p.lt.ref_p(im) ) then`
		795	`i0 = im`
		796	`else`
		797	`i1 = im`
		798	`endif`
		799	`if ( (i1-i0).gt.1 ) goto 100`
		800
		801	`c Linear interpolation in between`
		802	`z = ref_z(i0) + ( p - ref_p(i0) ) / ( ref_p(i1) - ref_p(i0) ) *`
		803	`> ( ref_z(i1) - ref_z(i0) )`
		804
3	michaesp	805	`end`
5	michaesp	806
		807
		808
		809
		810

Subversion Repositories lagranto.ecmwf

(root)/trunk/lib/iotra.f @ 15 – Rev 5