Context Navigation

source: codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90 @ 842

Last change on this file since 842 was 842, checked in by jisesh, 5 years ago
devel: separate module for compute_caldyn_kv
File size: 39.3 KB

Rev	Line
[362]	1	MODULE caldyn_kernels_hevi_mod
	2	USE icosa
[369]	3	USE trace
	4	USE omp_para
	5	USE disvert_mod
[362]	6	USE transfert_mod
[731]	7	USE caldyn_vars_mod
[362]	8	IMPLICIT NONE
	9	PRIVATE
	10
[562]	11	REAL(rstd), PARAMETER :: pbot=1e5, rho_bot=1e6
[368]	12
[538]	13	LOGICAL, SAVE :: debug_hevi_solver = .FALSE.
[368]	14
[842]	15	PUBLIC :: compute_caldyn_Coriolis, &
[369]	16	compute_caldyn_slow_hydro, compute_caldyn_slow_NH, &
[366]	17	compute_caldyn_solver, compute_caldyn_fast
[362]	18
	19	CONTAINS
	20
[562]	21	SUBROUTINE compute_NH_geopot(tau, phis, m_ik, m_il, theta, W_il, Phi_il)
[368]	22	REAL(rstd),INTENT(IN) :: tau ! solve Phi-tau*dPhi/dt = Phi_rhs
[562]	23	REAL(rstd),INTENT(IN) :: phis(iim*jjm)
[368]	24	REAL(rstd),INTENT(IN) :: m_ik(iim*jjm,llm)
	25	REAL(rstd),INTENT(IN) :: m_il(iim*jjm,llm+1)
	26	REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm)
	27	REAL(rstd),INTENT(IN) :: W_il(iim*jjm,llm+1) ! vertical momentum
	28	REAL(rstd),INTENT(INOUT) :: Phi_il(iim*jjm,llm+1) ! geopotential
	29
	30	REAL(rstd) :: Phi_star_il(iim*jjm,llm+1)
	31	REAL(rstd) :: p_ik(iim*jjm,llm) ! pressure
	32	REAL(rstd) :: R_il(iim*jjm,llm+1) ! rhs of tridiag problem
	33	REAL(rstd) :: x_il(iim*jjm,llm+1) ! solution of tridiag problem
	34	REAL(rstd) :: A_ik(iim*jjm,llm) ! off-diagonal coefficients of tridiag problem
	35	REAL(rstd) :: B_il(iim*jjm,llm+1) ! diagonal coefficients of tridiag problem
	36	REAL(rstd) :: C_ik(iim*jjm,llm) ! Thomas algorithm
	37	REAL(rstd) :: D_il(iim*jjm,llm+1) ! Thomas algorithm
	38	REAL(rstd) :: gamma, rho_ij, X_ij, Y_ij
[657]	39	REAL(rstd) :: wil, tau2_g, g2, gm2, ml_g2, c2_mik, vreff
[368]	40
[538]	41	INTEGER :: iter, ij, l, ij_omp_begin_ext, ij_omp_end_ext
[368]	42
[603]	43	CALL distrib_level(ij_begin_ext,ij_end_ext, ij_omp_begin_ext,ij_omp_end_ext)
[538]	44
[573]	45	IF(dysl) THEN
[562]	46	#define PHI_BOT(ij) phis(ij)
[612]	47	#include "../kernels_hex/compute_NH_geopot.k90"
[573]	48	#undef PHI_BOT
	49	ELSE
[368]	50	! FIXME : vertical OpenMP parallelism will not work
	51
	52	tau2_g=tau*tau/g
	53	g2=g*g
	54	gm2 = g**-2
	55	gamma = 1./(1.-kappa)
	56
	57	! compute Phi_star
	58	DO l=1,llm+1
	59	!DIR$ SIMD
	60	DO ij=ij_begin_ext,ij_end_ext
	61	Phi_star_il(ij,l) = Phi_il(ij,l) + taug2(W_il(ij,l)/m_il(ij,l)-tau)
	62	ENDDO
	63	ENDDO
	64
	65	! Newton-Raphson iteration : Phi_il contains current guess value
[377]	66	DO iter=1,5 ! 2 iterations should be enough
[368]	67
	68	! Compute pressure, A_ik
	69	DO l=1,llm
	70	!DIR$ SIMD
	71	DO ij=ij_begin_ext,ij_end_ext
	72	rho_ij = (g*m_ik(ij,l))/(Phi_il(ij,l+1)-Phi_il(ij,l))
	73	X_ij = (cpp/preff)kappatheta(ij,l)*rho_ij
	74	p_ik(ij,l) = preff(X_ij*gamma)
	75	c2_mik = gammap_ik(ij,l)/(rho_ijm_ik(ij,l)) ! c^2 = gammaRT = gamma*p/rho
	76	A_ik(ij,l) = c2_mik(tau/grho_ij)**2
	77	ENDDO
	78	ENDDO
	79
	80	! Compute residual, B_il
	81	! bottom interface l=1
	82	!DIR$ SIMD
	83	DO ij=ij_begin_ext,ij_end_ext
	84	ml_g2 = gm2*m_il(ij,1)
	85	B_il(ij,1) = A_ik(ij,1) + ml_g2 + tau2_g*rho_bot
	86	R_il(ij,1) = ml_g2*( Phi_il(ij,1)-Phi_star_il(ij,1)) &
[565]	87	+ tau2_g( p_ik(ij,1)-pbot+rho_bot(Phi_il(ij,1)-phis(ij)) )
[368]	88	ENDDO
	89	! inner interfaces
	90	DO l=2,llm
	91	!DIR$ SIMD
	92	DO ij=ij_begin_ext,ij_end_ext
	93	ml_g2 = gm2*m_il(ij,l)
	94	B_il(ij,l) = A_ik(ij,l)+A_ik(ij,l-1) + ml_g2
	95	R_il(ij,l) = ml_g2*( Phi_il(ij,l)-Phi_star_il(ij,l)) &
	96	+ tau2_g*(p_ik(ij,l)-p_ik(ij,l-1))
	97	! consistency check : if Wil=0 and initial state is in hydrostatic balance
	98	! then Phi_star_il(ij,l) = Phi_il(ij,l) - tau^2*g^2
	99	! and residual = tau^2*(ml+(1/g)dl_pi)=0
	100	ENDDO
	101	ENDDO
	102	! top interface l=llm+1
	103	!DIR$ SIMD
	104	DO ij=ij_begin_ext,ij_end_ext
	105	ml_g2 = gm2*m_il(ij,llm+1)
	106	B_il(ij,llm+1) = A_ik(ij,llm) + ml_g2
	107	R_il(ij,llm+1) = ml_g2*( Phi_il(ij,llm+1)-Phi_star_il(ij,llm+1)) &
	108	+ tau2_g*( ptop-p_ik(ij,llm) )
	109	ENDDO
	110
	111	! FIXME later
	112	! the lines below modify the tridiag problem
	113	! for flat, rigid boundary conditions at top and bottom :
	114	! zero out A(1), A(llm), R(1), R(llm+1)
	115	! => x(l)=0 at l=1,llm+1
	116	DO ij=ij_begin_ext,ij_end_ext
	117	A_ik(ij,1) = 0.
	118	A_ik(ij,llm) = 0.
	119	R_il(ij,1) = 0.
	120	R_il(ij,llm+1) = 0.
	121	ENDDO
	122
	123	IF(debug_hevi_solver) THEN ! print Linf(residual)
	124	PRINT *, '[hevi_solver] R,p', iter, MAXVAL(ABS(R_il)), MAXVAL(p_ik)
	125	END IF
	126
	127	! Solve -A(l-1)x(l-1) + B(l)x(l) - A(l)x(l+1) = R(l) using Thomas algorithm
	128	! Forward sweep :
	129	! C(0)=0, C(l) = -A(l) / (B(l)+A(l-1)C(l-1)),
	130	! D(0)=0, D(l) = (R(l)+A(l-1)D(l-1)) / (B(l)+A(l-1)C(l-1))
	131	! bottom interface l=1
	132	!DIR$ SIMD
	133	DO ij=ij_begin_ext,ij_end_ext
	134	X_ij = 1./B_il(ij,1)
	135	C_ik(ij,1) = -A_ik(ij,1) * X_ij
	136	D_il(ij,1) = R_il(ij,1) * X_ij
	137	ENDDO
	138	! inner interfaces/layers
	139	DO l=2,llm
	140	!DIR$ SIMD
	141	DO ij=ij_begin_ext,ij_end_ext
	142	X_ij = 1./(B_il(ij,l) + A_ik(ij,l-1)*C_ik(ij,l-1))
	143	C_ik(ij,l) = -A_ik(ij,l) * X_ij
	144	D_il(ij,l) = (R_il(ij,l)+A_ik(ij,l-1)D_il(ij,l-1)) X_ij
	145	ENDDO
	146	ENDDO
	147	! top interface l=llm+1
	148	!DIR$ SIMD
	149	DO ij=ij_begin_ext,ij_end_ext
	150	X_ij = 1./(B_il(ij,llm+1) + A_ik(ij,llm)*C_ik(ij,llm))
	151	D_il(ij,llm+1) = (R_il(ij,llm+1)+A_ik(ij,llm)D_il(ij,llm)) X_ij
	152	ENDDO
	153
	154	! Back substitution :
	155	! x(i) = D(i)-C(i)x(i+1), x(N+1)=0
	156	! + Newton-Raphson update
	157	x_il=0. ! FIXME
	158	! top interface l=llm+1
	159	!DIR$ SIMD
	160	DO ij=ij_begin_ext,ij_end_ext
	161	x_il(ij,llm+1) = D_il(ij,llm+1)
	162	Phi_il(ij,llm+1) = Phi_il(ij,llm+1) - x_il(ij,llm+1)
	163	ENDDO
	164	! lower interfaces
	165	DO l=llm,1,-1
	166	!DIR$ SIMD
	167	DO ij=ij_begin_ext,ij_end_ext
	168	x_il(ij,l) = D_il(ij,l) - C_ik(ij,l)*x_il(ij,l+1)
	169	Phi_il(ij,l) = Phi_il(ij,l) - x_il(ij,l)
	170	ENDDO
	171	ENDDO
	172
	173	IF(debug_hevi_solver) THEN
	174	PRINT *, '[hevi_solver] A,B', iter, MAXVAL(ABS(A_ik)),MAXVAL(ABS(B_il))
	175	PRINT *, '[hevi_solver] C,D', iter, MAXVAL(ABS(C_ik)),MAXVAL(ABS(D_il))
	176	DO l=1,llm+1
[821]	177	WRITE(*,'(A,I2.1,I3.2,E9.2)') '[hevi_solver] x', iter,l, MAXVAL(ABS(x_il(:,l)))
[368]	178	END DO
	179	END IF
	180
	181	END DO ! Newton-Raphson
[538]	182
[573]	183	END IF ! dysl
[368]	184
	185	END SUBROUTINE compute_NH_geopot
	186
[562]	187	SUBROUTINE compute_caldyn_solver(tau,phis, rhodz,theta,pk, geopot,W, m_il,pres, dPhi,dW,du)
[366]	188	REAL(rstd),INTENT(IN) :: tau ! "solve" Phi-tau*dPhi/dt = Phi_rhs
[562]	189	REAL(rstd),INTENT(IN) :: phis(iim*jjm)
[366]	190	REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm)
[538]	191	REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm,nqdyn)
[366]	192	REAL(rstd),INTENT(OUT) :: pk(iim*jjm,llm)
	193	REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)
	194	REAL(rstd),INTENT(INOUT) :: W(iim*jjm,llm+1) ! OUT if tau>0
[558]	195	REAL(rstd),INTENT(OUT) :: m_il(iim*jjm,llm+1) ! rhodz averaged to interfaces
	196	REAL(rstd),INTENT(OUT) :: pres(iim*jjm,llm) ! pressure
[366]	197	REAL(rstd),INTENT(OUT) :: dW(iim*jjm,llm+1)
	198	REAL(rstd),INTENT(OUT) :: dPhi(iim*jjm,llm+1)
[369]	199	REAL(rstd),INTENT(OUT) :: du(3iimjjm,llm)
[366]	200
[558]	201	REAL(rstd) :: berni(iim*jjm,llm) ! (W/m_il)^2
[573]	202	REAL(rstd) :: berni1(iim*jjm) ! (W/m_il)^2
[657]	203	REAL(rstd) :: gamma, rho_ij, T_ij, X_ij, Y_ij, vreff, Rd, Cvd, Rd_preff
[368]	204	INTEGER :: ij, l
[366]	205
	206	CALL trace_start("compute_caldyn_solver")
	207
[538]	208	Rd=cpp*kappa
	209
[573]	210	IF(dysl) THEN
	211
[558]	212	!$OMP BARRIER
[657]	213
	214	#include "../kernels_hex/caldyn_mil.k90"
	215	IF(tau>0) THEN ! solve implicit problem for geopotential
	216	CALL compute_NH_geopot(tau,phis, rhodz, m_il, theta, W, geopot)
	217	END IF
[562]	218	#define PHI_BOT(ij) phis(ij)
[612]	219	#include "../kernels_hex/caldyn_solver.k90"
[573]	220	#undef PHI_BOT
[558]	221	!$OMP BARRIER
[573]	222
	223	ELSE
	224
	225	#define BERNI(ij) berni1(ij)
[368]	226	! FIXME : vertical OpenMP parallelism will not work
[366]	227
[368]	228	! average m_ik to interfaces => m_il
	229	!DIR$ SIMD
	230	DO ij=ij_begin_ext,ij_end_ext
	231	m_il(ij,1) = .5*rhodz(ij,1)
	232	ENDDO
	233	DO l=2,llm
	234	!DIR$ SIMD
	235	DO ij=ij_begin_ext,ij_end_ext
	236	m_il(ij,l) = .5*(rhodz(ij,l-1)+rhodz(ij,l))
	237	ENDDO
	238	ENDDO
	239	!DIR$ SIMD
	240	DO ij=ij_begin_ext,ij_end_ext
	241	m_il(ij,llm+1) = .5*rhodz(ij,llm)
	242	ENDDO
	243
	244	IF(tau>0) THEN ! solve implicit problem for geopotential
[565]	245	CALL compute_NH_geopot(tau, phis, rhodz, m_il, theta, W, geopot)
[366]	246	END IF
	247
	248	! Compute pressure, stored temporarily in pk
	249	! kappa = R/Cp
	250	! 1-kappa = Cv/Cp
	251	! Cp/Cv = 1/(1-kappa)
	252	gamma = 1./(1.-kappa)
[368]	253	DO l=1,llm
[366]	254	!DIR$ SIMD
[368]	255	DO ij=ij_begin_ext,ij_end_ext
[366]	256	rho_ij = (g*rhodz(ij,l))/(geopot(ij,l+1)-geopot(ij,l))
[538]	257	X_ij = (cpp/preff)kappatheta(ij,l,1)*rho_ij
[366]	258	! kappa.theta.rho = p/exner
	259	! => X = (p/p0)/(exner/Cp)
	260	! = (p/p0)^(1-kappa)
	261	pk(ij,l) = preff(X_ij*gamma)
	262	ENDDO
	263	ENDDO
	264
[369]	265	! Update W, compute tendencies
[368]	266	DO l=2,llm
[366]	267	!DIR$ SIMD
[368]	268	DO ij=ij_begin_ext,ij_end_ext
	269	dW(ij,l) = (1./g)*(pk(ij,l-1)-pk(ij,l)) - m_il(ij,l)
	270	W(ij,l) = W(ij,l)+tau*dW(ij,l) ! update W
	271	dPhi(ij,l) = ggW(ij,l)/m_il(ij,l)
[366]	272	ENDDO
	273	! PRINT *,'Max dPhi', l,ij_begin,ij_end, MAXVAL(abs(dPhi(ij_begin:ij_end,l)))
	274	! PRINT *,'Max dW', l,ij_begin,ij_end, MAXVAL(abs(dW(ij_begin:ij_end,l)))
	275	ENDDO
	276	! Lower BC (FIXME : no orography yet !)
	277	DO ij=ij_begin,ij_end
	278	dPhi(ij,1)=0
	279	W(ij,1)=0
	280	dW(ij,1)=0
	281	dPhi(ij,llm+1)=0 ! rigid lid
	282	W(ij,llm+1)=0
	283	dW(ij,llm+1)=0
	284	ENDDO
	285	! Upper BC p=ptop
[368]	286	! DO ij=ij_omp_begin_ext,ij_omp_end_ext
	287	! dPhi(ij,llm+1) = W(ij,llm+1)/rhodz(ij,llm)
	288	! dW(ij,llm+1) = (1./g)(pk(ij,llm)-ptop) - .5rhodz(ij,llm)
	289	! ENDDO
[366]	290
[375]	291	! Compute Exner function (needed by compute_caldyn_fast) and du=-g^2.grad(w^2)
[368]	292	DO l=1,llm
[366]	293	!DIR$ SIMD
[368]	294	DO ij=ij_begin_ext,ij_end_ext
[366]	295	pk(ij,l) = cpp((pk(ij,l)/preff)*kappa) ! other formulae possible if exponentiation is slow
[538]	296	BERNI(ij) = (-.25gg)*( &
[375]	297	(W(ij,l)/m_il(ij,l))**2 &
[369]	298	+ (W(ij,l+1)/m_il(ij,l+1))**2 )
[366]	299	ENDDO
[369]	300	DO ij=ij_begin,ij_end
[538]	301	du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right))
	302	du(ij+u_lup,l) = ne_lup *(BERNI(ij)-BERNI(ij+t_lup))
	303	du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))
[369]	304	ENDDO
[366]	305	ENDDO
[538]	306	#undef BERNI
	307
[573]	308	END IF ! dysl
	309
[366]	310	CALL trace_end("compute_caldyn_solver")
	311
	312	END SUBROUTINE compute_caldyn_solver
	313
	314	SUBROUTINE compute_caldyn_fast(tau,u,rhodz,theta,pk,geopot,du)
	315	REAL(rstd),INTENT(IN) :: tau ! "solve" u-tau*du/dt = rhs
	316	REAL(rstd),INTENT(INOUT) :: u(iim3jjm,llm) ! OUT if tau>0
	317	REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm)
[405]	318	REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm,nqdyn)
[366]	319	REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm)
	320	REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)
[369]	321	REAL(rstd),INTENT(INOUT) :: du(iim3jjm,llm)
[362]	322	REAL(rstd) :: berni(iim*jjm,llm) ! Bernoulli function
[405]	323	REAL(rstd) :: berniv(iim*jjm,llm) ! moist Bernoulli function
[362]	324
	325	INTEGER :: i,j,ij,l
[837]	326	REAL(rstd) :: cp_ik, qv, temp, chi, nu, due, due_right, due_lup, due_ldown
[362]	327
	328	CALL trace_start("compute_caldyn_fast")
[366]	329
[562]	330	IF(dysl_caldyn_fast) THEN
[612]	331	#include "../kernels_hex/caldyn_fast.k90"
[562]	332	ELSE
	333
[366]	334	! Compute Bernoulli term
[362]	335	IF(boussinesq) THEN
	336	DO l=ll_begin,ll_end
	337	!DIR$ SIMD
	338	DO ij=ij_begin,ij_end
	339	berni(ij,l) = pk(ij,l)
	340	! from now on pk contains the vertically-averaged geopotential
	341	pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
[401]	342	END DO
	343	END DO
[362]	344	ELSE ! compressible
	345
	346	DO l=ll_begin,ll_end
[401]	347	SELECT CASE(caldyn_thermo)
	348	CASE(thermo_theta) ! vdp = theta.dpi => B = Phi
	349	!DIR$ SIMD
	350	DO ij=ij_begin,ij_end
	351	berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
	352	END DO
	353	CASE(thermo_entropy) ! vdp = dG + sdT => B = Phi + G, G=h-Ts=T*(cpp-s)
	354	!DIR$ SIMD
	355	DO ij=ij_begin,ij_end
	356	berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
[405]	357	+ pk(ij,l)*(cpp-theta(ij,l,1)) ! pk=temperature, theta=entropy
[401]	358	END DO
[405]	359	CASE(thermo_moist)
	360	!DIR$ SIMD
	361	DO ij=ij_begin,ij_end
	362	! du/dt = grad(Bd)+rv.grad(Bv)+s.grad(T)
	363	! Bd = Phi + gibbs_d
	364	! Bv = Phi + gibbs_v
	365	! pk=temperature, theta=entropy
	366	qv = theta(ij,l,2)
	367	temp = pk(ij,l)
	368	chi = log(temp/Treff)
	369	nu = (chi(cpp+qvcppv)-theta(ij,l,1))/(Rd+qv*Rv) ! log(p/preff)
	370	berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
	371	+ temp(cpp(1.-chi)+Rd*nu)
	372	berniv(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
	373	+ temp(cppv(1.-chi)+Rv*nu)
	374	END DO
[401]	375	END SELECT
	376	END DO
[362]	377
	378	END IF ! Boussinesq/compressible
	379
[369]	380	!!! u:=u+tau*du, du = -grad(B)-theta.grad(pi)
[362]	381	DO l=ll_begin,ll_end
[405]	382	IF(caldyn_thermo == thermo_moist) THEN
	383	!DIR$ SIMD
	384	DO ij=ij_begin,ij_end
	385	due_right = berni(ij+t_right,l)-berni(ij,l) &
	386	+ 0.5*(theta(ij,l,1)+theta(ij+t_right,l,1)) &
	387	*(pk(ij+t_right,l)-pk(ij,l)) &
	388	+ 0.5*(theta(ij,l,2)+theta(ij+t_right,l,2)) &
	389	*(berniv(ij+t_right,l)-berniv(ij,l))
	390
	391	due_lup = berni(ij+t_lup,l)-berni(ij,l) &
	392	+ 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1)) &
	393	*(pk(ij+t_lup,l)-pk(ij,l)) &
	394	+ 0.5*(theta(ij,l,2)+theta(ij+t_lup,l,2)) &
	395	*(berniv(ij+t_lup,l)-berniv(ij,l))
	396
	397	due_ldown = berni(ij+t_ldown,l)-berni(ij,l) &
	398	+ 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1)) &
	399	*(pk(ij+t_ldown,l)-pk(ij,l)) &
	400	+ 0.5*(theta(ij,l,2)+theta(ij+t_ldown,l,2)) &
	401	*(berniv(ij+t_ldown,l)-berniv(ij,l))
	402
	403	du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
	404	du(ij+u_lup,l) = du(ij+u_lup,l) - ne_lup*due_lup
	405	du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
	406	u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
	407	u(ij+u_lup,l) = u(ij+u_lup,l) + tau*du(ij+u_lup,l)
	408	u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
	409	END DO
	410	ELSE
	411	!DIR$ SIMD
	412	DO ij=ij_begin,ij_end
	413	due_right = 0.5*(theta(ij,l,1)+theta(ij+t_right,l,1)) &
	414	*(pk(ij+t_right,l)-pk(ij,l)) &
	415	+ berni(ij+t_right,l)-berni(ij,l)
	416	due_lup = 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1)) &
	417	*(pk(ij+t_lup,l)-pk(ij,l)) &
	418	+ berni(ij+t_lup,l)-berni(ij,l)
	419	due_ldown = 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1)) &
	420	*(pk(ij+t_ldown,l)-pk(ij,l)) &
	421	+ berni(ij+t_ldown,l)-berni(ij,l)
	422	du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
	423	du(ij+u_lup,l) = du(ij+u_lup,l) - ne_lup*due_lup
	424	du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
	425	u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
	426	u(ij+u_lup,l) = u(ij+u_lup,l) + tau*du(ij+u_lup,l)
	427	u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
	428	END DO
	429	END IF
	430	END DO
[562]	431
	432	END IF ! dysl
[362]	433	CALL trace_end("compute_caldyn_fast")
	434
	435	END SUBROUTINE compute_caldyn_fast
	436
[369]	437	SUBROUTINE compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du)
	438	REAL(rstd),INTENT(IN) :: hflux(3iimjjm,llm) ! hflux in kg/s
[404]	439	REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm,nqdyn) ! active scalars
[369]	440	REAL(rstd),INTENT(IN) :: qu(3iimjjm,llm)
[362]	441	REAL(rstd),INTENT(OUT) :: convm(iim*jjm,llm) ! mass flux convergence
[404]	442	REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm,nqdyn)
[369]	443	REAL(rstd),INTENT(INOUT) :: du(3iimjjm,llm)
	444
[538]	445	REAL(rstd) :: Ftheta(3iimjjm,llm) ! potential temperature flux
	446	REAL(rstd) :: uu_right, uu_lup, uu_ldown, du_trisk, divF
[404]	447	INTEGER :: ij,iq,l,kdown
[362]	448
[369]	449	CALL trace_start("compute_caldyn_Coriolis")
[362]	450
[562]	451	IF(dysl_caldyn_coriolis) THEN
[573]	452
[612]	453	#include "../kernels_hex/coriolis.k90"
[562]	454
	455	ELSE
[538]	456	#define FTHETA(ij) Ftheta(ij,1)
	457
[369]	458	DO l=ll_begin, ll_end
	459	! compute theta flux
[426]	460	DO iq=1,nqdyn
[362]	461	!DIR$ SIMD
[404]	462	DO ij=ij_begin_ext,ij_end_ext
[538]	463	FTHETA(ij+u_right) = 0.5*(theta(ij,l,iq)+theta(ij+t_right,l,iq)) &
[369]	464	* hflux(ij+u_right,l)
[538]	465	FTHETA(ij+u_lup) = 0.5*(theta(ij,l,iq)+theta(ij+t_lup,l,iq)) &
[404]	466	* hflux(ij+u_lup,l)
[538]	467	FTHETA(ij+u_ldown) = 0.5*(theta(ij,l,iq)+theta(ij+t_ldown,l,iq)) &
[404]	468	* hflux(ij+u_ldown,l)
	469	END DO
	470	! horizontal divergence of fluxes
[426]	471	!DIR$ SIMD
[404]	472	DO ij=ij_begin,ij_end
	473	! dtheta_rhodz = -div(flux.theta)
	474	dtheta_rhodz(ij,l,iq)= &
[538]	475	-1./Ai(ij)(ne_rightFTHETA(ij+u_right) + &
	476	ne_rup*FTHETA(ij+u_rup) + &
	477	ne_lup*FTHETA(ij+u_lup) + &
	478	ne_left*FTHETA(ij+u_left) + &
	479	ne_ldown*FTHETA(ij+u_ldown) + &
	480	ne_rdown*FTHETA(ij+u_rdown) )
[404]	481	END DO
	482	END DO
	483
[426]	484	!DIR$ SIMD
[362]	485	DO ij=ij_begin,ij_end
	486	! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
	487	convm(ij,l)= -1./Ai(ij)(ne_righthflux(ij+u_right,l) + &
	488	ne_rup*hflux(ij+u_rup,l) + &
	489	ne_lup*hflux(ij+u_lup,l) + &
	490	ne_left*hflux(ij+u_left,l) + &
	491	ne_ldown*hflux(ij+u_ldown,l) + &
[404]	492	ne_rdown*hflux(ij+u_rdown,l))
	493	END DO ! ij
	494	END DO ! llm
[362]	495
	496	!!! Compute potential vorticity (Coriolis) contribution to du
[369]	497	SELECT CASE(caldyn_conserv)
[362]	498
[733]	499	CASE(conserv_energy) ! energy-conserving TRiSK
[362]	500
	501	DO l=ll_begin,ll_end
	502	!DIR$ SIMD
	503	DO ij=ij_begin,ij_end
	504	uu_right = &
	505	wee(ij+u_right,1,1)hflux(ij+u_rup,l)(qu(ij+u_right,l)+qu(ij+u_rup,l))+ &
	506	wee(ij+u_right,2,1)hflux(ij+u_lup,l)(qu(ij+u_right,l)+qu(ij+u_lup,l))+ &
	507	wee(ij+u_right,3,1)hflux(ij+u_left,l)(qu(ij+u_right,l)+qu(ij+u_left,l))+ &
	508	wee(ij+u_right,4,1)hflux(ij+u_ldown,l)(qu(ij+u_right,l)+qu(ij+u_ldown,l))+ &
	509	wee(ij+u_right,5,1)hflux(ij+u_rdown,l)(qu(ij+u_right,l)+qu(ij+u_rdown,l))+ &
	510	wee(ij+u_right,1,2)hflux(ij+t_right+u_ldown,l)(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l))+ &
	511	wee(ij+u_right,2,2)hflux(ij+t_right+u_rdown,l)(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l))+ &
	512	wee(ij+u_right,3,2)hflux(ij+t_right+u_right,l)(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+ &
	513	wee(ij+u_right,4,2)hflux(ij+t_right+u_rup,l)(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l))+ &
	514	wee(ij+u_right,5,2)hflux(ij+t_right+u_lup,l)(qu(ij+u_right,l)+qu(ij+t_right+u_lup,l))
	515	uu_lup = &
	516	wee(ij+u_lup,1,1)hflux(ij+u_left,l)(qu(ij+u_lup,l)+qu(ij+u_left,l)) + &
	517	wee(ij+u_lup,2,1)hflux(ij+u_ldown,l)(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) + &
	518	wee(ij+u_lup,3,1)hflux(ij+u_rdown,l)(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) + &
	519	wee(ij+u_lup,4,1)hflux(ij+u_right,l)(qu(ij+u_lup,l)+qu(ij+u_right,l)) + &
	520	wee(ij+u_lup,5,1)hflux(ij+u_rup,l)(qu(ij+u_lup,l)+qu(ij+u_rup,l)) + &
	521	wee(ij+u_lup,1,2)hflux(ij+t_lup+u_right,l)(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l)) + &
	522	wee(ij+u_lup,2,2)hflux(ij+t_lup+u_rup,l)(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l)) + &
	523	wee(ij+u_lup,3,2)hflux(ij+t_lup+u_lup,l)(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) + &
	524	wee(ij+u_lup,4,2)hflux(ij+t_lup+u_left,l)(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l)) + &
	525	wee(ij+u_lup,5,2)hflux(ij+t_lup+u_ldown,l)(qu(ij+u_lup,l)+qu(ij+t_lup+u_ldown,l))
	526	uu_ldown = &
	527	wee(ij+u_ldown,1,1)hflux(ij+u_rdown,l)(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) + &
	528	wee(ij+u_ldown,2,1)hflux(ij+u_right,l)(qu(ij+u_ldown,l)+qu(ij+u_right,l)) + &
	529	wee(ij+u_ldown,3,1)hflux(ij+u_rup,l)(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) + &
	530	wee(ij+u_ldown,4,1)hflux(ij+u_lup,l)(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) + &
	531	wee(ij+u_ldown,5,1)hflux(ij+u_left,l)(qu(ij+u_ldown,l)+qu(ij+u_left,l)) + &
	532	wee(ij+u_ldown,1,2)hflux(ij+t_ldown+u_lup,l)(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l)) + &
	533	wee(ij+u_ldown,2,2)hflux(ij+t_ldown+u_left,l)(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l)) + &
	534	wee(ij+u_ldown,3,2)hflux(ij+t_ldown+u_ldown,l)(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) + &
	535	wee(ij+u_ldown,4,2)hflux(ij+t_ldown+u_rdown,l)(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l)) + &
	536	wee(ij+u_ldown,5,2)hflux(ij+t_ldown+u_right,l)(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_right,l))
[369]	537	du(ij+u_right,l) = du(ij+u_right,l) + .5*uu_right
	538	du(ij+u_lup,l) = du(ij+u_lup,l) + .5*uu_lup
	539	du(ij+u_ldown,l) = du(ij+u_ldown,l) + .5*uu_ldown
[362]	540	ENDDO
	541	ENDDO
	542
[735]	543	CASE(conserv_gassmann) ! energy-conserving TRiSK modified by Gassmann (2018)
	544
	545	DO l=ll_begin,ll_end
	546	!DIR$ SIMD
	547	DO ij=ij_begin,ij_end
	548	uu_right = &
	549	wee(ij+u_right,1,1)hflux(ij+u_rup,l) qu(ij+t_right+u_lup,l)+ &
	550	wee(ij+u_right,2,1)hflux(ij+u_lup,l) qu(ij+u_rup,l)+ &
	551	.5wee(ij+u_right,3,1)hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+ &
	552	wee(ij+u_right,4,1)hflux(ij+u_ldown,l)qu(ij+u_rdown,l)+ &
	553	wee(ij+u_right,5,1)hflux(ij+u_rdown,l)qu(ij+t_right+u_ldown,l)+ &
	554	wee(ij+u_right,1,2)hflux(ij+t_right+u_ldown,l)qu(ij+u_rdown,l)+ &
	555	wee(ij+u_right,2,2)hflux(ij+t_right+u_rdown,l)qu(ij+t_right+u_ldown,l)+ &
	556	.5wee(ij+u_right,3,2)hflux(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+ &
	557	wee(ij+u_right,4,2)hflux(ij+t_right+u_rup,l)qu(ij+t_right+u_lup,l)+ &
	558	wee(ij+u_right,5,2)hflux(ij+t_right+u_lup,l)qu(ij+u_rup,l)
	559	uu_lup = &
	560	wee(ij+u_lup,1,1)hflux(ij+u_left,l)qu(ij+t_lup+u_ldown,l) + &
	561	wee(ij+u_lup,2,1)hflux(ij+u_ldown,l)qu(ij+u_left,l) + &
	562	.5wee(ij+u_lup,3,1)hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) + &
	563	wee(ij+u_lup,4,1)hflux(ij+u_right,l)qu(ij+u_rup,l) + &
	564	wee(ij+u_lup,5,1)hflux(ij+u_rup,l)qu(ij+t_lup+u_right,l)+ &
	565	wee(ij+u_lup,1,2)hflux(ij+t_lup+u_right,l)qu(ij+u_rup,l) + &
	566	wee(ij+u_lup,2,2)hflux(ij+t_lup+u_rup,l)qu(ij+t_lup+u_right,l) + &
	567	.5wee(ij+u_lup,3,2)hflux(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) + &
	568	wee(ij+u_lup,4,2)hflux(ij+t_lup+u_left,l)qu(ij+t_lup+u_ldown,l) + &
	569	wee(ij+u_lup,5,2)hflux(ij+t_lup+u_ldown,l)qu(ij+u_left,l)
	570	uu_ldown = &
	571	wee(ij+u_ldown,1,1)hflux(ij+u_rdown,l)qu(ij+t_ldown,l+u_right) + &
	572	wee(ij+u_ldown,2,1)hflux(ij+u_right,l)qu(ij+u_rdown,l) + &
	573	.5wee(ij+u_ldown,3,1)hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) + &
	574	wee(ij+u_ldown,4,1)hflux(ij+u_lup,l)qu(ij+u_left,l) + &
	575	wee(ij+u_ldown,5,1)hflux(ij+u_left,l)qu(ij+t_ldown+u_lup,l) + &
	576	wee(ij+u_ldown,1,2)hflux(ij+t_ldown+u_lup,l)qu(ij+u_left,l) + &
	577	wee(ij+u_ldown,2,2)hflux(ij+t_ldown+u_left,l)qu(ij+t_ldown+u_lup,l) + &
	578	.5wee(ij+u_ldown,3,2)hflux(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) + &
	579	wee(ij+u_ldown,4,2)hflux(ij+t_ldown+u_rdown,l)qu(ij+t_ldown+u_right,l) + &
	580	wee(ij+u_ldown,5,2)hflux(ij+t_ldown+u_right,l)qu(ij+u_rdown,l)
	581	du(ij+u_right,l) = du(ij+u_right,l) + uu_right
	582	du(ij+u_lup,l) = du(ij+u_lup,l) + uu_lup
	583	du(ij+u_ldown,l) = du(ij+u_ldown,l) + uu_ldown
	584	ENDDO
	585	ENDDO
	586
[733]	587	CASE(conserv_enstrophy) ! enstrophy-conserving TRiSK
[362]	588
	589	DO l=ll_begin,ll_end
	590	!DIR$ SIMD
	591	DO ij=ij_begin,ij_end
	592	uu_right = &
	593	wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+ &
	594	wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+ &
	595	wee(ij+u_right,3,1)*hflux(ij+u_left,l)+ &
	596	wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+ &
	597	wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+ &
	598	wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+ &
	599	wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+ &
	600	wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+ &
	601	wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+ &
	602	wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)
	603	uu_lup = &
	604	wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+ &
	605	wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+ &
	606	wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+ &
	607	wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+ &
	608	wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+ &
	609	wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+ &
	610	wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+ &
	611	wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+ &
	612	wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+ &
	613	wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)
	614	uu_ldown = &
	615	wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+ &
	616	wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+ &
	617	wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+ &
	618	wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+ &
	619	wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+ &
	620	wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+ &
	621	wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+ &
	622	wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+ &
	623	wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+ &
	624	wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)
	625
[734]	626	du(ij+u_right,l) = du(ij+u_right,l) + uu_right*qu(ij+u_right,l)
	627	du(ij+u_lup,l) = du(ij+u_lup,l) + uu_lup*qu(ij+u_lup,l)
	628	du(ij+u_ldown,l) = du(ij+u_ldown,l) + uu_ldown*qu(ij+u_ldown,l)
[369]	629	END DO
	630	END DO
[362]	631	CASE DEFAULT
	632	STOP
	633	END SELECT
[538]	634	#undef FTHETA
[362]	635
[562]	636	END IF ! dysl
	637
[369]	638	CALL trace_end("compute_caldyn_Coriolis")
	639
	640	END SUBROUTINE compute_caldyn_Coriolis
	641
[735]	642	SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hv, hflux,Kv,du, zero)
[529]	643	LOGICAL, INTENT(IN) :: zero
[369]	644	REAL(rstd),INTENT(IN) :: u(3iimjjm,llm) ! prognostic "velocity"
[734]	645	REAL(rstd),INTENT(IN) :: Kv(2iimjjm,llm) ! kinetic energy at vertices
[735]	646	REAL(rstd),INTENT(IN) :: hv(2iimjjm,llm) ! height/mass averaged to vertices
[369]	647	REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm)
	648	REAL(rstd),INTENT(OUT) :: hflux(3iimjjm,llm) ! hflux in kg/s
[529]	649	REAL(rstd),INTENT(INOUT) :: du(3iimjjm,llm)
[369]	650
[537]	651	REAL(rstd) :: berni(iim*jjm,llm) ! Bernoulli function
[573]	652	REAL(rstd) :: berni1(iim*jjm) ! Bernoulli function
[537]	653	REAL(rstd) :: uu_right, uu_lup, uu_ldown, ke, uu
[369]	654	INTEGER :: ij,l
	655
	656	CALL trace_start("compute_caldyn_slow_hydro")
	657
[562]	658	IF(dysl_slow_hydro) THEN
[573]	659
[537]	660	#define BERNI(ij,l) berni(ij,l)
[612]	661	#include "../kernels_hex/caldyn_slow_hydro.k90"
[538]	662	#undef BERNI
[562]	663
	664	ELSE
	665
[573]	666	#define BERNI(ij) berni1(ij)
[537]	667
[369]	668	DO l = ll_begin, ll_end
	669	! Compute mass fluxes
[733]	670	IF (caldyn_conserv==conserv_energy) CALL test_message(req_qu)
[735]	671
	672	IF(caldyn_kinetic==kinetic_trisk) THEN
	673	!DIR$ SIMD
	674	DO ij=ij_begin_ext,ij_end_ext
	675	uu_right=0.5(rhodz(ij,l)+rhodz(ij+t_right,l))u(ij+u_right,l)
	676	uu_lup=0.5(rhodz(ij,l)+rhodz(ij+t_lup,l))u(ij+u_lup,l)
	677	uu_ldown=0.5(rhodz(ij,l)+rhodz(ij+t_ldown,l))u(ij+u_ldown,l)
	678	uu_right= uu_right*le_de(ij+u_right)
	679	uu_lup = uu_lup *le_de(ij+u_lup)
	680	uu_ldown= uu_ldown*le_de(ij+u_ldown)
	681	hflux(ij+u_right,l)=uu_right
	682	hflux(ij+u_lup,l) =uu_lup
	683	hflux(ij+u_ldown,l)=uu_ldown
	684	ENDDO
	685	ELSE ! mass flux deriving from consistent kinetic energy
	686	!DIR$ SIMD
	687	DO ij=ij_begin_ext,ij_end_ext
	688	uu_right=0.5(hv(ij+z_rup,l)+hv(ij+z_rdown,l))u(ij+u_right,l)
	689	uu_lup=0.5(hv(ij+z_up,l)+hv(ij+z_lup,l))u(ij+u_lup,l)
	690	uu_ldown=0.5(hv(ij+z_ldown,l)+hv(ij+z_down,l))u(ij+u_ldown,l)
	691	uu_right= uu_right*le_de(ij+u_right)
	692	uu_lup = uu_lup *le_de(ij+u_lup)
	693	uu_ldown= uu_ldown*le_de(ij+u_ldown)
	694	hflux(ij+u_right,l)=uu_right
	695	hflux(ij+u_lup,l) =uu_lup
	696	hflux(ij+u_ldown,l)=uu_ldown
	697	ENDDO
	698	END IF
	699
[369]	700	! Compute Bernoulli=kinetic energy
[734]	701	IF(caldyn_kinetic==kinetic_trisk) THEN
	702	!DIR$ SIMD
	703	DO ij=ij_begin,ij_end
	704	BERNI(ij) = &
	705	1/(4Ai(ij))(le_de(ij+u_right)u(ij+u_right,l)*2 + &
	706	le_de(ij+u_rup)u(ij+u_rup,l)*2 + &
	707	le_de(ij+u_lup)u(ij+u_lup,l)*2 + &
	708	le_de(ij+u_left)u(ij+u_left,l)*2 + &
	709	le_de(ij+u_ldown)u(ij+u_ldown,l)*2 + &
	710	le_de(ij+u_rdown)u(ij+u_rdown,l)*2 )
	711	ENDDO
	712	ELSE
	713	!DIR$ SIMD
	714	DO ij=ij_begin,ij_end
	715	BERNI(ij) = Riv(ij,vup) *Kv(ij+z_up,l) + &
	716	Riv(ij,vlup) *Kv(ij+z_lup,l) + &
	717	Riv(ij,vldown)*Kv(ij+z_ldown,l) + &
	718	Riv(ij,vdown) *Kv(ij+z_down,l) + &
	719	Riv(ij,vrdown)*Kv(ij+z_rdown,l) + &
	720	Riv(ij,vrup) *Kv(ij+z_rup,l)
	721	END DO
	722	END IF
[375]	723	! Compute du=-grad(Bernoulli)
[529]	724	IF(zero) THEN
	725	!DIR$ SIMD
	726	DO ij=ij_begin,ij_end
[537]	727	du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right))
	728	du(ij+u_lup,l) = ne_lup*(BERNI(ij)-BERNI(ij+t_lup))
	729	du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))
[529]	730	END DO
	731	ELSE
	732	!DIR$ SIMD
	733	DO ij=ij_begin,ij_end
	734	du(ij+u_right,l) = du(ij+u_right,l) + &
[537]	735	ne_right*(BERNI(ij)-BERNI(ij+t_right))
[529]	736	du(ij+u_lup,l) = du(ij+u_lup,l) + &
[537]	737	ne_lup*(BERNI(ij)-BERNI(ij+t_lup))
[529]	738	du(ij+u_ldown,l) = du(ij+u_ldown,l) + &
[537]	739	ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))
[529]	740	END DO
	741	END IF
[369]	742	END DO
[573]	743
[538]	744	#undef BERNI
[573]	745
[562]	746	END IF ! dysl
[369]	747	CALL trace_end("compute_caldyn_slow_hydro")
	748	END SUBROUTINE compute_caldyn_slow_hydro
[362]	749
[558]	750	SUBROUTINE compute_caldyn_slow_NH(u,rhodz,Phi,W, F_el,gradPhi2,w_il, hflux,du,dPhi,dW)
[369]	751	REAL(rstd),INTENT(IN) :: u(3iimjjm,llm) ! prognostic "velocity"
	752	REAL(rstd),INTENT(IN) :: rhodz(iimjjm,llm) ! rhodz
	753	REAL(rstd),INTENT(IN) :: Phi(iim*jjm,llm+1) ! prognostic geopotential
	754	REAL(rstd),INTENT(IN) :: W(iim*jjm,llm+1) ! prognostic vertical momentum
[362]	755
[369]	756	REAL(rstd),INTENT(OUT) :: hflux(3iimjjm,llm) ! hflux in kg/s
	757	REAL(rstd),INTENT(OUT) :: du(3iimjjm,llm)
	758	REAL(rstd),INTENT(OUT) :: dW(iim*jjm,llm+1)
	759	REAL(rstd),INTENT(OUT) :: dPhi(iim*jjm,llm+1)
	760
[558]	761	REAL(rstd) :: w_il(iim*jjm,llm+1) ! Wil/mil
[369]	762	REAL(rstd) :: F_el(3iimjjm,llm+1) ! NH mass flux
[558]	763	REAL(rstd) :: gradPhi2(iimjjm,llm+1) ! grad_Phi*2
[369]	764	REAL(rstd) :: DePhil(3iimjjm,llm+1) ! grad(Phi)
[539]	765
	766	INTEGER :: ij,l,kdown,kup
	767	REAL(rstd) :: W_el, W2_el, uu_right, uu_lup, uu_ldown, gPhi2, dP, divG, u2, uu
	768
	769	REAL(rstd) :: berni(iim*jjm,llm) ! Bernoulli function
	770	REAL(rstd) :: G_el(3iimjjm,llm+1) ! horizontal flux of W
	771	REAL(rstd) :: v_el(3iimjjm,llm+1)
[369]	772
[573]	773	REAL(rstd) :: berni1(iim*jjm) ! Bernoulli function
	774	REAL(rstd) :: G_el1(3iimjjm) ! horizontal flux of W
	775	REAL(rstd) :: v_el1(3iimjjm)
	776
[369]	777	CALL trace_start("compute_caldyn_slow_NH")
	778
[573]	779	IF(dysl) THEN
	780
[558]	781	!$OMP BARRIER
[612]	782	#include "../kernels_hex/caldyn_slow_NH.k90"
[558]	783	!$OMP BARRIER
[573]	784
	785	ELSE
	786
	787	#define BERNI(ij) berni1(ij)
	788	#define G_EL(ij) G_el1(ij)
	789	#define V_EL(ij) v_el1(ij)
	790
[369]	791	DO l=ll_begin, ll_endp1 ! compute on l levels (interfaces)
	792	IF(l==1) THEN
	793	kdown=1
	794	ELSE
	795	kdown=l-1
	796	END IF
	797	IF(l==llm+1) THEN
	798	kup=llm
	799	ELSE
	800	kup=l
	801	END IF
[377]	802	! below : "checked" means "formula also valid when kup=kdown (top/bottom)"
[369]	803	! compute mil, wil=Wil/mil
	804	DO ij=ij_begin_ext, ij_end_ext
[377]	805	w_il(ij,l) = 2.*W(ij,l)/(rhodz(ij,kdown)+rhodz(ij,kup)) ! checked
[369]	806	END DO
	807	! compute DePhi, v_el, G_el, F_el
	808	! v_el, W2_el and therefore G_el incorporate metric factor le_de
	809	! while DePhil, W_el and F_el don't
	810	DO ij=ij_begin_ext, ij_end_ext
	811	! Compute on edge 'right'
	812	W_el = .5*( W(ij,l)+W(ij+t_right,l) )
	813	DePhil(ij+u_right,l) = ne_right*(Phi(ij+t_right,l)-Phi(ij,l))
	814	F_el(ij+u_right,l) = DePhil(ij+u_right,l)*W_el
	815	W2_el = .5le_de(ij+u_right) &
	816	( W(ij,l)w_il(ij,l) + W(ij+t_right,l)w_il(ij+t_right,l) )
[573]	817	V_EL(ij+u_right) = .5le_de(ij+u_right)(u(ij+u_right,kup)+u(ij+u_right,kdown)) ! checked
	818	G_EL(ij+u_right) = V_EL(ij+u_right)W_el - DePhil(ij+u_right,l)W2_el
[369]	819	! Compute on edge 'lup'
	820	W_el = .5*( W(ij,l)+W(ij+t_lup,l) )
	821	DePhil(ij+u_lup,l) = ne_lup*(Phi(ij+t_lup,l)-Phi(ij,l))
	822	F_el(ij+u_lup,l) = DePhil(ij+u_lup,l)*W_el
	823	W2_el = .5le_de(ij+u_lup) &
	824	( W(ij,l)w_il(ij,l) + W(ij+t_lup,l)w_il(ij+t_lup,l) )
[573]	825	V_EL(ij+u_lup) = .5le_de(ij+u_lup)( u(ij+u_lup,kup) + u(ij+u_lup,kdown)) ! checked
	826	G_EL(ij+u_lup) = V_EL(ij+u_lup)W_el - DePhil(ij+u_lup,l)W2_el
[369]	827	! Compute on edge 'ldown'
	828	W_el = .5*( W(ij,l)+W(ij+t_ldown,l) )
	829	DePhil(ij+u_ldown,l) = ne_ldown*(Phi(ij+t_ldown,l)-Phi(ij,l))
	830	F_el(ij+u_ldown,l) = DePhil(ij+u_ldown,l)*W_el
	831	W2_el = .5le_de(ij+u_ldown) &
	832	( W(ij,l)w_il(ij,l) + W(ij+t_ldown,l)w_il(ij+t_ldown,l) )
[573]	833	V_EL(ij+u_ldown) = .5le_de(ij+u_ldown)( u(ij+u_ldown,kup) + u(ij+u_ldown,kdown)) ! checked
	834	G_EL(ij+u_ldown) = V_EL(ij+u_ldown)W_el - DePhil(ij+u_ldown,l)W2_el
[369]	835	END DO
	836	! compute GradPhi2, dPhi, dW
	837	DO ij=ij_begin_ext, ij_end_ext
	838	gradPhi2(ij,l) = &
	839	1/(2Ai(ij))(le_de(ij+u_right)DePhil(ij+u_right,l)*2 + &
	840	le_de(ij+u_rup)DePhil(ij+u_rup,l)*2 + &
	841	le_de(ij+u_lup)DePhil(ij+u_lup,l)*2 + &
	842	le_de(ij+u_left)DePhil(ij+u_left,l)*2 + &
	843	le_de(ij+u_ldown)DePhil(ij+u_ldown,l)*2 + &
	844	le_de(ij+u_rdown)DePhil(ij+u_rdown,l)*2 )
[377]	845
	846	dPhi(ij,l) = gradPhi2(ij,l)w_il(ij,l) -1/(2Ai(ij))* &
[573]	847	( DePhil(ij+u_right,l)*V_EL(ij+u_right) + & ! -v.gradPhi,
	848	DePhil(ij+u_rup,l)*V_EL(ij+u_rup) + & ! v_el already has le_de
	849	DePhil(ij+u_lup,l)*V_EL(ij+u_lup) + &
	850	DePhil(ij+u_left,l)*V_EL(ij+u_left) + &
	851	DePhil(ij+u_ldown,l)*V_EL(ij+u_ldown) + &
	852	DePhil(ij+u_rdown,l)*V_EL(ij+u_rdown) )
[377]	853
[369]	854	dW(ij,l) = -1./Ai(ij)*( & ! -div(G_el),
[573]	855	ne_right*G_EL(ij+u_right) + & ! G_el already has le_de
	856	ne_rup*G_EL(ij+u_rup) + &
	857	ne_lup*G_EL(ij+u_lup) + &
	858	ne_left*G_EL(ij+u_left) + &
	859	ne_ldown*G_EL(ij+u_ldown) + &
	860	ne_rdown*G_EL(ij+u_rdown))
[369]	861	END DO
	862	END DO
[377]	863
[369]	864	DO l=ll_begin, ll_end ! compute on k levels (layers)
	865	! Compute berni at scalar points
	866	DO ij=ij_begin_ext, ij_end_ext
[573]	867	BERNI(ij) = &
[369]	868	1/(4Ai(ij))( &
	869	le_de(ij+u_right)u(ij+u_right,l)*2 + &
	870	le_de(ij+u_rup)u(ij+u_rup,l)*2 + &
	871	le_de(ij+u_lup)u(ij+u_lup,l)*2 + &
	872	le_de(ij+u_left)u(ij+u_left,l)*2 + &
	873	le_de(ij+u_ldown)u(ij+u_ldown,l)*2 + &
	874	le_de(ij+u_rdown)u(ij+u_rdown,l)*2 ) &
[499]	875	- .25( gradPhi2(ij,l) w_il(ij,l)**2 + &
[369]	876	gradPhi2(ij,l+1)w_il(ij,l+1)*2 )
	877	END DO
	878	! Compute mass flux and grad(berni) at edges
	879	DO ij=ij_begin_ext, ij_end_ext
	880	! Compute on edge 'right'
	881	uu_right = 0.5(rhodz(ij,l)+rhodz(ij+t_right,l))u(ij+u_right,l) &
	882	-0.5*(F_el(ij+u_right,l)+F_el(ij+u_right,l+1))
	883	hflux(ij+u_right,l) = uu_right*le_de(ij+u_right)
[573]	884	du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right))
[369]	885	! Compute on edge 'lup'
	886	uu_lup = 0.5(rhodz(ij,l)+rhodz(ij+t_lup,l))u(ij+u_lup,l) &
	887	-0.5*(F_el(ij+u_lup,l)+F_el(ij+u_lup,l+1))
	888	hflux(ij+u_lup,l) = uu_lup*le_de(ij+u_lup)
[573]	889	du(ij+u_lup,l) = ne_lup*(BERNI(ij)-BERNI(ij+t_lup))
[369]	890	! Compute on edge 'ldown'
	891	uu_ldown = 0.5(rhodz(ij,l)+rhodz(ij+t_ldown,l))u(ij+u_ldown,l) &
	892	-0.5*(F_el(ij+u_ldown,l)+F_el(ij+u_ldown,l+1))
	893	hflux(ij+u_ldown,l) = uu_ldown*le_de(ij+u_ldown)
[573]	894	du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))
[369]	895	END DO
	896	END DO
	897
[573]	898	#undef V_EL
	899	#undef G_EL
	900	#undef BERNI
	901
	902	END IF ! dysl
	903
[369]	904	CALL trace_end("compute_caldyn_slow_NH")
	905
	906	END SUBROUTINE compute_caldyn_slow_NH
	907
[362]	908	END MODULE caldyn_kernels_hevi_mod

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

Download in other formats: