mirror of https://github.com/Askill/claude.git
Remove dead library files
The content of these files were moved to the .pyx files of the same name, so they could be statically compiled by cython for improved performance.
This commit is contained in:
anothersimulacrum
parent
f265687986
commit
d3de7b14c9
|
|
@ -1,84 +0,0 @@
|
|||
# claude low level library
|
||||
|
||||
import numpy as np
|
||||
cimport numpy as np
|
||||
|
||||
ctypedef np.float64_t 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
|
||||
def scalar_gradient_x(a,dx,nlon,i,j,k):
|
||||
return (a[i,(j+1)%nlon,k]-a[i,(j-1)%nlon,k])/dx[i]
|
||||
|
||||
def scalar_gradient_x_2D(a,dx,nlon,i,j)
|
||||
return (a[i,(j+1)%nlon]-a[i,(j-1)%nlon])/dx[i]
|
||||
|
||||
# gradient of scalar field a in the local y direction at point i,j
|
||||
def scalar_gradient_y(a,dy,nlat,i,j,k):
|
||||
if i == 0:
|
||||
return 2*(a[i+1,j]-a[i,j])/dy
|
||||
elif i == nlat-1:
|
||||
return 2*(a[i,j]-a[i-1,j])/dy
|
||||
else:
|
||||
return (a[i+1,j]-a[i-1,j])/dy
|
||||
|
||||
def scalar_gradient_y_2d(a,dy,nlat,i,j)
|
||||
if i == 0:
|
||||
return 2*(a[i+1,j]-a[i,j])/dy
|
||||
elif i == nlat-1:
|
||||
return 2*(a[i,j]-a[i-1,j])/dy
|
||||
else:
|
||||
return (a[i+1,j]-a[i-1,j])/dy
|
||||
|
||||
def scalar_gradient_z(a,dz,i,j,k):
|
||||
output = np.zeros_like(a)
|
||||
nlevels = len(dz)
|
||||
if output.ndim == 1:
|
||||
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])
|
||||
else:
|
||||
if k == 0:
|
||||
return (a[i,j,k+1]-a[i,j,k])/dz[k]
|
||||
elif k == nlevels-1:
|
||||
return (a[i,j,k]-a[i,j,k-1])/dz[k]
|
||||
else:
|
||||
return (a[i,j,k+1]-a[i,j,k-1])/(2*dz[k])
|
||||
|
||||
def surface_optical_depth(lat):
|
||||
return 4 + np.cos(lat*np.pi/90)*2.5/2
|
||||
|
||||
def thermal_radiation(a):
|
||||
return sigma*(a**4)
|
||||
|
||||
# power incident on (lat,lon) at time t
|
||||
def solar(insolation, lat, lon, t, day, year, axial_tilt):
|
||||
sun_longitude = -t % day
|
||||
sun_longitude *= 360/day
|
||||
sun_latitude = axial_tilt*np.cos(t*2*np.pi/year)
|
||||
|
||||
value = insolation*np.cos((lat-sun_latitude)*np.pi/180)
|
||||
|
||||
if value < 0:
|
||||
return 0
|
||||
else:
|
||||
|
||||
lon_diff = lon-sun_longitude
|
||||
value *= np.cos(lon_diff*np.pi/180)
|
||||
|
||||
if value < 0:
|
||||
if lat + sun_latitude > 90:
|
||||
return insolation*np.cos((lat+sun_latitude)*np.pi/180)*np.cos(lon_diff*np.pi/180)
|
||||
elif lat + sun_latitude < -90:
|
||||
return insolation*np.cos((lat+sun_latitude)*np.pi/180)*np.cos(lon_diff*np.pi/180)
|
||||
else:
|
||||
return 0
|
||||
else:
|
||||
return value
|
||||
|
||||
def profile(a):
|
||||
return np.mean(np.mean(a,axis=0),axis=0)
|
||||
|
|
@ -1,101 +0,0 @@
|
|||
# claude_top_level_library
|
||||
|
||||
import claude_low_level_library as low_level
|
||||
import numpy as np
|
||||
|
||||
# laplacian of scalar field a
|
||||
def laplacian(a):
|
||||
output = np.zeros_like(a)
|
||||
if output.ndim == 2:
|
||||
for i in np.arange(1,nlat-1):
|
||||
for j in range(nlon):
|
||||
output[i,j] = (scalar_gradient_x_2D(a,i,(j+1)%nlon) - scalar_gradient_x_2D(a,i,(j-1)%nlon))/dx[i] + (scalar_gradient_y_2D(a,i+1,j) - scalar_gradient_y_2D(a,i-1,j))/dy
|
||||
return output
|
||||
if output.ndim == 3:
|
||||
for i in np.arange(1,nlat-1):
|
||||
for j in range(nlon):
|
||||
for k in range(nlevels-1):
|
||||
output[i,j,k] = (scalar_gradient_x(a,i,(j+1)%nlon,k) - scalar_gradient_x(a,i,(j-1)%nlon,k))/dx[i] + (scalar_gradient_y(a,i+1,j,k) - scalar_gradient_y(a,i-1,j,k))/dy + (scalar_gradient_z(a,i,j,k+1)-scalar_gradient_z(a,i,j,k-1))/(2*dz[k])
|
||||
return output
|
||||
|
||||
# divergence of (a*u) where a is a scalar field and u is the atmospheric velocity field
|
||||
def divergence_with_scalar(a,u,v,dx,dy):
|
||||
output = np.zeros_like(a)
|
||||
nlat, nlon, nlevels = output.shape[:]
|
||||
|
||||
au = a*u
|
||||
av = a*v
|
||||
|
||||
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.1*scalar_gradient_z(a*w,i,j,k)
|
||||
return output
|
||||
|
||||
def 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):
|
||||
# calculate change in temperature of ground and atmosphere due to radiative imbalance
|
||||
nlat, nlon, nlevels = temperature_atmos.shape[:]
|
||||
|
||||
upward_radiation = np.zeros(nlevels)
|
||||
downward_radiation = np.zeros(nlevels)
|
||||
optical_depth = np.zeros(nlevels)
|
||||
Q = np.zeros(nlevels)
|
||||
|
||||
for i in range(nlat):
|
||||
for j in range(nlon):
|
||||
# calculate optical depth
|
||||
pressure_profile = air_pressure[i,j,:]
|
||||
density_profile = air_density[i,j,:]
|
||||
fl = 0.1
|
||||
optical_depth = low_level.surface_optical_depth(lat[i])*(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])
|
||||
|
||||
# calculate downward longwave flux, bc is zero at TOA (in model)
|
||||
downward_radiation[-1] = 0
|
||||
for k in np.arange(0,nlevels-1)[::-1]:
|
||||
downward_radiation[k] = (downward_radiation[k+1] - low_level.thermal_radiation(temperature_atmos[i,j,k])*(optical_depth[k+1]-optical_depth[k]))/(1 + optical_depth[k] - optical_depth[k+1])
|
||||
|
||||
# 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,0,0,k)/(1E3*density_profile[k])
|
||||
# make sure model does not have a higher top than 50km!!
|
||||
# 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)
|
||||
|
||||
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
|
||||
|
||||
def velocity_calculation(u,v,air_pressure,old_pressure,air_density,coriolis,gravity,dx,dy,dt):
|
||||
# introduce temporary arrays to update velocity in the atmosphere
|
||||
u_temp = np.zeros_like(u)
|
||||
v_temp = np.zeros_like(v)
|
||||
w_temp = np.zeros_like(u)
|
||||
|
||||
nlat,nlon,nlevels = air_pressure.shape[:]
|
||||
|
||||
# calculate acceleration of atmosphere using primitive equations on beta-plane
|
||||
for i in np.arange(1,nlat-1):
|
||||
for j in range(nlon):
|
||||
for k in range(nlevels):
|
||||
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) + coriolis[i]*v[i,j,k] - low_level.scalar_gradient_x(air_pressure,dx,nlon,i,j,k)/air_density[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) - coriolis[i]*u[i,j,k] - low_level.scalar_gradient_y(air_pressure,dy,nlat,i,j,k)/air_density[i,j,k] )
|
||||
w_temp[i,j,k] += -(air_pressure[i,j,k]-old_pressure[i,j,k])/(dt*air_density[i,j,k]*gravity)
|
||||
|
||||
u += u_temp
|
||||
v += v_temp
|
||||
w = w_temp
|
||||
|
||||
# approximate friction
|
||||
u *= 0.95
|
||||
v *= 0.95
|
||||
|
||||
return u,v,w
|
||||
Loading…
Reference in New Issue