Merge bug fixes on convective adjustment and SW water vapor absorption.

climlab/convection/convadj.py

@@ -102,57 +102,124 @@ def Akamaev_adjustment_multidim(theta, q, beta, n_k, theta_k, s_k, t_k):
 def Akamaev_adjustment(theta, q, beta, n_k, theta_k, s_k, t_k):
     '''Single column only.'''
     L = q.size  # number of vertical levels
-    k = 0
-    n_k[0] = 1
-    theta_k[0] = theta[0]
-
-    for l in range(1, L):
+    # Akamaev step 1
+    k = 1
+    n_k[k-1] = 1
+    theta_k[k-1] = theta[k-1]
+    l = 2
+    while True:
+        # Akamaev step 2
         n = 1
-        thistheta = theta[l]
-        done = False
-        while not done:
-            if (theta_k[k] > thistheta):
-                s = 0.
-                t = 0.
-                # stratification is unstable
-                if n == 1:
-                    # current layer is not an earlier-formed neutral layer
-                    s = q[l]
-                    t = s * thistheta
-                if (n_k[k] < 2):
-                    # lower adjacent level is not an earlier-formed neutral layer
-                    s_k[k] = q[l-n]
-                    t_k[k] = s_k[k] * theta_k[k]
-                #  join current and underlying layers
-                n += n_k[k]
-                s += s_k[k]
-                s_k[k] = s
-                t += t_k[k]
-                t_k[k] = t
-                thistheta = t / s
-                if k == 0:
-                    # joint neutral layer in the first one, done checking lower layers
-                    done = True
-                else:
-                    k -= 1
-                    # go back and check stability of the lower adjacent layer
+        thistheta = theta[l-1]
+        while True:
+            # Akamaev step 3
+            if theta_k[k-1] <= thistheta:
+                # Akamaev step 6
+                k += 1
+                break  # to step 7
             else:
-                k += 1  # statification is stable
-                done = True
-        # if l < L-1:
-        n_k[k] = n
-        theta_k[k] = thistheta
-    #  finished looping through to check stability
+                if n <= 1:
+                    s = q[l-1]
+                    t = s*thistheta
+                # Akamaev step 4
+                if n_k[k-1] <= 1:
+                # lower adjacent level is not an earlier-formed neutral layer
+                    s_k[k-1] = q[l-n-1]
+                    t_k[k-1] = s_k[k-1] * theta_k[k-1]
+                # Akamaev step 5
+                    #  join current and underlying layers
+                n += n_k[k-1]
+                s += s_k[k-1]
+                t += t_k[k-1]
+                s_k[k-1] = s
+                t_k[k-1] = t
+                thistheta = t/s
+                if k==1:
+                    #  joint neutral layer is the first one
+                    break  # to step 7
+                k -= 1
+                # back to step 3
+        # Akamaev step 7
+        if l == L:  # the scan is over
+            break  # to step 8
+        l += 1
+        n_k[k-1] = n
+        theta_k[k-1] = thistheta
+        # back to step 2
+
+
+#
+#    for l in range(2, L+1):  # l = 2,3,4,...,L
+#        n = 1
+#        thistheta = theta[l-1]
+#        done = False
+#        while not done:
+#            if (theta_k[k-1] > thistheta):
+#                s = 0.
+#                t = 0.
+#                # stratification is unstable
+#                if n == 1:
+#                    # current layer is not an earlier-formed neutral layer
+#                    s = q[l-1]
+#                    t = s * thistheta
+#                if (n_k[k-1] < 2):
+#                    # lower adjacent level is not an earlier-formed neutral layer
+#                    s_k[k-1] = q[l-n-1]
+#                    t_k[k-1] = s_k[k-1] * theta_k[k-1]
+#                #  join current and underlying layers
+#                n += n_k[k-1]
+#                s += s_k[k-1]
+#                s_k[k-1] = s
+#                t += t_k[k-1]
+#                t_k[k-1] = t
+#                thistheta = t / s
+#                if k == 1-1:
+#                    # joint neutral layer in the first one, done checking lower layers
+#                    done = True
+#                else:
+#                    k -= 1
+#                    # go back and check stability of the lower adjacent layer
+#            else:
+#                k += 1  # statification is stable
+#                done = True
+#        if l < L:
+#            n_k[k-1] = n
+#            theta_k[k-1] = thistheta
+#    #  finished looping through to check stability
 
     # update the potential temperatures
-    for l in range(L-1, -1, -1):
-        theta[l] = thistheta
-        if n > 1:
-            n -= 1
-        else:
-            k -= 1
-            n = n_k[k]
-            thistheta = theta_k[k]
+    while True:
+        while True:
+            # Akamaev step 8
+            if n==1: # current model level was not included in any neutral layer
+                break  # to step 11
+            while True:
+                #  Akamaev step 9
+                theta[l-1] = thistheta
+                if n==1:
+                    break
+                #  Akamaev step 10
+                l -= 1
+                n -= 1
+                #  back to step 9
+        # Akamaev step 11
+        if k==1:
+            break
+        k -= 1
+        l -= 1
+        n = n_k[k-1]
+        thistheta = theta_k[k-1]
+        # back to step 8
+
+#    for l in range(L, 0, -1):
+#        if n > 1:
+#        while n>1:
+#            theta[l-1] = thistheta
+#            n -= 1
+#        else:
+#            k -= 1
+#            n = n_k[k]
+#            thistheta = theta_k[k]
     return theta
 
 #  Attempt to use numba to compile the Akamaev_adjustment function
@@ -162,5 +229,6 @@ def Akamaev_adjustment(theta, q, beta, n_k, theta_k, s_k, t_k):
 try:
     from numba import jit
     Akamaev_adjustment = jit(signature_or_function=Akamaev_adjustment)
+    #print 'Compiling Akamaev_adjustment() with numba.'
 except:
     pass

climlab/model/column.py

@@ -14,7 +14,7 @@
 from climlab.radiation.radiation import Radiation, RadiationSW
 from climlab.convection.convadj import ConvectiveAdjustment
 from climlab.surface.surface_radiation import SurfaceRadiation
-from climlab.radiation.nband import ThreeBandSW, FourBandLW
+from climlab.radiation.nband import ThreeBandSW, FourBandLW, FourBandSW
 from climlab.radiation.water_vapor import ManabeWaterVapor
 
 class GreyRadiationModel(TimeDependentProcess):

climlab/radiation/nband.py

@@ -151,10 +151,6 @@ def __init__(self, **kwargs):
             channel 2 is remaining radiation (72%)
         '''
         super(ThreeBandSW, self).__init__(**kwargs)
-        #  Three SW channels:
-        # channel 0 is Hartley and Huggins band (UV, 1%, 200 - 340 nm)
-        # channel 1 is Chappuis band (27%, 450 - 800 nm)
-        # channel 2 is remaining radiation (72%)
         #   fraction of the total solar flux in each band:
         self.band_fraction = np.array([0.01, 0.27, 0.72])
         if 'CO2' not in self.absorber_vmr:
@@ -167,7 +163,8 @@ def __init__(self, **kwargs):
         O3 = np.array([200.E-24, 0.285E-24, 0.]) * const.Rd / const.kBoltzmann
         #self.absorption_cross_section['O3'] = np.reshape(O3,
         #    (self.num_channels, 1))
-        H2O = np.array([0.002, 0.002, 0.002])
+        #H2O = np.array([0.002, 0.002, 0.002])
+        H2O = np.array([0., 0., 0.001])
         for n in range(self.Tatm.domain.numdims):
             H2O = H2O[:, np.newaxis]
             O3 = O3[:, np.newaxis]
@@ -183,6 +180,57 @@ def __init__(self, **kwargs):
     def emissivity(self):
         # This ensures that emissivity is always zero for shortwave classes
         return np.zeros_like(self.absorptivity)
+
+
+class FourBandSW(NbandRadiation):
+    #  The most recent RMCM program uses a four-channel SW model
+    #  this is probably a better way to do it... distinguishes between
+    #  visible band with no absorption and near-infrared with weak H2O absorption
+    #   But this needs some tuning and better documentation
+    def __init__(self, **kwargs):
+        '''A four-band mdoel for shortwave radiation.
+    
+        The spectral decomposition used here is largely based on the
+        "Moist Radiative-Convective Model" by Aarnout van Delden, Utrecht University
+        a.j.vandelden@uu.nl
+        http://www.staff.science.uu.nl/~delde102/RCM.htm
+
+        Four SW channels:
+            channel 0 is Hartley and Huggins band (UV, 6%, <340 nm)
+            channel 1 is part of visible with no O3 absorption (14%, 340 - 500 nm)
+            channel 2 is Chappuis band (27%, 500 - 700 nm)
+            channel 3 is near-infrared (53%, > 700 nm)
+        '''
+        super(FourBandSW, self).__init__(**kwargs)
+        #   fraction of the total solar flux in each band:
+        self.band_fraction = np.array([0.06, 0.14, 0.27, 0.53])
+        if 'CO2' not in self.absorber_vmr:
+            self.absorber_vmr['CO2'] = 380.E-6 * np.ones_like(self.Tatm)
+        if 'O3' not in self.absorber_vmr:
+            self.absorber_vmr['O3'] = np.zeros_like(self.Tatm)
+        if 'H2O' not in self.absorber_vmr:
+            self.absorber_vmr['H2O'] = self.q
+        ##  absorption cross-sections in m**2 / kg
+        O3 = np.array([200.E-24, 0., 0.285E-24, 0.]) * const.Rd / const.kBoltzmann
+        #self.absorption_cross_section['O3'] = np.reshape(O3,
+        #    (self.num_channels, 1))
+        H2O = np.array([0., 0., 0., 0.0012])
+        for n in range(self.Tatm.domain.numdims):
+            H2O = H2O[:, np.newaxis]
+            O3 = O3[:, np.newaxis]
+        self.absorption_cross_section['O3'] = O3
+        self.absorption_cross_section['H2O'] = H2O
+        #self.absorption_cross_section['H2O'] = np.reshape(H2O,
+        #    (self.num_channels, 1))
+        self.absorption_cross_section['CO2'] = \
+            np.zeros_like(self.absorption_cross_section['O3'])
+        self.cosZen = 0.5  # cosine of the average solar zenith angle
+
+    @property
+    def emissivity(self):
+        # This ensures that emissivity is always zero for shortwave classes
+        return np.zeros_like(self.absorptivity)
+
 
 
 class FourBandLW(NbandRadiation):

climlab/solar/insolation.py

@@ -105,12 +105,17 @@ def daily_insolation(lat, day, orb=const.orb_present, S0=None, day_type=1):
     oldsettings = np.seterr(invalid='ignore')
     # Compute Ho, the hour angle at sunrise / sunset
     #  Check for no sunrise or no sunset: Berger 1978 eqn (8),(9)
-    Ho = np.where( abs( delta_big ) - np.math.pi / 2. + abs( phi_big ) < 0., np.arccos( -np.tan( phi_big ) * np.tan( delta_big ) ), 
+    Ho = np.where( abs( delta_big ) - np.math.pi / 2. + abs( phi_big ) < 0., 
+                  np.arccos( -np.tan( phi_big ) * np.tan( delta_big ) ), 
             np.where( phi_big * delta_big > 0. , np.math.pi, 0. ) )
+    # this is not really the daily average cosine of the zenith angle...
+    #  it's the integral from sunrise to sunset of that quantity...
+    coszen = (Ho*np.sin(phi_big)*np.sin(delta_big) + 
+              np.cos(phi_big)*np.cos(delta_big)*np.sin(Ho))
     # Compute insolation: Berger 1978 eq (10)
-    Fsw = S0 / np.math.pi * ( ( 1 + ecc * np.cos( lambda_long_big - np.deg2rad(long_peri) ))**2 / (1 - ecc**2)**2 * 
-        ( Ho * np.sin( phi_big ) * np.sin( delta_big ) + np.cos( phi_big ) * np.cos( delta_big ) * np.sin( Ho ) ) )
-
+    Fsw = S0/np.math.pi*( (1. + ecc*np.cos(lambda_long_big - 
+          np.deg2rad(long_peri)))**2 / (1. - ecc**2)**2 * coszen)
+            
     #  Remove singleton dimensions and return
     return np.squeeze( Fsw )