Added EKE calculation and saving

Author: evanmason

haversine_distmat.f90

@@ -27,7 +27,7 @@
 !
 !    If you have ifort on your system, change 'gfortran' to 'intelem'
 !
-!   Version 1.2.1
+!   Version 1.3.0
 !
 !===========================================================================
 

make_eddy_track_AVISO.py

@@ -24,7 +24,7 @@
 
 make_eddy_track_AVISO.py
 
-Version 1.2.1
+Version 1.3.0
 
 
 Scroll down to line ~640 to get started
@@ -401,9 +401,30 @@ def __init__(self, AVISO_file, lonmin, lonmax, latmin, latmax, with_pad=True, us
         self.set_index_padding()
         self.make_gridmask(with_pad, use_maskoceans).uvmask()
         self.get_AVISO_f_pm_pn()
+        self.set_u_v_eke()
 
 
 
+    def set_u_v_eke(self, pad=2):
+        '''
+        '''
+        double_pad = pad * 2
+        if self.zero_crossing:
+            u1 = np.empty((self.ip0, self.jp1 - self.jp0))
+            u0 = np.empty((self._lon.shape[0] - self.ip1, self.jp1 - self.jp0))
+            self.u = np.ma.concatenate((u0, u1), axis=1)
+        else:
+            self.u = np.empty((self.jp1 - self.jp0, self.ip1 - self.ip0))
+        self.v = np.empty_like(self.u)
+        self.eke = np.empty((self.u.shape[0] - double_pad,
+                             self.u.shape[1] - double_pad))
+        return self
+
+
+
+
+
+
     def get_AVISO_data(self, AVISO_file):
         '''
         Read nc data from AVISO file
@@ -492,13 +513,21 @@ def getSurfGeostrVel(self, zeta):
         by Evan Mason
         Changed to gv2, 14 May 07...
         '''
-        def gv2(zeta, gof, pm, pn, umask, vmask): # Pierrick's version
-            ugv = -gof * self.v2rho_2d(vmask * (zeta[1:] - zeta[:-1]) \
-                                        * 0.5 * (pn[1:] + pn[:-1]))
-            vgv =  gof * self.u2rho_2d(umask * (zeta[:, 1:] - zeta[:, :-1]) \
-                                        * 0.5 * (pm[:, 1:] + pm[:, :-1]))
-            return ugv, vgv
-        return gv2(zeta, self.gof(), self.pm(), self.pn(), self.umask(), self.vmask())
+        #def gv2(zeta, gof, pm, pn, umask, vmask): # Pierrick's version
+        self.u[:] = -self.gof() * self.v2rho_2d(self.vmask() * (zeta.data[1:] - zeta.data[:-1]) \
+                                        * 0.5 * (self.pn()[1:] + self.pn()[:-1]))
+        self.v[:] =  self.gof() * self.u2rho_2d(self.umask() * (zeta.data[:, 1:] - zeta.data[:, :-1]) \
+                                        * 0.5 * (self.pm()[:, 1:] + self.pm()[:, :-1]))
+            #return ugv, vgv
+        #return gv2(zeta, self.gof(), self.pm(), self.pn(), self.umask(), self.vmask())
+        return self
+
+    def getEKE(self, u, v):
+        '''
+        '''
+        np.add(u**2, v**2, out=self.eke)
+        self.eke *= 0.5
+        return self
 
 
     def lon(self):
@@ -782,15 +811,15 @@ def pcol_2dxy(self, x, y):
     subdomain = True
     if the_domain in 'Global':
 
-        #lonmin = -36.     # Canary
-        #lonmax = -7.5
-        #latmin = 18.
-        #latmax = 33.2
+        lonmin = -36.     # Canary
+        lonmax = -7.5
+        latmin = 18.
+        latmax = 33.2
 
-        lonmin = -65.     # Corridor
-        lonmax = -5.5
-        latmin = 11.5
-        latmax = 38.5
+        #lonmin = -65.     # Corridor
+        #lonmax = -5.5
+        #latmin = 11.5
+        #latmax = 38.5
 
         #lonmin = -179.     # SEP
         #lonmax = -65
@@ -1174,14 +1203,16 @@ def pcol_2dxy(self, x, y):
                 #else:
                     #umask, vmask, junk = sla_grd.uvpmask(np.zeros_like(sla))
 
-                # Multiply by 0.01 for m/s
-                u, v = sla_grd.getSurfGeostrVel(sla * 0.01)
+                # Multiply by 0.01 for m
+                sla_grd.getSurfGeostrVel(sla * 0.01)
 
                 # Remove padded boundary
                 sla = sla[sla_grd.jup0:sla_grd.jup1, sla_grd.iup0:sla_grd.iup1]
-                u = u[sla_grd.jup0:sla_grd.jup1, sla_grd.iup0:sla_grd.iup1]
-                v = v[sla_grd.jup0:sla_grd.jup1, sla_grd.iup0:sla_grd.iup1]
-
+                u = sla_grd.u[sla_grd.jup0:sla_grd.jup1, sla_grd.iup0:sla_grd.iup1]
+                v = sla_grd.v[sla_grd.jup0:sla_grd.jup1, sla_grd.iup0:sla_grd.iup1]
+                
+                sla_grd.getEKE(u, v)
+                
                 A_eddy.sla = np.ma.copy(sla)
                 C_eddy.sla = np.ma.copy(sla)
                 A_eddy.slacopy = np.ma.copy(sla)
@@ -1207,7 +1238,9 @@ def pcol_2dxy(self, x, y):
                 if anim_figs:
                    animfig = plt.figure(999)
                    animax = animfig.add_subplot(111)
-
+                   # Colorbar axis
+                   animax_cbar = get_cax(animax, dx=0.03, width=.05, position='b')
+                   
             if 'Q' in diag_type:
                 CS = ax.contour(sla_grd.lon(),
                                  sla_grd.lat(), qparam, qparameter)
@@ -1325,7 +1358,7 @@ def pcol_2dxy(self, x, y):
                     #print 'figure saving'
                     #tt = time.time()
                     anim_figure(A_eddy, C_eddy, Mx, My, pMx, pMy, plt.cm.RdBu_r, rtime, diag_type, 
-                                savedir, 'SLA ' + tit, animax)
+                                savedir, 'SLA ' + tit, animax, animax_cbar)
                     #print 'figure saving done in %s seconds\n' %(time.time() - tt)
                 #anim_fig.start()
 

make_eddy_track_ROMS.py

@@ -24,7 +24,7 @@
 
 make_eddy_track_ROMS.py
 
-Version 1.1.0
+Version 1.3.0
 
 
 

make_eddy_tracker_list_obj.py

@@ -25,7 +25,7 @@
 
 make_eddy_tracker_list_obj.py
 
-Version 1.2.1
+Version 1.3.0
 
 
 ===========================================================================
@@ -123,7 +123,7 @@ class track (object):
         saved2nc - Becomes True once saved to netcdf file
         dayzero - True at first appearance of eddy
     '''
-    def __init__(self, eddy_index, datatype, lon, lat, time, Uavg,
+    def __init__(self, eddy_index, datatype, lon, lat, time, Uavg, teke,
                        radius_s, radius_e, amplitude, bounds,
                        temp=None,
                        salt=None):
@@ -134,6 +134,7 @@ def __init__(self, eddy_index, datatype, lon, lat, time, Uavg,
         self.lat        = np.atleast_1d(lat)
         self.ocean_time = np.atleast_1d(time)
         self.Uavg       = np.atleast_1d(Uavg)
+        self.teke       = np.atleast_1d(teke)
         self.radius_s   = np.atleast_1d(radius_s) # speed-based eddy radius
         self.radius_e   = np.atleast_1d(radius_e) # effective eddy radius
         self.amplitude  = np.atleast_1d(amplitude)
@@ -146,7 +147,7 @@ def __init__(self, eddy_index, datatype, lon, lat, time, Uavg,
         self.saved2nc   = False
 
 
-    def append_pos(self, lon, lat, time, Uavg,
+    def append_pos(self, lon, lat, time, Uavg, teke,
                    radius_s, radius_e, amplitude, bounds, temp, salt):
         '''
         Append track updates
@@ -155,6 +156,7 @@ def append_pos(self, lon, lat, time, Uavg,
         self.lat = np.r_[self.lat, lat]
         self.ocean_time = np.r_[self.ocean_time, time]
         self.Uavg = np.r_[self.Uavg, Uavg]
+        self.teke = np.r_[self.teke, teke]
         self.radius_s = np.r_[self.radius_s, radius_s]
         self.radius_e = np.r_[self.radius_e, radius_e]
         self.amplitude = np.r_[self.amplitude, amplitude]
@@ -210,6 +212,7 @@ def __init__(self, datatype, track_duration_min):
         self.new_radii_e = np.array([])
         self.new_amp    = np.array([])
         self.new_Uavg   = np.array([])
+        self.new_teke   = np.array([])
         if 'ROMS' in self.datatype:
             self.new_temp = np.array([])
             self.new_salt = np.array([])
@@ -221,6 +224,7 @@ def __init__(self, datatype, track_duration_min):
         self.old_radii_e  = np.array([])
         self.old_amp    = np.array([])
         self.old_Uavg   = np.array([])
+        self.old_teke   = np.array([])
         if 'ROMS' in self.datatype:
             self.old_temp = np.array([])
             self.old_salt = np.array([])
@@ -235,24 +239,24 @@ def __init__(self, datatype, track_duration_min):
         assert datatype in ('ROMS', 'AVISO'), "Unknown string in 'datatype' parameter"
 
 
-    def append_list(self, lon, lat, time, Uavg,
+    def append_list(self, lon, lat, time, Uavg, teke,
             radius_s, radius_e, amplitude, bounds, temp=None, salt=None):
         '''
         Append a new 'track' object to the list
         '''
         self.tracklist.append(track(self.index, self.datatype,
-                                    lon, lat, time, Uavg,
+                                    lon, lat, time, Uavg, teke,
                                     radius_s, radius_e,
                                     amplitude, bounds,
                                     temp, salt))
 
 
-    def update_track(self, index, lon, lat, time, Uavg,
+    def update_track(self, index, lon, lat, time, Uavg, teke,
              radius_s, radius_e, amplitude, bounds, temp=None, salt=None):
         '''
         Update a track at index
         '''
-        self.tracklist[index].append_pos(lon, lat, time, Uavg,
+        self.tracklist[index].append_pos(lon, lat, time, Uavg, teke,
                                          radius_s, radius_e,
                                          amplitude, bounds,
                                          temp, salt)
@@ -268,6 +272,7 @@ def reset_holding_variables(self):
         self.new_radii_tmp_e = np.array([])
         self.new_amp_tmp = np.array([])
         self.new_Uavg_tmp = np.array([])
+        self.new_teke_tmp = np.array([])
         self.new_time_tmp = np.array([])
         self.new_temp_tmp = np.array([])
         self.new_salt_tmp = np.array([])
@@ -285,6 +290,7 @@ def set_old_variables(self):
         self.old_radii_e = np.copy(self.new_radii_tmp_e)
         self.old_amp = np.copy(self.new_amp_tmp)
         self.old_Uavg = np.copy(self.new_Uavg_tmp)
+        self.old_teke = np.copy(self.new_teke_tmp)
         self.old_temp = np.copy(self.new_temp_tmp)
         self.old_salt = np.copy(self.new_salt_tmp)
 
@@ -391,6 +397,7 @@ def create_netcdf(self, directory, savedir, title,
         nc.createVariable('A', 'f8', ('Nobs'), fill_value=self.fillval)
         nc.createVariable('L', 'f8', ('Nobs'), fill_value=self.fillval)
         nc.createVariable('U', 'f8', ('Nobs'), fill_value=self.fillval)
+        nc.createVariable('Teke', 'f8', ('Nobs'), fill_value=self.fillval)
         nc.createVariable('radius_e', 'f8', ('Nobs'), fill_value=self.fillval)
         if 'Q' in self.diag_type:
             nc.createVariable('qparameter', 'f8', ('Nobs'), fill_value=self.fillval)
@@ -446,10 +453,11 @@ def create_netcdf(self, directory, savedir, title,
                                         'between the extremum of SSH within ' + \
                                         'the eddy and the SSH around the contour ' + \
                                         'defining the eddy perimeter'
+        
         nc.variables['L'].units = 'km'
         nc.variables['L'].min_val = self.radmin / 1000.
         nc.variables['L'].max_val = self.radmax / 1000.
-        nc.variables['L'].long_name = 'radius scale'
+        nc.variables['L'].long_name = 'speed radius scale'
         nc.variables['L'].description = 'radius of a circle whose area is equal ' + \
                                         'to that enclosed by the contour of ' + \
                                         'maximum circum-average speed'
@@ -460,6 +468,13 @@ def create_netcdf(self, directory, savedir, title,
         nc.variables['U'].long_name = 'maximum circum-averaged speed'
         nc.variables['U'].description = 'average speed of the contour defining ' + \
                                         'the radius scale L'
+        nc.variables['Teke'].units = 'm^2/sec^2'
+        #nc.variables['Teke'].min = 0.
+        #nc.variables['Teke'].max = 376.6
+        nc.variables['Teke'].long_name = 'sum EKE within contour Ceff'
+        nc.variables['Teke'].description = 'sum of eddy kinetic energy within contour defining ' + \
+                                        'the effective radius'
+        
         nc.variables['radius_e'].units = 'km'
         nc.variables['radius_e'].min_val = self.radmin / 1000.
         nc.variables['radius_e'].max_val = self.radmax / 1000.
@@ -495,6 +510,7 @@ def _reduce_inactive_tracks(self):
                 track.qparameter = False
                 track.amplitude = False
                 track.Uavg = False
+                track.teke = False
                 track.radius_s = False
                 track.radius_e = False
                 track.bounds = False
@@ -549,6 +565,7 @@ def write2chelton_nc(self, savedir, rtime):
                         nc.variables['lat'][self.ncind:tend] = np.array([self.tracklist[i].lat])
                         nc.variables['A'][self.ncind:tend] = np.array([self.tracklist[i].amplitude])
                         nc.variables['U'][self.ncind:tend] = np.array([self.tracklist[i].Uavg]) * 100. # to cm/s
+                        nc.variables['Teke'][self.ncind:tend] = np.array([self.tracklist[i].teke])
                         nc.variables['L'][self.ncind:tend] = np.array([self.tracklist[i].radius_s]) * np.array(1e-3)
                         nc.variables['radius_e'][self.ncind:tend] = np.array([self.tracklist[i].radius_e]) * np.array(1e-3)
                         if 'ROMS' in self.datatype:
@@ -594,6 +611,7 @@ def write2chelton_nc(self, savedir, rtime):
         self.old_radii_e = self.new_radii_e[lasti:]
         self.old_amp     = self.new_amp[lasti:]
         self.old_Uavg    = self.new_Uavg[lasti:]
+        self.old_teke    = self.new_teke[lasti:]
         if 'ROMS' in self.datatype:
             self.old_temp = self.new_temp[lasti:]
             self.old_salt = self.new_salt[lasti:]
@@ -604,6 +622,7 @@ def write2chelton_nc(self, savedir, rtime):
         self.new_radii_e = np.array([])
         self.new_amp     = np.array([])
         self.new_Uavg    = np.array([])
+        self.new_teke    = np.array([])
         self.new_time    = np.array([])
         if 'ROMS' in self.datatype:
             self.new_temp = np.array([])
@@ -669,6 +688,14 @@ def insert_at_index(self, xarr, ind, x):
             self.new_amp = np.zeros((newsize))
             self.new_amp[:tmp.size] = tmp
             self.new_amp[ind] = x
+        elif 'new_Uavg' in xarr:
+            self.new_Uavg = np.zeros((newsize))
+            self.new_Uavg[:tmp.size] = tmp
+            self.new_Uavg[ind] = x
+        elif 'new_teke' in xarr:
+            self.new_teke = np.zeros((newsize))
+            self.new_teke[:tmp.size] = tmp
+            self.new_teke[ind] = x
         elif 'new_temp' in xarr:
             self.new_temp = np.zeros((newsize))
             self.new_temp[:tmp.size] = tmp
@@ -683,7 +710,7 @@ def insert_at_index(self, xarr, ind, x):
 
 
     def update_eddy_properties(self, centlon, centlat, eddy_radius_s, eddy_radius_e,
-                               amplitude, Uavg, rtime, bounds, cent_temp=None, cent_salt=None):
+                               amplitude, Uavg, teke, rtime, bounds, cent_temp=None, cent_salt=None):
         '''
         Append new variable values to track arrays
         '''
@@ -693,6 +720,7 @@ def update_eddy_properties(self, centlon, centlat, eddy_radius_s, eddy_radius_e,
         self.new_radii_tmp_e = np.r_[self.new_radii_tmp_e, eddy_radius_e]
         self.new_amp_tmp = np.r_[self.new_amp_tmp, amplitude]
         self.new_Uavg_tmp = np.r_[self.new_Uavg_tmp, Uavg]
+        self.new_teke_tmp = np.r_[self.new_teke_tmp, teke]
         self.new_time_tmp = np.r_[self.new_time_tmp, rtime]
         if 'ROMS' in self.datatype:
             self.new_temp_tmp = np.r_[self.new_temp_tmp, cent_temp]