Rewrite EmanuelConvection wrapper so we can use the dry adjustment.

Not working yet!

climlab/__init__.py

@@ -7,7 +7,7 @@
 has been documented for a comprehensive understanding and traceability.
 '''
 
-__version__ = '0.6.4'
+__version__ = '0.6.5.dev0'
 
 # this should ensure that we can still import constants.py as climlab.constants
 from .utils import constants

climlab/convection/_emanuel_convection/Driver.f90

@@ -22,12 +22,13 @@
 
 
 subroutine emanuel_convection(T, Q, QS, U, V, TRA, P, PH, &
-                NCOL,  ND,  NL,   NTRA,   DELT, CBMFold, &
+                NCOL,  ND,  NL,   NTRA,   DELT, IPBL, CBMFold, &
                 CPD, CPV, CL, RV, RD, LV0, G, ROWL, MINORIG, &
                 ELCRIT, TLCRIT, ENTP, SIGD, SIGS, OMTRAIN, OMTSNOW, &
                 COEFFR, COEFFS, CU, BETA, DTMAX, ALPHA, DAMP, &
                 IFLAG,  FT,     FQ,   FU,     FV,  FTRA, &
-                PRECIP, WD,   TPRIME, QPRIME, CBMFnew    )
+                PRECIP, WD,   TPRIME, QPRIME, CBMFnew, &
+                Tout, Qout, QSout, Uout, Vout, TRAout)
 
 ! INPUT
     integer, intent(in) :: ND, NTRA ! number of layers, number of tracers
@@ -36,6 +37,7 @@ subroutine emanuel_convection(T, Q, QS, U, V, TRA, P, PH, &
     real, intent(in) :: T(NCOL,ND),Q(NCOL,ND),QS(NCOL,ND),U(NCOL,ND),V(NCOL,ND)
     real, intent(in) :: TRA(NCOL,ND,NTRA),P(ND),PH(ND+1)
     real, intent(in) :: DELT ! The model time step (sec) between calls to CONVECT
+    real, intent(in) :: IPBL ! switch to bypass the dry adiabatic adjustment (bypass if IPBL=0)
     real, intent(in) :: CBMFold(NCOL) ! The cloud base mass flux ((kg/m**2)/s)
     real, intent(in) :: CPD, CPV, CL, RV, RD, LV0, G, ROWL  ! thermodynamic constants
     integer, intent(in) :: MINORIG  ! index of lowest level from which convection may originate
@@ -47,17 +49,35 @@ subroutine emanuel_convection(T, Q, QS, U, V, TRA, P, PH, &
     real, intent(out) :: FTRA(NCOL,ND,NTRA)
     real, intent(out) :: PRECIP(NCOL), WD(NCOL), TPRIME(NCOL), QPRIME(NCOL)
     real, intent(out) :: CBMFnew(NCOL)
+    real, intent(out) :: Tout(NCOL,ND),Qout(NCOL,ND),QSout(NCOL,ND)
+    real, intent(out) :: Uout(NCOL,ND),Vout(NCOL,ND),TRAout(NCOL,ND,NTRA)
 !  These are not comments! Necessary directives to f2py to handle array dimensions
 !f2py depend(ND) P,PH
 !f2py depend(NCOL,ND) T,Q,QS,U,V
 !f2py depend(NCOL,ND,NTRA) TRA
 !f2py depend(NCOL,ND) FT,FQ,FU,FV
 !f2py depend(NCOL,ND,NTRA) FTRA
 !f2py depend(NCOL) IFLAG, CBMFold, CBMFnew, PRECIP, WD, TPRIME, QPRIME
+!f2py depend(NCOL,ND) Tout,Qout,QSout,Uout,Vout
+!f2py depend(NCOL,ND,NTRA) TRAout
 
     do j = 1, NCOL
+      Tout(j,:) = T(j,:)
+      Qout(j,:) = Q(j,:)
+      QSout(j,:) = QS(j,:)
+      Uout(j,:) = U(j,:)
+      Vout(j,:) = V(j,:)
+      TRAout(j,:,:) = TRA(j,:,:)
+      PRECIP(j) = 0.
+      if (IPBL.NE.0) then
+        !  Local arrays Tj, Qj etc will be altered by this subroutine
+        call DRYADJUSTMENT(Tout(j,:),Qout(j,:),QSout(j,:),Uout(j,:), &
+             Vout(j,:),TRAout(j,:,:),P,PH, &
+             ND, NL, NTRA, DELT, CPD, CPV, CL, RV, RD, LV0, G, ROWL, PRECIP(j))
+      end if
       CBMFnew(j) = 0. + CBMFold(j)  ! will be updated during call to CONVECT
-      call CONVECT(T(j,:),Q(j,:),QS(j,:),U(j,:),V(j,:),TRA(j,:,:),P,PH, &
+      call CONVECT(Tout(j,:),Qout(j,:),QSout(j,:),Uout(j,:),Vout(j,:), &
+             TRAout(j,:,:),P,PH, &
              ND,NL,NTRA,DELT,MINORIG, &
              CPD, CPV, CL, RV, RD, LV0, G, ROWL, &
              ELCRIT, TLCRIT, ENTP, SIGD, SIGS, OMTRAIN, OMTSNOW, &

climlab/convection/_emanuel_convection/convect.f

@@ -265,7 +265,8 @@ SUBROUTINE CONVECT
 C     1    (CPD*(1.-Q(I))+Q(I)*CPV)
 C         TH(I)=T(I)*(1000.0/P(I))**RDCP
 C    7   CONTINUE
-        PRECIP=0.0
+C BRIAN  --  PRECIP is passed in already initialized -- it may be non-zero from dry adjustment
+C        PRECIP=0.0
         WD=0.0
         TPRIME=0.0
         QPRIME=0.0

climlab/convection/_emanuel_convection/dry_adjustment.f

@@ -15,111 +15,115 @@
 *****                          Kerry Emanuel                          *****
 ***************************************************************************
 C
-      SUBROUTINE DRYADJUSTMENT
+        SUBROUTINE DRYADJUSTMENT
      *    (T,   Q,    QS,     U,    V,      TRA,    P,    PH,
      *     ND,  NL,   NTRA,   DELT,
-     *     CPD, CPV, CL, RV, RD, LV0, G, ROWL,
+     *     CPD, CPV, CL, RV, RD, LV0, G, ROWL, PRECIP
      *     )
 C    ***  THE PARAMETER NA SHOULD IN GENERAL BE GREATER THAN   ***
 C    ***                OR EQUAL TO  ND + 1                    ***
 C
-        PARAMETER (NA=70)
-        REAL TH(NA), PRECIP
-
-        PRECIP=0.0
+      PARAMETER (NA=70)
+C
+      REAL T(ND),Q(ND),QS(ND),U(ND),V(ND),TRA(ND,NTRA),P(ND),PH(ND)
+      REAL TRATM(NA)
+      REAL TH(NA),TOLD(NA)
+      REAL LV(NA),LV0,H(NA),HP(NA),GZ(NA),HM(NA)
 
 C  Compute potential temperature
         DO 7 I=1,NL+1
-          RDCP=(RD*(1.-Q(I))+Q(I)*RV)/
-     1     (CPD*(1.-Q(I))+Q(I)*CPV)
-          TH(I)=T(I)*(1000.0/P(I))**RDCP
+         RDCP=(RD*(1.-Q(I))+Q(I)*RV)/
+     1    (CPD*(1.-Q(I))+Q(I)*CPV)
+         TH(I)=T(I)*(1000.0/P(I))**RDCP
     7   CONTINUE
+
 C
 C     ***            PERFORM DRY ADIABATIC ADJUSTMENT            ***
 C
-      JC=0
-      DO 30 I=NL-1,1,-1
-       JN=0
-        SUM=TH(I)*(1.+Q(I)*EPSI-Q(I))
-       DO 10 J=I+1,NL
-        SUM=SUM+TH(J)*(1.+Q(J)*EPSI-Q(J))
-        THBAR=SUM/FLOAT(J+1-I)
-        IF((TH(J)*(1.+Q(J)*EPSI-Q(J))).LT.THBAR)JN=J
-      10    CONTINUE
-       IF(I.EQ.1)JN=MAX(JN,2)
-       IF(JN.EQ.0)GOTO 30
-      12    CONTINUE
-       AHM=0.0
-       RM=0.0
-       UM=0.0
-       VM=0.0
-       DO K=1,NTRA
-        TRATM(K)=0.0
-       END DO
-       DO 15 J=I,JN
-        AHM=AHM+(CPD*(1.-Q(J))+Q(J)*CPV)*T(J)*(PH(J)-PH(J+1))
-        RM=RM+Q(J)*(PH(J)-PH(J+1))
-        UM=UM+U(J)*(PH(J)-PH(J+1))
-        VM=VM+V(J)*(PH(J)-PH(J+1))
-        DO K=1,NTRA
-         TRATM(K)=TRATM(K)+TRA(J,K)*(PH(J)-PH(J+1))
-        END DO
-      15    CONTINUE
-       DPHINV=1./(PH(I)-PH(JN+1))
-       RM=RM*DPHINV
-       UM=UM*DPHINV
-       VM=VM*DPHINV
-       DO K=1,NTRA
-        TRATM(K)=TRATM(K)*DPHINV
-       END DO
-       A2=0.0
-       DO 20 J=I,JN
-        Q(J)=RM
-        U(J)=UM
-        V(J)=VM
-        DO K=1,NTRA
-         TRA(J,K)=TRATM(K)
-        END DO
-        RDCP=(RD*(1.-Q(J))+Q(J)*RV)/
-      1     (CPD*(1.-Q(J))+Q(J)*CPV)
-        X=(0.001*P(J))**RDCP
-        TOLD(J)=T(J)
-        T(J)=X
-        A2=A2+(CPD*(1.-Q(J))+Q(J)*CPV)*X*(PH(J)-PH(J+1))
-      20    CONTINUE
-       DO 25 J=I,JN
-        TH(J)=AHM/A2
-        T(J)=T(J)*TH(J)
-        TC=TOLD(J)-273.15
-        ALV=LV0-CPVMCL*TC
-        QS(J)=QS(J)+QS(J)*(1.+QS(J)*(EPSI-1.))*ALV*(T(J)-
-      1     TOLD(J))/(RV*TOLD(J)*TOLD(J))
-      25    CONTINUE
-       IF((TH(JN+1)*(1.+Q(JN+1)*EPSI-Q(JN+1))).LT.
-      1    (TH(JN)*(1.+Q(JN)*EPSI-Q(JN))))THEN
-        JN=JN+1
-        GOTO 12
-       END IF
-       IF(I.EQ.1)JC=JN
-      30   CONTINUE
+        JC=0
+        DO 30 I=NL-1,1,-1
+         JN=0
+          SUM=TH(I)*(1.+Q(I)*EPSI-Q(I))
+         DO 10 J=I+1,NL
+          SUM=SUM+TH(J)*(1.+Q(J)*EPSI-Q(J))
+          THBAR=SUM/FLOAT(J+1-I)
+          IF((TH(J)*(1.+Q(J)*EPSI-Q(J))).LT.THBAR)JN=J
+   10    CONTINUE
+         IF(I.EQ.1)JN=MAX(JN,2)
+         IF(JN.EQ.0)GOTO 30
+   12    CONTINUE
+         AHM=0.0
+         RM=0.0
+         UM=0.0
+         VM=0.0
+         DO K=1,NTRA
+          TRATM(K)=0.0
+         END DO
+         DO 15 J=I,JN
+          AHM=AHM+(CPD*(1.-Q(J))+Q(J)*CPV)*T(J)*(PH(J)-PH(J+1))
+          RM=RM+Q(J)*(PH(J)-PH(J+1))
+          UM=UM+U(J)*(PH(J)-PH(J+1))
+          VM=VM+V(J)*(PH(J)-PH(J+1))
+          DO K=1,NTRA
+           TRATM(K)=TRATM(K)+TRA(J,K)*(PH(J)-PH(J+1))
+          END DO
+   15    CONTINUE
+         DPHINV=1./(PH(I)-PH(JN+1))
+         RM=RM*DPHINV
+         UM=UM*DPHINV
+         VM=VM*DPHINV
+         DO K=1,NTRA
+          TRATM(K)=TRATM(K)*DPHINV
+         END DO
+         A2=0.0
+         DO 20 J=I,JN
+          Q(J)=RM
+          U(J)=UM
+          V(J)=VM
+          DO K=1,NTRA
+           TRA(J,K)=TRATM(K)
+          END DO
+          RDCP=(RD*(1.-Q(J))+Q(J)*RV)/
+     1     (CPD*(1.-Q(J))+Q(J)*CPV)
+          X=(0.001*P(J))**RDCP
+          TOLD(J)=T(J)
+          T(J)=X
+          A2=A2+(CPD*(1.-Q(J))+Q(J)*CPV)*X*(PH(J)-PH(J+1))
+   20    CONTINUE
+         DO 25 J=I,JN
+          TH(J)=AHM/A2
+          T(J)=T(J)*TH(J)
+          TC=TOLD(J)-273.15
+          ALV=LV0-CPVMCL*TC
+          QS(J)=QS(J)+QS(J)*(1.+QS(J)*(EPSI-1.))*ALV*(T(J)-
+     1     TOLD(J))/(RV*TOLD(J)*TOLD(J))
+   25    CONTINUE
+         IF((TH(JN+1)*(1.+Q(JN+1)*EPSI-Q(JN+1))).LT.
+     1    (TH(JN)*(1.+Q(JN)*EPSI-Q(JN))))THEN
+          JN=JN+1
+          GOTO 12
+         END IF
+         IF(I.EQ.1)JC=JN
+   30   CONTINUE
 C
 C   ***   Remove any supersaturation that results from adjustment ***
 C
       IF(JC.GT.1)THEN
-      DO 38 J=1,JC
-        IF(QS(J).LT.Q(J))THEN
-         ALV=LV0-CPVMCL*(T(J)-273.15)
-         TNEW=T(J)+ALV*(Q(J)-QS(J))/(CPD*(1.-Q(J))+
-      1      CL*Q(J)+QS(J)*(CPV-CL+ALV*ALV/(RV*T(J)*T(J))))
-         ALVNEW=LV0-CPVMCL*(TNEW-273.15)
-         QNEW=(ALV*Q(J)-(TNEW-T(J))*(CPD*(1.-Q(J))+CL*Q(J)))/ALVNEW
-         PRECIP=PRECIP+24.*3600.*1.0E5*(PH(J)-PH(J+1))*
-      1      (Q(J)-QNEW)/(G*DELT*ROWL)
-         T(J)=TNEW
-         Q(J)=QNEW
-         QS(J)=QNEW
-        END IF
-      38  CONTINUE
+       DO 38 J=1,JC
+          IF(QS(J).LT.Q(J))THEN
+           ALV=LV0-CPVMCL*(T(J)-273.15)
+           TNEW=T(J)+ALV*(Q(J)-QS(J))/(CPD*(1.-Q(J))+
+     1      CL*Q(J)+QS(J)*(CPV-CL+ALV*ALV/(RV*T(J)*T(J))))
+           ALVNEW=LV0-CPVMCL*(TNEW-273.15)
+           QNEW=(ALV*Q(J)-(TNEW-T(J))*(CPD*(1.-Q(J))+CL*Q(J)))/ALVNEW
+           PRECIP=PRECIP+24.*3600.*1.0E5*(PH(J)-PH(J+1))*
+     1      (Q(J)-QNEW)/(G*DELT*ROWL)
+           T(J)=TNEW
+           Q(J)=QNEW
+           QS(J)=QNEW
+          END IF
+   38  CONTINUE
       END IF
+C
       RETURN
       END

climlab/convection/_emanuel_convection/setup.py

@@ -38,6 +38,7 @@ def configuration(parent_package='', top_path=None):
 def gen_source(ext, build_dir):
     thispath = config.local_path
     sourcelist = []
+    sourcelist.append(join(thispath,'dry_adjustment.f'))
     sourcelist.append(join(thispath,'convect.f'))
     sourcelist.append(join(thispath,'Driver.f90'))
     try:

climlab/convection/emanuel_convection.py

@@ -85,6 +85,7 @@ class EmanuelConvection(TimeDependentProcess):
         - DTMAX = 0.9, maximum negative temperature perturbation a lifted parcel is allowed to have below its LFC
         - ALPHA = 0.2, first parameter that controls the rate of approach to quasi-equilibrium
         - DAMP = 0.1, second parameter that controls the rate of approach to quasi-equilibrium (DAMP must be less than 1)
+        - IPBL = 0, switch to bypass the dry convective adjustment (bypass if IPBL==0)
 
     Tendencies computed:
 
@@ -162,6 +163,7 @@ def __init__(self,
             DTMAX=0.9,
             ALPHA=0.2,
             DAMP=0.1,
+            IPBL=0,
             **kwargs):
         super(EmanuelConvection, self).__init__(**kwargs)
         self.time_type = 'explicit'
@@ -186,6 +188,7 @@ def __init__(self,
         self.add_input('DTMAX', DTMAX)
         self.add_input('ALPHA', ALPHA)
         self.add_input('DAMP', DAMP)
+        self.add_input('IPBL', IPBL)
 
     def _compute(self):
         #  Invert arrays so the first element is the bottom of column
@@ -210,13 +213,17 @@ def _compute(self):
         TRA = np.zeros((NCOL,ND,NTRA), order='F')  # tracers ignored
         DELT = float(self.timestep)
         CBMF = self.CBMF
-        (IFLAG, FT, FQ, FU, FV, FTRA, PRECIP, WD, TPRIME, QPRIME, CBMFnew) = \
-            convect(T, Q, QS, U, V, TRA, P, PH, NCOL, ND, NL, NTRA, DELT, CBMF,
+        (IFLAG, FT, FQ, FU, FV, FTRA, PRECIP, WD, TPRIME, QPRIME, CBMFnew,
+        Tout, Qout, QSout, Uout, Vout, TRAout) = \
+            convect(T, Q, QS, U, V, TRA, P, PH, NCOL, ND, NL, NTRA, DELT, self.IPBL, CBMF,
                     CPD, CPV, CL, RV, RD, LV0, G, ROWL, self.MINORIG,
                     self.ELCRIT, self.TLCRIT, self.ENTP, self.SIGD, self.SIGS,
                     self.OMTRAIN, self.OMTSNOW, self.COEFFR, self.COEFFS,
                     self.CU, self.BETA, self.DTMAX, self.ALPHA, self.DAMP
                     )
+        # If dry adjustment is being used then the tendencies need to be adjusted
+        FT += (Tout - T) / DELT
+        FQ += (Qout - Q) / DELT
         tendencies = {'Tatm': _convect_to_climlab(FT)*np.ones_like(self.state['Tatm']),
                       'q': _convect_to_climlab(FQ)*np.ones_like(self.state['q'])}
         if 'Ts' in self.state:

conda-recipe/meta.yaml

@@ -1,5 +1,5 @@
 {% set name = "climlab" %}
-{% set version = "0.6.4" %}
+{% set version = "0.6.5.dev0" %}
 
 package:
   name: {{ name|lower }}
@@ -9,7 +9,7 @@ source:
   path: ../
 
 build:
-  number: 1
+  number: 0
 
 requirements:
   build:

setup.py

@@ -1,7 +1,7 @@
 import os, sys
 import textwrap
 
-VERSION = '0.6.4'
+VERSION = '0.6.5.dev0'
 
 # BEFORE importing setuptools, remove MANIFEST. Otherwise it may not be
 # properly updated when the contents of directories change (true for distutils,