da_moist_phys_adj.inc
References to this file elsewhere.
1 subroutine da_moist_phys_adj(xb, xa, xp)
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 adjoint 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 (xb_type), intent(in) :: xb ! First guess structure.
17 type (x_type), intent(inout) :: xa ! increment structure.
18 type (xpose_type), intent(inout) :: xp ! Dimensions and xpose buffers.
19
20 real, dimension(ims:ime,jms:jme,kms:kme) :: T_OLD,T_NEW
21 real, dimension(ims:ime,jms:jme,kms:kme) :: Q_OLD,Q_NEW
22 real, dimension(ims:ime,jms:jme,kms:kme) :: QCW_OLD,QCW_NEW
23 real, dimension(ims:ime,jms:jme,kms:kme) :: QRN_OLD,QRN_NEW
24
25 real, dimension(kms:kme) :: EES, QVSS
26 real, dimension(kms:kme) :: EES9, QVSS9
27
28 real, dimension(kms:kme) :: DT
29 real, dimension(kms:kme) :: QVT,QCT,QRT,TTT
30 real, dimension(kms:kme) :: QVT9,QCT9,QRT9,TTT9
31 real, dimension(kms:kme) :: SCR2,SCR3,SCR4,SCR5,SCR6
32 real, dimension(kms:kme) :: DUM31
33 real, dimension(kms:kme) :: PRA,PRC,PRD,PRE
34 real, dimension(kms:kme) :: SCR31,SCR42,SCR71
35 real, dimension(kms:kme) :: DUM112,DUM113,DUM211,DUM411
36 real, dimension(kms:kme) :: PRA2,PRC2
37
38 real, dimension(kms:kme) :: SCR29,SCR39,SCR49,SCR59,SCR69
39 real, dimension(kms:kme) :: DUM319
40 real, dimension(kms:kme) :: PRA9,PRC9,PRD9,PRE9
41 real, dimension(kms:kme) :: SCR319,SCR429,SCR719
42 real, dimension(kms:kme) :: DUM1129,DUM1139,DUM2119,DUM4119
43 real, dimension(kms:kme) :: TMP
44
45 integer :: i, j, k
46
47 if (trace_use) call da_trace_entry("da_moist_phys_adj")
48
49 do k=kts,kte
50 do j=jts,jte
51 do i=its,ite
52 QRN_NEW(i,j,k) = xa % qrn (i,j,k)
53 QRN_OLD(i,j,k) = xa % qrn (i,j,k)
54 QCW_NEW(i,j,k) = xa % qcw (i,j,k)
55 QCW_OLD(i,j,k) = xa % qcw (i,j,k)
56 Q_NEW(i,j,k) = xa % q (i,j,k)
57 Q_OLD(i,j,k) = xa % q (i,j,k)
58 T_NEW(i,j,k) = xa % t (i,j,k)
59 T_OLD(i,j,k) = xa % t (i,j,k)
60 end do
61 end do
62 end do
63
64 call da_filter_adj(t_new, xp)
65 call da_filter_adj(q_new, xp)
66 call da_filter_adj(qcw_new, xp)
67 call da_filter_adj(qrn_new, xp)
68
69 do k=kts,kte
70 do j=jts,jte
71 do i=its,ite
72 xa % qrn (i,j,k) = QRN_NEW(i,j,k)
73 QRN_OLD(i,j,k) = QRN_OLD(i,j,k) - QRN_NEW(i,j,k)
74 xa % qcw (i,j,k) = QCW_NEW(i,j,k)
75 QCW_OLD(i,j,k) = QCW_OLD(i,j,k) - QCW_NEW(i,j,k)
76 xa % q (i,j,k) = Q_NEW(i,j,k)
77 Q_OLD(i,j,k) = Q_OLD(i,j,k) - Q_NEW(i,j,k)
78 xa % t (i,j,k) = T_NEW(i,j,k)
79 T_OLD(i,j,k) = T_OLD(i,j,k) - T_NEW(i,j,k)
80 end do
81 end do
82 end do
83
84 do j=jts,jte
85 do i=its,ite
86
87 ! Preparation
88
89 do K=kts,kte
90 DT(k) = xb%delt(i,j,k)
91 end do
92
93 do K=kts,kte
94
95 if (DT(k) <= 0.) cycle
96
97 if( xb%t(I,J,K) > TO )then
98 EES(K)=SVP1*EXP(SVP2*(xb%t(I,J,K)-SVPT0)/(xb%t(I,J,K)-SVP3))
99 else
100 EES(K)=.611*EXP(22.514-6.15E3/xb%t(I,J,K))
101 end if
102
103 QVSS(K)=622.*EES(K)/(xb%p(I,J,K)-EES(K))
104
105 SCR4(K)=xb%q(I,J,K)/QVSS(K)
106
107 end do
108
109 do K=kts,kte
110
111 if (DT(k) <= 0.) cycle
112
113 if(xb%qcw(I,J,K) > 0.) then
114 SCR2(K)=xb%qcw(I,J,K)
115 else
116 SCR2(K)=0.
117 end if
118 if(xb%qrn(I,J,K) > 1.E-25) then
119 SCR3(K)=xb%qrn(I,J,K)
120 else
121 SCR3(K)=1.E-25
122 end if
123 SCR5(K)=xb%q(I,J,K)/SCR4(K)
124
125 SCR6(K)=xb%p(I,J,K)/(gas_constant*xb%t(I,J,K))
126
127 DUM31(K)=3.1484E6-XLV1*xb%t(I,J,K)
128
129 end do
130
131 call da_autoc(DT,SCR2,PRC,kts,kte,kms,kme)
132
133 do k=kts,kte
134 PRC2(K)=PRC(K)
135 end do
136
137 call da_accre(DT,SCR2,SCR3,PRA,kts,kte,kms,kme)
138
139 do k=kts,kte
140 PRA2(K)=PRA(K)
141 end do
142
143 call da_evapo(DT,SCR3,SCR5,xb%q(I,J,:),PRE,SCR6,kts,kte,kms,kme)
144
145 do K=kts, kte
146
147 if (DT(k) <= 0.) cycle
148
149 ! Readjust
150
151 DUM112(K)=(PRC(k)+PRA(k))*DT(k)
152 if (DUM112(K) > SCR2(k)) then
153 PRC(k)=SCR2(K)*PRC(K)/DUM112(K)
154 PRA(k)=SCR2(K)*PRA(K)/DUM112(K)
155 end if
156 QVT(K)=-PRE(K)
157 QCT(K)=-PRC(K)-PRA(K)
158 QRT(K)=PRC(K)+PRA(K)+PRE(K)
159 if(xb%t(I,J,K).GT.TO)then
160 DUM411(K)=DUM31(K)
161 else
162 DUM411(K)=XLS
163 end if
164 PRD(K)=cp*(1.+0.887*xb%q(I,J,K))
165 TTT(K)=-DUM411(K)*QVT(K)/PRD(K)
166
167 DUM113(K)=xb%q(I,J,K)+DT(K)*QVT(K)
168 if(DUM113(K) > 1.E-12 ) then
169 SCR42(K)=DUM113(K)
170 else
171 SCR42(K)=1.E-12
172 end if
173 DUM211(K)=xb%qcw(I,J,K)+QCT(K)*DT(k)
174 if(DUM211(K) > 0.) then
175 SCR31(K)=DUM211(K)
176 else
177 SCR31(K)=0.
178 end if
179 SCR71(K)=xb%t(I,J,K)+TTT(K)*DT(k)
180 end do
181
182 call da_condens_adj(DT,SCR31,SCR42,SCR71,DUM31,PRD, &
183 QVT,QCT,QRT,TTT, &
184 xb%p(I,J,:),xb%t(I,J,:),xb%q(I,J,:), &
185 xb%qcw(I,J,:),xb%qrn(I,J,:), &
186 SCR319,SCR429,SCR719,DUM319,PRD9, &
187 QVT9,QCT9,QRT9,TTT9, &
188 xa%p(I,J,:),xa%t(I,J,:),xa%q(I,J,:), &
189 xa%qcw(I,J,:),xa%qrn(I,J,:),kts,kte)
190
191 do K=kts, kte
192 if (DT(k) <= 0.) cycle
193
194 ! Readjust
195
196 xa%t(I,J,K)=xa%t(I,J,K)+SCR719(K)
197 TTT9(K)=TTT9(K)+DT(k)*SCR719(K)
198 DUM2119(K)=0.
199 if(DUM211(K) > 0.) then
200 DUM2119(K)=SCR319(K)
201 end if
202 xa%qcw(I,J,K)=xa%qcw(I,J,K)+DUM2119(K)
203 QCT9(K)=QCT9(K)+DT(K)*DUM2119(K)
204 DUM1139(K)=0.
205 if(DUM113(K) > 1.E-12 ) then
206 DUM1139(K)=SCR429(K)
207 end if
208 xa%q(I,J,K)=xa%q(I,J,K)+DUM1139(K)
209 QVT9(K)=QVT9(K)+DT(K)*DUM1139(K)
210 PRD9(K)=PRD9(K)+DUM411(K)*QVT(K)/(PRD(K)*PRD(K))*TTT9(K)
211 QVT9(K)=QVT9(K)-DUM411(K)/PRD(K)*TTT9(K)
212 DUM4119(K)=-QVT(K)/PRD(K)*TTT9(K)
213 xa%q(I,J,K)=xa%q(I,J,K)+cp*0.887*PRD9(K)
214
215 if(xb%t(I,J,K).GT.TO)then
216 DUM319(K)=DUM319(K)+DUM4119(K)
217 end if
218 PRC9(K)=QRT9(K)
219 PRA9(K)=QRT9(K)
220 PRE9(K)=QRT9(K)
221 PRC9(K)=PRC9(K)-QCT9(K)
222 PRA9(K)=PRA9(K)-QCT9(K)
223 PRE9(K)=PRE9(K)-QVT9(K)
224
225 SCR29(K)=0.
226 DUM1129(K)=0.
227 if (DUM112(K) > SCR2(k)) then
228 DUM1129(K)=-SCR2(K)*PRA2(K)/(DUM112(K)*DUM112(K))*PRA9(K)
229 SCR29(K)=PRA2(K)/DUM112(K)*PRA9(K)
230 PRA9(K)=PRA9(K)*SCR2(K)/DUM112(K)
231 DUM1129(K)=DUM1129(K)-SCR2(K)*PRC2(K)/(DUM112(K)*DUM112(K))*PRC9(K)
232 SCR29(K)=SCR29(K)+PRC2(K)/DUM112(K)*PRC9(K)
233 PRC9(K)=PRC9(K)*SCR2(K)/DUM112(K)
234 end if
235 PRC9(K)=PRC9(K)+DT(K)*DUM1129(K)
236 PRA9(K)=PRA9(K)+DT(K)*DUM1129(K)
237
238 end do
239
240 call da_evapo_adj(DT,SCR3,SCR5,xb%q(I,J,:),PRE,SCR6, &
241 SCR39,SCR59,xa%q(I,J,:),PRE9,SCR69, &
242 kts,kte,kms,kme)
243
244 call da_accre_adj(DT,SCR2,SCR3,PRA,SCR29,SCR39,PRA9,kts,kte,kms,kme)
245
246 call da_autoc_adj(DT,SCR2,PRC,SCR29,PRC9,kts,kte,kms,kme)
247
248 ! Preparation
249
250 do K=kts,kte
251
252 if (DT(k) <= 0.) cycle
253
254 xa%t(I,J,K)=xa%t(I,J,K)-XLV1*DUM319(K)
255
256 xa%p(I,J,K)=xa%p(I,J,K)+SCR69(K)/(gas_constant*xb%t(I,J,K))
257 xa%t(I,J,K)=xa%t(I,J,K)-xb%p(I,J,K)/ &
258 (gas_constant*xb%t(I,J,K)**2)*SCR69(K)
259 xa%q(I,J,K)=xa%q(I,J,K)+SCR59(K)/SCR4(K)
260 SCR49(K)=-xb%q(I,J,K)/SCR4(K)**2*SCR59(K)
261
262 if(xb%qrn(I,J,K) > 1.E-25) then
263 xa%qrn(I,J,K)=xa%qrn(I,J,K)+SCR39(K)
264 end if
265 if(xb%qcw(I,J,K) > 0.) then
266 xa%qcw(I,J,K)=xa%qcw(I,J,K)+SCR29(K)
267 end if
268
269 end do
270
271 do K=kts,kte
272
273 if (DT(k) <= 0.) cycle
274
275 xa%q(I,J,K)=xa%q(I,J,K)+SCR49(K)/QVSS(K)
276 QVSS9(K)=-xb%q(I,J,K)/QVSS(K)**2*SCR49(K)
277 TMP(K)=622./((xb%p(I,J,K)-EES(K))**2)
278 EES9(K)=TMP(K)*xb%p(I,J,K)*QVSS9(K)
279 xa%p(I,J,K)=xa%p(I,J,K)-TMP(K)*EES(K)*QVSS9(K)
280 if( xb%t(I,J,K) > TO )then
281 xa%t(I,J,K)=xa%t(I,J,K)+EES(K)*SVP2*(SVPT0-SVP3)/ &
282 ((xb%t(I,J,K)-SVP3)*(xb%t(I,J,K)-SVP3))*EES9(K)
283 else
284 xa%t(I,J,K)=xa%t(I,J,K)+EES(K)*6.15E3/(xb%t(I,J,K)* &
285 xb%t(I,J,K))*EES9(K)
286 end if
287
288 end do
289
290 do K=kts,kte
291 xa % qt (i,j,k) = xa % qt (i,j,k) + xa % q(i,j,k)
292 xa % qcw(i,j,k) = xa % qcw(i,j,k) - xa % q(i,j,k)
293 xa % qrn(i,j,k) = xa % qrn(i,j,k) - xa % q(i,j,k)
294 end do
295 end do
296 end do
297
298 do k=kts,kte
299 do j=jts,jte
300 do i=its,ite
301 xa % qrn (i,j,k) = xa % qrn (i,j,k) + QRN_OLD(i,j,k)
302 xa % qcw (i,j,k) = xa % qcw (i,j,k) + QCW_OLD(i,j,k)
303 xa % q (i,j,k) = xa % q (i,j,k) + Q_OLD(i,j,k)
304 xa % t (i,j,k) = xa % t (i,j,k) + T_OLD(i,j,k)
305 end do
306 end do
307 end do
308
309 if (trace_use) call da_trace_exit("da_moist_phys_adj")
310
311 end subroutine da_moist_phys_adj
312
313