Calculate tendency with PBL turbulence models.
65 implicit none
66 class(LocalMesh3D), intent(in), target :: lmesh
67 class(ElementBase3D), intent(in) :: elem
68 class(ElementBase1D), intent(in) :: elem1D
69 real(RP), intent(out) :: RHOU_tp(elem%Np,lmesh%NeA)
70 real(RP), intent(out) :: RHOV_tp(elem%Np,lmesh%NeA)
71 real(RP), intent(out) :: DRHOT_tp(elem%Np,lmesh%NeA)
72 type(LocalMeshFieldBaseList), intent(inout) :: RHOQ_tp_list(QA)
73 real(RP), intent(in) :: DDENS_(elem%Np,lmesh%NeA)
74 real(RP), intent(in) :: MOMX_(elem%Np,lmesh%NeA)
75 real(RP), intent(in) :: MOMY_(elem%Np,lmesh%NeA)
76 real(RP), intent(in) :: DRHOT_(elem%Np,lmesh%NeA)
77 type(LocalMeshFieldBaseList), intent(in) :: QTRC_list(QA)
78 real(RP), intent(in) :: PT_(elem%Np,lmesh%NeA)
79 real(RP), intent(in) :: DENS_hyd(elem%Np,lmesh%NeA)
80 real(RP), intent(in) :: PRES_hyd(elem%Np,lmesh%NeA)
81 real(RP), intent(in) :: NU(elem%Np,lmesh%NeA)
82 real(RP), intent(in) :: KH(elem%Np,lmesh%NeA)
83 class(ElementOperationBase3D), intent(in) :: element3D_operation
84 real(RP), intent(in) :: C_IP
85 real(RP), intent(in) :: dtsec
86 logical, intent(in) :: is_bound(elem%NfpTot,lmesh%Ne)
87 logical, intent(in) :: use_delta_form
88
89 class(LocalMesh2D), pointer :: lmesh2D
90 class(ElementBase2D), pointer :: elem2D
91
92 integer :: iq
93 real(RP) :: QTRC00_(elem%Np,QA,lmesh%Ne)
94
95 real(RP) :: PROG_VARS (elem%Np,lmesh%NeX*lmesh%NeY,lmesh%NeZ,3+QA)
96 real(RP) :: alph_M(elem%NfpTot,lmesh%Ne)
97 real(RP) :: alph_H(elem%NfpTot,lmesh%Ne)
98 real(RP) :: GsqrtV(elem%Np,lmesh%Ne)
99
100 integer :: vmapM(elem%NfpTot,lmesh%Ne)
101 integer :: vmapP(elem%NfpTot,lmesh%Ne)
102 integer :: ke_xy, ke_z, ke, ke2d
103 integer :: p
104
105 integer :: im, jm
106
107 real(RP) :: DENS(elem%Np,lmesh%Ne)
108 real(RP), allocatable :: b1D_ij(:,:,:,:,:)
109 real(RP), allocatable :: BndMatL(:,:,:,:,:)
110 real(RP), allocatable :: BndMatD(:,:,:,:,:)
111 real(RP), allocatable :: G(:,:,:,:,:,:)
112
113 real(RP) :: impl_fac, r_impl_fac
114
115
116 lmesh2d => lmesh%lcmesh2D
117 elem2d => lmesh2d%refElem2D
118 impl_fac = 1.0_rp * dtsec
119
120 call lmesh%GetVmapZ3D( vmapm, vmapp )
122 elem%Nnode_v )
123
124 allocate( b1d_ij(im*elem%Nnode_v,3+qa,jm,lmesh%Ne2D,lmesh%NeZ) )
125 allocate( bndmatd(im*elem%Nnode_v,elem%Nnode_v,jm,lmesh%Ne2D,2) )
126 allocate( bndmatl(im*elem%Nnode_v,elem%Nnode_v,jm,lmesh%Ne2D,2) )
127 allocate( g(im*elem%Nnode_v,elem%Nnode_v,jm,lmesh%Ne2D,lmesh%NeZ,2) )
128
129
130
131 do iq = 1, qa
132 do ke=lmesh%NeS, lmesh%NeE
133 qtrc00_(:,iq,ke) = qtrc_list(iq)%ptr%val(:,ke)
134 end do
135 end do
136
137
138
139 do ke_z =1, lmesh%NeZ
140 do ke_xy=1, lmesh%NeX * lmesh%NeY
141 ke = ke_xy + (ke_z-1)*lmesh%Ne2D
142 ke2d = lmesh%EMap3Dto2D(ke)
143
144 dens(:,ke) = dens_hyd(:,ke) + ddens_(:,ke)
145
146 prog_vars(:,ke_xy,ke_z,1) = momx_(:,ke)
147 prog_vars(:,ke_xy,ke_z,2) = momy_(:,ke)
148 prog_vars(:,ke_xy,ke_z,3) = dens(:,ke) * pt_(:,ke)
149 do iq = 1, qa
150 prog_vars(:,ke_xy,ke_z,3+iq) = dens(:,ke) * qtrc00_(:,iq,ke)
151 end do
152
153 do p=1, elem%Np
154 gsqrtv(p,ke) = lmesh%Gsqrt(p,ke) / lmesh%GsqrtH(elem%IndexH2Dto3D(p),ke2d)
155 end do
156 end do
157 end do
158
159
160
161 call eval_ax( rhou_tp, rhov_tp, drhot_tp, rhoq_tp_list, alph_m, alph_h, &
162 prog_vars, momx_, momy_, pt_, qtrc00_, nu, kh, dens, gsqrtv, &
163 impl_fac, dtsec, lmesh, elem, vmapm, vmapp, is_bound, &
164 element3d_operation, c_ip, im, jm, b1d_ij, use_delta_form )
165
166 call vi_solve( prog_vars, &
167 bndmatl, bndmatd, g, b1d_ij, &
168 dens, nu, kh, gsqrtv, c_ip, dtsec, impl_fac, &
169 im, jm, lmesh, elem, elem1d, use_delta_form )
170
171
172 r_impl_fac = 1.0_rp / impl_fac
173
174
175 do ke_z =1, lmesh%NeZ
176 do ke_xy=1, lmesh%NeX * lmesh%NeY
177 ke = ke_xy + (ke_z-1)*lmesh%Ne2D
178 rhou_tp(:,ke) = ( prog_vars(:,ke_xy,ke_z,1) - momx_(:,ke) ) * r_impl_fac
179 rhov_tp(:,ke) = ( prog_vars(:,ke_xy,ke_z,2) - momy_(:,ke) ) * r_impl_fac
180 drhot_tp(:,ke) = ( prog_vars(:,ke_xy,ke_z,3) - dens(:,ke) * pt_(:,ke) ) * r_impl_fac
181
182 do iq=1, qa
183 rhoq_tp_list(iq)%ptr%val(:,ke) = ( prog_vars(:,ke_xy,ke_z,3+iq) - dens(:,ke) * qtrc00_(:,iq,ke) ) * r_impl_fac
184 end do
185 end do
186 end do
187
188 return
module FElib / Fluid dyn solver / Atmosphere / HEVI / Common
subroutine, public atm_dyn_dgm_hevi_common_linalgebra_get_param(im, jm, nnode_h1d)