c c 3/10/94: c Taken from doptuv.f (this dir), and changed name from redline to doppler. c This routine takes ht,t,u,v, and ver as input (HM,TM,ZM,VM,VER, c dimensioned np (which is usually kmx in calling routines)), and c returns scalars doppler v,u,t (VELV,VELZ,ATTPV) c Note input VER is not necessarilly from 6300A (was assumed to be in c doptuv.f) c Keeping in quadrat at end for use in integrating vol emission rate (ver) c to obtain brightness (not used in doppler routine) c Note: for redline emissions pass in sigma=15867.852 c for greenline emissions pass sigma=17924. c Example call: c call doppler(kmx,z,t,u,v,emis,sigma,dopv,dopu,dopt, c + dattpv,attpz,dattpz,spval) c CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE doppler(NP,HM,TM,ZM,VM,VER,sigma, + VELV,VELZ,ATTPV,DATTPV,ATTPZ,DATTPZ,spval) DIMENSION HT(NP), TT(NP), 1ZT(NP), VT(NP), BT(NP), DL(NP), DOP(NP), AS(61), AV(65), AZ(61), 2AZC(61), AVC(61), WE(61), AWS(61),ASW(61), WEV(61), WEZ(61) DIMENSION HM(NP),TM(NP),ZM(NP),VM(NP),VER(NP) DATA PI/3.141592654/, BB/3.5811516E-7/, CSPEED/3.0E10/, 1 AMASS/16.0/, JMAX/61/ c R = 6370.0 ALPHA = 20.0 PI180 = PI/180.0 B = ( BB * SIGMA)**2/(ALOG(2.0) * AMASS) NPN = NP -1 DO 2 I = 1 , NPN II = I + 1 HT(I) = (HM(I) + HM(II))/2.0 TT(I) = (TM(I) + TM(II))/2.0 ZT(I) = (ZM(I) + ZM(II))/2.0 VT(I) = (VM(I) + VM(II))/2.0 BT(I) = (VER(I) + VER(II))/2.0 2 CONTINUE c 1 RADALPH = ALPHA * PI180 RADALPH = ALPHA * PI180 RSINA = R * SIN(RADALPH) RSINASQ = RSINA**2 RCOSA = R * COS(RADALPH) TWOR = 2.0 * R SUMB = 0.0 SUMBT = 0.0 SUMBZ = 0.0 SUMBV = 0.0 SUMBL = 0.0 SUMBH = 0.0 SUMZ = 0.0 SUMZB = 0.0 DO 3 I = 1 ,NPN II = I + 1 DH = HM(II) HH = HM(I) DPH = SQRT( RSINASQ + TWOR*DH + DH**2) DOH = SQRT( RSINASQ + TWOR*HH + HH**2) DLM = (DPH + DOH)/2.0 - RSINA DL(I) = DPH - DOH DOP(I) = RCOSA/(R + HT(I)) BD = BT(I) * DL (I) SUMB = SUMB + BD SUMBT = SUMBT + BD*TT(I) SUMBZ = SUMBZ + BD*DOP(I)*ZT(I) SUMBV = SUMBV + BD*DOP(I)*VT(I) SUMBL = SUMBL + BD*DLM SUMBH = SUMBH + BD*HT(I) BD = ( HM(II) - HM(I)) * BT(I) SUMZ = SUMZ + BD*HT(I) SUMZB = SUMZB + BD 3 CONTINUE ABT = SUMBT/SUMB ABZ = SUMBZ/SUMB ABV = SUMBV/SUMB ABL = SUMBL/SUMB c ABH = SUMBH/SUMB ABZH = SUMZ/SUMZB c AHL = R*(SQRT(1.0 + (ABL/R)**2 + 2.0*ABL*RSINA/R**2) - 1.0) c ABZC = ABZ*(R+ABZH)/RCOSA c ABVC = ABV*(R+ABZH)/RCOSA SUMB = SUMB/10. SUMZB = SUMZB/10. SIGMAV = SIGMA * ( ABV/CSPEED) SIGMAZ = SIGMA * ( ABZ/CSPEED) DDG = SQRT(B*ABT*ALOG(2.0))/10.0 DGS = DDG*30.0 JC = JMAX/2 DGSS = -DGS DO 11 I = 1, JMAX AV(I) = 0.0 AZ(I) = 0.0 AS(I) = DGSS DGSS = DGSS + DDG 11 CONTINUE BSUM = 0.0 DO 12 K = 1, NPN TBT = B * TT(K) OBT = -1.0/TBT ATOP = SQRT(-OBT/PI)* BT(K) * DL(K) BSUM = BSUM + ATOP SIGMAHZ = SIGMA * (ZT(K) * DOP(K)/CSPEED) SIGMAHV = SIGMA * (VT(K) * DOP(K)/CSPEED) SIGMA0Z = SIGMAHZ - SIGMAZ SIGMA0V = SIGMAHV - SIGMAV DGSS = DGS J = JMAX + 1 DO 13 I = 1, JC JI = J - I SIGV1 = EXP(OBT*( -DGSS - SIGMA0V)**2) SIGV2 = EXP(OBT*( DGSS - SIGMA0V)**2) SIGZ1 = EXP(OBT*( -DGSS - SIGMA0Z)**2) SIGZ2 = EXP(OBT*( DGSS - SIGMA0Z)**2) AV(I) = AV(I) + ATOP * SIGV1 AV(JI) = AV(JI) + ATOP * SIGV2 AZ(I) = AZ(I) + ATOP * SIGZ1 AZ(JI) = AZ(JI) + ATOP * SIGZ2 DGSS = DGSS - DDG 13 CONTINUE AV( JC + 1 ) = AV( JC + 1 ) + ATOP * EXP(OBT*(SIGMA0V**2)) AZ( JC + 1 ) = AZ( JC + 1 ) + ATOP * EXP(OBT*(SIGMA0Z**2)) 12 CONTINUE AVZ = -1.0E30 AVV = AVZ IV = JC + 1 IZ = IV POSZ = 0.0 c POSV = POSZ DO 15 I = 1, JMAX WEV(I) = AV(I)*0.129 AV(I) = ALOG(WEV(I)) WEZ(I) = AZ(I)*0.129 AZ(I) = ALOG(WEZ(I)) WE(I) = 1.0 IF(AV(I) .LE. AVV) GO TO 14 AVV = AV(I) c POSV = AS(I) IV = I 14 IF ( AZ(I) .LE. AVZ ) GO TO 15 AVZ = AZ(I) POSZ = AS(I) IZ = I 15 CONTINUE IVS = IV - 5 DDG2 = DDG/2.0 DO 19 I = 1, 10 AWS(I) = ( AV(IVS + 1 ) - AV(IVS))/DDG ASW(I) = AS(IVS) + DDG2 IVS = IVS + 1 19 CONTINUE CALL ONEPAR(ASW,AWS,WE,10,CROV,CRB,SLOPE,SB,CR1,CRB1,SL1,SB1,R2) CROV = -CROV/SLOPE VELV = (( SIGMAV+CROV)*CSPEED/SIGMA)*(R+ABZH)/RCOSA IZS = IZ - 5 DO 20 I = 1, 10 AWS(I) = ( AZ(IZS + 1 ) - AZ(IZS))/DDG ASW(I) = AS(IZS) + DDG2 IZS = IZS + 1 20 CONTINUE CALL ONEPAR(ASW,AWS,WE,10,CROZ,CRB,SLOPE,SB,CR1,CRB1,SL1,SB1,R2) CROZ = -CROZ/SLOPE VELZ = (( SIGMAZ+CROZ)*CSPEED/SIGMA)*(R+ABZH)/RCOSA DO 18 I = 1, JMAX AZC(I) = (AS(I) + CROZ)**2 AVC(I) = (AS(I) + CROV)**2 18 CONTINUE CALL ONEPAR(AZC,AZ,WEZ,JMAX,CR,CRB,SLOPE,SB,CR1,CRB1,SL1,SB1,R2) ATTPZ = -1.0/(SLOPE * B ) DATTPZ = B*SB*ATTPZ**2 CALL ONEPAR(AVC,AV,WEV,JMAX,CR,CRB,SLOPE,SB,CR1,CRB1,SL1,SB1,R2) ATTPV = -1.0/(SLOPE * B ) DATTPV = B*SB*ATTPV**2 c 10 CONTINUE END c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c SUBROUTINE ONEPAR(X,Y,N,NP,CR,CRB,SL,SB,CR1,CRB1,SL1,SB1,R2) DIMENSION X (NP), Y (NP), N (NP) REAL N SUMX = 0.0 SUMXQ = 0.0 SUMY = 0.0 SUMYQ = 0.0 SUMXY = 0.0 SUMN = 0.0 SUMX2 = 0.0 SUMY2 = 0.0 SUMXTY = 0.0 SL = 0.0 CR = 0.0 c SX2 = 0.0 c SY2 = 0.0 c SYX2 = 0.0 R2 = 0.0 SB = 0.0 c T = 0.0 DO 1 I = 1, NP SUMX = SUMX + X (I) * N (I) SUMY = SUMY + Y (I) * N (I) SUMXY = SUMXY + X (I) * Y (I) * N (I) SUMXQ = SUMXQ + X (I) * * 2 * N (I) SUMYQ = SUMYQ + Y (I) * * 2 * N (I) SUMN = SUMN + N (I) 1 CONTINUE SUMXY = SUMXY * SUMN SUMXQ = SUMXQ * SUMN SUMXTY = SUMX * SUMY SUMYQ = SUMYQ * SUMN SUMX2 = SUMX * * 2 SUMY2 = SUMY * * 2 XL = SUMXQ - SUMX2 YL = SUMYQ - SUMY2 IF(ABS(XL).LT.1.0E-25.OR.ABS(YL).LT.1.0E-25)RETURN SL = (SUMXY - SUMXTY) / XL SB = SQRT (ABS (SUMN / XL)) SL1 = (SUMXY - SUMXTY) / YL SB1 = SQRT( ABS(SUMN/YL)) XL = XL * SUMN YL = YL * SUMN c SYX2 = SQRT (ABS (SUMX2 / XL)) CR = (SUMXQ * SUMY - SUMX * SUMXY) / XL CRB = SUMXQ/XL CR1 = (SUMYQ*SUMX - SUMY*SUMXY)/ YL CRB1 = SUMYQ/YL CRB1 = SQRT(ABS(CRB1)) CRB = SQRT( CRB) R2 = SQRT (ABS (SL * SL1)) RETURN END