Context Navigation

geometry.f90 @ 1023

Last change on this file since 1023 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

Rev	Line
[12]	1	MODULE geometry
	2	USE field_mod
[356]	3	IMPLICIT NONE
	4
[12]	5	TYPE t_geometry
[286]	6	TYPE(t_field),POINTER :: centroid(:)
[12]	7	TYPE(t_field),POINTER :: xyz_i(:)
	8	TYPE(t_field),POINTER :: xyz_e(:)
	9	TYPE(t_field),POINTER :: xyz_v(:)
[286]	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(:)
[15]	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(:)
[12]	20	TYPE(t_field),POINTER :: Ai(:)
	21	TYPE(t_field),POINTER :: Av(:)
	22	TYPE(t_field),POINTER :: de(:)
	23	TYPE(t_field),POINTER :: le(:)
[362]	24	TYPE(t_field),POINTER :: le_de(:) ! le/de, 0. if de=0.
[12]	25	TYPE(t_field),POINTER :: Riv(:)
[427]	26	TYPE(t_field),POINTER :: S1(:)
	27	TYPE(t_field),POINTER :: S2(:)
[12]	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
[186]	35	TYPE(t_geometry),SAVE,TARGET :: geom
	36
[12]	37
[15]	38	REAL(rstd),POINTER :: Ai(:) ! area of a cell
[186]	39	!$OMP THREADPRIVATE(Ai)
[286]	40	REAL(rstd),POINTER :: centroid(:,:) ! coordinate of the centroid of the cell
	41	!$OMP THREADPRIVATE(centroid)
[15]	42	REAL(rstd),POINTER :: xyz_i(:,:) ! coordinate of the center of the cell (voronoi)
[186]	43	!$OMP THREADPRIVATE(xyz_i)
[15]	44	REAL(rstd),POINTER :: xyz_e(:,:) ! coordinate of a wind point on the cell on a edge
[186]	45	!$OMP THREADPRIVATE(xyz_e)
[286]	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)
[15]	56	REAL(rstd),POINTER :: ep_e(:,:) ! perpendicular unit vector of a edge (outsider)
[186]	57	!$OMP THREADPRIVATE(ep_e)
[15]	58	REAL(rstd),POINTER :: et_e(:,:) ! tangeantial unit vector of a edge
[186]	59	!$OMP THREADPRIVATE(et_e)
[15]	60	REAL(rstd),POINTER :: elon_i(:,:) ! unit longitude vector on the center
[186]	61	!$OMP THREADPRIVATE(elon_i)
[15]	62	REAL(rstd),POINTER :: elat_i(:,:) ! unit latitude vector on the center
[186]	63	!$OMP THREADPRIVATE(elat_i)
[15]	64	REAL(rstd),POINTER :: elon_e(:,:) ! unit longitude vector on a wind point
[186]	65	!$OMP THREADPRIVATE(elon_e)
[15]	66	REAL(rstd),POINTER :: elat_e(:,:) ! unit latitude vector on a wind point
[186]	67	!$OMP THREADPRIVATE(elat_e)
[15]	68	REAL(rstd),POINTER :: Av(:) ! area of dual mesk cell
[186]	69	!$OMP THREADPRIVATE(Av)
[15]	70	REAL(rstd),POINTER :: de(:) ! distance from a neighbour == lenght of an edge of the dual mesh
[186]	71	!$OMP THREADPRIVATE(de)
[15]	72	REAL(rstd),POINTER :: le(:) ! lenght of a edge
[186]	73	!$OMP THREADPRIVATE(le)
[427]	74	REAL(rstd),POINTER :: le_de(:) ! le/de
[362]	75	!$OMP THREADPRIVATE(le_de)
[427]	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)
[15]	80	REAL(rstd),POINTER :: Riv(:,:) ! weight
[186]	81	!$OMP THREADPRIVATE(Riv)
[15]	82	REAL(rstd),POINTER :: Riv2(:,:) ! weight
[186]	83	!$OMP THREADPRIVATE(Riv2)
[15]	84	INTEGER,POINTER :: ne(:,:) ! convention for the way on the normal wind on an edge
[186]	85	!$OMP THREADPRIVATE(ne)
[15]	86	REAL(rstd),POINTER :: Wee(:,:,:) ! weight
[186]	87	!$OMP THREADPRIVATE(Wee)
[15]	88	REAL(rstd),POINTER :: bi(:) ! orographie
[186]	89	!$OMP THREADPRIVATE(bi)
[15]	90	REAL(rstd),POINTER :: fv(:) ! coriolis (evaluted on a vertex)
[186]	91	!$OMP THREADPRIVATE(fv)
[12]	92
	93
[146]	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
[153]	100
[12]	101	CONTAINS
	102
	103	SUBROUTINE allocate_geometry
	104	USE field_mod
	105	IMPLICIT NONE
	106
[266]	107	CALL allocate_field(geom%Ai,field_t,type_real,name='Ai')
[286]	108
[12]	109	CALL allocate_field(geom%xyz_i,field_t,type_real,3)
[286]	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)
[15]	114	CALL allocate_field(geom%centroid,field_t,type_real,3)
[286]	115
[12]	116	CALL allocate_field(geom%xyz_e,field_u,type_real,3)
[286]	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)
[15]	121	CALL allocate_field(geom%ep_e,field_u,type_real,3)
	122	CALL allocate_field(geom%et_e,field_u,type_real,3)
[286]	123
[12]	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)
[362]	127	CALL allocate_field(geom%le_de,field_u,type_real)
[12]	128	CALL allocate_field(geom%bi,field_t,type_real)
	129	CALL allocate_field(geom%Av,field_z,type_real)
[427]	130	CALL allocate_field(geom%S1,field_t,type_real,6)
	131	CALL allocate_field(geom%S2,field_t,type_real,6)
[12]	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
[186]	145	!!$OMP MASTER
[12]	146	Ai=geom%Ai(ind)
	147	xyz_i=geom%xyz_i(ind)
[15]	148	centroid=geom%centroid(ind)
[12]	149	xyz_e=geom%xyz_e(ind)
[15]	150	ep_e=geom%ep_e(ind)
	151	et_e=geom%et_e(ind)
[286]	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)
[15]	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)
[12]	160	xyz_v=geom%xyz_v(ind)
	161	de=geom%de(ind)
	162	le=geom%le(ind)
[362]	163	le_de=geom%le_de(ind)
[12]	164	Av=geom%Av(ind)
[427]	165	S1=geom%S1(ind)
	166	S2=geom%S2(ind)
[12]	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)
[186]	173	!!$OMP END MASTER
	174	!!$OMP BARRIER
[12]	175	END SUBROUTINE swap_geometry
[153]	176
	177	SUBROUTINE update_circumcenters
	178	USE domain_mod
	179	USE dimensions
	180	USE spherical_geom_mod
	181	USE vector
	182	USE transfert_mod
[295]	183	USE omp_para
[153]	184
	185	IMPLICIT NONE
	186	REAL(rstd) :: x1(3),x2(3)
	187	INTEGER :: ind,i,j,n,k
[294]	188	TYPE(t_message),SAVE :: message0, message1
[186]	189	LOGICAL, SAVE :: first=.TRUE.
	190	!$OMP THREADPRIVATE(first)
	191
	192	IF (first) THEN
[294]	193	CALL init_message(geom%xyz_i, req_i0 ,message0)
	194	CALL init_message(geom%xyz_i, req_i1 ,message1)
[186]	195	first=.FALSE.
	196	ENDIF
	197
[294]	198	CALL transfert_message(geom%xyz_i,message0)
	199	CALL transfert_message(geom%xyz_i,message1)
[186]	200
[153]	201	DO ind=1,ndomain
[295]	202	IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
[153]	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
[356]	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
[1014]	229	LOGICAL :: read_metric=.FALSE.
	230
	231	CALL getin('read_metric', read_metric)
[356]	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
[1014]	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
[356]	256	ENDDO
	257	ENDDO
[1014]	258	ENDIF
	259
[356]	260	END SUBROUTINE remap_schmidt_loc
	261
[15]	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
[186]	269	USE getin_mod
[295]	270	USE omp_para
[1014]	271	USE checksum_mod
[15]	272	IMPLICIT NONE
[17]	273	INTEGER :: nb_it=0
[15]	274	TYPE(t_domain),POINTER :: d
[953]	275	INTEGER :: ind,it,i,j,n
[15]	276	REAL(rstd) :: vect(3,6)
	277	REAL(rstd) :: centr(3)
	278	REAL(rstd) :: sum
	279	LOGICAL :: check
[1014]	280	LOGICAL :: read_metric=.FALSE.
[12]	281
[17]	282	CALL getin('optim_it',nb_it)
[1014]	283	CALL getin('read_metric', read_metric)
	284	IF (read_metric) nb_it=0
[17]	285
[15]	286	DO ind=1,ndomain
[295]	287	IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
[15]	288	d=>domain(ind)
	289	CALL swap_dimensions(ind)
	290	CALL swap_geometry(ind)
[388]	291	xyz_i(:,1) = 0 ; xyz_i(:,2) = 0 ; xyz_i(:,3) = 1
	292
[15]	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
[153]	301	CALL update_circumcenters
	302
[882]	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
[15]	324
	325
	326	DO it=1,nb_it
	327	IF (MOD(it,100)==0) THEN
[324]	328	check=is_master
[15]	329	ELSE
	330	check=.FALSE.
	331	ENDIF
	332
	333	sum=0
	334	DO ind=1,ndomain
[295]	335	IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master ) CYCLE
[15]	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
[324]	357	IF (check) PRINT *,"it = ",it," diff centroid circumcenter ",sum
[15]	358
[153]	359	CALL update_circumcenters
[15]	360
	361	ENDDO
	362
	363	END SUBROUTINE optimize_geometry
[153]	364
[356]	365	SUBROUTINE update_domain
	366	! copy position of generators and vertices back into domain(:)%xyz/vertex
[370]	367	! so that XIOS/create_header_gen gets the right values
[356]	368	USE omp_para
	369	USE dimensions
[370]	370	USE domain_mod
[963]	371	USE transfert_mod
[370]	372
[356]	373	INTEGER :: ind,i,j,k,n
	374	TYPE(t_domain),POINTER :: d
[469]	375	TYPE(t_field),POINTER,SAVE :: xyz_glo(:), xyz_loc(:), vertex_glo(:), vertex_loc(:)
[370]	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
[356]	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)
[370]	385	xyz = xyz_loc(ind)
	386	xyz(:,:) = xyz_i(:,:)
	387	vertex = vertex_loc(ind)
	388	vertex(:,:) = xyz_v(:,:)
	389	END DO
[469]	390
	391	!$OMP BARRIER
	392	!$OMP MASTER
[370]	393	CALL allocate_field_glo(xyz_glo, field_t, type_real, 3)
	394	CALL allocate_field_glo(vertex_glo, field_z, type_real, 3)
[469]	395
[963]	396
[370]	397	CALL gather_field(xyz_loc, xyz_glo)
[469]	398	CALL gather_field(vertex_loc, vertex_glo)
	399
[393]	400	CALL bcast_field(xyz_glo)
	401	CALL bcast_field(vertex_glo)
[370]	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,:)
[356]	411	DO k=0,5
[370]	412	d%vertex(:,k,i,j) = vertex(n+d%z_pos(k+1),:)
[356]	413	END DO
	414	END DO
	415	END DO
	416	END DO
[370]	417
	418	CALL deallocate_field_glo(vertex_glo)
	419	CALL deallocate_field_glo(xyz_glo)
[469]	420	!$OMP END MASTER
	421	!$OMP BARRIER
	422	CALL deallocate_field(vertex_loc)
	423	CALL deallocate_field(xyz_loc)
[370]	424
[356]	425	END SUBROUTINE update_domain
	426
[12]	427	SUBROUTINE set_geometry
	428	USE metric
[15]	429	USE vector
[12]	430	USE spherical_geom_mod
	431	USE domain_mod
	432	USE dimensions
[15]	433	USE transfert_mod
[186]	434	USE getin_mod
[295]	435	USE omp_para
[1014]	436	USE checksum_mod
[12]	437	IMPLICIT NONE
	438
	439	REAL(rstd) :: surf(6)
	440	REAL(rstd) :: surf_v(6)
[15]	441	REAL(rstd) :: vect(3,6)
[21]	442	REAL(rstd) :: centr(3)
[899]	443	REAL(rstd) :: vet(3),vep(3)
[12]	444	INTEGER :: ind,i,j,k,n
	445	TYPE(t_domain),POINTER :: d
[427]	446	REAL(rstd) :: S12
[12]	447	REAL(rstd) :: w(6)
[15]	448	REAL(rstd) :: lon,lat
[12]	449	INTEGER :: ii_glo,jj_glo
[15]	450
[12]	451
[15]	452	CALL optimize_geometry
[356]	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
[1014]	458
[153]	459	DO ind=1,ndomain
[295]	460	IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master ) CYCLE
[12]	461	d=>domain(ind)
	462	CALL swap_dimensions(ind)
	463	CALL swap_geometry(ind)
[388]	464	lon_i(:)=0 ; lat_i(:)=0
	465	lon_e(:)=0 ; lat_e(:)=0
[12]	466	DO j=jj_begin-1,jj_end+1
	467	DO i=ii_begin-1,ii_end+1
	468	n=(j-1)*iim+i
[15]	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
[12]	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
[15]	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))
[12]	494
[15]	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))
[362]	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
[12]	507	Ai(n)=0
	508	DO k=0,5
[15]	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))
[12]	511	Ai(n)=Ai(n)+surf(k+1)
[15]	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
[12]	519	ENDDO
	520
[47]	521	! Sign convention : Ringler et al., JCP 2010, eq. 21 p. 3071
	522	! Normal component is along outgoing normal vector if ne=1
[12]	523
[15]	524	CALL cross_product2(xyz_v(n+z_rdown,:),xyz_v(n+z_rup,:),vep)
	525	IF (norm(vep)>1e-30) THEN
[47]	526	vep(:)=vep(:)/norm(vep) ! Inward normal vector
	527	CALL cross_product2(vep,xyz_e(n+u_right,:),vet) ! Counter-clockwise tangent vector
[15]	528	vet(:)=vet(:)/norm(vet)
[47]	529	ep_e(n+u_right,:)=-vep(:)*ne(n,right)
[15]	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)
[47]	538	ep_e(n+u_lup,:)=-vep(:)*ne(n,lup)
[15]	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)
[47]	547	ep_e(n+u_ldown,:)=-vep(:)*ne(n,ldown)
[15]	548	et_e(n+u_ldown,:)=vet(:)*ne(n,ldown)
	549	ENDIF
	550
	551	CALL xyz2lonlat(xyz_i(n,:),lon,lat)
[286]	552	lon_i(n)=lon
	553	lat_i(n)=lat
[15]	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)
[286]	563	lon_e(n+u_right)=lon
	564	lat_e(n+u_right)=lat
[15]	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)
[286]	573	lon_e(n+u_lup)=lon
	574	lat_e(n+u_lup)=lat
[15]	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)
[286]	583	lon_e(n+u_ldown)=lon
	584	lat_e(n+u_ldown)=lat
[15]	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
[12]	593	DO k=0,5
[427]	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)
[12]	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
[604]	648	le_de(n+u_ldown)=0
[12]	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
[604]	655	le_de(n+u_right)=0
[12]	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
[604]	661	le_de(n+u_rup)=0
[12]	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
[604]	667	le_de(n+u_lup)=0
[12]	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
[186]	696
[15]	697	CALL transfert_request(geom%Ai,req_i1)
	698	CALL transfert_request(geom%centroid,req_i1)
[186]	699
[295]	700	! CALL surf_triangle(d%xyz(:,ii_begin,jj_begin),d%xyz(:,ii_begin,jj_end),d%xyz(:,ii_end,jj_begin),S)
[151]	701
[953]	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
[12]	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 @ 1023

Download in other formats: