Merge pull request #267 from xcompact3d/remove_dbg_schemes

Remove dbg schemes

src/CMakeLists.txt

@@ -49,6 +49,7 @@ add_executable(xcompact3d
                     tools.f90
                     transeq.f90
                     turbine.f90
+		    utilities.f90
                     variables.f90
                     visu.f90
                     xcompact3d.f90)

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

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
 

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

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