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