source: codes/icosagcm/trunk/src/kernels/caldyn_slow_NH.k90 @ 580

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

trunk : upgrading to devel
File size: 9.9 KB
Line 
1   !--------------------------------------------------------------------------
2   !---------------------------- caldyn_slow_NH ----------------------------------
3   IF (ll_begin==1) THEN
4      !DIR$ SIMD
5      DO ij=ij_begin_ext, ij_end_ext
6         w_il(ij,1) = 2.*W(ij,1)/rhodz(ij,1)
7      END DO
8   END IF
9   DO l = ll_beginp1, ll_end
10      !DIR$ SIMD
11      DO ij=ij_begin_ext, ij_end_ext
12         w_il(ij,l) = 2.*W(ij,l)/(rhodz(ij,l-1)+rhodz(ij,l))
13      END DO
14   END DO
15   IF(ll_endp1==llm+1) THEN
16      !DIR$ SIMD
17      DO ij=ij_begin_ext, ij_end_ext
18         w_il(ij,llm+1) = 2.*W(ij,llm+1)/rhodz(ij,llm)
19      END DO
20   END IF
21   IF (ll_begin==1) THEN
22      !DIR$ SIMD
23      DO ij=ij_begin_ext, ij_end_ext
24         ! compute DePhi, v_el, G_el, F_el
25         ! v_el, W2_el and therefore G_el incorporate metric factor le_de
26         ! while DePhil, W_el and F_el do not
27         W_el = .5*( W(ij+t_right,1)+W(ij,1) )
28         DePhil(ij+u_right,1) = ne_right*(Phi(ij+t_right,1)-Phi(ij,1))
29         F_el(ij+u_right,1) = DePhil(ij+u_right,1)*W_el
30         W2_el = .5*le_de(ij+u_right) * ( W(ij,1)*w_il(ij,1) + W(ij+t_right,1)*w_il(ij+t_right,1) )
31         v_el(ij+u_right,1) = .5*le_de(ij+u_right)*(u(ij+u_right,1)+u(ij+u_right,1)) ! checked
32         G_el(ij+u_right,1) = v_el(ij+u_right,1)*W_el - DePhil(ij+u_right,1)*W2_el
33         ! compute DePhi, v_el, G_el, F_el
34         ! v_el, W2_el and therefore G_el incorporate metric factor le_de
35         ! while DePhil, W_el and F_el do not
36         W_el = .5*( W(ij+t_lup,1)+W(ij,1) )
37         DePhil(ij+u_lup,1) = ne_lup*(Phi(ij+t_lup,1)-Phi(ij,1))
38         F_el(ij+u_lup,1) = DePhil(ij+u_lup,1)*W_el
39         W2_el = .5*le_de(ij+u_lup) * ( W(ij,1)*w_il(ij,1) + W(ij+t_lup,1)*w_il(ij+t_lup,1) )
40         v_el(ij+u_lup,1) = .5*le_de(ij+u_lup)*(u(ij+u_lup,1)+u(ij+u_lup,1)) ! checked
41         G_el(ij+u_lup,1) = v_el(ij+u_lup,1)*W_el - DePhil(ij+u_lup,1)*W2_el
42         ! compute DePhi, v_el, G_el, F_el
43         ! v_el, W2_el and therefore G_el incorporate metric factor le_de
44         ! while DePhil, W_el and F_el do not
45         W_el = .5*( W(ij+t_ldown,1)+W(ij,1) )
46         DePhil(ij+u_ldown,1) = ne_ldown*(Phi(ij+t_ldown,1)-Phi(ij,1))
47         F_el(ij+u_ldown,1) = DePhil(ij+u_ldown,1)*W_el
48         W2_el = .5*le_de(ij+u_ldown) * ( W(ij,1)*w_il(ij,1) + W(ij+t_ldown,1)*w_il(ij+t_ldown,1) )
49         v_el(ij+u_ldown,1) = .5*le_de(ij+u_ldown)*(u(ij+u_ldown,1)+u(ij+u_ldown,1)) ! checked
50         G_el(ij+u_ldown,1) = v_el(ij+u_ldown,1)*W_el - DePhil(ij+u_ldown,1)*W2_el
51      END DO
52   END IF
53   DO l = ll_beginp1, ll_end
54      !DIR$ SIMD
55      DO ij=ij_begin_ext, ij_end_ext
56         ! compute DePhi, v_el, G_el, F_el
57         ! v_el, W2_el and therefore G_el incorporate metric factor le_de
58         ! while DePhil, W_el and F_el do not
59         W_el = .5*( W(ij+t_right,l)+W(ij,l) )
60         DePhil(ij+u_right,l) = ne_right*(Phi(ij+t_right,l)-Phi(ij,l))
61         F_el(ij+u_right,l) = DePhil(ij+u_right,l)*W_el
62         W2_el = .5*le_de(ij+u_right) * ( W(ij,l)*w_il(ij,l) + W(ij+t_right,l)*w_il(ij+t_right,l) )
63         v_el(ij+u_right,l) = .5*le_de(ij+u_right)*(u(ij+u_right,l)+u(ij+u_right,l-1)) ! checked
64         G_el(ij+u_right,l) = v_el(ij+u_right,l)*W_el - DePhil(ij+u_right,l)*W2_el
65         ! compute DePhi, v_el, G_el, F_el
66         ! v_el, W2_el and therefore G_el incorporate metric factor le_de
67         ! while DePhil, W_el and F_el do not
68         W_el = .5*( W(ij+t_lup,l)+W(ij,l) )
69         DePhil(ij+u_lup,l) = ne_lup*(Phi(ij+t_lup,l)-Phi(ij,l))
70         F_el(ij+u_lup,l) = DePhil(ij+u_lup,l)*W_el
71         W2_el = .5*le_de(ij+u_lup) * ( W(ij,l)*w_il(ij,l) + W(ij+t_lup,l)*w_il(ij+t_lup,l) )
72         v_el(ij+u_lup,l) = .5*le_de(ij+u_lup)*(u(ij+u_lup,l)+u(ij+u_lup,l-1)) ! checked
73         G_el(ij+u_lup,l) = v_el(ij+u_lup,l)*W_el - DePhil(ij+u_lup,l)*W2_el
74         ! compute DePhi, v_el, G_el, F_el
75         ! v_el, W2_el and therefore G_el incorporate metric factor le_de
76         ! while DePhil, W_el and F_el do not
77         W_el = .5*( W(ij+t_ldown,l)+W(ij,l) )
78         DePhil(ij+u_ldown,l) = ne_ldown*(Phi(ij+t_ldown,l)-Phi(ij,l))
79         F_el(ij+u_ldown,l) = DePhil(ij+u_ldown,l)*W_el
80         W2_el = .5*le_de(ij+u_ldown) * ( W(ij,l)*w_il(ij,l) + W(ij+t_ldown,l)*w_il(ij+t_ldown,l) )
81         v_el(ij+u_ldown,l) = .5*le_de(ij+u_ldown)*(u(ij+u_ldown,l)+u(ij+u_ldown,l-1)) ! checked
82         G_el(ij+u_ldown,l) = v_el(ij+u_ldown,l)*W_el - DePhil(ij+u_ldown,l)*W2_el
83      END DO
84   END DO
85   IF(ll_endp1==llm+1) THEN
86      !DIR$ SIMD
87      DO ij=ij_begin_ext, ij_end_ext
88         ! compute DePhi, v_el, G_el, F_el
89         ! v_el, W2_el and therefore G_el incorporate metric factor le_de
90         ! while DePhil, W_el and F_el do not
91         W_el = .5*( W(ij+t_right,llm+1)+W(ij,llm+1) )
92         DePhil(ij+u_right,llm+1) = ne_right*(Phi(ij+t_right,llm+1)-Phi(ij,llm+1))
93         F_el(ij+u_right,llm+1) = DePhil(ij+u_right,llm+1)*W_el
94         W2_el = .5*le_de(ij+u_right) * ( W(ij,llm+1)*w_il(ij,llm+1) + W(ij+t_right,llm+1)*w_il(ij+t_right,llm+1) )
95         v_el(ij+u_right,llm+1) = .5*le_de(ij+u_right)*(u(ij+u_right,llm)+u(ij+u_right,llm)) ! checked
96         G_el(ij+u_right,llm+1) = v_el(ij+u_right,llm+1)*W_el - DePhil(ij+u_right,llm+1)*W2_el
97         ! compute DePhi, v_el, G_el, F_el
98         ! v_el, W2_el and therefore G_el incorporate metric factor le_de
99         ! while DePhil, W_el and F_el do not
100         W_el = .5*( W(ij+t_lup,llm+1)+W(ij,llm+1) )
101         DePhil(ij+u_lup,llm+1) = ne_lup*(Phi(ij+t_lup,llm+1)-Phi(ij,llm+1))
102         F_el(ij+u_lup,llm+1) = DePhil(ij+u_lup,llm+1)*W_el
103         W2_el = .5*le_de(ij+u_lup) * ( W(ij,llm+1)*w_il(ij,llm+1) + W(ij+t_lup,llm+1)*w_il(ij+t_lup,llm+1) )
104         v_el(ij+u_lup,llm+1) = .5*le_de(ij+u_lup)*(u(ij+u_lup,llm)+u(ij+u_lup,llm)) ! checked
105         G_el(ij+u_lup,llm+1) = v_el(ij+u_lup,llm+1)*W_el - DePhil(ij+u_lup,llm+1)*W2_el
106         ! compute DePhi, v_el, G_el, F_el
107         ! v_el, W2_el and therefore G_el incorporate metric factor le_de
108         ! while DePhil, W_el and F_el do not
109         W_el = .5*( W(ij+t_ldown,llm+1)+W(ij,llm+1) )
110         DePhil(ij+u_ldown,llm+1) = ne_ldown*(Phi(ij+t_ldown,llm+1)-Phi(ij,llm+1))
111         F_el(ij+u_ldown,llm+1) = DePhil(ij+u_ldown,llm+1)*W_el
112         W2_el = .5*le_de(ij+u_ldown) * ( W(ij,llm+1)*w_il(ij,llm+1) + W(ij+t_ldown,llm+1)*w_il(ij+t_ldown,llm+1) )
113         v_el(ij+u_ldown,llm+1) = .5*le_de(ij+u_ldown)*(u(ij+u_ldown,llm)+u(ij+u_ldown,llm)) ! checked
114         G_el(ij+u_ldown,llm+1) = v_el(ij+u_ldown,llm+1)*W_el - DePhil(ij+u_ldown,llm+1)*W2_el
115      END DO
116   END IF
117   DO l = ll_begin, ll_endp1
118      !DIR$ SIMD
119      DO ij=ij_begin_ext, ij_end_ext
120         ! compute GradPhi2, dPhi, dW
121         gPhi2=0.
122         dP=0.
123         divG=0
124         gPhi2 = gPhi2 + le_de(ij+u_rup)*DePhil(ij+u_rup,l)**2
125         dP = dP + le_de(ij+u_rup)*DePhil(ij+u_rup,l)*v_el(ij+u_rup,l)
126         divG = divG + ne_rup*G_el(ij+u_rup,l) ! -div(G_el), G_el already has le_de
127         gPhi2 = gPhi2 + le_de(ij+u_lup)*DePhil(ij+u_lup,l)**2
128         dP = dP + le_de(ij+u_lup)*DePhil(ij+u_lup,l)*v_el(ij+u_lup,l)
129         divG = divG + ne_lup*G_el(ij+u_lup,l) ! -div(G_el), G_el already has le_de
130         gPhi2 = gPhi2 + le_de(ij+u_left)*DePhil(ij+u_left,l)**2
131         dP = dP + le_de(ij+u_left)*DePhil(ij+u_left,l)*v_el(ij+u_left,l)
132         divG = divG + ne_left*G_el(ij+u_left,l) ! -div(G_el), G_el already has le_de
133         gPhi2 = gPhi2 + le_de(ij+u_ldown)*DePhil(ij+u_ldown,l)**2
134         dP = dP + le_de(ij+u_ldown)*DePhil(ij+u_ldown,l)*v_el(ij+u_ldown,l)
135         divG = divG + ne_ldown*G_el(ij+u_ldown,l) ! -div(G_el), G_el already has le_de
136         gPhi2 = gPhi2 + le_de(ij+u_rdown)*DePhil(ij+u_rdown,l)**2
137         dP = dP + le_de(ij+u_rdown)*DePhil(ij+u_rdown,l)*v_el(ij+u_rdown,l)
138         divG = divG + ne_rdown*G_el(ij+u_rdown,l) ! -div(G_el), G_el already has le_de
139         gPhi2 = gPhi2 + le_de(ij+u_right)*DePhil(ij+u_right,l)**2
140         dP = dP + le_de(ij+u_right)*DePhil(ij+u_right,l)*v_el(ij+u_right,l)
141         divG = divG + ne_right*G_el(ij+u_right,l) ! -div(G_el), G_el already has le_de
142         gradPhi2(ij,l) = 1./(2.*Ai(ij)) * gPhi2
143         dPhi(ij,l) = gradPhi2(ij,l)*w_il(ij,l) - 1./(2.*Ai(ij))*dP
144         dW(ij,l) = (-1./Ai(ij))*divG
145      END DO
146   END DO
147   ! We need a barrier here because we compute gradPhi2, F_el and w_il above and do a vertical average below
148   !$OMP BARRIER
149   DO l = ll_begin, ll_end
150      !DIR$ SIMD
151      DO ij=ij_begin_ext, ij_end_ext
152         ! Compute berni at scalar points
153         u2=0.
154         u2 = u2 + le_de(ij+u_rup)*u(ij+u_rup,l)**2
155         u2 = u2 + le_de(ij+u_lup)*u(ij+u_lup,l)**2
156         u2 = u2 + le_de(ij+u_left)*u(ij+u_left,l)**2
157         u2 = u2 + le_de(ij+u_ldown)*u(ij+u_ldown,l)**2
158         u2 = u2 + le_de(ij+u_rdown)*u(ij+u_rdown,l)**2
159         u2 = u2 + le_de(ij+u_right)*u(ij+u_right,l)**2
160         berni(ij,l) = 1./(4.*Ai(ij)) * u2 - .25*( gradPhi2(ij,l)*w_il(ij,l)**2 + gradPhi2(ij,l+1)*w_il(ij,l+1)**2 )
161      END DO
162   END DO
163   DO l = ll_begin, ll_end
164      !DIR$ SIMD
165      DO ij=ij_begin_ext, ij_end_ext
166         ! Compute mass flux and grad(berni)
167         uu = .5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l) - .5*( F_el(ij+u_right,l)+F_el(ij+u_right,l+1) )
168         hflux(ij+u_right,l) = le_de(ij+u_right)*uu
169         du(ij+u_right,l) = ne_right*(berni(ij,l)-berni(ij+t_right,l))
170         ! Compute mass flux and grad(berni)
171         uu = .5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l) - .5*( F_el(ij+u_lup,l)+F_el(ij+u_lup,l+1) )
172         hflux(ij+u_lup,l) = le_de(ij+u_lup)*uu
173         du(ij+u_lup,l) = ne_lup*(berni(ij,l)-berni(ij+t_lup,l))
174         ! Compute mass flux and grad(berni)
175         uu = .5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l) - .5*( F_el(ij+u_ldown,l)+F_el(ij+u_ldown,l+1) )
176         hflux(ij+u_ldown,l) = le_de(ij+u_ldown)*uu
177         du(ij+u_ldown,l) = ne_ldown*(berni(ij,l)-berni(ij+t_ldown,l))
178      END DO
179   END DO
180   !---------------------------- caldyn_slow_NH ----------------------------------
181   !--------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.