diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 1859857b8..5a70f957f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -49,6 +49,7 @@ add_executable(xcompact3d tools.f90 transeq.f90 turbine.f90 + utilities.f90 variables.f90 visu.f90 xcompact3d.f90) diff --git a/src/Case-ABL.f90 b/src/Case-ABL.f90 index 906111236..fc8c9d777 100644 --- a/src/Case-ABL.f90 +++ b/src/Case-ABL.f90 @@ -383,7 +383,6 @@ subroutine wall_sgs_slip(ux,uy,uz,phi,nut1,wallfluxx,wallfluxy,wallfluxz) use var, only: di1, di3 use var, only: sxy1, syz1, heatflux, ta2, tb2, ta3, tb3 use ibm_param, only : ubcx, ubcz - use dbg_schemes, only: log_prec, tanh_prec, sqrt_prec, abs_prec, atan_prec implicit none real(mytype),dimension(xsize(1),xsize(2),xsize(3)),intent(in) :: ux,uy,uz, nut1 @@ -512,19 +511,19 @@ subroutine wall_sgs_slip(ux,uy,uz,phi,nut1,wallfluxx,wallfluxy,wallfluxz) if (itherm==0) then Q_HAve = TempFlux else if (itherm==1) then - Q_HAve =-k_roughness**two*S_HAve*(Phi_HAve-(T_wall+TempRate*t))/((log_prec(delta/z_zero)-PsiM_HAve)*(log_prec(delta/z_zero)-PsiH_HAve)) + Q_HAve =-k_roughness**two*S_HAve*(Phi_HAve-(T_wall+TempRate*t))/((log(delta/z_zero)-PsiM_HAve)*(log(delta/z_zero)-PsiH_HAve)) endif - L_HAve=-(k_roughness*S_HAve/(log_prec(delta/z_zero)-PsiM_HAve))**three*Phi_HAve/(k_roughness*gravv*Q_HAve) + L_HAve=-(k_roughness*S_HAve/(log(delta/z_zero)-PsiM_HAve))**three*Phi_HAve/(k_roughness*gravv*Q_HAve) if (istrat==0) then PsiM_HAve=-4.8_mytype*delta/L_HAve PsiH_HAve=-7.8_mytype*delta/L_HAve else if (istrat==1) then zeta_HAve=(one-sixteen*delta/L_HAve)**zptwofive - PsiM_HAve=two*log_prec(half*(one+zeta_HAve))+log_prec(zpfive*(one+zeta_HAve**two))-two*atan_prec(zeta_HAve)+pi/two - PsiH_HAve=two*log_prec(half*(one+zeta_HAve**two)) + PsiM_HAve=two*log(half*(one+zeta_HAve))+log(zpfive*(one+zeta_HAve**two))-two*atan(zeta_HAve)+pi/two + PsiH_HAve=two*log(half*(one+zeta_HAve**two)) endif ii = ii + 1 - OL_diff = abs_prec((L_HAve - Lold)/Lold) + OL_diff = abs((L_HAve - Lold)/Lold) Lold = L_HAve if (ii==50) exit enddo @@ -548,30 +547,30 @@ subroutine wall_sgs_slip(ux,uy,uz,phi,nut1,wallfluxx,wallfluxy,wallfluxz) do i=1,xsize(1) ! Horizontally-averaged formulation if(iwallmodel==1) then - tauwallxy(i,k)=-(k_roughness/(log_prec(delta/z_zero)-PsiM_HAve))**two*ux_HAve*S_HAve - tauwallzy(i,k)=-(k_roughness/(log_prec(delta/z_zero)-PsiM_HAve))**two*uz_HAve*S_HAve + tauwallxy(i,k)=-(k_roughness/(log(delta/z_zero)-PsiM_HAve))**two*ux_HAve*S_HAve + tauwallzy(i,k)=-(k_roughness/(log(delta/z_zero)-PsiM_HAve))**two*uz_HAve*S_HAve ! Local formulation else ux12=half*(uxf1(i,1,k)+uxf1(i,2,k)) uz12=half*(uzf1(i,1,k)+uzf1(i,2,k)) - S12=sqrt_prec(ux12**2.+uz12**2.) + S12=sqrt(ux12**2.+uz12**2.) if (iscalar==1) then Phi12= half*(phif1(i,1,k)+ phif1(i,2,k)) + Tstat12 do ii=1,10 - if (itherm==1) heatflux(i,k)=-k_roughness**two*S12*(Phi12-(T_wall+TempRate*t))/((log_prec(delta/z_zero)-PsiM(i,k))*(log_prec(delta/z_zero)-PsiH(i,k))) - Obukhov(i,k)=-(k_roughness*S12/(log_prec(delta/z_zero)-PsiM(i,k)))**three*Phi12/(k_roughness*gravv*heatflux(i,k)) + if (itherm==1) heatflux(i,k)=-k_roughness**two*S12*(Phi12-(T_wall+TempRate*t))/((log(delta/z_zero)-PsiM(i,k))*(log(delta/z_zero)-PsiH(i,k))) + Obukhov(i,k)=-(k_roughness*S12/(log(delta/z_zero)-PsiM(i,k)))**three*Phi12/(k_roughness*gravv*heatflux(i,k)) if (istrat==0) then PsiM(i,k)=-4.8_mytype*delta/Obukhov(i,k) PsiH(i,k)=-7.8_mytype*delta/Obukhov(i,k) else if (istrat==1) then zeta(i,k)=(one-sixteen*delta/Obukhov(i,k))**zptwofive - PsiM(i,k)=two*log_prec(half*(one+zeta(i,k)))+log_prec(zpfive*(one+zeta(i,k)**2.))-two*atan_prec(zeta(i,k))+pi/two - PsiH(i,k)=two*log_prec(half*(one+zeta(i,k)**two)) + PsiM(i,k)=two*log(half*(one+zeta(i,k)))+log(zpfive*(one+zeta(i,k)**2.))-two*atan(zeta(i,k))+pi/two + PsiH(i,k)=two*log(half*(one+zeta(i,k)**two)) endif enddo endif - tauwallxy(i,k)=-(k_roughness/(log_prec(delta/z_zero)-PsiM(i,k)))**two*ux12*S12 - tauwallzy(i,k)=-(k_roughness/(log_prec(delta/z_zero)-PsiM(i,k)))**two*uz12*S12 + tauwallxy(i,k)=-(k_roughness/(log(delta/z_zero)-PsiM(i,k)))**two*ux12*S12 + tauwallzy(i,k)=-(k_roughness/(log(delta/z_zero)-PsiM(i,k)))**two*uz12*S12 endif ! Apply second-order upwind scheme for the near wall ! Below should change for non-uniform grids, same for wall_sgs_slip_scalar @@ -617,7 +616,7 @@ subroutine wall_sgs_slip(ux,uy,uz,phi,nut1,wallfluxx,wallfluxy,wallfluxz) endif ! Compute friction velocity u_shear and boundary layer height - u_shear=k_roughness*S_HAve/(log_prec(delta/z_zero)-PsiM_HAve) + u_shear=k_roughness*S_HAve/(log(delta/z_zero)-PsiM_HAve) if (iheight==1) call boundary_height(ux,uy,uz,dBL) if (iscalar==1) zL=dBL/L_HAve @@ -692,7 +691,6 @@ subroutine wall_sgs_noslip(ux1,uy1,uz1,nut1,wallsgsx1,wallsgsy1,wallsgsz1) use variables use var, only: di1, di2, di3 use var, only: sxy1, syz1, tb1, ta2, tb2 - use dbg_schemes, only: log_prec, sqrt_prec use ibm_param, only : ubcx, ubcy, ubcz implicit none @@ -738,9 +736,9 @@ subroutine wall_sgs_noslip(ux1,uy1,uz1,nut1,wallsgsx1,wallsgsy1,wallsgsz1) y_sampling=delta-real(j0,mytype)*dy ux_delta=(1-y_sampling/dy)*ta2(i,j0+1,k)+(y_sampling/dy)*ta2(i,j0+2,k) uz_delta=(1-y_sampling/dy)*tb2(i,j0+1,k)+(y_sampling/dy)*tb2(i,j0+2,k) - S_delta=sqrt_prec(ux_delta**2.+uz_delta**2.) - tauwallxy2(i,k)=-(k_roughness/(log_prec(delta/z_zero)))**two*ux_delta*S_delta - tauwallzy2(i,k)=-(k_roughness/(log_prec(delta/z_zero)))**two*uz_delta*S_delta + S_delta=sqrt(ux_delta**2.+uz_delta**2.) + tauwallxy2(i,k)=-(k_roughness/(log(delta/z_zero)))**two*ux_delta*S_delta + tauwallzy2(i,k)=-(k_roughness/(log(delta/z_zero)))**two*uz_delta*S_delta txy2(i,2,k) = tauwallxy2(i,k) tyz2(i,2,k) = tauwallzy2(i,k) enddo @@ -809,8 +807,8 @@ subroutine wall_sgs_noslip(ux1,uy1,uz1,nut1,wallsgsx1,wallsgsy1,wallsgsz1) call MPI_ALLREDUCE(uz_HAve_local,uz_HAve,1,real_type,MPI_SUM,MPI_COMM_WORLD,code) ux_HAve=ux_HAve/(nxc*nzc) uz_HAve=uz_HAve/(nxc*nzc) - S_HAve=sqrt_prec(ux_HAve**2.+uz_HAve**2.) - u_shear=k_roughness*S_HAve/log_prec(5*dy/z_zero) + S_HAve=sqrt(ux_HAve**2.+uz_HAve**2.) + u_shear=k_roughness*S_HAve/log(5*dy/z_zero) if (mod(itime,ilist)==0.and.nrank==0) then ! Write u_shear in file @@ -835,7 +833,6 @@ subroutine forceabl (ux) ! Routine to force constant flow rate use param use var use MPI - use dbg_schemes, only: log_prec implicit none @@ -995,7 +992,6 @@ subroutine damping_zone (dux1,duy1,duz1,ux1,uy1,uz1) ! Damping zone for ABL use param use var, only: yp - use dbg_schemes, only: sin_prec, cos_prec, log_prec implicit none @@ -1022,7 +1018,7 @@ subroutine damping_zone (dux1,duy1,duz1,ux1,uy1,uz1) ! Damping zone for ABL ! Damping for non-neutral ABL if (ibuoyancy==1) then if (y>=damp_lo) then - lambda=sin_prec(half*pi*(y-damp_lo)/(yly-damp_lo))**two + lambda=sin(half*pi*(y-damp_lo)/(yly-damp_lo))**two else lambda=zero endif @@ -1034,13 +1030,13 @@ subroutine damping_zone (dux1,duy1,duz1,ux1,uy1,uz1) ! Damping zone for ABL if (y>=(dBL+half*dheight)) then lambda=one elseif (y>=(dBL-half*dheight).and.y<(dBL+half*dheight)) then - lambda=half*(one-cos_prec(pi*(y-(dBL-half*dheight))/dheight)) + lambda=half*(one-cos(pi*(y-(dBL-half*dheight))/dheight)) else lambda=zero endif xloc=real(i-1,mytype)*dx if (iconcprec.eq.1.and.xloc.ge.pdl) lambda=0. - dux1(i,j,k,1)=dux1(i,j,k,1)-coeff*lambda*(ux1(i,j,k)-ustar/k_roughness*log_prec(dBL/z_zero)) + dux1(i,j,k,1)=dux1(i,j,k,1)-coeff*lambda*(ux1(i,j,k)-ustar/k_roughness*log(dBL/z_zero)) duy1(i,j,k,1)=duy1(i,j,k,1)-coeff*lambda*(uy1(i,j,k)-UG(2)) duz1(i,j,k,1)=duz1(i,j,k,1)-coeff*lambda*(uz1(i,j,k)-UG(3)) endif @@ -1058,7 +1054,6 @@ subroutine damping_zone_scalar (dphi1,phi1) use param use var, only: yp - use dbg_schemes, only: sin_prec implicit none @@ -1077,7 +1072,7 @@ subroutine damping_zone_scalar (dphi1,phi1) if (istret==0) y=real(j+xstart(2)-1-1,mytype)*dy if (istret/=0) y=yp(j+xstart(2)-1) if (y>=damp_lo) then - lambda=sin_prec(half*pi*(y-damp_lo)/(yly-damp_lo))**two + lambda=sin(half*pi*(y-damp_lo)/(yly-damp_lo))**two else lambda=zero endif diff --git a/src/Case-Cylinder-wake.f90 b/src/Case-Cylinder-wake.f90 index 1e01e4b1b..c4b8b8197 100644 --- a/src/Case-Cylinder-wake.f90 +++ b/src/Case-Cylinder-wake.f90 @@ -26,7 +26,6 @@ subroutine geomcomplex_cyl(epsi,nxi,nxf,ny,nyi,nyf,nzi,nzf,dx,yp,remp) use param, only : one, two, ten use ibm_param - use dbg_schemes, only: sqrt_prec implicit none @@ -69,7 +68,7 @@ subroutine geomcomplex_cyl(epsi,nxi,nxf,ny,nyi,nyf,nzi,nzf,dx,yp,remp) ym=yp(j) do i=nxi,nxf xm=real(i-1,mytype)*dx - r=sqrt_prec((xm-cexx)**two+(ym-ceyy)**two) + r=sqrt((xm-cexx)**two+(ym-ceyy)**two) if (r-ra.gt.zeromach) then cycle endif @@ -209,7 +208,6 @@ subroutine init_cyl (ux1,uy1,uz1,phi1) USE variables USE param USE MPI - use dbg_schemes, only: exp_prec implicit none @@ -252,7 +250,7 @@ subroutine init_cyl (ux1,uy1,uz1,phi1) do j=1,xsize(2) if (istret.eq.0) y=(j+xstart(2)-1-1)*dy-yly/2. if (istret.ne.0) y=yp(j+xstart(2)-1)-yly/2. - um=exp_prec(-zptwo*y*y) + um=exp(-zptwo*y*y) do i=1,xsize(1) ux1(i,j,k)=um*ux1(i,j,k) uy1(i,j,k)=um*uy1(i,j,k) @@ -290,7 +288,6 @@ subroutine postprocess_cyl(ux1,uy1,uz1,ep1) USE var, only : ta1,tb1,tc1,td1,te1,tf1,tg1,th1,ti1,di1 USE var, only : ta2,tb2,tc2,td2,te2,tf2,di2,ta3,tb3,tc3,td3,te3,tf3,di3 USE ibm_param - use dbg_schemes, only: sqrt_prec real(mytype),intent(in),dimension(xsize(1),xsize(2),xsize(3)) :: ux1, uy1, uz1, ep1 diff --git a/src/Case-Mixing-layer.f90 b/src/Case-Mixing-layer.f90 index 90066fcb0..810167be5 100644 --- a/src/Case-Mixing-layer.f90 +++ b/src/Case-Mixing-layer.f90 @@ -22,7 +22,6 @@ subroutine init_mixlayer (rho1,ux1,uy1,uz1) use param, only : u1, u2, dens1, dens2 use param, only : half, one, two, four, eight, sixteen use param, only : ntime, nrhotime - use dbg_schemes, only: sin_prec, cos_prec, exp_prec, tanh_prec, sqrt_prec use MPI implicit none @@ -57,8 +56,8 @@ subroutine init_mixlayer (rho1,ux1,uy1,uz1) rhomax = max(dens1, dens2) T1 = pressure0 / dens1 T2 = pressure0 / dens2 - u1 = sqrt_prec(dens2 / dens1) / (sqrt_prec(dens2 / dens1) + one) - u2 = -sqrt_prec(dens1 / dens2) / (one + sqrt_prec(dens1 / dens2)) + u1 = sqrt(dens2 / dens1) / (sqrt(dens2 / dens1) + one) + u2 = -sqrt(dens1 / dens2) / (one + sqrt(dens1 / dens2)) M = 0.2_mytype rspech = 1.4_mytype heatcap = (one / (T2 * (rspech - one))) * ((u1 - u2) / M)**2 @@ -71,7 +70,7 @@ subroutine init_mixlayer (rho1,ux1,uy1,uz1) !! Set mean field ux1(i, j, k) = ux1(i, j, k) + half * (u1 + u2) & - + half * (u1 - u2) * tanh_prec(two * y) + + half * (u1 - u2) * tanh(two * y) uy1(i, j, k) = zero uz1(i, j, k) = zero @@ -84,14 +83,14 @@ subroutine init_mixlayer (rho1,ux1,uy1,uz1) rho1(i, j, k, 1) = MIN(rho1(i, j, k, 1), rhomax) ! Calculate disturbance field (as given in Fortune2004) - disturb_decay = 0.025_mytype * (u1 - u2) * exp_prec(-0.05_mytype * (y**2)) - u_disturb = disturb_decay * (sin_prec(eight * PI * x / xlx) & - + sin_prec(four * PI * x / xlx) / eight & - + sin_prec(two * PI * x / xlx) / sixteen) + disturb_decay = 0.025_mytype * (u1 - u2) * exp(-0.05_mytype * (y**2)) + u_disturb = disturb_decay * (sin(eight * PI * x / xlx) & + + sin(four * PI * x / xlx) / eight & + + sin(two * PI * x / xlx) / sixteen) u_disturb = (0.05_mytype * y * xlx / PI) * u_disturb - v_disturb = disturb_decay * (cos_prec(eight * PI * x / xlx) & - + cos_prec(four * PI * x / xlx) / eight & - + cos_prec(two * PI * x / xlx) / sixteen) + v_disturb = disturb_decay * (cos(eight * PI * x / xlx) & + + cos(four * PI * x / xlx) / eight & + + cos(two * PI * x / xlx) / sixteen) ux1(i, j, k) = ux1(i, j, k) + u_disturb uy1(i, j, k) = uy1(i, j, k) + v_disturb diff --git a/src/Case-TBL.f90 b/src/Case-TBL.f90 index c74a90f8e..465557679 100644 --- a/src/Case-TBL.f90 +++ b/src/Case-TBL.f90 @@ -175,7 +175,6 @@ subroutine blasius() use decomp_2d_io use MPI use param, only : zero, zptwo, zpeight, one, nine - use dbg_schemes, only: exp_prec, sqrt_prec implicit none @@ -213,7 +212,7 @@ subroutine blasius() -17.018112827391400_mytype*eta_bl**4 +0.819582894357566_mytype*eta_bl**3 & -0.0601348202321954_mytype*eta_bl**2 +2.989739912704050_mytype*eta_bl**1 - f_bl=f_bl+(1-exp_prec(-delta_eta/delta_int))*eps_eta + f_bl=f_bl+(1-exp(-delta_eta/delta_int))*eps_eta @@ -232,14 +231,14 @@ subroutine blasius() +130.109496961069_mytype*eta_bl**4 -7.64507811014497000_mytype*eta_bl**3 & +6.94303207046209_mytype*eta_bl**2 -0.00209716712558639_mytype*eta_bl**1 ! & - g_bl=g_bl+(1-exp_prec(-delta_eta/delta_int))*eps_eta + g_bl=g_bl+(1-exp(-delta_eta/delta_int))*eps_eta x_bl=one/(4.91_mytype**2*xnu) bxx1(j,k)=f_bl/1.0002014996204402_mytype/1.0000000359138641_mytype !To assure 1.0 in infinity - bxy1(j,k)=g_bl*sqrt_prec(xnu/x_bl)/1.000546554_mytype + bxy1(j,k)=g_bl*sqrt(xnu/x_bl)/1.000546554_mytype bxz1(j,k)=zero enddo enddo @@ -269,7 +268,7 @@ subroutine blasius() -0.0601348202321954_mytype*eta_bl**2 +2.989739912704050_mytype*eta_bl**1 - f_bl_inf=f_bl_inf+(1-exp_prec(-delta_eta/delta_int))*eps_eta + f_bl_inf=f_bl_inf+(1-exp(-delta_eta/delta_int))*eps_eta f_bl_inf=f_bl_inf/1.0002014996204402_mytype/1.0000000359138641_mytype !To assure 1.0 in infinity #ifdef DEBG @@ -292,7 +291,7 @@ subroutine blasius() +6.94303207046209_mytype*eta_bl**2 -0.00209716712558639_mytype*eta_bl**1 - g_bl_inf=g_bl_inf+(1-exp_prec(-delta_eta/delta_int))*eps_eta + g_bl_inf=g_bl_inf+(1-exp(-delta_eta/delta_int))*eps_eta g_bl_inf=g_bl_inf/1.000546554_mytype #ifdef DEBG if (nrank == 0) write(*,*)'g_bl_inf ', g_bl_inf diff --git a/src/Case-TGV.f90 b/src/Case-TGV.f90 index d42997f8c..ce552115f 100644 --- a/src/Case-TGV.f90 +++ b/src/Case-TGV.f90 @@ -28,7 +28,6 @@ subroutine init_tgv (ux1,uy1,uz1,ep1,phi1) use variables use param use MPI - use dbg_schemes, only: sin_prec, cos_prec implicit none @@ -67,10 +66,10 @@ subroutine init_tgv (ux1,uy1,uz1,ep1,phi1) do i=1,xsize(1) x=real(i-1,mytype)*dx - ux1(i,j,k)=+sin_prec(x)*cos_prec(y)*cos_prec(z) - uy1(i,j,k)=-cos_prec(x)*sin_prec(y)*cos_prec(z) + ux1(i,j,k)=+sin(x)*cos(y)*cos(z) + uy1(i,j,k)=-cos(x)*sin(y)*cos(z) if (iscalar == 1) then - phi1(i,j,k,1:numscalar)=sin_prec(x)*sin_prec(y)*cos_prec(z) + phi1(i,j,k,1:numscalar)=sin(x)*sin(y)*cos(z) endif uz1(i,j,k)=zero enddo @@ -575,7 +574,6 @@ subroutine error_tgv2D(ux1,uy1,phi1) use MPI use param, only : one, two, xnu, ifirst, itime use variables, only : numscalar, sc - use dbg_schemes, only: sin_prec, cos_prec implicit none @@ -610,8 +608,8 @@ subroutine error_tgv2D(ux1,uy1,phi1) do i = 1,xsize(1) x = real(i+xstart(1)-1-1,mytype)*dx ! Initial solution - solx0 = sin_prec(x) * cos_prec(y) - soly0 = - cos_prec(x) * sin_prec(y) + solx0 = sin(x) * cos(y) + soly0 = - cos(x) * sin(y) ! Analytical solution solxt = solx0 * xdamping(1) solyt = soly0 * ydamping(1) @@ -639,7 +637,7 @@ subroutine error_tgv2D(ux1,uy1,phi1) do i = 1,xsize(1) x = real(i+xstart(1)-1-1,mytype)*dx ! Initial solution - sols0 = sin_prec(x) * sin_prec(y) + sols0 = sin(x) * sin(y) ! Analytical solution solst = sols0 * sdamping(1,l) ! Time discrete solution @@ -736,7 +734,6 @@ subroutine compute_tgv2D_errors(xdamping, ydamping, sdamping) use param, only : one, two, xnu, ifirst, itime, itimescheme, iimplicit use variables, only : numscalar, sc - use dbg_schemes, only: exp_prec implicit none @@ -755,12 +752,12 @@ subroutine compute_tgv2D_errors(xdamping, ydamping, sdamping) sdamping(:,:) = one ! Compute analytical damping - coef(1) = exp_prec(-two*dt*xnu) + coef(1) = exp(-two*dt*xnu) do it = ifirst, itime xdamping(1) = xdamping(1) * coef(1) ydamping(1) = ydamping(1) * coef(1) do l = 1, numscalar - sdamping(1,l) = sdamping(1,l) * exp_prec(-two*dt*xnu/sc(l)) + sdamping(1,l) = sdamping(1,l) * exp(-two*dt*xnu/sc(l)) enddo enddo @@ -838,7 +835,6 @@ subroutine compute_k2(kin, k2out) use param use derivx, only : alsaix, asix, bsix, csix, dsix - use dbg_schemes, only: cos_prec implicit none @@ -851,11 +847,11 @@ subroutine compute_k2(kin, k2out) endif endif - k2out = asix * two * (one - cos_prec(kin*dx)) & - + four * bsix * half * (one - cos_prec(two*kin*dx)) & - + nine * csix * (two / nine) * (one - cos_prec(three*kin*dx)) & - + sixteen * dsix * (one / eight) * (one - cos_prec(four*kin*dx)) - k2out = k2out / (one + two * alsaix * cos_prec(kin*dx)) + k2out = asix * two * (one - cos(kin*dx)) & + + four * bsix * half * (one - cos(two*kin*dx)) & + + nine * csix * (two / nine) * (one - cos(three*kin*dx)) & + + sixteen * dsix * (one / eight) * (one - cos(four*kin*dx)) + k2out = k2out / (one + two * alsaix * cos(kin*dx)) end subroutine compute_k2 diff --git a/src/Case-Uniform.f90 b/src/Case-Uniform.f90 index 90b5b6d79..531ed3912 100644 --- a/src/Case-Uniform.f90 +++ b/src/Case-Uniform.f90 @@ -237,7 +237,6 @@ subroutine postprocess_uniform(ux1,uy1,uz1,ep1) USE var, only : ta1,tb1,tc1,td1,te1,tf1,tg1,th1,ti1,di1 USE var, only : ta2,tb2,tc2,td2,te2,tf2,di2,ta3,tb3,tc3,td3,te3,tf3,di3 USE ibm_param - use dbg_schemes, only: sqrt_prec real(mytype),intent(in),dimension(xsize(1),xsize(2),xsize(3)) :: ux1, uy1, uz1, ep1 diff --git a/src/acl_elem.f90 b/src/acl_elem.f90 index 70faab818..65110b31f 100644 --- a/src/acl_elem.f90 +++ b/src/acl_elem.f90 @@ -149,7 +149,6 @@ subroutine set_actuatorline_geometry(actuatorline) !******************************************************************************* use param, only: zero, one - use dbg_schemes, only: cos_prec, sin_prec implicit none type(ActuatorLineType), intent(inout) :: actuatorline @@ -181,8 +180,8 @@ subroutine set_actuatorline_geometry(actuatorline) actuatorline%pitch(istation)=pitch(istation)*conrad endif - actuatorline%tx(istation)= cos_prec(actuatorline%pitch(istation)) - actuatorline%ty(istation)=-sin_prec(actuatorline%pitch(istation)) + actuatorline%tx(istation)= cos(actuatorline%pitch(istation)) + actuatorline%ty(istation)=-sin(actuatorline%pitch(istation)) actuatorline%tz(istation)= zero actuatorline%C(istation)=ctoR(istation)*length actuatorline%thick(istation)=thick(istation) @@ -233,7 +232,6 @@ subroutine Compute_ActuatorLine_Forces(act_line,rho_air,visc,dt,time) ! !******************************************************************************* - use dbg_schemes, only: sqrt_prec, sin_prec, cos_prec use param, only: zpone, zpfive, zero, one, two, four, eight implicit none @@ -280,7 +278,7 @@ subroutine Compute_ActuatorLine_Forces(act_line,rho_air,visc,dt,time) urdn=nxe*(u-ub)+nye*(v-vb)+nze*(w-wb) ! Normal urdc=txe*(u-ub)+tye*(v-vb)+tze*(w-wb) ! Tangential - ur=sqrt_prec(urdn**2.0+urdc**2.0) + ur=sqrt(urdn**2.0+urdc**2.0) act_line%EUr(ielem)=ur ! This is the dynamic angle of attack @@ -329,8 +327,8 @@ subroutine Compute_ActuatorLine_Forces(act_line,rho_air,visc,dt,time) CMAM=-CNAM/four-urdn*urdc/(eight*ur**2) CT=CT+CTAM CN=CN+CNAM - CL=CN*cos_prec(alpha)-CT*sin_prec(alpha) - CD=CN*sin_prec(alpha)+CT*cos_prec(alpha) + CL=CN*cos(alpha)-CT*sin(alpha) + CD=CN*sin(alpha)+CT*cos(alpha) CM25=CM25+CMAM endif @@ -338,7 +336,7 @@ subroutine Compute_ActuatorLine_Forces(act_line,rho_air,visc,dt,time) if (act_line%do_random_walk_forcing) then Strouhal=0.17_mytype freq=Strouhal*ur/max(ElemChord,0.0001_mytype) - CL=CL*(one+zpone*sin_prec(two*pi*freq*time)+0.05_mytype*(-one+two*rand(ielem))) + CL=CL*(one+zpone*sin(two*pi*freq*time)+0.05_mytype*(-one+two*rand(ielem))) CD=CD*(one+0.05_mytype*(-one+two*rand(ielem))) endif @@ -349,8 +347,8 @@ subroutine Compute_ActuatorLine_Forces(act_line,rho_air,visc,dt,time) endif ! Tangential and normal coeffs - CN=CL*cos_prec(alpha)+CD*sin_prec(alpha) - CT=-CL*sin_prec(alpha)+CD*cos_prec(alpha) + CN=CL*cos(alpha)+CD*sin(alpha) + CT=-CL*sin(alpha)+CD*cos(alpha) ! Apply coeffs to calculate tangential and normal Forces FN=zpfive*rho_air*CN*ElemArea*ur**2.0 @@ -402,7 +400,6 @@ subroutine compute_Tower_Forces(tower,rho_air,visc,time,CL,CD,Str) !******************************************************************************* use param, only: zero, zpfive, one, two - use dbg_schemes, only: sin_prec, cos_prec implicit none type(ActuatorLineType), intent(inout) :: tower @@ -449,13 +446,13 @@ subroutine compute_Tower_Forces(tower,rho_air,visc,time,CL,CD,Str) ! Apply random walk on the lift and drag forces freq=Str*ur/max(Diameter,0.0001_mytype) - tower%ECL(ielem)=CL*sin_prec(two*freq*pi*time) + tower%ECL(ielem)=CL*sin(two*freq*pi*time) tower%ECL(ielem)=tower%ECL(ielem)*(one+0.05_mytype*(-one+two*rand(ielem))) tower%ECD(ielem)=CD ! Tangential and normal coeffs - CN=tower%ECL(ielem)*cos_prec(alpha)+tower%ECD(ielem)*sin_prec(alpha) - CT=-tower%ECL(ielem)*sin_prec(alpha)+tower%ECD(ielem)*cos_prec(alpha) + CN=tower%ECL(ielem)*cos(alpha)+tower%ECD(ielem)*sin(alpha) + CT=-tower%ECL(ielem)*sin(alpha)+tower%ECD(ielem)*cos(alpha) ! Apply Coeffs to calculate tangential and normal Forces FN=zpfive*rho_air*CN*ElemArea*ur**2.0 @@ -485,8 +482,6 @@ subroutine pitch_actuator_line(act_line,pitch_angle) ! !******************************************************************************* - use dbg_schemes, only: sqrt_prec - implicit none type(ActuatorLineType), intent(inout) :: act_line real(mytype), intent(inout) :: pitch_angle ! Pitch in degrees @@ -497,7 +492,7 @@ subroutine pitch_actuator_line(act_line,pitch_angle) Nstation=act_line%Nelem+1 ! Compute the blade-wise unit vector - S=sqrt_prec((act_line%QCX(NStation)-act_line%COR(1))**2. + & + S=sqrt((act_line%QCX(NStation)-act_line%COR(1))**2. + & (act_line%QCY(NStation)-act_line%COR(2))**2. + & (act_line%QCZ(NStation)-act_line%COR(3))**2.) @@ -532,7 +527,6 @@ subroutine populate_blade_airfoils(NElem,NData,EAirfoil,AirfoilData,ETtoC) !******************************************************************************* use param, only: zero - use dbg_schemes, only: abs_prec implicit none integer, intent(in) :: NElem,NData @@ -559,7 +553,7 @@ subroutine populate_blade_airfoils(NElem,NData,EAirfoil,AirfoilData,ETtoC) do idata=1,NData Thicks(idata)=AirfoilData(idata)%tc - diffthicks(idata)=abs_prec(AirfoilData(idata)%tc-ETtoC(ielem)) + diffthicks(idata)=abs(AirfoilData(idata)%tc-ETtoC(ielem)) enddo imin=minloc(Thicks,1) imax=maxloc(Thicks,1) @@ -582,7 +576,6 @@ subroutine rotate_actuatorline(actuatorline,origin,rotN,theta) ! !******************************************************************************* - use dbg_schemes, only: sqrt_prec use param, only : zero implicit none @@ -619,7 +612,7 @@ subroutine rotate_actuatorline(actuatorline,origin,rotN,theta) tztmp=actuatorline%tz(ielem) call QuatRot(txtmp,tytmp,tztmp,theta,nrx,nry,nrz,zero,zero,zero,vrx,vry,vrz) - VMag=sqrt_prec(vrx**2+vry**2+vrz**2) + VMag=sqrt(vrx**2+vry**2+vrz**2) actuatorline%tx(ielem)=vrx/VMag actuatorline%ty(ielem)=vry/VMag actuatorline%tz(ielem)=vrz/VMag @@ -673,7 +666,6 @@ subroutine make_actuatorline_geometry(blade) !******************************************************************************* use param, only: two - use dbg_schemes, only: sqrt_prec implicit none type(ActuatorLineType), intent(inout) :: blade ! For simplity. In fact this is an actuator line. diff --git a/src/acl_model.f90 b/src/acl_model.f90 index 9568bf379..9c56354f5 100644 --- a/src/acl_model.f90 +++ b/src/acl_model.f90 @@ -471,7 +471,6 @@ subroutine actuator_line_model_update(current_time,dt) use decomp_2d_constants, only : mytype use param, only: zpfive, two, sixty, onehundredeighty - use dbg_schemes, only: sin_prec implicit none real(mytype), intent(inout) :: current_time, dt @@ -576,7 +575,7 @@ subroutine actuator_line_model_update(current_time,dt) ! Do harmonic pitch control for all elements (stations) of the actuator line Nstation=ActuatorLine(i)%NElem+1 do j=1,Nstation - ActuatorLine(i)%pitch(j)=actuatorline(i)%pitchAmp*sin_prec(actuatorline(i)%angular_pitch_freq*(ctime-ActuatorLine(i)%pitch_start_time)) + ActuatorLine(i)%pitch(j)=actuatorline(i)%pitchAmp*sin(actuatorline(i)%angular_pitch_freq*(ctime-ActuatorLine(i)%pitch_start_time)) enddo if (nrank==0.and.mod(itime,ilist)==0) then write(*,*) 'Harmonic pitch :' diff --git a/src/acl_source.f90 b/src/acl_source.f90 index 1cf18a85c..6219396b0 100644 --- a/src/acl_source.f90 +++ b/src/acl_source.f90 @@ -8,7 +8,6 @@ module actuator_line_source use decomp_2d_constants, only: real_type use variables, only: ilist use param, only: itime, zero, half, one - use dbg_schemes, only: sin_prec, sqrt_prec use actuator_line_model_utils use actuator_line_model @@ -320,7 +319,7 @@ subroutine Compute_Momentum_Source_Term_pointwise if (istret.ne.0) ymesh=yp(j) do i=xstart(1),xend(1) xmesh=real(i-1,mytype)*dx - dist = sqrt_prec((Sx(isource)-xmesh)**2.+(Sy(isource)-ymesh)**2.+(Sz(isource)-zmesh)**2.) + dist = sqrt((Sx(isource)-xmesh)**2.+(Sy(isource)-ymesh)**2.+(Sz(isource)-zmesh)**2.) if (distone) then + Ftip=two/pi*acos(exp(-g1*turbine%Nblades*half*(one/rroot-one)/sin(phi))) + if (abs(exp(-g1*turbine%Nblades*half*(one/rroot-one)/sin(phi)))>one) then if (nrank==0) write(*,*) 'Something went wrong with the tip correction model -- phi =', phi Ftip=one endif @@ -378,13 +377,13 @@ subroutine Compute_Turbine_EndEffects(turbine) if (turbine%EndEffectModel_is_Glauert) then g1=one else if (turbine%EndEffectModel_is_Shen) then - g1=exp_prec(-turbine%ShenCoeff_c1*(turbine%NBlades*turbine%TSR-turbine%ShenCoeff_c2))+zpone + g1=exp(-turbine%ShenCoeff_c1*(turbine%NBlades*turbine%TSR-turbine%ShenCoeff_c2))+zpone else write(*,*) 'Only Glauert and Shen et al. are available at the moment' stop endif - Froot=two/pi*acos_prec(exp_prec(-g1*turbine%Nblades*half*(one/rtip-one)/sin_prec(phi))) - if (abs_prec(exp_prec(-g1*turbine%Nblades*half*(one/rtip-one)/sin_prec(phi)))>one) then + Froot=two/pi*acos(exp(-g1*turbine%Nblades*half*(one/rtip-one)/sin(phi))) + if (abs(exp(-g1*turbine%Nblades*half*(one/rtip-one)/sin(phi)))>one) then if (nrank==0) write(*,*) 'Something went wrong with the root correction model -- phi =', phi Froot=one endif @@ -647,7 +646,7 @@ subroutine compute_rotor_upstream_velocity(turbine) Turbine%Ux_upstream=Ux Turbine%Uy_upstream=Uy Turbine%Uz_upstream=Uz - Turbine%Uref=sqrt_prec(Ux**2.0+Uy**2.0+Uz**2.0) + Turbine%Uref=sqrt(Ux**2.0+Uy**2.0+Uz**2.0) return diff --git a/src/acl_utils.f90 b/src/acl_utils.f90 index 1c449b66e..d3f9e6947 100644 --- a/src/acl_utils.f90 +++ b/src/acl_utils.f90 @@ -6,7 +6,6 @@ module actuator_line_model_utils use decomp_2d_constants, only: mytype use param, only: zero, one, two - use dbg_schemes, only: sqrt_prec, cos_prec, exp_prec, sin_prec implicit none public QuatRot, cross, IsoKernel, AnIsoKernel, int2str @@ -97,7 +96,7 @@ subroutine QuatRot(vx,vy,vz,Theta,Rx,Ry,Rz,Ox,Oy,Oz,vRx,vRy,vRz) ! vR: Rotated vector ! Force normalize nR - RMag=sqrt_prec(Rx**2.0+Ry**2.0+Rz**2.0) + RMag=sqrt(Rx**2.0+Ry**2.0+Rz**2.0) nRx=Rx/RMag nRy=Ry/RMag nRz=Rz/RMag @@ -109,7 +108,7 @@ subroutine QuatRot(vx,vy,vz,Theta,Rx,Ry,Rz,Ox,Oy,Oz,vRx,vRy,vRz) p=reshape([zero,vOx,vOy,vOz],[4,1]) ! Rotation quaternion and conjugate - q=(/cos_prec(Theta/2),nRx*sin_prec(Theta/2),nRy*sin_prec(Theta/2),nRz*sin_prec(Theta/2)/) + q=(/cos(Theta/2),nRx*sin(Theta/2),nRy*sin(Theta/2),nRz*sin(Theta/2)/) qbar=(/q(1),-q(2),-q(3),-q(4)/) QL=transpose(reshape((/q(1), -q(2), -q(3), -q(4), & @@ -148,7 +147,7 @@ subroutine IDW(Ncol,Xcol,Ycol,Zcol,Fxcol,Fycol,Fzcol,p,Xmesh,Ymesh,Zmesh,Fxmesh, wsum=zero do i=1,Ncol - d(i)=sqrt_prec((Xcol(i)-Xmesh)**2+(Ycol(i)-Ymesh)**2+(Zcol(i)-Zmesh)**2) + d(i)=sqrt((Xcol(i)-Xmesh)**2+(Ycol(i)-Ymesh)**2+(Zcol(i)-Zmesh)**2) w(i)=one/d(i)**p wsum=wsum+w(i) enddo @@ -184,9 +183,9 @@ real(mytype) function IsoKernel(dr,epsilon_par,dim) real(mytype), intent(in) :: dr, epsilon_par if (dim==2) then - IsoKernel = one/(epsilon_par**2*pi)*exp_prec(-(dr/epsilon_par)**2.0) + IsoKernel = one/(epsilon_par**2*pi)*exp(-(dr/epsilon_par)**2.0) else if (dim==3) then - IsoKernel = one/(epsilon_par**3.0*pi**1.5)*exp_prec(-(dr/epsilon_par)**2.0) + IsoKernel = one/(epsilon_par**3.0*pi**1.5)*exp(-(dr/epsilon_par)**2.0) else write(*,*) "1D source not implemented" stop @@ -211,7 +210,7 @@ real(mytype) function AnIsoKernel(dx,dy,dz,nx,ny,nz,tx,ty,tz,sx,sy,sz,ec,et,es) s=dx*sx+dy*sy+dz*sz ! Spanwise projection if (abs(s)<=es) then - AnIsoKernel = exp_prec(-((n/et)**2.0+(t/ec)**2.0))/(ec*et*pi) + AnIsoKernel = exp(-((n/et)**2.0+(t/ec)**2.0))/(ec*et*pi) else AnIsoKernel = zero endif diff --git a/src/case.f90 b/src/case.f90 index 239ef41cf..96a20e21b 100644 --- a/src/case.f90 +++ b/src/case.f90 @@ -12,7 +12,6 @@ module case use tgv use cyl use hill - use dbg_schemes use channel use mixlayer use gravitycur diff --git a/src/dynstall_legacy.f90 b/src/dynstall_legacy.f90 index 1683e958c..919975ea8 100644 --- a/src/dynstall_legacy.f90 +++ b/src/dynstall_legacy.f90 @@ -84,7 +84,6 @@ subroutine LB_EvalIdealCL(AOA,AOA0,CLa,RefFlag,CLID) use decomp_2d_constants, only : mytype use param, only: one, two, thirty - use dbg_schemes, only: sin_prec implicit none real(mytype) :: AOA, AOA0, CLa @@ -120,7 +119,7 @@ subroutine LB_EvalIdealCL(AOA,AOA0,CLa,RefFlag,CLID) if (abs(aID) 1.0e-4) then + ! if (abs(out(i,j,k)) > 1.0e-4) then ! write(*,*) 'AFTER',i,j,k,out(i,j,k),xyzk ! end if @@ -416,8 +411,6 @@ end subroutine poisson_000 subroutine poisson_100(rhs) - use dbg_schemes, only: abs_prec - implicit none real(mytype), dimension(:,:,:), intent(INOUT) :: rhs @@ -428,10 +421,6 @@ subroutine poisson_100(rhs) integer :: nx,ny,nz, i,j,k, itmp - complex(mytype) :: cx - real(mytype) :: rl, iy - external cx, rl, iy - 100 format(1x,a8,3I4,2F12.6) nx = nx_global - 1 @@ -471,7 +460,7 @@ subroutine poisson_100(rhs) do k = sp%xst(3), sp%xen(3) do j = sp%xst(2), sp%xen(2) do i = sp%xst(1), sp%xen(1) - if (abs_prec(cw1(i,j,k)) > 1.0e-4) then + if (abs(cw1(i,j,k)) > 1.0e-4) then write(*,100) 'START', i, j, k, cw1(i,j,k) end if end do @@ -490,7 +479,7 @@ subroutine poisson_100(rhs) cw1(i,j,k) = cx(tmp1 * bz(k) + tmp2 * az(k), & tmp2 * bz(k) - tmp1 * az(k)) #ifdef DEBG - if (abs_prec(cw1(i,j,k)) > 1.0e-4) & + if (abs(cw1(i,j,k)) > 1.0e-4) & write(*,100) 'after z',i,j,k,cw1(i,j,k) #endif end do @@ -507,7 +496,7 @@ subroutine poisson_100(rhs) tmp2 * by(j) - tmp1 * ay(j)) if (j > (ny/2+1)) cw1(i,j,k) = -cw1(i,j,k) #ifdef DEBG - if (abs_prec(cw1(i,j,k)) > 1.0e-4) & + if (abs(cw1(i,j,k)) > 1.0e-4) & write(*,100) 'after y',i,j,k,cw1(i,j,k) #endif end do @@ -540,7 +529,7 @@ subroutine poisson_100(rhs) do k = sp%xst(3), sp%xen(3) do j = sp%xst(2), sp%xen(2) do i = sp%xst(1), sp%xen(1) - if (abs_prec(cw1b(i,j,k)) > 1.0e-4) then + if (abs(cw1b(i,j,k)) > 1.0e-4) then write(*,100) 'after x',i,j,k,cw1b(i,j,k) end if end do @@ -555,20 +544,20 @@ subroutine poisson_100(rhs) tmp1 = rl(kxyz(i,j,k)) tmp2 = iy(kxyz(i,j,k)) ! CANNOT DO A DIVISION BY ZERO - if ((abs_prec(tmp1) < epsilon).and.(abs_prec(tmp2) < epsilon)) then + if ((abs(tmp1) < epsilon).and.(abs(tmp2) < epsilon)) then cw1b(i,j,k)=cx(zero, zero) end if - if ((abs_prec(tmp1) < epsilon).and.(abs_prec(tmp2) >= epsilon)) then + if ((abs(tmp1) < epsilon).and.(abs(tmp2) >= epsilon)) then cw1b(i,j,k)=cx(zero, iy(cw1b(i,j,k)) / (-tmp2)) end if - if ((abs_prec(tmp1) >= epsilon).and.(abs_prec(tmp2) < epsilon)) then + if ((abs(tmp1) >= epsilon).and.(abs(tmp2) < epsilon)) then cw1b(i,j,k)=cx(rl(cw1b(i,j,k)) / (-tmp1), zero) end if - if ((abs_prec(tmp1) >= epsilon).and.(abs_prec(tmp2) >= epsilon)) then + if ((abs(tmp1) >= epsilon).and.(abs(tmp2) >= epsilon)) then cw1b(i,j,k)=cx(rl(cw1b(i,j,k)) / (-tmp1), iy(cw1b(i,j,k)) / (-tmp2)) end if #ifdef DEBG - if (abs_prec(cw1b(i,j,k)) > 1.0e-4) & + if (abs(cw1b(i,j,k)) > 1.0e-4) & write(*,100) 'AFTER',i,j,k,cw1b(i,j,k) #endif end do @@ -603,7 +592,7 @@ subroutine poisson_100(rhs) do k = sp%xst(3),sp%xen(3) do j = sp%xst(2),sp%xen(2) do i = sp%xst(1),sp%xen(1) - if (abs_prec(cw1(i,j,k)) > 1.0e-4) then + if (abs(cw1(i,j,k)) > 1.0e-4) then write(*,100) 'AFTER X',i,j,k,cw1(i,j,k) end if end do @@ -621,7 +610,7 @@ subroutine poisson_100(rhs) tmp2 * by(j) + tmp1 * ay(j)) if (j > (ny/2+1)) cw1(i,j,k) = -cw1(i,j,k) #ifdef DEBG - if (abs_prec(cw1(i,j,k)) > 1.0e-4) & + if (abs(cw1(i,j,k)) > 1.0e-4) & write(*,100) 'AFTER Y',i,j,k,cw1(i,j,k) #endif end do @@ -637,7 +626,7 @@ subroutine poisson_100(rhs) cw1(i,j,k) = cx(tmp1 * bz(k) - tmp2 * az(k), & tmp2 * bz(k) + tmp1 * az(k)) #ifdef DEBG - if (abs_prec(cw1(i,j,k)) > 1.0e-4) & + if (abs(cw1(i,j,k)) > 1.0e-4) & write(*,100) 'END',i,j,k,cw1(i,j,k) #endif end do @@ -674,8 +663,6 @@ end subroutine poisson_100 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine poisson_010(rhs) - use dbg_schemes, only: abs_prec - implicit none real(mytype), dimension(:,:,:), intent(INOUT) :: rhs @@ -686,10 +673,6 @@ subroutine poisson_010(rhs) integer :: nx,ny,nz, i,j,k - complex(mytype) :: cx - real(mytype) :: rl, iy - external cx, rl, iy - 100 format(1x,a8,3I4,2F12.6) nx = nx_global @@ -727,7 +710,7 @@ subroutine poisson_010(rhs) do k = sp%xst(3), sp%xen(3) do j = sp%xst(2), sp%xen(2) do i = sp%xst(1), sp%xen(1) - if (abs_prec(cw1(i,j,k)) > 1.0e-4) then + if (abs(cw1(i,j,k)) > 1.0e-4) then write(*,100) 'START',i,j,k,cw1(i,j,k) end if end do @@ -746,7 +729,7 @@ subroutine poisson_010(rhs) cw1(i,j,k) = cx(tmp1 * bz(k) + tmp2 * az(k), & tmp2 * bz(k) - tmp1 * az(k)) #ifdef DEBG - if (abs_prec(cw1(i,j,k)) > 1.0e-4) & + if (abs(cw1(i,j,k)) > 1.0e-4) & write(*,100) 'after z',i,j,k,cw1(i,j,k) #endif end do @@ -763,7 +746,7 @@ subroutine poisson_010(rhs) tmp2 * bx(i) - tmp1 * ax(i)) if (i.gt.(nx/2+1)) cw1(i,j,k)=-cw1(i,j,k) #ifdef DEBG - if (abs_prec(cw1(i,j,k)) > 1.0e-4) & + if (abs(cw1(i,j,k)) > 1.0e-4) & write(*,100) 'after x',i,j,k,cw1(i,j,k) #endif end do @@ -799,7 +782,7 @@ subroutine poisson_010(rhs) do k = sp%yst(3), sp%yen(3) do j = sp%yst(2), sp%yen(2) do i = sp%yst(1), sp%yen(1) - if (abs_prec(cw2b(i,j,k)) > 1.0e-4) then + if (abs(cw2b(i,j,k)) > 1.0e-4) then write(*,100) 'after y',i,j,k,cw2b(i,j,k) write(*,*)kxyz(i,j,k) end if @@ -818,16 +801,16 @@ subroutine poisson_010(rhs) tmp1 = rl(kxyz(i,j,k)) tmp2 = iy(kxyz(i,j,k)) !CANNOT DO A DIVISION BY ZERO - if ((abs_prec(tmp1) < epsilon).and.(abs_prec(tmp2) < epsilon)) then + if ((abs(tmp1) < epsilon).and.(abs(tmp2) < epsilon)) then cw2b(i,j,k) = cx(zero, zero) end if - if ((abs_prec(tmp1) < epsilon).and.(abs_prec(tmp2) >= epsilon)) then + if ((abs(tmp1) < epsilon).and.(abs(tmp2) >= epsilon)) then cw2b(i,j,k) = cx(zero, iy(cw2b(i,j,k)) / (-tmp2)) end if - if ((abs_prec(tmp1) >= epsilon).and.(abs_prec(tmp2) < epsilon)) then + if ((abs(tmp1) >= epsilon).and.(abs(tmp2) < epsilon)) then cw2b(i,j,k) = cx(rl(cw2b(i,j,k)) / (-tmp1), zero) end if - if ((abs_prec(tmp1) >= epsilon).and.(abs_prec(tmp2) >= epsilon)) then + if ((abs(tmp1) >= epsilon).and.(abs(tmp2) >= epsilon)) then cw2b(i,j,k) = cx(rl(cw2b(i,j,k)) / (-tmp1), iy(cw2b(i,j,k)) / (-tmp2)) end if end do @@ -926,7 +909,7 @@ subroutine poisson_010(rhs) do k = sp%yst(3), sp%yen(3) do j = sp%yst(2), sp%yen(2) do i = sp%yst(1), sp%yen(1) - if (abs_prec(cw2b(i,j,k)) > 1.0e-4) then + if (abs(cw2b(i,j,k)) > 1.0e-4) then write(*,100) 'AFTER',i,j,k,cw2b(i,j,k) write(*,*)kxyz(i,j,k) end if @@ -965,7 +948,7 @@ subroutine poisson_010(rhs) do k = sp%xst(3),sp%xen(3) do j = sp%xst(2),sp%xen(2) do i = sp%xst(1),sp%xen(1) - if (abs_prec(cw1(i,j,k)) > 1.0e-4) then + if (abs(cw1(i,j,k)) > 1.0e-4) then write(*,100) 'AFTER Y',i,j,k,cw1(i,j,k) end if end do @@ -983,7 +966,7 @@ subroutine poisson_010(rhs) tmp2 * bx(i) + tmp1 * ax(i)) if (i > (nx/2 + 1)) cw1(i,j,k) = -cw1(i,j,k) #ifdef DEBG - if (abs_prec(cw1(i,j,k)) > 1.0e-4) & + if (abs(cw1(i,j,k)) > 1.0e-4) & write(*,100) 'AFTER X',i,j,k,cw1(i,j,k) #endif end do @@ -999,7 +982,7 @@ subroutine poisson_010(rhs) cw1(i,j,k) = cx(tmp1 * bz(k) - tmp2 * az(k), & tmp2 * bz(k) + tmp1 * az(k)) #ifdef DEBG - if (abs_prec(cw1(i,j,k)) > 1.0e-4) & + if (abs(cw1(i,j,k)) > 1.0e-4) & write(*,100) 'END',i,j,k,cw1(i,j,k) #endif end do @@ -1034,10 +1017,6 @@ end subroutine poisson_010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine poisson_11x(rhs) - use dbg_schemes, only: abs_prec - use mpi - - implicit none real(mytype), dimension(:,:,:), intent(INOUT) :: rhs @@ -1048,9 +1027,6 @@ subroutine poisson_11x(rhs) integer :: nx,ny,nz, i,j,k - complex(mytype) :: cx - real(mytype) :: rl, iy - external cx, rl, iy #ifdef DEBG real(mytype) :: dep, dep1 integer :: code @@ -1230,16 +1206,16 @@ subroutine poisson_11x(rhs) tmp1 = rl(kxyz(i,j,k)) tmp2 = iy(kxyz(i,j,k)) !CANNOT DO A DIVISION BY ZERO - if ((abs_prec(tmp1) < epsilon).and.(abs_prec(tmp2) < epsilon)) then + if ((abs(tmp1) < epsilon).and.(abs(tmp2) < epsilon)) then cw1b(i,j,k) = cx(zero, zero) end if - if ((abs_prec(tmp1) < epsilon).and.(abs_prec(tmp2) >= epsilon)) then + if ((abs(tmp1) < epsilon).and.(abs(tmp2) >= epsilon)) then cw1b(i,j,k) = cx(zero, iy(cw1b(i,j,k)) / (-tmp2)) end if - if ((abs_prec(tmp1) >= epsilon).and.(abs_prec(tmp2) < epsilon)) then + if ((abs(tmp1) >= epsilon).and.(abs(tmp2) < epsilon)) then cw1b(i,j,k) = cx(rl(cw1b(i,j,k)) / (-tmp1), zero) end if - if ((abs_prec(tmp1) >= epsilon).and.(abs_prec(tmp2) >= epsilon)) then + if ((abs(tmp1) >= epsilon).and.(abs(tmp2) >= epsilon)) then cw1b(i,j,k) = cx(real(cw1b(i,j,k)) / (-tmp1), iy(cw1b(i,j,k)) / (-tmp2)) end if end do @@ -1492,7 +1468,6 @@ end subroutine poisson_11x subroutine abxyz(ax,ay,az,bx,by,bz,nx,ny,nz,bcx,bcy,bcz) use param - use dbg_schemes, only: sin_prec, cos_prec implicit none @@ -1506,42 +1481,42 @@ subroutine abxyz(ax,ay,az,bx,by,bz,nx,ny,nz,bcx,bcy,bcz) if (bcx == 0) then do i = 1, nx - ax(i) = sin_prec(real(i-1, kind=mytype)*PI/real(nx, kind=mytype)) - bx(i) = cos_prec(real(i-1, kind=mytype)*PI/real(nx, kind=mytype)) + ax(i) = sin(real(i-1, kind=mytype)*PI/real(nx, kind=mytype)) + bx(i) = cos(real(i-1, kind=mytype)*PI/real(nx, kind=mytype)) end do elseif (bcx == 1) then do i = 1, nx - ax(i) = sin_prec(real(i-1, kind=mytype)*PI*half/ & + ax(i) = sin(real(i-1, kind=mytype)*PI*half/ & real(nx, kind=mytype)) - bx(i) = cos_prec(real(i-1, kind=mytype)*PI*half/ & + bx(i) = cos(real(i-1, kind=mytype)*PI*half/ & real(nx, kind=mytype)) end do end if if (bcy == 0) then do j = 1, ny - ay(j) = sin_prec(real(j-1, kind=mytype)*PI/real(ny, kind=mytype)) - by(j) = cos_prec(real(j-1, kind=mytype)*PI/real(ny, kind=mytype)) + ay(j) = sin(real(j-1, kind=mytype)*PI/real(ny, kind=mytype)) + by(j) = cos(real(j-1, kind=mytype)*PI/real(ny, kind=mytype)) end do elseif (bcy == 1) then do j = 1, ny - ay(j) = sin_prec(real(j-1, kind=mytype)*PI*half/ & + ay(j) = sin(real(j-1, kind=mytype)*PI*half/ & real(ny, kind=mytype)) - by(j) = cos_prec(real(j-1, kind=mytype)*PI*half/ & + by(j) = cos(real(j-1, kind=mytype)*PI*half/ & real(ny, kind=mytype)) end do end if if (bcz == 0) then do k = 1, nz - az(k) = sin_prec(real(k-1, kind=mytype)*PI/real(nz, kind=mytype)) - bz(k) = cos_prec(real(k-1, kind=mytype)*PI/real(nz, kind=mytype)) + az(k) = sin(real(k-1, kind=mytype)*PI/real(nz, kind=mytype)) + bz(k) = cos(real(k-1, kind=mytype)*PI/real(nz, kind=mytype)) end do elseif (bcz == 1) then do k = 1, nz - az(k) = sin_prec(real(k-1, kind=mytype)*PI*half/ & + az(k) = sin(real(k-1, kind=mytype)*PI*half/ & real(nz, kind=mytype)) - bz(k) = cos_prec(real(k-1, kind=mytype)*PI*half/ & + bz(k) = cos(real(k-1, kind=mytype)*PI*half/ & real(nz, kind=mytype)) end do end if @@ -1562,7 +1537,6 @@ subroutine waves () use decomp_2d use variables use decomp_2d_fft - use dbg_schemes, only: sin_prec, cos_prec implicit none @@ -1580,10 +1554,6 @@ subroutine waves () real(mytype) :: ytt_rl,xtt_rl,ztt_rl,yt1_rl,xt1_rl,zt1_rl real(mytype) :: xtt1_rl,ytt1_rl,ztt1_rl - complex(mytype) :: cx - real(mytype) :: rl, iy - external cx, rl, iy - xkx = zero xk2 = zero yky = zero @@ -1595,8 +1565,8 @@ subroutine waves () if (bcx == 0) then do i = 1, nx/2 + 1 w = twopi * (i-1) / nx - wp = acix6 * two * dx * sin_prec(w * half) + bcix6 * two * dx * sin_prec(three * half * w) - wp = wp / (one + two * alcaix6 * cos_prec(w)) + wp = acix6 * two * dx * sin(w * half) + bcix6 * two * dx * sin(three * half * w) + wp = wp / (one + two * alcaix6 * cos(w)) ! xkx(i) = cx_one_one * (nx * wp / xlx) exs(i) = cx_one_one * (nx * w / xlx) @@ -1611,8 +1581,8 @@ subroutine waves () else do i = 1, nx w = twopi * half * (i-1) / nxm - wp = acix6 * two * dx * sin_prec(w * half) +(bcix6 * two * dx) * sin_prec(three * half * w) - wp = wp / (one + two * alcaix6 * cos_prec(w)) + wp = acix6 * two * dx * sin(w * half) +(bcix6 * two * dx) * sin(three * half * w) + wp = wp / (one + two * alcaix6 * cos(w)) ! xkx(i) = cx_one_one * nxm * wp / xlx exs(i) = cx_one_one * nxm * w / xlx @@ -1628,8 +1598,8 @@ subroutine waves () if (bcy == 0) then do j = 1, ny/2 + 1 w = twopi * (j-1) / ny - wp = aciy6 * two * dy * sin_prec(w * half) + bciy6 * two * dy * sin_prec(three * half * w) - wp = wp / (one + two * alcaiy6 * cos_prec(w)) + wp = aciy6 * two * dy * sin(w * half) + bciy6 * two * dy * sin(three * half * w) + wp = wp / (one + two * alcaiy6 * cos(w)) ! if (istret == 0) yky(j) = cx_one_one * (ny * wp / yly) if (istret /= 0) yky(j) = cx_one_one * (ny * wp) @@ -1645,8 +1615,8 @@ subroutine waves () else do j = 1, ny w = twopi * half * (j-1) / nym - wp = aciy6 * two * dy * sin_prec(w * half) +(bciy6 * two *dy) * sin_prec(three * half * w) - wp = wp / (one + two * alcaiy6 * cos_prec(w)) + wp = aciy6 * two * dy * sin(w * half) +(bciy6 * two *dy) * sin(three * half * w) + wp = wp / (one + two * alcaiy6 * cos(w)) ! if (istret == 0) yky(j) = cx_one_one * (nym * wp / yly) if (istret /= 0) yky(j) = cx_one_one * (nym * wp) @@ -1663,8 +1633,8 @@ subroutine waves () if (bcz == 0) then do k = 1, nz/2 + 1 w = twopi * (k-1) / nz - wp = aciz6 * two * dz * sin_prec(w * half) + (bciz6 * two * dz) * sin_prec(three * half * w) - wp = wp / (one + two * alcaiz6 * cos_prec(w)) + wp = aciz6 * two * dz * sin(w * half) + (bciz6 * two * dz) * sin(three * half * w) + wp = wp / (one + two * alcaiz6 * cos(w)) ! zkz(k) = cx_one_one * (nz * wp / zlz) ezs(k) = cx_one_one * (nz * w / zlz) @@ -1675,10 +1645,10 @@ subroutine waves () do k= 1, nz/2 + 1 w = pi * (k-1) / nzm w1 = pi * (nzm-k+1) / nzm - wp = aciz6 * two * dz * sin_prec(w * half)+(bciz6 * two * dz) * sin_prec(three * half * w) - wp = wp / (one + two * alcaiz6 * cos_prec(w)) - w1p = aciz6 * two * dz * sin_prec(w1 * half) + (bciz6 * two * dz) * sin_prec(three * half * w1) - w1p = w1p / (one + two * alcaiz6 * cos_prec(w1)) + wp = aciz6 * two * dz * sin(w * half)+(bciz6 * two * dz) * sin(three * half * w) + wp = wp / (one + two * alcaiz6 * cos(w)) + w1p = aciz6 * two * dz * sin(w1 * half) + (bciz6 * two * dz) * sin(three * half * w1) + w1p = w1p / (one + two * alcaiz6 * cos(w1)) ! zkz(k) = cx(nzm * wp / zlz, -nzm * w1p / zlz) ezs(k) = cx(nzm * w / zlz, nzm * w1 / zlz) @@ -1701,21 +1671,21 @@ subroutine waves () rlexs = rl(exs(i)) * dx ! xtt_rl = two * & - (bicix6 * cos_prec(rlexs * onepfive) + cicix6 * cos_prec(rlexs * twopfive) + dicix6 * cos_prec(rlexs * threepfive)) + (bicix6 * cos(rlexs * onepfive) + cicix6 * cos(rlexs * twopfive) + dicix6 * cos(rlexs * threepfive)) ! ytt_rl = two * & - (biciy6 * cos_prec(rleys * onepfive) + ciciy6 * cos_prec(rleys * twopfive) + diciy6 * cos_prec(rleys * threepfive)) + (biciy6 * cos(rleys * onepfive) + ciciy6 * cos(rleys * twopfive) + diciy6 * cos(rleys * threepfive)) ! ztt_rl = two * & - (biciz6 * cos_prec(rlezs * onepfive) + ciciz6 * cos_prec(rlezs * twopfive) + diciz6 * cos_prec(rlezs * threepfive)) + (biciz6 * cos(rlezs * onepfive) + ciciz6 * cos(rlezs * twopfive) + diciz6 * cos(rlezs * threepfive)) ! - xtt1_rl = two * aicix6 * cos_prec(rlexs * half) - ytt1_rl = two * aiciy6 * cos_prec(rleys * half) - ztt1_rl = two * aiciz6 * cos_prec(rlezs * half) + xtt1_rl = two * aicix6 * cos(rlexs * half) + ytt1_rl = two * aiciy6 * cos(rleys * half) + ztt1_rl = two * aiciz6 * cos(rlezs * half) ! - xt1_rl = one + two * ailcaix6 * cos_prec(rlexs) - yt1_rl = one + two * ailcaiy6 * cos_prec(rleys) - zt1_rl = one + two * ailcaiz6 * cos_prec(rlezs) + xt1_rl = one + two * ailcaix6 * cos(rlexs) + yt1_rl = one + two * ailcaiy6 * cos(rleys) + zt1_rl = one + two * ailcaiz6 * cos(rlezs) ! xt2 = xk2(i) * ((((ytt1_rl + ytt_rl) / yt1_rl) * ((ztt1_rl + ztt_rl) / zt1_rl))**2) yt2 = yk2(j) * ((((xtt1_rl + xtt_rl) / xt1_rl) * ((ztt1_rl + ztt_rl) / zt1_rl))**2) @@ -1743,21 +1713,21 @@ subroutine waves () rlexs = rl(exs(i)) * dx ! xtt_rl = two * & - (bicix6 * cos_prec(rlexs * onepfive) + cicix6 * cos_prec(rlexs * twopfive) + dicix6 * cos_prec(rlexs * threepfive)) + (bicix6 * cos(rlexs * onepfive) + cicix6 * cos(rlexs * twopfive) + dicix6 * cos(rlexs * threepfive)) ! ytt_rl = two * & - (biciy6 * cos_prec(rleys * onepfive) + ciciy6 * cos_prec(rleys * twopfive) + diciy6 * cos_prec(rleys * threepfive)) + (biciy6 * cos(rleys * onepfive) + ciciy6 * cos(rleys * twopfive) + diciy6 * cos(rleys * threepfive)) ! ztt_rl = two * & - (biciz6 * cos_prec(rlezs * onepfive) + ciciz6 * cos_prec(rlezs * twopfive) + diciz6 * cos_prec(rlezs * threepfive)) + (biciz6 * cos(rlezs * onepfive) + ciciz6 * cos(rlezs * twopfive) + diciz6 * cos(rlezs * threepfive)) ! - xtt1_rl = two * aicix6 * cos_prec(rlexs * half) - ytt1_rl = two * aiciy6 * cos_prec(rleys * half) - ztt1_rl = two * aiciz6 * cos_prec(rlezs * half) + xtt1_rl = two * aicix6 * cos(rlexs * half) + ytt1_rl = two * aiciy6 * cos(rleys * half) + ztt1_rl = two * aiciz6 * cos(rlezs * half) ! - xt1_rl = one + two * ailcaix6 * cos_prec(rlexs) - yt1_rl = one + two * ailcaiy6 * cos_prec(rleys) - zt1_rl = one + two * ailcaiz6 * cos_prec(rlezs) + xt1_rl = one + two * ailcaix6 * cos(rlexs) + yt1_rl = one + two * ailcaiy6 * cos(rleys) + zt1_rl = one + two * ailcaiz6 * cos(rlezs) ! xt2 = xk2(i) * ((((ytt1_rl + ytt_rl) / yt1_rl) * ((ztt1_rl + ztt_rl) / zt1_rl))**2) yt2 = yk2(j) * ((((xtt1_rl + xtt_rl) / xt1_rl) * ((ztt1_rl + ztt_rl) / zt1_rl))**2) @@ -1785,26 +1755,26 @@ subroutine waves () rlexs = rl(exs(i)) * dx ! xtt_rl = two * & - (bicix6 * cos_prec(rlexs * onepfive) + cicix6 * cos_prec(rlexs * twopfive) + dicix6 * cos_prec(rlexs * threepfive)) + (bicix6 * cos(rlexs * onepfive) + cicix6 * cos(rlexs * twopfive) + dicix6 * cos(rlexs * threepfive)) ! ytt_rl = two * & - (biciy6 * cos_prec(rleys * onepfive) + ciciy6 * cos_prec(rleys * twopfive) + diciy6 * cos_prec(rleys * threepfive)) + (biciy6 * cos(rleys * onepfive) + ciciy6 * cos(rleys * twopfive) + diciy6 * cos(rleys * threepfive)) ! ztt = two * cx( & - biciz6 * cos_prec(rlezs * onepfive) + ciciz6 * cos_prec(rlezs * twopfive) + diciz6 * cos_prec(rlezs * threepfive),& - biciz6 * cos_prec(iyezs * onepfive) + ciciz6 * cos_prec(iyezs * twopfive) + diciz6 * cos_prec(iyezs * threepfive)) + biciz6 * cos(rlezs * onepfive) + ciciz6 * cos(rlezs * twopfive) + diciz6 * cos(rlezs * threepfive),& + biciz6 * cos(iyezs * onepfive) + ciciz6 * cos(iyezs * twopfive) + diciz6 * cos(iyezs * threepfive)) ! - xtt1_rl = two * aicix6 * cos_prec(rlexs * half) - ytt1_rl = two * aiciy6 * cos_prec(rleys * half) + xtt1_rl = two * aicix6 * cos(rlexs * half) + ytt1_rl = two * aiciy6 * cos(rleys * half) ! - ztt1 = two * cx(aiciz6 * cos_prec(rlezs * half),& - aiciz6 * cos_prec(iyezs * half)) + ztt1 = two * cx(aiciz6 * cos(rlezs * half),& + aiciz6 * cos(iyezs * half)) ! - xt1_rl = one + two * ailcaix6 * cos_prec(rlexs) - yt1_rl = one + two * ailcaiy6 * cos_prec(rleys) + xt1_rl = one + two * ailcaix6 * cos(rlexs) + yt1_rl = one + two * ailcaiy6 * cos(rleys) ! - zt1 = cx((one + two * ailcaiz6 * cos_prec(rlezs)),& - (one + two * ailcaiz6 * cos_prec(iyezs))) + zt1 = cx((one + two * ailcaiz6 * cos(rlezs)),& + (one + two * ailcaiz6 * cos(iyezs))) ! tmp1 = cx(rl(ztt1 + ztt) / rl(zt1),& iy(ztt1 + ztt) / iy(zt1)) @@ -1848,11 +1818,9 @@ subroutine matrice_refinement() use variables use param use var - use MPI use derivX use derivY - use derivZ - use dbg_schemes, only: cos_prec + use derivZ implicit none @@ -1870,11 +1838,6 @@ subroutine matrice_refinement() complex(mytype) :: ytt,xtt,ztt,yt1,xt1,yt2,xt2 complex(mytype) :: xtt1,ytt1,ztt1,zt1,zt2,tmp1,tmp2,tmp3 - - complex(mytype) :: cx - real(mytype) :: rl, iy - external cx, rl, iy - real(mytype) :: xtt_rl, xtt1_rl, xt1_rl real(mytype) :: rlexs @@ -1889,9 +1852,9 @@ subroutine matrice_refinement() do i = sp%yst(1),sp%yen(1) ! rlexs = rl(exs(i)) * dx - xtt_rl=two * (bicix6 * cos_prec(rlexs * onepfive) + cicix6 * cos_prec(rlexs * twopfive) + dicix6 * cos_prec(rlexs * threepfive)) - xtt1_rl=two * aicix6 * cos_prec(rlexs * half) - xt1_rl= one + two * ailcaix6 * cos_prec(rlexs) + xtt_rl=two * (bicix6 * cos(rlexs * onepfive) + cicix6 * cos(rlexs * twopfive) + dicix6 * cos(rlexs * threepfive)) + xtt1_rl=two * aicix6 * cos(rlexs * half) + xt1_rl= one + two * ailcaix6 * cos(rlexs) ! transx_rl(i) = (xtt1_rl + xtt_rl) / xt1_rl transx_rl2(i) = transx_rl(i)**2 @@ -1903,9 +1866,9 @@ subroutine matrice_refinement() do j = sp%yst(2),sp%yen(2) ! rleys = rl(eys(j)) * dy - ytt_rl=two * (biciy6 * cos_prec(rleys * onepfive) + ciciy6 * cos_prec(rleys * twopfive) + diciy6 * cos_prec(rleys * threepfive)) - ytt1_rl=two * aiciy6 * cos_prec(rleys * half) - yt1_rl=one + two * ailcaiy6 * cos_prec(rleys) + ytt_rl=two * (biciy6 * cos(rleys * onepfive) + ciciy6 * cos(rleys * twopfive) + diciy6 * cos(rleys * threepfive)) + ytt1_rl=two * aiciy6 * cos(rleys * half) + yt1_rl=one + two * ailcaiy6 * cos(rleys) ! transy_rl(j) = (ytt1_rl + ytt_rl) / yt1_rl transy_rl2(j) = transy_rl(j)**2 @@ -1918,9 +1881,9 @@ subroutine matrice_refinement() do k = sp%yst(3),sp%yen(3) ! rlezs = rl(ezs(k)) * dz - ztt_rl=two * (biciz6 * cos_prec(rlezs * onepfive) + ciciz6 * cos_prec(rlezs * twopfive) + diciz6 * cos_prec(rlezs * threepfive)) - ztt1_rl=two * aiciz6 * cos_prec(rlezs * half) - zt1_rl=one + two * ailcaiz6 * cos_prec(rlezs) + ztt_rl=two * (biciz6 * cos(rlezs * onepfive) + ciciz6 * cos(rlezs * twopfive) + diciz6 * cos(rlezs * threepfive)) + ztt1_rl=two * aiciz6 * cos(rlezs * half) + zt1_rl=one + two * ailcaiz6 * cos(rlezs) ! transz_rl(k) = (ztt1_rl + ztt_rl) / zt1_rl transz_rl2(k) = transz_rl(k)**2 @@ -1936,12 +1899,12 @@ subroutine matrice_refinement() ! rlezs = rl(ezs(k)) * dz iyezs = iy(ezs(k)) * dz - ztt = two * cx(biciz6 * cos_prec(rlezs * onepfive) + ciciz6 * cos_prec(rlezs * twopfive), & - biciz6 * cos_prec(iyezs * onepfive) + ciciz6 * cos_prec(iyezs * twopfive)) - ztt1 = two * cx(aiciz6 * cos_prec(rlezs * half),& - aiciz6 * cos_prec(iyezs * half)) - zt1 = cx(one + two * ailcaiz6 * cos_prec(rlezs),& - one + two * ailcaiz6 * cos_prec(iyezs)) + ztt = two * cx(biciz6 * cos(rlezs * onepfive) + ciciz6 * cos(rlezs * twopfive), & + biciz6 * cos(iyezs * onepfive) + ciciz6 * cos(iyezs * twopfive)) + ztt1 = two * cx(aiciz6 * cos(rlezs * half),& + aiciz6 * cos(iyezs * half)) + zt1 = cx(one + two * ailcaiz6 * cos(rlezs),& + one + two * ailcaiz6 * cos(iyezs)) ! transz_rl(k) = rl(ztt1 + ztt) / rl(zt1) transz_rl2(k) = transz_rl(k)**2 diff --git a/src/tools.f90 b/src/tools.f90 index f7bed7e33..66055dc16 100644 --- a/src/tools.f90 +++ b/src/tools.f90 @@ -7,6 +7,7 @@ module tools use decomp_2d_constants use decomp_2d_mpi use decomp_2d + use utilities implicit none @@ -35,7 +36,6 @@ subroutine test_scalar_min_max(phi) use param use var use mpi - use dbg_schemes, only: abs_prec implicit none @@ -68,7 +68,7 @@ subroutine test_scalar_min_max(phi) write(*,*) 'Phi'//char(48+is)//' min max=', real(phimin1,4), real(phimax1,4) - if (abs_prec(phimax1) > 100._mytype) then !if phi control turned off + if (abs(phimax1) > 100._mytype) then !if phi control turned off write(*,*) 'Scalar diverged! SIMULATION IS STOPPED!' call MPI_ABORT(MPI_COMM_WORLD,code,ierror); stop endif @@ -85,7 +85,6 @@ subroutine test_speed_min_max(ux,uy,uz) use param use var use mpi - use dbg_schemes, only: abs_prec implicit none @@ -127,7 +126,7 @@ subroutine test_speed_min_max(ux,uy,uz) write(*,*) 'U,V,W max=',real(uxmax1,4),real(uymax1,4),real(uzmax1,4) !print *,'CFL=',real(abs(max(uxmax1,uymax1,uzmax1)*dt)/min(dx,dy,dz),4) - if((abs_prec(uxmax1)>=onehundred).or.(abs_prec(uymax1)>=onehundred).OR.(abs_prec(uzmax1)>=onehundred)) then + if((abs(uxmax1)>=onehundred).or.(abs(uymax1)>=onehundred).OR.(abs(uzmax1)>=onehundred)) then write(*,*) 'Velocity diverged! SIMULATION IS STOPPED!' call MPI_ABORT(MPI_COMM_WORLD,code,ierror) stop @@ -1092,7 +1091,6 @@ subroutine cfl_compute(uxmax,uymax,uzmax) use param use variables use var - use dbg_schemes, only: sqrt_prec, abs_prec implicit none @@ -1103,7 +1101,7 @@ subroutine cfl_compute(uxmax,uymax,uzmax) real(mytype) :: visc ! Set the constants (this is true for periodic boundaries) - sigma_conv=[zero, sqrt_prec(three), 2.85_mytype] + sigma_conv=[zero, sqrt(three), 2.85_mytype] sigma_diff=[two, 2.5_mytype, 2.9_mytype] if(jles==0) then @@ -1114,13 +1112,13 @@ subroutine cfl_compute(uxmax,uymax,uzmax) ! This is considering 1D peridic boundaries ! Do x-direction - cfl_x_adv =abs_prec(uxmax) * dt / dx + cfl_x_adv =abs(uxmax) * dt / dx cfl_x_diff = visc * dt / dx**2 ! Do y-direction - cfl_y_adv = abs_prec(uymax) * dt / dy + cfl_y_adv = abs(uymax) * dt / dy cfl_y_diff = visc * dt / dy**2 ! Do z-direction - cfl_z_adv = abs_prec(uzmax) * dt / dz + cfl_z_adv = abs(uzmax) * dt / dz cfl_z_diff = visc * dt / dz**2 ! So far we will focus on uniform grids @@ -1132,7 +1130,7 @@ subroutine cfl_compute(uxmax,uymax,uzmax) 1003 format('CFL y-direction (Adv and Diff) =',F9.4,',',F9.4) write(*,1004) cfl_z_adv, cfl_z_diff 1004 format('CFL z-direction (Adv and Diff) =',F9.4,',',F9.4) - cfl_conv_lim = sigma_conv(itimescheme) / sqrt_prec(three) + cfl_conv_lim = sigma_conv(itimescheme) / sqrt(three) cfl_diff_lim = sigma_diff(itimescheme) / six write(*,1005) cfl_conv_lim, cfl_diff_lim write(*,*) ' ' @@ -1151,7 +1149,6 @@ subroutine stretching() use param use var use mpi - use dbg_schemes, only: abs_prec, sqrt_prec, sin_prec, cos_prec, tan_prec, atan_prec implicit none @@ -1340,7 +1337,7 @@ subroutine inversion5_v1(aaa_in,eee,spI) use param use var use mpi - use dbg_schemes, only: abs_prec + use utilities implicit none @@ -1362,10 +1359,6 @@ subroutine inversion5_v1(aaa_in,eee,spI) real(mytype) :: tmp1,tmp2,tmp3,tmp4 - complex(mytype) :: cx - real(mytype) :: rl, iy - external cx, rl, iy - aaa = aaa_in do i = 1, 2 @@ -1397,12 +1390,12 @@ subroutine inversion5_v1(aaa_in,eee,spI) do k = spI%yst(3), spI%yen(3) do j = spI%yst(1), spI%yen(1) - if (abs_prec(rl(aaa(j,ny/2-1,k,3))) > epsilon) then + if (abs(rl(aaa(j,ny/2-1,k,3))) > epsilon) then tmp1 = rl(aaa(j,ny/2,k,2)) / rl(aaa(j,ny/2-1,k,3)) else tmp1 = zero endif - if (abs_prec(iy(aaa(j,ny/2-1,k,3))) > epsilon) then + if (abs(iy(aaa(j,ny/2-1,k,3))) > epsilon) then tmp2 = iy(aaa(j,ny/2,k,2)) / iy(aaa(j,ny/2-1,k,3)) else tmp2 = zero @@ -1411,14 +1404,14 @@ subroutine inversion5_v1(aaa_in,eee,spI) b1(j,k) = cx(rl(aaa(j,ny/2,k,3)) - tmp1 * rl(aaa(j,ny/2-1,k,4)),& iy(aaa(j,ny/2,k,3)) - tmp2 * iy(aaa(j,ny/2-1,k,4))) - if (abs_prec(rl(b1(j,k))) > epsilon) then + if (abs(rl(b1(j,k))) > epsilon) then tmp1 = rl(sr(j,k)) / rl(b1(j,k)) tmp3 = rl(eee(j,ny/2,k)) / rl(b1(j,k)) - tmp1 * rl(eee(j,ny/2-1,k)) else tmp1 = zero tmp3 = zero endif - if (abs_prec(iy(b1(j,k))) > epsilon) then + if (abs(iy(b1(j,k))) > epsilon) then tmp2 = iy(sr(j,k)) / iy(b1(j,k)) tmp4 = iy(eee(j,ny/2,k)) / iy(b1(j,k)) - tmp2 * iy(eee(j,ny/2-1,k)) else @@ -1428,12 +1421,12 @@ subroutine inversion5_v1(aaa_in,eee,spI) a1(j,k) = cx(tmp1,tmp2) eee(j,ny/2,k) = cx(tmp3,tmp4) - if (abs_prec(rl(aaa(j,ny/2-1,k,3))) > epsilon) then + if (abs(rl(aaa(j,ny/2-1,k,3))) > epsilon) then tmp1 = one / rl(aaa(j,ny/2-1,k,3)) else tmp1 = zero endif - if (abs_prec(iy(aaa(j,ny/2-1,k,3))) > epsilon) then + if (abs(iy(aaa(j,ny/2-1,k,3))) > epsilon) then tmp2 = one / iy(aaa(j,ny/2-1,k,3)) else tmp2 = zero @@ -1449,12 +1442,12 @@ subroutine inversion5_v1(aaa_in,eee,spI) do i = ny/2 - 2, 1, -1 do k = spI%yst(3), spI%yen(3) do j = spI%yst(1), spI%yen(1) - if (abs_prec(rl(aaa(j,i,k,3))) > epsilon) then + if (abs(rl(aaa(j,i,k,3))) > epsilon) then tmp1 = one / rl(aaa(j,i,k,3)) else tmp1 = zero endif - if (abs_prec(iy(aaa(j,i,k,3))) > epsilon) then + if (abs(iy(aaa(j,i,k,3))) > epsilon) then tmp2 = one/iy(aaa(j,i,k,3)) else tmp2 = zero @@ -1484,7 +1477,7 @@ subroutine inversion5_v2(aaa,eee,spI) use param use var use MPI - use dbg_schemes, only: abs_prec + use utilities implicit none @@ -1506,10 +1499,6 @@ subroutine inversion5_v2(aaa,eee,spI) real(mytype) :: tmp1,tmp2,tmp3,tmp4 - complex(mytype) :: cx - real(mytype) :: rl, iy - external cx, rl, iy - do i = 1, 2 ja(i) = 4 - i jb(i) = 5 - i @@ -1538,12 +1527,12 @@ subroutine inversion5_v2(aaa,eee,spI) enddo do k = spI%yst(3), spI%yen(3) do j = spI%yst(1), spI%yen(1) - if (abs_prec(rl(aaa(j,nym-1,k,3))) > epsilon) then + if (abs(rl(aaa(j,nym-1,k,3))) > epsilon) then tmp1 = rl(aaa(j,nym,k,2)) / rl(aaa(j,nym-1,k,3)) else tmp1 = zero endif - if (abs_prec(iy(aaa(j,nym-1,k,3))) > epsilon) then + if (abs(iy(aaa(j,nym-1,k,3))) > epsilon) then tmp2 = iy(aaa(j,nym,k,2)) / iy(aaa(j,nym-1,k,3)) else tmp2 = zero @@ -1551,14 +1540,14 @@ subroutine inversion5_v2(aaa,eee,spI) sr(j,k) = cx(tmp1,tmp2) b1(j,k) = cx(rl(aaa(j,nym,k,3)) - tmp1 * rl(aaa(j,nym-1,k,4)),& iy(aaa(j,nym,k,3)) - tmp2 * iy(aaa(j,nym-1,k,4))) - if (abs_prec(rl(b1(j,k))) > epsilon) then + if (abs(rl(b1(j,k))) > epsilon) then tmp1 = rl(sr(j,k)) / rl(b1(j,k)) tmp3 = rl(eee(j,nym,k)) / rl(b1(j,k)) - tmp1 * rl(eee(j,nym-1,k)) else tmp1 = zero tmp3 = zero endif - if (abs_prec(iy(b1(j,k))) > epsilon) then + if (abs(iy(b1(j,k))) > epsilon) then tmp2 = iy(sr(j,k)) / iy(b1(j,k)) tmp4 = iy(eee(j,nym,k)) / iy(b1(j,k)) - tmp2 * iy(eee(j,nym-1,k)) else @@ -1568,12 +1557,12 @@ subroutine inversion5_v2(aaa,eee,spI) a1(j,k) = cx(tmp1, tmp2) eee(j,nym,k) = cx(tmp3, tmp4) - if (abs_prec(rl(aaa(j,nym-1,k,3))) > epsilon) then + if (abs(rl(aaa(j,nym-1,k,3))) > epsilon) then tmp1 = one / rl(aaa(j,nym-1,k,3)) else tmp1 = zero endif - if (abs_prec(iy(aaa(j,nym-1,k,3))) > epsilon) then + if (abs(iy(aaa(j,nym-1,k,3))) > epsilon) then tmp2 = one / iy(aaa(j,nym-1,k,3)) else tmp2 = zero @@ -1589,12 +1578,12 @@ subroutine inversion5_v2(aaa,eee,spI) do i = nym - 2, 1, -1 do k = spI%yst(3), spI%yen(3) do j = spI%yst(1), spI%yen(1) - if (abs_prec(rl(aaa(j,i,k,3))) > epsilon) then + if (abs(rl(aaa(j,i,k,3))) > epsilon) then tmp1 = one / rl(aaa(j,i,k,3)) else tmp1 = zero endif - if (abs_prec(iy(aaa(j,i,k,3))) > epsilon) then + if (abs(iy(aaa(j,i,k,3))) > epsilon) then tmp2 = one / iy(aaa(j,i,k,3)) else tmp2 = zero @@ -1623,7 +1612,6 @@ subroutine tripping(tb,ta) use decomp_2d_constants use decomp_2d_mpi use mpi - use dbg_schemes, only: sqrt_prec, sin_prec, exp_prec implicit none @@ -1657,7 +1645,7 @@ subroutine tripping(tb,ta) call random_number(randx) h_coeff(j)=one*(randx-zpfive) enddo - h_coeff=h_coeff/sqrt_prec(real(z_modes,mytype)) + h_coeff=h_coeff/sqrt(real(z_modes,mytype)) endif !Initialization h_nxt (always bounded by xsize(3)^2 operations) @@ -1668,7 +1656,7 @@ subroutine tripping(tb,ta) h_nxt(k)=zero z_pos=-zlz*zpfive+(xstart(3)+(k-1)-1)*dz do j=1,z_modes - h_nxt(k)= h_nxt(k)+h_coeff(j)*sin_prec(two*pi*j*z_pos/zlz) + h_nxt(k)= h_nxt(k)+h_coeff(j)*sin(two*pi*j*z_pos/zlz) enddo enddo end if @@ -1687,7 +1675,7 @@ subroutine tripping(tb,ta) call random_number(randx) h_coeff(j)=one*(randx-zpfive) enddo - h_coeff=h_coeff/sqrt_prec(real(z_modes,mytype)) !Non-dimensionalization + h_coeff=h_coeff/sqrt(real(z_modes,mytype)) !Non-dimensionalization end if call MPI_BCAST(h_coeff,z_modes,real_type,0,MPI_COMM_WORLD,code) @@ -1698,7 +1686,7 @@ subroutine tripping(tb,ta) h_nxt(k)=zero z_pos=-zlz*zpfive+(xstart(3)+(k-1)-1)*dz do j=1,z_modes - h_nxt(k)= h_nxt(k)+h_coeff(j)*sin_prec(two*pi*j*z_pos/zlz) + h_nxt(k)= h_nxt(k)+h_coeff(j)*sin(two*pi*j*z_pos/zlz) enddo enddo endif @@ -1716,8 +1704,8 @@ subroutine tripping(tb,ta) do k=1,xsize(3) !g(z)*EXP_F(X,Y) ta(i,j,k)=((one-b_tr)*h_i(k)+b_tr*h_nxt(k)) - !ta(i,j,k)=A_tr*exp_prec(-((x_pos-x0_tr)/xs_tr)**2-(y_pos/ys_tr)**2)*ta(i,j,k) - ta(i,j,k)=A_tr*exp_prec(-((x_pos-x0_tr)/xs_tr)**2-((y_pos-zpfive)/ys_tr)**2)*ta(i,j,k) + !ta(i,j,k)=A_tr*exp(-((x_pos-x0_tr)/xs_tr)**2-(y_pos/ys_tr)**2)*ta(i,j,k) + ta(i,j,k)=A_tr*exp(-((x_pos-x0_tr)/xs_tr)**2-((y_pos-zpfive)/ys_tr)**2)*ta(i,j,k) tb(i,j,k)=tb(i,j,k)+ta(i,j,k) z_pos=-zlz*zpfive+(xstart(3)+(k-1)-1)*dz @@ -1745,7 +1733,6 @@ subroutine tbl_tripping(tb,ta) use decomp_2d_constants use decomp_2d_mpi use mpi - use dbg_schemes, only: sqrt_prec, exp_prec, sin_prec implicit none @@ -1782,11 +1769,11 @@ subroutine tbl_tripping(tb,ta) nxt_itr=1 do j=1,z_modes call random_number(randx) - h_coeff1(j)=one*(randx-zpfive)/sqrt_prec(real(z_modes,mytype)) + h_coeff1(j)=one*(randx-zpfive)/sqrt(real(z_modes,mytype)) call random_number(randx) phase1(j) = two*pi*randx call random_number(randx) - h_coeff2(j)=one*(randx-zpfive)/sqrt_prec(real(z_modes,mytype)) + h_coeff2(j)=one*(randx-zpfive)/sqrt(real(z_modes,mytype)) call random_number(randx) phase2(j) = two*pi*randx enddo @@ -1806,8 +1793,8 @@ subroutine tbl_tripping(tb,ta) h_2(k)=zero z_pos=-zlz*zpfive+real(xstart(3)+(k-1)-1,mytype)*dz do j=1,z_modes - h_1(k)= h_1(k)+h_coeff1(j)*sin_prec(two*pi*real(j,mytype)*z_pos/zlz+phase1(j)) - h_2(k)= h_2(k)+h_coeff2(j)*sin_prec(two*pi*real(j,mytype)*z_pos/zlz+phase2(j)) + h_1(k)= h_1(k)+h_coeff1(j)*sin(two*pi*real(j,mytype)*z_pos/zlz+phase1(j)) + h_2(k)= h_2(k)+h_coeff2(j)*sin(two*pi*real(j,mytype)*z_pos/zlz+phase2(j)) enddo enddo end if @@ -1823,7 +1810,7 @@ subroutine tbl_tripping(tb,ta) if (nrank == 0) then do j=1,z_modes call random_number(randx) - h_coeff1(j)=one*(randx-zpfive)/sqrt_prec(real(z_modes,mytype)) + h_coeff1(j)=one*(randx-zpfive)/sqrt(real(z_modes,mytype)) call random_number(randx) phase1(j) = two*pi*randx enddo @@ -1837,7 +1824,7 @@ subroutine tbl_tripping(tb,ta) h_1(k)=zero z_pos=-zlz*zpfive+real(xstart(3)+(k-1)-1,mytype)*dz do j=1,z_modes - h_1(k)= h_1(k)+h_coeff1(j)*sin_prec(two*pi*real(j,mytype)*z_pos/zlz+phase1(j)) + h_1(k)= h_1(k)+h_coeff1(j)*sin(two*pi*real(j,mytype)*z_pos/zlz+phase1(j)) enddo enddo endif @@ -1853,7 +1840,7 @@ subroutine tbl_tripping(tb,ta) y_pos=yp(xstart(2)+(j-1)) do k=1,xsize(3) ta(i,j,k)=((one-b_tr)*h_1(k)+b_tr*h_2(k)) - ta(i,j,k)=A_tr*exp_prec(-((x_pos-x0_tr_tbl)/xs_tr_tbl)**2-((y_pos-0.05_mytype)/ys_tr_tbl)**2)*ta(i,j,k) + ta(i,j,k)=A_tr*exp(-((x_pos-x0_tr_tbl)/xs_tr_tbl)**2-((y_pos-0.05_mytype)/ys_tr_tbl)**2)*ta(i,j,k) tb(i,j,k)=tb(i,j,k)+ta(i,j,k) z_pos=-zlz*zpfive+real(xstart(3)+(k-1)-1,mytype)*dz @@ -1869,48 +1856,6 @@ subroutine tbl_tripping(tb,ta) end subroutine tbl_tripping !################################################################## !################################################################## -function rl(complexnumber) - - use param - - implicit none - - real(mytype) :: rl - complex(mytype) :: complexnumber - - rl = real(complexnumber, kind=mytype) - -end function rl -!################################################################## -!################################################################## -function iy(complexnumber) - - use param - - implicit none - - real(mytype) :: iy - complex(mytype) :: complexnumber - - iy = aimag(complexnumber) - -end function iy -!################################################################## -!################################################################## -function cx(realpart,imaginarypart) - - use param - - implicit none - - complex(mytype) :: cx - real(mytype) :: realpart, imaginarypart - - cx = cmplx(realpart, imaginarypart, kind=mytype) - -end function cx -!################################################################## -!################################################################## subroutine calc_temp_eos(temp, rho, phi, mweight, xlen, ylen, zlen) use decomp_2d diff --git a/src/utilities.f90 b/src/utilities.f90 new file mode 100644 index 000000000..d2bee80d4 --- /dev/null +++ b/src/utilities.f90 @@ -0,0 +1,54 @@ +!Copyright (c) 2012-2022, Xcompact3d +!This file is part of Xcompact3d (xcompact3d.com) +!SPDX-License-Identifier: BSD 3-Clause + +module utilities + + use decomp_2d_constants, only : mytype + + implicit none + + private + + public :: rl, iy, cx + +contains + + !################################################################## + function rl(complexnumber) + + implicit none + + real(mytype) :: rl + complex(mytype) :: complexnumber + + rl = real(complexnumber, kind=mytype) + +end function rl +!################################################################## +!################################################################## +function iy(complexnumber) + + implicit none + + real(mytype) :: iy + complex(mytype) :: complexnumber + + iy = aimag(complexnumber) + +end function iy +!################################################################## +!################################################################## +function cx(realpart,imaginarypart) + + implicit none + + complex(mytype) :: cx + real(mytype) :: realpart, imaginarypart + + cx = cmplx(realpart, imaginarypart, kind=mytype) + +end function cx +!################################################################## + +end module utilities