use humid.in and correct unsat profile

tools/drizzle.py

@@ -1,22 +1,24 @@
 import h5py
 from scipy.special import lambertw
-import os
 import afidtools as afid
 import numpy as np
 
-fld = os.environ['WORK']+'/RainyBenard/drizzle_IC'
+fld = '../examples/RainyBenard'
 inputs = afid.InputParams(fld)
 nxm = inputs.nxm
 print(nxm)
 
-saturated_top = False
+mpar = afid.MoistParams(fld)
+alpha = mpar.alpha
+beta = mpar.beta
+gamma = mpar.gamma
 
-alpha = 3.0
-beta = 2.0
-if saturated_top:
-    gamma = beta/(1 - np.exp(-alpha))
-else:
-    gamma = beta
+saturated_top = (mpar.qfixN==1)
+if gamma < 0:
+    if saturated_top:
+        gamma = beta/(1 - np.exp(-alpha))
+    else:
+        gamma = beta
 print(alpha, beta, gamma)
 
 def W(z):
@@ -46,17 +48,18 @@ def W(z):
 
 if not saturated_top:
     # Calculate height above which humidity not at saturation
-    D = alpha*(beta + 1)*np.exp(1-alpha)/(1 + alpha*beta*np.exp(1-alpha))
-    B = beta*D - 1
-    zc = 1 + 1/alpha/(B-beta)
+    A = beta - (1 + gamma)/(1 + gamma*alpha*np.exp(1-alpha))
+    B = alpha*np.exp(1-alpha)*(1 + gamma)/(1 + gamma*alpha*np.exp(1-alpha))
+    zc = 1 - (1 + gamma*alpha*np.exp(1-alpha))/alpha/(1+gamma)
     # Prescribe linear profiles above this level
-    b[z > zc] = beta - 1 + B*(z[z > zc] - 1)
-    q[z > zc] = D*(1 - z[z > zc])
+    b[z > zc] = beta - 1 + A*(z[z > zc] - 1)
+    q[z > zc] = B*(1 - z[z > zc])
 
 import matplotlib.pyplot as plt
 fig, ax = plt.subplots(layout='constrained')
 ax.plot(b, z)
 ax.plot(q, z)
+ax.grid(True)
 fig.savefig('drizzle.png')
 
 with h5py.File(fld+'/drizzle.h5','w') as f: