| blob | 246a8134f1ab6d06e2918cd27f449804f8592ecd |
1 !WRF:model_layer:physics
2 !
3 !
4 !
5 !
6 module a_module_bl_gwdo
7 contains
8 !
9 ! Generated by TAPENADE (INRIA, Tropics team)
10 ! Tapenade 3.10 (r5498) - 20 Jan 2015 09:48
11 !
12 ! Differentiation of gwdo in reverse (adjoint) mode:
13 ! gradient of useful results: rublten p3di u3d dusfcg z dvsfcg
14 ! dtauy3d rvblten t3d qv3d pi3d v3d dtaux3d
15 ! with respect to varying inputs: rublten p3di u3d dusfcg z dvsfcg
16 ! dtauy3d rvblten t3d qv3d pi3d v3d dtaux3d
17 ! RW status of diff variables: rublten:in-out p3di:incr u3d:incr
18 ! dusfcg:in-out z:incr dvsfcg:in-out dtauy3d:in-out
19 ! rvblten:in-out t3d:incr qv3d:incr pi3d:incr v3d:incr
20 ! dtaux3d:in-out
21 !-------------------------------------------------------------------------------
22 SUBROUTINE GWDO_B(u3d, u3db, v3d, v3db, t3d, t3db, qv3d, qv3db, p3d, &
23 & p3di, p3dib, pi3d, pi3db, z, zb, rublten, rubltenb, rvblten, rvbltenb&
24 & , dtaux3d, dtaux3db, dtauy3d, dtauy3db, dusfcg, dusfcgb, dvsfcg, &
25 & dvsfcgb, var2d, oc12d, oa2d1, oa2d2, oa2d3, oa2d4, ol2d1, ol2d2, ol2d3&
26 & , ol2d4, znu, znw, p_top, cp, g, rd, rv, ep1, pi, dt, dx, &
27 & kpbl2d, itimestep, ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, &
28 & kms, kme, its, ite, jts, jte, kts, kte)
29 IMPLICIT NONE
30 !
31 !-------------------------------------------------------------------------------
32 !
33 !-- u3d 3d u-velocity interpolated to theta points (m/s)
34 !-- v3d 3d v-velocity interpolated to theta points (m/s)
35 !-- t3d temperature (k)
36 !-- qv3d 3d water vapor mixing ratio (kg/kg)
37 !-- p3d 3d pressure (pa)
38 !-- p3di 3d pressure (pa) at interface level
39 !-- pi3d 3d exner function (dimensionless)
40 !-- rublten u tendency due to pbl parameterization (m/s/s)
41 !-- rvblten v tendency due to pbl parameterization (m/s/s)
42 !-- znu eta values (sigma values)
43 !-- cp heat capacity at constant pressure for dry air (j/kg/k)
44 !-- g acceleration due to gravity (m/s^2)
45 !-- rd gas constant for dry air (j/kg/k)
46 !-- z height above sea level (m)
47 !-- rv gas constant for water vapor (j/kg/k)
48 !-- dt time step (s)
49 !-- dx model grid interval (m)
50 !-- ep1 constant for virtual temperature (r_v/r_d - 1) (dimensionless)
51 !-- ids start index for i in domain
52 !-- ide end index for i in domain
53 !-- jds start index for j in domain
54 !-- jde end index for j in domain
55 !-- kds start index for k in domain
56 !-- kde end index for k in domain
57 !-- ims start index for i in memory
58 !-- ime end index for i in memory
59 !-- jms start index for j in memory
60 !-- jme end index for j in memory
61 !-- kms start index for k in memory
62 !-- kme end index for k in memory
63 !-- its start index for i in tile
64 !-- ite end index for i in tile
65 !-- jts start index for j in tile
66 !-- jte end index for j in tile
67 !-- kts start index for k in tile
68 !-- kte end index for k in tile
69 !
70 !-------------------------------------------------------------------------------
71 INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, ims, ime, jms, &
72 & jme, kms, kme, its, ite, jts, jte, kts, kte
73 INTEGER, INTENT(IN) :: itimestep
74 !
75 REAL, INTENT(IN) :: dt, dx, cp, g, rd, rv, ep1, pi
76 !
77 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: qv3d, p3d, &
78 & pi3d, t3d, z
79 REAL, DIMENSION(ims:ime, kms:kme, jms:jme) :: qv3db, pi3db, t3db, zb
80 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: p3di
81 REAL, DIMENSION(ims:ime, kms:kme, jms:jme) :: p3dib
82 !
83 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: rublten, &
84 & rvblten
85 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: rubltenb
86 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: dtaux3d, &
87 & dtauy3d
88 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: dtaux3db
89 !
90 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: u3d, v3d
91 REAL, DIMENSION(ims:ime, kms:kme, jms:jme) :: u3db, v3db
92 !
93 INTEGER, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: kpbl2d
94 REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: dusfcg, dvsfcg
95 REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: dusfcgb
96 !
97 REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: var2d, oc12d, oa2d1, &
98 & oa2d2, oa2d3, oa2d4, ol2d1, ol2d2, ol2d3, ol2d4
99 !
100 REAL, DIMENSION(kms:kme), OPTIONAL, INTENT(IN) :: znu, znw
101 !
102 REAL, OPTIONAL, INTENT(IN) :: p_top
103 !
104 !local
105 !
106 REAL, DIMENSION(its:ite, kts:kte) :: delprsi, pdh
107 REAL, DIMENSION(its:ite, kts:kte) :: delprsib, pdhb
108 REAL, DIMENSION(its:ite, kts:kte + 1) :: pdhi
109 REAL, DIMENSION(its:ite, kts:kte+1) :: pdhib
110 REAL, DIMENSION(its:ite, 4) :: oa4, ol4
111 INTEGER :: i, j, k, kdt, kpblmax
112 INTEGER :: branch
113 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: rvbltenb
114 REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: dvsfcgb
115 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: dtauy3db
116 !
117 DO k=kts,kte
118 IF (znu(k) .GT. 0.6) kpblmax = k + 1
119 END DO
120 !
121 DO j=jts,jte
122 DO k=kts,kte+1
123 DO i=its,ite
124 IF (k .LE. kte) THEN
125 CALL PUSHREAL8(pdh(i, k))
126 pdh(i, k) = p3d(i, k, j)
127 CALL PUSHCONTROL1B(0)
128 ELSE
129 CALL PUSHCONTROL1B(1)
130 END IF
131 CALL PUSHREAL8(pdhi(i, k))
132 pdhi(i, k) = p3di(i, k, j)
133 END DO
134 END DO
135 !
136 DO k=kts,kte
137 DO i=its,ite
138 CALL PUSHREAL8(delprsi(i, k))
139 delprsi(i, k) = pdhi(i, k) - pdhi(i, k+1)
140 END DO
141 END DO
142 DO i=its,ite
143 CALL PUSHREAL8(oa4(i, 1))
144 oa4(i, 1) = oa2d1(i, j)
145 CALL PUSHREAL8(oa4(i, 2))
146 oa4(i, 2) = oa2d2(i, j)
147 CALL PUSHREAL8(oa4(i, 3))
148 oa4(i, 3) = oa2d3(i, j)
149 CALL PUSHREAL8(oa4(i, 4))
150 oa4(i, 4) = oa2d4(i, j)
151 CALL PUSHREAL8(ol4(i, 1))
152 ol4(i, 1) = ol2d1(i, j)
153 CALL PUSHREAL8(ol4(i, 2))
154 ol4(i, 2) = ol2d2(i, j)
155 CALL PUSHREAL8(ol4(i, 3))
156 ol4(i, 3) = ol2d3(i, j)
157 CALL PUSHREAL8(ol4(i, 4))
158 ol4(i, 4) = ol2d4(i, j)
159 END DO
160 END DO
161 delprsib = 0.0
162 pdhib = 0.0
163 pdhb = 0.0
164 DO j=jte,jts,-1
165 CALL GWDO2D_B(rublten(ims, kms, j), rubltenb(ims, kms, j), rvblten(&
166 & ims, kms, j), rvbltenb(ims, kms, j), dtaux3d(ims, kms, j), &
167 & dtaux3db(ims, kms, j), dtauy3d(ims, kms, j), dtauy3db(ims, &
168 & kms, j), u3d(ims, kms, j), u3db(ims, kms, j), v3d(ims, kms, &
169 & j), v3db(ims, kms, j), t3d(ims, kms, j), t3db(ims, kms, j), &
170 & qv3d(ims, kms, j), qv3db(ims, kms, j), delprsi(its, kts), &
171 & delprsib(its, kts), pdhi(its, kts), pdhib(its, kts), pdh(its&
172 & , kts), pdhb(its, kts), pi3d(ims, kms, j), pi3db(ims, kms, j&
173 & ), z(ims, kms, j), zb(ims, kms, j), rcl=1.0, kpblmax=kpblmax&
174 & , var=var2d(ims, j), oc1=oc12d(ims, j), oa4=oa4, ol4=ol4, &
175 & dusfc=dusfcg(ims, j), dusfcb=dusfcgb(ims, j), dvsfc=dvsfcg(&
176 & ims, j), dvsfcb=dvsfcgb(ims, j), g=g, cp=cp, rd=rd, rv=rv, &
177 & fv=ep1, pi=pi, dxmeter=dx, deltim=dt, kpbl=kpbl2d(ims, j), &
178 & kdt=itimestep, lat=j, ids=ids, ide=ide, jds=jds, jde=jde, &
179 & kds=kds, kde=kde, ims=ims, ime=ime, jms=jms, jme=jme, kms=&
180 & kms, kme=kme, its=its, ite=ite, jts=jts, jte=jte, kts=kts, &
181 & kte=kte)
182 DO i=ite,its,-1
183 CALL POPREAL8(ol4(i, 4))
184 CALL POPREAL8(ol4(i, 3))
185 CALL POPREAL8(ol4(i, 2))
186 CALL POPREAL8(ol4(i, 1))
187 CALL POPREAL8(oa4(i, 4))
188 CALL POPREAL8(oa4(i, 3))
189 CALL POPREAL8(oa4(i, 2))
190 CALL POPREAL8(oa4(i, 1))
191 END DO
192 DO k=kte,kts,-1
193 DO i=ite,its,-1
194 CALL POPREAL8(delprsi(i, k))
195 pdhib(i, k) = pdhib(i, k) + delprsib(i, k)
196 pdhib(i, k+1) = pdhib(i, k+1) - delprsib(i, k)
197 delprsib(i, k) = 0.0
198 END DO
199 END DO
200 DO k=kte+1,kts,-1
201 DO i=ite,its,-1
202 CALL POPREAL8(pdhi(i, k))
203 p3dib(i, k, j) = p3dib(i, k, j) + pdhib(i, k)
204 pdhib(i, k) = 0.0
205 CALL POPCONTROL1B(branch)
206 IF (branch .EQ. 0) THEN
207 CALL POPREAL8(pdh(i, k))
208 pdhb(i, k) = 0.0
209 END IF
210 END DO
211 END DO
212 END DO
213 END SUBROUTINE GWDO_B
215 ! Differentiation of gwdo2d in reverse (adjoint) mode:
216 ! gradient of useful results: v1 dvsfc dvdt dtauy2d prsi
217 ! prsl dusfc del t1 q1 dudt dtaux2d u1 zl prslk
218 ! with respect to varying inputs: v1 dvsfc dvdt dtauy2d prsi
219 ! prsl dusfc del t1 q1 dudt dtaux2d u1 zl prslk
220 !-------------------------------------------------------------------------------
221 !
222 !-------------------------------------------------------------------------------
223 SUBROUTINE GWDO2D_B(dudt, dudtb, dvdt, dvdtb, dtaux2d, dtaux2db, dtauy2d&
224 & , dtauy2db, u1, u1b, v1, v1b, t1, t1b, q1, q1b, del, delb, prsi, prsib&
225 & , prsl, prslb, prslk, prslkb, zl, zlb, rcl, kpblmax, var, oc1, oa4, &
226 & ol4, dusfc, dusfcb, dvsfc, dvsfcb, g, cp, rd, rv, fv, pi, dxmeter, &
227 & deltim, kpbl, kdt, lat, ids, ide, jds, jde, kds, kde, ims, ime, jms, &
228 & jme, kms, kme, its, ite, jts, jte, kts, kte)
229 IMPLICIT NONE
230 !-------------------------------------------------------------------------------
231 INTEGER :: kdt, lat, latd, lond, kpblmax, ids, ide, jds, jde, kds, kde&
232 & , ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts, kte
233 !
234 REAL :: g, rd, rv, fv, cp, pi, dxmeter, deltim, rcl
235 REAL :: dudt(ims:ime, kms:kme), dvdt(ims:ime, kms:kme), dtaux2d(ims:&
236 & ime, kms:kme), dtauy2d(ims:ime, kms:kme), u1(ims:ime, kms:kme), v1(ims&
237 & :ime, kms:kme), t1(ims:ime, kms:kme), q1(ims:ime, kms:kme), zl(ims:ime&
238 & , kms:kme), prsl(its:ite, kts:kte), prslk(ims:ime, kms:kme)
239 REAL :: dudtb(ims:ime, kms:kme), dvdtb(ims:ime, kms:kme), dtaux2db(ims&
240 & :ime, kms:kme), dtauy2db(ims:ime, kms:kme), u1b(ims:ime, kms:kme), v1b&
241 & (ims:ime, kms:kme), t1b(ims:ime, kms:kme), q1b(ims:ime, kms:kme), zlb(&
242 & ims:ime, kms:kme), prslb(its:ite, kts:kte), prslkb(ims:ime, kms:kme)
243 REAL :: prsi(its:ite, kts:kte+1), del(its:ite, kts:kte)
244 REAL :: prsib(its:ite, kts:kte+1), delb(its:ite, kts:kte)
245 REAL :: oa4(its:ite, 4), ol4(its:ite, 4)
246 !
247 INTEGER :: kpbl(ims:ime)
248 REAL :: var(ims:ime), oc1(ims:ime), dusfc(ims:ime), dvsfc(ims:ime)
249 REAL :: dusfcb(ims:ime), dvsfcb(ims:ime)
250 !
251 ! critical richardson number for wave breaking : ! larger drag with larger value
252 !
253 REAL, PARAMETER :: ric=0.25
254 !
255 REAL, PARAMETER :: dw2min=1.
256 REAL, PARAMETER :: rimin=-100.
257 REAL, PARAMETER :: bnv2min=1.0e-5
258 REAL, PARAMETER :: efmin=0.0
259 REAL, PARAMETER :: efmax=10.0
260 REAL, PARAMETER :: xl=4.0e4
261 REAL, PARAMETER :: critac=1.0e-5
262 REAL, PARAMETER :: gmax=1.
263 REAL, PARAMETER :: veleps=1.0
264 REAL, PARAMETER :: factop=0.5
265 REAL, PARAMETER :: frc=1.0
266 REAL, PARAMETER :: ce=0.8
267 REAL, PARAMETER :: cg=0.5
268 INTEGER, PARAMETER :: kpblmin=2
269 !
270 ! local variables
271 !
272 INTEGER :: i, k, lcap, lcapp1, nwd, idir, klcap, kp1, ikount, kk
273 !
274 REAL :: rcs, rclcs, csg, fdir, cleff, cs, rcsks, wdir, ti, rdz, temp, &
275 & tem2, dw2, shr2, bvf2, rdelks, wtkbj, tem, gfobnv, hd, fro, rim, temc&
276 & , tem1, efact, temv, dtaux, dtauy
277 REAL :: rcsksb, tib, rdzb, tem2b, dw2b, shr2b, bvf2b, rdelksb, wtkbjb&
278 & , temb, gfobnvb, hdb, temcb, tem1b, efactb, dtauxb, dtauyb
279 !
280 LOGICAL :: ldrag(its:ite), icrilv(its:ite), flag(its:ite), kloop1(its:&
281 & ite)
282 !
283 REAL :: taub(its:ite), taup(its:ite, kts:kte+1), xn(its:ite), yn(its:&
284 & ite), ubar(its:ite), vbar(its:ite), fr(its:ite), ulow(its:ite), rulow(&
285 & its:ite), bnv(its:ite), oa(its:ite), ol(its:ite), roll(its:ite), dtfac&
286 & (its:ite), brvf(its:ite), xlinv(its:ite), delks(its:ite), delks1(its:&
287 & ite), bnv2(its:ite, kts:kte), usqj(its:ite, kts:kte), taud(its:ite, &
288 & kts:kte), ro(its:ite, kts:kte), vtk(its:ite, kts:kte), vtj(its:ite, &
289 & kts:kte), zlowtop(its:ite), velco(its:ite, kts:kte-1), coefm(its:ite)
290 REAL :: taubb(its:ite), taupb(its:ite, kts:kte+1), xnb(its:ite), ynb(&
291 & its:ite), ubarb(its:ite), vbarb(its:ite), frb(its:ite), ulowb(its:ite)&
292 & , rulowb(its:ite), bnvb(its:ite), rollb(its:ite), dtfacb(its:ite), &
293 & brvfb(its:ite), delksb(its:ite), delks1b(its:ite), bnv2b(its:ite, kts:&
294 & kte), usqjb(its:ite, kts:kte), taudb(its:ite, kts:kte), rob(its:ite, &
295 & kts:kte), vtkb(its:ite, kts:kte), vtjb(its:ite, kts:kte), velcob(its:&
296 & ite, kts:kte-1)
297 !
298 INTEGER :: kbl(its:ite), klowtop(its:ite)
299 !
300 LOGICAL :: iope
301 INTEGER, PARAMETER :: mdir=8
302 INTEGER :: nwdir(mdir)
303 !
304 ! variables for flow-blocking drag
305 !
306 REAL, PARAMETER :: frmax=10.
307 REAL, PARAMETER :: olmin=1.0e-5
308 REAL, PARAMETER :: odmin=0.1
309 REAL, PARAMETER :: odmax=10.
310 REAL, PARAMETER :: erad=6371.315e+3
311 INTEGER :: komax(its:ite)
312 INTEGER :: kblk
313 REAL :: cd
314 REAL :: zblk, tautem
315 REAL :: zblkb, tautemb
316 REAL :: pe, ke
317 REAL :: delx, dely, dxy4(4), dxy4p(4)
318 REAL :: dxy(its:ite), dxyp(its:ite)
319 REAL :: ol4p(4), olp(its:ite), od(its:ite)
320 REAL :: taufb(its:ite, kts:kte+1)
321 REAL :: taufbb(its:ite, kts:kte+1)
322 REAL :: tmp
323 REAL :: tmp0
324 REAL :: tmp1
325 REAL :: tmp2
326 REAL :: tmp3
327 INTEGER :: branch
328 INTEGER :: ad_to
329 REAL :: temp3
330 REAL :: x3b
331 REAL :: temp2
332 REAL :: y1b
333 REAL :: temp1
334 REAL :: temp0
335 REAL :: tempb9
336 REAL :: tempb8
337 REAL :: tempb7
338 REAL :: max2b
339 REAL :: tempb6
340 REAL :: tempb5
341 REAL :: tempb4
342 REAL :: tempb19
343 REAL :: tempb3
344 REAL :: tempb18
345 REAL :: tempb2
346 REAL :: tempb17
347 REAL :: tempb1
348 REAL :: tempb16
349 REAL :: tempb0
350 REAL :: tempb15
351 REAL :: tempb14
352 REAL :: tempb13
353 REAL :: tmpb3
354 REAL :: tmpb
355 REAL :: tempb12
356 REAL :: tmpb2
357 REAL :: tmpb1
358 REAL :: tempb11
359 REAL :: tmpb0
360 REAL :: tempb10
361 REAL :: x4
362 REAL :: x3
363 REAL :: x2
364 INTEGER :: x1
365 REAL :: x2b
366 REAL :: temp12
367 REAL :: temp11
368 REAL :: temp10
369 REAL :: tempb
370 REAL :: x4b
371 REAL :: max3
372 REAL :: max2
373 REAL :: max1
374 REAL :: temp9
375 REAL :: temp8
376 REAL :: tempb20
377 REAL :: temp7
378 REAL :: temp6
379 REAL :: temp5
380 REAL :: y1
381 REAL :: temp4
382 DATA nwdir /6, 7, 5, 8, 2, 3, 1, 4/
383 !
384 !---- constants
385 !
386 rcs = SQRT(rcl)
387 cs = 1./SQRT(rcl)
388 csg = cs*g
389 lcap = kte
390 lcapp1 = lcap + 1
391 fdir = mdir/(2.0*pi)
392 !
393 !--- calculate length of grid for flow-blocking drag
394 !
395 delx = dxmeter
396 dely = dxmeter
397 dxy4(1) = delx
398 dxy4(2) = dely
399 dxy4(3) = SQRT(delx*delx + dely*dely)
400 dxy4(4) = dxy4(3)
401 dxy4p(1) = dxy4(2)
402 dxy4p(2) = dxy4(1)
403 dxy4p(3) = dxy4(4)
404 dxy4p(4) = dxy4(3)
405 !
406 !
407 !-----initialize arrays
408 !
409 DO i=its,ite
410 klowtop(i) = 0
411 kbl(i) = 0
412 END DO
413 !
414 DO i=its,ite
415 xn(i) = 0.0
416 yn(i) = 0.0
417 ubar(i) = 0.0
418 vbar(i) = 0.0
419 roll(i) = 0.0
420 taub(i) = 0.0
421 taup(i, 1) = 0.0
422 oa(i) = 0.0
423 ol(i) = 0.0
424 ulow(i) = 0.0
425 dtfac(i) = 1.0
426 ldrag(i) = .false.
427 icrilv(i) = .false.
428 END DO
429 !
430 DO k=kts,kte
431 DO i=its,ite
432 usqj(i, k) = 0.0
433 bnv2(i, k) = 0.0
434 vtj(i, k) = 0.0
435 vtk(i, k) = 0.0
436 taup(i, k) = 0.0
437 taud(i, k) = 0.0
438 END DO
439 END DO
440 !
441 DO i=its,ite
442 taup(i, kte+1) = 0.0
443 END DO
444 !
445 ! initialize array for flow-blocking drag
446 !
447 taufb(its:ite, kts:kte+1) = 0.0
448 !
449 DO k=kts,kte
450 DO i=its,ite
451 vtj(i, k) = t1(i, k)*(1.+fv*q1(i, k))
452 vtk(i, k) = vtj(i, k)/prslk(i, k)
453 ! density kg/m**3
454 ro(i, k) = 1./rd*prsl(i, k)/vtj(i, k)
455 END DO
456 END DO
457 !
458 ! determine reference level: maximum of 2*var and pbl heights
459 !
460 DO i=its,ite
461 zlowtop(i) = 2.*var(i)
462 END DO
463 !
464 DO i=its,ite
465 kloop1(i) = .true.
466 END DO
467 !
468 DO k=kts+1,kte
469 DO i=its,ite
470 IF (kloop1(i) .AND. zl(i, k) - zl(i, 1) .GE. zlowtop(i)) THEN
471 klowtop(i) = k + 1
472 kloop1(i) = .false.
473 END IF
474 END DO
475 END DO
476 !
477 DO i=its,ite
478 IF (kpbl(i) .LT. klowtop(i)) THEN
479 kbl(i) = klowtop(i)
480 ELSE
481 kbl(i) = kpbl(i)
482 END IF
483 IF (kbl(i) .GT. kpblmax) THEN
484 x1 = kpblmax
485 ELSE
486 x1 = kbl(i)
487 END IF
488 IF (x1 .LT. kpblmin) THEN
489 kbl(i) = kpblmin
490 ELSE
491 kbl(i) = x1
492 END IF
493 END DO
494 !
495 ! determine the level of maximum orographic height
496 !
497 komax(:) = kbl(:)
498 !
499 DO i=its,ite
500 delks(i) = 1.0/(prsi(i, 1)-prsi(i, kbl(i)))
501 delks1(i) = 1.0/(prsl(i, 1)-prsl(i, kbl(i)))
502 END DO
503 !
504 ! compute low level averages within pbl
505 !
506 DO k=kts,kpblmax
507 DO i=its,ite
508 IF (k .LT. kbl(i)) THEN
509 rcsks = rcs*del(i, k)*delks(i)
510 rdelks = del(i, k)*delks(i)
511 ! pbl u mean
512 ubar(i) = ubar(i) + rcsks*u1(i, k)
513 ! pbl v mean
514 vbar(i) = vbar(i) + rcsks*v1(i, k)
515 ! ro mean
516 roll(i) = roll(i) + rdelks*ro(i, k)
517 CALL PUSHCONTROL1B(1)
518 ELSE
519 CALL PUSHCONTROL1B(0)
520 END IF
521 END DO
522 END DO
523 !
524 ! figure out low-level horizontal wind direction
525 !
526 ! nwd 1 2 3 4 5 6 7 8
527 ! wd w s sw nw e n ne se
528 !
529 DO i=its,ite
530 wdir = ATAN2(ubar(i), vbar(i)) + pi
531 idir = MOD(NINT(fdir*wdir), mdir) + 1
532 nwd = nwdir(idir)
533 oa(i) = (1-2*INT((nwd-1)/4))*oa4(i, MOD(nwd-1, 4)+1)
534 ol(i) = ol4(i, MOD(nwd-1, 4)+1)
535 !
536 !----- compute orographic width along (ol) and perpendicular (olp)
537 !----- the direction of wind
538 !
539 ol4p(1) = ol4(i, 2)
540 ol4p(2) = ol4(i, 1)
541 ol4p(3) = ol4(i, 4)
542 ol4p(4) = ol4(i, 3)
543 olp(i) = ol4p(MOD(nwd-1, 4)+1)
544 IF (ol(i) .LT. olmin) THEN
545 max1 = olmin
546 ELSE
547 max1 = ol(i)
548 END IF
549 !
550 !----- compute orographic direction (horizontal orographic aspect ratio)
551 !
552 od(i) = olp(i)/max1
553 IF (od(i) .GT. odmax) THEN
554 od(i) = odmax
555 ELSE
556 od(i) = od(i)
557 END IF
558 IF (od(i) .LT. odmin) THEN
559 od(i) = odmin
560 ELSE
561 od(i) = od(i)
562 END IF
563 !
564 !----- compute length of grid in the along(dxy) and cross(dxyp) wind directions
565 !
566 dxy(i) = dxy4(MOD(nwd-1, 4)+1)
567 dxyp(i) = dxy4p(MOD(nwd-1, 4)+1)
568 END DO
569 !
570 !--- saving richardson number in usqj for migwdi
571 !
572 DO k=kts,kte-1
573 DO i=its,ite
574 ti = 2.0/(t1(i, k)+t1(i, k+1))
575 rdz = 1./(zl(i, k+1)-zl(i, k))
576 CALL PUSHREAL8(tem1)
577 tem1 = u1(i, k) - u1(i, k+1)
578 CALL PUSHREAL8(tem2)
579 tem2 = v1(i, k) - v1(i, k+1)
580 dw2 = rcl*(tem1*tem1+tem2*tem2)
581 IF (dw2 .LT. dw2min) THEN
582 CALL PUSHREAL8(max2)
583 max2 = dw2min
584 CALL PUSHCONTROL1B(0)
585 ELSE
586 CALL PUSHREAL8(max2)
587 max2 = dw2
588 CALL PUSHCONTROL1B(1)
589 END IF
590 shr2 = max2*rdz*rdz
591 CALL PUSHREAL8(bvf2)
592 bvf2 = g*(g/cp+rdz*(vtj(i, k+1)-vtj(i, k)))*ti
593 IF (bvf2/shr2 .LT. rimin) THEN
594 usqj(i, k) = rimin
595 CALL PUSHCONTROL1B(0)
596 ELSE
597 usqj(i, k) = bvf2/shr2
598 CALL PUSHCONTROL1B(1)
599 END IF
600 bnv2(i, k) = 2.0*g*rdz*(vtk(i, k+1)-vtk(i, k))/(vtk(i, k+1)+vtk(i&
601 & , k))
602 IF (bnv2(i, k) .LT. bnv2min) THEN
603 bnv2(i, k) = bnv2min
604 CALL PUSHCONTROL1B(0)
605 ELSE
606 CALL PUSHCONTROL1B(1)
607 bnv2(i, k) = bnv2(i, k)
608 END IF
609 END DO
610 END DO
611 !
612 !----compute the "low level" or 1/3 wind magnitude (m/s)
613 !
614 DO i=its,ite
615 x2 = SQRT(ubar(i)*ubar(i) + vbar(i)*vbar(i))
616 IF (x2 .LT. 1.0) THEN
617 ulow(i) = 1.0
618 CALL PUSHCONTROL1B(0)
619 ELSE
620 ulow(i) = x2
621 CALL PUSHCONTROL1B(1)
622 END IF
623 rulow(i) = 1./ulow(i)
624 END DO
625 !
626 DO k=kts,kte-1
627 DO i=its,ite
628 velco(i, k) = 0.5*rcs*((u1(i, k)+u1(i, k+1))*ubar(i)+(v1(i, k)+v1(&
629 & i, k+1))*vbar(i))
630 CALL PUSHREAL8(velco(i, k))
631 velco(i, k) = velco(i, k)*rulow(i)
632 IF (velco(i, k) .LT. veleps .AND. velco(i, k) .GT. 0.) THEN
633 velco(i, k) = veleps
634 CALL PUSHCONTROL1B(1)
635 ELSE
636 CALL PUSHCONTROL1B(0)
637 END IF
638 END DO
639 END DO
640 !
641 ! no drag when critical level in the base layer
642 !
643 DO i=its,ite
644 ldrag(i) = velco(i, 1) .LE. 0.
645 END DO
646 !
647 ! no drag when velco.lt.0
648 !
649 DO k=kpblmin,kpblmax
650 DO i=its,ite
651 IF (k .LT. kbl(i)) ldrag(i) = ldrag(i) .OR. velco(i, k) .LE. 0.
652 END DO
653 END DO
654 !
655 ! no drag when bnv2.lt.0
656 !
657 DO k=kts,kpblmax
658 DO i=its,ite
659 IF (k .LT. kbl(i)) ldrag(i) = ldrag(i) .OR. bnv2(i, k) .LT. 0.
660 END DO
661 END DO
662 !
663 !-----the low level weighted average ri is stored in usqj(1,1; im)
664 !-----the low level weighted average n**2 is stored in bnv2(1,1; im)
665 !---- this is called bnvl2 in phys_gwd_alpert_sub not bnv2
666 !---- rdelks (del(k)/delks) vert ave factor so we can * instead of /
667 !
668 DO i=its,ite
669 wtkbj = (prsl(i, 1)-prsl(i, 2))*delks1(i)
670 CALL PUSHREAL8(bnv2(i, 1))
671 bnv2(i, 1) = wtkbj*bnv2(i, 1)
672 CALL PUSHREAL8(usqj(i, 1))
673 usqj(i, 1) = wtkbj*usqj(i, 1)
674 END DO
675 !
676 DO k=kpblmin,kpblmax
677 DO i=its,ite
678 IF (k .LT. kbl(i)) THEN
679 rdelks = (prsl(i, k)-prsl(i, k+1))*delks1(i)
680 tmp = bnv2(i, 1) + bnv2(i, k)*rdelks
681 CALL PUSHREAL8(bnv2(i, 1))
682 bnv2(i, 1) = tmp
683 tmp0 = usqj(i, 1) + usqj(i, k)*rdelks
684 CALL PUSHREAL8(usqj(i, 1))
685 usqj(i, 1) = tmp0
686 CALL PUSHCONTROL1B(1)
687 ELSE
688 CALL PUSHCONTROL1B(0)
689 END IF
690 END DO
691 END DO
692 !
693 DO i=its,ite
694 ldrag(i) = ldrag(i) .OR. bnv2(i, 1) .LE. 0.0
695 ldrag(i) = ldrag(i) .OR. ulow(i) .EQ. 1.0
696 ldrag(i) = ldrag(i) .OR. var(i) .LE. 0.0
697 END DO
698 !
699 ! set all ri low level values to the low level value
700 !
701 DO k=kpblmin,kpblmax
702 DO i=its,ite
703 IF (k .LT. kbl(i)) THEN
704 tmp1 = usqj(i, 1)
705 CALL PUSHREAL8(usqj(i, k))
706 usqj(i, k) = tmp1
707 CALL PUSHCONTROL1B(1)
708 ELSE
709 CALL PUSHCONTROL1B(0)
710 END IF
711 END DO
712 END DO
713 !
714 DO i=its,ite
715 IF (.NOT.ldrag(i)) THEN
716 bnv(i) = SQRT(bnv2(i, 1))
717 fr(i) = bnv(i)*rulow(i)*2.*var(i)*od(i)
718 IF (fr(i) .GT. frmax) THEN
719 fr(i) = frmax
720 CALL PUSHCONTROL1B(0)
721 ELSE
722 CALL PUSHCONTROL1B(1)
723 fr(i) = fr(i)
724 END IF
725 xn(i) = ubar(i)*rulow(i)
726 yn(i) = vbar(i)*rulow(i)
727 CALL PUSHCONTROL1B(1)
728 ELSE
729 CALL PUSHCONTROL1B(0)
730 END IF
731 END DO
732 !
733 ! compute the base level stress and store it in taub
734 ! calculate enhancement factor, number of mountains & aspect
735 ! ratio const. use simplified relationship between standard
736 ! deviation & critical hgt
737 !
738 DO i=its,ite
739 IF (.NOT.ldrag(i)) THEN
740 CALL PUSHREAL8(efact)
741 efact = (oa(i)+2.)**(ce*fr(i)/frc)
742 IF (efact .LT. efmin) THEN
743 x3 = efmin
744 CALL PUSHCONTROL1B(0)
745 ELSE
746 x3 = efact
747 CALL PUSHCONTROL1B(1)
748 END IF
749 IF (x3 .GT. efmax) THEN
750 efact = efmax
751 CALL PUSHCONTROL1B(0)
752 ELSE
753 efact = x3
754 CALL PUSHCONTROL1B(1)
755 END IF
756 !!!!!!! cleff (effective grid length) is highly tunable parameter
757 !!!!!!! the bigger (smaller) value produce weaker (stronger) wave drag
758 cleff = SQRT(dxy(i)**2. + dxyp(i)**2.)
759 IF (dxmeter .LT. cleff) THEN
760 max3 = cleff
761 ELSE
762 max3 = dxmeter
763 END IF
764 cleff = 3.*max3
765 coefm(i) = (1.+ol(i))**(oa(i)+1.)
766 xlinv(i) = coefm(i)/cleff
767 tem = fr(i)*fr(i)*oc1(i)
768 gfobnv = gmax*tem/((tem+cg)*bnv(i))
769 taub(i) = xlinv(i)*roll(i)*ulow(i)*ulow(i)*ulow(i)*gfobnv*efact
770 CALL PUSHCONTROL1B(1)
771 ELSE
772 taub(i) = 0.0
773 xn(i) = 0.0
774 yn(i) = 0.0
775 CALL PUSHCONTROL1B(0)
776 END IF
777 END DO
778 !
779 ! now compute vertical structure of the stress.
780 !
781 DO k=kts,kpblmax
782 DO i=its,ite
783 IF (k .LE. kbl(i)) THEN
784 taup(i, k) = taub(i)
785 CALL PUSHCONTROL1B(1)
786 ELSE
787 CALL PUSHCONTROL1B(0)
788 END IF
789 END DO
790 END DO
791 !
792 ! vertical level k loop!
793 DO k=kpblmin,kte-1
794 CALL PUSHINTEGER4(kp1)
795 kp1 = k + 1
796 DO i=its,ite
797 !
798 ! unstablelayer if ri < ric
799 ! unstable layer if upper air vel comp along surf vel <=0 (crit lay)
800 ! at (u-c)=0. crit layer exists and bit vector should be set (.le.)
801 !
802 IF (k .GE. kbl(i)) THEN
803 icrilv(i) = (icrilv(i) .OR. usqj(i, k) .LT. ric) .OR. velco(i, k&
804 & ) .LE. 0.0
805 IF (bnv2(i, k) .LT. bnv2min) THEN
806 CALL PUSHREAL8(brvf(i))
807 brvf(i) = bnv2min
808 CALL PUSHCONTROL1B(0)
809 ELSE
810 CALL PUSHREAL8(brvf(i))
811 brvf(i) = bnv2(i, k)
812 CALL PUSHCONTROL1B(1)
813 END IF
814 ! brunt-vaisala frequency
815 CALL PUSHREAL8(brvf(i))
816 brvf(i) = SQRT(brvf(i))
817 CALL PUSHCONTROL1B(1)
818 ELSE
819 CALL PUSHCONTROL1B(0)
820 END IF
821 END DO
822 !
823 DO i=its,ite
824 IF (k .GE. kbl(i) .AND. (.NOT.ldrag(i))) THEN
825 IF (.NOT.icrilv(i) .AND. taup(i, k) .GT. 0.0) THEN
826 temv = 1.0/velco(i, k)
827 CALL PUSHREAL8(tem1)
828 tem1 = coefm(i)/dxy(i)*(ro(i, kp1)+ro(i, k))*brvf(i)*velco(i, &
829 & k)*0.5
830 CALL PUSHREAL8(hd)
831 hd = SQRT(taup(i, k)/tem1)
832 fro = brvf(i)*hd*temv
833 !
834 ! rim is the minimum-richardson number by shutts (1985)
835 !
836 CALL PUSHREAL8(tem2)
837 tem2 = SQRT(usqj(i, k))
838 tem = 1. + tem2*fro
839 rim = usqj(i, k)*(1.-fro)/(tem*tem)
840 !
841 ! check stability to employ the 'saturation hypothesis'
842 ! of lindzen (1981) except at tropospheric downstream regions
843 !
844 IF (rim .LE. ric) THEN
845 ! saturation hypothesis!
846 IF (oa(i) .LE. 0. .OR. kp1 .GE. kpblmin) THEN
847 temc = 2.0 + 1.0/tem2
848 hd = velco(i, k)*(2.*SQRT(temc)-temc)/brvf(i)
849 taup(i, kp1) = tem1*hd*hd
850 CALL PUSHCONTROL3B(4)
851 ELSE
852 CALL PUSHCONTROL3B(3)
853 END IF
854 ELSE
855 ! no wavebreaking!
856 tmp2 = taup(i, k)
857 taup(i, kp1) = tmp2
858 CALL PUSHCONTROL3B(2)
859 END IF
860 ELSE
861 CALL PUSHCONTROL3B(1)
862 END IF
863 ELSE
864 CALL PUSHCONTROL3B(0)
865 END IF
866 END DO
867 END DO
868 !
869 IF (lcap .LT. kte) THEN
870 DO klcap=lcapp1,kte
871 DO i=its,ite
872 tmp3 = prsi(i, klcap)/prsi(i, lcap)*taup(i, lcap)
873 CALL PUSHREAL8(taup(i, klcap))
874 taup(i, klcap) = tmp3
875 END DO
876 END DO
877 CALL PUSHCONTROL1B(1)
878 ELSE
879 CALL PUSHCONTROL1B(0)
880 END IF
881 DO i=its,ite
882 IF (.NOT.ldrag(i)) THEN
883 !
884 !------- determine the height of flow-blocking layer
885 !
886 CALL PUSHINTEGER4(kblk)
887 kblk = 0
888 pe = 0.0
889 DO k=kte,kpblmin,-1
890 IF (kblk .EQ. 0 .AND. k .LE. komax(i)) THEN
891 pe = pe + bnv2(i, k)*(zl(i, komax(i))-zl(i, k))*del(i, k)/g/ro&
892 & (i, k)
893 ke = 0.5*((rcs*u1(i, k))**2.+(rcs*v1(i, k))**2.)
894 !
895 !---------- apply flow-blocking drag when pe >= ke
896 !
897 IF (pe .GE. ke) THEN
898 CALL PUSHINTEGER4(kblk)
899 kblk = k
900 IF (kblk .GT. kbl(i)) THEN
901 kblk = kbl(i)
902 ELSE
903 kblk = kblk
904 END IF
905 CALL PUSHREAL8(zblk)
906 zblk = zl(i, kblk) - zl(i, kts)
907 CALL PUSHCONTROL2B(2)
908 ELSE
909 CALL PUSHCONTROL2B(1)
910 END IF
911 ELSE
912 CALL PUSHCONTROL2B(0)
913 END IF
914 END DO
915 IF (kblk .NE. 0) THEN
916 !
917 !--------- compute flow-blocking stress
918 !
919 IF (2.0 - 1.0/od(i) .LT. 0.0) THEN
920 CALL PUSHREAL8(cd)
921 cd = 0.0
922 CALL PUSHCONTROL1B(0)
923 ELSE
924 CALL PUSHREAL8(cd)
925 cd = 2.0 - 1.0/od(i)
926 CALL PUSHCONTROL1B(1)
927 END IF
928 taufb(i, kts) = 0.5*roll(i)*coefm(i)/dxy(i)**2*cd*dxyp(i)*olp(i)&
929 & *zblk*ulow(i)**2
930 tautem = taufb(i, kts)/FLOAT(kblk-kts)
931 DO k=kts+1,kblk
932 taufb(i, k) = taufb(i, k-1) - tautem
933 END DO
934 CALL PUSHINTEGER4(k - 1)
935 !
936 !----------sum orographic GW stress and flow-blocking stress
937 !
938 CALL PUSHREAL8ARRAY(taup(i, :), kte - kts + 2)
939 taup(i, :) = taup(i, :) + taufb(i, :)
940 CALL PUSHCONTROL2B(2)
941 ELSE
942 CALL PUSHCONTROL2B(1)
943 END IF
944 ELSE
945 CALL PUSHCONTROL2B(0)
946 END IF
947 END DO
948 !
949 ! calculate - (g)*d(tau)/d(pressure) and deceleration terms dtaux, dtauy
950 !
951 DO k=kts,kte
952 DO i=its,ite
953 taud(i, k) = 1.*(taup(i, k+1)-taup(i, k))*csg/del(i, k)
954 END DO
955 END DO
956 !
957 ! limit de-acceleration (momentum deposition ) at top to 1/2 value
958 ! the idea is some stuff must go out the 'top'
959 !
960 DO klcap=lcap,kte
961 DO i=its,ite
962 taud(i, klcap) = taud(i, klcap)*factop
963 END DO
964 END DO
965 !
966 ! if the gravity wave drag would force a critical line
967 ! in the lower ksmm1 layers during the next deltim timestep,
968 ! then only apply drag until that critical line is reached.
969 !
970 DO k=kts,kpblmax-1
971 DO i=its,ite
972 IF (k .LE. kbl(i)) THEN
973 IF (taud(i, k) .NE. 0.) THEN
974 x4 = velco(i, k)/(deltim*rcs*taud(i, k))
975 IF (x4 .GE. 0.) THEN
976 y1 = x4
977 CALL PUSHCONTROL1B(0)
978 ELSE
979 y1 = -x4
980 CALL PUSHCONTROL1B(1)
981 END IF
982 IF (dtfac(i) .GT. y1) THEN
983 dtfac(i) = y1
984 CALL PUSHCONTROL2B(2)
985 ELSE
986 dtfac(i) = dtfac(i)
987 CALL PUSHCONTROL2B(3)
988 END IF
989 ELSE
990 CALL PUSHCONTROL2B(1)
991 END IF
992 ELSE
993 CALL PUSHCONTROL2B(0)
994 END IF
995 END DO
996 END DO
997 !
998 DO k=kts,kte
999 DO i=its,ite
1000 CALL PUSHREAL8(taud(i, k))
1001 taud(i, k) = taud(i, k)*dtfac(i)
1002 END DO
1003 END DO
1004 DO i=ite,its,-1
1005 dvsfcb(i) = -(rcs*dvsfcb(i)/g)
1006 dusfcb(i) = -(rcs*dusfcb(i)/g)
1007 END DO
1008 taudb = 0.0
1009 xnb = 0.0
1010 ynb = 0.0
1011 dtfacb = 0.0
1012 DO k=kte,kts,-1
1013 DO i=ite,its,-1
1014 dtaux = taud(i, k)*xn(i)
1015 dtauy = taud(i, k)*yn(i)
1016 dtauyb = dvdtb(i, k) + dtauy2db(i, k) + del(i, k)*dvsfcb(i)
1017 delb(i, k) = delb(i, k) + dtaux*dusfcb(i) + dtauy*dvsfcb(i)
1018 dtauxb = dudtb(i, k) + dtaux2db(i, k) + del(i, k)*dusfcb(i)
1019 dtauy2db(i, k) = 0.0
1020 dtaux2db(i, k) = 0.0
1021 taudb(i, k) = taudb(i, k) + xn(i)*dtauxb + yn(i)*dtauyb
1022 ynb(i) = ynb(i) + taud(i, k)*dtauyb
1023 xnb(i) = xnb(i) + taud(i, k)*dtauxb
1024 CALL POPREAL8(taud(i, k))
1025 dtfacb(i) = dtfacb(i) + taud(i, k)*taudb(i, k)
1026 taudb(i, k) = dtfac(i)*taudb(i, k)
1027 END DO
1028 END DO
1029 DO i=ite,its,-1
1030 dvsfcb(i) = 0.0
1031 dusfcb(i) = 0.0
1032 END DO
1033 velcob = 0.0
1034 DO k=kpblmax-1,kts,-1
1035 DO i=ite,its,-1
1036 CALL POPCONTROL2B(branch)
1037 IF (branch .GE. 2) THEN
1038 IF (branch .EQ. 2) THEN
1039 y1b = dtfacb(i)
1040 dtfacb(i) = 0.0
1041 ELSE
1042 y1b = 0.0
1043 END IF
1044 CALL POPCONTROL1B(branch)
1045 IF (branch .EQ. 0) THEN
1046 x4b = y1b
1047 ELSE
1048 x4b = -y1b
1049 END IF
1050 temp12 = deltim*rcs*taud(i, k)
1051 velcob(i, k) = velcob(i, k) + x4b/temp12
1052 taudb(i, k) = taudb(i, k) - velco(i, k)*deltim*rcs*x4b/temp12**2
1053 END IF
1054 END DO
1055 END DO
1056 DO klcap=kte,lcap,-1
1057 DO i=ite,its,-1
1058 taudb(i, klcap) = factop*taudb(i, klcap)
1059 END DO
1060 END DO
1061 taupb = 0.0
1062 DO k=kte,kts,-1
1063 DO i=ite,its,-1
1064 tempb20 = csg*taudb(i, k)/del(i, k)
1065 taupb(i, k+1) = taupb(i, k+1) + tempb20
1066 taupb(i, k) = taupb(i, k) - tempb20
1067 delb(i, k) = delb(i, k) - (taup(i, k+1)-taup(i, k))*tempb20/del(i&
1068 & , k)
1069 taudb(i, k) = 0.0
1070 END DO
1071 END DO
1072 rollb = 0.0
1073 zblkb = 0.0
1074 taufbb = 0.0
1075 ulowb = 0.0
1076 DO i=ite,its,-1
1077 CALL POPCONTROL2B(branch)
1078 IF (branch .NE. 0) THEN
1079 IF (branch .NE. 1) THEN
1080 CALL POPREAL8ARRAY(taup(i, :), kte - kts + 2)
1081 taufbb(i, :) = taufbb(i, :) + taupb(i, :)
1082 tautemb = 0.0
1083 CALL POPINTEGER4(ad_to)
1084 DO k=ad_to,kts+1,-1
1085 taufbb(i, k-1) = taufbb(i, k-1) + taufbb(i, k)
1086 tautemb = tautemb - taufbb(i, k)
1087 taufbb(i, k) = 0.0
1088 END DO
1089 taufbb(i, kts) = taufbb(i, kts) + tautemb/FLOAT(kblk-kts)
1090 tempb18 = cd*0.5*olp(i)*coefm(i)*dxyp(i)*taufbb(i, kts)
1091 tempb19 = ulow(i)**2*tempb18/dxy(i)**2
1092 rollb(i) = rollb(i) + zblk*tempb19
1093 zblkb = zblkb + roll(i)*tempb19
1094 ulowb(i) = ulowb(i) + roll(i)*zblk*2*ulow(i)*tempb18/dxy(i)**2
1095 taufbb(i, kts) = 0.0
1096 CALL POPCONTROL1B(branch)
1097 IF (branch .EQ. 0) THEN
1098 CALL POPREAL8(cd)
1099 ELSE
1100 CALL POPREAL8(cd)
1101 END IF
1102 END IF
1103 DO k=kpblmin,kte,1
1104 CALL POPCONTROL2B(branch)
1105 IF (branch .NE. 0) THEN
1106 IF (branch .NE. 1) THEN
1107 CALL POPREAL8(zblk)
1108 zlb(i, kblk) = zlb(i, kblk) + zblkb
1109 zlb(i, kts) = zlb(i, kts) - zblkb
1110 CALL POPINTEGER4(kblk)
1111 zblkb = 0.0
1112 END IF
1113 END IF
1114 END DO
1115 CALL POPINTEGER4(kblk)
1116 END IF
1117 END DO
1118 CALL POPCONTROL1B(branch)
1119 IF (branch .NE. 0) THEN
1120 DO klcap=kte,lcapp1,-1
1121 DO i=ite,its,-1
1122 CALL POPREAL8(taup(i, klcap))
1123 tmpb3 = taupb(i, klcap)
1124 taupb(i, klcap) = 0.0
1125 tempb17 = tmpb3/prsi(i, lcap)
1126 prsib(i, klcap) = prsib(i, klcap) + taup(i, lcap)*tempb17
1127 taupb(i, lcap) = taupb(i, lcap) + prsi(i, klcap)*tempb17
1128 prsib(i, lcap) = prsib(i, lcap) - prsi(i, klcap)*taup(i, lcap)*&
1129 & tempb17/prsi(i, lcap)
1130 END DO
1131 END DO
1132 END IF
1133 brvfb = 0.0
1134 bnv2b = 0.0
1135 rob = 0.0
1136 usqjb = 0.0
1137 DO k=kte-1,kpblmin,-1
1138 DO i=ite,its,-1
1139 CALL POPCONTROL3B(branch)
1140 IF (branch .GE. 2) THEN
1141 IF (branch .EQ. 2) THEN
1142 tmpb2 = taupb(i, kp1)
1143 taupb(i, kp1) = 0.0
1144 taupb(i, k) = taupb(i, k) + tmpb2
1145 tem1b = 0.0
1146 tem2b = 0.0
1147 ELSE IF (branch .EQ. 3) THEN
1148 tem1b = 0.0
1149 tem2b = 0.0
1150 ELSE
1151 tem1b = hd**2*taupb(i, kp1)
1152 hdb = tem1*2*hd*taupb(i, kp1)
1153 taupb(i, kp1) = 0.0
1154 temc = 2.0 + 1.0/tem2
1155 temp11 = SQRT(temc)
1156 tempb16 = (2.*temp11-temc)*hdb/brvf(i)
1157 temp10 = velco(i, k)/brvf(i)
1158 velcob(i, k) = velcob(i, k) + tempb16
1159 brvfb(i) = brvfb(i) - temp10*tempb16
1160 IF (temc .EQ. 0.0) THEN
1161 temcb = -(temp10*hdb)
1162 ELSE
1163 temcb = (2.*temp10/(2.0*temp11)-temp10)*hdb
1164 END IF
1165 tem2b = -(temcb/tem2**2)
1166 END IF
1167 CALL POPREAL8(tem2)
1168 IF (.NOT.usqj(i, k) .EQ. 0.0) usqjb(i, k) = usqjb(i, k) + tem2b/&
1169 & (2.0*SQRT(usqj(i, k)))
1170 CALL POPREAL8(hd)
1171 CALL POPREAL8(tem1)
1172 temp9 = brvf(i)/dxy(i)
1173 tempb13 = coefm(i)*0.5*tem1b
1174 tempb14 = velco(i, k)*temp9*tempb13
1175 tempb15 = (ro(i, kp1)+ro(i, k))*tempb13
1176 rob(i, kp1) = rob(i, kp1) + tempb14
1177 rob(i, k) = rob(i, k) + tempb14
1178 velcob(i, k) = velcob(i, k) + temp9*tempb15
1179 brvfb(i) = brvfb(i) + velco(i, k)*tempb15/dxy(i)
1180 END IF
1181 END DO
1182 DO i=ite,its,-1
1183 CALL POPCONTROL1B(branch)
1184 IF (branch .NE. 0) THEN
1185 CALL POPREAL8(brvf(i))
1186 IF (brvf(i) .EQ. 0.0) THEN
1187 brvfb(i) = 0.0
1188 ELSE
1189 brvfb(i) = brvfb(i)/(2.0*SQRT(brvf(i)))
1190 END IF
1191 CALL POPCONTROL1B(branch)
1192 IF (branch .EQ. 0) THEN
1193 CALL POPREAL8(brvf(i))
1194 brvfb(i) = 0.0
1195 ELSE
1196 CALL POPREAL8(brvf(i))
1197 bnv2b(i, k) = bnv2b(i, k) + brvfb(i)
1198 brvfb(i) = 0.0
1199 END IF
1200 END IF
1201 END DO
1202 CALL POPINTEGER4(kp1)
1203 END DO
1204 taubb = 0.0
1205 DO k=kpblmax,kts,-1
1206 DO i=ite,its,-1
1207 CALL POPCONTROL1B(branch)
1208 IF (branch .NE. 0) THEN
1209 taubb(i) = taubb(i) + taupb(i, k)
1210 taupb(i, k) = 0.0
1211 END IF
1212 END DO
1213 END DO
1214 bnvb = 0.0
1215 frb = 0.0
1216 DO i=ite,its,-1
1217 CALL POPCONTROL1B(branch)
1218 IF (branch .EQ. 0) THEN
1219 ynb(i) = 0.0
1220 xnb(i) = 0.0
1221 taubb(i) = 0.0
1222 ELSE
1223 tem = fr(i)*fr(i)*oc1(i)
1224 gfobnv = gmax*tem/((tem+cg)*bnv(i))
1225 tempb10 = xlinv(i)*ulow(i)**3*taubb(i)
1226 ulowb(i) = ulowb(i) + roll(i)*gfobnv*efact*xlinv(i)*3*ulow(i)**2*&
1227 & taubb(i)
1228 rollb(i) = rollb(i) + gfobnv*efact*tempb10
1229 gfobnvb = roll(i)*efact*tempb10
1230 efactb = roll(i)*gfobnv*tempb10
1231 taubb(i) = 0.0
1232 temp8 = (cg+tem)*bnv(i)
1233 tempb11 = gmax*gfobnvb/temp8
1234 tempb12 = -(tem*tempb11/temp8)
1235 temb = bnv(i)*tempb12 + tempb11
1236 bnvb(i) = bnvb(i) + (cg+tem)*tempb12
1237 frb(i) = frb(i) + oc1(i)*2*fr(i)*temb
1238 CALL POPCONTROL1B(branch)
1239 IF (branch .EQ. 0) THEN
1240 x3b = 0.0
1241 ELSE
1242 x3b = efactb
1243 END IF
1244 CALL POPCONTROL1B(branch)
1245 IF (branch .EQ. 0) THEN
1246 efactb = 0.0
1247 ELSE
1248 efactb = x3b
1249 END IF
1250 CALL POPREAL8(efact)
1251 IF (.NOT.oa(i) + 2. .LE. 0.0) frb(i) = frb(i) + ce*(oa(i)+2.)**(ce&
1252 & *(fr(i)/frc))*LOG(oa(i)+2.)*efactb/frc
1253 END IF
1254 END DO
1255 vbarb = 0.0
1256 ubarb = 0.0
1257 rulowb = 0.0
1258 DO i=ite,its,-1
1259 CALL POPCONTROL1B(branch)
1260 IF (branch .NE. 0) THEN
1261 vbarb(i) = vbarb(i) + rulow(i)*ynb(i)
1262 rulowb(i) = rulowb(i) + ubar(i)*xnb(i) + vbar(i)*ynb(i)
1263 ynb(i) = 0.0
1264 ubarb(i) = ubarb(i) + rulow(i)*xnb(i)
1265 xnb(i) = 0.0
1266 CALL POPCONTROL1B(branch)
1267 IF (branch .EQ. 0) frb(i) = 0.0
1268 tempb9 = var(i)*2.*od(i)*frb(i)
1269 bnvb(i) = bnvb(i) + rulow(i)*tempb9
1270 rulowb(i) = rulowb(i) + bnv(i)*tempb9
1271 frb(i) = 0.0
1272 IF (.NOT.bnv2(i, 1) .EQ. 0.0) bnv2b(i, 1) = bnv2b(i, 1) + bnvb(i)/&
1273 & (2.0*SQRT(bnv2(i, 1)))
1274 bnvb(i) = 0.0
1275 END IF
1276 END DO
1277 DO k=kpblmax,kpblmin,-1
1278 DO i=ite,its,-1
1279 CALL POPCONTROL1B(branch)
1280 IF (branch .NE. 0) THEN
1281 CALL POPREAL8(usqj(i, k))
1282 tmpb1 = usqjb(i, k)
1283 usqjb(i, k) = 0.0
1284 usqjb(i, 1) = usqjb(i, 1) + tmpb1
1285 END IF
1286 END DO
1287 END DO
1288 delks1b = 0.0
1289 DO k=kpblmax,kpblmin,-1
1290 DO i=ite,its,-1
1291 CALL POPCONTROL1B(branch)
1292 IF (branch .NE. 0) THEN
1293 tmpb0 = bnv2b(i, 1)
1294 rdelks = (prsl(i, k)-prsl(i, k+1))*delks1(i)
1295 CALL POPREAL8(usqj(i, 1))
1296 tmpb = usqjb(i, 1)
1297 usqjb(i, 1) = tmpb
1298 usqjb(i, k) = usqjb(i, k) + rdelks*tmpb
1299 CALL POPREAL8(bnv2(i, 1))
1300 rdelksb = bnv2(i, k)*tmpb0 + usqj(i, k)*tmpb
1301 bnv2b(i, 1) = tmpb0
1302 bnv2b(i, k) = bnv2b(i, k) + rdelks*tmpb0
1303 prslb(i, k) = prslb(i, k) + delks1(i)*rdelksb
1304 prslb(i, k+1) = prslb(i, k+1) - delks1(i)*rdelksb
1305 delks1b(i) = delks1b(i) + (prsl(i, k)-prsl(i, k+1))*rdelksb
1306 END IF
1307 END DO
1308 END DO
1309 DO i=ite,its,-1
1310 wtkbj = (prsl(i, 1)-prsl(i, 2))*delks1(i)
1311 CALL POPREAL8(usqj(i, 1))
1312 CALL POPREAL8(bnv2(i, 1))
1313 wtkbjb = bnv2(i, 1)*bnv2b(i, 1) + usqj(i, 1)*usqjb(i, 1)
1314 usqjb(i, 1) = wtkbj*usqjb(i, 1)
1315 bnv2b(i, 1) = wtkbj*bnv2b(i, 1)
1316 prslb(i, 1) = prslb(i, 1) + delks1(i)*wtkbjb
1317 prslb(i, 2) = prslb(i, 2) - delks1(i)*wtkbjb
1318 delks1b(i) = delks1b(i) + (prsl(i, 1)-prsl(i, 2))*wtkbjb
1319 END DO
1320 DO k=kte-1,kts,-1
1321 DO i=ite,its,-1
1322 CALL POPCONTROL1B(branch)
1323 IF (branch .NE. 0) velcob(i, k) = 0.0
1324 CALL POPREAL8(velco(i, k))
1325 rulowb(i) = rulowb(i) + velco(i, k)*velcob(i, k)
1326 velcob(i, k) = rulow(i)*velcob(i, k)
1327 tempb8 = rcs*0.5*velcob(i, k)
1328 u1b(i, k) = u1b(i, k) + ubar(i)*tempb8
1329 u1b(i, k+1) = u1b(i, k+1) + ubar(i)*tempb8
1330 ubarb(i) = ubarb(i) + (u1(i, k)+u1(i, k+1))*tempb8
1331 v1b(i, k) = v1b(i, k) + vbar(i)*tempb8
1332 v1b(i, k+1) = v1b(i, k+1) + vbar(i)*tempb8
1333 vbarb(i) = vbarb(i) + (v1(i, k)+v1(i, k+1))*tempb8
1334 velcob(i, k) = 0.0
1335 END DO
1336 END DO
1337 DO i=ite,its,-1
1338 ulowb(i) = ulowb(i) - rulowb(i)/ulow(i)**2
1339 rulowb(i) = 0.0
1340 CALL POPCONTROL1B(branch)
1341 IF (branch .EQ. 0) THEN
1342 ulowb(i) = 0.0
1343 x2b = 0.0
1344 ELSE
1345 x2b = ulowb(i)
1346 ulowb(i) = 0.0
1347 END IF
1348 IF (ubar(i)**2 + vbar(i)**2 .EQ. 0.0) THEN
1349 tempb7 = 0.0
1350 ELSE
1351 tempb7 = x2b/(2.0*SQRT(ubar(i)**2+vbar(i)**2))
1352 END IF
1353 ubarb(i) = ubarb(i) + 2*ubar(i)*tempb7
1354 vbarb(i) = vbarb(i) + 2*vbar(i)*tempb7
1355 END DO
1356 vtjb = 0.0
1357 vtkb = 0.0
1358 DO k=kte-1,kts,-1
1359 DO i=ite,its,-1
1360 CALL POPCONTROL1B(branch)
1361 IF (branch .EQ. 0) bnv2b(i, k) = 0.0
1362 rdz = 1./(zl(i, k+1)-zl(i, k))
1363 temp6 = vtk(i, k+1) + vtk(i, k)
1364 temp7 = vtk(i, k+1) - vtk(i, k)
1365 tempb5 = g*2.0*bnv2b(i, k)/temp6
1366 tempb6 = -(rdz*temp7*tempb5/temp6)
1367 rdzb = temp7*tempb5
1368 vtkb(i, k+1) = vtkb(i, k+1) + tempb6 + rdz*tempb5
1369 vtkb(i, k) = vtkb(i, k) + tempb6 - rdz*tempb5
1370 bnv2b(i, k) = 0.0
1371 CALL POPCONTROL1B(branch)
1372 IF (branch .EQ. 0) THEN
1373 usqjb(i, k) = 0.0
1374 bvf2b = 0.0
1375 shr2b = 0.0
1376 ELSE
1377 shr2 = max2*rdz*rdz
1378 bvf2b = usqjb(i, k)/shr2
1379 shr2b = -(bvf2*usqjb(i, k)/shr2**2)
1380 usqjb(i, k) = 0.0
1381 END IF
1382 ti = 2.0/(t1(i, k)+t1(i, k+1))
1383 CALL POPREAL8(bvf2)
1384 temp5 = vtj(i, k+1) - vtj(i, k)
1385 tempb4 = g*ti*bvf2b
1386 rdzb = rdzb + max2*2*rdz*shr2b + temp5*tempb4
1387 vtjb(i, k+1) = vtjb(i, k+1) + rdz*tempb4
1388 vtjb(i, k) = vtjb(i, k) - rdz*tempb4
1389 tib = g*(g/cp+rdz*temp5)*bvf2b
1390 max2b = rdz**2*shr2b
1391 CALL POPCONTROL1B(branch)
1392 IF (branch .EQ. 0) THEN
1393 CALL POPREAL8(max2)
1394 dw2b = 0.0
1395 ELSE
1396 CALL POPREAL8(max2)
1397 dw2b = max2b
1398 END IF
1399 tem1b = rcl*2*tem1*dw2b
1400 tem2b = rcl*2*tem2*dw2b
1401 CALL POPREAL8(tem2)
1402 v1b(i, k) = v1b(i, k) + tem2b
1403 v1b(i, k+1) = v1b(i, k+1) - tem2b
1404 CALL POPREAL8(tem1)
1405 u1b(i, k) = u1b(i, k) + tem1b
1406 u1b(i, k+1) = u1b(i, k+1) - tem1b
1407 temp4 = zl(i, k+1) - zl(i, k)
1408 tempb2 = -(rdzb/temp4**2)
1409 zlb(i, k+1) = zlb(i, k+1) + tempb2
1410 zlb(i, k) = zlb(i, k) - tempb2
1411 temp3 = t1(i, k) + t1(i, k+1)
1412 tempb3 = -(2.0*tib/temp3**2)
1413 t1b(i, k) = t1b(i, k) + tempb3
1414 t1b(i, k+1) = t1b(i, k+1) + tempb3
1415 END DO
1416 END DO
1417 delksb = 0.0
1418 DO k=kpblmax,kts,-1
1419 DO i=ite,its,-1
1420 CALL POPCONTROL1B(branch)
1421 IF (branch .NE. 0) THEN
1422 rdelks = del(i, k)*delks(i)
1423 rdelksb = ro(i, k)*rollb(i)
1424 rob(i, k) = rob(i, k) + rdelks*rollb(i)
1425 rcsks = rcs*del(i, k)*delks(i)
1426 rcsksb = u1(i, k)*ubarb(i) + v1(i, k)*vbarb(i)
1427 v1b(i, k) = v1b(i, k) + rcsks*vbarb(i)
1428 u1b(i, k) = u1b(i, k) + rcsks*ubarb(i)
1429 delb(i, k) = delb(i, k) + rcs*delks(i)*rcsksb + delks(i)*rdelksb
1430 delksb(i) = delksb(i) + rcs*del(i, k)*rcsksb + del(i, k)*rdelksb
1431 END IF
1432 END DO
1433 END DO
1434 DO i=ite,its,-1
1435 temp2 = prsl(i, 1) - prsl(i, kbl(i))
1436 tempb0 = -(delks1b(i)/temp2**2)
1437 prslb(i, 1) = prslb(i, 1) + tempb0
1438 prslb(i, kbl(i)) = prslb(i, kbl(i)) - tempb0
1439 delks1b(i) = 0.0
1440 temp1 = prsi(i, 1) - prsi(i, kbl(i))
1441 tempb1 = -(delksb(i)/temp1**2)
1442 prsib(i, 1) = prsib(i, 1) + tempb1
1443 prsib(i, kbl(i)) = prsib(i, kbl(i)) - tempb1
1444 delksb(i) = 0.0
1445 END DO
1446 DO k=kte,kts,-1
1447 DO i=ite,its,-1
1448 tempb = vtkb(i, k)/prslk(i, k)
1449 temp0 = rd*vtj(i, k)
1450 prslb(i, k) = prslb(i, k) + rob(i, k)/temp0
1451 vtjb(i, k) = vtjb(i, k) + tempb - prsl(i, k)*rd*rob(i, k)/temp0**2
1452 rob(i, k) = 0.0
1453 prslkb(i, k) = prslkb(i, k) - vtj(i, k)*tempb/prslk(i, k)
1454 vtkb(i, k) = 0.0
1455 t1b(i, k) = t1b(i, k) + (fv*q1(i, k)+1.)*vtjb(i, k)
1456 q1b(i, k) = q1b(i, k) + t1(i, k)*fv*vtjb(i, k)
1457 vtjb(i, k) = 0.0
1458 END DO
1459 END DO
1460 DO k=kte,kts,-1
1461 DO i=ite,its,-1
1462 dtauy2db(i, k) = 0.0
1463 dtaux2db(i, k) = 0.0
1464 END DO
1465 END DO
1466 END SUBROUTINE GWDO2D_B
1468 end module a_module_bl_gwdo