var/da/da_setup_structures/da_cumulus.inc

   1 subroutine da_cumulus (zcb, tcb, qcb, pcb, pk, te, z, t, q, lcb, lct, pct, zct, kts, kte)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: TBD
   5    !-----------------------------------------------------------------------
   6
   7    implicit none
   8
   9    integer, intent(in)    :: kts, kte
  10    real,    intent(inout) :: zcb, tcb, qcb, pcb
  11    real,    intent(in)    :: pk(kts:kte)
  12    real,    intent(in)    :: te(kts:kte)
  13    real,    intent(out)   :: z(kts:kte)
  14    real,    intent(out)   :: t(kts:kte)
  15    real,    intent(out)   :: q(kts:kte)
  16    integer, intent(out)   :: lcb, lct
  17    real,    intent(out)   :: pct, zct
  18
  19    integer   :: k, ia, l, ncb
  20    real      :: cp, r, hl, em, et, p
  21    real      :: tll, qll, pll, zll, tbar, pbar, qbar
  22    real      :: dp, dz, ddt, dt
  23
  24    if (trace_use) call da_trace_entry("da_cumulus")
  25
  26    cp=1004.0
  27    r=2000.0/7.0
  28    hl=2.49e06
  29    dt=0.1
  30    ia=1000
  31
  32    do k = kts, kte
  33       z(k) = 0.0
  34       t(k) = 0.0
  35       q(k) = 0.0
  36    end do
  37
  38    em=gravity*zcb+cp*tcb+hl*qcb
  39
  40    ncb=kts
  41
  42    if (pk(kte) > pcb) then
  43       ncb=kte
  44    end if
  45
  46    do l=kte-1,kts,-1
  47       if (pk(l) > pcb) then
  48          ncb=l+1
  49          exit
  50       end if
  51    end do
  52
  53    do l=ncb,kte
  54       p=pk(l)
  55       do k=1,ia
  56          if (l == ncb) then
  57             tll=tcb
  58             qll=qcb
  59             pll=pcb
  60             zll=zcb
  61          else
  62             tll=t(l-1)
  63             qll=q(l-1)
  64             pll=pk(l-1)
  65             zll=z(l-1)
  66          end if
  67
  68          t(l)=tll-(k*dt)
  69
  70          call da_qfrmrh(p, t(l), 100.0, q(l))
  71
  72          tbar=0.5*(t(l)+tll)
  73          qbar=0.5*(q(l)+qll)
  74          pbar=0.5*(p+pll)
  75          dp=pll-p
  76          dz=(r*tbar*(1.0+0.61*qbar)*dp)/(gravity*pbar)
  77          z(l)=zll+dz
  78          et=gravity*z(l)+cp*t(l)+hl*q(l)
  79          if ((et-em) <= 0.0) exit
  80       end do
  81    end do
  82
  83    lct=ncb
  84
  85    do k=kte,ncb+1,-1
  86       ddt=t(k)-te(k)
  87
  88       if (ddt >= 0.0) then
  89          lct=k
  90          exit
  91       end if
  92    end do
  93
  94    lcb=lct
  95
  96    do k=ncb,kte
  97       ddt=t(k)-te(k)
  98       if (ddt >= 0.0) then
  99          lcb=k
 100          exit
 101       end if
 102    end do
 103
 104    pct=pk(lct)
 105    zct=z(lct)
 106    pcb=pk(lcb)
 107    zcb=z(lcb)
 108
 109    if (trace_use) call da_trace_exit("da_cumulus")
 110
 111 end subroutine da_cumulus
 112
 113