C+ C NAME: C Extractd C PURPOSE: C Extracts interpolated velocity and density IPS deconvolved maps at C a given time and position. C CALLING SEQUENCE: C call Extractd(bExtract,RR,dRR,XCbe,XCtbeg,XCtend,nLng,nLat,nMap,nT,FALLOFF, C & VV,DD,XCshift,Rconst,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar) C INPUTS: C bExtract(12) logical .TRUE.: extract; .FALSE.: do not extract C 1<=ITP<=9 Planets (1=Mercury, 2=Venus, 3=Earth, etc.) C 10=Ulysses, 11=Helios 1; 12=Helios 2 C RR real height of deconvolution C dRR real distance between heights C XCbeg real start Carrington variable of arrays C XCend real end Carrington variable of arrays C nLng integer # longitudes C nLat integer # latitudes C nMap integer # heights C nT integer # times C FALLOFF real Density falloff with solar distance C VV(nLng,nLat,nT)real velocity map at height RR C DD(nLng,nLat,nT)real density map at height RR C XCshift(nLng,nLat,nMap,nT,3) real shift amount from array C Rconst real Tomography time constant ~25.38/27.275, but specific to the interval. C DVfact(nLng,nLat,nMap,nT) real velocity difference at height C DDfact(nLng,nLat,nMap,nT) real density difference at height C NCoff integer C XCstrt real Start rotation value C nCar integer # Carrington rotations C JDCar real*8 Julian Date at beginning of rotations C OUTPUTS: C FUNCTIONS/SUBROUTINES: C SunNewcomb, ECLIP_EQUATOR,WR2DARR C- subroutine Extractd(bExtract,RR,dRR,XCbe,XCtbeg,XCtend,nLng,nLat,nMap,nT,FALLOFF, & VV,DD,XCshift,Rconst,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar) parameter (NTP = 12) logical bExtract(NTP) character cPrefix(NTP)*2 /'ME','VE','EA','MA','JU','SA','UR','NE','PL','UL','H1','H2'/ real VV(nLng,nLat,nT), & DD(nLng,nLat,nT), & XCshift(nLng,nLat,nMap,nT,3), & DDfact (nLng,nLat,nMap,nT), & DVfact (nLng,nLat,nMap,nT) , & dXC(3) real*8 JDCar(nCar),JD,JEpoch,JDCntr real RVect(3),XCbe(2,nT) equivalence (ObsLng,RVect(1)),(ObsLat,RVect(2)),(ObsDis,RVect(3)) character cFile*14, & cStr*90 /';Year DOY Hour Height Lat Long Density Velocity'/, & cDate*23 logical bGetLun parameter (NTT = 60) ! 2*NT+1=# times extracted (covers 30 days) real*8 dT /0.25d0/! Time increment for extracted time series include 'MAPCOORDINATES.H' Bad = BadR4() dTinc = dT XCcntr = XCstrt+0.5 ! Carrington variable at center of rotation I = XCcntr ! Time when Earth was at center of rotation JDCntr = JDCar(I)+(JDCar(I+1)-JDCar(I))*(XCcntr-I) do ITP=1,NTP if (bExtract(ITP)) then ! Heliographic pos of extraction point at JDcntr call Julian(1,iYr,Doy,JDCntr,JEpoch) call ExtractPositionn(ITP,iYr,Doy,RVect) JD = JDcntr+(XCcntr-XMAP_OBS_POS(XCcntr,ObsLng))*SUN__SiderealP I = JD JD = I+nint((JD-I)/dT)*dT write (cFile,'(A,A,F8.3)') cPrefix(ITP),'_',NCoff+XCstrt if (.not. bGetLun(iU,cFile)) call Say('Extract','E','NOUNIT', & 'iGetLun unable to assign logical unit number') open (iU,file=cFile,status='NEW') call OSGetDateAndTime(cDate) write (iU,'(4A)') '; ',cFile(:itrim(cFile)),' created on ',cDate write (iU,'(A)' ) cStr(:itrim(cStr)) XCinc = 0.0 do I=-NTT,NTT call Julian(1,iYr,Doy,JD+I*dT,JEpoch) call ExtractPositionn(ITP,iYr,Doy,RVect) ObsD = Bad ObsV = Bad ObsXC = XMAP_OBS_POS(XCstrt,ObsLng) if(I.eq.-NTT) Xtbeg = ObsXC - XCstrt + 1.0 ! 1'st time if(I.gt.-NTT) then ! Other times dObs = ObsXC-ObsXCL else dObs = 0. end if ObsXCL = ObsXC if(abs(dObs).gt.0.5) XCinc = XCinc + 1.0 ObsXC = ObsXC + XCinc C print *, 'ObsXC', ObsXC ObsTim = Xtbeg + dTinc*(I+NTT)/(25.38/Rconst) C ObsTim = Xtbeg + dTinc*(I+NTT)/25.38 C print *, 'ObsTim',ObsTim C print *, 'XCstrt ObsLng ObsLat ObsXCL ObsXC', XCstrt, ObsLng, ObsLat, OBsXCL, ObsXC ObsXCN = ObsXC ObsXCN = min(XCbe(2,nT),max(XCbe(1,1),ObsXCN)) ObsTimN = ObsTim ObsTimN = min(XCtend,max(XCtbeg,ObstimN)) C print *, 'ObsTimN',ObsTimN C print *, 'XCstrt ObsLng ObsLat ObsXCL ObsXCN', XCstrt, ObsLng, ObsLat, OBsXCL, ObsXCN call Get4Dval(3,XCbe,XCtbeg,XCtend,RR,dRR,nLng,nLat, & nMap,nT,XCshift,1,ObsTimN,ObsXCN,ObsLat,ObsDis,dXC) if (dXC(1) .eq. Bad) stop 'Extract: oops 1' if (dXC(2) .eq. Bad) stop 'Extract: oops 2' if (dXC(3) .eq. Bad) stop 'Extract: oops 3' call Get4Dval(1,XCbe,XCtbeg,XCtend,RR,dRR,nLng,nLat, & nMap,nT,DDfact,1,ObsTimN,ObsXCN,ObsLat,ObsDis,DF) call Get4Dval(1,XCbe,XCtbeg,XCtend,RR,dRR,nLng,nLat, & nMap,nt,DVfact,1,ObsTimN,ObsXCN,ObsLat,ObsDis,VF) ObsXCN = ObsXC+dXC (1) ObsXCN = min(XCbe(2,nT),max(XCbe(1,1),ObsXCN)) ObsLa = ObsLat + dXC(2) ObsTimN = ObsTim+dXC (3) ObsTimN = min(XCtend,max(XCtbeg,ObstimN)) call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,DD,1, & ObsXCN,ObsLa,ObsTimN,ObsD) C print *, 'DF ObsD', DF, ObsD if (DF .ne. Bad .and. ObsD .ne. Bad) ObsD = ObsD*DF*((RR/ObsDis)**FALLOFF) call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,VV,1, & ObsXCN,ObsLa,ObsTimN,ObsV) C print *, 'VF ObsV', VF, ObsV if (VF .ne. Bad .and. ObsV .ne. Bad) ObsV = ObsV*VF ObsD = max(ObsD,-9999.999) ! Should fit output format ObsV = max(ObsV,-9999.999) write (iU,'(I5,I4,3F7.3,3F9.3)') iYr,int(Doy),(Doy-int(Doy))*24,ObsDis,ObsLat,ObsLng,ObsD,ObsV end do iU = iFreeLun(iU) end if end do return end subroutine ExtractPositionn(ITP,iYr,Doy,RVect) real RVect(3),VVect(5) real*8 dLngSun,dLatSun,dDisSun if (ITP .eq. 3) then ! Earth call SunNewcomb(0,iYr,Doy,dLngSun,dLatSun,dDisSun) RVect(1) = dLngSun+180 RVect(2) = -dLatSun/3600 ! arcsec -> degrees RVect(3) = dDisSun else if (ITP .eq. 10) then ! Ulysses call UlyssesOrbit(iYr,Doy,RVect,VVect) else if (ITP .eq. 11 .or. ! Helios 1 & ITP .eq. 12) then ! Helios 2 call HOSOrbit(ITP-10,iYr,Doy,RVect(3),RVect(1),VR,VT) RVect(2) = 0 else call PlanetOrbit(ITP,iYr,Doy,RVect,VVect) end if ! Convert to heliographic coordinates call ECLIPTIC_HELIOGRAPHIC(0,iYr,Doy,RVect(1),RVect(2)) return end