diff --git a/Pressure.f90 b/Pressure.f90
index 853e24f..fa5d916 100644
--- a/Pressure.f90
+++ b/Pressure.f90
@@ -12,8 +12,6 @@ Subroutine Pressure(HydroParam,MeshParam,dt)
     type(MeshGridParam) :: MeshParam
     type(HydrodynamicParam) :: HydroParam
 
-    !HydroParam%iBarot = 0
-    HydroParam%PBarc = 0.d0
     ! 1. Define the Cell centered Horizontal Eddy-Viscosity:
     If (HydroParam%iBarot == 1) Then         ! User-defined Horizontal Diffusion
         !!$OMP parallel do default(none) shared(MeshParam,HydroParam) private(iLayer,iEdge,l,r,iLayer_bar,NearZero,dt)
@@ -30,24 +28,20 @@ Subroutine Pressure(HydroParam,MeshParam,dt)
             
             Do iLayer = HydroParam%Smallm(iEdge), HydroParam%CapitalM(iEdge)
                 If ( r == 0 .or. HydroParam%H(iEdge) <= HydroParam%PCRI+NearZero) Then
-                    HydroParam%PBarc(iLayer,iEdge) = 0.d0
+                    HydroParam%PBarc(iLayer,iEdge) = 0.
                 Else
-                    HydroParam%PBarc(iLayer,iEdge) = 0.d0
-                    If (iLayer>=HydroParam%ElSmallm(l).and.iLayer>=HydroParam%ElSmallm(r)) Then 
+                    HydroParam%PBarc(iLayer,iEdge) = 0.
+                    If (iLayer>=HydroParam%ElSmallm(l).and.iLayer>=HydroParam%ElSmallm(r)) Then
                         Do iLayer_bar=iLayer,HydroParam%CapitalM(iEdge)
                             If (iLayer_bar==iLayer) Then 
-                                HydroParam%PBarc(iLayer,iEdge) = HydroParam%PBarc(iLayer,iEdge) + HydroParam%g*dt/(MeshParam%CirDistance(iEdge)*HydroParam%rho0)*0.5d0*((HydroParam%Theta*HydroParam%sDRhoW(iLayer_bar,r)+(1.-HydroParam%Theta)*HydroParam%sDRhoWt(iLayer_bar,r))-(HydroParam%Theta*HydroParam%sDRhoW(iLayer_bar,l)+(1.-HydroParam%Theta)*HydroParam%sDRhoWt(iLayer_bar,l)))!(HydroParam%sDRhoW(iLayer_bar,r)-HydroParam%sDRhoW(iLayer_bar,l))*HydroParam%DZj(iLayer_bar,iEdge) !HydroParam%g*dt/(MeshParam%CirDistance(iEdge)*HydroParam%rho0)*0.5*(HydroParam%sDRhoW(iLayer_bar,r)-HydroParam%sDRhoW(iLayer_bar,l))*HydroParam%DZj(iLayer_bar,iEdge)
+                                HydroParam%PBarc(iLayer,iEdge) = HydroParam%PBarc(iLayer,iEdge) + HydroParam%g*dt/(MeshParam%CirDistance(iEdge)*HydroParam%rho0)*0.5*(HydroParam%sDRhoW(iLayer_bar,r)-HydroParam%sDRhoW(iLayer_bar,l))*HydroParam%DZj(iLayer_bar,iEdge) !HydroParam%g*dt/(MeshParam%CirDistance(iEdge)*HydroParam%rho0)*0.5*(HydroParam%sDRhoW(iLayer_bar,r)-HydroParam%sDRhoW(iLayer_bar,l))*HydroParam%DZj(iLayer_bar,iEdge)
                             Else
-                                HydroParam%PBarc(iLayer,iEdge) = HydroParam%PBarc(iLayer,iEdge) + HydroParam%g*dt/(MeshParam%CirDistance(iEdge)*HydroParam%rho0)*1.d0*((HydroParam%Theta*HydroParam%sDRhoW(iLayer_bar,r)+(1.-HydroParam%Theta)*HydroParam%sDRhoWt(iLayer_bar,r))-(HydroParam%Theta*HydroParam%sDRhoW(iLayer_bar,l)+(1.-HydroParam%Theta)*HydroParam%sDRhoWt(iLayer_bar,l)))!(HydroParam%sDRhoW(iLayer_bar,r)-HydroParam%sDRhoW(iLayer_bar,l))*HydroParam%DZj(iLayer_bar,iEdge) !HydroParam%g*dt/(MeshParam%CirDistance(iEdge)*HydroParam%rho0)*1.0*(HydroParam%sDRhoW(iLayer_bar,r)-HydroParam%sDRhoW(iLayer_bar,l))*HydroParam%DZj(iLayer_bar,iEdge)
+                                HydroParam%PBarc(iLayer,iEdge) = HydroParam%PBarc(iLayer,iEdge) + HydroParam%g*dt/(MeshParam%CirDistance(iEdge)*HydroParam%rho0)*1.*(HydroParam%sDRhoW(iLayer_bar,r)-HydroParam%sDRhoW(iLayer_bar,l))*HydroParam%DZj(iLayer_bar,iEdge) !HydroParam%g*dt/(MeshParam%CirDistance(iEdge)*HydroParam%rho0)*1.0*(HydroParam%sDRhoW(iLayer_bar,r)-HydroParam%sDRhoW(iLayer_bar,l))*HydroParam%DZj(iLayer_bar,iEdge)
                             EndIf
                         EndDo
                     EndIf
                 EndIf
                 
-                !if (abs( HydroParam%PBarc(iLayer,iEdge))<1.e-8) then
-                !     HydroParam%PBarc(iLayer,iEdge)=0.d0
-                !endif
-                
                 HydroParam%Fu(iLayer,iEdge) = HydroParam%Fu(iLayer,iEdge) - HydroParam%PBarc(iLayer,iEdge)
             EndDo       ! Layer Loop
         EndDo       ! Edge Loop
diff --git a/UpWindCWC_HR_GATS.F90 b/UpWindCWC_HR_GATS.F90
index 37d3f2a..0c52427 100644
--- a/UpWindCWC_HR_GATS.F90
+++ b/UpWindCWC_HR_GATS.F90
@@ -289,7 +289,7 @@ Subroutine UpWindCWC_HR_GATS(index,uLoadVar,dVar,VarEst,VarEstP,dt,dtday,HydroPa
                 EndIf
                 
                 If (HydroParam%Theta*HydroParam%w(iLayer+1,iElem)+(1.-HydroParam%Theta)*HydroParam%wt(iLayer+1,iElem)>0) Then !top layer outflow 
-                    !if (iLayer < HydroParam%ElCapitalM(iElem)) then
+                    if (iLayer < HydroParam%ElCapitalM(iElem)) then
                         !set boundary condition for psi_cel top layer outflow
                         b = VarEstP(iLayer,iElem)
                         bd = HydroParam%DZit(iLayer,iElem)
@@ -313,11 +313,11 @@ Subroutine UpWindCWC_HR_GATS(index,uLoadVar,dVar,VarEst,VarEstP,dt,dtday,HydroPa
                         dendist = (0.5d0*(bd+cd))  !(0.5d0*(HydroParam%DZit(iLayer,iElem)+HydroParam%DZit(iLayer-1,iElem))+NearZero))  
                     
                         if (abs(a-b)<=NearZero.or.abs(b-c)<=NearZero) then
-                            HydroParam%psi_cell(iLayer,iElem) = 0
+                            HydroParam%psi_cell(iLayer+1,iElem) = 0
                             cycle
                         endif
                         if (abs(numdist)<=NearZero.or.abs(dendist)<=NearZero) then
-                            HydroParam%psi_cell(iLayer,iElem) = 0
+                            HydroParam%psi_cell(iLayer+1,iElem) = 0
                             cycle
                         endif
                     
@@ -326,9 +326,9 @@ Subroutine UpWindCWC_HR_GATS(index,uLoadVar,dVar,VarEst,VarEstP,dt,dtday,HydroPa
                         r_face =  ((VarEstP(iLayer,iElem)-VarEstP(iLayer-1,iElem))/(VarEstP(iLayer+1,iElem)-VarEstP(iLayer,iElem)+NearZero))*(0.5d0*(HydroParam%DZit(iLayer,iElem)+HydroParam%DZit(iLayer+1,iElem))/(0.5d0*(HydroParam%DZit(iLayer,iElem)+HydroParam%DZit(iLayer-1,iElem))+NearZero)) !r_face = (r_num_plus(iLayer)/(abs(VarEstP(iLayer+1,iElem)-VarEstP(iLayer,iElem))*(r_den_plus(iLayer)) + NearZero))
                         Courant = 0.5*(HydroParam%Theta*HydroParam%u(iLayer,Face)*+(1.-HydroParam%Theta)*HydroParam%ut(iLayer,Face))*dt/MeshParam%CirDistance(Face)
                         Call Psi_value(psi_flag,r_face,Courant,fi_small,psi_face)
-                        HydroParam%psi_cell(iLayer,iElem) = psi_face
+                        HydroParam%psi_cell(iLayer+1,iElem) = psi_face
                     
-                        !endif
+                        endif
                     EndIf  
                 
                 !Do iEdge = 1,4