diff --git a/claude_top_level_library.pyx b/claude_top_level_library.pyx index e6f70f0..77151e7 100644 --- a/claude_top_level_library.pyx +++ b/claude_top_level_library.pyx @@ -53,7 +53,7 @@ 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) + 1*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(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 @@ -74,7 +74,7 @@ def thermal_advection(np.ndarray temperature_atmos,np.ndarray air_pressure,np.nd 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[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) @@ -129,7 +129,7 @@ def radiation_calculation(np.ndarray temperature_world, np.ndarray temperature_a return temperature_world, temperature_atmos -def velocity_calculation(np.ndarray u,np.ndarray v,np.ndarray air_pressure,np.ndarray old_pressure,np.ndarray air_density,np.ndarray coriolis,DTYPE_f gravity,np.ndarray dx,DTYPE_f dy,DTYPE_f dt): +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): # 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) @@ -144,22 +144,20 @@ def velocity_calculation(np.ndarray u,np.ndarray v,np.ndarray air_pressure,np.nd inv_dt_grav = 1/(dt*gravity) # calculate acceleration of atmosphere using primitive equations on beta-plane - for i in np.arange(3,nlat-3).tolist(): + 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) + coriolis[i]*v[i,j,k] - low_level.scalar_gradient_x(air_pressure,dx,nlon,i,j,k)*inv_density ) - 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)*inv_density ) - w_temp[i,j,k] += -(air_pressure[i,j,k]-old_pressure[i,j,k])*inv_density*inv_dt_grav + 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 += u_temp v += v_temp - w = w_temp - - # approximate friction - u *= 0.95 - v *= 0.95 + w += w_temp + # approximate surface friction u[:,:,0] *= 0.8 v[:,:,0] *= 0.8 diff --git a/toy_model.py b/toy_model.py index c9715a1..0a609d4 100644 --- a/toy_model.py +++ b/toy_model.py @@ -22,8 +22,8 @@ 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 = 7*day +dt_main = 60*0.1 +spinup_length = 10*day ### @@ -32,7 +32,7 @@ advection_boundary = 3 # how many gridpoints away from poles to apply advectio smoothing_parameter_t = 0.6 smoothing_parameter_u = 0.6 smoothing_parameter_v = 0.6 -smoothing_parameter_w = 0.3 +smoothing_parameter_w = 0.2 save = False # save current state to file? load = True # load initial state from file? @@ -103,6 +103,7 @@ 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 @@ -201,10 +202,10 @@ while True: if velocity: 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,v,w = top_level.velocity_calculation(u,v,w,air_pressure,ref_pressure_profile,air_density,coriolis,gravity,dx,dy,dz,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.3) + w = top_level.smoothing_3D(w,smoothing_parameter_w,0.1) u[(advection_boundary,-advection_boundary-1),:,:] *= 0.5 v[(advection_boundary,-advection_boundary-1),:,:] *= 0.5 @@ -241,7 +242,6 @@ while True: # time_taken = float(round(time.time() - before_advection,3)) # print('Advection: ',str(time_taken),'s') - # before_plot = time.time() if plot: @@ -268,8 +268,8 @@ while True: 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].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)