C+ C NAME: C MkShiftd n 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 underwent C considerable modification on about 4/1/02 when the way of determining ratios C were not accumulated with each height, but were accumulated relative to the base C at the source surface. This version, which is where the majority of the C time is spent in the tomography program, is more than twice as slow as the C old version of the mkshiftd subroutine. C CATEGORY: C Data processing C CALLING SEQUENCE: C call MkShiftdn(XCbe,XCtbeg,XCtend,ALng,nLng,nLat,nMap,nT, C VV,DD,XCshift,Rconst,Vratio,Dratio,RR,dRR,CLRV,CLRD,VLL,VUL, C NSIDE,Vtmp,Dtmp,VZtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt) C INPUTS: C XCbe(2,nT) real Boundary values of time maps C XCtbeg real Beginning time interval C XCtend real 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 VV(nLng,nLat,nT)real Map of velocity values at a given height RR C DD(nLng,nLat,nT)real Map of density values at a given height RR C Rconst real ~25.38/27.275, but specific to this interval C RR real Distance above sun for reference velocity map C dRR real Distance between individual maps (in AU) 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 C Vtmp(nLng,nLat,nT,3) real Internal scratch array C Dtmp(nLng,nLat,nT) real Internal scratch array C VZtmp(nLng,nLat,nT) real Internal scratch array C XLTtmp(nLng,nT,2) real Internal scratch array C XLtmp(nLng,nT) real Internal scratch array C XLtmpt(nLng,nT) real Internal scratch array C Dtmpt(nT) real Internal scratch array C OUTPUTS: C XCshift(nLng,nLat,nMap,nT,3) real XCshift contains the final shifts at all heights C in terms of fractions of a Carrington rotation C Vratio (nLng,nLat,nMap,nT) real velocity ratios C Dratio (nLng,nLat,nMap,nT) real density ratios 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 MkShiftdn(XCbe,XCtbeg,XCtend,ALng,nLng,nLat,nMap,nT, & VV,DD,XCshift,Rconst,Vratio,Dratio,RR,dRR,CLRV,CLRD,VLL,VUL, & NSIDE,Vtmp,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt) real VV(nLng,nLat,nT), ! Input velocity map at height RR & DD(nLng,nLat,nT), ! Input density map at height RR & XCshift(nLng,nLat,nMap,nT,3), ! Map of accumulated shifts at heights (from 1 RR) & Vratio (nLng,nLat,nMap,nT), ! Map of accumulated velocity ratios at heights (from 1 RR) & Dratio (nLng,nLat,nMap,nT), ! Map of accumulated density ratios at heights (from 1 RR) & XCbe (2,nT), ! Carrington time intervals & XLT(2), ! Internal array & Vtmp (nLng,nLat,nT,3), ! Scratch arrays & Dtmp (nLng,nLat,nT), & XLTtmp(nLng,nT,2), & XLtmp (nLng,nT), & XLtmpt(nLng,nT), & Dtmpt (nLng,nT) integer NDLT(2) ! Internal array byte NSIDE(*) ! (9*nLng*nLat) include 'MAPCOORDINATES.H' ! Contains executable C statements, so must follow last declaration statement if (nRR .lt. 1 .or. nRR .gt. nMap) call Say('MkShiftd','F','Illegal','reference distance') Bad = BadR4() VEL1AU = 425.0 C print*, 'Into this mkshiftdn' NR13 = 0 NR77 = 0 KPVOLD = 0 NDLT(1) = nLng NDLT(2) = nT ANt = XCtend-XCtbeg C print *, ' In MKshiftd - XCtbeg XCtend, Rconst', XCtbeg, XCtend , Rconst C------- C Store input arrays in Vtmp and Dtmp call ArrR4Copy2 ( nLngLatnt,VV,DD,Vtmp,Dtmp) ! Copy input arrays into scratch arrays call ArrR4SetMinMax(-nLngLatnt,Vtmp,VLL,VUL) ! Constrain Vtmp to [VLL,VUL] (exclude bad values) 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. do N=1,nT call GridSphere2D(ALng,nLng,nLat,1,Vtmp(1,1,N,1),CLRV,4,0.0,0.0)! Fill in bad values call ArrR4Copy (nLngLat,Vtmp(1,1,N,1),Vtmp(1,1,N,3)) call GridSphere2D(ALng,nLng,nLat,1,Dtmp(1,1,N),CLRD,4,0.0,0.0) call ArrR4Zero (nLngLat, XCshift(1,1,nRR,N,1)) ! Shifts at RR are 0 call ArrR4Zero (nLngLat, XCshift(1,1,nRR,N,2)) ! Shifts at RR are 0 call ArrR4Zero (nLngLat, XCshift(1,1,nRR,N,3)) ! Shifts at RR are 0 call ArrR4Constant(nLngLat,1.,Vratio (1,1,nRR,N)) ! V ratios at RR are 1. call ArrR4Constant(nLngLat,1.,Dratio (1,1,nRR,N)) ! D ratios at RR are 1. end do 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) KPV = 1 do KP=nRR+1,nMap KM = KP-1 ! Height index of previous map KMV = KPV KPV = 3-KMV dXCL = SUN__SPIRAL*(RRhght(KP)-RRhght(KM)) ! Shift per (km/s)^(-1) from map KM to map KP dXCT = dXCL*Rconst ! Shift per(km/s)^(-1) from map KM to map KP in time Lmin = dXCL/VUL*(nLng1/ALng) ! Minimum longitude shift (based on maximum velocity) Lmax = 1 + dXCL/VLL*(nLng1/ALng) ! Maximum longitude shift (based on minimum velocity) Lmint = dXCT/VUL*((nT-1)/ANt) ! Minimum time shift (based on maximum velocity) Lmaxt = 1 + dXCT/VLL*((nT-1)/ANt) ! Maximum time shift (based on minimum velocity) C print *, 'Lmin, Lmax, Lmint, Lmaxt =',Lmin, Lmax, Lmint, Lmaxt do J=1,nLat do N=1,nt XCbeg = XCbe(1,N) XCend = XCbe(2,N) do L=1,nLng C C Longitude at longitude L at time N on level KM shifts to C longitude XLtmp(L,N) at level KP C XLtmp(L,N) = XCvar(XCfull(L)+dXCL/Vtmp(L,J,N,KMV)) ! Shifted longitude on scale [1,nLng] C C Time at longitude L at time N on level KM shifts to time C XLtmpt(L,N) at level KP C XLtmpt(L,N) = XCtvar(XCtfull(N)+dXCT/Vtmp(L,J,N,KMV)) ! Shifted time on scale [1,nT] Dtmpt(L,N) = Dtmp(L,J,N) end do end do do I=1,nLng ! Loop over all longitudes on level KP do N=1,nT ! Loop over all times on level KP VW = 0. ! Accumulate longitude weighted velocities W = 0. ! Accumulate longitude weights C C Calculate the weighted velocities and their weights at each longitude C do L=I+Lmin,min(nLng,I+1+Lmax) ! Search to the right of I on level KM at time N do M=N+Lmint,min(nT,N+1+Lmaxt) ! Search to lower times of N on level KM at longitude I F = XLtmp(L,M)-I ! Distance between I and XLtmp(L,N) F = abs(F) if (F .lt. 1.) then ! Collect points within one grid spacing Ft = XLtmpt(L,M)-N ! Distance between N and XLtmpt(I,M) Ft = abs(Ft) if (Ft .lt. 1.) then ! Collect points within one grid spacing XCbeg = XCbe(1,M) XCend = XCbe(2,M) XLT(1) = XCvar(XCfull(L)+XCshift(L,J,KM,M,1)) XLT(2) = XCtvar(XCtfull(M)+XCshift(L,J,KM,M,3)) C print *, ' ' C print *, ' Before FLINT - XLtmpt XLT(1) XLT(2) XCfull XCtFull XCshift C & L J KM M', XLtmpt(L,M), XLT(1), XLT(2), XCfull(L), XCtfull(M), XCshift(L,J,KM,M), L, J, KM, M XCV = FLINT(2,NDLT,Dtmpt,XLT,0.00002) if (XCV .ne. Bad) then F = (1.-Ft)*(1.-F)*Dratio(L,J,KM,M)*XCV ! Weight with 1-F and with density at M on level KM F = max(0.,F) ! Redundant if input densities are all >=0 VW = VW + F*Vtmp(L,J,M,KMV) W = W + F else ! Density beyond time (or longitude) domain F = (1.-Ft)*(1.-F)*1.0 VW = VW + F*Vtmp(L,J,M,KMV) W = W + F end if C print *, ' After FLINTt - XCV F I J KM M VW W', XCV, F, I, J, KM, M, VW, W end if end if end do end do if (W .eq. 0) then Vtmp(I,J,N,KPV) = Bad else Vtmp(I,J,N,KPV) = VW/W end if C print *, 'V, I J N KPV', Vtmp(I,J,N,KPV), I, J, N, KPV end do end do end do C------- C Fill holes and smooth V maps using half the binwidth as window. C Then add the new shift to the old shift closer to the reference surface. C Amount of new shift needed to get to the last location is C print *, 'Through first loop of nMap', KP C C Put in in 2000 (BVJ, UCSD/CASS) C NTT = nt NTT2 = NTT/2 do N=NTT2,NTT-1 NNM = NTT - N NNP = N + 1 call ArrR4getminmax(nLngLat,Vtmp(1,1,NNM,KPV),aminm,amaxm) call ArrR4getminmax(nLngLat,Vtmp(1,1,NNP,KPV),aminp,amaxp) if(amaxp.eq.Bad.and.N.eq.NTT2) then print *, 'Bad middle Vtmp in MkShift, Replace with constant velocity' call ArrR4Constant(nLngLat,VEL1AU,Vtmp(1,1,NNM+1,KPV)) end if if(amaxm.eq.Bad) then print *, 'All bad values in Vtmp in MkShift NNM =',NNM,' at Ht. =',KP call ArrR4Copy(nLngLat,Vtmp(1,1,NNM+1,KPV),Vtmp(1,1,NNM,KPV)) end if if(amaxp.eq.Bad) then print *, 'All bad values in Vtmp in MkShift NNP =',NNP,' at Ht. =',KP call ArrR4Copy(nLngLat,Vtmp(1,1,NNP-1,KPV),Vtmp(1,1,NNP,KPV)) end if end do do N=1,nt call GridSphere2D(ALng,nLng,nLat,1,Vtmp(1,1,N,KPV),90./(nLat-1),4,0.0,90.0) ! Simply fill the holes end do do J=1,nLat do N=1,nt do I=1,nLng XLTtmp(I,N,1) = XCShift(I,J,KM,N,3) XLTtmp(I,N,2) = Vtmp ( I,J,N,3) end do end do do N=1,nT XCbeg = XCbe(1,N) XCend = XCbe(2,N) do I=1,nLng ! Loop over locations at level KP dXC0L = dXCL/Vtmp(I,J,N,KPV) ! Shift needed to go back to level KM in longitude dXC0T = dXCT/Vtmp(I,J,N,KPV) ! Shift needed to go back to level KM in time XLT(1) = min(1.*nLng,XCvar(XCfull(I)-dXC0L)) ! Longitude on level KM on scale [1,nLng] ! Interpolate shifts at level KM XLT(2) = max(1.,XCtvar(XCtfull(N)-dXC0T)) ! Time on level KM on scale [1,nT] ! Interpolate shifts at level KM C if(J.eq.1) print *, ' Begin' XCshift(I,J,KP,N,1) = FLINT(2,NDLT,XLTtmp(1,1,1),XLT,0.00002)-dXC0L ! Accumulated longitude shift XCshift(I,J,KP,N,2) = 0.0 ! Accumulated latitude shift XCshift(I,J,KP,N,3) = FLINT(2,NDLT,XLTtmp(1,1,1),XLT,0.00002)-dXC0T ! Accumulated time shift dXCshiftL = XCshift(I,J,KP,N,1) ! Shift needed to go back to base dXCshiftT = XCshift(I,J,KP,N,3) ! Shift needed to go back to base in time XLT(1) = min(1.*nLng,XCvar(XCfull(I)+dXCshiftL)) ! Longitude on base level on scale [1,nLng] XLT(2) = max(1.,XCtvar(XCtfull(N)+dXCshiftT)) ! Time on base level on scale [1,nT] R = Vtmp(I,J,N,KPV)/FLINT(2,NDLT,XLTtmp(1,1,2),XLT,0.00002) ! Velocity ratio between this level and that ! projected from the base level at this height Vratio(I,J,KP,N) = R ! Accumulated V ratio on this level Dratio(I,J,KP,N) = 1.0/R ! Accumulated D ratio this level if mv = constant end do end do end do end do return end