diff --git a/claude_low_level_library.py b/claude_low_level_library.py deleted file mode 100644 index e81d5f8..0000000 --- a/claude_low_level_library.py +++ /dev/null @@ -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) \ No newline at end of file diff --git a/claude_top_level_library.py b/claude_top_level_library.py deleted file mode 100644 index 79f6745..0000000 --- a/claude_top_level_library.py +++ /dev/null @@ -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 \ No newline at end of file