Closed HenryDane closed 10 months ago
After a lot of updates, I've gotten to this:
I removed adv/diff, switched to A+BT, and switched to P2 albedo.
I'm genuinely at a loss as for why the results are different.
Comparisons of relevant quantities are shown below:
Ts.min(), Ts.max()
# (219.03, 320.48)
ebm.Ts.min(), ebm.Ts.max()
# (Field(6.78257524), Field(106.71048209))
# ...
Qs.max(), Qs.min()
# (1413.4, -0.0)
ebm.insolation.max(), ebm.insolation.min()
# (Field(1405.65663901), Field(0.))
# (array(1413.44308192), array(0.))
# ^ (climlab insol evaluated at higher res)
# ...
alphas.max(), alphas.min()
# (0.57999, 0.20501)
ebm.albedo.max(), ebm.albedo.min()
# (Field(0.57639724), Field(0.20860276))
# ...
((1 - alphas) * Qs).max(), ((1 - alphas) * Qs).min()
# (1079.900157, -0.0)
ebm.ASR.max(), ebm.ASR.min()
# (Field(1070.57541206), Field(0.))
# ...
OLRs.max(), OLRs.min()
# (294.6600000000001, 91.76000000000005)
ebm.OLR.max(), ebm.OLR.min()
# (Field(423.42774292), Field(223.57798349))
# ...
dTdts.min(), dTdts.max()
# (-2.4288e-06, 6.2807e-06)
ebm.tendencies['Ts'].min(), ebm.tendencies['Ts'].max()
# (Field(-3.37556695e-06), Field(5.29673904e-06))
The full output is attached below. It seems that the temperature slowly decreases over the model's runtime. I do not yet know why.
As of 2e72a39
(on ebm_work
) I get same results as:
ebm = climlab.EBM(num_lat=16, num_lon=16, timestep=climlab.constants.seconds_per_hour, water_depth=30)
iins = climlab.radiation.insolation.InstantInsolation(domains=ebm.Ts.domain,timestep=ebm.time['timestep'])
ebm.add_subprocess('insolation', iins)
ebm.remove_subprocess('diffusion')
ebm.subprocess('SW').insolation = iins.insolation
ebm.remove_subprocess('albedo')
albpr = climlab.surface.albedo.P2Albedo(state=ebm.state)
ebm.add_subprocess('albedo', albpr)
ebm.subprocess('SW').albedo = albpr.albedo
# same initial condition
ebm.Ts[:] = 0.0
# run
ebm.integrate_years(3.0)
A timestep of 5 mins is shown below. Using a timestep of 1 hour produces oscillations in the right-hand plot such that it has 24 peaks.
As of 1b430d8
(on ebm_work
) I get the same results as:
ebm = climlab.EBM(num_lat=16, num_lon=16, timestep=climlab.constants.seconds_per_hour, water_depth=30)
iins = climlab.radiation.insolation.InstantInsolation(domains=ebm.Ts.domain,timestep=ebm.time['timestep'])
ebm.add_subprocess('insolation', iins)
ebm.remove_subprocess('diffusion')
ebm.subprocess('SW').insolation = iins.insolation
# same initial condition
ebm.Ts[:] = 0.0
# run
ebm.integrate_years(3)
Note that this requires 4c4a942
on ebm_work
in HenryDane/climlab
.
Resulting plot comparison is:
The problem with the adv-diff schema is the boundary conditions. climlab
uses a tri-diagonal solver which uses fluxes to solve the adv-diff equation and also assigns the fluxes at the boundaries to zero. I need to implement this.
Ref: https://climlab.readthedocs.io/en/latest/api/climlab.dynamics.adv_diff_numerics.html
As of 109c6de
(on ebm_work
), model now matches the default climlab ebm (with instant insolation).
ebm = climlab.EBM(num_lat=16, num_lon=16, timestep=climlab.constants.seconds_per_hour, water_depth=30)
iins = climlab.radiation.insolation.InstantInsolation(domains=ebm.Ts.domain,timestep=ebm.time['timestep'])
ebm.add_subprocess('insolation', iins)
ebm.subprocess('SW').insolation = iins.insolation
It is quite noticeable that the temperatures produced by this model are too large at the equator and too small at the poles. I've attached a plot below.
This needs to get fixed before this is useable.