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)`
		58	`endif`
		59
		60	`c Read header information`
		61	`call read_hea(fid,time,vars,ntra,ntim,ncol,mode)`
		62
		63	`end`
		64
		65
		66	`c ----------------------------------------------------------------`
		67	`c Open a trajectory file for wrinting`
		68	`c ----------------------------------------------------------------`
		69
		70	`subroutine wopen_tra(fid,filename,ntra,ntim,ncol,time,vars,mode)`
		71
		72	`implicit none`
		73
		74	`c Declaration of subroutine parameters`
		75	`integer fid`
		76	`character*80 filename`
		77	`integer mode`
		78	`integer ntra,ntim,ncol`
		79	`integer time(6)`
		80	`character*80 vars(ncol)`
		81
		82	`c Auxiliary variables`
		83	`integer vardim(4)`
		84	`real varmin(4),varmax(4),stag(4)`
		85	`real mdv`
		86	`character*80 cfn`
		87	`integer ierr`
		88	`integer i`
		89	`character*80 varname`
		90	`real rtime(6)`
		91
		92	`c Open file`
		93	`if (mode.eq.1) then`
		94	`fid = 10`
		95	`open(fid,file=filename)`
		96	`elseif (mode.eq.2) then`
		97	`fid = 10`
		98	`open(fid,file=filename)`
		99	`elseif (mode.eq.3) then`
		100	`open(fid,file=filename,form='unformatted')`
		101	`elseif (mode.eq.4) then`
		102	`vardim(1)=ntra`
		103	`vardim(2)=1`
		104	`vardim(3)=1`
		105	`vardim(4)=1`
		106	`cfn =trim(filename)//'_cst'`
		107	`mdv =-999.98999`
		108	`call crecdf(filename,fid,varmin,varmax,3,cfn,ierr)`
		109	`endif`
		110
		111	`c Write header information`
		112	`call write_hea(fid,time,vars,ntra,ntim,ncol,mode)`
		113
		114	`end`
		115
		116
		117	`c ----------------------------------------------------------------`
		118	`c Read a trajectory`
		119	`c ----------------------------------------------------------------`
		120
		121	`subroutine read_tra(fid,tra,ntra,ntim,ncol,mode)`
		122
		123	`implicit none`
		124
		125	`c Declaration of subroutine parameters`
		126	`integer fid`
		127	`integer ntim`
		128	`integer ncol`
		129	`integer ntra`
		130	`real tra(ntra,ntim,ncol)`
		131	`integer mode`
		132
		133	`c Auxiliary variables`
		134	`integer i,j,n`
		135	`real arr(ntra)`
		136	`integer ntimes`
		137	`real times(1000)`
		138	`integer ierr`
		139	`character*80 vars(ncol+2)`
		140	`integer nvars`
		141
		142	`c Read ascii mode, sorted by trajectory (mode=1)`
		143	`if (mode.eq.1) then`
		144	`read(fid,*,end=100)`
		145	`do n=1,ntra`
		146	`do i=1,ntim`
		147	`read(fid,*,end=110) (tra(n,i,j),j=1,ncol)`
		148	`enddo`
		149	`enddo`
		150
		151	`c Read ascii mode, sorted by time (mode=2)`
		152	`elseif (mode.eq.2) then`
		153	`read(fid,*,end=100)`
		154	`do i=1,ntim`
		155	`do n=1,ntra`
		156	`read(fid,*,end=100) (tra(n,i,j),j=1,ncol)`
		157	`enddo`
		158	`enddo`
		159
		160	`c Read fortran mode (mode=3)`
		161	`elseif (mode.eq.3) then`
		162	`read(fid) tra`
		163
		164	`c Read IVE netcdf mode (mode=4)`
		165	`elseif (mode.eq.4) then`
		166	`call gettimes(fid,times,ntimes,ierr)`
		167	`call getvars(fid,nvars,vars,ierr)`
		168	`do i=1,ntim`
		169	`do j=1,ncol`
		170	`if (j.eq.1) then`
		171	`do n=1,ntra`
		172	`tra(n,i,1)=times(i)`
		173	`enddo`
		174	`else`
		175	`call getdat(fid,vars(j),times(i),0,arr,ierr)`
		176	`do n=1,ntra`
		177	`tra(n,i,j)=arr(n)`
		178	`enddo`
		179	`endif`
		180	`enddo`
		181	`enddo`
		182
		183	`endif`
		184
		185	`return`
		186
		187	`c End of file has been reached: set negative <fid>`
		188	`100 fid=-fid`
		189	`return`
		190
		191	`c Error: incomplete trajectory`
		192	`110 print*,'<read_tra>: Incomplete trajectory... Stop'`
		193	`stop`
		194
		195	`end`
		196
		197
		198	`c ----------------------------------------------------------------`
		199	`c Write a trajectory`
		200	`c ----------------------------------------------------------------`
		201
		202	`subroutine write_tra(fid,tra,ntra,ntim,ncol,mode)`
		203
		204	`implicit none`
		205
		206	`c Declaration of subroutine parameters`
		207	`integer fid`
		208	`integer ntim`
		209	`integer ncol`
		210	`integer ntra`
		211	`real tra(ntra,ntim,ncol)`
		212	`integer mode`
		213
		214	`c Auxiliary variables`
		215	`integer i,j,n`
		216	`real arr(ntra)`
		217	`integer ierr`
		218	`real time`
		219	`character*80 vars(ncol+2)`
		220	`integer nvars`
		221
		222	`c Write ascii mode, sorted by trajectory (mode=1)`
		223	`if (mode.eq.1) then`
		224	`do n=1,ntra`
		225	`write(fid,*)`
		226	`do i=1,ntim`
		227
		228	`c Avoid ugly *s or missing space in output`
		229	`do j=5,ncol`
		230	`if ( abs(tra(n,i,j)).gt.9999.) then`
		231	`print*,'Format problem : ',tra(n,i,j),' -> -999.99'`
		232	`tra(n,i,j) = -999.99`
		233	`endif`
		234	`enddo`
		235
		236	`write(fid,'(1f7.2,f9.2,f8.2,i6,100f10.3)')`
		237	`> (tra(n,i,j),j=1,3), ! time, lon, lat`
		238	`> nint(tra(n,i,4)), ! p`
		239	`> (tra(n,i,j),j=5,ncol) ! fields`
		240	`enddo`
		241	`enddo`
		242
		243	`c Write ascii mode, sorted by time (mode=2)`
		244	`elseif (mode.eq.2) then`
		245	`do i=1,ntim`
		246	`write(fid,*)`
		247	`do n=1,ntra`
		248
		249	`c Avoid ugly *s or missing space in output`
		250	`do j=5,ncol`
		251	`if ( abs(tra(n,i,j)).gt.9999.) then`
		252	`print*,'Format problem : ',tra(n,i,j),' -> -999.99'`
		253	`tra(n,i,j) = -999.99`
		254	`endif`
		255	`enddo`
		256
		257	`write(fid,'(1f7.2,f9.2,f8.2,i6,100f10.3)')`
		258	`> (tra(n,i,j),j=1,3), ! time, lon, lat`
		259	`> nint(tra(n,i,4)), ! p`
		260	`> (tra(n,i,j),j=5,ncol) ! fields`
		261	`enddo`
		262	`enddo`
		263
		264	`c Write fortran mode (mode=3)`
		265	`elseif (mode.eq.3) then`
		266	`write(fid) tra`
		267
		268	`c Write netcdf mode (mode=4)`
		269	`elseif (mode.eq.4) then`
		270	`call getvars(fid,nvars,vars,ierr)`
		271	`do i=1,ntim`
		272	`time=tra(1,i,1)`
		273	`do j=2,ncol`
		274	`do n=1,ntra`
		275	`arr(n)=tra(n,i,j)`
		276	`enddo`
		277	`call putdat(fid,vars(j),time,0,arr,ierr)`
		278	`enddo`
		279	`enddo`
		280	`endif`
		281
		282	`end`
		283
		284
		285	`c ----------------------------------------------------------------`
		286	`c Read header from trajectory file`
		287	`c ----------------------------------------------------------------`
		288
		289	`subroutine read_hea(fid,time,vars,ntra,ntim,ncol,mode)`
		290
		291	`implicit none`
		292
		293	`c Declaration of subroutine parameters`
		294	`integer fid`
		295	`integer time(6)`
		296	`integer ntra,ntim,ncol`
		297	`character*80 vars(ncol)`
		298	`integer mode`
		299
		300	`c Auxiliary variables`
		301	`integer i`
		302	`character ch(500)`
		303	`character*500 str`
		304	`integer ich(500)`
		305	`integer isstr,ileft,iright`
		306	`character*80 varname`
		307	`real rtime(6)`
		308	`integer ierr`
		309	`integer nvars`
		310	`character*15 str1`
		311	`character str2`
		312	`character*13 str3`
		313	`character*4 str4`
		314	`character*80 linestr`
		315	`integer itmp1,itmp2`
		316
		317	`c Read ascii format (mode=1,2)`
		318	`if ( (mode.eq.1).or.(mode.eq.2) ) then`
		319
		320	`c Read the time specification (old and new format)`
		321	`read(fid,'(a80)') linestr`
		322
		323	`if ( linestr(1:14).eq.'Reference date' ) then`
		324	`read(linestr,'(a15,i4,i2,i2,a1,i2,i2,a13,i8,a4)')`
		325	`> str1,`
		326	`> time(1),time(2),time(3),str2,time(4),time(5),`
		327	`> str3,time(6),str4`
		328
		329	`elseif ( linestr(1:11).eq.'time period' ) then`
		330	`read(linestr,'(a12,i4,i2,i2,a1,i2,a4,i6,a1,i2,a2)')`
		331	`> str1,`
		332	`> time(1),time(2),time(3),str2,time(4),`
		333	`> str3,itmp1,str3,itmp2,str4`
		334	`time(5) = 0`
		335	`time(6) = itmp1 * 60 + itmp2`
		336
		337	`endif`
		338
		339	`c Skip the empty line and read field names`
		340	`read(fid,*)`
		341	`read(fid,'(a500)',end=100) str`
		342	`do i=1,500`
		343	`ch(i)=str(i:i)`
		344	`enddo`
		345
		346	`c Split the input string`
		347	`isstr=0`
		348	`nvars=0`
		349	`do i=1,500`
		350	`if ( (isstr.eq.0).and.(ch(i).ne.' ') ) then`
		351	`isstr=1`
		352	`ileft=i`
		353	`elseif ( (isstr.eq.1).and.(ch(i).eq.' ') ) then`
		354	`isstr=0`
		355	`iright=i-1`
		356	`nvars=nvars+1`
		357	`vars(nvars)=str(ileft:iright)`
		358	`endif`
		359	`enddo`
		360
		361	`c Skip the empty line`
		362	`read(fid,*,end=100)`
		363
		364	`c Read fortran mode (mode=3)`
		365	`elseif (mode.eq.3) then`
		366	`read(fid) ntra,ntim,ncol`
		367	`read(fid) time`
		368	`read(fid) vars`
		369
		370	`c Read IVE netcdf mode (mode=4)`
		371	`elseif (mode.eq.4) then`
		372	`call getvars(fid,nvars,vars,ierr)`
		373	`varname='BASEDATE'`
		374	`call getdat(fid,varname,0.,0,rtime,ierr)`
		375	`do i=1,6`
		376	`time(i)=nint(rtime(i))`
		377	`enddo`
		378	`endif`
		379
		380	`return`
		381
		382	`c End of file has been reached`
		383	`100 fid=-fid`
		384	`return`
		385
		386	`c Excetion handling`
		387	`110 print*,'<read_hea>: Unexspected time format.... Stop'`
		388	`stop`
		389
		390	`end`
		391
		392
		393	`c ----------------------------------------------------------------`
		394	`c Write header to trajectory file (in ascii mode)`
		395	`c ----------------------------------------------------------------`
		396
		397	`subroutine write_hea(fid,time,vars,ntra,ntim,ncol,mode)`
		398
		399	`implicit none`
		400
		401	`c Declaration of subroutine parameters`
		402	`integer fid`
		403	`integer time(6)`
		404	`integer ntra,ntim,ncol`
		405	`character*80 vars(ncol)`
		406	`integer mode`
		407
		408	`c Auxiliary variables`
		409	`integer i`
		410	`character*500 str`
		411	`character*4 str1`
		412	`character*2 str2,str3,str4,str5,str6`
		413	`integer vardim(4)`
		414	`real varmin(4),varmax(4),stag(4)`
		415	`real mdv`
		416	`integer ierr`
		417	`character*80 varname`
		418	`real rtime(6)`
		419	`integer nvars`
		420
		421	`c Write ascii format (mode=1,2)`
		422	`if ( (mode.eq.1).or.(mode.eq.2) ) then`
		423
		424	`c Get the strings for output`
		425	`write(str1,'(i4)') time(1)`
		426	`write(str2,'(i2)') time(2)`
		427	`write(str3,'(i2)') time(3)`
		428	`write(str4,'(i2)') time(4)`
		429	`write(str5,'(i2)') time(5)`
		430	`if (time(2).eq. 0) str2(1:1)='0'`
		431	`if (time(3).eq. 0) str3(1:1)='0'`
		432	`if (time(4).eq. 0) str4(1:1)='0'`
		433	`if (time(5).eq. 0) str5(1:1)='0'`
		434	`if (time(2).lt.10) str2(1:1)='0'`
		435	`if (time(3).lt.10) str3(1:1)='0'`
		436	`if (time(4).lt.10) str4(1:1)='0'`
		437	`if (time(5).lt.10) str5(1:1)='0'`
		438
		439	`c Write the time specification`
		440	`write(fid,'(a15,a4,a2,a2,a1,a2,a2,a13,i8,a4)')`
		441	`> 'Reference date ',`
		442	`> str1,str2,str3,'_',str4,str5,`
		443	`> ' / Time range',time(6), ' min'`
		444	`write(fid,*)`
		445
		446	`c Write variable names`
		447	`str=''`
		448	`do i=1,ncol`
		449	`str=trim(str)//trim(vars(i))`
		450	`enddo`
		451	`write(fid,'(a6,a9,a8,a6,100a10)') (trim(vars(i)),i=1,ncol)`
		452	`write(fid,'(a6,a9,a8,a6,100a10)')`
		453	`> '------','---------','--------','------',`
		454	`> ('----------',i=5,ncol)`
		455
		456	`c Write fortran mode (mode=3)`
		457	`elseif (mode.eq.3) then`
		458	`write(fid) ntra,ntim,ncol`
		459	`write(fid) time`
		460	`write(fid) vars`
		461
		462	`c Write IVE netcdf format (mode=4)`
		463	`elseif (mode.eq.4) then`
		464	`vardim(1)=ntra`
		465	`vardim(2)=1`
		466	`vardim(3)=1`
		467	`vardim(4)=1`
		468	`mdv =-999.98999`
		469	`do i=2,ncol`
		470	`call putdef(fid,vars(i),4,mdv,vardim,`
		471	`> varmin,varmax,stag,ierr)`
		472	`enddo`
		473	`varname='BASEDATE'`
		474	`vardim(1)=6`
		475	`call putdef(fid,varname,4,mdv,vardim,`
		476	`> varmin,varmax,stag,ierr)`
		477	`do i=1,6`
		478	`rtime(i)=real(time(i))`
		479	`enddo`
		480	`call putdat(fid,varname,0.,0,rtime,ierr)`
		481
		482	`endif`
		483
		484	`end`
		485
		486
		487	`c ----------------------------------------------------------------`
		488	`c Close a trajectory file`
		489	`c ----------------------------------------------------------------`
		490
		491	`subroutine close_tra(fid,mode)`
		492
		493	`implicit none`
		494
		495	`c Declaration of subroutine parameters`
		496	`integer fid`
		497	`integer mode`
		498
		499	`c Auxiliary variables`
		500	`integer ierr`
		501
		502	`c Close file`
		503	`if (mode.eq.1) then`
		504	`close(abs(fid))`
		505	`elseif (mode.eq.2) then`
		506	`close(abs(fid))`
		507	`elseif (mode.eq.3) then`
		508	`close(fid)`
		509	`elseif (mode.eq.4) then`
		510	`call clscdf(fid,ierr)`
		511	`endif`
		512
		513	`end`
		514
		515
		516	`c ----------------------------------------------------------------`
		517	`c Determine the mode of a trajectory file`
		518	`c ----------------------------------------------------------------`
		519
		520	`subroutine mode_tra(mode,filename)`
		521
		522	`implicit none`
		523
		524	`c Declaration of subroutine parameters`
		525	`integer mode`
		526	`character*80 filename`
		527
		528	`c Auxiliary variables`
		529	`integer len`
		530	`character char0,char1`
		531
		532	`c Get mode`
		533	`mode=-1`
		534
		535	`len = len_trim(filename)`
		536
		537	`char0 = filename((len-1):(len-1))`
		538	`char1 = filename(len:len)`
		539
		540	`if ( (char0.eq.'.').and.(char1.eq.'1') ) mode=1`
		541	`if ( (char0.eq.'.').and.(char1.eq.'2') ) mode=2`
		542	`if ( (char0.eq.'.').and.(char1.eq.'3') ) mode=3`
		543	`if ( (char0.eq.'.').and.(char1.eq.'4') ) mode=4`
		544
		545	`end`
		546
		547
		548	`c ----------------------------------------------------------------`
		549	`c Get dimension of a trajectory file`
		550	`c ----------------------------------------------------------------`
		551
		552	`subroutine info_tra(filename,ntra,ntim,ncol,mode)`
		553
		554	`implicit none`
		555
		556	`c Declaration of subroutine parameters`
		557	`integer fid`
		558	`character*80 filename`
		559	`integer mode`
		560	`integer ntra,ntim,ncol`
		561
		562	`c Auxiliary variables`
		563	`integer vardim(4)`
		564	`real varmin(4),varmax(4),stag(4)`
		565	`real mdv`
		566	`character*80 cfn`
		567	`integer ierr`
		568	`integer i,ndim`
		569	`character*80 vars(100)`
		570	`integer nvars`
		571	`integer ntimes`
		572	`real times(100)`
		573	`character*500 str`
		574	`integer nline0,nline1,nline2`
		575	`integer isstr,isok`
		576	`character ch`
		577
		578	`c Open file`
		579	`if (mode.eq.1) then`
		580	`fid=10`
		581	`open(fid,file=filename)`
		582	`elseif (mode.eq.2) then`
		583	`fid=10`
		584	`open(fid,file=filename)`
		585	`elseif (mode.eq.3) then`
		586	`fid=10`
		587	`open(fid,file=filename,form='unformatted')`
		588	`elseif (mode.eq.4) then`
		589	`call cdfopn(filename,fid,ierr)`
		590	`endif`
		591
		592	`c Get dimension information`
		593	`if ( (mode.eq.1).or.(mode.eq.2) ) then`
		594	`read(fid,*)`
		595	`read(fid,*)`
		596	`read(fid,'(a500)') str`
		597	`read(fid,*)`
		598
		599	`c Get the number of columns`
		600	`isstr=0`
		601	`ncol =0`
		602	`do i=1,500`
		603	`ch = str(i:i)`
		604
		605	`if ( (isstr.eq.0).and.(ch.ne.' ') ) then`
		606	`isstr=1`
		607	`elseif ( (isstr.eq.1).and.(ch.eq.' ') ) then`
		608	`isstr=0`
		609	`ncol=ncol+1`
		610	`endif`
		611	`enddo`
		612
		613	`c Get the first data block`
		614	`nline0 = 5`
		615	`nline1 = 5`
		616	`read(fid,*)`
		617	`100 read(fid,'(a500)',end=110) str`
		618	`if (str.ne.'') then`
		619	`nline1 = nline1 + 1`
		620	`goto 100`
		621	`endif`
		622	`110 continue`
		623
		624	`c Get the total numbers of lines in the data block`
		625	`nline2 = nline1`
		626	`120 read(fid,*,end=130)`
		627	`nline2 = nline2 + 1`
		628	`goto 120`
		629	`130 nline2 = nline2 + 1`
		630
		631	`c Set the dimensions`
		632	`if (mode.eq.1) then`
		633	`ntim = nline1 - nline0`
		634	`ntra = (nline2-nline0+1)/(ntim+1)`
		635	`else`
		636	`ntra = nline1 - nline0`
		637	`ntim = (nline2-nline0+1)/(ntra+1)`
		638	`endif`
		639
		640	`elseif (mode.eq.3) then`
		641	`read(fid) ntra,ntim,ncol`
		642
		643	`elseif (mode.eq.4) then`
		644	`call gettimes(fid,times,ntimes,ierr)`
		645	`call getvars(fid,nvars,vars,ierr)`
		646	`call getdef(fid,trim(vars(2)),ndim,mdv,vardim,`
		647	`> varmin,varmax,stag,ierr)`
		648	`ntra = vardim(1)`
		649	`ntim = ntimes`
		650	`ncol = nvars-1`
		651
		652	`endif`
		653
		654	`c Close file`
		655	`if (mode.eq.1) then`
		656	`close(fid)`
		657	`elseif (mode.eq.2) then`
		658	`close(fid)`
		659	`elseif (mode.eq.3) then`
		660	`close(fid)`
		661	`elseif (mode.eq.4) then`
		662	`call clscdf(fid,ierr)`
		663	`endif`
		664
		665
		666	`end`
		667

Subversion Repositories lagranto.ecmwf

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