1 | !-------------------------------------------------------------------------- |
---|
2 | !---------------------------- caldyn_slow_NH ---------------------------------- |
---|
3 | !$OMP DO SCHEDULE(STATIC) |
---|
4 | DO ij = 1, primal_num |
---|
5 | kdown = 1 |
---|
6 | kup = 1 |
---|
7 | l=1 |
---|
8 | w_il(l,ij) = 2.*W(l,ij)/rhodz(kup,ij) |
---|
9 | DO l = 2, llm |
---|
10 | kdown = l-1 |
---|
11 | kup = l |
---|
12 | w_il(l,ij) = 2.*W(l,ij)/(rhodz(kdown,ij)+rhodz(kup,ij)) |
---|
13 | END DO |
---|
14 | kdown = llm |
---|
15 | kup = llm |
---|
16 | l=llm+1 |
---|
17 | w_il(l,ij) = 2.*W(l,ij)/rhodz(kdown,ij) |
---|
18 | END DO |
---|
19 | !$OMP END DO |
---|
20 | !$OMP DO SCHEDULE(STATIC) |
---|
21 | DO edge = 1, edge_num |
---|
22 | ij_left = left(edge) |
---|
23 | ij_right = right(edge) |
---|
24 | kdown = 1 |
---|
25 | kup = 1 |
---|
26 | l=1 |
---|
27 | ! compute DePhi, v_el, G_el, F_el |
---|
28 | ! v_el, W2_el and therefore G_el incorporate metric factor le_de |
---|
29 | ! while DePhil, W_el and F_el do not |
---|
30 | W_el = .5*( W(l,ij_right)+W(l,ij_left) ) |
---|
31 | DePhil(l,edge) = 1.*(Phi(l,ij_right)-Phi(l,ij_left)) |
---|
32 | F_el(l,edge) = DePhil(l,edge)*W_el |
---|
33 | W2_el = .5*le_de(edge) * ( W(l,ij_left)*w_il(l,ij_left) + W(l,ij_right)*w_il(l,ij_right) ) |
---|
34 | v_el(l,edge) = .5*le_de(edge)*(u(kup,edge)+u(kdown,edge)) ! checked |
---|
35 | G_el(l,edge) = v_el(l,edge)*W_el - DePhil(l,edge)*W2_el |
---|
36 | DO l = 2, llm |
---|
37 | kdown = l-1 |
---|
38 | kup = l |
---|
39 | ! compute DePhi, v_el, G_el, F_el |
---|
40 | ! v_el, W2_el and therefore G_el incorporate metric factor le_de |
---|
41 | ! while DePhil, W_el and F_el do not |
---|
42 | W_el = .5*( W(l,ij_right)+W(l,ij_left) ) |
---|
43 | DePhil(l,edge) = 1.*(Phi(l,ij_right)-Phi(l,ij_left)) |
---|
44 | F_el(l,edge) = DePhil(l,edge)*W_el |
---|
45 | W2_el = .5*le_de(edge) * ( W(l,ij_left)*w_il(l,ij_left) + W(l,ij_right)*w_il(l,ij_right) ) |
---|
46 | v_el(l,edge) = .5*le_de(edge)*(u(kup,edge)+u(kdown,edge)) ! checked |
---|
47 | G_el(l,edge) = v_el(l,edge)*W_el - DePhil(l,edge)*W2_el |
---|
48 | END DO |
---|
49 | kdown = llm |
---|
50 | kup = llm |
---|
51 | l=llm+1 |
---|
52 | ! compute DePhi, v_el, G_el, F_el |
---|
53 | ! v_el, W2_el and therefore G_el incorporate metric factor le_de |
---|
54 | ! while DePhil, W_el and F_el do not |
---|
55 | W_el = .5*( W(l,ij_right)+W(l,ij_left) ) |
---|
56 | DePhil(l,edge) = 1.*(Phi(l,ij_right)-Phi(l,ij_left)) |
---|
57 | F_el(l,edge) = DePhil(l,edge)*W_el |
---|
58 | W2_el = .5*le_de(edge) * ( W(l,ij_left)*w_il(l,ij_left) + W(l,ij_right)*w_il(l,ij_right) ) |
---|
59 | v_el(l,edge) = .5*le_de(edge)*(u(kup,edge)+u(kdown,edge)) ! checked |
---|
60 | G_el(l,edge) = v_el(l,edge)*W_el - DePhil(l,edge)*W2_el |
---|
61 | END DO |
---|
62 | !$OMP END DO |
---|
63 | ! compute GradPhi2, dPhi, dW |
---|
64 | !$OMP DO SCHEDULE(STATIC) |
---|
65 | DO ij = 1, primal_num |
---|
66 | DO l = 1, llm+1 |
---|
67 | gPhi2=0. |
---|
68 | dP=0. |
---|
69 | divG=0 |
---|
70 | DO iedge = 1, primal_deg(ij) |
---|
71 | edge = primal_edge(iedge,ij) |
---|
72 | gPhi2 = gPhi2 + le_de(edge)*DePhil(l,edge)**2 |
---|
73 | dP = dP + le_de(edge)*DePhil(l,edge)*v_el(l,edge) |
---|
74 | divG = divG + primal_ne(iedge,ij)*G_el(l,edge) ! -div(G_el), G_el already has le_de |
---|
75 | END DO |
---|
76 | gradPhi2(l,ij) = 1./(2.*Ai(ij)) * gPhi2 |
---|
77 | dPhi(l,ij) = gradPhi2(l,ij)*w_il(l,ij) - 1./(2.*Ai(ij))*dP |
---|
78 | dW(l,ij) = (-1./Ai(ij))*divG |
---|
79 | END DO |
---|
80 | END DO |
---|
81 | !$OMP END DO |
---|
82 | ! We need a barrier here because we compute gradPhi2, F_el and w_il above and do a vertical average below |
---|
83 | !$OMP BARRIER |
---|
84 | ! Compute berni at scalar points |
---|
85 | !$OMP DO SCHEDULE(STATIC) |
---|
86 | DO ij = 1, primal_num |
---|
87 | DO l = 1, llm |
---|
88 | u2=0. |
---|
89 | DO iedge = 1, primal_deg(ij) |
---|
90 | edge = primal_edge(iedge,ij) |
---|
91 | u2 = u2 + le_de(edge)*u(l,edge)**2 |
---|
92 | END DO |
---|
93 | 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 ) |
---|
94 | END DO |
---|
95 | END DO |
---|
96 | !$OMP END DO |
---|
97 | !$OMP DO SCHEDULE(STATIC) |
---|
98 | DO edge = 1, edge_num |
---|
99 | ij_left = left(edge) |
---|
100 | ij_right = right(edge) |
---|
101 | DO l = 1, llm |
---|
102 | ! Compute mass flux and grad(berni) |
---|
103 | uu = .5*(rhodz(l,ij_left)+rhodz(l,ij_right))*u(l,edge) - .5*( F_el(l,edge)+F_el(l+1,edge) ) |
---|
104 | hflux(l,edge) = le_de(edge)*uu |
---|
105 | du(l,edge) = 1.*(berni(l,ij_left)-berni(l,ij_right)) |
---|
106 | END DO |
---|
107 | END DO |
---|
108 | !$OMP END DO |
---|
109 | !---------------------------- caldyn_slow_NH ---------------------------------- |
---|
110 | !-------------------------------------------------------------------------- |
---|