! 3 and 2 dimensional Euler solver - written in general coordinates ! Spatial Schemes- MP5, WENO5, WENO-Z, IG4, BVD algorithm ! Riemann flux estimate - HLLC ! Written by Hemanth Chandra Vamsi K and Ameraswara Sainadh Ch.[29th April 2021] PROGRAM MAIN_3D_MB_Gen_Coords use declare_variables ! Pre_processing.f90 routines call INPUTS() !--> nbocks = 1, NI(1) = 400, NJ(1) = 400, NK(1) = 1 call ALLOCATE_VARIABLES() call COEFFICENTS() call GENERATE_GRID() call CONSERVATIVE_METRICS() call INTERPOLATE_METRICS_ph() call INITIALIZE_NON_DIMENSIONALIZE() ! Main time loop [where parallelization is needed] call CPU_TIME(start_time) !$acc data copyin(NI,NJ,NK,nblocks,nprims,nconserv), & !$acc copyin(CFL,f2D,fViscous), & !$acc copyin(RK_factor1,RK_factor2,RK_factor3), & !$acc copyin(gamma,Pr,Re), & !$acc copyin(Ix,Iy,Iz,Jx,Jy,Jz,Kx,Ky,Kz,Jac), & !$acc copyin(Ix_iph,Iy_iph,Iz_iph,Jx_iph,Jy_iph,Jz_iph,Kx_iph,Ky_iph,Kz_iph,Jac_iph), & !$acc copyin(Ix_c,Iy_c,Iz_c), & !$acc copyin(Ix_jph,Iy_jph,Iz_jph,Jx_jph,Jy_jph,Jz_jph,Kx_jph,Ky_jph,Kz_jph,Jac_jph), & !$acc copyin(Jx_c,Jy_c,Jz_c), & !$acc copyin(Ix_kph,Iy_kph,Iz_kph,Jx_kph,Jy_kph,Jz_kph,Kx_kph,Ky_kph,Kz_kph,Jac_kph), & !$acc copyin(Kx_c,Ky_c,Kz_c), & !$acc copyin(Qp,Qc,Qc_initial), & !$acc copyin(Qp_iphL,Qp_iphR,Qp_jphL,Qp_jphR,Qp_kphL,Qp_kphR), & !$acc copyin(Residual_RHS), & !$acc copyin(aa_tdmaI,bb_tdmaI,cc_tdmaI,aa_tdmaJ,bb_tdmaJ,cc_tdmaJ,aa_tdmaK,bb_tdmaK,cc_tdmaK), & !$acc copyin(aa_tdmaI2,bb_tdmaI2,cc_tdmaI2,aa_tdmaJ2,bb_tdmaJ2,cc_tdmaJ2,aa_tdmaK2,bb_tdmaK2,cc_tdmaK2), & !$acc copyin(Qp_I,Qp_J,Qp_K),& !$acc copyin(X_CC,Y_CC) ! DO WHILE (time < t_end) DO WHILE (iter < 2) call RUNTIME_ROUTINE() Do step = 1, rk_steps call BOUNDARY_CONDITIONS() call MP5_FACE_INTERPOLATION_WITH_PRIMS_CHARS3() call HLL_INTERCELL_FFLUX_ESTIMATE_I() call HLL_INTERCELL_GFLUX_ESTIMATE_J() IF (f2D .ne. 1) call HLL_INTERCELL_HFLUX_ESTIMATE_K() IF (fViscous == 1) call VISCOUS_FLUXES() !$acc parallel loop gang collapse(4) default(present) DO n_cons = 1,nconserv DO nbl = 1,nblocks DO k = 1, NKmax DO j = 1, NJmax !$acc loop vector DO i = 1, NImax if (k.le.NK(nbl).and.j.le.NJ(nbl).and.i.le.NI(nbl)) then Qc(i,j,k,nbl,n_cons) = RK_factor1(step)*Qc_initial(i,j,k,nbl,n_cons) & + RK_factor2(step)*Qc(i,j,k,nbl,n_cons) & + RK_factor3(step)*Residual_RHS(i,j,k,nbl,n_cons)*time_step*Jac(i,j,k,nbl) endif ENDDO ENDDO ENDDO ENDDO ENDDO IF (fPositivity == 1) call POSITIVITY_PRESERVING() call SET_PRIMITIVES() ENDDO ENDDO !$acc end data call CPU_TIME(end_time) Print*, 'Total wall clock time = ', (end_time-start_time), 'secs' ! Post-processing call OUTPUT(1) print*, 'SOLUTION - (at domain center)' print*, Qp(NImax/2, NJmax/2, 1,1,1) print*, Qp(NImax/2, NJmax/2, 1,1,2) print*, Qp(NImax/2, NJmax/2, 1,1,3) print*, Qp(NImax/2, NJmax/2, 1,1,4) print*, Qp(NImax/2, NJmax/2, 1,1,5) END SUBROUTINE RUNTIME_ROUTINE() use declare_variables implicit none double precision :: U_contra, V_contra, W_contra, Cx, Cy, Cz, mu_L, vx, vy, vz, T !================================================ IF(MOD(iter,5000) == 0) call OUTPUT(1) !$acc parallel loop gang collapse(4) default(present) DO n_cons = 1,nconserv DO nbl = 1,nblocks DO k = 1, NKmax DO j = 1, NJmax !$acc loop vector DO i = 1, NImax if (k.le.NK(nbl).and.j.le.NJ(nbl).and.i.le.NI(nbl)) then Qc_initial(i,j,k,nbl,n_cons) = Qc(i,j,k,nbl,n_cons) endif ENDDO ENDDO ENDDO ENDDO ENDDO time = time + time_step !================================================ IF (iter == 0) time_step = 1000.d0 IF (fViscous == 1) THEN mu_L = 1.d0/Re !$acc parallel loop gang collapse(3) reduction(min:time_step) default(present) DO nbl = 1,nblocks DO k = 1, NKmax DO j = 1, NJmax !$acc loop vector DO i = 1, NImax if (k.le.NK(nbl).and.j.le.NJ(nbl).and.i.le.NI(nbl)) then ! Temperature and viscosity (n0n -dimensional) ! T = Gamma*Mach**2*Qp(i,j,k,nbl,5)/Qp(i,j,k,nbl,1) ! mu_L = (T**1.5d0)*(1.d0+110.5d0/T_inf)/(T+110.5d0/T_inf)/Re vx = mu_L*(Ix(i,j,k,nbl)**2+Iy(i,j,k,nbl)**2+Iz(i,j,k,nbl)**2) vy = mu_L*(Jx(i,j,k,nbl)**2+Jy(i,j,k,nbl)**2+Jz(i,j,k,nbl)**2) vz = mu_L*(Kx(i,j,k,nbl)**2+Ky(i,j,k,nbl)**2+Kz(i,j,k,nbl)**2) Cx = SQRT(Gamma*Qp(i,j,k,nbl,5)/Qp(i,j,k,nbl,1)*(Ix(i,j,k,nbl)**2+Iy(i,j,k,nbl)**2+Iz(i,j,k,nbl)**2)) Cy = SQRT(Gamma*Qp(i,j,k,nbl,5)/Qp(i,j,k,nbl,1)*(Jx(i,j,k,nbl)**2+Jy(i,j,k,nbl)**2+Jz(i,j,k,nbl)**2)) Cz = SQRT(Gamma*Qp(i,j,k,nbl,5)/Qp(i,j,k,nbl,1)*(Kx(i,j,k,nbl)**2+Ky(i,j,k,nbl)**2+Kz(i,j,k,nbl)**2)) U_contra = Qp(i,j,k,nbl,2)*Ix(i,j,k,nbl) + Qp(i,j,k,nbl,3)*Iy(i,j,k,nbl) + Qp(i,j,k,nbl,4)*Iz(i,j,k,nbl) V_contra = Qp(i,j,k,nbl,2)*Jx(i,j,k,nbl) + Qp(i,j,k,nbl,3)*Jy(i,j,k,nbl) + Qp(i,j,k,nbl,4)*Jz(i,j,k,nbl) W_contra = Qp(i,j,k,nbl,2)*Kx(i,j,k,nbl) + Qp(i,j,k,nbl,3)*Ky(i,j,k,nbl) + Qp(i,j,k,nbl,4)*Kz(i,j,k,nbl) time_step = MIN(time_step, CFL/MAX( ABS(U_contra)+Cx, ABS(V_contra)+Cy, ABS(W_contra)+Cz), CFL/MAX(vx,vy,vz)) endif ENDDO ENDDO ENDDO ENDDO ELSE !$acc parallel loop gang collapse(3) reduction(min:time_step) default(present) DO nbl = 1,nblocks DO k = 1, NKmax DO j = 1, NJmax !$acc loop vector DO i = 1, NImax if (k.le.NK(nbl).and.j.le.NJ(nbl).and.i.le.NI(nbl)) then Cx = SQRT(Gamma*Qp(i,j,k,nbl,5)/Qp(i,j,k,nbl,1)*(Ix(i,j,k,nbl)**2+Iy(i,j,k,nbl)**2+Iz(i,j,k,nbl)**2)) Cy = SQRT(Gamma*Qp(i,j,k,nbl,5)/Qp(i,j,k,nbl,1)*(Jx(i,j,k,nbl)**2+Jy(i,j,k,nbl)**2+Jz(i,j,k,nbl)**2)) Cz = SQRT(Gamma*Qp(i,j,k,nbl,5)/Qp(i,j,k,nbl,1)*(Kx(i,j,k,nbl)**2+Ky(i,j,k,nbl)**2+Kz(i,j,k,nbl)**2)) U_contra = Qp(i,j,k,nbl,2)*Ix(i,j,k,nbl) + Qp(i,j,k,nbl,3)*Iy(i,j,k,nbl) + Qp(i,j,k,nbl,4)*Iz(i,j,k,nbl) V_contra = Qp(i,j,k,nbl,2)*Jx(i,j,k,nbl) + Qp(i,j,k,nbl,3)*Jy(i,j,k,nbl) + Qp(i,j,k,nbl,4)*Jz(i,j,k,nbl) W_contra = Qp(i,j,k,nbl,2)*Kx(i,j,k,nbl) + Qp(i,j,k,nbl,3)*Ky(i,j,k,nbl) + Qp(i,j,k,nbl,4)*Kz(i,j,k,nbl) time_step = MIN(time_step, CFL/MAX( ABS(U_contra)+Cx, ABS(V_contra)+Cy, ABS(W_contra)+Cz)) endif ENDDO ENDDO ENDDO ENDDO ENDIF IF (time + time_step > t_end) time_step = t_end-time iter = iter +1 WRITE(*,1) iter, time, time_step 1 format(i8,f8.4,f12.8) IF (time + time_step=0.d0) THEN FFlux_iph(i,j,k,nbl,1) = (rho_L*U_contra_L) /Jac_iph(i,j,k,nbl) FFlux_iph(i,j,k,nbl,2) = (rho_L*U_contra_L*Ux_L + P_L*Ix_iph(i,j,k,nbl)) /Jac_iph(i,j,k,nbl) FFlux_iph(i,j,k,nbl,3) = (rho_L*U_contra_L*Uy_L + P_L*Iy_iph(i,j,k,nbl)) /Jac_iph(i,j,k,nbl) FFlux_iph(i,j,k,nbl,4) = ((rho_L*U_contra_L*Uz_L + P_L*Iz_iph(i,j,k,nbl)) /Jac_iph(i,j,k,nbl))*(1-f2D) FFlux_iph(i,j,k,nbl,5) = (E_L + P_L) * U_contra_L /Jac_iph(i,j,k,nbl) ELSEIF(S_L <= 0.d0 .and. S_R>=0.d0) THEN FFlux_iph(i,j,k,nbl,1) = (S_R*(rho_L*U_contra_L) /Jac_iph(i,j,k,nbl) - S_L*(rho_R*U_contra_R) /Jac_iph(i,j,k,nbl) & + S_L*S_R*(Rho_R - Rho_L)/Jac_iph(i,j,k,nbl))/(S_R - S_L) FFlux_iph(i,j,k,nbl,2) = (S_R*(rho_L*U_contra_L*Ux_L + P_L*Ix_iph(i,j,k,nbl)) /Jac_iph(i,j,k,nbl) - & S_L*(rho_R*U_contra_R*Ux_R + P_R*Ix_iph(i,j,k,nbl)) /Jac_iph(i,j,k,nbl) + & S_L*S_R*(Rho_R*Ux_R - Rho_L*Ux_L)/Jac_iph(i,j,k,nbl))/(S_R - S_L) FFlux_iph(i,j,k,nbl,3) = (S_R*(rho_L*U_contra_L*Uy_L + P_L*Iy_iph(i,j,k,nbl)) /Jac_iph(i,j,k,nbl) & - S_L*(rho_R*U_contra_R*Uy_R + P_R*Iy_iph(i,j,k,nbl)) /Jac_iph(i,j,k,nbl) + S_L*S_R*(Rho_R*Uy_R - & Rho_L*Uy_L)/Jac_iph(i,j,k,nbl))/(S_R - S_L) FFlux_iph(i,j,k,nbl,4) =((S_R*((rho_L*U_contra_L*Uz_L + P_L*Iz_iph(i,j,k,nbl)) /Jac_iph(i,j,k,nbl))*(1-f2D)& - S_L*((rho_R*U_contra_R*Uz_R + P_R*Iz_iph(i,j,k,nbl)) /Jac_iph(i,j,k,nbl))*(1-f2D) + S_L*S_R*(Rho_R*Uz_R - & Rho_L*Uz_L)/Jac_iph(i,j,k,nbl))/(S_R - S_L))*(1-f2D) FFlux_iph(i,j,k,nbl,5) = (S_R*(E_L + P_L) * U_contra_L /Jac_iph(i,j,k,nbl) - S_L*(E_R + P_R) * U_contra_R & /Jac_iph(i,j,k,nbl) + S_L*S_R*(E_R - E_L)/Jac_iph(i,j,k,nbl))/(S_R - S_L) ELSEIF(S_R<=0.d0) THEN FFlux_iph(i,j,k,nbl,1) = (rho_R*U_contra_R) /Jac_iph(i,j,k,nbl) FFlux_iph(i,j,k,nbl,2) = (rho_R*U_contra_R*Ux_R + P_R*Ix_iph(i,j,k,nbl)) /Jac_iph(i,j,k,nbl) FFlux_iph(i,j,k,nbl,3) = (rho_R*U_contra_R*Uy_R + P_R*Iy_iph(i,j,k,nbl)) /Jac_iph(i,j,k,nbl) FFlux_iph(i,j,k,nbl,4) = ((rho_R*U_contra_R*Uz_R + P_R*Iz_iph(i,j,k,nbl)) /Jac_iph(i,j,k,nbl))*(1-f2D) FFlux_iph(i,j,k,nbl,5) = (E_R + P_R) * U_contra_R /Jac_iph(i,j,k,nbl) ENDIF endif ENDDO ENDDO ENDDO ENDDO !$acc parallel loop gang collapse(4) default(present) DO n_cons = 1,nconserv DO nbl = 1,nblocks DO k = 1, NKmax DO j = 1, NJmax !$acc loop vector DO i = 1, NImax if (k.le.NK(nbl).and.j.le.NJ(nbl).and.i.le.NI(nbl)) then Residual_RHS(i,j,k,nbl,n_cons) = -(Fflux_iph(i,j,k,nbl,n_cons) - Fflux_iph(i-1,j,k,nbl,n_cons)) endif ENDDO ENDDO ENDDO ENDDO ENDDO !$acc end data END SUBROUTINE SUBROUTINE HLL_INTERCELL_GFLUX_ESTIMATE_J() ! Task: Given primitave vars on the cell faces (Qp_iphR, Qp_iphL), find "FLUX ON EACH FACE" use declare_variables implicit none !Variables for HLL double precision, dimension(NImax, 0:NJmax, NKmax, nblocks, nconserv) :: Gflux_jph double precision :: Rho_L, Rho_R, Ux_L, Ux_R, Uy_L, Uy_R, Uz_L, Uz_R, P_L, P_R, C_L, C_R, E_L, E_R, h_L, h_R, E_local double precision :: V_contra_L, V_contra_R, S_L, S_R ! double precision, dimension(nconserv) :: F_L, F_R ! double precision :: F_L1, F_L2, F_L3, F_L4, F_L5, F_R1, F_R2, F_R3, F_R4, F_R5 double precision :: sqrt_rho_L,sqrt_rho_R,Rho,divisor,u_average,v_average,w_average,h_average,vcontra_ave,uvw_average,C_average !$acc data create(Gflux_jph) !$acc parallel loop gang collapse(3) default(present) DO nbl = 1,nblocks DO k = 1,NKmax DO j = 0,NJmax !$acc loop vector DO i = 1,NImax if (k.le.NK(nbl).and.j.le.NJ(nbl).and.i.le.NI(nbl)) then !Local variables Rho_L = Qp_jphL(i,j,k,nbl,1) Rho_R = Qp_jphR(i,j,k,nbl,1) Ux_L = Qp_jphL(i,j,k,nbl,2) Ux_R = Qp_jphR(i,j,k,nbl,2) Uy_L = Qp_jphL(i,j,k,nbl,3) Uy_R = Qp_jphR(i,j,k,nbl,3) Uz_L = Qp_jphL(i,j,k,nbl,4) Uz_R = Qp_jphR(i,j,k,nbl,4) P_L = Qp_jphL(i,j,k,nbl,5) P_R = Qp_jphR(i,j,k,nbl,5) C_L = Sqrt(Gamma*P_L/Rho_L*(Jx_jph(i,j,k,nbl)**2+Jy_jph(i,j,k,nbl)**2+Jz_jph(i,j,k,nbl)**2)) C_R = Sqrt(Gamma*P_R/Rho_R*(Jx_jph(i,j,k,nbl)**2+Jy_jph(i,j,k,nbl)**2+Jz_jph(i,j,k,nbl)**2)) E_L = P_L/(Gamma-1.d0) + Rho_L*(Ux_L**2+Uy_L**2+Uz_L**2)/2.d0 E_R = P_R/(Gamma-1.d0) + Rho_R*(Ux_R**2+Uy_R**2+Uz_R**2)/2.d0 h_L = (E_L + P_L)/Rho_L h_R = (E_R + P_R)/Rho_R V_contra_L = Ux_L*Jx_jph(i,j,k,nbl) + Uy_L*Jy_jph(i,j,k,nbl) + Uz_L*Jz_jph(i,j,k,nbl) V_contra_R = Ux_R*Jx_jph(i,j,k,nbl) + Uy_R*Jy_jph(i,j,k,nbl) + Uz_R*Jz_jph(i,j,k,nbl) ! Roe averaging sqrt_rho_L = sqrt(Rho_L) sqrt_rho_R = sqrt(Rho_R) Rho = sqrt(Rho_R/Rho_L)*Rho_L divisor = 1.d0/(sqrt_rho_R+sqrt_rho_L) h_average = ((h_L*sqrt_rho_L) + (h_R*sqrt_rho_R))*divisor u_average = ((Ux_L*sqrt_rho_L) + (Ux_R*sqrt_rho_R))*divisor v_average = ((Uy_L*sqrt_rho_L) + (Uy_R*sqrt_rho_R))*divisor w_average = ((Uz_L*sqrt_rho_L) + (Uz_R*sqrt_rho_R))*divisor vcontra_ave = u_average*Jx_jph(i,j,k,nbl) + v_average*Jy_jph(i,j,k,nbl) + w_average*Jz_jph(i,j,k,nbl) uvw_average = 0.5d0 * (u_average**2.d0 + v_average**2.d0 + w_average**2.d0) C_average = sqrt((gamma - 1.d0) * (h_average - uvw_average))*SQRT(Jx_jph(i,j,k,nbl)**2+Jy_jph(i,j,k,nbl)**2+Jz_jph(i,j,k,nbl)**2) ! Wave speeds S_L = MIN(V_contra_L - C_L, vcontra_ave - C_average) S_R = MAX(V_contra_R + C_R, vcontra_ave + C_average) !S_L = MIN(V_contra_L - C_L, V_contra_R - C_R) !S_R = MAX(V_contra_L + C_L, V_contra_R + C_R) IF (S_L>=0.d0) THEN GFlux_jph(i,j,k,nbl,1) = (rho_L*V_contra_L) /Jac_jph(i,j,k,nbl) GFlux_jph(i,j,k,nbl,2) = (rho_L*V_contra_L*Ux_L + P_L*Jx_jph(i,j,k,nbl)) /Jac_jph(i,j,k,nbl) GFlux_jph(i,j,k,nbl,3) = (rho_L*V_contra_L*Uy_L + P_L*Jy_jph(i,j,k,nbl)) /Jac_jph(i,j,k,nbl) GFlux_jph(i,j,k,nbl,4) = ((rho_L*V_contra_L*Uz_L + P_L*Jz_jph(i,j,k,nbl)) /Jac_jph(i,j,k,nbl))*(1-f2D) GFlux_jph(i,j,k,nbl,5) = (E_L + P_L) * V_contra_L /Jac_jph(i,j,k,nbl) ELSEIF(S_L <= 0.d0 .and. S_R>=0.d0) THEN GFlux_jph(i,j,k,nbl,1) = (S_R*(rho_L*V_contra_L) /Jac_jph(i,j,k,nbl) - S_L*(rho_R*V_contra_R) & /Jac_jph(i,j,k,nbl) + S_L*S_R*(Rho_R - Rho_L)/Jac_jph(i,j,k,nbl))/(S_R - S_L) GFlux_jph(i,j,k,nbl,2) = (S_R*(rho_L*V_contra_L*Ux_L + P_L*Jx_jph(i,j,k,nbl)) /Jac_jph(i,j,k,nbl) & - S_L*(rho_R*V_contra_R*Ux_R + P_R*Jx_jph(i,j,k,nbl)) /Jac_jph(i,j,k,nbl) + S_L*S_R*(Rho_R*Ux_R - & Rho_L*Ux_L)/Jac_jph(i,j,k,nbl))/(S_R - S_L) GFlux_jph(i,j,k,nbl,3) = (S_R*(rho_L*V_contra_L*Uy_L + P_L*Jy_jph(i,j,k,nbl)) /Jac_jph(i,j,k,nbl) & - S_L*(rho_R*V_contra_R*Uy_R + P_R*Jy_jph(i,j,k,nbl)) /Jac_jph(i,j,k,nbl) + S_L*S_R*(Rho_R*Uy_R - & Rho_L*Uy_L)/Jac_jph(i,j,k,nbl))/(S_R - S_L) GFlux_jph(i,j,k,nbl,4) =((S_R*((rho_L*V_contra_L*Uz_L + P_L*Jz_jph(i,j,k,nbl)) & /Jac_jph(i,j,k,nbl))*(1-f2D) - S_L*((rho_R*V_contra_R*Uz_R + P_R*Jz_jph(i,j,k,nbl)) & /Jac_jph(i,j,k,nbl))*(1-f2D) + S_L*S_R*(Rho_R*Uz_R - Rho_L*Uz_L)& /Jac_jph(i,j,k,nbl))/(S_R - S_L))*(1-f2D) GFlux_jph(i,j,k,nbl,5) = (S_R*(E_L + P_L) * V_contra_L /Jac_jph(i,j,k,nbl) - S_L*(E_R + P_R)& * V_contra_R /Jac_jph(i,j,k,nbl) + S_L*S_R*(E_R - E_L)/Jac_jph(i,j,k,nbl))/(S_R - S_L) ELSEIF(S_R<=0.d0) THEN GFlux_jph(i,j,k,nbl,1) = (rho_R*V_contra_R) /Jac_jph(i,j,k,nbl) GFlux_jph(i,j,k,nbl,2) = (rho_R*V_contra_R*Ux_R + P_R*Jx_jph(i,j,k,nbl)) /Jac_jph(i,j,k,nbl) GFlux_jph(i,j,k,nbl,3) = (rho_R*V_contra_R*Uy_R + P_R*Jy_jph(i,j,k,nbl)) /Jac_jph(i,j,k,nbl) GFlux_jph(i,j,k,nbl,4) = ((rho_R*V_contra_R*Uz_R + P_R*Jz_jph(i,j,k,nbl)) /Jac_jph(i,j,k,nbl))*(1-f2D) GFlux_jph(i,j,k,nbl,5) = (E_R + P_R) * V_contra_R /Jac_jph(i,j,k,nbl) ENDIF endif ENDDO ENDDO ENDDO ENDDO !$acc parallel loop gang collapse(4) default(present) DO n_cons = 1,nconserv DO nbl = 1,nblocks DO k = 1,NKmax DO j = 1,NJmax !$acc loop vector DO i = 1,NImax if (k.le.NK(nbl).and.j.le.NJ(nbl).and.i.le.NI(nbl)) then Residual_RHS(i,j,k,nbl,n_cons) = Residual_RHS(i,j,k,nbl,n_cons) - & (Gflux_jph(i,j,k,nbl,n_cons) - Gflux_jph(i,j-1,k,nbl,n_cons)) endif ENDDO ENDDO ENDDO ENDDO ENDDO !$acc end data END SUBROUTINE SUBROUTINE HLL_INTERCELL_HFLUX_ESTIMATE_K() ! Task: Given primitave vars on the cell faces (Qp_iphR, Qp_iphL), find "FLUX ON EACH FACE" use declare_variables implicit none double precision, dimension(NImax, NJmax, 0:NKmax, nblocks, nconserv) :: Hflux_kph double precision :: Rho_L, Rho_R, Ux_L, Ux_R, Uy_L, Uy_R, Uz_L, Uz_R, P_L, P_R, C_L, C_R, E_L, E_R, h_L, h_R, E_local double precision :: W_contra_L, W_contra_R, S_L, S_R ! double precision, dimension(nconserv) :: F_L, F_R ! double precision :: F_L1, F_L2, F_L3, F_L4, F_L5, F_R1, F_R2, F_R3, F_R4, F_R5 double precision :: sqrt_rho_L,sqrt_rho_R,Rho,divisor,u_average,v_average,w_average,h_average,wcontra_ave,uvw_average,C_average !$acc data create(Hflux_kph) !$acc parallel loop gang collapse(3) default(present) DO nbl = 1,nblocks DO k = 0,NKmax DO j = 1,NJmax !$acc loop vector DO i = 1,NImax if (k.le.NK(nbl).and.j.le.NJ(nbl).and.i.le.NI(nbl)) then !Local variables Rho_L = Qp_kphL(i,j,k,nbl,1) Rho_R = Qp_kphR(i,j,k,nbl,1) Ux_L = Qp_kphL(i,j,k,nbl,2) Ux_R = Qp_kphR(i,j,k,nbl,2) Uy_L = Qp_kphL(i,j,k,nbl,3) Uy_R = Qp_kphR(i,j,k,nbl,3) Uz_L = Qp_kphL(i,j,k,nbl,4) Uz_R = Qp_kphR(i,j,k,nbl,4) P_L = Qp_kphL(i,j,k,nbl,5) P_R = Qp_kphR(i,j,k,nbl,5) C_L = Sqrt(Gamma*P_L/Rho_L*(Kx_kph(i,j,k,nbl)**2+Ky_kph(i,j,k,nbl)**2+Kz_kph(i,j,k,nbl)**2)) C_R = Sqrt(Gamma*P_R/Rho_R*(Kx_kph(i,j,k,nbl)**2+Ky_kph(i,j,k,nbl)**2+Kz_kph(i,j,k,nbl)**2)) E_L = P_L/(Gamma-1.d0) + Rho_L*(Ux_L**2+Uy_L**2+Uz_L**2)/2.d0 E_R = P_R/(Gamma-1.d0) + Rho_R*(Ux_R**2+Uy_R**2+Uz_R**2)/2.d0 h_L = (E_L + P_L)/Rho_L h_R = (E_R + P_R)/Rho_R W_contra_L = Ux_L*Kx_kph(i,j,k,nbl) + Uy_L*Ky_kph(i,j,k,nbl) + Uz_L*Kz_kph(i,j,k,nbl) W_contra_R = Ux_R*Kx_kph(i,j,k,nbl) + Uy_R*Ky_kph(i,j,k,nbl) + Uz_R*Kz_kph(i,j,k,nbl) ! Roe averaging sqrt_rho_L = sqrt(Rho_L) sqrt_rho_R = sqrt(Rho_R) Rho = sqrt(Rho_R/Rho_L)*Rho_L divisor = 1.d0/(sqrt_rho_R+sqrt_rho_L) h_average = ((h_L*sqrt_rho_L) + (h_R*sqrt_rho_R))*divisor u_average = ((Ux_L*sqrt_rho_L) + (Ux_R*sqrt_rho_R))*divisor v_average = ((Uy_L*sqrt_rho_L) + (Uy_R*sqrt_rho_R))*divisor w_average = ((Uz_L*sqrt_rho_L) + (Uz_R*sqrt_rho_R))*divisor wcontra_ave = u_average*Kx_kph(i,j,k,nbl) + v_average*Ky_kph(i,j,k,nbl) + w_average*Kz_kph(i,j,k,nbl) uvw_average = 0.5d0 * (u_average**2.d0 + v_average**2.d0 + w_average**2.d0) C_average = sqrt((gamma-1.d0) * (h_average-uvw_average))*SQRT(Kx_kph(i,j,k,nbl)**2+Ky_kph(i,j,k,nbl)**2+Kz_kph(i,j,k,nbl)**2) ! Wave speeds S_L = MIN(W_contra_L - C_L, wcontra_ave - C_average) S_R = MAX(W_contra_R + C_R, wcontra_ave + C_average) !S_L = MIN(W_contra_L - C_L, W_contra_R - C_R) !S_R = MAX(W_contra_L + C_L, W_contra_R + C_R)1 IF (S_L>=0.d0) THEN Hflux_kph(i,j,k,nbl,1) = (rho_L*W_contra_L) /Jac_kph(i,j,k,nbl) Hflux_kph(i,j,k,nbl,2) = (rho_L*W_contra_L*Ux_L + P_L*Kx_kph(i,j,k,nbl)) /Jac_kph(i,j,k,nbl) Hflux_kph(i,j,k,nbl,3) = (rho_L*W_contra_L*Uy_L + P_L*Ky_kph(i,j,k,nbl)) /Jac_kph(i,j,k,nbl) Hflux_kph(i,j,k,nbl,4) = ((rho_L*W_contra_L*Uz_L + P_L*Kz_kph(i,j,k,nbl)) /Jac_kph(i,j,k,nbl))*(1-f2D) Hflux_kph(i,j,k,nbl,5) = (E_L + P_L) * W_contra_L /Jac_kph(i,j,k,nbl) ELSEIF(S_L <= 0.d0 .and. S_R>=0.d0) THEN Hflux_kph(i,j,k,nbl,1) = (S_R*(rho_L*W_contra_L) /Jac_kph(i,j,k,nbl) - S_L*(rho_R*W_contra_R) & /Jac_kph(i,j,k,nbl) + S_L*S_R*(Rho_R - Rho_L))/(S_R - S_L)/Jac_kph(i,j,k,nbl) Hflux_kph(i,j,k,nbl,2) = (S_R*(rho_L*W_contra_L*Ux_L + P_L*Kx_kph(i,j,k,nbl)) /Jac_kph(i,j,k,nbl)& - S_L*(rho_R*W_contra_R*Ux_R + P_R*Kx_kph(i,j,k,nbl)) /Jac_kph(i,j,k,nbl) + S_L*S_R*(Rho_R*Ux_R -& Rho_L*Ux_L))/(S_R - S_L)/Jac_kph(i,j,k,nbl) Hflux_kph(i,j,k,nbl,3) = (S_R*(rho_L*W_contra_L*Uy_L + P_L*Ky_kph(i,j,k,nbl)) /Jac_kph(i,j,k,nbl) & - S_L*(rho_R*W_contra_R*Uy_R + P_R*Ky_kph(i,j,k,nbl)) /Jac_kph(i,j,k,nbl) + S_L*S_R*(Rho_R*Uy_R - & Rho_L*Uy_L))/(S_R - S_L)/Jac_kph(i,j,k,nbl) Hflux_kph(i,j,k,nbl,4) =((S_R*((rho_L*W_contra_L*Uz_L + P_L*Kz_kph(i,j,k,nbl)) /Jac_kph(i,j,k,nbl))& *(1-f2D) - S_L*((rho_R*W_contra_R*Uz_R + P_R*Kz_kph(i,j,k,nbl)) /Jac_kph(i,j,k,nbl))*(1-f2D) + & S_L*S_R*(Rho_R*Uz_R - Rho_L*Uz_L))/(S_R - S_L)/Jac_kph(i,j,k,nbl))*(1-f2D) Hflux_kph(i,j,k,nbl,5) = (S_R*(E_L + P_L) * W_contra_L /Jac_kph(i,j,k,nbl) - S_L*(E_R + P_R) * & W_contra_R /Jac_kph(i,j,k,nbl) + S_L*S_R*(E_R - E_L))/(S_R - S_L)/Jac_kph(i,j,k,nbl) ELSEIF(S_R<=0.d0) THEN Hflux_kph(i,j,k,nbl,1) = (rho_R*W_contra_R) /Jac_kph(i,j,k,nbl) Hflux_kph(i,j,k,nbl,2) = (rho_R*W_contra_R*Ux_R + P_R*Kx_kph(i,j,k,nbl)) /Jac_kph(i,j,k,nbl) Hflux_kph(i,j,k,nbl,3) = (rho_R*W_contra_R*Uy_R + P_R*Ky_kph(i,j,k,nbl)) /Jac_kph(i,j,k,nbl) Hflux_kph(i,j,k,nbl,4) = ((rho_R*W_contra_R*Uz_R + P_R*Kz_kph(i,j,k,nbl)) /Jac_kph(i,j,k,nbl))*(1-f2D) Hflux_kph(i,j,k,nbl,5) = (E_R + P_R) * W_contra_R /Jac_kph(i,j,k,nbl) ENDIF endif ENDDO ENDDO ENDDO ENDDO !$acc parallel loop gang collapse(4) default(present) DO n_cons = 1,nconserv DO nbl = 1,nblocks DO k = 1,NKmax DO j = 1,NJmax !$acc loop vector DO i = 1,NImax if (k.le.NK(nbl).and.j.le.NJ(nbl).and.i.le.NI(nbl)) then Residual_RHS(i,j,k,nbl,n_cons) = Residual_RHS(i,j,k,nbl,n_cons) - & (Hflux_kph(i,j,k,nbl,n_cons) - Hflux_kph(i,j,k-1,nbl,n_cons)) endif ENDDO ENDDO ENDDO ENDDO ENDDO !$acc end data END SUBROUTINE !=============== ! VISCOUS FLUXES !=============== SUBROUTINE VISCOUS_FLUXES() use declare_variables implicit none double precision, dimension(1:NImax, 1:NJmax, 1:NKmax, nblocks, 1:4) :: Visc_Fflux double precision, dimension(1:NImax, 1:NJmax, 1:NKmax, nblocks, 1:4) :: Visc_Gflux double precision, dimension(1:NImax, 1:NJmax, 1:NKmax, nblocks, 1:4) :: Visc_Hflux double precision, dimension(1:NImax, 1:NJmax, 1:NKmax, nblocks, 1:4) :: ViscFlux_der double precision, parameter :: funGamma = 1.4d0/((1.4d0-1.d0)*0.71d0) double precision :: rx, ry, rz, ux, uy, uz, vx, vy, vz, wx, wy, wz, px, py, pz, bx, by, bz!, Tx, Ty, Tz double precision :: T,Blk_Visc_Term, tau_xx, tau_yy, tau_zz, tau_xy, tau_yz, tau_zx, mu_L!, Fourier_factor double precision :: var1, var2, var3, var4, var5, var6, var7, var8, var9, var10, var11, var12 double precision Ixx,Iyy,Izz,Jac1 double precision Jxx,Jyy,Jzz double precision Kxx,Kyy,Kzz double precision QpI1,QpI2,QpI3,QPI4,QPI5 double precision QpJ1,QpJ2,QpJ3,QPJ4,QPJ5 double precision QpK1,QpK2,QpK3,QPK4,QPK5 !$acc parallel loop gang collapse(4) default(present) DO n_prim = 1, nprims DO nbl = 1,nblocks DO k = 1,NKmax DO j = 1,NJmax !$acc loop vector DO i = 1,NImax if (k.le.(NK(nbl)).and.j.le.(NJ(nbl)) .and.i.le.(NI(nbl))) then Qp_I(i,j,k,nbl,n_prim) = (1.0d0/12.0d0)*Qp(i-2,j,k,nbl,n_prim)- & (2.0d0/3.0d0)*Qp(i-1,j,k,nbl,n_prim)+ & (0.d0)*Qp(i+0,j,k,nbl,n_prim)+ & (2.0d0/3.0d0)*Qp(i+1,j,k,nbl,n_prim)- & (1.0d0/12.0d0)*Qp(i+2,j,k,nbl,n_prim) Qp_J(i,j,k,nbl,n_prim) = (1.0d0/12.0d0)*Qp(i,j-2,k,nbl,n_prim)-& (2.0d0/3.0d0)*Qp(i,j-1,k,nbl,n_prim)+ & (0.d0)*Qp(i,j+0,k,nbl,n_prim)+ & (2.0d0/3.0d0)*Qp(i,j+1,k,nbl,n_prim)- & (1.0d0/12.0d0)*Qp(i,j+2,k,nbl,n_prim) Qp_K(i,j,k,nbl,n_prim) = (1.0d0/12.0d0)*Qp(i,j,k-2,nbl,n_prim)-& (2.0d0/3.0d0)*Qp(i,j,k-1,nbl,n_prim)+ & (0.d0)*Qp(i,j,k+0,nbl,n_prim)+ & (2.0d0/3.0d0)*Qp(i,j,k+1,nbl,n_prim)- & (1.0d0/12.0d0)*Qp(i,j,k+2,nbl,n_prim) endif ENDDO ENDDO ENDDO ENDDO ENDDO ! call DERIVATIVE_I5D_COMPACT_E6_CENTRAL(Qp,Qp_I,nprims) ! call DERIVATIVE_J5D_COMPACT_E6_CENTRAL(Qp,Qp_J,nprims) ! IF( f2D .ne. 1) call DERIVATIVE_J5D_COMPACT_E6_CENTRAL(Qp,Qp_K,nprims) !$acc data create(Visc_Fflux,Visc_Gflux,Visc_Hflux,ViscFlux_der) !acc loop seq DO nbl = 1,nblocks !$acc parallel loop gang collapse(2) default(present) DO k = 1, NK(nbl) DO j = 1, NJ(nbl) !$acc loop vector DO i = 1, NI(nbl) QpI1 = Qp_I(i,j,k,nbl,1) QpI2 = Qp_I(i,j,k,nbl,2) QpI3 = Qp_I(i,j,k,nbl,3) QpI4 = Qp_I(i,j,k,nbl,4) QpI5 = Qp_I(i,j,k,nbl,5) QpJ1 = Qp_J(i,j,k,nbl,1) QpJ2 = Qp_J(i,j,k,nbl,2) QpJ3 = Qp_J(i,j,k,nbl,3) QpJ4 = Qp_J(i,j,k,nbl,4) QpJ5 = Qp_J(i,j,k,nbl,5) QpK1 = Qp_K(i,j,k,nbl,1) QpK2 = Qp_K(i,j,k,nbl,2) QpK3 = Qp_K(i,j,k,nbl,3) QpK4 = Qp_K(i,j,k,nbl,4) QpK5 = Qp_K(i,j,k,nbl,5) Ixx = Ix(i,j,k,nbl) Jxx = Jx(i,j,k,nbl) Kxx = Kx(i,j,k,nbl) !chain rule :: to calculate local u_x, v_x,...... rx = Ixx*QpI1 + Jxx*QpJ1 + Kxx*QpK1 ux = Ixx*QpI2 + Jxx*QpJ2 + Kxx*QpK2 vx = Ixx*QpI3 + Jxx*QpJ3 + Kxx*QpK3 wx = Ixx*QpI4 + Jxx*QpJ4 + Kxx*QpK4 px = Ixx*QpI5 + Jxx*QpJ5 + Kxx*QpK5 Iyy = Iy(i,j,k,nbl) Jyy = Jy(i,j,k,nbl) Kyy = Ky(i,j,k,nbl) ry = Iyy*QpI1 + Jyy*QpJ1 + Kyy*QpK1 uy = Iyy*QpI2 + Jyy*QpJ2 + Kyy*QpK2 vy = Iyy*QpI3 + Jyy*QpJ3 + Kyy*QpK3 wy = Iyy*QpI4 + Jyy*QpJ4 + Kyy*QpK4 py = Iyy*QpI5 + Jyy*QpJ5 + Kyy*QpK5 Izz = Iz(i,j,k,nbl) Jzz = Jz(i,j,k,nbl) Kzz = Kz(i,j,k,nbl) ry = Izz*QpI1 + Jzz*QpJ1 + Kzz*QpK1 uy = Izz*QpI2 + Jzz*QpJ2 + Kzz*QpK2 vy = Izz*QpI3 + Jzz*QpJ3 + Kzz*QpK3 wy = Izz*QpI4 + Jzz*QpJ4 + Kzz*QpK4 py = Izz*QpI5 + Jzz*QpJ5 + Kzz*QpK5 ! =============================== !Temperature and viscosity at iph face (if velocity non-dimensionalized with 'U_inf-freestream Vel') T = Gamma*Mach**2*Qp(i,j,k,nbl,5)/Qp(i,j,k,nbl,1) mu_L = (T**1.5d0)*(1.4045988d0)/(T + 0.4045988d0)/Re !Temperature and viscosity at iph face (if velocity non-dimensionalized with 'a- freestream sound speed') ! T = Gamma*Qp_iph(5)/Qp_iph(1) ! mu_L = (T**1.5d0)*(1.4045988d0)/(T + 0.4045988d0)*Mach/Re ! For dimensional solver ! T = Qp_iph(5)/Qp_iph(1)/Rgas ! mu_L = 1.d0/Re ! =============================== !since shear stress tensor is a symmetric matrix, we only need to define 6 terms Blk_Visc_Term = (-2.d0/3.d0)*(ux + vy + wz) tau_xx = mu_L*(2.d0*ux + Blk_Visc_Term) ; tau_xy = mu_L*(uy + vx) tau_yy = mu_L*(2.d0*vy + Blk_Visc_Term) ; tau_yz = mu_L*(vz + wy) tau_zz = mu_L*(2.d0*wz + Blk_Visc_Term) ; tau_zx = mu_L*(wx + uz) ! Fourier_factor = mu_L*funGamma ! Tx = px/Qp(i,j,k,nbl,1) - Qp(i,j,k,nbl,5)/(Qp(i,j,k,nbl,1))**2.d0*rx ! Ty = py/Qp(i,j,k,nbl,1) - Qp(i,j,k,nbl,5)/(Qp(i,j,k,nbl,1))**2.d0*ry ! Tz = pz/Qp(i,j,k,nbl,1) - Qp(i,j,k,nbl,5)/(Qp(i,j,k,nbl,1))**2.d0*rz bx = (tau_xx*Qp(i,j,k,nbl,2) + tau_xy*Qp(i,j,k,nbl,3) + tau_zx*Qp(i,j,k,nbl,4)) & + mu_L*funGamma*(px/Qp(i,j,k,nbl,1) - Qp(i,j,k,nbl,5)/(Qp(i,j,k,nbl,1))**2.d0*rx) by = (tau_xy*Qp(i,j,k,nbl,2) + tau_yy*Qp(i,j,k,nbl,3) + tau_yz*Qp(i,j,k,nbl,4)) & + mu_L*funGamma*(py/Qp(i,j,k,nbl,1) - Qp(i,j,k,nbl,5)/(Qp(i,j,k,nbl,1))**2.d0*ry) bz = (tau_zx*Qp(i,j,k,nbl,2) + tau_yz*Qp(i,j,k,nbl,3) + tau_zz*Qp(i,j,k,nbl,4)) & + mu_L*funGamma*(pz/Qp(i,j,k,nbl,1) - Qp(i,j,k,nbl,5)/(Qp(i,j,k,nbl,1))**2.d0*rz) Jac1 = Jac(i,j,k,nbl) Visc_Fflux(i,j,k,nbl,1) =(Ixx*tau_xx+Iyy*tau_xy + Izz*tau_zx)/Jac1 Visc_Fflux(i,j,k,nbl,2) =(Ixx*tau_xy+Iyy*tau_yy + Izz*tau_yz)/Jac1 Visc_Fflux(i,j,k,nbl,3) =(Ixx*tau_zx+Iyy*tau_yz + Izz*tau_zz)/Jac1*(1-f2D) Visc_Fflux(i,j,k,nbl,4) =(Ixx*bx +Iyy*by + Izz*bz )/Jac1 Visc_Gflux(i,j,k,nbl,1) =(Jxx*tau_xx+Jyy*tau_xy + Jzz*tau_zx)/Jac1 Visc_Gflux(i,j,k,nbl,2) =(Jxx*tau_xy+Jyy*tau_yy + Jzz*tau_yz)/Jac1 Visc_Gflux(i,j,k,nbl,3) =(Jxx*tau_zx+Jyy*tau_yz + Jzz*tau_zz)/Jac1*(1-f2D) Visc_Gflux(i,j,k,nbl,4) =(Jxx*bx +Jyy*by + Jzz*bz )/Jac1 Visc_Hflux(i,j,k,nbl,1) =(Kxx*tau_xx+ Kyy*tau_xy+ Kzz*tau_zx)/Jac1 Visc_Hflux(i,j,k,nbl,2) =(Kxx*tau_xy+ Kyy*tau_yy+ Kzz*tau_yz)/Jac1 Visc_Hflux(i,j,k,nbl,3) =(Kxx*tau_zx+ Kyy*tau_yz+ Kzz*tau_zz)/Jac1 Visc_Hflux(i,j,k,nbl,4) =(Kxx*bx + Kyy*by + Kzz*bz )/Jac1 ENDDO ENDDO ENDDO ENDDO call DERIVATIVE_I5D_COMPACT(Visc_Fflux,ViscFlux_der,4,1) !$acc loop seq DO nbl = 1,nblocks !$acc parallel loop gang collapse(3) default(present) DO n_cons = 2,nconserv DO k = 1,NK(nbl) DO j = 1,NJ(nbl) !$acc loop vector DO i = 1,NI(nbl) Residual_RHS(i,j,k,nbl,n_cons) = Residual_RHS(i,j,k,nbl,n_cons) + ViscFlux_der(i,j,k,nbl,n_cons-1) ENDDO ENDDO ENDDO ENDDO ENDDO call DERIVATIVE_J5D_COMPACT(Visc_Gflux,ViscFlux_der,4,1) !$acc loop seq DO nbl = 1,nblocks !$acc parallel loop gang collapse(3) default(present) DO n_cons = 2,nconserv DO k = 1,NK(nbl) DO j = 1,NJ(nbl) !$acc loop vector DO i = 1,NI(nbl) Residual_RHS(i,j,k,nbl,n_cons) = Residual_RHS(i,j,k,nbl,n_cons) + ViscFlux_der(i,j,k,nbl,n_cons-1) ENDDO ENDDO ENDDO ENDDO ENDDO IF (f2D .ne. 1) THEN call DERIVATIVE_K5D_COMPACT(Visc_Hflux,ViscFlux_der,4,1) !$acc loop seq DO nbl = 1,nblocks !$acc parallel loop gang collapse(3) default(present) DO n_cons = 2,nconserv DO k = 1,NK(nbl) DO j = 1,NJ(nbl) !$acc loop vector DO i = 1,NI(nbl) Residual_RHS(i,j,k,nbl,n_cons) = Residual_RHS(i,j,k,nbl,n_cons) + ViscFlux_der(i,j,k,nbl,n_cons-1) ENDDO ENDDO ENDDO ENDDO ENDDO ENDIF !$acc end data END !================ ! SET PRIMITIVES !================ SUBROUTINE SET_PRIMITIVES() use declare_variables implicit none !$acc parallel loop gang collapse(3) default(present) DO nbl = 1,nblocks DO k = 1, NKmax DO j = 1, NJmax !$acc loop vector DO i = 1, NImax if (k.le.NK(nbl).and.j.le.NJ(nbl).and.i.le.NI(nbl)) then Qp(i,j,k,nbl,1) = Qc(i,j,k,nbl,1) !rho Qp(i,j,k,nbl,2) = Qc(i,j,k,nbl,2)/Qc(i,j,k,nbl,1) !u Qp(i,j,k,nbl,3) = Qc(i,j,k,nbl,3)/Qc(i,j,k,nbl,1) !v Qp(i,j,k,nbl,4) = Qc(i,j,k,nbl,4)/Qc(i,j,k,nbl,1) !w Qp(i,j,k,nbl,5) = (Qc(i,j,k,nbl,5) - 0.5d0*Qc(i,j,k,nbl,1)* & ((Qc(i,j,k,nbl,2)/Qc(i,j,k,nbl,1))**2 + (Qc(i,j,k,nbl,3)/Qc(i,j,k,nbl,1))**2 & + (Qc(i,j,k,nbl,4)/Qc(i,j,k,nbl,1))**2))*(gamma-1.d0) !P !(Qp(i,j,k,nbl,2)**2 + Qp(i,j,k,nbl,3)**2 + Qp(i,j,k,nbl,4)**2))*(gamma-1.d0) !P ! Clipping shit IF (Qp(i,j,k,nbl,1) < 0.d0) Qp(i,j,k,nbl,1) = 0.01d0 IF (Qp(i,j,k,nbl,5) < 0.d0) Qp(i,j,k,nbl,5) = 0.01d0 endif ENDDO ENDDO ENDDO ENDDO END SUBROUTINE SUBROUTINE POSITIVITY_PRESERVING() use declare_variables implicit none ! also make a do loop for prims here (its possible) !$acc loop seq DO nbl = 1,nblocks !$acc parallel loop gang collapse(3) default(present) DO n_prim = 1,nprims DO k = 1,NK(nbl) DO j = 1,NJ(nbl) !$acc loop vector DO i = 1,NI(nbl) IF (Qp_iphL(i,j,k,nbl,1) .lt. 0.d0 .or. Qp_iphL(i,j,k,nbl,5) .lt. 0.d0 .or. & Qp_iphR(i,j,k,nbl,1) .lt. 0.d0 .or. Qp_iphR(i,j,k,nbl,5) .lt. 0.d0) THEN Qp_iphL(i,j,k,nbl,n_prim) = Qp(i,j,k,nbl,n_prim) Qp_iphR(i,j,k,nbl,n_prim) = Qp(i+1,j,k,nbl,n_prim) ENDIF IF (Qp_jphL(i,j,k,nbl,1) .lt. 0.d0 .or. Qp_jphL(i,j,k,nbl,5) .lt. 0.d0 .or. & Qp_jphR(i,j,k,nbl,1) .lt. 0.d0 .or. Qp_jphR(i,j,k,nbl,5) .lt. 0.d0) THEN Qp_jphL(i,j,k,nbl,n_prim) = Qp(i,j,k,nbl,n_prim) Qp_jphR(i,j,k,nbl,n_prim) = Qp(i,j+1,k,nbl,n_prim) ENDIF IF (Qp_kphL(i,j,k,nbl,1) .lt. 0.d0 .or. Qp_kphL(i,j,k,nbl,5) .lt. 0.d0 .or. & Qp_kphR(i,j,k,nbl,1) .lt. 0.d0 .or. Qp_kphR(i,j,k,nbl,5) .lt. 0.d0) THEN Qp_kphL(i,j,k,nbl,n_prim) = Qp(i,j,k,nbl,n_prim) Qp_kphR(i,j,k,nbl,n_prim) = Qp(i,j,k+1,nbl,n_prim) ENDIF ENDDO ENDDO ENDDO ENDDO ENDDO END