Context Navigation

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

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

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

Download in other formats: