Enforce BC for the normal at walls

src/bound.f90

@@ -80,7 +80,7 @@ end subroutine bounduvw
   !
   subroutine boundp(cbc,n,bc,nb,is_bound,dl,dzc,p)
     !
-    ! imposes pressure boundary conditions
+    ! imposes pressure and scalar boundary conditions
     !
     implicit none
     character(len=1), intent(in), dimension(0:1,3) :: cbc

src/main.f90

@@ -49,6 +49,7 @@ program cans
   use mod_output         , only: out0d,gen_alias,out1d,out1d_chan,out2d,out3d,write_log_output,write_visu_2d,write_visu_3d
   use mod_param          , only: rkcoeff,small, &
                                  nb,is_bound,cbcvel,bcvel,cbcpre,bcpre,cbcsca,bcsca,cbcpsi,bcpsi, &
+                                 cbcnox,bcnox,cbcnoy,bcnoy,cbcnoz,bcnoz, &
                                  icheck,iout0d,iout1d,iout2d,iout3d,isave, &
                                  nstep,time_max,tw_max,stop_type,restart,is_overwrite_save,nsaves_max, &
                                  datadir,   &
@@ -327,9 +328,9 @@ program cans
   call acdi_cmpt_phi(n,seps,psi,phi)
   call cmpt_norm_curv(n,dli,dzci,dzfi,phi,normx,normy,normz,kappa)
   call boundp(cbcpsi,n,bcpsi,nb,is_bound,dl,dzc,kappa)
-  call boundp(cbcpsi,n,bcpsi,nb,is_bound,dl,dzc,normx)
-  call boundp(cbcpsi,n,bcpsi,nb,is_bound,dl,dzc,normy)
-  call boundp(cbcpsi,n,bcpsi,nb,is_bound,dl,dzc,normz)
+  call boundp(cbcnox,n,bcnox,nb,is_bound,dl,dzc,normx)
+  call boundp(cbcnoy,n,bcnoy,nb,is_bound,dl,dzc,normy)
+  call boundp(cbcnoz,n,bcnoz,nb,is_bound,dl,dzc,normz)
   !
 #if defined(_CONSERVATIVE_MOMENTUM)
   !$acc kernels default(present) async(1)
@@ -391,9 +392,9 @@ program cans
         call acdi_cmpt_phi(n,seps,psi,phi)
         call cmpt_norm_curv(n,dli,dzci,dzfi,phi,normx,normy,normz,kappa)
         call boundp(cbcpsi,n,bcpre,nb,is_bound,dl,dzc,kappa)
-        call boundp(cbcpsi,n,bcpsi,nb,is_bound,dl,dzc,normx)
-        call boundp(cbcpsi,n,bcpsi,nb,is_bound,dl,dzc,normy)
-        call boundp(cbcpsi,n,bcpsi,nb,is_bound,dl,dzc,normz)
+        call boundp(cbcnox,n,bcnox,nb,is_bound,dl,dzc,normx)
+        call boundp(cbcnoy,n,bcnoy,nb,is_bound,dl,dzc,normy)
+        call boundp(cbcnoz,n,bcnoz,nb,is_bound,dl,dzc,normz)
       end if
 #if defined(_SCALAR)
       call tm_scal(tm_coeff,n,dli,dzci,dzfi,dt,ssource,rho12,ka12,cp12,psi,u,v,w,s)

src/param.f90

@@ -55,6 +55,12 @@ module mod_param
 real(rp)        , protected, dimension(0:1,3)   ::   bcpsi
 character(len=1), protected, dimension(0:1,3)   ::  cbcsca
 real(rp)        , protected, dimension(0:1,3)   ::   bcsca
+character(len=1), protected, dimension(0:1,3)   ::  cbcnox
+real(rp)        , protected, dimension(0:1,3)   ::   bcnox
+character(len=1), protected, dimension(0:1,3)   ::  cbcnoy
+real(rp)        , protected, dimension(0:1,3)   ::   bcnoy
+character(len=1), protected, dimension(0:1,3)   ::  cbcnoz
+real(rp)        , protected, dimension(0:1,3)   ::   bcnoz
 !
 real(rp), protected, dimension(3) :: bforce,gacc
 real(rp), protected               :: ssource
@@ -167,6 +173,44 @@ subroutine read_input(myid)
 #if defined(_CONSTANT_COEFFS_POISSON)
     rho0 = minval(rho12(:))
 #endif
+    !
+    ! set normals BC based on cbcpsi
+    !
+    block
+      integer :: b,c
+      logical :: is_wall
+      cbcnox = cbcpsi
+      cbcnoy = cbcpsi
+      cbcnoz = cbcpsi
+      bcnox  = bcpsi
+      bcnoy  = bcpsi
+      bcnoz  = bcpsi
+      do b = 0,1
+        do c = 1,3
+          is_wall = cbcpsi(b,c)=='N'
+          select case(c)
+          case(1)
+            if (is_wall) then
+              cbcnox(b,1) = 'D'
+              cbcnoy(b,1) = 'I'
+              cbcnoz(b,1) = 'I'
+            end if
+          case(2)
+            if (is_wall) then
+              cbcnox(b,2) = 'I'
+              cbcnoy(b,2) = 'D'
+              cbcnoz(b,2) = 'I'
+            end if
+          case(3)
+            if (is_wall) then
+              cbcnox(b,3) = 'I'
+              cbcnoy(b,3) = 'I'
+              cbcnoz(b,3) = 'D'
+            end if
+          end select
+        end do
+      end do
+    end block
 #if defined(_OPENACC)
     !
     ! read cuDecomp parameter file cudecomp.in, if it exists