pkg/atm_phys/dargan_bettsmiller_mod.F90

b2ea1d2979 Jean*0001 module dargan_bettsmiller_mod
                0002
                0003 ! modified by pog:
                0004 !
                0005 ! - changed comments to say r is mixing ratio in capecalcnew
                0006 ! - removed unused avg_bl code
                0007 ! - removed second order in p code that was commented out
                0008 ! - add fatal error if LCL table lookup out of range
                0009 ! - changed so mixing ratio has conventional definition
                0010 ! - changed lcl table routine call so uses conventional mixing ratio
                0011 ! - didn't change approximate saturated adiabatic ascent integration or
                0012 !   approximate 'saturate out at constant p if r0>rs' calculation to be
                0013 !   consistent with conventional definition of mixing ratio
                0014 ! - CAPE/CIN calculation uses virtual temperature if option set
                0015
                0016 !----------------------------------------------------------------------
                0017 !use            fms_mod, only:  file_exist, error_mesg, open_file,  &
                0018 !                               check_nml_error, mpp_pe, FATAL,  &
                0019 !                               close_file
                0020
                0021 use simple_sat_vapor_pres_mod, only:  escomp, descomp
                0022 use     gcm_params_mod, only: gcm_LEN_MBUF, gcm_SQZ_R, gcm_stdMsgUnit
                0023 use      constants_mod, only:  HLv,HLs,Cp_air,Grav,rdgas,rvgas, &
                0024                                kappa
                0025
                0026 implicit none
                0027 private
                0028 !---------------------------------------------------------------------
                0029 !  ---- public interfaces ----
                0030
                0031    public  dargan_bettsmiller, dargan_bettsmiller_init
                0032
                0033 !-----------------------------------------------------------------------
                0034 !   ---- version number ----
                0035
15388b052e Jean*0036  character(len=128) :: version = '$Id: dargan_bettsmiller_mod.F90,v 1.3 2017/08/11 20:48:50 jmc Exp $'
b2ea1d2979 Jean*0037  character(len=128) :: tag = '$Name:  $'
                0038
                0039 !-----------------------------------------------------------------------
                0040 !   ---- local/private data ----
                0041
                0042     logical :: do_init=.true.
                0043
                0044 !-----------------------------------------------------------------------
                0045 !   --- namelist ----
                0046
                0047 real    :: tau_bm=7200.
                0048 real    :: rhbm = .8
                0049 logical :: do_virtual = .false.
                0050 logical :: do_shallower = .false.
                0051 logical :: do_changeqref = .false.
                0052 logical :: do_envsat = .false.
                0053 logical :: do_taucape = .false.
                0054 logical :: do_bm_shift = .false.
                0055 real    :: capetaubm = 900.
                0056 real    :: tau_min = 2400.
                0057
                0058 namelist /dargan_bettsmiller_nml/  tau_bm, rhbm, &
                0059                             do_shallower, do_changeqref, &
                0060                             do_envsat, do_taucape, capetaubm, tau_min, &
                0061                             do_virtual, do_bm_shift
                0062
                0063 !-----------------------------------------------------------------------
                0064 !           description of namelist variables
                0065 !
                0066 !  tau_bm    =  betts-miller relaxation timescale (seconds)
                0067 !
                0068 !  rhbm      = relative humidity that you're relaxing towards
                0069 !
                0070 !  do_shallower = do the shallow convection scheme where it chooses a smaller
                0071 !                 depth such that precipitation is zero
                0072 !
                0073 !  do_changeqref = do the shallow convection scheme where if changes the
                0074 !                  profile of both q and T in order make precip zero
                0075 !
                0076 !  do_envsat = reference profile is rhbm times saturated wrt environment
                0077 !              (if false, it's rhbm times parcel)
                0078 !
                0079 !  do_taucape = scheme where taubm is proportional to CAPE**-1/2
                0080 !
                0081 !  capetaubm = for the above scheme, the value of CAPE for which
                0082 !              tau = tau_bm
                0083 !
                0084 !  tau_min   = minimum relaxation time allowed for the above scheme
                0085 !
                0086 !  do_virtual = use virtual temperature for CAPE/CIN calculation
                0087 !
                0088 !  do_bm_shift = always conserve enthalpy by shifting temperature
                0089 !-----------------------------------------------------------------------
                0090
                0091 contains
                0092
                0093 !#######################################################################
                0094
                0095    subroutine dargan_bettsmiller (dt, tin, qin, pfull, phalf, coldT, &
                0096                                   rain, snow, tdel, qdel, q_ref, bmflag, &
                0097                                   klzbs, cape, cin, t_ref,invtau_bm_t,invtau_bm_q, &
                0098                                   capeflag, bi,bj,myIter, myThid, mask, conv)
                0099
                0100 !-----------------------------------------------------------------------
                0101 !
                0102 !                     Betts-Miller Convection Scheme
                0103 !
                0104 !-----------------------------------------------------------------------
                0105 !
                0106 !   input:  dt       time step in seconds
                0107 !           tin      temperature at full model levels
                0108 !           qin      specific humidity of water vapor at full
                0109 !                      model levels
                0110 !           pfull    pressure at full model levels
                0111 !           phalf    pressure at half (interface) model levels
                0112 !           coldT    should precipitation be snow at this point?
                0113 !   optional:
                0114 !           mask     optional mask (0 or 1.)
                0115 !           conv     logical flag; if true then no betts-miller
                0116 !                       adjustment is performed at that grid-point or
                0117 !                       model level
                0118 !
                0119 !  output:  rain     liquid precipitation (kg/m2)
                0120 !           snow     frozen precipitation (kg/m2)
                0121 !           tdel     temperature tendency at full model levels
                0122 !           qdel     specific humidity tendency (of water vapor) at
                0123 !                      full model levels
                0124 !           bmflag   flag for which routines you're calling
599643daec Jean*0125 !           klzbs    stored (integer part) klzb and ktop (decimal part) values
b2ea1d2979 Jean*0126 !           cape     convectively available potential energy
                0127 !           cin      convective inhibition (this and the above are before the
                0128 !                    adjustment)
                0129 !           invtau_bm_t temperature relaxation timescale
                0130 !           invtau_bm_q humidity relaxation timescale
                0131 !           capeflag a flag that says why cape=0
                0132
                0133 !-----------------------------------------------------------------------
                0134 !--------------------- interface arguments -----------------------------
                0135
                0136    real   , intent(in) , dimension(:,:,:) :: tin, qin, pfull, phalf
                0137    real   , intent(in)                    :: dt
                0138    logical   , intent(in) , dimension(:,:):: coldT
                0139    real   , intent(out), dimension(:,:)   :: rain,snow, bmflag, klzbs, cape, &
                0140        cin, invtau_bm_t, invtau_bm_q, capeflag
                0141    real   , intent(out), dimension(:,:,:) :: tdel, qdel, q_ref, t_ref
                0142    integer, intent(in)                    :: bi,bj, myIter
                0143    integer, intent(in)                    :: myThid
                0144    real   , intent(in) , dimension(:,:,:), optional :: mask
                0145    logical, intent(in) , dimension(:,:,:), optional :: conv
                0146 !-----------------------------------------------------------------------
                0147 !---------------------- local data -------------------------------------
                0148
599643daec Jean*0149 !logical,dimension(size(tin,1),size(tin,2),size(tin,3)) :: do_adjust
                0150    real,dimension(size(tin,1),size(tin,2),size(tin,3)) :: rin
                0151    real,dimension(size(tin,1),size(tin,2))             :: precip, precip_t
                0152    real,dimension(size(tin,3))                         :: eref, rpc, tpc
b2ea1d2979 Jean*0153
                0154    real                                                ::  &
                0155        cape1, cin1, tot, deltak, deltaq, qrefint, deltaqfrac, deltaqfrac2, &
599643daec Jean*0156        ptopfrac, es, capeflag1, small
                0157  integer  i, j, k, ix, jx, kx, klzb, ktop
b2ea1d2979 Jean*0158 !-----------------------------------------------------------------------
                0159 !     computation of precipitation by betts-miller scheme
                0160 !-----------------------------------------------------------------------
                0161       capeflag1 = 0.
                0162
                0163 !     if (do_init) call error_mesg ('dargan_bettsmiller',  &
                0164 !                        'dargan_bettsmiller_init has not been called.', FATAL)
                0165
                0166       ix=size(tin,1)
                0167       jx=size(tin,2)
                0168       kx=size(tin,3)
                0169       small = 1.e-10
                0170
599643daec Jean*0171 ! initialise output:
                0172       precip = 0.
                0173       tdel   = 0.
                0174       qdel   = 0.
                0175       t_ref  = tin
                0176       q_ref  = qin
                0177       bmflag = 0.
                0178       klzbs  = 0.
                0179       cape   = 0.
                0180       cin    = 0.
                0181       invtau_bm_t = 0.
                0182       invtau_bm_q = 0.
                0183       precip_t = 0.
                0184
b2ea1d2979 Jean*0185 ! calculate r (where r is the mixing ratio)
                0186        rin = qin/(1.0 - qin)
                0187
599643daec Jean*0188        do j=1,jx
                0189           do i=1,ix
b2ea1d2979 Jean*0190              cape1 = 0.
                0191              cin1 = 0.
                0192              tot = 0.
                0193              klzb=0
                0194 ! the bmflag is written out to show what aspects of the bm scheme is called
                0195 ! bmflag = 0 is no cape, no convection
                0196 ! bmflag = 1 is shallow conv, the predicted precip is less than zero
                0197 ! bmflag = 2 is deep convection
                0198              bmflag(i,j) = 0.
                0199              tpc = tin(i,j,:)
                0200              rpc = rin(i,j,:)
                0201 ! calculate cape, cin, level of zero buoyancy, and parcel properties
                0202 ! new code (second order in delta ln p and exact LCL calculation)
                0203              call capecalcnew( kx,  pfull(i,j,:),  phalf(i,j,:),&
                0204                             cp_air, rdgas, rvgas, hlv, kappa, tin(i,j,:), &
                0205                             rin(i,j,:), cape1, cin1, tpc, &
                0206                             rpc, klzb, i,j,bi,bj,myIter,myThid)
                0207
                0208 ! set values for storage
                0209              capeflag(i,j) = capeflag1
                0210              cape(i,j) = cape1
                0211              cin(i,j) = cin1
                0212              klzbs(i,j) = klzb
                0213              if(cape1.gt.0.) then
                0214 !             if((tot.gt.0.).and.(cape1.gt.0.)) then
                0215                 bmflag(i,j) = 1.
                0216 ! reference temperature is just that of the parcel all the way up
599643daec Jean*0217 !               t_ref(i,j,:) = tpc
                0218                 t_ref(i,j,klzb:kx) = tpc(klzb:kx)
b2ea1d2979 Jean*0219                 do k=klzb,kx
                0220 ! sets reference spec hum to a certain relative hum (change to vapor pressure,
                0221 ! multiply by rhbm, then back to spec humid)
                0222                    if(do_envsat) then
                0223                       call escomp(tin(i,j,k),es)
                0224                       es = es*rhbm
                0225                       rpc(k) = mixing_ratio(es, pfull(i,j,k))
                0226                       q_ref(i,j,k) = rpc(k)/(1 + rpc(k))
                0227                    else
                0228                       eref(k) = rhbm*pfull(i,j,k)*rpc(k)/(rdgas/rvgas + rpc(k))
                0229                       rpc(k) = mixing_ratio(eref(k),pfull(i,j,k))
                0230                       q_ref(i,j,k) = rpc(k)/(1 + rpc(k))
                0231                    endif
                0232                 end do
                0233 ! set the tendencies to zero where you don't adjust
                0234 ! set the reference profiles to be the original profiles (for diagnostic
                0235 ! purposes only --  you can think of this as what you're relaxing to in
                0236 ! areas above the actual convection
599643daec Jean*0237 !               do k=1,max(klzb-1,1)
                0238 !                  qdel(i,j,k) = 0.0
                0239 !                  tdel(i,j,k) = 0.0
                0240 !                  q_ref(i,j,k) = qin(i,j,k)
                0241 !                  t_ref(i,j,k) = tin(i,j,k)
                0242 !               end do
b2ea1d2979 Jean*0243 ! initialize p to zero for the loop
                0244                 precip(i,j) = 0.
                0245                 precip_t(i,j) = 0.
                0246 ! makes t_bm prop to (CAPE)**-.5.  Gives a relaxation time of tau_bm when
                0247 ! CAPE = sqrt(capetaubm)
                0248                 if(do_taucape) then
                0249                    tau_bm = sqrt(capetaubm)*tau_bm/sqrt(cape1)
                0250                    if(tau_bm.lt.tau_min) tau_bm = tau_min
                0251                 endif
                0252                 do k=klzb, kx
                0253 ! relax to reference profiles
                0254                    tdel(i,j,k) = - (tin(i,j,k) - t_ref(i,j,k))/tau_bm*dt
                0255                    qdel(i,j,k) = - (qin(i,j,k) - q_ref(i,j,k))/tau_bm*dt
                0256 ! Precipitation can be calculated already, based on the change in q on the
                0257 ! way up (this isn't altered in the energy conservation scheme).
599643daec Jean*0258                    precip(i,j)  = precip(i,j) - qdel(i,j,k)           &
                0259                                 *(phalf(i,j,k+1)- phalf(i,j,k))/grav
                0260                    precip_t(i,j)= precip_t(i,j) + cp_air/(hlv+small)*tdel(i,j,k) &
                0261                                 *(phalf(i,j,k+1)-phalf(i,j,k))/grav
b2ea1d2979 Jean*0262                 end do
                0263                 if ((precip(i,j).gt.0.).and.(precip_t(i,j).gt.0.)) then
                0264 ! If precip > 0, then correct energy.
                0265                    bmflag(i,j) = 2.
                0266 ! not simple scheme: shift the reference profile of temperature
                0267 ! deltak is the energy correction that you make to the temperature reference
                0268 ! profile
                0269                    if(precip(i,j).gt.precip_t(i,j) .and. (.not. do_bm_shift)) then
                0270 ! if the q precip is greater, then lengthen the relaxation timescale on q to
                0271 ! conserve energy.  qdel is therefore changed.
                0272                       invtau_bm_q(i,j) = precip_t(i,j)/precip(i,j)/tau_bm
                0273                       qdel(i,j,klzb:kx) = tau_bm*invtau_bm_q(i,j)* &
                0274                          qdel(i,j,klzb:kx)
                0275                       precip(i,j) = precip_t(i,j)
                0276                       invtau_bm_t(i,j) = 1./tau_bm
                0277                    else
                0278
                0279                          deltak = 0.
                0280                          do k=klzb, kx
                0281 ! Calculate the integrated difference in energy change within each level.
                0282                             deltak = deltak - (tdel(i,j,k) + hlv/cp_air*&
                0283                                      qdel(i,j,k))* &
                0284                                      (phalf(i,j,k+1) - phalf(i,j,k))
                0285                          end do
                0286 ! Divide by total pressure.
                0287                          deltak = deltak/(phalf(i,j,kx+1) - phalf(i,j,klzb))
                0288 ! Adjust the reference profile (uniformly with height), and correspondingly
                0289 ! the temperature change.
                0290                          t_ref(i,j,klzb:kx) = t_ref(i,j,klzb:kx)+ &
                0291                               deltak*tau_bm/dt
                0292                          tdel(i,j,klzb:kx) = tdel(i,j,klzb:kx) + deltak
                0293                    endif
                0294                 else if(precip_t(i,j).gt.0.) then
                0295 ! If precip < 0, then do the shallow conv routine.
                0296 ! First option: do_shallower = true
                0297 ! This chooses the depth of convection based on choosing the height that
                0298 ! it can make precip zero, i.e., subtract off heights until that precip
                0299 ! becomes positive.
                0300                    if (do_shallower) then
                0301 ! ktop is the new top of convection.  set this initially to klzb.
                0302                       ktop = klzb
                0303 ! Work your way down until precip is positive again.
                0304                       do while ((precip(i,j).lt.0.).and.(ktop.le.kx))
                0305                          precip(i,j) = precip(i,j) - qdel(i,j,ktop)* &
                0306                                   (phalf(i,j,ktop) - phalf(i,j,ktop+1))/grav
                0307                          ktop = ktop + 1
                0308                       end do
                0309 ! since there will be an overshoot (precip is going to be greater than zero
                0310 ! once we finish this), the actual new top of convection is somewhere between
                0311 ! the current ktop, and one level above this.  set ktop to the level above.
                0312                       ktop = ktop - 1
599643daec Jean*0313                       klzbs(i,j) = klzbs(i,j) + float(ktop)/( kx + 1.d0 )
b2ea1d2979 Jean*0314 ! Adjust the tendencies in the places above back to zero, and the reference
                0315 ! profiles back to the original t,q.
                0316                       if (ktop.gt.klzb) then
                0317                          qdel(i,j,klzb:ktop-1) = 0.
599643daec Jean*0318 !                        q_ref(i,j,klzb:ktop-1) = qin(i,j,klzb:ktop-1)
b2ea1d2979 Jean*0319                          tdel(i,j,klzb:ktop-1) = 0.
599643daec Jean*0320 !                        t_ref(i,j,klzb:ktop-1) = tin(i,j,klzb:ktop-1)
b2ea1d2979 Jean*0321                       end if
                0322 ! Then make the change only a fraction of the new top layer so the precip is
                0323 ! identically zero.
                0324 ! Calculate the fractional penetration of convection through that top layer.
                0325 ! This is the amount necessary to make precip identically zero.
                0326                       if (precip(i,j).gt.0.) then
                0327                          ptopfrac = precip(i,j)/(qdel(i,j,ktop)* &
                0328                             (phalf(i,j,ktop+1) - phalf(i,j,ktop)))*grav
                0329 ! Reduce qdel in the top layer by this fraction.
                0330                          qdel(i,j,ktop) = ptopfrac*qdel(i,j,ktop)
                0331 ! Set precip to zero
                0332                          precip(i,j) = 0.
                0333 ! Now change the reference temperature in such a way to make the net
                0334 ! heating zero.
                0335
                0336 !! modification: pog
                0337 !! Reduce tdel in the top layer
                0338                          tdel(i,j,ktop) = ptopfrac*tdel(i,j,ktop)
                0339 !! end modification: pog
                0340
                0341                          deltak = 0.
                0342                          if (ktop.lt.kx) then
                0343 !! modification: pog
                0344 !! include the full mass of the top layer when calculating delta_k
                0345 !! also use the entire mass of the column when normalizing
                0346 !
                0347 ! Integrate temperature tendency
                0348                             do k=ktop,kx
                0349                                deltak = deltak + tdel(i,j,k)* &
                0350                                    (phalf(i,j,k) - phalf(i,j,k+1))
                0351                             end do
                0352
                0353 ! Normalize by the pressure difference.
                0354                             deltak = deltak/(phalf(i,j,kx+1) - phalf(i,j,ktop))
                0355
                0356 !! end modification: pog
                0357
                0358 ! Subtract this value uniformly from tdel, and make the according change to
                0359 ! t_ref.
                0360                             do k=ktop,kx
                0361                                tdel(i,j,k) = tdel(i,j,k) + deltak
                0362                                t_ref(i,j,k) = t_ref(i,j,k) + deltak*tau_bm/dt
                0363                             end do
                0364                          end if
                0365                       else
                0366                          precip(i,j) = 0.
                0367                          qdel(i,j,kx) = 0.
599643daec Jean*0368 !                        q_ref(i,j,kx) = qin(i,j,kx)
b2ea1d2979 Jean*0369                          tdel(i,j,kx) = 0.
599643daec Jean*0370 !                        t_ref(i,j,kx) = tin(i,j,kx)
                0371 !                        invtau_bm_t(i,j) = 0.
                0372 !                        invtau_bm_q(i,j) = 0.
b2ea1d2979 Jean*0373                       end if
                0374                    else if(do_changeqref) then
                0375 ! Change the reference profile of q by a certain fraction so that precip is
                0376 ! zero.  This involves calculating the total integrated q_ref dp (this is the
                0377 ! quantity intqref), as well as the necessary change in q_ref (this is the
                0378 ! quantity deltaq).  Then the fractional change in q_ref at each level (the
                0379 ! quantity deltaqfrac) is 1-deltaq/intqref.  (have to multiply q_ref by
                0380 ! 1-deltaq/intqref at every level)  Then the change in qdel is
                0381 ! -deltaq/intqref*q_ref*dt/tau_bm.
                0382 ! Change the reference profile of T by a uniform amount so that precip is zero.
                0383                       deltak = 0.
                0384                       deltaq = 0.
                0385                       qrefint = 0.
                0386                       do k=klzb,kx
                0387 ! deltaq = a positive quantity (since int qdel is positive).  It's how
                0388 ! much q_ref must be changed by, in an integrated sense.  The requisite
                0389 ! change in qdel is this without the factors of tau_bm and dt.
                0390                          deltaq = deltaq - qdel(i,j,k)*tau_bm/dt* &
                0391                                    (phalf(i,j,k) - phalf(i,j,k+1))
                0392 ! deltak = the amount tdel needs to be changed
                0393                          deltak  = deltak  + tdel(i,j,k)* &
                0394                                    (phalf(i,j,k) - phalf(i,j,k+1))
                0395 ! qrefint = integrated value of qref
                0396                          qrefint = qrefint - q_ref(i,j,k)* &
                0397                                    (phalf(i,j,k) - phalf(i,j,k+1))
                0398                       end do
                0399 ! Normalize deltak by total pressure.
                0400                       deltak  = deltak /(phalf(i,j,kx+1) - phalf(i,j,klzb))
                0401 ! multiplying factor for q_ref is 1 + the ratio
                0402                       deltaqfrac = 1. - deltaq/qrefint
                0403 ! multiplying factor for qdel adds dt/tau_bm
                0404                       deltaqfrac2 = - deltaq/qrefint*dt/tau_bm
                0405                       precip(i,j) = 0.0
                0406                       do k=klzb,kx
                0407                          qdel(i,j,k) = qdel(i,j,k) + deltaqfrac2*q_ref(i,j,k)
                0408                          q_ref(i,j,k) = deltaqfrac*q_ref(i,j,k)
                0409                          tdel(i,j,k) = tdel(i,j,k) + deltak
                0410                          t_ref(i,j,k) = t_ref(i,j,k) + deltak*tau_bm/dt
                0411                       end do
                0412                    else
                0413                       precip(i,j) = 0.
                0414                       tdel(i,j,:) = 0.
                0415                       qdel(i,j,:) = 0.
599643daec Jean*0416 !                     invtau_bm_t(i,j) = 0.
                0417 !                     invtau_bm_q(i,j) = 0.
b2ea1d2979 Jean*0418                    end if
                0419 ! for cases where virtual temp predicts CAPE but precip_t < 0.
                0420 ! - also note cape and precip_t are different because cape
                0421 ! involves integration wrt dlogp and precip_t wrt dp
                0422                 else
                0423                    tdel(i,j,:) = 0.0
                0424                    qdel(i,j,:) = 0.0
                0425                    precip(i,j) = 0.0
599643daec Jean*0426 !                  q_ref(i,j,:) = qin(i,j,:)
                0427 !                  t_ref(i,j,:) = tin(i,j,:)
                0428 !                  invtau_bm_t(i,j) = 0.
                0429 !                  invtau_bm_q(i,j) = 0.
b2ea1d2979 Jean*0430                 end if
                0431 ! if no CAPE, set tendencies to zero.
599643daec Jean*0432 !            else
                0433 !               tdel(i,j,:) = 0.0
                0434 !               qdel(i,j,:) = 0.0
                0435 !               precip(i,j) = 0.0
                0436 !               q_ref(i,j,:) = qin(i,j,:)
                0437 !               t_ref(i,j,:) = tin(i,j,:)
                0438 !               invtau_bm_t(i,j) = 0.
                0439 !               invtau_bm_q(i,j) = 0.
b2ea1d2979 Jean*0440              end if
                0441           end do
                0442        end do
                0443
                0444        rain = precip
                0445        snow = 0.
599643daec Jean*0446 !      snow = precip_t  !to output value of precip_t (debug)
b2ea1d2979 Jean*0447
                0448    end subroutine dargan_bettsmiller
                0449
                0450 !#######################################################################
                0451
                0452 !all new cape calculation.
                0453
                0454       subroutine capecalcnew(kx,p,phalf,cp_air,rdgas,rvgas,hlv,kappa,tin,rin,&
                0455                              cape,cin,tp,rp,klzb, i,j,bi,bj,myIter,myThid)
                0456
                0457 !
                0458 !    Input:
                0459 !
                0460 !    kx          number of levels
                0461 !    p           pressure (index 1 refers to TOA, index kx refers to surface)
                0462 !    phalf       pressure at half levels
                0463 !    cp_air      specific heat of dry air
                0464 !    rdgas       gas constant for dry air
                0465 !    rvgas       gas constant for water vapor
                0466 !    hlv         latent heat of vaporization
                0467 !    kappa       the constant kappa
                0468 !    tin         temperature of the environment
                0469 !    rin         mixing ratio of the environment
                0470 !
                0471 !    Output:
                0472 !    cape        Convective available potential energy
                0473 !    cin         Convective inhibition (if there's no LFC, then this is set
                0474 !                to zero)
                0475 !    tp          Parcel temperature (set to the environmental temperature
                0476 !                where no adjustment)
                0477 !    rp          Parcel mixing ratio (set to the environmental humidity
                0478 !                where no adjustment, and set to the saturation humidity at
                0479 !                the parcel temperature below the LCL)
                0480 !    klzb        Level of zero buoyancy
                0481 !
                0482 !    Algorithm:
                0483 !    Start with surface parcel.
                0484 !    Calculate the lifting condensation level (uses an analytic formula and a
                0485 !       lookup table).
                0486 !    Calculate parcel ascent up to LZB.
                0487 !    Calculate CAPE and CIN.
                0488       implicit none
                0489       integer, intent(in)                    :: kx
                0490       real, intent(in), dimension(:)         :: p, phalf, tin, rin
                0491       real, intent(in)                       :: rdgas, rvgas, hlv, kappa, cp_air
                0492       integer, intent(out)                   :: klzb
                0493       real, intent(out), dimension(:)        :: tp, rp
                0494       real, intent(out)                      :: cape, cin
                0495       integer, intent(in)                    :: i,j,bi,bj,myIter
                0496       integer, intent(in)                    :: myThid
599643daec Jean*0497       integer            :: k, klcl, klfc
b2ea1d2979 Jean*0498       logical            :: nocape
599643daec Jean*0499       real, dimension(kx)   :: tin_virtual
b2ea1d2979 Jean*0500       real                  :: t0, r0, es, rs, theta0, pstar, value, tlcl, &
599643daec Jean*0501                                a, b, dtdlnp, &
                0502                                plcl, plzb, small
b2ea1d2979 Jean*0503
                0504       pstar = 1.e5
                0505 ! so we can run dry limit (one expression involves 1/hlv)
                0506       small = 1.e-10
                0507
                0508       nocape = .true.
                0509       cape = 0.
                0510       cin = 0.
                0511       plcl = 0.
                0512       plzb = 0.
                0513       klfc = 0
                0514       klcl = 0
                0515       klzb = 0
                0516       tp(1:kx) = tin(1:kx)
                0517       rp(1:kx) = rin(1:kx)
                0518
                0519       ! calculate the virtual temperature
                0520       do k=1,kx
                0521        tin_virtual(k) = virtual_temp(tin(k), rin(k))
                0522       enddo
                0523
                0524 ! start with surface parcel
                0525       t0 = tin(kx)
                0526       r0 = rin(kx)
                0527 ! calculate the lifting condensation level by the following:
                0528 ! are you saturated to begin with?
                0529       call escomp(t0,es)
                0530       rs = mixing_ratio(es, p(kx))
                0531       if (r0.ge.rs) then
                0532 ! if you�re already saturated, set lcl to be the surface value.
                0533          plcl = p(kx)
                0534 ! the first level where you�re completely saturated.
                0535          klcl = kx
                0536 ! saturate out to get the parcel temp and humidity at this level
                0537 ! first order (in delta T) accurate expression for change in temp
                0538 ! note this is assumed to happen at constant pressure, also the resulting
                0539 ! saturation mixing ratio is different from rs and this is accounted for in
                0540 ! the formula, but the formula is approximate and should be replaced
                0541          tp(kx) = t0 + (r0 - rs)/(cp_air/(hlv+small) + hlv*rs/rvgas/t0**2.)
                0542          call escomp(tp(kx),es)
                0543          rp(kx) = mixing_ratio(es, p(kx))
                0544       else
                0545 ! if not saturated to begin with, use the analytic expression to calculate the
                0546 ! exact pressure and temperature where you�re saturated.
                0547          theta0 = tin(kx)*(pstar/p(kx))**kappa
                0548 ! the expression that we utilize is log(r/(Rd/Rv+r)/theta**(1/kappa)*pstar) =
                0549 ! log(es/T**(1/kappa))
                0550 ! The right hand side of this is only a function of temperature, therefore
                0551 ! this is put into a lookup table to solve for temperature.
                0552          if (r0.gt.0.) then
                0553             value = log(theta0**(-1/kappa)*pstar*r0/(rdgas/rvgas+r0))
                0554 !           call lcltabl(value,tlcl,myThid)
                0555             call lcltabl(value,tlcl,i,j,bi,bj,myIter,myThid)
                0556             plcl = pstar*(tlcl/theta0)**(1/kappa)
                0557 ! just in case plcl is very high up
                0558             if (plcl.lt.p(1)) then
                0559                plcl = p(1)
                0560                tlcl = theta0*(plcl/pstar)**kappa
                0561                write (*,*) 'hi lcl'
                0562             end if
                0563             k = kx
                0564          else
                0565 ! if the parcel mixing ratio is zero or negative, set lcl to top level
                0566             plcl = p(1)
                0567             tlcl = theta0*(plcl/pstar)**kappa
                0568 !            write (*,*) 'zero r0', r0
                0569             go to 11
                0570          end if
                0571 ! calculate the parcel temperature (adiabatic ascent) below the LCL.
                0572 ! the mixing ratio stays the same
                0573          do while (p(k).gt.plcl)
                0574             tp(k) = theta0*(p(k)/pstar)**kappa
                0575             call escomp(tp(k),es)
                0576 ! note rp is not the actual parcel mixing ratio here, but rather will be used
                0577 ! to calculate the reference moisture profile using rhbm
                0578             rp(k) = mixing_ratio(es,p(k))
                0579 ! this definition of CIN contains everything below the LCL
                0580 ! we use the actual parcel mixing ratio for the virtual temperature effect
                0581             cin = cin + rdgas*(tin_virtual(k)-virtual_temp(tp(k),r0))*log(phalf(k+1)/phalf(k))
                0582             k = k-1
                0583          end do
                0584 ! first level where you�re saturated at the level
                0585          klcl = k
                0586 ! do a saturated ascent to get the parcel temp at the LCL.
                0587 ! use your 2nd order equation up to the pressure above.
                0588 ! moist adaibat derivatives: (use the lcl values for temp, humid, and
                0589 ! pressure)
                0590 ! note moist adiabat derivatives are approximate and should be replaced
                0591 ! with more accurate expressions (see e.g. Holton pg 503 for derivation)
                0592          a = kappa*tlcl + hlv/cp_air*r0
                0593          b = hlv**2.*r0/cp_air/rvgas/tlcl**2.
                0594          dtdlnp = a/(1. + b)
                0595 ! first order in p
                0596 !         tp(klcl) = tlcl + dtdlnp*log(p(klcl)/plcl)
                0597 ! second order in p (RK2)
                0598 ! first get temp halfway up
                0599          tp(klcl) = tlcl + dtdlnp*log(p(klcl)/plcl)/2.
                0600          if ((tp(klcl).lt.173.16).and.nocape) go to 11
                0601          call escomp(tp(klcl),es)
                0602          rp(klcl) = mixing_ratio(es,(p(klcl) + plcl)/2)
                0603          a = kappa*tp(klcl) + hlv/cp_air*rp(klcl)
                0604          b = hlv**2./cp_air/rvgas*rp(klcl)/tp(klcl)**2.
                0605          dtdlnp = a/(1. + b)
                0606 ! second half of RK2
                0607          tp(klcl) = tlcl + dtdlnp*log(p(klcl)/plcl)
                0608          if ((tp(klcl).lt.173.16).and.nocape) go to 11
                0609          call escomp(tp(klcl),es)
                0610          rp(klcl) = mixing_ratio(es,p(klcl))
                0611 !         write (*,*) 'tp, rp klcl:kx, new', tp(klcl:kx), rp(klcl:kx)
                0612 ! CAPE/CIN stuff
                0613          if ((virtual_temp(tp(klcl),rp(klcl)).lt.tin_virtual(klcl)).and.nocape) then
                0614 ! if you�re not yet buoyant, then add to the CIN and continue
                0615             cin = cin + rdgas*(tin_virtual(klcl) - &
                0616                  virtual_temp(tp(klcl),rp(klcl)))*log(phalf(klcl+1)/phalf(klcl))
                0617          else
                0618 ! if you�re buoyant, then add to cape
                0619             cape = cape + rdgas*(virtual_temp(tp(klcl),rp(klcl)) - &
                0620                   tin_virtual(klcl))*log(phalf(klcl+1)/phalf(klcl))
                0621 ! if it�s the first time buoyant, then set the level of free convection to k
                0622             if (nocape) then
                0623                nocape = .false.
                0624                klfc = klcl
                0625             endif
                0626          end if
                0627       end if
                0628 ! then average the properties over the boundary layer if so desired.  to give
                0629 ! a new "parcel".  this may not be saturated at the LCL, so make sure you get
                0630 ! to a level where it is before moist adiabatic ascent!
                0631 !
                0632 ! then, start at the LCL, and do moist adiabatic ascent by the first order
                0633 ! scheme -- 2nd order as well
                0634       do k=klcl-1,1,-1
                0635 ! note moist adiabat derivatives are approximate and should be replaced
                0636 ! with more accurate expressions (see e.g. Holton pg 503 for derivation)
                0637          a = kappa*tp(k+1) + hlv/cp_air*rp(k+1)
                0638          b = hlv**2./cp_air/rvgas*rp(k+1)/tp(k+1)**2.
                0639          dtdlnp = a/(1. + b)
                0640 ! first order in p
                0641 !         tp(k) = tp(k+1) + dtdlnp*log(p(k)/p(k+1))
                0642 ! second order in p (RK2)
                0643 ! first get temp halfway up
                0644          tp(k) = tp(k+1) + dtdlnp*log(p(k)/p(k+1))/2.
                0645          if ((tp(k).lt.173.16).and.nocape) go to 11
                0646          call escomp(tp(k),es)
                0647          rp(k) = mixing_ratio(es,(p(k) + p(k+1))/2)
                0648          a = kappa*tp(k) + hlv/cp_air*rp(k)
                0649          b = hlv**2./cp_air/rvgas*rp(k)/tp(k)**2.
                0650          dtdlnp = a/(1. + b)
                0651 ! second half of RK2
                0652          tp(k) = tp(k+1) + dtdlnp*log(p(k)/p(k+1))
                0653 ! if you're below the lookup table value, just presume that there's no way
                0654 ! you could have cape and call it quits
                0655          if ((tp(k).lt.173.16).and.nocape) go to 11
                0656          call escomp(tp(k),es)
                0657          rp(k) = mixing_ratio(es,p(k))
                0658          if ((virtual_temp(tp(k),rp(k)).lt.tin_virtual(k)).and.nocape) then
                0659 ! if you�re not yet buoyant, then add to the CIN and continue
                0660             cin = cin + rdgas*(tin_virtual(k)-virtual_temp(tp(k),rp(k)))*log(phalf(k+1)/phalf(k))
                0661          elseif((virtual_temp(tp(k),rp(k)).lt.tin_virtual(k)).and.(.not.nocape)) then
                0662 ! if you have CAPE, and it�s your first time being negatively buoyant,
                0663 ! then set the level of zero buoyancy to k+1, and stop the moist ascent
                0664             klzb = k+1
                0665             go to 11
                0666          else
                0667 ! if you�re buoyant, then add to cape
                0668             cape = cape + rdgas*(virtual_temp(tp(k),rp(k))-tin_virtual(k))*log(phalf(k+1)/phalf(k))
                0669 ! if it�s the first time buoyant, then set the level of free convection to k
                0670             if (nocape) then
                0671                nocape = .false.
                0672                klfc = k
                0673             endif
                0674          end if
                0675       end do
                0676  11   if(nocape) then
                0677 ! this is if you made it through without having a LZB
                0678 ! set LZB to be the top level.
                0679          plzb = p(1)
                0680          klzb = 0
                0681          klfc = 0
                0682          cin = 0.
                0683          tp(1:kx) = tin(1:kx)
                0684          rp(1:kx) = rin(1:kx)
                0685       end if
                0686 !      write (*,*) 'plcl, klcl, tlcl, r0 new', plcl, klcl, tlcl, r0
                0687 !      write (*,*) 'tp, rp new', tp, rp
                0688 !       write (*,*) 'tp, new', tp
                0689 !       write (*,*) 'tin new', tin
                0690 !       write (*,*) 'klcl, klfc, klzb new', klcl, klfc, klzb
                0691       end subroutine capecalcnew
                0692
                0693 ! lookup table for the analytic evaluation of LCL
                0694 !     subroutine lcltabl(value,tlcl,myThid)
                0695       subroutine lcltabl(value,tlcl,i,j,bi,bj,myIter,myThid)
                0696 !
                0697 ! Table of values used to compute the temperature of the lifting condensation
                0698 ! level.
                0699 !
                0700 ! the expression that we utilize is log(r/(Rd/Rv+r)/theta**(1/kappa)*pstar) =
                0701 ! log(es/T**(1/kappa))
                0702 !
                0703 ! Gives the values of the temperature for the following range:
                0704 !   starts with -23, is uniformly distributed up to -10.4.  There are a
                0705 ! total of 127 values, and the increment is .1.
                0706 !
                0707       implicit none
                0708       real, intent(in)     :: value
                0709       real, intent(out)    :: tlcl
                0710       integer, intent(in)  :: i,j,bi,bj,myIter
                0711       integer, intent(in)  :: myThid
                0712
                0713       integer              :: ival
                0714       real, dimension(127) :: lcltable
                0715       real                 :: v1, v2
                0716
                0717       integer              :: errorMessageUnit
                0718       DATA errorMessageUnit / 15 /
                0719
                0720       data lcltable/  1.7364512e+02,   1.7427449e+02,   1.7490874e+02, &
                0721       1.7554791e+02,   1.7619208e+02,   1.7684130e+02,   1.7749563e+02, &
                0722       1.7815514e+02,   1.7881989e+02,   1.7948995e+02,   1.8016539e+02, &
                0723       1.8084626e+02,   1.8153265e+02,   1.8222461e+02,   1.8292223e+02, &
                0724       1.8362557e+02,   1.8433471e+02,   1.8504972e+02,   1.8577068e+02, &
                0725       1.8649767e+02,   1.8723077e+02,   1.8797006e+02,   1.8871561e+02, &
                0726       1.8946752e+02,   1.9022587e+02,   1.9099074e+02,   1.9176222e+02, &
                0727       1.9254042e+02,   1.9332540e+02,   1.9411728e+02,   1.9491614e+02, &
                0728       1.9572209e+02,   1.9653521e+02,   1.9735562e+02,   1.9818341e+02, &
                0729       1.9901870e+02,   1.9986158e+02,   2.0071216e+02,   2.0157057e+02, &
                0730       2.0243690e+02,   2.0331128e+02,   2.0419383e+02,   2.0508466e+02, &
                0731       2.0598391e+02,   2.0689168e+02,   2.0780812e+02,   2.0873335e+02, &
                0732       2.0966751e+02,   2.1061074e+02,   2.1156316e+02,   2.1252493e+02, &
                0733       2.1349619e+02,   2.1447709e+02,   2.1546778e+02,   2.1646842e+02, &
                0734       2.1747916e+02,   2.1850016e+02,   2.1953160e+02,   2.2057364e+02, &
                0735       2.2162645e+02,   2.2269022e+02,   2.2376511e+02,   2.2485133e+02, &
                0736       2.2594905e+02,   2.2705847e+02,   2.2817979e+02,   2.2931322e+02, &
                0737       2.3045895e+02,   2.3161721e+02,   2.3278821e+02,   2.3397218e+02, &
                0738       2.3516935e+02,   2.3637994e+02,   2.3760420e+02,   2.3884238e+02, &
                0739       2.4009473e+02,   2.4136150e+02,   2.4264297e+02,   2.4393941e+02, &
                0740       2.4525110e+02,   2.4657831e+02,   2.4792136e+02,   2.4928053e+02, &
                0741       2.5065615e+02,   2.5204853e+02,   2.5345799e+02,   2.5488487e+02, &
                0742       2.5632953e+02,   2.5779231e+02,   2.5927358e+02,   2.6077372e+02, &
                0743       2.6229310e+02,   2.6383214e+02,   2.6539124e+02,   2.6697081e+02, &
                0744       2.6857130e+02,   2.7019315e+02,   2.7183682e+02,   2.7350278e+02, &
                0745       2.7519152e+02,   2.7690354e+02,   2.7863937e+02,   2.8039954e+02, &
                0746       2.8218459e+02,   2.8399511e+02,   2.8583167e+02,   2.8769489e+02, &
                0747       2.8958539e+02,   2.9150383e+02,   2.9345086e+02,   2.9542719e+02, &
                0748       2.9743353e+02,   2.9947061e+02,   3.0153922e+02,   3.0364014e+02, &
                0749       3.0577420e+02,   3.0794224e+02,   3.1014515e+02,   3.1238386e+02, &
                0750       3.1465930e+02,   3.1697246e+02,   3.1932437e+02,   3.2171609e+02, &
                0751       3.2414873e+02,   3.2662343e+02,   3.2914139e+02,   3.3170385e+02 /
                0752
                0753 !     v1 = value
                0754       v1 = MIN( MAX( value, -23.d0 ), -10.4d0 )
                0755 !     if (value.lt.-23.0) call error_mesg ('dargan_bettsmiller',  &
                0756 !                        'lcltable: value too low', FATAL)
                0757
                0758 !     if (value.gt.-10.4) call error_mesg ('dargan_bettsmiller',  &
                0759 !                        'lcltable: value too high', FATAL)
                0760       if (value.lt.-23.0) then
                0761           CALL PRINT_ERROR( 'In: dargan_bettsmiller '//  &
                0762                             'lcltable: value too low', myThid)
892dca4bcd Jean*0763           WRITE(errorMessageUnit,'(A,4I4,I6,1PE14.6)')   &
                0764             'i,j,bi,bj,myIter,value=',i,j,bi,bj,myIter,value
b2ea1d2979 Jean*0765 !         STOP 'ABNORMAL END: S/R LCLTABL (dargan_bettsmiller_mod)'
                0766       endif
                0767       if (value.gt.-10.4) then
                0768           CALL PRINT_ERROR( 'In: dargan_bettsmiller '//  &
                0769                             'lcltable: value too high', myThid)
892dca4bcd Jean*0770           WRITE(errorMessageUnit,'(A,4I4,I6,1PE14.6)')   &
                0771             'i,j,bi,bj,myIter,value=',i,j,bi,bj,myIter,value
b2ea1d2979 Jean*0772 !         STOP 'ABNORMAL END: S/R LCLTABL (dargan_bettsmiller_mod)'
                0773       endif
                0774       ival = floor(10.*(v1 + 23.0))
                0775       v2 = -230. + ival
                0776       v1 = 10.*v1
                0777       tlcl = (v2 + 1.0 - v1)*lcltable(ival+1) + (v1 - v2)*lcltable(ival+2)
                0778
                0779       end subroutine lcltabl
                0780
                0781 !#######################################################################
                0782
                0783    subroutine dargan_bettsmiller_init (myThid)
                0784
                0785 !-----------------------------------------------------------------------
                0786 !
                0787 !        initialization for bettsmiller
                0788 !
                0789 !-----------------------------------------------------------------------
                0790
599643daec Jean*0791 ! integer  unit,io,ierr
b2ea1d2979 Jean*0792   integer, intent(in) ::myThid
                0793 !-------------------------------------------------------------------------------------
599643daec Jean*0794 !integer, dimension(3) :: half = (/1,2,4/)
b2ea1d2979 Jean*0795 !integer :: ierr, io
                0796 integer         :: iUnit
                0797 CHARACTER*(gcm_LEN_MBUF) :: msgBuf
                0798 !-------------------------------------------------------------------------------
                0799
                0800 !----------- read namelist ---------------------------------------------
                0801
                0802 !      if (file_exist('input.nml')) then
                0803 !         unit = open_file (file='input.nml', action='read')
                0804 !         ierr=1; do while (ierr /= 0)
                0805 !            read  (unit, nml=dargan_bettsmiller_nml, iostat=io, end=10)
                0806 !            ierr = check_nml_error (io,'dargan_bettsmiller_nml')
                0807 !         enddo
                0808 !  10     call close_file (unit)
                0809 !      endif
                0810
                0811 !---------- output namelist --------------------------------------------
                0812
                0813 !      unit = open_file (file='logfile.out', action='append')
                0814 !      if ( mpp_pe() == 0 ) then
                0815 !           write (unit,'(/,80("="),/(a))') trim(version), trim(tag)
                0816 !           write (unit,nml=dargan_bettsmiller_nml)
                0817 !      endif
                0818 !      call close_file (unit)
                0819
                0820 !      do_init=.false.
                0821
                0822 !    _BARRIER
                0823 !    _BEGIN_MASTER(myThid)
                0824      CALL BARRIER(myThid)
                0825      IF ( myThid.EQ.1 ) THEN
                0826
                0827      WRITE(msgBuf,'(A)') 'DARGAN_BETTSMILLER_init: opening data.atm_gray'
                0828      CALL PRINT_MESSAGE( msgBuf, gcm_stdMsgUnit, gcm_SQZ_R, myThid )
                0829      CALL OPEN_COPY_DATA_FILE(                                      &
                0830                            'data.atm_gray', 'DARGAN_BETTSMILLER_init',       &
                0831                            iUnit,                                   &
                0832                            myThid )
                0833 !    Read parameters from open data file
                0834      READ(UNIT=iUnit,NML=dargan_bettsmiller_nml)
                0835      WRITE(msgBuf,'(A)')                                            &
                0836           'DARGAB_BETTSMILLER_INIT: finished reading data.atm_gray'
                0837      CALL PRINT_MESSAGE( msgBuf, gcm_stdMsgUnit, gcm_SQZ_R, myThid )
                0838 !    Close the open data file
15388b052e Jean*0839 #ifdef SINGLE_DISK_IO
b2ea1d2979 Jean*0840      CLOSE(iUnit)
15388b052e Jean*0841 #else
                0842      CLOSE(iUnit,STATUS='DELETE')
                0843 #endif /* SINGLE_DISK_IO */
b2ea1d2979 Jean*0844
                0845      ENDIF
                0846      CALL BARRIER(myThid)
                0847
                0848    return
                0849    end subroutine dargan_bettsmiller_init
                0850
                0851 !#######################################################################
                0852
                0853  real function mixing_ratio(vapor_pressure, pressure)
                0854
                0855  ! calculates the mixing ratio from the vapor pressure and pressure
                0856
                0857       implicit none
                0858       real, intent(in)     :: vapor_pressure, pressure
                0859
                0860       mixing_ratio = rdgas*vapor_pressure/rvgas/(pressure-vapor_pressure)
                0861
                0862  end function mixing_ratio
                0863
                0864 !#######################################################################
                0865
                0866  real function virtual_temp(temp, r)
                0867
                0868  ! calculates the virtual temperature from the temperature and mixing ratio
                0869  ! consistent with the approximation used in the fms code
                0870
                0871       implicit none
                0872       real, intent(in)     :: temp         ! temperature
                0873       real, intent(in)     :: r            ! mixing ratio
                0874
                0875       real                 :: q            ! specific humidity
                0876
                0877       if (do_virtual) then
                0878        q = r/(1.0+r)
                0879        virtual_temp = temp*(1.0+q*(rvgas/rdgas-1.0))
                0880       else
                0881        virtual_temp = temp
                0882       endif
                0883
                0884  end function virtual_temp
                0885
                0886 !#######################################################################
                0887
                0888 end module dargan_bettsmiller_mod
                0889