askill
/
claude
mirror of https://github.com/Askill/claude.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Added adiabatic motion to convection 

- Converted temperature to potential temperature for thermal advection
- Added optional argument for vertical smoothing in 3D FFT function
- Small changes to plotting
This commit is contained in:
Simon Clark 2020-09-16 20:57:03 +01:00
parent a3ebf0605c
commit 663573f911
4 changed files with 72 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
claude_low_level_library.pyx
View File

@ -52,7 +52,7 @@ def scalar_gradient_z_1D(np.ndarray a,np.ndarray dz,np.int_t k):
return (a[k+1]-a[k-1])/(2*dz[k])
def surface_optical_depth(DTYPE_f lat):
cdef DTYPE_f inv_90
# cdef DTYPE_f inv_90
return 4 + np.cos(lat*inv_90)*2
def thermal_radiation(DTYPE_f a):
@ -85,4 +85,26 @@ def solar(DTYPE_f insolation,DTYPE_f lat,DTYPE_f lon,np.int_t t,DTYPE_f day,DT
return value
def profile(np.ndarray a):
return np.mean(np.mean(a,axis=0),axis=0)
return np.mean(np.mean(a,axis=0),axis=0)
def t_to_theta(np.ndarray temperature_atmos, np.ndarray air_pressure):
cdef np.ndarray output = np.zeros_like(temperature_atmos)
cdef np.int_t i,j,k
cdef DTYPE_f inv_p0
for i in range(output.shape[0]):
for j in range(output.shape[1]):
inv_p0 = 1/air_pressure[i,j,0]
for k in range(output.shape[2]):
output[i,j,k] = temperature_atmos[i,j,k]*(air_pressure[i,j,k]*inv_p0)**0.286
return output
def theta_to_t(np.ndarray theta, np.ndarray air_pressure):
cdef np.ndarray output = np.zeros_like(theta)
cdef np.int_t i,j,k
cdef DTYPE_f inv_p0
for i in range(output.shape[0]):
for j in range(output.shape[1]):
inv_p0 = 1/air_pressure[i,j,0]
for k in range(output.shape[2]):
output[i,j,k] = theta[i,j,k]*(air_pressure[i,j,k]*inv_p0)**-0.286
return output

38
claude_top_level_library.pyx
View File

@ -53,7 +53,31 @@ def divergence_with_scalar(np.ndarray a,np.ndarray u,np.ndarray v,np.ndarray w,n
for i in range(nlat):
for j in range(nlon):
for k in range(nlevels):
output[i,j,k] = low_level.scalar_gradient_x(au,dx,nlon,i,j,k) + low_level.scalar_gradient_y(av,dy,nlat,i,j,k) + 0.05*low_level.scalar_gradient_z(aw,dz,i,j,k)
output[i,j,k] = low_level.scalar_gradient_x(au,dx,nlon,i,j,k) + low_level.scalar_gradient_y(av,dy,nlat,i,j,k) + 1*low_level.scalar_gradient_z(aw,dz,i,j,k)
return output
# specifically thermal advection, converts temperature field to potential temperature field, then advects (adiabatically), then converts back to temperature
def thermal_advection(np.ndarray temperature_atmos,np.ndarray air_pressure,np.ndarray u,np.ndarray v,np.ndarray w,np.ndarray dx,DTYPE_f dy,np.ndarray dz):
cdef np.ndarray output = np.zeros_like(temperature_atmos)
cdef np.ndarray au, av, aw
cdef np.int_t nlat, nlon, nlevels, i, j, k
cdef np.ndarray theta = low_level.t_to_theta(temperature_atmos,air_pressure)
nlat = output.shape[0]
nlon = output.shape[1]
nlevels = output.shape[2]
au = theta*u
av = theta*v
aw = theta*w
for i in range(nlat):
for j in range(nlon):
for k in range(nlevels):
output[i,j,k] = low_level.scalar_gradient_x(au,dx,nlon,i,j,k) + low_level.scalar_gradient_y(av,dy,nlat,i,j,k) + 1*low_level.scalar_gradient_z(aw,dz,i,j,k)
output = low_level.theta_to_t(output,air_pressure)
return output
def radiation_calculation(np.ndarray temperature_world, np.ndarray temperature_atmos, np.ndarray air_pressure, np.ndarray air_density, np.ndarray heat_capacity_earth, np.ndarray albedo, DTYPE_f insolation, np.ndarray lat, np.ndarray lon, np.ndarray heights, np.ndarray dz, np.int_t t, np.int_t dt, DTYPE_f day, DTYPE_f year, DTYPE_f axial_tilt):
@ -141,13 +165,13 @@ def velocity_calculation(np.ndarray u,np.ndarray v,np.ndarray air_pressure,np.nd
return u,v,w
def smoothing_3D(np.ndarray a,DTYPE_f smooth_parameter):
nlat = a.shape[0]
nlon = a.shape[1]
nlevels = a.shape[2]
def smoothing_3D(np.ndarray a,DTYPE_f smooth_parameter, DTYPE_f vert_smooth_parameter=0.5):
cdef np.int_t nlat = a.shape[0]
cdef np.int_t nlon = a.shape[1]
cdef np.int_t nlevels = a.shape[2]
smooth_parameter *= 0.5
test = np.fft.fftn(a)
cdef np.ndarray test = np.fft.fftn(a)
test[int(nlat*smooth_parameter):int(nlat*(1-smooth_parameter)),:,:] = 0
test[:,int(nlon*smooth_parameter):int(nlon*(1-smooth_parameter)),:] = 0
# test[:,:int(nlevels*smooth_parameter):int(nlevels*(1-smooth_parameter))] = 0
test[:,:,int(nlevels*vert_smooth_parameter):int(nlevels*(1-vert_smooth_parameter))] = 0
return np.fft.ifftn(test).real

BIN
save_file.p
View File

Binary file not shown.

34
toy_model.py
View File

@ -18,28 +18,28 @@ top = 40E3 # top of atmosphere (m)
planet_radius = 6.4E6 # define the planet's radius (m)
insolation = 1370 # TOA radiation from star (W m^-2)
gravity = 9.81 # define surface gravity for planet (m s^-2)
axial_tilt = -23.5/2 # tilt of rotational axis w.r.t. solar plane
axial_tilt = -23.5/2 # tilt of rotational axis w.r.t. solar plane
year = 365*day # length of year (s)
dt_spinup = 60*137
dt_main = 60*7
spinup_length = 3*day
spinup_length = 7*day
###
advection = True # if you want to include advection set this to be True
advection_boundary = 3 # how many gridpoints away from poles to apply advection
smoothing_parameter_t = 0.9
smoothing_parameter_u = 0.8
smoothing_parameter_v = 0.8
smoothing_parameter_t = 0.6
smoothing_parameter_u = 0.6
smoothing_parameter_v = 0.6
smoothing_parameter_w = 0.3
save = True # save current state to file?
save = False # save current state to file?
load = True # load initial state from file?
plot = False # display plots of output?
level_plots = False # display plots of output on vertical levels?
nplots = 3 # how many levels you want to see plots of (evenly distributed through column)
plot = True # display plots of output?
level_plots = True # display plots of output on vertical levels?
nplots = 4 # how many levels you want to see plots of (evenly distributed through column)
###########################
@ -203,8 +203,8 @@ while True:
before_velocity = time.time()
u,v,w = top_level.velocity_calculation(u,v,air_pressure,old_pressure,air_density,coriolis,gravity,dx,dy,dt)
u = top_level.smoothing_3D(u,smoothing_parameter_u)
# v = top_level.smoothing_3D(v,smoothing_parameter_v)
# w = top_level.smoothing_3D(w,smoothing_parameter_w)
v = top_level.smoothing_3D(v,smoothing_parameter_v)
w = top_level.smoothing_3D(w,smoothing_parameter_w,0.3)
u[(advection_boundary,-advection_boundary-1),:,:] *= 0.5
v[(advection_boundary,-advection_boundary-1),:,:] *= 0.5
@ -221,12 +221,12 @@ while True:
# time_taken = float(round(time.time() - before_velocity,3))
# print('Velocity: ',str(time_taken),'s')
# before_advection = time.time()
if advection:
# allow for thermal advection in the atmosphere, and heat diffusion in the atmosphere and the ground
# atmosp_addition = dt*(thermal_diffusivity_air*laplacian(temperature_atmos))
# before_advection = time.time()
atmosp_addition = dt*top_level.divergence_with_scalar(temperature_atmos,u,v,w,dx,dy,dz)
# allow for thermal advection in the atmosphere, and heat diffusion in the atmosphere and the ground
atmosp_addition = dt*top_level.thermal_advection(temperature_atmos,air_pressure,u,v,w,dx,dy,dz)
temperature_atmos[advection_boundary:-advection_boundary,:,:] -= atmosp_addition[advection_boundary:-advection_boundary,:,:]
temperature_atmos[advection_boundary-1,:,:] -= 0.5*atmosp_addition[advection_boundary-1,:,:]
temperature_atmos[-advection_boundary,:,:] -= 0.5*atmosp_addition[-advection_boundary,:,:]
@ -239,7 +239,6 @@ while True:
# temperature_world -= dt*(thermal_diffusivity_roc*top_level.laplacian_2D(temperature_world,dx,dy))
# time_taken = float(round(time.time() - before_advection,3))
# print('Advection: ',str(time_taken),'s')
@ -270,6 +269,7 @@ while True:
ax[0].set_xlabel('Longitude')
ax[1].contourf(heights_plot, lat_z_plot, np.transpose(np.mean(temperature_atmos,axis=1)), cmap='seismic',levels=15)
# ax[1].contourf(heights_plot, lat_z_plot, np.transpose(np.mean(w,axis=1)), cmap='seismic',levels=15)
ax[1].plot(lat_plot,tropopause_height,color='black',linestyle='--',linewidth=3,alpha=0.5)
if velocity:
ax[1].contour(heights_plot,lat_z_plot, np.transpose(np.mean(u,axis=1)), colors='white',levels=20,linewidths=1,alpha=0.8)
@ -282,7 +282,7 @@ while True:
f.colorbar(test, cax=cbar_ax)
cbar_ax.set_title('Temperature (K)')
f.suptitle( 'Time ' + str(round(24*t/day,2)) + ' hours' )
f.suptitle( 'Time ' + str(round(t/day,2)) + ' days' )
if level_plots:
quiver_padding = int(50/resolution)