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