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