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

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

dynamico tree creation

YM
File size: 8.0 KB
Line 
1MODULE geometry
2  USE field_mod
3 
4  TYPE t_geometry
5    TYPE(t_field),POINTER :: xyz_i(:)
6    TYPE(t_field),POINTER :: xyz_e(:)
7    TYPE(t_field),POINTER :: xyz_v(:)
8    TYPE(t_field),POINTER :: Ai(:)
9    TYPE(t_field),POINTER :: Av(:)
10    TYPE(t_field),POINTER :: de(:)
11    TYPE(t_field),POINTER :: le(:)
12    TYPE(t_field),POINTER :: Riv(:)
13    TYPE(t_field),POINTER :: Riv2(:)
14    TYPE(t_field),POINTER :: ne(:)
15    TYPE(t_field),POINTER :: Wee(:)
16    TYPE(t_field),POINTER :: bi(:)
17    TYPE(t_field),POINTER :: fv(:)
18  END TYPE t_geometry   
19
20  TYPE(t_geometry),TARGET :: geom
21 
22  REAL(rstd),POINTER :: Ai(:)
23  REAL(rstd),POINTER :: xyz_i(:,:)
24  REAL(rstd),POINTER :: xyz_e(:,:)
25  REAL(rstd),POINTER :: xyz_v(:,:)
26  REAL(rstd),POINTER :: Av(:)
27  REAL(rstd),POINTER :: de(:)
28  REAL(rstd),POINTER :: le(:)
29  REAL(rstd),POINTER :: Riv(:,:)
30  REAL(rstd),POINTER :: Riv2(:,:)
31  INTEGER,POINTER :: ne(:,:)
32  REAL(rstd),POINTER :: Wee(:,:,:)
33  REAL(rstd),POINTER :: bi(:)
34  REAL(rstd),POINTER :: fv(:)
35     
36
37
38CONTAINS
39
40  SUBROUTINE allocate_geometry
41  USE field_mod
42  IMPLICIT NONE
43 
44    CALL allocate_field(geom%Ai,field_t,type_real)
45    CALL allocate_field(geom%xyz_i,field_t,type_real,3)
46    CALL allocate_field(geom%xyz_e,field_u,type_real,3)
47    CALL allocate_field(geom%xyz_v,field_z,type_real,3)
48    CALL allocate_field(geom%de,field_u,type_real)
49    CALL allocate_field(geom%le,field_u,type_real)
50    CALL allocate_field(geom%bi,field_t,type_real)
51    CALL allocate_field(geom%Av,field_z,type_real)
52    CALL allocate_field(geom%Riv,field_t,type_real,6)
53    CALL allocate_field(geom%Riv2,field_t,type_real,6)
54    CALL allocate_field(geom%ne,field_t,type_integer,6)
55    CALL allocate_field(geom%Wee,field_u,type_real,5,2)
56    CALL allocate_field(geom%bi,field_t,type_real)
57    CALL allocate_field(geom%fv,field_z,type_real)
58
59  END SUBROUTINE allocate_geometry
60 
61 
62  SUBROUTINE swap_geometry(ind)
63  USE field_mod
64  IMPLICIT NONE
65    INTEGER,INTENT(IN) :: ind
66    Ai=geom%Ai(ind)
67    xyz_i=geom%xyz_i(ind)
68    xyz_e=geom%xyz_e(ind)
69    xyz_v=geom%xyz_v(ind)
70    de=geom%de(ind)
71    le=geom%le(ind)
72    Av=geom%Av(ind)
73    Riv=geom%Riv(ind)
74    Riv2=geom%Riv2(ind)
75    ne=geom%ne(ind)
76    Wee=geom%Wee(ind)
77    bi=geom%bi(ind)
78    fv=geom%fv(ind)
79   
80  END SUBROUTINE swap_geometry
81   
82  SUBROUTINE set_geometry
83  USE metric
84  USE spherical_geom_mod
85  USE domain_mod
86  USE dimensions
87  IMPLICIT NONE
88
89    REAL(rstd) :: surf(6)
90    REAL(rstd) :: surf_v(6)
91    INTEGER :: ind,i,j,k,n
92    TYPE(t_domain),POINTER :: d
93    REAL(rstd) ::  S
94    REAL(rstd) :: w(6)
95    REAl(rstd) :: lon,lat
96    INTEGER :: ii_glo,jj_glo
97    REAL(rstd) :: S1,S2
98     
99    DO ind=1,ndomain
100      d=>domain(ind)
101      CALL swap_dimensions(ind)
102      CALL swap_geometry(ind)
103      DO j=jj_begin-1,jj_end+1
104        DO i=ii_begin-1,ii_end+1
105          n=(j-1)*iim+i
106       
107          xyz_i(n,:)=d%xyz(:,i,j) 
108          xyz_v(n+z_up,:)=d%vertex(:,vup-1,i,j) 
109          xyz_v(n+z_down,:)=d%vertex(:,vdown-1,i,j) 
110         
111          CALL xyz2lonlat(xyz_v(n+z_up,:),lon,lat)
112          fv(n+z_up)=2*sin(lat)*omega
113          CALL xyz2lonlat(xyz_v(n+z_down,:),lon,lat)
114          fv(n+z_down)=2*sin(lat)*omega
115         
116          bi(n)=0. 
117       
118          CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,right-1,i,j),de(n+u_right))
119          CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,lup-1,i,j),de(n+u_lup))
120          CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,ldown-1,i,j),de(n+u_ldown))
121         
122          CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,right-1,i,j),0.5,xyz_e(n+u_right,:))
123          CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,lup-1,i,j),0.5,xyz_e(n+u_lup,:))
124          CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,ldown-1,i,j),0.5,xyz_e(n+u_ldown,:))
125         
126          CALL dist_cart(d%vertex(:,vrdown-1,i,j),d%vertex(:,vrup-1,i,j),le(n+u_right))
127          CALL dist_cart(d%vertex(:,vup-1,i,j),d%vertex(:,vlup-1,i,j),le(n+u_lup))
128          CALL dist_cart(d%vertex(:,vldown-1,i,j),d%vertex(:,vdown-1,i,j),le(n+u_ldown))
129
130         
131          Ai(n)=0
132          DO k=0,5
133            CALL surf_triangle(d%xyz(:,i,j),d%neighbour(:,k,i,j),d%neighbour(:,MOD((k+1+6),6),i,j),surf_v(k+1))
134            CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,MOD((k-1+6),6),i,j),d%vertex(:,k,i,j),surf(k+1))
135            ne(n,k+1)=d%ne(k,i,j)
136            Ai(n)=Ai(n)+surf(k+1)
137          ENDDO
138
139
140          DO k=0,5
141            CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,k,i,j),d%neighbour(:,k,i,j),S1)
142            CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,k,i,j),d%neighbour(:,MOD(k+1+6,6),i,j),S2)
143            Riv(n,k+1)=0.5*(S1+S2)/Ai(n)
144            Riv2(n,k+1)=0.5*(S1+S2)/surf_v(k+1)
145          ENDDO
146
147          DO k=1,6
148            IF (ABS(surf_v(k))<1e-30) THEN
149              Riv(n,k)=0.
150            ENDIF
151          ENDDO
152         
153          Av(n+z_up)=surf_v(vup)+1e-100
154          Av(n+z_down)=surf_v(vdown)+1e-100
155       
156        ENDDO
157      ENDDO
158     
159      DO j=jj_begin,jj_end
160        DO i=ii_begin,ii_end
161          n=(j-1)*iim+i
162   
163          CALL compute_wee(n,right,w)
164          Wee(n+u_right,:,1)=w(1:5)
165
166          CALL compute_wee(n+t_right,left,w)
167          Wee(n+u_right,:,2)=w(1:5)
168
169
170          CALL compute_wee(n,lup,w)
171          Wee(n+u_lup,:,1)=w(1:5)
172
173          CALL compute_wee(n+t_lup,rdown,w)
174          Wee(n+u_lup,:,2)=w(1:5)
175
176
177          CALL compute_wee(n,ldown,w)
178          Wee(n+u_ldown,:,1)=w(1:5)
179
180          CALL compute_wee(n+t_ldown,rup,w)
181          Wee(n+u_ldown,:,2)=w(1:5)
182
183        ENDDO
184      ENDDO
185     
186      DO j=jj_begin,jj_end
187        DO i=ii_begin,ii_end
188          n=(j-1)*iim+i
189          ii_glo=d%ii_begin_glo-d%ii_begin+i
190          jj_glo=d%jj_begin_glo-d%jj_begin+j
191         
192          IF (ii_glo==1 .AND. jj_glo==1) THEN
193            le(n+u_ldown)=0
194            xyz_v(n+z_ldown,:)=xyz_v(n+z_down,:)
195                       
196          ENDIF
197
198          IF (ii_glo==iim_glo .AND. jj_glo==1) THEN
199            le(n+u_right)=0
200            xyz_v(n+z_rdown,:)=xyz_v(n+z_rup,:)
201          ENDIF
202
203          IF (ii_glo==iim_glo .AND. jj_glo==jjm_glo) THEN
204            le(n+u_rup)=0
205            xyz_v(n+z_rup,:)=xyz_v(n+z_up,:)
206          ENDIF
207
208          IF (ii_glo==1 .AND. jj_glo==jjm_glo) THEN
209            le(n+u_lup)=0
210            xyz_v(n+z_up,:)=xyz_v(n+z_lup,:)
211          ENDIF
212         
213        ENDDO
214      ENDDO
215
216      DO j=jj_begin-1,jj_end+1
217        DO i=ii_begin-1,ii_end+1
218          n=(j-1)*iim+i
219          xyz_i(n,:)=xyz_i(n,:) * radius
220          xyz_v(n+z_up,:)=xyz_v(n+z_up,:) * radius
221          xyz_v(n+z_down,:)=xyz_v(n+z_down,:) *radius
222          de(n+u_right)=de(n+u_right) * radius
223          de(n+u_lup)=de(n+u_lup)*radius
224          de(n+u_ldown)=de(n+u_ldown)*radius
225          xyz_e(n+u_right,:)=xyz_e(n+u_right,:)*radius
226          xyz_e(n+u_lup,:)=xyz_e(n+u_lup,:)*radius
227          xyz_e(n+u_ldown,:)=xyz_e(n+u_ldown,:)*radius
228          le(n+u_right)=le(n+u_right)*radius
229          le(n+u_lup)=le(n+u_lup)*radius
230          le(n+u_ldown)=le(n+u_ldown)*radius
231          Ai(n)=Ai(n)*radius**2
232          Av(n+z_up)=Av(n+z_up)*radius**2
233          Av(n+z_down)=Av(n+z_down)*radius**2
234        ENDDO
235      ENDDO
236                 
237    ENDDO
238    CALL surf_triangle(d%xyz(:,ii_begin,jj_begin),d%xyz(:,ii_begin,jj_end),d%xyz(:,ii_end,jj_begin),S)
239!    PRINT *,"Surf triangle : ",S*20/(4*Pi)
240 
241  END SUBROUTINE set_geometry
242 
243  SUBROUTINE compute_wee(n,pos,w)
244  IMPLICIT NONE
245    INTEGER,INTENT(IN) :: n
246    INTEGER,INTENT(IN) :: pos
247    REAL(rstd),INTENT(OUT) ::w(6)
248
249    REAL(rstd) :: ne_(0:5)
250    REAL(rstd) :: Riv_(6)
251    INTEGER    :: k
252   
253   
254    DO k=0,5
255      ne_(k)=ne(n,MOD(pos-1+k+6,6)+1) 
256      Riv_(k+1)=Riv(n,MOD(pos-1+k+6,6)+1)
257    ENDDO
258         
259    w(1)=-ne_(0)*ne_(1)*(Riv_(1)-0.5)
260    w(2)=-ne_(2)*(ne_(0)*Riv_(2)-w(1)*ne_(1))
261    w(3)=-ne_(3)*(ne_(0)*Riv_(3)-w(2)*ne_(2))
262    w(4)=-ne_(4)*(ne_(0)*Riv_(4)-w(3)*ne_(3))
263    w(5)=-ne_(5)*(ne_(0)*Riv_(5)-w(4)*ne_(4))
264    w(6)=ne_(0)*ne_(5)*(Riv_(6)-0.5)
265   
266!    IF ( ABS(w(5)-w(6))>1e-20) PRINT *, "pb pour wee : w(5)!=w(6)",sum(Riv_(:))
267
268   END SUBROUTINE compute_wee
269           
270
271 
272  SUBROUTINE compute_geometry
273  IMPLICIT NONE
274    CALL allocate_geometry
275    CALL set_geometry
276   
277  END SUBROUTINE compute_geometry
278 
279END MODULE geometry
Note: See TracBrowser for help on using the repository browser.