!-------------------------------------------------------------------------- !---------------------------- caldyn_slow_NH ---------------------------------- !$OMP DO SCHEDULE(STATIC) DO ij = 1, primal_num kdown = 1 kup = 1 l=1 w_il(l,ij) = 2.*W(l,ij)/rhodz(kup,ij) DO l = 2, llm kdown = l-1 kup = l w_il(l,ij) = 2.*W(l,ij)/(rhodz(kdown,ij)+rhodz(kup,ij)) END DO kdown = llm kup = llm l=llm+1 w_il(l,ij) = 2.*W(l,ij)/rhodz(kdown,ij) END DO !$OMP END DO !$OMP DO SCHEDULE(STATIC) DO edge = 1, edge_num ij_left = left(edge) ij_right = right(edge) kdown = 1 kup = 1 l=1 ! compute DePhi, v_el, G_el, F_el ! v_el, W2_el and therefore G_el incorporate metric factor le_de ! while DePhil, W_el and F_el do not W_el = .5*( W(l,ij_right)+W(l,ij_left) ) DePhil(l,edge) = 1.*(Phi(l,ij_right)-Phi(l,ij_left)) F_el(l,edge) = DePhil(l,edge)*W_el W2_el = .5*le_de(edge) * ( W(l,ij_left)*w_il(l,ij_left) + W(l,ij_right)*w_il(l,ij_right) ) v_el(l,edge) = .5*le_de(edge)*(u(kup,edge)+u(kdown,edge)) ! checked G_el(l,edge) = v_el(l,edge)*W_el - DePhil(l,edge)*W2_el DO l = 2, llm kdown = l-1 kup = l ! compute DePhi, v_el, G_el, F_el ! v_el, W2_el and therefore G_el incorporate metric factor le_de ! while DePhil, W_el and F_el do not W_el = .5*( W(l,ij_right)+W(l,ij_left) ) DePhil(l,edge) = 1.*(Phi(l,ij_right)-Phi(l,ij_left)) F_el(l,edge) = DePhil(l,edge)*W_el W2_el = .5*le_de(edge) * ( W(l,ij_left)*w_il(l,ij_left) + W(l,ij_right)*w_il(l,ij_right) ) v_el(l,edge) = .5*le_de(edge)*(u(kup,edge)+u(kdown,edge)) ! checked G_el(l,edge) = v_el(l,edge)*W_el - DePhil(l,edge)*W2_el END DO kdown = llm kup = llm l=llm+1 ! compute DePhi, v_el, G_el, F_el ! v_el, W2_el and therefore G_el incorporate metric factor le_de ! while DePhil, W_el and F_el do not W_el = .5*( W(l,ij_right)+W(l,ij_left) ) DePhil(l,edge) = 1.*(Phi(l,ij_right)-Phi(l,ij_left)) F_el(l,edge) = DePhil(l,edge)*W_el W2_el = .5*le_de(edge) * ( W(l,ij_left)*w_il(l,ij_left) + W(l,ij_right)*w_il(l,ij_right) ) v_el(l,edge) = .5*le_de(edge)*(u(kup,edge)+u(kdown,edge)) ! checked G_el(l,edge) = v_el(l,edge)*W_el - DePhil(l,edge)*W2_el END DO !$OMP END DO ! compute GradPhi2, dPhi, dW !$OMP DO SCHEDULE(STATIC) DO ij = 1, primal_num DO l = 1, llm+1 gPhi2=0. dP=0. divG=0 DO iedge = 1, primal_deg(ij) edge = primal_edge(iedge,ij) gPhi2 = gPhi2 + le_de(edge)*DePhil(l,edge)**2 dP = dP + le_de(edge)*DePhil(l,edge)*v_el(l,edge) divG = divG + primal_ne(iedge,ij)*G_el(l,edge) ! -div(G_el), G_el already has le_de END DO gradPhi2(l,ij) = 1./(2.*Ai(ij)) * gPhi2 dPhi(l,ij) = gradPhi2(l,ij)*w_il(l,ij) - 1./(2.*Ai(ij))*dP dW(l,ij) = (-1./Ai(ij))*divG END DO END DO !$OMP END DO ! We need a barrier here because we compute gradPhi2, F_el and w_il above and do a vertical average below !$OMP BARRIER ! Compute berni at scalar points !$OMP DO SCHEDULE(STATIC) DO ij = 1, primal_num DO l = 1, llm u2=0. DO iedge = 1, primal_deg(ij) edge = primal_edge(iedge,ij) u2 = u2 + le_de(edge)*u(l,edge)**2 END DO berni(l,ij) = 1./(4.*Ai(ij)) * u2 - .25*( gradPhi2(l,ij)*w_il(l,ij)**2 + gradPhi2(l+1,ij)*w_il(l+1,ij)**2 ) END DO END DO !$OMP END DO !$OMP DO SCHEDULE(STATIC) DO edge = 1, edge_num ij_left = left(edge) ij_right = right(edge) DO l = 1, llm ! Compute mass flux and grad(berni) uu = .5*(rhodz(l,ij_left)+rhodz(l,ij_right))*u(l,edge) - .5*( F_el(l,edge)+F_el(l+1,edge) ) hflux(l,edge) = le_de(edge)*uu du(l,edge) = 1.*(berni(l,ij_left)-berni(l,ij_right)) END DO END DO !$OMP END DO !---------------------------- caldyn_slow_NH ---------------------------------- !--------------------------------------------------------------------------