program ampl5_grads PARAMETER (NF=5,m=2,ng=(NF-1)/2,mm=2) parameter (ntime=248,ilat=36,kgit=10) implicit real*4(a-h,o-z) real amph(ilat,kgit,mm),amplt(ntime,ilat,kgit,mm), $ AMPAVE(kgit),PHSAVE(kgit) real ghr(ntime),ghrinp(ntime),ghi(ntime),ghiinp(ntime), $ tHour(ntime) real ghrDiff(ntime), ghiDiff(ntime) real FUN(NF,ntime),SQ(NF),A(NF,NF),AM1(NF,NF),X(NF),Y(NF) real Xstw(ilat,kgit), Ystw(ilat,kgit) real SIGMA(NF), WORK(2*NF), errx(nf),erry(nf),ierr real period(nf), amplr(nf), ampli(nf) real Aw(ng,ilat,kgit),Ae(ng,ilat,kgit) real Pw(ng,ilat,kgit),Pe(ng,ilat,kgit) real apwe(ng+1,ilat,kgit,4) pi = 4. * atan(1.0) dolg=0. c------------------------------------ M=1 Period1 = 24. Period2 = 12. c Period3 = 8. ! c Period4 = 6. ! c---------------------------------------- om1 = 2. * pi / Period1 om2 = 2. * pi / Period2 c om3 = 2. * pi / Period3 c om4 = 2. * pi / Period4 Period(1)=period1 Period(2)=period2 c Period(3)=period3 c Period(4)=period4 c---- read data open (99, FILE='ERR.dat') open (11, FILE='lh_m2_La.dx' , $ form='unformatted', * access='direct', status='old',recl=4*ILAT*kgit*mm) open (16, FILE='lh_La_m2.dx' , $ form='unformatted', * access='direct', status='unknown',recl=4*(ng+1)*ILAT*kgit*4) print*,'ilat=',ilat,' kgit=',kgit,' ntime=',ntime ntime1=1 ntime2=ntime1+247 shour=0. do 112 jt=ntime1,ntime2 read(11, rec=jt) amph do kz=1,kgit do ii=1,ilat amplt(jt,ii,kz,1)=amph(ii,kz,1) amplt(jt,ii,kz,2)=amph(ii,kz,2) end do end do tHour(jt-ntime1+1) = float(jt-ntime1)*3. shour=shour+tHour(jt-ntime1+1) 112 continue hourav=shour/float(ntime2-ntime1) do 100 kz=1,kgit SAVE=0. CAVE=0. do 200 ii=1,ilat c....................................calculating Cos and Sin components do jt=ntime1,ntime2 ghrinp(jt-ntime1+1)=amplt(jt,ii,kz,1)*cos(float(m)* $ (dolg-amplt(jt,ii,kz,2))*pi/180.) ghiinp(jt-ntime1+1)=amplt(jt,ii,kz,1)*sin(float(m)* $ (dolg-amplt(jt,ii,kz,2))*pi/180.) c if(kz.eq.10.and.ii.eq.25) then c write(10,110)tHour(jt-ntime1+1), c $ ghrinp(jt-ntime1+1),ghiinp(jt-ntime1+1) c end if 110 format(1x,3e13.4) end do c-------------------------------------------------------------------- j0=1 j1=1 DO j=1,ntime2-ntime1 ghr(j)=ghrinp(j+j1-1) ghi(j)=ghiinp(j+j1-1) FUN(1,j)=1. FUN(2,j)=sin(om1 * (tHour(j+j1-1)-tHour(j1))) FUN(3,j)=COS(om1 * (tHour(j+j1-1)-tHour(j1))) FUN(4,j)=sin(om2 * (tHour(j+j1-1)-tHour(j1))) FUN(5,j)=COS(om2 * (tHour(j+j1-1)-tHour(j1))) c FUN(6,j)=sin(om3 * (tHour(j+j1-1)-tHour(j1))) c FUN(7,j)=COS(om3 * (tHour(j+j1-1)-tHour(j1))) c FUN(8,j)=sin(om4 * (tHour(j+j1-1)-tHour(j1))) c FUN(9,j)=COS(om4 * (tHour(j+j1-1)-tHour(j1))) ENDDO DO I=1,NF S=0. DO J=1,ntime2-ntime1 S=S+FUN(I,J)*FUN(I,J) ENDDO SQ(I)=SQRT(S) DO J=1,ntime2-ntime1 FUN(I,J)=FUN(I,J)/SQ(I) ENDDO ENDDO DO I=1,NF DO K=1,NF S=0. DO J=1,ntime2-ntime1 S=S+FUN(I,J)*FUN(K,J) ENDDO AM1(I,K)=S ENDDO ENDDO tol=1.e-5 call LGINF(AM1,nf,nf,nf,tol,A,nf,SIGMA,WORK,ierr) DO I=1,NF S=0. DO J=1,ntime2-ntime1 S=S+ghr(J)*FUN(I,J) ENDDO WORK(I)=S ENDDO call PRMV(nf,nf,nf,A,WORK,X) DO I=1,NF X(I)=X(I)/SQ(I) if(I.eq.1) Xstw(ii,kz)=X(I) ENDDO DO I=1,NF S=0. DO J=1,ntime2-ntime1 S=S+ghi(J)*FUN(I,J) ENDDO WORK(I)=S ENDDO call PRMV(nf,nf,nf,A,WORK,Y) DO I=1,NF Y(I)=Y(I)/SQ(I) if(I.eq.1) Ystw(ii,kz)=Y(I) ENDDO C--------------- CALCULATIONS OF ERRORS IN COEFFICIENTS DO I=1,NF SSu=0. SSv=0. DO J=1,ntime2-ntime1 Su=0. Sv=0. DO L=1,NF Su=Su+A(I,L)*FUN(L,J)*0.1 Sv=Sv+A(I,L)*FUN(L,J)*0.1 END DO SSu=SSu+Su*Su SSv=SSv+Sv*Sv END DO ERRX(I)= SQRT(SSu)/SQ(I) ERRY(I)= SQRT(SSv)/SQ(I) END DO do 300 i=1,ng is=2+2*(i-1) ic=3+2*(i-1) ys=y(is) yc=y(ic) xs=x(is) xc=x(ic) amplI(i)=sqrt(ys*ys+yc*yc) amplR(i)=sqrt(xs*xs+xc*xc) ERRrS=ERRX(IS) ERRrC=ERRX(IC) ERRiS=ERRY(IS) ERRiC=ERRY(IC) ERRi=SQRT(ERRiS*ERRiS*YS*YS+ERRiC*ERRiC*YC*YC)/AMPLI(i) ERRr=SQRT(ERRrS*ERRrS*XS*XS+ERRrC*ERRrC*XC*XC)/AMPLR(i) Frequ= 1./period(i) phaseI=180.*ATAN2(ys,yc)/pi if(phaseI.lt.0.)phaseI=phaseI+360. phaseR=180.*ATAN2(xs,xc)/pi if(phaseR.lt.0.)phaseR=phaseR+360. c---------------------------------------------------------------- phaseR=phaseR+(thour(j1)-thour(j0))*360./period(i) phaseI=phaseI+(thour(j1)-thour(j0))*360./period(i) 601 continue if(phaser.gt.360.)then phaser=phaser-360. Go to 601 End if 602 continue if(phasei.gt.360.)then phasei=phasei-360. Go to 602 End if errphR=180./(xs*xs+xc*xc)*SQRT(xc*xc*errx(is)*errx(is)+ $ xs*xs*errx(ic)*errx(ic))/pi errphI=180./(ys*ys+yc*yc)*SQRT(yc*yc*erry(is)*erry(is)+ $ ys*ys*erry(ic)*erry(ic))/pi c------------------------- CALCULATION OF THE WESTWARD AND EASTWARD c PROPAGATING WAVES Aws= 0.5*(amplR(i)*cos(pi/2.-phaseR*pi/180.)- $ amplI(i)*cos( phaseI*pi/180.)) Awc= 0.5*(amplR(i)*cos( phaseR*pi/180.)+ $ amplI(i)*cos(pi/2.-phaseI*pi/180.)) Aes=-0.5*(amplR(i)*cos(pi/2.-phaseR*pi/180.)+ $ amplI(i)*cos( phaseI*pi/180.)) Aec= 0.5*(amplR(i)*cos( phaseR*pi/180.)- $ amplI(i)*cos(pi/2.-phaseI*pi/180.)) Awest=SQRT(Aws*Aws+Awc*Awc) Aeast=SQRT(Aes*Aes+Aec*Aec) phaseW=180.*ATAN2(Aws,Awc)/pi if(phaseW.lt.0.)phaseW=phaseW+360. PhaseW=PhaseW/float(m) phaseE=180.*ATAN2(Aes,Aec)/pi if(phaseE.lt.0.)phaseE=phaseE+360. PhaseE=PhaseE/float(m) Aw(i,ii,kz)=Awest Ae(i,ii,kz)=Aeast Pw(i,ii,kz)= phaseW Pe(i,ii,kz)= phaseE apwe(i,ii,kz,1)=Aw(i,ii,kz) apwe(i,ii,kz,2)=Ae(i,ii,kz) apwe(i,ii,kz,3)=Pw(i,ii,kz) apwe(i,ii,kz,4)=Pe(i,ii,kz) c---- real and imaginary parts for 2-nd harmonic above the equator IF(I.NE.2) GO TO 400 IF(II.LT.18.OR.II.GT.19) GO TO 400 SAVE = SAVE+Aws/2. CAVE = CAVE+Awc/2. 400 CONTINUE c--------------- end harmonics 300 continue ampl= sqrt(Xstw(ii,kz)*Xstw(ii,kz) & +Ystw(ii,kz)*Ystw(ii,kz)) phas= 360.-180./pi & *ATAN2(Ystw(ii,kz),Xstw(ii,kz)) if( phas.gt.360) phas = phas -360. phas = phas/float(m) apwe(ng+1,ii,kz,1) = ampl apwe(ng+1,ii,kz,2) = 0. apwe(ng+1,ii,kz,3) = phas apwe(ng+1,ii,kz,4) = 0. jout=ntime2-ntime1 DO j=1,Jout tDout=thour(j0+j-1) ghrAPPROX=0. ghiAPPROX=0. DO I=1,nf ghrApprox = ghrAPPROX+X(I)*SQ(I)*FUN(I,J) ghiApprox = ghiAPPROX+Y(I)*SQ(I)*FUN(I,J) ENDDO ghrDiff(j) = ghr(j) - ghrApprox ghiDiff(j) = ghi(j) - ghiApprox Enddo c-------------- end latitude 200 continue AMPAVE(kz) = SQRT(SAVE*SAVE+CAVE*CAVE) PHSAVE(kz) = 180./pi*ATAN2(SAVE,CAVE) IF(PHSAVE(kz).LT.0.) PHSAVE(kz) = PHSAVE(kz)+360. PHSAVE(kz) = PHSAVE(kz)/float(m) zoben =135. dz = zoben/(float(kgit)-0.5) z = dz*0.5+float(kz-1)*dz c Write(10,110)z,AMPAVE(kz),PHSAVE(kz) c-------------- end altitude 100 continue print*,'AW8=',Aw(1,20,20),' Pw8=',Pw(1,20,20) print*,'AWtdw=',Aw(2,20,20),' Pwtdw=',Pw(2,20,20) print*,'X(1)=',X(1),' Y(1)=',Y(1) write(16,rec=1) apwe end