C+ C NAME: C GridSphere3D C PURPOSE: C Fills 'holes' and/or smooths a 3D array using a Gaussian weighting function. C Two of the array dimensions represent a regular grid in spherical C coordinates (longitude and latitude). The third represents a regular C grid in radial distance. C CALLING SEQUENCE: subroutine GridSphere3D(dXC,nLng,nLat,nR,Z,WidDeg,R1,iBadZ,WThreshold,ClipLng) C INPUTS: C dXC real range of Carrington variable covered by nLng C i.e. number of 360 deg rotations. C nLng integer # longitudes C nLat integer # latitudes C nR integer # radial distances C Z(nLng,nLat,nR) real 3D array of function values C WidDeg real Width of Gaussian used for angular C smoothing (in degrees) (WidDeg<0 has C special meaning; see PROCEDURE) C R1 real (3D version only) C R1 > 0:Radial distance for inner boundary (in C grid spacings), i.e. location of Z(*,*,1) C R1=0 IS INVALID; MUST BE R1>0 C R1 < 0:(value does not matter) C suppresses the radial smoothing, i.e. C each subarray Z(*,*,i),i=1,nR is C 2D smoothed over longitude and latitude only. C iBadZ integer 0: Invalid elements remain invalid; C valid elements are replaced by a C smoothed value C 1: All elements (valid and invalid) are C replaced by smoothed values. C 2: Invalid elements are replaced by C smoothed values; valid elements remain C untouched. C C The following two options are useful if all C bad bins have to be filled in with something. C Use them at your own risk (see PROCEDURE). C C 3: Same as iBadZ=1, but if any invalid C elements remain at the end, these are C all filled in by a call to GridFill C 4: Same as iBadZ=2, but if any invalid C elements remain at the end, these are C all filled in by a call to GridFill C C >=10: Use open grid for longitude and latitude (?) C >=100: Use lookup table for Gaussian C C WThreshold real The replacement by a smoothed value is C made only if total weight is C larger than/equal to the threshold value. C (see PROCEDURE). C ClipLng real Longitude difference (see PROCEDURE) C (usually set to zero). C OUTPUTS: C Z(nLng,nLat,nR) real smoothed array of function values C INCLUDE: include 'math.h' C CALLS: C iArrR4ValuePresent, Say, ArrR4Copy, GridFill, GaussLookup C RESTRICTIONS: C > GridSphere3D uses internal scratch arrays, which are C dimensioned using 3 parameters, LNG,LAT and LNR. If the scratch C arrays are too small, program execution stops. C > iBadZ=3 and iBadZ=4: C GridFill uses a 2D algorithm to fill in holes. If nR>1 then GridFill C is applied to each level k=1,nR separately. C PROCEDURE: C > Invalid or 'Bad' array elements are indicated by the value BadR4() C > The 1st dimension represents a regular grid in longitude covering C [0,360] degrees, i.e. Lng(i) = 360*(i-1)/(nLng-1), i=1,nLng. C The 2nd dimension represents a regular grid in latitude covering C [-90,+90] degrees, i.e. Lat(j) = -90*[-1+2*(j-1)/(nLat-1)], j=1,nLat. C The 3rd dimension represents a regular grid in radial distance, C i.e. R(k) = R1+k-1 C > Obviously, 2D arrays can be smoothed by setting nR=1. C > GridSphere2D is just a special call to GridSphere3D with R1=-1 (the C value does not matter as long as it's negative). C > The weighting is Gaussian with weight=1 at the origin. To avoid loosing C valid elements for iBadZ=0 or 1, set the threshold weight to a C value less than 1. C > The angular weighting factor involves the angular distance between C two locations on a sphere. The points are connected by two arcs, C one is the minimum angular distance A (0 180+3*WidDeg is probably not a bad place to start. C > If WidDeg > 0, the smoothing procedure is speeded up by looping C over only part of the array when calculating the smoothed function C value. Excluded are all elements which are more than 3 times the C Gaussian width away from the element for which the smoothed value is C calculated (the weight would be less than exp(-9)). The exclusion C is based on an approximate calculation of the distance between points. C Near the poles there may be a problem (see the explanation in the code), C but I think it's under control now. C In case you don't trust the result set WidDeg to a negative value. C If WidDeg <0, then the absolute value is used as Gaussian width, and C the whole array is used in the calculation of the smoothed average. C The procedure will be a LOT slower though. C > If ClipLng=0 then the averaging algorithm is mod(XCend-XCbeg), i.e. if C in the input array column 1 (XC = XCend) and column nLng (XC = XCbeg) C are the same this is also true for the output array. C HOWEVER this symmetry is broken if iBadZ=3 or 4 because of the extra C call to GridFill. C MODIFICATION HISTORY: C FEB-1997, Paul Hick (UCSD/CASS; pphick@ucsd.edu) C- real dXC integer nLng integer nLat integer nR real Z(nLng,nLat,nR) real WidDeg real R1 integer iBadZ real WThreshold real ClipLng character cSay*12 /'GridSphere3D'/ parameter (LNG =1800) ! 3*108+1 Dimensions scratch arrays parameter (LAT = 900) parameter (LNR = 1) parameter (WGauss = 3.0) real cosLng(0:LNG-1) ! Scratch arrays real cosLat(LAT) real sinLat(LAT) real Ztmp (LNG*LAT*LNR) integer*1 NSIDE(LNG*LAT*9) equivalence (Ztmp,NSIDE) logical bLimitScan logical bInclPole logical bTestPole logical bFillAll logical bOpenGrid logical bLookup logical bRange integer nLngOld /0/ save nLngOld integer nLatOld /0/ save nLatOld real dLng, dLat, rLat0 save dLng, dLat, cosLng, sinLat, cosLat, rLat0 real range_lo /0.0/ save range_lo real range_hi /0.0/ save range_hi logical bWeightsOn save bWeightsOn real Zweight(LNG*LAT*LNR) save Zweight real ZThreshold save ZThreshold fncLng(i) = i*dLng ! Func for longitude difference (radians) fncLat(j) = rLat0+j*dLat ! Func for latitude (radians) fncRad(k) = R0+k ! Func for radial distance (grid spacings) indZ(i,j,k) = ((k-1)*nLat+j-1)*nLng+i ! 1-Dim index for addressing Ztmp rBad = BadR4() iFillBadZ = iBadZ ! Pick up input value range = range_hi-range_lo bRange = range .ne. 0.0 bLookup = iFillBadZ .ge. 100 iFillBadZ = mod(iFillBadZ,100) bOpenGrid = iFillBadZ .ge. 10 iFillBadZ = mod(iFillBadZ,10) bFillAll = iFillBadZ .ge. 3 ! Test for 3,4 if (bFillAll) iFillBadZ = iFillBadZ-2 ! Change 3,4 to 1,2 j = nLng*nLat*nR i = iArrR4ValuePresent(j,Z,rBad) if (iFillBadZ .eq. 2 .and. i .eq. 0) return ! No bad values to fill in if (i .eq. j) then ! Return without doing anything call Say(cSay,'W','Invalid','no valid array elements') return end if if (LNG .lt. nLng .or. LAT .lt. nLat .or. LNG*LAT*LNR .lt. j) & call Say(cSay,'E','Parameter','too small: LAT, LNG or LNR') !------- ! A new set of cosines and sines are calculated if the input values of nLng ! or nLat are changed. The arrays are explicitly included in a SAVE statement. if (nLng .ne. nLngOld) then call Say(cSay,'I','Refreshing','cosLng array') if (bOpenGrid) then dLng = 2.0*MATH__PI/nLng*dXC ! Longitudinal grid spacing (radians) else dLng = 2.0*MATH__PI/(nLng-1)*dXC! Longitudinal grid spacing (radians) end if do i=0,nLng-1 cosLng(i) = cos(fncLng(i)) end do nLngOld = nLng end if if (nLat .ne. nLatOld) then call Say(cSay,'I','Refreshing','cosLat and sinLat array') if (bOpenGrid) then dLat = MATH__PI/nLat ! Latitudinal grid spacing (radians) rLat0 = -0.5*MATH__PI+0.5*dLat-dLat else dLat = MATH__PI/(nLat-1) ! Latitudinal grid spacing (radians) rLat0 = -0.5*MATH__PI-dLat end if do j=1,nLat Tmp = fncLat(j) sinLat(j) = sin(Tmp) cosLat(j) = max(0.0,cos(Tmp)) ! NOT redundant !!! end do nLatOld = nLat end if WidRad = MATH__RPD*abs(WidDeg) ! Gaussian width in radians bLimitScan = WidDeg .gt. 0.0 ! Activates fast version R0 = abs(R1)-1.0 ! Used by statement func fncRad i3DSmooth = 1 ! Full 3D smoothing if (R1 .lt. 0.0) i3DSmooth = 0 ! 2D smoothing only iClip = ClipLng*MATH__RPD/dLng ! Limit on longitude difference if (iClip .le. 0) iClip = nLng-1 !------- ! Smoothed values are calculated for each grid point (i,j,k) in the map. This is ! implemented by the outer loops (loop indices k,j,i). !$omp parallel private(i,j,k,Rk,dRad,kLimit,k1lo,k1hi,sLat,cLat, !$omp& jLimit,j1lo,j1hi,bInclPole,jSouth,jNorth,iLimit,i1lo,i1hi,ii, !$omp& bTestPole,W,F,k1,j1,i1,dR,sLat1,cLat1,i1lo1,i1hi1,i1s,dA,dW,dZ) do k=1,nR Rk = fncRad(k) ! Radial distance (grid spacings) dRad = Rk*WidRad ! Radial width of Gaussian (grid spacings) !------- ! Limit range of k1-values to be included to several times the Gaussian width ! Note 1: if bLimitScan=.FALSE. then [k1lo,k1hi]=[1,nR]. ! Note 2: if i3DSmooth=0 (2D smoothing only) then k1lo=k1hi=k, i.e. only ! elements with k1=k will be used in the smoothing. kLimit = nR-1 if (bLimitScan) kLimit = WGauss*dRad kLimit = i3DSmooth*kLimit k1lo = max( 1,k-kLimit) ! Keep inside [1,nR] k1hi = min(nR,k+kLimit) Rk = Rk/dRad ! = 1./WidRad (independent of k!!) ! This is why R1 must be >0 (divide by zero) !$omp do schedule(static) do j=1,nLat sLat = sinLat(j) cLat = cosLat(j) !------- ! Limit range of j1-values to be included to several times the ! Gaussian width. Note that this limit on j1 is strict: all ! elements with j1 outside [j1lo,j1hi] do have an angular distance ! bigger than the imposed limit. Note that if bLimitScan=.FALSE. ! then [j1lo,j1hi]=[1,nLat]. jLimit = nLat-1 if (bLimitScan) jLimit = WGauss*WidRad/dLat j1lo = max( 1,j-jLimit) ! Keep inside [1,nLat] j1hi = min(nLat,j+jLimit) bInclPole = bLimitScan .and. (j1lo .eq. 1 .or. j1hi .eq. nLat) if (bOpenGrid) then jSouth = j-1 jNorth = nLat+nLat+1-j else jSouth = j-2 jNorth = nLat+nLat-j end if do i=1,nLng if (iFillBadZ .eq. 0 .and. Z(i,j,k) .eq. rBad) then !-------- ! iFillBadZ=0: DO NOT fill bad values; only smooth the valid ! function values. ! Do not calculate new function values for bad elements. Ztmp(indZ(i,j,k)) = rBad ! Keep bad elements else if (iFillBadZ .eq. 2 .and. Z(i,j,k) .ne. rBad .and. & (.not. bWeightsOn .or. Zweight(indZ(i,j,k)) .ge. ZThreshold)) then !------- ! iFillBadZ=2: NO smoothing; only fill in bad values) ! If the function value is valid, carry the value unmodified into Ztmp. Ztmp(indZ(i,j,k)) = Z(i,j,k) ! Don't change good elements else ! Calculate a smoothed value !------- ! Limit range of i1-values to be included to several times the ! Gaussian width. Note that this limit on i1 is NOT strict: NOT all ! elements with i1 outside [i1lo,i1hi] necessarily have an angular ! distance bigger than the imposed limit. iLimit is measured along a ! small circle at constant latitude through the point (i,j,k). Near ! the poles this is not a good measure for the angular distance between ! points (i,j,k) and (i1,j1,k1). In particular, if (i,j,k) is close the ! pole, then points (i1,j1,k1) close to the same pole but near the ! opposite longitude may be inappropriately discarded (I think). This ! is corrected for once we know more about the points (i1,j1,k1) (in ! particular the value of j1 is important). bTestPole is used later to ! decide whether a correction is necessary. Note that if bLimitScan= ! .FALSE. the range [i1lo,i1hi] is a range of nLng ! elements centered on element i, i.e. all longitudes are included. iLimit = nLng/2 if (bLimitScan .and. cosLat(j) .ne. 0.) & iLimit = min(iLimit,int(WGauss*WidRad/cosLat(j)/dLng)) i1lo = i-iLimit ! Could be outside [1,nLng] !! i1hi = i+iLimit if (nLng/2*2 .eq. nLng) i1lo = i1lo+1 ! Don't remove !! bTestPole = bInclPole .and. i1hi-i1lo+1 .ne. nLng !-------- ! For each grid point (i,j) the smoothed value is calculated ! by averaging over the entire array (if WidDeg .lt. 0) or over ! a limited part of the map near point (i,j) (if WidDeg .gt. 0). ! This is implemented by the inner loops (loop indices j1,i1,k1). ! ! There are two ways to do the weighting: ! ! 1. Use the real distance between points P(i,j,k) and P1(i1,j1,k1) to ! calculate the Gaussian weight factor: ! Use the cosine rule for a spherical triangle to get the angular distance ! A = angle(PSP1), between P and P1 (S=Sun center): ! cos(A) = sLat*sLat1+cLat*cLat1*cosLng(ii) ! Then use the cosine rule in plane triagle SPP1 to get the distance, d, ! between P1 and P: ! d^2 = Rk*Rk+Rk1*Rk1-2*Rk*Rk1*cos(A) ! With Rk and Rk1 expressed in units of dRad (Gaussian width in units ! of the grid spacing, the weight factor becomes dW = exp(-d^2). ! This is the most elegant way to do the smoothing for a 3D array, but ! has the disadvantage that it does not reduce to an averaging over angular ! coordinates if nR=1. Since I want to use the same procedure for the ! smoothing of a 2D map in angular coordinates, the following compromise ! is used: ! ! 2. Calculate the angular distance, A, between P and P1 as above. ! Convert A to a linear distance, dA=Rk*A. The radial distance between P ! and P1 is dR=Rk-Rk1 (Rk,Rk1 in units of dRad). The weight factor is ! calculated by combining these two independent distance measures: ! dW = exp(-dA^2-dR^2). Since dA^2+dR^2 is NOT the true distance between ! P and P1, this is NOT a true spatial Gaussian smoothing procedure. ! This way if a 2D array is averaged (nR=1) then dR=0, and the above dW ! reduced to a weight factor dW = exp(-dA^2)=exp(-[A/WidRad]^2. W = 0.0 ! Initialize sum of weights F = 0.0 ! Initialize sum of weighted function values do k1=k1lo,k1hi dR = Rk-fncRad(k1)/dRad do j1=j1lo,j1hi sLat1 = sLat*sinLat(j1) cLat1 = cLat*cosLat(j1) !------- ! The following is kinda kludgy, and probably unnecessary. ! If j1lo=1 or j1hi=nLat, then there may be points (i1,j1) across the pole ! from (i,j), which are excluded from the range [i1lo,i1hi], while still ! being close to (i,j). The distance between these points across the pole ! will be j-1+j1-1=jSouth+j1 across the south pole and nLat-j+nLat-j1=jNorth-j1 ! across north pole. If this distance is smaller than jLimit, then make sure ! to include them by including all longitudes [1,nLng] at latitude j1. i1lo1 = i1lo i1hi1 = i1hi if (bTestPole .and. min(jSouth+j1,jNorth-j1) .lt. jLimit) then i1lo1 = 1 i1hi1 = nLng end if !------- ! Note that i1lo1 and i1hi1 could be outside the range [1,nLng] do i1s=i1lo1,i1hi1 ! Don't touch the loop variable i1 = 1+mod(i1s-1+nLng,nLng) ! Put back in range [1,nLng] ii = abs(i1-i) !------- ! The restriction ii<=iClip is used to avoid inclusion of points ! almost 360 degree away from ii if (Z(i1,j1,k1) .ne. rBad .and. ii .le. iClip) then dA = sLat1+cLat1*cosLng(ii) dA = max(-1.0,min(1.0,dA))! NOT redundant !! dA = Rk*acos(dA) ! Rk=1./WidRad if (bLookup) then dW = GaussLookup(dA,dR) else dW = exp(-dA*dA-dR*dR)! Weight factor end if dZ = Z(i1,j1,k1) if (bRange .and. W .ne. 0.0) then dA = dZ-F/W if (abs(dA) .gt. 0.5*range) dZ = dZ-sign(range,dA) end if W = W+dW ! Sum weight factors F = F+dW*dZ ! Sum weighted function values end if end do end do end do if (W .gt. WThreshold) then ! Check threshold dZ = F/W if (bRange) then if (dZ .le. range_lo) then dZ = dZ+range else if (dZ .gt. range_hi) then dZ = dZ-range end if end if Ztmp(indZ(i,j,k)) = dZ else Ztmp(indZ(i,j,k)) = rBad end if end if end do end do !$omp end do nowait end do !$omp end parallel call ArrR4Copy(nLng*nLat*nR,Ztmp,Z) ! Copy into output array if (bFillAll) then do k=1,nR call GridFill(0,nLng,nLat,Z(1,1,k),NSIDE,Zmin,Zmax) end do end if range_lo = 0.0 range_hi = 0.0 bWeightsOn = .FALSE. return C+ C NAME: C GridSphereRange C PURPOSE: C CALLING SEQUENCE: entry GridSphereRange(range_lo_,range_hi_) C MODIFICATION HISTORY: C MAY-2007, Paul Hick (pphick@ucsd.edu) C- range_lo = range_lo_ range_hi = range_hi_ return C+ C NAME: C GridSphereWeight C PURPOSE: C CALLING SEQUENCE: entry GridSphereWeight(nLng,nLat,nR,Z,ZThreshold_) C MODIFICATION HISTORY: C MAY-2007, Paul Hick (pphick@ucsd.edu) C- bWeightsOn = .TRUE. !$omp parallel private(i,j,k) do k=1,nR !$omp do schedule(static) do j=1,nLat do i=1,nLng Zweight(indZ(i,j,k)) = Z(i,j,k) end do end do !$omp end do nowait end do !$omp end parallel ZThreshold = ZThreshold_ return end