C+ C NAME: C MkShift C PURPOSE: C Makes shifts at even heights to use in determining projections. C The shifts 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. C CATEGORY: C Data processing C CALLING SEQUENCE: subroutine MkShift(XCbeg,XCend,nLng,nLat,nMap,VV,DD,XC3D,Vratio, & Dratio,RR,dRR,CLRV,CLRD,VLL,VUL,NSIDE,Vtmp,Dtmp,XLtmp) C INPUTS: C nLng integer # longitude points C nLat integer # latitude points C nMap integer # heights (height begins at SUN__RAU AU) C VV(nLng,nLat) real Map of velocity values at a given height RR C DD(nLng,nLat) real Map of density values at a given height RR 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,2) real Internal scratch array C Dtmp(nLng,nLat) real Internal scratch array C XLtmp(nLng) real Internal scratch array C OUTPUTS: C XC3D(nLng,nLat,nMap) real XC3D contains the final shifts at all heights C in terms of fractions of a Carrington rotation C Vratio (nLng,nLat,nMap) real velocity ratios C Dratio (nLng,nLat,nMap) real density ratios (normalized density ratio???) C CALLS: C Say, ArrR4Copy2, ArrR4Zero, ArrR4Constant, ArrR4SetMinMax, FLINT, GridSphere2D C INCLUDE: C include 'mapcoordinates.h' C PROCEDURE: C > XC3D(I,J,N) is the shift needed to trace location (I,J,N) back to the reference C surface, i.e. XC = XCfrac(I)+XC3D(I,J,N) defines the origin of point (I,J,N) C at the reference surface in terms of a modified Carrington variable. C XC3D(I,J,N) will be negative above the reference surface (N>nRR), positive C below the reference surface (N Bad grid points in VV and DD are filled with 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/CASS); revision C- real XCbeg real XCend integer nLng integer nLat integer nMap real VV(nLng,nLat) ! Input velocity map at height RR real DD(nLng,nLat) ! Input density map at height RR real XC3D(nLng,nLat,nMap) ! Map of accumulated shifts at heights (from 1 RR) real Vratio (nLng,nLat,nMap) ! Map of accumulated velocity ratios at heights (from 1 RR) real Dratio (nLng,nLat,nMap) ! Map of accumulated density ratios at heights (from 1 RR) real RR real dRR real CLRV real CLRD real VLL real VUL integer*1 NSIDE(*)! 9*nLng*nLat) real Vtmp (nLng,nLat,2) ! Scratch arrays real Dtmp (nLng,nLat) real XLtmp(nLng) ! Include file contains statement functions, so must follow last declaration statement include 'mapcoordinates.h' Bad = BadR4() ALng = XCEnd-XCBeg !------- ! Store input arrays in Vtmp and Dtmp call ArrR4Copy2 ( nLngLat,VV,DD,Vtmp,Dtmp) ! Copy input arrays into scratch arrays call ArrR4SetMinMax(-nLngLat,Vtmp,VLL,VUL) ! Constrain Vtmp to [VLL,VUL] (exclude bad values) !------- ! in the main program calls to this subroutine are preceeded by INTERP calls. Hence there are ! no bad values at this point. The GridSphere2D calls should be redundant. call GridSphere2D(ALng,nLng,nLat,1,Vtmp,CLRV,4,0.0,0.0)! Fill in bad values call GridSphere2D(ALng,nLng,nLat,1,Dtmp,CLRD,4,0.0,0.0) call ArrR4Zero (nLngLat, XC3D(1,1,nRR)) ! Shifts at RR are 0 call ArrR4Constant(nLngLat,1.,Vratio (1,1,nRR)) ! V ratios at RR are 1. call ArrR4Constant(nLngLat,1.,Dratio (1,1,nRR)) ! D ratios at RR are 1. !------- ! The known velocities for the previous level (KM) are stored in Vtmp(*,*,KMV) ! 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 and KMV exchange values 1 and 2 KPV = 3-KMV dXC = SUN__SPIRAL*(RRhght(KP)-RRhght(KM)) ! Shift per (km/s)^(-1) from map KM to map KP Lmin = dXC/VUL*(ALng*nLng1+1) ! Minimum shift (based on maximum velocity) Lmax = dXC/VLL*(ALng*nLng1+2) ! Maximum shift (based on minimum velocity do J=1,nLat ! print *, 'J =', J do L=1,nLng ! Longitude L on level KM shifts to longitude XLtmp(L) at level KP XLtmp(L) = XCvar(XCfull(L)+dXC/Vtmp(L,J,KMV)) ! Shifted longitude on scale [1,nLng] end do do I=1,nLng ! Loop over all longitudes on level KP VW = 0. ! Accumulate weighted velocities W = 0. ! Accumulate weights do LL=I+Lmin,min(nLng,I+1+Lmax) ! Search to the right of I on level KM L = LL F = XLtmp(L)-I ! Distance between I and XLtmp(L) F = abs(F) if (F .lt. 1.) then ! Collect points within one grid spacing tmp = XCvar(XCfull(L)+XC3D(L,J,KM)) tmp = FLINT(1,nLng,Dtmp(1,J),tmp,0.0) ! Density at source surface if (tmp .ne. Bad) then ! Weigh with 1-F and with density at L on level KM F = (1.-F)*Dratio(L,J,KM)*tmp F = max(0.,F) ! Redundant if input densities are all >=0 VW = VW+F*Vtmp(L,J,KMV) ! Momentum flux: nvr^2 ! This assumes that Dratio*Dtmp is a normalized density !tmp = SW_Model_DevFromR2(Dratio(L,J,KM)*tmp,Vtmp(L,J,KMV)) !if (tmp .eq. 0.0) then W = W+F ! Mass: nr^2 = constant !else ! W = W+F*(RRhght(KM)/RRhght(KP))**tmp ! Mass nr^(2+delta) = constant !end if end if end if end do if (W .eq. 0) then Vtmp(I,J,KPV) = Bad else Vtmp(I,J,KPV) = VW/W end if end do end do !------- ! Fill holes and smooth V maps using half the binwidth as window. ! Then add the new shift to the old shift closer to the reference surface. ! Amount of new shift needed to get to the last location is ! call GridSphere2D(ALng,nLng,nLat,1,Vtmp(1,1,KPV),90./(nLat-1),4,0.0,90.0) ! Original call GridSphere2D(ALng,nLng,nLat,1,Vtmp(1,1,KPV),90./(nLat-1),3,0.0,90.0) ! Smoothing now used do J=1,nLat do I=1,nLng ! Loop over locations at level KP dXC0 = dXC/Vtmp(I,J,KPV) ! Shift needed to go back to level KM XL = XCvar(XCfull(I)-dXC0) ! Longitude on level KM on scale [1,nLng] ! Interpolate shifts at level KM if (XL .gt. nLng) XL = nLng XC3D(I,J,KP) = FLINT(1,nLng,XC3D(1,J,KM ),XL,0.0)-dXC0! Accumulated shift R = Vtmp(I,J,KPV)/FLINT(1,nLng,Vtmp (1,J,KMV),XL,0.0) ! Velocity ratio between levels Vratio(I,J,KP) = FLINT(1,nLng,Vratio (1,J,KM ),XL,0.0)*R ! Accumulated V ratio Dratio(I,J,KP) = FLINT(1,nLng,Dratio (1,J,KM ),XL,0.0)/R ! Accumulated D ratio end do end do end do return end