Changeset 727 for configs/testing/python

Timestamp:

08/22/18 17:30:37 (6 years ago)

Author:

dubos

Message:

testing : compare output to reference outputs

File:

• configs/testing/python/common.py (modified) (4 diffs)

configs/testing/python/common.py

@@ -44 +44 @@
         #        plt.yticks(np.arange(0, 10001, 1000))
         plt.savefig('%s.png'%varname)
+        plt.close()
     lon, lat, Omega,T,u,Phi = getvars(nc, 'lon','lat','OMEGA', 'T', 'U', 'PHI')
     plot_var(nlon,nlat,llm, lon,Phi,Omega,'Omega')
@@ -64 +65 @@
         axis_latitude()
         plt.savefig('T850.png')
+        plt.close()
         
-    def plot_dT(nlon,nlat,llm, lon,T,p,Phi): # perturbation temp, final time
+    def plot_dT(nlon,nlat,llm, lon,T_ref,T,p,Phi): # perturbation temp, final time
         # vertical slice at final time
         print 'Reading data ...'
-        T,p, Phi = T[-1,:,nlat/2,:], p[-1,:,nlat/2,:], Phi[-1,:,nlat/2,:]
+        T_ref,T,p,Phi = [ x[-1,:,nlat/2,:] for x in T_ref,T,p,Phi]
         print '... done.'
         N, Teq, peq = 0.01, 300., 1e5
@@ -85 +87 @@
         plt.yticks(np.arange(0, 10001, 1000))
         plt.savefig('dT.png')
+        plt.close()
+
+        plt.figure(figsize=(12,6))
+        plt.contourf(lon2,z,T-T_ref)
+        plt.colorbar() 
+        plt.title('$T-T_{ref}$')
+        axis_longitude()
+        plt.ylabel('z (m)')
+        plt.yticks(np.arange(0, 10001, 1000))
+        plt.savefig('dT_ref.png')
+        plt.close()
 
     lon, lat, T850, T, Phi, p = getvars(nc, 'lon','lat','T850', 'T', 'PHI','P')
-    plot_dT(nlon,nlat,llm, lon,T,p,Phi)
+    T_ref, = getvars(nc_ref, 'T')
+    plot_dT(nlon,nlat,llm, lon,T_ref,T,p,Phi)
+#    plot_dT_ref(nlon,nlat,llm, lon,T_ref,T,p,Phi)
     plot_T850(lon,lat,T850)
 
@@ -116 +131 @@
 #--------------------------- MAIN ----------------------------
 
-gridfile = 'netcdf/output_dcmip2016_regular.nc'
-nc = cdf.Dataset(gridfile, "r")
+filename = 'output_dcmip2016_regular.nc'
+nc = cdf.Dataset('netcdf/%s'%filename, "r")
+nc_ref = cdf.Dataset('netcdf_ref/%s'%filename, "r")
+
 llm, nlon, nlat, ntime = getdims(nc, 'lev','lon','lat','time_counter')
 Cpd, kappa, g = 1004.5, 0.2857143, 9.80616