Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

geometry.f90

Last change on this file was 1014, checked in by ymipsl, 4 years ago

1j+1j=2j : when restart from a file, in order to do not modify metric, it must

do not optimize the metric even if the optim_it getin parameter is set to a value > 0
do not call the schimtt tranform, even if the zoom procedure return identity, due to rounding error (probably due to finite decimal of pi)

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

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

Context Navigation

source: codes/icosagcm/trunk/src/sphere/geometry.f90

Download in other formats: