Context Navigation

source: codes/icosagcm/trunk/src/advect.f90 @ 238

Last change on this file since 238 was 186, checked in by ymipsl, 10 years ago

Add new openMP parallelism based on distribution of domains on threads. There is no more limitation of number of threads by MPI process.

YM

File size: 21.9 KB

Rev	Line
[22]	1	MODULE advect_mod
	2	USE icosa
	3	IMPLICIT NONE
[17]	4
[138]	5	PRIVATE
	6
	7	PUBLIC :: init_advect, compute_backward_traj, compute_gradq3d, compute_advect_horiz
	8
[17]	9	CONTAINS
	10
[22]	11	!==========================================================================
[17]	12
[148]	13	SUBROUTINE init_advect(normal,tangent,one_over_sqrt_leng)
[22]	14	IMPLICIT NONE
	15	REAL(rstd),INTENT(OUT) :: normal(3iimjjm,3)
	16	REAL(rstd),INTENT(OUT) :: tangent(3iimjjm,3)
[148]	17	REAL(rstd),INTENT(OUT) :: one_over_sqrt_leng(iim*jjm)
[17]	18
[174]	19	INTEGER :: ij
[22]	20
[174]	21	!$SIMD
	22	DO ij=ij_begin,ij_end
[22]	23
[174]	24	CALL cross_product2(xyz_v(ij+z_rdown,:),xyz_v(ij+z_rup,:),normal(ij+u_right,:))
	25	normal(ij+u_right,:)=normal(ij+u_right,:)/sqrt(sum(normal(ij+u_right,:)*2)+1e-50)ne(ij,right)
[22]	26
[174]	27	CALL cross_product2(xyz_v(ij+z_up,:),xyz_v(ij+z_lup,:),normal(ij+u_lup,:))
	28	normal(ij+u_lup,:)=normal(ij+u_lup,:)/sqrt(sum(normal(ij+u_lup,:)*2)+1e-50)ne(ij,lup)
[22]	29
[174]	30	CALL cross_product2(xyz_v(ij+z_ldown,:),xyz_v(ij+z_down,:),normal(ij+u_ldown,:))
	31	normal(ij+u_ldown,:)=normal(ij+u_ldown,:)/sqrt(sum(normal(ij+u_ldown,:)*2)+1e-50)ne(ij,ldown)
[22]	32
[174]	33	tangent(ij+u_right,:)=xyz_v(ij+z_rup,:)-xyz_v(ij+z_rdown,:)
	34	tangent(ij+u_right,:)=tangent(ij+u_right,:)/sqrt(sum(tangent(ij+u_right,:)**2)+1e-50)
[22]	35
[174]	36	tangent(ij+u_lup,:)=xyz_v(ij+z_lup,:)-xyz_v(ij+z_up,:)
	37	tangent(ij+u_lup,:)=tangent(ij+u_lup,:)/sqrt(sum(tangent(ij+u_lup,:)**2)+1e-50)
[22]	38
[174]	39	tangent(ij+u_ldown,:)=xyz_v(ij+z_down,:)-xyz_v(ij+z_ldown,:)
	40	tangent(ij+u_ldown,:)=tangent(ij+u_ldown,:)/sqrt(sum(tangent(ij+u_ldown,:)**2)+1e-50)
[148]	41
[174]	42	one_over_sqrt_leng(ij) = 1./sqrt(max(sum((xyz_v(ij+z_up,:) - xyz_i(ij,:))2),sum((xyz_v(ij+z_down,:) - xyz_i(ij,:))2), &
	43	sum((xyz_v(ij+z_rup,:) - xyz_i(ij,:))2),sum((xyz_v(ij+z_rdown,:) - xyz_i(ij,:))2), &
	44	sum((xyz_v(ij+z_lup,:) - xyz_i(ij,:))2),sum((xyz_v(ij+z_ldown,:) - xyz_i(ij,:))2)) )
[22]	45	ENDDO
	46
[17]	47	END SUBROUTINE init_advect
	48
[22]	49	!=======================================================================================
[17]	50
[186]	51	SUBROUTINE compute_gradq3d(qi_in,one_over_sqrt_leng_in,gradq3d_out,xyz_i,xyz_v)
[148]	52	USE trace
[151]	53	USE omp_para
[22]	54	IMPLICIT NONE
[186]	55	REAL(rstd),INTENT(IN) :: qi_in(iim*jjm,llm)
	56	REAL(rstd),INTENT(IN) :: one_over_sqrt_leng_in(iim*jjm)
	57	REAL(rstd),INTENT(IN) :: xyz_i(iim*jjm,3)
	58	REAL(rstd),INTENT(IN) :: xyz_v(2iimjjm,3)
	59	REAL(rstd),INTENT(OUT) :: gradq3d_out(iim*jjm,llm,3)
[17]	60	REAL(rstd) :: maxq,minq,minq_c,maxq_c
[148]	61	REAL(rstd) :: alphamx,alphami,alpha ,maggrd
[17]	62	REAL(rstd) :: leng1,leng2
	63	REAL(rstd) :: arr(2iimjjm)
[148]	64	REAL(rstd) :: ar(iim*jjm)
[17]	65	REAL(rstd) :: gradtri(2iimjjm,llm,3)
[174]	66	INTEGER :: ij,k,ind,l
[186]	67	REAL(rstd) :: qi(iim*jjm,llm)
	68	REAL(rstd) :: one_over_sqrt_leng(iim*jjm)
	69	REAL(rstd) :: gradq3d(iim*jjm,llm,3)
	70	REAL(rstd) :: detx,dety,detz,det
	71	REAL(rstd) :: A(3,3), a11,a12,a13,a21,a22,a23,a31,a32,a33
	72	REAL(rstd) :: x1,x2,x3
	73	REAL(rstd) :: dq(3)
[138]	74
[186]	75	qi=qi_in
	76	one_over_sqrt_leng=one_over_sqrt_leng_in
	77
	78	CALL trace_start("compute_gradq3d1")
[148]	79
[138]	80	! TODO : precompute ar, drop arr as output argument of gradq ?
	81
[22]	82	!========================================================================================== GRADIENT
[138]	83	! Compute gradient at triangles solving a linear system
	84	! arr = area of triangle joining centroids of hexagons
[186]	85	! DO l = ll_begin,ll_end
	86	!!$SIMD
	87	! DO ij=ij_begin_ext,ij_end_ext
	88	!! CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,:),arr(ij+z_up))
	89	!! CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,:),arr(ij+z_down))
	90	! CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,1),gradtri(ij+z_up,l,2),gradtri(ij+z_up,l,3),arr(ij+z_up))
	91	! CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,1),gradtri(ij+z_down,l,2),gradtri(ij+z_down,l,3),arr(ij+z_down))
	92	! END DO
	93	! END DO
	94
[151]	95	DO l = ll_begin,ll_end
[174]	96	!$SIMD
	97	DO ij=ij_begin_ext,ij_end_ext
[186]	98	! CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,1),gradtri(ij+z_up,l,2),gradtri(ij+z_up,l,3),arr(ij+z_up))
	99
	100
	101	A(1,1)=xyz_i(ij+t_rup,1)-xyz_i(ij,1); A(1,2)=xyz_i(ij+t_rup,2)-xyz_i(ij,2); A(1,3)=xyz_i(ij+t_rup,3)-xyz_i(ij,3)
	102	A(2,1)=xyz_i(ij+t_lup,1)-xyz_i(ij,1); A(2,2)=xyz_i(ij+t_lup,2)-xyz_i(ij,2); A(2,3)=xyz_i(ij+t_lup,3)-xyz_i(ij,3)
	103	A(3,1)=xyz_v(ij+z_up,1); A(3,2)= xyz_v(ij+z_up,2); A(3,3)=xyz_v(ij+z_up,3)
	104
	105	dq(1) = qi(ij+t_rup,l)-qi(ij,l)
	106	dq(2) = qi(ij+t_lup,l)-qi(ij,l)
	107	dq(3) = 0.0
	108
	109
	110	! CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),det)
	111
	112	a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
	113	a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
	114	a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
	115
	116	x1 = a11 * (a22 * a33 - a23 * a32)
	117	x2 = a12 * (a21 * a33 - a23 * a31)
	118	x3 = a13 * (a21 * a32 - a22 * a31)
	119	det = x1 - x2 + x3
	120
	121	! CALL determinant(dq(1),dq(2),dq(3),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),detx)
	122
	123	a11=dq(1) ; a12=dq(2) ; a13=dq(3)
	124	a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
	125	a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
	126
	127	x1 = a11 * (a22 * a33 - a23 * a32)
	128	x2 = a12 * (a21 * a33 - a23 * a31)
	129	x3 = a13 * (a21 * a32 - a22 * a31)
	130	detx = x1 - x2 + x3
	131
	132	! CALL determinant(A(1,1),A(2,1),A(3,1),dq(1),dq(2),dq(3),A(1,3),A(2,3),A(3,3),dety)
	133
	134	a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
	135	a21=dq(1) ; a22=dq(2) ; a23=dq(3)
	136	a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
	137
	138	x1 = a11 * (a22 * a33 - a23 * a32)
	139	x2 = a12 * (a21 * a33 - a23 * a31)
	140	x3 = a13 * (a21 * a32 - a22 * a31)
	141	dety = x1 - x2 + x3
	142
	143	! CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),dq(1),dq(2),dq(3),detz)
	144
	145	a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
	146	a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
	147	a31=dq(1) ; a32=dq(2) ; a33=dq(3)
	148
	149	x1 = a11 * (a22 * a33 - a23 * a32)
	150	x2 = a12 * (a21 * a33 - a23 * a31)
	151	x3 = a13 * (a21 * a32 - a22 * a31)
	152	detz = x1 - x2 + x3
	153
	154	gradtri(ij+z_up,l,1) = detx
	155	gradtri(ij+z_up,l,2) = dety
	156	gradtri(ij+z_up,l,3) = detz
	157	arr(ij+z_up) = det
	158
	159	ENDDO
	160
	161	DO ij=ij_begin_ext,ij_end_ext
	162
	163
	164	! CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,1),gradtri(ij+z_down,l,2),gradtri(ij+z_down,l,3),arr(ij+z_down))
	165
	166	A(1,1)=xyz_i(ij+t_ldown,1)-xyz_i(ij,1); A(1,2)=xyz_i(ij+t_ldown,2)-xyz_i(ij,2); A(1,3)=xyz_i(ij+t_ldown,3)-xyz_i(ij,3)
	167	A(2,1)=xyz_i(ij+t_rdown,1)-xyz_i(ij,1); A(2,2)=xyz_i(ij+t_rdown,2)-xyz_i(ij,2); A(2,3)=xyz_i(ij+t_rdown,3)-xyz_i(ij,3)
	168	A(3,1)=xyz_v(ij+z_down,1); A(3,2)= xyz_v(ij+z_down,2); A(3,3)=xyz_v(ij+z_down,3)
	169
	170	dq(1) = qi(ij+t_ldown,l)-qi(ij,l)
	171	dq(2) = qi(ij+t_rdown,l)-qi(ij,l)
	172	dq(3) = 0.0
	173
	174
	175	! CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),det)
	176
	177	a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
	178	a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
	179	a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
	180
	181	x1 = a11 * (a22 * a33 - a23 * a32)
	182	x2 = a12 * (a21 * a33 - a23 * a31)
	183	x3 = a13 * (a21 * a32 - a22 * a31)
	184	det = x1 - x2 + x3
	185
	186	! CALL determinant(dq(1),dq(2),dq(3),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),detx)
	187
	188	a11=dq(1) ; a12=dq(2) ; a13=dq(3)
	189	a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
	190	a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
	191
	192	x1 = a11 * (a22 * a33 - a23 * a32)
	193	x2 = a12 * (a21 * a33 - a23 * a31)
	194	x3 = a13 * (a21 * a32 - a22 * a31)
	195	detx = x1 - x2 + x3
	196
	197	! CALL determinant(A(1,1),A(2,1),A(3,1),dq(1),dq(2),dq(3),A(1,3),A(2,3),A(3,3),dety)
	198
	199	a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
	200	a21=dq(1) ; a22=dq(2) ; a23=dq(3)
	201	a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
	202
	203	x1 = a11 * (a22 * a33 - a23 * a32)
	204	x2 = a12 * (a21 * a33 - a23 * a31)
	205	x3 = a13 * (a21 * a32 - a22 * a31)
	206	dety = x1 - x2 + x3
	207
	208	! CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),dq(1),dq(2),dq(3),detz)
	209
	210	a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
	211	a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
	212	a31=dq(1) ; a32=dq(2) ; a33=dq(3)
	213
	214	x1 = a11 * (a22 * a33 - a23 * a32)
	215	x2 = a12 * (a21 * a33 - a23 * a31)
	216	x3 = a13 * (a21 * a32 - a22 * a31)
	217	detz = x1 - x2 + x3
	218
	219	gradtri(ij+z_down,l,1) = detx
	220	gradtri(ij+z_down,l,2) = dety
	221	gradtri(ij+z_down,l,3) = detz
	222	arr(ij+z_down) = det
	223
	224	END DO
[22]	225	END DO
[17]	226
[174]	227	!$SIMD
[186]	228	DO ij=ij_begin,ij_end
	229	ar(ij) = arr(ij+z_up)+arr(ij+z_lup)+arr(ij+z_ldown)+arr(ij+z_down)+arr(ij+z_rdown)+arr(ij+z_rup)+1.e-50
[148]	230	ENDDO
[186]	231	CALL trace_end("compute_gradq3d1")
	232	CALL trace_start2("compute_gradq3d2")
[148]	233
	234	DO k=1,3
[151]	235	DO l =ll_begin,ll_end
[174]	236	!$SIMD
	237	DO ij=ij_begin,ij_end
	238	gradq3d(ij,l,k) = ( gradtri(ij+z_up,l,k) + gradtri(ij+z_down,l,k) + &
	239	gradtri(ij+z_rup,l,k) + gradtri(ij+z_ldown,l,k) + &
	240	gradtri(ij+z_lup,l,k)+ gradtri(ij+z_rdown,l,k) ) / ar(ij)
[148]	241	END DO
	242	END DO
	243	ENDDO
[186]	244	CALL trace_end2("compute_gradq3d2")
	245	CALL trace_start("compute_gradq3d3")
[22]	246
	247	!============================================================================================= LIMITING
[151]	248	DO l =ll_begin,ll_end
[174]	249	!$SIMD
	250	DO ij=ij_begin,ij_end
[186]	251	! maggrd = dot_product(gradq3d(ij,l,:),gradq3d(ij,l,:))
	252	maggrd = gradq3d(ij,l,1)gradq3d(ij,l,1) + gradq3d(ij,l,2)gradq3d(ij,l,2) + gradq3d(ij,l,3)*gradq3d(ij,l,3)
[22]	253	maggrd = sqrt(maggrd)
[174]	254	maxq_c = qi(ij,l) + maggrd*one_over_sqrt_leng(ij)
	255	minq_c = qi(ij,l) - maggrd*one_over_sqrt_leng(ij)
	256	maxq = max(qi(ij,l),qi(ij+t_right,l),qi(ij+t_lup,l),qi(ij+t_rup,l),qi(ij+t_left,l), &
	257	qi(ij+t_rdown,l),qi(ij+t_ldown,l))
	258	minq = min(qi(ij,l),qi(ij+t_right,l),qi(ij+t_lup,l),qi(ij+t_rup,l),qi(ij+t_left,l), &
	259	qi(ij+t_rdown,l),qi(ij+t_ldown,l))
	260	alphamx = (maxq - qi(ij,l)) ; alphamx = alphamx/(maxq_c - qi(ij,l) )
[22]	261	alphamx = max(alphamx,0.0)
[174]	262	alphami = (minq - qi(ij,l)); alphami = alphami/(minq_c - qi(ij,l))
[22]	263	alphami = max(alphami,0.0)
	264	alpha = min(alphamx,alphami,1.0)
[186]	265	! gradq3d(ij,l,:) = alpha*gradq3d(ij,l,:)
	266	gradq3d(ij,l,1) = alpha*gradq3d(ij,l,1)
	267	gradq3d(ij,l,2) = alpha*gradq3d(ij,l,2)
	268	gradq3d(ij,l,3) = alpha*gradq3d(ij,l,3)
[22]	269	END DO
[17]	270	END DO
[148]	271
[186]	272	CALL trace_end("compute_gradq3d3")
	273
	274	gradq3d_out=gradq3d
	275
	276	CONTAINS
	277
	278	SUBROUTINE gradq(n0,l,n1,n2,n3,q,dq1,dq2,dq3,det)
	279	IMPLICIT NONE
	280	INTEGER, INTENT(IN) :: n0,l,n1,n2,n3
	281	REAL(rstd), INTENT(IN) :: q(iim*jjm,llm)
	282	! REAL(rstd), INTENT(OUT) :: dq(3), det
	283	REAL(rstd), INTENT(OUT) :: dq1,dq2,dq3,det
	284	REAL(rstd) :: dq(3)
	285
	286	REAL(rstd) ::detx,dety,detz
	287	INTEGER :: info
	288	INTEGER :: IPIV(3)
	289
	290	REAL(rstd) :: A(3,3)
	291	REAL(rstd) :: B(3)
	292
	293	! TODO : replace A by A1,A2,A3
	294	A(1,1)=xyz_i(n1,1)-xyz_i(n0,1); A(1,2)=xyz_i(n1,2)-xyz_i(n0,2); A(1,3)=xyz_i(n1,3)-xyz_i(n0,3)
	295	A(2,1)=xyz_i(n2,1)-xyz_i(n0,1); A(2,2)=xyz_i(n2,2)-xyz_i(n0,2); A(2,3)=xyz_i(n2,3)-xyz_i(n0,3)
	296	A(3,1)=xyz_v(n3,1); A(3,2)= xyz_v(n3,2); A(3,3)=xyz_v(n3,3)
	297
	298	dq(1) = q(n1,l)-q(n0,l)
	299	dq(2) = q(n2,l)-q(n0,l)
	300	dq(3) = 0.0
	301
	302	! CALL DGESV(3,1,A,3,IPIV,dq(:),3,info)
	303	! CALL determinant(A(:,1),A(:,2),A(:,3),det)
	304	! CALL determinant(dq,A(:,2),A(:,3),detx)
	305	! CALL determinant(A(:,1),dq,A(:,3),dety)
	306	! CALL determinant(A(:,1),A(:,2),dq,detz)
	307	! dq(1) = detx
	308	! dq(2) = dety
	309	! dq(3) = detz
	310
	311	CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),det)
	312	CALL determinant(dq(1),dq(2),dq(3),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),dq1)
	313	CALL determinant(A(1,1),A(2,1),A(3,1),dq(1),dq(2),dq(3),A(1,3),A(2,3),A(3,3),dq2)
	314	CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),dq(1),dq(2),dq(3),dq3)
	315
	316	END SUBROUTINE gradq
	317
	318	!==========================================================================
	319	! PURE SUBROUTINE determinant(a1,a2,a3,det)
	320	! IMPLICIT NONE
	321	! REAL(rstd), DIMENSION(3), INTENT(IN) :: a1,a2,a3
	322	! REAL(rstd), INTENT(OUT) :: det
	323	! REAL(rstd) :: x1,x2,x3
	324	! x1 = a1(1) * (a2(2) * a3(3) - a2(3) * a3(2))
	325	! x2 = a1(2) * (a2(1) * a3(3) - a2(3) * a3(1))
	326	! x3 = a1(3) * (a2(1) * a3(2) - a2(2) * a3(1))
	327	! det = x1 - x2 + x3
	328	! END SUBROUTINE determinant
	329
	330	SUBROUTINE determinant(a11,a12,a13,a21,a22,a23,a31,a32,a33,det)
	331	IMPLICIT NONE
	332	REAL(rstd), INTENT(IN) :: a11,a12,a13,a21,a22,a23,a31,a32,a33
	333	REAL(rstd), INTENT(OUT) :: det
	334	REAL(rstd) :: x1,x2,x3
	335	x1 = a11 * (a22 * a33 - a23 * a32)
	336	x2 = a12 * (a21 * a33 - a23 * a31)
	337	x3 = a13 * (a21 * a32 - a22 * a31)
	338	det = x1 - x2 + x3
	339	END SUBROUTINE determinant
	340
	341
[22]	342	END SUBROUTINE compute_gradq3d
	343
[138]	344	! Backward trajectories, for use with Miura approach
[137]	345	SUBROUTINE compute_backward_traj(normal,tangent,ue,tau, cc)
[148]	346	USE trace
[151]	347	USE omp_para
[22]	348	IMPLICIT NONE
[137]	349	REAL(rstd),INTENT(IN) :: normal(3iimjjm,3)
	350	REAL(rstd),INTENT(IN) :: tangent(3iimjjm,3)
	351	REAL(rstd),INTENT(IN) :: ue(iim3jjm,llm)
	352	REAL(rstd),INTENT(OUT) :: cc(3iimjjm,llm,3) ! start of backward trajectory
	353	REAL(rstd),INTENT(IN) :: tau
	354
[138]	355	REAL(rstd) :: v_e(3), up_e, qe, ed(3)
[174]	356	INTEGER :: ij,l
[137]	357
[148]	358	CALL trace_start("compute_backward_traj")
	359
[138]	360	! TODO : compute normal displacement ue*tau as hfluxt / mass(upwind) then reconstruct tangential displacement
	361
[137]	362	! reconstruct tangential wind then 3D wind at edge then cc = edge midpoint - u*tau
[151]	363	DO l = ll_begin,ll_end
[174]	364	!$SIMD
	365	DO ij=ij_begin,ij_end
	366	up_e =1/de(ij+u_right)*( &
	367	wee(ij+u_right,1,1)le(ij+u_rup)ue(ij+u_rup,l)+ &
	368	wee(ij+u_right,2,1)le(ij+u_lup)ue(ij+u_lup,l)+ &
	369	wee(ij+u_right,3,1)le(ij+u_left)ue(ij+u_left,l)+ &
	370	wee(ij+u_right,4,1)le(ij+u_ldown)ue(ij+u_ldown,l)+ &
	371	wee(ij+u_right,5,1)le(ij+u_rdown)ue(ij+u_rdown,l)+ &
	372	wee(ij+u_right,1,2)le(ij+t_right+u_ldown)ue(ij+t_right+u_ldown,l)+ &
	373	wee(ij+u_right,2,2)le(ij+t_right+u_rdown)ue(ij+t_right+u_rdown,l)+ &
	374	wee(ij+u_right,3,2)le(ij+t_right+u_right)ue(ij+t_right+u_right,l)+ &
	375	wee(ij+u_right,4,2)le(ij+t_right+u_rup)ue(ij+t_right+u_rup,l)+ &
	376	wee(ij+u_right,5,2)le(ij+t_right+u_lup)ue(ij+t_right+u_lup,l) &
[137]	377	)
[174]	378	v_e = ue(ij+u_right,l)normal(ij+u_right,:) + up_etangent(ij+u_right,:)
	379	cc(ij+u_right,l,:) = xyz_e(ij+u_right,:) - v_e*tau
[137]	380
[174]	381	up_e=1/de(ij+u_lup)*( &
	382	wee(ij+u_lup,1,1)le(ij+u_left)ue(ij+u_left,l)+ &
	383	wee(ij+u_lup,2,1)le(ij+u_ldown)ue(ij+u_ldown,l)+ &
	384	wee(ij+u_lup,3,1)le(ij+u_rdown)ue(ij+u_rdown,l)+ &
	385	wee(ij+u_lup,4,1)le(ij+u_right)ue(ij+u_right,l)+ &
	386	wee(ij+u_lup,5,1)le(ij+u_rup)ue(ij+u_rup,l)+ &
	387	wee(ij+u_lup,1,2)le(ij+t_lup+u_right)ue(ij+t_lup+u_right,l)+ &
	388	wee(ij+u_lup,2,2)le(ij+t_lup+u_rup)ue(ij+t_lup+u_rup,l)+ &
	389	wee(ij+u_lup,3,2)le(ij+t_lup+u_lup)ue(ij+t_lup+u_lup,l)+ &
	390	wee(ij+u_lup,4,2)le(ij+t_lup+u_left)ue(ij+t_lup+u_left,l)+ &
	391	wee(ij+u_lup,5,2)le(ij+t_lup+u_ldown)ue(ij+t_lup+u_ldown,l) &
[137]	392	)
[174]	393	v_e = ue(ij+u_lup,l)normal(ij+u_lup,:) + up_etangent(ij+u_lup,:)
	394	cc(ij+u_lup,l,:) = xyz_e(ij+u_lup,:) - v_e*tau
	395
[137]	396
[174]	397	up_e=1/de(ij+u_ldown)*( &
	398	wee(ij+u_ldown,1,1)le(ij+u_rdown)ue(ij+u_rdown,l)+ &
	399	wee(ij+u_ldown,2,1)le(ij+u_right)ue(ij+u_right,l)+ &
	400	wee(ij+u_ldown,3,1)le(ij+u_rup)ue(ij+u_rup,l)+ &
	401	wee(ij+u_ldown,4,1)le(ij+u_lup)ue(ij+u_lup,l)+ &
	402	wee(ij+u_ldown,5,1)le(ij+u_left)ue(ij+u_left,l)+ &
	403	wee(ij+u_ldown,1,2)le(ij+t_ldown+u_lup)ue(ij+t_ldown+u_lup,l)+ &
	404	wee(ij+u_ldown,2,2)le(ij+t_ldown+u_left)ue(ij+t_ldown+u_left,l)+ &
	405	wee(ij+u_ldown,3,2)le(ij+t_ldown+u_ldown)ue(ij+t_ldown+u_ldown,l)+ &
	406	wee(ij+u_ldown,4,2)le(ij+t_ldown+u_rdown)ue(ij+t_ldown+u_rdown,l)+ &
	407	wee(ij+u_ldown,5,2)le(ij+t_ldown+u_right)ue(ij+t_ldown+u_right,l) &
[137]	408	)
[174]	409	v_e = ue(ij+u_ldown,l)normal(ij+u_ldown,:) + up_etangent(ij+u_ldown,:)
	410	cc(ij+u_ldown,l,:) = xyz_e(ij+u_ldown,:) - v_e*tau
[137]	411	ENDDO
	412	END DO
[148]	413
	414	CALL trace_end("compute_backward_traj")
	415
[137]	416	END SUBROUTINE compute_backward_traj
	417
	418	! Horizontal transport (S. Dubey, T. Dubos)
	419	! Slope-limited van Leer approach with hexagons
[138]	420	SUBROUTINE compute_advect_horiz(update_mass,hfluxt,cc,gradq3d, mass,qi)
[148]	421	USE trace
[151]	422	USE omp_para
[137]	423	IMPLICIT NONE
[138]	424	LOGICAL, INTENT(IN) :: update_mass
	425	REAL(rstd), INTENT(IN) :: gradq3d(iim*jjm,llm,3)
	426	REAL(rstd), INTENT(IN) :: hfluxt(3iimjjm,llm) ! mass flux
	427	REAL(rstd), INTENT(IN) :: cc(3iimjjm,llm,3) ! barycenter of quadrilateral, where q is evaluated (1-point quadrature)
	428	REAL(rstd), INTENT(INOUT) :: mass(iim*jjm,llm)
	429	REAL(rstd), INTENT(INOUT) :: qi(iim*jjm,llm)
[22]	430
[138]	431	REAL(rstd) :: dq,dmass,qe,ed(3), newmass
	432	REAL(rstd) :: qflux(3iimjjm,llm)
[174]	433	INTEGER :: ij,k,l
[17]	434
[148]	435	CALL trace_start("compute_advect_horiz")
	436
[138]	437	! evaluate tracer field at point cc using piecewise linear reconstruction
[136]	438	! q(cc)= q0 + gradq.(cc-xyz_i) with xi centroid of hexagon
[138]	439	! ne*hfluxt>0 iff outgoing
[151]	440	DO l = ll_begin,ll_end
[174]	441
	442	!$SIMD
	443	DO ij=ij_begin_ext,ij_end_ext
	444
	445	IF (ne(ij,right)*hfluxt(ij+u_right,l)>0) THEN
	446	ed = cc(ij+u_right,l,:) - xyz_i(ij,:)
[186]	447	! qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed)
	448	qe = qi(ij,l)+gradq3d(ij,l,1)ed(1)+gradq3d(ij,l,2)ed(2)+gradq3d(ij,l,3)*ed(3)
[174]	449	ELSE
	450	ed = cc(ij+u_right,l,:) - xyz_i(ij+t_right,:)
[186]	451	! qe = qi(ij+t_right,l)+sum2(gradq3d(ij+t_right,l,:),ed)
	452	qe = qi(ij+t_right,l) + gradq3d(ij+t_right,l,1)ed(1)+gradq3d(ij+t_right,l,2)ed(2)+gradq3d(ij+t_right,l,3)*ed(3)
[174]	453	ENDIF
	454	qflux(ij+u_right,l) = hfluxt(ij+u_right,l)*qe
[137]	455
[174]	456	IF (ne(ij,lup)*hfluxt(ij+u_lup,l)>0) THEN
	457	ed = cc(ij+u_lup,l,:) - xyz_i(ij,:)
[186]	458	! qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed)
	459	qe = qi(ij,l) + gradq3d(ij,l,1)ed(1)+gradq3d(ij,l,2)ed(2)+gradq3d(ij,l,3)*ed(3)
[174]	460	ELSE
	461	ed = cc(ij+u_lup,l,:) - xyz_i(ij+t_lup,:)
[186]	462	! qe = qi(ij+t_lup,l)+sum2(gradq3d(ij+t_lup,l,:),ed)
	463	qe = qi(ij+t_lup,l) + gradq3d(ij+t_lup,l,1)ed(1)+gradq3d(ij+t_lup,l,2)ed(2)+gradq3d(ij+t_lup,l,3)*ed(3)
[174]	464	ENDIF
	465	qflux(ij+u_lup,l) = hfluxt(ij+u_lup,l)*qe
	466
	467	IF (ne(ij,ldown)*hfluxt(ij+u_ldown,l)>0) THEN
	468	ed = cc(ij+u_ldown,l,:) - xyz_i(ij,:)
[186]	469	! qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed)
	470	qe = qi(ij,l) + gradq3d(ij,l,1)ed(1)+gradq3d(ij,l,2)ed(2)+gradq3d(ij,l,3)*ed(3)
[174]	471	ELSE
	472	ed = cc(ij+u_ldown,l,:) - xyz_i(ij+t_ldown,:)
[186]	473	! qe = qi(ij+t_ldown,l)+sum2(gradq3d(ij+t_ldown,l,:),ed)
	474	qe = qi(ij+t_ldown,l)+gradq3d(ij+t_ldown,l,1)ed(1)+gradq3d(ij+t_ldown,l,2)ed(2)+gradq3d(ij+t_ldown,l,3)*ed(3)
[174]	475	ENDIF
	476	qflux(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)*qe
	477	END DO
[22]	478	END DO
[17]	479
[138]	480	! update q and, if update_mass, update mass
[151]	481	DO l = ll_begin,ll_end
[174]	482	!$SIMD
	483	DO ij=ij_begin,ij_end
[138]	484	! sign convention as Ringler et al. (2010) eq. 21
[174]	485	dmass = hfluxt(ij+u_right,l)*ne(ij,right) &
	486	+ hfluxt(ij+u_lup,l) *ne(ij,lup) &
	487	+ hfluxt(ij+u_ldown,l)*ne(ij,ldown) &
	488	+ hfluxt(ij+u_rup,l) *ne(ij,rup) &
	489	+ hfluxt(ij+u_left,l) *ne(ij,left) &
	490	+ hfluxt(ij+u_rdown,l)*ne(ij,rdown)
[22]	491
[174]	492	dq = qflux(ij+u_right,l) *ne(ij,right) &
	493	+ qflux(ij+u_lup,l) *ne(ij,lup) &
	494	+ qflux(ij+u_ldown,l) *ne(ij,ldown) &
	495	+ qflux(ij+u_rup,l) *ne(ij,rup) &
	496	+ qflux(ij+u_left,l) *ne(ij,left) &
	497	+ qflux(ij+u_rdown,l) *ne(ij,rdown)
[17]	498
[174]	499
	500	newmass = mass(ij,l) - dmass/Ai(ij)
	501	qi(ij,l) = (qi(ij,l)*mass(ij,l) - dq/Ai(ij) ) / newmass
	502	IF(update_mass) mass(ij,l) = newmass
	503
[17]	504	END DO
[22]	505	END DO
[17]	506
[148]	507	CALL trace_end("compute_advect_horiz")
[22]	508	CONTAINS
[138]	509
[22]	510	!====================================================================================
[138]	511	PURE REAL(rstd) FUNCTION sum2(a1,a2)
[22]	512	IMPLICIT NONE
[138]	513	REAL(rstd),INTENT(IN):: a1(3), a2(3)
[22]	514	sum2 = a1(1)a2(1)+a1(2)a2(2)+a1(3)*a2(3)
[141]	515	! sum2 = 0. ! Godunov scheme
[22]	516	END FUNCTION sum2
	517
	518	END SUBROUTINE compute_advect_horiz
[138]	519
[17]	520
[22]	521	END MODULE advect_mod

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

Download in other formats: