Add production limiter in sst model and fix openmp issue in jacobi pr…

…econditioner

lib/Condiff.f03

@@ -12,9 +12,11 @@ Subroutine Condiff(scalar)
 real(8) St(Ic,Jc),Sm(Ic,Jc),fw1(Ic,Jc),alphastar(Ic,Jc),betastar(Ic,Jc),alpha(Ic,Jc),beta(Ic,Jc),Dwt(Ic,Jc),C3e(Ic,Jc)
 character(*) scalar
 character(6) wallfunutype,wallfunktype
+logical(1) productlimit
 
 wallfunutype='parvel'
 wallfunktype='genlaw'
+productlimit=.true.
 
 !$OMP PARALLEL
 if(scalar=='U'.or.scalar=='V') then
@@ -374,13 +376,19 @@ Subroutine Condiff(scalar)
        end if
       end if
       !b(i,j)=min(b(i,j),10*Dampk)
-     else if(Walltreat=='lr') then   
-      b(i,j)=min(mut(i,j)*St(i,j)**2*Jg(i,j)*dx*dy,10*Dampk)-Dampk
-      !b(i,j)=mut(i,j)*St(i,j)**2*Jg(i,j)*dx*dy-Dampk
+     else if(Walltreat=='lr') then
+      if(productlimit) then
+       b(i,j)=min(mut(i,j)*St(i,j)**2*Jg(i,j)*dx*dy,10*Dampk)-Dampk
+      else
+       b(i,j)=mut(i,j)*St(i,j)**2*Jg(i,j)*dx*dy-Dampk
+      end if
      end if
     else
-     b(i,j)=min(mut(i,j)*St(i,j)**2*Jg(i,j)*dx*dy,10*Dampk)-Dampk
-     !b(i,j)=mut(i,j)*St(i,j)**2*Jg(i,j)*dx*dy-Dampk
+     if(productlimit) then
+      b(i,j)=min(mut(i,j)*St(i,j)**2*Jg(i,j)*dx*dy,10*Dampk)-Dampk
+     else
+      b(i,j)=mut(i,j)*St(i,j)**2*Jg(i,j)*dx*dy-Dampk
+     end if
     end if
    else if(scalar=='Te') then
     if(j==1.and.(i>=Ib1.and.i<=Ib2)) then

lib/Linearsolver.f03

@@ -134,12 +134,16 @@ Subroutine CGSTAB(aM,b,F,F0,a,Ic,Jc,Ib1,Ib2,scalar)
  beta=(rho*alpha)/(rho0*omiga)
  !$OMP WORKSHARE
  p=rms+beta*(p-omiga*v)
+ !$OMP END WORKSHARE
  if(pretype=='JAC') then
+  !$OMP WORKSHARE
   y=a*p/aM(1,:,:)
+  !$OMP END WORKSHARE
  else
+  !$OMP WORKSHARE
   y=p
+  !$OMP END WORKSHARE
  end if
- !$OMP END WORKSHARE
  if(pretype=='ILU'.or.pretype=='SSOR') then
   Call Sorprecond(aM,aD,y,a,Ic,Jc,Ib1,Ib2,scalar,pretype)
  end if
@@ -167,12 +171,16 @@ Subroutine CGSTAB(aM,b,F,F0,a,Ic,Jc,Ib1,Ib2,scalar)
  alpha=rho/sumvrms
  !$OMP WORKSHARE
  s=rms-alpha*v
+ !$OMP END WORKSHARE
  if(pretype=='JAC') then
+  !$OMP WORKSHARE
   z=a*s/aM(1,:,:)
+  !$OMP END WORKSHARE
  else
+  !$OMP WORKSHARE
   z=s
+  !$OMP END WORKSHARE
  end if
- !$OMP END WORKSHARE
  if(pretype=='ILU'.or.pretype=='SSOR') then
   Call Sorprecond(aM,aD,z,a,Ic,Jc,Ib1,Ib2,scalar,pretype)
  end if
@@ -194,13 +202,15 @@ Subroutine CGSTAB(aM,b,F,F0,a,Ic,Jc,Ib1,Ib2,scalar)
    end DO
  end DO
  !$OMP END DO
- !$OMP WORKSHARE
  if(pretype=='JAC') then
+  !$OMP WORKSHARE
   pt=a*t/aM(1,:,:)
+  !$OMP END WORKSHARE
  else
+  !$OMP WORKSHARE
   pt=t
+  !$OMP END WORKSHARE
  end if
- !$OMP END WORKSHARE
  if(pretype=='ILU'.or.pretype=='SSOR') then
   Call Sorprecond(aM,aD,pt,a,Ic,Jc,Ib1,Ib2,scalar,pretype)
  end if