module_diffusion_em.F
References to this file elsewhere.
1 ! WRF:MODEL_LAYER:PHYSICS
2
3 MODULE module_diffusion_em
4
5 USE module_configure
6 USE module_bc
7 USE module_state_description
8 USE module_big_step_utilities_em
9 USE module_model_constants
10 USE module_wrf_error
11
12 CONTAINS
13
14 !=======================================================================
15 !=======================================================================
16
17 SUBROUTINE cal_deform_and_div( config_flags, u, v, w, div, &
18 defor11, defor22, defor33, &
19 defor12, defor13, defor23, &
20 u_base, v_base, msfux, msfuy, &
21 msfvx, msfvy, msftx, msfty, &
22 rdx, rdy, dn, dnw, rdz, rdzw, &
23 fnm, fnp, cf1, cf2, cf3, zx, zy, &
24 ids, ide, jds, jde, kds, kde, &
25 ims, ime, jms, jme, kms, kme, &
26 its, ite, jts, jte, kts, kte )
27
28 ! History: Sep 2003 Changes by Jason Knievel and George Bryan, NCAR
29 ! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
30 ! ... ...
31
32 ! Purpose: This routine calculates deformation and 3-d divergence.
33
34 ! References: Klemp and Wilhelmson (JAS 1978)
35 ! Chen and Dudhia (NCAR WRF physics report 2000)
36
37 !-----------------------------------------------------------------------
38 ! Comments 10-MAR-05
39 ! Equations 13a-f, Chen and Dudhia 2000, Appendix A:
40 ! Eqn 13a: D11=defor11= 2m^2 * (partial du^/dX + partial dpsi/dx * partial du^/dpsi)
41 ! Eqn 13b: D22=defor22= 2m^2 * (partial dv^/dY + partial dpsi/dy * partial dv^/dpsi)
42 ! Eqn 13c: D33=defor33= 2 * partial dw/dz [SIMPLER FORM]
43 ! Eqn 13d: D12=defor12= m^2 * (partial dv^/dX + partial du^/dY +
44 ! partial dpsi/dx * partial dv^/dpsi +
45 ! partial dpsi/dy * partial du^/dpsi)
46 ! Eqn 13e: D13=defor13= m^2 * (partial dw^/dX + partial dpsi/dx * partial dw^/dpsi)
47 ! + partial du/dz [SIMPLER FORM]
48 ! Eqn 13f: D23=defor23= m^2 * (partial dw^/dY + partial dpsi/dy * partial dw^/dpsi)
49 ! + partial dv/dz [SIMPLER FORM]
50 !-----------------------------------------------------------------------
51 ! Begin declarations.
52
53 IMPLICIT NONE
54
55 TYPE( grid_config_rec_type ), INTENT( IN ) &
56 :: config_flags
57
58 INTEGER, INTENT( IN ) &
59 :: ids, ide, jds, jde, kds, kde, &
60 ims, ime, jms, jme, kms, kme, &
61 its, ite, jts, jte, kts, kte
62
63 REAL, INTENT( IN ) &
64 :: rdx, rdy, cf1, cf2, cf3
65
66 REAL, DIMENSION( kms:kme ), INTENT( IN ) &
67 :: fnm, fnp, dn, dnw, u_base, v_base
68
69 REAL, DIMENSION( ims:ime , jms:jme ), INTENT( IN ) &
70 :: msfux, msfuy, msfvx, msfvy, msftx, msfty
71
72 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
73 :: u, v, w, zx, zy, rdz, rdzw
74
75 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
76 :: defor11, defor22, defor33, defor12, defor13, defor23, div
77
78 ! Local variables.
79
80 INTEGER &
81 :: i, j, k, ktf, ktes1, ktes2, i_start, i_end, j_start, j_end
82
83 REAL &
84 :: tmp, tmpzx, tmpzy, tmpzeta_z, cft1, cft2
85
86 REAL, DIMENSION( its:ite, jts:jte ) &
87 :: mm, zzavg, zeta_zd12
88
89 REAL, DIMENSION( its-2:ite+2, kts:kte, jts-2:jte+2 ) &
90 :: tmp1, hat, hatavg
91
92 ! End declarations.
93 !-----------------------------------------------------------------------
94
95 ! Comments 10-MAR-2005
96 ! Treat all differentials as 'div-style' [or 'curl-style'],
97 ! i.e., du/dX becomes (in map coordinate space) mx*my * d(u/my)/dx,
98 ! NB - all equations referred to here are from Chen and Dudhia 2002, from the
99 ! WRF physics documents web pages:
100 ! http://www.mmm.ucar.edu/wrf/users/docs/wrf-doc-physics.pdf
101
102 !=======================================================================
103 ! In the following section, calculate 3-d divergence and the first three
104 ! (defor11, defor22, defor33) of six deformation terms.
105
106 ktes1 = kte-1
107 ktes2 = kte-2
108
109 cft2 = - 0.5 * dnw(ktes1) / dn(ktes1)
110 cft1 = 1.0 - cft2
111
112 ktf = MIN( kte, kde-1 )
113
114 i_start = its
115 i_end = MIN( ite, ide-1 )
116 j_start = jts
117 j_end = MIN( jte, jde-1 )
118
119 ! Square the map scale factor.
120
121 DO j = j_start, j_end
122 DO i = i_start, i_end
123 mm(i,j) = msftx(i,j) * msfty(i,j)
124 END DO
125 END DO
126
127 !-----------------------------------------------------------------------
128 ! Calculate du/dx.
129
130 ! Apply a coordinate transformation to zonal velocity, u.
131
132 DO j = j_start, j_end
133 DO k = kts, ktf
134 DO i = i_start, i_end+1
135 hat(i,k,j) = u(i,k,j) / msfuy(i,j)
136 END DO
137 END DO
138 END DO
139
140 ! Average in x and z.
141
142 DO j=j_start,j_end
143 DO k=kts+1,ktf
144 DO i=i_start,i_end
145 hatavg(i,k,j) = 0.5 * &
146 ( fnm(k) * ( hat(i,k ,j) + hat(i+1, k,j) ) + &
147 fnp(k) * ( hat(i,k-1,j) + hat(i+1,k-1,j) ) )
148 END DO
149 END DO
150 END DO
151
152 ! Extrapolate to top and bottom of domain (to w levels).
153
154 DO j = j_start, j_end
155 DO i = i_start, i_end
156 hatavg(i,1,j) = 0.5 * ( &
157 cf1 * hat(i ,1,j) + &
158 cf2 * hat(i ,2,j) + &
159 cf3 * hat(i ,3,j) + &
160 cf1 * hat(i+1,1,j) + &
161 cf2 * hat(i+1,2,j) + &
162 cf3 * hat(i+1,3,j) )
163 hatavg(i,kte,j) = 0.5 * ( &
164 cft1 * ( hat(i,ktes1,j) + hat(i+1,ktes1,j) ) + &
165 cft2 * ( hat(i,ktes2,j) + hat(i+1,ktes2,j) ) )
166 END DO
167 END DO
168
169 ! Comments 10-MAR-05
170 ! Eqn 13a: D11=defor11= 2m^2 * (partial du^/dX + partial dpsi/dx * partial du^/dpsi)
171 ! Below, D11 is set = 2*tmp1
172 ! => tmp1 = m^2 * (partial du^/dX + partial dpsi/dx * partial du^/dpsi)
173 ! tmpzx = averaged value of dpsi/dx (=zx)
174
175 DO j = j_start, j_end
176 DO k = kts, ktf
177 DO i = i_start, i_end
178 tmpzx = 0.25 * ( &
179 zx(i,k ,j) + zx(i+1,k ,j) + &
180 zx(i,k+1,j) + zx(i+1,k+1,j) )
181 tmp1(i,k,j) = ( hatavg(i,k+1,j) - hatavg(i,k,j) ) *tmpzx * rdzw(i,k,j)
182 ! tmp1 to here = partial dpsi/dx * partial du^/dpsi:
183 END DO
184 END DO
185 END DO
186
187 DO j = j_start, j_end
188 DO k = kts, ktf
189 DO i = i_start, i_end
190 tmp1(i,k,j) = mm(i,j) * ( rdx * ( hat(i+1,k,j) - hat(i,k,j) ) - &
191 tmp1(i,k,j))
192 END DO
193 END DO
194 END DO
195
196 ! End calculation of du/dx.
197 !-----------------------------------------------------------------------
198
199 !-----------------------------------------------------------------------
200 ! Calculate defor11 (2*du/dx).
201 ! Comments 10-MAR-05
202 ! Eqn 13a: D11=defor11= 2 m^2 * (partial du^/dX + partial dpsi/dx * partial du^/dpsi)
203 ! = 2*tmp1
204
205 DO j = j_start, j_end
206 DO k = kts, ktf
207 DO i = i_start, i_end
208 defor11(i,k,j) = 2.0 * tmp1(i,k,j)
209 END DO
210 END DO
211 END DO
212
213 ! End calculation of defor11.
214 !-----------------------------------------------------------------------
215
216 !-----------------------------------------------------------------------
217 ! Calculate zonal divergence (du/dx) and add it to the divergence array.
218
219 DO j = j_start, j_end
220 DO k = kts, ktf
221 DO i = i_start, i_end
222 div(i,k,j) = tmp1(i,k,j)
223 END DO
224 END DO
225 END DO
226
227 ! End calculation of zonal divergence.
228 !-----------------------------------------------------------------------
229
230 !-----------------------------------------------------------------------
231 ! Calculate dv/dy.
232
233 ! Apply a coordinate transformation to meridional velocity, v.
234
235 DO j = j_start, j_end+1
236 DO k = kts, ktf
237 DO i = i_start, i_end
238 ! Because msfvx at the poles will be undefined (1./0.), we will have
239 ! trouble. But we are OK since v at the poles is 0., and that takes
240 ! precedence in this case.
241 IF ((config_flags%polar) .AND. ((j == jds) .OR. (j == jde))) THEN
242 hat(i,k,j) = 0.
243 ELSE ! normal code
244 hat(i,k,j) = v(i,k,j) / msfvx(i,j)
245 ENDIF
246 END DO
247 END DO
248 END DO
249
250 ! Account for the slope in y of eta surfaces.
251
252 DO j=j_start,j_end
253 DO k=kts+1,ktf
254 DO i=i_start,i_end
255 hatavg(i,k,j) = 0.5 * ( &
256 fnm(k) * ( hat(i,k ,j) + hat(i,k ,j+1) ) + &
257 fnp(k) * ( hat(i,k-1,j) + hat(i,k-1,j+1) ) )
258 END DO
259 END DO
260 END DO
261
262 ! Extrapolate to top and bottom of domain (to w levels).
263
264 DO j = j_start, j_end
265 DO i = i_start, i_end
266 hatavg(i,1,j) = 0.5 * ( &
267 cf1 * hat(i,1,j ) + &
268 cf2 * hat(i,2,j ) + &
269 cf3 * hat(i,3,j ) + &
270 cf1 * hat(i,1,j+1) + &
271 cf2 * hat(i,2,j+1) + &
272 cf3 * hat(i,3,j+1) )
273 hatavg(i,kte,j) = 0.5 * ( &
274 cft1 * ( hat(i,ktes1,j) + hat(i,ktes1,j+1) ) + &
275 cft2 * ( hat(i,ktes2,j) + hat(i,ktes2,j+1) ) )
276 END DO
277 END DO
278
279 ! Comments 10-MAR-05
280 ! Eqn 13b: D22=defor22= 2m^2 * (partial dv^/dY + partial dpsi/dy * partial dv^/dpsi)
281 ! Below, D22 is set = 2*tmp1
282 ! => tmp1 = m^2 * (partial dv^/dY + partial dpsi/dy * partial dv^/dpsi)
283 ! tmpzy = averaged value of dpsi/dy (=zy)
284
285 DO j = j_start, j_end
286 DO k = kts, ktf
287 DO i = i_start, i_end
288 tmpzy = 0.25 * ( &
289 zy(i,k ,j) + zy(i,k ,j+1) + &
290 zy(i,k+1,j) + zy(i,k+1,j+1) )
291 tmp1(i,k,j) = ( hatavg(i,k+1,j) - hatavg(i,k,j) ) * tmpzy * rdzw(i,k,j)
292 ! tmp1 to here = partial dpsi/dy * partial dv^/dpsi:
293 END DO
294 END DO
295 END DO
296
297 DO j = j_start, j_end
298 DO k = kts, ktf
299 DO i = i_start, i_end
300 tmp1(i,k,j) = mm(i,j) * ( &
301 rdy * ( hat(i,k,j+1) - hat(i,k,j) ) - tmp1(i,k,j) )
302 END DO
303 END DO
304 END DO
305
306 ! End calculation of dv/dy.
307 !-----------------------------------------------------------------------
308
309 !-----------------------------------------------------------------------
310 ! Calculate defor22 (2*dv/dy).
311 ! Comments 10-MAR-05
312 ! Eqn 13b: D22=defor22= 2 m^2 * (partial dv^/dY + partial dpsi/dy * partial dv^/dpsi)
313 ! = 2*tmp1
314
315 DO j = j_start, j_end
316 DO k = kts, ktf
317 DO i = i_start, i_end
318 defor22(i,k,j) = 2.0 * tmp1(i,k,j)
319 END DO
320 END DO
321 END DO
322
323 ! End calculation of defor22.
324 !-----------------------------------------------------------------------
325
326 !-----------------------------------------------------------------------
327 ! Calculate meridional divergence (dv/dy) and add it to the divergence
328 ! array.
329
330 DO j = j_start, j_end
331 DO k = kts, ktf
332 DO i = i_start, i_end
333 div(i,k,j) = div(i,k,j) + tmp1(i,k,j)
334 END DO
335 END DO
336 END DO
337
338 ! End calculation of meridional divergence.
339 !-----------------------------------------------------------------------
340
341 !-----------------------------------------------------------------------
342 ! Comments 10-MAR-05
343 ! Eqn 13c: D33=defor33= 2 * partial dw/dz
344 ! Below, D33 is set = 2*tmp1
345 ! => tmp1 = partial dw/dz
346
347 ! Calculate dw/dz.
348
349 DO j = j_start, j_end
350 DO k = kts, ktf
351 DO i = i_start, i_end
352 tmp1(i,k,j) = ( w(i,k+1,j) - w(i,k,j) ) * rdzw(i,k,j)
353 END DO
354 END DO
355 END DO
356
357 ! End calculation of dw/dz.
358 !-----------------------------------------------------------------------
359
360 !-----------------------------------------------------------------------
361 ! Calculate defor33 (2*dw/dz).
362
363 DO j = j_start, j_end
364 DO k = kts, ktf
365 DO i = i_start, i_end
366 defor33(i,k,j) = 2.0 * tmp1(i,k,j)
367 END DO
368 END DO
369 END DO
370
371 ! End calculation of defor33.
372 !-----------------------------------------------------------------------
373
374 !-----------------------------------------------------------------------
375 ! Calculate vertical divergence (dw/dz) and add it to the divergence
376 ! array.
377
378 DO j = j_start, j_end
379 DO k = kts, ktf
380 DO i = i_start, i_end
381 div(i,k,j) = div(i,k,j) + tmp1(i,k,j)
382 END DO
383 END DO
384 END DO
385
386 ! End calculation of vertical divergence.
387 !-----------------------------------------------------------------------
388
389 ! Three-dimensional divergence is now finished and values are in array
390 ! "div." Also, the first three (defor11, defor22, defor33) of six
391 ! deformation terms are now calculated at pressure points.
392 !=======================================================================
393
394 ! Comments 10-MAR-2005
395 ! Treat all differentials as 'div-style' [or 'curl-style'],
396 ! i.e., du/dY becomes (in map coordinate space) mx*my * d(u/mx)/dy,
397 ! dv/dX becomes (in map coordinate space) mx*my * d(v/my)/dx,
398 ! (see e.g. Haltiner and Williams p. 441)
399
400 !=======================================================================
401 ! Calculate the final three deformations (defor12, defor13, defor23) at
402 ! vorticity points.
403
404 i_start = its
405 i_end = ite
406 j_start = jts
407 j_end = jte
408
409 IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
410 config_flags%nested) i_start = MAX( ids+1, its )
411 IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
412 config_flags%nested) i_end = MIN( ide-1, ite )
413 IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
414 config_flags%nested) j_start = MAX( jds+1, jts )
415 IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
416 config_flags%nested) j_end = MIN( jde-1, jte )
417 IF ( config_flags%periodic_x ) i_start = its
418 IF ( config_flags%periodic_x ) i_end = ite
419
420
421 !-----------------------------------------------------------------------
422 ! Calculate du/dy.
423
424 ! First, calculate an average mapscale factor.
425
426 ! Comments 10-MAR-05
427 ! du/dy => need u map scale factor in x (which is defined at u points)
428 ! averaged over j and j-1
429 ! dv/dx => need v map scale factor in y (which is defined at v points)
430 ! averaged over i and i-1
431
432 DO j = j_start, j_end
433 DO i = i_start, i_end
434 mm(i,j) = 0.25 * ( msfux(i,j-1) + msfux(i,j) ) * ( msfvy(i-1,j) + msfvy(i,j) )
435 END DO
436 END DO
437
438 ! Apply a coordinate transformation to zonal velocity, u.
439
440 DO j =j_start-1, j_end
441 DO k =kts, ktf
442 DO i =i_start, i_end
443 ! Fixes to set_physical_bc2/3d for polar boundary conditions
444 ! remove issues with loop over j
445 hat(i,k,j) = u(i,k,j) / msfux(i,j)
446 END DO
447 END DO
448 END DO
449
450 ! Average in y and z.
451
452 DO j=j_start,j_end
453 DO k=kts+1,ktf
454 DO i=i_start,i_end
455 hatavg(i,k,j) = 0.5 * ( &
456 fnm(k) * ( hat(i,k ,j-1) + hat(i,k ,j) ) + &
457 fnp(k) * ( hat(i,k-1,j-1) + hat(i,k-1,j) ) )
458 END DO
459 END DO
460 END DO
461
462 ! Extrapolate to top and bottom of domain (to w levels).
463
464 DO j = j_start, j_end
465 DO i = i_start, i_end
466 hatavg(i,1,j) = 0.5 * ( &
467 cf1 * hat(i,1,j-1) + &
468 cf2 * hat(i,2,j-1) + &
469 cf3 * hat(i,3,j-1) + &
470 cf1 * hat(i,1,j ) + &
471 cf2 * hat(i,2,j ) + &
472 cf3 * hat(i,3,j ) )
473 hatavg(i,kte,j) = 0.5 * ( &
474 cft1 * ( hat(i,ktes1,j-1) + hat(i,ktes1,j) ) + &
475 cft2 * ( hat(i,ktes2,j-1) + hat(i,ktes2,j) ) )
476 END DO
477 END DO
478
479 ! tmpzy = averaged value of dpsi/dy (=zy) on vorticity grid
480 ! tmp1 = partial dpsi/dy * partial du^/dpsi
481 DO j = j_start, j_end
482 DO k = kts, ktf
483 DO i = i_start, i_end
484 tmpzy = 0.25 * ( &
485 zy(i-1,k ,j) + zy(i,k ,j) + &
486 zy(i-1,k+1,j) + zy(i,k+1,j) )
487 tmp1(i,k,j) = ( hatavg(i,k+1,j) - hatavg(i,k,j) ) * &
488 0.25 * tmpzy * ( rdzw(i,k,j) + rdzw(i-1,k,j) + &
489 rdzw(i-1,k,j-1) + rdzw(i,k,j-1) )
490 END DO
491 END DO
492 END DO
493
494 ! End calculation of du/dy.
495 !----------------------------------------------------------------------
496
497 !-----------------------------------------------------------------------
498 ! Add the first term to defor12 (du/dy+dv/dx) at vorticity points.
499
500 ! Comments 10-MAR-05
501 ! Eqn 13d: D12=defor12= m^2 * (partial dv^/dX + partial du^/dY +
502 ! partial dpsi/dx * partial dv^/dpsi +
503 ! partial dpsi/dy * partial du^/dpsi)
504 ! Here deal with m^2 * (partial du^/dY + partial dpsi/dy * partial du^/dpsi)
505 ! Still need to add v^ terms:
506 ! m^2 * (partial dv^/dX + partial dpsi/dx * partial dv^/dpsi)
507
508 DO j = j_start, j_end
509 DO k = kts, ktf
510 DO i = i_start, i_end
511 defor12(i,k,j) = mm(i,j) * ( &
512 rdy * ( hat(i,k,j) - hat(i,k,j-1) ) - tmp1(i,k,j) )
513 END DO
514 END DO
515 END DO
516
517 ! End addition of the first term to defor12.
518 !-----------------------------------------------------------------------
519
520 !-----------------------------------------------------------------------
521 ! Calculate dv/dx.
522
523 ! Apply a coordinate transformation to meridional velocity, v.
524
525 DO j = j_start, j_end
526 DO k = kts, ktf
527 DO i = i_start-1, i_end
528 hat(i,k,j) = v(i,k,j) / msfvy(i,j)
529 END DO
530 END DO
531 END DO
532
533 ! Account for the slope in x of eta surfaces.
534
535 DO j = j_start, j_end
536 DO k = kts+1, ktf
537 DO i = i_start, i_end
538 hatavg(i,k,j) = 0.5 * ( &
539 fnm(k) * ( hat(i-1,k ,j) + hat(i,k ,j) ) + &
540 fnp(k) * ( hat(i-1,k-1,j) + hat(i,k-1,j) ) )
541 END DO
542 END DO
543 END DO
544
545 ! Extrapolate to top and bottom of domain (to w levels).
546
547 DO j = j_start, j_end
548 DO i = i_start, i_end
549 hatavg(i,1,j) = 0.5 * ( &
550 cf1 * hat(i-1,1,j) + &
551 cf2 * hat(i-1,2,j) + &
552 cf3 * hat(i-1,3,j) + &
553 cf1 * hat(i ,1,j) + &
554 cf2 * hat(i ,2,j) + &
555 cf3 * hat(i ,3,j) )
556 hatavg(i,kte,j) = 0.5 * ( &
557 cft1 * ( hat(i,ktes1,j) + hat(i-1,ktes1,j) ) + &
558 cft2 * ( hat(i,ktes2,j) + hat(i-1,ktes2,j) ) )
559 END DO
560 END DO
561
562 ! Fixes to set_physical_bc2/3d have made any check for polar B.C.'s
563 ! unnecessary in this place. zx, rdzw, and hatavg are all defined
564 ! in places they need to be and the values at the poles are replications
565 ! of the values one grid point in, so the averaging over j and j-1 works
566 ! to act as just using the value at j or j-1 (with out extra code).
567 !
568 ! tmpzx = averaged value of dpsi/dx (=zx) on vorticity grid
569 ! tmp1 = partial dpsi/dx * partial dv^/dpsi
570 DO j = j_start, j_end
571 DO k = kts, ktf
572 DO i = i_start, i_end
573 tmpzx = 0.25 * ( &
574 zx(i,k ,j-1) + zx(i,k ,j) + &
575 zx(i,k+1,j-1) + zx(i,k+1,j) )
576 tmp1(i,k,j) = ( hatavg(i,k+1,j) - hatavg(i,k,j) ) * &
577 0.25 * tmpzx * ( rdzw(i,k,j) + rdzw(i,k,j-1) + &
578 rdzw(i-1,k,j-1) + rdzw(i-1,k,j) )
579 END DO
580 END DO
581 END DO
582
583 ! End calculation of dv/dx.
584 !-----------------------------------------------------------------------
585
586 !-----------------------------------------------------------------------
587 ! Add the second term to defor12 (du/dy+dv/dx) at vorticity points.
588
589 ! Comments 10-MAR-05
590 ! Eqn 13d: D12=defor12= m^2 * (partial dv^/dX + partial du^/dY +
591 ! partial dpsi/dx * partial dv^/dpsi +
592 ! partial dpsi/dy * partial du^/dpsi)
593 ! Here adding v^ terms:
594 ! m^2 * (partial dv^/dX + partial dpsi/dx * partial dv^/dpsi)
595
596 DO j = j_start, j_end
597 DO k = kts, ktf
598 DO i = i_start, i_end
599 defor12(i,k,j) = defor12(i,k,j) + &
600 mm(i,j) * ( &
601 rdx * ( hat(i,k,j) - hat(i-1,k,j) ) - tmp1(i,k,j) )
602 END DO
603 END DO
604 END DO
605
606 ! End addition of the second term to defor12.
607 !-----------------------------------------------------------------------
608
609 !-----------------------------------------------------------------------
610 ! Update the boundary for defor12 (might need to change later).
611
612 IF ( .NOT. config_flags%periodic_x .AND. i_start .EQ. ids+1 ) THEN
613 DO j = jts, jte
614 DO k = kts, kte
615 defor12(ids,k,j) = defor12(ids+1,k,j)
616 END DO
617 END DO
618 END IF
619
620 IF ( .NOT. config_flags%periodic_y .AND. j_start .EQ. jds+1) THEN
621 DO k = kts, kte
622 DO i = its, ite
623 defor12(i,k,jds) = defor12(i,k,jds+1)
624 END DO
625 END DO
626 END IF
627
628 IF ( .NOT. config_flags%periodic_x .AND. i_end .EQ. ide-1) THEN
629 DO j = jts, jte
630 DO k = kts, kte
631 defor12(ide,k,j) = defor12(ide-1,k,j)
632 END DO
633 END DO
634 END IF
635
636 IF ( .NOT. config_flags%periodic_y .AND. j_end .EQ. jde-1) THEN
637 DO k = kts, kte
638 DO i = its, ite
639 defor12(i,k,jde) = defor12(i,k,jde-1)
640 END DO
641 END DO
642 END IF
643
644 ! End update of boundary for defor12.
645 !-----------------------------------------------------------------------
646
647 ! Comments 10-MAR-05
648 ! Further deformation terms not needed for 2-dimensional Smagorinsky diffusion,
649 ! so those terms have not been dealt with yet.
650 ! A "y" has simply been added to all map scale factors to allow the model to
651 ! compile without errors.
652
653 !-----------------------------------------------------------------------
654 ! Calculate dw/dx.
655
656 if ((config_flags%diff_opt.eq.2).or.(config_flags%km_opt.eq.2).or.(config_flags%km_opt.eq.3)) then
657
658 i_start = its
659 i_end = MIN( ite, ide-1 )
660 j_start = jts
661 j_end = MIN( jte, jde-1 )
662
663 IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
664 config_flags%nested) i_start = MAX( ids+1, its )
665 IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
666 config_flags%nested) j_start = MAX( jds+1, jts )
667
668 IF ( config_flags%periodic_x ) i_start = its
669 IF ( config_flags%periodic_x ) i_end = MIN( ite, ide )
670 IF ( config_flags%periodic_y ) j_end = MIN( jte, jde )
671
672 ! Square the mapscale factor.
673
674 DO j = jts, jte
675 DO i = its, ite
676 mm(i,j) = msfuy(i,j) * msfuy(i,j)
677 END DO
678 END DO
679
680 ! Apply a coordinate transformation to vertical velocity, w. This is for both
681 ! defor13 and defor23.
682
683 DO j = j_start, j_end
684 DO k = kts, kte
685 DO i = i_start, i_end
686 hat(i,k,j) = w(i,k,j) / msfty(i,j)
687 END DO
688 END DO
689 END DO
690
691 i = i_start-1
692 DO j = j_start, MIN( jte, jde-1 )
693 DO k = kts, kte
694 hat(i,k,j) = w(i,k,j) / msfty(i,j)
695 END DO
696 END DO
697
698 j = j_start-1
699 DO k = kts, kte
700 DO i = i_start, MIN( ite, ide-1 )
701 hat(i,k,j) = w(i,k,j) / msfty(i,j)
702 END DO
703 END DO
704
705 ! QUESTION: What is this for?
706
707 DO j = j_start, j_end
708 DO k = kts, ktf
709 DO i = i_start, i_end
710 hatavg(i,k,j) = 0.25 * ( &
711 hat(i ,k ,j) + &
712 hat(i ,k+1,j) + &
713 hat(i-1,k ,j) + &
714 hat(i-1,k+1,j) )
715 END DO
716 END DO
717 END DO
718
719 ! Calculate dw/dx.
720
721 DO j = j_start, j_end
722 DO k = kts+1, ktf
723 DO i = i_start, i_end
724 tmp1(i,k,j) = ( hatavg(i,k,j) - hatavg(i,k-1,j) ) * zx(i,k,j) * &
725 0.5 * ( rdz(i,k,j) + rdz(i-1,k,j) )
726 END DO
727 END DO
728 END DO
729
730 ! End calculation of dw/dx.
731 !-----------------------------------------------------------------------
732
733 !-----------------------------------------------------------------------
734 ! Add the first term (dw/dx) to defor13 (dw/dx+du/dz) at vorticity
735 ! points.
736
737 DO j = j_start, j_end
738 DO k = kts+1, ktf
739 DO i = i_start, i_end
740 defor13(i,k,j) = mm(i,j) * ( &
741 rdx * ( hat(i,k,j) - hat(i-1,k,j) ) - tmp1(i,k,j) )
742 END DO
743 END DO
744 END DO
745
746 DO j = j_start, j_end
747 DO i = i_start, i_end
748 defor13(i,kts,j ) = 0.0
749 defor13(i,ktf+1,j) = 0.0
750 END DO
751 END DO
752
753 ! End addition of the first term to defor13.
754 !-----------------------------------------------------------------------
755
756 !-----------------------------------------------------------------------
757 ! Calculate du/dz.
758
759 IF ( config_flags%mix_full_fields ) THEN
760
761 DO j = j_start, j_end
762 DO k = kts+1, ktf
763 DO i = i_start, i_end
764 tmp1(i,k,j) = ( u(i,k,j) - u(i,k-1,j) ) * &
765 0.5 * ( rdz(i,k,j) + rdz(i-1,k,j) )
766 END DO
767 END DO
768 END DO
769
770 ELSE
771
772 DO j = j_start, j_end
773 DO k = kts+1, ktf
774 DO i = i_start, i_end
775 tmp1(i,k,j) = ( u(i,k,j) - u_base(k) - u(i,k-1,j) + u_base(k-1) ) * &
776 0.5 * ( rdz(i,k,j) + rdz(i-1,k,j) )
777 END DO
778 END DO
779 END DO
780
781 END IF
782
783 !-----------------------------------------------------------------------
784 ! Add the second term (du/dz) to defor13 (dw/dx+du/dz) at vorticity
785 ! points.
786
787 DO j = j_start, j_end
788 DO k = kts+1, ktf
789 DO i = i_start, i_end
790 defor13(i,k,j) = defor13(i,k,j) + tmp1(i,k,j)
791 END DO
792 END DO
793 END DO
794
795 ! End addition of the second term to defor13.
796 !-----------------------------------------------------------------------
797
798 !-----------------------------------------------------------------------
799 ! Calculate dw/dy.
800
801 i_start = its
802 i_end = MIN( ite, ide-1 )
803 j_start = jts
804 j_end = MIN( jte, jde-1 )
805
806 IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
807 config_flags%nested) i_start = MAX( ids+1, its )
808 IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
809 config_flags%nested) j_start = MAX( jds+1, jts )
810 IF ( config_flags%periodic_y ) j_end = MIN( jte, jde )
811 IF ( config_flags%periodic_x ) i_start = its
812 IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
813
814 ! Square mapscale factor.
815
816 DO j = jts, jte
817 DO i = its, ite
818 mm(i,j) = msfvy(i,j) * msfvy(i,j)
819 END DO
820 END DO
821
822 ! QUESTION: What is this for?
823
824 DO j = j_start, j_end
825 DO k = kts, ktf
826 DO i = i_start, i_end
827 hatavg(i,k,j) = 0.25 * ( &
828 hat(i,k ,j ) + &
829 hat(i,k+1,j ) + &
830 hat(i,k ,j-1) + &
831 hat(i,k+1,j-1) )
832 END DO
833 END DO
834 END DO
835
836 ! Calculate dw/dy and store in tmp1.
837
838 DO j = j_start, j_end
839 DO k = kts+1, ktf
840 DO i = i_start, i_end
841 tmp1(i,k,j) = ( hatavg(i,k,j) - hatavg(i,k-1,j) ) * zy(i,k,j) * &
842 0.5 * ( rdz(i,k,j) + rdz(i,k,j-1) )
843 END DO
844 END DO
845 END DO
846
847 ! End calculation of dw/dy.
848 !-----------------------------------------------------------------------
849
850 !-----------------------------------------------------------------------
851 ! Add the first term (dw/dy) to defor23 (dw/dy+dv/dz) at vorticity
852 ! points.
853
854 DO j = j_start, j_end
855 DO k = kts+1, ktf
856 DO i = i_start, i_end
857 defor23(i,k,j) = mm(i,j) * ( &
858 rdy * ( hat(i,k,j) - hat(i,k,j-1) ) - tmp1(i,k,j) )
859 END DO
860 END DO
861 END DO
862
863 DO j = j_start, j_end
864 DO i = i_start, i_end
865 defor23(i,kts,j ) = 0.0
866 defor23(i,ktf+1,j) = 0.0
867 END DO
868 END DO
869
870 ! End addition of the first term to defor23.
871 !-----------------------------------------------------------------------
872
873 !-----------------------------------------------------------------------
874 ! Calculate dv/dz.
875
876 IF ( config_flags%mix_full_fields ) THEN
877
878 DO j = j_start, j_end
879 DO k = kts+1, ktf
880 DO i = i_start, i_end
881 tmp1(i,k,j) = ( v(i,k,j) - v(i,k-1,j) ) * &
882 0.5 * ( rdz(i,k,j) + rdz(i,k,j-1) )
883 END DO
884 END DO
885 END DO
886
887 ELSE
888
889 DO j = j_start, j_end
890 DO k = kts+1, ktf
891 DO i = i_start, i_end
892 tmp1(i,k,j) = ( v(i,k,j) - v_base(k) - v(i,k-1,j) + v_base(k-1) ) * &
893 0.5 * ( rdz(i,k,j) + rdz(i,k,j-1) )
894 END DO
895 END DO
896 END DO
897
898 END IF
899
900 ! End calculation of dv/dz.
901 !-----------------------------------------------------------------------
902
903 !-----------------------------------------------------------------------
904 ! Add the second term (dv/dz) to defor23 (dw/dy+dv/dz) at vorticity
905 ! points.
906
907 ! Add tmp1 to defor23.
908
909 DO j = j_start, j_end
910 DO k = kts+1, ktf
911 DO i = i_start, i_end
912 defor23(i,k,j) = defor23(i,k,j) + tmp1(i,k,j)
913 END DO
914 END DO
915 END DO
916
917 ! End addition of the second term to defor23.
918
919 else ! defor13 and defor23 not needed
920
921 DO j = j_start, j_end
922 DO k = kts+1, ktf
923 DO i = i_start, i_end
924 defor13(i,k,j)= 0.
925 defor23(i,k,j)= 0.
926 enddo
927 enddo
928 enddo
929 endif
930 !-----------------------------------------------------------------------
931
932 !-----------------------------------------------------------------------
933 ! Update the boundary for defor13 and defor23 (might need to change
934 ! later).
935
936 IF ( .NOT. config_flags%periodic_x .AND. i_start .EQ. ids+1) THEN
937 DO j = jts, jte
938 DO k = kts, kte
939 defor13(ids,k,j) = defor13(ids+1,k,j)
940 defor23(ids,k,j) = defor23(ids+1,k,j)
941 END DO
942 END DO
943 END IF
944
945 IF ( .NOT. config_flags%periodic_y .AND. j_start .EQ. jds+1) THEN
946 DO k = kts, kte
947 DO i = its, ite
948 defor13(i,k,jds) = defor13(i,k,jds+1)
949 defor23(i,k,jds) = defor23(i,k,jds+1)
950 END DO
951 END DO
952 END IF
953
954 IF ( .NOT. config_flags%periodic_x .AND. i_end .EQ. ide-1) THEN
955 DO j = jts, jte
956 DO k = kts, kte
957 defor13(ide,k,j) = defor13(ide-1,k,j)
958 defor23(ide,k,j) = defor23(ide-1,k,j)
959 END DO
960 END DO
961 END IF
962
963 IF ( .NOT. config_flags%periodic_y .AND. j_end .EQ. jde-1) THEN
964 DO k = kts, kte
965 DO i = its, ite
966 defor13(i,k,jde) = defor13(i,k,jde-1)
967 defor23(i,k,jde) = defor23(i,k,jde-1)
968 END DO
969 END DO
970 END IF
971
972 ! End update of boundary for defor13 and defor23.
973 !-----------------------------------------------------------------------
974
975 ! The second three (defor12, defor13, defor23) of six deformation terms
976 ! are now calculated at vorticity points.
977 !=======================================================================
978
979 END SUBROUTINE cal_deform_and_div
980
981 !=======================================================================
982 !=======================================================================
983
984 SUBROUTINE calculate_km_kh( config_flags, dt, &
985 dampcoef, zdamp, damp_opt, &
986 xkmh, xkmv, xkhh, xkhv, &
987 BN2, khdif, kvdif, div, &
988 defor11, defor22, defor33, &
989 defor12, defor13, defor23, &
990 tke, p8w, t8w, theta, t, p, moist, &
991 dn, dnw, dx, dy, rdz, rdzw, cr_len, &
992 n_moist, cf1, cf2, cf3, warm_rain, &
993 kh_tke_upper_bound, kv_tke_upper_bound, &
994 msftx, msfty, &
995 ids, ide, jds, jde, kds, kde, &
996 ims, ime, jms, jme, kms, kme, &
997 its, ite, jts, jte, kts, kte )
998
999 ! History: Sep 2003 Changes by George Bryan and Jason Knievel, NCAR
1000 ! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
1001 ! ... ...
1002
1003 ! Purpose: This routine calculates exchange coefficients for the TKE
1004 ! scheme.
1005
1006 ! References: Klemp and Wilhelmson (JAS 1978)
1007 ! Deardorff (B-L Meteor 1980)
1008 ! Chen and Dudhia (NCAR WRF physics report 2000)
1009
1010 !-----------------------------------------------------------------------
1011 ! Begin declarations.
1012
1013 IMPLICIT NONE
1014
1015 TYPE( grid_config_rec_type ), INTENT( IN ) &
1016 :: config_flags
1017
1018 INTEGER, INTENT( IN ) &
1019 :: n_moist, damp_opt, &
1020 ids, ide, jds, jde, kds, kde, &
1021 ims, ime, jms, jme, kms, kme, &
1022 its, ite, jts, jte, kts, kte
1023
1024 LOGICAL, INTENT( IN ) &
1025 :: warm_rain
1026
1027 REAL, INTENT( IN ) &
1028 :: cr_len, dx, dy, zdamp, dt, dampcoef, cf1, cf2, cf3, khdif, kvdif
1029
1030 REAL, DIMENSION( kms:kme ), INTENT( IN ) &
1031 :: dnw, dn
1032
1033 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, n_moist ), INTENT( INOUT ) &
1034 :: moist
1035
1036 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
1037 :: xkmv, xkmh, xkhv, xkhh, BN2
1038
1039 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), INTENT( IN ) &
1040 :: defor11, defor22, defor33, defor12, defor13, defor23, &
1041 div, rdz, rdzw, p8w, t8w, theta, t, p
1042
1043 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
1044 :: tke
1045
1046 REAL, INTENT( IN ) &
1047 :: kh_tke_upper_bound, kv_tke_upper_bound
1048
1049 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
1050 :: msftx, msfty
1051
1052 ! Local variables.
1053
1054 INTEGER &
1055 :: i_start, i_end, j_start, j_end, ktf, i, j, k
1056
1057 ! End declarations.
1058 !-----------------------------------------------------------------------
1059
1060 ktf = MIN( kte, kde-1 )
1061 i_start = its
1062 i_end = MIN( ite, ide-1 )
1063 j_start = jts
1064 j_end = MIN( jte, jde-1 )
1065
1066 ! Select a scheme for calculating diffusion coefficients.
1067
1068 km_coef: SELECT CASE( config_flags%km_opt )
1069
1070 CASE (1)
1071 CALL isotropic_km( config_flags, xkmh, xkmv, &
1072 xkhh, xkhv, khdif, kvdif, &
1073 ids, ide, jds, jde, kds, kde, &
1074 ims, ime, jms, jme, kms, kme, &
1075 its, ite, jts, jte, kts, kte )
1076 CASE (2)
1077 CALL calculate_N2( config_flags, BN2, moist, &
1078 theta, t, p, p8w, t8w, &
1079 dnw, dn, rdz, rdzw, &
1080 n_moist, cf1, cf2, cf3, warm_rain, &
1081 ids, ide, jds, jde, kds, kde, &
1082 ims, ime, jms, jme, kms, kme, &
1083 its, ite, jts, jte, kts, kte )
1084
1085 CALL tke_km( config_flags, xkmh, xkmv, &
1086 xkhh, xkhv, BN2, tke, p8w, t8w, theta, &
1087 rdz, rdzw, dx, dy, cr_len, &
1088 kh_tke_upper_bound, kv_tke_upper_bound, &
1089 ids, ide, jds, jde, kds, kde, &
1090 ims, ime, jms, jme, kms, kme, &
1091 its, ite, jts, jte, kts, kte )
1092 CASE (3)
1093 CALL calculate_N2( config_flags, BN2, moist, &
1094 theta, t, p, p8w, t8w, &
1095 dnw, dn, rdz, rdzw, &
1096 n_moist, cf1, cf2, cf3, warm_rain, &
1097 ids, ide, jds, jde, kds, kde, &
1098 ims, ime, jms, jme, kms, kme, &
1099 its, ite, jts, jte, kts, kte )
1100
1101 CALL smag_km( config_flags, xkmh, xkmv, &
1102 xkhh, xkhv, BN2, div, &
1103 defor11, defor22, defor33, &
1104 defor12, defor13, defor23, &
1105 rdzw, dx, dy, cr_len, &
1106 ids, ide, jds, jde, kds, kde, &
1107 ims, ime, jms, jme, kms, kme, &
1108 its, ite, jts, jte, kts, kte )
1109 CASE (4)
1110 CALL smag2d_km( config_flags, xkmh, xkmv, &
1111 xkhh, xkhv, defor11, defor22, defor12, &
1112 rdzw, dx, dy, &
1113 ids, ide, jds, jde, kds, kde, &
1114 ims, ime, jms, jme, kms, kme, &
1115 its, ite, jts, jte, kts, kte )
1116 CASE DEFAULT
1117 CALL wrf_error_fatal( 'Please choose diffusion coefficient scheme' )
1118
1119 END SELECT km_coef
1120
1121 IF ( damp_opt .eq. 1 ) THEN
1122 CALL cal_dampkm( config_flags, xkmh, xkhh, xkmv, xkhv, &
1123 dx, dy, dt, dampcoef, rdz, rdzw, zdamp, &
1124 msftx, msfty, &
1125 ids, ide, jds, jde, kds, kde, &
1126 ims, ime, jms, jme, kms, kme, &
1127 its, ite, jts, jte, kts, kte )
1128 END IF
1129
1130 END SUBROUTINE calculate_km_kh
1131
1132 !=======================================================================
1133
1134 SUBROUTINE cal_dampkm( config_flags,xkmh,xkhh,xkmv,xkhv, &
1135 dx,dy,dt,dampcoef, &
1136 rdz, rdzw ,zdamp, &
1137 msftx, msfty, &
1138 ids,ide, jds,jde, kds,kde, &
1139 ims,ime, jms,jme, kms,kme, &
1140 its,ite, jts,jte, kts,kte )
1141
1142 !-----------------------------------------------------------------------
1143 ! Begin declarations.
1144
1145 IMPLICIT NONE
1146
1147 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
1148
1149 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1150 ims, ime, jms, jme, kms, kme, &
1151 its, ite, jts, jte, kts, kte
1152
1153 REAL , INTENT(IN ) :: zdamp,dx,dy,dt,dampcoef
1154
1155
1156 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: xkmh , &
1157 xkhh , &
1158 xkmv , &
1159 xkhv
1160
1161 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: rdz, &
1162 rdzw
1163
1164 REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: msftx, &
1165 msfty
1166 ! LOCAL VARS
1167
1168 INTEGER :: i_start, i_end, j_start, j_end, ktf, ktfm1, i, j, k
1169 REAL :: kmmax,kmmvmax,degrad90,dz,tmp
1170 REAL :: ds
1171 REAL , DIMENSION( its:ite ) :: deltaz
1172 REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: dampk,dampkv
1173
1174 ! End declarations.
1175 !-----------------------------------------------------------------------
1176
1177 ktf = min(kte,kde-1)
1178 ktfm1 = ktf-1
1179
1180 i_start = its
1181 i_end = MIN(ite,ide-1)
1182 j_start = jts
1183 j_end = MIN(jte,jde-1)
1184
1185 ! keep upper damping diffusion away from relaxation zones at boundaries if used
1186 IF(config_flags%specified .OR. config_flags%nested)THEN
1187 i_start = MAX(i_start,ids+config_flags%spec_bdy_width-1)
1188 i_end = MIN(i_end,ide-config_flags%spec_bdy_width)
1189 j_start = MAX(j_start,jds+config_flags%spec_bdy_width-1)
1190 j_end = MIN(j_end,jde-config_flags%spec_bdy_width)
1191 ENDIF
1192
1193 kmmax=dx*dx/dt
1194 degrad90=DEGRAD*90.
1195 DO j = j_start, j_end
1196
1197 k=ktf
1198 DO i = i_start, i_end
1199 ! Unmodified dx used above may produce very large diffusivities
1200 ! when msftx is very large. And the above formula ignores the fact
1201 ! that dy may now be different from dx as well. Let's fix that by
1202 ! defining a "ds" as the minimum of the "real-space" (physical
1203 ! distance) values of dx and dy, and then using that smallest value
1204 ! to calculate a point-by-point kmmax
1205 ds = MIN(dx/msftx(i,j),dy/msfty(i,j))
1206 kmmax=ds*ds/dt
1207
1208 ! deltaz(i)=0.5*dnw(k)/zeta_z(i,j)
1209 ! dz=dnw(k)/zeta_z(i,j)
1210 dz = 1./rdzw(i,k,j)
1211 deltaz(i) = 0.5*dz
1212
1213 kmmvmax=dz*dz/dt
1214 tmp=min(deltaz(i)/zdamp,1.)
1215 dampk(i,k,j)=cos(degrad90*tmp)*kmmax*dampcoef
1216 dampkv(i,k,j)=cos(degrad90*tmp)*kmmvmax*dampcoef
1217 ! set upper limit on vertical K (based on horizontal K)
1218 dampkv(i,k,j)=min(dampkv(i,k,j),dampk(i,k,j))
1219
1220 ENDDO
1221
1222 DO k = ktfm1,kts,-1
1223 DO i = i_start, i_end
1224 ! Unmodified dx used above may produce very large diffusivities
1225 ! when msftx is very large. And the above formula ignores the fact
1226 ! that dy may now be different from dx as well. Let's fix that by
1227 ! defining a "ds" as the minimum of the "real-space" (physical
1228 ! distance) values of dx and dy, and then using that smallest value
1229 ! to calculate a point-by-point kmmax
1230 ds = MIN(dx/msftx(i,j),dy/msfty(i,j))
1231 kmmax=ds*ds/dt
1232
1233 ! deltaz(i)=deltaz(i)+dn(k)/zeta_z(i,j)
1234 ! dz=dnw(k)/zeta_z(i,j)
1235 dz = 1./rdz(i,k,j)
1236 deltaz(i) = deltaz(i) + dz
1237 dz = 1./rdzw(i,k,j)
1238
1239 kmmvmax=dz*dz/dt
1240 tmp=min(deltaz(i)/zdamp,1.)
1241 dampk(i,k,j)=cos(degrad90*tmp)*kmmax*dampcoef
1242 dampkv(i,k,j)=cos(degrad90*tmp)*kmmvmax*dampcoef
1243 ! set upper limit on vertical K (based on horizontal K)
1244 dampkv(i,k,j)=min(dampkv(i,k,j),dampk(i,k,j))
1245 ENDDO
1246 ENDDO
1247
1248 ENDDO
1249
1250 DO j = j_start, j_end
1251 DO k = kts,ktf
1252 DO i = i_start, i_end
1253 xkmh(i,k,j)=max(xkmh(i,k,j),dampk(i,k,j))
1254 xkhh(i,k,j)=max(xkhh(i,k,j),dampk(i,k,j))
1255 xkmv(i,k,j)=max(xkmv(i,k,j),dampkv(i,k,j))
1256 xkhv(i,k,j)=max(xkhv(i,k,j),dampkv(i,k,j))
1257 ENDDO
1258 ENDDO
1259 ENDDO
1260
1261 END SUBROUTINE cal_dampkm
1262
1263 !=======================================================================
1264 !=======================================================================
1265
1266 SUBROUTINE calculate_N2( config_flags, BN2, moist, &
1267 theta, t, p, p8w, t8w, &
1268 dnw, dn, rdz, rdzw, &
1269 n_moist, cf1, cf2, cf3, warm_rain, &
1270 ids, ide, jds, jde, kds, kde, &
1271 ims, ime, jms, jme, kms, kme, &
1272 its, ite, jts, jte, kts, kte )
1273
1274 !-----------------------------------------------------------------------
1275 ! Begin declarations.
1276
1277 IMPLICIT NONE
1278
1279 TYPE( grid_config_rec_type ), INTENT( IN ) &
1280 :: config_flags
1281
1282 INTEGER, INTENT( IN ) &
1283 :: n_moist, &
1284 ids, ide, jds, jde, kds, kde, &
1285 ims, ime, jms, jme, kms, kme, &
1286 its, ite, jts, jte, kts, kte
1287
1288 LOGICAL, INTENT( IN ) &
1289 :: warm_rain
1290
1291 REAL, INTENT( IN ) &
1292 :: cf1, cf2, cf3
1293
1294 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
1295 :: BN2
1296
1297 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
1298 :: rdz, rdzw, theta, t, p, p8w, t8w
1299
1300 REAL, DIMENSION( kms:kme ), INTENT( IN ) &
1301 :: dnw, dn
1302
1303 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), INTENT( INOUT ) &
1304 :: moist
1305
1306 ! Local variables.
1307
1308 INTEGER &
1309 :: i, j, k, ktf, ispe, ktes1, ktes2, &
1310 i_start, i_end, j_start, j_end
1311
1312 REAL &
1313 :: coefa, thetaep1, thetaem1, qc_cr, es, tc, qlpqi, qsw, qsi, &
1314 tmpdz, xlvqv, thetaesfc, thetasfc, qvtop, qvsfc, thetatop, thetaetop
1315
1316 REAL, DIMENSION( its:ite, jts:jte ) &
1317 :: tmp1sfc, tmp1top
1318
1319 REAL, DIMENSION( its:ite, kts:kte, jts:jte ) &
1320 :: tmp1, qvs, qctmp
1321
1322 ! End declarations.
1323 !-----------------------------------------------------------------------
1324
1325 qc_cr = 0.00001 ! in Kg/Kg
1326
1327 ktf = MIN( kte, kde-1 )
1328 ktes1 = kte-1
1329 ktes2 = kte-2
1330
1331 i_start = its
1332 i_end = MIN( ite, ide-1 )
1333 j_start = jts
1334 j_end = MIN( jte, jde-1 )
1335
1336 IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
1337 config_flags%nested) i_start = MAX( ids+1, its )
1338 IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
1339 config_flags%nested) i_end = MIN( ide-2, ite )
1340 IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
1341 config_flags%nested) j_start = MAX( jds+1, jts )
1342 IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
1343 config_flags%nested) j_end = MIN( jde-2 ,jte )
1344 IF ( config_flags%periodic_x ) i_start = its
1345 IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
1346
1347 IF ( P_QC .GT. PARAM_FIRST_SCALAR) THEN
1348 DO j = j_start, j_end
1349 DO k = kts, ktf
1350 DO i = i_start, i_end
1351 qctmp(i,k,j) = moist(i,k,j,P_QC)
1352 END DO
1353 END DO
1354 END DO
1355 ELSE
1356 DO j = j_start, j_end
1357 DO k = kts, ktf
1358 DO i = i_start, i_end
1359 qctmp(i,k,j) = 0.0
1360 END DO
1361 END DO
1362 END DO
1363 END IF
1364
1365 DO j = jts, jte
1366 DO k = kts, kte
1367 DO i = its, ite
1368 tmp1(i,k,j) = 0.0
1369 END DO
1370 END DO
1371 END DO
1372
1373 DO j = jts,jte
1374 DO i = its,ite
1375 tmp1sfc(i,j) = 0.0
1376 tmp1top(i,j) = 0.0
1377 END DO
1378 END DO
1379
1380 DO ispe = PARAM_FIRST_SCALAR, n_moist
1381 IF ( ispe .EQ. P_QV .OR. ispe .EQ. P_QC .OR. ispe .EQ. P_QI) THEN
1382 DO j = j_start, j_end
1383 DO k = kts, ktf
1384 DO i = i_start, i_end
1385 tmp1(i,k,j) = tmp1(i,k,j) + moist(i,k,j,ispe)
1386 END DO
1387 END DO
1388 END DO
1389
1390 DO j = j_start, j_end
1391 DO i = i_start, i_end
1392 tmp1sfc(i,j) = tmp1sfc(i,j) + &
1393 cf1 * moist(i,1,j,ispe) + &
1394 cf2 * moist(i,2,j,ispe) + &
1395 cf3 * moist(i,3,j,ispe)
1396 tmp1top(i,j) = tmp1top(i,j) + &
1397 moist(i,ktes1,j,ispe) + &
1398 ( moist(i,ktes1,j,ispe) - moist(i,ktes2,j,ispe) ) * &
1399 0.5 * dnw(ktes1) / dn(ktes1)
1400 END DO
1401 END DO
1402 END IF
1403 END DO
1404
1405 ! Calculate saturation mixing ratio.
1406
1407 DO j = j_start, j_end
1408 DO k = kts, ktf
1409 DO i = i_start, i_end
1410 tc = t(i,k,j) - SVPT0
1411 es = 1000.0 * SVP1 * EXP( SVP2 * tc / ( t(i,k,j) - SVP3 ) )
1412 qvs(i,k,j) = EP_2 * es / ( p(i,k,j) - es )
1413 END DO
1414 END DO
1415 END DO
1416
1417 DO j = j_start, j_end
1418 DO k = kts+1, ktf-1
1419 DO i = i_start, i_end
1420 tmpdz = 1.0 / rdz(i,k,j) + 1.0 / rdz(i,k+1,j)
1421 IF ( moist(i,k,j,P_QV) .GE. qvs(i,k,j) .OR. qctmp(i,k,j) .GE. qc_cr) THEN
1422 xlvqv = XLV * moist(i,k,j,P_QV)
1423 coefa = ( 1.0 + xlvqv / R_d / t(i,k,j) ) / &
1424 ( 1.0 + XLV * xlvqv / Cp / R_v / t(i,k,j) / t(i,k,j) ) / &
1425 theta(i,k,j)
1426 thetaep1 = theta(i,k+1,j) * &
1427 ( 1.0 + XLV * qvs(i,k+1,j) / Cp / t(i,k+1,j) )
1428 thetaem1 = theta(i,k-1,j) * &
1429 ( 1.0 + XLV * qvs(i,k-1,j) / Cp / t(i,k-1,j) )
1430 BN2(i,k,j) = g * ( coefa * ( thetaep1 - thetaem1 ) / tmpdz - &
1431 ( tmp1(i,k+1,j) - tmp1(i,k-1,j) ) / tmpdz )
1432 ELSE
1433 BN2(i,k,j) = g * ( (theta(i,k+1,j) - theta(i,k-1,j) ) / &
1434 theta(i,k,j) / tmpdz + &
1435 1.61 * ( moist(i,k+1,j,P_QV) - moist(i,k-1,j,P_QV) ) / &
1436 tmpdz - &
1437 ( tmp1(i,k+1,j) - tmp1(i,k-1,j) ) / tmpdz )
1438 ENDIF
1439 END DO
1440 END DO
1441 END DO
1442
1443 k = kts
1444 DO j = j_start, j_end
1445 DO i = i_start, i_end
1446 tmpdz = 1.0 / rdz(i,k+1,j) + 0.5 / rdzw(i,k,j)
1447 thetasfc = T8w(i,kts,j) / ( p8w(i,k,j) / p1000mb )**( R_d / Cp )
1448 IF ( moist(i,k,j,P_QV) .GE. qvs(i,k,j) .OR. qctmp(i,k,j) .GE. qc_cr) THEN
1449 qvsfc = cf1 * qvs(i,1,j) + &
1450 cf2 * qvs(i,2,j) + &
1451 cf3 * qvs(i,3,j)
1452 xlvqv = XLV * moist(i,k,j,P_QV)
1453 coefa = ( 1.0 + xlvqv / R_d / t(i,k,j) ) / &
1454 ( 1.0 + XLV * xlvqv / Cp / R_v / t(i,k,j) / t(i,k,j) ) / &
1455 theta(i,k,j)
1456 thetaep1 = theta(i,k+1,j) * &
1457 ( 1.0 + XLV * qvs(i,k+1,j) / Cp / t(i,k+1,j) )
1458 thetaesfc = thetasfc * &
1459 ( 1.0 + XLV * qvsfc / Cp / t8w(i,kts,j) )
1460 BN2(i,k,j) = g * ( coefa * ( thetaep1 - thetaesfc ) / tmpdz - &
1461 ( tmp1(i,k+1,j) - tmp1sfc(i,j) ) / tmpdz )
1462 ELSE
1463 qvsfc = cf1 * moist(i,1,j,P_QV) + &
1464 cf2 * moist(i,2,j,P_QV) + &
1465 cf3 * moist(i,3,j,P_QV)
1466 ! BN2(i,k,j) = g * ( ( theta(i,k+1,j) - thetasfc ) / &
1467 ! theta(i,k,j) / tmpdz + &
1468 ! 1.61 * ( moist(i,k+1,j,P_QV) - qvsfc ) / &
1469 ! tmpdz - &
1470 ! ( tmp1(i,k+1,j) - tmp1sfc(i,j) ) / tmpdz )
1471 !...... MARTA: change in computation of BN2 at the surface, WCS 040331
1472
1473 tmpdz= 1./rdzw(i,k,j) ! controlare come calcola rdzw
1474 BN2(i,k,j) = g * ( ( theta(i,k+1,j) - theta(i,k,j)) / &
1475 theta(i,k,j) / tmpdz + &
1476 1.61 * ( moist(i,k+1,j,P_QV) - qvsfc ) / &
1477 tmpdz - &
1478 ( tmp1(i,k+1,j) - tmp1sfc(i,j) ) / tmpdz )
1479 ! end of MARTA/WCS change
1480
1481 ENDIF
1482 END DO
1483 END DO
1484
1485
1486 !...... MARTA: change in computation of BN2 at the top, WCS 040331
1487 DO j = j_start, j_end
1488 DO i = i_start, i_end
1489 BN2(i,ktf,j)=BN2(i,ktf-1,j)
1490 END DO
1491 END DO
1492 ! end of MARTA/WCS change
1493
1494 END SUBROUTINE calculate_N2
1495
1496 !=======================================================================
1497 !=======================================================================
1498
1499 SUBROUTINE isotropic_km( config_flags, &
1500 xkmh,xkmv,xkhh,xkhv,khdif,kvdif, &
1501 ids,ide, jds,jde, kds,kde, &
1502 ims,ime, jms,jme, kms,kme, &
1503 its,ite, jts,jte, kts,kte )
1504
1505 !-----------------------------------------------------------------------
1506 ! Begin declarations.
1507
1508 IMPLICIT NONE
1509
1510 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
1511
1512 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1513 ims, ime, jms, jme, kms, kme, &
1514 its, ite, jts, jte, kts, kte
1515
1516 REAL , INTENT(IN ) :: khdif,kvdif
1517
1518 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: xkmh, &
1519 xkmv, &
1520 xkhh, &
1521 xkhv
1522 ! LOCAL VARS
1523
1524 INTEGER :: i_start, i_end, j_start, j_end, ktf, i, j, k
1525 REAL :: khdif3,kvdif3
1526
1527 ! End declarations.
1528 !-----------------------------------------------------------------------
1529
1530 ktf = kte
1531
1532 i_start = its
1533 i_end = MIN(ite,ide-1)
1534 j_start = jts
1535 j_end = MIN(jte,jde-1)
1536
1537 ! khdif3=khdif*3.
1538 ! kvdif3=kvdif*3.
1539 khdif3=khdif/prandtl
1540 kvdif3=kvdif/prandtl
1541
1542 DO j = j_start, j_end
1543 DO k = kts, ktf
1544 DO i = i_start, i_end
1545 xkmh(i,k,j)=khdif
1546 xkmv(i,k,j)=kvdif
1547 xkhh(i,k,j)=khdif3
1548 xkhv(i,k,j)=kvdif3
1549 ENDDO
1550 ENDDO
1551 ENDDO
1552
1553 END SUBROUTINE isotropic_km
1554
1555 !=======================================================================
1556 !=======================================================================
1557
1558 SUBROUTINE smag_km( config_flags,xkmh,xkmv,xkhh,xkhv,BN2, &
1559 div,defor11,defor22,defor33,defor12, &
1560 defor13,defor23, &
1561 rdzw,dx,dy,cr_len_in, &
1562 ids,ide, jds,jde, kds,kde, &
1563 ims,ime, jms,jme, kms,kme, &
1564 its,ite, jts,jte, kts,kte )
1565
1566 !-----------------------------------------------------------------------
1567 ! Begin declarations.
1568
1569 IMPLICIT NONE
1570
1571 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
1572
1573 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1574 ims, ime, jms, jme, kms, kme, &
1575 its, ite, jts, jte, kts, kte
1576
1577 REAL , INTENT(IN ) :: cr_len_in, dx, dy
1578
1579
1580 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: BN2, &
1581 rdzw
1582
1583 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: xkmh, &
1584 xkmv, &
1585 xkhh, &
1586 xkhv
1587
1588 REAL , DIMENSION( ims:ime , kms:kme, jms:jme ), INTENT(IN ) :: &
1589 defor11, &
1590 defor22, &
1591 defor33, &
1592 defor12, &
1593 defor13, &
1594 defor23, &
1595 div
1596
1597 ! LOCAL VARS
1598
1599 INTEGER :: i_start, i_end, j_start, j_end, ktf, i, j, k
1600 REAL :: deltas, tmp, pr, mlen_h, mlen_v, cr_len
1601
1602 REAL, DIMENSION( its:ite , kts:kte , jts:jte ) :: def2
1603
1604 ! End declarations.
1605 !-----------------------------------------------------------------------
1606
1607 ktf = min(kte,kde-1)
1608
1609 i_start = its
1610 i_end = MIN(ite,ide-1)
1611 j_start = jts
1612 j_end = MIN(jte,jde-1)
1613
1614 IF ( config_flags%open_xs .or. config_flags%specified .or. &
1615 config_flags%nested) i_start = MAX(ids+1,its)
1616 IF ( config_flags%open_xe .or. config_flags%specified .or. &
1617 config_flags%nested) i_end = MIN(ide-2,ite)
1618 IF ( config_flags%open_ys .or. config_flags%specified .or. &
1619 config_flags%nested) j_start = MAX(jds+1,jts)
1620 IF ( config_flags%open_ye .or. config_flags%specified .or. &
1621 config_flags%nested) j_end = MIN(jde-2,jte)
1622 IF ( config_flags%periodic_x ) i_start = its
1623 IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
1624
1625 pr = prandtl
1626 cr_len = cr_len_in
1627
1628 do j=j_start,j_end
1629 do k=kts,ktf
1630 do i=i_start,i_end
1631 def2(i,k,j)=0.5*(defor11(i,k,j)*defor11(i,k,j) + &
1632 defor22(i,k,j)*defor22(i,k,j) + &
1633 defor33(i,k,j)*defor33(i,k,j))
1634 enddo
1635 enddo
1636 enddo
1637
1638 do j=j_start,j_end
1639 do k=kts,ktf
1640 do i=i_start,i_end
1641 tmp=0.25*(defor12(i ,k,j)+defor12(i ,k,j+1)+ &
1642 defor12(i+1,k,j)+defor12(i+1,k,j+1))
1643 def2(i,k,j)=def2(i,k,j)+tmp*tmp
1644 enddo
1645 enddo
1646 enddo
1647
1648 do j=j_start,j_end
1649 do k=kts,ktf
1650 do i=i_start,i_end
1651 tmp=0.25*(defor13(i ,k+1,j)+defor13(i ,k,j)+ &
1652 defor13(i+1,k+1,j)+defor13(i+1,k,j))
1653 def2(i,k,j)=def2(i,k,j)+tmp*tmp
1654 enddo
1655 enddo
1656 enddo
1657
1658 do j=j_start,j_end
1659 do k=kts,ktf
1660 do i=i_start,i_end
1661 tmp=0.25*(defor23(i,k+1,j )+defor23(i,k,j )+ &
1662 defor23(i,k+1,j+1)+defor23(i,k,j+1))
1663 def2(i,k,j)=def2(i,k,j)+tmp*tmp
1664 enddo
1665 enddo
1666 enddo
1667 !
1668 cr_len = dx + 1. ! hardwire for mixing length = (dx*dy*dz)**(1/3).
1669 ! remove this for the alternate formulation
1670
1671 IF (dx .gt. cr_len) THEN
1672 mlen_h=sqrt(dx*dy)
1673 DO j = j_start, j_end
1674 DO k = kts, ktf
1675 DO i = i_start, i_end
1676 mlen_v= 1./rdzw(i,k,j)
1677 tmp=max(0.,def2(i,k,j)-BN2(i,k,j)/pr)
1678 tmp=tmp**0.5
1679 xkmh(i,k,j)=max(c_s*c_s*mlen_h*mlen_h*tmp, 1.0E-6*mlen_h*mlen_h )
1680 xkmh(i,k,j)=min(xkmh(i,k,j), 10.*mlen_h )
1681 xkmv(i,k,j)=max(c_s*c_s*mlen_v*mlen_v*tmp, 1.0E-6*mlen_v*mlen_v )
1682 xkmv(i,k,j)=min(xkmv(i,k,j), 10.*mlen_v )
1683 xkhh(i,k,j)=xkmh(i,k,j)/pr
1684 xkhv(i,k,j)=xkmv(i,k,j)/pr
1685 ENDDO
1686 ENDDO
1687 ENDDO
1688 ELSE
1689 DO j = j_start, j_end
1690 DO k = kts, ktf
1691 DO i = i_start, i_end
1692 deltas=(dx*dy/rdzw(i,k,j))**0.33333333
1693 tmp=max(0.,def2(i,k,j)-BN2(i,k,j)/pr)
1694 tmp=tmp**0.5
1695 xkmh(i,k,j)=max(c_s*c_s*deltas*deltas*tmp, 1.0E-6*deltas*deltas )
1696 xkmh(i,k,j)=min(xkmh(i,k,j), 10.*deltas )
1697 xkmv(i,k,j)=xkmh(i,k,j)
1698 xkhh(i,k,j)=xkmh(i,k,j)/pr
1699 xkhv(i,k,j)=xkmv(i,k,j)/pr
1700 ENDDO
1701 ENDDO
1702 ENDDO
1703 ENDIF
1704
1705 END SUBROUTINE smag_km
1706
1707 !=======================================================================
1708 !=======================================================================
1709
1710 SUBROUTINE smag2d_km( config_flags,xkmh,xkmv,xkhh,xkhv, &
1711 defor11,defor22,defor12, &
1712 rdzw,dx,dy, &
1713 ids,ide, jds,jde, kds,kde, &
1714 ims,ime, jms,jme, kms,kme, &
1715 its,ite, jts,jte, kts,kte )
1716
1717 !-----------------------------------------------------------------------
1718 ! Begin declarations.
1719
1720 IMPLICIT NONE
1721
1722 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
1723
1724 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1725 ims, ime, jms, jme, kms, kme, &
1726 its, ite, jts, jte, kts, kte
1727
1728 REAL , INTENT(IN ) :: dx, dy
1729
1730
1731 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rdzw
1732
1733 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: xkmh, &
1734 xkmv, &
1735 xkhh, &
1736 xkhv
1737
1738 REAL , DIMENSION( ims:ime , kms:kme, jms:jme ), INTENT(IN ) :: &
1739 defor11, &
1740 defor22, &
1741 defor12
1742
1743 ! LOCAL VARS
1744
1745 INTEGER :: i_start, i_end, j_start, j_end, ktf, i, j, k
1746 REAL :: deltas, tmp, pr, mlen_h
1747
1748 REAL, DIMENSION( its:ite , kts:kte , jts:jte ) :: def2
1749
1750 ! End declarations.
1751 !-----------------------------------------------------------------------
1752
1753 ktf = min(kte,kde-1)
1754
1755 i_start = its
1756 i_end = MIN(ite,ide-1)
1757 j_start = jts
1758 j_end = MIN(jte,jde-1)
1759
1760 IF ( config_flags%open_xs .or. config_flags%specified .or. &
1761 config_flags%nested) i_start = MAX(ids+1,its)
1762 IF ( config_flags%open_xe .or. config_flags%specified .or. &
1763 config_flags%nested) i_end = MIN(ide-2,ite)
1764 IF ( config_flags%open_ys .or. config_flags%specified .or. &
1765 config_flags%nested) j_start = MAX(jds+1,jts)
1766 IF ( config_flags%open_ye .or. config_flags%specified .or. &
1767 config_flags%nested) j_end = MIN(jde-2,jte)
1768 IF ( config_flags%periodic_x ) i_start = its
1769 IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
1770
1771 pr=prandtl
1772
1773 do j=j_start,j_end
1774 do k=kts,ktf
1775 do i=i_start,i_end
1776 def2(i,k,j)=0.25*((defor11(i,k,j)-defor22(i,k,j))*(defor11(i,k,j)-defor22(i,k,j)))
1777 tmp=0.25*(defor12(i ,k,j)+defor12(i ,k,j+1)+ &
1778 defor12(i+1,k,j)+defor12(i+1,k,j+1))
1779 def2(i,k,j)=def2(i,k,j)+tmp*tmp
1780 enddo
1781 enddo
1782 enddo
1783 !
1784 mlen_h=sqrt(dx*dy)
1785 DO j = j_start, j_end
1786 DO k = kts, ktf
1787 DO i = i_start, i_end
1788 tmp=sqrt(def2(i,k,j))
1789 ! xkmh(i,k,j)=max(c_s*c_s*mlen_h*mlen_h*tmp, 1.0E-6*mlen_h*mlen_h )
1790 xkmh(i,k,j)=c_s*c_s*mlen_h*mlen_h*tmp
1791 xkmh(i,k,j)=min(xkmh(i,k,j), 10.*mlen_h )
1792 xkmv(i,k,j)=0.
1793 xkhh(i,k,j)=xkmh(i,k,j)/pr
1794 xkhv(i,k,j)=0.
1795 ENDDO
1796 ENDDO
1797 ENDDO
1798
1799 END SUBROUTINE smag2d_km
1800
1801 !=======================================================================
1802 !=======================================================================
1803
1804 SUBROUTINE tke_km( config_flags, xkmh, xkmv, xkhh, xkhv, &
1805 bn2, tke, p8w, t8w, theta, &
1806 rdz, rdzw, dx,dy, cr_len_in, &
1807 kh_tke_upper_bound, kv_tke_upper_bound, &
1808 ids, ide, jds, jde, kds, kde, &
1809 ims, ime, jms, jme, kms, kme, &
1810 its, ite, jts, jte, kts, kte )
1811
1812 ! History: Sep 2003 Changes by Jason Knievel and George Bryan, NCAR
1813 ! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
1814 ! ... ...
1815
1816 ! Purpose: This routine calculates the exchange coefficients for the
1817 ! TKE turbulence parameterization.
1818
1819 ! References: Klemp and Wilhelmson (JAS 1978)
1820 ! Chen and Dudhia (NCAR WRF physics report 2000)
1821
1822 !-----------------------------------------------------------------------
1823 ! Begin declarations.
1824
1825 IMPLICIT NONE
1826
1827 TYPE( grid_config_rec_type ), INTENT( IN ) &
1828 :: config_flags
1829
1830 INTEGER, INTENT( IN ) &
1831 :: ids, ide, jds, jde, kds, kde, &
1832 ims, ime, jms, jme, kms, kme, &
1833 its, ite, jts, jte, kts, kte
1834
1835 REAL, INTENT( IN ) &
1836 :: cr_len_in, dx, dy
1837
1838 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
1839 :: tke, p8w, t8w, theta, rdz, rdzw, bn2
1840
1841 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
1842 :: xkmh, xkmv, xkhh, xkhv
1843
1844 REAL, INTENT( IN ) &
1845 :: kh_tke_upper_bound, kv_tke_upper_bound
1846
1847 ! Local variables.
1848
1849 REAL, DIMENSION( its:ite, kts:kte, jts:jte ) &
1850 :: l_scale
1851
1852 REAL, DIMENSION( its:ite, kts:kte, jts:jte ) &
1853 :: dthrdn
1854
1855 REAL &
1856 :: deltas, tmp, mlen_s, mlen_h, mlen_v, tmpdz, &
1857 thetasfc, thetatop, minkx, pr_inv, pr_inv_h, pr_inv_v, cr_len
1858
1859 INTEGER &
1860 :: i_start, i_end, j_start, j_end, ktf, i, j, k
1861
1862 REAL, PARAMETER :: tke_seed_value = 1.e-06
1863 REAL :: tke_seed
1864 REAL, PARAMETER :: epsilon = 1.e-10
1865
1866 ! End declarations.
1867 !-----------------------------------------------------------------------
1868
1869 ktf = MIN( kte, kde-1 )
1870 i_start = its
1871 i_end = MIN( ite, ide-1 )
1872 j_start = jts
1873 j_end = MIN( jte, jde-1 )
1874
1875 IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
1876 config_flags%nested) i_start = MAX( ids+1, its )
1877 IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
1878 config_flags%nested) i_end = MIN( ide-2, ite )
1879 IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
1880 config_flags%nested) j_start = MAX( jds+1, jts )
1881 IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
1882 config_flags%nested) j_end = MIN( jde-2, jte)
1883 IF ( config_flags%periodic_x ) i_start = its
1884 IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
1885
1886 ! in the absence of surface drag or a surface heat flux, there
1887 ! is no way to generate tke without pre-existing tke. Use
1888 ! tke_seed if the drag and flux are off.
1889
1890 cr_len = cr_len_in
1891 tke_seed = tke_seed_value
1892 if( (config_flags%tke_drag_coefficient .gt. epsilon) .or. &
1893 (config_flags%tke_heat_flux .gt. epsilon) ) tke_seed = 0.
1894
1895 DO j = j_start, j_end
1896 DO k = kts+1, ktf-1
1897 DO i = i_start, i_end
1898 tmpdz = 1.0 / ( rdz(i,k+1,j) + rdz(i,k,j) )
1899 dthrdn(i,k,j) = ( theta(i,k+1,j) - theta(i,k-1,j) ) / tmpdz
1900 END DO
1901 END DO
1902 END DO
1903
1904 k = kts
1905 DO j = j_start, j_end
1906 DO i = i_start, i_end
1907 tmpdz = 1.0 / ( rdzw(i,k+1,j) + rdzw(i,k,j) )
1908 thetasfc = T8w(i,kts,j) / ( p8w(i,k,j) / p1000mb )**( R_d / Cp )
1909 dthrdn(i,k,j) = ( theta(i,k+1,j) - thetasfc ) / tmpdz
1910 END DO
1911 END DO
1912
1913 k = ktf
1914 DO j = j_start, j_end
1915 DO i = i_start, i_end
1916 tmpdz = 1.0 / rdz(i,k,j) + 0.5 / rdzw(i,k,j)
1917 thetatop = T8w(i,kde,j) / ( p8w(i,kde,j) / p1000mb )**( R_d / Cp )
1918 dthrdn(i,k,j) = ( thetatop - theta(i,k-1,j) ) / tmpdz
1919 END DO
1920 END DO
1921
1922 cr_len = dx + 1.0 ! hardwire for mixing length = (dx*dy*dz)**(1/3).
1923 ! remove this for the alternate formulation
1924
1925 IF ( dx .gt. cr_len ) THEN
1926 mlen_h = SQRT( dx * dy )
1927 DO j = j_start, j_end
1928 DO k = kts, ktf
1929 DO i = i_start, i_end
1930 tmp = SQRT( MAX( tke(i,k,j), tke_seed ) )
1931 deltas = 1.0 / rdzw(i,k,j)
1932 mlen_v = deltas
1933 IF ( dthrdn(i,k,j) .GT. 0.) THEN
1934 mlen_s = 0.76 * tmp / ( ABS( g / theta(i,k,j) * dthrdn(i,k,j) ) )**0.5
1935 mlen_v = MIN( mlen_v, mlen_s )
1936 END IF
1937 xkmh(i,k,j) = MAX( c_k * tmp * mlen_h, 1.0E-6 * mlen_h * mlen_h )
1938 xkmh(i,k,j) = MIN( xkmh(i,k,j), 10.0 * mlen_h )
1939 xkmv(i,k,j) = MAX( c_k * tmp * mlen_v, 1.0E-6 * deltas * deltas )
1940 xkmv(i,k,j) = MIN( xkmv(i,k,j), 10.0 * deltas )
1941 pr_inv_h = 1./prandtl
1942 pr_inv_v = 1.0 + 2.0 * mlen_v / deltas
1943 xkhh(i,k,j) = xkmh(i,k,j) * pr_inv_h
1944 xkhv(i,k,j) = xkmv(i,k,j) * pr_inv_v
1945 END DO
1946 END DO
1947 END DO
1948 ELSE
1949 CALL calc_l_scale( config_flags, tke, BN2, l_scale, &
1950 i_start, i_end, ktf, j_start, j_end, &
1951 dx, dy, rdzw, &
1952 ids, ide, jds, jde, kds, kde, &
1953 ims, ime, jms, jme, kms, kme, &
1954 its, ite, jts, jte, kts, kte )
1955 DO j = j_start, j_end
1956 DO k = kts, ktf
1957 DO i = i_start, i_end
1958 tmp = SQRT( MAX( tke(i,k,j), tke_seed ) )
1959 deltas = ( dx * dy / rdzw(i,k,j) )**0.33333333
1960 xkmh(i,k,j) = c_k * tmp * l_scale(i,k,j)
1961 xkmh(i,k,j) = MIN( kh_tke_upper_bound, xkmh(i,k,j) )
1962 xkmv(i,k,j) = c_k * tmp * l_scale(i,k,j)
1963 xkmv(i,k,j) = MIN( kv_tke_upper_bound, xkmv(i,k,j) )
1964 pr_inv = 1.0 + 2.0 * l_scale(i,k,j) / deltas
1965 xkhh(i,k,j) = MIN( kh_tke_upper_bound*pr_inv, xkmh(i,k,j) * pr_inv )
1966 xkhv(i,k,j) = MIN( kv_tke_upper_bound*pr_inv, xkmv(i,k,j) * pr_inv )
1967 END DO
1968 END DO
1969 END DO
1970 END IF
1971
1972 END SUBROUTINE tke_km
1973
1974 !=======================================================================
1975 !=======================================================================
1976
1977 SUBROUTINE calc_l_scale( config_flags, tke, BN2, l_scale, &
1978 i_start, i_end, ktf, j_start, j_end, &
1979 dx, dy, rdzw, &
1980 ids, ide, jds, jde, kds, kde, &
1981 ims, ime, jms, jme, kms, kme, &
1982 its, ite, jts, jte, kts, kte )
1983
1984 ! History: Sep 2003 Written by Bryan and Knievel, NCAR
1985
1986 ! Purpose: This routine calculates the length scale, based on stability,
1987 ! for TKE parameterization of subgrid-scale turbulence.
1988
1989 !-----------------------------------------------------------------------
1990 ! Begin declarations.
1991
1992 IMPLICIT NONE
1993
1994 TYPE( grid_config_rec_type ), INTENT( IN ) &
1995 :: config_flags
1996
1997 INTEGER, INTENT( IN ) &
1998 :: i_start, i_end, ktf, j_start, j_end, &
1999 ids, ide, jds, jde, kds, kde, &
2000 ims, ime, jms, jme, kms, kme, &
2001 its, ite, jts, jte, kts, kte
2002
2003 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
2004 :: BN2, tke, rdzw
2005
2006 REAL, INTENT( IN ) &
2007 :: dx, dy
2008
2009 REAL, DIMENSION( its:ite, kts:kte, jts:jte ), INTENT( OUT ) &
2010 :: l_scale
2011
2012 ! Local variables.
2013
2014 INTEGER &
2015 :: i, j, k
2016
2017 REAL &
2018 :: deltas, tmp
2019
2020 ! End declarations.
2021 !-----------------------------------------------------------------------
2022
2023 DO j = j_start, j_end
2024 DO k = kts, ktf
2025 DO i = i_start, i_end
2026 deltas = ( dx * dy / rdzw(i,k,j) )**0.33333333
2027 l_scale(i,k,j) = deltas
2028
2029 IF ( BN2(i,k,j) .gt. 1.0e-6 ) THEN
2030 tmp = SQRT( MAX( tke(i,k,j), 1.0e-6 ) )
2031 l_scale(i,k,j) = 0.76 * tmp / SQRT( BN2(i,k,j) )
2032 l_scale(i,k,j) = MIN( l_scale(i,k,j), deltas)
2033 l_scale(i,k,j) = MAX( l_scale(i,k,j), 0.001 * deltas )
2034 END IF
2035
2036 END DO
2037 END DO
2038 END DO
2039
2040 END SUBROUTINE calc_l_scale
2041
2042 !=======================================================================
2043 !=======================================================================
2044
2045 SUBROUTINE horizontal_diffusion_2 ( rt_tendf, ru_tendf, rv_tendf, rw_tendf, &
2046 tke_tendf, &
2047 moist_tendf, n_moist, &
2048 chem_tendf, n_chem, &
2049 scalar_tendf, n_scalar, &
2050 thp, theta, mu, tke, config_flags, &
2051 defor11, defor22, defor12, &
2052 defor13, defor23, div, &
2053 moist, chem, scalar, &
2054 msfux, msfuy, msfvx, msfvy, &
2055 msftx, msfty, xkmh, xkhh,km_opt, &
2056 rdx, rdy, rdz, rdzw, fnm, fnp, &
2057 cf1, cf2, cf3, &
2058 ht, dx, dy, dt, zx, zy, dn, dnw, &
2059 ids, ide, jds, jde, kds, kde, &
2060 ims, ime, jms, jme, kms, kme, &
2061 its, ite, jts, jte, kts, kte )
2062
2063 !-----------------------------------------------------------------------
2064 ! Begin declarations.
2065
2066 IMPLICIT NONE
2067
2068 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2069
2070 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2071 ims, ime, jms, jme, kms, kme, &
2072 its, ite, jts, jte, kts, kte
2073
2074 INTEGER , INTENT(IN ) :: n_moist, n_chem, n_scalar, km_opt
2075
2076 REAL , INTENT(IN ) :: cf1, cf2, cf3, dx, dy, dt
2077
2078 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
2079 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
2080 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
2081 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
2082
2083 REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfux, &
2084 msfuy, &
2085 msfvx, &
2086 msfvy, &
2087 msftx, &
2088 msfty, &
2089 mu, ht
2090
2091 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::rt_tendf,&
2092 ru_tendf,&
2093 rv_tendf,&
2094 rw_tendf,&
2095 tke_tendf
2096
2097 REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), &
2098 INTENT(INOUT) :: moist_tendf
2099
2100 REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_chem), &
2101 INTENT(INOUT) :: chem_tendf
2102
2103 REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_scalar), &
2104 INTENT(INOUT) :: scalar_tendf
2105
2106 REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), &
2107 INTENT(IN ) :: moist
2108
2109 REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_chem), &
2110 INTENT(IN ) :: chem
2111
2112 REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_scalar) , &
2113 INTENT(IN ) :: scalar
2114
2115 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor11, &
2116 defor22, &
2117 defor12, &
2118 defor13, &
2119 defor23, &
2120 div, &
2121 xkmh, &
2122 xkhh, &
2123 zx, &
2124 zy, &
2125 theta, &
2126 thp, &
2127 tke, &
2128 rdz, &
2129 rdzw
2130
2131
2132 REAL , INTENT(IN ) :: rdx, &
2133 rdy
2134
2135 ! LOCAL VARS
2136
2137 INTEGER :: im, ic, is
2138
2139 ! REAL , DIMENSION(its-1:ite+1, kts:kte, jts-1:jte+1) :: xkhh
2140
2141 ! End declarations.
2142 !-----------------------------------------------------------------------
2143
2144 !-----------------------------------------------------------------------
2145 ! Call diffusion subroutines.
2146
2147 CALL horizontal_diffusion_u_2( ru_tendf, mu, config_flags, &
2148 defor11, defor12, div, &
2149 tke(ims,kms,jms), &
2150 msfux, msfuy, xkmh, rdx, rdy, fnm, fnp, &
2151 zx, zy, rdzw, dx, dy, dt, &
2152 ids, ide, jds, jde, kds, kde, &
2153 ims, ime, jms, jme, kms, kme, &
2154 its, ite, jts, jte, kts, kte )
2155
2156 CALL horizontal_diffusion_v_2( rv_tendf, mu, config_flags, &
2157 defor12, defor22, div, &
2158 tke(ims,kms,jms), &
2159 msfvx, msfvy, xkmh, rdx, rdy, fnm, fnp, &
2160 zx, zy, rdzw, dx, dy, dt, &
2161 ids, ide, jds, jde, kds, kde, &
2162 ims, ime, jms, jme, kms, kme, &
2163 its, ite, jts, jte, kts, kte )
2164
2165 CALL horizontal_diffusion_w_2( rw_tendf, mu, config_flags, &
2166 defor13, defor23, div, &
2167 tke(ims,kms,jms), &
2168 msftx, msfty, xkmh, rdx, rdy, fnm, fnp, &
2169 zx, zy, rdz, dx, dy, dt, &
2170 ids, ide, jds, jde, kds, kde, &
2171 ims, ime, jms, jme, kms, kme, &
2172 its, ite, jts, jte, kts, kte )
2173
2174 CALL horizontal_diffusion_s ( rt_tendf, mu, config_flags, thp, &
2175 msftx, msfty, msfux, msfuy, &
2176 msfvx, msfvy, xkhh, rdx, rdy, &
2177 fnm, fnp, cf1, cf2, cf3, &
2178 ht, dx, dy, dt, zx, zy, rdz, rdzw, dnw, dn, &
2179 .false., 'T', &
2180 ids, ide, jds, jde, kds, kde, &
2181 ims, ime, jms, jme, kms, kme, &
2182 its, ite, jts, jte, kts, kte )
2183
2184 IF (km_opt .eq. 2) &
2185 CALL horizontal_diffusion_s ( tke_tendf(ims,kms,jms), &
2186 mu, config_flags, &
2187 tke(ims,kms,jms), &
2188 msftx, msfty, msfux, msfuy, &
2189 msfvx, msfvy, xkhh, rdx, rdy, &
2190 fnm, fnp, cf1, cf2, cf3, &
2191 ht, dx, dy, dt, zx, zy, rdz, rdzw, dnw, dn, &
2192 .true., 'E', &
2193 ids, ide, jds, jde, kds, kde, &
2194 ims, ime, jms, jme, kms, kme, &
2195 its, ite, jts, jte, kts, kte )
2196
2197 IF (n_moist .ge. PARAM_FIRST_SCALAR) THEN
2198
2199 moist_loop: do im = PARAM_FIRST_SCALAR, n_moist
2200
2201 CALL horizontal_diffusion_s( moist_tendf(ims,kms,jms,im), &
2202 mu, config_flags, &
2203 moist(ims,kms,jms,im), &
2204 msftx, msfty, msfux, msfuy, &
2205 msfvx, msfvy, xkhh, rdx, rdy, &
2206 fnm, fnp, cf1, cf2, cf3, &
2207 ht, dx, dy, dt, zx, zy, rdz, rdzw, dnw, dn, &
2208 .false., 'M', &
2209 ids, ide, jds, jde, kds, kde, &
2210 ims, ime, jms, jme, kms, kme, &
2211 its, ite, jts, jte, kts, kte )
2212
2213 ENDDO moist_loop
2214
2215 ENDIF
2216
2217 IF (n_chem .ge. PARAM_FIRST_SCALAR) THEN
2218
2219 chem_loop: do ic = PARAM_FIRST_SCALAR, n_chem
2220
2221 CALL horizontal_diffusion_s( chem_tendf(ims,kms,jms,ic), &
2222 mu, config_flags, &
2223 chem(ims,kms,jms,ic), &
2224 msftx, msfty, msfux, msfuy, &
2225 msfvx, msfvy, xkhh, rdx, rdy, &
2226 fnm, fnp, cf1, cf2, cf3, &
2227 ht, dx, dy, dt, zx, zy, rdz, rdzw, dnw, dn, &
2228 .false., 'S', &
2229 ids, ide, jds, jde, kds, kde, &
2230 ims, ime, jms, jme, kms, kme, &
2231 its, ite, jts, jte, kts, kte )
2232
2233 ENDDO chem_loop
2234
2235 ENDIF
2236 IF (n_scalar .ge. PARAM_FIRST_SCALAR) THEN
2237
2238 scalar_loop: do is = PARAM_FIRST_SCALAR, n_scalar
2239
2240 CALL horizontal_diffusion_s( scalar_tendf(ims,kms,jms,is), &
2241 mu, config_flags, &
2242 scalar(ims,kms,jms,is), &
2243 msftx, msfty, msfux, msfuy, &
2244 msfvx, msfvy, xkhh, rdx, rdy, &
2245 fnm, fnp, cf1, cf2, cf3, &
2246 ht, dx, dy, dt, zx, zy, rdz, rdzw, dnw, dn, &
2247 .false., 'S', &
2248 ids, ide, jds, jde, kds, kde, &
2249 ims, ime, jms, jme, kms, kme, &
2250 its, ite, jts, jte, kts, kte )
2251
2252 ENDDO scalar_loop
2253
2254 ENDIF
2255
2256 END SUBROUTINE horizontal_diffusion_2
2257
2258 !=======================================================================
2259 !=======================================================================
2260
2261 SUBROUTINE horizontal_diffusion_u_2( tendency, mu, config_flags, &
2262 defor11, defor12, div, tke, &
2263 msfux, msfuy, &
2264 xkmh, rdx, rdy, fnm, fnp, &
2265 zx, zy, rdzw, dx, dy, dt, &
2266 ids, ide, jds, jde, kds, kde, &
2267 ims, ime, jms, jme, kms, kme, &
2268 its, ite, jts, jte, kts, kte )
2269
2270 !-----------------------------------------------------------------------
2271 ! Begin declarations.
2272
2273 IMPLICIT NONE
2274
2275 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2276
2277 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2278 ims, ime, jms, jme, kms, kme, &
2279 its, ite, jts, jte, kts, kte
2280
2281 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
2282 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
2283
2284 REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfux, &
2285 msfuy, &
2286 mu
2287
2288 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
2289
2290 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: rdzw
2291
2292
2293 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor11, &
2294 defor12, &
2295 div, &
2296 tke, &
2297 xkmh, &
2298 zx, &
2299 zy
2300
2301 REAL , INTENT(IN ) :: rdx, &
2302 rdy,dx, dy, dt
2303 ! Local data
2304
2305 INTEGER :: i, j, k, ktf
2306
2307 INTEGER :: i_start, i_end, j_start, j_end
2308 INTEGER :: is_ext,ie_ext,js_ext,je_ext
2309
2310 REAL , DIMENSION( ims:ime, kms:kme, jms:jme) :: xkxloc
2311
2312 REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau1avg, &
2313 titau2avg, &
2314 titau1, &
2315 titau2
2316 ! new
2317 ! zxavg, &
2318 ! zyavg
2319 REAL :: mrdx, mrdy, rcoup
2320
2321 REAL :: tmpzy, tmpzeta_z
2322
2323 REAL :: term1, term2, term3
2324
2325 ! End declarations.
2326 !-----------------------------------------------------------------------
2327
2328 ktf=MIN(kte,kde-1)
2329
2330 !-----------------------------------------------------------------------
2331 ! u : p (.), u(|), w(-)
2332 !
2333 ! p u p u u u
2334 !
2335 ! p | . | . | . | k+1 | . | . | . | k+1
2336 !
2337 ! w - 13 - - k+1 13 k+1
2338 !
2339 ! p | 11 O 11 | . | k | 12 O 12 | . | k
2340 !
2341 ! w - 13 - - k 13 k
2342 !
2343 ! p | . | . | . | k-1 | . | . | . | k-1
2344 !
2345 ! i-1 i i i+1 j-1 j j j+1 j+1
2346 !
2347
2348 i_start = its
2349 i_end = ite
2350 j_start = jts
2351 j_end = MIN(jte,jde-1)
2352
2353 IF ( config_flags%open_xs .or. config_flags%specified .or. &
2354 config_flags%nested) i_start = MAX(ids+1,its)
2355 IF ( config_flags%open_xe .or. config_flags%specified .or. &
2356 config_flags%nested) i_end = MIN(ide-1,ite)
2357 IF ( config_flags%open_ys .or. config_flags%specified .or. &
2358 config_flags%nested) j_start = MAX(jds+1,jts)
2359 IF ( config_flags%open_ye .or. config_flags%specified .or. &
2360 config_flags%nested) j_end = MIN(jde-2,jte)
2361 IF ( config_flags%periodic_x ) i_start = its
2362 IF ( config_flags%periodic_x ) i_end = ite
2363
2364 ! titau1 = titau11
2365 is_ext=1
2366 ie_ext=0
2367 js_ext=0
2368 je_ext=0
2369
2370 ! Limit the diffusion coefficient to 40% of the theoretical stability limit
2371 DO j = j_start-1, j_end+1
2372 DO k = kts,ktf
2373 DO i = i_start-1, i_end+1
2374 xkxloc(i,k,j) = min(xkmh(i,k,j),0.05*dx*dy/dt*(1./rdzw(i,k,j)/max(1./rdzw(i,k,j),&
2375 abs(zx(i,k,j))*dx,abs(zx(i+1,k,j))*dx,abs(zx(i,k,j-1))*dx,abs(zx(i-1,k,j))*dx,abs(zx(i,k,j+1))*dx,&
2376 abs(zy(i,k,j))*dy,abs(zy(i-1,k,j))*dy,abs(zy(i,k,j+1))*dy,abs(zy(i-1,k,j+1))*dy))**2)
2377 ENDDO
2378 ENDDO
2379 ENDDO
2380
2381 CALL cal_titau_11_22_33( config_flags, titau1, &
2382 mu, tke, xkxloc, defor11, &
2383 is_ext, ie_ext, js_ext, je_ext, &
2384 ids, ide, jds, jde, kds, kde, &
2385 ims, ime, jms, jme, kms, kme, &
2386 its, ite, jts, jte, kts, kte )
2387
2388 ! titau2 = titau12
2389 is_ext=0
2390 ie_ext=0
2391 js_ext=0
2392 je_ext=1
2393
2394
2395 CALL cal_titau_12_21( config_flags, titau2, &
2396 mu, xkxloc, defor12, &
2397 is_ext, ie_ext, js_ext, je_ext, &
2398 ids, ide, jds, jde, kds, kde, &
2399 ims, ime, jms, jme, kms, kme, &
2400 its, ite, jts, jte, kts, kte )
2401
2402 ! titau1avg = titau11avg
2403 ! titau2avg = titau12avg
2404
2405 DO j = j_start, j_end
2406 DO k = kts+1,ktf
2407 DO i = i_start, i_end
2408 titau1avg(i,k,j)=0.5*(fnm(k)*(titau1(i-1,k ,j)+titau1(i,k ,j))+ &
2409 fnp(k)*(titau1(i-1,k-1,j)+titau1(i,k-1,j)))
2410 titau2avg(i,k,j)=0.5*(fnm(k)*(titau2(i,k ,j+1)+titau2(i,k ,j))+ &
2411 fnp(k)*(titau2(i,k-1,j+1)+titau2(i,k-1,j)))
2412 tmpzy = 0.25*( zy(i-1,k,j )+zy(i,k,j )+ &
2413 zy(i-1,k,j+1)+zy(i,k,j+1) )
2414 ! tmpzeta_z = 0.5*(zeta_z(i,j)+zeta_z(i-1,j))
2415 ! titau1avg(i,k,j)=titau1avg(i,k,j)*zx(i,k,j)*tmpzeta_z
2416 ! titau2avg(i,k,j)=titau2avg(i,k,j)*tmpzy *tmpzeta_z
2417
2418 titau1avg(i,k,j)=titau1avg(i,k,j)*zx(i,k,j)
2419 titau2avg(i,k,j)=titau2avg(i,k,j)*tmpzy
2420
2421 ENDDO
2422 ENDDO
2423 ENDDO
2424 !
2425 DO j = j_start, j_end
2426 DO i = i_start, i_end
2427 titau1avg(i,kts,j)=0.
2428 titau1avg(i,ktf+1,j)=0.
2429 titau2avg(i,kts,j)=0.
2430 titau2avg(i,ktf+1,j)=0.
2431 ENDDO
2432 ENDDO
2433 !
2434 DO j = j_start, j_end
2435 DO k = kts,ktf
2436 DO i = i_start, i_end
2437
2438 mrdx=msfux(i,j)*rdx
2439 mrdy=msfuy(i,j)*rdy
2440 tendency(i,k,j)=tendency(i,k,j)- &
2441 (mrdx*(titau1(i,k,j )-titau1(i-1,k,j))+ &
2442 mrdy*(titau2(i,k,j+1)-titau2(i,k,j ))- &
2443 msfuy(i,j)*rdzw(i,k,j)*((titau1avg(i,k+1,j)-titau1avg(i,k,j))+ &
2444 (titau2avg(i,k+1,j)-titau2avg(i,k,j)) &
2445 ) )
2446 ENDDO
2447 ENDDO
2448 ENDDO
2449
2450 END SUBROUTINE horizontal_diffusion_u_2
2451
2452 !=======================================================================
2453 !=======================================================================
2454
2455 SUBROUTINE horizontal_diffusion_v_2( tendency, mu, config_flags, &
2456 defor12, defor22, div, tke, &
2457 msfvx, msfvy, &
2458 xkmh, rdx, rdy, fnm, fnp, &
2459 zx, zy, rdzw, dx, dy, dt, &
2460 ids, ide, jds, jde, kds, kde, &
2461 ims, ime, jms, jme, kms, kme, &
2462 its, ite, jts, jte, kts, kte )
2463
2464 !-----------------------------------------------------------------------
2465 ! Begin declarations.
2466
2467 IMPLICIT NONE
2468
2469 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2470
2471 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2472 ims, ime, jms, jme, kms, kme, &
2473 its, ite, jts, jte, kts, kte
2474
2475 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
2476 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
2477
2478 REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfvx, &
2479 msfvy, &
2480 mu
2481
2482 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: tendency
2483
2484 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor12, &
2485 defor22, &
2486 div, &
2487 tke, &
2488 xkmh, &
2489 zx, &
2490 zy, &
2491 rdzw
2492
2493 REAL , INTENT(IN ) :: rdx, &
2494 rdy,dx, dy, dt
2495
2496 ! Local data
2497
2498 INTEGER :: i, j, k, ktf
2499
2500 INTEGER :: i_start, i_end, j_start, j_end
2501 INTEGER :: is_ext,ie_ext,js_ext,je_ext
2502
2503 REAL , DIMENSION( ims:ime, kms:kme, jms:jme) :: xkxloc
2504
2505 REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau1avg, &
2506 titau2avg, &
2507 titau1, &
2508 titau2
2509 ! new
2510 ! zxavg, &
2511 ! zyavg
2512
2513 REAL :: mrdx, mrdy, rcoup
2514
2515 REAL :: tmpzx, tmpzeta_z
2516
2517 ! End declarations.
2518 !-----------------------------------------------------------------------
2519
2520 ktf=MIN(kte,kde-1)
2521
2522 !-----------------------------------------------------------------------
2523 ! v : p (.), v(+), w(-)
2524 !
2525 ! p v p v v v
2526 !
2527 ! p + . + . + . + k+1 + . + . + . + k+1
2528 !
2529 ! w - 23 - - k+1 23 k+1
2530 !
2531 ! p + 22 O 22 + . + k + 21 O 21 + . + k
2532 !
2533 ! w - 23 - - k 23 k
2534 !
2535 ! p + . + . + . + k-1 + . + . + . + k-1
2536 !
2537 ! j-1 j j j+1 i-1 i i i+1 i+1
2538 !
2539
2540 i_start = its
2541 i_end = MIN(ite,ide-1)
2542 j_start = jts
2543 j_end = jte
2544
2545 IF ( config_flags%open_xs .or. config_flags%specified .or. &
2546 config_flags%nested) i_start = MAX(ids+1,its)
2547 IF ( config_flags%open_xe .or. config_flags%specified .or. &
2548 config_flags%nested) i_end = MIN(ide-2,ite)
2549 IF ( config_flags%open_ys .or. config_flags%specified .or. &
2550 config_flags%nested) j_start = MAX(jds+1,jts)
2551 IF ( config_flags%open_ye .or. config_flags%specified .or. &
2552 config_flags%nested) j_end = MIN(jde-1,jte)
2553 IF ( config_flags%periodic_x ) i_start = its
2554 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2555
2556 ! titau1 = titau21
2557 is_ext=0
2558 ie_ext=1
2559 js_ext=0
2560 je_ext=0
2561
2562 ! Limit the diffusion coefficient to 40% of the theoretical stability limit
2563 DO j = j_start-1, j_end+1
2564 DO k = kts,ktf
2565 DO i = i_start-1, i_end+1
2566 xkxloc(i,k,j) = min(xkmh(i,k,j),0.05*dx*dy/dt*(1./rdzw(i,k,j)/max(1./rdzw(i,k,j),&
2567 abs(zx(i,k,j))*dx,abs(zx(i+1,k,j))*dx,abs(zx(i,k,j-1))*dx,abs(zx(i+1,k,j-1))*dx,&
2568 abs(zy(i,k,j))*dy,abs(zy(i-1,k,j))*dy,abs(zy(i,k,j+1))*dy,abs(zy(i,k,j-1))*dy,abs(zy(i+1,k,j))*dy))**2)
2569 ENDDO
2570 ENDDO
2571 ENDDO
2572
2573 CALL cal_titau_12_21( config_flags, titau1, &
2574 mu, xkxloc, defor12, &
2575 is_ext,ie_ext,js_ext,je_ext, &
2576 ids, ide, jds, jde, kds, kde, &
2577 ims, ime, jms, jme, kms, kme, &
2578 its, ite, jts, jte, kts, kte )
2579
2580 ! titau2 = titau22
2581 is_ext=0
2582 ie_ext=0
2583 js_ext=1
2584 je_ext=0
2585
2586 CALL cal_titau_11_22_33( config_flags, titau2, &
2587 mu, tke, xkxloc, defor22, &
2588 is_ext, ie_ext, js_ext, je_ext, &
2589 ids, ide, jds, jde, kds, kde, &
2590 ims, ime, jms, jme, kms, kme, &
2591 its, ite, jts, jte, kts, kte )
2592
2593 DO j = j_start, j_end
2594 DO k = kts+1,ktf
2595 DO i = i_start, i_end
2596 titau1avg(i,k,j)=0.5*(fnm(k)*(titau1(i+1,k ,j)+titau1(i,k ,j))+ &
2597 fnp(k)*(titau1(i+1,k-1,j)+titau1(i,k-1,j)))
2598 titau2avg(i,k,j)=0.5*(fnm(k)*(titau2(i,k ,j-1)+titau2(i,k ,j))+ &
2599 fnp(k)*(titau2(i,k-1,j-1)+titau2(i,k-1,j)))
2600
2601 tmpzx = 0.25*( zx(i,k,j )+zx(i+1,k,j )+ &
2602 zx(i,k,j-1)+zx(i+1,k,j-1) )
2603
2604
2605 titau1avg(i,k,j)=titau1avg(i,k,j)*tmpzx
2606 titau2avg(i,k,j)=titau2avg(i,k,j)*zy(i,k,j)
2607
2608
2609 ENDDO
2610 ENDDO
2611 ENDDO
2612
2613 DO j = j_start, j_end
2614 DO i = i_start, i_end
2615 titau1avg(i,kts,j)=0.
2616 titau1avg(i,ktf+1,j)=0.
2617 titau2avg(i,kts,j)=0.
2618 titau2avg(i,ktf+1,j)=0.
2619 ENDDO
2620 ENDDO
2621 !
2622 DO j = j_start, j_end
2623 DO k = kts,ktf
2624 DO i = i_start, i_end
2625
2626 mrdx=msfvx(i,j)*rdx
2627 mrdy=msfvy(i,j)*rdy
2628 tendency(i,k,j)=tendency(i,k,j)- &
2629 (mrdy*(titau2(i ,k,j)-titau2(i,k,j-1))+ &
2630 mrdx*(titau1(i+1,k,j)-titau1(i,k,j ))- &
2631 msfvy(i,j)*rdzw(i,k,j)*((titau1avg(i,k+1,j)-titau1avg(i,k,j))+ &
2632 (titau2avg(i,k+1,j)-titau2avg(i,k,j)) &
2633 ) &
2634 )
2635
2636 ENDDO
2637 ENDDO
2638 ENDDO
2639
2640 END SUBROUTINE horizontal_diffusion_v_2
2641
2642 !=======================================================================
2643 !=======================================================================
2644
2645 SUBROUTINE horizontal_diffusion_w_2( tendency, mu, config_flags, &
2646 defor13, defor23, div, tke, &
2647 msftx, msfty, &
2648 xkmh, rdx, rdy, fnm, fnp, &
2649 zx, zy, rdz, dx, dy, dt, &
2650 ids, ide, jds, jde, kds, kde, &
2651 ims, ime, jms, jme, kms, kme, &
2652 its, ite, jts, jte, kts, kte )
2653
2654 !-----------------------------------------------------------------------
2655 ! Begin declarations.
2656
2657 IMPLICIT NONE
2658
2659 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2660
2661 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2662 ims, ime, jms, jme, kms, kme, &
2663 its, ite, jts, jte, kts, kte
2664
2665 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
2666 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
2667
2668 REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msftx, &
2669 msfty, &
2670 mu
2671
2672 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: tendency
2673
2674 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor13, &
2675 defor23, &
2676 div, &
2677 tke, &
2678 xkmh, &
2679 zx, &
2680 zy, &
2681 rdz
2682
2683 REAL , INTENT(IN ) :: rdx, &
2684 rdy,dx, dy, dt
2685
2686 ! Local data
2687
2688 INTEGER :: i, j, k, ktf
2689
2690 INTEGER :: i_start, i_end, j_start, j_end
2691 INTEGER :: is_ext,ie_ext,js_ext,je_ext
2692
2693 REAL , DIMENSION( ims:ime, kms:kme, jms:jme) :: xkxloc
2694
2695 REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau1avg, &
2696 titau2avg, &
2697 titau1, &
2698 titau2
2699 ! new
2700 ! zxavg, &
2701 ! zyavg
2702
2703 REAL :: mrdx, mrdy, rcoup
2704
2705 REAL :: tmpzx, tmpzy, tmpzeta_z
2706
2707 ! End declarations.
2708 !-----------------------------------------------------------------------
2709
2710 ktf=MIN(kte,kde-1)
2711
2712 !-----------------------------------------------------------------------
2713 ! w : p (.), u(|), v(+), w(-)
2714 !
2715 ! p u p u p v p v
2716 !
2717 ! w - - - k+1 w - - - k+1
2718 !
2719 ! p . | 33 | . k p . + 33 + . k
2720 !
2721 ! w - 31 O 31 - k w - 32 O 32 - k
2722 !
2723 ! p . | 33 | . k-1 p . | 33 | . k-1
2724 !
2725 ! w - - - k-1 w - - - k-1
2726 !
2727 ! i-1 i i i+1 j-1 j j j+1
2728 !
2729 i_start = its
2730 i_end = MIN(ite,ide-1)
2731 j_start = jts
2732 j_end = MIN(jte,jde-1)
2733
2734 IF ( config_flags%open_xs .or. config_flags%specified .or. &
2735 config_flags%nested) i_start = MAX(ids+1,its)
2736 IF ( config_flags%open_xe .or. config_flags%specified .or. &
2737 config_flags%nested) i_end = MIN(ide-2,ite)
2738 IF ( config_flags%open_ys .or. config_flags%specified .or. &
2739 config_flags%nested) j_start = MAX(jds+1,jts)
2740 IF ( config_flags%open_ye .or. config_flags%specified .or. &
2741 config_flags%nested) j_end = MIN(jde-2,jte)
2742 IF ( config_flags%periodic_x ) i_start = its
2743 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2744
2745 ! titau1 = titau31
2746 is_ext=0
2747 ie_ext=1
2748 js_ext=0
2749 je_ext=0
2750
2751 ! Limit the diffusion coefficient to 40% of the theoretical stability limit
2752 DO j = j_start-1, j_end+1
2753 DO k = kts,ktf
2754 DO i = i_start-1, i_end+1
2755 xkxloc(i,k,j) = min(xkmh(i,k,j),0.05*dx*dy/dt*(1./rdz(i,k,j))/max(1./rdz(i,k,j),&
2756 abs(zx(i,k,j))*dx,abs(zx(i+1,k,j))*dx,&
2757 abs(zy(i,k,j))*dy,abs(zy(i,k,j+1))*dy)**2 )
2758 ENDDO
2759 ENDDO
2760 ENDDO
2761
2762 CALL cal_titau_13_31( config_flags, titau1, defor13, &
2763 mu, xkxloc, fnm, fnp, &
2764 is_ext, ie_ext, js_ext, je_ext, &
2765 ids, ide, jds, jde, kds, kde, &
2766 ims, ime, jms, jme, kms, kme, &
2767 its, ite, jts, jte, kts, kte )
2768
2769 ! titau2 = titau32
2770 is_ext=0
2771 ie_ext=0
2772 js_ext=0
2773 je_ext=1
2774
2775 CALL cal_titau_23_32( config_flags, titau2, defor23, &
2776 mu, xkxloc, fnm, fnp, &
2777 is_ext, ie_ext, js_ext, je_ext, &
2778 ids, ide, jds, jde, kds, kde, &
2779 ims, ime, jms, jme, kms, kme, &
2780 its, ite, jts, jte, kts, kte )
2781
2782 ! titau1avg = titau31avg * zx * zeta_z = titau13avg * zx * zeta_z
2783 ! titau2avg = titau32avg * zy * zeta_z = titau23avg * zy * zeta_z
2784
2785 DO j = j_start, j_end
2786 DO k = kts,ktf
2787 DO i = i_start, i_end
2788 titau1avg(i,k,j)=0.25*(titau1(i+1,k+1,j)+titau1(i,k+1,j)+ &
2789 titau1(i+1,k ,j)+titau1(i,k ,j))
2790 titau2avg(i,k,j)=0.25*(titau2(i,k+1,j+1)+titau2(i,k+1,j)+ &
2791 titau2(i,k ,j+1)+titau2(i,k ,j))
2792 ! new
2793 tmpzx =0.25*( zx(i,k ,j)+zx(i+1,k ,j)+ &
2794 zx(i,k+1,j)+zx(i+1,k+1,j) )
2795 tmpzy =0.25*( zy(i,k ,j)+zy(i,k ,j+1)+ &
2796 zy(i,k+1,j)+zy(i,k+1,j+1) )
2797
2798 titau1avg(i,k,j)=titau1avg(i,k,j)*tmpzx
2799 titau2avg(i,k,j)=titau2avg(i,k,j)*tmpzy
2800 ! titau1avg(i,k,j)=titau1avg(i,k,j)*tmpzx*zeta_z(i,j)
2801 ! titau2avg(i,k,j)=titau2avg(i,k,j)*tmpzy*zeta_z(i,j)
2802 ENDDO
2803 ENDDO
2804 ENDDO
2805
2806 DO j = j_start, j_end
2807 DO i = i_start, i_end
2808 titau1avg(i,ktf+1,j)=0.
2809 titau2avg(i,ktf+1,j)=0.
2810 ENDDO
2811 ENDDO
2812
2813 DO j = j_start, j_end
2814 DO k = kts+1,ktf
2815 DO i = i_start, i_end
2816
2817 mrdx=msftx(i,j)*rdx
2818 mrdy=msfty(i,j)*rdy
2819
2820 tendency(i,k,j)=tendency(i,k,j)- &
2821 (mrdx*(titau1(i+1,k,j)-titau1(i,k,j))+ &
2822 mrdy*(titau2(i,k,j+1)-titau2(i,k,j))- &
2823 msfty(i,j)*rdz(i,k,j)*(titau1avg(i,k,j)-titau1avg(i,k-1,j)+ &
2824 titau2avg(i,k,j)-titau2avg(i,k-1,j) &
2825 ) &
2826 )
2827 ! msft(i,j)/dn(k)*(titau1avg(i,k,j)-titau1avg(i,k-1,j)+ &
2828 ! titau2avg(i,k,j)-titau2avg(i,k-1,j) &
2829 ! ) &
2830 ! )
2831 ENDDO
2832 ENDDO
2833 ENDDO
2834
2835 END SUBROUTINE horizontal_diffusion_w_2
2836
2837 !=======================================================================
2838 !=======================================================================
2839
2840 SUBROUTINE horizontal_diffusion_s (tendency, mu, config_flags, var, &
2841 msftx, msfty, msfux, msfuy, &
2842 msfvx, msfvy, xkhh, rdx, rdy, &
2843 fnm, fnp, cf1, cf2, cf3, &
2844 ht, dx, dy, dt, zx, zy, rdz, rdzw, dn, dnw, &
2845 doing_tke, charvar, &
2846 ids, ide, jds, jde, kds, kde, &
2847 ims, ime, jms, jme, kms, kme, &
2848 its, ite, jts, jte, kts, kte )
2849
2850 !-----------------------------------------------------------------------
2851 ! Begin declarations.
2852
2853 IMPLICIT NONE
2854
2855 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2856
2857 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2858 ims, ime, jms, jme, kms, kme, &
2859 its, ite, jts, jte, kts, kte
2860
2861 LOGICAL, INTENT(IN ) :: doing_tke
2862
2863 REAL , INTENT(IN ) :: cf1, cf2, cf3
2864
2865 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
2866 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
2867 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
2868 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
2869
2870 REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfux
2871 REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfuy
2872 REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfvx
2873 REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfvy
2874 REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msftx
2875 REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfty
2876
2877 REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: mu, ht
2878
2879 ! REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1), &
2880 ! INTENT(IN ) :: xkhh
2881
2882 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: tendency
2883
2884 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: &
2885 xkhh, &
2886 rdz, &
2887 rdzw
2888
2889 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: var, &
2890 zx, &
2891 zy
2892
2893 REAL , INTENT(IN ) :: rdx, &
2894 rdy,dx, dy, dt
2895 character*1 charvar
2896
2897 ! Local data
2898
2899 INTEGER :: i, j, k, ktf
2900
2901 INTEGER :: i_start, i_end, j_start, j_end
2902
2903 REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: H1avg, &
2904 H2avg, &
2905 H1, &
2906 H2, &
2907 xkxavg
2908 ! new
2909 ! zxavg, &
2910 ! zyavg
2911
2912 REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: tmptendf
2913
2914 REAL :: mrdx, mrdy, rcoup
2915 REAL :: tmpzx, tmpzy, tmpzeta_z,tgrmin,tgrmax
2916 INTEGER :: ktes1,ktes2,ii,jj
2917 logical :: limtgrad
2918
2919 ! End declarations.
2920 !-----------------------------------------------------------------------
2921
2922 ktf=MIN(kte,kde-1)
2923
2924 !-----------------------------------------------------------------------
2925 ! scalars: t (.), u(|), v(+), w(-)
2926 !
2927 ! t u t u t v t v
2928 !
2929 ! w - 3 - k+1 w - 3 - k+1
2930 !
2931 ! t . 1 O 1 . k t . 2 O 2 . k
2932 !
2933 ! w - 3 - k w - 3 - k
2934 !
2935 ! t . | . | . k-1 t . + . + . k-1
2936 !
2937 ! w - - - k-1 w - - - k-1
2938 !
2939 ! t i-1 i i i+1 j-1 j j j+1
2940 !
2941
2942 ktes1=kte-1
2943 ktes2=kte-2
2944
2945 i_start = its
2946 i_end = MIN(ite,ide-1)
2947 j_start = jts
2948 j_end = MIN(jte,jde-1)
2949
2950 limtgrad = .true.
2951
2952 IF ( config_flags%open_xs .or. config_flags%specified .or. &
2953 config_flags%nested) i_start = MAX(ids+1,its)
2954 IF ( config_flags%open_xe .or. config_flags%specified .or. &
2955 config_flags%nested) i_end = MIN(ide-2,ite)
2956 IF ( config_flags%open_ys .or. config_flags%specified .or. &
2957 config_flags%nested) j_start = MAX(jds+1,jts)
2958 IF ( config_flags%open_ye .or. config_flags%specified .or. &
2959 config_flags%nested) j_end = MIN(jde-2,jte)
2960 IF ( config_flags%periodic_x ) i_start = its
2961 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2962
2963 if (limtgrad) then ! limit vertical temperature gradient in metric term to reduce possible overcorrection
2964 tgrmin = -0.005 ! -0.5K/100m
2965 tgrmax = 0.08 ! +8K/100m
2966 else
2967 tgrmin = -100.
2968 tgrmax = 100.
2969 endif
2970
2971 ! diffusion of the TKE needs mutiple 2
2972
2973 IF ( doing_tke ) THEN
2974 DO j = j_start, j_end
2975 DO k = kts,ktf
2976 DO i = i_start, i_end
2977 tmptendf(i,k,j)=tendency(i,k,j)
2978 ENDDO
2979 ENDDO
2980 ENDDO
2981 ENDIF
2982
2983 ! H1 = partial var over partial x
2984
2985 ! Limit the diffusion coefficient to 40% of the theoretical stability limit
2986 DO j = j_start, j_end
2987 DO k = kts, ktf
2988 DO i = i_start, i_end + 1
2989 ! new
2990 ! zxavg(i,k,j) =0.5*( zx(i-1,k,j)+ zx(i,k,j))
2991 ! xkxavg(i,k,j)=0.5*(xkhh(i-1,k,j)+xkhh(i,k,j))
2992 xkxavg(i,k,j)=min(0.5*(xkhh(i-1,k,j)+xkhh(i,k,j)),0.05*dx*dx/dt*0.5*(1./rdzw(i,k,j)+1./rdzw(i-1,k,j))/ &
2993 max(0.5*(1./rdzw(i,k,j)+1./rdzw(i-1,k,j)),abs(zx(i,k,j))*dx,abs(zx(i+1,k,j))*dx )**2)
2994
2995 ENDDO
2996 ENDDO
2997 ENDDO
2998
2999 DO j = j_start, j_end
3000 DO k = kts+1, ktf
3001 DO i = i_start, i_end + 1
3002 H1avg(i,k,j)=0.5*(fnm(k)*(var(i-1,k ,j)+var(i,k ,j))+ &
3003 fnp(k)*(var(i-1,k-1,j)+var(i,k-1,j)))
3004 ENDDO
3005 ENDDO
3006 ENDDO
3007
3008
3009 DO j = j_start, j_end
3010 DO i = i_start, i_end + 1
3011 H1avg(i,kts ,j)=0.5*(cf1*var(i ,1,j)+cf2*var(i ,2,j)+ &
3012 cf3*var(i ,3,j)+cf1*var(i-1,1,j)+ &
3013 cf2*var(i-1,2,j)+cf3*var(i-1,3,j))
3014 H1avg(i,ktf+1,j)=0.5*(var(i,ktes1,j)+(var(i,ktes1,j)- &
3015 var(i,ktes2,j))*0.5*dnw(ktes1)/dn(ktes1)+ &
3016 var(i-1,ktes1,j)+(var(i-1,ktes1,j)- &
3017 var(i-1,ktes2,j))*0.5*dnw(ktes1)/dn(ktes1))
3018 ENDDO
3019 ENDDO
3020
3021 if (charvar.eq.'T') then ! use temperature gradient limiter to avoid instability
3022
3023 DO j = j_start, j_end
3024 DO k = kts, ktf
3025 DO i = i_start, i_end + 1
3026 tmpzx = 0.5*( zx(i,k,j)+ zx(i,k+1,j))
3027 H1(i,k,j)=-msfuy(i,j)*xkxavg(i,k,j)*( &
3028 rdx*(var(i,k,j)-var(i-1,k,j)) - tmpzx* &
3029 max(tgrmin,min(tgrmax,(H1avg(i,k+1,j)-H1avg(i,k,j))*rdzw(i,k,j))) )
3030 ENDDO
3031 ENDDO
3032 ENDDO
3033
3034 else
3035 DO j = j_start, j_end
3036 DO k = kts, ktf
3037 DO i = i_start, i_end + 1
3038 ! new
3039 tmpzx = 0.5*( zx(i,k,j)+ zx(i,k+1,j))
3040 H1(i,k,j)=-msfuy(i,j)*xkxavg(i,k,j)*( &
3041 rdx*(var(i,k,j)-var(i-1,k,j)) - tmpzx* &
3042 (H1avg(i,k+1,j)-H1avg(i,k,j))*rdzw(i,k,j) )
3043
3044 ! tmpzeta_z = 0.5*(zeta_z(i,j)+zeta_z(i-1,j))
3045 ! H1(i,k,j)=-msfuy(i,j)*xkxavg(i,k,j)*( &
3046 ! rdx*(var(i,k,j)-var(i-1,k,j)) - tmpzx*tmpzeta_z* &
3047 ! (H1avg(i,k+1,j)-H1avg(i,k,j))/dnw(k))
3048 ENDDO
3049 ENDDO
3050 ENDDO
3051 endif
3052
3053 ! H2 = partial var over partial y
3054
3055
3056 ! Limit the diffusion coefficient to 40% of the theoretical stability limit
3057 DO j = j_start, j_end + 1
3058 DO k = kts, ktf
3059 DO i = i_start, i_end
3060 ! new
3061 ! zyavg(i,k,j) =0.5*( zy(i,k,j-1)+ zy(i,k,j))
3062 ! xkxavg(i,k,j)=0.5*(xkhh(i,k,j-1)+xkhh(i,k,j))
3063 xkxavg(i,k,j)=min(0.5*(xkhh(i,k,j-1)+xkhh(i,k,j)),0.05*dy*dy/dt*0.5*(1./rdzw(i,k,j)+1./rdzw(i,k,j-1))/ &
3064 max(0.5*(1./rdzw(i,k,j)+1./rdzw(i,k,j-1)),abs(zy(i,k,j))*dy,abs(zy(i,k,j+1))*dy )**2)
3065
3066 ENDDO
3067 ENDDO
3068 ENDDO
3069
3070 DO j = j_start, j_end + 1
3071 DO k = kts+1, ktf
3072 DO i = i_start, i_end
3073 ! new
3074 H2avg(i,k,j)=0.5*(fnm(k)*(var(i,k ,j-1)+var(i,k ,j))+ &
3075 fnp(k)*(var(i,k-1,j-1)+var(i,k-1,j)))
3076 ENDDO
3077 ENDDO
3078 ENDDO
3079
3080
3081 DO j = j_start, j_end + 1
3082 DO i = i_start, i_end
3083 H2avg(i,kts ,j)=0.5*(cf1*var(i,1,j )+cf2*var(i ,2,j)+ &
3084 cf3*var(i,3,j )+cf1*var(i,1,j-1)+ &
3085 cf2*var(i,2,j-1)+cf3*var(i,3,j-1))
3086 H2avg(i,ktf+1,j)=0.5*(var(i,ktes1,j)+(var(i,ktes1,j)- &
3087 var(i,ktes2,j))*0.5*dnw(ktes1)/dn(ktes1)+ &
3088 var(i,ktes1,j-1)+(var(i,ktes1,j-1)- &
3089 var(i,ktes2,j-1))*0.5*dnw(ktes1)/dn(ktes1))
3090 ENDDO
3091 ENDDO
3092
3093 if (charvar.eq.'T') then ! use temperature gradient limiter to avoid instability
3094
3095 DO j = j_start, j_end + 1
3096 DO k = kts, ktf
3097 DO i = i_start, i_end
3098 tmpzy = 0.5*( zy(i,k,j)+ zy(i,k+1,j))
3099
3100 H2(i,k,j)=-msfvy(i,j)*xkxavg(i,k,j)*( &
3101 rdy*(var(i,k,j)-var(i,k,j-1)) - tmpzy* &
3102 max(tgrmin,min(tgrmax,(H2avg(i ,k+1,j)-H2avg(i,k,j))*rdzw(i,k,j))) )
3103 ENDDO
3104 ENDDO
3105 ENDDO
3106
3107 else
3108 DO j = j_start, j_end + 1
3109 DO k = kts, ktf
3110 DO i = i_start, i_end
3111 ! new
3112 tmpzy = 0.5*( zy(i,k,j)+ zy(i,k+1,j))
3113
3114 H2(i,k,j)=-msfvy(i,j)*xkxavg(i,k,j)*( &
3115 rdy*(var(i,k,j)-var(i,k,j-1)) - tmpzy* &
3116 (H2avg(i ,k+1,j)-H2avg(i,k,j))*rdzw(i,k,j))
3117
3118 ! tmpzeta_z = 0.5*(zeta_z(i,j)+zeta_z(i,j-1))
3119 ! H2(i,k,j)=-msfvy(i,j)*xkxavg(i,k,j)*( &
3120 ! rdy*(var(i,k,j)-var(i,k,j-1)) - tmpzy*tmpzeta_z* &
3121 ! (H2avg(i ,k+1,j)-H2avg(i,k,j))/dnw(k))
3122 ENDDO
3123 ENDDO
3124 ENDDO
3125 endif
3126
3127 DO j = j_start, j_end
3128 DO k = kts+1, ktf
3129 DO i = i_start, i_end
3130 H1avg(i,k,j)=0.5*(fnm(k)*(H1(i+1,k ,j)+H1(i,k ,j))+ &
3131 fnp(k)*(H1(i+1,k-1,j)+H1(i,k-1,j)))
3132 H2avg(i,k,j)=0.5*(fnm(k)*(H2(i,k ,j+1)+H2(i,k ,j))+ &
3133 fnp(k)*(H2(i,k-1,j+1)+H2(i,k-1,j)))
3134 ! new
3135 ! zxavg(i,k,j)=fnm(k)*zx(i,k,j)+fnp(k)*zx(i,k-1,j)
3136 ! zyavg(i,k,j)=fnm(k)*zy(i,k,j)+fnp(k)*zy(i,k-1,j)
3137
3138 ! H1avg(i,k,j)=zx*H1avg*zeta_z
3139 ! H2avg(i,k,j)=zy*H2avg*zeta_z
3140
3141 tmpzx = 0.5*( zx(i,k,j)+ zx(i+1,k,j ))
3142 tmpzy = 0.5*( zy(i,k,j)+ zy(i ,k,j+1))
3143
3144 H1avg(i,k,j)=H1avg(i,k,j)*tmpzx
3145 H2avg(i,k,j)=H2avg(i,k,j)*tmpzy
3146
3147 ! H1avg(i,k,j)=H1avg(i,k,j)*tmpzx*zeta_z(i,j)
3148 ! H2avg(i,k,j)=H2avg(i,k,j)*tmpzy*zeta_z(i,j)
3149 ENDDO
3150 ENDDO
3151 ENDDO
3152
3153 DO j = j_start, j_end
3154 DO i = i_start, i_end
3155 H1avg(i,kts ,j)=0.
3156 H1avg(i,ktf+1,j)=0.
3157 H2avg(i,kts ,j)=0.
3158 H2avg(i,ktf+1,j)=0.
3159 ENDDO
3160 ENDDO
3161
3162 DO j = j_start, j_end
3163 DO k = kts,ktf
3164 DO i = i_start, i_end
3165
3166 mrdx=msftx(i,j)*rdx
3167 mrdy=msfty(i,j)*rdy
3168
3169 tendency(i,k,j)=tendency(i,k,j)- &
3170 (mrdx*0.5*((mu(i+1,j)+mu(i,j))*H1(i+1,k,j)- &
3171 (mu(i-1,j)+mu(i,j))*H1(i ,k,j))+ &
3172 mrdy*0.5*((mu(i,j+1)+mu(i,j))*H2(i,k,j+1)- &
3173 (mu(i,j-1)+mu(i,j))*H2(i,k,j ))- &
3174 msfty(i,j)*mu(i,j)*(H1avg(i,k+1,j)-H1avg(i,k,j)+ &
3175 H2avg(i,k+1,j)-H2avg(i,k,j) &
3176 )*rdzw(i,k,j) &
3177 )
3178
3179 ENDDO
3180 ENDDO
3181 ENDDO
3182
3183 IF ( doing_tke ) THEN
3184 DO j = j_start, j_end
3185 DO k = kts,ktf
3186 DO i = i_start, i_end
3187 tendency(i,k,j)=tmptendf(i,k,j)+2.* &
3188 (tendency(i,k,j)-tmptendf(i,k,j))
3189 ENDDO
3190 ENDDO
3191 ENDDO
3192 ENDIF
3193
3194 END SUBROUTINE horizontal_diffusion_s
3195
3196 !=======================================================================
3197 !=======================================================================
3198
3199 SUBROUTINE vertical_diffusion_2 ( ru_tendf, rv_tendf, rw_tendf, rt_tendf, &
3200 tke_tendf, moist_tendf, n_moist, &
3201 chem_tendf, n_chem, &
3202 scalar_tendf, n_scalar, &
3203 u_2, v_2, &
3204 thp,u_base,v_base,t_base,qv_base,mu,tke, &
3205 config_flags,defor13,defor23,defor33,div, &
3206 moist, chem, scalar, xkmv, xkhv,km_opt, &
3207 fnm, fnp, dn, dnw, rdz, rdzw, &
3208 hfx, qfx, ust, rho, &
3209 ids, ide, jds, jde, kds, kde, &
3210 ims, ime, jms, jme, kms, kme, &
3211 its, ite, jts, jte, kts, kte )
3212
3213 !-----------------------------------------------------------------------
3214 ! Begin declarations.
3215
3216 IMPLICIT NONE
3217
3218 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3219
3220 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3221 ims, ime, jms, jme, kms, kme, &
3222 its, ite, jts, jte, kts, kte
3223
3224 INTEGER , INTENT(IN ) :: n_moist, n_chem, n_scalar, km_opt
3225
3226 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
3227 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
3228 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
3229 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
3230 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
3231
3232 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: qv_base
3233 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: u_base
3234 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: v_base
3235 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: t_base
3236
3237 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::ru_tendf,&
3238 rv_tendf,&
3239 rw_tendf,&
3240 tke_tendf,&
3241 rt_tendf
3242
3243 REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), &
3244 INTENT(INOUT) :: moist_tendf
3245
3246 REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_chem), &
3247 INTENT(INOUT) :: chem_tendf
3248
3249 REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_scalar) , &
3250 INTENT(INOUT) :: scalar_tendf
3251
3252 REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), &
3253 INTENT(INOUT) :: moist
3254
3255 REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_chem), &
3256 INTENT(INOUT) :: chem
3257
3258 REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_scalar) , &
3259 INTENT(IN ) :: scalar
3260
3261 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor13, &
3262 defor23, &
3263 defor33, &
3264 div, &
3265 xkmv, &
3266 xkhv, &
3267 tke, &
3268 rdz, &
3269 u_2, &
3270 v_2, &
3271 rdzw
3272
3273 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: rho
3274 REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: hfx, &
3275 qfx
3276 REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: ust
3277 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: thp
3278
3279 ! LOCAL VAR
3280
3281 REAL , DIMENSION( ims:ime, kms:kme, jms:jme) :: var_mix
3282
3283 INTEGER :: im, i,j,k
3284 INTEGER :: i_start, i_end, j_start, j_end
3285
3286 ! REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: xkhv
3287
3288 !***************************************************************************
3289 !***************************************************************************
3290 !MODIFICA VARIABILI PER I FLUSSI
3291 !
3292 REAL :: V0_u,V0_v,tao_xz,tao_yz,ustar,cd0
3293 REAL :: xsfc,psi1,vk2,zrough,lnz
3294 REAL :: heat_flux, moist_flux, heat_flux0
3295 REAL :: cpm
3296 !
3297 !FINE MODIFICA VARIABILI PER I FLUSSI
3298 !***************************************************************************
3299 !
3300
3301 ! End declarations.
3302 !-----------------------------------------------------------------------
3303
3304 i_start = its
3305 i_end = MIN(ite,ide-1)
3306 j_start = jts
3307 j_end = MIN(jte,jde-1)
3308 !
3309 !-----------------------------------------------------------------------
3310
3311 CALL vertical_diffusion_u_2( ru_tendf, config_flags, mu, &
3312 defor13, xkmv, &
3313 dnw, rdzw, fnm, fnp, &
3314 ids, ide, jds, jde, kds, kde, &
3315 ims, ime, jms, jme, kms, kme, &
3316 its, ite, jts, jte, kts, kte )
3317
3318
3319 CALL vertical_diffusion_v_2( rv_tendf, config_flags, mu, &
3320 defor23, xkmv, &
3321 dnw, rdzw, fnm, fnp, &
3322 ids, ide, jds, jde, kds, kde, &
3323 ims, ime, jms, jme, kms, kme, &
3324 its, ite, jts, jte, kts, kte )
3325
3326 CALL vertical_diffusion_w_2( rw_tendf, config_flags, mu, &
3327 defor33, tke(ims,kms,jms), &
3328 div, xkmv, &
3329 dn, rdz, &
3330 ids, ide, jds, jde, kds, kde, &
3331 ims, ime, jms, jme, kms, kme, &
3332 its, ite, jts, jte, kts, kte )
3333
3334 !*****************************************
3335 !*****************************************
3336 ! MODIFICA al flusso di momento alla parete
3337 !
3338 vflux: SELECT CASE( config_flags%isfflx )
3339 CASE (0) ! Assume cd a constant, specified in namelist
3340 cd0 = config_flags%tke_drag_coefficient ! constant drag coefficient
3341 ! set in namelist.input
3342 !
3343 !calcolo del modulo della velocita
3344 DO j = j_start, j_end
3345 DO i = i_start, i_end+1
3346 V0_u=0.
3347 tao_xz=0.
3348 V0_u= sqrt((u_2(i,kts,j)**2) + &
3349 (((v_2(i ,kts,j )+ &
3350 v_2(i ,kts,j+1)+ &
3351 v_2(i-1,kts,j )+ &
3352 v_2(i-1,kts,j+1))/4)**2))+epsilon
3353 tao_xz=cd0*V0_u*u_2(i,kts,j)
3354 ru_tendf(i,kts,j)=ru_tendf(i,kts,j) &
3355 -0.25*(mu(i,j)+mu(i-1,j))*tao_xz*(rdzw(i,kts,j)+rdzw(i-1,kts,j))
3356 ENDDO
3357 ENDDO
3358
3359 !
3360 DO j = j_start, j_end+1
3361 DO i = i_start, i_end
3362 V0_v=0.
3363 tao_yz=0.
3364 V0_v= sqrt((v_2(i,kts,j)**2) + &
3365 (((u_2(i ,kts,j )+ &
3366 u_2(i ,kts,j-1)+ &
3367 u_2(i+1,kts,j )+ &
3368 u_2(i+1,kts,j-1))/4)**2))+epsilon
3369 tao_yz=cd0*V0_v*v_2(i,kts,j)
3370 rv_tendf(i,kts,j)=rv_tendf(i,kts,j) &
3371 -0.25*(mu(i,j)+mu(i,j-1))*tao_yz*(rdzw(i,kts,j)+rdzw(i,kts,j-1))
3372 ENDDO
3373 ENDDO
3374
3375 CASE (1) ! ustar computed from surface routine
3376 DO j = j_start, j_end
3377 DO i = i_start, i_end+1
3378 V0_u=0.
3379 tao_xz=0.
3380 V0_u= sqrt((u_2(i,kts,j)**2) + &
3381 (((v_2(i ,kts,j )+ &
3382 v_2(i ,kts,j+1)+ &
3383 v_2(i-1,kts,j )+ &
3384 v_2(i-1,kts,j+1))/4)**2))+epsilon
3385 tao_xz=ust(i,j)*ust(i,j)*u_2(i,kts,j)/V0_u
3386 ru_tendf(i,kts,j)=ru_tendf(i,kts,j) &
3387 -0.25*(mu(i,j)+mu(i-1,j))*tao_xz*(rdzw(i,kts,j)+rdzw(i-1,kts,j))
3388 ENDDO
3389 ENDDO
3390
3391 !
3392 DO j = j_start, j_end+1
3393 DO i = i_start, i_end
3394 V0_v=0.
3395 tao_yz=0.
3396 V0_v= sqrt((v_2(i,kts,j)**2) + &
3397 (((u_2(i ,kts,j )+ &
3398 u_2(i ,kts,j-1)+ &
3399 u_2(i+1,kts,j )+ &
3400 u_2(i+1,kts,j-1))/4)**2))+epsilon
3401 tao_yz=ust(i,j)*ust(i,j)*v_2(i,kts,j)/V0_v
3402 rv_tendf(i,kts,j)=rv_tendf(i,kts,j) &
3403 -0.25*(mu(i,j)+mu(i,j-1))*tao_yz*(rdzw(i,kts,j)+rdzw(i,kts,j-1))
3404 ENDDO
3405 ENDDO
3406
3407 CASE DEFAULT
3408 CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' )
3409 END SELECT vflux
3410
3411 !
3412 ! FINE MODIFICA al flusso di momento alla parete
3413 !*****************************************
3414 !*****************************************
3415
3416 IF ( config_flags%mix_full_fields ) THEN
3417
3418 DO j=jts,min(jte,jde-1)
3419 DO k=kts,kte-1
3420 DO i=its,min(ite,ide-1)
3421 var_mix(i,k,j) = thp(i,k,j)
3422 ENDDO
3423 ENDDO
3424 ENDDO
3425
3426 ELSE
3427
3428 DO j=jts,min(jte,jde-1)
3429 DO k=kts,kte-1
3430 DO i=its,min(ite,ide-1)
3431 var_mix(i,k,j) = thp(i,k,j) - t_base(k)
3432 ENDDO
3433 ENDDO
3434 ENDDO
3435
3436 END IF
3437
3438 CALL vertical_diffusion_s( rt_tendf, config_flags, var_mix, mu, xkhv, &
3439 dn, dnw, rdz, rdzw, fnm, fnp, &
3440 .false., &
3441 ids, ide, jds, jde, kds, kde, &
3442 ims, ime, jms, jme, kms, kme, &
3443 its, ite, jts, jte, kts, kte )
3444
3445
3446 !*****************************************
3447 !*****************************************
3448 !MODIFICA al flusso di calore
3449 !
3450 !
3451 hflux: SELECT CASE( config_flags%isfflx )
3452 CASE (0) ! with fixed surface heat flux given in the namelist
3453 heat_flux = config_flags%tke_heat_flux ! constant heat flux value
3454 ! set in namelist.input
3455 DO j = j_start, j_end
3456 DO i = i_start, i_end
3457 rt_tendf(i,kts,j)=rt_tendf(i,kts,j) &
3458 +mu(i,j)*heat_flux*rdzw(i,kts,j)
3459 ENDDO
3460 ENDDO
3461
3462 CASE (1) ! use surface heat flux computed from surface routine
3463 DO j = j_start, j_end
3464 DO i = i_start, i_end
3465
3466 cpm = cp * (1. + 0.8 * moist(i,kts,j,P_QV))
3467 heat_flux = hfx(i,j)/cpm/rho(i,1,j)
3468 rt_tendf(i,kts,j)=rt_tendf(i,kts,j) &
3469 +mu(i,j)*heat_flux*rdzw(i,kts,j)
3470
3471 ENDDO
3472 ENDDO
3473
3474 CASE DEFAULT
3475 CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' )
3476 END SELECT hflux
3477
3478 !
3479 ! FINE MODIFICA al flusso di calore
3480 !*****************************************
3481 !*****************************************
3482
3483 If (km_opt .eq. 2) then
3484 CALL vertical_diffusion_s( tke_tendf(ims,kms,jms), &
3485 config_flags, tke(ims,kms,jms), &
3486 mu, xkhv, &
3487 dn, dnw, rdz, rdzw, fnm, fnp, &
3488 .true., &
3489 ids, ide, jds, jde, kds, kde, &
3490 ims, ime, jms, jme, kms, kme, &
3491 its, ite, jts, jte, kts, kte )
3492 endif
3493
3494 IF (n_moist .ge. PARAM_FIRST_SCALAR) THEN
3495
3496 moist_loop: do im = PARAM_FIRST_SCALAR, n_moist
3497
3498 IF ( (.not. config_flags%mix_full_fields) .and. (im == P_QV) ) THEN
3499
3500 DO j=jts,min(jte,jde-1)
3501 DO k=kts,kte-1
3502 DO i=its,min(ite,ide-1)
3503 var_mix(i,k,j) = moist(i,k,j,im) - qv_base(k)
3504 ENDDO
3505 ENDDO
3506 ENDDO
3507
3508 ELSE
3509
3510 DO j=jts,min(jte,jde-1)
3511 DO k=kts,kte-1
3512 DO i=its,min(ite,ide-1)
3513 var_mix(i,k,j) = moist(i,k,j,im)
3514 ENDDO
3515 ENDDO
3516 ENDDO
3517
3518 END IF
3519
3520
3521 CALL vertical_diffusion_s( moist_tendf(ims,kms,jms,im), &
3522 config_flags, var_mix, &
3523 mu, xkhv, &
3524 dn, dnw, rdz, rdzw, fnm, fnp, &
3525 .false., &
3526 ids, ide, jds, jde, kds, kde, &
3527 ims, ime, jms, jme, kms, kme, &
3528 its, ite, jts, jte, kts, kte )
3529
3530 !*****************************************
3531 !*****************************************
3532 !MODIFICATIONS for water vapor flux
3533 !
3534 !
3535 qflux: SELECT CASE( config_flags%isfflx )
3536 CASE (0)
3537 ! do nothing
3538 CASE (1) ! with surface moisture flux
3539 IF ( im == P_QV ) THEN
3540 DO j = j_start, j_end
3541 DO i = i_start, i_end
3542 moist_flux = qfx(i,j)
3543 moist_tendf(i,kts,j,im)=moist_tendf(i,kts,j,im) &
3544 +mu(i,j)*moist_flux*rdzw(i,kts,j)
3545 ENDDO
3546 ENDDO
3547 ENDIF
3548 CASE DEFAULT
3549 CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' )
3550 END SELECT qflux
3551 !
3552 ! END MODIFICATIONS for water vapor flux
3553 !*****************************************
3554 !*****************************************
3555 ENDDO moist_loop
3556
3557 ENDIF
3558
3559 IF (n_chem .ge. PARAM_FIRST_SCALAR) THEN
3560
3561 chem_loop: do im = PARAM_FIRST_SCALAR, n_chem
3562
3563 CALL vertical_diffusion_s( chem_tendf(ims,kms,jms,im), &
3564 config_flags, chem(ims,kms,jms,im), &
3565 mu, xkhv, &
3566 dn, dnw, rdz, rdzw, fnm, fnp, &
3567 .false., &
3568 ids, ide, jds, jde, kds, kde, &
3569 ims, ime, jms, jme, kms, kme, &
3570 its, ite, jts, jte, kts, kte )
3571 ENDDO chem_loop
3572
3573 ENDIF
3574
3575
3576 IF (n_scalar .ge. PARAM_FIRST_SCALAR) THEN
3577
3578 scalar_loop: do im = PARAM_FIRST_SCALAR, n_scalar
3579
3580 CALL vertical_diffusion_s( scalar_tendf(ims,kms,jms,im), &
3581 config_flags, scalar(ims,kms,jms,im), &
3582 mu, xkhv, &
3583 dn, dnw, rdz, rdzw, fnm, fnp, &
3584 .false., &
3585 ids, ide, jds, jde, kds, kde, &
3586 ims, ime, jms, jme, kms, kme, &
3587 its, ite, jts, jte, kts, kte )
3588 ENDDO scalar_loop
3589
3590 ENDIF
3591
3592 END SUBROUTINE vertical_diffusion_2
3593
3594 !=======================================================================
3595 !=======================================================================
3596
3597 SUBROUTINE vertical_diffusion_u_2( tendency, config_flags, mu, &
3598 defor13, xkmv, &
3599 dnw, rdzw, fnm, fnp, &
3600 ids, ide, jds, jde, kds, kde, &
3601 ims, ime, jms, jme, kms, kme, &
3602 its, ite, jts, jte, kts, kte )
3603
3604 !-----------------------------------------------------------------------
3605 ! Begin declarations.
3606
3607 IMPLICIT NONE
3608
3609 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3610
3611 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3612 ims, ime, jms, jme, kms, kme, &
3613 its, ite, jts, jte, kts, kte
3614
3615 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
3616 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
3617 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
3618 ! REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: zeta_z
3619
3620 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
3621
3622 REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
3623 INTENT(IN ) ::defor13, &
3624 xkmv, &
3625 rdzw
3626 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
3627
3628 ! LOCAL VARS
3629
3630 INTEGER :: i, j, k, ktf
3631
3632 INTEGER :: i_start, i_end, j_start, j_end
3633 INTEGER :: is_ext,ie_ext,js_ext,je_ext
3634
3635 REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau3
3636
3637 REAL , DIMENSION( its:ite, jts:jte) :: zzavg
3638
3639 REAL :: rdzu
3640
3641 ! End declarations.
3642 !-----------------------------------------------------------------------
3643
3644 ktf=MIN(kte,kde-1)
3645
3646 i_start = its
3647 i_end = ite
3648 j_start = jts
3649 j_end = MIN(jte,jde-1)
3650
3651 IF ( config_flags%open_xs .or. config_flags%specified .or. &
3652 config_flags%nested) i_start = MAX(ids+1,its)
3653 IF ( config_flags%open_xe .or. config_flags%specified .or. &
3654 config_flags%nested) i_end = MIN(ide-1,ite)
3655 IF ( config_flags%open_ys .or. config_flags%specified .or. &
3656 config_flags%nested) j_start = MAX(jds+1,jts)
3657 IF ( config_flags%open_ye .or. config_flags%specified .or. &
3658 config_flags%nested) j_end = MIN(jde-2,jte)
3659 IF ( config_flags%periodic_x ) i_start = its
3660 IF ( config_flags%periodic_x ) i_end = ite
3661
3662 ! titau3 = titau13
3663 is_ext=0
3664 ie_ext=0
3665 js_ext=0
3666 je_ext=0
3667 CALL cal_titau_13_31( config_flags, titau3, defor13, &
3668 mu, xkmv, fnm, fnp, &
3669 is_ext, ie_ext, js_ext, je_ext, &
3670 ids, ide, jds, jde, kds, kde, &
3671 ims, ime, jms, jme, kms, kme, &
3672 its, ite, jts, jte, kts, kte )
3673 !
3674 DO j = j_start, j_end
3675 DO k=kts+1,ktf
3676 DO i = i_start, i_end
3677
3678 rdzu = 2./(1./rdzw(i,k,j) + 1./rdzw(i-1,k,j))
3679 tendency(i,k,j)=tendency(i,k,j)-rdzu*(titau3(i,k+1,j)-titau3(i,k,j))
3680
3681 ENDDO
3682 ENDDO
3683 ENDDO
3684
3685 ! ******** MODIF...
3686 ! we will pick up the surface drag (titau3(i,kts,j)) later
3687 !
3688 DO j = j_start, j_end
3689 k=kts
3690 DO i = i_start, i_end
3691
3692 rdzu = 2./(1./rdzw(i,k,j) + 1./rdzw(i-1,k,j))
3693 tendency(i,k,j)=tendency(i,k,j)-rdzu*(titau3(i,k+1,j))
3694 ENDDO
3695 ENDDO
3696 ! ******** MODIF...
3697
3698 END SUBROUTINE vertical_diffusion_u_2
3699
3700 !=======================================================================
3701 !=======================================================================
3702
3703 SUBROUTINE vertical_diffusion_v_2( tendency, config_flags, mu, &
3704 defor23, xkmv, &
3705 dnw, rdzw, fnm, fnp, &
3706 ids, ide, jds, jde, kds, kde, &
3707 ims, ime, jms, jme, kms, kme, &
3708 its, ite, jts, jte, kts, kte )
3709 !-----------------------------------------------------------------------
3710 ! Begin declarations.
3711
3712 IMPLICIT NONE
3713
3714 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3715
3716 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3717 ims, ime, jms, jme, kms, kme, &
3718 its, ite, jts, jte, kts, kte
3719 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
3720 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
3721 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
3722 ! REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: zeta_z
3723
3724 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
3725
3726 REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
3727 INTENT(IN ) ::defor23, &
3728 xkmv, &
3729 rdzw
3730
3731 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
3732
3733 ! LOCAL VARS
3734
3735 INTEGER :: i, j, k, ktf
3736
3737 INTEGER :: i_start, i_end, j_start, j_end
3738 INTEGER :: is_ext,ie_ext,js_ext,je_ext
3739
3740 REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau3
3741
3742 REAL , DIMENSION( its:ite, jts:jte) :: zzavg
3743
3744 REAL :: rdzv
3745
3746 ! End declarations.
3747 !-----------------------------------------------------------------------
3748
3749 ktf=MIN(kte,kde-1)
3750
3751 i_start = its
3752 i_end = MIN(ite,ide-1)
3753 j_start = jts
3754 j_end = jte
3755
3756 IF ( config_flags%open_xs .or. config_flags%specified .or. &
3757 config_flags%nested) i_start = MAX(ids+1,its)
3758 IF ( config_flags%open_xe .or. config_flags%specified .or. &
3759 config_flags%nested) i_end = MIN(ide-2,ite)
3760 IF ( config_flags%open_ys .or. config_flags%specified .or. &
3761 config_flags%nested) j_start = MAX(jds+1,jts)
3762 IF ( config_flags%open_ye .or. config_flags%specified .or. &
3763 config_flags%nested) j_end = MIN(jde-1,jte)
3764 IF ( config_flags%periodic_x ) i_start = its
3765 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
3766
3767 ! titau3 = titau23
3768 is_ext=0
3769 ie_ext=0
3770 js_ext=0
3771 je_ext=0
3772 CALL cal_titau_23_32( config_flags, titau3, defor23, &
3773 mu, xkmv, fnm, fnp, &
3774 is_ext, ie_ext, js_ext, je_ext, &
3775 ids, ide, jds, jde, kds, kde, &
3776 ims, ime, jms, jme, kms, kme, &
3777 its, ite, jts, jte, kts, kte )
3778
3779 DO j = j_start, j_end
3780 DO k = kts+1,ktf
3781 DO i = i_start, i_end
3782
3783 rdzv = 2./(1./rdzw(i,k,j) + 1./rdzw(i,k,j-1))
3784 tendency(i,k,j)=tendency(i,k,j)-rdzv*(titau3(i,k+1,j)-titau3(i,k,j))
3785
3786 ENDDO
3787 ENDDO
3788 ENDDO
3789
3790 ! ******** MODIF...
3791 ! we will pick up the surface drag (titau3(i,kts,j)) later
3792 !
3793 DO j = j_start, j_end
3794 k=kts
3795 DO i = i_start, i_end
3796
3797 rdzv = 2./(1./rdzw(i,k,j) + 1./rdzw(i,k,j-1))
3798 tendency(i,k,j)=tendency(i,k,j)-rdzv*(titau3(i,k+1,j))
3799
3800 ENDDO
3801 ENDDO
3802 ! ******** MODIF...
3803
3804 END SUBROUTINE vertical_diffusion_v_2
3805
3806 !=======================================================================
3807 !=======================================================================
3808
3809 SUBROUTINE vertical_diffusion_w_2(tendency, config_flags, mu, &
3810 defor33, tke, div, xkmv, &
3811 dn, rdz, &
3812 ids, ide, jds, jde, kds, kde, &
3813 ims, ime, jms, jme, kms, kme, &
3814 its, ite, jts, jte, kts, kte )
3815
3816 !-----------------------------------------------------------------------
3817 ! Begin declarations.
3818
3819 IMPLICIT NONE
3820
3821 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3822
3823 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3824 ims, ime, jms, jme, kms, kme, &
3825 its, ite, jts, jte, kts, kte
3826
3827 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
3828
3829 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
3830
3831 REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
3832 INTENT(IN ) ::defor33, &
3833 tke, &
3834 div, &
3835 xkmv, &
3836 rdz
3837
3838 REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: mu
3839
3840 ! LOCAL VARS
3841
3842 INTEGER :: i, j, k, ktf
3843
3844 INTEGER :: i_start, i_end, j_start, j_end
3845 INTEGER :: is_ext,ie_ext,js_ext,je_ext
3846
3847 REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau3
3848
3849 ! End declarations.
3850 !-----------------------------------------------------------------------
3851
3852 ktf=MIN(kte,kde-1)
3853
3854 i_start = its
3855 i_end = MIN(ite,ide-1)
3856 j_start = jts
3857 j_end = MIN(jte,jde-1)
3858
3859 IF ( config_flags%open_xs .or. config_flags%specified .or. &
3860 config_flags%nested) i_start = MAX(ids+1,its)
3861 IF ( config_flags%open_xe .or. config_flags%specified .or. &
3862 config_flags%nested) i_end = MIN(ide-2,ite)
3863 IF ( config_flags%open_ys .or. config_flags%specified .or. &
3864 config_flags%nested) j_start = MAX(jds+1,jts)
3865 IF ( config_flags%open_ye .or. config_flags%specified .or. &
3866 config_flags%nested) j_end = MIN(jde-2,jte)
3867 IF ( config_flags%periodic_x ) i_start = its
3868 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
3869
3870 ! titau3 = titau33
3871 is_ext=0
3872 ie_ext=0
3873 js_ext=0
3874 je_ext=0
3875 CALL cal_titau_11_22_33( config_flags, titau3, &
3876 mu, tke, xkmv, defor33, &
3877 is_ext, ie_ext, js_ext, je_ext, &
3878 ids, ide, jds, jde, kds, kde, &
3879 ims, ime, jms, jme, kms, kme, &
3880 its, ite, jts, jte, kts, kte )
3881
3882 ! DO j = j_start, j_end
3883 ! DO k = kts+1, ktf
3884 ! DO i = i_start, i_end
3885 ! titau3(i,k,j)=titau3(i,k,j)*zeta_z(i,j)
3886 ! ENDDO
3887 ! ENDDO
3888 ! ENDDO
3889
3890 DO j = j_start, j_end
3891 DO k = kts+1, ktf
3892 DO i = i_start, i_end
3893 tendency(i,k,j)=tendency(i,k,j)-rdz(i,k,j)*(titau3(i,k,j)-titau3(i,k-1,j))
3894 ENDDO
3895 ENDDO
3896 ENDDO
3897
3898 END SUBROUTINE vertical_diffusion_w_2
3899
3900 !=======================================================================
3901 !=======================================================================
3902
3903 SUBROUTINE vertical_diffusion_s( tendency, config_flags, var, mu, xkhv, &
3904 dn, dnw, rdz, rdzw, fnm, fnp, &
3905 doing_tke, &
3906 ids, ide, jds, jde, kds, kde, &
3907 ims, ime, jms, jme, kms, kme, &
3908 its, ite, jts, jte, kts, kte )
3909
3910 !-----------------------------------------------------------------------
3911 ! Begin declarations.
3912
3913 IMPLICIT NONE
3914
3915 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3916
3917 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3918 ims, ime, jms, jme, kms, kme, &
3919 its, ite, jts, jte, kts, kte
3920
3921 LOGICAL, INTENT(IN ) :: doing_tke
3922
3923 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
3924 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
3925 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
3926 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
3927
3928 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
3929
3930 REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(IN) :: xkhv
3931
3932 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: mu
3933
3934 REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
3935 INTENT(IN ) :: var, &
3936 rdz, &
3937 rdzw
3938 ! LOCAL VARS
3939
3940 INTEGER :: i, j, k, ktf
3941
3942 INTEGER :: i_start, i_end, j_start, j_end
3943
3944 REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: H3, &
3945 xkxavg, &
3946 rravg
3947
3948 REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: tmptendf
3949
3950 ! End declarations.
3951 !-----------------------------------------------------------------------
3952
3953 ktf=MIN(kte,kde-1)
3954
3955 i_start = its
3956 i_end = MIN(ite,ide-1)
3957 j_start = jts
3958 j_end = MIN(jte,jde-1)
3959
3960 IF ( config_flags%open_xs .or. config_flags%specified .or. &
3961 config_flags%nested) i_start = MAX(ids+1,its)
3962 IF ( config_flags%open_xe .or. config_flags%specified .or. &
3963 config_flags%nested) i_end = MIN(ide-2,ite)
3964 IF ( config_flags%open_ys .or. config_flags%specified .or. &
3965 config_flags%nested) j_start = MAX(jds+1,jts)
3966 IF ( config_flags%open_ye .or. config_flags%specified .or. &
3967 config_flags%nested) j_end = MIN(jde-2,jte)
3968 IF ( config_flags%periodic_x ) i_start = its
3969 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
3970
3971 IF (doing_tke) THEN
3972 DO j = j_start, j_end
3973 DO k = kts,ktf
3974 DO i = i_start, i_end
3975 tmptendf(i,k,j)=tendency(i,k,j)
3976 ENDDO
3977 ENDDO
3978 ENDDO
3979 ENDIF
3980
3981 ! H3
3982
3983 xkxavg = 0.
3984
3985 DO j = j_start, j_end
3986 DO k = kts+1,ktf
3987 DO i = i_start, i_end
3988 xkxavg(i,k,j)=fnm(k)*xkhv(i,k,j)+fnp(k)*xkhv(i,k-1,j)
3989 H3(i,k,j)=-xkxavg(i,k,j)*(var(i,k,j)-var(i,k-1,j))*rdz(i,k,j)
3990 ! H3(i,k,j)=-xkxavg(i,k,j)*zeta_z(i,j)* &
3991 ! (var(i,k,j)-var(i,k-1,j))/dn(k)
3992 ENDDO
3993 ENDDO
3994 ENDDO
3995
3996 DO j = j_start, j_end
3997 DO i = i_start, i_end
3998 H3(i,kts,j)=0.
3999 H3(i,ktf+1,j)=0.
4000 ! H3(i,kts,j)=H3(i,kts+1,j)
4001 ! H3(i,ktf+1,j)=H3(i,ktf,j)
4002 ENDDO
4003 ENDDO
4004
4005 DO j = j_start, j_end
4006 DO k = kts,ktf
4007 DO i = i_start, i_end
4008 tendency(i,k,j)=tendency(i,k,j) &
4009 -mu(i,j)*(H3(i,k+1,j)-H3(i,k,j))*rdzw(i,k,j)
4010 ENDDO
4011 ENDDO
4012 ENDDO
4013
4014 IF (doing_tke) THEN
4015 DO j = j_start, j_end
4016 DO k = kts,ktf
4017 DO i = i_start, i_end
4018 tendency(i,k,j)=tmptendf(i,k,j)+2.* &
4019 (tendency(i,k,j)-tmptendf(i,k,j))
4020 ENDDO
4021 ENDDO
4022 ENDDO
4023 ENDIF
4024
4025 END SUBROUTINE vertical_diffusion_s
4026
4027 !=======================================================================
4028 !=======================================================================
4029
4030 SUBROUTINE cal_titau_11_22_33( config_flags, titau, &
4031 mu, tke, xkx, defor, &
4032 is_ext, ie_ext, js_ext, je_ext, &
4033 ids, ide, jds, jde, kds, kde, &
4034 ims, ime, jms, jme, kms, kme, &
4035 its, ite, jts, jte, kts, kte )
4036
4037 ! History: Sep 2003 Changes by George Bryan and Jason Knievel, NCAR
4038 ! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
4039 ! Aug 2000 Original code by Shu-Hua Chen, UC-Davis
4040
4041 ! Purpose: This routine calculates stress terms (taus) for use in
4042 ! the calculation of production of TKE by sheared wind
4043
4044 ! References: Klemp and Wilhelmson (JAS 1978)
4045 ! Deardorff (B-L Meteor 1980)
4046 ! Chen and Dudhia (NCAR WRF physics report 2000)
4047
4048 ! Key:
4049
4050 !-----------------------------------------------------------------------
4051 ! Begin declarations.
4052
4053 IMPLICIT NONE
4054
4055 TYPE( grid_config_rec_type ), INTENT( IN ) &
4056 :: config_flags
4057
4058 INTEGER, INTENT( IN ) &
4059 :: ids, ide, jds, jde, kds, kde, &
4060 ims, ime, jms, jme, kms, kme, &
4061 its, ite, jts, jte, kts, kte
4062
4063 INTEGER, INTENT( IN ) &
4064 :: is_ext, ie_ext, js_ext, je_ext
4065
4066 REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ), INTENT( INOUT ) &
4067 :: titau
4068
4069 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
4070 :: defor, xkx, tke
4071
4072 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
4073 :: mu
4074
4075 ! Local variables.
4076
4077 INTEGER &
4078 :: i, j, k, ktf, i_start, i_end, j_start, j_end
4079
4080 ! End declarations.
4081 !-----------------------------------------------------------------------
4082
4083 ktf = MIN( kte, kde-1 )
4084
4085 i_start = its
4086 i_end = ite
4087 j_start = jts
4088 j_end = jte
4089
4090 IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
4091 config_flags%nested) i_start = MAX( ids+1, its )
4092 IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
4093 config_flags%nested) i_end = MIN( ide-1, ite )
4094 IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
4095 config_flags%nested) j_start = MAX( jds+1, jts )
4096 IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
4097 config_flags%nested) j_end = MIN( jde-1, jte )
4098 IF ( config_flags%periodic_x ) i_start = its
4099 IF ( config_flags%periodic_x ) i_end = ite
4100
4101 i_start = i_start - is_ext
4102 i_end = i_end + ie_ext
4103 j_start = j_start - js_ext
4104 j_end = j_end + je_ext
4105
4106 IF ( config_flags%km_opt .EQ. 2) THEN
4107 DO j = j_start,j_end
4108 DO k = kts,ktf
4109 DO i = i_start,i_end
4110 titau(i,k,j) = mu(i,j) * ( - xkx(i,k,j) * ( defor(i,k,j) ) )
4111 END DO
4112 END DO
4113 END DO
4114 ELSE
4115 DO j = j_start, j_end
4116 DO k = kts, ktf
4117 DO i = i_start, i_end
4118 titau(i,k,j) = - mu(i,j) * xkx(i,k,j) * defor(i,k,j)
4119 END DO
4120 END DO
4121 END DO
4122 END IF
4123
4124 END SUBROUTINE cal_titau_11_22_33
4125
4126 !=======================================================================
4127 !=======================================================================
4128
4129 SUBROUTINE cal_titau_12_21( config_flags, titau, &
4130 mu, xkx, defor, &
4131 is_ext, ie_ext, js_ext, je_ext, &
4132 ids, ide, jds, jde, kds, kde, &
4133 ims, ime, jms, jme, kms, kme, &
4134 its, ite, jts, jte, kts, kte )
4135
4136 ! History: Sep 2003 Modifications by George Bryan and Jason Knievel, NCAR
4137 ! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
4138 ! Aug 2000 Original code by Shu-Hua Chen, UC-Davis
4139
4140 ! Pusrpose This routine calculates the stress terms (taus) for use in
4141 ! the calculation of production of TKE by sheared wind
4142
4143 ! References: Klemp and Wilhelmson (JAS 1978)
4144 ! Deardorff (B-L Meteor 1980)
4145 ! Chen and Dudhia (NCAR WRF physics report 2000)
4146
4147 ! Key:
4148
4149 !-----------------------------------------------------------------------
4150 ! Begin declarations.
4151
4152 IMPLICIT NONE
4153
4154 TYPE( grid_config_rec_type), INTENT( IN ) &
4155 :: config_flags
4156
4157 INTEGER, INTENT( IN ) &
4158 :: ids, ide, jds, jde, kds, kde, &
4159 ims, ime, jms, jme, kms, kme, &
4160 its, ite, jts, jte, kts, kte
4161
4162 INTEGER, INTENT( IN ) &
4163 :: is_ext, ie_ext, js_ext, je_ext
4164
4165 REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ), INTENT( INOUT ) &
4166 :: titau
4167
4168 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
4169 :: defor, xkx
4170
4171 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
4172 :: mu
4173
4174 ! Local variables.
4175
4176 INTEGER &
4177 :: i, j, k, ktf, i_start, i_end, j_start, j_end
4178
4179 REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ) &
4180 :: xkxavg
4181
4182 REAL, DIMENSION( its-1:ite+1, jts-1:jte+1 ) &
4183 :: muavg
4184
4185 ! End declarations.
4186 !-----------------------------------------------------------------------
4187
4188 ktf = MIN( kte, kde-1 )
4189
4190 ! Needs one more point in the x and y directions.
4191
4192 i_start = its
4193 i_end = ite
4194 j_start = jts
4195 j_end = jte
4196
4197 IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
4198 config_flags%nested ) i_start = MAX( ids+1, its )
4199 IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
4200 config_flags%nested ) i_end = MIN( ide-1, ite )
4201 IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
4202 config_flags%nested ) j_start = MAX( jds+1, jts )
4203 IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
4204 config_flags%nested ) j_end = MIN( jde-1, jte )
4205 IF ( config_flags%periodic_x ) i_start = its
4206 IF ( config_flags%periodic_x ) i_end = ite
4207
4208 i_start = i_start - is_ext
4209 i_end = i_end + ie_ext
4210 j_start = j_start - js_ext
4211 j_end = j_end + je_ext
4212
4213 DO j = j_start, j_end
4214 DO k = kts, ktf
4215 DO i = i_start, i_end
4216 xkxavg(i,k,j) = 0.25 * ( xkx(i-1,k,j ) + xkx(i,k,j ) + &
4217 xkx(i-1,k,j-1) + xkx(i,k,j-1) )
4218 END DO
4219 END DO
4220 END DO
4221
4222 DO j = j_start, j_end
4223 DO i = i_start, i_end
4224 muavg(i,j) = 0.25 * ( mu(i-1,j ) + mu(i,j ) + &
4225 mu(i-1,j-1) + mu(i,j-1) )
4226 END DO
4227 END DO
4228
4229 ! titau12 or titau21
4230
4231 DO j = j_start, j_end
4232 DO k = kts, ktf
4233 DO i = i_start, i_end
4234 titau(i,k,j) = - muavg(i,j) * xkxavg(i,k,j) * defor(i,k,j)
4235 END DO
4236 END DO
4237 END DO
4238
4239 END SUBROUTINE cal_titau_12_21
4240
4241 !=======================================================================
4242
4243 SUBROUTINE cal_titau_13_31( config_flags, titau, &
4244 defor, mu, xkx, fnm, fnp, &
4245 is_ext, ie_ext, js_ext, je_ext, &
4246 ids, ide, jds, jde, kds, kde, &
4247 ims, ime, jms, jme, kms, kme, &
4248 its, ite, jts, jte, kts, kte )
4249
4250 ! History: Sep 2003 Modifications by George Bryan and Jason Knievel, NCAR
4251 ! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
4252 ! Aug 2000 Original code by Shu-Hua Chen, UC-Davis
4253
4254 ! Purpose: This routine calculates the stress terms (taus) for use in
4255 ! the calculation of production of TKE by sheared wind
4256
4257 ! References: Klemp and Wilhelmson (JAS 1978)
4258 ! Deardorff (B-L Meteor 1980)
4259 ! Chen and Dudhia (NCAR WRF physics report 2000)
4260
4261 ! Key:
4262
4263 !-----------------------------------------------------------------------
4264 ! Begin declarations.
4265
4266 IMPLICIT NONE
4267
4268 TYPE( grid_config_rec_type), INTENT( IN ) &
4269 :: config_flags
4270
4271 INTEGER, INTENT( IN ) &
4272 :: ids, ide, jds, jde, kds, kde, &
4273 ims, ime, jms, jme, kms, kme, &
4274 its, ite, jts, jte, kts, kte
4275
4276 INTEGER, INTENT( IN ) &
4277 :: is_ext, ie_ext, js_ext, je_ext
4278
4279 REAL, DIMENSION( kms:kme ), INTENT( IN ) &
4280 :: fnm, fnp
4281
4282 REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ), INTENT( INOUT ) &
4283 :: titau
4284
4285 REAL, DIMENSION( ims:ime, kms:kme, jms:jme), INTENT( IN ) &
4286 :: defor, xkx
4287
4288 REAL, DIMENSION( ims:ime, jms:jme), INTENT( IN ) &
4289 :: mu
4290
4291 ! Local variables.
4292
4293 INTEGER &
4294 :: i, j, k, ktf, i_start, i_end, j_start, j_end
4295
4296 REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ) &
4297 :: xkxavg
4298
4299 REAL, DIMENSION( its-1:ite+1, jts-1:jte+1 ) &
4300 :: muavg
4301
4302 ! End declarations.
4303 !-----------------------------------------------------------------------
4304
4305 ktf = MIN( kte, kde-1 )
4306
4307 ! Find ide-1 and jde-1 for averaging to p point.
4308
4309 i_start = its
4310 i_end = ite
4311 j_start = jts
4312 j_end = MIN( jte, jde-1 )
4313
4314 IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
4315 config_flags%nested) i_start = MAX( ids+1, its )
4316 IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
4317 config_flags%nested) i_end = MIN( ide-1, ite )
4318 IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
4319 config_flags%nested) j_start = MAX( jds+1, jts )
4320 IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
4321 config_flags%nested) j_end = MIN( jde-2, jte )
4322 IF ( config_flags%periodic_x ) i_start = its
4323 IF ( config_flags%periodic_x ) i_end = ite
4324
4325 i_start = i_start - is_ext
4326 i_end = i_end + ie_ext
4327 j_start = j_start - js_ext
4328 j_end = j_end + je_ext
4329
4330 DO j = j_start, j_end
4331 DO k = kts+1, ktf
4332 DO i = i_start, i_end
4333 xkxavg(i,k,j) = 0.5 * ( fnm(k) * ( xkx(i,k ,j) + xkx(i-1,k ,j) ) + &
4334 fnp(k) * ( xkx(i,k-1,j) + xkx(i-1,k-1,j) ) )
4335 END DO
4336 END DO
4337 END DO
4338
4339 DO j = j_start, j_end
4340 DO i = i_start, i_end
4341 muavg(i,j) = 0.5 * ( mu(i,j) + mu(i-1,j) )
4342 END DO
4343 END DO
4344
4345 DO j = j_start, j_end
4346 DO k = kts+1, ktf
4347 DO i = i_start, i_end
4348 titau(i,k,j) = - muavg(i,j) * xkxavg(i,k,j) * defor(i,k,j)
4349 ENDDO
4350 ENDDO
4351 ENDDO
4352
4353 DO j = j_start, j_end
4354 DO i = i_start, i_end
4355 titau(i,kts ,j) = 0.0
4356 titau(i,ktf+1,j) = 0.0
4357 ENDDO
4358 ENDDO
4359
4360 END SUBROUTINE cal_titau_13_31
4361
4362 !=======================================================================
4363 !=======================================================================
4364
4365 SUBROUTINE cal_titau_23_32( config_flags, titau, defor, &
4366 mu, xkx, fnm, fnp, &
4367 is_ext, ie_ext, js_ext, je_ext, &
4368 ids, ide, jds, jde, kds, kde, &
4369 ims, ime, jms, jme, kms, kme, &
4370 its, ite, jts, jte, kts, kte )
4371
4372 ! History: Sep 2003 Changes by George Bryan and Jason Knievel, NCAR
4373 ! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
4374 ! Aug 2000 Original code by Shu-Hua Chen, UC-Davis
4375
4376 ! Purpose: This routine calculates stress terms (taus) for use in
4377 ! the calculation of production of TKE by sheared wind
4378
4379 ! References: Klemp and Wilhelmson (JAS 1978)
4380 ! Deardorff (B-L Meteor 1980)
4381 ! Chen and Dudhia (NCAR WRF physics report 2000)
4382
4383 ! Key:
4384
4385 !-----------------------------------------------------------------------
4386 ! Begin declarations.
4387
4388 IMPLICIT NONE
4389
4390 TYPE( grid_config_rec_type ), INTENT( IN ) &
4391 :: config_flags
4392
4393 INTEGER, INTENT( IN ) &
4394 :: ids, ide, jds, jde, kds, kde, &
4395 ims, ime, jms, jme, kms, kme, &
4396 its, ite, jts, jte, kts, kte
4397
4398 INTEGER, INTENT( IN ) &
4399 :: is_ext,ie_ext,js_ext,je_ext
4400
4401 REAL, DIMENSION( kms:kme ), INTENT( IN ) &
4402 :: fnm, fnp
4403
4404 REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ), INTENT( INOUT ) &
4405 :: titau
4406
4407 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
4408 :: defor, xkx
4409
4410 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
4411 :: mu
4412
4413 ! Local variables.
4414
4415 INTEGER &
4416 :: i, j, k, ktf, i_start, i_end, j_start, j_end
4417
4418 REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ) &
4419 :: xkxavg
4420
4421 REAL, DIMENSION( its-1:ite+1, jts-1:jte+1 ) &
4422 :: muavg
4423
4424 ! End declarations.
4425 !-----------------------------------------------------------------------
4426
4427 ktf = MIN( kte, kde-1 )
4428
4429 ! Find ide-1 and jde-1 for averaging to p point.
4430
4431 i_start = its
4432 i_end = MIN( ite, ide-1 )
4433 j_start = jts
4434 j_end = jte
4435
4436 IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
4437 config_flags%nested) i_start = MAX( ids+1, its )
4438 IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
4439 config_flags%nested) i_end = MIN( ide-2, ite )
4440 IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
4441 config_flags%nested) j_start = MAX( jds+1, jts )
4442 IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
4443 config_flags%nested) j_end = MIN( jde-1, jte )
4444 IF ( config_flags%periodic_x ) i_start = its
4445 IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
4446
4447 i_start = i_start - is_ext
4448 i_end = i_end + ie_ext
4449 j_start = j_start - js_ext
4450 j_end = j_end + je_ext
4451
4452 DO j = j_start, j_end
4453 DO k = kts+1, ktf
4454 DO i = i_start, i_end
4455 xkxavg(i,k,j) = 0.5 * ( fnm(k) * ( xkx(i,k ,j) + xkx(i,k ,j-1) ) + &
4456 fnp(k) * ( xkx(i,k-1,j) + xkx(i,k-1,j-1) ) )
4457 END DO
4458 END DO
4459 END DO
4460
4461 DO j = j_start, j_end
4462 DO i = i_start, i_end
4463 muavg(i,j) = 0.5 * ( mu(i,j) + mu(i,j-1) )
4464 END DO
4465 END DO
4466
4467 DO j = j_start, j_end
4468 DO k = kts+1, ktf
4469 DO i = i_start, i_end
4470 titau(i,k,j) = - muavg(i,j) * xkxavg(i,k,j) * defor(i,k,j)
4471 END DO
4472 END DO
4473 END DO
4474
4475 DO j = j_start, j_end
4476 DO i = i_start, i_end
4477 titau(i,kts ,j) = 0.0
4478 titau(i,ktf+1,j) = 0.0
4479 END DO
4480 END DO
4481
4482 END SUBROUTINE cal_titau_23_32
4483
4484 !=======================================================================
4485 !=======================================================================
4486
4487 SUBROUTINE phy_bc ( config_flags,div,defor11,defor22,defor33, &
4488 defor12,defor13,defor23,xkmh,xkmv,xkhh,xkhv,tke, &
4489 RUBLTEN, RVBLTEN, &
4490 ids, ide, jds, jde, kds, kde, &
4491 ims, ime, jms, jme, kms, kme, &
4492 ips, ipe, jps, jpe, kps, kpe, &
4493 its, ite, jts, jte, kts, kte )
4494
4495 !------------------------------------------------------------------------------
4496 ! Begin declarations.
4497
4498 IMPLICIT NONE
4499
4500 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
4501
4502 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4503 ims, ime, jms, jme, kms, kme, &
4504 ips, ipe, jps, jpe, kps, kpe, &
4505 its, ite, jts, jte, kts, kte
4506
4507 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::RUBLTEN, &
4508 RVBLTEN, &
4509 defor11, &
4510 defor22, &
4511 defor33, &
4512 defor12, &
4513 defor13, &
4514 defor23, &
4515 xkmh, &
4516 xkmv, &
4517 xkhh, &
4518 xkhv, &
4519 tke, &
4520 div
4521
4522 ! End declarations.
4523 !-----------------------------------------------------------------------
4524
4525 IF(config_flags%bl_pbl_physics .GT. 0) THEN
4526
4527 CALL set_physical_bc3d( RUBLTEN , 't', config_flags, &
4528 ids, ide, jds, jde, kds, kde, &
4529 ims, ime, jms, jme, kms, kme, &
4530 ips, ipe, jps, jpe, kps, kpe, &
4531 its, ite, jts, jte, kts, kte )
4532
4533 CALL set_physical_bc3d( RVBLTEN , 't', config_flags, &
4534 ids, ide, jds, jde, kds, kde, &
4535 ims, ime, jms, jme, kms, kme, &
4536 ips, ipe, jps, jpe, kps, kpe, &
4537 its, ite, jts, jte, kts, kte )
4538
4539 ENDIF
4540
4541 ! move out of the conditional, below; this one is needed for
4542 ! all diff_opt cases. JM
4543 CALL set_physical_bc3d( xkmh , 't', config_flags, &
4544 ids, ide, jds, jde, kds, kde, &
4545 ims, ime, jms, jme, kms, kme, &
4546 ips, ipe, jps, jpe, kps, kpe, &
4547 its, ite, jts, jte, kts, kte )
4548
4549 IF(config_flags%diff_opt .eq. 2) THEN
4550
4551 CALL set_physical_bc3d( xkmv , 't', config_flags, &
4552 ids, ide, jds, jde, kds, kde, &
4553 ims, ime, jms, jme, kms, kme, &
4554 ips, ipe, jps, jpe, kps, kpe, &
4555 its, ite, jts, jte, kts, kte )
4556
4557 CALL set_physical_bc3d( xkhh , 't', config_flags, &
4558 ids, ide, jds, jde, kds, kde, &
4559 ims, ime, jms, jme, kms, kme, &
4560 ips, ipe, jps, jpe, kps, kpe, &
4561 its, ite, jts, jte, kts, kte )
4562
4563 CALL set_physical_bc3d( xkhv , 't', config_flags, &
4564 ids, ide, jds, jde, kds, kde, &
4565 ims, ime, jms, jme, kms, kme, &
4566 ips, ipe, jps, jpe, kps, kpe, &
4567 its, ite, jts, jte, kts, kte )
4568
4569 CALL set_physical_bc3d( tke , 't', config_flags, &
4570 ids, ide, jds, jde, kds, kde, &
4571 ims, ime, jms, jme, kms, kme, &
4572 ips, ipe, jps, jpe, kps, kpe, &
4573 its, ite, jts, jte, kts, kte )
4574
4575 CALL set_physical_bc3d( div , 't', config_flags, &
4576 ids, ide, jds, jde, kds, kde, &
4577 ims, ime, jms, jme, kms, kme, &
4578 ips, ipe, jps, jpe, kps, kpe, &
4579 its, ite, jts, jte, kts, kte )
4580
4581 CALL set_physical_bc3d( defor11 , 't', config_flags, &
4582 ids, ide, jds, jde, kds, kde, &
4583 ims, ime, jms, jme, kms, kme, &
4584 ips, ipe, jps, jpe, kps, kpe, &
4585 its, ite, jts, jte, kts, kte )
4586
4587 CALL set_physical_bc3d( defor22 , 't', config_flags, &
4588 ids, ide, jds, jde, kds, kde, &
4589 ims, ime, jms, jme, kms, kme, &
4590 ips, ipe, jps, jpe, kps, kpe, &
4591 its, ite, jts, jte, kts, kte )
4592
4593 CALL set_physical_bc3d( defor33 , 't', config_flags, &
4594 ids, ide, jds, jde, kds, kde, &
4595 ims, ime, jms, jme, kms, kme, &
4596 ips, ipe, jps, jpe, kps, kpe, &
4597 its, ite, jts, jte, kts, kte )
4598
4599 CALL set_physical_bc3d( defor12 , 'd', config_flags, &
4600 ids, ide, jds, jde, kds, kde, &
4601 ims, ime, jms, jme, kms, kme, &
4602 ips, ipe, jps, jpe, kps, kpe, &
4603 its, ite, jts, jte, kts, kte )
4604
4605 CALL set_physical_bc3d( defor13 , 'e', config_flags, &
4606 ids, ide, jds, jde, kds, kde, &
4607 ims, ime, jms, jme, kms, kme, &
4608 ips, ipe, jps, jpe, kps, kpe, &
4609 its, ite, jts, jte, kts, kte )
4610
4611 CALL set_physical_bc3d( defor23 , 'f', config_flags, &
4612 ids, ide, jds, jde, kds, kde, &
4613 ims, ime, jms, jme, kms, kme, &
4614 ips, ipe, jps, jpe, kps, kpe, &
4615 its, ite, jts, jte, kts, kte )
4616
4617 ENDIF
4618
4619 END SUBROUTINE phy_bc
4620
4621 !=======================================================================
4622 !=======================================================================
4623
4624 SUBROUTINE tke_rhs( tendency, BN2, config_flags, &
4625 defor11, defor22, defor33, &
4626 defor12, defor13, defor23, &
4627 u, v, w, div, tke, mu, &
4628 theta, p, p8w, t8w, z, fnm, fnp, &
4629 cf1, cf2, cf3, msftx, msfty, &
4630 xkmh, xkmv, xkhv, &
4631 rdx, rdy, dx, dy, dt, zx, zy, &
4632 rdz, rdzw, dn, dnw, cr_len, &
4633 hfx, qfx, qv, ust, rho, &
4634 ids, ide, jds, jde, kds, kde, &
4635 ims, ime, jms, jme, kms, kme, &
4636 its, ite, jts, jte, kts, kte )
4637
4638 !-----------------------------------------------------------------------
4639 ! Begin declarations.
4640
4641 IMPLICIT NONE
4642
4643 TYPE( grid_config_rec_type ), INTENT( IN ) &
4644 :: config_flags
4645
4646 INTEGER, INTENT( IN ) &
4647 :: ids, ide, jds, jde, kds, kde, &
4648 ims, ime, jms, jme, kms, kme, &
4649 its, ite, jts, jte, kts, kte
4650
4651 REAL, INTENT( IN ) &
4652 :: cf1, cf2, cf3, dt, rdx, rdy, dx, dy, cr_len
4653
4654 REAL, DIMENSION( kms:kme ), INTENT( IN ) &
4655 :: fnm, fnp, dnw, dn
4656
4657 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
4658 :: msftx, msfty
4659
4660 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
4661 :: tendency
4662
4663 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
4664 :: defor11, defor22, defor33, defor12, defor13, defor23, &
4665 div, BN2, tke, xkmh, xkmv, xkhv, zx, zy, u, v, w, theta, &
4666 p, p8w, t8w, z, rdz, rdzw
4667
4668 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
4669 :: mu
4670
4671 REAL, DIMENSION ( ims:ime, jms:jme ), INTENT( IN ) &
4672 :: hfx, ust, qfx
4673 REAL, DIMENSION ( ims:ime, kms:kme, jms:jme ), INTENT ( IN ) &
4674 :: qv, rho
4675
4676 ! Local variables.
4677
4678 INTEGER &
4679 :: i, j, k, ktf, i_start, i_end, j_start, j_end
4680
4681 ! End declarations.
4682 !-----------------------------------------------------------------------
4683
4684 CALL tke_shear( tendency, config_flags, &
4685 defor11, defor22, defor33, &
4686 defor12, defor13, defor23, &
4687 u, v, w, tke, ust, mu, fnm, fnp, &
4688 cf1, cf2, cf3, msftx, msfty, &
4689 xkmh, xkmv, &
4690 rdx, rdy, zx, zy, rdz, rdzw, dnw, dn, &
4691 ids, ide, jds, jde, kds, kde, &
4692 ims, ime, jms, jme, kms, kme, &
4693 its, ite, jts, jte, kts, kte )
4694
4695 CALL tke_buoyancy( tendency, config_flags, mu, &
4696 tke, xkhv, BN2, theta, dt, &
4697 hfx, qfx, qv, rho, &
4698 ids, ide, jds, jde, kds, kde, &
4699 ims, ime, jms, jme, kms, kme, &
4700 its, ite, jts, jte, kts, kte )
4701
4702 CALL tke_dissip( tendency, config_flags, &
4703 mu, tke, bn2, theta, p8w, t8w, z, &
4704 dx, dy,rdz, rdzw, cr_len, &
4705 ids, ide, jds, jde, kds, kde, &
4706 ims, ime, jms, jme, kms, kme, &
4707 its, ite, jts, jte, kts, kte )
4708
4709 ! Set a lower limit on TKE.
4710
4711 ktf = MIN( kte, kde-1 )
4712 i_start = its
4713 i_end = MIN( ite, ide-1 )
4714 j_start = jts
4715 j_end = MIN( jte, jde-1 )
4716
4717 IF ( config_flags%open_xs .or. config_flags%specified .or. &
4718 config_flags%nested) i_start = MAX(ids+1,its)
4719 IF ( config_flags%open_xe .or. config_flags%specified .or. &
4720 config_flags%nested) i_end = MIN(ide-2,ite)
4721 IF ( config_flags%open_ys .or. config_flags%specified .or. &
4722 config_flags%nested) j_start = MAX(jds+1,jts)
4723 IF ( config_flags%open_ye .or. config_flags%specified .or. &
4724 config_flags%nested) j_end = MIN(jde-2,jte)
4725 IF ( config_flags%periodic_x ) i_start = its
4726 IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
4727
4728 DO j = j_start, j_end
4729 DO k = kts, ktf
4730 DO i = i_start, i_end
4731 tendency(i,k,j) = max( tendency(i,k,j), -mu(i,j) * max( 0.0 , tke(i,k,j) ) / dt )
4732 END DO
4733 END DO
4734 END DO
4735
4736 END SUBROUTINE tke_rhs
4737
4738 !=======================================================================
4739 !=======================================================================
4740
4741 SUBROUTINE tke_buoyancy( tendency, config_flags, mu, &
4742 tke, xkhv, BN2, theta, dt, &
4743 hfx, qfx, qv, rho, &
4744 ids, ide, jds, jde, kds, kde, &
4745 ims, ime, jms, jme, kms, kme, &
4746 its, ite, jts, jte, kts, kte )
4747
4748 !-----------------------------------------------------------------------
4749 ! Begin declarations.
4750
4751 IMPLICIT NONE
4752
4753 TYPE( grid_config_rec_type ), INTENT( IN ) &
4754 :: config_flags
4755
4756 INTEGER, INTENT( IN ) &
4757 :: ids, ide, jds, jde, kds, kde, &
4758 ims, ime, jms, jme, kms, kme, &
4759 its, ite, jts, jte, kts, kte
4760
4761 REAL, INTENT( IN ) &
4762 :: dt
4763
4764 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
4765 :: tendency
4766
4767 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
4768 :: xkhv, tke, BN2, theta
4769
4770 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
4771 :: mu
4772
4773 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT ( IN ) &
4774 :: qv, rho
4775
4776 REAL, DIMENSION(ims:ime, jms:jme ), INTENT ( IN ) :: hfx, qfx
4777
4778 ! Local variables.
4779
4780 INTEGER &
4781 :: i, j, k, ktf
4782
4783 INTEGER &
4784 :: i_start, i_end, j_start, j_end
4785
4786 REAL :: heat_flux, heat_flux0
4787
4788 REAL :: cpm
4789
4790 ! End declarations.
4791 !-----------------------------------------------------------------------
4792
4793 !-----------------------------------------------------------------------
4794 ! Add to the TKE tendency the term that accounts for production of TKE
4795 ! due to buoyant motions.
4796
4797 ktf = MIN( kte, kde-1 )
4798 i_start = its
4799 i_end = MIN( ite, ide-1 )
4800 j_start = jts
4801 j_end = MIN( jte, jde-1 )
4802
4803 IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
4804 config_flags%nested ) i_start = MAX( ids+1, its )
4805 IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
4806 config_flags%nested ) i_end = MIN( ide-2, ite )
4807 IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
4808 config_flags%nested ) j_start = MAX( jds+1, jts )
4809 IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
4810 config_flags%nested ) j_end = MIN( jde-2, jte )
4811 IF ( config_flags%periodic_x ) i_start = its
4812 IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
4813
4814 DO j = j_start, j_end
4815 DO k = kts+1, ktf
4816 DO i = i_start, i_end
4817 tendency(i,k,j) = tendency(i,k,j) - mu(i,j) * xkhv(i,k,j) * BN2(i,k,j)
4818 END DO
4819 END DO
4820 END DO
4821
4822 ! MARTA: change in the computation of the tke's tendency at the surface.
4823 ! the buoyancy flux is the average of the surface heat flux (0.06) and the
4824 ! flux at the first w level
4825 !
4826 ! WCS 040331
4827
4828 hflux: SELECT CASE( config_flags%isfflx )
4829 CASE (0) ! with fixed surface heat flux given in the namelist
4830 heat_flux0 = config_flags%tke_heat_flux ! constant heat flux value
4831 ! set in namelist.input
4832 ! LES mods
4833 K=KTS
4834 DO j = j_start, j_end
4835 DO i = i_start, i_end
4836 heat_flux = heat_flux0
4837 tendency(i,k,j)= tendency(i,k,j) - &
4838 mu(i,j)*((xkhv(i,k,j)*BN2(i,k,j))- (g/theta(i,k,j))*heat_flux)/2.
4839
4840 ENDDO
4841 ENDDO
4842 CASE (1) ! use surface heat flux computed from surface routine
4843 K=KTS
4844 DO j = j_start, j_end
4845 DO i = i_start, i_end
4846 cpm = cp * (1. + 0.8*qv(i,k,j))
4847 heat_flux = (hfx(i,j)/cpm)/rho(i,k,j)
4848
4849 tendency(i,k,j)= tendency(i,k,j) - &
4850 mu(i,j)*((xkhv(i,k,j)*BN2(i,k,j))- (g/theta(i,k,j))*heat_flux)/2.
4851
4852 ENDDO
4853 ENDDO
4854 CASE DEFAULT
4855 CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' )
4856 END SELECT hflux
4857
4858 ! end of MARTA/WCS change
4859
4860 ! The tendency array now includes production of TKE from buoyant
4861 ! motions.
4862 !-----------------------------------------------------------------------
4863
4864 END SUBROUTINE tke_buoyancy
4865
4866 !=======================================================================
4867 !=======================================================================
4868
4869 SUBROUTINE tke_dissip( tendency, config_flags, &
4870 mu, tke, bn2, theta, p8w, t8w, z, &
4871 dx, dy, rdz, rdzw, cr_len_in, &
4872 ids, ide, jds, jde, kds, kde, &
4873 ims, ime, jms, jme, kms, kme, &
4874 its, ite, jts, jte, kts, kte )
4875
4876 ! History: Sep 2003 Changes by George Bryan and Jason Knievel, NCAR
4877 ! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
4878 ! Aug 2000 Original code by Shu-Hua Chen, UC-Davis
4879
4880 ! Purpose: This routine calculates dissipation of turbulent kinetic
4881 ! energy.
4882
4883 ! References: Klemp and Wilhelmson (JAS 1978)
4884 ! Deardorff (B-L Meteor 1980)
4885 ! Chen and Dudhia (NCAR WRF physics report 2000)
4886
4887 !-----------------------------------------------------------------------
4888 ! Begin declarations.
4889
4890 IMPLICIT NONE
4891
4892 TYPE( grid_config_rec_type ), INTENT( IN ) &
4893 :: config_flags
4894
4895 INTEGER, INTENT( IN ) &
4896 :: ids, ide, jds, jde, kds, kde, &
4897 ims, ime, jms, jme, kms, kme, &
4898 its, ite, jts, jte, kts, kte
4899
4900 REAL, INTENT( IN ) &
4901 :: dx, dy, cr_len_in
4902
4903 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
4904 :: tendency
4905
4906 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
4907 :: tke, bn2, theta, p8w, t8w, z, rdz, rdzw
4908
4909 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
4910 :: mu
4911
4912 ! Local variables.
4913
4914 REAL, DIMENSION( its:ite, kts:kte, jts:jte ) &
4915 :: dthrdn
4916
4917 REAL, DIMENSION( its:ite, kts:kte, jts:jte ) &
4918 :: l_scale
4919
4920 REAL, DIMENSION( its:ite ) &
4921 :: sumtke, sumtkez
4922
4923 INTEGER &
4924 :: i, j, k, ktf, i_start, i_end, j_start, j_end
4925
4926 REAL &
4927 :: disp_len, deltas, coefc, tmpdz, len_s, thetasfc, &
4928 thetatop, len_0, tketmp, tmp, cr_len, ce1, ce2
4929
4930 ! End declarations.
4931 !-----------------------------------------------------------------------
4932
4933 ce1 = ( c_k / 0.10 ) * 0.19
4934 ce2 = max( 0.0 , 0.93 - ce1 )
4935
4936 ktf = MIN( kte, kde-1 )
4937 i_start = its
4938 i_end = MIN(ite,ide-1)
4939 j_start = jts
4940 j_end = MIN(jte,jde-1)
4941
4942 IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
4943 config_flags%nested) i_start = MAX( ids+1, its )
4944 IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
4945 config_flags%nested) i_end = MIN( ide-2, ite )
4946 IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
4947 config_flags%nested) j_start = MAX( jds+1, jts )
4948 IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
4949 config_flags%nested) j_end = MIN( jde-2, jte )
4950 IF ( config_flags%periodic_x ) i_start = its
4951 IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
4952
4953 cr_len = cr_len_in
4954 cr_len = dx + 1.0 ! hardwire for mixing length = (dx*dy*dz)**(1/3).
4955 ! remove this for the alternate formulation
4956
4957 IF (dx .gt. cr_len) THEN
4958
4959 DO j = j_start, j_end
4960 DO i = i_start, i_end
4961 sumtke(i) = 0.0
4962 sumtkez(i) = 0.0
4963 END DO
4964 DO k = kts, ktf
4965 DO i = i_start, i_end
4966 tketmp = MAX( tke(i,k,j), 0.0 )
4967 sumtke(i) = sumtke(i) + SQRT(tketmp) / rdzw(i,k,j)
4968 sumtkez(i) = sumtkez(i)+ sumtke(i) * z(i,k,j)
4969 IF ( ABS( sumtke(i) ) .gt. 0.01 ) THEN
4970 len_0 = 0.2 * sumtkez(i) / sumtke(i)
4971 ELSE
4972 len_0 = 80.0
4973 ENDIF
4974 len_0 = MIN( 80.0, len_0)
4975 l_scale(i,k,j) = KARMAN * z(i,k,j) / &
4976 ( 1.0 + KARMAN * z(i,k,j) / len_0 )
4977 tendency(i,k,j) = tendency(i,k,j) - &
4978 mu(i,j) * 2.0 * SQRT( 2.0 ) / 15.0 * &
4979 tketmp**1.5 / l_scale(i,k,j)
4980 END DO
4981 END DO
4982 END DO
4983 ELSE
4984 CALL calc_l_scale( config_flags, tke, BN2, l_scale, &
4985 i_start, i_end, ktf, j_start, j_end, &
4986 dx, dy, rdzw, &
4987 ids, ide, jds, jde, kds, kde, &
4988 ims, ime, jms, jme, kms, kme, &
4989 its, ite, jts, jte, kts, kte )
4990 DO j = j_start, j_end
4991 DO k = kts, ktf
4992 DO i = i_start, i_end
4993 deltas = ( dx * dy / rdzw(i,k,j) )**0.33333333
4994 tketmp = MAX( tke(i,k,j), 1.0e-6 )
4995
4996 ! Apply Deardorff's (1980) "wall effect" at the bottom of the domain.
4997
4998 IF ( k .eq. kts .or. k .eq. ktf ) then
4999 coefc = 3.9
5000 ELSE
5001 coefc = ce1 + ce2 * l_scale(i,k,j) / deltas
5002 END IF
5003
5004 tendency(i,k,j) = tendency(i,k,j) - &
5005 mu(i,j) * coefc * tketmp**1.5 / l_scale(i,k,j)
5006 END DO
5007 END DO
5008 END DO
5009 ENDIF
5010
5011 END SUBROUTINE tke_dissip
5012
5013 !=======================================================================
5014 !=======================================================================
5015
5016 SUBROUTINE tke_shear( tendency, config_flags, &
5017 defor11, defor22, defor33, &
5018 defor12, defor13, defor23, &
5019 u, v, w, tke, ust, mu, fnm, fnp, &
5020 cf1, cf2, cf3, msftx, msfty, &
5021 xkmh, xkmv, &
5022 rdx, rdy, zx, zy, rdz, rdzw, dn, dnw, &
5023 ids, ide, jds, jde, kds, kde, &
5024 ims, ime, jms, jme, kms, kme, &
5025 its, ite, jts, jte, kts, kte )
5026
5027 ! History: Sep 2003 Rewritten by George Bryan and Jason Knievel,
5028 ! NCAR
5029 ! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
5030 ! Aug 2000 Original code by Shu-Hua Chen, UC-Davis
5031
5032 ! Purpose: This routine calculates the production of turbulent
5033 ! kinetic energy by stresses due to sheared wind.
5034
5035 ! References: Klemp and Wilhelmson (JAS 1978)
5036 ! Deardorff (B-L Meteor 1980)
5037 ! Chen and Dudhia (NCAR WRF physics report 2000)
5038
5039 ! Key:
5040
5041 ! avg temporary working array
5042 ! cf1
5043 ! cf2
5044 ! cf3
5045 ! defor11 deformation term ( du/dx + du/dx )
5046 ! defor12 deformation term ( dv/dx + du/dy ); same as defor21
5047 ! defor13 deformation term ( dw/dx + du/dz ); same as defor31
5048 ! defor22 deformation term ( dv/dy + dv/dy )
5049 ! defor23 deformation term ( dw/dy + dv/dz ); same as defor32
5050 ! defor33 deformation term ( dw/dz + dw/dz )
5051 ! div 3-d divergence
5052 ! dn
5053 ! dnw
5054 ! fnm
5055 ! fnp
5056 ! msftx
5057 ! msfty
5058 ! rdx
5059 ! rdy
5060 ! tendency
5061 ! titau tau (stress tensor) with a tilde, indicating division by
5062 ! a map-scale factor and the fraction of the total modeled
5063 ! atmosphere beneath a given altitude (titau = tau/m/zeta)
5064 ! tke turbulent kinetic energy
5065
5066 !-----------------------------------------------------------------------
5067 ! Begin declarations.
5068
5069 IMPLICIT NONE
5070
5071 TYPE( grid_config_rec_type ), INTENT( IN ) &
5072 :: config_flags
5073
5074 INTEGER, INTENT( IN ) &
5075 :: ids, ide, jds, jde, kds, kde, &
5076 ims, ime, jms, jme, kms, kme, &
5077 its, ite, jts, jte, kts, kte
5078
5079 REAL, INTENT( IN ) &
5080 :: cf1, cf2, cf3, rdx, rdy
5081
5082 REAL, DIMENSION( kms:kme ), INTENT( IN ) &
5083 :: fnm, fnp, dn, dnw
5084
5085 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
5086 :: msftx, msfty
5087
5088 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
5089 :: tendency
5090
5091 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
5092 :: defor11, defor22, defor33, defor12, defor13, defor23, &
5093 tke, xkmh, xkmv, zx, zy, u, v, w, rdz, rdzw
5094
5095 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
5096 :: mu
5097
5098 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
5099 :: ust
5100
5101 ! Local variables.
5102
5103 INTEGER &
5104 :: i, j, k, ktf, ktes1, ktes2, &
5105 i_start, i_end, j_start, j_end, &
5106 is_ext, ie_ext, js_ext, je_ext
5107
5108 REAL &
5109 :: mtau
5110
5111 REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ) &
5112 :: avg, titau, tmp2
5113
5114 REAL, DIMENSION( its:ite, kts:kte, jts:jte ) &
5115 :: titau12, tmp1, zxavg, zyavg
5116
5117 REAL :: absU, cd0, Cd
5118
5119 ! End declarations.
5120 !-----------------------------------------------------------------------
5121
5122 ktf = MIN( kte, kde-1 )
5123 ktes1 = kte-1
5124 ktes2 = kte-2
5125
5126 i_start = its
5127 i_end = MIN( ite, ide-1 )
5128 j_start = jts
5129 j_end = MIN( jte, jde-1 )
5130
5131 IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
5132 config_flags%nested ) i_start = MAX( ids+1, its )
5133 IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
5134 config_flags%nested ) i_end = MIN( ide-2, ite )
5135 IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
5136 config_flags%nested ) j_start = MAX( jds+1, jts )
5137 IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
5138 config_flags%nested ) j_end = MIN( jde-2, jte )
5139 IF ( config_flags%periodic_x ) i_start = its
5140 IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
5141
5142 DO j = j_start, j_end
5143 DO k = kts, ktf
5144 DO i = i_start, i_end
5145 zxavg(i,k,j) = 0.25 * ( zx(i,k ,j) + zx(i+1,k ,j) + &
5146 zx(i,k+1,j) + zx(i+1,k+1,j) )
5147 zyavg(i,k,j) = 0.25 * ( zy(i,k ,j) + zy(i,k ,j+1) + &
5148 zy(i,k+1,j) + zy(i,k+1,j+1) )
5149 END DO
5150 END DO
5151 END DO
5152
5153 ! Begin calculating production of turbulence due to shear. The approach
5154 ! is to add together contributions from six terms, each of which is the
5155 ! square of a deformation that is then multiplied by an exchange
5156 ! coefficiant. The same exchange coefficient is assumed for horizontal
5157 ! and vertical coefficients for some of the terms (the vertical value is
5158 ! the one used).
5159
5160 ! For defor11.
5161
5162 DO j = j_start, j_end
5163 DO k = kts, ktf
5164 DO i = i_start, i_end
5165 tendency(i,k,j) = tendency(i,k,j) + 0.5 * &
5166 mu(i,j) * xkmh(i,k,j) * ( ( defor11(i,k,j) )**2 )
5167 END DO
5168 END DO
5169 END DO
5170
5171 ! For defor22.
5172
5173 DO j = j_start, j_end
5174 DO k = kts, ktf
5175 DO i = i_start, i_end
5176 tendency(i,k,j) = tendency(i,k,j) + 0.5 * &
5177 mu(i,j) * xkmh(i,k,j) * ( ( defor22(i,k,j) )**2 )
5178 END DO
5179 END DO
5180 END DO
5181
5182 ! For defor33.
5183
5184 DO j = j_start, j_end
5185 DO k = kts, ktf
5186 DO i = i_start, i_end
5187 tendency(i,k,j) = tendency(i,k,j) + 0.5 * &
5188 mu(i,j) * xkmv(i,k,j) * ( ( defor33(i,k,j) )**2 )
5189 END DO
5190 END DO
5191 END DO
5192
5193 ! For defor12.
5194
5195 DO j = j_start, j_end
5196 DO k = kts, ktf
5197 DO i = i_start, i_end
5198 avg(i,k,j) = 0.25 * &
5199 ( ( defor12(i ,k,j)**2 ) + ( defor12(i ,k,j+1)**2 ) + &
5200 ( defor12(i+1,k,j)**2 ) + ( defor12(i+1,k,j+1)**2 ) )
5201 END DO
5202 END DO
5203 END DO
5204
5205 DO j = j_start, j_end
5206 DO k = kts, ktf
5207 DO i = i_start, i_end
5208 tendency(i,k,j) = tendency(i,k,j) + mu(i,j) * xkmh(i,k,j) * avg(i,k,j)
5209 END DO
5210 END DO
5211 END DO
5212
5213 ! For defor13.
5214
5215 DO j = j_start, j_end
5216 DO k = kts+1, ktf
5217 DO i = i_start, i_end+1
5218 tmp2(i,k,j) = defor13(i,k,j)
5219 END DO
5220 END DO
5221 END DO
5222
5223 DO j = j_start, j_end
5224 DO i = i_start, i_end+1
5225 tmp2(i,kts ,j) = 0.0
5226 tmp2(i,ktf+1,j) = 0.0
5227 END DO
5228 END DO
5229
5230 DO j = j_start, j_end
5231 DO k = kts, ktf
5232 DO i = i_start, i_end
5233 avg(i,k,j) = 0.25 * &
5234 ( ( tmp2(i ,k+1,j)**2 ) + ( tmp2(i ,k,j)**2 ) + &
5235 ( tmp2(i+1,k+1,j)**2 ) + ( tmp2(i+1,k,j)**2 ) )
5236 END DO
5237 END DO
5238 END DO
5239
5240 DO j = j_start, j_end
5241 DO k = kts, ktf
5242 DO i = i_start, i_end
5243 tendency(i,k,j) = tendency(i,k,j) + mu(i,j) * xkmv(i,k,j) * avg(i,k,j)
5244 END DO
5245 END DO
5246 END DO
5247
5248 !MARTA: add the drag at the surface; WCS 040331
5249 K=KTS
5250
5251 uflux: SELECT CASE( config_flags%isfflx )
5252 CASE (0) ! Assume cd a constant, specified in namelist
5253
5254 cd0 = config_flags%tke_drag_coefficient ! drag coefficient set
5255 ! in namelist.input
5256 DO j = j_start, j_end
5257 DO i = i_start, i_end
5258
5259 absU=0.5*sqrt((u(i,k,j)+u(i+1,k,j))**2+(v(i,k,j)+v(i,k,j+1))**2)
5260 Cd = cd0
5261 tendency(i,k,j) = tendency(i,k,j) + &
5262 mu(i,j)*( (u(i,k,j)+u(i+1,k,j))*0.5* &
5263 Cd*absU*(defor13(i,kts+1,j)+defor13(i+1,kts+1,j))*0.5 )
5264
5265 END DO
5266 END DO
5267
5268 CASE (1) ! ustar computed from surface routine
5269
5270 DO j = j_start, j_end
5271 DO i = i_start, i_end
5272
5273 absU=0.5*sqrt((u(i,k,j)+u(i+1,k,j))**2+(v(i,k,j)+v(i,k,j+1))**2)+epsilon
5274 Cd = (ust(i,j)**2)/(absU**2)
5275 tendency(i,k,j) = tendency(i,k,j) + &
5276 mu(i,j)*( (u(i,k,j)+u(i+1,k,j))*0.5* &
5277 Cd*absU*(defor13(i,kts+1,j)+defor13(i+1,kts+1,j))*0.5 )
5278
5279 END DO
5280 END DO
5281
5282 CASE DEFAULT
5283 CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' )
5284 END SELECT uflux
5285 ! end of MARTA/WCS change
5286
5287 ! For defor23.
5288
5289 DO j = j_start, j_end+1
5290 DO k = kts+1, ktf
5291 DO i = i_start, i_end
5292 tmp2(i,k,j) = defor23(i,k,j)
5293 END DO
5294 END DO
5295 END DO
5296
5297 DO j = j_start, j_end+1
5298 DO i = i_start, i_end
5299 tmp2(i,kts, j) = 0.0
5300 tmp2(i,ktf+1,j) = 0.0
5301 END DO
5302 END DO
5303
5304 DO j = j_start, j_end
5305 DO k = kts, ktf
5306 DO i = i_start, i_end
5307 avg(i,k,j) = 0.25 * &
5308 ( ( tmp2(i,k+1,j )**2 ) + ( tmp2(i,k,j )**2) + &
5309 ( tmp2(i,k+1,j+1)**2 ) + ( tmp2(i,k,j+1)**2) )
5310 END DO
5311 END DO
5312 END DO
5313
5314 DO j = j_start, j_end
5315 DO k = kts, ktf
5316 DO i = i_start, i_end
5317 tendency(i,k,j) = tendency(i,k,j) + mu(i,j) * xkmv(i,k,j) * avg(i,k,j)
5318 END DO
5319 END DO
5320 END DO
5321
5322 !MARTA: add the drag at the surface; WCS 040331
5323 K=KTS
5324
5325 vflux: SELECT CASE( config_flags%isfflx )
5326 CASE (0) ! Assume cd a constant, specified in namelist
5327
5328 cd0 = config_flags%tke_drag_coefficient ! constant drag coefficient
5329 ! set in namelist.input
5330 DO j = j_start, j_end
5331 DO i = i_start, i_end
5332
5333 absU=0.5*sqrt((u(i,k,j)+u(i+1,k,j))**2+(v(i,k,j)+v(i,k,j+1))**2)
5334 Cd = cd0
5335 tendency(i,k,j) = tendency(i,k,j) + &
5336 mu(i,j)*( (v(i,k,j)+v(i,k,j+1))*0.5* &
5337 Cd*absU*(defor23(i,kts+1,j)+defor23(i,kts+1,j+1))*0.5 )
5338
5339 END DO
5340 END DO
5341
5342 CASE (1) ! ustar computed from surface routine
5343
5344 DO j = j_start, j_end
5345 DO i = i_start, i_end
5346
5347 absU=0.5*sqrt((u(i,k,j)+u(i+1,k,j))**2+(v(i,k,j)+v(i,k,j+1))**2)+epsilon
5348 Cd = (ust(i,j)**2)/(absU**2)
5349 tendency(i,k,j) = tendency(i,k,j) + &
5350 mu(i,j)*( (v(i,k,j)+v(i,k,j+1))*0.5* &
5351 Cd*absU*(defor23(i,kts+1,j)+defor23(i,kts+1,j+1))*0.5 )
5352
5353 END DO
5354 END DO
5355
5356 CASE DEFAULT
5357 CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' )
5358 END SELECT vflux
5359 ! end of MARTA/WCS change
5360
5361 END SUBROUTINE tke_shear
5362
5363 !=======================================================================
5364 !=======================================================================
5365
5366 SUBROUTINE compute_diff_metrics( config_flags, ph, phb, z, rdz, rdzw, &
5367 zx, zy, rdx, rdy, &
5368 ids, ide, jds, jde, kds, kde, &
5369 ims, ime, jms, jme, kms, kme, &
5370 its, ite, jts, jte, kts, kte )
5371
5372 !-----------------------------------------------------------------------
5373 ! Begin declarations.
5374
5375 IMPLICIT NONE
5376
5377 TYPE( grid_config_rec_type ), INTENT( IN ) &
5378 :: config_flags
5379
5380 INTEGER, INTENT( IN ) &
5381 :: ids, ide, jds, jde, kds, kde, &
5382 ims, ime, jms, jme, kms, kme, &
5383 its, ite, jts, jte, kts, kte
5384
5385 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
5386 :: ph, phb
5387
5388 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( OUT ) &
5389 :: rdz, rdzw, zx, zy, z
5390
5391 REAL, INTENT( IN ) &
5392 :: rdx, rdy
5393
5394 ! Local variables.
5395
5396 INTEGER &
5397 :: i, j, k, i_start, i_end, j_start, j_end, ktf
5398
5399 ! End declarations.
5400 !-----------------------------------------------------------------------
5401
5402 ktf = MIN( kte, kde-1 )
5403
5404 ! Bug fix, WCS, 22 april 2002.
5405
5406 ! We need rdzw in halo for average to u and v points.
5407
5408 j_start = jts-1
5409 j_end = jte
5410
5411 ! Begin with dz computations.
5412
5413 DO j = j_start, j_end
5414
5415 IF ( ( j_start >= jts ) .AND. ( j_end <= MIN( jte, jde-1 ) ) ) THEN
5416 i_start = its-1
5417 i_end = ite
5418 ELSE
5419 i_start = its
5420 i_end = MIN( ite, ide-1 )
5421 END IF
5422
5423 ! Compute z at w points for rdz and rdzw computations. We'll switch z
5424 ! to z at p points before returning
5425
5426 DO k = 1, kte
5427
5428 ! Bug fix, WCS, 22 april 2002
5429
5430 DO i = i_start, i_end
5431 z(i,k,j) = ( ph(i,k,j) + phb(i,k,j) ) / g
5432 END DO
5433 END DO
5434
5435 DO k = 1, ktf
5436 DO i = i_start, i_end
5437 rdzw(i,k,j) = 1.0 / ( z(i,k+1,j) - z(i,k,j) )
5438 END DO
5439 END DO
5440
5441 DO k = 2, ktf
5442 DO i = i_start, i_end
5443 rdz(i,k,j) = 2.0 / ( z(i,k+1,j) - z(i,k-1,j) )
5444 END DO
5445 END DO
5446
5447 ! Bug fix, WCS, 22 april 2002; added the following code
5448
5449 DO i = i_start, i_end
5450 rdz(i,1,j) = 2./(z(i,2,j)-z(i,1,j))
5451 END DO
5452
5453 END DO
5454
5455 ! End bug fix.
5456
5457 ! Now compute zx and zy; we'll assume that the halo for ph and phb is
5458 ! properly filled.
5459
5460 i_start = its
5461 i_end = MIN( ite, ide-1 )
5462 j_start = jts
5463 j_end = MIN( jte, jde-1 )
5464
5465 DO j = j_start, j_end
5466 DO k = 1, kte
5467 DO i = MAX( ids+1, its ), i_end
5468 zx(i,k,j) = rdx * ( phb(i,k,j) - phb(i-1,k,j) ) / g
5469 END DO
5470 END DO
5471 END DO
5472
5473 DO j = j_start, j_end
5474 DO k = 1, kte
5475 DO i = MAX( ids+1, its ), i_end
5476 zx(i,k,j) = zx(i,k,j) + rdx * ( ph(i,k,j) - ph(i-1,k,j) ) / g
5477 END DO
5478 END DO
5479 END DO
5480
5481 DO j = MAX( jds+1, jts ), j_end
5482 DO k = 1, kte
5483 DO i = i_start, i_end
5484 zy(i,k,j) = rdy * ( phb(i,k,j) - phb(i,k,j-1) ) / g
5485 END DO
5486 END DO
5487 END DO
5488
5489 DO j = MAX( jds+1, jts ), j_end
5490 DO k = 1, kte
5491 DO i = i_start, i_end
5492 zy(i,k,j) = zy(i,k,j) + rdy * ( ph(i,k,j) - ph(i,k,j-1) ) / g
5493 END DO
5494 END DO
5495 END DO
5496
5497 ! Some b.c. on zx and zy.
5498
5499 IF ( .NOT. config_flags%periodic_x ) THEN
5500
5501 IF ( ite == ide ) THEN
5502 DO j = j_start, j_end
5503 DO k = 1, ktf
5504 zx(ide,k,j) = 0.0
5505 END DO
5506 END DO
5507 END IF
5508
5509 IF ( its == ids ) THEN
5510 DO j = j_start, j_end
5511 DO k = 1, ktf
5512 zx(ids,k,j) = 0.0
5513 END DO
5514 END DO
5515 END IF
5516
5517 ELSE
5518
5519 IF ( ite == ide ) THEN
5520 DO j=j_start,j_end
5521 DO k=1,ktf
5522 zx(ide,k,j) = rdx * ( phb(ide,k,j) - phb(ide-1,k,j) ) / g
5523 END DO
5524 END DO
5525
5526 DO j = j_start, j_end
5527 DO k = 1, ktf
5528 zx(ide,k,j) = zx(ide,k,j) + rdx * ( ph(ide,k,j) - ph(ide-1,k,j) ) / g
5529 END DO
5530 END DO
5531 END IF
5532
5533 IF ( its == ids ) THEN
5534 DO j = j_start, j_end
5535 DO k = 1, ktf
5536 zx(ids,k,j) = rdx * ( phb(ids,k,j) - phb(ids-1,k,j) ) / g
5537 END DO
5538 END DO
5539
5540 DO j =j_start,j_end
5541 DO k =1,ktf
5542 zx(ids,k,j) = zx(ids,k,j) + rdx * ( ph(ids,k,j) - ph(ids-1,k,j) ) / g
5543 END DO
5544 END DO
5545 END IF
5546
5547 END IF
5548
5549 IF ( .NOT. config_flags%periodic_y ) THEN
5550
5551 IF ( jte == jde ) THEN
5552 DO k =1, ktf
5553 DO i =i_start, i_end
5554 zy(i,k,jde) = 0.0
5555 END DO
5556 END DO
5557 END IF
5558
5559 IF ( jts == jds ) THEN
5560 DO k =1, ktf
5561 DO i =i_start, i_end
5562 zy(i,k,jds) = 0.0
5563 END DO
5564 END DO
5565 END IF
5566
5567 ELSE
5568
5569 IF ( jte == jde ) THEN
5570 DO j=j_start, j_end
5571 DO k=1, ktf
5572 zy(i,k,jde) = rdy * ( phb(i,k,jde) - phb(i,k,jde-1) ) / g
5573 END DO
5574 END DO
5575
5576 DO j = j_start, j_end
5577 DO k = 1, ktf
5578 zy(i,k,jde) = zy(i,k,jde) + rdy * ( ph(i,k,jde) - ph(i,k,jde-1) ) / g
5579 END DO
5580 END DO
5581 END IF
5582
5583 IF ( jts == jds ) THEN
5584 DO j = j_start, j_end
5585 DO k = 1, ktf
5586 zy(i,k,jds) = rdy * ( phb(i,k,jds) - phb(i,k,jds-1) ) / g
5587 END DO
5588 END DO
5589
5590 DO j = j_start, j_end
5591 DO k = 1, ktf
5592 zy(i,k,jds) = zy(i,k,jds) + rdy * ( ph(i,k,jds) - ph(i,k,jds-1) ) / g
5593 END DO
5594 END DO
5595 END IF
5596
5597 END IF
5598
5599 ! Calculate z at p points.
5600
5601 DO j = j_start, j_end
5602 DO k = 1, ktf
5603 DO i = i_start, i_end
5604 z(i,k,j) = 0.5 * &
5605 ( ph(i,k,j) + phb(i,k,j) + ph(i,k+1,j) + phb(i,k+1,j) ) / g
5606 END DO
5607 END DO
5608 END DO
5609
5610 END SUBROUTINE compute_diff_metrics
5611
5612 !=======================================================================
5613 !=======================================================================
5614
5615 END MODULE module_diffusion_em
5616
5617 !=======================================================================
5618 !=======================================================================
5619