MODULE DWeighting ! Initial letter "D" indicates a DOUBLE PRECISION (REAL*8) version, ! created 2015.07 as part of a systematic upgrade of several of my codes ! (Shells, NeoKinema, FiniteMap, & NeoKineMap) to 64-bit precision. ! It is intended that the "D" FUNCTIONs and SUBROUTINEs here should be ! logically equivalent to those in MODULE Weighting, ! just more precise. Naturally, they now all take REAL*8 arguments ! in place of the old REAL arguments. ! Computes "weights" (with mean value 1.0) ! for a list of "number_of_data" datum positions, ! each of which is specified by "theta_radians" ! (colatitude, measured from North pole) and ! "phi_radians" (longitude, measured East from Greenwich). ! Exploits pre-existing code (USE DIcosahedron) ! to create a global finite-element grid of equal-sized ! spherical triangles. ! Then, each datum position is assigned to a finite element, ! and counts of total data in each finite element are kept. ! The surface area (actually, solid angle, in steradians) ! of the element is equally divided among these data. ! Finally, surface areas of empty elements are added to ! the area associated with the nearest datum. ! The result is that weights are proportional to surface areas !"associated with" each datum, where "association" is procedurally ! defined by this algorithm. ! by Peter Bird, UCLA, November 2006. ! (c) Copyright 2006 by G. Peter Bird and the ! Regents of the University of California; ! REAL*8 revision copyright 2015. USE DSphere ! Peter Bird's Fortran 90 MODULE file DSphere.f90 PRIVATE ! <=== Very Important! We only desire access through ! CALL DInitialize_Weighting(subdivision) ! and ! CALL DPerform_Weighting(...) ! ; ! we don't want to make all internal variables and arrays visible! PUBLIC DInitialize_Weighting PUBLIC DIs_Initialization_Needed PUBLIC DPerform_Weighting ! The following values need to have a "lifetime" that extends ! beyond the end of Initialize_Weighting, ! so they are available for later calls to Perform_Weighting. SAVE INTEGER :: n_slice, numNod, numEl INTEGER, DIMENSION(:,:), ALLOCATABLE :: nodes INTEGER, DIMENSION(:,:), ALLOCATABLE :: neighbor ! neighbors of each element LOGICAL :: initialized = .FALSE. REAL*8, DIMENSION(:,:), ALLOCATABLE :: center ! center uvecs of elements REAL*8, DIMENSION(:), ALLOCATABLE :: a_ ! (plane) areas of elements REAL*8, DIMENSION(:,:), ALLOCATABLE :: node_uvec CONTAINS SUBROUTINE DInitialize_Weighting (subdivision) !Create global finite-element grid and store it. !Suggested values of subdivision are 2 ... 5. USE DIcosahedron ! Peter Bird's Fortran 90 MODULE file DIcosahedron.f90 IMPLICIT NONE SAVE LOGICAL :: chatty INTEGER, INTENT(IN) :: subdivision !Note that its value is captured and kept here, as n_slice: n_slice = subdivision numNod = 2 + 10 * (4**subdivision) IF (ALLOCATED(node_uvec)) DEALLOCATE(node_uvec) ALLOCATE ( node_uvec(3, numNod) ) numEl = 20 * (4**subdivision) IF (ALLOCATED(nodes)) DEALLOCATE(nodes) ALLOCATE ( nodes(3, numEl) ) CALL DMake_Global_Grid (n_slice = subdivision, & ! only input(!) & numNod = numNod, node_uvec = node_uvec, & ! output: number of nodes, unit vectors of nodes, & numEl = numEl, nodes = nodes) ! number of elements, element definitions IF (ALLOCATED(a_)) DEALLOCATE(a_) ALLOCATE ( a_(numEl) ) IF (ALLOCATED(center)) DEALLOCATE(center) ALLOCATE ( center(3, numEl) ) IF (ALLOCATED(neighbor)) DEALLOCATE(neighbor) ALLOCATE ( neighbor(3, numEl) ) chatty = .FALSE. CALL DLearn_Spherical_Triangles (numEl, nodes, node_uvec, chatty, & & a_, center, neighbor) initialized = .TRUE. END SUBROUTINE DInitialize_Weighting LOGICAL FUNCTION DIs_Initialization_Needed() IMPLICIT NONE DIs_Initialization_Needed = (.NOT.initialized) ! from private to public END FUNCTION DIs_Initialization_Needed SUBROUTINE DPerform_Weighting (number_of_data, theta_radians, phi_radians, & ! INTENT(IN) & weights) ! INTENT(OUT) USE DSphere ! Peter Bird's Fortran 90 MODULE file DSphere.f90 IMPLICIT NONE INTEGER :: best INTEGER, INTENT(IN) :: number_of_data REAL*8 :: distance, trial_distance REAL*8, DIMENSION(:), INTENT(IN) :: theta_radians, phi_radians REAL*8, DIMENSION(:), INTENT(OUT) :: weights INTEGER :: i, iEle, j, n_allocated INTEGER, DIMENSION(:), ALLOCATABLE :: data_in_element ! one entry for each element INTEGER, DIMENSION(:), ALLOCATABLE :: element_assignment ! one entry for each datum LOGICAL :: cold_start, success REAL*8 :: element_steradians, & & s1, s2, s3, sphere_steradians REAL*8, DIMENSION(3) :: uvec REAL*8, DIMENSION(:,:), ALLOCATABLE :: uvecs IF (.NOT.initialized) THEN WRITE (*, "(' ERROR: You must CALL Initialize_Weighting(subdivision)'& &/' before you are allowed to CALL Perform_Weighting()!')") CALL DPause() STOP END IF sphere_steradians = 4.0D0 * 3.141592653589793D0 element_steradians = sphere_steradians / numEl ALLOCATE ( uvecs(3, number_of_data) ) DO j = 1, number_of_data CALL DThetaPhi_2_Uvec (theta_radians(j), phi_radians(j), uvec) uvecs(1:3, j) = uvec(1:3) END DO ALLOCATE ( data_in_element(numEl) ) data_in_element = 0 ! whole array (just initializing) ALLOCATE ( element_assignment(number_of_data) ) element_assignment = 0 ! whole array, just for tidiness ! Assign data to elements: DO j = 1, number_of_data uvec(1:3) = uvecs(1:3, j) cold_start = .TRUE. CALL DWhich_Spherical_Triangle (uvec, cold_start, & & numEl, nodes, node_uvec, center, a_, neighbor, & & success, iEle, s1, s2, s3) IF (success) THEN element_assignment(j) = iEle data_in_element(iEle) = data_in_element(iEle) + 1 ELSE WRITE (*, "(' ERROR: Datum #', I6, ' did not fit into any element of global icosahedral grid.')") j CALL DPause() STOP END IF END DO weights = 0.0D0 ! whole array; it will be used to accumulate steradians (until the very last line) !Assign element_steradians to one or more data "weights": DO i = 1, numEl IF (data_in_element(i) > 0) THEN n_allocated = 0 allocating: DO j = 1, number_of_data IF (element_assignment(j) == i) THEN weights(j) = weights(j) + element_steradians / data_in_element(i) n_allocated = n_allocated + 1 IF (n_allocated == data_in_element(i)) EXIT allocating END IF END DO allocating ELSE ! data_in_element(i) == 0 !Find nearest datum, and give it the whole solid angle. !First, find its center: uvec(1) = ( node_uvec(1, nodes(1, i)) + node_uvec(1, nodes(2, i)) + node_uvec(1, nodes(3, i)) ) / 3.0 uvec(2) = ( node_uvec(2, nodes(1, i)) + node_uvec(2, nodes(2, i)) + node_uvec(2, nodes(3, i)) ) / 3.0 uvec(3) = ( node_uvec(3, nodes(1, i)) + node_uvec(3, nodes(2, i)) + node_uvec(3, nodes(3, i)) ) / 3.0 CALL DMake_Uvec(uvec, uvec) !Second, search for datum nearest to element center: distance = 2.01D0 ! just initializing best = 0 ! ditto DO j = 1, number_of_data trial_distance = DSQRT( (uvec(1) - uvecs(1, j))**2 + & & (uvec(2) - uvecs(2, j))**2 + & & (uvec(3) - uvecs(3, j))**2 ) IF (trial_distance < distance) THEN distance = trial_distance best = j END IF END DO ! j = 1, numEl !Note: "best" is the only important information retained from this search. weights(best) = weights(best) + element_steradians END IF END DO ! i = 1, numEl weights(1:number_of_data) = weights(1:number_of_data) * (number_of_data / sphere_steradians) !Now, weights are renormalized so that their sum is number_of_data, and their mean value is 1.0D0 DEALLOCATE ( element_assignment ) DEALLOCATE ( data_in_element ) DEALLOCATE ( uvecs ) END SUBROUTINE DPerform_Weighting SUBROUTINE DDumb_s123 (element, vector, node, xyz_nod, center, a_, & & s1, s2, s3) ! Finds s1, s2, s3 coordinates of position vector "in" element ! (whether the point is actually in the element or NOT). IMPLICIT NONE INTEGER, INTENT(IN) :: element ! element # REAL*8, DIMENSION(3), INTENT(IN) :: vector ! uvec to point INTEGER, DIMENSION(:,:), INTENT(IN) :: node ! element definitions REAL*8, DIMENSION(:,:), INTENT(IN) :: xyz_nod ! uvecs of nodes REAL*8, DIMENSION(:,:), INTENT(IN) :: center ! uvec of each element (uvec) REAL*8, DIMENSION(:), INTENT(IN) :: a_ ! element areas (plane; R == 1.0) REAL*8, INTENT(OUT) :: s1, s2, s3 ! results !- - - - - - - - - - - - - - - - - - - - - - - - INTEGER :: i1, i2, i3 REAL*8, DIMENSION(3) :: tv, tvi, tvo, tv1, tv2, v1 REAL*8 :: d1, dc, t IF (element == 0) THEN WRITE (*,"(' ERROR: element = 0 in Dumb_s123')") CALL DTraceback END IF i1 = node(1, element) i2 = node(2, element) i3 = node(3, element) !shorten(?) vector to just touch plane element -> v1 tv1 = center(1:3, element) dc = DOT_PRODUCT(vector, tv1) IF (dc <= 0.0D0) THEN WRITE (*,"(' ERROR: Internal vector >= 90 deg. from element in Dumb_s123')") CALL DTraceback END IF tv2 = xyz_nod(1:3, i1) d1 = DOT_PRODUCT(tv2, tv1) t = d1 / dc v1 = t * vector tvi = xyz_nod(1:3,i3) - xyz_nod(1:3,i2) tvo = v1(1:3) - xyz_nod(1:3,i3) CALL DCross(tvi, tvo, tv) s1 = 0.5D0 * DOT_PRODUCT(tv1, tv) / a_(element) tvi = xyz_nod(1:3,i1) - xyz_nod(1:3,i3) tvo = v1(1:3) - xyz_nod(1:3,i1) CALL DCross(tvi, tvo, tv) s2 = 0.5D0 * DOT_PRODUCT(tv1, tv) / a_(element) s3 = 1.0D0 - s1 - s2 END SUBROUTINE DDumb_s123 SUBROUTINE DLearn_Spherical_Triangles (numEl, nodes, node_uvec, chatty, & & a_, center, neighbor) !Creates arrays needed by lookup subr. Which_Spherical_Triangle: ! a_ = area of plane triangle below element, when radius == 1.0 ! center = uvec pointing to center of element ! neighbor = neighboring spherical triangular elements on each ! side (or zero if none); note that algorithm ! depends on node-location match, not on node-number ! match, and therefore ignores intevening faults. !These arrays are only meaningful for finite element grids used ! with SHELLS and/or RESTORE. IMPLICIT NONE INTEGER, INTENT(IN) :: numEl ! number of spherical triangle elements INTEGER, DIMENSION(:,:), INTENT(IN) :: nodes ! element definitions REAL*8, DIMENSION(:,:), INTENT(IN) :: node_uvec ! uvecs of nodes LOGICAL, INTENT(IN) :: chatty REAL*8, DIMENSION(:), INTENT(OUT) :: a_ REAL*8, DIMENSION(:,:), INTENT(OUT) :: center INTEGER, DIMENSION(:,:), INTENT(OUT) :: neighbor INTEGER :: furthest, ia, ib, i1, i2, i3, j, j1, j2, j3, k, l_, m, step_aside REAL*8, DIMENSION(3) :: a, b, c, t, u IF (chatty) WRITE (*,"(' Learning the spherical triangles...')") furthest = (numEl + 1) / 2 ! INTEGER division is intentional neighbor = 0 ! whole array, initialized to "no neighbor on this side" homes: DO l_ = 1, numEl !first, a_ i1 = nodes(1,l_) i2 = nodes(2,l_) i3 = nodes(3,l_) a = node_uvec(1:3,i2) - node_uvec(1:3,i1) b = node_uvec(1:3,i3) - node_uvec(1:3,i2) CALL DCross (a, b, c) a_(l_) = 0.5D0 * DSQRT(c(1)**2 + c(2)**2 + c(3)**2) !second, compute center t(1:3) = (node_uvec(1:3,i1)+node_uvec(1:3,i2)+node_uvec(1:3,i3)) / 3.0D0 CALL DMake_Uvec(t, u) center(1:3, l_) = u(1:3) !third, find neighbor(?) for each side of element sides: DO j = 1, 3 ! 3 sides k = 1 + MOD (j, 3) ia = nodes(k, l_) ! 1st node along side ib = nodes(1 + MOD (k, 3), l_) ! 2nd node along side strangers: DO step_aside = 1, furthest m = l_ + step_aside ! I also try -step_aside, below m = 1 + MOD(m-1, numEl) ! wraps around j1 = nodes(1, m) j2 = nodes(2, m) j3 = nodes(3, m) IF (DSame(j1, ib) .AND. DSame(j2, ia)) THEN neighbor(j, l_) = m EXIT strangers ELSE IF (DSame(j2, ib) .AND. DSame(j3, ia)) THEN neighbor(j, l_) = m EXIT strangers ELSE IF (DSame(j3, ib) .AND. DSame(j1, ia)) THEN neighbor(j, l_) = m EXIT strangers END IF m = l_ - step_aside ! I also try +step_aside, above m = 1 + MOD(m-1+numEl, numEl) ! wraps around j1 = nodes(1, m) j2 = nodes(2, m) j3 = nodes(3, m) IF (DSame(j1, ib) .AND. DSame(j2, ia)) THEN neighbor(j, l_) = m EXIT strangers ELSE IF (DSame(j2, ib) .AND. DSame(j3, ia)) THEN neighbor(j, l_) = m EXIT strangers ELSE IF (DSame(j3, ib) .AND. DSame(j1, ia)) THEN neighbor(j, l_) = m EXIT strangers END IF END DO strangers END DO sides IF (chatty) WRITE (*,"('+Learning the spherical triangles...',I6)") l_ END DO homes IF (chatty) WRITE (*,"('+Learning the spherical triangles...DONE ')") CONTAINS LOGICAL FUNCTION DSame(i,j) ! Are node_uvec #i and #j the same vector? INTEGER, INTENT(IN) :: i, j REAL*8, PARAMETER :: epsilon = 1.0D-5 ! 64 m on Earth !the logic is: !Same = (node_uvec(1,i) == node_uvec(1,j)).AND. & ! & (node_uvec(2,i) == node_uvec(2,j)).AND. & ! & (node_uvec(3,i) == node_uvec(3,j)) !But, it is written this way for speed: IF (ABS(node_uvec(1,i) - node_uvec(1,j)) <= epsilon) THEN IF (ABS(node_uvec(2,i) - node_uvec(2,j)) <= epsilon) THEN IF (ABS(node_uvec(3,i) - node_uvec(3,j)) <= epsilon) THEN DSame = .TRUE. ELSE DSame = .FALSE. END IF ELSE DSame = .FALSE. END IF ELSE DSame = .FALSE. END IF END FUNCTION DSame END SUBROUTINE DLearn_Spherical_Triangles SUBROUTINE DPause() ! This routine is called after any error message, ! and also after successful program completion, ! in order to let the user have time to read the message ! when the job is run interactively. ! To run the job as a batch job (e.g., from a Windows .bat file), ! comment out the READ so that this routine does not stop ! the flow of the program. IMPLICIT none WRITE(*,"(' Press [Enter]...')") ! ------------------------------------- !CCCC Following line could be commented out so that jobs won't pause !CCCC on batch-processing computer systems. !CCCC Remove the CCCCC for interactive computer systems! !CCCC READ(*, *) ! ------------------------------------- END SUBROUTINE DPause SUBROUTINE DPull_in(s) ! If necessary, adjusts internal coordinates s(1..3) so ! that none is negative. IMPLICIT NONE REAL*8, DIMENSION(3), INTENT(INOUT) :: s INTEGER, DIMENSION(1) :: array ! stupid, to satisfy MINLOC REAL*8 :: factor, lowest, highest, medium INTEGER :: side, sidea, sideb lowest = MINVAL(s) IF (lowest < 0.0D0) THEN highest = MAXVAL(s) medium = 1.00D0 - lowest - highest IF (medium > 0.0D0) THEN ! correct to nearest edge array = MINLOC(s) side = array(1) s(side) = 0.0D0 sidea = 1 + MOD(side, 3) sideb = 1 + MOD(sidea, 3) factor = 1.00D0 / (1.00D0 - lowest) s(sidea) = factor * s(sidea) ! s(sideb) = factor * s(sideb) would be logical s(sideb) = 1.00D0 - s(sidea) ! is safer ELSE ! correct to nearest vertex array = MAXLOC(s) side = array(1) s(side) = 1.00D0 sidea = 1 + MOD(side, 3) sideb = 1 + MOD(sidea, 3) s(sidea) = 0.0D0 s(sideb) = 0.0D0 END IF END IF END SUBROUTINE DPull_in SUBROUTINE DWhich_Spherical_Triangle (b_, cold_start, & & num_ele, node, xyz_node, center, a_, neighbor, & & success, iEle, s1, s2, s3) !Locates a point (b_, a uvec) in element iEle with internal !coordinates (s1, s2, s3) in a SHELLS or RESTORE .feg. !and reports success. !If (cold_start), makes no use of input iEle, s1, s2, s3 !If not, uses these values to initialize the search. ! !Note that Learn_Spherical_Triangles can be used to initialize !necessary arrays a_, center, and neighbor. ! !Beware of variable name changes (numEl = num_ele, nodes = node, etc.). IMPLICIT NONE REAL*8, DIMENSION(3), INTENT(IN) :: b_ ! uvec of unknown point LOGICAL, INTENT(IN) :: cold_start ! mode switch INTEGER, INTENT(IN) :: num_ele ! count of elements INTEGER, DIMENSION(:,:),INTENT(IN) :: node ! element definitions REAL*8, DIMENSION(:,:),INTENT(IN) :: xyz_node ! uvecs of nodes REAL*8, DIMENSION(:,:),INTENT(IN) :: center ! center uvecs of elements REAL*8, DIMENSION(:), INTENT(IN) :: a_ ! (plane) areas of elements INTEGER, DIMENSION(:,:),INTENT(IN) :: neighbor ! neighbors of each element LOGICAL, INTENT(OUT) :: success ! OUTPUT INTEGER, INTENT(INOUT) :: iEle ! OUTPUT REAL*8, INTENT(INOUT) :: s1, s2, s3 ! OUTPUT !- - - - - - - - - - - - - - - - - - - - - - - - - - - - INTEGER :: back1, back2, back3, i, iet, l_ REAL*8 :: r2, r2min, s1t, s2t, s3t REAL*8, DIMENSION(3) :: s_temp, tv ! establish defaults (not found) in case of quick exit success = .FALSE. IF (cold_start) THEN iEle = 0 ! default s1 = 0.0D0; s2 = 0.0D0; s3 = 0.0D0 ! default !find closest element center to initialize search r2min = 4.01D0 ! radians DO l_ = 1, num_ele r2 = (b_(1) - center(1,l_))**2 +(b_(2) - center(2,l_))**2 +(b_(3) - center(3,l_))**2 IF (r2 < r2min) THEN r2min = r2 iet = l_ END IF END DO ! If closest element center is more than 1 radian away, give uP. tv = center(1:3, iet) IF (DOT_PRODUCT(b_, tv) < 0.540D0) RETURN END IF ! cold_start ! initialize search memory (with impossible numbers) back1 = -1 back2 = -2 back3 = -3 is_it_here: DO ! first, check for infinite loop between 2 elements! IF (iet == back2) THEN ! in loop; force location in one or the other! CALL DDumb_s123 (iet, b_, node, xyz_node, center, a_, & & s1t, s2t, s3t) s_temp(1) = s1t; s_temp(2) = s2t; s_temp(3) = s3t CALL DPull_in(s_temp) s1t = s_temp(1); s2t = s_temp(2); s3t = s_temp(3) EXIT is_it_here ! then, check for infinite loop between 3 elements! ELSE IF (iet == back3) THEN ! in loop; force location in one or the other! CALL DDumb_s123 (iet, b_, node, xyz_node, center, a_, & & s1t, s2t, s3t) s_temp(1) = s1t; s_temp(2) = s2t; s_temp(3) = s3t CALL DPull_in(s_temp) s1t = s_temp(1); s2t = s_temp(2); s3t = s_temp(3) EXIT is_it_here ELSE ! normal operation CALL DDumb_s123 (iet, b_, node, xyz_node, center, a_, & & s1t, s2t, s3t) IF ((s1t < s2t) .AND. (s1t < s3t)) THEN ! s1 is most negative; most critical IF (s1t >= 0.0D0) THEN EXIT is_it_here ! success ELSE i = neighbor(1, iet) IF (i > 0) THEN back3 = back2 back2 = back1 back1 = iet iet = i CYCLE is_it_here ELSE RETURN ! fell off edge of grid ENDIF ENDIF ELSE IF ((s2t < s1t) .AND. (s2t < s3t)) THEN ! s2 is most negative; most critical IF (s2t >= 0.0D0) THEN EXIT is_it_here ! success ELSE i = neighbor(2, iet) IF (i > 0) THEN back3 = back2 back2 = back1 back1 = iet iet = i CYCLE is_it_here ELSE RETURN ! fell off edge of grid ENDIF ENDIF ELSE ! s3 is most negative; most critical IF (s3t >= 0.0D0) THEN EXIT is_it_here ! success ELSE i = neighbor(3, iet) IF (i > 0) THEN back3 = back2 back2 = back1 back1 = iet iet = i CYCLE is_it_here ELSE RETURN ! fell off edge of grid ENDIF ENDIF END IF END IF ! in/not in a loop END DO is_it_here ! successful completion iEle = iet s1 = s1t s2 = s2t s3 = s3t success = .TRUE. END SUBROUTINE DWhich_Spherical_Triangle END MODULE DWeighting