| blob | 8d66a24fc493faf2aadac94309a007ef66a1b6d9 |
1 !WRF:model_layer:physics
2 !
3 !
4 !
5 !
6 module g_module_bl_gwdo
7 contains
9 ! Generated by TAPENADE (INRIA, Tropics team)
10 ! Tapenade 3.10 (r5498) - 20 Jan 2015 09:48
11 !
12 ! Differentiation of gwdo in forward (tangent) mode:
13 ! variations of useful results: rublten dusfcg dvsfcg dtauy3d
14 ! rvblten 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:in u3d:in
18 ! dusfcg:in-out z:in dvsfcg:in-out dtauy3d:in-out
19 ! rvblten:in-out t3d:in qv3d:in pi3d:in v3d:in dtaux3d:in-out
20 !-------------------------------------------------------------------------------
21 SUBROUTINE G_GWDO(u3d, u3dd, v3d, v3dd, t3d, t3dd, qv3d, qv3dd, p3d, &
22 & p3di, p3did, pi3d, pi3dd, z, zd, rublten, rubltend, rvblten, rvbltend&
23 & , dtaux3d, dtaux3dd, dtauy3d, dtauy3dd, dusfcg, dusfcgd, dvsfcg, &
24 & dvsfcgd, var2d, oc12d, oa2d1, oa2d2, oa2d3, oa2d4, ol2d1, ol2d2, ol2d3&
25 & , ol2d4, znu, znw, p_top, cp, g, rd, rv, ep1, pi, dt, dx, &
26 & kpbl2d, itimestep, ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, &
27 & kms, kme, its, ite, jts, jte, kts, kte)
28 IMPLICIT NONE
29 !
30 !-------------------------------------------------------------------------------
31 !
32 !-- u3d 3d u-velocity interpolated to theta points (m/s)
33 !-- v3d 3d v-velocity interpolated to theta points (m/s)
34 !-- t3d temperature (k)
35 !-- qv3d 3d water vapor mixing ratio (kg/kg)
36 !-- p3d 3d pressure (pa)
37 !-- p3di 3d pressure (pa) at interface level
38 !-- pi3d 3d exner function (dimensionless)
39 !-- rublten u tendency due to pbl parameterization (m/s/s)
40 !-- rvblten v tendency due to pbl parameterization (m/s/s)
41 !-- znu eta values (sigma values)
42 !-- cp heat capacity at constant pressure for dry air (j/kg/k)
43 !-- g acceleration due to gravity (m/s^2)
44 !-- rd gas constant for dry air (j/kg/k)
45 !-- z height above sea level (m)
46 !-- rv gas constant for water vapor (j/kg/k)
47 !-- dt time step (s)
48 !-- dx model grid interval (m)
49 !-- ep1 constant for virtual temperature (r_v/r_d - 1) (dimensionless)
50 !-- ids start index for i in domain
51 !-- ide end index for i in domain
52 !-- jds start index for j in domain
53 !-- jde end index for j in domain
54 !-- kds start index for k in domain
55 !-- kde end index for k in domain
56 !-- ims start index for i in memory
57 !-- ime end index for i in memory
58 !-- jms start index for j in memory
59 !-- jme end index for j in memory
60 !-- kms start index for k in memory
61 !-- kme end index for k in memory
62 !-- its start index for i in tile
63 !-- ite end index for i in tile
64 !-- jts start index for j in tile
65 !-- jte end index for j in tile
66 !-- kts start index for k in tile
67 !-- kte end index for k in tile
68 !
69 !-------------------------------------------------------------------------------
70 INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, ims, ime, jms, &
71 & jme, kms, kme, its, ite, jts, jte, kts, kte
72 INTEGER, INTENT(IN) :: itimestep
73 !
74 REAL, INTENT(IN) :: dt, dx, cp, g, rd, rv, ep1, pi
75 !
76 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: qv3d, p3d, &
77 & pi3d, t3d, z
78 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: qv3dd, pi3dd&
79 & , t3dd, zd
80 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: p3di
81 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: p3did
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) :: rubltend&
86 & , rvbltend
87 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: dtaux3d, &
88 & dtauy3d
89 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: dtaux3dd&
90 & , dtauy3dd
91 !
92 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: u3d, v3d
93 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: u3dd, v3dd
94 !
95 INTEGER, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: kpbl2d
96 REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: dusfcg, dvsfcg
97 REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: dusfcgd, dvsfcgd
98 !
99 REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: var2d, oc12d, oa2d1, &
100 & oa2d2, oa2d3, oa2d4, ol2d1, ol2d2, ol2d3, ol2d4
101 !
102 REAL, DIMENSION(kms:kme), OPTIONAL, INTENT(IN) :: znu, znw
103 !
104 REAL, OPTIONAL, INTENT(IN) :: p_top
105 !
106 !local
107 !
108 REAL, DIMENSION(its:ite, kts:kte) :: delprsi, pdh
109 REAL, DIMENSION(its:ite, kts:kte) :: delprsid, pdhd
110 REAL, DIMENSION(its:ite, kts:kte + 1) :: pdhi
111 REAL, DIMENSION(its:ite, kts:kte+1) :: pdhid
112 REAL, DIMENSION(its:ite, 4) :: oa4, ol4
113 INTEGER :: i, j, k, kdt, kpblmax
114 !
115 DO k=kts,kte
116 IF (znu(k) .GT. 0.6) kpblmax = k + 1
117 END DO
118 delprsid = 0.0
119 pdhid = 0.0
120 pdhd = 0.0
121 !
122 DO j=jts,jte
123 DO k=kts,kte+1
124 DO i=its,ite
125 IF (k .LE. kte) THEN
126 pdhd(i, k) = 0.0
127 pdh(i, k) = p3d(i, k, j)
128 END IF
129 pdhid(i, k) = p3did(i, k, j)
130 pdhi(i, k) = p3di(i, k, j)
131 END DO
132 END DO
133 !
134 DO k=kts,kte
135 DO i=its,ite
136 delprsid(i, k) = pdhid(i, k) - pdhid(i, k+1)
137 delprsi(i, k) = pdhi(i, k) - pdhi(i, k+1)
138 END DO
139 END DO
140 DO i=its,ite
141 oa4(i, 1) = oa2d1(i, j)
142 oa4(i, 2) = oa2d2(i, j)
143 oa4(i, 3) = oa2d3(i, j)
144 oa4(i, 4) = oa2d4(i, j)
145 ol4(i, 1) = ol2d1(i, j)
146 ol4(i, 2) = ol2d2(i, j)
147 ol4(i, 3) = ol2d3(i, j)
148 ol4(i, 4) = ol2d4(i, j)
149 END DO
150 CALL GWDO2D_D(dudt=rublten(ims, kms, j), dudtd=rubltend(ims, kms, j)&
151 & , dvdt=rvblten(ims, kms, j), dvdtd=rvbltend(ims, kms, j), &
152 & dtaux2d=dtaux3d(ims, kms, j), dtaux2dd=dtaux3dd(ims, kms, j)&
153 & , dtauy2d=dtauy3d(ims, kms, j), dtauy2dd=dtauy3dd(ims, kms, &
154 & j), u1=u3d(ims, kms, j), u1d=u3dd(ims, kms, j), v1=v3d(ims, &
155 & kms, j), v1d=v3dd(ims, kms, j), t1=t3d(ims, kms, j), t1d=&
156 & t3dd(ims, kms, j), q1=qv3d(ims, kms, j), q1d=qv3dd(ims, kms&
157 & , j), del=delprsi(its, kts), deld=delprsid(its, kts), prsi=&
158 & pdhi(its, kts), prsid=pdhid(its, kts), prsl=pdh(its, kts), &
159 & prsld=pdhd(its, kts), prslk=pi3d(ims, kms, j), prslkd=pi3dd(&
160 & ims, kms, j), zl=z(ims, kms, j), zld=zd(ims, kms, j), rcl=&
161 & 1.0, kpblmax=kpblmax, var=var2d(ims, j), oc1=oc12d(ims, j), &
162 & oa4=oa4, ol4=ol4, dusfc=dusfcg(ims, j), dusfcd=dusfcgd(ims, &
163 & j), dvsfc=dvsfcg(ims, j), dvsfcd=dvsfcgd(ims, j), g=g, cp=cp&
164 & , rd=rd, rv=rv, fv=ep1, pi=pi, dxmeter=dx, deltim=dt, kpbl=&
165 & kpbl2d(ims, j), kdt=itimestep, lat=j, ids=ids, ide=ide, jds=&
166 & jds, jde=jde, kds=kds, kde=kde, ims=ims, ime=ime, jms=jms, &
167 & jme=jme, kms=kms, kme=kme, its=its, ite=ite, jts=jts, jte=&
168 & jte, kts=kts, kte=kte)
169 END DO
170 END SUBROUTINE G_GWDO
172 ! Differentiation of gwdo2d in forward (tangent) mode:
173 ! variations of useful results: dvsfc dvdt dtauy2d dusfc dudt
174 ! dtaux2d
175 ! with respect to varying inputs: v1 dvsfc dvdt dtauy2d prsi
176 ! prsl dusfc del t1 q1 dudt dtaux2d u1 zl prslk
177 !-------------------------------------------------------------------------------
178 !
179 !-------------------------------------------------------------------------------
180 SUBROUTINE GWDO2D_D(dudt, dudtd, dvdt, dvdtd, dtaux2d, dtaux2dd, dtauy2d&
181 & , dtauy2dd, u1, u1d, v1, v1d, t1, t1d, q1, q1d, del, deld, prsi, prsid&
182 & , prsl, prsld, prslk, prslkd, zl, zld, rcl, kpblmax, var, oc1, oa4, &
183 & ol4, dusfc, dusfcd, dvsfc, dvsfcd, g, cp, rd, rv, fv, pi, dxmeter, &
184 & deltim, kpbl, kdt, lat, ids, ide, jds, jde, kds, kde, ims, ime, jms, &
185 & jme, kms, kme, its, ite, jts, jte, kts, kte)
186 IMPLICIT NONE
187 !-------------------------------------------------------------------------------
188 INTEGER :: kdt, lat, latd, lond, kpblmax, ids, ide, jds, jde, kds, kde&
189 & , ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts, kte
190 !
191 REAL :: g, rd, rv, fv, cp, pi, dxmeter, deltim, rcl
192 REAL :: dudt(ims:ime, kms:kme), dvdt(ims:ime, kms:kme), dtaux2d(ims:&
193 & ime, kms:kme), dtauy2d(ims:ime, kms:kme), u1(ims:ime, kms:kme), v1(ims&
194 & :ime, kms:kme), t1(ims:ime, kms:kme), q1(ims:ime, kms:kme), zl(ims:ime&
195 & , kms:kme), prsl(its:ite, kts:kte), prslk(ims:ime, kms:kme)
196 REAL :: dudtd(ims:ime, kms:kme), dvdtd(ims:ime, kms:kme), dtaux2dd(ims&
197 & :ime, kms:kme), dtauy2dd(ims:ime, kms:kme), u1d(ims:ime, kms:kme), v1d&
198 & (ims:ime, kms:kme), t1d(ims:ime, kms:kme), q1d(ims:ime, kms:kme), zld(&
199 & ims:ime, kms:kme), prsld(its:ite, kts:kte), prslkd(ims:ime, kms:kme)
200 REAL :: prsi(its:ite, kts:kte+1), del(its:ite, kts:kte)
201 REAL :: prsid(its:ite, kts:kte+1), deld(its:ite, kts:kte)
202 REAL :: oa4(its:ite, 4), ol4(its:ite, 4)
203 !
204 INTEGER :: kpbl(ims:ime)
205 REAL :: var(ims:ime), oc1(ims:ime), dusfc(ims:ime), dvsfc(ims:ime)
206 REAL :: dusfcd(ims:ime), dvsfcd(ims:ime)
207 !
208 ! critical richardson number for wave breaking : ! larger drag with larger value
209 !
210 REAL, PARAMETER :: ric=0.25
211 !
212 REAL, PARAMETER :: dw2min=1.
213 REAL, PARAMETER :: rimin=-100.
214 REAL, PARAMETER :: bnv2min=1.0e-5
215 REAL, PARAMETER :: efmin=0.0
216 REAL, PARAMETER :: efmax=10.0
217 REAL, PARAMETER :: xl=4.0e4
218 REAL, PARAMETER :: critac=1.0e-5
219 REAL, PARAMETER :: gmax=1.
220 REAL, PARAMETER :: veleps=1.0
221 REAL, PARAMETER :: factop=0.5
222 REAL, PARAMETER :: frc=1.0
223 REAL, PARAMETER :: ce=0.8
224 REAL, PARAMETER :: cg=0.5
225 INTEGER, PARAMETER :: kpblmin=2
226 !
227 ! local variables
228 !
229 INTEGER :: i, k, lcap, lcapp1, nwd, idir, klcap, kp1, ikount, kk
230 !
231 REAL :: rcs, rclcs, csg, fdir, cleff, cs, rcsks, wdir, ti, rdz, temp, &
232 & tem2, dw2, shr2, bvf2, rdelks, wtkbj, tem, gfobnv, hd, fro, rim, temc&
233 & , tem1, efact, temv, dtaux, dtauy
234 REAL :: rcsksd, tid, rdzd, tem2d, dw2d, shr2d, bvf2d, rdelksd, wtkbjd&
235 & , temd, gfobnvd, hdd, temcd, tem1d, efactd, dtauxd, dtauyd
236 !
237 LOGICAL :: ldrag(its:ite), icrilv(its:ite), flag(its:ite), kloop1(its:&
238 & ite)
239 !
240 REAL :: taub(its:ite), taup(its:ite, kts:kte+1), xn(its:ite), yn(its:&
241 & ite), ubar(its:ite), vbar(its:ite), fr(its:ite), ulow(its:ite), rulow(&
242 & its:ite), bnv(its:ite), oa(its:ite), ol(its:ite), roll(its:ite), dtfac&
243 & (its:ite), brvf(its:ite), xlinv(its:ite), delks(its:ite), delks1(its:&
244 & ite), bnv2(its:ite, kts:kte), usqj(its:ite, kts:kte), taud(its:ite, &
245 & kts:kte), ro(its:ite, kts:kte), vtk(its:ite, kts:kte), vtj(its:ite, &
246 & kts:kte), zlowtop(its:ite), velco(its:ite, kts:kte-1), coefm(its:ite)
247 REAL :: taubd(its:ite), taupd(its:ite, kts:kte+1), xnd(its:ite), ynd(&
248 & its:ite), ubard(its:ite), vbard(its:ite), frd(its:ite), ulowd(its:ite)&
249 & , rulowd(its:ite), bnvd(its:ite), rolld(its:ite), dtfacd(its:ite), &
250 & brvfd(its:ite), delksd(its:ite), delks1d(its:ite), bnv2d(its:ite, kts:&
251 & kte), usqjd(its:ite, kts:kte), taudd(its:ite, kts:kte), rod(its:ite, &
252 & kts:kte), vtkd(its:ite, kts:kte), vtjd(its:ite, kts:kte), velcod(its:&
253 & ite, kts:kte-1)
254 !
255 INTEGER :: kbl(its:ite), klowtop(its:ite)
256 !
257 LOGICAL :: iope
258 INTEGER, PARAMETER :: mdir=8
259 INTEGER :: nwdir(mdir)
260 !
261 ! variables for flow-blocking drag
262 !
263 REAL, PARAMETER :: frmax=10.
264 REAL, PARAMETER :: olmin=1.0e-5
265 REAL, PARAMETER :: odmin=0.1
266 REAL, PARAMETER :: odmax=10.
267 REAL, PARAMETER :: erad=6371.315e+3
268 INTEGER :: komax(its:ite)
269 INTEGER :: kblk
270 REAL :: cd
271 REAL :: zblk, tautem
272 REAL :: zblkd, tautemd
273 REAL :: pe, ke
274 REAL :: delx, dely, dxy4(4), dxy4p(4)
275 REAL :: dxy(its:ite), dxyp(its:ite)
276 REAL :: ol4p(4), olp(its:ite), od(its:ite)
277 REAL :: taufb(its:ite, kts:kte+1)
278 REAL :: taufbd(its:ite, kts:kte+1)
279 REAL :: result1
280 REAL :: result1d
281 REAL :: arg1
282 REAL :: arg1d
283 REAL :: pwx1
284 REAL :: pwy1
285 REAL :: pwy1d
286 REAL :: arg10
287 REAL :: pwy10
288 REAL :: max2d
289 REAL :: x4
290 REAL :: x3
291 REAL :: x2
292 REAL :: x2d
293 INTEGER :: x1
294 REAL :: x4d
295 REAL :: max3
296 REAL :: max2
297 REAL :: max1
298 REAL :: x3d
299 REAL :: y1
300 REAL :: y1d
301 DATA nwdir /6, 7, 5, 8, 2, 3, 1, 4/
302 !
303 !---- constants
304 !
305 rcs = SQRT(rcl)
306 result1 = SQRT(rcl)
307 cs = 1./result1
308 csg = cs*g
309 lcap = kte
310 lcapp1 = lcap + 1
311 fdir = mdir/(2.0*pi)
312 !
313 !--- calculate length of grid for flow-blocking drag
314 !
315 delx = dxmeter
316 dely = dxmeter
317 dxy4(1) = delx
318 dxy4(2) = dely
319 arg1 = delx*delx + dely*dely
320 dxy4(3) = SQRT(arg1)
321 dxy4(4) = dxy4(3)
322 dxy4p(1) = dxy4(2)
323 dxy4p(2) = dxy4(1)
324 dxy4p(3) = dxy4(4)
325 dxy4p(4) = dxy4(3)
326 !
327 !
328 !-----initialize arrays
329 !
330 dtaux = 0.0
331 dtauy = 0.0
332 DO i=its,ite
333 klowtop(i) = 0
334 kbl(i) = 0
335 END DO
336 !
337 DO i=its,ite
338 xn(i) = 0.0
339 yn(i) = 0.0
340 ubar(i) = 0.0
341 vbar(i) = 0.0
342 roll(i) = 0.0
343 taub(i) = 0.0
344 taup(i, 1) = 0.0
345 oa(i) = 0.0
346 ol(i) = 0.0
347 ulow(i) = 0.0
348 dtfac(i) = 1.0
349 ldrag(i) = .false.
350 icrilv(i) = .false.
351 flag(i) = .true.
352 END DO
353 !
354 DO k=kts,kte
355 DO i=its,ite
356 usqj(i, k) = 0.0
357 bnv2(i, k) = 0.0
358 vtj(i, k) = 0.0
359 vtk(i, k) = 0.0
360 taup(i, k) = 0.0
361 taud(i, k) = 0.0
362 dtaux2dd(i, k) = 0.0
363 dtaux2d(i, k) = 0.0
364 dtauy2dd(i, k) = 0.0
365 dtauy2d(i, k) = 0.0
366 END DO
367 END DO
368 !
369 DO i=its,ite
370 taup(i, kte+1) = 0.0
371 xlinv(i) = 1.0/xl
372 END DO
373 !
374 ! initialize array for flow-blocking drag
375 !
376 taufb(its:ite, kts:kte+1) = 0.0
377 komax(its:ite) = 0
378 vtjd = 0.0
379 vtkd = 0.0
380 rod = 0.0
381 !
382 DO k=kts,kte
383 DO i=its,ite
384 vtjd(i, k) = t1d(i, k)*(1.+fv*q1(i, k)) + t1(i, k)*fv*q1d(i, k)
385 vtj(i, k) = t1(i, k)*(1.+fv*q1(i, k))
386 vtkd(i, k) = (vtjd(i, k)*prslk(i, k)-vtj(i, k)*prslkd(i, k))/prslk&
387 & (i, k)**2
388 vtk(i, k) = vtj(i, k)/prslk(i, k)
389 ! density kg/m**3
390 rod(i, k) = (prsld(i, k)*vtj(i, k)/rd-prsl(i, k)*vtjd(i, k)/rd)/&
391 & vtj(i, k)**2
392 ro(i, k) = 1./rd*prsl(i, k)/vtj(i, k)
393 END DO
394 END DO
395 !
396 ! determine reference level: maximum of 2*var and pbl heights
397 !
398 DO i=its,ite
399 zlowtop(i) = 2.*var(i)
400 END DO
401 !
402 DO i=its,ite
403 kloop1(i) = .true.
404 END DO
405 !
406 DO k=kts+1,kte
407 DO i=its,ite
408 IF (kloop1(i) .AND. zl(i, k) - zl(i, 1) .GE. zlowtop(i)) THEN
409 klowtop(i) = k + 1
410 kloop1(i) = .false.
411 END IF
412 END DO
413 END DO
414 !
415 DO i=its,ite
416 IF (kpbl(i) .LT. klowtop(i)) THEN
417 kbl(i) = klowtop(i)
418 ELSE
419 kbl(i) = kpbl(i)
420 END IF
421 IF (kbl(i) .GT. kpblmax) THEN
422 x1 = kpblmax
423 ELSE
424 x1 = kbl(i)
425 END IF
426 IF (x1 .LT. kpblmin) THEN
427 kbl(i) = kpblmin
428 ELSE
429 kbl(i) = x1
430 END IF
431 END DO
432 !
433 ! determine the level of maximum orographic height
434 !
435 komax(:) = kbl(:)
436 delksd = 0.0
437 delks1d = 0.0
438 !
439 DO i=its,ite
440 delksd(i) = -((prsid(i, 1)-prsid(i, kbl(i)))/(prsi(i, 1)-prsi(i, kbl&
441 & (i)))**2)
442 delks(i) = 1.0/(prsi(i, 1)-prsi(i, kbl(i)))
443 delks1d(i) = -((prsld(i, 1)-prsld(i, kbl(i)))/(prsl(i, 1)-prsl(i, &
444 & kbl(i)))**2)
445 delks1(i) = 1.0/(prsl(i, 1)-prsl(i, kbl(i)))
446 END DO
447 rolld = 0.0
448 vbard = 0.0
449 ubard = 0.0
450 !
451 ! compute low level averages within pbl
452 !
453 DO k=kts,kpblmax
454 DO i=its,ite
455 IF (k .LT. kbl(i)) THEN
456 rcsksd = rcs*(deld(i, k)*delks(i)+del(i, k)*delksd(i))
457 rcsks = rcs*del(i, k)*delks(i)
458 rdelksd = deld(i, k)*delks(i) + del(i, k)*delksd(i)
459 rdelks = del(i, k)*delks(i)
460 ! pbl u mean
461 ubard(i) = ubard(i) + rcsksd*u1(i, k) + rcsks*u1d(i, k)
462 ubar(i) = ubar(i) + rcsks*u1(i, k)
463 ! pbl v mean
464 vbard(i) = vbard(i) + rcsksd*v1(i, k) + rcsks*v1d(i, k)
465 vbar(i) = vbar(i) + rcsks*v1(i, k)
466 ! ro mean
467 rolld(i) = rolld(i) + rdelksd*ro(i, k) + rdelks*rod(i, k)
468 roll(i) = roll(i) + rdelks*ro(i, k)
469 END IF
470 END DO
471 END DO
472 !
473 ! figure out low-level horizontal wind direction
474 !
475 ! nwd 1 2 3 4 5 6 7 8
476 ! wd w s sw nw e n ne se
477 !
478 DO i=its,ite
479 wdir = ATAN2(ubar(i), vbar(i)) + pi
480 idir = MOD(NINT(fdir*wdir), mdir) + 1
481 nwd = nwdir(idir)
482 oa(i) = (1-2*INT((nwd-1)/4))*oa4(i, MOD(nwd-1, 4)+1)
483 ol(i) = ol4(i, MOD(nwd-1, 4)+1)
484 !
485 !----- compute orographic width along (ol) and perpendicular (olp)
486 !----- the direction of wind
487 !
488 ol4p(1) = ol4(i, 2)
489 ol4p(2) = ol4(i, 1)
490 ol4p(3) = ol4(i, 4)
491 ol4p(4) = ol4(i, 3)
492 olp(i) = ol4p(MOD(nwd-1, 4)+1)
493 IF (ol(i) .LT. olmin) THEN
494 max1 = olmin
495 ELSE
496 max1 = ol(i)
497 END IF
498 !
499 !----- compute orographic direction (horizontal orographic aspect ratio)
500 !
501 od(i) = olp(i)/max1
502 IF (od(i) .GT. odmax) THEN
503 od(i) = odmax
504 ELSE
505 od(i) = od(i)
506 END IF
507 IF (od(i) .LT. odmin) THEN
508 od(i) = odmin
509 ELSE
510 od(i) = od(i)
511 END IF
512 !
513 !----- compute length of grid in the along(dxy) and cross(dxyp) wind directions
514 !
515 dxy(i) = dxy4(MOD(nwd-1, 4)+1)
516 dxyp(i) = dxy4p(MOD(nwd-1, 4)+1)
517 END DO
518 bnv2d = 0.0
519 usqjd = 0.0
520 !
521 !--- saving richardson number in usqj for migwdi
522 !
523 DO k=kts,kte-1
524 DO i=its,ite
525 tid = -(2.0*(t1d(i, k)+t1d(i, k+1))/(t1(i, k)+t1(i, k+1))**2)
526 ti = 2.0/(t1(i, k)+t1(i, k+1))
527 rdzd = -((zld(i, k+1)-zld(i, k))/(zl(i, k+1)-zl(i, k))**2)
528 rdz = 1./(zl(i, k+1)-zl(i, k))
529 tem1d = u1d(i, k) - u1d(i, k+1)
530 tem1 = u1(i, k) - u1(i, k+1)
531 tem2d = v1d(i, k) - v1d(i, k+1)
532 tem2 = v1(i, k) - v1(i, k+1)
533 dw2d = rcl*(tem1d*tem1+tem1*tem1d+tem2d*tem2+tem2*tem2d)
534 dw2 = rcl*(tem1*tem1+tem2*tem2)
535 IF (dw2 .LT. dw2min) THEN
536 max2 = dw2min
537 max2d = 0.0
538 ELSE
539 max2d = dw2d
540 max2 = dw2
541 END IF
542 shr2d = (max2d*rdz+max2*rdzd)*rdz + max2*rdz*rdzd
543 shr2 = max2*rdz*rdz
544 bvf2d = g*((rdzd*(vtj(i, k+1)-vtj(i, k))+rdz*(vtjd(i, k+1)-vtjd(i&
545 & , k)))*ti+(g/cp+rdz*(vtj(i, k+1)-vtj(i, k)))*tid)
546 bvf2 = g*(g/cp+rdz*(vtj(i, k+1)-vtj(i, k)))*ti
547 IF (bvf2/shr2 .LT. rimin) THEN
548 usqjd(i, k) = 0.0
549 usqj(i, k) = rimin
550 ELSE
551 usqjd(i, k) = (bvf2d*shr2-bvf2*shr2d)/shr2**2
552 usqj(i, k) = bvf2/shr2
553 END IF
554 bnv2d(i, k) = (2.0*g*(rdzd*(vtk(i, k+1)-vtk(i, k))+rdz*(vtkd(i, k+&
555 & 1)-vtkd(i, k)))*(vtk(i, k+1)+vtk(i, k))-2.0*g*rdz*(vtk(i, k+1)-&
556 & vtk(i, k))*(vtkd(i, k+1)+vtkd(i, k)))/(vtk(i, k+1)+vtk(i, k))**2
557 bnv2(i, k) = 2.0*g*rdz*(vtk(i, k+1)-vtk(i, k))/(vtk(i, k+1)+vtk(i&
558 & , k))
559 IF (bnv2(i, k) .LT. bnv2min) THEN
560 bnv2d(i, k) = 0.0
561 bnv2(i, k) = bnv2min
562 ELSE
563 bnv2(i, k) = bnv2(i, k)
564 END IF
565 END DO
566 END DO
567 rulowd = 0.0
568 ulowd = 0.0
569 !
570 !----compute the "low level" or 1/3 wind magnitude (m/s)
571 !
572 DO i=its,ite
573 arg1d = ubard(i)*ubar(i) + ubar(i)*ubard(i) + vbard(i)*vbar(i) + &
574 & vbar(i)*vbard(i)
575 arg1 = ubar(i)*ubar(i) + vbar(i)*vbar(i)
576 IF (arg1 .EQ. 0.0) THEN
577 x2d = 0.0
578 ELSE
579 x2d = arg1d/(2.0*SQRT(arg1))
580 END IF
581 x2 = SQRT(arg1)
582 IF (x2 .LT. 1.0) THEN
583 ulowd(i) = 0.0
584 ulow(i) = 1.0
585 ELSE
586 ulowd(i) = x2d
587 ulow(i) = x2
588 END IF
589 rulowd(i) = -(ulowd(i)/ulow(i)**2)
590 rulow(i) = 1./ulow(i)
591 END DO
592 velcod = 0.0
593 !
594 DO k=kts,kte-1
595 DO i=its,ite
596 velcod(i, k) = 0.5*rcs*((u1d(i, k)+u1d(i, k+1))*ubar(i)+(u1(i, k)+&
597 & u1(i, k+1))*ubard(i)+(v1d(i, k)+v1d(i, k+1))*vbar(i)+(v1(i, k)+&
598 & v1(i, k+1))*vbard(i))
599 velco(i, k) = 0.5*rcs*((u1(i, k)+u1(i, k+1))*ubar(i)+(v1(i, k)+v1(&
600 & i, k+1))*vbar(i))
601 velcod(i, k) = velcod(i, k)*rulow(i) + velco(i, k)*rulowd(i)
602 velco(i, k) = velco(i, k)*rulow(i)
603 IF (velco(i, k) .LT. veleps .AND. velco(i, k) .GT. 0.) THEN
604 velcod(i, k) = 0.0
605 velco(i, k) = veleps
606 END IF
607 END DO
608 END DO
609 !
610 ! no drag when critical level in the base layer
611 !
612 DO i=its,ite
613 ldrag(i) = velco(i, 1) .LE. 0.
614 END DO
615 !
616 ! no drag when velco.lt.0
617 !
618 DO k=kpblmin,kpblmax
619 DO i=its,ite
620 IF (k .LT. kbl(i)) ldrag(i) = ldrag(i) .OR. velco(i, k) .LE. 0.
621 END DO
622 END DO
623 !
624 ! no drag when bnv2.lt.0
625 !
626 DO k=kts,kpblmax
627 DO i=its,ite
628 IF (k .LT. kbl(i)) ldrag(i) = ldrag(i) .OR. bnv2(i, k) .LT. 0.
629 END DO
630 END DO
631 !
632 !-----the low level weighted average ri is stored in usqj(1,1; im)
633 !-----the low level weighted average n**2 is stored in bnv2(1,1; im)
634 !---- this is called bnvl2 in phys_gwd_alpert_sub not bnv2
635 !---- rdelks (del(k)/delks) vert ave factor so we can * instead of /
636 !
637 DO i=its,ite
638 wtkbjd = (prsld(i, 1)-prsld(i, 2))*delks1(i) + (prsl(i, 1)-prsl(i, 2&
639 & ))*delks1d(i)
640 wtkbj = (prsl(i, 1)-prsl(i, 2))*delks1(i)
641 bnv2d(i, 1) = wtkbjd*bnv2(i, 1) + wtkbj*bnv2d(i, 1)
642 bnv2(i, 1) = wtkbj*bnv2(i, 1)
643 usqjd(i, 1) = wtkbjd*usqj(i, 1) + wtkbj*usqjd(i, 1)
644 usqj(i, 1) = wtkbj*usqj(i, 1)
645 END DO
646 !
647 DO k=kpblmin,kpblmax
648 DO i=its,ite
649 IF (k .LT. kbl(i)) THEN
650 rdelksd = (prsld(i, k)-prsld(i, k+1))*delks1(i) + (prsl(i, k)-&
651 & prsl(i, k+1))*delks1d(i)
652 rdelks = (prsl(i, k)-prsl(i, k+1))*delks1(i)
653 bnv2d(i, 1) = bnv2d(i, 1) + bnv2d(i, k)*rdelks + bnv2(i, k)*&
654 & rdelksd
655 bnv2(i, 1) = bnv2(i, 1) + bnv2(i, k)*rdelks
656 usqjd(i, 1) = usqjd(i, 1) + usqjd(i, k)*rdelks + usqj(i, k)*&
657 & rdelksd
658 usqj(i, 1) = usqj(i, 1) + usqj(i, k)*rdelks
659 END IF
660 END DO
661 END DO
662 !
663 DO i=its,ite
664 ldrag(i) = ldrag(i) .OR. bnv2(i, 1) .LE. 0.0
665 ldrag(i) = ldrag(i) .OR. ulow(i) .EQ. 1.0
666 ldrag(i) = ldrag(i) .OR. var(i) .LE. 0.0
667 END DO
668 !
669 ! set all ri low level values to the low level value
670 !
671 DO k=kpblmin,kpblmax
672 DO i=its,ite
673 IF (k .LT. kbl(i)) THEN
674 usqjd(i, k) = usqjd(i, 1)
675 usqj(i, k) = usqj(i, 1)
676 END IF
677 END DO
678 END DO
679 bnvd = 0.0
680 xnd = 0.0
681 ynd = 0.0
682 frd = 0.0
683 !
684 DO i=its,ite
685 IF (.NOT.ldrag(i)) THEN
686 IF (bnv2(i, 1) .EQ. 0.0) THEN
687 bnvd(i) = 0.0
688 ELSE
689 bnvd(i) = bnv2d(i, 1)/(2.0*SQRT(bnv2(i, 1)))
690 END IF
691 bnv(i) = SQRT(bnv2(i, 1))
692 frd(i) = 2.*var(i)*od(i)*(bnvd(i)*rulow(i)+bnv(i)*rulowd(i))
693 fr(i) = bnv(i)*rulow(i)*2.*var(i)*od(i)
694 IF (fr(i) .GT. frmax) THEN
695 frd(i) = 0.0
696 fr(i) = frmax
697 ELSE
698 fr(i) = fr(i)
699 END IF
700 xnd(i) = ubard(i)*rulow(i) + ubar(i)*rulowd(i)
701 xn(i) = ubar(i)*rulow(i)
702 ynd(i) = vbard(i)*rulow(i) + vbar(i)*rulowd(i)
703 yn(i) = vbar(i)*rulow(i)
704 END IF
705 END DO
706 taubd = 0.0
707 !
708 ! compute the base level stress and store it in taub
709 ! calculate enhancement factor, number of mountains & aspect
710 ! ratio const. use simplified relationship between standard
711 ! deviation & critical hgt
712 !
713 DO i=its,ite
714 IF (.NOT.ldrag(i)) THEN
715 pwx1 = oa(i) + 2.
716 pwy1d = ce*frd(i)/frc
717 pwy1 = ce*fr(i)/frc
718 IF (pwx1 .GT. 0.0) THEN
719 efactd = LOG(pwx1)*pwx1**pwy1*pwy1d
720 ELSE
721 efactd = 0.0
722 END IF
723 efact = pwx1**pwy1
724 IF (efact .LT. efmin) THEN
725 x3 = efmin
726 x3d = 0.0
727 ELSE
728 x3d = efactd
729 x3 = efact
730 END IF
731 IF (x3 .GT. efmax) THEN
732 efact = efmax
733 efactd = 0.0
734 ELSE
735 efactd = x3d
736 efact = x3
737 END IF
738 !!!!!!! cleff (effective grid length) is highly tunable parameter
739 !!!!!!! the bigger (smaller) value produce weaker (stronger) wave drag
740 arg10 = dxy(i)**2. + dxyp(i)**2.
741 cleff = SQRT(arg10)
742 IF (dxmeter .LT. cleff) THEN
743 max3 = cleff
744 ELSE
745 max3 = dxmeter
746 END IF
747 cleff = 3.*max3
748 pwx1 = 1. + ol(i)
749 pwy10 = oa(i) + 1.
750 coefm(i) = pwx1**pwy10
751 xlinv(i) = coefm(i)/cleff
752 temd = oc1(i)*(frd(i)*fr(i)+fr(i)*frd(i))
753 tem = fr(i)*fr(i)*oc1(i)
754 gfobnvd = (gmax*temd*(tem+cg)*bnv(i)-gmax*tem*(temd*bnv(i)+(tem+cg&
755 & )*bnvd(i)))/((tem+cg)*bnv(i))**2
756 gfobnv = gmax*tem/((tem+cg)*bnv(i))
757 taubd(i) = xlinv(i)*(((rolld(i)*ulow(i)+roll(i)*ulowd(i))*gfobnv*&
758 & efact+roll(i)*ulow(i)*(gfobnvd*efact+gfobnv*efactd))*ulow(i)**2+&
759 & roll(i)*ulow(i)*gfobnv*efact*(ulowd(i)*ulow(i)+ulow(i)*ulowd(i))&
760 & )
761 taub(i) = xlinv(i)*roll(i)*ulow(i)*ulow(i)*ulow(i)*gfobnv*efact
762 ELSE
763 taubd(i) = 0.0
764 taub(i) = 0.0
765 xnd(i) = 0.0
766 xn(i) = 0.0
767 ynd(i) = 0.0
768 yn(i) = 0.0
769 END IF
770 END DO
771 taupd = 0.0
772 !
773 ! now compute vertical structure of the stress.
774 !
775 DO k=kts,kpblmax
776 DO i=its,ite
777 IF (k .LE. kbl(i)) THEN
778 taupd(i, k) = taubd(i)
779 taup(i, k) = taub(i)
780 END IF
781 END DO
782 END DO
783 brvfd = 0.0
784 !
785 ! vertical level k loop!
786 DO k=kpblmin,kte-1
787 kp1 = k + 1
788 DO i=its,ite
789 !
790 ! unstablelayer if ri < ric
791 ! unstable layer if upper air vel comp along surf vel <=0 (crit lay)
792 ! at (u-c)=0. crit layer exists and bit vector should be set (.le.)
793 !
794 IF (k .GE. kbl(i)) THEN
795 icrilv(i) = (icrilv(i) .OR. usqj(i, k) .LT. ric) .OR. velco(i, k&
796 & ) .LE. 0.0
797 IF (bnv2(i, k) .LT. bnv2min) THEN
798 brvfd(i) = 0.0
799 brvf(i) = bnv2min
800 ELSE
801 brvfd(i) = bnv2d(i, k)
802 brvf(i) = bnv2(i, k)
803 END IF
804 ! brunt-vaisala frequency
805 IF (brvf(i) .EQ. 0.0) THEN
806 brvfd(i) = 0.0
807 ELSE
808 brvfd(i) = brvfd(i)/(2.0*SQRT(brvf(i)))
809 END IF
810 brvf(i) = SQRT(brvf(i))
811 END IF
812 END DO
813 !
814 DO i=its,ite
815 IF (k .GE. kbl(i) .AND. (.NOT.ldrag(i))) THEN
816 IF (.NOT.icrilv(i) .AND. taup(i, k) .GT. 0.0) THEN
817 temv = 1.0/velco(i, k)
818 tem1d = coefm(i)*0.5*((rod(i, kp1)+rod(i, k))*brvf(i)*velco(i&
819 & , k)+(ro(i, kp1)+ro(i, k))*(brvfd(i)*velco(i, k)+brvf(i)*&
820 & velcod(i, k)))/dxy(i)
821 tem1 = coefm(i)/dxy(i)*(ro(i, kp1)+ro(i, k))*brvf(i)*velco(i, &
822 & k)*0.5
823 hd = SQRT(taup(i, k)/tem1)
824 fro = brvf(i)*hd*temv
825 !
826 ! rim is the minimum-richardson number by shutts (1985)
827 !
828 IF (usqj(i, k) .EQ. 0.0) THEN
829 tem2d = 0.0
830 ELSE
831 tem2d = usqjd(i, k)/(2.0*SQRT(usqj(i, k)))
832 END IF
833 tem2 = SQRT(usqj(i, k))
834 tem = 1. + tem2*fro
835 rim = usqj(i, k)*(1.-fro)/(tem*tem)
836 !
837 ! check stability to employ the 'saturation hypothesis'
838 ! of lindzen (1981) except at tropospheric downstream regions
839 !
840 IF (rim .LE. ric) THEN
841 ! saturation hypothesis!
842 IF (oa(i) .LE. 0. .OR. kp1 .GE. kpblmin) THEN
843 temcd = -(tem2d/tem2**2)
844 temc = 2.0 + 1.0/tem2
845 IF (temc .EQ. 0.0) THEN
846 result1d = 0.0
847 ELSE
848 result1d = temcd/(2.0*SQRT(temc))
849 END IF
850 result1 = SQRT(temc)
851 hdd = ((velcod(i, k)*(2.*result1-temc)+velco(i, k)*(2.*&
852 & result1d-temcd))*brvf(i)-velco(i, k)*(2.*result1-temc)*&
853 & brvfd(i))/brvf(i)**2
854 hd = velco(i, k)*(2.*result1-temc)/brvf(i)
855 taupd(i, kp1) = (tem1d*hd+tem1*hdd)*hd + tem1*hd*hdd
856 taup(i, kp1) = tem1*hd*hd
857 END IF
858 ELSE
859 ! no wavebreaking!
860 taupd(i, kp1) = taupd(i, k)
861 taup(i, kp1) = taup(i, k)
862 END IF
863 END IF
864 END IF
865 END DO
866 END DO
867 !
868 IF (lcap .LT. kte) THEN
869 DO klcap=lcapp1,kte
870 DO i=its,ite
871 taupd(i, klcap) = (prsid(i, klcap)*prsi(i, lcap)-prsi(i, klcap)*&
872 & prsid(i, lcap))*taup(i, lcap)/prsi(i, lcap)**2 + prsi(i, klcap&
873 & )*taupd(i, lcap)/prsi(i, lcap)
874 taup(i, klcap) = prsi(i, klcap)/prsi(i, lcap)*taup(i, lcap)
875 END DO
876 END DO
877 zblkd = 0.0
878 taufbd = 0.0
879 ELSE
880 zblkd = 0.0
881 taufbd = 0.0
882 END IF
883 DO i=its,ite
884 IF (.NOT.ldrag(i)) THEN
885 !
886 !------- determine the height of flow-blocking layer
887 !
888 kblk = 0
889 pe = 0.0
890 DO k=kte,kpblmin,-1
891 IF (kblk .EQ. 0 .AND. k .LE. komax(i)) THEN
892 pe = pe + bnv2(i, k)*(zl(i, komax(i))-zl(i, k))*del(i, k)/g/ro&
893 & (i, k)
894 ke = 0.5*((rcs*u1(i, k))**2.+(rcs*v1(i, k))**2.)
895 !
896 !---------- apply flow-blocking drag when pe >= ke
897 !
898 IF (pe .GE. ke) THEN
899 kblk = k
900 IF (kblk .GT. kbl(i)) THEN
901 kblk = kbl(i)
902 ELSE
903 kblk = kblk
904 END IF
905 zblkd = zld(i, kblk) - zld(i, kts)
906 zblk = zl(i, kblk) - zl(i, kts)
907 END IF
908 END IF
909 END DO
910 IF (kblk .NE. 0) THEN
911 !
912 !--------- compute flow-blocking stress
913 !
914 IF (2.0 - 1.0/od(i) .LT. 0.0) THEN
915 cd = 0.0
916 ELSE
917 cd = 2.0 - 1.0/od(i)
918 END IF
919 taufbd(i, kts) = cd*dxyp(i)*olp(i)*(0.5*coefm(i)*rolld(i)*zblk*&
920 & ulow(i)**2/dxy(i)**2+0.5*roll(i)*coefm(i)*(zblkd*ulow(i)**2+&
921 & zblk*2*ulow(i)*ulowd(i))/dxy(i)**2)
922 taufb(i, kts) = 0.5*roll(i)*coefm(i)/dxy(i)**2*cd*dxyp(i)*olp(i)&
923 & *zblk*ulow(i)**2
924 tautemd = taufbd(i, kts)/FLOAT(kblk-kts)
925 tautem = taufb(i, kts)/FLOAT(kblk-kts)
926 DO k=kts+1,kblk
927 taufbd(i, k) = taufbd(i, k-1) - tautemd
928 taufb(i, k) = taufb(i, k-1) - tautem
929 END DO
930 !
931 !----------sum orographic GW stress and flow-blocking stress
932 !
933 taupd(i, :) = taupd(i, :) + taufbd(i, :)
934 taup(i, :) = taup(i, :) + taufb(i, :)
935 END IF
936 END IF
937 END DO
938 taudd = 0.0
939 !
940 ! calculate - (g)*d(tau)/d(pressure) and deceleration terms dtaux, dtauy
941 !
942 DO k=kts,kte
943 DO i=its,ite
944 taudd(i, k) = (csg*(taupd(i, k+1)-taupd(i, k))*del(i, k)-(taup(i, &
945 & k+1)-taup(i, k))*csg*deld(i, k))/del(i, k)**2
946 taud(i, k) = 1.*(taup(i, k+1)-taup(i, k))*csg/del(i, k)
947 END DO
948 END DO
949 !
950 ! limit de-acceleration (momentum deposition ) at top to 1/2 value
951 ! the idea is some stuff must go out the 'top'
952 !
953 DO klcap=lcap,kte
954 DO i=its,ite
955 taudd(i, klcap) = factop*taudd(i, klcap)
956 taud(i, klcap) = taud(i, klcap)*factop
957 END DO
958 END DO
959 dtfacd = 0.0
960 !
961 ! if the gravity wave drag would force a critical line
962 ! in the lower ksmm1 layers during the next deltim timestep,
963 ! then only apply drag until that critical line is reached.
964 !
965 DO k=kts,kpblmax-1
966 DO i=its,ite
967 IF (k .LE. kbl(i)) THEN
968 IF (taud(i, k) .NE. 0.) THEN
969 x4d = (velcod(i, k)*deltim*rcs*taud(i, k)-velco(i, k)*deltim*&
970 & rcs*taudd(i, k))/(deltim*rcs*taud(i, k))**2
971 x4 = velco(i, k)/(deltim*rcs*taud(i, k))
972 IF (x4 .GE. 0.) THEN
973 y1d = x4d
974 y1 = x4
975 ELSE
976 y1d = -x4d
977 y1 = -x4
978 END IF
979 IF (dtfac(i) .GT. y1) THEN
980 dtfacd(i) = y1d
981 dtfac(i) = y1
982 ELSE
983 dtfac(i) = dtfac(i)
984 END IF
985 END IF
986 END IF
987 END DO
988 END DO
989 !
990 DO i=its,ite
991 dusfcd(i) = 0.0
992 dusfc(i) = 0.
993 dvsfcd(i) = 0.0
994 dvsfc(i) = 0.
995 END DO
996 !
997 DO k=kts,kte
998 DO i=its,ite
999 taudd(i, k) = taudd(i, k)*dtfac(i) + taud(i, k)*dtfacd(i)
1000 taud(i, k) = taud(i, k)*dtfac(i)
1001 dtauxd = taudd(i, k)*xn(i) + taud(i, k)*xnd(i)
1002 dtaux = taud(i, k)*xn(i)
1003 dtauyd = taudd(i, k)*yn(i) + taud(i, k)*ynd(i)
1004 dtauy = taud(i, k)*yn(i)
1005 dtaux2dd(i, k) = dtauxd
1006 dtaux2d(i, k) = dtaux
1007 dtauy2dd(i, k) = dtauyd
1008 dtauy2d(i, k) = dtauy
1009 dudtd(i, k) = dtauxd + dudtd(i, k)
1010 dudt(i, k) = dtaux + dudt(i, k)
1011 dvdtd(i, k) = dtauyd + dvdtd(i, k)
1012 dvdt(i, k) = dtauy + dvdt(i, k)
1013 dusfcd(i) = dusfcd(i) + dtauxd*del(i, k) + dtaux*deld(i, k)
1014 dusfc(i) = dusfc(i) + dtaux*del(i, k)
1015 dvsfcd(i) = dvsfcd(i) + dtauyd*del(i, k) + dtauy*deld(i, k)
1016 dvsfc(i) = dvsfc(i) + dtauy*del(i, k)
1017 END DO
1018 END DO
1019 !
1020 DO i=its,ite
1021 dusfcd(i) = -(rcs*dusfcd(i)/g)
1022 dusfc(i) = (-(1./g*rcs))*dusfc(i)
1023 dvsfcd(i) = -(rcs*dvsfcd(i)/g)
1024 dvsfc(i) = (-(1./g*rcs))*dvsfc(i)
1025 END DO
1026 !
1027 RETURN
1028 END SUBROUTINE GWDO2D_D
1030 end module g_module_bl_gwdo