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,ktf+1,j)=0.
2616 titau2avg(i,ktf+1,j)=0.
2617 ENDDO
2618 ENDDO
2619
2620 DO j = j_start, j_end
2621 DO k = kts+1,ktf
2622 DO i = i_start, i_end
2623
2624 mrdx=msft(i,j)*rdx
2625 mrdy=msft(i,j)*rdy
2626
2627 tendency(i,k,j)=tendency(i,k,j)- &
2628 (mrdx*(titau1(i+1,k,j)-titau1(i,k,j))+ &
2629 mrdy*(titau2(i,k,j+1)-titau2(i,k,j))- &
2630 msft(i,j)*rdz(i,k,j)*(titau1avg(i,k,j)-titau1avg(i,k-1,j)+ &
2631 titau2avg(i,k,j)-titau2avg(i,k-1,j) &
2632 ) &
2633 )
2634 ! msft(i,j)/dn(k)*(titau1avg(i,k,j)-titau1avg(i,k-1,j)+ &
2635 ! titau2avg(i,k,j)-titau2avg(i,k-1,j) &
2636 ! ) &
2637 ! )
2638 ENDDO
2639 ENDDO
2640 ENDDO
2641
2642 END SUBROUTINE horizontal_diffusion_w_2
2643
2644 !=======================================================================
2645 !=======================================================================
2646
2647 SUBROUTINE horizontal_diffusion_s (tendency, mu, config_flags, var, &
2648 msft, msfu, msfv, xkhh, rdx, rdy, &
2649 fnm, fnp, cf1, cf2, cf3, &
2650 zx, zy, rdz, rdzw, dn, dnw, &
2651 doing_tke, &
2652 ids, ide, jds, jde, kds, kde, &
2653 ims, ime, jms, jme, kms, kme, &
2654 its, ite, jts, jte, kts, kte )
2655
2656 !-----------------------------------------------------------------------
2657 ! Begin declarations.
2658
2659 IMPLICIT NONE
2660
2661 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2662
2663 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2664 ims, ime, jms, jme, kms, kme, &
2665 its, ite, jts, jte, kts, kte
2666
2667 LOGICAL, INTENT(IN ) :: doing_tke
2668
2669 REAL , INTENT(IN ) :: cf1, cf2, cf3
2670
2671 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
2672 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
2673 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
2674 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
2675
2676 REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfu
2677 REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfv
2678 REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msft
2679
2680 REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: mu
2681
2682 ! REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1), &
2683 ! INTENT(IN ) :: xkhh
2684
2685 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: tendency
2686
2687 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: &
2688 xkhh, &
2689 rdz, &
2690 rdzw
2691
2692 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: var, &
2693 zx, &
2694 zy
2695
2696 REAL , INTENT(IN ) :: rdx, &
2697 rdy
2698
2699 ! Local data
2700
2701 INTEGER :: i, j, k, ktf
2702
2703 INTEGER :: i_start, i_end, j_start, j_end
2704
2705 REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: H1avg, &
2706 H2avg, &
2707 H1, &
2708 H2, &
2709 xkxavg
2710 ! new
2711 ! zxavg, &
2712 ! zyavg
2713
2714 REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: tmptendf
2715
2716 REAL :: mrdx, mrdy, rcoup
2717 REAL :: tmpzx, tmpzy, tmpzeta_z
2718 INTEGER :: ktes1,ktes2
2719
2720 ! End declarations.
2721 !-----------------------------------------------------------------------
2722
2723 ktf=MIN(kte,kde-1)
2724
2725 !-----------------------------------------------------------------------
2726 ! scalars: t (.), u(|), v(+), w(-)
2727 !
2728 ! t u t u t v t v
2729 !
2730 ! w - 3 - k+1 w - 3 - k+1
2731 !
2732 ! t . 1 O 1 . k t . 2 O 2 . k
2733 !
2734 ! w - 3 - k w - 3 - k
2735 !
2736 ! t . | . | . k-1 t . + . + . k-1
2737 !
2738 ! w - - - k-1 w - - - k-1
2739 !
2740 ! t i-1 i i i+1 j-1 j j j+1
2741 !
2742
2743 ktes1=kte-1
2744 ktes2=kte-2
2745
2746 i_start = its
2747 i_end = MIN(ite,ide-1)
2748 j_start = jts
2749 j_end = MIN(jte,jde-1)
2750
2751 IF ( config_flags%open_xs .or. config_flags%specified .or. &
2752 config_flags%nested) i_start = MAX(ids+1,its)
2753 IF ( config_flags%open_xe .or. config_flags%specified .or. &
2754 config_flags%nested) i_end = MIN(ide-2,ite)
2755 IF ( config_flags%open_ys .or. config_flags%specified .or. &
2756 config_flags%nested) j_start = MAX(jds+1,jts)
2757 IF ( config_flags%open_ye .or. config_flags%specified .or. &
2758 config_flags%nested) j_end = MIN(jde-2,jte)
2759 IF ( config_flags%periodic_x ) i_start = its
2760 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2761
2762 ! diffusion of the TKE needs mutiple 2
2763
2764 IF ( doing_tke ) THEN
2765 DO j = j_start, j_end
2766 DO k = kts,ktf
2767 DO i = i_start, i_end
2768 tmptendf(i,k,j)=tendency(i,k,j)
2769 ENDDO
2770 ENDDO
2771 ENDDO
2772 ENDIF
2773
2774 ! H1 = partial var over partial x
2775
2776 DO j = j_start, j_end
2777 DO k = kts, ktf
2778 DO i = i_start, i_end + 1
2779 ! new
2780 ! zxavg(i,k,j) =0.5*( zx(i-1,k,j)+ zx(i,k,j))
2781 xkxavg(i,k,j)=0.5*(xkhh(i-1,k,j)+xkhh(i,k,j))
2782 ENDDO
2783 ENDDO
2784 ENDDO
2785
2786 DO j = j_start, j_end
2787 DO k = kts+1, ktf
2788 DO i = i_start, i_end + 1
2789 H1avg(i,k,j)=0.5*(fnm(k)*(var(i-1,k ,j)+var(i,k ,j))+ &
2790 fnp(k)*(var(i-1,k-1,j)+var(i,k-1,j)))
2791 ENDDO
2792 ENDDO
2793 ENDDO
2794
2795 DO j = j_start, j_end
2796 DO i = i_start, i_end + 1
2797 H1avg(i,kts ,j)=0.5*(cf1*var(i ,1,j)+cf2*var(i ,2,j)+ &
2798 cf3*var(i ,3,j)+cf1*var(i-1,1,j)+ &
2799 cf2*var(i-1,2,j)+cf3*var(i-1,3,j))
2800 H1avg(i,ktf+1,j)=0.5*(var(i,ktes1,j)+(var(i,ktes1,j)- &
2801 var(i,ktes2,j))*0.5*dnw(ktes1)/dn(ktes1)+ &
2802 var(i-1,ktes1,j)+(var(i-1,ktes1,j)- &
2803 var(i-1,ktes2,j))*0.5*dnw(ktes1)/dn(ktes1))
2804 ENDDO
2805 ENDDO
2806
2807 DO j = j_start, j_end
2808 DO k = kts, ktf
2809 DO i = i_start, i_end + 1
2810 ! new
2811 tmpzx = 0.5*( zx(i,k,j)+ zx(i,k+1,j))
2812 H1(i,k,j)=-msfu(i,j)*xkxavg(i,k,j)*( &
2813 rdx*(var(i,k,j)-var(i-1,k,j)) - tmpzx* &
2814 (H1avg(i,k+1,j)-H1avg(i,k,j))*rdzw(i,k,j) )
2815
2816 ! tmpzeta_z = 0.5*(zeta_z(i,j)+zeta_z(i-1,j))
2817 ! H1(i,k,j)=-msfu(i,j)*xkxavg(i,k,j)*( &
2818 ! rdx*(var(i,k,j)-var(i-1,k,j)) - tmpzx*tmpzeta_z* &
2819 ! (H1avg(i,k+1,j)-H1avg(i,k,j))/dnw(k))
2820 ENDDO
2821 ENDDO
2822 ENDDO
2823
2824 ! H2 = partial var over partial y
2825
2826 DO j = j_start, j_end + 1
2827 DO k = kts, ktf
2828 DO i = i_start, i_end
2829 ! new
2830 ! zyavg(i,k,j) =0.5*( zy(i,k,j-1)+ zy(i,k,j))
2831 xkxavg(i,k,j)=0.5*(xkhh(i,k,j-1)+xkhh(i,k,j))
2832 ENDDO
2833 ENDDO
2834 ENDDO
2835
2836 DO j = j_start, j_end + 1
2837 DO k = kts+1, ktf
2838 DO i = i_start, i_end
2839 ! new
2840 H2avg(i,k,j)=0.5*(fnm(k)*(var(i,k ,j-1)+var(i,k ,j))+ &
2841 fnp(k)*(var(i,k-1,j-1)+var(i,k-1,j)))
2842 ENDDO
2843 ENDDO
2844 ENDDO
2845
2846 DO j = j_start, j_end + 1
2847 DO i = i_start, i_end
2848 H2avg(i,kts ,j)=0.5*(cf1*var(i,1,j )+cf2*var(i ,2,j)+ &
2849 cf3*var(i,3,j )+cf1*var(i,1,j-1)+ &
2850 cf2*var(i,2,j-1)+cf3*var(i,3,j-1))
2851 H2avg(i,ktf+1,j)=0.5*(var(i,ktes1,j)+(var(i,ktes1,j)- &
2852 var(i,ktes2,j))*0.5*dnw(ktes1)/dn(ktes1)+ &
2853 var(i,ktes1,j-1)+(var(i,ktes1,j-1)- &
2854 var(i,ktes2,j-1))*0.5*dnw(ktes1)/dn(ktes1))
2855 ENDDO
2856 ENDDO
2857
2858 DO j = j_start, j_end + 1
2859 DO k = kts, ktf
2860 DO i = i_start, i_end
2861 ! new
2862 tmpzy = 0.5*( zy(i,k,j)+ zy(i,k+1,j))
2863
2864 H2(i,k,j)=-msfv(i,j)*xkxavg(i,k,j)*( &
2865 rdy*(var(i,k,j)-var(i,k,j-1)) - tmpzy* &
2866 (H2avg(i ,k+1,j)-H2avg(i,k,j))*rdzw(i,k,j))
2867
2868 ! tmpzeta_z = 0.5*(zeta_z(i,j)+zeta_z(i,j-1))
2869 ! H2(i,k,j)=-msfv(i,j)*xkxavg(i,k,j)*( &
2870 ! rdy*(var(i,k,j)-var(i,k,j-1)) - tmpzy*tmpzeta_z* &
2871 ! (H2avg(i ,k+1,j)-H2avg(i,k,j))/dnw(k))
2872 ENDDO
2873 ENDDO
2874 ENDDO
2875
2876 DO j = j_start, j_end
2877 DO k = kts+1, ktf
2878 DO i = i_start, i_end
2879 H1avg(i,k,j)=0.5*(fnm(k)*(H1(i+1,k ,j)+H1(i,k ,j))+ &
2880 fnp(k)*(H1(i+1,k-1,j)+H1(i,k-1,j)))
2881 H2avg(i,k,j)=0.5*(fnm(k)*(H2(i,k ,j+1)+H2(i,k ,j))+ &
2882 fnp(k)*(H2(i,k-1,j+1)+H2(i,k-1,j)))
2883 ! new
2884 ! zxavg(i,k,j)=fnm(k)*zx(i,k,j)+fnp(k)*zx(i,k-1,j)
2885 ! zyavg(i,k,j)=fnm(k)*zy(i,k,j)+fnp(k)*zy(i,k-1,j)
2886
2887 ! H1avg(i,k,j)=zx*H1avg*zeta_z
2888 ! H2avg(i,k,j)=zy*H2avg*zeta_z
2889
2890 tmpzx = 0.5*( zx(i,k,j)+ zx(i+1,k,j ))
2891 tmpzy = 0.5*( zy(i,k,j)+ zy(i ,k,j+1))
2892
2893 H1avg(i,k,j)=H1avg(i,k,j)*tmpzx
2894 H2avg(i,k,j)=H2avg(i,k,j)*tmpzy
2895
2896 ! H1avg(i,k,j)=H1avg(i,k,j)*tmpzx*zeta_z(i,j)
2897 ! H2avg(i,k,j)=H2avg(i,k,j)*tmpzy*zeta_z(i,j)
2898 ENDDO
2899 ENDDO
2900 ENDDO
2901
2902 DO j = j_start, j_end
2903 DO i = i_start, i_end
2904 H1avg(i,kts ,j)=0.
2905 H1avg(i,ktf+1,j)=0.
2906 H2avg(i,kts ,j)=0.
2907 H2avg(i,ktf+1,j)=0.
2908 ENDDO
2909 ENDDO
2910
2911 DO j = j_start, j_end
2912 DO k = kts,ktf
2913 DO i = i_start, i_end
2914
2915 mrdx=msft(i,j)*rdx
2916 mrdy=msft(i,j)*rdy
2917
2918 tendency(i,k,j)=tendency(i,k,j)- &
2919 (mrdx*0.5*((mu(i+1,j)+mu(i,j))*H1(i+1,k,j)- &
2920 (mu(i-1,j)+mu(i,j))*H1(i ,k,j))+ &
2921 mrdy*0.5*((mu(i,j+1)+mu(i,j))*H2(i,k,j+1)- &
2922 (mu(i,j-1)+mu(i,j))*H2(i,k,j ))- &
2923 msft(i,j)*mu(i,j)*(H1avg(i,k+1,j)-H1avg(i,k,j)+ &
2924 H2avg(i,k+1,j)-H2avg(i,k,j) &
2925 )*rdzw(i,k,j) &
2926 )
2927
2928 ENDDO
2929 ENDDO
2930 ENDDO
2931
2932 IF ( doing_tke ) THEN
2933 DO j = j_start, j_end
2934 DO k = kts,ktf
2935 DO i = i_start, i_end
2936 tendency(i,k,j)=tmptendf(i,k,j)+2.* &
2937 (tendency(i,k,j)-tmptendf(i,k,j))
2938 ENDDO
2939 ENDDO
2940 ENDDO
2941 ENDIF
2942
2943 END SUBROUTINE horizontal_diffusion_s
2944
2945 !=======================================================================
2946 !=======================================================================
2947
2948 SUBROUTINE vertical_diffusion_2 ( ru_tendf, rv_tendf, rw_tendf, rt_tendf, &
2949 tke_tendf, moist_tendf, n_moist, &
2950 chem_tendf, n_chem, &
2951 scalar_tendf, n_scalar, &
2952 u_2, v_2, &
2953 thp,u_base,v_base,t_base,qv_base,mu,tke, &
2954 config_flags,defor13,defor23,defor33,div, &
2955 moist, chem, scalar, xkmv, xkhv,km_opt, &
2956 fnm, fnp, dn, dnw, rdz, rdzw, &
2957 ids, ide, jds, jde, kds, kde, &
2958 ims, ime, jms, jme, kms, kme, &
2959 its, ite, jts, jte, kts, kte )
2960
2961 !-----------------------------------------------------------------------
2962 ! Begin declarations.
2963
2964 IMPLICIT NONE
2965
2966 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2967
2968 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2969 ims, ime, jms, jme, kms, kme, &
2970 its, ite, jts, jte, kts, kte
2971
2972 INTEGER , INTENT(IN ) :: n_moist, n_chem, n_scalar, km_opt
2973
2974 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
2975 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
2976 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
2977 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
2978 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
2979
2980 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: qv_base
2981 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: u_base
2982 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: v_base
2983 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: t_base
2984
2985 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::ru_tendf,&
2986 rv_tendf,&
2987 rw_tendf,&
2988 tke_tendf,&
2989 rt_tendf
2990
2991 REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), &
2992 INTENT(INOUT) :: moist_tendf
2993
2994 REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_chem), &
2995 INTENT(INOUT) :: chem_tendf
2996
2997 REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_scalar) , &
2998 INTENT(INOUT) :: scalar_tendf
2999
3000 REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), &
3001 INTENT(INOUT) :: moist
3002
3003 REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_chem), &
3004 INTENT(INOUT) :: chem
3005
3006 REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_scalar) , &
3007 INTENT(IN ) :: scalar
3008
3009 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor13, &
3010 defor23, &
3011 defor33, &
3012 div, &
3013 xkmv, &
3014 xkhv, &
3015 tke, &
3016 rdz, &
3017 u_2, &
3018 v_2, &
3019 rdzw
3020
3021 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: thp
3022
3023 ! LOCAL VAR
3024
3025 REAL , DIMENSION( ims:ime, kms:kme, jms:jme) :: var_mix
3026
3027 INTEGER :: im, i,j,k
3028 INTEGER :: i_start, i_end, j_start, j_end
3029
3030 ! REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: xkhv
3031
3032 !***************************************************************************
3033 !***************************************************************************
3034 !MODIFICA VARIABILI PER I FLUSSI
3035 !
3036 REAL , DIMENSION( ims:ime, jms:jme) :: Cd
3037 REAL :: V0_u,V0_v,tao_xz,tao_yz,ustar,cd0
3038 REAL :: xsfc,psi1,vk2,zrough,lnz
3039 REAL :: heat_flux
3040 !
3041 !FINE MODIFICA VARIABILI PER I FLUSSI
3042 !***************************************************************************
3043 !
3044
3045 ! End declarations.
3046 !-----------------------------------------------------------------------
3047
3048 i_start = its
3049 i_end = MIN(ite,ide-1)
3050 j_start = jts
3051 j_end = MIN(jte,jde-1)
3052 !
3053 !-----------------------------------------------------------------------
3054
3055 CALL vertical_diffusion_u_2( ru_tendf, config_flags, mu, &
3056 defor13, xkmv, &
3057 dnw, rdzw, fnm, fnp, &
3058 ids, ide, jds, jde, kds, kde, &
3059 ims, ime, jms, jme, kms, kme, &
3060 its, ite, jts, jte, kts, kte )
3061
3062
3063 CALL vertical_diffusion_v_2( rv_tendf, config_flags, mu, &
3064 defor23, xkmv, &
3065 dnw, rdzw, fnm, fnp, &
3066 ids, ide, jds, jde, kds, kde, &
3067 ims, ime, jms, jme, kms, kme, &
3068 its, ite, jts, jte, kts, kte )
3069
3070 CALL vertical_diffusion_w_2( rw_tendf, config_flags, mu, &
3071 defor33, tke(ims,kms,jms), &
3072 div, xkmv, &
3073 dn, rdz, &
3074 ids, ide, jds, jde, kds, kde, &
3075 ims, ime, jms, jme, kms, kme, &
3076 its, ite, jts, jte, kts, kte )
3077
3078 !*****************************************
3079 !*****************************************
3080 ! MODIFICA al flusso di momento alla parete
3081 !
3082 cd0 = config_flags%tke_drag_coefficient ! constant drag coefficient
3083 ! set in namelist.input
3084 DO j = j_start, j_end+1
3085 DO i = i_start, i_end+1
3086 Cd(i,j)= cd0
3087 ENDDO
3088 ENDDO
3089 !
3090 !calcolo del modulo della velocita
3091 DO j = j_start, j_end
3092 DO i = i_start, i_end+1
3093 V0_u=0.
3094 tao_xz=0.
3095 V0_u= sqrt((u_2(i,kts,j)**2) + &
3096 (((v_2(i ,kts,j )+ &
3097 v_2(i ,kts,j+1)+ &
3098 v_2(i-1,kts,j )+ &
3099 v_2(i-1,kts,j+1))/4)**2))+epsilon
3100 tao_xz=Cd(i,j)*V0_u*u_2(i,kts,j)
3101 ru_tendf(i,kts,j)=ru_tendf(i,kts,j) &
3102 -0.25*(mu(i,j)+mu(i-1,j))*tao_xz*(rdzw(i,kts,j)+rdzw(i-1,kts,j))
3103 ENDDO
3104 ENDDO
3105
3106 !
3107 DO j = j_start, j_end+1
3108 DO i = i_start, i_end
3109 V0_v=0.
3110 tao_yz=0.
3111 V0_v= sqrt((v_2(i,kts,j)**2) + &
3112 (((u_2(i ,kts,j )+ &
3113 u_2(i ,kts,j-1)+ &
3114 u_2(i+1,kts,j )+ &
3115 u_2(i+1,kts,j-1))/4)**2))+epsilon
3116 tao_yz=Cd(i,j)*V0_v*v_2(i,kts,j)
3117 rv_tendf(i,kts,j)=rv_tendf(i,kts,j) &
3118 -0.25*(mu(i,j)+mu(i,j-1))*tao_yz*(rdzw(i,kts,j)+rdzw(i,kts,j-1))
3119 ENDDO
3120 ENDDO
3121 !
3122 ! FINE MODIFICA al flusso di momento alla parete
3123 !*****************************************
3124 !*****************************************
3125
3126 IF ( config_flags%mix_full_fields ) THEN
3127
3128 DO j=jts,min(jte,jde-1)
3129 DO k=kts,kte-1
3130 DO i=its,min(ite,ide-1)
3131 var_mix(i,k,j) = thp(i,k,j)
3132 ENDDO
3133 ENDDO
3134 ENDDO
3135
3136 ELSE
3137
3138 DO j=jts,min(jte,jde-1)
3139 DO k=kts,kte-1
3140 DO i=its,min(ite,ide-1)
3141 var_mix(i,k,j) = thp(i,k,j) - t_base(k)
3142 ENDDO
3143 ENDDO
3144 ENDDO
3145
3146 END IF
3147
3148 CALL vertical_diffusion_s( rt_tendf, config_flags, var_mix, mu, xkhv, &
3149 dn, dnw, rdz, rdzw, fnm, fnp, &
3150 .false., &
3151 ids, ide, jds, jde, kds, kde, &
3152 ims, ime, jms, jme, kms, kme, &
3153 its, ite, jts, jte, kts, kte )
3154
3155
3156 !*****************************************
3157 !*****************************************
3158 !MODIFICA al flusso di calore
3159 !
3160 !
3161 heat_flux = config_flags%tke_heat_flux ! constant heat flux value
3162 ! set in namelist.input
3163 DO j = j_start, j_end
3164 DO i = i_start, i_end
3165
3166 rt_tendf(i,kts,j)=rt_tendf(i,kts,j) &
3167 +mu(i,j)*heat_flux*rdzw(i,kts,j)
3168
3169 ENDDO
3170 ENDDO
3171 !
3172 ! FINE MODIFICA al flusso di calore
3173 !*****************************************
3174 !*****************************************
3175
3176 If (km_opt .eq. 2) then
3177 CALL vertical_diffusion_s( tke_tendf(ims,kms,jms), &
3178 config_flags, tke(ims,kms,jms), &
3179 mu, xkhv, &
3180 dn, dnw, rdz, rdzw, fnm, fnp, &
3181 .true., &
3182 ids, ide, jds, jde, kds, kde, &
3183 ims, ime, jms, jme, kms, kme, &
3184 its, ite, jts, jte, kts, kte )
3185 endif
3186
3187 IF (n_moist .ge. PARAM_FIRST_SCALAR) THEN
3188
3189 moist_loop: do im = PARAM_FIRST_SCALAR, n_moist
3190
3191 IF ( (.not. config_flags%mix_full_fields) .and. (im == P_QV) ) THEN
3192
3193 DO j=jts,min(jte,jde-1)
3194 DO k=kts,kte-1
3195 DO i=its,min(ite,ide-1)
3196 var_mix(i,k,j) = moist(i,k,j,im) - qv_base(k)
3197 ENDDO
3198 ENDDO
3199 ENDDO
3200
3201 ELSE
3202
3203 DO j=jts,min(jte,jde-1)
3204 DO k=kts,kte-1
3205 DO i=its,min(ite,ide-1)
3206 var_mix(i,k,j) = moist(i,k,j,im)
3207 ENDDO
3208 ENDDO
3209 ENDDO
3210
3211 END IF
3212
3213
3214 CALL vertical_diffusion_s( moist_tendf(ims,kms,jms,im), &
3215 config_flags, var_mix, &
3216 mu, xkhv, &
3217 dn, dnw, rdz, rdzw, fnm, fnp, &
3218 .false., &
3219 ids, ide, jds, jde, kds, kde, &
3220 ims, ime, jms, jme, kms, kme, &
3221 its, ite, jts, jte, kts, kte )
3222
3223 ENDDO moist_loop
3224
3225 ENDIF
3226
3227 IF (n_chem .ge. PARAM_FIRST_SCALAR) THEN
3228
3229 chem_loop: do im = PARAM_FIRST_SCALAR, n_chem
3230
3231 CALL vertical_diffusion_s( chem_tendf(ims,kms,jms,im), &
3232 config_flags, chem(ims,kms,jms,im), &
3233 mu, xkhv, &
3234 dn, dnw, rdz, rdzw, fnm, fnp, &
3235 .false., &
3236 ids, ide, jds, jde, kds, kde, &
3237 ims, ime, jms, jme, kms, kme, &
3238 its, ite, jts, jte, kts, kte )
3239 ENDDO chem_loop
3240
3241 ENDIF
3242
3243
3244 IF (n_scalar .ge. PARAM_FIRST_SCALAR) THEN
3245
3246 scalar_loop: do im = PARAM_FIRST_SCALAR, n_scalar
3247
3248 CALL vertical_diffusion_s( scalar_tendf(ims,kms,jms,im), &
3249 config_flags, scalar(ims,kms,jms,im), &
3250 mu, xkhv, &
3251 dn, dnw, rdz, rdzw, fnm, fnp, &
3252 .false., &
3253 ids, ide, jds, jde, kds, kde, &
3254 ims, ime, jms, jme, kms, kme, &
3255 its, ite, jts, jte, kts, kte )
3256 ENDDO scalar_loop
3257
3258 ENDIF
3259
3260 END SUBROUTINE vertical_diffusion_2
3261
3262 !=======================================================================
3263 !=======================================================================
3264
3265 SUBROUTINE vertical_diffusion_u_2( tendency, config_flags, mu, &
3266 defor13, xkmv, &
3267 dnw, rdzw, fnm, fnp, &
3268 ids, ide, jds, jde, kds, kde, &
3269 ims, ime, jms, jme, kms, kme, &
3270 its, ite, jts, jte, kts, kte )
3271
3272 !-----------------------------------------------------------------------
3273 ! Begin declarations.
3274
3275 IMPLICIT NONE
3276
3277 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3278
3279 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3280 ims, ime, jms, jme, kms, kme, &
3281 its, ite, jts, jte, kts, kte
3282
3283 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
3284 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
3285 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
3286 ! REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: zeta_z
3287
3288 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
3289
3290 REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
3291 INTENT(IN ) ::defor13, &
3292 xkmv, &
3293 rdzw
3294 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
3295
3296 ! LOCAL VARS
3297
3298 INTEGER :: i, j, k, ktf
3299
3300 INTEGER :: i_start, i_end, j_start, j_end
3301 INTEGER :: is_ext,ie_ext,js_ext,je_ext
3302
3303 REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau3
3304
3305 REAL , DIMENSION( its:ite, jts:jte) :: zzavg
3306
3307 REAL :: rdzu
3308
3309 ! End declarations.
3310 !-----------------------------------------------------------------------
3311
3312 ktf=MIN(kte,kde-1)
3313
3314 i_start = its
3315 i_end = ite
3316 j_start = jts
3317 j_end = MIN(jte,jde-1)
3318
3319 IF ( config_flags%open_xs .or. config_flags%specified .or. &
3320 config_flags%nested) i_start = MAX(ids+1,its)
3321 IF ( config_flags%open_xe .or. config_flags%specified .or. &
3322 config_flags%nested) i_end = MIN(ide-1,ite)
3323 IF ( config_flags%open_ys .or. config_flags%specified .or. &
3324 config_flags%nested) j_start = MAX(jds+1,jts)
3325 IF ( config_flags%open_ye .or. config_flags%specified .or. &
3326 config_flags%nested) j_end = MIN(jde-2,jte)
3327 IF ( config_flags%periodic_x ) i_start = its
3328 IF ( config_flags%periodic_x ) i_end = ite
3329
3330 ! titau3 = titau13
3331 is_ext=0
3332 ie_ext=0
3333 js_ext=0
3334 je_ext=0
3335 CALL cal_titau_13_31( config_flags, titau3, defor13, &
3336 mu, xkmv, fnm, fnp, &
3337 is_ext, ie_ext, js_ext, je_ext, &
3338 ids, ide, jds, jde, kds, kde, &
3339 ims, ime, jms, jme, kms, kme, &
3340 its, ite, jts, jte, kts, kte )
3341 !
3342 DO j = j_start, j_end
3343 DO k=kts+1,ktf
3344 DO i = i_start, i_end
3345
3346 rdzu = 2./(1./rdzw(i,k,j) + 1./rdzw(i-1,k,j))
3347 tendency(i,k,j)=tendency(i,k,j)-rdzu*(titau3(i,k+1,j)-titau3(i,k,j))
3348
3349 ENDDO
3350 ENDDO
3351 ENDDO
3352
3353 ! ******** MODIF...
3354 ! we will pick up the surface drag (titau3(i,kts,j)) later
3355 !
3356 DO j = j_start, j_end
3357 k=kts
3358 DO i = i_start, i_end
3359
3360 rdzu = 2./(1./rdzw(i,k,j) + 1./rdzw(i-1,k,j))
3361 tendency(i,k,j)=tendency(i,k,j)-rdzu*(titau3(i,k+1,j))
3362 ENDDO
3363 ENDDO
3364 ! ******** MODIF...
3365
3366 END SUBROUTINE vertical_diffusion_u_2
3367
3368 !=======================================================================
3369 !=======================================================================
3370
3371 SUBROUTINE vertical_diffusion_v_2( tendency, config_flags, mu, &
3372 defor23, xkmv, &
3373 dnw, rdzw, fnm, fnp, &
3374 ids, ide, jds, jde, kds, kde, &
3375 ims, ime, jms, jme, kms, kme, &
3376 its, ite, jts, jte, kts, kte )
3377 !-----------------------------------------------------------------------
3378 ! Begin declarations.
3379
3380 IMPLICIT NONE
3381
3382 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3383
3384 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3385 ims, ime, jms, jme, kms, kme, &
3386 its, ite, jts, jte, kts, kte
3387 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
3388 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
3389 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
3390 ! REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: zeta_z
3391
3392 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
3393
3394 REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
3395 INTENT(IN ) ::defor23, &
3396 xkmv, &
3397 rdzw
3398
3399 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
3400
3401 ! LOCAL VARS
3402
3403 INTEGER :: i, j, k, ktf
3404
3405 INTEGER :: i_start, i_end, j_start, j_end
3406 INTEGER :: is_ext,ie_ext,js_ext,je_ext
3407
3408 REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau3
3409
3410 REAL , DIMENSION( its:ite, jts:jte) :: zzavg
3411
3412 REAL :: rdzv
3413
3414 ! End declarations.
3415 !-----------------------------------------------------------------------
3416
3417 ktf=MIN(kte,kde-1)
3418
3419 i_start = its
3420 i_end = MIN(ite,ide-1)
3421 j_start = jts
3422 j_end = jte
3423
3424 IF ( config_flags%open_xs .or. config_flags%specified .or. &
3425 config_flags%nested) i_start = MAX(ids+1,its)
3426 IF ( config_flags%open_xe .or. config_flags%specified .or. &
3427 config_flags%nested) i_end = MIN(ide-2,ite)
3428 IF ( config_flags%open_ys .or. config_flags%specified .or. &
3429 config_flags%nested) j_start = MAX(jds+1,jts)
3430 IF ( config_flags%open_ye .or. config_flags%specified .or. &
3431 config_flags%nested) j_end = MIN(jde-1,jte)
3432 IF ( config_flags%periodic_x ) i_start = its
3433 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
3434
3435 ! titau3 = titau23
3436 is_ext=0
3437 ie_ext=0
3438 js_ext=0
3439 je_ext=0
3440 CALL cal_titau_23_32( config_flags, titau3, defor23, &
3441 mu, xkmv, fnm, fnp, &
3442 is_ext, ie_ext, js_ext, je_ext, &
3443 ids, ide, jds, jde, kds, kde, &
3444 ims, ime, jms, jme, kms, kme, &
3445 its, ite, jts, jte, kts, kte )
3446
3447 DO j = j_start, j_end
3448 DO k = kts+1,ktf
3449 DO i = i_start, i_end
3450
3451 rdzv = 2./(1./rdzw(i,k,j) + 1./rdzw(i,k,j-1))
3452 tendency(i,k,j)=tendency(i,k,j)-rdzv*(titau3(i,k+1,j)-titau3(i,k,j))
3453
3454 ENDDO
3455 ENDDO
3456 ENDDO
3457
3458 ! ******** MODIF...
3459 ! we will pick up the surface drag (titau3(i,kts,j)) later
3460 !
3461 DO j = j_start, j_end
3462 k=kts
3463 DO i = i_start, i_end
3464
3465 rdzv = 2./(1./rdzw(i,k,j) + 1./rdzw(i,k,j-1))
3466 tendency(i,k,j)=tendency(i,k,j)-rdzv*(titau3(i,k+1,j))
3467
3468 ENDDO
3469 ENDDO
3470 ! ******** MODIF...
3471
3472 END SUBROUTINE vertical_diffusion_v_2
3473
3474 !=======================================================================
3475 !=======================================================================
3476
3477 SUBROUTINE vertical_diffusion_w_2(tendency, config_flags, mu, &
3478 defor33, tke, div, xkmv, &
3479 dn, rdz, &
3480 ids, ide, jds, jde, kds, kde, &
3481 ims, ime, jms, jme, kms, kme, &
3482 its, ite, jts, jte, kts, kte )
3483
3484 !-----------------------------------------------------------------------
3485 ! Begin declarations.
3486
3487 IMPLICIT NONE
3488
3489 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3490
3491 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3492 ims, ime, jms, jme, kms, kme, &
3493 its, ite, jts, jte, kts, kte
3494
3495 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
3496
3497 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
3498
3499 REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
3500 INTENT(IN ) ::defor33, &
3501 tke, &
3502 div, &
3503 xkmv, &
3504 rdz
3505
3506 REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: mu
3507
3508 ! LOCAL VARS
3509
3510 INTEGER :: i, j, k, ktf
3511
3512 INTEGER :: i_start, i_end, j_start, j_end
3513 INTEGER :: is_ext,ie_ext,js_ext,je_ext
3514
3515 REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau3
3516
3517 ! End declarations.
3518 !-----------------------------------------------------------------------
3519
3520 ktf=MIN(kte,kde-1)
3521
3522 i_start = its
3523 i_end = MIN(ite,ide-1)
3524 j_start = jts
3525 j_end = MIN(jte,jde-1)
3526
3527 IF ( config_flags%open_xs .or. config_flags%specified .or. &
3528 config_flags%nested) i_start = MAX(ids+1,its)
3529 IF ( config_flags%open_xe .or. config_flags%specified .or. &
3530 config_flags%nested) i_end = MIN(ide-2,ite)
3531 IF ( config_flags%open_ys .or. config_flags%specified .or. &
3532 config_flags%nested) j_start = MAX(jds+1,jts)
3533 IF ( config_flags%open_ye .or. config_flags%specified .or. &
3534 config_flags%nested) j_end = MIN(jde-2,jte)
3535 IF ( config_flags%periodic_x ) i_start = its
3536 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
3537
3538 ! titau3 = titau33
3539 is_ext=0
3540 ie_ext=0
3541 js_ext=0
3542 je_ext=0
3543 CALL cal_titau_11_22_33( config_flags, titau3, &
3544 mu, tke, xkmv, defor33, &
3545 is_ext, ie_ext, js_ext, je_ext, &
3546 ids, ide, jds, jde, kds, kde, &
3547 ims, ime, jms, jme, kms, kme, &
3548 its, ite, jts, jte, kts, kte )
3549
3550 ! DO j = j_start, j_end
3551 ! DO k = kts+1, ktf
3552 ! DO i = i_start, i_end
3553 ! titau3(i,k,j)=titau3(i,k,j)*zeta_z(i,j)
3554 ! ENDDO
3555 ! ENDDO
3556 ! ENDDO
3557
3558 DO j = j_start, j_end
3559 DO k = kts+1, ktf
3560 DO i = i_start, i_end
3561 tendency(i,k,j)=tendency(i,k,j)-rdz(i,k,j)*(titau3(i,k,j)-titau3(i,k-1,j))
3562 ENDDO
3563 ENDDO
3564 ENDDO
3565
3566 END SUBROUTINE vertical_diffusion_w_2
3567
3568 !=======================================================================
3569 !=======================================================================
3570
3571 SUBROUTINE vertical_diffusion_s( tendency, config_flags, var, mu, xkhv, &
3572 dn, dnw, rdz, rdzw, fnm, fnp, &
3573 doing_tke, &
3574 ids, ide, jds, jde, kds, kde, &
3575 ims, ime, jms, jme, kms, kme, &
3576 its, ite, jts, jte, kts, kte )
3577
3578 !-----------------------------------------------------------------------
3579 ! Begin declarations.
3580
3581 IMPLICIT NONE
3582
3583 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3584
3585 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3586 ims, ime, jms, jme, kms, kme, &
3587 its, ite, jts, jte, kts, kte
3588
3589 LOGICAL, INTENT(IN ) :: doing_tke
3590
3591 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
3592 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
3593 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
3594 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
3595
3596 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
3597
3598 REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(IN) :: xkhv
3599
3600 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: mu
3601
3602 REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
3603 INTENT(IN ) :: var, &
3604 rdz, &
3605 rdzw
3606 ! LOCAL VARS
3607
3608 INTEGER :: i, j, k, ktf
3609
3610 INTEGER :: i_start, i_end, j_start, j_end
3611
3612 REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: H3, &
3613 xkxavg, &
3614 rravg
3615
3616 REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: tmptendf
3617
3618 ! End declarations.
3619 !-----------------------------------------------------------------------
3620
3621 ktf=MIN(kte,kde-1)
3622
3623 i_start = its
3624 i_end = MIN(ite,ide-1)
3625 j_start = jts
3626 j_end = MIN(jte,jde-1)
3627
3628 IF ( config_flags%open_xs .or. config_flags%specified .or. &
3629 config_flags%nested) i_start = MAX(ids+1,its)
3630 IF ( config_flags%open_xe .or. config_flags%specified .or. &
3631 config_flags%nested) i_end = MIN(ide-2,ite)
3632 IF ( config_flags%open_ys .or. config_flags%specified .or. &
3633 config_flags%nested) j_start = MAX(jds+1,jts)
3634 IF ( config_flags%open_ye .or. config_flags%specified .or. &
3635 config_flags%nested) j_end = MIN(jde-2,jte)
3636 IF ( config_flags%periodic_x ) i_start = its
3637 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
3638
3639 IF (doing_tke) THEN
3640 DO j = j_start, j_end
3641 DO k = kts,ktf
3642 DO i = i_start, i_end
3643 tmptendf(i,k,j)=tendency(i,k,j)
3644 ENDDO
3645 ENDDO
3646 ENDDO
3647 ENDIF
3648
3649 ! H3
3650
3651 xkxavg = 0.
3652
3653 DO j = j_start, j_end
3654 DO k = kts+1,ktf
3655 DO i = i_start, i_end
3656 xkxavg(i,k,j)=fnm(k)*xkhv(i,k,j)+fnp(k)*xkhv(i,k-1,j)
3657 H3(i,k,j)=-xkxavg(i,k,j)*(var(i,k,j)-var(i,k-1,j))*rdz(i,k,j)
3658 ! H3(i,k,j)=-xkxavg(i,k,j)*zeta_z(i,j)* &
3659 ! (var(i,k,j)-var(i,k-1,j))/dn(k)
3660 ENDDO
3661 ENDDO
3662 ENDDO
3663
3664 DO j = j_start, j_end
3665 DO i = i_start, i_end
3666 H3(i,kts,j)=0.
3667 H3(i,ktf+1,j)=0.
3668 ! H3(i,kts,j)=H3(i,kts+1,j)
3669 ! H3(i,ktf+1,j)=H3(i,ktf,j)
3670 ENDDO
3671 ENDDO
3672
3673 DO j = j_start, j_end
3674 DO k = kts,ktf
3675 DO i = i_start, i_end
3676 tendency(i,k,j)=tendency(i,k,j) &
3677 -mu(i,j)*(H3(i,k+1,j)-H3(i,k,j))*rdzw(i,k,j)
3678 ENDDO
3679 ENDDO
3680 ENDDO
3681
3682 IF (doing_tke) THEN
3683 DO j = j_start, j_end
3684 DO k = kts,ktf
3685 DO i = i_start, i_end
3686 tendency(i,k,j)=tmptendf(i,k,j)+2.* &
3687 (tendency(i,k,j)-tmptendf(i,k,j))
3688 ENDDO
3689 ENDDO
3690 ENDDO
3691 ENDIF
3692
3693 END SUBROUTINE vertical_diffusion_s
3694
3695 !=======================================================================
3696 !=======================================================================
3697
3698 SUBROUTINE cal_titau_11_22_33( config_flags, titau, &
3699 mu, tke, xkx, defor, &
3700 is_ext, ie_ext, js_ext, je_ext, &
3701 ids, ide, jds, jde, kds, kde, &
3702 ims, ime, jms, jme, kms, kme, &
3703 its, ite, jts, jte, kts, kte )
3704
3705 ! History: Sep 2003 Changes by George Bryan and Jason Knievel, NCAR
3706 ! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
3707 ! Aug 2000 Original code by Shu-Hua Chen, UC-Davis
3708
3709 ! Purpose: This routine calculates stress terms (taus) for use in
3710 ! the calculation of production of TKE by sheared wind
3711
3712 ! References: Klemp and Wilhelmson (JAS 1978)
3713 ! Deardorff (B-L Meteor 1980)
3714 ! Chen and Dudhia (NCAR WRF physics report 2000)
3715
3716 ! Key:
3717
3718 !-----------------------------------------------------------------------
3719 ! Begin declarations.
3720
3721 IMPLICIT NONE
3722
3723 TYPE( grid_config_rec_type ), INTENT( IN ) &
3724 :: config_flags
3725
3726 INTEGER, INTENT( IN ) &
3727 :: ids, ide, jds, jde, kds, kde, &
3728 ims, ime, jms, jme, kms, kme, &
3729 its, ite, jts, jte, kts, kte
3730
3731 INTEGER, INTENT( IN ) &
3732 :: is_ext, ie_ext, js_ext, je_ext
3733
3734 REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ), INTENT( INOUT ) &
3735 :: titau
3736
3737 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
3738 :: defor, xkx, tke
3739
3740 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
3741 :: mu
3742
3743 ! Local variables.
3744
3745 INTEGER &
3746 :: i, j, k, ktf, i_start, i_end, j_start, j_end
3747
3748 ! End declarations.
3749 !-----------------------------------------------------------------------
3750
3751 ktf = MIN( kte, kde-1 )
3752
3753 i_start = its
3754 i_end = ite
3755 j_start = jts
3756 j_end = jte
3757
3758 IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
3759 config_flags%nested) i_start = MAX( ids+1, its )
3760 IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
3761 config_flags%nested) i_end = MIN( ide-1, ite )
3762 IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
3763 config_flags%nested) j_start = MAX( jds+1, jts )
3764 IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
3765 config_flags%nested) j_end = MIN( jde-1, jte )
3766 IF ( config_flags%periodic_x ) i_start = its
3767 IF ( config_flags%periodic_x ) i_end = ite
3768
3769 i_start = i_start - is_ext
3770 i_end = i_end + ie_ext
3771 j_start = j_start - js_ext
3772 j_end = j_end + je_ext
3773
3774 IF ( config_flags%km_opt .EQ. 2) THEN
3775 DO j = j_start,j_end
3776 DO k = kts,ktf
3777 DO i = i_start,i_end
3778 titau(i,k,j) = mu(i,j) * ( - xkx(i,k,j) * ( defor(i,k,j) ) )
3779 END DO
3780 END DO
3781 END DO
3782 ELSE
3783 DO j = j_start, j_end
3784 DO k = kts, ktf
3785 DO i = i_start, i_end
3786 titau(i,k,j) = - mu(i,j) * xkx(i,k,j) * defor(i,k,j)
3787 END DO
3788 END DO
3789 END DO
3790 END IF
3791
3792 END SUBROUTINE cal_titau_11_22_33
3793
3794 !=======================================================================
3795 !=======================================================================
3796
3797 SUBROUTINE cal_titau_12_21( config_flags, titau, &
3798 mu, xkx, defor, &
3799 is_ext, ie_ext, js_ext, je_ext, &
3800 ids, ide, jds, jde, kds, kde, &
3801 ims, ime, jms, jme, kms, kme, &
3802 its, ite, jts, jte, kts, kte )
3803
3804 ! History: Sep 2003 Modifications by George Bryan and Jason Knievel, NCAR
3805 ! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
3806 ! Aug 2000 Original code by Shu-Hua Chen, UC-Davis
3807
3808 ! Pusrpose This routine calculates the stress terms (taus) for use in
3809 ! the calculation of production of TKE by sheared wind
3810
3811 ! References: Klemp and Wilhelmson (JAS 1978)
3812 ! Deardorff (B-L Meteor 1980)
3813 ! Chen and Dudhia (NCAR WRF physics report 2000)
3814
3815 ! Key:
3816
3817 !-----------------------------------------------------------------------
3818 ! Begin declarations.
3819
3820 IMPLICIT NONE
3821
3822 TYPE( grid_config_rec_type), INTENT( IN ) &
3823 :: config_flags
3824
3825 INTEGER, INTENT( IN ) &
3826 :: ids, ide, jds, jde, kds, kde, &
3827 ims, ime, jms, jme, kms, kme, &
3828 its, ite, jts, jte, kts, kte
3829
3830 INTEGER, INTENT( IN ) &
3831 :: is_ext, ie_ext, js_ext, je_ext
3832
3833 REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ), INTENT( INOUT ) &
3834 :: titau
3835
3836 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
3837 :: defor, xkx
3838
3839 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
3840 :: mu
3841
3842 ! Local variables.
3843
3844 INTEGER &
3845 :: i, j, k, ktf, i_start, i_end, j_start, j_end
3846
3847 REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ) &
3848 :: xkxavg
3849
3850 REAL, DIMENSION( its-1:ite+1, jts-1:jte+1 ) &
3851 :: muavg
3852
3853 ! End declarations.
3854 !-----------------------------------------------------------------------
3855
3856 ktf = MIN( kte, kde-1 )
3857
3858 ! Needs one more point in the x and y directions.
3859
3860 i_start = its
3861 i_end = ite
3862 j_start = jts
3863 j_end = jte
3864
3865 IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
3866 config_flags%nested ) i_start = MAX( ids+1, its )
3867 IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
3868 config_flags%nested ) i_end = MIN( ide-1, ite )
3869 IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
3870 config_flags%nested ) j_start = MAX( jds+1, jts )
3871 IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
3872 config_flags%nested ) j_end = MIN( jde-1, jte )
3873 IF ( config_flags%periodic_x ) i_start = its
3874 IF ( config_flags%periodic_x ) i_end = ite
3875
3876 i_start = i_start - is_ext
3877 i_end = i_end + ie_ext
3878 j_start = j_start - js_ext
3879 j_end = j_end + je_ext
3880
3881 DO j = j_start, j_end
3882 DO k = kts, ktf
3883 DO i = i_start, i_end
3884 xkxavg(i,k,j) = 0.25 * ( xkx(i-1,k,j ) + xkx(i,k,j ) + &
3885 xkx(i-1,k,j-1) + xkx(i,k,j-1) )
3886 END DO
3887 END DO
3888 END DO
3889
3890 DO j = j_start, j_end
3891 DO i = i_start, i_end
3892 muavg(i,j) = 0.25 * ( mu(i-1,j ) + mu(i,j ) + &
3893 mu(i-1,j-1) + mu(i,j-1) )
3894 END DO
3895 END DO
3896
3897 ! titau12 or titau21
3898
3899 DO j = j_start, j_end
3900 DO k = kts, ktf
3901 DO i = i_start, i_end
3902 titau(i,k,j) = - muavg(i,j) * xkxavg(i,k,j) * defor(i,k,j)
3903 END DO
3904 END DO
3905 END DO
3906
3907 END SUBROUTINE cal_titau_12_21
3908
3909 !=======================================================================
3910
3911 SUBROUTINE cal_titau_13_31( config_flags, titau, &
3912 defor, mu, xkx, fnm, fnp, &
3913 is_ext, ie_ext, js_ext, je_ext, &
3914 ids, ide, jds, jde, kds, kde, &
3915 ims, ime, jms, jme, kms, kme, &
3916 its, ite, jts, jte, kts, kte )
3917
3918 ! History: Sep 2003 Modifications by George Bryan and Jason Knievel, NCAR
3919 ! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
3920 ! Aug 2000 Original code by Shu-Hua Chen, UC-Davis
3921
3922 ! Purpose: This routine calculates the stress terms (taus) for use in
3923 ! the calculation of production of TKE by sheared wind
3924
3925 ! References: Klemp and Wilhelmson (JAS 1978)
3926 ! Deardorff (B-L Meteor 1980)
3927 ! Chen and Dudhia (NCAR WRF physics report 2000)
3928
3929 ! Key:
3930
3931 !-----------------------------------------------------------------------
3932 ! Begin declarations.
3933
3934 IMPLICIT NONE
3935
3936 TYPE( grid_config_rec_type), INTENT( IN ) &
3937 :: config_flags
3938
3939 INTEGER, INTENT( IN ) &
3940 :: ids, ide, jds, jde, kds, kde, &
3941 ims, ime, jms, jme, kms, kme, &
3942 its, ite, jts, jte, kts, kte
3943
3944 INTEGER, INTENT( IN ) &
3945 :: is_ext, ie_ext, js_ext, je_ext
3946
3947 REAL, DIMENSION( kms:kme ), INTENT( IN ) &
3948 :: fnm, fnp
3949
3950 REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ), INTENT( INOUT ) &
3951 :: titau
3952
3953 REAL, DIMENSION( ims:ime, kms:kme, jms:jme), INTENT( IN ) &
3954 :: defor, xkx
3955
3956 REAL, DIMENSION( ims:ime, jms:jme), INTENT( IN ) &
3957 :: mu
3958
3959 ! Local variables.
3960
3961 INTEGER &
3962 :: i, j, k, ktf, i_start, i_end, j_start, j_end
3963
3964 REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ) &
3965 :: xkxavg
3966
3967 REAL, DIMENSION( its-1:ite+1, jts-1:jte+1 ) &
3968 :: muavg
3969
3970 ! End declarations.
3971 !-----------------------------------------------------------------------
3972
3973 ktf = MIN( kte, kde-1 )
3974
3975 ! Find ide-1 and jde-1 for averaging to p point.
3976
3977 i_start = its
3978 i_end = ite
3979 j_start = jts
3980 j_end = MIN( jte, jde-1 )
3981
3982 IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
3983 config_flags%nested) i_start = MAX( ids+1, its )
3984 IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
3985 config_flags%nested) i_end = MIN( ide-1, ite )
3986 IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
3987 config_flags%nested) j_start = MAX( jds+1, jts )
3988 IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
3989 config_flags%nested) j_end = MIN( jde-2, jte )
3990 IF ( config_flags%periodic_x ) i_start = its
3991 IF ( config_flags%periodic_x ) i_end = ite
3992
3993 i_start = i_start - is_ext
3994 i_end = i_end + ie_ext
3995 j_start = j_start - js_ext
3996 j_end = j_end + je_ext
3997
3998 DO j = j_start, j_end
3999 DO k = kts+1, ktf
4000 DO i = i_start, i_end
4001 xkxavg(i,k,j) = 0.5 * ( fnm(k) * ( xkx(i,k ,j) + xkx(i-1,k ,j) ) + &
4002 fnp(k) * ( xkx(i,k-1,j) + xkx(i-1,k-1,j) ) )
4003 END DO
4004 END DO
4005 END DO
4006
4007 DO j = j_start, j_end
4008 DO i = i_start, i_end
4009 muavg(i,j) = 0.5 * ( mu(i,j) + mu(i-1,j) )
4010 END DO
4011 END DO
4012
4013 DO j = j_start, j_end
4014 DO k = kts+1, ktf
4015 DO i = i_start, i_end
4016 titau(i,k,j) = - muavg(i,j) * xkxavg(i,k,j) * defor(i,k,j)
4017 ENDDO
4018 ENDDO
4019 ENDDO
4020
4021 DO j = j_start, j_end
4022 DO i = i_start, i_end
4023 titau(i,kts ,j) = 0.0
4024 titau(i,ktf+1,j) = 0.0
4025 ENDDO
4026 ENDDO
4027
4028 END SUBROUTINE cal_titau_13_31
4029
4030 !=======================================================================
4031 !=======================================================================
4032
4033 SUBROUTINE cal_titau_23_32( config_flags, titau, defor, &
4034 mu, xkx, fnm, fnp, &
4035 is_ext, ie_ext, js_ext, je_ext, &
4036 ids, ide, jds, jde, kds, kde, &
4037 ims, ime, jms, jme, kms, kme, &
4038 its, ite, jts, jte, kts, kte )
4039
4040 ! History: Sep 2003 Changes by George Bryan and Jason Knievel, NCAR
4041 ! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
4042 ! Aug 2000 Original code by Shu-Hua Chen, UC-Davis
4043
4044 ! Purpose: This routine calculates stress terms (taus) for use in
4045 ! the calculation of production of TKE by sheared wind
4046
4047 ! References: Klemp and Wilhelmson (JAS 1978)
4048 ! Deardorff (B-L Meteor 1980)
4049 ! Chen and Dudhia (NCAR WRF physics report 2000)
4050
4051 ! Key:
4052
4053 !-----------------------------------------------------------------------
4054 ! Begin declarations.
4055
4056 IMPLICIT NONE
4057
4058 TYPE( grid_config_rec_type ), INTENT( IN ) &
4059 :: config_flags
4060
4061 INTEGER, INTENT( IN ) &
4062 :: ids, ide, jds, jde, kds, kde, &
4063 ims, ime, jms, jme, kms, kme, &
4064 its, ite, jts, jte, kts, kte
4065
4066 INTEGER, INTENT( IN ) &
4067 :: is_ext,ie_ext,js_ext,je_ext
4068
4069 REAL, DIMENSION( kms:kme ), INTENT( IN ) &
4070 :: fnm, fnp
4071
4072 REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ), INTENT( INOUT ) &
4073 :: titau
4074
4075 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
4076 :: defor, xkx
4077
4078 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
4079 :: mu
4080
4081 ! Local variables.
4082
4083 INTEGER &
4084 :: i, j, k, ktf, i_start, i_end, j_start, j_end
4085
4086 REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ) &
4087 :: xkxavg
4088
4089 REAL, DIMENSION( its-1:ite+1, jts-1:jte+1 ) &
4090 :: muavg
4091
4092 ! End declarations.
4093 !-----------------------------------------------------------------------
4094
4095 ktf = MIN( kte, kde-1 )
4096
4097 ! Find ide-1 and jde-1 for averaging to p point.
4098
4099 i_start = its
4100 i_end = MIN( ite, ide-1 )
4101 j_start = jts
4102 j_end = jte
4103
4104 IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
4105 config_flags%nested) i_start = MAX( ids+1, its )
4106 IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
4107 config_flags%nested) i_end = MIN( ide-2, ite )
4108 IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
4109 config_flags%nested) j_start = MAX( jds+1, jts )
4110 IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
4111 config_flags%nested) j_end = MIN( jde-1, jte )
4112 IF ( config_flags%periodic_x ) i_start = its
4113 IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
4114
4115 i_start = i_start - is_ext
4116 i_end = i_end + ie_ext
4117 j_start = j_start - js_ext
4118 j_end = j_end + je_ext
4119
4120 DO j = j_start, j_end
4121 DO k = kts+1, ktf
4122 DO i = i_start, i_end
4123 xkxavg(i,k,j) = 0.5 * ( fnm(k) * ( xkx(i,k ,j) + xkx(i,k ,j-1) ) + &
4124 fnp(k) * ( xkx(i,k-1,j) + xkx(i,k-1,j-1) ) )
4125 END DO
4126 END DO
4127 END DO
4128
4129 DO j = j_start, j_end
4130 DO i = i_start, i_end
4131 muavg(i,j) = 0.5 * ( mu(i,j) + mu(i,j-1) )
4132 END DO
4133 END DO
4134
4135 DO j = j_start, j_end
4136 DO k = kts+1, ktf
4137 DO i = i_start, i_end
4138 titau(i,k,j) = - muavg(i,j) * xkxavg(i,k,j) * defor(i,k,j)
4139 END DO
4140 END DO
4141 END DO
4142
4143 DO j = j_start, j_end
4144 DO i = i_start, i_end
4145 titau(i,kts ,j) = 0.0
4146 titau(i,ktf+1,j) = 0.0
4147 END DO
4148 END DO
4149
4150 END SUBROUTINE cal_titau_23_32
4151
4152 !=======================================================================
4153 !=======================================================================
4154
4155 SUBROUTINE phy_bc ( config_flags,div,defor11,defor22,defor33, &
4156 defor12,defor13,defor23,xkmh,xkmhd,xkmv,xkhh,xkhv,tke, &
4157 RUBLTEN, RVBLTEN, &
4158 ids, ide, jds, jde, kds, kde, &
4159 ims, ime, jms, jme, kms, kme, &
4160 ips, ipe, jps, jpe, kps, kpe, &
4161 its, ite, jts, jte, kts, kte )
4162
4163 !------------------------------------------------------------------------------
4164 ! Begin declarations.
4165
4166 IMPLICIT NONE
4167
4168 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
4169
4170 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4171 ims, ime, jms, jme, kms, kme, &
4172 ips, ipe, jps, jpe, kps, kpe, &
4173 its, ite, jts, jte, kts, kte
4174
4175 REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::RUBLTEN, &
4176 RVBLTEN, &
4177 defor11, &
4178 defor22, &
4179 defor33, &
4180 defor12, &
4181 defor13, &
4182 defor23, &
4183 xkmh, &
4184 xkmhd, &
4185 xkmv, &
4186 xkhh, &
4187 xkhv, &
4188 tke, &
4189 div
4190
4191 ! End declarations.
4192 !-----------------------------------------------------------------------
4193
4194 IF(config_flags%bl_pbl_physics .GT. 0) THEN
4195
4196 CALL set_physical_bc3d( RUBLTEN , 't', config_flags, &
4197 ids, ide, jds, jde, kds, kde, &
4198 ims, ime, jms, jme, kms, kme, &
4199 ips, ipe, jps, jpe, kps, kpe, &
4200 its, ite, jts, jte, kts, kte )
4201
4202 CALL set_physical_bc3d( RVBLTEN , 't', config_flags, &
4203 ids, ide, jds, jde, kds, kde, &
4204 ims, ime, jms, jme, kms, kme, &
4205 ips, ipe, jps, jpe, kps, kpe, &
4206 its, ite, jts, jte, kts, kte )
4207
4208 ENDIF
4209
4210 ! move out of the conditional, below; this one is needed for
4211 ! all diff_opt cases. JM
4212 CALL set_physical_bc3d( xkmhd , 't', config_flags, &
4213 ids, ide, jds, jde, kds, kde, &
4214 ims, ime, jms, jme, kms, kme, &
4215 ips, ipe, jps, jpe, kps, kpe, &
4216 its, ite, jts, jte, kts, kte )
4217
4218 IF(config_flags%diff_opt .eq. 2) THEN
4219
4220 CALL set_physical_bc3d( xkmh , 't', config_flags, &
4221 ids, ide, jds, jde, kds, kde, &
4222 ims, ime, jms, jme, kms, kme, &
4223 ips, ipe, jps, jpe, kps, kpe, &
4224 its, ite, jts, jte, kts, kte )
4225
4226
4227 CALL set_physical_bc3d( xkmv , 't', config_flags, &
4228 ids, ide, jds, jde, kds, kde, &
4229 ims, ime, jms, jme, kms, kme, &
4230 ips, ipe, jps, jpe, kps, kpe, &
4231 its, ite, jts, jte, kts, kte )
4232
4233 CALL set_physical_bc3d( xkhh , 't', config_flags, &
4234 ids, ide, jds, jde, kds, kde, &
4235 ims, ime, jms, jme, kms, kme, &
4236 ips, ipe, jps, jpe, kps, kpe, &
4237 its, ite, jts, jte, kts, kte )
4238
4239 CALL set_physical_bc3d( xkhv , 't', config_flags, &
4240 ids, ide, jds, jde, kds, kde, &
4241 ims, ime, jms, jme, kms, kme, &
4242 ips, ipe, jps, jpe, kps, kpe, &
4243 its, ite, jts, jte, kts, kte )
4244
4245 CALL set_physical_bc3d( tke , 't', config_flags, &
4246 ids, ide, jds, jde, kds, kde, &
4247 ims, ime, jms, jme, kms, kme, &
4248 ips, ipe, jps, jpe, kps, kpe, &
4249 its, ite, jts, jte, kts, kte )
4250
4251 CALL set_physical_bc3d( div , 't', config_flags, &
4252 ids, ide, jds, jde, kds, kde, &
4253 ims, ime, jms, jme, kms, kme, &
4254 ips, ipe, jps, jpe, kps, kpe, &
4255 its, ite, jts, jte, kts, kte )
4256
4257 CALL set_physical_bc3d( defor11 , 't', config_flags, &
4258 ids, ide, jds, jde, kds, kde, &
4259 ims, ime, jms, jme, kms, kme, &
4260 ips, ipe, jps, jpe, kps, kpe, &
4261 its, ite, jts, jte, kts, kte )
4262
4263 CALL set_physical_bc3d( defor22 , 't', config_flags, &
4264 ids, ide, jds, jde, kds, kde, &
4265 ims, ime, jms, jme, kms, kme, &
4266 ips, ipe, jps, jpe, kps, kpe, &
4267 its, ite, jts, jte, kts, kte )
4268
4269 CALL set_physical_bc3d( defor33 , 't', config_flags, &
4270 ids, ide, jds, jde, kds, kde, &
4271 ims, ime, jms, jme, kms, kme, &
4272 ips, ipe, jps, jpe, kps, kpe, &
4273 its, ite, jts, jte, kts, kte )
4274
4275 CALL set_physical_bc3d( defor12 , 'd', config_flags, &
4276 ids, ide, jds, jde, kds, kde, &
4277 ims, ime, jms, jme, kms, kme, &
4278 ips, ipe, jps, jpe, kps, kpe, &
4279 its, ite, jts, jte, kts, kte )
4280
4281 CALL set_physical_bc3d( defor13 , 'e', config_flags, &
4282 ids, ide, jds, jde, kds, kde, &
4283 ims, ime, jms, jme, kms, kme, &
4284 ips, ipe, jps, jpe, kps, kpe, &
4285 its, ite, jts, jte, kts, kte )
4286
4287 CALL set_physical_bc3d( defor23 , 'f', config_flags, &
4288 ids, ide, jds, jde, kds, kde, &
4289 ims, ime, jms, jme, kms, kme, &
4290 ips, ipe, jps, jpe, kps, kpe, &
4291 its, ite, jts, jte, kts, kte )
4292
4293 ENDIF
4294
4295 END SUBROUTINE phy_bc
4296
4297 !=======================================================================
4298 !=======================================================================
4299
4300 SUBROUTINE tke_rhs( tendency, BN2, config_flags, &
4301 defor11, defor22, defor33, &
4302 defor12, defor13, defor23, &
4303 u, v, w, div, tke, mu, &
4304 theta, p, p8w, t8w, z, fnm, fnp, &
4305 cf1, cf2, cf3, msft, xkmh, xkmv, xkhv, &
4306 rdx, rdy, dx, dy, dt, zx, zy, &
4307 rdz, rdzw, dn, dnw, cr_len, &
4308 ids, ide, jds, jde, kds, kde, &
4309 ims, ime, jms, jme, kms, kme, &
4310 its, ite, jts, jte, kts, kte )
4311
4312 !-----------------------------------------------------------------------
4313 ! Begin declarations.
4314
4315 IMPLICIT NONE
4316
4317 TYPE( grid_config_rec_type ), INTENT( IN ) &
4318 :: config_flags
4319
4320 INTEGER, INTENT( IN ) &
4321 :: ids, ide, jds, jde, kds, kde, &
4322 ims, ime, jms, jme, kms, kme, &
4323 its, ite, jts, jte, kts, kte
4324
4325 REAL, INTENT( IN ) &
4326 :: cf1, cf2, cf3, dt, rdx, rdy, dx, dy, cr_len
4327
4328 REAL, DIMENSION( kms:kme ), INTENT( IN ) &
4329 :: fnm, fnp, dnw, dn
4330
4331 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
4332 :: msft
4333
4334 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
4335 :: tendency
4336
4337 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
4338 :: defor11, defor22, defor33, defor12, defor13, defor23, &
4339 div, BN2, tke, xkmh, xkmv, xkhv, zx, zy, u, v, w, theta, &
4340 p, p8w, t8w, z, rdz, rdzw
4341
4342 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
4343 :: mu
4344
4345 ! Local variables.
4346
4347 INTEGER &
4348 :: i, j, k, ktf, i_start, i_end, j_start, j_end
4349
4350 ! End declarations.
4351 !-----------------------------------------------------------------------
4352
4353 CALL tke_shear( tendency, config_flags, &
4354 defor11, defor22, defor33, &
4355 defor12, defor13, defor23, &
4356 u, v, w, tke, mu, fnm, fnp, &
4357 cf1, cf2, cf3, msft, xkmh, xkmv, &
4358 rdx, rdy, zx, zy, rdz, rdzw, dnw, dn, &
4359 ids, ide, jds, jde, kds, kde, &
4360 ims, ime, jms, jme, kms, kme, &
4361 its, ite, jts, jte, kts, kte )
4362
4363 CALL tke_buoyancy( tendency, config_flags, mu, &
4364 tke, xkhv, BN2, theta, dt, &
4365 ids, ide, jds, jde, kds, kde, &
4366 ims, ime, jms, jme, kms, kme, &
4367 its, ite, jts, jte, kts, kte )
4368
4369 CALL tke_dissip( tendency, config_flags, &
4370 mu, tke, bn2, theta, p8w, t8w, z, &
4371 dx, dy,rdz, rdzw, cr_len, &
4372 ids, ide, jds, jde, kds, kde, &
4373 ims, ime, jms, jme, kms, kme, &
4374 its, ite, jts, jte, kts, kte )
4375
4376 ! Set a lower limit on TKE.
4377
4378 ktf = MIN( kte, kde-1 )
4379 i_start = its
4380 i_end = MIN( ite, ide-1 )
4381 j_start = jts
4382 j_end = MIN( jte, jde-1 )
4383
4384 IF ( config_flags%open_xs .or. config_flags%specified .or. &
4385 config_flags%nested) i_start = MAX(ids+1,its)
4386 IF ( config_flags%open_xe .or. config_flags%specified .or. &
4387 config_flags%nested) i_end = MIN(ide-2,ite)
4388 IF ( config_flags%open_ys .or. config_flags%specified .or. &
4389 config_flags%nested) j_start = MAX(jds+1,jts)
4390 IF ( config_flags%open_ye .or. config_flags%specified .or. &
4391 config_flags%nested) j_end = MIN(jde-2,jte)
4392 IF ( config_flags%periodic_x ) i_start = its
4393 IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
4394
4395 DO j = j_start, j_end
4396 DO k = kts, ktf
4397 DO i = i_start, i_end
4398 tendency(i,k,j) = max( tendency(i,k,j), -mu(i,j) * max( 0.0 , tke(i,k,j) ) / dt )
4399 END DO
4400 END DO
4401 END DO
4402
4403 END SUBROUTINE tke_rhs
4404
4405 !=======================================================================
4406 !=======================================================================
4407
4408 SUBROUTINE tke_buoyancy( tendency, config_flags, mu, &
4409 tke, xkhv, BN2, theta, dt, &
4410 ids, ide, jds, jde, kds, kde, &
4411 ims, ime, jms, jme, kms, kme, &
4412 its, ite, jts, jte, kts, kte )
4413
4414 !-----------------------------------------------------------------------
4415 ! Begin declarations.
4416
4417 IMPLICIT NONE
4418
4419 TYPE( grid_config_rec_type ), INTENT( IN ) &
4420 :: config_flags
4421
4422 INTEGER, INTENT( IN ) &
4423 :: ids, ide, jds, jde, kds, kde, &
4424 ims, ime, jms, jme, kms, kme, &
4425 its, ite, jts, jte, kts, kte
4426
4427 REAL, INTENT( IN ) &
4428 :: dt
4429
4430 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
4431 :: tendency
4432
4433 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
4434 :: xkhv, tke, BN2, theta
4435
4436 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
4437 :: mu
4438
4439 ! Local variables.
4440
4441 INTEGER &
4442 :: i, j, k, ktf
4443
4444 INTEGER &
4445 :: i_start, i_end, j_start, j_end
4446
4447 REAL :: heat_flux
4448
4449 ! End declarations.
4450 !-----------------------------------------------------------------------
4451
4452 !-----------------------------------------------------------------------
4453 ! Add to the TKE tendency the term that accounts for production of TKE
4454 ! due to buoyant motions.
4455
4456 ktf = MIN( kte, kde-1 )
4457 i_start = its
4458 i_end = MIN( ite, ide-1 )
4459 j_start = jts
4460 j_end = MIN( jte, jde-1 )
4461
4462 IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
4463 config_flags%nested ) i_start = MAX( ids+1, its )
4464 IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
4465 config_flags%nested ) i_end = MIN( ide-2, ite )
4466 IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
4467 config_flags%nested ) j_start = MAX( jds+1, jts )
4468 IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
4469 config_flags%nested ) j_end = MIN( jde-2, jte )
4470 IF ( config_flags%periodic_x ) i_start = its
4471 IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
4472
4473 DO j = j_start, j_end
4474 DO k = kts+1, ktf
4475 DO i = i_start, i_end
4476 tendency(i,k,j) = tendency(i,k,j) - mu(i,j) * xkhv(i,k,j) * BN2(i,k,j)
4477 END DO
4478 END DO
4479 END DO
4480
4481 ! MARTA: change in the computation of the tke's tendency at the surface.
4482 ! the buoyancy flux is the average of the surface heat flux (0.06) and the
4483 ! flux at the first w level
4484 !
4485 ! WCS 040331
4486
4487 heat_flux = config_flags%tke_heat_flux ! constant heat flux value
4488 ! set in namelist.input
4489 IF(abs(heat_flux).lt.1.0e-6)THEN
4490 K=KTS
4491 DO j = j_start, j_end
4492 DO i = i_start, i_end
4493 tendency(i,k,j)= tendency(i,k,j) - &
4494 mu(i,j)*xkhv(i,k,j)*BN2(i,k,j)
4495 ENDDO
4496 ENDDO
4497 ELSE
4498 K=KTS
4499 DO j = j_start, j_end
4500 DO i = i_start, i_end
4501 tendency(i,k,j)= tendency(i,k,j) - &
4502 mu(i,j)*((xkhv(i,k,j)*BN2(i,k,j))- (g/theta(i,k,j))*heat_flux)/2.
4503
4504 ENDDO
4505 ENDDO
4506 ENDIF
4507 ! end of MARTA/WCS change
4508
4509 ! The tendency array now includes production of TKE from buoyant
4510 ! motions.
4511 !-----------------------------------------------------------------------
4512
4513 END SUBROUTINE tke_buoyancy
4514
4515 !=======================================================================
4516 !=======================================================================
4517
4518 SUBROUTINE tke_dissip( tendency, config_flags, &
4519 mu, tke, bn2, theta, p8w, t8w, z, &
4520 dx, dy, rdz, rdzw, cr_len_in, &
4521 ids, ide, jds, jde, kds, kde, &
4522 ims, ime, jms, jme, kms, kme, &
4523 its, ite, jts, jte, kts, kte )
4524
4525 ! History: Sep 2003 Changes by George Bryan and Jason Knievel, NCAR
4526 ! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
4527 ! Aug 2000 Original code by Shu-Hua Chen, UC-Davis
4528
4529 ! Purpose: This routine calculates dissipation of turbulent kinetic
4530 ! energy.
4531
4532 ! References: Klemp and Wilhelmson (JAS 1978)
4533 ! Deardorff (B-L Meteor 1980)
4534 ! Chen and Dudhia (NCAR WRF physics report 2000)
4535
4536 !-----------------------------------------------------------------------
4537 ! Begin declarations.
4538
4539 IMPLICIT NONE
4540
4541 TYPE( grid_config_rec_type ), INTENT( IN ) &
4542 :: config_flags
4543
4544 INTEGER, INTENT( IN ) &
4545 :: ids, ide, jds, jde, kds, kde, &
4546 ims, ime, jms, jme, kms, kme, &
4547 its, ite, jts, jte, kts, kte
4548
4549 REAL, INTENT( IN ) &
4550 :: dx, dy, cr_len_in
4551
4552 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
4553 :: tendency
4554
4555 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
4556 :: tke, bn2, theta, p8w, t8w, z, rdz, rdzw
4557
4558 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
4559 :: mu
4560
4561 ! Local variables.
4562
4563 REAL, DIMENSION( its:ite, kts:kte, jts:jte ) &
4564 :: dthrdn
4565
4566 REAL, DIMENSION( its:ite, kts:kte, jts:jte ) &
4567 :: l_scale
4568
4569 REAL, DIMENSION( its:ite ) &
4570 :: sumtke, sumtkez
4571
4572 INTEGER &
4573 :: i, j, k, ktf, i_start, i_end, j_start, j_end
4574
4575 REAL &
4576 :: disp_len, deltas, coefc, tmpdz, len_s, thetasfc, &
4577 thetatop, len_0, tketmp, tmp, cr_len, ce1, ce2
4578
4579 ! End declarations.
4580 !-----------------------------------------------------------------------
4581
4582 ce1 = ( c_k / 0.10 ) * 0.19
4583 ce2 = max( 0.0 , 0.93 - ce1 )
4584
4585 ktf = MIN( kte, kde-1 )
4586 i_start = its
4587 i_end = MIN(ite,ide-1)
4588 j_start = jts
4589 j_end = MIN(jte,jde-1)
4590
4591 IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
4592 config_flags%nested) i_start = MAX( ids+1, its )
4593 IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
4594 config_flags%nested) i_end = MIN( ide-2, ite )
4595 IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
4596 config_flags%nested) j_start = MAX( jds+1, jts )
4597 IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
4598 config_flags%nested) j_end = MIN( jde-2, jte )
4599 IF ( config_flags%periodic_x ) i_start = its
4600 IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
4601
4602 cr_len = cr_len_in
4603 cr_len = dx + 1.0 ! hardwire for mixing length = (dx*dy*dz)**(1/3).
4604 ! remove this for the alternate formulation
4605
4606 IF (dx .gt. cr_len) THEN
4607
4608 DO j = j_start, j_end
4609 DO i = i_start, i_end
4610 sumtke(i) = 0.0
4611 sumtkez(i) = 0.0
4612 END DO
4613 DO k = kts, ktf
4614 DO i = i_start, i_end
4615 tketmp = MAX( tke(i,k,j), 0.0 )
4616 sumtke(i) = sumtke(i) + SQRT(tketmp) / rdzw(i,k,j)
4617 sumtkez(i) = sumtkez(i)+ sumtke(i) * z(i,k,j)
4618 IF ( ABS( sumtke(i) ) .gt. 0.01 ) THEN
4619 len_0 = 0.2 * sumtkez(i) / sumtke(i)
4620 ELSE
4621 len_0 = 80.0
4622 ENDIF
4623 len_0 = MIN( 80.0, len_0)
4624 l_scale(i,k,j) = KARMAN * z(i,k,j) / &
4625 ( 1.0 + KARMAN * z(i,k,j) / len_0 )
4626 tendency(i,k,j) = tendency(i,k,j) - &
4627 mu(i,j) * 2.0 * SQRT( 2.0 ) / 15.0 * &
4628 tketmp**1.5 / l_scale(i,k,j)
4629 END DO
4630 END DO
4631 END DO
4632 ELSE
4633 CALL calc_l_scale( config_flags, tke, BN2, l_scale, &
4634 i_start, i_end, ktf, j_start, j_end, &
4635 dx, dy, rdzw, &
4636 ids, ide, jds, jde, kds, kde, &
4637 ims, ime, jms, jme, kms, kme, &
4638 its, ite, jts, jte, kts, kte )
4639 DO j = j_start, j_end
4640 DO k = kts, ktf
4641 DO i = i_start, i_end
4642 deltas = ( dx * dy / rdzw(i,k,j) )**0.33333333
4643 tketmp = MAX( tke(i,k,j), 1.0e-6 )
4644
4645 ! Apply Deardorff's (1980) "wall effect" at the bottom of the domain.
4646
4647 IF ( k .eq. kts .or. k .eq. ktf ) then
4648 coefc = 3.9
4649 ELSE
4650 coefc = ce1 + ce2 * l_scale(i,k,j) / deltas
4651 END IF
4652
4653 tendency(i,k,j) = tendency(i,k,j) - &
4654 mu(i,j) * coefc * tketmp**1.5 / l_scale(i,k,j)
4655 END DO
4656 END DO
4657 END DO
4658 ENDIF
4659
4660 END SUBROUTINE tke_dissip
4661
4662 !=======================================================================
4663 !=======================================================================
4664
4665 SUBROUTINE tke_shear( tendency, config_flags, &
4666 defor11, defor22, defor33, &
4667 defor12, defor13, defor23, &
4668 u, v, w, tke, mu, fnm, fnp, &
4669 cf1, cf2, cf3, msft, xkmh, xkmv, &
4670 rdx, rdy, zx, zy, rdz, rdzw, dn, dnw, &
4671 ids, ide, jds, jde, kds, kde, &
4672 ims, ime, jms, jme, kms, kme, &
4673 its, ite, jts, jte, kts, kte )
4674
4675 ! History: Sep 2003 Rewritten by George Bryan and Jason Knievel,
4676 ! NCAR
4677 ! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
4678 ! Aug 2000 Original code by Shu-Hua Chen, UC-Davis
4679
4680 ! Purpose: This routine calculates the production of turbulent
4681 ! kinetic energy by stresses due to sheared wind.
4682
4683 ! References: Klemp and Wilhelmson (JAS 1978)
4684 ! Deardorff (B-L Meteor 1980)
4685 ! Chen and Dudhia (NCAR WRF physics report 2000)
4686
4687 ! Key:
4688
4689 ! avg temporary working array
4690 ! cf1
4691 ! cf2
4692 ! cf3
4693 ! defor11 deformation term ( du/dx + du/dx )
4694 ! defor12 deformation term ( dv/dx + du/dy ); same as defor21
4695 ! defor13 deformation term ( dw/dx + du/dz ); same as defor31
4696 ! defor22 deformation term ( dv/dy + dv/dy )
4697 ! defor23 deformation term ( dw/dy + dv/dz ); same as defor32
4698 ! defor33 deformation term ( dw/dz + dw/dz )
4699 ! div 3-d divergence
4700 ! dn
4701 ! dnw
4702 ! fnm
4703 ! fnp
4704 ! msft
4705 ! rdx
4706 ! rdy
4707 ! tendency
4708 ! titau tau (stress tensor) with a tilde, indicating division by
4709 ! a map-scale factor and the fraction of the total modeled
4710 ! atmosphere beneath a given altitude (titau = tau/m/zeta)
4711 ! tke turbulent kinetic energy
4712
4713 !-----------------------------------------------------------------------
4714 ! Begin declarations.
4715
4716 IMPLICIT NONE
4717
4718 TYPE( grid_config_rec_type ), INTENT( IN ) &
4719 :: config_flags
4720
4721 INTEGER, INTENT( IN ) &
4722 :: ids, ide, jds, jde, kds, kde, &
4723 ims, ime, jms, jme, kms, kme, &
4724 its, ite, jts, jte, kts, kte
4725
4726 REAL, INTENT( IN ) &
4727 :: cf1, cf2, cf3, rdx, rdy
4728
4729 REAL, DIMENSION( kms:kme ), INTENT( IN ) &
4730 :: fnm, fnp, dn, dnw
4731
4732 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
4733 :: msft
4734
4735 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
4736 :: tendency
4737
4738 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
4739 :: defor11, defor22, defor33, defor12, defor13, defor23, &
4740 tke, xkmh, xkmv, zx, zy, u, v, w, rdz, rdzw
4741
4742 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
4743 :: mu
4744
4745 ! Local variables.
4746
4747 INTEGER &
4748 :: i, j, k, ktf, ktes1, ktes2, &
4749 i_start, i_end, j_start, j_end, &
4750 is_ext, ie_ext, js_ext, je_ext
4751
4752 REAL &
4753 :: mtau
4754
4755 REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ) &
4756 :: avg, titau, tmp2
4757
4758 REAL, DIMENSION( its:ite, kts:kte, jts:jte ) &
4759 :: titau12, tmp1, zxavg, zyavg
4760
4761 REAL :: absU, cd0
4762
4763 ! End declarations.
4764 !-----------------------------------------------------------------------
4765
4766 ktf = MIN( kte, kde-1 )
4767 ktes1 = kte-1
4768 ktes2 = kte-2
4769
4770 i_start = its
4771 i_end = MIN( ite, ide-1 )
4772 j_start = jts
4773 j_end = MIN( jte, jde-1 )
4774
4775 IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
4776 config_flags%nested ) i_start = MAX( ids+1, its )
4777 IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
4778 config_flags%nested ) i_end = MIN( ide-2, ite )
4779 IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
4780 config_flags%nested ) j_start = MAX( jds+1, jts )
4781 IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
4782 config_flags%nested ) j_end = MIN( jde-2, jte )
4783 IF ( config_flags%periodic_x ) i_start = its
4784 IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
4785
4786 DO j = j_start, j_end
4787 DO k = kts, ktf
4788 DO i = i_start, i_end
4789 zxavg(i,k,j) = 0.25 * ( zx(i,k ,j) + zx(i+1,k ,j) + &
4790 zx(i,k+1,j) + zx(i+1,k+1,j) )
4791 zyavg(i,k,j) = 0.25 * ( zy(i,k ,j) + zy(i,k ,j+1) + &
4792 zy(i,k+1,j) + zy(i,k+1,j+1) )
4793 END DO
4794 END DO
4795 END DO
4796
4797 ! Begin calculating production of turbulence due to shear. The approach
4798 ! is to add together contributions from six terms, each of which is the
4799 ! square of a deformation that is then multiplied by an exchange
4800 ! coefficiant. The same exchange coefficient is assumed for horizontal
4801 ! and vertical coefficients for some of the terms (the vertical value is
4802 ! the one used).
4803
4804 ! For defor11.
4805
4806 DO j = j_start, j_end
4807 DO k = kts, ktf
4808 DO i = i_start, i_end
4809 tendency(i,k,j) = tendency(i,k,j) + 0.5 * &
4810 mu(i,j) * xkmh(i,k,j) * ( ( defor11(i,k,j) )**2 )
4811 END DO
4812 END DO
4813 END DO
4814
4815 ! For defor22.
4816
4817 DO j = j_start, j_end
4818 DO k = kts, ktf
4819 DO i = i_start, i_end
4820 tendency(i,k,j) = tendency(i,k,j) + 0.5 * &
4821 mu(i,j) * xkmh(i,k,j) * ( ( defor22(i,k,j) )**2 )
4822 END DO
4823 END DO
4824 END DO
4825
4826 ! For defor33.
4827
4828 DO j = j_start, j_end
4829 DO k = kts, ktf
4830 DO i = i_start, i_end
4831 tendency(i,k,j) = tendency(i,k,j) + 0.5 * &
4832 mu(i,j) * xkmv(i,k,j) * ( ( defor33(i,k,j) )**2 )
4833 END DO
4834 END DO
4835 END DO
4836
4837 ! For defor12.
4838
4839 DO j = j_start, j_end
4840 DO k = kts, ktf
4841 DO i = i_start, i_end
4842 avg(i,k,j) = 0.25 * &
4843 ( ( defor12(i ,k,j)**2 ) + ( defor12(i ,k,j+1)**2 ) + &
4844 ( defor12(i+1,k,j)**2 ) + ( defor12(i+1,k,j+1)**2 ) )
4845 END DO
4846 END DO
4847 END DO
4848
4849 DO j = j_start, j_end
4850 DO k = kts, ktf
4851 DO i = i_start, i_end
4852 tendency(i,k,j) = tendency(i,k,j) + mu(i,j) * xkmh(i,k,j) * avg(i,k,j)
4853 END DO
4854 END DO
4855 END DO
4856
4857 ! For defor13.
4858
4859 DO j = j_start, j_end
4860 DO k = kts+1, ktf
4861 DO i = i_start, i_end+1
4862 tmp2(i,k,j) = defor13(i,k,j)
4863 END DO
4864 END DO
4865 END DO
4866
4867 DO j = j_start, j_end
4868 DO i = i_start, i_end+1
4869 tmp2(i,kts ,j) = 0.0
4870 tmp2(i,ktf+1,j) = 0.0
4871 END DO
4872 END DO
4873
4874 DO j = j_start, j_end
4875 DO k = kts, ktf
4876 DO i = i_start, i_end
4877 avg(i,k,j) = 0.25 * &
4878 ( ( tmp2(i ,k+1,j)**2 ) + ( tmp2(i ,k,j)**2 ) + &
4879 ( tmp2(i+1,k+1,j)**2 ) + ( tmp2(i+1,k,j)**2 ) )
4880 END DO
4881 END DO
4882 END DO
4883
4884 DO j = j_start, j_end
4885 DO k = kts, ktf
4886 DO i = i_start, i_end
4887 tendency(i,k,j) = tendency(i,k,j) + mu(i,j) * xkmv(i,k,j) * avg(i,k,j)
4888 END DO
4889 END DO
4890 END DO
4891
4892 !MARTA: add the drag at the surface; WCS 040331
4893 K=KTS
4894
4895 cd0 = config_flags%tke_drag_coefficient ! drag coefficient set
4896 ! in namelist.input
4897 DO j = j_start, j_end
4898 DO i = i_start, i_end
4899
4900 absU=0.5*sqrt((u(i,k,j)+u(i+1,k,j))**2+(v(i,k,j)+v(i,k,j+1))**2)
4901
4902 tendency(i,k,j) = tendency(i,k,j) + &
4903 mu(i,j)*( (u(i,k,j)+u(i+1,k,j))*0.5* &
4904 cd0*absU*defor13(i,kts+1,j))
4905
4906 END DO
4907 END DO
4908 ! end of MARTA/WCS change
4909
4910 ! For defor23.
4911
4912 DO j = j_start, j_end+1
4913 DO k = kts+1, ktf
4914 DO i = i_start, i_end
4915 tmp2(i,k,j) = defor23(i,k,j)
4916 END DO
4917 END DO
4918 END DO
4919
4920 DO j = j_start, j_end+1
4921 DO i = i_start, i_end
4922 tmp2(i,kts, j) = 0.0
4923 tmp2(i,ktf+1,j) = 0.0
4924 END DO
4925 END DO
4926
4927 DO j = j_start, j_end
4928 DO k = kts, ktf
4929 DO i = i_start, i_end
4930 avg(i,k,j) = 0.25 * &
4931 ( ( tmp2(i,k+1,j )**2 ) + ( tmp2(i,k,j )**2) + &
4932 ( tmp2(i,k+1,j+1)**2 ) + ( tmp2(i,k,j+1)**2) )
4933 END DO
4934 END DO
4935 END DO
4936
4937 DO j = j_start, j_end
4938 DO k = kts, ktf
4939 DO i = i_start, i_end
4940 tendency(i,k,j) = tendency(i,k,j) + mu(i,j) * xkmv(i,k,j) * avg(i,k,j)
4941 END DO
4942 END DO
4943 END DO
4944
4945 !MARTA: add the drag at the surface; WCS 040331
4946 K=KTS
4947
4948 cd0 = config_flags%tke_drag_coefficient ! constant drag coefficient
4949 ! set in namelist.input
4950 DO j = j_start, j_end
4951 DO i = i_start, i_end
4952
4953 absU=0.5*sqrt((u(i,k,j)+u(i+1,k,j))**2+(v(i,k,j)+v(i,k,j+1))**2)
4954
4955 tendency(i,k,j) = tendency(i,k,j) + &
4956 mu(i,j)*( (v(i,k,j)+v(i,k,j+1))*0.5* &
4957 cd0*absU*defor23(i,kts+1,j))
4958
4959 END DO
4960 END DO
4961 ! end of MARTA/WCS change
4962
4963 END SUBROUTINE tke_shear
4964
4965 !=======================================================================
4966 !=======================================================================
4967
4968 SUBROUTINE compute_diff_metrics( config_flags, ph, phb, z, rdz, rdzw, &
4969 zx, zy, rdx, rdy, &
4970 ids, ide, jds, jde, kds, kde, &
4971 ims, ime, jms, jme, kms, kme, &
4972 its, ite, jts, jte, kts, kte )
4973
4974 !-----------------------------------------------------------------------
4975 ! Begin declarations.
4976
4977 IMPLICIT NONE
4978
4979 TYPE( grid_config_rec_type ), INTENT( IN ) &
4980 :: config_flags
4981
4982 INTEGER, INTENT( IN ) &
4983 :: ids, ide, jds, jde, kds, kde, &
4984 ims, ime, jms, jme, kms, kme, &
4985 its, ite, jts, jte, kts, kte
4986
4987 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
4988 :: ph, phb
4989
4990 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( OUT ) &
4991 :: rdz, rdzw, zx, zy, z
4992
4993 REAL, INTENT( IN ) &
4994 :: rdx, rdy
4995
4996 ! Local variables.
4997
4998 INTEGER &
4999 :: i, j, k, i_start, i_end, j_start, j_end, ktf
5000
5001 ! End declarations.
5002 !-----------------------------------------------------------------------
5003
5004 ktf = MIN( kte, kde-1 )
5005
5006 ! Bug fix, WCS, 22 april 2002.
5007
5008 ! We need rdzw in halo for average to u and v points.
5009
5010 j_start = jts-1
5011 j_end = jte
5012
5013 ! Begin with dz computations.
5014
5015 DO j = j_start, j_end
5016
5017 IF ( ( j_start >= jts ) .AND. ( j_end <= MIN( jte, jde-1 ) ) ) THEN
5018 i_start = its-1
5019 i_end = ite
5020 ELSE
5021 i_start = its
5022 i_end = MIN( ite, ide-1 )
5023 END IF
5024
5025 ! Compute z at w points for rdz and rdzw computations. We'll switch z
5026 ! to z at p points before returning
5027
5028 DO k = 1, kte
5029
5030 ! Bug fix, WCS, 22 april 2002
5031
5032 DO i = i_start, i_end
5033 z(i,k,j) = ( ph(i,k,j) + phb(i,k,j) ) / g
5034 END DO
5035 END DO
5036
5037 DO k = 1, ktf
5038 DO i = i_start, i_end
5039 rdzw(i,k,j) = 1.0 / ( z(i,k+1,j) - z(i,k,j) )
5040 END DO
5041 END DO
5042
5043 DO k = 2, ktf
5044 DO i = i_start, i_end
5045 rdz(i,k,j) = 2.0 / ( z(i,k+1,j) - z(i,k-1,j) )
5046 END DO
5047 END DO
5048
5049 ! Bug fix, WCS, 22 april 2002; added the following code
5050
5051 DO i = i_start, i_end
5052 rdz(i,1,j) = 2./(z(i,2,j)-z(i,1,j))
5053 END DO
5054
5055 END DO
5056
5057 ! End bug fix.
5058
5059 ! Now compute zx and zy; we'll assume that the halo for ph and phb is
5060 ! properly filled.
5061
5062 i_start = its
5063 i_end = MIN( ite, ide-1 )
5064 j_start = jts
5065 j_end = MIN( jte, jde-1 )
5066
5067 DO j = j_start, j_end
5068 DO k = 1, kte
5069 DO i = MAX( ids+1, its ), i_end
5070 zx(i,k,j) = rdx * ( phb(i,k,j) - phb(i-1,k,j) ) / g
5071 END DO
5072 END DO
5073 END DO
5074
5075 DO j = j_start, j_end
5076 DO k = 1, kte
5077 DO i = MAX( ids+1, its ), i_end
5078 zx(i,k,j) = zx(i,k,j) + rdx * ( ph(i,k,j) - ph(i-1,k,j) ) / g
5079 END DO
5080 END DO
5081 END DO
5082
5083 DO j = MAX( jds+1, jts ), j_end
5084 DO k = 1, kte
5085 DO i = i_start, i_end
5086 zy(i,k,j) = rdy * ( phb(i,k,j) - phb(i,k,j-1) ) / g
5087 END DO
5088 END DO
5089 END DO
5090
5091 DO j = MAX( jds+1, jts ), j_end
5092 DO k = 1, kte
5093 DO i = i_start, i_end
5094 zy(i,k,j) = zy(i,k,j) + rdy * ( ph(i,k,j) - ph(i,k,j-1) ) / g
5095 END DO
5096 END DO
5097 END DO
5098
5099 ! Some b.c. on zx and zy.
5100
5101 IF ( .NOT. config_flags%periodic_x ) THEN
5102
5103 IF ( ite == ide ) THEN
5104 DO j = j_start, j_end
5105 DO k = 1, ktf
5106 zx(ide,k,j) = 0.0
5107 END DO
5108 END DO
5109 END IF
5110
5111 IF ( its == ids ) THEN
5112 DO j = j_start, j_end
5113 DO k = 1, ktf
5114 zx(ids,k,j) = 0.0
5115 END DO
5116 END DO
5117 END IF
5118
5119 ELSE
5120
5121 IF ( ite == ide ) THEN
5122 DO j=j_start,j_end
5123 DO k=1,ktf
5124 zx(ide,k,j) = rdx * ( phb(ide,k,j) - phb(ide-1,k,j) ) / g
5125 END DO
5126 END DO
5127
5128 DO j = j_start, j_end
5129 DO k = 1, ktf
5130 zx(ide,k,j) = zx(ide,k,j) + rdx * ( ph(ide,k,j) - ph(ide-1,k,j) ) / g
5131 END DO
5132 END DO
5133 END IF
5134
5135 IF ( its == ids ) THEN
5136 DO j = j_start, j_end
5137 DO k = 1, ktf
5138 zx(ids,k,j) = rdx * ( phb(ids,k,j) - phb(ids-1,k,j) ) / g
5139 END DO
5140 END DO
5141
5142 DO j =j_start,j_end
5143 DO k =1,ktf
5144 zx(ids,k,j) = zx(ids,k,j) + rdx * ( ph(ids,k,j) - ph(ids-1,k,j) ) / g
5145 END DO
5146 END DO
5147 END IF
5148
5149 END IF
5150
5151 IF ( .NOT. config_flags%periodic_y ) THEN
5152
5153 IF ( jte == jde ) THEN
5154 DO k =1, ktf
5155 DO i =i_start, i_end
5156 zy(i,k,jde) = 0.0
5157 END DO
5158 END DO
5159 END IF
5160
5161 IF ( jts == jds ) THEN
5162 DO k =1, ktf
5163 DO i =i_start, i_end
5164 zy(i,k,jds) = 0.0
5165 END DO
5166 END DO
5167 END IF
5168
5169 ELSE
5170
5171 IF ( jte == jde ) THEN
5172 DO j=j_start, j_end
5173 DO k=1, ktf
5174 zy(i,k,jde) = rdy * ( phb(i,k,jde) - phb(i,k,jde-1) ) / g
5175 END DO
5176 END DO
5177
5178 DO j = j_start, j_end
5179 DO k = 1, ktf
5180 zy(i,k,jde) = zy(i,k,jde) + rdy * ( ph(i,k,jde) - ph(i,k,jde-1) ) / g
5181 END DO
5182 END DO
5183 END IF
5184
5185 IF ( jts == jds ) THEN
5186 DO j = j_start, j_end
5187 DO k = 1, ktf
5188 zy(i,k,jds) = rdy * ( phb(i,k,jds) - phb(i,k,jds-1) ) / g
5189 END DO
5190 END DO
5191
5192 DO j = j_start, j_end
5193 DO k = 1, ktf
5194 zy(i,k,jds) = zy(i,k,jds) + rdy * ( ph(i,k,jds) - ph(i,k,jds-1) ) / g
5195 END DO
5196 END DO
5197 END IF
5198
5199 END IF
5200
5201 ! Calculate z at p points.
5202
5203 DO j = j_start, j_end
5204 DO k = 1, ktf
5205 DO i = i_start, i_end
5206 z(i,k,j) = 0.5 * &
5207 ( ph(i,k,j) + phb(i,k,j) + ph(i,k+1,j) + phb(i,k+1,j) ) / g
5208 END DO
5209 END DO
5210 END DO
5211
5212 END SUBROUTINE compute_diff_metrics
5213
5214 !=======================================================================
5215 !=======================================================================
5216
5217 END MODULE module_diffusion_em
5218
5219 !=======================================================================
5220 !=======================================================================
5221
5222