source: codes/icosagcm/devel/src/kernels_unst/caldyn_vert_NH.k90 @ 624

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

devel : added unstructured-mesh kernels
File size: 2.1 KB
Line 
1   !--------------------------------------------------------------------------
2   !---------------------------- caldyn_vert_NH ----------------------------------
3   !$OMP DO SCHEDULE(STATIC)
4   DO ij = 1, primal_num
5      l=1
6      w_ij = ( W(l,ij)+.5*W(l+1,ij) )/mass(l,ij)
7      wflux_ij = .5*(wflux(l,ij)+wflux(l+1,ij))
8      W_etadot(l,ij) = wflux_ij*w_ij
9      eta_dot(l,ij) = wflux_ij / mass(l,ij)
10      wcov(l,ij) = w_ij*(geopot(l+1,ij)-geopot(l,ij))
11      DO l = 2, llm-1
12         w_ij = .5*( W(l,ij)+W(l+1,ij) )/mass(l,ij)
13         wflux_ij = .5*(wflux(l,ij)+wflux(l+1,ij))
14         W_etadot(l,ij) = wflux_ij*w_ij
15         eta_dot(l,ij) = wflux_ij / mass(l,ij)
16         wcov(l,ij) = w_ij*(geopot(l+1,ij)-geopot(l,ij))
17      END DO
18      l=llm
19      w_ij = ( .5*W(l,ij)+W(l+1,ij) )/mass(l,ij)
20      wflux_ij = .5*(wflux(l,ij)+wflux(l+1,ij))
21      W_etadot(l,ij) = wflux_ij*w_ij
22      eta_dot(l,ij) = wflux_ij / mass(l,ij)
23      wcov(l,ij) = w_ij*(geopot(l+1,ij)-geopot(l,ij))
24   END DO
25   !$OMP END DO
26   ! add NH term to du
27   !$OMP DO SCHEDULE(STATIC)
28   DO edge = 1, edge_num
29      ij_left = left(edge)
30      ij_right = right(edge)
31      DO l = 1, llm
32         du(l,edge) = du(l,edge) - .5*(wcov(l,ij_left)+wcov(l,ij_right))*1.*(eta_dot(l,ij_right)-eta_dot(l,ij_left))
33      END DO
34   END DO
35   !$OMP END DO
36   ! add NH terms to dW, dPhi
37   ! FIXME : TODO top and bottom
38   !$OMP DO SCHEDULE(STATIC)
39   DO ij = 1, primal_num
40      DO l = 2, llm
41         dPhi(l,ij)=dPhi(l,ij)-wflux(l,ij)*(geopot(l+1,ij)-geopot(l-1,ij))/(mass(l-1,ij)+mass(l,ij))
42      END DO
43   END DO
44   !$OMP END DO
45   ! We need a barrier here because we compute W_etadot above and do a vertical difference below
46   !$OMP BARRIER
47   !$OMP DO SCHEDULE(STATIC)
48   DO ij = 1, primal_num
49      l=1
50      dW(l,ij) = dW(l,ij) - W_etadot(l,ij)
51      DO l = 2, llm
52         dW(l,ij) = dW(l,ij) + W_etadot(l-1,ij) - W_etadot(l,ij)
53      END DO
54      l=llm+1
55      dW(l,ij) = dW(l,ij) + W_etadot(l-1,ij)
56   END DO
57   !$OMP END DO
58   !---------------------------- caldyn_vert_NH ----------------------------------
59   !--------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.