How to implement a system-dependent upper boundary condition (switching between Neumann and Dirichlet) for 1D Richards' equation using implicit finite difference?

I am solving the 1D mixed-form Richards' equation in MATLAB using a fully implicit finite difference scheme (Picard iteration).

My boundary conditions are:

Bottom Boundary (z=80 cm): the lower boundary was set as a variable pressure head with a low negative pressure head of 51 cm.

Top Boundary (z=0 cm): Precipitation flux q = -P(t).

My question:

Specifically, how should I write the MATLAB code to handle the top boundary condition, where there is a time-dependent precipitation flux P(t), subject to a maximum allowable surface pressure head of h = 20 cm.

% 1D Integer-order Richards' Equation (Picard Iteration Scheme)
clc;
clear; close all;

L = 80;                % Column length [cm]
T = 816;               % Total time [h]

Ks = 3.2;              % cm/h (Saturated hydraulic conductivity)
alpha_v = 0.076;       % 1/m  (van Genuchten parameter; if using cm throughout, strictly it should be converted to 1/cm)
nv = 1.86;
mv = 1 - 1/nv;
theta_s = 0.320;
theta_r = 0.032;

h_bot = -51;           % cm (Constant pressure head at the bottom boundary)

% Initial water content -> Initial pressure head
theta_init = 0.08;
Se = (theta_init-theta_r)/(theta_s-theta_r);
Psi_init = -( (Se^(-1/mv)-1)^(1/nv) )/alpha_v;  % cm (Negative)

dz = 0.1;              % cm
dt = 0.05;             % h

z = 0:dz:L;
t = 0:dt:T;
M = length(z);
N = length(t);

% Apply upper boundary pulse infiltration (Flux, positive downwards, unit: cm/h)
qTop = zeros(1,N);
q_pulse = 15/0.5;                  % 30 cm/h
events = [0 0.5; 240 240.5; 480 480.5];  % hours
for k = 1:size(events,1)
    qTop(t>=events(k,1) & t<=events(k,2)) = q_pulse;
end

% van Genuchten-Mualem Model
Kfun = @(h) (Ks .* ((1 - (alpha_v * abs(h)).^(nv-1) .* (1 + (alpha_v * abs(h)).^nv).^(-mv)).^2) ./ ...
               (1 + (alpha_v * abs(h)).^nv).^(mv/2));
Cfun = @(h) (theta_s - theta_r) .* alpha_v .* nv .* mv .* ((alpha_v .*abs(h)).^(nv-1)) .* ((1 + ( alpha_v *abs(h)).^nv).^(-mv-1));
Qfun = @(h) (theta_r + (theta_s - theta_r) .* (1 + (alpha_v *abs(h)).^nv).^(-mv));

% Result arrays
h = zeros(M, N);
theta = zeros(M, N);
q = zeros(M, N);

h(:,1) = Psi_init;
theta(:,1) = Qfun(Psi_init);

% Picard parameters
maxIter = 200;
tol = 1e-6;                  
iterate = zeros(N,1);

% Workspace for iterations
KK = zeros(M,1);
CC = zeros(M,1);

% Main time-stepping loop
for n = 1:N-1

    % Previous time step (converged solution)
    h_old = h(:,n);

    % Picard initial guess: use previous time step
    h_it = h_old;

    for it = 1:maxIter

        % Update nonlinear coefficients using current iteration solution (Picard: fix coefficients, solve linear equation for new h)
        for j = 1:M
            if h_it(j) > 0
                KK(j) = Ks;
                CC(j) = 0;
            else
                KK(j) = Kfun(h_it(j));
                CC(j) = Cfun(h_it(j));
            end
        end

        % Assemble tridiagonal linear system B*h_new = d
        B = zeros(M,M);
        d = zeros(M,1);

        % q = -K * (dh/dz + 1)
        h_max = 20;   % Maximum allowable pressure head (ponding depth)
        
        K12 = 0.5*(KK(1)+KK(2));
        qT = qTop(n+1);
        
          
         if qT > 0  % Water is being applied
            if h_it(1) < h_max
                B(1,1) = -1; B(1,2) = 1;
                d(1)   = dz * (qT/max(K12,1e-12) + 1); % [Correction] z is positive downwards, gravity term is -1
            else
                B(1,1) = 1; B(1,2) = 0; 
                d(1)   = h_max; % Ponding limit reached, force constant head
            end
        else       % Water application stopped
            if h_it(1) > 0
                B(1,1) = 1; B(1,2) = 0; 
                d(1)   = 0; % Surface is ponded, let it infiltrate under potential energy, set head to 0
            else
                B(1,1) = -1; B(1,2) = 1;
                d(1)   = dz * (0/max(K12,1e-12) + 1); % No ponding, zero flux boundary (qT=0)
            end
        end


        for j = 2:M-1
            Kjp = 0.5*(KK(j)+KK(j+1));   % K_{j+1/2}
            Kjm = 0.5*(KK(j)+KK(j-1));   % K_{j-1/2}

            B(j,j-1) = -(dt/(dz^2))*Kjm;
            B(j,j)   =  CC(j) + (dt/(dz^2))*(Kjm + Kjp);
            B(j,j+1) = -(dt/(dz^2))*Kjp;

            d(j) = CC(j)*h_old(j) - (dt/dz)*(Kjp - Kjm);
        end

        % Bottom boundary: Constant head
         B(M, M-1) = 0;
        B(M,M) = 1;
        d(M)   = h_bot;
        h_new = ThomasAlgorithm(B, d);

        % Convergence criterion
        err = max(abs(h_new - h_it));
        h_it = h_new;

        if err < tol
            iterate(n+1) = it;
            break;
        end
    end

    % Save converged solution for the current time step
    h(:,n+1) = h_it;

    % Calculate flux q
    for j = 1:M
        h_j = h(j, n+1);
        if h_j > 0
            K_j = Ks;
        else
            K_j = Kfun(h_j);
        end

        if j == 1
            dh_dz = (h(j+1, n+1) - h_j) / (-dz);
        elseif j == M
            dh_dz = (h_j - h(j-1, n+1)) / (-dz);
        else
            dh_dz = (h(j+1, n+1) - h(j-1, n+1)) / (-2*dz);
        end

        q(j, n+1) = -K_j * (dh_dz + 1);
    end
    
  for j = 1:M
       if h(j,n+1)<0
        theta(j,n+1)=Qfun(h(j,n+1));
       else
        theta(j,n+1)=theta_s;   
       
       end
  end

    fprintf('Picard-Richards: Time step %d / %d, Iterations: %d\n', n, N-1, max(1,iterate(n+1)));
end

Unrecognized function or variable 'ThomasAlgorithm'.