102
105 use scale_const, only: &
106 ohm => const_ohm
107 implicit none
108
109 class(LocalMesh3D), intent(in) :: lmesh
110 class(ElementBase3D), intent(in) :: elem
111 class(LocalMesh2D), intent(in) :: lmesh2D
112 class(ElementBase2D), intent(in) :: elem2D
113 class(ElementOperationBase3D), intent(in) :: element3D_operation
114 type(SparseMat), intent(in) :: Dx, Dy, Dz, Sx, Sy, Sz, Lift
115 real(RP), intent(out) :: DENS_dt(elem%Np,lmesh%NeA)
116 real(RP), intent(out) :: MOMX_dt(elem%Np,lmesh%NeA)
117 real(RP), intent(out) :: MOMY_dt(elem%Np,lmesh%NeA)
118 real(RP), intent(out) :: MOMZ_dt(elem%Np,lmesh%NeA)
119 real(RP), intent(out) :: EnTot_dt(elem%Np,lmesh%NeA)
120 real(RP), intent(in) :: DDENS_(elem%Np,lmesh%NeA)
121 real(RP), intent(in) :: MOMX_(elem%Np,lmesh%NeA)
122 real(RP), intent(in) :: MOMY_(elem%Np,lmesh%NeA)
123 real(RP), intent(in) :: MOMZ_(elem%Np,lmesh%NeA)
124 real(RP), intent(in) :: ETOT_(elem%Np,lmesh%NeA)
125 real(RP), intent(in) :: DPRES_(elem%Np,lmesh%NeA)
126 real(RP), intent(in) :: DENS_hyd(elem%Np,lmesh%NeA)
127 real(RP), intent(in) :: PRES_hyd(elem%Np,lmesh%NeA)
128 real(RP), intent(in) :: PRES_hyd_ref(elem%Np,lmesh%NeA)
129 real(RP), intent(in) :: THERM_hyd(elem%Np,lmesh%NeA)
130 real(RP), intent(in) :: CORIOLIS(elem2D%Np,lmesh2D%NeA)
131 real(RP), intent(in) :: Rtot (elem%Np,lmesh%NeA)
132 real(RP), intent(in) :: CVtot(elem%Np,lmesh%NeA)
133 real(RP), intent(in) :: CPtot(elem%Np,lmesh%NeA)
134 real(RP), intent(in) :: DPhydDx(elem%Np,lmesh%NeA)
135 real(RP), intent(in) :: DPhydDy(elem%Np,lmesh%NeA)
136
137 real(RP) :: Fx(elem%Np), Fy(elem%Np), Fz(elem%Np), LiftDelFlx(elem%Np)
138 real(RP) :: DPRES_hyd(elem%Np), GradPhyd_x(elem%Np), GradPhyd_y(elem%Np)
139 real(RP) :: del_flux(elem%NfpTot,lmesh%Ne,PRGVAR_NUM)
140 real(RP) :: del_flux_hyd(elem%NfpTot,lmesh%Ne,2)
141 real(RP) :: rdens_(elem%Np), u_(elem%Np), v_(elem%Np), w_(elem%Np), wt_(elem%np), drho(elem%Np)
142
143 real(RP) :: G11(elem%Np), G12(elem%Np), G22(elem%Np)
144 real(RP) :: GsqrtV(elem%Np), RGsqrtV(elem%Np)
145 real(RP) :: X2D(elem2D%Np,lmesh2D%Ne), Y2D(elem2D%Np,lmesh2D%Ne)
146 real(RP) :: X(elem%Np), Y(elem%Np), twoOVdel2(elem%Np)
147 real(RP) :: CORI(elem%Np,2)
148 logical :: is_panel1to4
149 real(RP) :: s
150
151 integer :: ke, ke2d
152 integer :: p, p12, p3
153
154 real(RP) :: gamm, rgamm
155 real(RP) :: rP0
156 real(RP) :: RovP0, P0ovR
157
158 real(RP) :: enthalpy_(elem%Np)
159 real(RP) :: u1_(elem%Np), u2_(elem%Np)
160
161
162 call prof_rapstart('cal_dyn_tend_bndflux', 3)
163 call get_ebnd_flux( &
164 del_flux, del_flux_hyd, &
165 ddens_, momx_, momy_, momz_, etot_, dpres_, dens_hyd, pres_hyd, &
166 rtot, cvtot, cptot, &
167 lmesh%Gsqrt, lmesh%GIJ(:,:,1,1), lmesh%GIJ(:,:,1,2), lmesh%GIJ(:,:,2,2), &
168 lmesh%G_ij(:,:,1,1), lmesh%G_ij(:,:,1,2), lmesh%G_ij(:,:,2,2), &
169 lmesh%GsqrtH, lmesh%GI3(:,:,1), lmesh%GI3(:,:,2), lmesh%zlev(:,:), &
170 lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), &
171 lmesh%vmapM, lmesh%vmapP, elem%IndexH2Dto3D_bnd, &
172 lmesh, elem, lmesh2d, elem2d )
173 call prof_rapend('cal_dyn_tend_bndflux', 3)
174
175
176 call prof_rapstart('cal_dyn_tend_interior', 3)
177 gamm = cpdry / cvdry
178 rgamm = cvdry / cpdry
179 rp0 = 1.0_rp / pres00
180 rovp0 = rdry * rp0
181 p0ovr = pres00 / rdry
182
183 s = 1.0_rp
184 is_panel1to4 = .true.
185 if ( lmesh%panelID == 5 ) then
186 is_panel1to4 = .false.
187 else if ( lmesh%panelID == 6 ) then
188 is_panel1to4 = .false.
189 s = - 1.0_rp
190 end if
191
192
193
194
195
196
197
198
199
200
201
202 do ke2d = lmesh2d%NeS, lmesh2d%NeE
203 x2d(:,ke2d) = tan(lmesh2d%pos_en(:,ke2d,1))
204 y2d(:,ke2d) = tan(lmesh2d%pos_en(:,ke2d,2))
205 end do
206
207
208 do ke = lmesh%NeS, lmesh%NeE
209
210 ke2d = lmesh%EMap3Dto2D(ke)
211 g11(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,1,1)
212 g12(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,1,2)
213 g22(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,2,2)
214 gsqrtv(:) = lmesh%Gsqrt(:,ke) / lmesh%GsqrtH(elem%IndexH2Dto3D,ke2d)
215 rgsqrtv(:) = 1.0_rp / gsqrtv(:)
216
217
218 rdens_(:) = 1.0_rp / ( ddens_(:,ke) + dens_hyd(:,ke) )
219 u_(:) = momx_(:,ke) * rdens_(:)
220 v_(:) = momy_(:,ke) * rdens_(:)
221 w_(:) = momz_(:,ke) * rdens_(:)
222 wt_(:) = w_(:) * rgsqrtv(:) + lmesh%GI3(:,ke,1) * u_(:) + lmesh%GI3(:,ke,2) * v_(:)
223 u1_(:) = lmesh%G_ij(elem%IndexH2Dto3D(:),ke2d,1,1) * u_(:) + lmesh%G_ij(elem%IndexH2Dto3D(:),ke2d,2,1) * v_(:)
224 u2_(:) = lmesh%G_ij(elem%IndexH2Dto3D(:),ke2d,2,1) * u_(:) + lmesh%G_ij(elem%IndexH2Dto3D(:),ke2d,2,2) * v_(:)
225
226
227
228
229
230
231 enthalpy_(:) = etot_(:,ke) + pres_hyd(:,ke) + dpres_(:,ke)
232
233 x(:) = x2d(elem%IndexH2Dto3D,ke2d)
234 y(:) = y2d(elem%IndexH2Dto3D,ke2d)
235 twoovdel2(:) = 2.0_rp / ( 1.0_rp + x(:)**2 + y(:)**2 )
236
237 cori(:,1) = s * ohm * twoovdel2(:) * ( - x(:) * y(:) * momx_(:,ke) + (1.0_rp + y(:)**2) * momy_(:,ke) )
238 cori(:,2) = s * ohm * twoovdel2(:) * ( - (1.0_rp + x(:)**2) * momx_(:,ke) + x(:) * y(:) * momy_(:,ke) )
239
240 if ( is_panel1to4 ) then
241 cori(:,1) = s * y(:) * cori(:,1)
242 cori(:,2) = s * y(:) * cori(:,2)
243 end if
244
245 drho(:) = matmul(intrpmat_vpordm1, ddens_(:,ke))
246
247
248
249 dpres_hyd(:) = pres_hyd(:,ke) - pres_hyd_ref(:,ke)
250
251 call sparsemat_matmul(dx, gsqrtv(:) * dpres_hyd(:), fx)
252 call sparsemat_matmul(dz, gsqrtv(:) * lmesh%GI3(:,ke,1) * dpres_hyd(:), fz)
253 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux_hyd(:,ke,1), liftdelflx)
254 gradphyd_x(:) = lmesh%Escale(:,ke,1,1) * fx(:) &
255 + lmesh%Escale(:,ke,3,3) * fz(:) &
256 + liftdelflx(:)
257
258 call sparsemat_matmul(dy, gsqrtv(:) * dpres_hyd(:), fy)
259 call sparsemat_matmul(dz, gsqrtv(:) * lmesh%GI3(:,ke,2) * dpres_hyd(:), fz)
260 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux_hyd(:,ke,2), liftdelflx)
261 gradphyd_y(:) = lmesh%Escale(:,ke,2,2) * fy(:) &
262 + lmesh%Escale(:,ke,3,3) * fz(:) &
263 + liftdelflx(:)
264
265
266 call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * momx_(:,ke), fx)
267 call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * momy_(:,ke), fy)
268 call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( ddens_(:,ke) + dens_hyd(:,ke) ) * wt_(:), fz)
269 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,dens_vid), liftdelflx)
270
271 dens_dt(:,ke) = - ( &
272 lmesh%Escale(:,ke,1,1) * fx(:) &
273 + lmesh%Escale(:,ke,2,2) * fy(:) &
274 + lmesh%Escale(:,ke,3,3) * fz(:) &
275 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke)
276
277
278 call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * ( u_(:) * momx_(:,ke) + g11(:) * dpres_(:,ke) ), fx)
279 call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * ( v_(:) * momx_(:,ke) + g12(:) * dpres_(:,ke) ), fy)
280 call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( wt_(:) * momx_(:,ke) &
281 + ( lmesh%GI3(:,ke,1) * g11(:) + lmesh%GI3(:,ke,2) * g12(:) ) * dpres_(:,ke) ), fz)
282 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momx_vid), liftdelflx)
283
284 momx_dt(:,ke) = &
285 - ( lmesh%Escale(:,ke,1,1) * fx(:) &
286 + lmesh%Escale(:,ke,2,2) * fy(:) &
287 + lmesh%Escale(:,ke,3,3) * fz(:) &
288 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke) &
289 + twoovdel2(:) * y(:) * &
290 ( - x(:) * y(:) * u_(:) + (1.0_rp + y(:)**2) * v_(:) ) * momx_(:,ke) &
291 - ( g11(:) * gradphyd_x(:) + g12(:) * gradphyd_y(:) ) * rgsqrtv(:) &
292 + cori(:,1)
293
294
295 call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * ( u_(:) * momy_(:,ke) + g12(:) * dpres_(:,ke) ), fx)
296 call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * ( v_(:) * momy_(:,ke) + g22(:) * dpres_(:,ke) ), fy)
297 call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( wt_(:) * momy_(:,ke) &
298 + ( lmesh%GI3(:,ke,1) * g12(:) + lmesh%GI3(:,ke,2) * g22(:) ) * dpres_(:,ke) ), fz)
299 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momy_vid), liftdelflx)
300
301 momy_dt(:,ke) = &
302 - ( lmesh%Escale(:,ke,1,1) * fx(:) &
303 + lmesh%Escale(:,ke,2,2) * fy(:) &
304 + lmesh%Escale(:,ke,3,3) * fz(:) &
305 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke) &
306 + twoovdel2(:) * x(:) * &
307 ( (1.0_rp + x(:)**2) * u_(:) - x(:) * y(:) * v_(:) ) * momy_(:,ke) &
308 - ( g12(:) * gradphyd_x(:) + g22(:) * gradphyd_y(:) ) * rgsqrtv(:) &
309 + cori(:,2)
310
311
312 call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * u_(:) * momz_(:,ke), fx)
313 call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * v_(:) * momz_(:,ke), fy)
314 call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( wt_(:) * momz_(:,ke) + rgsqrtv(:) * dpres_(:,ke) ), fz)
315 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momz_vid), liftdelflx)
316
317 momz_dt(:,ke) = &
318 - ( lmesh%Escale(:,ke,1,1) * fx(:) &
319 + lmesh%Escale(:,ke,2,2) * fy(:) &
320 + lmesh%Escale(:,ke,3,3) * fz(:) &
321 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke) &
322 - grav * drho(:)
323
324
325 call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * u_(:) * enthalpy_(:), fx)
326 call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * v_(:) * enthalpy_(:), fy)
327 call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * wt_(:) * enthalpy_(:), fz)
328 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,etot_vid), liftdelflx)
329
330 entot_dt(:,ke) = &
331 - ( lmesh%Escale(:,ke,1,1) * fx(:) &
332 + lmesh%Escale(:,ke,2,2) * fy(:) &
333 + lmesh%Escale(:,ke,3,3) * fz(:) &
334 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke)
335 end do
336
337
338 call prof_rapend('cal_dyn_tend_interior', 3)
339
340 return
module FElib / Fluid dyn solver / Atmosphere / Nonhydrostatic model / Numflux
subroutine, public atm_dyn_dgm_nonhydro3d_etot_heve_numflux_get_generalhvc(del_flux, del_flux_hyd, ddens_, momx_, momy_, momz_, etot_, dpres_, dens_hyd, pres_hyd, rtot, cvtot, cptot, gsqrt, g11, g12, g22, g_11, g_12, g_22, gsqrth, g13, g23, zlev, nx, ny, nz, vmapm, vmapp, im2dto3d, lmesh, elem, lmesh2d, elem2d)