diff --git a/Code/Source/svFSI/DISTRIBUTE.f b/Code/Source/svFSI/DISTRIBUTE.f
index 756c9c72..3e3a6de7 100755
--- a/Code/Source/svFSI/DISTRIBUTE.f
+++ b/Code/Source/svFSI/DISTRIBUTE.f
@@ -186,6 +186,10 @@ SUBROUTINE DISTRIBUTE
          CALL cm%bcast(cmmVarWall)
          CALL cm%bcast(shlEq)
          CALL cm%bcast(pstEq)
+!        Varwall properties---------------------------------------------
+         CALL cm%bcast(useVarWall)
+         CALL cm%bcast(nvwp)
+!        ---------------------------------------------------------------
          CALL cm%bcast(sstEq)
          CALL cm%bcast(cepEq)
          IF (rmsh%isReqd) THEN
@@ -256,6 +260,24 @@ SUBROUTINE DISTRIBUTE
          DEALLOCATE(tmpX)
       END IF
 
+!     Varwall properties------------------------------------------------
+!     Distribute variable wall properties (vWP0) to processors
+      flag = ALLOCATED(vWP0)
+      CALL cm%bcast(flag)
+      IF (flag) THEN
+         IF (cm%mas()) THEN
+            ALLOCATE(tmpX(nvwp,gtnNo))
+            tmpX = vWP0
+            DEALLOCATE(vWP0)
+         ELSE
+            ALLOCATE(tmpX(0,0))
+         END IF
+         ALLOCATE(vWP0(nvwp,tnNo))
+         vWP0 = LOCAL(tmpX)
+         DEALLOCATE(tmpX)
+      END IF
+!     ------------------------------------------------------------------
+
 !     Distribute initial flow quantities to processors
       flag = ALLOCATED(Pinit)
       CALL cm%bcast(flag)
diff --git a/Code/Source/svFSI/FSI.f b/Code/Source/svFSI/FSI.f
index 234e4ad8..56a88f8b 100755
--- a/Code/Source/svFSI/FSI.f
+++ b/Code/Source/svFSI/FSI.f
@@ -52,7 +52,7 @@ SUBROUTINE CONSTRUCT_FSI(lM, Ag, Yg, Dg)
       INTEGER(KIND=IKIND), ALLOCATABLE :: ptr(:)
       REAL(KIND=RKIND), ALLOCATABLE :: xl(:,:), al(:,:), yl(:,:),
      2   dl(:,:), bfl(:,:), fN(:,:), pS0l(:,:), pSl(:), ya_l(:),
-     3   lR(:,:), lK(:,:,:), lKd(:,:,:)
+     3   lR(:,:), lK(:,:,:), lKd(:,:,:), lVWP(:,:)
       REAL(KIND=RKIND), ALLOCATABLE :: xwl(:,:), xql(:,:), Nwx(:,:),
      2   Nwxx(:,:), Nqx(:,:)
 
@@ -72,7 +72,7 @@ SUBROUTINE CONSTRUCT_FSI(lM, Ag, Yg, Dg)
       ALLOCATE(ptr(eNoN), xl(nsd,eNoN), al(tDof,eNoN), yl(tDof,eNoN),
      2   dl(tDof,eNoN), bfl(nsd,eNoN), fN(nsd,nFn), pS0l(nsymd,eNoN),
      3   pSl(nsymd), ya_l(eNoN), lR(dof,eNoN), lK(dof*dof,eNoN,eNoN),
-     4   lKd(dof*nsd,eNoN,eNoN))
+     4   lKd(dof*nsd,eNoN,eNoN), lVWP(nvwp,eNoN))
 
 !     Loop over all elements of mesh
       DO e=1, lM%nEl
@@ -98,6 +98,12 @@ SUBROUTINE CONSTRUCT_FSI(lM, Ag, Yg, Dg)
             yl(:,a)  = Yg(:,Ac)
             dl(:,a)  = Dg(:,Ac)
             bfl(:,a) = Bf(:,Ac)
+
+!           Varwall properties -----------------------------------------
+!           Calculate local wall property
+            IF (useVarWall) lVWP(:,a) = vWP0(:,Ac)
+!           ------------------------------------------------------------
+
             IF (ALLOCATED(lM%fN)) THEN
                DO iFn=1, nFn
                   fN(:,iFn) = lM%fN((iFn-1)*nsd+1:iFn*nsd,e)
@@ -159,11 +165,11 @@ SUBROUTINE CONSTRUCT_FSI(lM, Ag, Yg, Dg)
 
                CASE (phys_lElas)
                   CALL LELAS3D(fs(1)%eNoN, w, fs(1)%N(:,g), Nwx, al, dl,
-     2               bfl, pS0l, pSl, lR, lK)
+     2               bfl, pS0l, pSl, lR, lK, lVWP)
 
                CASE (phys_struct)
                   CALL STRUCT3D(fs(1)%eNoN, nFn, w, fs(1)%N(:,g), Nwx,
-     2               al, yl, dl, bfl, fN, pS0l, pSl, ya_l, lR, lK)
+     2               al, yl, dl, bfl, fN, pS0l, pSl, ya_l, lR, lK, lVWP)
 
                CASE (phys_ustruct)
                   CALL USTRUCT3D_M(vmsStab, fs(1)%eNoN, fs(2)%eNoN, nFn,
diff --git a/Code/Source/svFSI/INITIALIZE.f b/Code/Source/svFSI/INITIALIZE.f
index 7cfe05c2..9379cfe3 100755
--- a/Code/Source/svFSI/INITIALIZE.f
+++ b/Code/Source/svFSI/INITIALIZE.f
@@ -154,6 +154,7 @@ SUBROUTINE INITIALIZE(timeP)
       IF (dFlag) i = 3*tDof
       IF (pstEq) i = i + nsymd
       IF (sstEq) i = i + nsd
+!      IF (useVarWall) i = i + nvwp
       IF (cepEq) THEN
          i = i + nXion
          IF (cem%cpld) i = i + 1
@@ -663,6 +664,10 @@ SUBROUTINE FINALIZE
       IF (ALLOCATED(pSn))      DEALLOCATE(pSn)
       IF (ALLOCATED(pSa))      DEALLOCATE(pSa)
 
+!     Varwall properties -----------------------------------------------
+      IF (ALLOCATED(vWP0))      DEALLOCATE(vWP0)
+!     ------------------------------------------------------------------
+
       IF (ALLOCATED(Pinit))    DEALLOCATE(Pinit)
       IF (ALLOCATED(Vinit))    DEALLOCATE(Vinit)
       IF (ALLOCATED(Dinit))    DEALLOCATE(Dinit)
diff --git a/Code/Source/svFSI/LELAS.f b/Code/Source/svFSI/LELAS.f
index c2a15013..940003db 100755
--- a/Code/Source/svFSI/LELAS.f
+++ b/Code/Source/svFSI/LELAS.f
@@ -47,14 +47,17 @@ SUBROUTINE CONSTRUCT_LELAS(lM, Ag, Dg)
 
       INTEGER(KIND=IKIND), ALLOCATABLE :: ptr(:)
       REAL(KIND=RKIND), ALLOCATABLE :: xl(:,:), al(:,:), dl(:,:),
-     2   bfl(:,:), pS0l(:,:), pSl(:), N(:), Nx(:,:), lR(:,:), lK(:,:,:)
+     2   bfl(:,:), pS0l(:,:), pSl(:), N(:), Nx(:,:), lR(:,:), lK(:,:,:),
+     3   lVWP(:,:)
 
       eNoN = lM%eNoN
 
 !     LELAS: dof = nsd
       ALLOCATE(ptr(eNoN), xl(nsd,eNoN), al(tDof,eNoN), dl(tDof,eNoN),
      2   bfl(nsd,eNoN), pS0l(nsymd,eNoN), pSl(nsymd), N(eNoN),
-     3   Nx(nsd,eNoN), lR(dof,eNoN), lK(dof*dof,eNoN,eNoN))
+     3   Nx(nsd,eNoN), lR(dof,eNoN), lK(dof*dof,eNoN,eNoN),
+     4   lVWP(nvwp,eNoN))
+
 
 !     Loop over all elements of mesh
       DO e=1, lM%nEl
@@ -76,8 +79,15 @@ SUBROUTINE CONSTRUCT_LELAS(lM, Ag, Dg)
             dl(:,a)  = Dg(:,Ac)
             bfl(:,a) = Bf(:,Ac)
             IF (ALLOCATED(pS0)) pS0l(:,a) = pS0(:,Ac)
+!           Varwall properties------------------------------------------
+!           Calculate local wall property
+            IF (useVarWall) THEN
+               lVWP(:,a) = vWP0(:,Ac)
+            END IF
+!        ---------------------------------------------------------------
          END DO
 
+
 !        Gauss integration
          lR = 0._RKIND
          lK = 0._RKIND
@@ -92,7 +102,7 @@ SUBROUTINE CONSTRUCT_LELAS(lM, Ag, Dg)
             pSl = 0._RKIND
             IF (nsd .EQ. 3) THEN
                CALL LELAS3D(eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl, lR,
-     2            lK)
+     2            lK, lVWP)
 
             ELSE IF (nsd .EQ. 2) THEN
                CALL LELAS2D(eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl, lR,
@@ -128,22 +138,21 @@ SUBROUTINE CONSTRUCT_LELAS(lM, Ag, Dg)
       END SUBROUTINE CONSTRUCT_LELAS
 !####################################################################
       PURE SUBROUTINE LELAS3D (eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl,
-     2   lR, lK)
+     2   lR, lK, lVWP)
       USE COMMOD
       IMPLICIT NONE
       INTEGER(KIND=IKIND), INTENT(IN) :: eNoN
       REAL(KIND=RKIND), INTENT(IN) :: w, N(eNoN), Nx(3,eNoN),
-     2   al(tDof,eNoN), dl(tDof,eNoN), bfl(3,eNoN), pS0l(6,eNoN)
+     2   al(tDof,eNoN), dl(tDof,eNoN), bfl(3,eNoN), pS0l(6,eNoN),
+     3   lVWP(nvwp,eNoN)
       REAL(KIND=RKIND), INTENT(INOUT) :: pSl(6), lR(dof,eNoN),
      2   lK(dof*dof,eNoN,eNoN)
 
       INTEGER(KIND=IKIND) a, b, i, j, k
       REAL(KIND=RKIND) NxdNx, rho, elM, nu, lambda, mu, divD, T1, amd,
-     2   wl, lDm, ed(6), ud(3), f(3), S0(6), S(6)
+     2   wl, lDm, ed(6), ud(3), f(3), S0(6), S(6), eVWP(nvwp), Cst(6,6)
 
       rho  = eq(cEq)%dmn(cDmn)%prop(solid_density)
-      elM  = eq(cEq)%dmn(cDmn)%prop(elasticity_modulus)
-      nu   = eq(cEq)%dmn(cDmn)%prop(poisson_ratio)
       f(1) = eq(cEq)%dmn(cDmn)%prop(f_x)
       f(2) = eq(cEq)%dmn(cDmn)%prop(f_y)
       f(3) = eq(cEq)%dmn(cDmn)%prop(f_z)
@@ -151,16 +160,11 @@ PURE SUBROUTINE LELAS3D (eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl,
       j    = i + 1
       k    = j + 1
 
-      lambda = elM*nu / (1._RKIND+nu) / (1._RKIND-2._RKIND*nu)
-      mu     = elM * 0.5_RKIND / (1._RKIND+nu)
-      lDm    = lambda/mu
-      T1     = eq(cEq)%af*eq(cEq)%beta*dt*dt
-      amd    = eq(cEq)%am/T1*rho
-      wl     = w*T1*mu
-
       ed = 0._RKIND
       ud = -f
       S0 = 0._RKIND
+      eVWP = 0._RKIND
+
       DO a=1, eNoN
          ud(1) = ud(1) + N(a)*(al(i,a)-bfl(1,a))
          ud(2) = ud(2) + N(a)*(al(j,a)-bfl(2,a))
@@ -173,6 +177,14 @@ PURE SUBROUTINE LELAS3D (eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl,
          ed(5) = ed(5) + Nx(3,a)*dl(j,a) + Nx(2,a)*dl(k,a)
          ed(6) = ed(6) + Nx(1,a)*dl(k,a) + Nx(3,a)*dl(i,a)
 
+!     ------------------------------------------------------------------
+!     Calculate local wall property
+!     Don't use if calculating mesh
+         IF (useVarWall .AND. (phys_mesh .NE. eq(cEq)%phys)) THEN
+            eVWP(:) = eVWP(:) + N(a)*lVWP(:,a)
+         END IF
+!     ------------------------------------------------------------------
+
          S0(1) = S0(1) + N(a)*pS0l(1,a)
          S0(2) = S0(2) + N(a)*pS0l(2,a)
          S0(3) = S0(3) + N(a)*pS0l(3,a)
@@ -180,8 +192,26 @@ PURE SUBROUTINE LELAS3D (eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl,
          S0(5) = S0(5) + N(a)*pS0l(5,a)
          S0(6) = S0(6) + N(a)*pS0l(6,a)
       END DO
+
+
+      IF (useVarWall .AND. (phys_mesh .NE. eq(cEq)%phys)) THEN
+         elM  = eVWP(1)
+         nu   = evWP(2)
+      ELSE
+         elM  = eq(cEq)%dmn(cDmn)%prop(elasticity_modulus)
+         nu   = eq(cEq)%dmn(cDmn)%prop(poisson_ratio)
+      END IF
+      
+      lambda = elM*nu / (1._RKIND+nu) / (1._RKIND-2._RKIND*nu)
+      mu     = elM * 0.5_RKIND / (1._RKIND+nu)
+      lDm    = lambda/mu
+      T1     = eq(cEq)%af*eq(cEq)%beta*dt*dt
+      amd    = eq(cEq)%am/T1*rho
+      wl     = w*T1*mu
+
       divD = lambda*(ed(1) + ed(2) + ed(3))
 
+
 !     Stress in Voigt notation
       S(1) = divD + 2._RKIND*mu*ed(1)
       S(2) = divD + 2._RKIND*mu*ed(2)
@@ -189,6 +219,7 @@ PURE SUBROUTINE LELAS3D (eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl,
       S(4) = mu*ed(4)  ! 2*eps_12
       S(5) = mu*ed(5)  ! 2*eps_23
       S(6) = mu*ed(6)  ! 2*eps_13
+
       pSl  = S
 
 !     Add prestress contribution
@@ -253,7 +284,7 @@ PURE SUBROUTINE LELAS2D (eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl,
 
       INTEGER(KIND=IKIND) a, b, i, j
       REAL(KIND=RKIND) NxdNx, rho, elM, nu, lambda, mu, divD, T1, amd,
-     2   wl, lDm, ed(3), ud(2), f(2), S0(3), S(3)
+     2   T2, wl, lDm, ed(3), ud(2), f(2), S0(3), S(3)
 
       rho  = eq(cEq)%dmn(cDmn)%prop(solid_density)
       elM  = eq(cEq)%dmn(cDmn)%prop(elasticity_modulus)
@@ -296,6 +327,8 @@ PURE SUBROUTINE LELAS2D (eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl,
 !     Add prestress contribution
       S = pSl + S0
 
+!     Need to add variable wall tangent matrix
+
       DO a=1, eNoN
          lR(1,a) = lR(1,a) + w*(rho*N(a)*ud(1) + Nx(1,a)*S(1) +
      2      Nx(2,a)*S(3))
@@ -306,9 +339,9 @@ PURE SUBROUTINE LELAS2D (eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl,
          DO b=1, eNoN
             NxdNx = Nx(1,a)*Nx(1,b) + Nx(2,a)*Nx(2,b)
 
-            T1 = amd*N(a)*N(b)/mu + NxdNx
+            T2 = amd*N(a)*N(b)/mu + NxdNx
 
-            lK(1,a,b) = lK(1,a,b) + wl*(T1
+            lK(1,a,b) = lK(1,a,b) + wl*(T2
      2         + (1._RKIND + lDm)*Nx(1,a)*Nx(1,b))
 
             lK(2,a,b) = lK(2,a,b) + wl*(lDm*Nx(1,a)*Nx(2,b)
@@ -317,7 +350,7 @@ PURE SUBROUTINE LELAS2D (eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl,
             lK(dof+1,a,b) = lK(dof+1,a,b) + wl*(lDm*Nx(2,a)*Nx(1,b)
      2         + Nx(1,a)*Nx(2,b))
 
-            lK(dof+2,a,b) = lK(dof+2,a,b) + wl*(T1
+            lK(dof+2,a,b) = lK(dof+2,a,b) + wl*(T2
      2         + (1._RKIND + lDm)*Nx(2,a)*Nx(2,b))
          END DO
       END DO
diff --git a/Code/Source/svFSI/MATMODELS.f b/Code/Source/svFSI/MATMODELS.f
index 36733867..0842716a 100755
--- a/Code/Source/svFSI/MATMODELS.f
+++ b/Code/Source/svFSI/MATMODELS.f
@@ -38,7 +38,7 @@
 
 !     Compute 2nd Piola-Kirchhoff stress and material stiffness tensors
 !     including both dilational and isochoric components
-      SUBROUTINE GETPK2CC(lDmn, F, nfd, fl, ya, S, Dm)
+      SUBROUTINE GETPK2CC(lDmn, F, nfd, fl, ya, S, Dm, eVWP)
       USE MATFUN
       USE COMMOD
       IMPLICIT NONE
@@ -46,6 +46,7 @@ SUBROUTINE GETPK2CC(lDmn, F, nfd, fl, ya, S, Dm)
       INTEGER(KIND=IKIND), INTENT(IN) :: nfd
       REAL(KIND=RKIND), INTENT(IN) :: F(nsd,nsd), fl(nsd,nfd), ya
       REAL(KIND=RKIND), INTENT(OUT) :: S(nsd,nsd), Dm(nsymd,nsymd)
+      REAL(KIND=RKIND), INTENT(IN), OPTIONAL :: eVWP(nvwp)
 
       TYPE(stModelType) :: stM
       REAL(KIND=RKIND) :: nd, Kp, J, J2d, J4d, trE, p, pl, Inv1, Inv2,
@@ -129,7 +130,12 @@ SUBROUTINE GETPK2CC(lDmn, F, nfd, fl, ya, S, Dm)
 
 !     NeoHookean model
       CASE (stIso_nHook)
-         g1 = 2._RKIND * stM%C10
+         IF (useVarWall) THEN
+!           Converting elastic modulus and poisson ratio to g1
+            g1 = eVWP(1)* 0.5_RKIND/(1._RKIND+eVWP(2))
+         ELSE
+            g1 = 2._RKIND * stM%C10
+         END IF
          Sb = g1*IDm
 
 !        Fiber reinforcement/active stress
diff --git a/Code/Source/svFSI/MOD.f b/Code/Source/svFSI/MOD.f
index 743f25ad..ccb4c0a6 100755
--- a/Code/Source/svFSI/MOD.f
+++ b/Code/Source/svFSI/MOD.f
@@ -837,6 +837,8 @@ MODULE COMMOD
       LOGICAL cmmInit
 !     Whether variable wall properties are used for CMM
       LOGICAL cmmVarWall
+!     Whether variable wall properties should be read and used
+      LOGICAL useVarWall
 !     Whether shell equation is being solved
       LOGICAL shlEq
 !     Whether PRESTRESS is being solved
@@ -893,6 +895,8 @@ MODULE COMMOD
       INTEGER(KIND=IKIND) rsTS
 !     Number of stress values to be stored
       INTEGER(KIND=IKIND) nsymd
+!     Number of variable wall properties to read in from mesh
+      INTEGER(KIND=IKIND) nvwp
 
 !     REAL VARIABLES
 !     Time step size
@@ -933,7 +937,7 @@ MODULE COMMOD
       REAL(KIND=RKIND), ALLOCATABLE :: Ao(:,:)
 !     New time derivative of variables
       REAL(KIND=RKIND), ALLOCATABLE :: An(:,:)
-!     Old integrated variables (dissplacement)
+!     Old integrated variables (displacement)
       REAL(KIND=RKIND), ALLOCATABLE :: Do(:,:)
 !     New integrated variables
       REAL(KIND=RKIND), ALLOCATABLE :: Dn(:,:)
@@ -964,6 +968,9 @@ MODULE COMMOD
       REAL(KIND=RKIND), ALLOCATABLE :: pSn(:,:)
       REAL(KIND=RKIND), ALLOCATABLE :: pSa(:)
 
+!     Variables for variable wall properties
+      REAL(KIND=RKIND), ALLOCATABLE :: vWP0(:,:)
+
 !     Temporary storage for initializing state variables
       REAL(KIND=RKIND), ALLOCATABLE :: Pinit(:)
       REAL(KIND=RKIND), ALLOCATABLE :: Vinit(:,:)
diff --git a/Code/Source/svFSI/POST.f b/Code/Source/svFSI/POST.f
index cf1c7568..81477042 100755
--- a/Code/Source/svFSI/POST.f
+++ b/Code/Source/svFSI/POST.f
@@ -521,12 +521,14 @@ SUBROUTINE TPOST(lM, m, res, resE, lD, lY, iEq, outGrp)
       REAL(KIND=RKIND) w, Jac, detF, Je, ya, Ja, elM, nu, lambda, mu,
      2   p, trS, vmises, xi(nsd), xi0(nsd), xp(nsd), ed(nsymd),
      3   Im(nsd,nsd), F(nsd,nsd), C(nsd,nsd), Eg(nsd,nsd), P1(nsd,nsd),
-     4   S(nsd,nsd), sigma(nsd,nsd), Dm(nsymd,nsymd)
+     4   S(nsd,nsd), sigma(nsd,nsd), Dm(nsymd,nsymd), eVWP(nvwp),
+     5   sigma_temp(6), Cst(6,6)
       TYPE(fsType) :: fs
 
       INTEGER, ALLOCATABLE :: eNds(:)
       REAL(KIND=RKIND), ALLOCATABLE :: xl(:,:), dl(:,:), yl(:,:),
-     2   fN(:,:), resl(:), Nx(:,:), N(:), sA(:), sF(:,:), sE(:)
+     2   fN(:,:), resl(:), Nx(:,:), N(:), sA(:), sF(:,:),
+     3   sE(:), lVWP(:,:)
 
       dof  = eq(iEq)%dof
       i    = eq(iEq)%s
@@ -555,13 +557,15 @@ SUBROUTINE TPOST(lM, m, res, resE, lD, lY, iEq, outGrp)
 
       ALLOCATE (sA(tnNo), sF(m,tnNo), sE(lM%nEl), xl(nsd,fs%eNoN),
      2   dl(tDof,fs%eNoN), yl(tDof,fs%eNoN), fN(nsd,nFn), resl(m),
-     3   Nx(nsd,fs%eNoN), N(fs%eNoN))
+     3   Nx(nsd,fs%eNoN), N(fs%eNoN), lVWP(nvwp,fs%eNoN))
 
+      S    = 0._RKIND
       sA   = 0._RKIND
       sF   = 0._RKIND
       sE   = 0._RKIND
       insd = nsd
       ya   = 0._RKIND
+
       IF (lM%lFib) insd = 1
 
       DO e=1, lM%nEl
@@ -571,7 +575,7 @@ SUBROUTINE TPOST(lM, m, res, resE, lD, lY, iEq, outGrp)
      2       cPhys .NE. phys_ustruct .AND.
      3       cPhys .NE. phys_lElas) CYCLE
 
-         IF (cPhys .EQ. phys_lElas) THEN
+         IF (cPhys .EQ. phys_lElas .AND. .NOT. useVarWall) THEN
             elM = eq(iEq)%dmn(cDmn)%prop(elasticity_modulus)
             nu  = eq(iEq)%dmn(cDmn)%prop(poisson_ratio)
             lambda = elM*nu/(1._RKIND + nu)/(1._RKIND - 2._RKIND*nu)
@@ -594,6 +598,11 @@ SUBROUTINE TPOST(lM, m, res, resE, lD, lY, iEq, outGrp)
             xl(:,a) = x(:,Ac)
             dl(:,a) = lD(:,Ac)
             yl(:,a) = lY(:,Ac)
+!           Varwall properties -----------------------------------------
+!           Calculate local wall property
+            IF (useVarWall) THEN
+                lVWP(:,a) = vWP0(:,Ac)
+            END IF
          END DO
 
          Je = 0._RKIND
@@ -607,6 +616,7 @@ SUBROUTINE TPOST(lM, m, res, resE, lD, lY, iEq, outGrp)
 
             Im = MAT_ID(nsd)
             F  = Im
+            eVWP = 0._RKIND
             DO a=1, fs%eNoN
                IF (nsd .EQ. 3) THEN
                   F(1,1) = F(1,1) + Nx(1,a)*dl(i,a)
@@ -618,6 +628,12 @@ SUBROUTINE TPOST(lM, m, res, resE, lD, lY, iEq, outGrp)
                   F(3,1) = F(3,1) + Nx(1,a)*dl(k,a)
                   F(3,2) = F(3,2) + Nx(2,a)*dl(k,a)
                   F(3,3) = F(3,3) + Nx(3,a)*dl(k,a)
+!                 Variable wall ----------------------------------------
+!                 Calculate local wall property
+                  IF (useVarWall) THEN
+                     eVWP(:) = eVWP(:) + N(a)*lVWP(:,a)
+                  END IF
+
                ELSE
                   F(1,1) = F(1,1) + Nx(1,a)*dl(i,a)
                   F(1,2) = F(1,2) + Nx(2,a)*dl(i,a)
@@ -625,6 +641,12 @@ SUBROUTINE TPOST(lM, m, res, resE, lD, lY, iEq, outGrp)
                   F(2,2) = F(2,2) + Nx(2,a)*dl(j,a)
                END IF
             END DO
+            IF (cPhys .EQ. phys_lElas .AND. useVarWall) THEN
+               elM = eVWP(1)
+               nu  = eVWP(2)
+               lambda = elM*nu/(1._RKIND + nu)/(1._RKIND - 2._RKIND*nu)
+               mu     = 0.5_RKIND*elM/(1._RKIND + nu)
+            END IF
             detF = MAT_DET(F, nsd)
 
             ed = 0._RKIND
@@ -693,7 +715,6 @@ SUBROUTINE TPOST(lM, m, res, resE, lD, lY, iEq, outGrp)
                END IF
 
             CASE (outGrp_stress, outGrp_cauchy, outGrp_mises)
-               sigma = 0._RKIND
                IF (cPhys .EQ. phys_lElas) THEN
                   IF (nsd .EQ. 3) THEN
                      detF = lambda*(ed(1) + ed(2) + ed(3))
@@ -734,7 +755,8 @@ SUBROUTINE TPOST(lM, m, res, resE, lD, lY, iEq, outGrp)
                   IF (.NOT.ISZERO(detF)) sigma(:,:) = sigma(:,:) / detF
 
                ELSE IF (cPhys .EQ. phys_struct) THEN
-                  CALL GETPK2CC(eq(iEq)%dmn(cDmn), F, nFn, fN, ya, S,Dm)
+                  CALL GETPK2CC(eq(iEq)%dmn(cDmn), F, nFn, fN, ya,
+     2                     S,Dm, eVWP)
                   P1 = MATMUL(F, S)
                   sigma = MATMUL(P1, TRANSPOSE(F))
                   IF (.NOT.ISZERO(detF)) sigma(:,:) = sigma(:,:) / detF
diff --git a/Code/Source/svFSI/READFILES.f b/Code/Source/svFSI/READFILES.f
index c465aefe..77e1ddd5 100755
--- a/Code/Source/svFSI/READFILES.f
+++ b/Code/Source/svFSI/READFILES.f
@@ -97,6 +97,7 @@ SUBROUTINE READFILES
          pstEq        = .FALSE.
          sstEq        = .FALSE.
          ibFlag       = .FALSE.
+         useVarWall   = .FALSE.
 
          i = IARGC()
          IF (i .NE. 0) THEN
@@ -223,6 +224,10 @@ SUBROUTINE READFILES
 
 !--------------------------------------------------------------------
 !     Reading equations
+
+      IF ((nsd .NE. 3) .AND. useVarWall) err = "Variable wall "//
+     2      "properties can only be used in 3-dimensions"
+
       nEq = list%srch("Add equation",ll=1)
       std = " Number of equations: "//nEq
       ALLOCATE(eq(nEq))
@@ -414,6 +419,7 @@ SUBROUTINE READEQ(lEq, list, eqName)
 
       THflag = .FALSE.
       propL  = prop_NA
+
       SELECT CASE (eqName)
 !     FLUID Navier-Stokes solver ------------------------------------
       CASE ('fluid')
@@ -427,7 +433,7 @@ SUBROUTINE READEQ(lEq, list, eqName)
          IF (nsd .EQ. 3) propL(5,1) = f_z
          CALL READDOMAIN(lEq, propL, list)
 
-         nDOP = (/11,2,3,0/)
+         nDOP = (/12,2,3,0/)
          outPuts(1)  = out_velocity
          outPuts(2)  = out_pressure
          outPuts(3)  = out_WSS
@@ -440,6 +446,7 @@ SUBROUTINE READEQ(lEq, list, eqName)
          outPuts(9)  = out_viscosity
          outPuts(10) = out_divergence
          outPuts(11) = out_acceleration
+         outPuts(12) = out_displacement
 
          CALL READLS(lSolver_NS, lEq, list)
 
@@ -476,12 +483,24 @@ SUBROUTINE READEQ(lEq, list, eqName)
       CASE ('lElas')
          lEq%phys = phys_lElas
 
-         propL(1,1) = solid_density
-         propL(2,1) = elasticity_modulus
-         propL(3,1) = poisson_ratio
-         propL(4,1) = f_x
-         propL(5,1) = f_y
-         IF (nsd .EQ. 3) propL(6,1) = f_z
+         IF (useVarWall .AND. (nvwp .LT. 2._RKIND)) THEN
+            err = "Number of variable wall properties for linear "//
+     2         "elastic material must be at least 2."
+         END IF
+
+         IF (useVarWall) THEN
+            propL(1,1) = solid_density
+            propL(2,1) = f_x
+            propL(3,1) = f_y
+            IF (nsd .EQ. 3) propL(4,1) = f_z
+         ELSE
+            propL(1,1) = solid_density
+            propL(2,1) = elasticity_modulus
+            propL(3,1) = poisson_ratio
+            propL(4,1) = f_x
+            propL(5,1) = f_y
+            IF (nsd .EQ. 3) propL(6,1) = f_z
+         END IF
          CALL READDOMAIN(lEq, propL, list)
 
          lPtr => list%get(pstEq, "Prestress")
@@ -492,7 +511,7 @@ SUBROUTINE READEQ(lEq, list, eqName)
             outPuts(2)  = out_stress
             outPuts(3)  = out_strain
          ELSE
-            nDOP = (/8,2,0,0/)
+            nDOP = (/9,2,0,0/)
             outPuts(1) = out_displacement
             outPuts(2) = out_mises
             outPuts(3) = out_stress
@@ -501,6 +520,7 @@ SUBROUTINE READEQ(lEq, list, eqName)
             outPuts(6) = out_acceleration
             outPuts(7) = out_integ
             outPuts(8) = out_jacobian
+            outPuts(9) = out_cauchy
          END IF
 
          CALL READLS(lSolver_CG, lEq, list)
@@ -560,6 +580,8 @@ SUBROUTINE READEQ(lEq, list, eqName)
          IF (nsd .EQ. 3) propL(8,1) = f_z
          CALL READDOMAIN(lEq, propL, list)
 
+         IF (useVarWall) err = "Varible wall for USTRUCT is not "//
+     2      "implemented yet"
          lPtr => list%get(pstEq, "Prestress")
          IF (pstEq) err = "Prestress for USTRUCT is not implemented yet"
 
@@ -598,6 +620,9 @@ SUBROUTINE READEQ(lEq, list, eqName)
          propL(8,1) = f_z
          CALL READDOMAIN(lEq, propL, list)
 
+         IF (useVarWall) err = "Varible wall for SHELLS is not "//
+     2      "implemented yet"
+
          lPtr => list%get(pstEq, "Prestress")
          IF (pstEq) err = "Prestress for SHELLS is not implemented yet"
 
@@ -610,6 +635,7 @@ SUBROUTINE READEQ(lEq, list, eqName)
 
 !     COUPLED MOMENTUM FLUID STRUCTURE INTERACTION equation solver---
       CASE ('CMM')
+
          IF (nsd .NE. 3) err = "CMM eq. is not implemented for 2D "//
      2      "domains"
          lEq%phys = phys_CMM
@@ -701,6 +727,9 @@ SUBROUTINE READEQ(lEq, list, eqName)
          CALL READDOMAIN(lEq, propL, list)
          IF (cmmInit) lEq%dmn(:)%prop(solid_density) = 0._RKIND
 
+         IF (useVarWall) err = "Varible wall for CMM should not be "//
+     2      "defined in mesh data block"
+
          CALL READLS(lSolver_GMRES, lEq, list)
 
 !     FLUID STRUCTURE INTERACTION equation solver---------------------
@@ -2487,6 +2516,15 @@ SUBROUTINE READMATMODEL(lDmn, lPD)
          err = "Undefined constitutive model used"
       END SELECT
 
+      IF (useVarWall .AND. (lDmn%stM%isoType .NE. stIso_nHook)) THEN
+         err = "Variable wall properties currently only implemented "//
+     2         "for isotropic Neo-Hookean material in STRUCT."
+      END IF
+      IF (useVarWall .AND. (nvwp .LT. 2._RKIND)) THEN
+         err = "Number of variable wall properties for isotropic "//
+     2         "Neo-Hookean must be at least 2."
+      END IF
+
 !     Fiber reinforcement stress
       lFib => lPD%get(ctmp, "Fiber reinforcement stress")
       IF (ASSOCIATED(lFib)) THEN
diff --git a/Code/Source/svFSI/READMSH.f b/Code/Source/svFSI/READMSH.f
index e993668e..6220479a 100755
--- a/Code/Source/svFSI/READMSH.f
+++ b/Code/Source/svFSI/READMSH.f
@@ -50,7 +50,7 @@ SUBROUTINE READMSH(list)
       INTEGER(KIND=IKIND) :: i, j, iM, iFa, a, b, Ac, e, lDof, lnNo
       REAL(KIND=RKIND) :: maxX(nsd), minX(nsd), fibN(nsd), rtmp
       CHARACTER(LEN=stdL) :: ctmp, fExt
-      TYPE(listType), POINTER :: lPtr, lPM
+      TYPE(listType), POINTER :: lPtr, lPM, lPtr2
       TYPE(stackType) :: avNds
       TYPE(fileType) :: fTmp
 
@@ -408,6 +408,39 @@ SUBROUTINE READMSH(list)
          END DO
       END IF
 
+!     Read varwall properties ------------------------------------------
+      flag = .FALSE.
+      DO iM=1, nMsh
+         lPM => list%get(msh(iM)%name,"Add mesh",iM)
+         lPtr2 => lPM%get(nvwp,"Number of variable wall properties")
+         lPtr => lPM%get(cTmp, "Variable wall properties file path")
+         IF (ASSOCIATED(lPtr) .AND. ASSOCIATED(lPtr2)) THEN
+            IF (rmsh%isReqd) THEN
+               err = "Variable wall properties "//
+     2            "is not currently allowed with remeshing"
+            END IF
+            flag = .TRUE.
+            useVarWall = .TRUE.
+            std = "Setting varwall flag to true"
+            ALLOCATE(msh(iM)%x(nvwp,msh(iM)%gnNo))
+            msh(iM)%x = 0._RKIND
+            CALL READVTUPDATA(msh(iM), cTmp, "varWallProps", nvwp, 1)
+         END IF
+      END DO
+      IF (flag) THEN
+         ALLOCATE(vWP0(nvwp,gtnNo))
+         vWP0 = 0._RKIND
+         DO iM=1, nMsh
+            IF (.NOT.ALLOCATED(msh(iM)%x)) CYCLE
+            DO a=1, msh(iM)%gnNo
+               Ac = msh(iM)%gN(a)
+               vWP0(:,Ac) = msh(iM)%x(:,a)
+            END DO
+            DEALLOCATE(msh(iM)%x)
+         END DO
+      END IF
+!     ------------------------------------------------------------------
+
 !      Load any initial data (velocity, pressure, displacements)
       IF (.NOT.resetSim) CALL LOADVARINI(list)
 
diff --git a/Code/Source/svFSI/REMESH.f b/Code/Source/svFSI/REMESH.f
index 2911ab11..f693f0d7 100755
--- a/Code/Source/svFSI/REMESH.f
+++ b/Code/Source/svFSI/REMESH.f
@@ -394,6 +394,10 @@ END SUBROUTINE INTERP
       IF (ALLOCATED(pSn))      DEALLOCATE(pSn)
       IF (ALLOCATED(pSa))      DEALLOCATE(pSa)
 
+!     Varwall properties -----------------------------------------------
+      IF (ALLOCATED(vWP0))      DEALLOCATE(vWP0)
+!     ------------------------------------------------------------------
+
       t2 = CPUT()
 
       std = " Time taken for remeshing: "//STR(t2-t1)//" (s)"
diff --git a/Code/Source/svFSI/STRUCT.f b/Code/Source/svFSI/STRUCT.f
index 933abcf5..1afe9a67 100755
--- a/Code/Source/svFSI/STRUCT.f
+++ b/Code/Source/svFSI/STRUCT.f
@@ -50,7 +50,7 @@ SUBROUTINE CONSTRUCT_dSOLID(lM, Ag, Yg, Dg)
       INTEGER(KIND=IKIND), ALLOCATABLE :: ptr(:)
       REAL(KIND=RKIND), ALLOCATABLE :: xl(:,:), al(:,:), yl(:,:),
      2   dl(:,:), bfl(:,:), fN(:,:), pS0l(:,:), pSl(:), ya_l(:), N(:),
-     3   Nx(:,:), lR(:,:), lK(:,:,:)
+     3   Nx(:,:), lR(:,:), lK(:,:,:), lVWP(:,:)
 
       eNoN = lM%eNoN
       nFn  = lM%nFn
@@ -60,10 +60,12 @@ SUBROUTINE CONSTRUCT_dSOLID(lM, Ag, Yg, Dg)
       ALLOCATE(ptr(eNoN), xl(nsd,eNoN), al(tDof,eNoN), yl(tDof,eNoN),
      2   dl(tDof,eNoN), bfl(nsd,eNoN), fN(nsd,nFn), pS0l(nsymd,eNoN),
      3   pSl(nsymd), ya_l(eNoN), N(eNoN), Nx(nsd,eNoN), lR(dof,eNoN),
-     4   lK(dof*dof,eNoN,eNoN))
+     4   lK(dof*dof,eNoN,eNoN), lVWP(nvwp,eNoN))
+
 
 !     Loop over all elements of mesh
       DO e=1, lM%nEl
+
 !        Update domain and proceed if domain phys and eqn phys match
          cDmn  = DOMAIN(lM, cEq, e)
          cPhys = eq(cEq)%dmn(cDmn)%phys
@@ -84,6 +86,11 @@ SUBROUTINE CONSTRUCT_dSOLID(lM, Ag, Yg, Dg)
             yl(:,a)  = Yg(:,Ac)
             dl(:,a)  = Dg(:,Ac)
             bfl(:,a) = Bf(:,Ac)
+
+!           Varwall properties-----------------------------------------------------
+!           Calculate local wall property
+            IF (useVarWall) lVWP(:,a) = vWP0(:,Ac)
+!           ------------------------------------------------------------
             IF (ALLOCATED(lM%fN)) THEN
                DO iFn=1, nFn
                   fN(:,iFn) = lM%fN((iFn-1)*nsd+1:iFn*nsd,e)
@@ -107,7 +114,7 @@ SUBROUTINE CONSTRUCT_dSOLID(lM, Ag, Yg, Dg)
             pSl = 0._RKIND
             IF (nsd .EQ. 3) THEN
                CALL STRUCT3D(eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN,
-     2            pS0l, pSl, ya_l, lR, lK)
+     2                pS0l, pSl, ya_l, lR, lK, lVWP)
 
             ELSE IF (nsd .EQ. 2) THEN
                CALL STRUCT2D(eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN,
@@ -138,20 +145,20 @@ SUBROUTINE CONSTRUCT_dSOLID(lM, Ag, Yg, Dg)
       END DO ! e: loop
 
       DEALLOCATE(ptr, xl, al, yl, dl, bfl, fN, pS0l, pSl, ya_l, N, Nx,
-     2   lR, lK)
+     2   lR, lK, lVWP)
 
       RETURN
       END SUBROUTINE CONSTRUCT_dSOLID
 !####################################################################
       SUBROUTINE STRUCT3D(eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN,
-     2   pS0l, pSl, ya_l, lR, lK)
+     2   pS0l, pSl, ya_l, lR, lK, lVWP)
       USE COMMOD
       USE ALLFUN
       IMPLICIT NONE
       INTEGER(KIND=IKIND), INTENT(IN) :: eNoN, nFn
       REAL(KIND=RKIND), INTENT(IN) :: w, N(eNoN), Nx(3,eNoN),
      2   al(tDof,eNoN), yl(tDof,eNoN), dl(tDof,eNoN), bfl(3,eNoN),
-     3   fN(3,nFn), pS0l(6,eNoN), ya_l(eNoN)
+     3   fN(3,nFn), pS0l(6,eNoN), ya_l(eNoN), lVWP(nvwp,eNoN)
       REAL(KIND=RKIND), INTENT(OUT) :: pSl(6)
       REAL(KIND=RKIND), INTENT(INOUT) :: lR(dof,eNoN),
      2   lK(dof*dof,eNoN,eNoN)
@@ -159,7 +166,7 @@ SUBROUTINE STRUCT3D(eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN,
       INTEGER(KIND=IKIND) :: a, b, i, j, k
       REAL(KIND=RKIND) :: rho, dmp, T1, amd, afl, ya_g, fb(3), ud(3),
      2   NxSNx, BmDBm, F(3,3), S(3,3), P(3,3), Dm(6,6), DBm(6,3),
-     3   Bm(6,3,eNoN), S0(3,3)
+     3   Bm(6,3,eNoN), S0(3,3), eVWP(nvwp)
 
 !     Define parameters
       rho     = eq(cEq)%dmn(cDmn)%prop(solid_density)
@@ -181,6 +188,8 @@ SUBROUTINE STRUCT3D(eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN,
       F(3,3) = 1._RKIND
       S0     = 0._RKIND
       ya_g   = 0._RKIND
+      eVWP   = 0._RKIND
+
       DO a=1, eNoN
          ud(1) = ud(1) + N(a)*(rho*(al(i,a)-bfl(1,a)) + dmp*yl(i,a))
          ud(2) = ud(2) + N(a)*(rho*(al(j,a)-bfl(2,a)) + dmp*yl(j,a))
@@ -196,6 +205,11 @@ SUBROUTINE STRUCT3D(eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN,
          F(3,2) = F(3,2) + Nx(2,a)*dl(k,a)
          F(3,3) = F(3,3) + Nx(3,a)*dl(k,a)
 
+!     Varwall properits-------------------------------------------------
+!     Calculate local wall property
+         IF (useVarWall) eVWP(:) = eVWP(:) + N(a)*lVWP(:,a)
+!     ------------------------------------------------------------------
+
          S0(1,1) = S0(1,1) + N(a)*pS0l(1,a)
          S0(2,2) = S0(2,2) + N(a)*pS0l(2,a)
          S0(3,3) = S0(3,3) + N(a)*pS0l(3,a)
@@ -211,7 +225,7 @@ SUBROUTINE STRUCT3D(eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN,
 
 !     2nd Piola-Kirchhoff tensor (S) and material stiffness tensor in
 !     Voigt notationa (Dm)
-      CALL GETPK2CC(eq(cEq)%dmn(cDmn), F, nFn, fN, ya_g, S, Dm)
+      CALL GETPK2CC(eq(cEq)%dmn(cDmn), F, nFn, fN, ya_g, S, Dm, eVWP)
 
 !     Prestress
       pSl(1) = S(1,1)
diff --git a/Code/Source/svFSI/VTKXML.f b/Code/Source/svFSI/VTKXML.f
index 30317f43..7fbb80b0 100755
--- a/Code/Source/svFSI/VTKXML.f
+++ b/Code/Source/svFSI/VTKXML.f
@@ -255,7 +255,7 @@ SUBROUTINE READVTUPDATA(lM, fName, kwrd, m, idx)
       REAL(KIND=RKIND), ALLOCATABLE :: tmpR(:,:)
 
       iStat = 0
-      std = " <VTK XML Parser> Loading file <"//TRIM(fName)//">"
+      std = " <VTK XML Parser> Loading file <"//TRIM(fName)//"> "//m
       CALL loadVTK(vtu, fName, iStat)
       IF (iStat .LT. 0) err = "VTU file read error (init)"