I removed the with torch.no_grad() statement and just detached the mean and var tensors.
Also, I’ve added a condition for the save method as used in this example 2 and scripting seems to work.
Here is a reduced code snippet:
@torch.jit.ignore
def save_statisctics(x : torch.Tensor):
    print('save called')


class BatchNorm2d(nn.BatchNorm2d):
    def __init__(self, num_features, eps=1e-5, momentum=0.1,
                 affine=True, track_running_stats=True, save=False):
        super(BatchNorm2d, self).__init__(
            num_features, eps, momentum, affine, track_running_stats)
        self.save = save

    #@torch.jit.script_method
    def forward(self, input):
        self._check_input_dim(input)

        exponential_average_factor = 0.0

        if self.training and self.track_running_stats:
            if self.num_batches_tracked is not None:
                self.num_batches_tracked += 1
                if self.momentum is None:  # use cumulative moving average
                    exponential_average_factor = 1.0 / float(self.num_batches_tracked)
                else:  # use exponential moving average
                    exponential_average_factor = self.momentum

        # calculate running estimates
        if self.training:

            if self.save:
                save_statisctics(self.running_mean)

            mean = input.mean([0, 2, 3]).detach()
            # use biased var in train
            var = input.var([0, 2, 3], unbiased=False).detach()
            n = input.numel() / input.size(1)
            
            self.running_mean = exponential_average_factor * mean\
                + (1 - exponential_average_factor) * self.running_mean
            # update running_var with unbiased var
            self.running_var = exponential_average_factor * var * n / (n - 1)\
                + (1 - exponential_average_factor) * self.running_var
        else:
            mean = self.running_mean
            var = self.running_var

        input = (input - mean[None, :, None, None]) / (torch.sqrt(var[None, :, None, None] + self.eps))
        if self.affine:
            input = input * self.weight[None, :, None, None] + self.bias[None, :, None, None]

        return input


class CNNNet_CIFAR_MyBN(nn.Module):

    def __init__(self, conv1_dim=100, save=False):
        super(CNNNet_CIFAR_MyBN, self).__init__()
        self.conv1 = nn.Conv2d(3, 1, 5, stride=1, padding=2)
        self.normalize1 = BatchNorm2d(conv1_dim, save=save)

    def forward(self, x):
        x = F.relu(self.normalize1((self.conv1(x))))
        return x


model = CNNNet_CIFAR_MyBN(save=False)
x = torch.randn(1, 3, 32, 32)
out = model(x)
scripted_model = torch.jit.script(model)
out = scripted_model(x)

Let me know, if this would work for you.