program ocupa c---------------------------------------------------------------------- double precision timdif c CHARACTER NSTA*10 dimension ab(6,200),nsta(200),dist(200),az(200),azin(200) double precision az,dist common/dados/ab,nd,nsta,dist,az,azin,vp,vs c integer iocupa(21,24) c character*12 kar,cdata1*10,cdata2*10,ctime1*8,ctime2*8 integer*2 zstime(6),zetime(6),zatime(6),z1time(6),z2time(6) integer dias(12),diasb(12),diasn(12) c data diasn/0,31,59,90,120,151,181,212,243,273,304,334/ data diasb/0,31,60,91,121,152,182,213,244,274,305,335/ c...............................................First campaign day zatime(1)=02 zatime(2)=331 zatime(3)=0 zatime(4)=0 zatime(5)=0 zatime(6)=0 vp=6.0 vs=3.46 ab(2,1)=0.0 ab(3,1)=0.0 ab(5,1)=0.0 ab(6,1)=0.0 ab(2,2)=0.0 ab(3,2)=0.0 ab(5,2)=0.0 ab(6,2)=0.0 c do 5 ij=1,21 do 4 ih=1,24 iocupa(ij,ih)=0 4 continue 5 continue c open(unit=9,file='all-mcs.txt', status='old') open(unit=10,file='all-mcs.log', status='unknown') open(unit=11,file='all-mcs.ocu', status='unknown') c...............................................................leitura sumhh=0. sumdd=0. do 100 ii=1,100000 read(9,'(a)',end=999,err=999) kar write(*,*) kar read(9,1000,end=999,err=999) rlon1,rlat1,cdata1,ctime1 1000 format(f11.7,f12.7,t42,a10,t53,a8) read(cdata1,1001) iaa,imes,idd 1001 format(i4,1x,i2,1x,i2) iaa=iaa-2000 izz1=iaa/4 izz2=izz1*4 if (iaa.eq.izz2) then c write(*,*) ' *** Ano Bissexto ***',kano do 543 i=1,12 543 dias(i)=diasb(i) else do 654 i=1,12 654 dias(i)=diasn(i) endif ijul=dias(imes)+idd zstime(1)=iaa zstime(2)=ijul read(ctime1,1002) ihh,imm,iss 1002 format(i2,1x,i2,1x,i2) zstime(3)=ihh zstime(4)=imm zstime(5)=iss zstime(6)=00 read(9,1000,end=999,err=999) rlon2,rlat2,cdata2,ctime2 read(cdata2,1001) iaa,imes,idd iaa=iaa-2000 ijul=dias(imes)+idd zetime(1)=iaa zetime(2)=ijul read(ctime2,1002) ihh,imm,iss zetime(3)=ihh zetime(4)=imm zetime(5)=iss zetime(6)=00 strtim = sngl(timdif(zetime,zstime)) delhh=strtim/3600. sumhh=sumhh+delhh ab(1,1)=rlat1 ab(4,1)=rlon1 ab(1,2)=rlat2 ab(4,2)=rlon2 nd=2 call subdis() deldd=sngl(dist(2)) sumdd=sumdd+deldd write(10,*) kar,deldd,delhh c tbeg=sngl(timdif(zstime,zatime))/3600. tend=sngl(timdif(zetime,zatime))/3600. write(*,*) tbeg,tend do 50 ij=1,21 do 40 ih=1,24 if (iocupa(ij,ih).eq.1) goto 40 z1time(1)=02 z1time(2)=330+ij z1time(3)=ih-1 z1time(4)=0 z1time(5)=0 z1time(6)=0 t1=sngl(timdif(z1time,zatime))/3600. z2time(1)=02 z2time(2)=330+ij z2time(3)=ih z2time(4)=0 z2time(5)=0 z2time(6)=0 t2=sngl(timdif(z2time,zatime))/3600. if (t2.gt.tbeg .and. t1.lt.tend) iocupa(ij,ih)=1 40 continue 50 continue c 100 continue c 999 write(10,*) sumdd,sumhh c iplot=0 do 200 il=0,6 iplot=0 do 180 ic=1,72 ii=(ic-1)/24 ij=1+il*3+ii ih=ic-ii*24 c write(*,*) ij,ih c.....................................Dentro if (iocupa(ij,ih).eq.1) then if (iplot.eq.0) then x1=float(ic-1) x2=float(ic) iplot=1 else x2=float(ic) endif c.....................................Fora else if (iplot.eq.0) then else y1=float(il) y2=float(il)+0.20 write(11,2001) 2001 format('>') write(11,*) x1,y1 write(11,*) x1,y2 write(11,*) x2,y2 write(11,*) x2,y1 write(11,*) x1,y1 iplot=0 endif endif 180 continue if (iplot.eq.1) then y1=float(il) y2=float(il)+0.20 write(11,2001) write(11,*) x1,y1 write(11,*) x1,y2 write(11,*) x2,y2 write(11,*) x2,y1 write(11,*) x1,y1 iplot=0 endif 200 continue c stop end c c****************************************************** c double precision function timdif(time1,time2) c----------------------------------------------------------------------- c c Returns the time difference in seconds between two PC-DASI format c time arrays: timdif = time1 - time2 c c----------------------------------------------------------------------- implicit double precision (s,t) logical timset integer*2 time1(6), time2(6) timdif = -12345.0d0 c c Verify that times are defined c if (.not. timset(time1) .or. .not. timset(time2)) goto 20 if (time1(1) .ne. time2(1)) goto 30 c c Compute time difference c spd = 86400.0d0 sph = 3600.0d0 spm = 60.0d0 spms = dble(1.0/1000.0) tmp1 = spd*dble(time1(2)) + sph*dble(time1(3)) + & spm*dble(time1(4)) + dble(time1(5)) + & spms*dble(time1(6)) tmp2 = spd*dble(time2(2)) + sph*dble(time2(3)) + & spm*dble(time2(4)) + dble(time2(5)) + & spms*dble(time2(6)) timdif = tmp1 - tmp2 return c 20 write(*,*) 'TIMDIF undefined: times not completely defined!' return c 30 write(*,*) 'TIMDIF undefined: years are different.' return c end c c****************************************************** c logical function timset(time) c----------------------------------------------------------------------- c c .true. if entire time array is defined. c c----------------------------------------------------------------------- integer*2 time(6) c timset = .false. do 10 i = 1, 6 10 if (time(i) .lt. 0) return timset = .true. c return end c c************************************************************ c subroutine subdis() c INTEGER RR,SS,PD,RRR,S,TD,ALFA,CON CHARACTER NSTA*10,IFILEN*9,NEPI*9,NSTA1*10,idtiro*4 CHARACTER*1 ICAR1,fout*20 dimension ab(6,200),nsta(200),dist(200),az(200),azin(200) DIMENSION BL(200),BP(200),ALTI(200),ICW(200) DIMENSION TH(200),PHI(200),XDEG(200) DIMENSION DELTA(200),DELTE(200),T4(200),T5(200) DOUBLE PRECISION THK,PHK,SINTK,COSTK,TANTK,TH2,V1,THG,D,E,F,A, %B,C,THC,PHC,TANTC,TH2I,V2,AL,DL,A12,AZ,COSA12,DM,EO,EO2,EO3, %EOP,EOQ,U1,BO,CO,C2,C4,Z1,Z2,X2,Y2,U2,DIST,D1,E1,F1,C1,A1,B1, %SC,SD,XDEG,TEST,DELTA,T4,T5,T2,DELTE,RADI,RAD,EC,ELL,EL,EP, %BL,EEL,EEP,BP,TH,PHI,PAXIS,EAXIS common/dados/ab,nd,nsta,dist,az,azin,vp,vs c DO 113 I=1,ND BP(I)=AB(1,I)+AB(2,I)/60.+AB(3,I)/3600. ICW(I)=1 if(AB(6,I).LT.0.) ICW(I)=-1 if(AB(6,I).LT.0.) AB(6,I)=-AB(6,I) if(AB(5,I).LT.0.) ICW(I)=-1 IF(AB(5,I).LT.0.) AB(5,I)=-AB(5,I) IF(AB(4,I).LT.0) ICW(I)=-1 IF(ICW(I).LT.0) GO TO 232 BL(I)=AB(4,I)+AB(5,I)/60.+AB(6,I)/3600. GO TO 113 232 BL(I)=AB(4,I)-AB(5,I)/60.-AB(6,I)/3600. 113 CONTINUE C 107 DO 108 I=1,ND TH(I)=BP(I) PHI(I)=BL(I) 108 CONTINUE c EEL=BL(1) EEP=BP(1) NEPI=NSTA(1) c RADI=57.2957795 RAD=1.0/RADI ELL=.99327733 EC=.672267002E-02 C ELL=.993277 C EC=.672267E-02 EL=EC/ELL EP=1.0+EL PAXIS=6356.912 EAXIS=6378.388 THK=RAD*EEP PHK=RAD*EEL SINTK= DSIN(THK) COSTK= DCOS(THK) TANTK=SINTK/COSTK TH2=THK+THK V1=DEXP(.230258509E+01*(.380544268E+01-(.7323684E-3)*DCOS(TH2) 1+(.6175 E-6)*DCOS(TH2+TH2)-(.7E-9)*DCOS(TH2+TH2+TH2))) C V1=EXP(.230259E+01*(.380544E+01-(.732368E-3)*COS(TH2) C 1+(.6175 E-6)*COS(TH2+TH2)-(.7E-9)*COS(TH2+TH2+TH2))) THG=DATAN(ELL*TANTK) D=DSIN(PHK) E=-DCOS(PHK) F=-DCOS(THG) A=F*E B=-D*F C=DSIN(THG) DO 38 I=2,ND THC=RAD*TH(I) PHC=RAD*PHI(I) TANTC=DSIN(THC)/DCOS(THC) TH2I= THC+THC V2=DEXP(.230258509E+01*(.380544268E+01-(.7323684E-3)*DCOS(TH2I) 1+(.6175 E-6)*DCOS(TH2I+TH2I)-(.7E-9)*DCOS(TH2I+TH2I+TH2I))) C V2=EXP(.230259E+01*(.380544E+01-(.732368E-3)*COS(TH2I) C 1+(.6175 E-6)*COS(TH2I+TH2I)-(.7E-9)*COS(TH2I+TH2I+TH2I))) AL= TANTC/(EP*TANTK)+EC*DSQRT((EP+TANTC*TANTC)/(EP+TANTK*TANTK)) DL= PHC-PHK A12=DATAN2(DSIN(DL),((AL-DCOS(DL))*SINTK)) AZ(I)=RADI*A12 IF(AZ(I).LT.0.0) AZ(I)=AZ(I)+360. COSA12=DCOS(A12) DM=COSTK*COSA12 DM=DM*DM+SINTK*SINTK EO= EL*DM EO2=EO*EO EO3=EO*EO2 EOP=1.+EO EOQ= DSQRT(EOP) U1=DATAN2(SINTK,(COSTK*EOQ*COSA12)) BO=V1*DSQRT(1.+EL*COSTK*COSTK*COSA12*COSA12)/EOP CO=1.+0.25*EO-(3.*EO2-1.25*EO3)/64. C2=-0.125*EO+EO2/32.-15.0*EO3/1024. C4=(-EO2+0.75*EO3)/256. Z1=V1*(1.-EC)*SINTK Z2=V2*(1.-EC)*DSIN(THC) X2=V2*DCOS(THC)*DCOS(DL) Y2=V2*DCOS(THC)*DSIN(DL) U2=DATAN((V1*SINTK+EOP*(Z2-Z1))/(EOQ*(X2*COSA12-Y2*SINTK* & DSIN(A12)))) DIST(I)=BO*(CO*(U2-U1)+C2*(DSIN(U2+U2)-DSIN(U1+U1)) & +C4*(DSIN(U2+U2+U2+U2)-DSIN(U1+U1+U1+U1))) DIST(I)=DABS(DIST(I)) THG= DATAN(ELL*TANTC) D1=DSIN(PHC) D1=DSIN(PHC) E1=-DCOS(PHC) F1=-DCOS(THG) C1=DSIN(THG) A1=F1*E1 B1=-D1*F1 SC=A*A1+B*B1+C*C1 SD=DSQRT(0.25*((A-A1)*(A-A1)+(B-B1)*(B-B1)+(C-C1)*(C-C1))* & ((A+A1)*(A+A1)+(B+B1)*(B+B1)+(C+C1)*(C+C1))) XDEG(I)=RADI*DATAN2(SD,SC) TEST=DIST(I)-111.0*XDEG(I) IF(DABS(TEST)-100.)24,23,23 23 U2=U2+0.314159265E+01 DIST(I)=BO*(CO*(U2-U1)+C2*(DSIN(U2+U2)-DSIN(U1+U1)) & +C4*(DSIN(U2+U2+U2+U2)-DSIN(U1+U1+U1+U1))) 24 DIST(I)=DABS(DIST(I)) DELTA(I)=DIST(I)/vp T4(I)=T3+DELTA(I) T5(I)=T2 202 IF(T4(I).GT.60.)T5(I)=T5(I)+1. IF(T4(I).GT.60.)T4(I)=T4(I)-60. IF(T4(I).GT.60.)GO TO 202 DELTE(I)=DIST(I)/vs 25 AL=TANTK/(EP*TANTC)+EC*DSQRT((EP+TANTK*TANTK)/(EP+TANTC*TANTC)) AZIN(I)=RADI*DATAN2(DSIN(-DL),((AL-DCOS(DL))*DSIN(THC))) IF(AZIN(I).LT.0.0)AZIN(I)=AZIN(I)+360.0 38 CONTINUE c c return end