C+ C NAME: C Extractdnnm C PURPOSE: C Extracts interpolated velocity, density, and radial, tangential, and latitude magnetic field IPS deconvolved measurements at C a given time and position. C CALLING SEQUENCE: C call Extractdn(bForecast,bTrace,bExtract,RR,dRR,XCbe,XCtbeg,XCtend,nLng,nLat,nMap,nT,NinterD,FALLOFF, C & nMap,nT,NinterD,FALLOFF,VV,DD,XCshift,DVfact,DDfact,TT3D,BRss,NCoff,XCstrt,nCar,JDCar) C INPUTS: C bForecast logical If .TRUE. this is forecast mode and then add a few days to the output 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 XCbe(2,nT) real start and ending Carrington variable of arrays C XCtbeg real*8 beginning time of observations C XCtend real*8 ending time of observations C nLng integer # longitudes C nLat integer # latitudes C nMap integer # heights C nT integer # times C NinterD integer # of intermediate data times interpolated 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 DVfact(nLng,nLat,nMap,nT) real velocity difference at height C DDfact(nLng,nLat,nMap,nT) real density difference at height C BRss(nLng,nLat,jT3D-kT3D+1) real scratch array; source surface magnetic field C C BR3D(nLng,nLat,nRad) real scratch array; radial component written to file C BT3D(nLng,nLat,nRad) real scratch array; tangential component written to file 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 Extractdnnm(bForecast,bTrace,distTr,bExtract,cPrefix,NTP,RR,dRR,XCbe,XCtbeg,XCtend,nLng,nLat, & nMap,nT,NinterD,FALLOFF,VV,DD,XCshift,DVfact,DDfact,TT3D,BRss,NCoff,XCstrt,nCar,JDCar) logical bExtract(NTP) logical bForecast,bTrace character cPrefix(NTP)*2 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), & dXCC(3), & dXCCC(3,100), & dXCCCC(3), & ObsVV(100), & DistptC(100) real TT3D(*) real BRss(*) real*8 JDCar(nCar),JD,JEpoch,JDCntr real*8 XCtbeg,XCtend,XCtF,XCtfull,X,Doy8,ObsTim,ObsTimN,Doy8C real*8 dT ! Time increment for extracted time series real*8 DoypY /0.0d0/ real ALng /3.0/ 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 character cFile1*20 logical bGetLun external EARTH8 include 'MAPCOORDINATES.H' VT = 100.0*Sun__Omega*Sun__AU print *, ' Into extractdn.f nT, XCstrt, XCtbeg, XCtend, VT =',nT,XCstrt, XCtbeg, XCtend, VT XCtF = 0.0d0 if(bForecast) XCtF = 1.0d0 Bad = BadR4() dT = (XCtend - XCtbeg)/(dble((NinterD+1)*(nT-1))) ! extract interval length in days NTT = nint(20.0/dT) ! 2*20 =# times extracteint(Doy8),(Doy8-int(Doy8))*24d (potentially covers +-20 days from the center of the rotation) C print *, 'Extractdn Interval length in days', dT XCcntr = XCstrt + 0.5 ! Carrington variable at center of rotation I = XCcntr ! Carrington interval before Earth was at center of rotation JDCntr = JDCar(I)+(JDCar(I+1)-JDCar(I))*dble(XCcntr-I) do ITP=1,NTP if (bExtract(ITP)) then ! Heliographic pos of extraction point at JDcntr call Julian8(1,iYr,Doy8,JDCntr,JEpoch) call ExtractPositionn8(ITP,iYr,Doy8,RVect) JD = JDcntr+(XCcntr-XMAP_OBS_POS(XCcntr,ObsLng))*Sun__SynodicP 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') write (cFile1,'(A,A,F8.3,A)') cPrefix(ITP),'_',NCoff+XCstrt,'_trace' if (.not.bGetLun(iU1,cFile1)) call Say('Extract','E','NOUNIT', & 'iGetLun unable to assign logical unit number for Trace') open (iU1,file=cFile1,status='NEW') call OSGetDateAndTime(cDate) write (iU,'(4A)') '; ',cFile(:itrim(cFile)),' created on ',cDate write (iU1,'(4A)') '; ',cFile1(:itrim(cFile1)),' created on ',cDate XCinc = 0.0 ObsXCL = -100. iYrL = 0 xxinc = 0.0 DoypY = 0.0d0 do I=-NTT,NTT call Julian8(1,iYr,Doy8,dble(JD+I*dT),JEpoch) iYrr = iYr C if(Doy8 .lt. 180.0d0 .and. XCtend+XCtF .gt. 360.0d0) then ! Changed from previous on 24 Jan 201 BVJ C iYrr = iYrr + 1 C xxcinc = -1.0 C end if if(iYrr.gt.iYrL.and.I.gt.-NTT) then ! Year change. Changed to .ge. on 24 Jan 2012 BVJ Leap = 0 ! Suppose it's not a leap year if ((iYrL)/4*4 .eq. (iYrL) .and. ((iYrL) .lt. 1582 .or. ! It's a leap year & iYrL/100*100 .ne. iYrL .or. iYrL/400*400 .eq. iYrL)) Leap = 1 DoypY = dble(365+Leap) iYrL = iYrr end if iYrL = iYrr C print *, 'to here 1', iyr, iyrr, iYrL, Doy8, DoypV, XCtbeg, Doy8+DoypY, XCtend+XCtF if((Doy8+DoypY).ge.XCtbeg.and.(Doy8+DoypY).le.XCtend+XCtF) then C print *, 'to here 2', iyr, iyrr, iYrL, Doy8, DoypV, XCtbeg, Doy8+DoypY, XCtend+XCtF call ExtractPositionn8(ITP,iYr,Doy8,RVect) ObsD = Bad ObsV = Bad ObsXC = XMAP_OBS_POS(XCstrt,ObsLng) C if(I.gt.-NTT) then ! Other times if(ObsXCL.gt.-99.) then ! Changed from previous on 20 April 2011 BVJ dObs = ObsXC-ObsXCL else dObs = 0. end if ObsXCL = ObsXC if(abs(dObs).gt.0.5) XCinc = XCinc + 1.0 + xxcinc ObsXC = ObsXC + XCinc C print *, 'ObsXC', ObsXC ObsTim = Doy8 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 + DoypY ObsTimN = min(XCtend,max(XCtbeg,ObstimN)) C print *, 'ObsTimN',ObsTimN call Get4Dval(3,XCbe,XCtbeg,XCtend,RR,dRR,nLng,nLat, & nMap,nT,XCshift,1,ObsTimN,ObsXCN,ObsLat,ObsDis,dXC) C ******** C if (dXC(1) .eq. Bad) stop 'Extract: oops 1' C if (dXC(2) .eq. Bad) stop 'Extract: oops 2' C if (dXC(3) .eq. Bad) stop 'Extract: oops 3' C write(*,'(A,3F12.4)') 'XCstrt ObsXCN ObsLng', float(NCOFF)+XCstrt, float(NCOFF)+ObsXCN, ObsLng C write(*,'(A,3F12.4)') 'Doy8+DoypY, ObsTimN, dXC(1) ', Doy8+DoypY, ObsTimN, dXC(1) if (dXC(1).ne.Bad.and.dXC(2).ne.Bad.and.dXC(3).ne.Bad) then 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) call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,BRss,1, ! Br*R^2 & ObsXCN,ObsLa,ObsTimN,ObsBr) C print *, 'After = ', JD+I*DT, iYr, DF, VF, Doy8, ObsXC, ObsXCN, ObsLa, ObsTimN, ObsDis, ObsD, ObsV if (VF .ne. Bad .and. ObsV .ne. Bad) ObsV = ObsV*VF if (ObsBr .ne. Bad .and. ObsV .ne. Bad) then XL = cosd(ObsLa) OBsBt = -ObsBr*(VT*XL/ObsV) ! Bt*R else ObsBr = Bad ! Erase good ObsBr if Obsv is bad ObsBt = Bad end if ObsBL = 0.0 C Write(*,'(A,2F7.4,F8.3,F8.4,3F8.3)') 'DF VF, ObsD, ObsV, ObsBr, ObsBt, ObsBL', DF, VF, ObsD, ObsV, ObsBr, ObsBt, ObsBL ObsD = max(ObsD, -9999.999) ! Should fit output format ObsV = max(ObsV, -9999.999) ObsBr = max(ObsBr,-9999.999) ObsBt = max(ObsBt,-9999.999) ObsBL = max(ObsBL,-9999.999) else print *, cPrefix(ITP),' is out of range of the 3D matrix' ObsD = -9999.999 ObsV = -9999.999 ObsBr = -9999.999 ObsBt = -9999.999 ObsBL = -9999.999 end if C ********** Distpt = distTr ObsVV(1) = ObsV if(ObsVV(1).eq.-9999.999) ObsVV(1) = 400.0 iTracen = 0 iTrace = 1 errorf = 0.00001 CONST3 = (sngl(Sun__AU)*100000000.0)/86400.0 CONST1 = CONST3/sngl(Sun__SIDEREALP) if(cPrefix(ITP).eq.'e3') print *, Sun__AU, Sun__SIDEREALP, CONST1, CONST3 101 continue Distp = ObsDis htfrac1 = CONST1*(Distpt-RR)/ObsVV(iTrace) htfrac2 = (Distpt-RR)/(ObsDis-RR) htfrac3 = CONST3*(Distpt-RR)/ObsVV(iTrace) dXC1 = dXC(1) - htfrac1 dXC2 = dXC(2) - htfrac2 dXC3 = dXC(3) - htfrac3 if(cPrefix(ITP).eq.'e3') print *, 'First', iTrace, Distp, Distpt, htfrac1, htfrac2, htfrac3, & dXC(1), dXC(2), dXC(3), dXC1, dXC2, dXC3 call Get4Dval(3,XCbe,XCtbeg,XCtend,RR,dRR,nLng,nLat, & nMap,nT,XCshift,1,ObsTimN+dXC3,ObsXCN+dXC1,ObsLat+dXC2,Distpt,dXCCCC) if(dXCCCC(1).ne.Bad.and.dXCCCC(2).ne.Bad.and.dXCCCC(3).ne.Bad) then call Get4Dval(1,XCbe,XCtbeg,XCtend,RR,dRR,nLng,nLat, & nMap,nT,DVfact,1,ObsTimN+dXC3,ObsXCN+dXC1,ObsLat+dXC2-dXCCCC(2),Distpt,VFF) call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,VV,1, & ObsXCN,ObsLa+dXC2-dXCCCC(2),ObsTimN,ObsVV(iTrace+1)) if(ObsVV(iTrace+1).ne.Bad.and.VFF.ne.Bad.and.VFF.ne.0.0) then ObsVV(iTrace+1) = ObsVV(iTrace+1)*VFF htfract1 = CONST1*(Distpt-RR)/ObsVV(iTrace+1) htfract3 = CONST3*(Distpt-RR)/ObsVV(iTrace+1) else dXCC(1) = Bad dXCC(3) = Bad end if else dXCC(1) = Bad dXCC(3) = Bad end if dXC1f = dXC(1)*errorf dXC2f = dXC(2)*errorf dXC3f = dXC(3)*errorf do iii=1,10 call Get4Dval(3,XCbe,XCtbeg,XCtend,RR,dRR,nLng,nLat, & nMap,nT,XCshift,1,ObsTimN+dXC3,ObsXCN+dXC1,ObsLat+dXC2,Distp,dXCC) C if(cPrefix(ITP).eq.'e3') print *, iii, Distp, Distpt, dXC(1), dXC(2), dXC(3), dXCC(1), dXCC(2), dXCC(3) Distp = Distpt dxC1L = dxC1 dxC2L = dxC2 dxC3L = dxC3 if (dXC(1).ne.Bad.and.dXC(2).ne.Bad.and.dXC(3).ne.Bad.and.dXCC(1).ne.Bad.and.dXCC(2).ne.Bad.and.dXCC(3).ne.Bad) then dxC1 = -dXC(1) + htfrac1 dxC2 = dXCC(2) - dXC(2) + htfrac2 ! After iteration dXCC(2) contains the increment of the pt to the base dxC3 = -dXC(3) + htfrac3 else go to 100 end if if(abs(dxC2L-dxC2).le.abs(dXC2f)) go to 100 end do 100 continue if (dXC(1).ne.Bad.and.dXC(2).ne.Bad.and.dXC(3).ne.Bad.and.dXCC(1).ne.Bad.and.dXCC(2).ne.Bad.and.dXCC(3).ne.Bad) then dXCC(1) = dXC(1) + htfrac1 ! Now change dXCC to the location of the object to the point. dXCC(2) = dXC(2) - dXCC(2) dXCC(3) = dXC(3) + htfrac3 else dXCC(1) = Bad dXCC(2) = Bad dXCC(3) = Bad end if if(iTrace.eq.1) then dXCC1 = dXCC(1) dXCC2 = dXCC(2) dXCC3 = dXCC(3) Distpt1 = Distpt end if if(bTrace) then dXCCC(1,iTrace) = dXCC(1) dXCCC(2,iTrace) = dXCC(2) dXCCC(3,iTrace) = dXCC(3) DistptC(iTrace) = Distpt iTrace = iTrace + 1 if(iTrace.gt.100) go to 102 Distpt = Distpt + 0.05 iDistpt = int(Distpt*20.0 + errorf) Distpt = float(iDistpt)*0.05 if(Distpt.ge.ObsDis) then if(iTracen.eq.0) then Distpt = ObsDis iTracen = iTrace go to 101 else go to 102 end if end if go to 101 end if 102 continue C print *, 'XCstrt ObsLng ObsLat ObsXCL ObsXCN dXC ', XCstrt, ObsLng, ObsLat, OBsXCL, ObsXCN, dXC(1), dXC(2), dXC(3) C ********** intDoy8C = -999 hrC = -9999.999 if(dXCC3.ne.Bad) then intDoy8C = int(Doy8 + dble(dXCC3)) hrC = (Doy8 + dble(dXCC3) - intDoy8C)*24 end if ObsLatC = -9999.999 if(dXCC2.ne.Bad) ObsLatC = ObsLat+dXCC2 ObslngdX = -9999.999 if(dXCC1.ne.Bad) then ObslngdX = ObsLng-dXCC1*360.0 if(ObslngdX.gt.360.0) ObslngdX = ObslngdX - 360.0 end if write (iU,'(I5,I4,3F7.3,6F9.3,I4,4F9.3)') iYr,int(Doy8),(Doy8-int(Doy8))*24,ObsDis,ObsLat,ObsLng, & ObsD,ObsV,ObsBr,ObsBt,ObsBL,intDoy8C,hrC,Distpt1,ObsLatC,ObsLngdX if(bTrace) then do nnn = 1, iTracen intDoy8C = -999 hrC = -9999.999 if(dXCCC(3,nnn).ne.Bad) then intDoy8C = int(Doy8 + dble(dXCCC(3,nnn))) hrC = (Doy8 + dble(dXCCC(3,nnn)) - intDoy8C)*24 end if ObsLatC = -9999.999 if(dXCCC(2,nnn).ne.Bad) ObsLatC = ObsLat+dXCCC(2,nnn) ObslngdX = -9999.999 if(dXCCC(1,nnn).ne.Bad) then ObslngdX = ObsLng-dXCCC(1,nnn)*360.0 if(ObslngdX.gt.360.0) ObslngdX = ObslngdX - 360.0 end if write (iU1,'(I5,I4,3F7.3,F9.3,2I4,4F9.3)') iYr,int(Doy8),(Doy8-int(Doy8))*24,ObsDis,ObsLat,ObsLng, & nnn,intDoy8C,hrC,DistptC(nnn),ObsLatC,ObsLngdX end do end if end if end do iU = iFreeLun(iU) iU1 = iFreeLun(iU1) end if end do return end C call ArrR4Constant(nTim,XCbe(1),XCbegMAT) C call ArrR4Constant(nTim,XCbe(2),XCendMAT) C do iB=1,B3D__N ! Magnetic field from different solar sites - use only nso data C do iTim=iT3D,jT3D C TTi = TT3D(iTim) C XCbeg = XCbegMAT(iTim-kT3D+1) ! Needed for statement fnc XCvar C XCend = XCendMAT(iTim-kT3D+1)