C+ C NAME: C shift_MOD C PURPOSE: C To accumulate shift values at given heights and times to use in determining C model projections to a lower source surface given the velocities at each point C in time. The shifts give the average longitude, latitude and time shifts in C fractions of a Carrington rotation required to get to the input reference C map at each longitude, latitude, height and time point. C CATEGORY: C Data processing C CALLING SEQUENCE: C call shift_MOD(NLNGM,NLATM,NMAPM,NTM,V3,D3,V3L,D3L,R_MOD,IYR,DOY,IYRL,DOYL,ICTX,FICTX,ICTXL,FICTXL, C V3D,D3D,V3DL,D3DL,V3DTMP,D3DTMP, C nCar,JDCar,NCoff,xcbe,xctbeg,xctend,aLng,nLng,nLat,nMap,nT,rr,drr,Rconst, C cLrV,cLrD,vTmp,dTmp,vv,dd,xcshift,vRatio,dRatio) C INPUTS: C THESE ARE MODEL INPUTS C NLNGM integer # of model longitude points (inertial frame includes both 0 and 360 deg.) C NLATM integer # of model latitude points (inertial frame from -90 to + 90 deg.) C NMAPM integer # of model heights (heights must go from same as tomography source surface C up to at least 3.0 AU) C V3(NLNGM,NLATM,NMAPM,3) real 3D velocity (west +, north +, outward +) input from model in km/sec at time IYR + DOY C D3(NLNGM,NLATM,NMAPM) real 3D density input from model in particles/cm^3 at time IYR + DOY C V3L(NLNGM,NLATM,NMAPM,3) real Save space for last 3D velocity input from model in km/sec at time IYRL + DOYL C D3L(NLNGM,NLATM,NMAPM) real Save space for last 3D density input from model in particles/cm^3 at time IYRL + DOYL C R_MOD(NMAPM) real Heights of the model calculation C IYR integer Year of the current model calculation (e.g., 20XX) (e.g., 2000) C DOY real Day of year for the current model calculation (begins at 1.0 on the first day) C IYRL integer Save space for the last year of the model calculation. C DOYL real Save space for the last day of year for the model calculation. C ICTX integer Beginning Carrington variable source surface at this time. (e.g., 1965) C FICTX real Beginning Carrington variable source surface fraction at this time. (e.g., 0.34567) C ICTXL integer Save space for the beginning Carrington variable source surface at the last time. C FICTXL real Save space for the beginning Carrington variable source surface fraction the last time. C Note - The save space variables need not be set, just allocated a space. C V3D((NLNGM-1)*3+1,NLATM,NMAPM,3) real Scratch arrays for model from calling program C D3D((NLNGM-1)*3+1,NLATM,NMAPM) real C V3DL((NLNGM-1)*3+1,NLATM,NMAPM,3) real C D3DL((NLNGM-1)*3+1,NLATM,NMAPM) real C V3DTMP(nLngt,nLatt,2,2,3) real C D3DTMP(nLngt,nLatt,2,2,3) real C C THESE ARE TOMOGRAPHY INPUTS AND THEIR VALUES C nCar integer # Carrington rotations C JDCar real*8 Julian Date at beginning of rotations C NCoff integer Carrington variable offset C xcbe(2,nT) real Tomography boundary values of time maps (passed through) C xctbeg real Tomography beginning time interval (passed through) C xctend real Tomography ending time interval (passed through) C aLng real Rotations per nLng1 (passed through) C nLng integer # tomography model longitude points (passed through) C nLat integer # tomography model latitude points (passed through) C nMap integer # tomography model tomographic heights (passed through) C nT integer # tomography model time intervals (passed through) C rr real Height of tomographic source surface (in AU). (must be identical to model value) C drr real Distance between points in tomographic height array (passed through) C Rconst real Tomography time constant ~25.38/27.275, but specific to the interval. C cLrV real Velocity map radius filter distance (in interval fractions) (passed through) C cLrD real Density map radius filter distance (in interval fractions) (passed through) C vTmp(nLng,nLat) real Internal scratch array C dTmp(nLng,nLat) real Internal scratch array C OUTPUTS: C THESE ARE OUTPUTS TO THE TOMOGRAPHY PROGRAM C vv(nLng,nLat,nT) real Map of velocity values at a given source surface height RR C dd(nLng,nLat,nT) real Map of density values at a given source surface height RR C xcshift(nLng,nLat,nMap,nT,3) real xcshift contains the final shifts for the tomography program C at the times and heights it needs in the coordinates it needs. C vRatio (nLng,nLat,nMap,nT) real velocity ratios C dRatio (nLng,nLat,nMap,nT) real density ratios C C FUNCTIONS/SUBROUTINES: C ArrR4Copy2, GridSphere2D, ArrR4Zero, ArrR4Constant, FLINT, get3Dtval C PROCEDURE: C > xcshift(i,j,n,itim,l) is the shift needed to trace location (i,j,n,itim,3) back C to the reference surface at a given time, i.e.: C xcshift(i,j,n,itim,1) defines the location of the origin of the point in longitude C at the reference surface in terms of a modified Carrington variable. C xcshift(i,j,n,itim,1) (longitude) will be negative above the reference C surface (n>1). C xcshift(i,j,n,itim,2) defines the location of the origin of the point in latitude C at the reference surface in terms of degrees. C xcshift(i,j,n,itim,3) defines the location of the origin of the point in time C at the reference surface in terms of a modified Carrington variable. C The time shift in Carrington coordinates will point to an earlier or the C same time. C Bad grid points in xcshift(i,j,n,itim,l) are filled with GridSphere2D. C The values of vRatio and dRatio are used to determine the changes C in the velocities and densities relative to their values at the source surface C at the point at the height i,j,n,itim,l. C MODIFICATION HISTORY: C APR-2002, B. Jackson (UCSD) C- subroutine shift_MOD(NLNGM,NLATM,NMAPM,NTM,V3,D3,V3L,D3L,R_MOD,IYR,DOY,IYRL,DOYL,ICTX,FICTX,ICTXL,FICTXL, & V3D,D3D,V3DL,D3DL,V3DTMP,D3DTMP, & nCar,JDCar,NCoff,xcbe,xctbeg,xctend,aLng,nLng,nLat,nMap,nT,rr,drr,Rconst, & cLrV,cLrD,vv,dd,xcshift,vRatio,dRatio) real V3(NLNGM,NLATM,NMAPM,3), ! Input 3D velocity from model (V west, north and out assumed +) & D3(NLNGM,NLATM,NMAPM), ! Input 3D density from model & V3L(NLNGM,NLATM,NMAPM,3), ! Save space for last 3D velocity from model & D3L(NLNGM,NLATM,NMAPM), ! Save space for last 3D density from model & R_MOD(NMAPM), ! Input heights from the Model C & XCbe(2,nT), ! Carrington time intervals & vv(nLng,nLat,nT), ! Tomography velocity map at height rr & dd(nLng,nLat,nT), ! Tomography density map at height rr & xcshift(nLng,nLat,nMap,nT,3), ! Map of accumulated tomographic model shifts at heights (from rr) & vRatio (nLng,nLat,nMap,nT), ! Map of accumulated tomographic velocity ratios at heights (from rr) & dRatio (nLng,nLat,nMap,nT), ! Map of accumulated tomographic density ratios at heights (from rr) C & XLT(3), ! Internal arrays & XLTM(4), & XLTMP(2),XLTMPL(2), C & V3D((NLNGM)*3+1,NLATM+1,NMAPM,3), ! These for NOT including 0, 360 -90 and +90 C & D3D((NLNGM)*3+1,NLATM+1,NMAPM), C & V3DL((NLNGM)*3+1,NLATM+1,NMAPM,3), C & D3DL((NLNGM)*3+1,NLATM+1,NMAPM), & V3D((NLNGM-1)*3+1,NLATM,NMAPM,3), ! These for including 0, 360 -90 and +90 & D3D((NLNGM-1)*3+1,NLATM,NMAPM), & V3DL((NLNGM-1)*3+1,NLATM,NMAPM,3), & D3DL((NLNGM-1)*3+1,NLATM,NMAPM), & V3Dtmp(nLng,nLat,2,2,3), & D3Dtmp(nLng,nLat,2,2) integer NDLT(3), ! Internal arrays & MNDLT(4), & NDLTMP(2), & NDLTM(2), & NDL(3) real*8 JDCar(nCar) external EARTH include 'MAPCOORDINATES.H' T1 = 0.0 if(IYRL .ne. 0) T1 = XMAP_SC_POS(EARTH,IYRL,DOYL,nCAR,JDCar) T2 = XMAP_SC_POS(EARTH,IYR,DOY,nCAR,JDCar) C print *, 'Into shift_MOD Last', T1, ICTXL, FICTXL, FICTXL + ICTXL - NCoff, ' Now', T2, ICTX, FICTX, FICTX + ICTX - NCoff do itim=2,nT T = XCtfull(itim) + XCbe(1,1) if(T .ge. T1 .and. T .lt. T2 .and. IYRL .ne. 0) then C print *, 'Into inner shift_MOD - T2, T1, itim, T, IYRL ', T2, T1, itim, T, IYRL LOINC = 1 ! MODEL longitude coordinates include 0 and 360 degrees C LOINC = 0 LAINC = 1 ! MODEL latitude coordinates include 0 and 180 degrees C LAINC = 0 DLATM = 180./(NLATM - LAINC) DLATM2 = (1-LAINC)*DLATM/2.0 BLAT = DLATM2 - 90.0 ! Beginning latitude value (near -90 degrees). XLTM(4) = 1. + (T - T1)/(T2 - T1) XLT(3) = XLTM(4) RoCon = Rconst*Sun__Spiral ! ~(1.49589E+08)/(60.*60.*24.*27.275) ! AU/CR?sec RCO = 1.0/Rconst ! ~27.275/25.38 C print *, 'Sun__Spiral =', rocon,RCO dnt = XCtfull(itim) - XCtfull(itim-1) ! Time change (at Earth ~27.275 days/rotation) in Carrington coordinates XTIM = dnt ! Time shift (in CR's at Earth) since the earlier rotation value dnth = dnt*rocon ! dnth in sec AU/km if R is in AU, dnt in CR at Earth Bad = BadR4() EPS = 0.00002 NLNGMT3 = (NLNGM-LOINC)*3 + 1 MNDLT(1) = nLng MNDLT(2) = nLat MNDLT(3) = 2 MNDLT(4) = 2 NDLTM(1) = NLNGM + (1 - LOINC) NDLTM(2) = NLATM + (1 - LAINC) NDLTMP(1) = NLNGMT3 NDLTMP(2) = NLATM + (1 - LAINC) NDLT(1) = nLng NDLT(2) = nLat NDLT(3) = nMap NDL(1) = nLng NDL(2) = nLat NDL(3) = 2 XCbeg = XCbe(1,itim) XCend = XCbe(2,itim) DLAT = 180./(nLat-1) XLON = aLng/(nLng-1) nLngLat = nLng*nLat C C Initialize the lowest level tomography values C do n=1,nT call ArrR4Zero (nLngLat, xcshift(1,1,1,n,1)) ! Shifts at rr are 0 call ArrR4Zero (nLngLat, xcshift(1,1,1,n,2)) ! Shifts at rr are 0 call ArrR4Zero (nLngLat, xcshift(1,1,1,n,3)) ! Shifts at rr are 0 call ArrR4Constant(nLngLat,1.,vRatio (1,1,1,n)) ! V ratios at rr are 1. call ArrR4Constant(nLngLat,1.,dRatio (1,1,1,n)) ! D ratios at rr are 1. end do C------- C Make V3D and D3D the same extended longitudinal length as the tomography program. C XCTXE = FICTX + ICTX - NCoff XCTXEL = FICTXL + ICTXL - NCoff nLng1 = NLNGMT3 nLat1 = NLATM - LAINC NLATMP = nLat1 + 1 NLNGM1 = NLNGM - 1 XCBeg = XCbe(1,itim) XCend = XCbe(2,itim) XCVE = XCvar(XCTXE + 1.) XCVEL = XCvar(XCTXEL + 1.) do i = 1, NLNGMT3 XCTX = i - XCVE + 1 if (XCTX .gt. NLNGM) XCTX = XCTX - NLNGM1 if (XCTX .gt. NLNGM) XCTX = XCTX - NLNGM1 if (XCTX .lt. 1.) XCTX = XCTX + NLNGM1 if (XCTX .lt. 1.) XCTX = XCTX + NLNGM1 if (XCTX .lt. 1.) XCTX = XCTX + NLNGM1 XLTMP(1) = XCTX XCTXL = i - XCVEL + 1 if (XCTXL .gt. NLNGM) XCTXL = XCTXL - NLNGM1 if (XCTXL .gt. NLNGM) XCTXL = XCTXL - NLNGM1 if (XCTXL .lt. 1.) XCTXL = XCTXL + NLNGM1 if (XCTXL .lt. 1.) XCTXL = XCTXL + NLNGM1 if (XCTXL .lt. 1.) XCTXL = XCTXL + NLNGM1 XLTMPL(1) = XCTXL C print *, 'i XCVE, XCVEL, XCTX, XCTXL', i, XCVE, XCVEL, XCTX, XCTXL do j=1, NLATMP XLA = (J-1)*DLATM + BLAT XLTMP(2) = max(1.,min(1.*NLATM,XLindx(XLA))) XLTMPL(2) = XLTMP(2) do n=2, NMAPM V3D(i,j,n,1) = FLINT(-2,NDLTM,V3(1,1,n,1),XLTMP,EPS) V3D(i,j,n,2) = FLINT(-2,NDLTM,V3(1,1,n,2),XLTMP,EPS) V3D(i,j,n,3) = FLINT(-2,NDLTM,V3(1,1,n,3),XLTMP,EPS) D3D(i,j,n) = FLINT(-2,NDLTM,D3(1,1,n), XLTMP,EPS) V3DL(i,j,n,1) = FLINT(-2,NDLTM,V3L(1,1,n,1),XLTMPL,EPS) V3DL(i,j,n,2) = FLINT(-2,NDLTM,V3L(1,1,n,2),XLTMPL,EPS) V3DL(i,j,n,3) = FLINT(-2,NDLTM,V3L(1,1,n,3),XLTMPL,EPS) D3DL(i,j,n) = FLINT(-2,NDLTM,D3L(1,1,n), XLTMPL,EPS) C if(n .eq. 1 .and. j .eq. 19) print *, 'shift_MOD base V3D(i,j,n,3) V3DL(i,j,n,3)',i,j,N,V3D(i,j,n,3),V3DL(i,j,n,3) C if(n.eq.15.and.j.eq.19) print *, 'D3D V3D',i,j,n,itim,XLTMPL(1),XLTMP(1),XLA,XLTMP(2),D3D(i,j,n),V3D(i,j,n,3) end do end do end do C C If there are bad model values, fill them (The bad values in the model data need to be previously filled with BadR4()) C do n=2,NMAPM call ArrR4getminmax(NLNGMT3*NLATMP,V3D(1,1,n,1),aminm,amaxm) if(aminm .eq. Bad) then print *, 'Bad values in height',n,' V3D 1' call GridSphere2D(aLng,NLNGMT3,NLATMP,1,V3D(1,1,n,1),cLrV,4,0.0,0.0)! Fill in bad values end if call ArrR4getminmax(NLNGMT3*NLATMP,V3D(1,1,n,2),aminm,amaxm) if(aminm .eq. Bad) then print *, 'Bad values in height',n,' V3D 2' call GridSphere2D(aLng,NLNGMT3,NLATMP,1,V3D(1,1,n,2),cLrV,4,0.0,0.0)! Fill in bad values end if call ArrR4getminmax(NLNGMT3*NLATMP,V3D(1,1,n,3),aminm,amaxm) if(aminm .eq. Bad) then print *, 'Bad values in height',n,' V3D 3' call GridSphere2D(aLng,NLNGMT3,NLATMP,1,V3D(1,1,n,3),cLrV,4,0.0,0.0)! Fill in bad values end if call ArrR4getminmax(NLNGMT3*NLATMP,D3D(1,1,n),aminm,amaxm) if(aminm .eq. Bad) then print *, 'Bad values in height',n,' D3D' call GridSphere2D(aLng,NLNGMT3,NLATMP,1,D3D(1,1,n), cLrD,4,0.0,0.0)! Fill in bad values end if call ArrR4getminmax(NLNGMT3*NLATMP,V3DL(1,1,n,1),aminm,amaxm) if(aminm .eq. Bad) then print *, 'Bad values in height',n,' V3DL 1' call GridSphere2D(aLng,NLNGMT3,NLATMP,1,V3DL(1,1,n,1),cLrV,4,0.0,0.0)! Fill in bad values end if call ArrR4getminmax(NLNGMT3*NLATMP,V3DL(1,1,n,2),aminm,amaxm) if(aminm .eq. Bad) then print *, 'Bad values in height',n,' V3DL 2' call GridSphere2D(aLng,NLNGMT3,NLATMP,1,V3DL(1,1,n,2),cLrV,4,0.0,0.0)! Fill in bad values end if call ArrR4getminmax(NLNGMT3*NLATMP,V3DL(1,1,n,3),aminm,amaxm) if(aminm .eq. Bad) then print *, 'Bad values in height',n,' V3DL 3' call GridSphere2D(aLng,NLNGMT3,NLATMP,1,V3DL(1,1,n,3),cLrV,4,0.0,0.0)! Fill in bad values end if call ArrR4getminmax(NLNGMT3*NLATMP,D3DL(1,1,n),aminm,amaxm) if(aminm .eq. Bad) then print *, 'Bad values in height',n,' D3DL' call GridSphere2D(aLng,NLNGMT3,NLATMP,1,D3DL(1,1,n), cLrD,4,0.0,0.0)! Fill in bad values end if end do C do n=1,NmapM C call ArrR4getminmax(NLNGMT3*NLATM,V3D(1,1,n,3),aminm,amaxm) C print *, 'V3D(1,1,n), n, aminm, amaxm', n,itim, aminm,amaxm C call ArrR4getminmax(NLNGMT3*NLATM,D3D(1,1,n),aminm,amaxm) C print *, 'D3D(1,1,n), n, aminm, amaxm',n, itim, aminm,amaxm C call ArrR4getminmax(NLNGMT3*NLATM,V3DL(1,1,n,3),aminm,amaxm) C print *, 'V3DL(1,1,n), n, aminm, amaxm', n,itim, aminm,amaxm C call ArrR4getminmax(NLNGMT3*NLATM,D3DL(1,1,n),aminm,amaxm) C print *, 'D3DL(1,1,n), n, aminm, amaxm',n, aminm,amaxm C end do C nLat1 = NLATM - LAINC C do i=1,nLng C XLO = (1-i)*XLON + Alng + XCbeg C XLTMP(1) = max(1.,min(1.*NLNGMT3,XCvar(XLO))) C do j=1,nLat C XLA = (J - 1)*DLAT - 90. C XLTMP(2) = XLindx(XLA) C V3DtmpT(i,j,2) = FLINT(2,NDLTMP,V3D(1,1,1,3),XLTMP,EPS) C D3DtmpT(i,j,2) = FLINT(2,NDLTMP,D3D(1,1,1),XLTMP,EPS) C V3DtmpT(i,j,1) = FLINT(2,NDLTMP,V3DL(1,1,1,3),XLTMP,EPS) C D3DtmpT(i,j,1) = FLINT(2,NDLTMP,D3DL(1,1,1),XLTMP,EPS) C end do C end do C nLng1 = nLng-1 C nLat1 = nLat-1 C do i = 1, nLng ! Interpolate in time, longitude and latitude for the first height for check. C XLT(1) = XCvar(XCfull(i)) C do j=1, nLat C XLT(2) = XLindx(XLdeg(j)) C vv(i,j,itim) = FLINT(3,NDL,V3DtmpT,XLT,EPS) C dd(i,j,itim) = FLINT(3,NDL,D3DtmpT,XLT,EPS) C if(j .eq. 10) print *, 'shift_MOD base vv(i,j,itim)',i,j,itim,vv(i,j,itim) C if(itim .eq. 2) then C vv(i,j,1) = vv(i,j,2) C dd(i,j,1) = dd(i,j,2) C end if C end do C end do C------- C print *, 'To the do loops' do k=2, NMAPM do n=1,nMap Rn = RRhght(n) if(R_MOD(k) . ge. Rn .and. R_MOD(k-1) . lt. Rn) then C print *, 'Into the height-time interpolation', k, n, R_MOD(k), Rn, itim XLTM(3) = 1. + (Rn - R_MOD(k-1))/(R_MOD(k) - R_MOD(k-1)) R2 = (RRhght(1)/Rn)**2 nLng1 = NLNGMT3 - 1 nLat1 = NLATM - LAINC do j=1,nLat XLA = (J - 1)*DLAT - 90. XLTMP(2) = XLindx(XLA) do i=1,nLng XLO = (1-i)*XLON + ALng + XCbeg XLTMP(1) = max(1.,min(1.*NLNGMT3,XCvar(XLO)) ) V3Dtmp(i,j,1,1,1) = FLINT(2,NDLTMP,V3DL(1,1,k-1,1),XLTMP,EPS) V3Dtmp(i,j,2,1,1) = FLINT(2,NDLTMP,V3DL(1,1,k,1),XLTMP,EPS) V3Dtmp(i,j,1,1,2) = FLINT(2,NDLTMP,V3DL(1,1,k-1,2),XLTMP,EPS) V3Dtmp(i,j,2,1,2) = FLINT(2,NDLTMP,V3DL(1,1,k,2),XLTMP,EPS) V3Dtmp(i,j,1,1,3) = FLINT(2,NDLTMP,V3DL(1,1,k-1,3),XLTMP,EPS) V3Dtmp(i,j,2,1,3) = FLINT(2,NDLTMP,V3DL(1,1,k,3),XLTMP,EPS) V3Dtmp(i,j,1,2,1) = FLINT(2,NDLTMP,V3D(1,1,k-1,1),XLTMP,EPS) V3Dtmp(i,j,2,2,1) = FLINT(2,NDLTMP,V3D(1,1,k,1),XLTMP,EPS) V3Dtmp(i,j,1,2,2) = FLINT(2,NDLTMP,V3D(1,1,k-1,2),XLTMP,EPS) V3Dtmp(i,j,2,2,2) = FLINT(2,NDLTMP,V3D(1,1,k,2),XLTMP,EPS) V3Dtmp(i,j,1,2,3) = FLINT(2,NDLTMP,V3D(1,1,k-1,3),XLTMP,EPS) V3Dtmp(i,j,2,2,3) = FLINT(2,NDLTMP,V3D(1,1,k,3),XLTMP,EPS) D3Dtmp(i,j,1,1) = FLINT(2,NDLTMP,D3DL(1,1,k-1),XLTMP,EPS) D3Dtmp(i,j,2,1) = FLINT(2,NDLTMP,D3DL(1,1,k),XLTMP,EPS) D3Dtmp(i,j,1,2) = FLINT(2,NDLTMP,D3D(1,1,k-1),XLTMP,EPS) D3Dtmp(i,j,2,2) = FLINT(2,NDLTMP,D3D(1,1,k),XLTMP,EPS) C if(n.eq.15.and.j.eq.10) print *, 'Tmp',i,j,n,itim,XLO,XLTMP(1),XLA,XLTMP(2),D3Dtmp(i,j,1,2),V3Dtmp(i,j,1,2,3) end do end do C call ArrR4getminmax(nLngLat*4,V3Dtmp(1,1,1,1,1),aminm,amaxm) C print *, 'aminm,amaxm V3Dtmp(i,j,1-2,1-2,1)', aminm,amaxm C call ArrR4getminmax(nLngLat*4,V3Dtmp(1,1,1,1,2),aminm,amaxm) C print *, 'aminm,amaxm V3Dtmp(i,j,1-2,1-2,2)', aminm,amaxm C call ArrR4getminmax(nLngLat*4,V3Dtmp(1,1,1,1,3),aminm,amaxm) C print *, 'aminm,amaxm V3Dtmp(i,j,1-2,1-2,3), n, Rn', aminm,amaxm,n,Rn,R_MOD(k-1),R_MOD(k) C call ArrR4getminmax(nLngLat*4,D3Dtmp(1,1,1,1),aminm,amaxm) C print *, 'aminm,amaxm D3Dtmp(i,j,1-2,1-2)', aminm,amaxm nLng1 = nLng - 1 nLat1 = nLat - 1 do j=1, nLat aLatj = XLdeg(j) XLTM(2) = XLindx(XLdeg(j)) do i=1, nLng if( itim .eq. 2) then ! The second tomography time must prepare the first tomo. shifts XLTM4SAVE = XLTM(4) XLTM(4) = 1.0 XLTM(1) = max(1.,min(1.*nLng,XCvar(XCfull(i) - dnt*RCO)) ) V11 = FLINT(4,MNDLT,V3Dtmp(1,1,1,1,1),XLTM,EPS) V22 = FLINT(4,MNDLT,V3Dtmp(1,1,1,1,2),XLTM,EPS) V33 = FLINT(4,MNDLT,V3Dtmp(1,1,1,1,3),XLTM,EPS) D33 = FLINT(4,MNDLT,D3Dtmp(1,1,1,1), XLTM,EPS) if(V33.eq.0.) stop 'Zero preliminary V33 velocity in shift_MOD' dnthh = RoCon*(RRhght(n) - RRhght(1))/V33 ! total time (rot.) to travel from RRhght(1) to RRhght(n) ALOCH = atan2d(dnthh*V11,Rn)/360. ! change of longitude of point from velocity (+ to west) ALACH = atan2d(dnthh*V22,Rn) ! change of latitude of point from velocity (+ to north) xcshift(i,j,n,itim-1,1) = ALOCH ! first longitude shift xcshift(i,j,n,itim-1,2) = -ALACH ! first latitude shift xcshift(i,j,n,itim-1,3) = -dnthh ! first time shift XC = xcshift(i,j,n,itim-1,1) + XCfull(i) XC = max(XCbeg,min(XCend,XC)) XL = xcshift(i,j,n,itim-1,2) + ALatj XL = max(-90.,min(90.,XL)) XCT = XCtfull(xctbeg) XCT = max(XCtbeg,min(XCtend,XCT)) call get3DTval(xcbe,xctbeg,xctend,nLng,nLat,nT,dd,1, & XC,XL,XCT,Den) if(D33 .ne. Bad) then dRatio(i,j,n,itim-1) = D33/(Den*R2) else dRatio(i,j,n,itim-1) = Bad end if call get3DTval(xcbe,xctbeg,xctend,nLng,nLat,nT,vv,1, & XC,XL,XCT,Vel) if(V33 .ne. Bad) then vRatio(i,j,n,itim-1) = V33/Vel else vRatio(i,j,n,itim-1) = Bad end if C if(n.eq.15.and.j.eq.10) print *, 'PRE',i,j,n,itim,XC,XL,XCT,dRatio(i,j,n,itim-1),vRatio(i,j,n,itim-1) C if(n.eq.15.and.j.eq.10) print *, ALOCH,ALatj,ALACH,Rn,D33,Den*R2,V33,Vel C if(n.eq.15.and.j.eq.10) print *, i,j,n,itim,dnthh,xcshift(i,j,n,itim-1,1),xcshift(i,j,n,itim-1,2),xcshift(i,j,n,itim-1,3) XLTM(4) = XLTM4SAVE ! Restore the interpolate value for itim = 2 end if XLTM(1) = XCvar(XCfull(i)) C C interpolate in height and time C V11 = FLINT(4,MNDLT,V3Dtmp(1,1,1,1,1),XLTM,EPS) V22 = FLINT(4,MNDLT,V3Dtmp(1,1,1,1,2),XLTM,EPS) V33 = FLINT(4,MNDLT,V3Dtmp(1,1,1,1,3),XLTM,EPS) D33 = FLINT(4,MNDLT,D3Dtmp(1,1,1,1), XLTM,EPS) if(V33.eq.0.) stop 'Zero V33 velocity in shift_MOD' dnthh = RoCon*drr/V33 ! total time (rot.) to travel from RRhght(n-1) to RRhght(n) RCH = V33*dnt/RoCon ! change of height particle travels during the tomo. time interval if(RCH .ge. drr) RCH = drr ! don't let the particle start lower than the earlier lower level if(RCH .lt. drr) dnthh = dnt ! if the particle hasn't gotten that far, time is the tomo. one ALOCH = atan2d(dnthh*V11,Rn)/360. ! change of longitude of point from velocity (+ to west) XLT(1) = XCvar(XCfull(i) + ALOCH) ALACH = atan2d(dnthh*V22,Rn) ! change of latitude of point from velocity (+ to north) XLT(2) = XLindx(ALatj - ALACH) XLT(3) = RRindx(Rn - RCH) xcshift(i,j,n,itim,1) = FLINT(-3,NDLT,xcshift(1,1,1,itim-1,1),XLT,EPS) + ALOCH ! longitude shift xcshift(i,j,n,itim,2) = FLINT(-3,NDLT,xcshift(1,1,1,itim-1,2),XLT,EPS) + ALACH ! latitude shift xcshift(i,j,n,itim,3) = FLINT(-3,NDLT,xcshift(1,1,1,itim-1,3),XLT,EPS) - dnthh ! time shift XC = xcshift(i,j,n,itim,1) + XCfull(i) XC = max(XCbeg,min(XCend,XC)) XL = xcshift(i,j,n,itim,2) + ALatj XL = max(-90.,min(90.,XL)) XCT = xcshift(i,j,n,itim,3) + XCtfull(itim) XCT = max(XCtbeg,min(XCtend,XCT)) call get3DTval(xcbe,xctbeg,xctend,nLng,nLat,nT,dd,1, & XC,XL,XCT,Den) if(D33 .ne. Bad) then dRatio(i,j,n,itim) = D33/(Den*R2) else dRatio(i,j,n,itim) = Bad end if call get3DTval(xcbe,xctbeg,xctend,nLng,nLat,nT,vv,1, & XC,XL,XCT,Vel) if(V33 .ne. Bad) then vRatio(i,j,n,itim) = V33/Vel else vRatio(i,j,n,itim) = Bad end if C if(n.eq.15.and.j.eq.10) print *, '1',i,j,n,itim,XC,XL,XCT,dRatio(i,j,n,itim),vRatio(i,j,n,itim) C if(n.eq.15.and.j.eq.10) print *, '1',ALOCH,ALatj,ALACH,Rn,D33,Den*R2,V33,Vel C if(n.eq.15.and.j.eq.10) print *, '1',i,j,n,itim,RCH,XLT(1),XLT(2),XLT(3),xcshift(i,j,n,itim,1),xcshift(i,j,n,itim,2),xcshift(i,j,n,itim,3) end do end do end if end do end do NN = 1 if(itim .eq. nT) NN = 2 do it=1,NN ! Handle the final tomography time differently to fill it with shifts, too. itimm = itim + it - 1 do n=2,nMap ! Now after the shifts are all determined in the inertial frame, add longitude shift do j=1, nLat ! to the last tomography time. The last tomography time is used above. do i=1, nLng xcshift(i,j,n,itimm-1,1) = xcshift(i,j,n,itimm-1,1) + xcshift(i,j,n,itimm-1,3)*RCO ! True longitude shift C if(n.eq.15.and.j.eq.10) print *, '2',i,j,n,itim,xcshift(i,j,n,itimm-1,1),xcshift(i,j,n,itimm-1,2),xcshift(i,j,n,itimm-1,3) end do end do end do end do end if end do do i=1, NLNGM do j=1, NLATM do n=1, NMAPM V3L(i,j,n,1) = V3(i,j,n,1) V3L(i,j,n,2) = V3(i,j,n,2) V3L(i,j,n,3) = V3(i,j,n,3) D3L(i,j,n) = D3(i,j,n) end do end do end do IYRL = IYR DOYL = DOY XCTXEL = XCTXE FICTXL = FICTX ICTXL = ICTX return end