Context Navigation

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

Last change on this file since 169 was 151, checked in by ymipsl, 11 years ago

Implementation of mixte parallelism MPI/OpenMP into src directory

YM

File size: 14.7 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
[22]	19	INTEGER :: i,j,n
	20
[17]	21	DO j=jj_begin-1,jj_end+1
[22]	22	DO i=ii_begin-1,ii_end+1
	23	n=(j-1)*iim+i
	24
	25	CALL cross_product2(xyz_v(n+z_rdown,:),xyz_v(n+z_rup,:),normal(n+u_right,:))
	26	normal(n+u_right,:)=normal(n+u_right,:)/sqrt(sum(normal(n+u_right,:)*2)+1e-50)ne(n,right)
	27
	28	CALL cross_product2(xyz_v(n+z_up,:),xyz_v(n+z_lup,:),normal(n+u_lup,:))
	29	normal(n+u_lup,:)=normal(n+u_lup,:)/sqrt(sum(normal(n+u_lup,:)*2)+1e-50)ne(n,lup)
	30
	31	CALL cross_product2(xyz_v(n+z_ldown,:),xyz_v(n+z_down,:),normal(n+u_ldown,:))
	32	normal(n+u_ldown,:)=normal(n+u_ldown,:)/sqrt(sum(normal(n+u_ldown,:)*2)+1e-50)ne(n,ldown)
	33
	34	tangent(n+u_right,:)=xyz_v(n+z_rup,:)-xyz_v(n+z_rdown,:)
	35	tangent(n+u_right,:)=tangent(n+u_right,:)/sqrt(sum(tangent(n+u_right,:)**2)+1e-50)
	36
	37	tangent(n+u_lup,:)=xyz_v(n+z_lup,:)-xyz_v(n+z_up,:)
	38	tangent(n+u_lup,:)=tangent(n+u_lup,:)/sqrt(sum(tangent(n+u_lup,:)**2)+1e-50)
	39
	40	tangent(n+u_ldown,:)=xyz_v(n+z_down,:)-xyz_v(n+z_ldown,:)
	41	tangent(n+u_ldown,:)=tangent(n+u_ldown,:)/sqrt(sum(tangent(n+u_ldown,:)**2)+1e-50)
[148]	42
	43	one_over_sqrt_leng(n) = 1./sqrt(max(sum((xyz_v(n+z_up,:) - xyz_i(n,:))2),sum((xyz_v(n+z_down,:) - xyz_i(n,:))2), &
	44	sum((xyz_v(n+z_rup,:) - xyz_i(n,:))2),sum((xyz_v(n+z_rdown,:) - xyz_i(n,:))2), &
	45	sum((xyz_v(n+z_lup,:) - xyz_i(n,:))2),sum((xyz_v(n+z_ldown,:) - xyz_i(n,:))2)) )
	46	ENDDO
[22]	47	ENDDO
	48
[17]	49	END SUBROUTINE init_advect
	50
[22]	51	!=======================================================================================
[17]	52
[148]	53	SUBROUTINE compute_gradq3d(qi,one_over_sqrt_leng,gradq3d)
	54	USE trace
[151]	55	USE omp_para
[22]	56	IMPLICIT NONE
	57	REAL(rstd),INTENT(IN) :: qi(iim*jjm,llm)
[148]	58	REAL(rstd),INTENT(IN) :: one_over_sqrt_leng(iim*jjm)
[23]	59	REAL(rstd),INTENT(OUT) :: gradq3d(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)
	66	INTEGER :: i,j,n,k,ind,l
[138]	67
[148]	68	CALL trace_start("compute_gradq3d")
	69
[138]	70	! TODO : precompute ar, drop arr as output argument of gradq ?
	71
[22]	72	!========================================================================================== GRADIENT
[138]	73	! Compute gradient at triangles solving a linear system
	74	! arr = area of triangle joining centroids of hexagons
[151]	75	DO l = ll_begin,ll_end
[22]	76	DO j=jj_begin-1,jj_end+1
	77	DO i=ii_begin-1,ii_end+1
	78	n=(j-1)*iim+i
[148]	79	CALL gradq(n,l,n+t_rup,n+t_lup,n+z_up,qi,gradtri(n+z_up,l,:),arr(n+z_up))
	80	CALL gradq(n,l,n+t_ldown,n+t_rdown,n+z_down,qi,gradtri(n+z_down,l,:),arr(n+z_down))
[22]	81	END DO
[17]	82	END DO
[22]	83	END DO
[17]	84
[148]	85	DO j=jj_begin-1,jj_end+1
	86	DO i=ii_begin-1,ii_end+1
	87	n=(j-1)*iim+i
	88	ar(n) = arr(n+z_up)+arr(n+z_lup)+arr(n+z_ldown)+arr(n+z_down)+arr(n+z_rdown)+arr(n+z_rup)+1.e-50
	89	ENDDO
	90	ENDDO
	91
	92	DO k=1,3
[151]	93	DO l =ll_begin,ll_end
[148]	94	DO j=jj_begin,jj_end
	95	DO i=ii_begin,ii_end
	96	n=(j-1)*iim+i
	97	gradq3d(n,l,k) = ( gradtri(n+z_up,l,k) + gradtri(n+z_down,l,k) + &
	98	gradtri(n+z_rup,l,k) + gradtri(n+z_ldown,l,k) + &
	99	gradtri(n+z_lup,l,k)+ gradtri(n+z_rdown,l,k) ) / ar(n)
	100	END DO
	101	END DO
	102	END DO
	103	ENDDO
[22]	104
	105	!============================================================================================= LIMITING
[151]	106	DO l =ll_begin,ll_end
[22]	107	DO j=jj_begin,jj_end
[17]	108	DO i=ii_begin,ii_end
	109	n=(j-1)*iim+i
[22]	110	maggrd = dot_product(gradq3d(n,l,:),gradq3d(n,l,:))
	111	maggrd = sqrt(maggrd)
[148]	112	maxq_c = qi(n,l) + maggrd*one_over_sqrt_leng(n)
	113	minq_c = qi(n,l) - maggrd*one_over_sqrt_leng(n)
[22]	114	maxq = max(qi(n,l),qi(n+t_right,l),qi(n+t_lup,l),qi(n+t_rup,l),qi(n+t_left,l), &
	115	qi(n+t_rdown,l),qi(n+t_ldown,l))
	116	minq = min(qi(n,l),qi(n+t_right,l),qi(n+t_lup,l),qi(n+t_rup,l),qi(n+t_left,l), &
	117	qi(n+t_rdown,l),qi(n+t_ldown,l))
	118	alphamx = (maxq - qi(n,l)) ; alphamx = alphamx/(maxq_c - qi(n,l) )
	119	alphamx = max(alphamx,0.0)
	120	alphami = (minq - qi(n,l)); alphami = alphami/(minq_c - qi(n,l))
	121	alphami = max(alphami,0.0)
	122	alpha = min(alphamx,alphami,1.0)
	123	gradq3d(n,l,:) = alpha*gradq3d(n,l,:)
	124	END DO
	125	END DO
[17]	126	END DO
[148]	127
	128	CALL trace_end("compute_gradq3d")
[22]	129	END SUBROUTINE compute_gradq3d
	130
[138]	131	! Backward trajectories, for use with Miura approach
[137]	132	SUBROUTINE compute_backward_traj(normal,tangent,ue,tau, cc)
[148]	133	USE trace
[151]	134	USE omp_para
[22]	135	IMPLICIT NONE
[137]	136	REAL(rstd),INTENT(IN) :: normal(3iimjjm,3)
	137	REAL(rstd),INTENT(IN) :: tangent(3iimjjm,3)
	138	REAL(rstd),INTENT(IN) :: ue(iim3jjm,llm)
	139	REAL(rstd),INTENT(OUT) :: cc(3iimjjm,llm,3) ! start of backward trajectory
	140	REAL(rstd),INTENT(IN) :: tau
	141
[138]	142	REAL(rstd) :: v_e(3), up_e, qe, ed(3)
[137]	143	INTEGER :: i,j,n,l
	144
[148]	145	CALL trace_start("compute_backward_traj")
	146
[138]	147	! TODO : compute normal displacement ue*tau as hfluxt / mass(upwind) then reconstruct tangential displacement
	148
[137]	149	! reconstruct tangential wind then 3D wind at edge then cc = edge midpoint - u*tau
[151]	150	DO l = ll_begin,ll_end
[137]	151	DO j=jj_begin-1,jj_end+1
	152	DO i=ii_begin-1,ii_end+1
	153	n=(j-1)*iim+i
	154
	155	up_e =1/de(n+u_right)*( &
	156	wee(n+u_right,1,1)le(n+u_rup)ue(n+u_rup,l)+ &
	157	wee(n+u_right,2,1)le(n+u_lup)ue(n+u_lup,l)+ &
	158	wee(n+u_right,3,1)le(n+u_left)ue(n+u_left,l)+ &
	159	wee(n+u_right,4,1)le(n+u_ldown)ue(n+u_ldown,l)+ &
	160	wee(n+u_right,5,1)le(n+u_rdown)ue(n+u_rdown,l)+ &
	161	wee(n+u_right,1,2)le(n+t_right+u_ldown)ue(n+t_right+u_ldown,l)+ &
	162	wee(n+u_right,2,2)le(n+t_right+u_rdown)ue(n+t_right+u_rdown,l)+ &
	163	wee(n+u_right,3,2)le(n+t_right+u_right)ue(n+t_right+u_right,l)+ &
	164	wee(n+u_right,4,2)le(n+t_right+u_rup)ue(n+t_right+u_rup,l)+ &
	165	wee(n+u_right,5,2)le(n+t_right+u_lup)ue(n+t_right+u_lup,l) &
	166	)
	167	v_e = ue(n+u_right,l)normal(n+u_right,:) + up_etangent(n+u_right,:)
	168	cc(n+u_right,l,:) = xyz_e(n+u_right,:) - v_e*tau
	169
	170	up_e=1/de(n+u_lup)*( &
	171	wee(n+u_lup,1,1)le(n+u_left)ue(n+u_left,l)+ &
	172	wee(n+u_lup,2,1)le(n+u_ldown)ue(n+u_ldown,l)+ &
	173	wee(n+u_lup,3,1)le(n+u_rdown)ue(n+u_rdown,l)+ &
	174	wee(n+u_lup,4,1)le(n+u_right)ue(n+u_right,l)+ &
	175	wee(n+u_lup,5,1)le(n+u_rup)ue(n+u_rup,l)+ &
	176	wee(n+u_lup,1,2)le(n+t_lup+u_right)ue(n+t_lup+u_right,l)+ &
	177	wee(n+u_lup,2,2)le(n+t_lup+u_rup)ue(n+t_lup+u_rup,l)+ &
	178	wee(n+u_lup,3,2)le(n+t_lup+u_lup)ue(n+t_lup+u_lup,l)+ &
	179	wee(n+u_lup,4,2)le(n+t_lup+u_left)ue(n+t_lup+u_left,l)+ &
	180	wee(n+u_lup,5,2)le(n+t_lup+u_ldown)ue(n+t_lup+u_ldown,l) &
	181	)
	182	v_e = ue(n+u_lup,l)normal(n+u_lup,:) + up_etangent(n+u_lup,:)
	183	cc(n+u_lup,l,:) = xyz_e(n+u_lup,:) - v_e*tau
	184
	185	up_e=1/de(n+u_ldown)*( &
	186	wee(n+u_ldown,1,1)le(n+u_rdown)ue(n+u_rdown,l)+ &
	187	wee(n+u_ldown,2,1)le(n+u_right)ue(n+u_right,l)+ &
	188	wee(n+u_ldown,3,1)le(n+u_rup)ue(n+u_rup,l)+ &
	189	wee(n+u_ldown,4,1)le(n+u_lup)ue(n+u_lup,l)+ &
	190	wee(n+u_ldown,5,1)le(n+u_left)ue(n+u_left,l)+ &
	191	wee(n+u_ldown,1,2)le(n+t_ldown+u_lup)ue(n+t_ldown+u_lup,l)+ &
	192	wee(n+u_ldown,2,2)le(n+t_ldown+u_left)ue(n+t_ldown+u_left,l)+ &
	193	wee(n+u_ldown,3,2)le(n+t_ldown+u_ldown)ue(n+t_ldown+u_ldown,l)+ &
	194	wee(n+u_ldown,4,2)le(n+t_ldown+u_rdown)ue(n+t_ldown+u_rdown,l)+ &
	195	wee(n+u_ldown,5,2)le(n+t_ldown+u_right)ue(n+t_ldown+u_right,l) &
	196	)
	197	v_e = ue(n+u_ldown,l)normal(n+u_ldown,:) + up_etangent(n+u_ldown,:)
	198	cc(n+u_ldown,l,:) = xyz_e(n+u_ldown,:) - v_e*tau
	199	ENDDO
	200	ENDDO
	201	END DO
[148]	202
	203	CALL trace_end("compute_backward_traj")
	204
[137]	205	END SUBROUTINE compute_backward_traj
	206
	207	! Horizontal transport (S. Dubey, T. Dubos)
	208	! Slope-limited van Leer approach with hexagons
[138]	209	SUBROUTINE compute_advect_horiz(update_mass,hfluxt,cc,gradq3d, mass,qi)
[148]	210	USE trace
[151]	211	USE omp_para
[137]	212	IMPLICIT NONE
[138]	213	LOGICAL, INTENT(IN) :: update_mass
	214	REAL(rstd), INTENT(IN) :: gradq3d(iim*jjm,llm,3)
	215	REAL(rstd), INTENT(IN) :: hfluxt(3iimjjm,llm) ! mass flux
	216	REAL(rstd), INTENT(IN) :: cc(3iimjjm,llm,3) ! barycenter of quadrilateral, where q is evaluated (1-point quadrature)
	217	REAL(rstd), INTENT(INOUT) :: mass(iim*jjm,llm)
	218	REAL(rstd), INTENT(INOUT) :: qi(iim*jjm,llm)
[22]	219
[138]	220	REAL(rstd) :: dq,dmass,qe,ed(3), newmass
	221	REAL(rstd) :: qflux(3iimjjm,llm)
[17]	222	INTEGER :: i,j,n,k,l
	223
[148]	224	CALL trace_start("compute_advect_horiz")
	225
[138]	226	! evaluate tracer field at point cc using piecewise linear reconstruction
[136]	227	! q(cc)= q0 + gradq.(cc-xyz_i) with xi centroid of hexagon
[138]	228	! ne*hfluxt>0 iff outgoing
[151]	229	DO l = ll_begin,ll_end
[22]	230	DO j=jj_begin-1,jj_end+1
	231	DO i=ii_begin-1,ii_end+1
	232	n=(j-1)*iim+i
[138]	233	if (ne(n,right)*hfluxt(n+u_right,l)>0) then
[22]	234	ed = cc(n+u_right,l,:) - xyz_i(n,:)
[138]	235	qe = qi(n,l)+sum2(gradq3d(n,l,:),ed)
[22]	236	else
	237	ed = cc(n+u_right,l,:) - xyz_i(n+t_right,:)
[138]	238	qe = qi(n+t_right,l)+sum2(gradq3d(n+t_right,l,:),ed)
[22]	239	endif
[138]	240	qflux(n+u_right,l) = hfluxt(n+u_right,l)*qe
[137]	241
[138]	242	if (ne(n,lup)*hfluxt(n+u_lup,l)>0) then
[22]	243	ed = cc(n+u_lup,l,:) - xyz_i(n,:)
[138]	244	qe = qi(n,l)+sum2(gradq3d(n,l,:),ed)
[22]	245	else
	246	ed = cc(n+u_lup,l,:) - xyz_i(n+t_lup,:)
[138]	247	qe = qi(n+t_lup,l)+sum2(gradq3d(n+t_lup,l,:),ed)
[22]	248	endif
[138]	249	qflux(n+u_lup,l) = hfluxt(n+u_lup,l)*qe
[137]	250
[138]	251	if (ne(n,ldown)*hfluxt(n+u_ldown,l)>0) then
[22]	252	ed = cc(n+u_ldown,l,:) - xyz_i(n,:)
[138]	253	qe = qi(n,l)+sum2(gradq3d(n,l,:),ed)
[22]	254	else
	255	ed = cc(n+u_ldown,l,:) - xyz_i(n+t_ldown,:)
[138]	256	qe = qi(n+t_ldown,l)+sum2(gradq3d(n+t_ldown,l,:),ed)
[22]	257	endif
[138]	258	qflux(n+u_ldown,l) = hfluxt(n+u_ldown,l)*qe
[22]	259	end do
	260	end do
	261	END DO
[17]	262
[138]	263	! update q and, if update_mass, update mass
[151]	264	DO l = ll_begin,ll_end
[138]	265	DO j=jj_begin,jj_end
	266	DO i=ii_begin,ii_end
	267	n=(j-1)*iim+i
	268	! sign convention as Ringler et al. (2010) eq. 21
	269	dmass = hfluxt(n+u_right,l)*ne(n,right) &
	270	+ hfluxt(n+u_lup,l) *ne(n,lup) &
	271	+ hfluxt(n+u_ldown,l)*ne(n,ldown) &
	272	+ hfluxt(n+u_rup,l) *ne(n,rup) &
	273	+ hfluxt(n+u_left,l) *ne(n,left) &
	274	+ hfluxt(n+u_rdown,l)*ne(n,rdown)
[22]	275
[138]	276	dq = qflux(n+u_right,l) *ne(n,right) &
	277	+ qflux(n+u_lup,l) *ne(n,lup) &
	278	+ qflux(n+u_ldown,l) *ne(n,ldown) &
	279	+ qflux(n+u_rup,l) *ne(n,rup) &
	280	+ qflux(n+u_left,l) *ne(n,left) &
	281	+ qflux(n+u_rdown,l) *ne(n,rdown)
[17]	282
[138]	283	newmass = mass(n,l) - dmass/Ai(n)
	284	qi(n,l) = (qi(n,l)*mass(n,l) - dq/Ai(n) ) / newmass
	285	IF(update_mass) mass(n,l) = newmass
	286	END DO
[17]	287	END DO
[22]	288	END DO
[17]	289
[148]	290	CALL trace_end("compute_advect_horiz")
[22]	291	CONTAINS
[138]	292
[22]	293	!====================================================================================
[138]	294	PURE REAL(rstd) FUNCTION sum2(a1,a2)
[22]	295	IMPLICIT NONE
[138]	296	REAL(rstd),INTENT(IN):: a1(3), a2(3)
[22]	297	sum2 = a1(1)a2(1)+a1(2)a2(2)+a1(3)*a2(3)
[141]	298	! sum2 = 0. ! Godunov scheme
[22]	299	END FUNCTION sum2
	300
	301	END SUBROUTINE compute_advect_horiz
[138]	302
[22]	303	!==========================================================================
[148]	304	PURE SUBROUTINE gradq(n0,l,n1,n2,n3,q,dq,det)
[22]	305	IMPLICIT NONE
[148]	306	INTEGER, INTENT(IN) :: n0,l,n1,n2,n3
	307	REAL(rstd), INTENT(IN) :: q(iim*jjm,llm)
[138]	308	REAL(rstd), INTENT(OUT) :: dq(3), det
	309
	310	REAL(rstd) ::detx,dety,detz
[17]	311	INTEGER :: info
	312	INTEGER :: IPIV(3)
	313
	314	REAL(rstd) :: A(3,3)
	315	REAL(rstd) :: B(3)
	316
[138]	317	! TODO : replace A by A1,A2,A3
	318	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)
	319	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)
	320	A(3,1)=xyz_v(n3,1); A(3,2)= xyz_v(n3,2); A(3,3)=xyz_v(n3,3)
[22]	321
[148]	322	dq(1) = q(n1,l)-q(n0,l)
	323	dq(2) = q(n2,l)-q(n0,l)
[17]	324	dq(3) = 0.0
[22]	325	! CALL DGESV(3,1,A,3,IPIV,dq(:),3,info)
	326	CALL determinant(A(:,1),A(:,2),A(:,3),det)
	327	CALL determinant(dq,A(:,2),A(:,3),detx)
	328	CALL determinant(A(:,1),dq,A(:,3),dety)
	329	CALL determinant(A(:,1),A(:,2),dq,detz)
	330	dq(1) = detx
	331	dq(2) = dety
	332	dq(3) = detz
[17]	333	END SUBROUTINE gradq
[138]	334
[22]	335	!==========================================================================
[138]	336	PURE SUBROUTINE determinant(a1,a2,a3,det)
[22]	337	IMPLICIT NONE
[138]	338	REAL(rstd), DIMENSION(3), INTENT(IN) :: a1,a2,a3
	339	REAL(rstd), INTENT(OUT) :: det
	340	REAL(rstd) :: x1,x2,x3
[22]	341	x1 = a1(1) * (a2(2) * a3(3) - a2(3) * a3(2))
	342	x2 = a1(2) * (a2(1) * a3(3) - a2(3) * a3(1))
	343	x3 = a1(3) * (a2(1) * a3(2) - a2(2) * a3(1))
	344	det = x1 - x2 + x3
	345	END SUBROUTINE determinant
	346
	347	END MODULE advect_mod

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

Download in other formats: