# Climate model Debug matrix size issues

2 views (last 30 days)

Show older comments

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

##### 1 Comment

Mathieu NOE
on 26 Oct 2020

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

### Answers (2)

Saurav Chaudhary
on 27 Oct 2020

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-

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

##### 0 Comments

Landen Alexander
on 18 Oct 2021

Edited: Landen Alexander
on 28 Oct 2021

##### 0 Comments

### See Also

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!