Update to establish running version with creating meshgrid

calc_ref_state.py

@@ -2,7 +2,24 @@
 # Program name  : calc_ref_state.py
 # written by    : Gabriel Wolf, g.a.wolf@reading.ac.uk
 #                 adapted from compute_reference_state_woa2009_for_Juan.m of Remi Tailleux
-# last modified : 13.09.2018
+# last modified : 19.09.2018
+# ############################################################################################
+# More information in the README.txt file
+# Necessary user input:
+#     manual user input under Input 0
+#     reading of data under step a
+# -> everything else should work by using this Code
+# Table of content of this programme:
+# - Input 0 :
+# - Input 1 :
+# - step a  :
+# - step b  :
+# - step c  :
+# - step d  :
+# - step e  :
+# - step f  :
+# - step g  :
+# - step h  :
 # ############################################################################################
 
 # Load necessary modules #####################################################################
@@ -33,7 +50,7 @@ saveplot   = input('Do you want to save the plots? Enter 1 for yes or 0 for no')
 dir_plot   = '/glusterfs/inspect/users/xg911182/Code/Python/plots_test/'
 lonlat_reg = [25.0, 25.0+360.0, -90.0, 90.0] # plotting range for 2dmap
 plot_2d    = 1 # plot 2d fields (lat-lon)
-plot_cs    = 1 # plot crosssections
+plot_cs    = 0 # plot crosssections
 cs_plot    = {'lon':[-47.5,172.5],'lat':[0]} # crosssections
 # Input 1 ####################################################################################
 # Physical constants
@@ -43,9 +60,8 @@ cs_plot    = {'lon':[-47.5,172.5],'lat':[0]} # crosssections
 #     rho0      : reference density     [kg m**-3]
 #     grav      : gravity               [m s**-2]
 #     gbuo      : 
-# PHYSICAL CONSTATS MISSING
 # ############################################################################################
-# constants stored in myconstants.py
+from myconstants import rho0, grav
 
 # step a #####################################################################################
 # Read necessary data to calculate reference state
@@ -71,16 +87,30 @@ print '----------------------------------------------'
 # step b #####################################################################################
 # Define grid for data from step a
 # input         : lon, lat, zlev (all from step a) 
-# output        : plev
+# output        : lat_xy, lat_yz, lat_xyz, lon_xy, lon_xyz, p_xyz, plev, z_yz
-#     plev      : MISSING DEFINITION
+#     lat_*     : meshgrid for latitude in x/y/z-dir  [degrees]
+#     lon_*     : meshgrid for longitude in x/y/z-dir [degrees]
+#     p_*       : meshgrid for pressure in x/y/z-dir  [dbar]
+#     plev      : gauge pressure (p_abs-10.1325 dbar) [dbar]
+#     z_*       : meshgrid for depth in x/y/z-dir     [m]
 # ############################################################################################
 beg_time = time.time()
 print ''
 print '----------------------------------------------'
-print 'Insert code to define grid of data'
+print 'Create 2d and 3d meshgrids'
+try:
+    plev  = grd.p
+except:
+    plev   = (rho0*grav/1e4)*grd.z
+    grd.p  = plev
+    grd.Up = 'dbar'
+    grd.Np = len(plev)
+lat_xy, lon_xy          = np.meshgrid(grd.lat,grd.lon)
+p_yz, lat_yz            = np.meshgrid(grd.z, grd.lat)
+lat_xyz, lon_xyz, p_xyz = np.meshgrid(grd.lat, grd.lon, grd.p)
 print '----------------------------------------------'
 end_time = time.time()
-print 'Elapsed time to define grid:       '+ '{:.2f}'.format(end_time-beg_time)+'s'
+print 'Elapsed time to create meshgrids:  '+ '{:.2f}'.format(end_time-beg_time)+'s'
 print '----------------------------------------------'
 
 # step c #####################################################################################
@@ -155,9 +185,11 @@ print '----------------------------------------------'
 
 # step g #####################################################################################
 # plot of full 2d fields
-# input         : date_str, dir_plot, i_dep, i_t, lonlat_reg
+# input         : data, date_str, dir_plot, grd, i_dep, i_t, lonlat_reg
+#     data      : output of step a
 #     date_str  : MISSING - DONE LOCALLY HERE
 #     dir_plot  : output of Input 0
+#     grd       : output of step a
 #     i_dep     : MISSING - DONE LOCALLY HERE
 #     i_t       : MISSING - DONE LOCALLY HERE
 #     lonlat_reg: output of Input 0
@@ -169,8 +201,8 @@ print ''
 print '----------------------------------------------'
 print 'Plot 2d fields (lat-lon)'
 list_allkeys = data.keys() # all variables
-i_dep = 0
+i_dep    = 0
-i_t   = 0
+i_t      = 0
 date_str = '2009'
 for i_key in range(0,len(list_allkeys)):
     # create plotname
@@ -191,20 +223,25 @@ print '----------------------------------------------'
 
 # step h #####################################################################################
 # plot of crossection fields
-# input         : date_str, dir_plot, i_t, plot_cs, saveplot
+# input         : cs_plot, data, date_str, dir_plot, grd, i_t, plot_cs, saveplot
+#     cs_plot   : output of Input 0
+#     data      : output of step a
 #     date_str  : MISSING - DONE LOCALLY HERE
 #     dir_plot  : output of Input 0
+#     grd       : output of step a 
 #     i_t       : MISSING - DONE LOCALLY HERE
 #     plot_cs   : output of Input 0
 #     saveplot  : output of Input 0
 # output        : None
 # saved/created : plots saved in dir_plot
 # ############################################################################################
+date_str = '2009'
+i_t      = 0
 if plot_cs:
     beg_time = time.time()
     print ''
     print '----------------------------------------------'
-    print 'Insert code to plot crossections'
+    print 'Plot lon-p and lat-p crossections'
     list_allkeys = data.keys() # all variables
     for i_key in range(0,len(list_allkeys)):
         for i_cs in range(0,len(cs_plot['lat'])):