added to output of TKE SGS scheme

pressel · Sep 3, 2015 · da66f8e · da66f8e
1 parent 46ef133
commit da66f8e
Show file tree

Hide file tree

Showing 6 changed files with 28 additions and 16 deletions.
diff --git a/Csrc/kinematics.h b/Csrc/kinematics.h
@@ -141,7 +141,8 @@ void compute_strain_rate_mag(const struct DimStruct *dims, double* restrict stra
     return;
 }
 
-void compute_wind_speed_angle(const struct DimStruct *dims, double* restrict u, double* restrict v, double* restrict wind_speed, double* restrict wind_angle){
+void compute_wind_speed_angle(const struct DimStruct *dims, double* restrict u, double* restrict v,
+ double* restrict wind_speed, double* restrict wind_angle, double u0, double v0){
     const size_t istride = dims->nlg[1] * dims->nlg[2];
     const size_t jstride = dims->nlg[2];
 
@@ -160,8 +161,8 @@ void compute_wind_speed_angle(const struct DimStruct *dims, double* restrict u,
             const size_t jshift = j*jstride;
             for(size_t k=kmin;k<kmax;k++){
                 const size_t ijk = ishift + jshift + k ;
-                const double u_interp = interp_2(u[ijk + istride], u[ijk]);
-                const double v_interp = interp_2(v[ijk + jstride], v[ijk]);
+                const double u_interp = interp_2(u[ijk + istride], u[ijk])+u0;
+                const double v_interp = interp_2(v[ijk + jstride], v[ijk])+v0;
                 wind_speed[ijk] = sqrt(u_interp * u_interp + v_interp * v_interp);
                 wind_angle[ijk] = atan2(v_interp,u_interp+1.0e-20);
             }

diff --git a/Csrc/sgs.h b/Csrc/sgs.h
@@ -7,10 +7,9 @@ double* restrict buoy_freq, double* restrict strain_rate_mag, double cs, double
     double delta = cbrt(dims->dx[0]*dims->dx[1]*dims->dx[2]);
 
     for (size_t i=0; i<dims->npg; i++){
-        double s2 = strain_rate_mag[i]*strain_rate_mag[i];
-        visc[i] = cs*cs*delta*delta*s2;
+        visc[i] = cs*cs*delta*delta*strain_rate_mag[i];
         if(buoy_freq[i] > 0.0){
-            double fb = sqrt(fmax(1.0 - buoy_freq[i]/(prt*s2),0.0));
+            double fb = sqrt(fmax(1.0 - buoy_freq[i]/(prt*strain_rate_mag[i]*strain_rate_mag[i]),0.0));
             visc[i] = visc[i] * fb;
         }
         diff[i] = visc[i]/prt;
@@ -21,7 +20,7 @@ double* restrict buoy_freq, double* restrict strain_rate_mag, double cs, double
 double tke_ell(double cn, double e, double buoy_freq, double delta){
     double ell;
     if(buoy_freq> 1.0e-10){
-        ell = fmax(fmin(cn*sqrt(fmax(e,0.0)/buoy_freq),delta),1.0e-10);
+        ell = fmax(fmin(cn*sqrt(fmax(e,0.0)/buoy_freq),delta),1e-10);
     }
     else{
         ell = delta;

diff --git a/Kinematics.pxd b/Kinematics.pxd
@@ -1,6 +1,7 @@
 cimport Grid
 cimport PrognosticVariables
 cimport ParallelMPI
+cimport ReferenceState
 from NetCDFIO cimport NetCDFIO_Stats
 
 cdef class Kinematics:
@@ -13,4 +14,4 @@ cdef class Kinematics:
         Py_ssize_t get_grad_shift(self, Grid.Grid Gr, Py_ssize_t vel_i, Py_ssize_t dx_j)
     cpdef initialize(self, Grid.Grid Gr, NetCDFIO_Stats NS, ParallelMPI.ParallelMPI Pa)
     cpdef update(self, Grid.Grid Gr, PrognosticVariables.PrognosticVariables PV)
-    cpdef stats_io(self, Grid.Grid Gr, PrognosticVariables.PrognosticVariables PV, NetCDFIO_Stats NS, ParallelMPI.ParallelMPI Pa)
+    cpdef stats_io(self, Grid.Grid Gr,ReferenceState.ReferenceState Ref, PrognosticVariables.PrognosticVariables PV, NetCDFIO_Stats NS, ParallelMPI.ParallelMPI Pa)
diff --git a/Kinematics.pyx b/Kinematics.pyx
@@ -1,6 +1,7 @@
 cimport Grid
 cimport PrognosticVariables
 cimport ParallelMPI
+cimport ReferenceState
 from NetCDFIO cimport NetCDFIO_Stats
 
 import numpy as np
@@ -10,7 +11,7 @@ cdef extern from "kinematics.h":
     void compute_velocity_gradient(Grid.DimStruct *dims, double *v, double *vgrad, long d)
     void compute_strain_rate(Grid.DimStruct *dims, double *vgrad, double *strain_rate)
     void compute_strain_rate_mag(Grid.DimStruct *dims, double *vgrad, double *strain_rate)
-    void compute_wind_speed_angle(Grid.DimStruct *dims, double *u, double *v, double *wind_speed, double *wind_angle)
+    void compute_wind_speed_angle(Grid.DimStruct *dims, double *u, double *v, double *wind_speed, double *wind_angle, double u0, double v0)
 
 cdef class Kinematics:
     def __init__(self):
@@ -50,7 +51,7 @@ cdef class Kinematics:
 
         return
 
-    cpdef stats_io(self, Grid.Grid Gr, PrognosticVariables.PrognosticVariables PV, NetCDFIO_Stats NS, ParallelMPI.ParallelMPI Pa):
+    cpdef stats_io(self, Grid.Grid Gr, ReferenceState.ReferenceState RS, PrognosticVariables.PrognosticVariables PV, NetCDFIO_Stats NS, ParallelMPI.ParallelMPI Pa):
         cdef:
             double [:] mean = np.empty((Gr.dims.nlg[2],),dtype=np.double,order='c')
             Py_ssize_t u_shift = PV.get_varshift(Gr, 'u')
@@ -60,12 +61,12 @@ cdef class Kinematics:
         mean = Pa.HorizontalMean(Gr,&self.strain_rate_mag[0])
         NS.write_profile('strain_rate_magnitude',mean[Gr.dims.gw:-Gr.dims.gw],Pa)
 
-        compute_wind_speed_angle(&Gr.dims, &PV.values[u_shift], &PV.values[v_shift], &self.wind_speed[0], &self.wind_angle[0])
+        compute_wind_speed_angle(&Gr.dims, &PV.values[u_shift], &PV.values[v_shift], &self.wind_speed[0], &self.wind_angle[0],RS.u0, RS.v0)
 
         mean = Pa.HorizontalMean(Gr,&self.wind_speed[0])
         NS.write_profile('wind_speed',mean[Gr.dims.gw:-Gr.dims.gw],Pa)
 
-        mean = Pa.HorizontalMean(Gr,&self.wind_speed[0])
+        mean = Pa.HorizontalMean(Gr,&self.wind_angle[0])
         NS.write_profile('wind_angle',mean[Gr.dims.gw:-Gr.dims.gw],Pa)
 
         return

diff --git a/SGS.pyx b/SGS.pyx
@@ -17,6 +17,7 @@ cdef extern from "sgs.h":
     void tke_shear_production(Grid.DimStruct *dims,  double* e_tendency, double* visc, double* strain_rate_mag)
     void tke_buoyant_production(Grid.DimStruct *dims,  double* e_tendency, double* diff, double* buoy_freq)
     void tke_surface(Grid.DimStruct *dims, double* e, double* lmo, double* ustar, double h_bl, double zb) nogil
+    double tke_ell(double cn, double e, double buoy_freq, double delta) nogil
 
 cdef class SGS:
     def __init__(self,namelist):
@@ -151,6 +152,7 @@ cdef class TKE:
         NS.add_profile('tke_shear_tendency', Gr, Pa)
         NS.add_profile('tke_buoyancy_tendency', Gr, Pa)
         NS.add_profile('tke_prandtl_number', Gr, Pa)
+        NS.add_profile('tke_mixing_length', Gr, Pa)
 
         return
 
@@ -221,7 +223,7 @@ cdef class TKE:
             Py_ssize_t e_shift = PV.get_varshift(Gr,'e')
             double [:] tmp_tendency  = np.zeros((Gr.dims.npg),dtype=np.double,order='c')
             double [:] mean_tendency = np.empty((Gr.dims.nlg[2],),dtype=np.double,order='c')
-            double [:] prt  = np.zeros((Gr.dims.npg),dtype=np.double,order='c')
+
             double [:] mean = np.empty((Gr.dims.nlg[2],),dtype=np.double,order='c')
 
 
@@ -240,13 +242,21 @@ cdef class TKE:
         cdef:
             Py_ssize_t i
             Py_ssize_t npg = Gr.dims.npg
+            double delta = (Gr.dims.dx[0] * Gr.dims.dx[1] * Gr.dims.dx[2])**(1.0/3.0)
+            double [:] prt  = np.zeros((Gr.dims.npg),dtype=np.double,order='c')
+            double [:] mixing_length = np.zeros((Gr.dims.npg),dtype=np.double,order='c')
 
         with nogil:
             for i in xrange(npg):
-                prt[i] = DV.values[visc_shift + i]/(DV.values[diff_shift + i] + 1.0e-20)
+                mixing_length[i] = tke_ell(self.cn, PV.values[e_shift+i], DV.values[bf_shift+i], delta)
+                prt[i] = delta/(delta + 2.0*mixing_length[i])
 
 
         mean = Pa.HorizontalMean(Gr,&prt[0])
         NS.write_profile('tke_prandtl_number',mean[Gr.dims.gw:-Gr.dims.gw],Pa)
 
+        mean = Pa.HorizontalMean(Gr,&mixing_length[0])
+        NS.write_profile('tke_mixing_length',mean[Gr.dims.gw:-Gr.dims.gw],Pa)
+
+
         return
diff --git a/Simulation3d.pyx b/Simulation3d.pyx
@@ -158,7 +158,7 @@ class Simulation3d:
                 self.Th.stats_io(self.Gr, self.PV, self.StatsIO, self.Pa)
                 self.Sur.stats_io(self.Gr, self.StatsIO, self.Pa)
                 self.SGS.stats_io(self.Gr,self.DV,self.PV,self.Ke,self.StatsIO,self.Pa)
-                self.Ke.stats_io(self.Gr,self.PV,self.StatsIO,self.Pa)
+                self.Ke.stats_io(self.Gr,self.Ref,self.PV,self.StatsIO,self.Pa)
         return
 
     def force_io(self):
@@ -171,6 +171,6 @@ class Simulation3d:
         self.Th.stats_io(self.Gr, self.PV, self.StatsIO, self.Pa)
         self.Sur.stats_io(self.Gr, self.StatsIO, self.Pa)
         self.SGS.stats_io(self.Gr,self.DV,self.PV,self.Ke,self.StatsIO,self.Pa)
-        self.Ke.stats_io(self.Gr,self.PV,self.StatsIO,self.Pa)
+        self.Ke.stats_io(self.Gr,self.Ref,self.PV,self.StatsIO,self.Pa)
         return