Changeset 15 for codes/icosagcm/trunk/src/spherical_geom.f90

Timestamp:

07/09/12 15:23:38 (12 years ago)

Author:

ymipsl

Message:

Update on 3D dynamic

YM

File:

• codes/icosagcm/trunk/src/spherical_geom.f90 (modified) (4 diffs)

codes/icosagcm/trunk/src/spherical_geom.f90

@@ -78 +78 @@
 
   REAL(rstd)  :: AB,AC,BC
-  REAL(rstd)  :: s
+  REAL(rstd)  :: s,x
   
   CALL dist_cart(A,B,AB)
@@ -85 +85 @@
   
   s=(AB+AC+BC)/2
-  surf=4*atan(sqrt( tan(s/2) * tan((s-AB)/2)  * tan((s-AC)/2) * tan((s-BC)/2)))
+  x=tan(s/2) * tan((s-AB)/2)  * tan((s-AC)/2) * tan((s-BC)/2)
+  IF (x<0) x=0.
+  surf=4*atan(sqrt( x))
   
 END SUBROUTINE surf_triangle 
@@ -142 +144 @@
 
 
-  SUBROUTINE circonscrit(A0,B0,C0,Center)
+  SUBROUTINE circumcenter(A0,B0,C0,Center)
   USE vector
   IMPLICIT NONE
@@ -157 +159 @@
     center=center/sqrt(sum(center**2))
     
-  END SUBROUTINE circonscrit
+  END SUBROUTINE circumcenter
     
 
-  SUBROUTINE centroide(A,B,C,Center)
+  SUBROUTINE compute_centroid(points,n,centr)
+  USE vector
   IMPLICIT NONE
-    REAL(rstd), INTENT(IN)  :: A(3),B(3),C(3)
-    REAL(rstd), INTENT(OUT) :: Center(3)
+    INTEGER :: n
+    REAL(rstd), INTENT(IN)  :: points(3,n)
+    REAL(rstd), INTENT(OUT) :: Centr(3)
+    
+    REAL(rstd) :: p1(3),p2(3),cross(3)
+    REAL(rstd) :: norm_cross
+    INTEGER :: i,j
+      
+      Centr(:)=0
+      DO i=1,n
+        j=MOD(i,n)+1
+        p1=points(:,i)/norm(points(:,i))
+        p2=points(:,j)/norm(points(:,j))
+        CALL cross_product2(p1,p2,cross)
+        norm_cross=norm(cross)
+        if (norm_cross<1e-10) CYCLE 
+        
+        Centr(:)=centr(:)+asin(norm_cross)*cross(:)/norm_cross
+      ENDDO
+      
+      Centr(:)=centr(:)/norm(centr(:))
   
-    REAL(rstd)  :: AB05(3),AC05(3),BC05(3)
-    REAL(rstd)  :: Vab(3),Vbc(3),Vca(3)
-    REAL(rstd)  :: t
-    INTEGER :: n1,n2
-       
-    AB05=(A+B)/2.
-    AC05=(A+C)/2.
-    BC05=(B+C)/2.
-    
-    CALL director(A,B,C,Vab)
-!    PRINT *,'director',sum((B-A)*Vab)
-    CALL director(B,C,A,Vbc)
-!    PRINT *,'director',sum((C-B)*Vbc)
-    CALL director(C,A,B,Vca)
-!    PRINT *,'director',sum((A-C)*Vca)
-   
-    n1=1 ; n2=2
-    
-    t=(Vbc(n1)*(BC05(n2)-AB05(n2))+(AB05(n1)-BC05(n1))*Vbc(n2))/(Vbc(n1)*Vab(n2)-Vab(n1)*Vbc(n2))
-    Center=AB05+t*Vab
-    Center=Center/sqrt(sum(Center**2))
-!    PRINT*,"Center :",A
-!    PRINT*,"Center :",B
-!    PRINT*,"Center :",C
-!    PRINT*,"Center :",Center
-   
-  CONTAINS 
-  
-    SUBROUTINE director(A,B,C,V)
-    IMPLICIT NONE
-    REAL(rstd), INTENT(IN)  :: A(3),B(3),C(3)
-    REAL(rstd), INTENT(OUT) :: V(3)
-    REAL(rstd) :: AB(3),AC(3)
-    REAL(rstd) :: t
-    
-      AB=B-A
-      AC=C-A
-    
-      t=sum(AB*AC)/sum(AB**2)
-      V=AB*t-AC
-    END SUBROUTINE director
-       
-  END SUBROUTINE centroide
+  END SUBROUTINE compute_centroid
 
 END MODULE spherical_geom_mod