da_moist_phys_lin.inc
References to this file elsewhere.
1 subroutine da_moist_phys_lin(grid)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Partition of the hydrometeors via the moist explicit scheme.
5 ! A warm rain process is used in this subroutine.
6 ! This is the tangent linear code of the scheme.
7 !
8 ! Method: The warm rain process is according to Hsie and Anthes (1984)
9 ! and Dudhia (1989)
10 !
11 ! Assumptions: 1) Model level stored top down.
12 !---------------------------------------------------------------------------
13
14 implicit none
15
16 type(domain), intent(inout) :: grid
17
18 real :: T_OLD(ims:ime,jms:jme,kms:kme),T_NEW(ims:ime,jms:jme,kms:kme)
19 real :: Q_OLD(ims:ime,jms:jme,kms:kme),Q_NEW(ims:ime,jms:jme,kms:kme)
20 real :: QCW_OLD(ims:ime,jms:jme,kms:kme),QCW_NEW(ims:ime,jms:jme,kms:kme)
21 real :: QRN_OLD(ims:ime,jms:jme,kms:kme),QRN_NEW(ims:ime,jms:jme,kms:kme)
22
23 real :: EES(kms:kme)
24 real :: QVSS(kms:kme)
25 real :: EES9(kms:kme)
26 real :: QVSS9(kms:kme)
27 real :: DT(its:ite,jts:jte,kms:kme)
28 real :: QVT(kms:kme)
29 real :: QCT(kms:kme)
30 real :: QRT(kms:kme)
31 real :: TTT(kms:kme)
32 real :: QVT9(kms:kme)
33 real :: QCT9(kms:kme)
34 real :: QRT9(kms:kme)
35 real :: TTT9(kms:kme)
36 real :: SCR2(kms:kme)
37 real :: SCR3(kms:kme)
38 real :: SCR4(kms:kme)
39 real :: SCR5(kms:kme)
40 real :: SCR6(kms:kme)
41 real :: DUM31(kms:kme)
42 real :: PRA(kms:kme)
43 real :: PRC(kms:kme)
44 real :: PRD(kms:kme)
45 real :: PRE(kms:kme)
46 real :: SCR31(kms:kme)
47 real :: SCR42(kms:kme)
48 real :: SCR71(kms:kme)
49 real :: DUM112(kms:kme)
50 real :: DUM113(kms:kme)
51 real :: DUM211(kms:kme)
52 real :: DUM411(kms:kme)
53 real :: SCR29(kms:kme)
54 real :: SCR39(kms:kme)
55 real :: SCR49(kms:kme)
56 real :: SCR59(kms:kme)
57 real :: SCR69(kms:kme)
58 real :: DUM319(kms:kme)
59 real :: PRA9(kms:kme)
60 real :: PRC9(kms:kme)
61 real :: PRD9(kms:kme)
62 real :: PRE9(kms:kme)
63 real :: SCR319(kms:kme)
64 real :: SCR429(kms:kme)
65 real :: SCR719(kms:kme)
66 real :: DUM1129(kms:kme)
67 real :: DUM4119(kms:kme)
68 real :: TMP(kms:kme)
69
70 integer, parameter :: qcth = 0.5e-3
71 integer, parameter :: qrth = 1.0e-6
72 integer, parameter :: alpha = 0.001
73 integer, parameter :: beta = 0.0486
74 integer, parameter :: gamma = 0.002
75
76 integer :: i, j, k
77
78 real :: qrth1
79
80 if (trace_use) call da_trace_entry("da_moist_phys_lin")
81
82 qrth1 = (QRTH*1.0e3)**0.875
83
84 T_OLD (its:ite,jts:jte,kts:kte) = grid%xa%t (its:ite,jts:jte,kts:kte)
85 Q_OLD (its:ite,jts:jte,kts:kte) = grid%xa%q (its:ite,jts:jte,kts:kte)
86 QCW_OLD(its:ite,jts:jte,kts:kte) = grid%xa%qcw(its:ite,jts:jte,kts:kte)
87 QRN_OLD(its:ite,jts:jte,kts:kte) = grid%xa%qrn(its:ite,jts:jte,kts:kte)
88
89 ! Preparation
90
91 grid%xa%q(its:ite,jts:jte,kts:kte) =grid%xa%qt(its:ite,jts:jte,kts:kte) - grid%xa%qcw(its:ite,jts:jte,kts:kte) - grid%xa %qrn(its:ite,jts:jte,kts:kte)
92 DT(its:ite,jts:jte,kts:kte) = grid%xb%delt(its:ite,jts:jte,kts:kte)
93
94 do j=jts,jte
95 do i=its,ite
96 do K=kts,kte
97
98 if (dt(i,j,k) <= 0.0) cycle
99
100 if ( grid%xb%t(I,J,K) > TO )then
101 EES(K)=SVP1*EXP(SVP2*(grid%xb%t(I,J,K)-SVPT0)/(grid%xb%t(I,J,K)-SVP3))
102 EES9(K)=EES(K)*SVP2*(SVPT0-SVP3)/((grid%xb%t(I,J,K)-SVP3) * (grid%xb%t(I,J,K)-SVP3))*grid%xa%t(I,J,K)
103 else
104 EES(K)=.611*EXP(22.514-6.15E3/grid%xb%t(I,J,K))
105 EES9(K)=EES(K)*6.15E3/(grid%xb%t(I,J,K)*grid%xb%t(I,J,K))*grid%xa%t(I,J,K)
106 end if
107
108 TMP(K)=622.0/((grid%xb%p(I,J,K)-EES(K))**2)
109 QVSS9(K)=TMP(K)*grid%xb%p(I,J,K)*EES9(K) - TMP(K)*EES(K)*grid%xa%p(I,J,K)
110 QVSS(K)=622.0*EES(K)/(grid%xb%p(I,J,K)-EES(K))
111
112 SCR49(K)=grid%xa%q(I,J,K)/QVSS(K)-grid%xb%q(I,J,K)/QVSS(K)**2*QVSS9(K)
113 SCR4(K)=grid%xb%q(I,J,K)/QVSS(K)
114
115 if (grid%xb%qcw(I,J,K) > 0.0) then
116 SCR29(K)=grid%xa%qcw(I,J,K)
117 SCR2(K)=grid%xb%qcw(I,J,K)
118 else
119 SCR29(K)=0.0
120 SCR2(K)=0.0
121 end if
122 if (grid%xb%qrn(I,J,K) > 1.0e-25) then
123 SCR39(K)=grid%xa%qrn(I,J,K)
124 SCR3(K)=grid%xb%qrn(I,J,K)
125 else
126 SCR39(K)=0.0
127 SCR3(K)=1.0E-25
128 end if
129 SCR59(K)=grid%xa%q(I,J,K)/SCR4(K)-grid%xb%q(I,J,K)/SCR4(K)**2*SCR49(K)
130 SCR5(K)=grid%xb%q(I,J,K)/SCR4(K)
131
132 SCR69(K)=grid%xa%p(I,J,K)/(gas_constant*grid%xb%t(I,J,K))-grid%xb%p(I,J,K)/ &
133 (gas_constant*grid%xb%t(I,J,K)**2)*grid%xa%t(I,J,K)
134 SCR6(K)=grid%xb%p(I,J,K)/(gas_constant*grid%xb%t(I,J,K))
135
136 DUM319(K)=-XLV1*grid%xa%t(I,J,K)
137 DUM31(K)=3.1484E6-XLV1*grid%xb%t(I,J,K)
138
139 ! Auto conversion
140
141 if (scr2(k) >= qcth) then
142 prc9(k) = alpha * scr29(k)
143 prc(k) = alpha * (scr2(k) - qcth)
144 else
145 prc9(k) = 0.0
146 prc(k) = 0.0
147 end if
148
149 ! Accretion
150
151 if (SCR2(k) > 0.0 .and. SCR3(k) > QRTH ) then
152 PRA9(K) = gamma * 0.875 * SCR2(k) * (SCR3(K)*1.0e3)**(-0.125) * 1.0e3 * SCR39(K) &
153 + gamma * SCR29(k) * (SCR3(K)*1.0e3)**0.875
154 PRA(k) = gamma * SCR2(k) * (SCR3(k)*1.0e3)**0.875
155 else if (SCR2(k) > 0.0 .and. SCR3(k) <= QRTH ) then
156 PRA9(K) = gamma * SCR29(k) * qrth1
157 PRA(k) = gamma * SCR2(k) * qrth1
158 else
159 PRA9(K) = 0.0
160 PRA(k) = 0.0
161 end if
162
163 end do
164
165 call da_evapo_lin(DT(i,j,:),SCR3,SCR5,grid%xb%q(I,J,:),PRE,SCR6, &
166 SCR39,SCR59,grid%xa%q(I,J,:),PRE9,SCR69, &
167 kts,kte,kms,kme)
168
169 do K=kts, kte
170
171 if (dt(i,j,k) <= 0.0) cycle
172
173 ! Readjust
174
175 DUM112(K)=(PRC(k)+PRA(k))*dt(i,j,k)
176 if (DUM112(K) > SCR2(k)) then
177 DUM1129(K)=(PRC9(k)+PRA9(k))*dt(i,j,k)
178 PRC9(K)=SCR29(K)*PRC(K)/DUM112(K) &
179 +PRC9(K)*SCR2(K)/DUM112(K) &
180 -SCR2(K)*PRC(K)/(DUM112(K)*DUM112(K))*DUM1129(K)
181 PRC(k)=SCR2(K)*PRC(K)/DUM112(K)
182 PRA9(K)=SCR29(K)*PRA(K)/DUM112(K) &
183 +PRA9(K)*SCR2(K)/DUM112(K) &
184 -SCR2(K)*PRA(K)/(DUM112(K)*DUM112(K))*DUM1129(K)
185 PRA(k)=SCR2(K)*PRA(K)/DUM112(K)
186 end if
187 QVT9(K)=-PRE9(K)
188 QVT(K)=-PRE(K)
189 QCT9(K)=-PRC9(K)-PRA9(K)
190 QCT(K)=-PRC(K)-PRA(K)
191 QRT9(K)=PRC9(K)+PRA9(K)+PRE9(K)
192 QRT(K)=PRC(K)+PRA(K)+PRE(K)
193 if (grid%xb%t(I,J,K).GT.TO)then
194 DUM4119(K)=DUM319(K)
195 DUM411(K)=DUM31(K)
196 else
197 DUM4119(K)=0.0
198 DUM411(K)=XLS
199 end if
200 PRD9(K)=cp*0.887*grid%xa%q(I,J,K)
201 PRD(K)=cp*(1.0+0.887*grid%xb%q(I,J,K))
202 TTT9(K)=-DUM4119(K)*QVT(K)/PRD(K) &
203 -QVT9(K)*DUM411(K)/PRD(K) &
204 +DUM411(K)*QVT(K)/(PRD(K)*PRD(K))*PRD9(K)
205 TTT(K)=-DUM411(K)*QVT(K)/PRD(K)
206
207 DUM113(K)=grid%xb%q(I,J,K)+dt(i,j,k)*QVT(K)
208 if (DUM113(K) > 1.0e-12 ) then
209 SCR429(K)=grid%xa%q(I,J,K)+dt(i,j,k)*QVT9(K)
210 SCR42(K)=DUM113(K)
211 else
212 SCR429(K)=0.0
213 SCR42(K)=1.0e-12
214 end if
215 DUM211(K)=grid%xb%qcw(I,J,K)+QCT(K)*dt(i,j,k)
216 if (DUM211(K) > 0.0) then
217 SCR319(K)=grid%xa%qcw(I,J,K)+QCT9(K)*dt(i,j,k)
218 SCR31(K)=DUM211(K)
219 else
220 SCR319(K)=0.0
221 SCR31(K)=0.0
222 end if
223 SCR719(K)=grid%xa%t(I,J,K)+TTT9(K)*dt(i,j,k)
224 SCR71(K)=grid%xb%t(I,J,K)+TTT(K)*dt(i,j,k)
225 end do
226
227 call da_condens_lin(DT(i,j,:),SCR31,SCR42,SCR71,DUM31,PRD, &
228 QVT,QCT,QRT,TTT, &
229 grid%xb%p(I,J,:),grid%xb%t(I,J,:),grid%xb%q(I,J,:), &
230 grid%xb%qcw(I,J,:),grid%xb%qrn(I,J,:), &
231 SCR319,SCR429,SCR719,DUM319,PRD9, &
232 QVT9,QCT9,QRT9,TTT9, &
233 grid%xa%p(I,J,:),grid%xa%t(I,J,:),grid%xa%q(I,J,:), &
234 grid%xa%qcw(I,J,:),grid%xa%qrn(I,J,:),kts,kte)
235 end do
236 end do
237
238 T_NEW (its:ite,jts:jte,kds:kde) = grid%xa%t (its:ite,jts:jte,kds:kde) - T_OLD (its:ite,jts:jte,kds:kde)
239 Q_NEW (its:ite,jts:jte,kds:kde) = grid%xa%q (its:ite,jts:jte,kds:kde) - Q_OLD (its:ite,jts:jte,kds:kde)
240 QCW_NEW(its:ite,jts:jte,kds:kde) = grid%xa%qcw (its:ite,jts:jte,kds:kde) - QCW_OLD(its:ite,jts:jte,kds:kde)
241 QRN_NEW(its:ite,jts:jte,kds:kde) = grid%xa%qrn (its:ite,jts:jte,kds:kde) - QRN_OLD(its:ite,jts:jte,kds:kde)
242
243 call da_filter(grid, t_new)
244 call da_filter(grid, q_new)
245 call da_filter(grid, qcw_new)
246 call da_filter(grid, qrn_new)
247
248 grid%xa%t (its:ite,jts:jte,kds:kde) = T_NEW (its:ite,jts:jte,kds:kde) + T_OLD (its:ite,jts:jte,kds:kde)
249 grid%xa%q (its:ite,jts:jte,kds:kde) = Q_NEW (its:ite,jts:jte,kds:kde) + Q_OLD (its:ite,jts:jte,kds:kde)
250 grid%xa%qcw (its:ite,jts:jte,kds:kde) = QCW_NEW(its:ite,jts:jte,kds:kde) + QCW_OLD(its:ite,jts:jte,kds:kde)
251 grid%xa%qrn (its:ite,jts:jte,kds:kde) = QRN_NEW(its:ite,jts:jte,kds:kde) + QRN_OLD(its:ite,jts:jte,kds:kde)
252
253 if (trace_use) call da_trace_exit("da_moist_phys_lin")
254
255 end subroutine da_moist_phys_lin
256
257