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