C+ C NAME: C Extractdvdm_3dmhd 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 Extractdvdm_3dmhd(bForecast,bExtract,cPrefixm,NTP,RR,dRR,XCbe,XCtbeg,XCtend,iYr,nLng,nLat, C & nMap,nT,NinterD,FALLOFFD,FALLOFFV,FALLOFFBR,FALLOFFBT,FALLOFFBN,VV,DD,XCshift3,XCshiftM, C & DVfact,DDfact,BRss,BRCs,BTCs,BNCs,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 cPrefixm character Prefix for the outparameter file C 1<=ITP<=9 Planets (1=Mercury, 2=Venus, 3=Earth, etc.) C 10=Ulysses, 11=Helios 1; 12=Helios 2 C NTP integer Total number of prefixes 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 iYr integer year at the beginning time of the tomography 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 FALLOFFD real Density falloff with solar distance C FALLOFFV real Velocity falloff with solar distance C FALLOFFBR real BR falloff with solar distance C FALLOFFBT real BT falloff with solar distance C FALLOFFBN real BN 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 XCshift3(nLng,nLat,nMap,nT,3) real shift amount from array C XCshiftM(nLng,nLat,nT,3) real shift amount to the magnetic field map from the base of the tomography 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 input array; CSSS radial source surface magnetic field C BRCs(nLng,nLat,jT3D-kT3D+1) real input array; Closed radial source surface magnetic field C BTCs(nLng,nLat,jT3D-kT3D+1) real input array; Closed tangential source surface magnetic field C BNCs(nLng,nLat,jT3D-kT3D+1) real input array; Closed normal source surface magnetic field C NCoff integer Carrington rotation offset 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 Extractdvdm_3dmhd(bForecast,bExtract,cPrefixm,NTP,RR,dRR,RRMS,XCbe,XCtbeg,XCtend,iYr,nLng,nLat, & nMap,nT,NinterD,FALLOFFD,FALLOFFV,FALLOFFBR,FALLOFFBT,FALLOFFBN,VV,DD,XCshift3,XCshiftM, & DVfact,DDfact,BRss,BRCs,BTCs,BNCs,NCoff,XCstrt,nCar,JDCar) logical bExtract(NTP) logical bForecast character cPrefixm(NTP)*3 real VV(nLng,nLat,nT), & DD(nLng,nLat,nT), & XCshift3(nLng,nLat,nMap,nT,3), ! This is the shift to get back to the 3DMHD tomographic origin base & XCshiftM(nLng,nLat,nT,3), ! This is the shift to get back to the magnetic field base from dRRM above it & DDfact (nLng,nLat,nMap,nT), & DVfact (nLng,nLat,nMap,nT), & dXC(3), & dXCM(3) real BRss(*) real BRCs(*) real BTCs(*) real BNCs(*) real*8 JDCar(nCar),JD,JEpoch,JDCntr real*8 XCtbeg,XCtend,XCtF,XCtE,XCtfull,X,Doy8,ObsTim,ObsTimN,ObsTimNM real*8 dT ! Time increment for extracted time series real ALng /3.0/ real RVect(3),XCbe(2,4*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 external EARTH8 include 'MAPCOORDINATES.H' write(*,'(A,I3,F9.6,2F8.3)') 'Into extractdvdm_3dmhd, nT, XCstrt, XCtbeg, XCtend', nT,XCstrt,XCtbeg,XCtend VT = 100.0*Sun__Omega*Sun__AU XCtF = 0.0d0 ! number of days past the tomography end time XCtE = 0.5d0*dble(Sun__SynodicP) ! Number of days past the end Carrington rotation center time if(bForecast) then XCtF = 1.0d0 ! number of days past the tomography end time time XCtE = 4.0d0 ! Number of days past the forecast time end if Bad = BadR4() dT = (XCtend - XCtbeg)/(dble((NinterD+1)*(nT-1))) ! extract interval length in days - typically 0.25 days or 6 hours NTTT = (XCtend+XCtF - XCtbeg)/dT + 1 ! total number of time steps possible from beginning of the tomography C print *, 'Extractdvdm 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) call Julian8(0,iYr,XCtbeg,JD,JEpoch) ! Julian date at XCtbeg I = int(JD) JDStrt = I + nint((JD-I)/dT)*dT ! Julian date at XCtbeg to nearest value of dT I = int(JDCntr) JDBeg = I + nint((JDCntr-I)/dT)*dT - 15.0d0 ! Beginning Julian date 15 days before center of rotation do ITP=1,NTP if (bExtract(ITP)) then print *, 'The location extracted is: ', cPrefixm(ITP) ! Heliographic pos of extraction point at JDcntr call Julian8(1,iYrc,Doy8,JDCntr,JEpoch) call ExtractPositionn8(ITP,iYrc,Doy8,RVect) JD = JDcntr+(XCcntr-XMAP_OBS_POS(XCcntr,ObsLng))*Sun__SynodicP I = int(JD) JD = I+nint((JD-I)/dT)*dT ! Julian date in days at observer location to nearest value of dT write (cFile,'(A,A,F8.3)') cPrefixm(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 XCinc = 0.0 ObsXCL = -100. xxinc = 0.0 do I=1,NTTT dTT = (dT * dble(I-1)) JD = JDStrt + dTT Doy8 = XCtbeg + dTT C print *, 'to here 1', JD, JDBeg, Doy8, XCtbeg, XCtend+XCtF, JDCntr+XCtE if(JD.ge.JDBeg.and.Doy8.gt.XCtbeg.and.Doy8.le.XCtend+XCtF.and.JD.le.JDCntr+XCtE) then ! Limits on extracted data C print *, 'to here 2', JD, JDBeg, Doy8, XCtbeg, XCtend+XCtF, JDCntr+XCtE call ExtractPositionn8(ITP,iYr,Doy8,RVect) ObsD = Bad ObsV = Bad ObsXC = XMAP_OBS_POS(XCstrt,ObsLng) 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 ObsTimN = min(XCtend,max(XCtbeg,ObstimN)) C print *, 'ObsTimN',ObsTimN call Get4Dval(3,XCbe,XCtbeg,XCtend,RR,dRR,nLng,nLat, & nMap,nT,XCshift3,1,ObsTimN,ObsXCN,ObsLat,ObsDis,dXC) C C print *, 'XCstrt ObsLng ObsLat ObsXCL ObsXCN dXC ', XCstrt, ObsLng, ObsLat, OBsXCL, ObsXCN, dXC(1), dXC(2), dXC(3) 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,4F12.4)') 'Doy8, ObsTimN, dXC(1), dXC(3) ', Doy8, ObsTimN, dXC(1), dXC(3) 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) ! Usual ObsXCN = min(XCbe(2,nT),max(XCbe(1,1),ObsXCN)) ObsLa = ObsLat + dXC(2) ! Usual ObsTimN = ObsTim + dXC(3) ! Usual ObsTimN = min(XCtend,max(XCtbeg,ObstimN)) C Plasma density and velocity 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)**FALLOFFD) call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,VV,1, & ObsXCN,ObsLa,ObsTimN,ObsV) if (VF .ne. Bad .and. ObsV .ne. Bad) ObsV = ObsV*VF*((RR/ObsDis)**FALLOFFV) C C Here put in the main difference between the 3dmhd extract and the kinematic model extract. C if(RR.gt.RRMS) then C call Get3TDval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,F3D,N,XC,XL,TI,FF) call Get3Dtval(XCbe,XCtbeg,XCtend,nLng,nLat,nT, ! Get from xcshiftM, dRRM above & XCshiftM(1,1,1,1),1,ObsXCN,ObsLa,ObsTimN,dXCM(1)) ! from lowest level to add to dXC call Get3Dtval(XCbe,XCtbeg,XCtend,nLng,nLat,nT, ! Get from xcshiftM, dRRM above & XCshiftM(1,1,1,2),1,ObsXCN,ObsLa,ObsTimN,dXCM(2)) ! from lowest level to add to dXC call Get3Dtval(XCbe,XCtbeg,XCtend,nLng,nLat,nT, ! Get from xcshiftM, dRRM above & XCshiftM(1,1,1,3),1,ObsXCN,ObsLa,ObsTimN,dXCM(3)) ! from lowest level to add to dXC C write (*,'(A,6F10.6)') 'DXCMs =', dXC(1), dXCM(1), dXC(2), dXCM(2), dXC(3), dXCM(3) dXCM(1) = dXCM(1) + dXC(1) ! Used to shift to magnetic field location dXCM(2) = dXCM(2) + dXC(2) dXCM(3) = dXCM(3) + dXC(3) else if(RR.lt.RRMS) then print *, 'RR is less than RRMS in extractdvdm_3DMHD.f not supported - stop' stop end if dXCM(1) = dXC(1) ! Used to shift to magnetic field location dXCM(2) = dXC(2) ! Zero change from dXC if RR = RRMS dXCM(3) = dXC(3) end if ObsXCNM = ObsXC + dXCM(1) ! Magnetic ObsXCNM = min(XCbe(2,nT),max(XCbe(1,1),ObsXCNM)) ObsLaM = ObsLat + dXCM(2) ! Magnetic ObsTimNM = ObsTim + dXCM(3) ! Magnetic ObsTimNM = min(XCtend,max(XCtbeg,ObstimNM)) C Magnetic Fields call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,BRss,1, ! Br*R^2 (True CSSS Br at this height) & ObsXCNM,ObsLaM,ObsTimNM,ObsBrr) call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,BRCs,1, ! Br*R^2 (True closed Br at this height) & ObsXCNM,ObsLaM,ObsTimNM,ObsBrc) call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,BTCs,1, ! Br*R^1 (Closed Bt at this height) & ObsXCNM,ObsLaM,ObsTimNM,ObsBtc) call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,BNCs,1, ! Br*R^1 (Closed Bn at this height) & ObsXCNM,ObsLaM,ObsTimNM,ObsBnc) C print *, 'After = ', JD+I*DT, iYr, DF, VF, Doy8, ObsXC, ObsXCN, ObsLa, ObsTimNM, ObsDis, ObsD, ObsV if (ObsBrr .ne. Bad .and. ObsBrc .ne. Bad) then ! Bt*R C OBsBr = (ObsBrr + ObsBrc)*((RR/ObsDis)**(FALLOFFBR-2.0)) OBsBr = (ObsBrr + ObsBrc)*((1.0/ObsDis)**FALLOFFBR) end if if (ObsBrr .ne. Bad .and. ObsBrc .eq. Bad) then C OBsBr = (ObsBrr)*((RR/ObsDis)**(FALLOFFBR-2.0)) OBsBr = (ObsBrr)*((1.0/ObsDis)**FALLOFFBR) end if if (ObsBrr .eq. Bad .and. ObsBrc .ne. Bad) then C OBsBr = (ObsBrc)*((RR/ObsDis)**(FALLOFFBR-2.0)) OBsBr = (ObsBrc)*((1.0/ObsDis)**FALLOFFBR) end if if (ObsBrr .eq. Bad .and. ObsBrc .eq. Bad) then ObsBr = Bad end if ! (This is the true Br at this height) if (ObsBr .ne. Bad .and. ObsV .ne. Bad) then ! Bt*R XL = cosd(ObsLa) C OBsBt1 = -ObsBr*(VT*XL/ObsV)*((RR/ObsDis)**(FALLOFFBT-1.0)) OBsBt1 = -ObsBr*(VT*XL/ObsV)*((1.0/ObsDis)**FALLOFFBT) else OBsBT1 = Bad end if if(ObsBtc.ne.Bad.and.OBsBT1.ne.Bad) then OBsBt = ObsBtc + OBsBT1 end if if(ObsBtc.ne.Bad.and.OBsBT1.eq.Bad) then OBsBt = ObsBtc end if if(ObsBtc.eq.Bad.and.OBsBT1.ne.Bad) then OBsBt = OBsBT1 end if if(ObsBtc.eq.Bad.and.OBsBT1.eq.Bad) then OBsBt = Bad end if if (ObsBnc .ne. Bad) then OBsBn = ObsBnc*((1.0/ObsDis)**FALLOFFBN) ! Bn falls-off at 1.34, and is not corrected to R^1.0 else ObsBn = Bad end if C Write(*,'(A,2F6.3,F6.3,F8.3,3F6.2)') 'DF, VF, ObsD, ObsV, ObsBr, ObsBt, ObsBn', C & DF, VF, ObsD, ObsV, ObsBr, ObsBt, ObsBn ObsD = max(ObsD, -9999.999) ! Should fit output format ObsV = max(ObsV, -9999.999) ObsBr = max(ObsBr,-9999.999) ObsBt = max(ObsBt,-9999.999) ObsBn = max(ObsBn,-9999.999) else print *, cPrefixm(ITP),' is out of range of the 3D matrix' ObsD = -9999.999 ObsV = -9999.999 ObsBr = -9999.999 ObsBt = -9999.999 ObsBn = -9999.999 end if C write (*,'(I5,I4,2F7.3,F8.3,7F7.3)') iYr,int(Doy8),(Doy8-int(Doy8))*24,ObsD,ObsV,ObsBrr,ObsBrc,ObsBr, C & ObsBtc,-ObsBr*(VT*XL/ObsV),ObsBt,ObsBn write (iU,'(I5,I4,3F7.3,6F9.3)') iYr,int(Doy8),(Doy8-int(Doy8))*24,ObsDis,ObsLat,ObsLng,ObsD,ObsV,ObsBr,ObsBt,ObsBn end if end do iU = iFreeLun(iU) end if end do return end