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