source: codes/icosagcm/devel/src/kernels_unst/caldyn_slow_NH.k90 @ 614

Last change on this file since 614 was 614, checked in by dubos, 7 years ago

devel : added unstructured-mesh kernels
File size: 4.2 KB
Line 
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   !--------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.