99 DENS_dt, MOMX_dt, MOMY_dt, MOMZ_dt, RHOT_dt, & ! (out)
100 ddens_, momx_, momy_, momz_, drhot_, dpres_, &
101 dens_hyd, pres_hyd, pres_hyd_ref, therm_hyd, &
102 coriolis, rtot, cvtot, cptot, dphyddx, dphyddy, &
103 element3d_operation, dx, dy, dz, sx, sy, sz, lift, &
104 lmesh, elem, lmesh2d, elem2d )
108 use scale_const,
only: &
117 type(
sparsemat),
intent(in) :: dx, dy, dz, sx, sy, sz, lift
118 real(rp),
intent(out) :: dens_dt(elem%np,lmesh%nea)
119 real(rp),
intent(out) :: momx_dt(elem%np,lmesh%nea)
120 real(rp),
intent(out) :: momy_dt(elem%np,lmesh%nea)
121 real(rp),
intent(out) :: momz_dt(elem%np,lmesh%nea)
122 real(rp),
intent(out) :: rhot_dt(elem%np,lmesh%nea)
123 real(rp),
intent(in) :: ddens_(elem%np,lmesh%nea)
124 real(rp),
intent(in) :: momx_(elem%np,lmesh%nea)
125 real(rp),
intent(in) :: momy_(elem%np,lmesh%nea)
126 real(rp),
intent(in) :: momz_(elem%np,lmesh%nea)
127 real(rp),
intent(in) :: drhot_(elem%np,lmesh%nea)
128 real(rp),
intent(in) :: dpres_(elem%np,lmesh%nea)
129 real(rp),
intent(in) :: dens_hyd(elem%np,lmesh%nea)
130 real(rp),
intent(in) :: pres_hyd(elem%np,lmesh%nea)
131 real(rp),
intent(in) :: pres_hyd_ref(elem%np,lmesh%nea)
132 real(rp),
intent(in) :: therm_hyd(elem%np,lmesh%nea)
133 real(rp),
intent(in) :: coriolis(elem2d%np,lmesh2d%nea)
134 real(rp),
intent(in) :: rtot (elem%np,lmesh%nea)
135 real(rp),
intent(in) :: cvtot(elem%np,lmesh%nea)
136 real(rp),
intent(in) :: cptot(elem%np,lmesh%nea)
137 real(rp),
intent(in) :: dphyddx(elem%np,lmesh%nea)
138 real(rp),
intent(in) :: dphyddy(elem%np,lmesh%nea)
140 real(rp) :: fx(elem%np), fy(elem%np), fz(elem%np), liftdelflx(elem%np)
141 real(rp) :: dpres_hyd(elem%np), gradphyd_x(elem%np), gradphyd_y(elem%np)
142 real(rp) :: del_flux(elem%nfptot,lmesh%ne,
prgvar_num)
143 real(rp) :: del_flux_hyd(elem%nfptot,lmesh%ne,2)
144 real(rp) :: rhot_(elem%np)
145 real(rp) :: rdens_(elem%np), u_(elem%np), v_(elem%np), w_(elem%np), wt_(elem%np), drho(elem%np)
147 real(rp) :: g11(elem%np), g12(elem%np), g22(elem%np)
148 real(rp) :: gsqrtv(elem%np), rgsqrtv(elem%np)
149 real(rp) :: x2d(elem2d%np,lmesh2d%ne), y2d(elem2d%np,lmesh2d%ne)
150 real(rp) :: x(elem%np), y(elem%np), twoovdel2(elem%np)
151 real(rp) :: cori(elem%np,2)
152 logical :: is_panel1to4
156 integer :: p, p12, p3
158 real(rp) :: gamm, rgamm
160 real(rp) :: rovp0, p0ovr
163 call prof_rapstart(
'cal_dyn_tend_bndflux', 3)
164 call get_ebnd_flux( &
165 del_flux, del_flux_hyd, &
166 ddens_, momx_, momy_, momz_, drhot_, dpres_, dens_hyd, pres_hyd, &
167 rtot, cvtot, cptot, &
168 lmesh%Gsqrt, lmesh%GIJ(:,:,1,1), lmesh%GIJ(:,:,1,2), lmesh%GIJ(:,:,2,2), &
169 lmesh%GsqrtH, lmesh%gam, lmesh%GI3(:,:,1), lmesh%GI3(:,:,2), &
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)
176 call prof_rapstart(
'cal_dyn_tend_interior', 3)
178 rgamm = cvdry / cpdry
179 rp0 = 1.0_rp / pres00
181 p0ovr = pres00 / rdry
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.
201 do ke2d = lmesh2d%NeS, lmesh2d%NeE
202 x2d(:,ke2d) = tan(lmesh2d%pos_en(:,ke2d,1))
203 y2d(:,ke2d) = tan(lmesh2d%pos_en(:,ke2d,2))
207 do ke = lmesh%NeS, lmesh%NeE
209 ke2d = lmesh%EMap3Dto2D(ke)
210 g11(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,1,1)
211 g12(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,1,2)
212 g22(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,2,2)
213 gsqrtv(:) = lmesh%Gsqrt(:,ke) / lmesh%GsqrtH(elem%IndexH2Dto3D,ke2d)
214 rgsqrtv(:) = 1.0_rp / gsqrtv(:)
217 rhot_(:) = p0ovr * ( pres_hyd(:,ke) * rp0 )**rgamm + drhot_(:,ke)
221 rdens_(:) = 1.0_rp / ( ddens_(:,ke) + dens_hyd(:,ke) )
222 u_(:) = momx_(:,ke) * rdens_(:)
223 v_(:) = momy_(:,ke) * rdens_(:)
224 w_(:) = momz_(:,ke) * rdens_(:)
225 wt_(:) = w_(:) * rgsqrtv(:) + lmesh%GI3(:,ke,1) * u_(:) + lmesh%GI3(:,ke,2) * v_(:)
227 x(:) = x2d(elem%IndexH2Dto3D,ke2d)
228 y(:) = y2d(elem%IndexH2Dto3D,ke2d)
229 twoovdel2(:) = 2.0_rp / ( 1.0_rp + x(:)**2 + y(:)**2 )
231 cori(:,1) = s * ohm * twoovdel2(:) * ( - x(:) * y(:) * momx_(:,ke) + ( 1.0_rp + y(:)**2 ) * momy_(:,ke) )
232 cori(:,2) = s * ohm * twoovdel2(:) * ( - ( 1.0_rp + x(:)**2 ) * momx_(:,ke) + x(:) * y(:) * momy_(:,ke) )
233 if ( is_panel1to4 )
then
234 cori(:,1) = s * y(:) * cori(:,1)
235 cori(:,2) = s * y(:) * cori(:,2)
242 dpres_hyd(:) = pres_hyd(:,ke) - pres_hyd_ref(:,ke)
246 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux_hyd(:,ke,1), liftdelflx)
247 gradphyd_x(:) = lmesh%Escale(:,ke,1,1) * fx(:) &
248 + lmesh%Escale(:,ke,3,3) * fz(:) &
253 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux_hyd(:,ke,2), liftdelflx)
254 gradphyd_y(:) = lmesh%Escale(:,ke,2,2) * fy(:) &
255 + lmesh%Escale(:,ke,3,3) * fz(:) &
261 call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( ddens_(:,ke) + dens_hyd(:,ke) ) * wt_(:), fz)
262 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,dens_vid), liftdelflx)
264 dens_dt(:,ke) = - ( &
265 lmesh%Escale(:,ke,1,1) * fx(:) &
266 + lmesh%Escale(:,ke,2,2) * fy(:) &
267 + lmesh%Escale(:,ke,3,3) * fz(:) &
268 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke)
271 call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * ( u_(:) * momx_(:,ke) + g11(:) * dpres_(:,ke) ), fx)
272 call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * ( v_(:) * momx_(:,ke) + g12(:) * dpres_(:,ke) ), fy)
274 + ( lmesh%GI3(:,ke,1) * g11(:) + lmesh%GI3(:,ke,2) * g12(:) ) * dpres_(:,ke) ), fz)
275 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momx_vid), liftdelflx)
278 - ( lmesh%Escale(:,ke,1,1) * fx(:) &
279 + lmesh%Escale(:,ke,2,2) * fy(:) &
280 + lmesh%Escale(:,ke,3,3) * fz(:) &
281 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke) &
282 - twoovdel2(:) * y(:) * &
283 ( x(:) * y(:) * u_(:) - (1.0_rp + y(:)**2) * v_(:) ) * momx_(:,ke) &
284 - ( g11(:) * gradphyd_x(:) + g12(:) * gradphyd_y(:) ) * rgsqrtv(:) &
288 call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * ( u_(:) * momy_(:,ke) + g12(:) * dpres_(:,ke) ), fx)
289 call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * ( v_(:) * momy_(:,ke) + g22(:) * dpres_(:,ke) ), fy)
291 + ( lmesh%GI3(:,ke,1) * g12(:) + lmesh%GI3(:,ke,2) * g22(:) ) * dpres_(:,ke) ), fz)
292 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momy_vid), liftdelflx)
295 - ( lmesh%Escale(:,ke,1,1) * fx(:) &
296 + lmesh%Escale(:,ke,2,2) * fy(:) &
297 + lmesh%Escale(:,ke,3,3) * fz(:) &
298 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke) &
299 - twoovdel2(:) * x(:) * &
300 ( - (1.0_rp + x(:)**2) * u_(:) + x(:) * y(:) * v_(:) ) * momy_(:,ke) &
301 - ( g12(:) * gradphyd_x(:) + g22(:) * gradphyd_y(:) ) * rgsqrtv(:) &
307 call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( wt_(:) * momz_(:,ke) + rgsqrtv(:) * dpres_(:,ke) ), fz)
308 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momz_vid), liftdelflx)
311 - ( lmesh%Escale(:,ke,1,1) * fx(:) &
312 + lmesh%Escale(:,ke,2,2) * fy(:) &
313 + lmesh%Escale(:,ke,3,3) * fz(:) &
314 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke) &
321 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,rhot_vid), liftdelflx)
324 - ( lmesh%Escale(:,ke,1,1) * fx(:) &
325 + lmesh%Escale(:,ke,2,2) * fy(:) &
326 + lmesh%Escale(:,ke,3,3) * fz(:) &
327 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke)
332 call prof_rapend(
'cal_dyn_tend_interior', 3)
339 DENS_dt, MOMX_dt, MOMY_dt, MOMZ_dt, RHOT_dt, & ! (out)
340 ddens_, momx_, momy_, momz_, drhot_, dpres_, &
341 dens_hyd, pres_hyd, pres_hyd_ref, therm_hyd, &
342 coriolis, rtot, cvtot, cptot, dphyddx, dphyddy, &
343 element3d_operation, dx, dy, dz, sx, sy, sz, lift, &
344 lmesh, elem, lmesh2d, elem2d )
348 use scale_const,
only: &
357 type(
sparsemat),
intent(in) :: dx, dy, dz, sx, sy, sz, lift
358 real(rp),
intent(out) :: dens_dt(elem%np,lmesh%nea)
359 real(rp),
intent(out) :: momx_dt(elem%np,lmesh%nea)
360 real(rp),
intent(out) :: momy_dt(elem%np,lmesh%nea)
361 real(rp),
intent(out) :: momz_dt(elem%np,lmesh%nea)
362 real(rp),
intent(out) :: rhot_dt(elem%np,lmesh%nea)
363 real(rp),
intent(in) :: ddens_(elem%np,lmesh%nea)
364 real(rp),
intent(in) :: momx_(elem%np,lmesh%nea)
365 real(rp),
intent(in) :: momy_(elem%np,lmesh%nea)
366 real(rp),
intent(in) :: momz_(elem%np,lmesh%nea)
367 real(rp),
intent(in) :: drhot_(elem%np,lmesh%nea)
368 real(rp),
intent(in) :: dpres_(elem%np,lmesh%nea)
369 real(rp),
intent(in) :: dens_hyd(elem%np,lmesh%nea)
370 real(rp),
intent(in) :: pres_hyd(elem%np,lmesh%nea)
371 real(rp),
intent(in) :: pres_hyd_ref(elem%np,lmesh%nea)
372 real(rp),
intent(in) :: therm_hyd(elem%np,lmesh%nea)
373 real(rp),
intent(in) :: coriolis(elem2d%np,lmesh2d%nea)
374 real(rp),
intent(in) :: rtot (elem%np,lmesh%nea)
375 real(rp),
intent(in) :: cvtot(elem%np,lmesh%nea)
376 real(rp),
intent(in) :: cptot(elem%np,lmesh%nea)
377 real(rp),
intent(in) :: dphyddx(elem%np,lmesh%nea)
378 real(rp),
intent(in) :: dphyddy(elem%np,lmesh%nea)
380 real(rp) :: flux(elem%np,3,5), dflux(elem%np,4,5)
381 real(rp) :: del_flux(elem%nfptot,
prgvar_num,lmesh%ne)
382 real(rp) :: u_, v_, w_, pt_
383 real(rp) :: drho(elem%np)
384 real(rp) :: rdens_(elem%np), gsqrtv(elem%np), rgsqrtv(elem%np), rgsqrt(elem%np)
385 real(rp) :: gsqrt_, gsqrtdpres_, e11, e22, e33
387 real(rp) :: g11(elem%np), g12(elem%np), g22(elem%np)
388 real(rp) :: x2d(elem2d%np,lmesh2d%ne), y2d(elem2d%np,lmesh2d%ne)
389 real(rp) :: x(elem%np), y(elem%np), twoovdel2(elem%np)
390 real(rp) :: cori(elem%np,2)
391 logical :: is_panel1to4
397 real(rp) :: gamm, rgamm
399 real(rp) :: rovp0, p0ovr
402 call prof_rapstart(
'cal_dyn_tend_bndflux', 3)
403 call get_ebnd_flux( &
405 ddens_, momx_, momy_, momz_, drhot_, dpres_, &
406 dens_hyd, pres_hyd, therm_hyd, rtot, cvtot, cptot, &
407 lmesh%Gsqrt, lmesh%GIJ(:,:,1,1), lmesh%GIJ(:,:,1,2), lmesh%GIJ(:,:,2,2), &
408 lmesh%GsqrtH, lmesh%gam, lmesh%GI3(:,:,1), lmesh%GI3(:,:,2), &
409 lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), &
410 lmesh%vmapM, lmesh%vmapP, elem%IndexH2Dto3D_bnd, &
411 lmesh, elem, lmesh2d, elem2d )
412 call prof_rapend(
'cal_dyn_tend_bndflux', 3)
415 call prof_rapstart(
'cal_dyn_tend_interior', 3)
417 rgamm = cvdry / cpdry
418 rp0 = 1.0_rp / pres00
420 p0ovr = pres00 / rdry
423 is_panel1to4 = .true.
424 if ( lmesh%panelID == 5 )
then
425 is_panel1to4 = .false.
426 else if ( lmesh%panelID == 6 )
then
427 is_panel1to4 = .false.
439 do ke2d = lmesh2d%NeS, lmesh2d%NeE
440 x2d(:,ke2d) = tan(lmesh2d%pos_en(:,ke2d,1))
441 y2d(:,ke2d) = tan(lmesh2d%pos_en(:,ke2d,2))
445 do ke = lmesh%NeS, lmesh%NeE
447 ke2d = lmesh%EMap3Dto2D(ke)
450 g11(p) = lmesh%GIJ(elem%IndexH2Dto3D(p),ke2d,1,1)
451 g12(p) = lmesh%GIJ(elem%IndexH2Dto3D(p),ke2d,1,2)
452 g22(p) = lmesh%GIJ(elem%IndexH2Dto3D(p),ke2d,2,2)
455 gsqrtv(p) = lmesh%Gsqrt(p,ke) / lmesh%GsqrtH(elem%IndexH2Dto3D(p),ke2d)
456 rgsqrtv(p) = 1.0_rp / gsqrtv(p)
457 rgsqrt(p) = 1.0_rp / lmesh%Gsqrt(p,ke)
458 rdens_(p) = 1.0_rp / ( ddens_(p,ke) + dens_hyd(p,ke) )
464 gsqrt_ = lmesh%Gsqrt(p,ke)
465 flux(p,1,dens_vid) = gsqrt_ * momx_(p,ke)
466 flux(p,2,dens_vid) = gsqrt_ * momy_(p,ke)
467 flux(p,3,dens_vid) = gsqrt_ * ( &
468 momz_(p,ke) * rgsqrtv(p) &
469 + lmesh%GI3(p,ke,1) * momx_(p,ke) &
470 + lmesh%GI3(p,ke,2) * momy_(p,ke) )
473 pt_ = ( therm_hyd(p,ke) + drhot_(p,ke) ) * rdens_(p)
475 flux(p,1,rhot_vid) = flux(p,1,dens_vid) * pt_
476 flux(p,2,rhot_vid) = flux(p,2,dens_vid) * pt_
477 flux(p,3,rhot_vid) = flux(p,3,dens_vid) * pt_
479 w_ = momz_(p,ke) * rdens_(p)
480 flux(p,1,momz_vid) = flux(p,1,dens_vid) * w_
481 flux(p,2,momz_vid) = flux(p,2,dens_vid) * w_
482 flux(p,3,momz_vid) = flux(p,3,dens_vid) * w_ + lmesh%Gsqrt(p,ke) * rgsqrtv(p) * dpres_(p,ke)
485 gsqrtdpres_ = lmesh%Gsqrt(p,ke) * dpres_(p,ke)
487 u_ = momx_(p,ke) * rdens_(p)
488 v_ = momy_(p,ke) * rdens_(p)
490 flux(p,1,momx_vid) = flux(p,1,dens_vid) * u_ + g11(p) * gsqrtdpres_
491 flux(p,2,momx_vid) = flux(p,2,dens_vid) * u_ + g12(p) * gsqrtdpres_
492 flux(p,3,momx_vid) = flux(p,3,dens_vid) * u_ + gsqrtdpres_ * ( g11(p) * lmesh%GI3(p,ke,1) + g12(p) * lmesh%GI3(p,ke,2) )
494 flux(p,1,momy_vid) = flux(p,1,dens_vid) * v_ + g12(p) * gsqrtdpres_
495 flux(p,2,momy_vid) = flux(p,2,dens_vid) * v_ + g22(p) * gsqrtdpres_
496 flux(p,3,momy_vid) = flux(p,3,dens_vid) * v_ + gsqrtdpres_ * ( g12(p) * lmesh%GI3(p,ke,1) + g22(p) * lmesh%GI3(p,ke,2) )
499 call element3d_operation%Div_var5( &
500 flux, del_flux(:,:,ke), &
505 e11 = lmesh%Escale(p,ke,1,1)
506 e22 = lmesh%Escale(p,ke,2,2)
507 e33 = lmesh%Escale(p,ke,3,3)
509 dens_dt(p,ke) = - ( &
510 e11 * dflux(p,1,dens_vid) &
511 + e22 * dflux(p,2,dens_vid) &
512 + e33 * dflux(p,3,dens_vid) &
513 + dflux(p,4,dens_vid) ) * rgsqrt(p)
515 rhot_dt(p,ke) = - ( &
516 e11 * dflux(p,1,rhot_vid) &
517 + e22 * dflux(p,2,rhot_vid) &
518 + e33 * dflux(p,3,rhot_vid) &
519 + dflux(p,4,rhot_vid) ) * rgsqrt(p)
522 call element3d_operation%VFilterPM1( ddens_(:,ke), &
526 e11 = lmesh%Escale(p,ke,1,1)
527 e22 = lmesh%Escale(p,ke,2,2)
528 e33 = lmesh%Escale(p,ke,3,3)
530 momz_dt(p,ke) = - ( &
531 e11 * dflux(p,1,momz_vid) &
532 + e22 * dflux(p,2,momz_vid) &
533 + e33 * dflux(p,3,momz_vid) &
534 + dflux(p,4,momz_vid) ) * rgsqrt(p) &
540 x(p) = x2d(elem%IndexH2Dto3D(p),ke2d)
541 y(p) = y2d(elem%IndexH2Dto3D(p),ke2d)
542 twoovdel2(p) = 2.0_rp / ( 1.0_rp + x(p)**2 + y(p)**2 )
546 cori(p,1) = s * ohm * twoovdel2(p) * ( - x(p) * y(p) * momx_(p,ke) + ( 1.0_rp + y(p)**2 ) * momy_(p,ke) )
547 cori(p,2) = s * ohm * twoovdel2(p) * ( - ( 1.0_rp + x(p)**2 ) * momx_(p,ke) + x(p) * y(p) * momy_(p,ke) )
549 if ( is_panel1to4 )
then
551 cori(p,1) = s * y(p) * cori(p,1)
552 cori(p,2) = s * y(p) * cori(p,2)
557 u_ = momx_(p,ke) * rdens_(p)
558 v_ = momy_(p,ke) * rdens_(p)
560 momx_dt(p,ke) = - ( g11(p) * dphyddx(p,ke) + g12(p) * dphyddy(p,ke) ) &
561 - twoovdel2(p) * y(p) * &
562 ( x(p) * y(p) * u_ - (1.0_rp + y(p)**2) * v_ ) * momx_(p,ke) &
565 momy_dt(p,ke) = - ( g12(p) * dphyddx(p,ke) + g22(p) * dphyddy(p,ke) ) &
566 - twoovdel2(p) * x(p) * &
567 ( - (1.0_rp + x(p)**2) * u_ + x(p) * y(p) * v_ ) * momy_(p,ke) &
572 e11 = lmesh%Escale(p,ke,1,1)
573 e22 = lmesh%Escale(p,ke,2,2)
574 e33 = lmesh%Escale(p,ke,3,3)
576 momx_dt(p,ke) = momx_dt(p,ke) - ( &
577 e11 * dflux(p,1,momx_vid) &
578 + e22 * dflux(p,2,momx_vid) &
579 + e33 * dflux(p,3,momx_vid) &
580 + dflux(p,4,momx_vid) ) * rgsqrt(p)
582 momy_dt(p,ke) = momy_dt(p,ke) - ( &
583 e11 * dflux(p,1,momy_vid) &
584 + e22 * dflux(p,2,momy_vid) &
585 + e33 * dflux(p,3,momy_vid) &
586 + dflux(p,4,momy_vid) ) * rgsqrt(p)
591 call prof_rapend(
'cal_dyn_tend_interior', 3)
598 DENS_dt, MOMX_dt, MOMY_dt, MOMZ_dt, RHOT_dt, & ! (out)
599 ddens_, momx_, momy_, momz_, drhot_, dpres_, &
600 dens_hyd, pres_hyd, pres_hyd_ref, therm_hyd, &
601 coriolis, rtot, cvtot, cptot, dphyddx, dphyddy, &
602 element3d_operation, dx, dy, dz, sx, sy, sz, lift, &
603 lmesh, elem, lmesh2d, elem2d )
607 use scale_const,
only: &
616 type(
sparsemat),
intent(in) :: dx, dy, dz, sx, sy, sz, lift
617 real(rp),
intent(out) :: dens_dt(elem%np,lmesh%nea)
618 real(rp),
intent(out) :: momx_dt(elem%np,lmesh%nea)
619 real(rp),
intent(out) :: momy_dt(elem%np,lmesh%nea)
620 real(rp),
intent(out) :: momz_dt(elem%np,lmesh%nea)
621 real(rp),
intent(out) :: rhot_dt(elem%np,lmesh%nea)
622 real(rp),
intent(in) :: ddens_(elem%np,lmesh%nea)
623 real(rp),
intent(in) :: momx_(elem%np,lmesh%nea)
624 real(rp),
intent(in) :: momy_(elem%np,lmesh%nea)
625 real(rp),
intent(in) :: momz_(elem%np,lmesh%nea)
626 real(rp),
intent(in) :: drhot_(elem%np,lmesh%nea)
627 real(rp),
intent(in) :: dpres_(elem%np,lmesh%nea)
628 real(rp),
intent(in) :: dens_hyd(elem%np,lmesh%nea)
629 real(rp),
intent(in) :: pres_hyd(elem%np,lmesh%nea)
630 real(rp),
intent(in) :: pres_hyd_ref(elem%np,lmesh%nea)
631 real(rp),
intent(in) :: therm_hyd(elem%np,lmesh%nea)
632 real(rp),
intent(in) :: coriolis(elem2d%np,lmesh2d%nea)
633 real(rp),
intent(in) :: rtot (elem%np,lmesh%nea)
634 real(rp),
intent(in) :: cvtot(elem%np,lmesh%nea)
635 real(rp),
intent(in) :: cptot(elem%np,lmesh%nea)
636 real(rp),
intent(in) :: dphyddx(elem%np,lmesh%nea)
637 real(rp),
intent(in) :: dphyddy(elem%np,lmesh%nea)
639 real(rp) :: fx(elem%np), fy(elem%np), fz(elem%np), liftdelflx(elem%np)
640 real(rp) :: dpres_hyd(elem%np), gradphyd_x(elem%np), gradphyd_y(elem%np)
641 real(rp) :: del_flux(elem%nfptot,lmesh%ne,
prgvar_num)
642 real(rp) :: del_flux_hyd(elem%nfptot,lmesh%ne,2)
643 real(rp) :: rhot_(elem%np)
644 real(rp) :: rdens_(elem%np), u_(elem%np), v_(elem%np), w_(elem%np), wt_(elem%np), drho(elem%np)
646 real(rp) :: g11(elem%np), g12(elem%np), g22(elem%np)
647 real(rp) :: gsqrtv(elem%np), rgsqrtv(elem%np), rgam2(elem%np)
648 real(rp) :: x2d(elem2d%np,lmesh2d%ne), y2d(elem2d%np,lmesh2d%ne)
649 real(rp) :: x(elem%np), y(elem%np), twoovdel2(elem%np)
650 real(rp) :: om1(elem%np), om2(elem%np), om3(elem%np), del(elem%np), r(elem%np)
651 logical :: is_panel1to4
662 call prof_rapstart(
'cal_dyn_tend_bndflux', 3)
663 call get_ebnd_flux( &
664 del_flux, del_flux_hyd, &
665 ddens_, momx_, momy_, momz_, drhot_, dpres_, dens_hyd, pres_hyd, &
666 rtot, cvtot, cptot, &
667 lmesh%Gsqrt, lmesh%GIJ(:,:,1,1), lmesh%GIJ(:,:,1,2), lmesh%GIJ(:,:,2,2), &
668 lmesh%GsqrtH, lmesh%gam, lmesh%GI3(:,:,1), lmesh%GI3(:,:,2), &
669 lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), &
670 lmesh%vmapM, lmesh%vmapP, elem%IndexH2Dto3D_bnd, &
671 lmesh, elem, lmesh2d, elem2d )
672 call prof_rapend(
'cal_dyn_tend_bndflux', 3)
675 call prof_rapstart(
'cal_dyn_tend_interior', 3)
676 rgamm = cvdry / cpdry
677 rp0 = 1.0_rp / pres00
678 p0ovr = pres00 / rdry
681 is_panel1to4 = .true.
682 if ( lmesh%panelID == 5 )
then
683 is_panel1to4 = .false.
684 else if ( lmesh%panelID == 6 )
then
685 is_panel1to4 = .false.
698 do ke2d = lmesh2d%NeS, lmesh2d%NeE
699 x2d(:,ke2d) = tan(lmesh2d%pos_en(:,ke2d,1))
700 y2d(:,ke2d) = tan(lmesh2d%pos_en(:,ke2d,2))
704 do ke = lmesh%NeS, lmesh%NeE
706 ke2d = lmesh%EMap3Dto2D(ke)
707 rgam2(:) = 1.0_rp / lmesh%gam(:,ke)**2
708 g11(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,1,1) * rgam2(:)
709 g12(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,1,2) * rgam2(:)
710 g22(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,2,2) * rgam2(:)
711 gsqrtv(:) = lmesh%Gsqrt(:,ke) * rgam2(:) / lmesh%GsqrtH(elem%IndexH2Dto3D,ke2d)
712 rgsqrtv(:) = 1.0_rp / gsqrtv(:)
715 rhot_(:) = p0ovr * ( pres_hyd(:,ke) * rp0 )**rgamm + drhot_(:,ke)
719 rdens_(:) = 1.0_rp / ( ddens_(:,ke) + dens_hyd(:,ke) )
720 u_(:) = momx_(:,ke) * rdens_(:)
721 v_(:) = momy_(:,ke) * rdens_(:)
722 w_(:) = momz_(:,ke) * rdens_(:)
723 wt_(:) = w_(:) * rgsqrtv(:) + lmesh%GI3(:,ke,1) * u_(:) + lmesh%GI3(:,ke,2) * v_(:)
725 x(:) = x2d(elem%IndexH2Dto3D,ke2d)
726 y(:) = y2d(elem%IndexH2Dto3D,ke2d)
727 del(:) = sqrt( 1.0_rp + x(:)**2 + y(:)**2 )
728 twoovdel2(:) = 2.0_rp / ( 1.0_rp + x(:)**2 + y(:)**2 )
730 r(:) = rplanet * lmesh%gam(:,ke)
734 if ( is_panel1to4 )
then
736 om2(:) = s * del(:) / ( r(:) * ( 1.0_rp + y(:)**2 ) )
737 om3(:) = s * y(:) / del(:)
739 om1(:) = - s * x(:) * del(:) / ( r(:) * ( 1.0_rp + x(:)**2 ) )
740 om2(:) = - s * y(:) * del(:) / ( r(:) * ( 1.0_rp + y(:)**2 ) )
748 dpres_hyd(:) = pres_hyd(:,ke) - pres_hyd_ref(:,ke)
752 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux_hyd(:,ke,1), liftdelflx)
753 gradphyd_x(:) = lmesh%Escale(:,ke,1,1) * fx(:) &
754 + lmesh%Escale(:,ke,3,3) * fz(:) &
759 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux_hyd(:,ke,2), liftdelflx)
760 gradphyd_y(:) = lmesh%Escale(:,ke,2,2) * fy(:) &
761 + lmesh%Escale(:,ke,3,3) * fz(:) &
767 call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( ddens_(:,ke) + dens_hyd(:,ke) ) * wt_(:), fz)
768 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,dens_vid), liftdelflx)
770 dens_dt(:,ke) = - ( &
771 lmesh%Escale(:,ke,1,1) * fx(:) &
772 + lmesh%Escale(:,ke,2,2) * fy(:) &
773 + lmesh%Escale(:,ke,3,3) * fz(:) &
774 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke)
777 call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * ( u_(:) * momx_(:,ke) + g11(:) * dpres_(:,ke) ), fx)
778 call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * ( v_(:) * momx_(:,ke) + g12(:) * dpres_(:,ke) ), fy)
780 + ( lmesh%GI3(:,ke,1) * g11(:) + lmesh%GI3(:,ke,2) * g12(:) ) * dpres_(:,ke) ), fz)
781 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momx_vid), liftdelflx)
784 - ( lmesh%Escale(:,ke,1,1) * fx(:) &
785 + lmesh%Escale(:,ke,2,2) * fy(:) &
786 + lmesh%Escale(:,ke,3,3) * fz(:) &
787 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke) &
788 - twoovdel2(:) * y(:) * &
789 ( x(:) * y(:) * u_(:) - ( 1.0_rp + y(:)**2 ) * v_(:) ) * momx_(:,ke) &
790 - 2.0_rp * u_(:) * momz_(:,ke) / r(:) &
791 - ( g11(:) * gradphyd_x(:) + g12(:) * gradphyd_y(:) ) * rgsqrtv(:) &
792 - lmesh%Gsqrt(:,ke) * ( g11(:) * ( om2(:) * momz_(:,ke) - om3(:) * momy_(:,ke) ) &
793 - g12(:) * ( om1(:) * momz_(:,ke) - om3(:) * momx_(:,ke) ) )
796 call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * ( u_(:) * momy_(:,ke) + g12(:) * dpres_(:,ke) ), fx)
797 call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * ( v_(:) * momy_(:,ke) + g22(:) * dpres_(:,ke) ), fy)
799 + ( lmesh%GI3(:,ke,1) * g12(:) + lmesh%GI3(:,ke,2) * g22(:) ) * dpres_(:,ke) ), fz)
800 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momy_vid), liftdelflx)
803 - ( lmesh%Escale(:,ke,1,1) * fx(:) &
804 + lmesh%Escale(:,ke,2,2) * fy(:) &
805 + lmesh%Escale(:,ke,3,3) * fz(:) &
806 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke) &
807 - twoovdel2(:) * x(:) * &
808 ( - (1.0_rp + x(:)**2) * u_(:) + x(:) * y(:) * v_(:) ) * momy_(:,ke) &
809 - 2.0_rp * v_(:) * momz_(:,ke) / r(:) &
810 - ( g12(:) * gradphyd_x(:) + g22(:) * gradphyd_y(:) ) * rgsqrtv(:) &
811 - lmesh%Gsqrt(:,ke) * ( g12(:) * ( om2(:) * momz_(:,ke) - om3(:) * momy_(:,ke) ) &
812 - g22(:) * ( om1(:) * momz_(:,ke) - om3(:) * momx_(:,ke) ) )
817 call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( wt_(:) * momz_(:,ke) + rgsqrtv(:) * dpres_(:,ke) ), fz)
818 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momz_vid), liftdelflx)
821 - ( lmesh%Escale(:,ke,1,1) * fx(:) &
822 + lmesh%Escale(:,ke,2,2) * fy(:) &
823 + lmesh%Escale(:,ke,3,3) * fz(:) &
824 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke) &
825 - 0.25_rp * r(:) * twoovdel2(:)**2 * ( 1.0_rp * x(:)**2 ) * ( 1.0_rp * y(:)**2 ) &
826 * ( - ( 1.0_rp + x(:)**2 ) * momx_(:,ke) * u_(:) &
827 + 2.0_rp * x(:) * y(:) * momx_(:,ke) * v_(:) &
828 - ( 1.0_rp + y(:)**2 ) * momy_(:,ke) * v_(:) ) &
829 + 2.0_rp * dpres_(:,ke) / r(:) &
830 - lmesh%Gsqrt(:,ke) * ( om1(:) * momy_(:,ke) - om2(:) * momx_(:,ke) ) &
831 - grav * rgam2(:) * drho(:)
837 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,rhot_vid), liftdelflx)
840 - ( lmesh%Escale(:,ke,1,1) * fx(:) &
841 + lmesh%Escale(:,ke,2,2) * fy(:) &
842 + lmesh%Escale(:,ke,3,3) * fz(:) &
843 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke)
847 call prof_rapend(
'cal_dyn_tend_interior', 3)