Context Navigation

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

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

Update_geometry bug fix : global data field value were not broadcasted from root process to others.

YM

File size: 22.0 KB

Line
1	MODULE geometry
2	USE field_mod
3	IMPLICIT NONE
4
5	TYPE t_geometry
6	TYPE(t_field),POINTER :: centroid(:)
7	TYPE(t_field),POINTER :: xyz_i(:)
8	TYPE(t_field),POINTER :: xyz_e(:)
9	TYPE(t_field),POINTER :: xyz_v(:)
10	TYPE(t_field),POINTER :: lon_i(:)
11	TYPE(t_field),POINTER :: lon_e(:)
12	TYPE(t_field),POINTER :: lat_i(:)
13	TYPE(t_field),POINTER :: lat_e(:)
14	TYPE(t_field),POINTER :: ep_e(:)
15	TYPE(t_field),POINTER :: et_e(:)
16	TYPE(t_field),POINTER :: elon_i(:)
17	TYPE(t_field),POINTER :: elat_i(:)
18	TYPE(t_field),POINTER :: elon_e(:)
19	TYPE(t_field),POINTER :: elat_e(:)
20	TYPE(t_field),POINTER :: Ai(:)
21	TYPE(t_field),POINTER :: Av(:)
22	TYPE(t_field),POINTER :: de(:)
23	TYPE(t_field),POINTER :: le(:)
24	TYPE(t_field),POINTER :: le_de(:) ! le/de, 0. if de=0.
25	TYPE(t_field),POINTER :: Riv(:)
26	TYPE(t_field),POINTER :: Riv2(:)
27	TYPE(t_field),POINTER :: ne(:)
28	TYPE(t_field),POINTER :: Wee(:)
29	TYPE(t_field),POINTER :: bi(:)
30	TYPE(t_field),POINTER :: fv(:)
31	END TYPE t_geometry
32
33	TYPE(t_geometry),SAVE,TARGET :: geom
34
35
36	REAL(rstd),POINTER :: Ai(:) ! area of a cell
37	!$OMP THREADPRIVATE(Ai)
38	REAL(rstd),POINTER :: centroid(:,:) ! coordinate of the centroid of the cell
39	!$OMP THREADPRIVATE(centroid)
40	REAL(rstd),POINTER :: xyz_i(:,:) ! coordinate of the center of the cell (voronoi)
41	!$OMP THREADPRIVATE(xyz_i)
42	REAL(rstd),POINTER :: xyz_e(:,:) ! coordinate of a wind point on the cell on a edge
43	!$OMP THREADPRIVATE(xyz_e)
44	REAL(rstd),POINTER :: xyz_v(:,:) ! coordinate of a vertex (center of the dual mesh)
45	!$OMP THREADPRIVATE(xyz_v)
46	REAL(rstd),POINTER :: lon_i(:) ! longitude of the center of the cell (voronoi)
47	!$OMP THREADPRIVATE(lon_i)
48	REAL(rstd),POINTER :: lon_e(:) ! longitude of a wind point on the cell on a edge
49	!$OMP THREADPRIVATE(lon_e)
50	REAL(rstd),POINTER :: lat_i(:) ! latitude of the center of the cell (voronoi)
51	!$OMP THREADPRIVATE(lat_i)
52	REAL(rstd),POINTER :: lat_e(:) ! latitude of a wind point on the cell on a edge
53	!$OMP THREADPRIVATE(lat_e)
54	REAL(rstd),POINTER :: ep_e(:,:) ! perpendicular unit vector of a edge (outsider)
55	!$OMP THREADPRIVATE(ep_e)
56	REAL(rstd),POINTER :: et_e(:,:) ! tangeantial unit vector of a edge
57	!$OMP THREADPRIVATE(et_e)
58	REAL(rstd),POINTER :: elon_i(:,:) ! unit longitude vector on the center
59	!$OMP THREADPRIVATE(elon_i)
60	REAL(rstd),POINTER :: elat_i(:,:) ! unit latitude vector on the center
61	!$OMP THREADPRIVATE(elat_i)
62	REAL(rstd),POINTER :: elon_e(:,:) ! unit longitude vector on a wind point
63	!$OMP THREADPRIVATE(elon_e)
64	REAL(rstd),POINTER :: elat_e(:,:) ! unit latitude vector on a wind point
65	!$OMP THREADPRIVATE(elat_e)
66	REAL(rstd),POINTER :: Av(:) ! area of dual mesk cell
67	!$OMP THREADPRIVATE(Av)
68	REAL(rstd),POINTER :: de(:) ! distance from a neighbour == lenght of an edge of the dual mesh
69	!$OMP THREADPRIVATE(de)
70	REAL(rstd),POINTER :: le(:) ! lenght of a edge
71	!$OMP THREADPRIVATE(le)
72	REAL(rstd),POINTER :: le_de(:) ! le/de
73	!$OMP THREADPRIVATE(le_de)
74	REAL(rstd),POINTER :: Riv(:,:) ! weight
75	!$OMP THREADPRIVATE(Riv)
76	REAL(rstd),POINTER :: Riv2(:,:) ! weight
77	!$OMP THREADPRIVATE(Riv2)
78	INTEGER,POINTER :: ne(:,:) ! convention for the way on the normal wind on an edge
79	!$OMP THREADPRIVATE(ne)
80	REAL(rstd),POINTER :: Wee(:,:,:) ! weight
81	!$OMP THREADPRIVATE(Wee)
82	REAL(rstd),POINTER :: bi(:) ! orographie
83	!$OMP THREADPRIVATE(bi)
84	REAL(rstd),POINTER :: fv(:) ! coriolis (evaluted on a vertex)
85	!$OMP THREADPRIVATE(fv)
86
87
88	INTEGER, PARAMETER :: ne_right=1
89	INTEGER, PARAMETER :: ne_rup=-1
90	INTEGER, PARAMETER :: ne_lup=1
91	INTEGER, PARAMETER :: ne_left=-1
92	INTEGER, PARAMETER :: ne_ldown=1
93	INTEGER, PARAMETER :: ne_rdown=-1
94
95	CONTAINS
96
97	SUBROUTINE allocate_geometry
98	USE field_mod
99	IMPLICIT NONE
100
101	CALL allocate_field(geom%Ai,field_t,type_real,name='Ai')
102
103	CALL allocate_field(geom%xyz_i,field_t,type_real,3)
104	CALL allocate_field(geom%lon_i,field_t,type_real)
105	CALL allocate_field(geom%lat_i,field_t,type_real)
106	CALL allocate_field(geom%elon_i,field_t,type_real,3)
107	CALL allocate_field(geom%elat_i,field_t,type_real,3)
108	CALL allocate_field(geom%centroid,field_t,type_real,3)
109
110	CALL allocate_field(geom%xyz_e,field_u,type_real,3)
111	CALL allocate_field(geom%lon_e,field_u,type_real)
112	CALL allocate_field(geom%lat_e,field_u,type_real)
113	CALL allocate_field(geom%elon_e,field_u,type_real,3)
114	CALL allocate_field(geom%elat_e,field_u,type_real,3)
115	CALL allocate_field(geom%ep_e,field_u,type_real,3)
116	CALL allocate_field(geom%et_e,field_u,type_real,3)
117
118	CALL allocate_field(geom%xyz_v,field_z,type_real,3)
119	CALL allocate_field(geom%de,field_u,type_real)
120	CALL allocate_field(geom%le,field_u,type_real)
121	CALL allocate_field(geom%le_de,field_u,type_real)
122	CALL allocate_field(geom%bi,field_t,type_real)
123	CALL allocate_field(geom%Av,field_z,type_real)
124	CALL allocate_field(geom%Riv,field_t,type_real,6)
125	CALL allocate_field(geom%Riv2,field_t,type_real,6)
126	CALL allocate_field(geom%ne,field_t,type_integer,6)
127	CALL allocate_field(geom%Wee,field_u,type_real,5,2)
128	CALL allocate_field(geom%bi,field_t,type_real)
129	CALL allocate_field(geom%fv,field_z,type_real)
130
131	END SUBROUTINE allocate_geometry
132
133
134	SUBROUTINE swap_geometry(ind)
135	USE field_mod
136	IMPLICIT NONE
137	INTEGER,INTENT(IN) :: ind
138	!!$OMP MASTER
139	Ai=geom%Ai(ind)
140	xyz_i=geom%xyz_i(ind)
141	centroid=geom%centroid(ind)
142	xyz_e=geom%xyz_e(ind)
143	ep_e=geom%ep_e(ind)
144	et_e=geom%et_e(ind)
145	lon_i=geom%lon_i(ind)
146	lat_i=geom%lat_i(ind)
147	lon_e=geom%lon_e(ind)
148	lat_e=geom%lat_e(ind)
149	elon_i=geom%elon_i(ind)
150	elat_i=geom%elat_i(ind)
151	elon_e=geom%elon_e(ind)
152	elat_e=geom%elat_e(ind)
153	xyz_v=geom%xyz_v(ind)
154	de=geom%de(ind)
155	le=geom%le(ind)
156	le_de=geom%le_de(ind)
157	Av=geom%Av(ind)
158	Riv=geom%Riv(ind)
159	Riv2=geom%Riv2(ind)
160	ne=geom%ne(ind)
161	Wee=geom%Wee(ind)
162	bi=geom%bi(ind)
163	fv=geom%fv(ind)
164	!!$OMP END MASTER
165	!!$OMP BARRIER
166	END SUBROUTINE swap_geometry
167
168	SUBROUTINE update_circumcenters
169	USE domain_mod
170	USE dimensions
171	USE spherical_geom_mod
172	USE vector
173	USE transfert_mod
174	USE omp_para
175
176	IMPLICIT NONE
177	REAL(rstd) :: x1(3),x2(3)
178	REAL(rstd) :: vect(3,6)
179	REAL(rstd) :: centr(3)
180	INTEGER :: ind,i,j,n,k
181	TYPE(t_message),SAVE :: message0, message1
182	LOGICAL, SAVE :: first=.TRUE.
183	!$OMP THREADPRIVATE(first)
184
185	IF (first) THEN
186	CALL init_message(geom%xyz_i, req_i0 ,message0)
187	CALL init_message(geom%xyz_i, req_i1 ,message1)
188	first=.FALSE.
189	ENDIF
190
191	CALL transfert_message(geom%xyz_i,message0)
192	CALL transfert_message(geom%xyz_i,message1)
193
194	DO ind=1,ndomain
195	IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
196	CALL swap_dimensions(ind)
197	CALL swap_geometry(ind)
198	DO j=jj_begin,jj_end
199	DO i=ii_begin,ii_end
200	n=(j-1)*iim+i
201	DO k=0,5
202	x1(:) = xyz_i(n+t_pos(k+1),:)
203	x2(:) = xyz_i(n+t_pos(MOD(k+1,6)+1),:)
204	if (norm(x1-x2)<1e-16) x2(:) = xyz_i(n+t_pos(MOD(k+2,6)+1),:)
205	CALL circumcenter(xyz_i(n,:), x1, x2, xyz_v(n+z_pos(k+1),:))
206	ENDDO
207	ENDDO
208	ENDDO
209	ENDDO
210
211	END SUBROUTINE update_circumcenters
212
213	SUBROUTINE remap_schmidt_loc
214	USE spherical_geom_mod
215	USE getin_mod
216	USE omp_para
217	USE domain_mod
218	USE dimensions
219	IMPLICIT NONE
220	INTEGER :: ind,i,j,n
221	REAL(rstd) :: schmidt_factor, schmidt_lon, schmidt_lat
222
223	! Schmidt transform parameters
224	schmidt_factor = 1.
225	CALL getin('schmidt_factor', schmidt_factor)
226	schmidt_factor = schmidt_factor**2.
227
228	schmidt_lon = 0.
229	CALL getin('schmidt_lon', schmidt_lon)
230	schmidt_lon = schmidt_lon * pi/180.
231
232	schmidt_lat = 45.
233	CALL getin('schmidt_lat', schmidt_lat)
234	schmidt_lat = schmidt_lat * pi/180.
235
236	DO ind=1,ndomain
237	IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
238	CALL swap_dimensions(ind)
239	CALL swap_geometry(ind)
240	DO j=jj_begin,jj_end
241	DO i=ii_begin,ii_end
242	n=(j-1)*iim+i
243	CALL schmidt_transform(xyz_i(n,:), schmidt_factor, schmidt_lon, schmidt_lat)
244	ENDDO
245	ENDDO
246	ENDDO
247	END SUBROUTINE remap_schmidt_loc
248
249	SUBROUTINE optimize_geometry
250	USE metric
251	USE spherical_geom_mod
252	USE domain_mod
253	USE dimensions
254	USE transfert_mod
255	USE vector
256	USE getin_mod
257	USE omp_para
258	IMPLICIT NONE
259	INTEGER :: nb_it=0
260	TYPE(t_domain),POINTER :: d
261	INTEGER :: ind,it,i,j,n,k
262	REAL(rstd) :: x1(3),x2(3)
263	REAL(rstd) :: vect(3,6)
264	REAL(rstd) :: centr(3)
265	REAL(rstd) :: sum
266	LOGICAL :: check
267
268
269	CALL getin('optim_it',nb_it)
270
271	DO ind=1,ndomain
272	IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
273	d=>domain(ind)
274	CALL swap_dimensions(ind)
275	CALL swap_geometry(ind)
276	xyz_i(:,1) = 0 ; xyz_i(:,2) = 0 ; xyz_i(:,3) = 1
277
278	DO j=jj_begin,jj_end
279	DO i=ii_begin,ii_end
280	n=(j-1)*iim+i
281	xyz_i(n,:)=d%xyz(:,i,j)
282	ENDDO
283	ENDDO
284	ENDDO
285
286	CALL update_circumcenters
287
288	DO ind=1,ndomain
289	IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master ) CYCLE
290	d=>domain(ind)
291	CALL swap_dimensions(ind)
292	CALL swap_geometry(ind)
293	DO j=jj_begin,jj_end
294	DO i=ii_begin,ii_end
295	n=(j-1)*iim+i
296	DO k=0,5
297	x1(:) = xyz_v(n+z_pos(k+1),:)
298	x2(:) = d%vertex(:,k,i,j)
299	IF (norm(x1-x2)>1e-10) THEN
300	PRINT*,"vertex diff ",ind,i,j,k
301	PRINT*,x1
302	PRINT*,x2
303	ENDIF
304	ENDDO
305	ENDDO
306	ENDDO
307	ENDDO
308
309
310	DO it=1,nb_it
311	IF (MOD(it,100)==0) THEN
312	check=is_master
313	ELSE
314	check=.FALSE.
315	ENDIF
316
317	sum=0
318	DO ind=1,ndomain
319	IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master ) CYCLE
320	CALL swap_dimensions(ind)
321	CALL swap_geometry(ind)
322	DO j=jj_begin,jj_end
323	DO i=ii_begin,ii_end
324	n=(j-1)*iim+i
325	vect(:,1)=xyz_v(n+z_rup,:)
326	vect(:,2)=xyz_v(n+z_up,:)
327	vect(:,3)=xyz_v(n+z_lup,:)
328	vect(:,4)=xyz_v(n+z_ldown,:)
329	vect(:,5)=xyz_v(n+z_down,:)
330	vect(:,6)=xyz_v(n+z_rdown,:)
331	CALL compute_centroid(vect,6,centr)
332	IF (check) THEN
333	sum=MAX(sum,norm(xyz_i(n,:)-centr(:)))
334	ENDIF
335	xyz_i(n,:)=centr(:)
336	ENDDO
337	ENDDO
338
339	ENDDO
340
341	IF (check) PRINT *,"it = ",it," diff centroid circumcenter ",sum
342
343	CALL update_circumcenters
344
345	ENDDO
346
347	END SUBROUTINE optimize_geometry
348
349	SUBROUTINE update_domain
350	! copy position of generators and vertices back into domain(:)%xyz/vertex
351	! so that XIOS/create_header_gen gets the right values
352	USE omp_para
353	USE dimensions
354	USE domain_mod
355	USE transfert_mpi_mod
356
357	INTEGER :: ind,i,j,k,n
358	TYPE(t_domain),POINTER :: d
359	TYPE(t_field),POINTER :: xyz_glo(:), xyz_loc(:), vertex_glo(:), vertex_loc(:)
360	REAL(rstd), POINTER :: xyz(:,:), vertex(:,:)
361
362	CALL allocate_field(xyz_loc, field_t, type_real, 3)
363	CALL allocate_field(vertex_loc, field_z, type_real, 3)
364
365	DO ind=1,ndomain
366	IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master ) CYCLE
367	CALL swap_dimensions(ind)
368	CALL swap_geometry(ind)
369	xyz = xyz_loc(ind)
370	xyz(:,:) = xyz_i(:,:)
371	vertex = vertex_loc(ind)
372	vertex(:,:) = xyz_v(:,:)
373	END DO
374
375	CALL allocate_field_glo(xyz_glo, field_t, type_real, 3)
376	CALL allocate_field_glo(vertex_glo, field_z, type_real, 3)
377	CALL gather_field(xyz_loc, xyz_glo)
378	CALL bcast_field(xyz_glo)
379	CALL gather_field(vertex_loc, vertex_glo)
380	CALL bcast_field(vertex_glo)
381	CALL deallocate_field(vertex_loc)
382	CALL deallocate_field(xyz_loc)
383
384	DO ind=1,ndomain_glo
385	d=>domain_glo(ind)
386	xyz = xyz_glo(ind)
387	vertex = vertex_glo(ind)
388	DO j=d%jj_begin,d%jj_end
389	DO i=d%ii_begin,d%ii_end
390	n=(j-1)*d%iim+i
391	d%xyz(:,i,j)=xyz(n,:)
392	DO k=0,5
393	d%vertex(:,k,i,j) = vertex(n+d%z_pos(k+1),:)
394	END DO
395	END DO
396	END DO
397	END DO
398
399	CALL deallocate_field_glo(vertex_glo)
400	CALL deallocate_field_glo(xyz_glo)
401
402	END SUBROUTINE update_domain
403
404	SUBROUTINE set_geometry
405	USE metric
406	USE vector
407	USE spherical_geom_mod
408	USE domain_mod
409	USE dimensions
410	USE transfert_mod
411	USE getin_mod
412	USE omp_para
413	IMPLICIT NONE
414
415	REAL(rstd) :: surf(6)
416	REAL(rstd) :: surf_v(6)
417	REAL(rstd) :: vect(3,6)
418	REAL(rstd) :: centr(3)
419	REAL(rstd) :: vet(3),vep(3), vertex(3)
420	INTEGER :: ind,i,j,k,n
421	TYPE(t_domain),POINTER :: d
422	REAL(rstd) :: S
423	REAL(rstd) :: w(6)
424	REAL(rstd) :: lon,lat
425	INTEGER :: ii_glo,jj_glo
426	REAL(rstd) :: S1,S2
427
428
429	CALL optimize_geometry
430	CALL remap_schmidt_loc
431	CALL update_circumcenters
432	! copy position of generators and vertices back into domain(:)%xyz/vertex
433	! so that XIOS gets the right values
434	CALL update_domain
435
436	DO ind=1,ndomain
437	IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master ) CYCLE
438	d=>domain(ind)
439	CALL swap_dimensions(ind)
440	CALL swap_geometry(ind)
441	lon_i(:)=0 ; lat_i(:)=0
442	lon_e(:)=0 ; lat_e(:)=0
443	DO j=jj_begin-1,jj_end+1
444	DO i=ii_begin-1,ii_end+1
445	n=(j-1)*iim+i
446
447	DO k=0,5
448	ne(n,k+1)=d%ne(k,i,j)
449	ENDDO
450
451	vect(:,1)=xyz_v(n+z_rup,:)
452	vect(:,2)=xyz_v(n+z_up,:)
453	vect(:,3)=xyz_v(n+z_lup,:)
454	vect(:,4)=xyz_v(n+z_ldown,:)
455	vect(:,5)=xyz_v(n+z_down,:)
456	vect(:,6)=xyz_v(n+z_rdown,:)
457	CALL compute_centroid(vect,6,centr)
458	centroid(n,:)=centr(:)
459
460
461	CALL xyz2lonlat(xyz_v(n+z_up,:),lon,lat)
462	fv(n+z_up)=2sin(lat)omega
463	CALL xyz2lonlat(xyz_v(n+z_down,:),lon,lat)
464	fv(n+z_down)=2sin(lat)omega
465
466	bi(n)=0.
467
468	CALL dist_cart(xyz_i(n,:),xyz_i(n+t_right,:),de(n+u_right))
469	CALL dist_cart(xyz_i(n,:),xyz_i(n+t_lup,:),de(n+u_lup))
470	CALL dist_cart(xyz_i(n,:),xyz_i(n+t_ldown,:),de(n+u_ldown))
471
472	CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_right,:),0.5,xyz_e(n+u_right,:))
473	CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_lup,:),0.5,xyz_e(n+u_lup,:))
474	CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_ldown,:),0.5,xyz_e(n+u_ldown,:))
475
476	CALL dist_cart(xyz_v(n+z_rdown,:), xyz_v(n+z_rup,:),le(n+u_right))
477	CALL dist_cart(xyz_v(n+z_up,:), xyz_v(n+z_lup,:),le(n+u_lup))
478	CALL dist_cart(xyz_v(n+z_ldown,:), xyz_v(n+z_down,:),le(n+u_ldown))
479
480	le_de(n+u_right)=le(n+u_right)/de(n+u_right) ! NaN possible but should be harmless
481	le_de(n+u_lup) =le(n+u_lup) /de(n+u_lup)
482	le_de(n+u_ldown)=le(n+u_ldown)/de(n+u_ldown)
483
484	Ai(n)=0
485	DO k=0,5
486	CALL surf_triangle(xyz_i(n,:),xyz_i(n+t_pos(k+1),:),xyz_i(n+t_pos(MOD((k+1+6),6)+1),:),surf_v(k+1))
487	CALL surf_triangle(xyz_i(n,:),xyz_v(n+z_pos(MOD((k-1+6),6)+1),:),xyz_v(n+z_pos(k+1),:),surf(k+1))
488	Ai(n)=Ai(n)+surf(k+1)
489	IF (i==ii_end .AND. j==jj_begin) THEN
490	IF (Ai(n)<1e20) THEN
491	ELSE
492	PRINT *,"PB !!",Ai(n),k,surf(k+1)
493	PRINT*,xyz_i(n,:),xyz_v(n+z_pos(MOD((k-1+6),6)+1),:),xyz_v(n+z_pos(k+1),:)
494	ENDIF
495	ENDIF
496	ENDDO
497
498	! Sign convention : Ringler et al., JCP 2010, eq. 21 p. 3071
499	! Normal component is along outgoing normal vector if ne=1
500
501	CALL cross_product2(xyz_v(n+z_rdown,:),xyz_v(n+z_rup,:),vep)
502	IF (norm(vep)>1e-30) THEN
503	vep(:)=vep(:)/norm(vep) ! Inward normal vector
504	CALL cross_product2(vep,xyz_e(n+u_right,:),vet) ! Counter-clockwise tangent vector
505	vet(:)=vet(:)/norm(vet)
506	ep_e(n+u_right,:)=-vep(:)*ne(n,right)
507	et_e(n+u_right,:)=vet(:)*ne(n,right)
508	ENDIF
509
510	CALL cross_product2(xyz_v(n+z_up,:),xyz_v(n+z_lup,:),vep)
511	IF (norm(vep)>1e-30) THEN
512	vep(:)=vep(:)/norm(vep)
513	CALL cross_product2(vep,xyz_e(n+u_lup,:),vet)
514	vet(:)=vet(:)/norm(vet)
515	ep_e(n+u_lup,:)=-vep(:)*ne(n,lup)
516	et_e(n+u_lup,:)=vet(:)*ne(n,lup)
517	ENDIF
518
519	CALL cross_product2(xyz_v(n+z_ldown,:),xyz_v(n+z_down,:),vep)
520	IF (norm(vep)>1e-30) THEN
521	vep(:)=vep(:)/norm(vep)
522	CALL cross_product2(vep,xyz_e(n+u_ldown,:),vet)
523	vet(:)=vet(:)/norm(vet)
524	ep_e(n+u_ldown,:)=-vep(:)*ne(n,ldown)
525	et_e(n+u_ldown,:)=vet(:)*ne(n,ldown)
526	ENDIF
527
528	CALL xyz2lonlat(xyz_i(n,:),lon,lat)
529	lon_i(n)=lon
530	lat_i(n)=lat
531	elon_i(n,1) = -sin(lon)
532	elon_i(n,2) = cos(lon)
533	elon_i(n,3) = 0
534	elat_i(n,1) = -cos(lon)*sin(lat)
535	elat_i(n,2) = -sin(lon)*sin(lat)
536	elat_i(n,3) = cos(lat)
537
538
539	CALL xyz2lonlat(xyz_e(n+u_right,:),lon,lat)
540	lon_e(n+u_right)=lon
541	lat_e(n+u_right)=lat
542	elon_e(n+u_right,1) = -sin(lon)
543	elon_e(n+u_right,2) = cos(lon)
544	elon_e(n+u_right,3) = 0
545	elat_e(n+u_right,1) = -cos(lon)*sin(lat)
546	elat_e(n+u_right,2) = -sin(lon)*sin(lat)
547	elat_e(n+u_right,3) = cos(lat)
548
549	CALL xyz2lonlat(xyz_e(n+u_lup,:),lon,lat)
550	lon_e(n+u_lup)=lon
551	lat_e(n+u_lup)=lat
552	elon_e(n+u_lup,1) = -sin(lon)
553	elon_e(n+u_lup,2) = cos(lon)
554	elon_e(n+u_lup,3) = 0
555	elat_e(n+u_lup,1) = -cos(lon)*sin(lat)
556	elat_e(n+u_lup,2) = -sin(lon)*sin(lat)
557	elat_e(n+u_lup,3) = cos(lat)
558
559	CALL xyz2lonlat(xyz_e(n+u_ldown,:),lon,lat)
560	lon_e(n+u_ldown)=lon
561	lat_e(n+u_ldown)=lat
562	elon_e(n+u_ldown,1) = -sin(lon)
563	elon_e(n+u_ldown,2) = cos(lon)
564	elon_e(n+u_ldown,3) = 0
565	elat_e(n+u_ldown,1) = -cos(lon)*sin(lat)
566	elat_e(n+u_ldown,2) = -sin(lon)*sin(lat)
567	elat_e(n+u_ldown,3) = cos(lat)
568
569
570	DO k=0,5
571	CALL surf_triangle(xyz_i(n,:), xyz_v(n+z_pos(k+1),:), xyz_i(n+t_pos(k+1),:),S1)
572	CALL surf_triangle(xyz_i(n,:), xyz_v(n+z_pos(k+1),:), xyz_i(n+t_pos(MOD(k+1+6,6)+1),:),S2)
573	! Riv(n,k+1)=0.5(S1+S2)(radius**2) ! Definition modified for DEC
574	Riv(n,k+1)=0.5*(S1+S2)/Ai(n)
575	Riv2(n,k+1)=0.5*(S1+S2)/surf_v(k+1)
576	ENDDO
577
578	DO k=1,6
579	IF (ABS(surf_v(k))<1e-30) THEN
580	Riv(n,k)=0.
581	ENDIF
582	ENDDO
583
584	Av(n+z_up)=surf_v(vup)+1e-100
585	Av(n+z_down)=surf_v(vdown)+1e-100
586
587	ENDDO
588	ENDDO
589
590	DO j=jj_begin,jj_end
591	DO i=ii_begin,ii_end
592	n=(j-1)*iim+i
593
594	CALL compute_wee(n,right,w)
595	Wee(n+u_right,:,1)=w(1:5)
596
597	CALL compute_wee(n+t_right,left,w)
598	Wee(n+u_right,:,2)=w(1:5)
599
600
601	CALL compute_wee(n,lup,w)
602	Wee(n+u_lup,:,1)=w(1:5)
603
604	CALL compute_wee(n+t_lup,rdown,w)
605	Wee(n+u_lup,:,2)=w(1:5)
606
607
608	CALL compute_wee(n,ldown,w)
609	Wee(n+u_ldown,:,1)=w(1:5)
610
611	CALL compute_wee(n+t_ldown,rup,w)
612	Wee(n+u_ldown,:,2)=w(1:5)
613
614	ENDDO
615	ENDDO
616
617	DO j=jj_begin,jj_end
618	DO i=ii_begin,ii_end
619	n=(j-1)*iim+i
620	ii_glo=d%ii_begin_glo-d%ii_begin+i
621	jj_glo=d%jj_begin_glo-d%jj_begin+j
622
623	IF (ii_glo==1 .AND. jj_glo==1) THEN
624	le(n+u_ldown)=0
625	xyz_v(n+z_ldown,:)=xyz_v(n+z_down,:)
626
627	ENDIF
628
629	IF (ii_glo==iim_glo .AND. jj_glo==1) THEN
630	le(n+u_right)=0
631	xyz_v(n+z_rdown,:)=xyz_v(n+z_rup,:)
632	ENDIF
633
634	IF (ii_glo==iim_glo .AND. jj_glo==jjm_glo) THEN
635	le(n+u_rup)=0
636	xyz_v(n+z_rup,:)=xyz_v(n+z_up,:)
637	ENDIF
638
639	IF (ii_glo==1 .AND. jj_glo==jjm_glo) THEN
640	le(n+u_lup)=0
641	xyz_v(n+z_up,:)=xyz_v(n+z_lup,:)
642	ENDIF
643
644	ENDDO
645	ENDDO
646
647	DO j=jj_begin-1,jj_end+1
648	DO i=ii_begin-1,ii_end+1
649	n=(j-1)*iim+i
650	xyz_i(n,:)=xyz_i(n,:) * radius
651	xyz_v(n+z_up,:)=xyz_v(n+z_up,:) * radius
652	xyz_v(n+z_down,:)=xyz_v(n+z_down,:) *radius
653	de(n+u_right)=de(n+u_right) * radius
654	de(n+u_lup)=de(n+u_lup)*radius
655	de(n+u_ldown)=de(n+u_ldown)*radius
656	xyz_e(n+u_right,:)=xyz_e(n+u_right,:)*radius
657	xyz_e(n+u_lup,:)=xyz_e(n+u_lup,:)*radius
658	xyz_e(n+u_ldown,:)=xyz_e(n+u_ldown,:)*radius
659	le(n+u_right)=le(n+u_right)*radius
660	le(n+u_lup)=le(n+u_lup)*radius
661	le(n+u_ldown)=le(n+u_ldown)*radius
662	Ai(n)=Ai(n)radius*2
663	Av(n+z_up)=Av(n+z_up)radius*2
664	Av(n+z_down)=Av(n+z_down)radius*2
665	ENDDO
666	ENDDO
667
668	ENDDO
669
670	CALL transfert_request(geom%Ai,req_i1)
671	CALL transfert_request(geom%centroid,req_i1)
672
673	! CALL surf_triangle(d%xyz(:,ii_begin,jj_begin),d%xyz(:,ii_begin,jj_end),d%xyz(:,ii_end,jj_begin),S)
674
675	END SUBROUTINE set_geometry
676
677	SUBROUTINE compute_wee(n,pos,w)
678	IMPLICIT NONE
679	INTEGER,INTENT(IN) :: n
680	INTEGER,INTENT(IN) :: pos
681	REAL(rstd),INTENT(OUT) ::w(6)
682
683	REAL(rstd) :: ne_(0:5)
684	REAL(rstd) :: Riv_(6)
685	INTEGER :: k
686
687
688	DO k=0,5
689	ne_(k)=ne(n,MOD(pos-1+k+6,6)+1)
690	Riv_(k+1)=Riv(n,MOD(pos-1+k+6,6)+1)
691	ENDDO
692
693	w(1)=-ne_(0)ne_(1)(Riv_(1)-0.5)
694	w(2)=-ne_(2)(ne_(0)Riv_(2)-w(1)*ne_(1))
695	w(3)=-ne_(3)(ne_(0)Riv_(3)-w(2)*ne_(2))
696	w(4)=-ne_(4)(ne_(0)Riv_(4)-w(3)*ne_(3))
697	w(5)=-ne_(5)(ne_(0)Riv_(5)-w(4)*ne_(4))
698	w(6)=ne_(0)ne_(5)(Riv_(6)-0.5)
699
700	! IF ( ABS(w(5)-w(6))>1e-20) PRINT *, "pb pour wee : w(5)!=w(6)",sum(Riv_(:))
701
702	END SUBROUTINE compute_wee
703
704
705
706	SUBROUTINE compute_geometry
707	IMPLICIT NONE
708	CALL allocate_geometry
709	CALL set_geometry
710
711	END SUBROUTINE compute_geometry
712
713	END MODULE geometry

Note: See TracBrowser for help on using the repository browser.

Download in other formats: