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

Rev	Author	Line No.	Line
2	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	`c Sebastian: Why 9999 -> conflicts with Pressure if it is in Pa`
		230	`do j=5,ncol`
		231	`c if ( abs(tra(n,i,j)).gt.9999.) then`
		232	`if ( abs(tra(n,i,j)).gt.999999.) then`
		233	`print*,'Format problem : ',tra(n,i,j),' -> -999.99'`
		234	`tra(n,i,j) = -999.99`
		235	`endif`
		236	`enddo`
		237
21	michaesp	238	`write(fid,'(1f7.2,f10.3,f9.3,i6,100f10.3)')`
2	michaesp	239	`> (tra(n,i,j),j=1,3), ! time, lon, lat`
		240	`> nint(tra(n,i,4)), ! p`
		241	`> (tra(n,i,j),j=5,ncol) ! fields`
		242	`enddo`
		243	`enddo`
		244
		245	`c Write ascii mode, sorted by time (mode=2)`
		246	`elseif (mode.eq.2) then`
		247	`do i=1,ntim`
		248	`write(fid,*)`
		249	`do n=1,ntra`
		250
		251	`c Avoid ugly *s or missing space in output`
		252	`do j=5,ncol`
		253	`if ( abs(tra(n,i,j)).gt.9999.) then`
		254	`print*,'Format problem : ',tra(n,i,j),' -> -999.99'`
		255	`tra(n,i,j) = -999.99`
		256	`endif`
		257	`enddo`
		258
21	michaesp	259	`write(fid,'(1f7.2,f10.3,f9.3,i6,100f10.3)')`
2	michaesp	260	`> (tra(n,i,j),j=1,3), ! time, lon, lat`
		261	`> nint(tra(n,i,4)), ! p`
		262	`> (tra(n,i,j),j=5,ncol) ! fields`
		263	`enddo`
		264	`enddo`
		265
		266	`c Write fortran mode (mode=3)`
		267	`elseif (mode.eq.3) then`
		268	`write(fid) tra`
		269
		270	`c Write netcdf mode (mode=4)`
		271	`elseif (mode.eq.4) then`
		272	`call getvars(fid,nvars,vars,ierr)`
		273	`do i=1,ntim`
		274	`time=tra(1,i,1)`
		275	`do j=2,ncol`
		276	`do n=1,ntra`
		277	`arr(n)=tra(n,i,j)`
		278	`enddo`
		279	`call putdat(fid,vars(j),time,0,arr,ierr)`
		280	`enddo`
		281	`enddo`
		282	`endif`
		283
		284	`end`
		285
		286
		287	`c ----------------------------------------------------------------`
		288	`c Read header from trajectory file`
		289	`c ----------------------------------------------------------------`
		290
		291	`subroutine read_hea(fid,time,vars,ntra,ntim,ncol,mode)`
		292
		293	`implicit none`
		294
		295	`c Declaration of subroutine parameters`
		296	`integer fid`
		297	`integer time(6)`
		298	`integer ntra,ntim,ncol`
		299	`character*80 vars(ncol)`
		300	`integer mode`
		301
		302	`c Auxiliary variables`
		303	`integer i`
		304	`character ch(500)`
		305	`character*500 str`
		306	`integer ich(500)`
		307	`integer isstr,ileft,iright`
		308	`character*80 varname`
		309	`real rtime(6)`
		310	`integer ierr`
		311	`integer nvars`
		312	`character*15 str1`
		313	`character str2`
		314	`character*13 str3`
		315	`character*4 str4`
		316	`character*80 linestr`
		317	`integer itmp1,itmp2`
		318
		319	`c Read ascii format (mode=1,2)`
		320	`if ( (mode.eq.1).or.(mode.eq.2) ) then`
		321
		322	`c Read the time specification (old and new format)`
		323	`read(fid,'(a80)') linestr`
11	michaesp	324
2	michaesp	325	`if ( linestr(1:14).eq.'Reference date' ) then`
		326	`read(linestr,'(a15,i4,i2,i2,a1,i2,i2,a13,i8,a4)')`
		327	`> str1,`
		328	`> time(1),time(2),time(3),str2,time(4),time(5),`
		329	`> str3,time(6),str4`
		330
		331	`elseif ( linestr(1:11).eq.'time period' ) then`
		332	`read(linestr,'(a12,i4,i2,i2,a1,i2,a4,i6,a1,i2,a2)')`
		333	`> str1,`
		334	`> time(1),time(2),time(3),str2,time(4),`
		335	`> str3,itmp1,str3,itmp2,str4`
		336	`time(5) = 0`
		337	`time(6) = itmp1 * 60 + itmp2`
		338
		339	`endif`
		340
		341	`c Skip the empty line and read field names`
		342	`read(fid,*)`
		343	`read(fid,'(a500)',end=100) str`
		344	`do i=1,500`
		345	`ch(i)=str(i:i)`
		346	`enddo`
		347
		348	`c Split the input string`
		349	`isstr=0`
		350	`nvars=0`
		351	`do i=1,500`
		352	`if ( (isstr.eq.0).and.(ch(i).ne.' ') ) then`
		353	`isstr=1`
		354	`ileft=i`
		355	`elseif ( (isstr.eq.1).and.(ch(i).eq.' ') ) then`
		356	`isstr=0`
		357	`iright=i-1`
		358	`nvars=nvars+1`
		359	`vars(nvars)=str(ileft:iright)`
		360	`endif`
		361	`enddo`
		362
		363	`c Skip the empty line`
		364	`read(fid,*,end=100)`
		365
		366	`c Read fortran mode (mode=3)`
		367	`elseif (mode.eq.3) then`
		368	`read(fid) ntra,ntim,ncol`
		369	`read(fid) time`
		370	`read(fid) vars`
		371
		372	`c Read IVE netcdf mode (mode=4)`
		373	`elseif (mode.eq.4) then`
		374	`call getvars(fid,nvars,vars,ierr)`
		375	`varname='BASEDATE'`
		376	`call getdat(fid,varname,0.,0,rtime,ierr)`
		377	`do i=1,6`
		378	`time(i)=nint(rtime(i))`
		379	`enddo`
		380	`endif`
		381
		382	`return`
		383
		384	`c End of file has been reached`
		385	`100 fid=-fid`
		386	`return`
		387
		388	`c Excetion handling`
		389	`110 print*,'<read_hea>: Unexspected time format.... Stop'`
		390	`stop`
		391
		392	`end`
		393
		394
		395	`c ----------------------------------------------------------------`
		396	`c Write header to trajectory file (in ascii mode)`
		397	`c ----------------------------------------------------------------`
		398
		399	`subroutine write_hea(fid,time,vars,ntra,ntim,ncol,mode)`
		400
		401	`implicit none`
		402
		403	`c Declaration of subroutine parameters`
		404	`integer fid`
		405	`integer time(6)`
		406	`integer ntra,ntim,ncol`
		407	`character*80 vars(ncol)`
		408	`integer mode`
		409
		410	`c Auxiliary variables`
		411	`integer i`
		412	`character*500 str`
		413	`character*4 str1`
		414	`character*2 str2,str3,str4,str5,str6`
		415	`integer vardim(4)`
		416	`real varmin(4),varmax(4),stag(4)`
		417	`real mdv`
		418	`integer ierr`
		419	`character*80 varname`
		420	`real rtime(6)`
		421	`integer nvars`
		422
		423	`c Write ascii format (mode=1,2)`
		424	`if ( (mode.eq.1).or.(mode.eq.2) ) then`
		425
		426	`c Get the strings for output`
		427	`write(str1,'(i4)') time(1)`
		428	`write(str2,'(i2)') time(2)`
		429	`write(str3,'(i2)') time(3)`
		430	`write(str4,'(i2)') time(4)`
		431	`write(str5,'(i2)') time(5)`
		432	`if (time(2).eq. 0) str2(1:1)='0'`
		433	`if (time(3).eq. 0) str3(1:1)='0'`
		434	`if (time(4).eq. 0) str4(1:1)='0'`
		435	`if (time(5).eq. 0) str5(1:1)='0'`
		436	`if (time(2).lt.10) str2(1:1)='0'`
		437	`if (time(3).lt.10) str3(1:1)='0'`
		438	`if (time(4).lt.10) str4(1:1)='0'`
		439	`if (time(5).lt.10) str5(1:1)='0'`
		440
		441	`c Write the time specification`
		442	`write(fid,'(a15,a4,a2,a2,a1,a2,a2,a13,i8,a4)')`
		443	`> 'Reference date ',`
		444	`> str1,str2,str3,'_',str4,str5,`
		445	`> ' / Time range',time(6), ' min'`
		446	`write(fid,*)`
		447
		448	`c Write variable names`
		449	`str=''`
		450	`do i=1,ncol`
		451	`str=trim(str)//trim(vars(i))`
		452	`enddo`
21	michaesp	453	`write(fid,'(a6,a10,a9,a6,100a10)') (trim(vars(i)),i=1,ncol)`
		454	`write(fid,'(a6,a10,a9,a6,100a10)')`
		455	`> '------','----------','---------','------',`
2	michaesp	456	`> ('----------',i=5,ncol)`
		457
		458	`c Write fortran mode (mode=3)`
		459	`elseif (mode.eq.3) then`
		460	`write(fid) ntra,ntim,ncol`
		461	`write(fid) time`
		462	`write(fid) vars`
		463
		464	`c Write IVE netcdf format (mode=4)`
		465	`elseif (mode.eq.4) then`
		466	`vardim(1)=ntra`
		467	`vardim(2)=1`
		468	`vardim(3)=1`
		469	`vardim(4)=1`
		470	`mdv =-999.98999`
		471	`do i=2,ncol`
		472	`call putdef(fid,vars(i),4,mdv,vardim,`
		473	`> varmin,varmax,stag,ierr)`
		474	`enddo`
		475	`varname='BASEDATE'`
		476	`vardim(1)=6`
		477	`call putdef(fid,varname,4,mdv,vardim,`
		478	`> varmin,varmax,stag,ierr)`
		479	`do i=1,6`
		480	`rtime(i)=real(time(i))`
		481	`enddo`
		482	`call putdat(fid,varname,0.,0,rtime,ierr)`
		483
		484	`endif`
		485
		486	`end`
		487
		488
		489	`c ----------------------------------------------------------------`
		490	`c Close a trajectory file`
		491	`c ----------------------------------------------------------------`
		492
		493	`subroutine close_tra(fid,mode)`
		494
		495	`implicit none`
		496
		497	`c Declaration of subroutine parameters`
		498	`integer fid`
		499	`integer mode`
		500
		501	`c Auxiliary variables`
		502	`integer ierr`
		503
		504	`c Close file`
		505	`if (mode.eq.1) then`
		506	`close(abs(fid))`
		507	`elseif (mode.eq.2) then`
		508	`close(abs(fid))`
		509	`elseif (mode.eq.3) then`
		510	`close(fid)`
		511	`elseif (mode.eq.4) then`
		512	`call clscdf(fid,ierr)`
		513	`endif`
		514
		515	`end`
		516
		517
		518	`c ----------------------------------------------------------------`
		519	`c Determine the mode of a trajectory file`
		520	`c ----------------------------------------------------------------`
		521
		522	`subroutine mode_tra(mode,filename)`
		523
		524	`implicit none`
		525
		526	`c Declaration of subroutine parameters`
		527	`integer mode`
		528	`character*80 filename`
		529
		530	`c Auxiliary variables`
		531	`integer len`
		532	`character char0,char1`
		533
		534	`c Get mode`
		535	`mode=-1`
		536
		537	`len = len_trim(filename)`
		538
		539	`char0 = filename((len-1):(len-1))`
		540	`char1 = filename(len:len)`
		541
		542	`if ( (char0.eq.'.').and.(char1.eq.'1') ) mode=1`
		543	`if ( (char0.eq.'.').and.(char1.eq.'2') ) mode=2`
		544	`if ( (char0.eq.'.').and.(char1.eq.'3') ) mode=3`
		545	`if ( (char0.eq.'.').and.(char1.eq.'4') ) mode=4`
		546
		547	`end`
		548
		549
		550	`c ----------------------------------------------------------------`
		551	`c Get dimension of a trajectory file`
		552	`c ----------------------------------------------------------------`
		553
		554	`subroutine info_tra(filename,ntra,ntim,ncol,mode)`
		555
		556	`implicit none`
		557
		558	`c Declaration of subroutine parameters`
		559	`integer fid`
		560	`character*80 filename`
		561	`integer mode`
		562	`integer ntra,ntim,ncol`
		563
		564	`c Auxiliary variables`
		565	`integer vardim(4)`
		566	`real varmin(4),varmax(4),stag(4)`
		567	`real mdv`
		568	`character*80 cfn`
		569	`integer ierr`
		570	`integer i,ndim`
		571	`character*80 vars(100)`
		572	`integer nvars`
		573	`integer ntimes`
		574	`real times(100)`
		575	`character*500 str`
		576	`integer nline0,nline1,nline2`
		577	`integer isstr,isok`
		578	`character ch`
		579
		580	`c Open file`
		581	`if (mode.eq.1) then`
		582	`fid=10`
		583	`open(fid,file=filename)`
		584	`elseif (mode.eq.2) then`
		585	`fid=10`
		586	`open(fid,file=filename)`
		587	`elseif (mode.eq.3) then`
		588	`fid=10`
		589	`open(fid,file=filename,form='unformatted')`
		590	`elseif (mode.eq.4) then`
		591	`call cdfopn(filename,fid,ierr)`
		592	`endif`
		593
		594	`c Get dimension information`
		595	`if ( (mode.eq.1).or.(mode.eq.2) ) then`
		596	`read(fid,*)`
		597	`read(fid,*)`
		598	`read(fid,'(a500)') str`
		599	`read(fid,*)`
		600
		601	`c Get the number of columns`
		602	`isstr=0`
		603	`ncol =0`
		604	`do i=1,500`
		605	`ch = str(i:i)`
		606
		607	`if ( (isstr.eq.0).and.(ch.ne.' ') ) then`
		608	`isstr=1`
		609	`elseif ( (isstr.eq.1).and.(ch.eq.' ') ) then`
		610	`isstr=0`
		611	`ncol=ncol+1`
		612	`endif`
		613	`enddo`
11	michaesp	614	`c backspace(fid)`
2	michaesp	615	`c Get the first data block`
		616	`nline0 = 5`
		617	`nline1 = 5`
		618	`read(fid,*)`
		619	`100 read(fid,'(a500)',end=110) str`
		620	`if (str.ne.'') then`
		621	`nline1 = nline1 + 1`
		622	`goto 100`
		623	`endif`
		624	`110 continue`
11	michaesp	625	`backspace(fid)`
2	michaesp	626
		627	`c Get the total numbers of lines in the data block`
		628	`nline2 = nline1`
		629	`120 read(fid,*,end=130)`
		630	`nline2 = nline2 + 1`
		631	`goto 120`
		632	`130 nline2 = nline2 + 1`
		633
		634	`c Set the dimensions`
		635	`if (mode.eq.1) then`
		636	`ntim = nline1 - nline0`
		637	`ntra = (nline2-nline0+1)/(ntim+1)`
		638	`else`
		639	`ntra = nline1 - nline0`
		640	`ntim = (nline2-nline0+1)/(ntra+1)`
		641	`endif`
		642
		643	`elseif (mode.eq.3) then`
		644	`read(fid) ntra,ntim,ncol`
		645
		646	`elseif (mode.eq.4) then`
		647	`call gettimes(fid,times,ntimes,ierr)`
		648	`call getvars(fid,nvars,vars,ierr)`
		649	`call getdef(fid,trim(vars(2)),ndim,mdv,vardim,`
		650	`> varmin,varmax,stag,ierr)`
		651	`ntra = vardim(1)`
		652	`ntim = ntimes`
		653	`ncol = nvars-1`
		654
		655	`endif`
		656
		657	`c Close file`
		658	`if (mode.eq.1) then`
		659	`close(fid)`
		660	`elseif (mode.eq.2) then`
		661	`close(fid)`
		662	`elseif (mode.eq.3) then`
		663	`close(fid)`
		664	`elseif (mode.eq.4) then`
		665	`call clscdf(fid,ierr)`
		666	`endif`
		667
		668
		669	`end`
		670

Subversion Repositories lagranto.wrf

(root)/trunk/lib/iotra.f @ 24 – Rev 21