source: codes/icosagcm/trunk/src/spherical_geom.f90 @ 12

Last change on this file since 12 was 12, checked in by ymipsl, 12 years ago

dynamico tree creation

YM
File size: 4.4 KB
Line 
1MODULE spherical_geom_mod
2USE genmod
3
4
5
6CONTAINS
7
8
9
10
11SUBROUTINE lonlat2xyz(lon,lat,xyz)
12IMPLICIT NONE
13  REAL(rstd),INTENT(IN) :: lon
14  REAL(rstd),INTENT(IN) :: lat
15  REAL(rstd),INTENT(OUT) :: xyz(3)
16 
17  xyz(1)=cos(lon)*cos(lat)
18  xyz(2)=sin(lon)*cos(lat)
19  xyz(3)=sin(lat)
20
21END SUBROUTINE lonlat2xyz
22
23
24SUBROUTINE xyz2lonlat(xyz,lon,lat)
25IMPLICIT NONE
26  REAL(rstd),INTENT(IN) :: xyz(3)
27  REAL(rstd),INTENT(OUT) :: lon
28  REAL(rstd),INTENT(OUT) :: lat
29 
30  REAL(rstd) :: coslat
31  REAL(rstd) :: xyzn(3)
32 
33  xyzn(:)=xyz(:)/sqrt(sum(xyz(:)**2))
34 
35  lat=asin(xyzn(3))
36  coslat=cos(lat)
37  IF (abs(coslat)<1e-15) THEN
38    lon=0.
39  ELSE
40    lon=acos(MAX(MIN((xyzn(1)/coslat),1.),-1.))
41    IF (xyzn(2)/coslat<0) lon=-lon
42  ENDIF
43   
44END SUBROUTINE xyz2lonlat
45
46
47SUBROUTINE dist_cart(A,B,d)
48USE vector
49IMPLICIT NONE
50  REAL(rstd),INTENT(IN)  :: A(3)
51  REAL(rstd),INTENT(IN)  :: B(3)
52  REAL(rstd),INTENT(OUT) :: d
53 
54   REAL(rstd)  :: n(3)
55   CALL cross_product2(A,B,n)
56   d=asin(sqrt(sum(n**2)))
57
58END SUBROUTINE dist_cart
59
60
61SUBROUTINE dist_lonlat(lonA,latA,lonB,latB,d)
62IMPLICIT NONE
63  REAL(rstd),INTENT(IN)  :: lonA
64  REAL(rstd),INTENT(IN)  :: latA
65  REAL(rstd),INTENT(IN)  :: lonB
66  REAL(rstd),INTENT(IN)  :: latB
67  REAL(rstd),INTENT(OUT) :: d
68 
69  d=acos(MAX(MIN(sin(latA)*sin(latB)+cos(latA)*cos(latB)*cos(lonA-lonB),1.),-1.))
70 
71END SUBROUTINE dist_lonlat
72
73SUBROUTINE surf_triangle(A,B,C,surf)
74  REAL(rstd),INTENT(IN)  :: A(3)
75  REAL(rstd),INTENT(IN)  :: B(3)
76  REAL(rstd),INTENT(IN)  :: C(3)
77  REAL(rstd),INTENT(OUT) :: Surf
78
79  REAL(rstd)  :: AB,AC,BC
80  REAL(rstd)  :: s
81 
82  CALL dist_cart(A,B,AB)
83  CALL dist_cart(A,C,AC)
84  CALL dist_cart(B,C,BC)
85 
86  s=(AB+AC+BC)/2
87  surf=4*atan(sqrt( tan(s/2) * tan((s-AB)/2)  * tan((s-AC)/2) * tan((s-BC)/2)))
88 
89END SUBROUTINE surf_triangle 
90
91
92SUBROUTINE div_arc(A,B,frac,C)
93IMPLICIT NONE
94  REAL(rstd),INTENT(IN)  :: A(3)
95  REAL(rstd),INTENT(IN)  :: B(3)
96  REAL(rstd),INTENT(IN)  :: frac
97  REAL(rstd),INTENT(OUT)  :: C(3)
98 
99  REAL(rstd) :: d
100  REAL(rstd) :: M(3,3)
101  REAL(rstd) :: alpha(3,3)
102  INTEGER    :: IPIV(3)
103  INTEGER    :: info
104  REAL(rstd) :: xa,xb,xc
105  REAL(rstd) :: ya,yb,yc
106  REAL(rstd) :: za,zb,zc
107  REAL(rstd) :: alpha_A,alpha_B,alpha_C
108  REAL(rstd) :: x,y,z
109  REAL(rstd) :: a1,a2,a3
110  REAL(rstd) :: b1,b2,b3
111 
112 
113  xa=A(1) ; ya=A(2) ; za=A(3)
114  xb=B(1) ; yb=B(2) ; zb=B(3)
115
116  CALL dist_cart(A,B,d)
117
118  C(1)=cos(frac*d) 
119  C(2)=cos((1-frac)*d)
120  C(3)=0.
121
122  xc=ya*zb-yb*za ; yc=-(xa*zb-xb*za) ; zc=xa*yb-xb*ya
123 
124  M(1,1)=xa ; M(1,2)=ya ; M(1,3)=za
125  M(2,1)=xb ; M(2,2)=yb ; M(2,3)=zb
126  M(3,1)=xc ; M(3,2)=yc ; M(3,3)=zc
127 
128  CALL DGESV(3,1,M,3,IPIV,C,3,info)
129 
130END SUBROUTINE div_arc
131
132SUBROUTINE div_arc_bis(A,B,frac,C)
133IMPLICIT NONE
134  REAL(rstd),INTENT(IN)  :: A(3)
135  REAL(rstd),INTENT(IN)  :: B(3)
136  REAL(rstd),INTENT(IN)  :: frac
137  REAL(rstd),INTENT(OUT)  :: C(3)
138 
139   C=A*(1-frac)+B*frac 
140   C=C/sqrt(sum(C**2))
141END SUBROUTINE div_arc_bis
142
143
144  SUBROUTINE circonscrit(A0,B0,C0,Center)
145  USE vector
146  IMPLICIT NONE
147    REAL(rstd), INTENT(IN)  :: A0(3),B0(3),C0(3)
148    REAL(rstd), INTENT(OUT) :: Center(3)
149   
150    REAL(rstd)  :: a(3),b(3),c(3)
151   
152    a=A0/sqrt(sum(A0**2))
153    b=B0/sqrt(sum(B0**2))
154    c=C0/sqrt(sum(C0**2))
155   
156    CALL Cross_product2(b-a,c-b,center)
157    center=center/sqrt(sum(center**2))
158   
159  END SUBROUTINE circonscrit
160   
161
162  SUBROUTINE centroide(A,B,C,Center)
163  IMPLICIT NONE
164    REAL(rstd), INTENT(IN)  :: A(3),B(3),C(3)
165    REAL(rstd), INTENT(OUT) :: Center(3)
166 
167    REAL(rstd)  :: AB05(3),AC05(3),BC05(3)
168    REAL(rstd)  :: Vab(3),Vbc(3),Vca(3)
169    REAL(rstd)  :: t
170    INTEGER :: n1,n2
171       
172    AB05=(A+B)/2.
173    AC05=(A+C)/2.
174    BC05=(B+C)/2.
175   
176    CALL director(A,B,C,Vab)
177!    PRINT *,'director',sum((B-A)*Vab)
178    CALL director(B,C,A,Vbc)
179!    PRINT *,'director',sum((C-B)*Vbc)
180    CALL director(C,A,B,Vca)
181!    PRINT *,'director',sum((A-C)*Vca)
182   
183    n1=1 ; n2=2
184   
185    t=(Vbc(n1)*(BC05(n2)-AB05(n2))+(AB05(n1)-BC05(n1))*Vbc(n2))/(Vbc(n1)*Vab(n2)-Vab(n1)*Vbc(n2))
186    Center=AB05+t*Vab
187    Center=Center/sqrt(sum(Center**2))
188!    PRINT*,"Center :",A
189!    PRINT*,"Center :",B
190!    PRINT*,"Center :",C
191!    PRINT*,"Center :",Center
192   
193  CONTAINS
194 
195    SUBROUTINE director(A,B,C,V)
196    IMPLICIT NONE
197    REAL(rstd), INTENT(IN)  :: A(3),B(3),C(3)
198    REAL(rstd), INTENT(OUT) :: V(3)
199    REAL(rstd) :: AB(3),AC(3)
200    REAL(rstd) :: t
201   
202      AB=B-A
203      AC=C-A
204   
205      t=sum(AB*AC)/sum(AB**2)
206      V=AB*t-AC
207    END SUBROUTINE director
208       
209  END SUBROUTINE centroide
210
211END MODULE spherical_geom_mod
212
213
Note: See TracBrowser for help on using the repository browser.