C+ C NAME: C Write3D_bb_UT C PURPOSE: C Write magnetic field output to b3d* file C CATEGORY: C I/O C CALLING SEQUENCE: subroutine Write3D_bb_UT(cWild,NT,nCar,JDCar,XCbegMAT,XCendMAT,XCbegROI,XCendROI,XC3D,VV3D,BR3D,BT3D,BRSS,BTSS) C INPUTS: C cWild*(*) character wild card for locating magnetic field files C in the format processed by href=ForeignFile=. C NT integer # times written to separate files C C nCar integer C JDCar(nCar) real*8 C C XCbegMAT(nTim) real start Carrington variable of matrix (- NCoff) C XCendMAT(nTim) real end Carrington variable of matrix (- NCoff) C (range of matrix typically is three Carrington rotations) C XCbegROI(nTim) real start Carrington variable of region of interest C XCendROI(nTim) real end Carrington variable of region of interest C (the region of interest typically covers one whole C Carrington rotation C XC3D(3,nLng,nLat,nRad,nTim) C real shift from source surface (in Carrington variable) C Array should not have any bad values. C VV3D(nLng,nLat,nRad,nTim) C real 3D velocity array (may contain bad values, filled in C by e.g. Write3D_nv. C BR3D(nLng,nLat,nRad) real scratch array C BT3D(nLng,nLat,nRad) real scratch array C BRSS(nLng,nLat,nTim) real scratch array to read Br C BTSS(nLng,nLat,nTim) real scratch array to read Bt (not used) C OUTPUTS: C (to file) C RESTRICTIONS: C > The VV3D array can have bad values in it. Usually this routine is called after C href=Write3D_nv= has added bad values to the velocity array. C > nRad must be larger/equal 3. C CALLS: C T3D_set, T3D_iget, T3D_get_grid, BadR4, ArrR4TimesConstant, Str2Str, Int2Str C Flt2Str, WR2DARR, FLINT C INCLUDE: include 'sun.h' include 't3d_param.h' include 't3d_array.h' include 't3d_index.h' C include 't3d_grid_fnc.h' C include 't3d_loc_fnc.h' C EXTERNAL: external EARTH C PROCEDURE: 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 The Stanford source surface maps at RR = 2.5 Rsun give the radial magnetic field BB in muT. C We are mostly interested in magnetic field at Earth, i.e. about r0=1 AU. We expect magnetic C fields of order (RR/r0)^2*BB ~ 10^-4*BB ~ 0.1 nT. This extra factor of 10^4 is included in the C output i.e. output is in units of 0.1 nT. C (In situ magnetic fields at Earth are typically in the nT range, at least one order of C magnitude larger than predicted from the Stanford maps, probably indicating that the Stanford C potential field maps underestimate the magnetic field near the current sheet by that much.) C MODIFICATION HISTORY: C SEP-1999, Paul Hick (UCSD/CASS; pphick@ucsd.edu) C- integer NT integer nCar real*8 JDCar(nCar) real XCbegMAT(*) real XCendMAT(*) real XCbegROI(*) real XCendROI(*) real VV3D (*) real XC3D (*) real BR3D (*) real BT3D (*) real BRSS (*) real BTSS (*) real*8 JDb real*8 JDe real*8 dJD real*8 FLINT8 real Tiny /0.0/ real*8 Tiny8 /0d0/ integer NDSS(3) integer ND3D(4) integer NDCC(5) real PPSS(3) real PP3D(4) real PPCC(5) equivalence (ND3D, NDCC(2)) real Scale(4) /0.0,1E-4, 0.0,1E-4/ real R_pwr(2) /2.0,1.0/ character cFile*120 character cSay*13 /'Write3D_bb_UT'/ logical bMod logical BField_Get logical bOK include 't3d_grid_fnc.h' include 't3d_loc_fnc.h' if (NT .le. 0) return call T3D_set(T3D__SCALE,4,Scale) ! Conversion to units of 0.1 nT call T3D_set(T3D__R_PWR,2,R_pwr) !------- ! Get source surface map. The source surface must be at 2.5 solar radii. If this is ! not the case T3D_ReadB returns .FALSE. ! Source surface arrays have dimensions [nLng,nLat,nTim] ! BTSS is not used (yet). call T3D_iget(T3D__NRAD,0,nRad) call T3D_iset(T3D__NRAD,0, 1) !------- ! The magnetic field is returned in nano Tesla (nT). ! If source surface magnetic field found, normalize by multiplying ! the radial component BR3D by RR*RR. bOK = BField_Get(cWild,TT,XCbegMAT,XCendMAT,BRSS) call T3D_iset(T3D__NRAD,0,nRad) if (.not. bOK) return call T3D_iget(T3D__MODE,0,MODE) bMod = iand(MODE,TOM__MOD) .ne. 0 call T3D_get_grid(TT,dTT,RR,dRR, nLng,nLat,nRad,nTim,dTTi,nLng1,nLat1) NDSS(1) = nLng ! Dimensions of source surface arrays NDSS(2) = nLat NDSS(3) = nTim !------- ! This also fills ND3D (equivalenced to NDCC) NDCC(1) = PRJ__N ! Dimensions of XC3D array NDCC(2) = nLng NDCC(3) = nLat NDCC(4) = nRad NDCC(5) = nTim ! May be degenerate (=1) Bad = BadR4() JDb = TTtime( 1) ! First and last time in Carrington variable JDe = TTtime(nTim) JDb = FLINT8(1,nCar,JDCar,JDb,Tiny8) ! First and last time in Julian dates JDe = FLINT8(1,nCar,JDCar,JDe,Tiny8) dJD = 1d0/NT ! Time increment in days JDb = int(JDb)+( 1+int((JDb-int(JDb))/dJD) )*dJD! First time rounded to dJD JDe = int(JDe)+( int((JDe-int(JDe))/dJD) )*dJD! Last time rounded to dJD nTim0 = 1+nint((JDe-JDb)/dJD) ! # times in equal steps in UT !------- ! Convert radial component to r^2*Br. ! The division by Scale(2) is used to change magnetic units. call ArrR4TimesConstant(-nLng*nLat*nTim,BRSS,RR*RR/Scale(2),BRSS) !------- ! For the magnetic field extraction we assume that the tomography source surface is the same ! as the magnetic field source surface (i.e. 2.5 solar radii for the Stanford maps). VT = 100*SUN__OMEGA*SUN__AU do L=0,nTim0-1 call Julian(1,iYr,Doy,JDb+L*dJD,JDe) ! Time: Julian day --> iYr,Doy TL = XMAP_SC_POS(EARTH,iYr,Doy,nCar,JDCar) ! Time (Carrington time) call T3D_iset(T3D__ITIME ,0,0 ) call T3D_set (T3D__TIME ,0,TL) TL = TTindx(TL) ! Time index XCbegL = FLINT(1,nTim,XCbegMAT,TL,Tiny) ! Matrix limits at time TL XCendL = FLINT(1,nTim,XCendMAT,TL,Tiny) call T3D_set (T3D__XCMAT ,0,XCbegL) call T3D_set (T3D__XCMAT+1,0,XCendL) call T3D_set (T3D__XCROI ,0,FLINT(1,nTim,XCbegROI,TL,Tiny)) call T3D_set (T3D__XCROI+1,0,FLINT(1,nTim,XCendROI,TL,Tiny)) PPCC(5) = TL ! Time index do K=1,nRad PPCC(4) = float(K) do J=1,nLat PPCC(3) = float(J) ! Latitude index XL = cosd(XLdeg(J)) do I=1,nLng !------ ! PPCC contains the position for interpolation on the 5-dim XC3D array. ! The first four are set to integer array indices (converted to float). ! Only the last (time) dimension will be interpolated. ! The 4-dim location PPCC(2:5) is used to interpolate on the 4-dim V3D array ! (remember that ND3D is equivalenced to NDCC(2:5). ! Note that we allow bad values in V3D, but not in XC3D. PPCC(2) = float(I) ! Longitude index V3D = FLINT(-4,ND3D,VV3D,PPCC(2),Tiny)! Velocity in grid point I,J,K at time TL PPCC(1) = PRJ__XC ! Longitude shift index XCshift = FLINT(5,NDCC,XC3D,PPCC,Tiny)! Shift in XC from grid point I,J,K at time TL PPCC(1) = PRJ__TT ! Time shift index TTshift = FLINT(5,NDCC,XC3D,PPCC,Tiny)! Shift in TT from grid point I,J,K at time TL !------- ! Find the Carrington variable of origin at the source surface of grid ! point I,J,K at time TL. Note that XCbeg and XCrange are set to the ! values at time TL. XCbeg = XCbegL ! Needed for statement fnc XCvar XCend = XCendL XCrange = XCend-XCbeg PPSS(1) = XCvar(I)+XCshift ! Carrington variable at source surface !------- ! Find the time of origin at the source surface of grid point I,J,K,L, PPSS(3) = TL+TTshift ! Time at source surface PPSS(3) = TTindx(PPSS(3)) ! Index for time in source surface array !------- ! To convert the Carrington longitude to an array index, the values of ! XCbeg and XCrange need to be set to the proper values at time PP(4) ! (the time of origin at the source surface). This is done by ! interpolation on the XCbegMAT and XCendMAT arrays. XCbeg = FLINT(1,nTim,XCbegMAT,PPSS(3),Tiny)! Needed for statement fncs XCindx, XCmod XCend = FLINT(1,nTim,XCendMAT,PPSS(3),Tiny) if (XCbeg .eq. Bad .or. XCend .eq. Bad) then BR = Bad VV = Bad else XCrange = XCend-XCbeg if (bMod) PPSS(1) = XCmod(PPSS(1)) ! Map to [XCbeg, XCend] if bMod is set PPSS(1) = XCindx(PPSS(1)) ! Index for longitude in source surface array PPSS(2) = float(J) PP3D(1) = PPSS(1) PP3D(2) = PPSS(2) PP3D(3) = 1.0 ! Index for radius (= source surface) PP3D(4) = PPSS(3) !------ ! The velocity at the source surface is calculated by 4-dim linear ! interpolation at longitude P(1), latitude P(2)=J, radius PP(3)=1 ! and time P(4) VV = FLINT(-4,ND3D,VV3D,PP3D,Tiny) ! Velocity at source surface BR = FLINT(-3,NDSS,BRSS,PPSS,Tiny) ! r^2*Br at source surface if (BR .ne. Bad .and. VV .ne. Bad .and. V3D .ne. Bad) then VV = 0.5*(VV+V3D) VV = -BR*(VT*XL/VV) ! r*Bt end if end if iloc = locVOL(I,J,K,1) BR3D(iloc) = BR BT3D(iloc) = VV end do end do end do call T3D_write_bb('bbut',t3d, BR3D, BT3D, 0, ' ', cFile) ! Arg t3d is not used end do return end