Climate model Debug matrix size issues

Hey all! I have a really complex model that runs, but does not run very well. I get some errors where my matrix values exceed 1 one where it cant exceed one. It took a long time to get this to run, I was hoping you guys could give me some general feed back on how to clean up my code :) let me know if you have suggestions or can fix errors in my equations.

% function columnmodel = 

PLOTSON = 1; % change this to 0 to hide debugging plots

% physical constants
sigma = 5.67e-8;
Cp_air = 1e3; %j/kgK
Cp_water = 4060e3; % j/kgK
rho_water = 1000; %kg/m3

% set up the model space
N_levels = 25; % how many levels of atmosphere

Z_atm = [ logspace(-2, log10(30), N_levels-2) 50 100]; 
% in kilometers, vertical spacing of layers: logarithmic between 10m and
% 30km, then two extra at 50 and 100km

% temperature and emissivity of the atmosphere
T_atm(Z_atm < 12) = 288 - Z_atm(Z_atm < 12)*6.5; %Kelvin, starting with a 6.5°C/km lapse rate
T_atm(Z_atm >= 12) = 210; % set constant above 12km
T_surf = 288; % Kelvin

rho_air = 1.225*exp(-9.81*.029*Z_atm*1000/(8.31*288)); %kg/m^3 density of air (approximate)
% Constant conditions
I_toa = 1361; % W/m^2, solar radiation at the top of the atmosphere total radiation 
%five times globe average = daytime-ish
I_SW_toa = .9*I_toa; % 90% shortwave
I_LW_toa = .1*I_toa; %10% long wave

Emissivity_LW_atmos = .5*ones(N_levels, 1);% emissivity of the atmosphere, LW
absorptionsw = 5e-4*ones(N_levels, 1); %absorption in SW
absorptionlw = 4e-1*ones(N_levels, 1); % absorption in LW
scattering = 2e-4*ones(N_levels, 1); % layer by layer reflectivity (scattering) of the atmosphere
%(if you want to add clouds, change this)
Ka = absorptionsw;
Ks = scattering;
% shortcuts for the 
ksw = (scattering+absorptionsw);
klw = absorptionlw;

%surface conditions
surface_albedo = .2;
E_surf = 1;

% time and vertical space step settings
dt = 15*60; % time step in seconds = 15 minutes
dz = [.01 diff(Z_atm) ]*1e3;  % vertical layer thickness in m, starting with a 10m surface layer

% thiq = ((Ka+Ks)*rho_air*dz);% optical thickness
% thiqlw = (Ka*rho_air*dz);
% step through time steps
for t = 1:30
    % EM radiation moves faster than a time step. 
    %   Model all of this as if it happens instantaneously
    %   using v for the atmosphere level counter 
    v = N_levels; % very top of the atmosphere   
    % initiate some variables as column vectors
            SWd = zeros(N_levels, 1);  % downwelling SW
            SWu = zeros(N_levels, 1);   % upwelling SW
            SWa = zeros(N_levels, 1);   % absorbed SW
            SWe = zeros(N_levels, 1);   % emitted SW

            LWd = zeros(N_levels, 1);  % downwelling LW
            LWu = zeros(N_levels, 1);   % upwelling LW
            LWa = zeros(N_levels, 1);   % absorbed LW
            LWe = zeros(N_levels, 1);   % emitted LW
    % very top of the atmosphere   
    SWd(v) = I_SW_toa*exp(-ksw(v)*rho_air(v)*dz(v));
    SWa(v) = I_SW_toa*((Ka(v)*rho_air(v)*dz(v)));
%     SWe(v) = %atmalbedo?
    %lw
    LWe(v) = Emissivity_LW_atmos(v)*sigma*T_atm(v)^4;
    LWd(v) = I_LW_toa*exp(-klw(v)*rho_air(v)*dz(v));
    LWa(v) = I_LW_toa*(1-exp(-(Ka(v)*rho_air(v)*dz(v))));
    % downwelling
    for v = (N_levels-1):-1:1 % start one level down from the top, increment by -1
    SWd(v) = SWd(v+1)*exp(-ksw(v)*rho_air(v)*dz(v));%energy going down 
    SWa(v) = SWd(v+1)*(-ksw(v)*rho_air(v)*dz(v));%energy down * thiqness (essentially trasmitted thru a layer?) 
    SWe(v) = 0;
    SWu(v) = 0; % how much is reflected back up??? why do we need this?

%     LWe(v) = 
    LWd(v) = I_LW_toa*exp(-klw(v)*rho_air(v)*dz(v));
    LWa(v) = I_LW_toa*1-exp(-(Ka(v)*rho_air(v)*dz(v)));   
    % nothing is reflected back up in the LW
    end
    % at the surface level
    v  = 1;
    % surface reflection, absorption, emission 
    SWu(1) = SWu(1)+surface_albedo*SWd(1)*exp(-ksw(v)*rho_air(v)*dz(v));
    SWa(1) = surface_albedo*SWd(1)*(ksw(v)*rho_air(v)*dz(v))+SWa(1);
    SWe(1) = 0; % no short wave emission 
    SurfSWa =  SWd(1)*(ksw(v)*rho_air(v)*dz(v))*(1-surface_albedo); % absorbed at the surface
    SurfSWe =  SWd(1)*(ksw(v)*rho_air(v)*dz(v))*surface_albedo; % emitted at the surface
    
    LWe(1) = LWe(1)+Emissivity_LW_atmos(1)*sigma*T_atm(t,v)^4;   
    SurfLWe =  LWd(1)*(klw(v)*rho_air(v)*dz(v))%E_surf(1)*sigma*T_surf(1)^4 + LWe(1)*exp(-(Ka*rho_air*dz));
    LWu(1) =Emissivity_LW_atmos(v)*sigma*T_atm(t,v)^4 + (SurfLWe(v)*exp(-(Ka(v)*rho_air(v)*dz(v)))); % just the part emitted that is going up
    LWa(1) = LWd(v+1)*(klw(v)*rho_air(v)*dz(v)); 
    SurfLWa =  LWd(v+1)*(klw(v)*rho_air(v)*dz(v));
    % upwelling 
    for v = 2:N_levels
            % upwelling short wave
            SWu(v) = SWu(v-1)*exp(-((Ka(v)+Ks(v))*rho_air(v-1)*dz(v-1)))+SWu(v-1); % how much gets passed up
            SWa(v) = SWu(v-1)*Ka(v)*rho_air(v-1)*dz(v-1)+SWa(v-1); % how much gets absorbed here
            % make sure to add it to what you absorbed on the way down
            % upwelling long wave
            LWe(v) = LWe(v)+Emissivity_LW_atmos(v)*sigma*T_atm(t,v)^4; % how much is emitted here pointed up
            LWu(v) = Emissivity_LW_atmos(v)*sigma*T_atm(t,v)^4 + LWu(v-1)*exp(-(Ka(v-1)*rho_air(v-1)*dz(v-1))); % how much gets passed up
            LWa(v) = LWd(v-1)*(klw(v)*rho_air(v)*dz(v)); % how much gets absorbed here    
            % nothing gets reflected down in the LW
    end
    
    Rnet = LWa + SWa - LWe - SWe; % add up the components of the net radiative balance
    Rnet_surf = SurfLWa - SurfLWe + SurfSWa - SurfSWe;
    %% calculate the temperature at the next time step from the current
    % update temperatures
    T_atmos(t+1, :) = Rnet;
              T_atmos(t+1, N_levels) = 210; % fix the top two because 
              % numerical instabilities add up at the edge of space. 
              T_atmos(t+1, N_levels-1) = 210;
    % surface temp
    T_surf(t+1) = Rnet_surf;      
    %% show results for debugging
    if PLOTSON
       
        figure(1)
        clf;
        subplot(1,4,1)
        semilogy(LWa, Z_atm, 'r')
        hold on
        semilogy(LWe, Z_atm, 'b')
        semilogy(LWu, Z_atm, 'g')
        semilogy(LWd, Z_atm, 'k')
        hold off
        xlabel('LW ')
        legend('absorbed', 'emitted', 'upwelling', 'downwelling')
        
        subplot(1,4,2)
        semilogy(SWa, Z_atm, 'r')
        hold on
        semilogy(SWe, Z_atm, 'b')
        semilogy(SWu, Z_atm, 'g')
        semilogy(SWd, Z_atm, 'k')
        hold off
        xlabel('SW ')
        legend('absorbed', 'emitted', 'upwelling', 'downwelling')
        
        subplot(1,4,3)
             
        semilogy(Rnet, Z_atm)

        xlabel('Rnet')
        
                
        subplot(1,4,4)
             
        semilogy(T_atm(t, :), Z_atm)

        xlabel('Temp')
        
    end
        
end
%% plot the end results
figure(3)
subplot(3,4,[1 2 3 5 6 7])
pcolor(1:t+1, Z_atm, (T_atm)' - repmat(T_atm(:,1)',length(kl), 1 ))
%set(gca, 'Ydir', 'normal');
ylim([0 30])
shading flat
colorbar;
caxis([-50 50])

subplot(3,4,[4 8])
plot(T_atm(end, 1:end-3)-273, Z_atm(1:end-3))
ylim([0 30])
xlabel('Temp, C')

subplot(3,4,[9:12])
plot((1:length(T_surf))*dt/(60*60), T_surf - 273)
xlabel('Hour of simulation')
ylabel('Surface temp')


% for making this a function:        
% tsurf = T_surf(end);
% tatmos = T_atmos(end,:);
% zatmos = Z_atmos;

%end

[Mathieu NOE]

hi

for such a complex code , it would be good if you could put some explanations at the beginning of your file about variables and their size.

next, you should do preallocation for variables that grow inside a loop like T_atmos and T_surf

Answer 1

[Saurav Chaudhary]

Identify bottlenecks by using the Profiler tool within MATLAB.

To know more about profiling refer below points:

• Profiling is a way to measure the time it takes to run your code and identify where MATLAB spends the most time.
• After you identify which functions are consuming the most time, you can evaluate them for possible performance improvements.
• You also can profile your code to determine which lines of code do not run.

You can refer this documentation link to see how to run profiler on your code-

https://www.mathworks.com/help/matlab/matlab_prog/profiling-for-improving-performance.html

You can also refer to this documentation page for better techniques that can be followed for MATLAB programming-

https://www.mathworks.com/help/matlab/matlab_prog/techniques-for-improving-performance.html

Answer 2

[Landen Alexander]

The ecology and climate in general are now in a very dangerous state. The entire planet is affected by humans and there is very little territory that would not change under the influence of human actions. You can click here and read that natural land has been in the industrialization stage for quite some time. While everyone says they want the earth to remain in its natural state, the scientist White explains that there isn't really a single place left unaffected by human development. The overall well-being of the land and nature must ultimately be a joint effort. No power can solve a dilemma on its own.