Climate model Debug matrix size issues

2 views (last 30 days)
Theodore Anderson
Theodore Anderson on 24 Oct 2020
Edited: Landen Alexander on 28 Oct 2021
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?
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
% 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
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
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')
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')
semilogy(Rnet, Z_atm)
semilogy(T_atm(t, :), Z_atm)
%% plot the end results
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
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')
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;
  1 Comment
Mathieu NOE
Mathieu NOE on 26 Oct 2020
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

Sign in to comment.

Answers (2)

Saurav Chaudhary
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-

Landen  Alexander
Landen Alexander on 18 Oct 2021
Edited: Landen Alexander on 28 Oct 2021
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.

Community Treasure Hunt

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

Start Hunting!