C C+ C NAME: C Write3D_bbtm_3n C PURPOSE: C Write magnetic field output to bb3d* file from input magnetic source files at the number and cadence of the tomography program C CATEGORY: C I/O C CALLING SEQUENCE: subroutine Write3D_bbtm_3n(bMagexx,XCbeGG,iT3D,nT3D,BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC,FALLOFFBN, & TT3D,XC3D,VV3D,BRss,BRCs,BTCs,BNCs,BR3D,BT3D,BN3D,iBBFlag) C INPUTS: C bMagexx logical If .true. write out magnetic 3D files. C XCbeGG(2,nTim) real Somehow related to the tomography times C Only XCbeGG(1,1) and XCbeGG(2,1) are used C iT3D integer 1st time index in TT3D processed C TT3D(iT3D) must be a valid time. C (the last time index processed is determined C internally as the last valid time in TT3D) C nT3D integer size of time dimension in arrays C TT3D,XC3D,VV3D and BRss. C (essential to control access to array C XC3D because of the extra trailing dimension) C BFACTORRR real A multiplicative factor used to provide a magnetic field enhancement over the CSSS C as suggested by T. Hoeksema (Fall AGU 2013) C BFACTORRC real A multiplicative factor used to provide a closed r field enhancement over the usual C C BFACTORTC real A multiplicative factor used to provide a closed t field enhancement over the usual C C BFACTORNC real A multiplicative factor used to provide a closed n field enhancement C C FALLOFFBN real Normal Field fall-off. This was found from the Helios spacecraft to be 1.34 C C The following three arrays are set up by Write3D_info. C C TT3D(nT3D) real Carrington times C These most likely will be the times spaced at C regular intervals in UT converted to Carrington times C XC3D(nLng,nLat,nRad,nT3D,3) real Shift needed to go back to source surface C (in Carrington variable). Will be negative. C Only XC3D(*,*,*,iT3D:jT3D,*) is used C VV3D(nLng,nLat,nRad,nT3D) real Velocity array C Only VV3D(*,*,*,iT3D:jT3D) is used C BRss(nLng,nLat,jT3D-kT3D+1) real scratch array; CSSS radial source surface magnetic field C C BRCs(nLng,nLat,jT3D-kT3D+1) real scratch array; Closed radial source surface magnetic field C C BTCs(nLng,nLat,jT3D-kT3D+1) real scratch array; Closed tangential source surface magnetic field C C BNCs(nLng,nLat,jT3D-kT3D+1) real scratch array; Closed normal source surface magnetic field C C BR3D(nLng,nLat,nRad) real scratch array; radial component written to file C C BT3D(nLng,nLat,nRad) real scratch array; tangential component written to file C C BN3D(nLng,nLat,nRad) real scratch array; normal component written to file C OUTPUTS: C (to file) C INCLUDE: include 'sun.h' include 't3d_array_3.h' include 'b3d_param_3.h' C CALLS: C FLINT, T3D_write_bb_3, T3D_get, T3D_set, T3D_iget, T3D_iset, T3D_get_grid_3 C ArrR4Step, ArrR4Constant, ArrR4Interpolate, BadR4, BField_Get_3, Say C ArrR4TimesConstant C PROCEDURE: C > TT3D may be only partially filled with valid times, e.g TT3D(kT3D:jT3D). C kT3D and jT3D are determined internally as the 1st and last valid time C (.ne. BadR4()) in TT3D. The input value iT3D must be greater/equal kT3D C and less/equal jT3D). Output files are produced only for times iT3D to jT3D. C kT3D usually is smaller than jT3D. Magnetic source surface data are collected C between time indices kT3D and jT3D. This ensures that heliospheric grid elements C at early times (just above iT3D) after traceback will still have a valid C time associated with them (above kT3D) at the source surface. C C > The radial dependence of Br follows from flux conservation: C Br = BB*(RR/r)^2 C The tangential component follows from geometrical considerations C by assuming that the magnetic field lies along the Parker spiral: C (Omega is angular velocity of Sun): C Bt/Br = r*cos(Lat)*Omega/V C Bt = BB*(Omega*RR*cos(Lat)/V)*(RR/r) C Output are 'normalized' magnetic field components (making the C values independent of radial distance): C r^2*Br = RR^2*BB C r *Bt = RR^2*BB *Omega*cos(Lat)/V C Expressing radial distances in units of r0 = 1 AU): C (r/r0)^2*Br = (RR/r0)^2*BB C (r/r0) *Bt = (RR/r0)^2*BB * r0*Omega*cos(Lat)/V C In terms of quantities defined in include file sun.h: C r0*Omega/V = Sun__AU*10^13*Sun__Omega*10^-6/(V[km/s]*10^5) = 100*Sun__AU*Sun__Omega/V[km/s] C C We are mostly interested in magnetic field at Earth, i.e. about r0=1 AU. C Insitu magnetic field are usually expressed in nT, so we use the same units here. C We write 'normalized magnetic fields' in units of nT to file: (RR/r0)^2*Br for the radial C and (RR/r0)*Bt for the tangential component. These quantities should be more or less C constant across the heliosphere, since the radial dropoff has been taken out. At 1 AU C (RR/r0 = 1) it should reproduce the insitu values. C C *TJD*'s unit troubles C X's output is in mT or 10^-2G. to bring to nT: C (15/215)^2 * (10^-6 T / mT) (10^9 nT/ T) C C MODIFICATION HISTORY: C SEP-1999, Paul Hick (UCSD/CASS) C ???-2002, Tamsen Dunn (UCSD/CASS; tdunn@ucsd.edu) Made time dep., fleshed out. C JUL-2002, Paul Hick (UCSD/CASS; pphick@ucsd.edu) C Valid magnetic fields are returned only in grid points where the C source surface Br and V are valid. Previously Br would be valid C if only the source surface Br was valid (even if V was bad). C AUG-2002, Tamsen Dunn (UCSD/CASS; tdunn@ucsd.edu) C Altered Bt component to depend entirely on V3D. If V3D is bad, use VSS. C This revealed bug below: C SEP-2002, Paul Hick (UCSD/CASS; pphick@ucsd.edu) C Fixed big error in use of XC3D (did not take into account that the time dimension C was different from the V3D array). Now three pieces of information are expected C that are calculated by Write3D_info: TT3D, times at which magnetic fields are C required; XC3D and V3D: shifts and velocity at times TT3D. XC3D and V3D now have C the same dimensions. C JUN-2003, Paul Hick (UCSD/CASS) C Added processing of multiple magnetic field types C OCT-2003, Paul Hick (UCSD/CASS; pphick@ucsd.edu) C Added extra arguments iT3D and jT3D C- real XCbeGG(*) ! (2,nTim) integer iT3D integer nT3D real TT3D(*) real XC3D(*) real VV3D(*) real BRss(*) real BRCs(*) real BTCs(*) real BNCs(*) real BR3D(*) real BT3D(*) real BN3D(*) character cSay*14 /'Write3D_bbtm_3'/ real eps /0.00002/ integer iBBFlag(4) integer Nss(3) ! Array dimensions at source surface real Pss(3) ! Location at source surface equivalence (Nss(1),nLng) equivalence (Nss(2),nLat) equivalence (Nss(3),nTim) real Scale(4) /0.0,1.0,0.0,1.0/ real PwrR(3) /2.0,1.0,1.34/ character cFile*256 logical BField_Get_3 logical bMagexx ! The normal component is given at RR with a -RR**-1.34 fall-off integer iB parameter (nMaxT = 600) real XCbegMAT(nMaxT) real XCendMAT(nMaxT) include 't3d_grid_fnc.h' !------- ! loc3D converts from 3D location (iLng,iLat,iRad,iTim,iPos) to linear ! offset relative to start of array. loc3D(iLng,iLat,iRad,iTim,iPos) = ((((iPos-1)*nT3D+iTim-1)*nRad+iRad-1)*nLat+iLat-1)*nLng+iLng PwrR(3) = FALLOFFBN Bad = BadR4() VT = 100.0*Sun__Omega*Sun__AU if (TT3D(iT3D) .eq. Bad) then ! Make sure TT3D(iT3D) is valid call Say(cSay,'W','time','at index iT3D is bad. No magnetic field calculated') return end if kT3D = iT3D-1 ! Search for first good time preceding iT3D do while (kT3D .ge. 1 .and. TT3D(kT3D) .ne. Bad) kT3D = kT3D-1 end do kT3D = kT3D+1 ! Index of earliest good time, or 1 jT3D = iT3D+1 ! Search for last good time after iT3D do while (jT3D .le. nT3D .and. TT3D(jT3D) .ne. Bad) jT3D = jT3D+1 end do jT3D = jT3D-1 ! Index of last good time, or nT3D print *, 'In ',cSay, kT3D,iT3D,jT3D,nT3D !------- ! Setting the iteration counter to one more than the total number of ! iterations will suppress adding the iteration number to the output file name. call T3D_iget(T3D__ITERATIONS, 0, L ) call T3D_iset(T3D__ITERATION , 0, L+1) Scale(2) = Scale(2)*BFACTORR Scale(4) = Scale(4)*BFACTORR call T3D_set(T3D__SCALE,4,Scale) ! No scale change, but set anyways call T3D_set(T3D__R_PWR,3,PwrR ) ! Radial normalization powers !------- ! Update T3D structure. Storing nT3D is essential because BField_Get_3 will use this ! to determine the dimension of the TT3D array. call T3D_iset(T3D__NTIM , 0, jT3D-kT3D+1) ! Must stay before T3D_get_grid_3 call !! !------- ! Extract grid information. Note that nTim = jT3D. TT, dTT and dTTi (not used) are zero. call T3D_get_grid_3(TT,dTT,RR,dRR,nLng,nLat,nRad,nTim, dTTi,nLng1,nLat1) if (nTim .gt. nMaxT) call Say(cSay,'E','nMaxT','parameter too small') call ArrR4Constant(nTim,XCbeGG(1),XCbegMAT) call ArrR4Constant(nTim,XCbeGG(2),XCendMAT) FIXBRCs = ((1.0/RR)**(2.0-FALLOFFBN)) ! Fixes the components to comply with Parker falloff (BVJ 2/2/2016) FIXBTCs = ((1.0/RR)**(1.0-FALLOFFBN)) FIXBNCs = -((1.0/RR)**(-FALLOFFBN)) ! Fixes the normal component to comply with a -RR**-1.34 fall-off (BVJ 5/5/2016) do iB=1,B3D__N ! there are four types of magnetic field data possible so far do iBB=1,4 iBBFlag(iBB) = 0 end do do iBB=1,4 ! there are four possible component maps for each magnetic field !------- ! The magnetic field is returned in nano Tesla (nT). if (iBB .eq. 1) then if (BField_Get_3(iB,iBB,TT3D(kT3D),XCbegMAT,XCendMAT,BRss)) then iBBFlag(IBB) = 1 C print *, 'Out of BField_Get_3 iBB = ', IBB go to 22 ! rr field end if else if (iBB .eq. 2) then if (BField_Get_3(iB,iBB,TT3D(kT3D),XCbegMAT,XCendMAT,BRCs)) then iBBFlag(IBB) = 1 go to 22 ! rc field end if else if (iBB .eq. 3) then if (BField_Get_3(iB,iBB,TT3D(kT3D),XCbegMAT,XCendMAT,BTCs)) then iBBFlag(IBB) = 1 C print *, 'Out of BField_Get_3 iBB = ', IBB go to 22 ! rt field end if else if (iBB .eq. 4) then if (BField_Get_3(iB,iBB,TT3D(kT3D),XCbegMAT,XCendMAT,BNCs)) then iBBFlag(IBB) = 1 C print *, 'Out of BField_Get_3 iBB = ', IBB go to 22 ! rn field end if end if if((iBBFlag(1)+iBBFlag(2)+iBBFlag(3)+iBBFlag(4)).eq.0) go to 23 ! no rr, rc, rt, or rn field data 22 continue if(iBB.eq.1) Write(*,'(A,3I3)') 'Write3d_bbtm_3n got CSSS regular fields' if(iBB.eq.2) Write(*,'(A,3I3)') 'Write3d_bbtm_3n got radial closed fields' if(iBB.eq.3) Write(*,'(A,3I3)') 'Write3d_bbtm_3n got tangential closed fields' if(iBB.eq.4) Write(*,'(A,3I3)') 'Write3d_bbtm_3n got normal closed fields' !------- ! Distance normalization: multiply radial component BR3D by RR*RR. ! The resulting r^2*Br should be in the range 1-10 nT. if(iBB.eq.1) call ArrR4TimesConstant(-nLng*nLat*nTim,BRss,RR*RR*BFACTORRR,BRss) if(iBB.eq.2) call ArrR4TimesConstant(-nLng*nLat*nTim,BRCs,RR*RR*FIXBRCs*BFACTORRC,BRCs) !((RR/1.3)**(2.0-1.34)) if(iBB.eq.3) call ArrR4TimesConstant(-nLng*nLat*nTim,BTCs,RR*FIXBTCs*BFACTORTC,BTCs) if(iBB.eq.4) call ArrR4TimesConstant(-nLng*nLat*nTim,BNCs,FIXBNCs*BFACTORNC,BNCs) end do C This begins the Output of the Magnetic field maps. if(bMagexx) then do iTim=iT3D,jT3D TTi = TT3D(iTim) XCbeg = XCbegMAT(iTim-kT3D+1) ! Needed for statement fnc XCvar XCend = XCendMAT(iTim-kT3D+1) XCrange = XCend-XCbeg XCb = TTi-0.5 ! Define region of interest (ROI) XCe = TTi+0.5 if (XCb .lt. XCbeg) then ! Keep ROI inside [XCbeg,XCend] XCb = XCbeg XCe = XCb+1.0 end if if (XCe .gt. XCend) then XCe = XCend XCb = XCe-1.0 end if call T3D_set (T3D__TIME , 0, TTi ) ! Time call T3D_iset(T3D__ITIME , 0, iTim-kT3D+1) ! Time index call T3D_set (T3D__XCMAT , 0, XCbeg) ! Begin of matrix call T3D_set (T3D__XCMAT+1, 0, XCend) ! End of matrix call T3D_set (T3D__XCROI , 0, XCb ) ! Begin ROI call T3D_set (T3D__XCROI+1, 0, XCe ) ! End ROI do iBB=1,4 if(iBB.eq.1) then ! First through at any time, make all bad to start do iRad=1,nRad do iLat=1,nLat do iLng=1,nLng iPos = loc3D(iLng,iLat,iRad,1,1) BR3D(iPos) = Bad BT3D(iPos) = Bad BN3D(iPos) = Bad end do end do end do end if do iRad=1,nRad Rhgt = RRhght(iRad) do iLat=1,nLat XL = cosd(XLdeg(iLat)) do iLng=1,nLng iPos = loc3D(iLng,iLat,iRad,1,1) !------- ! Trace heliospheric location back to the source surface. ! There are no bad values in XC3D. Statement functions can be used ! for longitude and latitude but NOT for the time. Pss(1) = XCindx( XCvar(iLng) + XC3D( loc3D(iLng,iLat,iRad,iTim,1) ) ) !------- ! We assume that the latitude shift is stored in units of grid spacings. ! This may be wrong. Currently the latitude shifts are zero so it doesn't ! matter. Pss(2) = float(iLat) + XC3D( loc3D(iLng,iLat,iRad,iTim,2) ) !------- ! To get the time index at the source surface we use the ArrR4Interpolate ! function. DO NOT USE the statement functions TTindx and TTtime here! ! Pss(3) will be bad if the traced-back time is outside the range of TT3D. ! Setting it to -1.0 forces FLINT to return a bad BR value. C Pss(3) = ArrR4Interpolate(nTim, TT3D(kT3D), TTi+XC3D( loc3D(iLng,iLat,iRad,iTim,3) )) Pss(3) = ArrR4Interpolate(nTim, TT3D(kT3D), TTi+(XC3D( loc3D(iLng,iLat,iRad,iTim,3) )/sngl(SUN__SYNODICP))) if (Pss(3) .eq. Bad) Pss(3) = -1.0 !------ ! BR is the normalized radial magnetic field at source surface location Pss ! (interpolate on source surface array BRss) C BR = FLINT(-3,Nss,BRss,Pss,eps) ! r^2*Br at source surface if(iBB.eq.1) then BR = FLINT(-3,Nss,BRss,Pss,eps) ! r^2*Br at source surface end if if(iBB.eq.2) then BRC = FLINT(-3,Nss,BRCs,Pss,eps) end if if(iBB.eq.3) BTC = FLINT(-3,Nss,BTCs,Pss,eps) ! r^1*Btc at source surface if(iBB.eq.4) BNC = FLINT(-3,Nss,BNCs,Pss,eps) ! r^1.34*Bnc at source surface if(iBB .eq. 1) then if(BR .ne. Bad) then BR3D(iPos) = BR else BR3D(iPos) = Bad end if end if if(iBB.eq.2) then if (BR3D(iPos).eq.Bad) BR3D(iPos) = 0.0 if (BRC .ne. Bad .and. BR3D(iPos) .ne. Bad) then BR3D(iPos) = BR3D(iPos) + BRC ! r^2*Br else BR3D(iPos) = Bad end if end if if(iBB.eq.3) then if (BTC .ne. Bad) then BT3D(iPos) = BTC ! r*Btc else BT3D(iPos) = Bad end if end if if(iBB .eq. 4) VV = VV3D(loc3D(iLng,iLat,iRad,iTim,1)) ! Velocity in heliosphere if (iBB .eq. 4) then if(BT3D(iPos) .eq. Bad .and. BR3D(iPos).ne. Bad) BT3D(iPos) = 0.0 if(BT3D(iPos) .ne. Bad .and. BR3D(iPos).eq. Bad) BR3D(iPos) = 0.0 if(BR3D(iPos) .ne. Bad .and. BT3D(iPos).ne. Bad .and. VV .ne. Bad) then BT3D(iPos) = -BR3D(iPos)*(VT*XL/VV) + BT3D(iPos) ! r*Bt + r*Btc else BT3D(iPos) = Bad end if if(BNC .ne. Bad) then BN3D(iPos) = BNC else BN3D(iPos) = Bad end if end if end do end do end do if(iBB.eq.4) call T3D_write_bb_3(B3D__PREFIX(iB),t3d, BR3D, BT3D, BN3D, 0, ' ', cFile) end do end do end if 23 continue end do return end