Update Simulation3d.pyx

Csrc/kinematics.h

@@ -111,7 +111,7 @@ void compute_strain_rate_mag(const struct DimStruct *dims, double* restrict stra
     //Off-diagonal components
     //Here factor of 2 arises because we invoke symmetry of tensor
     for(ssize_t vi1=0;vi1<dims->dims-1;vi1++){
-        for (ssize_t d=vi1;d<dims->dims;d++){
+        for (ssize_t d=vi1+1;d<dims->dims;d++){
             const ssize_t shift_s = 3 * dims->npg * vi1 + dims->npg * d ;
             for(ssize_t i=imin+1;i<imax;i++){
                 const ssize_t ishift = i*istride;

Csrc/microphysics_arctic_1m.h

@@ -341,20 +341,27 @@ void accretion_all(double density, double p0, double temperature, double ccn, do
         factor_s = density*GD3_SNOW*nsnow*pi*ALPHA_ACC_SNOW*C_SNOW*0.25*pow(snow_lam, (-D_SNOW-3.0));
     }
 
-    double src_ri = -piacr/density;
-    double src_rl = factor_r*e_rl*ql/density;
-    double src_si = factor_s*e_si*qi/density;
+    double src_ri_r = 0.0;
+    double src_ri_i = 0.0;
+    double src_ri_s = 0.0;
+    double src_sl_r = 0.0;
+    double src_sl_s = 0.0;
+    double src_rl_r = factor_r*e_rl*ql/density;
+    double src_si_s = factor_s*e_si*qi/density;
     double rime_sl = factor_s*e_sl*ql/density;
-    double src_sl = 0.0;
 
     if( temperature > 273.16 ){
-        src_sl = -cvl/lf0*(temperature-273.16)*rime_sl;
-        src_rl = src_rl + rime_sl - src_sl;
+        src_sl_s = -cvl/lf0*(temperature-273.16)*rime_sl;
+        src_sl_r = (1.0 + cvl/lf0*(temperature-273.16))*rime_sl;
     }
     else{
-        src_sl = rime_sl + (factor_r*e_ri*qi + piacr)/density;
+        src_ri_i = factor_r*e_ri*qi/density;
+        src_ri_r = -piacr/density;
+        src_sl_s = rime_sl;
     }
 
+    src_ri_s = -src_ri_r + src_ri_i; //both rain-ice collision terms contribute to snow tendency
+
     /* Now precip-precip interactions */
     double src_r = 0.0;
     double src_s = 0.0;
@@ -376,10 +383,10 @@ void accretion_all(double density, double p0, double temperature, double ccn, do
         }
     }
 
-    *qrain_tendency = src_r + src_rl + src_ri;
-    *qsnow_tendency = src_s + src_sl + src_si;
-    *ql_tendency = -(src_rl + rime_sl);
-    *qi_tendency = -(src_ri + src_si);
+    *qrain_tendency = src_r + src_rl_r + src_ri_r + src_sl_r;
+    *qsnow_tendency = src_s + src_sl_s + src_ri_s + src_si_s;
+    *ql_tendency = -src_rl_r - src_sl_s - src_sl_r;
+    *qi_tendency = -src_ri_i - src_si_s;
 
     return;
 };
@@ -1103,3 +1110,90 @@ void melt_snow_wrapper(const struct DimStruct *dims, double* density, double* te
     return;
 
 };
+
+void get_n_ice(const struct DimStruct *dims, double* density, double* qi, double n0i, double* n_ice){
+
+    const ssize_t istride = dims->nlg[1] * dims->nlg[2];
+    const ssize_t jstride = dims->nlg[2];
+    const ssize_t imin = dims->gw;
+    const ssize_t jmin = dims->gw;
+    const ssize_t kmin = dims->gw;
+    const ssize_t imax = dims->nlg[0]-dims->gw;
+    const ssize_t jmax = dims->nlg[1]-dims->gw;
+    const ssize_t kmax = dims->nlg[2]-dims->gw;
+
+    for(ssize_t i=imin; i<imax; i++){
+        const ssize_t ishift = i * istride;
+        for(ssize_t j=jmin; j<jmax; j++){
+            const ssize_t jshift = j * jstride;
+            for(ssize_t k=kmin; k<kmax; k++){
+                const ssize_t ijk = ishift + jshift + k;
+
+                double lambda_ = ice_lambda(density[ijk], qi[ijk], n0i);
+
+                n_ice[ijk] = n0i / lambda_;
+
+            }
+        }
+    }
+    return;
+
+};
+
+void get_n_rain(const struct DimStruct *dims, double* density, double* qrain, double* n0rain, double* n_rain){
+
+    const ssize_t istride = dims->nlg[1] * dims->nlg[2];
+    const ssize_t jstride = dims->nlg[2];
+    const ssize_t imin = dims->gw;
+    const ssize_t jmin = dims->gw;
+    const ssize_t kmin = dims->gw;
+    const ssize_t imax = dims->nlg[0]-dims->gw;
+    const ssize_t jmax = dims->nlg[1]-dims->gw;
+    const ssize_t kmax = dims->nlg[2]-dims->gw;
+
+    for(ssize_t i=imin; i<imax; i++){
+        const ssize_t ishift = i * istride;
+        for(ssize_t j=jmin; j<jmax; j++){
+            const ssize_t jshift = j * jstride;
+            for(ssize_t k=kmin; k<kmax; k++){
+                const ssize_t ijk = ishift + jshift + k;
+
+                double lambda_ = rain_lambda(density[ijk], qrain[ijk], n0rain[ijk]);
+
+                n_rain[ijk] = n0rain[ijk] / lambda_;
+
+            }
+        }
+    }
+    return;
+
+};
+
+void get_n_snow(const struct DimStruct *dims, double* density, double* qsnow, double* n0snow, double* n_snow){
+
+    const ssize_t istride = dims->nlg[1] * dims->nlg[2];
+    const ssize_t jstride = dims->nlg[2];
+    const ssize_t imin = dims->gw;
+    const ssize_t jmin = dims->gw;
+    const ssize_t kmin = dims->gw;
+    const ssize_t imax = dims->nlg[0]-dims->gw;
+    const ssize_t jmax = dims->nlg[1]-dims->gw;
+    const ssize_t kmax = dims->nlg[2]-dims->gw;
+
+    for(ssize_t i=imin; i<imax; i++){
+        const ssize_t ishift = i * istride;
+        for(ssize_t j=jmin; j<jmax; j++){
+            const ssize_t jshift = j * jstride;
+            for(ssize_t k=kmin; k<kmax; k++){
+                const ssize_t ijk = ishift + jshift + k;
+
+                double lambda_ = snow_lambda(density[ijk], qsnow[ijk], n0snow[ijk]);
+
+                n_snow[ijk] = n0snow[ijk] / lambda_;
+
+            }
+        }
+    }
+    return;
+
+};

Forcing.pxd

@@ -129,6 +129,10 @@ cdef class ForcingIsdacCC:
         double [:] nudge_coeff_scalars
         double [:] w_half
         double z_top
+        double [:] source_u_nudge
+        double [:] source_v_nudge
+        double [:] source_s_nudge
+        double [:] source_qt_nudge
 
         double divergence
     cpdef initialize(self, Grid.Grid Gr, ReferenceState.ReferenceState RS, Th, NetCDFIO_Stats NS, ParallelMPI.ParallelMPI Pa)

Forcing.pyx

@@ -822,7 +822,10 @@ cdef class ForcingIsdacCC:
         self.initial_u = np.zeros(Gr.dims.nlg[2],dtype=np.double,order='c')
         self.initial_v = np.zeros(Gr.dims.nlg[2],dtype=np.double,order='c')
         self.w_half =  np.zeros(Gr.dims.nlg[2],dtype=np.double,order='c')
-
+        self.source_u_nudge = np.zeros(Gr.dims.nlg[2],dtype=np.double,order='c')
+        self.source_v_nudge = np.zeros(Gr.dims.nlg[2],dtype=np.double,order='c')
+        self.source_s_nudge = np.zeros(Gr.dims.nlg[2],dtype=np.double,order='c')
+        self.source_qt_nudge = np.zeros(Gr.dims.nlg[2],dtype=np.double,order='c')
         cdef:
             Py_ssize_t k
 
@@ -837,8 +840,8 @@ cdef class ForcingIsdacCC:
         for k in xrange(Gr.dims.nlg[2]):
 
             #Set velocity profile
-            self.initial_v[k] = -2.0 + 0.003 * Gr.zl_half[k]
-            self.initial_u[k] = -7.0
+            self.initial_v[k] = -2.0 + 0.003 * Gr.zl_half[k] - RS.v0
+            self.initial_u[k] = -7.0 - RS.u0
 
             if Gr.zl_half[k] <= self.z_top:
                 self.nudge_coeff_velocities[k] = (1/7200.0)*0.5*(1 - cos(pi*Gr.zl_half[k]/self.z_top))
@@ -877,17 +880,53 @@ cdef class ForcingIsdacCC:
             Py_ssize_t s_shift = PV.get_varshift(Gr, 's')
             Py_ssize_t qt_shift = PV.get_varshift(Gr,'qt')
             Py_ssize_t t_shift = DV.get_varshift(Gr, 'temperature')
+            Py_ssize_t imin = Gr.dims.gw
+            Py_ssize_t jmin = Gr.dims.gw
+            Py_ssize_t kmin = Gr.dims.gw
+            Py_ssize_t imax = Gr.dims.nlg[0] - Gr.dims.gw
+            Py_ssize_t jmax = Gr.dims.nlg[1] - Gr.dims.gw
+            Py_ssize_t kmax = Gr.dims.nlg[2] - Gr.dims.gw
+            Py_ssize_t istride = Gr.dims.nlg[1] * Gr.dims.nlg[2]
+            Py_ssize_t jstride = Gr.dims.nlg[2]
+            Py_ssize_t i, j, k, ishift, jshift, ijk
 
         apply_subsidence(&Gr.dims,&RS.rho0[0],&RS.rho0_half[0],&self.w_half[0],&PV.values[s_shift],&PV.tendencies[s_shift])
         apply_subsidence(&Gr.dims,&RS.rho0[0],&RS.rho0_half[0],&self.w_half[0],&PV.values[qt_shift],&PV.tendencies[qt_shift])
         apply_subsidence(&Gr.dims,&RS.rho0[0],&RS.rho0_half[0],&self.w_half[0],&PV.values[u_shift],&PV.tendencies[u_shift])
         apply_subsidence(&Gr.dims,&RS.rho0[0],&RS.rho0_half[0],&self.w_half[0],&PV.values[v_shift],&PV.tendencies[v_shift])
 
-        apply_nudging(&Gr.dims,&self.nudge_coeff_velocities[0],&self.initial_u[0],&PV.values[u_shift],&PV.tendencies[u_shift])
-        apply_nudging(&Gr.dims,&self.nudge_coeff_velocities[0],&self.initial_v[0],&PV.values[v_shift],&PV.tendencies[v_shift])
-        apply_nudging(&Gr.dims,&self.nudge_coeff_scalars[0],&self.initial_entropy[0],&PV.values[s_shift],&PV.tendencies[s_shift])
-        apply_nudging(&Gr.dims,&self.nudge_coeff_scalars[0],&self.initial_qt[0],&PV.values[qt_shift],&PV.tendencies[qt_shift])
+        # apply_nudging(&Gr.dims,&self.nudge_coeff_velocities[0],&self.initial_u[0],&PV.values[u_shift],&PV.tendencies[u_shift])
+        # apply_nudging(&Gr.dims,&self.nudge_coeff_velocities[0],&self.initial_v[0],&PV.values[v_shift],&PV.tendencies[v_shift])
+        # apply_nudging(&Gr.dims,&self.nudge_coeff_scalars[0],&self.initial_entropy[0],&PV.values[s_shift],&PV.tendencies[s_shift])
+        # apply_nudging(&Gr.dims,&self.nudge_coeff_scalars[0],&self.initial_qt[0],&PV.values[qt_shift],&PV.tendencies[qt_shift])
 
+        #Get domain mean profiles
+        cdef:
+            double [:] umean = Pa.HorizontalMean(Gr, &PV.values[u_shift])
+            double [:] vmean = Pa.HorizontalMean(Gr, &PV.values[v_shift])
+            double [:] smean = Pa.HorizontalMean(Gr, &PV.values[s_shift])
+            double [:] qtmean = Pa.HorizontalMean(Gr, &PV.values[qt_shift])
+
+        # Calculate nudging
+        with nogil:
+            for k in xrange(kmin, kmax):
+                # Nudge mean wind profiles through entire depth
+                self.source_u_nudge[k] = -(umean[k] - self.initial_u[k]) * self.nudge_coeff_velocities[k]
+                self.source_v_nudge[k] = -(vmean[k] - self.initial_v[k]) * self.nudge_coeff_velocities[k]
+                self.source_s_nudge[k] = -(smean[k] - self.initial_entropy[k]) * self.nudge_coeff_scalars[k]
+                self.source_qt_nudge[k] = -(qtmean[k] - self.initial_qt[k]) * self.nudge_coeff_scalars[k]
+
+        with nogil:
+            for i in xrange(imin,imax):
+                ishift = i * istride
+                for j in xrange(jmin,jmax):
+                    jshift = j * jstride
+                    for k in xrange(kmin,kmax):
+                        ijk = ishift + jshift + k
+                        PV.tendencies[u_shift + ijk] += self.source_u_nudge[k]
+                        PV.tendencies[v_shift + ijk] += self.source_v_nudge[k]
+                        PV.tendencies[s_shift + ijk] += self.source_s_nudge[k]
+                        PV.tendencies[qt_shift + ijk] += self.source_qt_nudge[k]
 
         return
 
@@ -928,29 +967,29 @@ cdef class ForcingIsdacCC:
         mean_tendency = Pa.HorizontalMean(Gr,&tmp_tendency[0])
         NS.write_profile('v_subsidence_tendency',mean_tendency[Gr.dims.gw:-Gr.dims.gw],Pa)
 
-        tmp_tendency[:] = 0.0
-        apply_nudging(&Gr.dims,&self.nudge_coeff_velocities[0],&self.initial_u[0],&PV.values[u_shift],
-                      &tmp_tendency[0])
-        mean_tendency = Pa.HorizontalMean(Gr,&tmp_tendency[0])
-        NS.write_profile('u_nudging_tendency',mean_tendency[Gr.dims.gw:-Gr.dims.gw],Pa)
-
-        tmp_tendency[:] = 0.0
-        apply_nudging(&Gr.dims,&self.nudge_coeff_velocities[0],&self.initial_v[0],&PV.values[v_shift],
-                      &tmp_tendency[0])
-        mean_tendency = Pa.HorizontalMean(Gr,&tmp_tendency[0])
-        NS.write_profile('v_nudging_tendency',mean_tendency[Gr.dims.gw:-Gr.dims.gw],Pa)
-
-        tmp_tendency[:] = 0.0
-        apply_nudging(&Gr.dims,&self.nudge_coeff_scalars[0],&self.initial_entropy[0],&PV.values[s_shift],
-                      &tmp_tendency[0])
-        mean_tendency = Pa.HorizontalMean(Gr,&tmp_tendency[0])
-        NS.write_profile('s_nudging_tendency',mean_tendency[Gr.dims.gw:-Gr.dims.gw],Pa)
-
-        tmp_tendency[:] = 0.0
-        apply_nudging(&Gr.dims,&self.nudge_coeff_scalars[0],&self.initial_qt[0],&PV.values[qt_shift],
-                      &tmp_tendency[0])
-        mean_tendency = Pa.HorizontalMean(Gr,&tmp_tendency[0])
-        NS.write_profile('qt_nudging_tendency',mean_tendency[Gr.dims.gw:-Gr.dims.gw],Pa)
+        # tmp_tendency[:] = 0.0
+        # apply_nudging(&Gr.dims,&self.nudge_coeff_velocities[0],&self.initial_u[0],&PV.values[u_shift],
+        #               &tmp_tendency[0])
+        # mean_tendency = Pa.HorizontalMean(Gr,&tmp_tendency[0])
+        NS.write_profile('u_nudging_tendency',self.source_u_nudge[Gr.dims.gw:-Gr.dims.gw],Pa)
+        #
+        # tmp_tendency[:] = 0.0
+        # apply_nudging(&Gr.dims,&self.nudge_coeff_velocities[0],&self.initial_v[0],&PV.values[v_shift],
+        #               &tmp_tendency[0])
+        # mean_tendency = Pa.HorizontalMean(Gr,&tmp_tendency[0])
+        NS.write_profile('v_nudging_tendency',self.source_v_nudge[Gr.dims.gw:-Gr.dims.gw],Pa)
+
+        # tmp_tendency[:] = 0.0
+        # apply_nudging(&Gr.dims,&self.nudge_coeff_scalars[0],&self.initial_entropy[0],&PV.values[s_shift],
+        #               &tmp_tendency[0])
+        # mean_tendency = Pa.HorizontalMean(Gr,&tmp_tendency[0])
+        NS.write_profile('s_nudging_tendency',self.source_s_nudge[Gr.dims.gw:-Gr.dims.gw],Pa)
+
+        # tmp_tendency[:] = 0.0
+        # apply_nudging(&Gr.dims,&self.nudge_coeff_scalars[0],&self.initial_qt[0],&PV.values[qt_shift],
+        #               &tmp_tendency[0])
+        # mean_tendency = Pa.HorizontalMean(Gr,&tmp_tendency[0])
+        NS.write_profile('qt_nudging_tendency',self.source_qt_nudge[Gr.dims.gw:-Gr.dims.gw],Pa)
 
         return
 

Initialization.pyx

@@ -1231,62 +1231,44 @@ def InitSheba(namelist, Grid.Grid Gr,PrognosticVariables.PrognosticVariables PV,
         Py_ssize_t qt_varshift = PV.get_varshift(Gr,'qt')
         double [:] thetal = np.zeros((Gr.dims.nlg[2],),dtype=np.double,order='c')
         double [:] qt = np.zeros((Gr.dims.nlg[2],),dtype=np.double,order='c')
-        double [:] wt = np.zeros((Gr.dims.nlg[2],),dtype=np.double,order='c')
         double p_inv = 95700.0
-        double [:] ps
-        double [:] us
-        double [:] vs
         double [:] v = np.zeros((Gr.dims.nlg[2],),dtype=np.double,order='c')
         double [:] u = np.zeros((Gr.dims.nlg[2],),dtype=np.double,order='c')
 
-    for k in xrange(Gr.dims.nlg[2]):
 
-        #Set thetal and qt profile
-        if RS.p0_half[k] > p_inv:
-            thetal[k] = 255.7 #257.0
-            wt[k] = 0.91 #0.915 #Mixing ratio in g/kg
-        else:
-            thetal[k] = 255.7+6.9 #263.9
-            wt[k] = 0.8
+    sheba_init = np.loadtxt('./SHEBAdata/sheba_init.txt', comments='!')
+    sheba_p = sheba_init[:, 0]
+    sheba_t = sheba_init[:, 1]
+    sheba_rh = sheba_init[:, 2]
+    sheba_rl = sheba_init[:, 3]
+    sheba_thetal = sheba_t*(p_tilde/sheba_p)**(Rd/cpd)*np.exp(-(sheba_rl*2.501e6) / (sheba_t*cpd))
+    sheba_u = sheba_init[:, 4]
+    sheba_v = sheba_init[:, 5]
 
-    for k in xrange(Gr.dims.nlg[2]-1):
-        if RS.p0_half[k] < p_inv:
-            thetal[k+1] = thetal[k] + ((1.0/exner_c(RS.p0_half[k+1]))*fmin(3.631e-8*(p_inv - RS.p0_half[k+1]), 5.7e-4))*(RS.p0_half[k]-RS.p0_half[k+1])
-            wt[k+1] = wt[k] - 1.4e-5*(RS.p0_half[k]-RS.p0_half[k+1])
-
-        #Convert mixing ratio to specific humidity
-        qt[k] = wt[k]/(1.0 + wt[k]/1000.0)/1000.0
-
-    #Horizontal wind profiles from Colleen's JPLLES code
-
-    ps = np.array([1017.0,1012.0,1007.0,1002.0,997.0,992.0,987.0,982.0,977.0,972.0,967.0,962.0,957.01,957.0,
-                             956.99,952.0,947.0,942.0,937.0,932.0,927.0,922.0,917.0,912.0,907.0,902.0,897.0,892.0,
-                             887.0,882.0,877.0,872.0,867.0,862.0,857.0,852.0,847.0,842.0,837.0,832.0,827.0,822.0,817.0,
-                             812.0,807.0,802.0,797.0,792.0,787.0,782.0,777.0,772.0,767.0,762.0,757.0,752.0,747.0,742.0,
-                             737.0,732.0,727.0,722.0,717.0,712.0,707.0,702.0,697.0,692.0,687.0,682.0,677.0,672.0,667.0,
-                             662.0,657.0,652.0,647.0,642.0,637.0,632.0,627.0,622.0,617.0,612.0,607.0],dtype=np.double,order='c')*100.0
-
-    us = np.array([2.916,2.8999,2.7895,2.679,2.5568,2.417,2.2772,2.1374,1.9976,1.8856,1.7926,1.6996,1.6066,1.6066,
-                   1.6066,1.5136,1.4205,1.3248,1.2203,1.1159,1.0114,0.90701,0.80257,0.69813,0.59368,0.48564,0.37645,
-                   0.26726,0.15807,0.048881,-0.06031,-0.1695,-0.27869,-0.42739,-0.61744,-0.80748,-0.99753,-1.1876,
-                   -1.3776,-1.5677,-1.7577,-1.9477,-2.0576,-2.1326,-2.2076,-2.2826,-2.3576,-2.4326,-2.5076,-2.5826,
-                   -2.6576,-2.6911,-2.6961,-2.701,-2.7059,-2.7109,-2.7158,-2.7207,-2.7257,-2.7306,-2.7137,-2.653,
-                   -2.5923,-2.5316,-2.4708,-2.4101,-2.3494,-2.2886,-2.2279,-2.1662,-2.08,-1.9938,-1.9076,-1.8214,
-                   -1.7352,-1.649,-1.5629,-1.4767,-1.3905,-1.302,-1.213,-1.124,-1.035,-0.94601,-0.85701],dtype=np.double,order='c')
-
-    vs = np.array([2.8497,2.9023,3.2622,3.6221,3.8701,3.9526,4.0351,4.1177,4.2002,4.2743,4.3427,4.4112,4.4796,4.4796,
-                   4.4796,4.548,4.6165,4.6831,4.744,4.8048,4.8657,4.9266,4.9874,5.0483,5.1092,5.1671,5.2241,5.2811,
-                   5.3381,5.3951,5.4521,5.5091,5.5662,5.6245,5.6844,5.7442,5.8041,5.8639,5.9238,5.9836,6.0434,6.1033,
-                   6.1444,6.1774,6.2105,6.2435,6.2765,6.3096,6.3426,6.3756,6.4086,6.4722,6.5567,6.6412,6.7258,6.8103,
-                   6.8948,6.9794,7.0639,7.1484,7.2388,7.3407,7.4426,7.5446,7.6465,7.7485,7.8504,7.9524,8.0543,8.1574,
-                   8.2893,8.4211,8.553,8.6848,8.8167,8.9485,9.0804,9.2122,9.3441,9.493,9.6463,9.7995,9.9527,10.106,
-                   10.259],dtype=np.double,order='c')
+    sheba_pv = np.zeros_like(sheba_u)
+    for n, t_ in enumerate(sheba_t):
+        sheba_pv[n] = Th.get_pv_star(t_)*sheba_rh[n]/100.0
+    sheba_wt = eps_v * sheba_pv/(sheba_p - sheba_pv) + sheba_rl
+    sheba_qt = sheba_wt/(1.0 + sheba_wt)
+    sheba_ql = sheba_rl/(1.0 + sheba_rl)
 
     #Interpolate to LES grid
     for k in xrange(Gr.dims.nlg[2]):
-        u[k] = interp_pchip(RS.p0[k], ps[::-1], us[::-1])
-        v[k] = interp_pchip(RS.p0[k], ps[::-1], vs[::-1])
-
+        u[k] = interp_pchip(RS.p0_half[k], sheba_p, sheba_u)
+        v[k] = interp_pchip(RS.p0_half[k], sheba_p, sheba_v)
+        thetal[k] = interp_pchip(RS.p0_half[k], sheba_p, sheba_thetal)
+        qt[k] = interp_pchip(RS.p0_half[k], sheba_p, sheba_qt)
+
+    # plt.figure()
+    # plt.subplot(121)
+    # plt.plot(thetal, RS.p0_half)
+    # plt.plot(sheba_thetal, sheba_p)
+    # plt.ylim([110000, 80000])
+    # plt.subplot(122)
+    # plt.plot(qt, RS.p0_half)
+    # plt.plot(sheba_qt, sheba_p)
+    # plt.ylim([110000, 80000])
+    # plt.show()
 
     RS.u0 = 0.5 * (np.amax(u)+np.amin(u))
     RS.v0 = 0.5 * (np.amax(v)+np.amin(v))

Microphysics_Arctic_1M.pxd

@@ -125,6 +125,8 @@ cdef class Microphysics_Arctic_1M:
                  TimeStepping.TimeStepping TS, ParallelMPI.ParallelMPI Pa)
     cpdef stats_io(self, Grid.Grid Gr, ReferenceState.ReferenceState Ref, Th, PrognosticVariables.PrognosticVariables PV,
                    DiagnosticVariables.DiagnosticVariables DV, NetCDFIO_Stats NS, ParallelMPI.ParallelMPI Pa)
+    cpdef micro_fields(self, Grid.Grid Gr, ReferenceState.ReferenceState Ref, PrognosticVariables.PrognosticVariables PV,
+                   DiagnosticVariables.DiagnosticVariables DV)
     cpdef ice_stats(self, Grid.Grid Gr, ReferenceState.ReferenceState Ref, PrognosticVariables.PrognosticVariables PV,
                     DiagnosticVariables.DiagnosticVariables DV, NetCDFIO_Stats NS, ParallelMPI.ParallelMPI Pa)