Context Navigation

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

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

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

Download in other formats: