Post stream fix

Minor plotting addition and improved efficiency of vertical derivative

claude_low_level_library.pyx

@@ -34,7 +34,6 @@ def scalar_gradient_y_2d(a,dy,nlat,i,j):
 		return (a[i+1,j]-a[i-1,j])/dy
 
 def scalar_gradient_z(a,dz,i,j,k):
-	cdef np.ndarray output = np.zeros_like(a)
 	cdef np.int_t nlevels = len(dz)
 	if k == 0:
 		return (a[i,j,k+1]-a[i,j,k])/dz[k]
@@ -43,8 +42,7 @@ def scalar_gradient_z(a,dz,i,j,k):
 	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 i,np.int_t j,np.int_t k):
+def scalar_gradient_z_1D(np.ndarray a,np.ndarray dz,np.int_t k):
-	cdef np.ndarray output = np.zeros_like(a)
 	cdef np.int_t nlevels = len(dz)
 	if k == 0:
 		return (a[k+1]-a[k])/dz[k]

claude_top_level_library.c

@@ -3638,7 +3638,7 @@ static PyObject *__pyx_pf_24claude_top_level_library_2radiation_calculation(CYTH
  * 
  * 			# 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])
+ * 				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!!
  */
       __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 72, __pyx_L1_error)
@@ -3710,7 +3710,7 @@ static PyObject *__pyx_pf_24claude_top_level_library_2radiation_calculation(CYTH
         /* "claude_top_level_library.pyx":73
  * 			# 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])             # <<<<<<<<<<<<<<
+ * 				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!!
  * 				# approximate SW heating of ozone
  */
@@ -3735,8 +3735,8 @@ static PyObject *__pyx_pf_24claude_top_level_library_2radiation_calculation(CYTH
         }
         #if CYTHON_FAST_PYCALL
         if (PyFunction_Check(__pyx_t_17)) {
-          PyObject *__pyx_temp[6] = {__pyx_t_15, __pyx_t_3, ((PyObject *)__pyx_v_dz), __pyx_int_0, __pyx_int_0, __pyx_v_k};
+          PyObject *__pyx_temp[4] = {__pyx_t_15, __pyx_t_3, ((PyObject *)__pyx_v_dz), __pyx_v_k};
-          __pyx_t_1 = __Pyx_PyFunction_FastCall(__pyx_t_17, __pyx_temp+1-__pyx_t_12, 5+__pyx_t_12); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 73, __pyx_L1_error)
+          __pyx_t_1 = __Pyx_PyFunction_FastCall(__pyx_t_17, __pyx_temp+1-__pyx_t_12, 3+__pyx_t_12); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 73, __pyx_L1_error)
           __Pyx_XDECREF(__pyx_t_15); __pyx_t_15 = 0;
           __Pyx_GOTREF(__pyx_t_1);
           __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
@@ -3744,15 +3744,15 @@ static PyObject *__pyx_pf_24claude_top_level_library_2radiation_calculation(CYTH
         #endif
         #if CYTHON_FAST_PYCCALL
         if (__Pyx_PyFastCFunction_Check(__pyx_t_17)) {
-          PyObject *__pyx_temp[6] = {__pyx_t_15, __pyx_t_3, ((PyObject *)__pyx_v_dz), __pyx_int_0, __pyx_int_0, __pyx_v_k};
+          PyObject *__pyx_temp[4] = {__pyx_t_15, __pyx_t_3, ((PyObject *)__pyx_v_dz), __pyx_v_k};
-          __pyx_t_1 = __Pyx_PyCFunction_FastCall(__pyx_t_17, __pyx_temp+1-__pyx_t_12, 5+__pyx_t_12); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 73, __pyx_L1_error)
+          __pyx_t_1 = __Pyx_PyCFunction_FastCall(__pyx_t_17, __pyx_temp+1-__pyx_t_12, 3+__pyx_t_12); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 73, __pyx_L1_error)
           __Pyx_XDECREF(__pyx_t_15); __pyx_t_15 = 0;
           __Pyx_GOTREF(__pyx_t_1);
           __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
         } else
         #endif
         {
-          __pyx_t_16 = PyTuple_New(5+__pyx_t_12); if (unlikely(!__pyx_t_16)) __PYX_ERR(0, 73, __pyx_L1_error)
+          __pyx_t_16 = PyTuple_New(3+__pyx_t_12); if (unlikely(!__pyx_t_16)) __PYX_ERR(0, 73, __pyx_L1_error)
           __Pyx_GOTREF(__pyx_t_16);
           if (__pyx_t_15) {
             __Pyx_GIVEREF(__pyx_t_15); PyTuple_SET_ITEM(__pyx_t_16, 0, __pyx_t_15); __pyx_t_15 = NULL;
@@ -3762,15 +3762,9 @@ static PyObject *__pyx_pf_24claude_top_level_library_2radiation_calculation(CYTH
           __Pyx_INCREF(((PyObject *)__pyx_v_dz));
           __Pyx_GIVEREF(((PyObject *)__pyx_v_dz));
           PyTuple_SET_ITEM(__pyx_t_16, 1+__pyx_t_12, ((PyObject *)__pyx_v_dz));
-          __Pyx_INCREF(__pyx_int_0);
-          __Pyx_GIVEREF(__pyx_int_0);
-          PyTuple_SET_ITEM(__pyx_t_16, 2+__pyx_t_12, __pyx_int_0);
-          __Pyx_INCREF(__pyx_int_0);
-          __Pyx_GIVEREF(__pyx_int_0);
-          PyTuple_SET_ITEM(__pyx_t_16, 3+__pyx_t_12, __pyx_int_0);
           __Pyx_INCREF(__pyx_v_k);
           __Pyx_GIVEREF(__pyx_v_k);
-          PyTuple_SET_ITEM(__pyx_t_16, 4+__pyx_t_12, __pyx_v_k);
+          PyTuple_SET_ITEM(__pyx_t_16, 2+__pyx_t_12, __pyx_v_k);
           __pyx_t_3 = 0;
           __pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_17, __pyx_t_16, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 73, __pyx_L1_error)
           __Pyx_GOTREF(__pyx_t_1);
@@ -3949,7 +3943,7 @@ static PyObject *__pyx_pf_24claude_top_level_library_2radiation_calculation(CYTH
  * 
  * 			# 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])
+ * 				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!!
  */
       }

claude_top_level_library.pyx

@@ -70,7 +70,7 @@ 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,0,0,k)/(1E3*density_profile[k])
+				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!!
 				# approximate SW heating of ozone
 				if heights[k] > 20E3:

toy_model.py

@@ -23,17 +23,17 @@ year = 365*day					# length of year (s)
 
 dt_spinup = 60*137
 dt_main = 60*9
-spinup_length = 0*day
+spinup_length = 5*day
 
 ###
 
 advection = True 				# if you want to include advection set this to be True
 advection_boundary = 8			# how many gridpoints away from poles to apply advection
 
-save = False 					# save current state to file?
+save = True 					# save current state to file?
-load = False  					# load initial state from file?
+load = True  					# load initial state from file?
 
-plot = True 					# display plots of output?
+plot = False 					# display plots of output?
 level_plots = True 				# display plots of output on vertical levels?
 nplots = 5						# how many levels you want to see plots of (evenly distributed through column)
 
@@ -215,6 +215,18 @@ while True:
 	
 	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
+
 		# update plot
 		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)
@@ -226,6 +238,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')
+		ax[1].plot(lat_plot,tropopause_height,color='black',linestyle='--',linewidth=3,alpha=0.5)
 		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}$')
@@ -270,4 +283,5 @@ while True:
 	if save:
 		pickle.dump((temperature_atmos,temperature_world,u,v,w,t,air_density,albedo), open("save_file.p","wb"))
 
-	# sys.exit()
+	if np.isnan(u.max()):
+		sys.exit()