Merge pull request #49 from chowland/pf-tweaks

A few minor fixes for the phase-field model: mainly index issues
chowland · May 3, 2024 · f538996 · f538996
2 parents 02fff05 + abb1b7c
commit f538996
Showing 1 changed file with 14 additions and 14 deletions.
diff --git a/src/phasefield.F90 b/src/phasefield.F90
@@ -507,8 +507,8 @@ subroutine InterpTempMultigrid
     ! Fill temporary array with temperature field and BCs
     do ic=xstart(3)-lvlhalo,xend(3)+lvlhalo
         do jc=xstart(2)-lvlhalo,xend(2)+lvlhalo
-            tpdv(0,jc,ic) = 2.0*tempbp(1,jc,ic) - temp(1,jc,ic)
-            tpdv(nx,jc,ic) = 2.0*temptp(1,jc,ic) - temp(nxm,jc,ic)
+            tpdv(0,jc,ic) = (1.0 - 2.0*TfixS)*temp(1,jc,ic) + 2.0*TfixS*tempbp(1,jc,ic)
+            tpdv(nx,jc,ic) = (1.0 - 2.0*TfixN)*temp(nxm,jc,ic) + 2.0*TfixN*temptp(1,jc,ic)
             do kc=1,nxm
                 tpdv(kc,jc,ic) = temp(kc,jc,ic)
             end do
@@ -689,12 +689,12 @@ subroutine AddLatentHeat
                     kcr = irangs(kc)
 
                     qv3 = tpdvr(kcr-2:kcr+1,jcr-2:jcr+1,icr-2:icr+1)
-                    qv2(:,:) = qv3(:,:,1)*czsalc(1,ic) + qv3(:,:,2)*czsalc(2,ic) &
-                                + qv3(:,:,3)*czsalc(3,ic) + qv3(:,:,4)*czsalc(4,ic)
-                    qv1(:) = qv2(:,1)*cysalc(1,jc) + qv2(:,2)*cysalc(2,jc) &
-                            + qv2(:,3)*cysalc(3,jc) + qv2(:,4)*cysalc(4,jc)
+                    qv2(:,:) = qv3(:,:,1)*czphic(1,ic) + qv3(:,:,2)*czphic(2,ic) &
+                                + qv3(:,:,3)*czphic(3,ic) + qv3(:,:,4)*czphic(4,ic)
+                    qv1(:) = qv2(:,1)*cyphic(1,jc) + qv2(:,2)*cyphic(2,jc) &
+                            + qv2(:,3)*cyphic(3,jc) + qv2(:,4)*cyphic(4,jc)
 
-                    phi_rhs = sum(qv1(1:4)*cxsalc(1:4,kc))
+                    phi_rhs = sum(qv1(1:4)*cxphic(1:4,kc))
                     hro(kc,jc,ic) = hro(kc,jc,ic) + pf_S*phi_rhs*aldt
                 end do
             end do
@@ -706,13 +706,13 @@ subroutine AddLatentHeat
 
                     qv3 = tpdvr(kcr-1:kcr+2,jcr-1:jcr+2,icr-1:icr+2)
                     do ic=max(krangr(icr),xstart(3)),min(krangr(icr+1)-1,xend(3))
-                        qv2(:,:) = qv3(:,:,1)*czsalc(1,ic) + qv3(:,:,2)*czsalc(2,ic)&
-                                +qv3(:,:,3)*czsalc(3,ic) + qv3(:,:,4)*czsalc(4,ic)
+                        qv2(:,:) = qv3(:,:,1)*czphic(1,ic) + qv3(:,:,2)*czphic(2,ic)&
+                                +qv3(:,:,3)*czphic(3,ic) + qv3(:,:,4)*czphic(4,ic)
                         do jc=max(jrangr(jcr),xstart(2)),min(jrangr(jcr+1)-1,xend(2))
-                            qv1(:) = qv2(:,1)*cysalc(1,jc) + qv2(:,2)*cysalc(2,jc) &
-                                    +qv2(:,3)*cysalc(3,jc) + qv2(:,4)*cysalc(4,jc)
+                            qv1(:) = qv2(:,1)*cyphic(1,jc) + qv2(:,2)*cyphic(2,jc) &
+                                    +qv2(:,3)*cyphic(3,jc) + qv2(:,4)*cyphic(4,jc)
                             do kc=max(irangr(kcr),1),min(irangr(kcr+1)-1,nxm)
-                                phi_rhs = sum(qv1(1:4)*cxsalc(1:4,kc))
+                                phi_rhs = sum(qv1(1:4)*cxphic(1:4,kc))
 
                                 hro(kc,jc,ic) = hro(kc,jc,ic) + pf_S*phi_rhs*aldt
                             end do
@@ -898,9 +898,9 @@ subroutine SolveImpEqnUpdate_YZ_pf(q, rhs, axis)
 
             do kc=1,nxm
                 if (axis=="y") then
-                    philoc = 0.5*(phic(kc,jc,ic) + phic(kc,jc+1,ic))
+                    philoc = 0.5*(phic(kc,jc,ic) + phic(kc,jc-1,ic))
                 elseif (axis=="z") then
-                    philoc = 0.5*(phic(kc,jc,ic) + phic(kc,jc,ic+1))
+                    philoc = 0.5*(phic(kc,jc,ic) + phic(kc,jc,ic-1))
                 else
                     philoc = phic(kc,jc,ic)
                 end if