Changeset 247 for codes/icosagcm/trunk
- Timestamp:
- 07/22/14 16:27:58 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/spherical_geom.f90
r155 r247 179 179 180 180 181 SUBROUTINE circumcenter(A0,B0,C0, Center)181 SUBROUTINE circumcenter(A0,B0,C0,center) 182 182 USE vector 183 183 IMPLICIT NONE … … 185 185 REAL(rstd), INTENT(OUT) :: Center(3) 186 186 187 REAL(rstd) :: a(3),b(3),c(3) 187 REAL(rstd) :: a(3),b(3),c(3), ac(3), ab(3), p1(3), q(3), p2(3) 188 188 189 189 a=A0/sqrt(sum(A0**2)) … … 191 191 c=C0/sqrt(sum(C0**2)) 192 192 193 CALL Cross_product2(b-a,c-b,center) 194 center=center/sqrt(sum(center**2)) 195 193 ab=b-a 194 ac=c-a 195 CALL Cross_product2(ab,ac,p1) 196 IF(.FALSE.) THEN ! Direct solution, round-off error 197 center=p1/norm(p1) 198 ELSE ! Two-step solution, stable 199 q = SUM(ac**2)*ab-SUM(ab**2)*ac 200 CALL Cross_product2(p1,q,p2) 201 p2 = a + p2/(2.*SUM(p1**2)) 202 center = p2/norm(p2) 203 END IF 196 204 END SUBROUTINE circumcenter 197 205 … … 204 212 REAL(rstd), INTENT(OUT) :: Centr(3) 205 213 206 REAL(rstd) :: p1(3),p2(3),cross(3) 207 REAL(rstd) :: norm_cross 214 REAL(rstd) :: p1(3),p2(3),cross(3), cc(3) 215 REAL(rstd) :: norm_cross, area 208 216 INTEGER :: i,j 209 210 Centr(:)=0 211 DO i=1,n 212 j=MOD(i,n)+1 213 p1=points(:,i)/norm(points(:,i)) 214 p2=points(:,j)/norm(points(:,j)) 215 CALL cross_product2(p1,p2,cross) 216 norm_cross=norm(cross) 217 if (norm_cross<1e-10) CYCLE 218 219 Centr(:)=centr(:)+asin(norm_cross)*cross(:)/norm_cross 220 ENDDO 221 222 Centr(:)=centr(:)/norm(centr(:)) 223 217 218 centr(:)=0 219 IF(.FALSE.) THEN 220 ! Gauss formula (subject to round-off error) 221 DO i=1,n 222 j=MOD(i,n)+1 223 p1=points(:,i)/norm(points(:,i)) 224 p2=points(:,j)/norm(points(:,j)) 225 CALL cross_product2(p1,p2,cross) 226 norm_cross=norm(cross) 227 if (norm_cross<1e-10) CYCLE 228 centr(:)=centr(:)+asin(norm_cross)*cross(:)/norm_cross 229 ENDDO 230 ELSE 231 ! Simple area-weighted average (second-order accurate, stable) 232 cc=SUM(points,2) ! arithmetic average used as center 233 cc=cc/norm(cc) 234 DO i=1,n 235 j=MOD(i,n)+1 236 p1=points(:,i)/norm(points(:,i)) 237 p2=points(:,j)/norm(points(:,j)) 238 CALL surf_triangle(cc,p1,p2,area) 239 centr(:)=centr(:)+area*(p1+p2+cc) 240 ENDDO 241 END IF 242 243 centr(:)=centr(:)/norm(centr(:)) 244 224 245 END SUBROUTINE compute_centroid 225 246
Note: See TracChangeset
for help on using the changeset viewer.