c c------------------------------------------------------------------ c Begin file /home/sting/foster/loc/doptuv.f c------------------------------------------------------------------ c SUBROUTINE my_DOPTUV(LUOUT,UTC,IT,IX,JX,LX,fields) C C Calculate doppler temp and winds (uses Gonzalo's program REDLINE) C C ON INPUT: C LUOUT = fortran unit on which to write output info (not written if <=0 C UTC = current ut, IT = index of current ut C IX JX LX = indices for current lon, lat, loc C C ON OUTPUT: C DOPTN DOPUN DOPVN = Doppler temp, zonal wind, meridional wind C VER = volume emission rate (photons cm-3 s-1) C RAY = brightness in Rayleighs (integration of VER) C include 'tigcmloc.h' include 'doppler.h' c dimension fields(imx,kmx,jmx,nfget) DIMENSION HT(kmx),TN(kmx),UN(kmx),VN(kmx),TE(kmx),XNO2(kmx), + XNO(kmx),XNN2(kmx),XNO2P(kmx),XNE(kmx),VERK(kmx), + HTCM(kmx),sr63(kmx),xn2d(kmx),qo2(kmx),f6300(kmx) CHARACTER*40 TOPLAB DATA B1D/1.2/, RK3/8.E-12/, A1D/0.0093/, A6300/0.0071/ data isr63/1/, ifrst/0/ ! if isr63=1 -> add sr63 from solred.f to ver real epeak(2), eflx(2) save ifrst C ifrst = ifrst+1 call getloc(fields,ht,jx,ix,0,ixz) call getloc(fields,tn,jx,ix,0,ixt) call getloc(fields,un,jx,ix,0,ixu) call getloc(fields,vn,jx,ix,0,ixv) call getloc(fields,te,jx,ix,0,ixte) call getloc(fields,xno2,jx,ix,0,ixo2) call getloc(fields,xno,jx,ix,0,ixo1) call getloc(fields,xnn2,jx,ix,0,ixn2) call getloc(fields,xno2p,jx,ix,0,ixo2p) call getloc(fields,xn2d,jx,ix,0,ixn2d) call getloc(fields,xne,jx,ix,0,ixne) qo2(:) = xno2p(:) htcm(:) = ht(:) * 1.e+5 c write(6,"('doptuv: ix=',i4,' jx=',i4)") ix,jx c write(6,"('doptuv: ixz=',i2,' ifget(ixz)=',i2,' ht (1st,last)=', c + 2e12.4)") ixz,ifget(ixz),ht(1),ht(kmx) c write(6,"('doptuv: ixt=',i2,' ifget(ixt)=',i2,' tn (1st,last)=', c + 2e12.4)") ixt,ifget(ixt),tn(1),tn(kmx) c write(6,"('doptuv: ixu=',i2,' ifget(ixu)=',i2,' un (1st,last)=', c + 2e12.4)") ixu,ifget(ixu),un(1),un(kmx) c write(6,"('doptuv: ixv=',i2,' ifget(ixv)=',i2,' vn (1st,last)=', c + 2e12.4)") ixv,ifget(ixv),vn(1),vn(kmx) c write(6,"('doptuv: ixte=',i2,' ifget(ixte)=',i2,' te (1st,last)=', c + 2e12.4)") ixte,ifget(ixte),te(1),te(kmx) c write(6,"('doptuv: ixo2=',i2,' ifget(ixo2)=',i2,' o2 (1st,last)=', c + 2e12.4)") ixo2,ifget(ixo2),xno2(1),xno2(kmx) c write(6,"('doptuv: ixo1=',i2,' ifget(ixo1)=',i2,' o (1st,last)=', c + 2e12.4)") ixo1,ifget(ixo1),xno(1),xno(kmx) c write(6,"('doptuv: ixn2=',i2,' ifget(ixn2)=',i2,' n2 (1st,last)=', c + 2e12.4)") ixn2,ifget(ixn2),xnn2(1),xnn2(kmx) c write(6,"('doptuv: ixo2p=',i2,' ifget(ixo2p)=',i2, c + ' o2p (1st,last)=',2e12.4)") ixo2p,ifget(ixo2p),xno2p(1), c + xno2p(kmx) c write(6,"('doptuv: ixn2d=',i2,' ifget(ixn2d)=',i2, c + ' n2d (1st,last)=',2e12.4)") ixn2d,ifget(ixn2d),xn2d(1), c + xn2d(kmx) c write(6,"('doptuv: ixne=',i2,' ifget(ixne)=',i2, c + ' ne (1st,last)=',2e12.4)") ixne,ifget(ixne),xne(1),xne(kmx) c C WRITE(6,"(' DOPTUV: XNO2=',/(5E12.4))") XNO2 C WRITE(6,"(' DOPTUV: XNO=',/(5E12.4))") XNO C WRITE(6,"(' DOPTUV: XNN2=',/(5E12.4))") XNN2 * do 100 k=1,kmx * RK1 = 2.E-11 * EXP(107.8/TN(K)) * RK2 = 2.9E-11 * EXP(67.5/TN(K)) * IF (TE(K).GE.1200.) RKDR = 1.6E-7 * (TE(K)/300.)**(-0.55) * IF (TE(K).LT.1200.) RKDR = 1.95E-7 * (TE(K)/300.)**(-0.7) * VER(IT,K,LX) = (B1D * RKDR * XNO2P(K) * XNE(K)) / * + (A1D + RK1*XNN2(K) + RK2*XNO2(K) + RK3*XNO(K)) * A6300 * 100 CONTINUE * Modifications that Barbara made to the TIGCM for the * soft aurora. This is an awkward way of dealing with * auroral precipitation, but without major changes, this * is the only way possible.... if(utc.lt.3.67 .or. utc.gt.5.42) then eshift=10. eflx(1)=10. epeak(1)=1500. eflx(2)=3. epeak(2)=150. if(utc.gt.6.) then eflx(1)=(8.-utc)*10./2. eflx(2)=(8.-utc)*3./2. end if if(utc.gt.8.) then eflx(1)=0. eflx(2)=0. end if else if(utc.lt.4.08) then eshift=10. epeak(1)=200. eflx(1)=24.*(utc-3.67)/0.41+13. epeak(2)=0. eflx(2)=0. else if(utc.gt.4.42) then eshift=10. epeak(1)=200. eflx(1)=24.*(5.42-utc)+13. epeak(2)=0. eflx(2)=0. else eshift=10. epeak(1)=200. eflx(1)=37. epeak(2)=0. eflx(2)=0. end if eflx=0. ! for background case (no aurora) * if(utc.le.8.00) then * eshift=10. * eflx(1)=10. * epeak(1)=1500. * end if print*,' Calling AURORA at ',utc,' UT with eflx=',eflx call aurora(6,utc,'crs','rdt',eflx,epeak,eshift,49.,kmx,htcm, . xnn2,xno2,xno,xn2d,qo2,xne,tn,te,f6300) do k=1,kmx ver(it,k,lx)=f6300(k) end do c c 2/7/91: Add sr63 to volume emmission rate: c c SUBROUTINE SOLRED(IFRST,STL,zpht,xno,xno2,tn,SR63,f107,f107a, c + day,glat,glong) c if (isr63.gt.0) then slt = utc + gcmlon(ix)/15. if (slt.lt.0.) slt = 24.+slt day = iyd(1) - iyd(1)/1000*1000 call solred(ifrst,slt,ht,xno,xno2,xnn2,tn,sr63, + f107ms,f107ams,day,gcmlat(jx),gcmlon(ix)) ver(it,:,lx) = ver(it,:,lx) + sr63(:) endif VERK(:) = VER(IT,:,LX) C C Integrate volume emmission rate to find brightness in Rayleigh's: C RAY(IT,LX) = QUADRAT(kmx,HTCM,VERK) * 1.E-6 C C Call Gonzalo's routine to find doppler tn,uv,vn: C CALL REDLINE(LUOUT,kmx,HT,TN,UN,VN,VERK, VELV,VELZ,ATTPV,DATTPV, + ATTPZ,DATTPZ) C WRITE(6,"(' AFTER REDLINE: UT=',F6.2,' VELV VELZ=',2F10.3, C + ' ATTPV DATTPV=',2F10.2,' ATTPZ,DATTPZ=',2F10.2)") C + UTC,VELV,VELZ,ATTPV,DATTPV,ATTPZ,DATTPZ DOPTN(IT,LX) = ATTPV DOPUN(IT,LX) = VELZ DOPVN(IT,LX) = VELV C RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE REDLINE(LUOUT,NP,HM,TM,ZM,VM,VER,VELV,VELZ,ATTPV,DATTPV + ,ATTPZ,DATTPZ) DIMENSION HT(25), TT(25), 1ZT(25), VT(25), BT(25), DL(25), DOP(25), 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/, 1SIGMA/15867.852/, AMASS/16.0/, JMAX/61/ R = 6370.0 ALPHA = 20.0 C NP = 25 C OPEN (UNIT=1, FILE = 'ray_data',STATUS = 'OLD' ) C OPEN (UNIT=2, STATUS = 'UNKNOWN' ) C C WRITE(6,"(' ENTER REDLINE: NP=',I3,' HM=',/(5E12.4))") NP, C + (HM(I),I=1,NP) C WRITE(6,"(' TM=',/(5E12.4))") TM C WRITE(6,"(' ZM=',/(5E12.4))") ZM C WRITE(6,"(' VM=',/(5E12.4))") VM C WRITE(6,"(' VER=',/(5E12.4))") VER C PI180 = PI/180.0 B = ( BB * SIGMA)**2/(ALOG(2.0) * AMASS) NPN = NP -1 4 CONTINUE C READ(1,100,END=10)( HM(I), TM(I), ZM(I), VM(I), VER(I), I =1, NP ) 100 FORMAT ( /, (11X, F6.0,3X,4E12.3 )) 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 IF (LUOUT.GT.0) + WRITE(LUOUT,102)(HM(I),TM(I),ZM(I),VM(I),VER(I), I=1,NP) IF (LUOUT.GT.0) + WRITE(LUOUT,102)(HT(I),TT(I),ZT(I),VT(I),BT(I), I=1,NPN) 102 FORMAT(/,(11X,F6.0,3X,3F12.2,F12.7)) 1 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 IF (LUOUT.GT.0) + WRITE(LUOUT,102)(HT(I), TT(I), BT(I),DL(I),DOP(I), I= 1,NPN) ABT = SUMBT/SUMB ABZ = SUMBZ/SUMB ABV = SUMBV/SUMB ABL = SUMBL/SUMB ABH = SUMBH/SUMB ABZH = SUMZ/SUMZB AHL = R*(SQRT(1.0 + (ABL/R)**2 + 2.0*ABL*RSINA/R**2) - 1.0) ABZC = ABZ*(R+ABZH)/RCOSA ABVC = ABV*(R+ABZH)/RCOSA SUMB = SUMB/10. SUMZB = SUMZB/10. IF (LUOUT.GT.0) +WRITE(LUOUT,101) ABZH,ABH,AHL,ABL,ABT,ABZ,ABZC,ABV,ABVC,SUMB,SUMZB 101 FORMAT(/,1X, F10.2,2X, 'ZENITH HEIGHT',/ 1 ,1X, F10.2,2X, 'SLANT HEIGHT' ,/ 2 ,1X, F10.2,2X, 'MEAN SLANT HEIGHT' ,/ 3 ,1X, F10.2,2X, 'MEAN SLANT LENGTH' ,/ 4 ,1X, F10.2,2X, 'SLANT TEMPERATURE' ,/ 5 ,1X, F10.2,2X, 'ZONAL LINE-OF-SIGHT' ,/ 6 ,1X, F10.2,2X, 'ZONAL WIND CORRECTED' ,/ 7 ,1X, F10.2,2X, 'MERIDIONAL LINE OF SIGHT' ,/ 8 ,1X, F10.2,2X, 'MERIDIONAL CORRECTED' ,/ 9 ,1X, F10.2,2X, 'SLANT EMISSION RATE' ,/ A ,1X, F10.2,2X, 'ZENITH EMISSION RATE' ,/) 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 POSV = POSZ IF (LUOUT.GT.0) + WRITE(LUOUT,107)(AS(I), AV(I), AZ(I), I = 1, JMAX) 107 FORMAT(/,10X, ' AS, AV, AND AZ',/,(10X,E14.4,2F14.4)) 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) POSV = AS(I) IV = I 14 IF ( AZ(I) .LE. AVZ ) GO TO 15 AVZ = AZ(I) POSZ = AS(I) IZ = I 15 CONTINUE IF (LUOUT.GT.0) + WRITE(LUOUT,107)(AS(I), AV(I), AZ(I), I = 1, JMAX) 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 IF (LUOUT.GT.0) + WRITE(LUOUT,107)(ASW(I), AWS(I), WE(I), I = 1, 10) 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 IF (LUOUT.GT.0) + WRITE(LUOUT,107)(ASW(I), AWS(I), WE(I), I = 1, 10) 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 IF (LUOUT.GT.0) + WRITE(LUOUT,103) VELV,VELZ 103 FORMAT( 1X, F10.2, 2X, 'CORRECTED MERIDIONAL VELOCITY', 1 /,1X, F10.2, 2X, 'CORRECTED ZONAL VELOCITY ' ) 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 IF (LUOUT.GT.0) + WRITE(LUOUT,104) ATTPV, DATTPV, ATTPZ, DATTPZ 104 FORMAT( 1X,2F10.2, 2X, 'CORRECT MERIDIONAL TEMPERATURE', 1 /,1X,2F10.2, 2X, 'CORRECT ZONAL TEMPERATURE') C GO TO 4 10 CONTINUE END 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 SX2 = 0.0 SY2 = 0.0 SYX2 = 0.0 R2 = 0.0 SB = 0.0 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 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 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC REAL FUNCTION QUADRAT(NPTS,HTS,PARM) DIMENSION HTS(NPTS),PARM(NPTS) QUADRAT = 0. DO 100 K=1,NPTS-1 100 QUADRAT = QUADRAT + 0.5*(PARM(K)+PARM(K+1)) * (HTS(K+1)-HTS(K)) RETURN END