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

Fixed vertical motion, added pressure coordinates 

- Got some flies in this house
- Got some flies in this house
This commit is contained in:
Simon Clark 2020-11-04 20:55:18 +00:00
parent c899d6bf7a
commit bdb7b5fa17
5 changed files with 306 additions and 279 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

52
claude_low_level_library.pyx
View File

@ -7,6 +7,7 @@ cimport cython
ctypedef np.float64_t DTYPE_f
cdef float inv_180 = np.pi/180
cdef float inv_90 = np.pi/90
cdef DTYPE_f sigma = 5.67E-8
# define various useful differential functions:
# gradient of scalar field a in the local x direction at point i,j
@ -33,30 +34,19 @@ def scalar_gradient_y_2D(np.ndarray a,DTYPE_f dy,np.int_t nlat,np.int_t i,np.int
else:
return (a[i+1,j]-a[i-1,j])/dy
def scalar_gradient_z(np.ndarray a,np.ndarray dz,np.int_t i,np.int_t j,np.int_t k):
cdef np.int_t nlevels = len(dz)
def scalar_gradient_z_1D(np.ndarray a,np.ndarray pressure_levels,np.int_t k):
cdef np.int_t nlevels = len(pressure_levels)
if k == 0:
return (a[i,j,k+1]-a[i,j,k])/dz[k]
return -(a[k+1]-a[k])/(pressure_levels[k+1]-pressure_levels[k])
elif k == nlevels-1:
return (a[i,j,k]-a[i,j,k-1])/dz[k]
return -(a[k]-a[k-1])/(pressure_levels[k]-pressure_levels[k-1])
else:
return (a[i,j,k+1]-a[i,j,k-1])/(2*dz[k])
def scalar_gradient_z_1D(np.ndarray a,np.ndarray dz,np.int_t k):
cdef np.int_t nlevels = len(dz)
if k == 0:
return (a[k+1]-a[k])/dz[k]
elif k == nlevels-1:
return (a[k]-a[k-1])/dz[k]
else:
return (a[k+1]-a[k-1])/(2*dz[k])
return -(a[k+1]-a[k-1])/(pressure_levels[k+1]-pressure_levels[k-1])
def surface_optical_depth(DTYPE_f lat):
# cdef DTYPE_f inv_90
return 4 + np.cos(lat*inv_90)*2
return 4# + np.cos(lat*inv_90)*2
def thermal_radiation(DTYPE_f a):
cdef DTYPE_f sigma = 5.67E-8
return sigma*(a**4)
# power incident on (lat,lon) at time t
@ -87,24 +77,24 @@ def solar(DTYPE_f insolation,DTYPE_f lat,DTYPE_f lon,np.int_t t,DTYPE_f day,DT
def profile(np.ndarray a):
return np.mean(np.mean(a,axis=0),axis=0)
def t_to_theta(np.ndarray temperature_atmos, np.ndarray air_pressure):
def t_to_theta(np.ndarray temperature_atmos, np.ndarray pressure_levels):
cdef np.ndarray output = np.zeros_like(temperature_atmos)
cdef np.int_t i,j,k
cdef np.int_t 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
inv_p0 = 1/pressure_levels[0]
for k in range(len(pressure_levels)):
output[:,:,k] = temperature_atmos[:,:,k]*(pressure_levels[k]*inv_p0)**(-0.286)
return output
def theta_to_t(np.ndarray theta, np.ndarray air_pressure):
def theta_to_t(np.ndarray theta, np.ndarray pressure_levels):
cdef np.ndarray output = np.zeros_like(theta)
cdef np.int_t i,j,k
cdef np.int_t 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
inv_p0 = 1/pressure_levels[0]
for k in range(len(pressure_levels)):
output[:,:,k] = theta[:,:,k]*(pressure_levels[k]*inv_p0)**(0.286)
return output

2
claude_setup.py
View File

@ -4,5 +4,5 @@ import numpy
setup(
include_dirs=[numpy.get_include()],
ext_modules = cythonize("*.pyx")
ext_modules = cythonize("*.pyx", compiler_directives={'language_level' : "3"})
)

90
claude_top_level_library.pyx
View File

@ -37,7 +37,7 @@ def laplacian_3D(np.ndarray a,np.ndarray dx,DTYPE_f dy,np.ndarray dz):
return output
# divergence of (a*u) where a is a scalar field and u is the atmospheric velocity field
def divergence_with_scalar(np.ndarray a,np.ndarray u,np.ndarray v,np.ndarray w,np.ndarray dx,DTYPE_f dy,np.ndarray dz):
def divergence_with_scalar(np.ndarray a,np.ndarray u,np.ndarray v,np.ndarray w,np.ndarray dx,DTYPE_f dy,np.ndarray pressure_levels):
cdef np.ndarray output = np.zeros_like(a)
cdef np.ndarray au, av, aw
cdef np.int_t nlat, nlon, nlevels, i, j, k
@ -53,61 +53,39 @@ 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) + 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) + low_level.scalar_gradient_z_1D(aw[i,j,:],pressure_levels,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) + 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):
def radiation_calculation(np.ndarray temperature_world, np.ndarray potential_temperature, np.ndarray pressure_levels, np.ndarray heat_capacity_earth, np.ndarray albedo, DTYPE_f insolation, np.ndarray lat, np.ndarray lon, np.int_t t, np.int_t dt, DTYPE_f day, DTYPE_f year, DTYPE_f axial_tilt):
# calculate change in temperature of ground and atmosphere due to radiative imbalance
cdef np.int_t nlat,nlon,nlevels,i,j
cdef np.int_t nlat,nlon,nlevels,i,j,k
cdef DTYPE_f fl = 0.1
cdef np.ndarray upward_radiation,downward_radiation,optical_depth,Q, pressure_profile, density_profile
cdef DTYPE_f sun_lat
cdef np.ndarray upward_radiation,downward_radiation,optical_depth,Q,temperature_atmos
cdef DTYPE_f sun_lat, inv_day
inv_day = 1/(24*60*60)
temperature_atmos = low_level.theta_to_t(potential_temperature,pressure_levels)
nlat = temperature_atmos.shape[0]
nlon = temperature_atmos.shape[1]
nlevels = temperature_atmos.shape[2]
upward_radiation = np.zeros(nlevels)
downward_radiation = np.zeros(nlevels)
optical_depth = np.zeros(nlevels)
Q = np.zeros(nlevels)
for i in range(nlat):
sun_lat = low_level.surface_optical_depth(lat[i])
optical_depth = sun_lat*(fl*(pressure_levels/pressure_levels[0]) + (1-fl)*(pressure_levels/pressure_levels[0])**4)
for j in range(nlon):
# calculate optical depth
pressure_profile = air_pressure[i,j,:]
density_profile = air_density[i,j,:]
optical_depth = sun_lat*(fl*(pressure_profile/pressure_profile[0]) + (1-fl)*(pressure_profile/pressure_profile[0])**4)
# calculate upward longwave flux, bc is thermal radiation at surface
upward_radiation[0] = low_level.thermal_radiation(temperature_world[i,j])
for k in np.arange(1,nlevels):
upward_radiation[k] = (upward_radiation[k-1] - (optical_depth[k]-optical_depth[k-1])*(low_level.thermal_radiation(temperature_atmos[i,j,k])))/(1+optical_depth[k-1]-optical_depth[k])
upward_radiation[k] = (upward_radiation[k-1] - (optical_depth[k]-optical_depth[k-1])*(low_level.thermal_radiation(temperature_atmos[i,j,k])))/(1 + optical_depth[k-1] - optical_depth[k])
# calculate downward longwave flux, bc is zero at TOA (in model)
downward_radiation[-1] = 0
@ -116,42 +94,46 @@ def radiation_calculation(np.ndarray temperature_world, np.ndarray temperature_a
# gradient of difference provides heating at each level
for k in np.arange(nlevels):
Q[k] = -low_level.scalar_gradient_z_1D(upward_radiation-downward_radiation,dz,k)/(1E3*density_profile[k])
# make sure model does not have a higher top than 50km!!
Q[k] = -287*temperature_atmos[i,j,k]*low_level.scalar_gradient_z_1D(upward_radiation-downward_radiation,pressure_levels,k)/(1000*pressure_levels[k])
# approximate SW heating of ozone
if heights[k] > 20E3:
Q[k] += low_level.solar(5,lat[i],lon[j],t,day, year, axial_tilt)*((((heights[k]-20E3)/1E3)**2)/(30**2))/(24*60*60)
if pressure_levels[k] < 400*100:
Q[k] += low_level.solar(75,lat[i],lon[j],t,day, year, axial_tilt)*inv_day*(100/pressure_levels[k])
temperature_atmos[i,j,:] += Q*dt
# update surface temperature with shortwave radiation flux
temperature_world[i,j] += dt*((1-albedo[i,j])*(low_level.solar(insolation,lat[i],lon[j],t, day, year, axial_tilt) + downward_radiation[0]) - upward_radiation[0])/heat_capacity_earth[i,j]
return temperature_world, temperature_atmos
return temperature_world, low_level.t_to_theta(temperature_atmos,pressure_levels)
def velocity_calculation(np.ndarray u,np.ndarray v,np.ndarray w,np.ndarray air_pressure,np.ndarray ref_pressure_profile,np.ndarray air_density,np.ndarray coriolis,DTYPE_f gravity,np.ndarray dx,DTYPE_f dy,np.ndarray dz,DTYPE_f dt):
def velocity_calculation(np.ndarray u,np.ndarray v,np.ndarray w,np.ndarray pressure_levels,np.ndarray geopotential,np.ndarray potential_temperature,np.ndarray coriolis,DTYPE_f gravity,np.ndarray dx,DTYPE_f dy,DTYPE_f dt):
# introduce temporary arrays to update velocity in the atmosphere
cdef np.ndarray u_temp = np.zeros_like(u)
cdef np.ndarray v_temp = np.zeros_like(v)
cdef np.ndarray w_temp = np.zeros_like(u)
cdef np.ndarray temperature_atmos
cdef np.int_t nlat,nlon,nlevels,i,j,k
cdef DTYPE_f inv_density,inv_dt_grav
nlat = air_pressure.shape[0]
nlon = air_pressure.shape[1]
nlevels = air_pressure.shape[2]
inv_dt_grav = 1/(dt*gravity)
nlat = geopotential.shape[0]
nlon = geopotential.shape[1]
nlevels = len(pressure_levels)
# calculate acceleration of atmosphere using primitive equations on beta-plane
for i in np.arange(2,nlat-2).tolist():
for j in range(nlon):
for k in range(nlevels):
inv_density = 1/air_density[i,j,k]
u_temp[i,j,k] += dt*( -u[i,j,k]*low_level.scalar_gradient_x(u,dx,nlon,i,j,k) - v[i,j,k]*low_level.scalar_gradient_y(u,dy,nlat,i,j,k) - w[i,j,k]*low_level.scalar_gradient_z(u,dz,i,j,k) + coriolis[i]*v[i,j,k] - low_level.scalar_gradient_x(air_pressure,dx,nlon,i,j,k)*inv_density - 1E-6*u[i,j,k])
v_temp[i,j,k] += dt*( -u[i,j,k]*low_level.scalar_gradient_x(v,dx,nlon,i,j,k) - v[i,j,k]*low_level.scalar_gradient_y(v,dy,nlat,i,j,k) - w[i,j,k]*low_level.scalar_gradient_z(v,dz,i,j,k) - coriolis[i]*u[i,j,k] - low_level.scalar_gradient_y(air_pressure,dy,nlat,i,j,k)*inv_density - 1E-6*v[i,j,k])
# w_temp[i,j,k] += dt*( -u[i,j,k]*low_level.scalar_gradient_x(w,dx,nlon,i,j,k) - v[i,j,k]*low_level.scalar_gradient_y(w,dy,nlat,i,j,k) - w[i,j,k]*low_level.scalar_gradient_z(w,dz,i,j,k) - inv_density*low_level.scalar_gradient_z(air_pressure,dz,i,j,k) - gravity - 1E-6*w[i,j,k])
w_temp[i,j,k] += dt*( - inv_density*low_level.scalar_gradient_z_1D((air_pressure[i,j,:]-ref_pressure_profile),dz,k) - gravity - 1E-6*w[i,j,k])
u_temp[i,j,k] += dt*( -u[i,j,k]*low_level.scalar_gradient_x(u,dx,nlon,i,j,k) - v[i,j,k]*low_level.scalar_gradient_y(u,dy,nlat,i,j,k) - w[i,j,k]*low_level.scalar_gradient_z_1D(u[i,j,:],pressure_levels,k) + coriolis[i]*v[i,j,k] - low_level.scalar_gradient_x(geopotential,dx,nlon,i,j,k) - 1E-5*u[i,j,k])
v_temp[i,j,k] += dt*( -u[i,j,k]*low_level.scalar_gradient_x(v,dx,nlon,i,j,k) - v[i,j,k]*low_level.scalar_gradient_y(v,dy,nlat,i,j,k) - w[i,j,k]*low_level.scalar_gradient_z_1D(v[i,j,:],pressure_levels,k) - coriolis[i]*u[i,j,k] - low_level.scalar_gradient_y(geopotential,dy,nlat,i,j,k) - 1E-5*v[i,j,k])
temperature_atmos = low_level.theta_to_t(potential_temperature,pressure_levels)
for i in np.arange(2,nlat-2).tolist():
for j in range(nlon):
for k in np.arange(1,nlevels).tolist():
w_temp[i,j,k] = w_temp[i,j,k-1] - (pressure_levels[k]-pressure_levels[k-1])*pressure_levels[k]*gravity*( low_level.scalar_gradient_x(u,dx,nlon,i,j,k) + low_level.scalar_gradient_y(v,dy,nlat,i,j,k) )/(287*temperature_atmos[i,j,k])
u += u_temp
v += v_temp
@ -161,6 +143,10 @@ def velocity_calculation(np.ndarray u,np.ndarray v,np.ndarray w,np.ndarray air_p
u[:,:,0] *= 0.8
v[:,:,0] *= 0.8
# try to eliminate problems at top boundary
u[:,:,-1] *= 0.5
v[:,:,-1] *= 0.5
return u,v,w
def smoothing_3D(np.ndarray a,DTYPE_f smooth_parameter, DTYPE_f vert_smooth_parameter=0.5):

54
standard_atmosphere.txt
View File

@ -1,41 +1,13 @@
0 288.1 1.225E+0
2 275.2 1.007E+0
4 262.2 8.193E-1
6 249.2 6.601E-1
8 236.2 5.258E-1
10 223.3 4.135E-1
12 216.6 3.119E-1
14 216.6 2.279E-1
16 216.6 1.665E-1
18 216.6 1.216E-1
20 216.6 8.891E-2
22 218.6 6.451E-2
24 220.6 4.694E-2
26 222.5 3.426E-2
28 224.5 2.508E-2
30 226.5 1.841E-2
32 228.5 1.355E-2
34 233.7 9.887E-3
36 239.3 7.257E-3
38 244.8 5.366E-3
40 250.4 3.995E-3
42 255.9 2.995E-3
44 261.4 2.259E-3
46 266.9 1.714E-3
48 270.6 1.317E-3
50 270.6 1.027E-3
52 269.0 8.055E-4
54 263.5 6.389E-4
56 258.0 5.044E-4
58 252.5 3.962E-4
60 247.0 3.096E-4
62 241.5 2.407E-4
64 236.0 1.860E-4
66 230.5 1.429E-4
68 225.1 1.091E-4
70 219.6 8.281E-5
72 214.3 6.236E-5
74 210.3 4.637E-5
76 206.4 3.430E-5
78 202.5 2.523E-5
80 198.6 1.845E-5
0 288.1 1.225E+0 101300
2 275.2 1.007E+0 79500
4 262.2 8.193E-1 61700
6 249.2 6.601E-1 47200
8 236.2 5.258E-1 35600
10 223.3 4.135E-1 26500
20 216.6 8.891E-2 5500
30 226.5 1.841E-2 1200
40 250.4 3.995E-3 290
50 270.6 1.027E-3 80
60 247.0 3.096E-4 20
70 219.6 8.281E-5 5
80 198.6 1.845E-5 0.1

387
toy_model.py
View File

@ -8,38 +8,45 @@ import matplotlib.pyplot as plt
import time, sys, pickle
import claude_low_level_library as low_level
import claude_top_level_library as top_level
from scipy.interpolate import interp2d
# from twitch import prime_sub
######## CONTROL ########
day = 60*60*24 # define length of day (used for calculating Coriolis as well) (s)
resolution = 5 # how many degrees between latitude and longitude gridpoints
nlevels = 20 # how many vertical layers in the atmosphere
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
year = 365*day # length of year (s)
pressure_levels = np.array([1000,950,900,800,700,600,500,400,350,300,250,200,150,100,75,50,25,10,5,2,1])
pressure_levels *= 100
nlevels = len(pressure_levels)
top = nlevels
dt_spinup = 60*137
dt_main = 60*0.1
spinup_length = 10*day
dt_main = 60*3.5
spinup_length = 1*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
advection_boundary = 5 # how many gridpoints away from poles to apply advection
smoothing_parameter_t = 0.6
smoothing_parameter_u = 0.6
smoothing_parameter_v = 0.6
smoothing_parameter_w = 0.2
smoothing_parameter_w = 0.4
save = False # save current state to file?
load = True # load initial state from file?
load = False # load initial state from file?
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)
diagnostic = True # display raw fields for diagnostic purposes
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)
###########################
@ -49,46 +56,40 @@ lon = np.arange(0,360,resolution)
nlat = len(lat)
nlon = len(lon)
lon_plot, lat_plot = np.meshgrid(lon, lat)
heights = np.arange(0,top,top/nlevels)
heights_plot, lat_z_plot = np.meshgrid(lat,heights)
heights_plot, lat_z_plot = np.meshgrid(lat,pressure_levels/100)
# initialise arrays for various physical fields
temperature_world = np.zeros((nlat,nlon)) + 290
temperature_atmos = np.zeros((nlat,nlon,nlevels))
air_pressure = np.zeros_like(temperature_atmos)
u = np.zeros_like(temperature_atmos)
v = np.zeros_like(temperature_atmos)
w = np.zeros_like(temperature_atmos)
air_density = np.zeros_like(temperature_atmos)
potential_temperature = np.zeros((nlat,nlon,nlevels))
u = np.zeros_like(potential_temperature)
v = np.zeros_like(potential_temperature)
w = np.zeros_like(potential_temperature)
atmosp_addition = np.zeros_like(potential_temperature)
# #######################
##########################
# read temperature and density in from standard atmosphere
f = open("standard_atmosphere.txt", "r")
standard_height = []
standard_temp = []
standard_density = []
standard_pressure = []
for x in f:
h, t, r = x.split()
standard_height.append(float(h))
h, t, r, p = x.split()
standard_temp.append(float(t))
standard_density.append(float(r))
standard_pressure.append(float(p))
f.close()
density_profile = np.interp(x=heights/1E3,xp=standard_height,fp=standard_density)
temp_profile = np.interp(x=heights/1E3,xp=standard_height,fp=standard_temp)
# density_profile = np.interp(x=heights/1E3,xp=standard_height,fp=standard_density)
temp_profile = np.interp(x=pressure_levels[::-1],xp=standard_pressure[::-1],fp=standard_temp[::-1])[::-1]
for k in range(nlevels):
air_density[:,:,k] = density_profile[k]
temperature_atmos[:,:,k] = temp_profile[k]
potential_temperature[:,:,k] = temp_profile[k]
###########################
potential_temperature = low_level.t_to_theta(potential_temperature,pressure_levels)
geopotential = np.zeros_like(potential_temperature)
# weight_above = np.interp(x=heights/1E3,xp=standard_height,fp=standard_density)
# top_index = np.argmax(np.array(standard_height) >= top/1E3)
# if standard_height[top_index] == top/1E3:
# weight_above = np.trapz(np.interp(x=standard_height[top_index:],xp=standard_height,fp=standard_density),standard_height[top_index:])*gravity*1E3
# else:
# weight_above = np.trapz(np.interp(x=np.insert(standard_height[top_index:], 0, top/1E3),xp=standard_height,fp=standard_density),np.insert(standard_height[top_index:], 0, top/1E3))*gravity*1E3
sigma = np.zeros_like(pressure_levels)
kappa = 287/1000
for i in range(len(sigma)):
sigma[i] = 1E3*(pressure_levels[i]/pressure_levels[0])**kappa
###########################
@ -103,9 +104,6 @@ for i in range(nlat):
specific_gas = 287
thermal_diffusivity_roc = 1.5E-6
ref_pressure_profile = density_profile*specific_gas*temp_profile
air_pressure = air_density*specific_gas*temperature_atmos
# define planet size and various geometric constants
circumference = 2*np.pi*planet_radius
circle = np.pi*planet_radius**2
@ -119,60 +117,130 @@ angular_speed = 2*np.pi/day
for i in range(nlat):
dx[i] = dy*np.cos(lat[i]*np.pi/180)
coriolis[i] = angular_speed*np.sin(lat[i]*np.pi/180)
dz = np.zeros(nlevels)
for k in range(nlevels-1): dz[k] = heights[k+1] - heights[k]
dz[-1] = dz[-2]
#################### SHOW TIME ####################
# pole_lower_limit = 4
# polar_grid_resolution = dx[-pole_lower_limit]
# size_of_grid = planet_radius*np.cos(lat[-pole_lower_limit]*np.pi/180)
# x_coords = np.arange(-size_of_grid,size_of_grid,polar_grid_resolution)
# y_coords = np.arange(-size_of_grid,size_of_grid,polar_grid_resolution)
# xx,yy = np.meshgrid(x_coords,y_coords)
# polar_coords = []
# for i in range(xx.shape[0]):
# for j in range(xx.shape[1]):
# lat_point = np.arcsin((xx[i,j]**2 + yy[i,j]**2)**0.5/planet_radius)
# lon_point = np.arctan2(yy[i,j],xx[i,j])
# polar_coords.append((lat_point,lon_point))
# polar_data = np.zeros((pole_lower_limit,nlon))
# for i in range(pole_lower_limit):
# polar_data[i,:] = lat[i]
# f = interp2d(lat[:pole_lower_limit]*np.pi/180, lon*np.pi/180, np.transpose(polar_data))
# print(lat[:pole_lower_limit]*np.pi/180)
# print(polar_coords)
# polar_plane = np.zeros_like(xx)
# index = 0
# for i in range(xx.shape[0]):
# for j in range(xx.shape[1]):
# polar_plane[i,j] = f(polar_coords[index][0],polar_coords[index][1])
# index += 1
# print(polar_plane)
# plt.imshow(polar_plane)
# plt.show()
# sys.exit()
#######################################################################################################################################################################################################################
# INITIATE TIME
t = 0
if load:
# load in previous save file
temperature_atmos,temperature_world,u,v,w,t,air_pressure,air_density,albedo = pickle.load(open("save_file.p","rb"))
potential_temperature,temperature_world,u,v,w,t,albedo = pickle.load(open("save_file.p","rb"))
if plot:
# set up plot
f, ax = plt.subplots(2,figsize=(9,9))
f.canvas.set_window_title('CLAuDE')
test = ax[0].contourf(lon_plot, lat_plot, temperature_world, cmap='seismic')
ax[0].streamplot(lon_plot, lat_plot, u[:,:,0], v[:,:,0], color='white',density=1)
ax[1].contourf(heights_plot, lat_z_plot, np.transpose(np.mean(temperature_atmos,axis=1)), cmap='seismic',levels=15)
ax[1].contour(heights_plot,lat_z_plot, np.transpose(np.mean(u,axis=1)), colors='white',levels=20,linewidths=1,alpha=0.8)
ax[1].streamplot(heights_plot, lat_z_plot, np.transpose(np.mean(v,axis=1)),np.transpose(np.mean(10*w,axis=1)),color='black',density=0.75)
plt.subplots_adjust(left=0.1, right=0.75)
ax[0].set_title('Ground temperature')
ax[0].set_xlim(lon.min(),lon.max())
ax[1].set_title('Atmosphere temperature')
ax[1].set_xlim(lat.min(),lat.max())
ax[1].set_ylim(heights.min(),heights.max())
cbar_ax = f.add_axes([0.85, 0.15, 0.05, 0.7])
if not diagnostic:
# set up plot
f, ax = plt.subplots(2,figsize=(9,9))
f.canvas.set_window_title('CLAuDE')
test = ax[0].contourf(lon_plot, lat_plot, temperature_world, cmap='seismic')
ax[0].streamplot(lon_plot, lat_plot, u[:,:,0], v[:,:,0], color='white',density=1)
ax[1].contourf(heights_plot, lat_z_plot, np.transpose(np.mean(low_level.theta_to_t(potential_temperature,pressure_levels),axis=1))[:top,:], cmap='seismic',levels=15)
ax[1].contour(heights_plot,lat_z_plot, np.transpose(np.mean(u,axis=1))[:top,:], colors='white',levels=20,linewidths=1,alpha=0.8)
ax[1].quiver(heights_plot, lat_z_plot, np.transpose(np.mean(v,axis=1))[:top,:],np.transpose(np.mean(10*w,axis=1))[:top,:],color='black')
plt.subplots_adjust(left=0.1, right=0.75)
ax[0].set_title('Surface temperature')
ax[0].set_xlim(lon.min(),lon.max())
ax[1].set_title('Atmosphere temperature')
ax[1].set_xlim(lat.min(),lat.max())
ax[1].set_ylim((pressure_levels.max()/100,pressure_levels[:top].min()/100))
ax[1].set_yscale('log')
ax[1].set_ylabel('Pressure (hPa)')
ax[1].set_xlabel('Latitude')
cbar_ax = f.add_axes([0.85, 0.15, 0.05, 0.7])
f.colorbar(test, cax=cbar_ax)
cbar_ax.set_title('Temperature (K)')
f.suptitle( 'Time ' + str(round(t/day,2)) + ' days' )
f.colorbar(test, cax=cbar_ax)
cbar_ax.set_title('Temperature (K)')
f.suptitle( 'Time ' + str(round(t/day,2)) + ' days' )
if level_plots:
# allow for live updating as calculations take place
if level_plots:
level_divisions = int(np.floor(nlevels/nplots))
level_plots_levels = range(nlevels)[::level_divisions]
g, bx = plt.subplots(nplots,figsize=(9,8),sharex=True)
g.canvas.set_window_title('CLAuDE atmospheric levels')
for k, z in zip(range(nplots), level_plots_levels):
z += 1
bx[k].contourf(lon_plot, lat_plot, temperature_atmos[:,:,z], cmap='seismic')
bx[k].set_title(str(heights[z])+' km')
bx[k].set_ylabel('Latitude')
bx[-1].set_xlabel('Longitude')
level_divisions = int(np.floor(nlevels/nplots))
level_plots_levels = range(nlevels)[::level_divisions][::-1]
g, bx = plt.subplots(nplots,figsize=(9,8),sharex=True)
g.canvas.set_window_title('CLAuDE pressure levels')
for k, z in zip(range(nplots), level_plots_levels):
z += 1
bx[k].contourf(lon_plot, lat_plot, potential_temperature[:,:,z], cmap='seismic')
bx[k].set_title(str(pressure_levels[z]/100)+' hPa')
bx[k].set_ylabel('Latitude')
bx[-1].set_xlabel('Longitude')
else:
# set up plot
f, ax = plt.subplots(2,2,figsize=(9,9))
f.canvas.set_window_title('CLAuDE')
ax[0,0].contourf(heights_plot, lat_z_plot, np.transpose(np.mean(u,axis=1))[:top,:], cmap='seismic')
ax[0,0].set_title('u')
ax[0,1].contourf(heights_plot, lat_z_plot, np.transpose(np.mean(v,axis=1))[:top,:], cmap='seismic')
ax[0,1].set_title('v')
ax[1,0].contourf(heights_plot, lat_z_plot, np.transpose(np.mean(w,axis=1))[:top,:], cmap='seismic')
ax[1,0].set_title('w')
ax[1,1].contourf(heights_plot, lat_z_plot, np.transpose(np.mean(atmosp_addition,axis=1))[:top,:], cmap='seismic')
ax[1,1].set_title('atmosp_addition')
for axis in ax.ravel():
axis.set_ylim((pressure_levels.max()/100,pressure_levels[:top].min()/100))
axis.set_yscale('log')
f.suptitle( 'Time ' + str(round(t/day,2)) + ' days' )
plt.ion()
plt.show()
plt.pause(2)
if not diagnostic:
ax[0].cla()
ax[1].cla()
if level_plots:
for k in range(nplots):
bx[k].cla()
else:
ax[0,0].cla()
ax[0,1].cla()
ax[1,0].cla()
ax[1,1].cla()
while True:
initial_time = time.time()
@ -189,24 +257,26 @@ while True:
print('T: ',round(temperature_world.max()-273.15,1),' - ',round(temperature_world.min()-273.15,1),' C')
print('U: ',round(u.max(),2),' - ',round(u.min(),2),' V: ',round(v.max(),2),' - ',round(v.min(),2),' W: ',round(w.max(),2),' - ',round(w.min(),2))
# before_radiation = time.time()
temperature_world, temperature_atmos = top_level.radiation_calculation(temperature_world, temperature_atmos, air_pressure, air_density, heat_capacity_earth, albedo, insolation, lat, lon, heights, dz, t, dt, day, year, axial_tilt)
temperature_atmos = top_level.smoothing_3D(temperature_atmos,smoothing_parameter_t)
# time_taken = float(round(time.time() - before_radiation,3))
# print('Radiation: ',str(time_taken),'s')
before_radiation = time.time()
temperature_world, potential_temperature = top_level.radiation_calculation(temperature_world, potential_temperature, pressure_levels, heat_capacity_earth, albedo, insolation, lat, lon, t, dt, day, year, axial_tilt)
potential_temperature = top_level.smoothing_3D(potential_temperature,smoothing_parameter_t)
time_taken = float(round(time.time() - before_radiation,3))
print('Radiation: ',str(time_taken),'s')
# update air pressure
old_pressure = np.copy(air_pressure)
air_pressure = air_density*specific_gas*temperature_atmos
# update geopotential field
for k in np.arange(1,nlevels):
geopotential[:,:,k] = geopotential[:,:,k-1] - potential_temperature[:,:,k]*(sigma[k]-sigma[k-1])
geopotential = top_level.smoothing_3D(geopotential,smoothing_parameter_t)
if velocity:
before_velocity = time.time()
u,v,w = top_level.velocity_calculation(u,v,w,air_pressure,ref_pressure_profile,air_density,coriolis,gravity,dx,dy,dz,dt)
u,v,w = top_level.velocity_calculation(u,v,w,pressure_levels,geopotential,potential_temperature,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,0.1)
# boundary shite
u[(advection_boundary,-advection_boundary-1),:,:] *= 0.5
v[(advection_boundary,-advection_boundary-1),:,:] *= 0.5
w[(advection_boundary,-advection_boundary-1),:,:] *= 0.5
@ -214,97 +284,106 @@ while True:
u[:advection_boundary,:,:] = 0
v[:advection_boundary,:,:] = 0
w[:advection_boundary,:,:] = 0
u[-advection_boundary:,:,:] = 0
v[-advection_boundary:,:,:] = 0
w[-advection_boundary:,:,:] = 0
# time_taken = float(round(time.time() - before_velocity,3))
# print('Velocity: ',str(time_taken),'s')
w[:,:,-1] *= 0
w[:,:,-2] *= 0.1
w[:,:,-3] *= 0.5
time_taken = float(round(time.time() - before_velocity,3))
print('Velocity: ',str(time_taken),'s')
if advection:
# before_advection = time.time()
before_advection = time.time()
# 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,:,:]
# as density is now variable, allow for mass advection
density_addition = dt*top_level.divergence_with_scalar(air_density,u,v,w,dx,dy,dz)
air_density[advection_boundary:-advection_boundary,:,:] -= density_addition[advection_boundary:-advection_boundary,:,:]
air_density[(advection_boundary-1),:,:] -= 0.5*density_addition[advection_boundary-1,:,:]
air_density[-advection_boundary,:,:] -= 0.5*density_addition[-advection_boundary,:,:]
atmosp_addition = dt*top_level.divergence_with_scalar(potential_temperature,u,v,w,dx,dy,pressure_levels)
atmosp_addition[(-advection_boundary,advection_boundary-1),:,:] *= 0.5
atmosp_addition[:advection_boundary,:,:] *= 0
atmosp_addition[-advection_boundary:,:,:] *= 0
atmosp_addition[:,:,-1] *= 0
atmosp_addition[:,:,-2] *= 0.5
potential_temperature -= atmosp_addition
# 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')
time_taken = float(round(time.time() - before_advection,3))
print('Advection: ',str(time_taken),'s')
# before_plot = time.time()
if plot:
tropopause_height = np.zeros(nlat)
for i in range(nlat):
k = 2
if low_level.scalar_gradient_z_1D(np.mean(temperature_atmos[i,:,:],axis=0),dz,k) > 0:
k += 1
while low_level.scalar_gradient_z_1D(np.mean(temperature_atmos[i,:,:],axis=0),dz,k) < 0:
k += 1
else:
while low_level.scalar_gradient_z_1D(np.mean(temperature_atmos[i,:,:],axis=0),dz,k) < 0:
k += 1
tropopause_height[i] = heights[k]
# update plot
test = ax[0].contourf(lon_plot, lat_plot, temperature_world, cmap='seismic',levels=15)
if velocity: ax[0].streamplot(lon_plot, lat_plot, u[:,:,0], v[:,:,0], color='white',density=1)
ax[0].set_title('$\it{Ground} \quad \it{temperature}$')
ax[0].set_xlim((lon.min(),lon.max()))
ax[0].set_ylim((lat.min(),lat.max()))
ax[0].set_ylabel('Latitude')
ax[0].axhline(y=0,color='black',alpha=0.3)
ax[0].set_xlabel('Longitude')
if not diagnostic:
ax[0].contourf(lon_plot, lat_plot, temperature_world, cmap='seismic',levels=15)
if velocity:
ax[0].streamplot(lon_plot, lat_plot, u[:,:,0], v[:,:,0], color='white',density=0.75)
ax[0].set_title('$\it{Ground} \quad \it{temperature}$')
ax[0].set_xlim((lon.min(),lon.max()))
ax[0].set_ylim((lat.min(),lat.max()))
ax[0].set_ylabel('Latitude')
ax[0].axhline(y=0,color='black',alpha=0.3)
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)
ax[1].streamplot(heights_plot, lat_z_plot, np.transpose(np.mean(v,axis=1)),np.transpose(np.mean(10*w,axis=1)),color='black',density=0.75)
ax[1].set_title('$\it{Atmospheric} \quad \it{temperature}$')
ax[1].set_xlim((-90,90))
ax[1].set_ylim((0,heights.max()))
ax[1].set_ylabel('Height (m)')
ax[1].set_xlabel('Latitude')
f.colorbar(test, cax=cbar_ax)
cbar_ax.set_title('Temperature (K)')
f.suptitle( 'Time ' + str(round(t/day,2)) + ' days' )
test = ax[1].contourf(heights_plot, lat_z_plot, np.transpose(np.mean(low_level.theta_to_t(potential_temperature,pressure_levels),axis=1))[:top,:], cmap='seismic',levels=15)
# test = ax[1].contourf(heights_plot, lat_z_plot, np.transpose(np.mean(potential_temperature,axis=1)), cmap='seismic',levels=15)
# test = ax[1].contourf(heights_plot, lat_z_plot, np.transpose(np.mean(v,axis=1))[:top,:], cmap='seismic',levels=15)
if velocity:
ax[1].contour(heights_plot,lat_z_plot, np.transpose(np.mean(u,axis=1))[:top,:], colors='white',levels=20,linewidths=1,alpha=0.8)
ax[1].quiver(heights_plot, lat_z_plot, np.transpose(np.mean(v,axis=1))[:top,:],np.transpose(np.mean(10*w,axis=1))[:top,:],color='black')
ax[1].set_title('$\it{Atmospheric} \quad \it{temperature}$')
ax[1].set_xlim((-90,90))
ax[1].set_ylim((pressure_levels.max()/100,pressure_levels[:top].min()/100))
ax[1].set_ylabel('Pressure (hPa)')
ax[1].set_xlabel('Latitude')
ax[1].set_yscale('log')
f.colorbar(test, cax=cbar_ax)
cbar_ax.set_title('Temperature (K)')
f.suptitle( 'Time ' + str(round(t/day,2)) + ' days' )
if level_plots:
quiver_padding = int(50/resolution)
skip=(slice(None,None,2),slice(None,None,2))
for k, z in zip(range(nplots), level_plots_levels):
z += 1
bx[k].contourf(lon_plot, lat_plot, potential_temperature[:,:,z], cmap='seismic',levels=15)
bx[k].streamplot(lon_plot, lat_plot, u[:,:,z], v[:,:,z], color='white',density=1.5)
bx[k].set_title(str(round(pressure_levels[z]/100))+' hPa')
bx[k].set_ylabel('Latitude')
bx[k].set_xlim((lon.min(),lon.max()))
bx[k].set_ylim((lat.min(),lat.max()))
bx[-1].set_xlabel('Longitude')
else:
ax[0,0].contourf(heights_plot, lat_z_plot, np.transpose(np.mean(u,axis=1))[:top,:], cmap='seismic')
ax[0,0].set_title('u')
ax[0,1].contourf(heights_plot, lat_z_plot, np.transpose(np.mean(v,axis=1))[:top,:], cmap='seismic')
ax[0,1].set_title('v')
ax[1,0].contourf(heights_plot, lat_z_plot, np.transpose(np.mean(w,axis=1))[:top,:], cmap='seismic')
ax[1,0].set_title('w')
ax[1,1].contourf(heights_plot, lat_z_plot, np.transpose(np.mean(atmosp_addition,axis=1))[:top,:], cmap='seismic')
ax[1,1].set_title('atmosp_addition')
for axis in ax.ravel():
axis.set_ylim((pressure_levels.max()/100,pressure_levels[:top].min()/100))
axis.set_yscale('log')
f.suptitle( 'Time ' + str(round(t/day,2)) + ' days' )
if level_plots:
quiver_padding = int(50/resolution)
skip=(slice(None,None,2),slice(None,None,2))
for k, z in zip(range(nplots), level_plots_levels):
z += 1
bx[k].contourf(lon_plot, lat_plot, temperature_atmos[:,:,z], cmap='seismic',levels=15)
bx[k].streamplot(lon_plot, lat_plot, u[:,:,z], v[:,:,z], color='white',density=1.5)
bx[k].set_title(str(round(heights[z]/1E3))+' km')
bx[k].set_ylabel('Latitude')
bx[k].set_xlim((lon.min(),lon.max()))
bx[k].set_ylim((lat.min(),lat.max()))
bx[-1].set_xlabel('Longitude')
plt.pause(0.01)
ax[0].cla()
ax[1].cla()
if not diagnostic:
ax[0].cla()
ax[1].cla()
if level_plots:
for k in range(nplots):
bx[k].cla()
else:
ax[0,0].cla()
ax[0,1].cla()
ax[1,0].cla()
ax[1,1].cla()
if level_plots:
for k in range(nplots):
bx[k].cla()
# time_taken = float(round(time.time() - before_plot,3))
# print('Plotting: ',str(time_taken),'s')
@ -316,7 +395,7 @@ while True:
print('Time: ',str(time_taken),'s')
if save:
pickle.dump((temperature_atmos,temperature_world,u,v,w,t,air_pressure,air_density,albedo), open("save_file.p","wb"))
pickle.dump((potential_temperature,temperature_world,u,v,w,t,albedo), open("save_file.p","wb"))
if np.isnan(u.max()):
sys.exit()