bug fix in BouLac scheme (parcel energy computation)

PE_CBL_models.py

@@ -899,7 +899,7 @@ class CBL:
                 lup[k] = 0
 
                 # Compute vertical profile of energy perturbations for lifted parcels (from the current level up)
-                energy_profile_upwards = g/theta_0*( (theta[k+1:]+theta[k:-1])*0.5 - theta[k] ) * (zf[k+1:]-zf[k:-1]) - tke[k]
+                energy_profile_upwards = g/theta_0*np.cumsum(( (theta[k+1:]+theta[k:-1])*0.5 - theta[k] ) * (zf[k+1:]-zf[k:-1])) - tke[k]
 
                 # Track the level at which the energy perturbation profile changes sign
                 # from negative (below) to positive (above)
@@ -918,7 +918,7 @@ class CBL:
                 ldown[k] = 0
 
                 # Compute vertical profile of energy perturbations for sinking parcels (from the current level down)
-                energy_profile_downwards = g/theta_0*( (theta[:k]+theta[1:k+1])*0.5 - theta[k] ) * (zf[:k]-zf[1:k+1]) - tke[k]
+                energy_profile_downwards = g/theta_0*np.cumsum(( (theta[:k]+theta[1:k+1])*0.5 - theta[k] ) * (zf[:k]-zf[1:k+1]) )- tke[k]
 
                 # Track the lowermost level at which the energy perturbation profile
                 # is negative (in a stable atmosphere, it's always positive).
@@ -1107,7 +1107,6 @@ class CBL:
 
             # Determine the countergradient correction
             # According to Deardorff 1972 and Therry-Lacarrere 1982 (Eq. 41)
-            is_cgrad = True
             if is_cgrad:
                 tstar = H0/wstar
                 is_pbl = zf<h[j]
@@ -1204,14 +1203,6 @@ class CBL:
                 up        = u[:,j]
                 vp        = v[:,j]
 
-            """
-            print(thetap)
-            print(tkep)
-            if is_bwind:
-                print(up)
-                print(vp)
-            """
-
         # Form state vector
         xb[:nz] = theta[:,-1]
         if is_bwind: