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