C+ C NAME: C MkShiftdnma C PURPOSE: C To make XC shift maps at even heights to use in determining projections. C The maps give the average longitude shift in fractions of a rotation C required to get to the input reference velocity map at each latitude C and longitude point on the map. This version of the subroutine allows C the shifts to be determined from input volumetric data, and goes from C the top down. This is unlike the shift maker in the kinematic tomography C program that needs to "invent" conservation of mass and mass flux as the C different height surfaces are stepped outward. This program assumes C an external program has run to provide the volumes that this program can use, C and that these volumes have conserved masss and mass flux, and are otherwise C complete with (if possible) accurate 3-component magnetic field data. C Magnetic field data does not need to begin at the tomographic source surface C height, which is constrained to be located at the lower level of the input C volumetric data, but it must begin below that height. C C CATEGORY: C Data processing C CALLING SEQUENCE: C call MkShiftdnma(iPrint,XCbe,XCtbeg,XCtend,ALng,nLng,nLat,nMap,nT, C VV3,BB3,DD1,TT1,XCshift3,XCshiftM,Vratio3,Bratio3,Dratio1,Tratio1,RR,RRMS,KPF, C FALLOFFN,FALLOFFT,FALLOFFBR,FALLOFFBT,FALLOFFBN,CLRV,CLRD,VLL,VUL,VLLLL,VLLUL, C Vtmp,Dtmp,XLLTtmp) C INPUTS: C iPrint integer if 1, print out lots of information C XCbe(2,nT) real Boundary values of time maps C XCtbeg real*8 Beginning time interval C XCtend real*8 Ending time interval C Alng real Rotations per nLng1 C nLng integer # longitude points C nLat integer # latitude points C nMap integer # heights (height begins at Sun__RAu AU) C nT integer # time intervals C VV3(nLng,nLat,NMap,nT,3) real Maps of velocity component values at all heights and times C BB3(nLng,nLat,NMap,nT,3) real Maps of magnetic field component values at all heights and times C DD1(nLng,nLat,NMap,nT) real Maps of density values at all heights and times C TT1(nLng,nLat,NMap,nT) real Maps of temperature values at all heights and times C RR real Distance above sun for reference velocity map C RRMS real Distance above sun for reference magnetic field map C KPF integer The height of the first tomographic inversion map. C If 1, the magnetic field map and the first traceback map have the same values. C FALLOFFN real Fall off in density C FALLOFFT real Fall off in temperature C FALLOFFBR real Fall off in radial field C FALLOFFBT real Fall off in tangential field C FALLOFFBN real Fall off in normal field C CLRV real Velocity map radius filter distance C CLRD real Density map radius filter distance C VLL real Lower limit on traceback velocity C VUL real Upper limit on traceback velocity C VLLLL real Lower limit on traceback longitude, latitude velocity C VLLUL real Upper limit on traceback longitude, latitude velocity C C Vtmp(nLng,nLat,nT,3) real Internal scratch array C Dtmp(nLng,nLat,nT) real Internal scratch array C XLLTtmp(nLng,nT,8) real Internal scratch array C OUTPUTS: C XCshift3(nLng,nLat,nMap,nT,3) real XCshift contains the final accumulated shifts (longitude at all heights C in terms of fractions of a Carrington rotation), latitude, and time C XCshiftM(nLng,nLat,nT,3) real shift amount increments from RRMS to RR C Vratio3 (nLng,nLat,nMap,nT,3) real velocity ratios C Bratio3 (nLng,nLat,nMap,nT,3) real mag ratios C Dratio1 (nLng,nLat,nMap,nT) real density ratios C Tratio1 (nLng,nLat,nMap,nT) real temperature change C C FUNCTIONS/SUBROUTINES: C PROCEDURE: C > XCshift(I,J,N,K,3) is the shift needed to trace location (I,J,N,K) back C to the reference surface at a given time, i.e. C XC = XCfrac(I)+XCshift(I,J,N,K,3) defines the origin of point (I,J,N,K,3) C at the reference surface in terms of a modified Carrington variable. C XCshift(I,J,N,K,3) will be negative above the reference surface (N>nRR), C zero at the reference surface (N=nRR) and will point to an earlier C or the same time. Bad grid points in VV and DD are filled with C GridSphere2D C MODIFICATION HISTORY: C MAR-1997, B. Jackson (UCSD) (original name was MAKE_TS) C MAY-1998, Paul Hick (UCSD; pphick@ucsd.edu); revision C APR-1999, Paul Hick, B. Jackson (UCSD); revision C APR-2001, B. Jackson (UCSD); revision C- subroutine MkShiftdnma(iPrint,XCbe,XCtbeg,XCtend,ALng,nLng,nLat,nMap,nT, & VV3,BB3,DD1,TT1,XCshift3,XCshiftM,Vratio3,Bratio3,Dratio1,Tratio1,RR,RRMS,dRR,dRRMS,KPF, & FALLOFFN,FALLOFFT,FALLOFFBR,FALLOFFBT,FALLOFFBN,CLRV,CLRD,VLL,VUL,VLLLL,VLLUL, & Vtmp,Dtmp,XLLTtmp) real VV3 (nLng,nLat,NMap,nT,3), ! Input volumetric velocities & BB3 (nLng,nLat,NMap,nT,3), ! Input volumetric magnetic fields & DD1 (nLng,nLat,NMap,nT), ! Input volumetric densities & TT1 (nLng,nLat,NMap,nT), ! Input volumetric temperatures & XCshift3(nLng,nLat,nMap,nT,3), ! Map of accumulated shifts at heights (from 1 RR) (in long, lat, time) & XCshiftM(nLng,nLat,nT,3), ! shift amount increments from RRMS to RR & Vratio3 (nLng,nLat,nMap,nT,3), ! Map of accumulated velocity ratios at heights (from 1 RR) & Bratio3 (nLng,nLat,nMap,nT,3), ! Map of accumulated magnetic field ratios at heights (from 1 RR) & Dratio1 (nLng,nLat,nMap,nT), ! Map of accumulated density ratios at heights (from 1 RR) & Tratio1 (nLng,nLat,nMap,nT), ! Map of accumulated temperature ratios at heights (from 1 RR) & XCbe (2,nT), ! Carrington time intervals & XLT (3), ! Internal array & BB33(3), ! Internal array & VV33(3), ! Internal array & Vtmp (nLng,nLat,nT,3), ! Scratch arrays & Dtmp (nLng,nLat,nT), & XLLTtmp (nLng,nLat,nT,8) real*8 XCtbeg,XCtend,XCtfull,X integer NDLT(3) ! Internal array include 'MAPCOORDINATES.H' ! Contains executable C statements, so must follow last declaration statement C iPrint = 0 C iPrint = 1 if(iPrint.eq.1) print *, ' ' print*, 'Into mkshiftdnma' print*, 'Here 0 ', nRR, iPrint if (nRR .lt. 1 .or. nRR .gt. nMap) call Say('MkShiftd','F','Illegal','reference distance') Bad = BadR4() C VEL1AU = 425.0 NR13 = 0 NR77 = 0 KPVOLD = 0 NDLT(1) = nLng NDLT(2) = nLat NDLT(3) = nT AKMPAUPD = 1.495979e8/86400. ! Km/(AU*Day) AKMPAUPI = AKMPAUPD ! (Km/(AU*Day)) C The statement above was changed on 9/2/2011 BVJ C C print *, 'float(nT-1), AKMPAUPI, Sun__Spiral ', float(nT-1), AKMPAUPI, Sun__Spiral C ERot__Fract = Sun__Spiral*Sun__SiderealP/Sun__SynodicP ! 68.21347545*25.38/27.2753 C (Rotates a little slower than inertial from Earth) ERot__Fract = AKMPAUPI/Sun__SiderealP ! What longitude has this material come from on the Sun? C print *, ' In MKshiftd - NT-1, XCtbeg, XCtend, AKMPAUPI, ERot__Fract', NT-1, XCtbeg, XCtend, C & AKMPAUPI, ERot__Fract C------- C------- C The known velocities for the previous level (KM) are stored in Vtmp(*,*,*,KMV) C The velocities to be computed for the current level (KP) are stored in Vtmp(*,*,*,KPV)\ C C The above was the old way. Now the deconvolution must be done at the first level present C in the input model which is the location the inversion needs to take place. This is at C the input model boundary, which must be defined as the source surface, and the traceback C for the inversion needs to refer to this location in each higher level of the traceback C matrix. This is now called the KPF or first KP level. This is the place all vratios are C determined to be 1.0. If this is not done the velocities and densities and other parameters C now derived from the input model are unknown, or will need to be extrapolated downward. C C Physically this surface must start below all of the IPS lines of sight, which if defined C to be 11.6 degrees elongation, about what is needed for IPS sources at 327 MHz, is at 0.2 AU. C C This scheme presents a problem for magnetic field data that are extrapolated outward from C the solar surface since the way to extrapolate these magnetic fields have no way to have C their evolution determined below the inversion surface. The approximate solution to this C is used here is to have a second portion of the matrix (between the magnetic field source C surface and the inversion surface) be given a traceback using velocities from the lowest C level of the modeled velocities. Magnetic fields are traced outward from the magnetic field C source surface approximately using the lowest modeled velocities, and then from there given C a traceback in all of the higher levels to the inversion surface. C C in the main program calls to this subroutine are preceeded by INTERP calls. Hence there are C no bad values at this point. The GridSphere2D calls should be redundant. C C The ratios on the first levels are only approximate C if(iPrint.eq.1) print *, ' ' if(iPrint.eq.1) write(*,'(A,I3)'),'The new shift tomographic traceback surface KPF =',KPF if(iPrint.eq.1) print *, ' ' do N=1,nT C call GridSphere2D(ALng,nLng,nLat,1,Vtmp(1,1,N,1),CLRV,4,0.0,0.0)! Fill in bad values C call GridSphere2D(ALng,nLng,nLat,1,Vtmp(1,1,N,2),CLRV,4,0.0,0.0) C call GridSphere2D(ALng,nLng,nLat,1,Vtmp(1,1,N,3),CLRV,4,0.0,0.0) C call GridSphere2D(ALng,nLng,nLat,1,Dtmp(1,1,N),CLRD,4,0.0,0.0) call ArrR4Zero (nLngLat, XCshift3(1,1,KPF,N,1)) ! Shifts at level KPF are 0 call ArrR4Zero (nLngLat, XCshift3(1,1,KPF,N,2)) ! Shifts at level KPF are 0 call ArrR4Zero (nLngLat, XCshift3(1,1,KPF,N,3)) ! Shifts at level KPF are 0 call ArrR4Constant(nLngLat,1.,Vratio3(1,1,KPF,N,1)) ! V ratios at KPF are 1. call ArrR4Constant(nLngLat,1.,Vratio3(1,1,KPF,N,2)) ! V ratios at KPF are 1. call ArrR4Constant(nLngLat,1.,Vratio3(1,1,KPF,N,3)) ! V ratios at KPF are 1. call ArrR4Constant(nLngLat,1.,Bratio3(1,1,KPF,N,1)) ! V ratios at KPF are 1. call ArrR4Constant(nLngLat,1.,Bratio3(1,1,KPF,N,2)) ! V ratios at KPF are 1. call ArrR4Constant(nLngLat,1.,Bratio3(1,1,KPF,N,3)) ! V ratios at KPF are 1. call ArrR4Constant(nLngLat,1.,Dratio1(1,1,KPF,N)) ! D ratios at KPF are 1. call ArrR4Constant(nLngLat,1.,Tratio1(1,1,KPF,N)) ! T ratios at KPF are 1. end do C C In this do loop, deal with the shifts from each layer. C if(iPrint.eq.1) write (*,'(A,2F15.5)') , 'ERot__Fract, AKMPAUPI',ERot__Fract, AKMPAUPI nlat1 = nLat-1 nLng1 = nLng-1 nMap1 = nMap-1 nT1 = nT1 -1 KP = 1 C do K=1,nMap+1 do K=1,nMap if(K.eq.1) then KP = 1 KKK = KPF end if if(K.eq.KPF+1) then KP=1 ! This begins the accumulation of shifts at the second level of the shift matrix end if KKK = KP + KPF - 1 KM = KP ! Height index of previous map KP = KP + 1 ! Height index of current map if(iPrint.eq.1) print *, ' ' if(iPrint.eq.1) write (*,'(A,5I5)') , 'KPF, K, KP, KM, KKK',KPF, K, KP, KM, KKK do N=1,nT do J=1,nLat do I=1,nLng C if(K.eq.1.or.K.eq.KPF.or.K.eq.KPF+1) XCShift3(I,J,KKK,N,1) = 0.0 if(K.le.KPF+1) XCShift3(I,J,KKK,N,1) = 0.0 XLLTtmp(I,J,N,1) = XCShift3(I,J,KKK,N,1) ! Longitude C if(K.eq.1.or.K.eq.KPF.or.K.eq.KPF+1) XCShift3(I,J,KKK,N,2) = 0.0 if(K.le.KPF+1) XCShift3(I,J,KKK,N,2) = 0.0 XLLTtmp(I,J,N,2) = XCShift3(I,J,KKK,N,2) ! Latitude C if(K.eq.1.or.K.eq.KPF.or.K.eq.KPF+1) XCShift3(I,J,KKK,N,3) = 0.0 if(K.le.KPF+1) XCShift3(I,J,KKK,N,3) = 0.0 XLLTtmp(I,J,N,3) = XCShift3(I,J,KKK,N,3) ! Time end do end do end do if(iPrint.eq.1) write(*,'(A,3F8.4))'), & 'N=22,I=20,J=6, XLtp ',XLLTtmp(20,6,22,1),XLLTtmp(20,6,22,2),XLLTtmp(20,6,22,3) if(K.le.KPF) then C dXCL = ERot__Fract*(RRhght(KP-KPF)-RRMS) ! RRMS is the magnetic field source surface C dXCT = AKMPAUPI*(RRhght(KP-KPF)-RRMS) dXCL = ERot__Fract*dRRMS ! dRRMS is the magnetic field source surface distance delta dXCT = AKMPAUPI*dRRMS else dXCL = ERot__Fract*(RRhght(KP)-RRhght(KM)) ! Long. shift const. /(km/s)^(-1) from map KM to map KP without VT dXCT = AKMPAUPI*(RRhght(KP)-RRhght(KM)) ! Time shift const. per(km/s)^(-1) from map KM to map KP end if ALATA = atan(dXCT/RRhght(KP)) if(iPrint.eq.1) write (*,'(A,2I3,2F8.4,2F8.3)') , 'KP,KM,RRh(KP),RRh(KM),dXCL,dXCT ',KP,KM,RRhght(KP),RRhght(KM),dXCL,dXCT do N=1,nT do J=1,nLat ALAT = 90.0*(-1.0+2.0*(J-1)/nLat1) ! Shift in longitude depends on latitude do I=1,nLng if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6) write (*,'(A,2I3,2F8.4,2F8.3)'), & 'KP,KM,RRh(KP),RRh(KM),ALAT (latitude) ',KP,KM,RRhght(KP),RRhght(KM),ALAT VV331 = VV3(I,J,N,KP,1) VV332 = VV3(I,J,N,KP,2) VV333 = VV3(I,J,N,KP,3) if(K.lt.KPF) then VV331 = VV3(I,J,N,K,1) VV332 = VV3(I,J,N,K,2) VV333 = VV3(I,J,N,K,3) end if if(VV332.eq.0.0) then ! Shift in Longitude for tangential velocity dXCTLOO = 0.0 else if(ALAT.ne.-90.0.or.ALAT.ne.90.0) then ALONGA = ((dXCT/VV332)/RRhght(KP))/cos(ALAT) else ALONGA = 0.0 end if dXCTLOO = ALONGA end if if(VV333.eq.0.0) then ! Shift in Latitude for normal velocity dXCTLAA = 0.0 else dXCTLAA = atan((dXCT/VV333)/RRhght(KP)) end if if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6) write (*,'(A,4I3,2F8.4,2F8.3)'), & 'KP,N,J,I,RRh(KP),RRh(KM),dXLO,dXLA',K,N,J,I,RRhght(KP),RRhght(KM),dXCTLOO,dXCTLAA dXC0L = dXCL/VV331 - dXCTLOO ! XC Shift needed to go back to level KM in longitude dXC0LA = dXCLAA ! Deg Shift needed to go back to level KM in latitude dXC0T = dXCT/VV331 ! Time Shift needed to go back to level KM in time if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6) write(*,'(A,4I4,3F8.1))'), & 'KP,N,J,I,VV331,VV332,VV333',KP,N,J,I,VV331,VV332,VV333 if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6) write(*,'(A,4I4,3F8.4))'), & 'KP, N, J, I, dL, dLA, dT ',KP,N,J,I,dXC0L,dXC0LA,dXC0T XLT(3) = max(1.,XCtvar(XCtfull(N) - dXC0T)) ! Time index on level KM on scale [1,nT] ! Interpolate shifts at level KM II = Itvar(XLT(3)) ! Nearest interval time index XCbeg = XCbe(1,II) XCend = XCbe(2,II) XLT(1) = min(1.*nLng,XCvar(XCfull(I) - dXC0L)) ! Carrington variable -> index [1,nLng] ! Interpolate shifts at level KM XLT(2) = XLindx(XLdeg(J)-dXC0LA) ! Latitude (deg) -> index in [1,nLat] XLT(2) = min(float(nLat),XLT(2)) XLT(2) = max(1.,XLT(2)) ! Interpolate shifts at level KM XCshift3(I,J,K,N,1) = FLINT(3,NDLT,XLLTtmp(1,1,1,1),XLT,0.00002) - dXC0L ! Accumulated longitude shift XCshift3(I,J,K,N,2) = FLINT(3,NDLT,XLLTtmp(1,1,1,2),XLT,0.00002) - dXC0LA ! Accumulated latitude shift XCshift3(I,J,K,N,3) = FLINT(3,NDLT,XLLTtmp(1,1,1,3),XLT,0.00002) - dXC0T ! Accumulated time shift if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6) then print *, 'After Flint' write(*,'(A,4I4,3F8.4))'), & 'KP, N, J, I, XL1, XL2, Xl3',KP,N,J,I,XLT(1),XLT(2),XLT(3) write(*,'(A,4I4,3F8.4))'), & 'KP, N, J, I, FL1, FL2, FL3',KP,N,J,I,FLINT(3,NDLT,XLLTtmp(1,1,1,1),XLT,0.00002), & FLINT(3,NDLT,XLLTtmp(1,1,1,2),XLT,0.00002),FLINT(3,NDLT,XLLTtmp(1,1,1,3),XLT,0.00002) end if if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6) then write(*,'(A,4I4,3F8.4))'), & 'K, N, J, I, XS1, XS2, XS3',K,N,J,I,XCshift3(I,J,K,N,1),XCshift3(I,J,K,N,2),XCshift3(I,J,K,N,3) print *, 'End T Loop' end if end do end do end do end do C C In this do loop, deal with the ratios for each parameter in the input model. Ratios at KPF are 1.0 C if(iPrint.eq.1) print *, ' ' if(iPrint.eq.1) print *, 'Now begin second loop to provide ratios' if(iPrint.eq.1) print *, ' ' do N=1,nT do J=1,nLat do I=1,nLng XLLTtmp(I,J,N,1) = VV3(I,J,KPF,N,1) ! Base Velocity XLLTtmp(I,J,N,2) = VV3(I,J,KPF,N,2) XLLTtmp(I,J,N,3) = VV3(I,J,KPF,N,3) XLLTtmp(I,J,N,4) = BB3(I,J,1,N,1) ! Base Magnetic Field XLLTtmp(I,J,N,5) = BB3(I,J,1,N,2) XLLTtmp(I,J,N,6) = BB3(I,J,1,N,3) XLLTtmp(I,J,N,7) = DD1(I,J,KPF,N) ! Base Density XLLTtmp(I,J,N,8) = TT1(I,J,KPF,N) ! Base Temperature end do end do end do KP = 1 C do K=1,nMap+1 do K=1,nMap if(K.lt.KPF) then KP = 0 end if if(K.eq.KPF) then KP = 0 end if KP = KP + 1 BFR = (RRhght(KP)/RRhght(1))**FALLOFFBR BFT = (RRhght(KP)/RRhght(1))**FALLOFFBT BFN = (RRhght(KP)/RRhght(1))**FALLOFFBN DF = (RRhght(KP)/RRhght(1))**FALLOFFN TF = (RRhght(KP)/RRhght(1))**FALLOFFT do N=1,nT do J=1,nLat do I=1,nLng dXC0L = XCshift3(I,J,KP+KPF-1,N,1) ! Shift needed to go back to level 1 in longitude dXC0LA = XCshift3(I,J,KP+KPF-1,N,2) ! Shift needed to go back to level 1 in latitude dXC0T = XCshift3(I,J,KP+KPF-1,N,3) ! Shift needed to go back to level 1 in time if(iPrint.eq.1.and.N.eq.22.and.J.eq.6.and.I.eq.20) then print *, ' ' write(*,'(A,5I4,3F8.4))'), & 'K, KP, N, J, I, XS1, XS2, XS3',K,KP,N,J,I,dXC0L,dXC0LA,dXC0T end if if(iPrint.eq.1.and.N.eq.22.and.J.eq.6.and.I.eq.20) then write(*,'(A,5I4,3F8.4))'), & 'K, KP, N, J, I, XS1, XS2, XS3',K,KP,N,J,I,XCshift3(I,J,KP,N,1),XCshift3(I,J,KP,N,2),XCshift3(I,J,KP,N,3) end if XLT(3) = max(1.,XCtvar(XCtfull(N) - dXC0T)) ! Time index on level KP on scale [1,nT] ! Interpolate shifts at level KP II = Itvar(XLT(3)) ! Nearest interval time index XCbeg = XCbe(1,II) XCend = XCbe(2,II) XLT(1) = min(1.*nLng,XCvar(XCfull(I) - dXC0L)) ! Carrington variable -> index [1,nLng] ! Interpolate shifts at level KP XLT(2) = XLindx(XLdeg(J)-dXC0LA) ! Latitude (deg) -> index in [1,nLat] XLT(2) = min(float(nLat),XLT(2)) XLT(2) = max(1.,XLT(2)) ! Interpolate shifts at level KP if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6) then write(*,'(A,5I4,3F8.4))'), & 'K, KP, N, J, I, XS1, XS2, XS3',K,KP,N,J,I,dXC0L,dXC0LA,dXC0T write(*,'(A,5I4,3F8.4))'), & 'K, KP, N, J, I, XL1, XL2, XL3',K,KP,N,J,I,XLT(1),XLT(2),XLT(3) write(*,'(A,5I4,2F8.1,F12.1))'), & 'K, KP, N, J, I, XT1, XT7, XT8',K,KP,N,J,I,XLLTtmp(I,J,N,1),XLLTtmp(I,J,N,7),XLLTtmp(I,J,N,8) end if VV33(1) = FLINT(3,NDLT,XLLTtmp(1,1,1,1),XLT,0.00002) VV33(2) = FLINT(3,NDLT,XLLTtmp(1,1,1,2),XLT,0.00002) VV33(3) = FLINT(3,NDLT,XLLTtmp(1,1,1,3),XLT,0.00002) DD11 = FLINT(3,NDLT,XLLTtmp(1,1,1,7),XLT,0.00002) TT11 = FLINT(3,NDLT,XLLTtmp(1,1,1,8),XLT,0.00002) C call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,XLLTtmp(1,1,1,4),3, ! Interpolated radial field from level 1 C & dXC0L,dXC0LA,dXC0T,BB33) C Bratio3(I,J,K,N,1) = BB3(I,J,K,N,1)*BFR/BB33(1) C call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,XLLTtmp(1,1,1,5),3, ! Interpolated tangential field from level 1 C & dXC0L,dXC0LA,dXC0T,BB33) C Bratio3(I,J,K,N,2) = BB3(I,J,K,N,2)*BFT/BB33(2) C C call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,XLLTtmp(1,1,1,6),3, ! Interpolated normal field from level 1 C & dXC0L,dXC0LA,dXC0T,BB33) C Bratio3(I,J,K,N,3) = BB3(I,J,K,N,3)*BFN/BB33(3) C call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,XLLTtmp(1,1,1,1),1, ! Interpolated radial velocity from level KPF C & dXC0L,dXC0LA,dXC0T,VV33) Vratio3(I,J,K,N,1) = VV3(I,J,K,N,1)/VV33(1) C call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,XLLTtmp(1,1,1,2),3, ! Interpolated tangential velocity from level KPF C & dXC0L,dXC0LA,dXC0T,VV33) C Vratio3(I,J,K,N,2) = VV3(I,J,K,N,2)/VV33(2) C call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,XLLTtmp(1,1,1,3),3, ! Interpolated normal velocity from level KPF C & dXC0L,dXC0LA,dXC0T,VV33) C Vratio3(I,J,K,N,3) = VV3(I,J,K,N,3)/VV33(3) C call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,XLLTtmp(1,1,1,7),1, ! Interpolated density from level KPF C & dXC0L,dXC0LA,dXC0T,DD11) Dratio1(I,J,K,N) = DD1(I,J,K,N)*DF/DD11 C call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,XLLTtmp(1,1,1,8),1, ! Interpolated temperature from level KPF C & dXC0L,dXC0LA,dXC0T,TT11) Tratio1(I,J,K,N) = TT1(I,J,K,N)*TF/TT11 if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6) write(*,'(A,I3,4I4,2F8.5))'), & 'K, KP, N, J, I, VRatio, DRatio',K,KP,N,J,I,Vratio3(I,J,K,N,1),Dratio1(I,J,K,N) end do end do end do end do C C Ratios and the values of magnetic field at the boundary need to be delt with seperately below KPF. C if(KPF.eq.2000) then ! The tomographic map is inverted at the same level as the first magnetic map do N=1,nT if(N.eq.KPF) then call ArrR4Constant(nLngLat,1.,Vratio3(1,1,KPF,N,1)) ! V ratios at KPF are 1. call ArrR4Constant(nLngLat,1.,Vratio3(1,1,KPF,N,2)) ! V ratios at KPF are 1. call ArrR4Constant(nLngLat,1.,Vratio3(1,1,KPF,N,3)) ! V ratios at KPF are 1. call ArrR4Constant(nLngLat,1.,Bratio3(1,1,KPF,N,1)) ! Mag R ratios at KPF are 1. call ArrR4Constant(nLngLat,1.,Bratio3(1,1,KPF,N,2)) ! Mag T ratios at KPF are 1. call ArrR4Constant(nLngLat,1.,Bratio3(1,1,KPF,N,3)) ! Mag N ratios at KPF are 1. call ArrR4Constant(nLngLat,1.,Dratio1(1,1,KPF,N)) ! D ratios at KPF are 1. call ArrR4Constant(nLngLat,1.,Tratio1(1,1,KPF,N)) ! T ratios at KPF are 1. end if do J=1,nLat do I=1,nLng VV3 (I,J,1,N,1) = VV3 (I,J,2,N,1) VV3 (I,J,1,N,2) = VV3 (I,J,2,N,2) VV3 (I,J,1,N,3) = VV3 (I,J,2,N,3) BB3 (I,J,1,N,1) = BB3 (I,J,2,N,1) BB3 (I,J,1,N,2) = BB3 (I,J,2,N,2) BB3 (I,J,1,N,3) = BB3 (I,J,2,N,3) DD1 (I,J,1,N) = DD1 (I,J,2,N) TT1 (I,J,1,N) = TT1 (I,J,2,N) if(iPrint.eq.1.and.I.eq.20.and.J.eq.6) write(*,'(A,3I4,2F12.1))'), & 'N , J, I, VV3 DD1 ',N,J,I,VV3(I,J,2,N,1),DD1 (I,J,2,N) end do end do end do end if return end