1 | MODULE compute_caldyn_slow_NH_mod |
---|
2 | USE grid_param |
---|
3 | IMPLICIT NONE |
---|
4 | PRIVATE |
---|
5 | |
---|
6 | #include "../unstructured/unstructured.h90" |
---|
7 | |
---|
8 | PUBLIC :: compute_caldyn_slow_NH |
---|
9 | |
---|
10 | CONTAINS |
---|
11 | |
---|
12 | #ifdef BEGIN_DYSL |
---|
13 | |
---|
14 | KERNEL(caldyn_slow_NH) |
---|
15 | FORALL_CELLS_EXT('1', 'llm+1') |
---|
16 | ON_PRIMAL |
---|
17 | CST_IF(IS_INNER_INTERFACE, w_il(CELL) = 2.*W(CELL)/(rhodz(KDOWN(CELL))+rhodz(KUP(CELL))) ) |
---|
18 | CST_IF(IS_BOTTOM_LEVEL, w_il(CELL) = 2.*W(CELL)/rhodz(KUP(CELL)) ) |
---|
19 | CST_IF(IS_TOP_INTERFACE, w_il(CELL) = 2.*W(CELL)/rhodz(KDOWN(CELL)) ) |
---|
20 | END_BLOCK |
---|
21 | END_BLOCK |
---|
22 | FORALL_CELLS_EXT('1', 'llm+1') |
---|
23 | ON_EDGES |
---|
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(CELL2)+W(CELL1) ) |
---|
28 | DePhil(EDGE) = SIGN*(Phi(CELL2)-Phi(CELL1)) |
---|
29 | F_el(EDGE) = DePhil(EDGE)*W_el |
---|
30 | W2_el = .5*LE_DE * ( W(CELL1)*w_il(CELL1) + W(CELL2)*w_il(CELL2) ) |
---|
31 | v_el(EDGE) = .5*LE_DE*(u(KUP(EDGE))+u(KDOWN(EDGE))) ! checked |
---|
32 | G_el(EDGE) = v_el(EDGE)*W_el - DePhil(EDGE)*W2_el |
---|
33 | END_BLOCK |
---|
34 | END_BLOCK |
---|
35 | |
---|
36 | FORALL_CELLS_EXT('1', 'llm+1') |
---|
37 | ! compute GradPhi2, dPhi, dW |
---|
38 | ON_PRIMAL |
---|
39 | gPhi2=0. |
---|
40 | dP=0. |
---|
41 | divG=0 |
---|
42 | FORALL_EDGES |
---|
43 | gPhi2 = gPhi2 + LE_DE*DePhil(EDGE)**2 |
---|
44 | dP = dP + LE_DE*DePhil(EDGE)*v_el(EDGE) |
---|
45 | divG = divG + SIGN*G_el(EDGE) ! -div(G_el), G_el already has le_de |
---|
46 | END_BLOCK |
---|
47 | gradPhi2(CELL) = 1./(2.*AI) * gPhi2 |
---|
48 | dPhi(CELL) = gradPhi2(CELL)*w_il(CELL) - 1./(2.*AI)*dP |
---|
49 | dW(CELL) = (-1./AI)*divG |
---|
50 | END_BLOCK |
---|
51 | END_BLOCK |
---|
52 | |
---|
53 | ! We need a barrier here because we compute gradPhi2, F_el and w_il above and do a vertical average below |
---|
54 | BARRIER |
---|
55 | |
---|
56 | FORALL_CELLS_EXT() |
---|
57 | ! Compute berni at scalar points |
---|
58 | ON_PRIMAL |
---|
59 | u2=0. |
---|
60 | FORALL_EDGES |
---|
61 | u2 = u2 + LE_DE*u(EDGE)**2 |
---|
62 | END_BLOCK |
---|
63 | berni(CELL) = 1./(4.*AI) * u2 - .25*( gradPhi2(CELL)*w_il(CELL)**2 + gradPhi2(UP(CELL))*w_il(UP(CELL))**2 ) |
---|
64 | END_BLOCK |
---|
65 | END_BLOCK |
---|
66 | |
---|
67 | FORALL_CELLS_EXT() |
---|
68 | ON_EDGES |
---|
69 | ! Compute mass flux and grad(berni) |
---|
70 | uu = .5*(rhodz(CELL1)+rhodz(CELL2))*u(EDGE) - .5*( F_el(EDGE)+F_el(UP(EDGE)) ) |
---|
71 | hflux(EDGE) = LE_DE*uu |
---|
72 | du(EDGE) = SIGN*(berni(CELL1)-berni(CELL2)) |
---|
73 | END_BLOCK |
---|
74 | END_BLOCK |
---|
75 | |
---|
76 | END_BLOCK |
---|
77 | |
---|
78 | #endif END_DYSL |
---|
79 | |
---|
80 | SUBROUTINE compute_caldyn_vert_NH_unst(mass,geopot,W,wflux, eta_dot,wcov,W_etadot, du,dPhi,dW) |
---|
81 | USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT |
---|
82 | USE data_unstructured_mod, ONLY : enter_trace, exit_trace, & |
---|
83 | id_vert_NH, left,right |
---|
84 | FIELD_MASS :: mass, eta_dot, wcov, W_etadot ! IN, BUF*3 |
---|
85 | FIELD_GEOPOT :: geopot,W,wflux,dPhi,dW ! IN*3, INOUT*2 |
---|
86 | FIELD_U :: du ! INOUT |
---|
87 | DECLARE_INDICES |
---|
88 | NUM :: w_ij, wflux_ij |
---|
89 | START_TRACE(id_vert_NH, 6,0,1) |
---|
90 | #include "../kernels_unst/caldyn_vert_NH.k90" |
---|
91 | STOP_TRACE |
---|
92 | END SUBROUTINE compute_caldyn_vert_NH_unst |
---|
93 | |
---|
94 | SUBROUTINE compute_caldyn_slow_NH(u,rhodz,Phi,W, F_el,gradPhi2,w_il, hflux,du,dPhi,dW) |
---|
95 | USE icosa |
---|
96 | USE trace |
---|
97 | USE caldyn_vars_mod |
---|
98 | USE omp_para, ONLY : ll_begin, ll_end,ll_beginp1,ll_endp1 |
---|
99 | REAL(rstd),INTENT(IN) :: u(3*iim*jjm,llm) ! prognostic "velocity" |
---|
100 | REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) ! rho*dz |
---|
101 | REAL(rstd),INTENT(IN) :: Phi(iim*jjm,llm+1) ! prognostic geopotential |
---|
102 | REAL(rstd),INTENT(IN) :: W(iim*jjm,llm+1) ! prognostic vertical momentum |
---|
103 | |
---|
104 | REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s |
---|
105 | REAL(rstd),INTENT(OUT) :: du(3*iim*jjm,llm) |
---|
106 | REAL(rstd),INTENT(OUT) :: dW(iim*jjm,llm+1) |
---|
107 | REAL(rstd),INTENT(OUT) :: dPhi(iim*jjm,llm+1) |
---|
108 | |
---|
109 | REAL(rstd) :: w_il(iim*jjm,llm+1) ! Wil/mil |
---|
110 | REAL(rstd) :: F_el(3*iim*jjm,llm+1) ! NH mass flux |
---|
111 | REAL(rstd) :: gradPhi2(iim*jjm,llm+1) ! grad_Phi**2 |
---|
112 | REAL(rstd) :: DePhil(3*iim*jjm,llm+1) ! grad(Phi) |
---|
113 | |
---|
114 | INTEGER :: ij,l,kdown,kup |
---|
115 | REAL(rstd) :: W_el, W2_el, uu_right, uu_lup, uu_ldown, gPhi2, dP, divG, u2, uu |
---|
116 | |
---|
117 | REAL(rstd) :: berni(iim*jjm,llm) ! Bernoulli function |
---|
118 | REAL(rstd) :: G_el(3*iim*jjm,llm+1) ! horizontal flux of W |
---|
119 | REAL(rstd) :: v_el(3*iim*jjm,llm+1) |
---|
120 | |
---|
121 | REAL(rstd) :: berni1(iim*jjm) ! Bernoulli function |
---|
122 | REAL(rstd) :: G_el1(3*iim*jjm) ! horizontal flux of W |
---|
123 | REAL(rstd) :: v_el1(3*iim*jjm) |
---|
124 | |
---|
125 | CALL trace_start("compute_caldyn_slow_NH") |
---|
126 | |
---|
127 | IF(dysl) THEN |
---|
128 | |
---|
129 | !$OMP BARRIER |
---|
130 | #include "../kernels_hex/caldyn_slow_NH.k90" |
---|
131 | !$OMP BARRIER |
---|
132 | |
---|
133 | ELSE |
---|
134 | |
---|
135 | #define BERNI(ij) berni1(ij) |
---|
136 | #define G_EL(ij) G_el1(ij) |
---|
137 | #define V_EL(ij) v_el1(ij) |
---|
138 | |
---|
139 | DO l=ll_begin, ll_endp1 ! compute on l levels (interfaces) |
---|
140 | IF(l==1) THEN |
---|
141 | kdown=1 |
---|
142 | ELSE |
---|
143 | kdown=l-1 |
---|
144 | END IF |
---|
145 | IF(l==llm+1) THEN |
---|
146 | kup=llm |
---|
147 | ELSE |
---|
148 | kup=l |
---|
149 | END IF |
---|
150 | ! below : "checked" means "formula also valid when kup=kdown (top/bottom)" |
---|
151 | ! compute mil, wil=Wil/mil |
---|
152 | DO ij=ij_begin_ext, ij_end_ext |
---|
153 | w_il(ij,l) = 2.*W(ij,l)/(rhodz(ij,kdown)+rhodz(ij,kup)) ! checked |
---|
154 | END DO |
---|
155 | ! compute DePhi, v_el, G_el, F_el |
---|
156 | ! v_el, W2_el and therefore G_el incorporate metric factor le_de |
---|
157 | ! while DePhil, W_el and F_el don't |
---|
158 | DO ij=ij_begin_ext, ij_end_ext |
---|
159 | ! Compute on edge 'right' |
---|
160 | W_el = .5*( W(ij,l)+W(ij+t_right,l) ) |
---|
161 | DePhil(ij+u_right,l) = ne_right*(Phi(ij+t_right,l)-Phi(ij,l)) |
---|
162 | F_el(ij+u_right,l) = DePhil(ij+u_right,l)*W_el |
---|
163 | W2_el = .5*le_de(ij+u_right) * & |
---|
164 | ( W(ij,l)*w_il(ij,l) + W(ij+t_right,l)*w_il(ij+t_right,l) ) |
---|
165 | V_EL(ij+u_right) = .5*le_de(ij+u_right)*(u(ij+u_right,kup)+u(ij+u_right,kdown)) ! checked |
---|
166 | G_EL(ij+u_right) = V_EL(ij+u_right)*W_el - DePhil(ij+u_right,l)*W2_el |
---|
167 | ! Compute on edge 'lup' |
---|
168 | W_el = .5*( W(ij,l)+W(ij+t_lup,l) ) |
---|
169 | DePhil(ij+u_lup,l) = ne_lup*(Phi(ij+t_lup,l)-Phi(ij,l)) |
---|
170 | F_el(ij+u_lup,l) = DePhil(ij+u_lup,l)*W_el |
---|
171 | W2_el = .5*le_de(ij+u_lup) * & |
---|
172 | ( W(ij,l)*w_il(ij,l) + W(ij+t_lup,l)*w_il(ij+t_lup,l) ) |
---|
173 | V_EL(ij+u_lup) = .5*le_de(ij+u_lup)*( u(ij+u_lup,kup) + u(ij+u_lup,kdown)) ! checked |
---|
174 | G_EL(ij+u_lup) = V_EL(ij+u_lup)*W_el - DePhil(ij+u_lup,l)*W2_el |
---|
175 | ! Compute on edge 'ldown' |
---|
176 | W_el = .5*( W(ij,l)+W(ij+t_ldown,l) ) |
---|
177 | DePhil(ij+u_ldown,l) = ne_ldown*(Phi(ij+t_ldown,l)-Phi(ij,l)) |
---|
178 | F_el(ij+u_ldown,l) = DePhil(ij+u_ldown,l)*W_el |
---|
179 | W2_el = .5*le_de(ij+u_ldown) * & |
---|
180 | ( W(ij,l)*w_il(ij,l) + W(ij+t_ldown,l)*w_il(ij+t_ldown,l) ) |
---|
181 | V_EL(ij+u_ldown) = .5*le_de(ij+u_ldown)*( u(ij+u_ldown,kup) + u(ij+u_ldown,kdown)) ! checked |
---|
182 | G_EL(ij+u_ldown) = V_EL(ij+u_ldown)*W_el - DePhil(ij+u_ldown,l)*W2_el |
---|
183 | END DO |
---|
184 | ! compute GradPhi2, dPhi, dW |
---|
185 | DO ij=ij_begin_ext, ij_end_ext |
---|
186 | gradPhi2(ij,l) = & |
---|
187 | 1/(2*Ai(ij))*(le_de(ij+u_right)*DePhil(ij+u_right,l)**2 + & |
---|
188 | le_de(ij+u_rup)*DePhil(ij+u_rup,l)**2 + & |
---|
189 | le_de(ij+u_lup)*DePhil(ij+u_lup,l)**2 + & |
---|
190 | le_de(ij+u_left)*DePhil(ij+u_left,l)**2 + & |
---|
191 | le_de(ij+u_ldown)*DePhil(ij+u_ldown,l)**2 + & |
---|
192 | le_de(ij+u_rdown)*DePhil(ij+u_rdown,l)**2 ) |
---|
193 | |
---|
194 | dPhi(ij,l) = gradPhi2(ij,l)*w_il(ij,l) -1/(2*Ai(ij))* & |
---|
195 | ( DePhil(ij+u_right,l)*V_EL(ij+u_right) + & ! -v.gradPhi, |
---|
196 | DePhil(ij+u_rup,l)*V_EL(ij+u_rup) + & ! v_el already has le_de |
---|
197 | DePhil(ij+u_lup,l)*V_EL(ij+u_lup) + & |
---|
198 | DePhil(ij+u_left,l)*V_EL(ij+u_left) + & |
---|
199 | DePhil(ij+u_ldown,l)*V_EL(ij+u_ldown) + & |
---|
200 | DePhil(ij+u_rdown,l)*V_EL(ij+u_rdown) ) |
---|
201 | |
---|
202 | dW(ij,l) = -1./Ai(ij)*( & ! -div(G_el), |
---|
203 | ne_right*G_EL(ij+u_right) + & ! G_el already has le_de |
---|
204 | ne_rup*G_EL(ij+u_rup) + & |
---|
205 | ne_lup*G_EL(ij+u_lup) + & |
---|
206 | ne_left*G_EL(ij+u_left) + & |
---|
207 | ne_ldown*G_EL(ij+u_ldown) + & |
---|
208 | ne_rdown*G_EL(ij+u_rdown)) |
---|
209 | END DO |
---|
210 | END DO |
---|
211 | |
---|
212 | DO l=ll_begin, ll_end ! compute on k levels (layers) |
---|
213 | ! Compute berni at scalar points |
---|
214 | DO ij=ij_begin_ext, ij_end_ext |
---|
215 | BERNI(ij) = & |
---|
216 | 1/(4*Ai(ij))*( & |
---|
217 | le_de(ij+u_right)*u(ij+u_right,l)**2 + & |
---|
218 | le_de(ij+u_rup)*u(ij+u_rup,l)**2 + & |
---|
219 | le_de(ij+u_lup)*u(ij+u_lup,l)**2 + & |
---|
220 | le_de(ij+u_left)*u(ij+u_left,l)**2 + & |
---|
221 | le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 + & |
---|
222 | le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) & |
---|
223 | - .25*( gradPhi2(ij,l) *w_il(ij,l)**2 + & |
---|
224 | gradPhi2(ij,l+1)*w_il(ij,l+1)**2 ) |
---|
225 | END DO |
---|
226 | ! Compute mass flux and grad(berni) at edges |
---|
227 | DO ij=ij_begin_ext, ij_end_ext |
---|
228 | ! Compute on edge 'right' |
---|
229 | uu_right = 0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l) & |
---|
230 | -0.5*(F_el(ij+u_right,l)+F_el(ij+u_right,l+1)) |
---|
231 | hflux(ij+u_right,l) = uu_right*le_de(ij+u_right) |
---|
232 | du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right)) |
---|
233 | ! Compute on edge 'lup' |
---|
234 | uu_lup = 0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l) & |
---|
235 | -0.5*(F_el(ij+u_lup,l)+F_el(ij+u_lup,l+1)) |
---|
236 | hflux(ij+u_lup,l) = uu_lup*le_de(ij+u_lup) |
---|
237 | du(ij+u_lup,l) = ne_lup*(BERNI(ij)-BERNI(ij+t_lup)) |
---|
238 | ! Compute on edge 'ldown' |
---|
239 | uu_ldown = 0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l) & |
---|
240 | -0.5*(F_el(ij+u_ldown,l)+F_el(ij+u_ldown,l+1)) |
---|
241 | hflux(ij+u_ldown,l) = uu_ldown*le_de(ij+u_ldown) |
---|
242 | du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown)) |
---|
243 | END DO |
---|
244 | END DO |
---|
245 | |
---|
246 | #undef V_EL |
---|
247 | #undef G_EL |
---|
248 | #undef BERNI |
---|
249 | |
---|
250 | END IF ! dysl |
---|
251 | |
---|
252 | CALL trace_end("compute_caldyn_slow_NH") |
---|
253 | |
---|
254 | END SUBROUTINE compute_caldyn_slow_NH |
---|
255 | |
---|
256 | END MODULE compute_caldyn_slow_NH_mod |
---|