Moderation and Mediation model

Question

I am trying to fit a moderation mediation model (please see the figure)

So I have a clasic mediation model and a moderator for the association between the mediator and the outcome.

I want to fit this model in lavaan and I know how to do the mediation model but I am struggling to add the moderator. Please see here the code for the mediation model but any guidance about how to include the moderator would be really appreciated:

set.seed(1234)
X <- rnorm(100)
M <- 0.5*X + rnorm(100)
Y <- 0.7*M + rnorm(100)
W <- 0.3*M + rnorm(100)
data <- data.frame(X = X, Y = Y, M = M, W = W)
data$MW <- data$M * data$W



model <- ' # direct effect
             Y ~ c*X
           # mediator
            Y ~ b*M
            M ~ a*X

           # indirect effect (a*b)
             ab := a*b
           # total effect
             total := c + (a*b)
             Y ~~ Y
             M ~~ M
         '
fit <- sem(model, data = data)
summary(fit, standardized=T)


Plot <- semPaths(fit, whatLabels = "est",
                 sizeMan = 10,
                 edge.label.cex = 1.15,
                 style = "ram",
                 nCharNodes = 0, nCharEdges = 0)

Thank you so much in advance.

UPDATE:

I have been working according to your comments/suggested paper and I have this code. Do you think it is correct?

set.seed(1234)
X <- rnorm(100)
M <- 0.5*X + rnorm(100)
Y <- 0.7*M + rnorm(100)
W <- 0.3*M + rnorm(100)
data <- data.frame(X = X, Y = Y, M = M, W = W)
data$MW <- data$M * data$W
data <-as.data.frame(scale(data, scale=FALSE))
describe(data)



model <- ' # direct effect
             Y ~ 1+ c*X + w*W +wb*MW
           # mediator
            Y ~ b*M 
            M ~ a*X

           # indirect effect (a*b)
             ab := a*b
             wab := a*wb
             
             
             
           # total effect
             total := c + w + wb+ (a*b)
            #Variances an correlations
             Y ~~ Y
             M ~~ M
             X ~~ X
             MW~~MW
             MW ~~ X
             MW ~~ W
            X ~~ W
         '
fit <- sem(model, data = data)
summary(fit, standardized=T)


Plot <- semPaths(fit, whatLabels = "est",
                 sizeMan = 10,
                 edge.label.cex = 1.15,
                 style = "ram",
                 nCharNodes = 0, nCharEdges = 0)