After adding the channel model to the code the BER vs SNR curve doesn't look realistic. What went wrong in the code?

clc ;
clear all ;
close all ;
m =512; % Total number of OFDM symbols
N =1024; % Length of each OFDM symbol
M =16; % Size of the Constellation ( M can be 4 , 8 , 16 ,32 , 64 , 128 , 256)
Ncp =256 ; % Length of the Cyclic Prefix

Data = randi ([0 M-1] ,m , N ); % Generation of Random bits Matrix of size m by N

DataMod = qammod( Data , M);

DataMod_serialtoparallel = DataMod .'; % Performing Serial to Parallel Conversion
datamat = DataMod_serialtoparallel ; % Assigning the total data to a variable called datamat % Computation of Hermitian Symmetry Criteria
datamat (1 ,:) =0; % Assigning the First subcarrier to Zero
datamat (513 ,:) =0; % Assigning the Middle Subcarrier to Zero
datamat (514:1024 ,:) = flipud ( conj ( datamat (2:512 ,:) ) ) ; % Illustrating that only half of the subcarriers are exploited for data transmission as the remaining half are flipped complex conjugate versions of the previous ones .
d_ifft = ifft (( datamat ) ) ; % Computation of IFFT operation
d_ifft_paralleltoserial = d_ifft .'; % Parallel to Serial Conversion
CP_part = d_ifft_paralleltoserial (: ,end - Ncp +1: end) ; % Addition of Cyclic Prefix
DCOOFDM_CP =[ CP_part d_ifft_paralleltoserial ]; %Transmissin of DCO - OFDM signal

%% underwater VLC channel modelling
theta=70;
% semi-angle at half power
ma=-log10(2)/log10(cosd(theta));
%Lambertian order of emission
P_total=20;
%transmitted optical power by individual LED
Adet=1e-2;
%assumed a custom made PD with a detector physical area of 1cm^2
Ts=1;
%assumed an ideal optical filter
index=1.5;
%refractive index of a lens at a PD, value supported by paper 3.
C = 0.1162;
%underwater attenuation coefficient, value was taken from paper 3.
FOV=60*pi/180;
%FOV of a receiver
G_Con=(index^2)/sin(FOV);
%gain of an optical concentrator; ignore if no lens is used
lx=5; ly=5; lz=10;
% Underwater zone dimension in metre
h=10;
%the distance between source and receiver plane
XT=0; YT=0;
% position of LED;
Nx=lx*10; Ny=ly*10;
% number of grid in the receiver plane
x=-lx/2:lx/Nx:lx/2;
y=-ly/2:ly/Ny:ly/2;
[XR,YR]=meshgrid(x,y);
% receiver plane grid

D1=sqrt((XR-XT(1,1)).^2+(YR-YT(1,1)).^2+h^2);
% distance vector from source 1
cosphi_A1=h./D1;
% angle vector
H_A1=(ma+1)*Adet.*cosphi_A1.^(ma)./(2*pi.*D1.^2);
% channel DC gain for source 1
H_uw=H_A1.*exp(-C*D1);
% underwater channel DC gain for source 1; equation was supported 1,2,3

%Idea from same ACO-OFDM code
H_A2 = H_uw ./ norm ( H_uw );
d_channell = filter ( H_A2  (1 ,1:2),1 , DCOOFDM_CP .') .';

%%
bdc =7; %DC BIAS
clip = sqrt ((10.^( bdc /10) ) -1) ; % clipping factor k
bdcc = clip * sqrt (d_channell .* d_channell ) ;%Computation of DC bias
DCOOFDM_BIAS = bdcc + d_channell ; % Addition of DC bias to the cyclic prefix added signal
count =0;
snr_vector =0:1:30; % size of signal to noise ratio (SNR ) vector
for snr = snr_vector
    SNR = snr + 10* log10 ( log2 ( M ) ) ;
    count = count +1 ;
    DCOOFDM_with_chann = awgn ( DCOOFDM_BIAS , SNR ,'measured'); % Addition of AWGN
    
    % Receiver of DCO - OFDM
    DCOOFDM_with_chann1 = DCOOFDM_with_chann - bdcc ; %Removal of DC bias
    DCOOFDM_removal_CP = DCOOFDM_with_chann1 (: , Ncp +1: N + Ncp ) ; % Removal of Cyclic Prefix
    DCOOFDM_serialtoparallel = DCOOFDM_removal_CP .'; %Serial to Parallel Conversion
    
    DCOOFDM_parallel_fft = fft ( DCOOFDM_serialtoparallel) ; % Computation of FFT operation
    
    DCOOFDM_Demodulation = qamdemod (DCOOFDM_parallel_fft .', M );
    [~ , s_e1(count)]= symerr ( Data (: ,2:512),DCOOFDM_Demodulation (: ,2:512)) ;
end

% Plotting the BER curves
semilogy ( snr_vector , s_e1 ,'rd -','LineWidth' ,2); 
legend('FFT - based DCO -OFDM -7 dB of Bias') ;
axis ([0 30 10^ -4 1]) ;
xlabel ('SNR in dB') ;
ylabel ('BER') ;
grid on ;