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

Subversion Repositories lagranto.ecmwf

(root)/trunk/lib/iotra.f – Rev 41