| blob | c8d4b4d62895bf0109274150d325a585b22936ad |
1 !wrf:MODEL_LAYER:DYNAMICS
2 !
4 #if (RWORDSIZE == 4)
5 # define VPOWX vspowx
6 # define VPOW vspow
7 #else
8 # define VPOWX vpowx
9 # define VPOW vpow
10 #endif
13 MODULE module_big_step_utilities_em
15 USE module_model_constants
16 USE module_state_description, only: p_qg, p_qs, p_qi, gdscheme, tiedtkescheme, kfetascheme, g3scheme, &
17 p_qv, param_first_scalar, p_qr, p_qc
18 USE module_configure, ONLY : grid_config_rec_type
20 CONTAINS
22 !-------------------------------------------------------------------------------
24 SUBROUTINE calc_mu_uv ( config_flags, &
25 mu, mub, muu, muv, &
26 ids, ide, jds, jde, kds, kde, &
27 ims, ime, jms, jme, kms, kme, &
28 its, ite, jts, jte, kts, kte )
30 IMPLICIT NONE
32 ! Input data
34 TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
36 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
37 ims, ime, jms, jme, kms, kme, &
38 its, ite, jts, jte, kts, kte
40 REAL, DIMENSION( ims:ime , jms:jme ) , INTENT( OUT) :: muu, muv
41 REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu, mub
43 ! local stuff
45 INTEGER :: i, j, itf, jtf, im, jm
47 !<DESCRIPTION>
48 !
49 ! calc_mu_uv calculates the full column dry-air mass at the staggered
50 ! horizontal velocity points (u,v) and places the results in muu and muv.
51 ! This routine uses the reference state (mub) and perturbation state (mu)
52 !
53 !</DESCRIPTION>
56 itf=ite
57 jtf=MIN(jte,jde-1)
59 IF ( ( its .NE. ids ) .AND. ( ite .NE. ide ) ) THEN
60 DO j=jts,jtf
61 DO i=its,itf
62 muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
63 ENDDO
64 ENDDO
65 ELSE IF ( ( its .EQ. ids ) .AND. ( ite .NE. ide ) ) THEN
66 DO j=jts,jtf
67 DO i=its+1,itf
68 muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
69 ENDDO
70 ENDDO
71 i=its
72 im = its
73 if(config_flags%periodic_x) im = its-1
74 DO j=jts,jtf
75 ! muu(i,j) = mu(i,j) +mub(i,j)
76 ! fix for periodic b.c., 13 march 2004, wcs
77 muu(i,j) = 0.5*(mu(i,j)+mu(im,j)+mub(i,j)+mub(im,j))
78 ENDDO
79 ELSE IF ( ( its .NE. ids ) .AND. ( ite .EQ. ide ) ) THEN
80 DO j=jts,jtf
81 DO i=its,itf-1
82 muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
83 ENDDO
84 ENDDO
85 i=ite
86 im = ite-1
87 if(config_flags%periodic_x) im = ite
88 DO j=jts,jtf
89 ! muu(i,j) = mu(i-1,j) +mub(i-1,j)
90 ! fix for periodic b.c., 13 march 2004, wcs
91 muu(i,j) = 0.5*(mu(i-1,j)+mu(im,j)+mub(i-1,j)+mub(im,j))
92 ENDDO
93 ELSE IF ( ( its .EQ. ids ) .AND. ( ite .EQ. ide ) ) THEN
94 DO j=jts,jtf
95 DO i=its+1,itf-1
96 muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
97 ENDDO
98 ENDDO
99 i=its
100 im = its
101 if(config_flags%periodic_x) im = its-1
102 DO j=jts,jtf
103 ! muu(i,j) = mu(i,j) +mub(i,j)
104 ! fix for periodic b.c., 13 march 2004, wcs
105 muu(i,j) = 0.5*(mu(i,j)+mu(im,j)+mub(i,j)+mub(im,j))
106 ENDDO
107 i=ite
108 im = ite-1
109 if(config_flags%periodic_x) im = ite
110 DO j=jts,jtf
111 ! muu(i,j) = mu(i-1,j) +mub(i-1,j)
112 ! fix for periodic b.c., 13 march 2004, wcs
113 muu(i,j) = 0.5*(mu(i-1,j)+mu(im,j)+mub(i-1,j)+mub(im,j))
114 ENDDO
115 END IF
117 itf=MIN(ite,ide-1)
118 jtf=jte
120 IF ( ( jts .NE. jds ) .AND. ( jte .NE. jde ) ) THEN
121 DO j=jts,jtf
122 DO i=its,itf
123 muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
124 ENDDO
125 ENDDO
126 ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .NE. jde ) ) THEN
127 DO j=jts+1,jtf
128 DO i=its,itf
129 muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
130 ENDDO
131 ENDDO
132 j=jts
133 jm = jts
134 if(config_flags%periodic_y) jm = jts-1
135 DO i=its,itf
136 ! muv(i,j) = mu(i,j) +mub(i,j)
137 ! fix for periodic b.c., 13 march 2004, wcs
138 muv(i,j) = 0.5*(mu(i,j)+mu(i,jm)+mub(i,j)+mub(i,jm))
139 ENDDO
140 ELSE IF ( ( jts .NE. jds ) .AND. ( jte .EQ. jde ) ) THEN
141 DO j=jts,jtf-1
142 DO i=its,itf
143 muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
144 ENDDO
145 ENDDO
146 j=jte
147 jm = jte-1
148 if(config_flags%periodic_y) jm = jte
149 DO i=its,itf
150 muv(i,j) = mu(i,j-1) +mub(i,j-1)
151 ! fix for periodic b.c., 13 march 2004, wcs
152 muv(i,j) = 0.5*(mu(i,j-1)+mu(i,jm)+mub(i,j-1)+mub(i,jm))
153 ENDDO
154 ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .EQ. jde ) ) THEN
155 DO j=jts+1,jtf-1
156 DO i=its,itf
157 muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
158 ENDDO
159 ENDDO
160 j=jts
161 jm = jts
162 if(config_flags%periodic_y) jm = jts-1
163 DO i=its,itf
164 ! muv(i,j) = mu(i,j) +mub(i,j)
165 ! fix for periodic b.c., 13 march 2004, wcs
166 muv(i,j) = 0.5*(mu(i,j)+mu(i,jm)+mub(i,j)+mub(i,jm))
167 ENDDO
168 j=jte
169 jm = jte-1
170 if(config_flags%periodic_y) jm = jte
171 DO i=its,itf
172 ! muv(i,j) = mu(i,j-1) +mub(i,j-1)
173 ! fix for periodic b.c., 13 march 2004, wcs
174 muv(i,j) = 0.5*(mu(i,j-1)+mu(i,jm)+mub(i,j-1)+mub(i,jm))
175 ENDDO
176 END IF
178 END SUBROUTINE calc_mu_uv
180 !-------------------------------------------------------------------------------
182 SUBROUTINE calc_mu_uv_1 ( config_flags, &
183 mu, muu, muv, &
184 ids, ide, jds, jde, kds, kde, &
185 ims, ime, jms, jme, kms, kme, &
186 its, ite, jts, jte, kts, kte )
188 IMPLICIT NONE
190 ! Input data
192 TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
194 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
195 ims, ime, jms, jme, kms, kme, &
196 its, ite, jts, jte, kts, kte
198 REAL, DIMENSION( ims:ime , jms:jme ) , INTENT( OUT) :: muu, muv
199 REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
201 ! local stuff
203 INTEGER :: i, j, itf, jtf, im, jm
205 !<DESCRIPTION>
206 !
207 ! calc_mu_uv calculates the full column dry-air mass at the staggered
208 ! horizontal velocity points (u,v) and places the results in muu and muv.
209 ! This routine uses the full state (mu)
210 !
211 !</DESCRIPTION>
213 itf=ite
214 jtf=MIN(jte,jde-1)
216 IF ( ( its .NE. ids ) .AND. ( ite .NE. ide ) ) THEN
217 DO j=jts,jtf
218 DO i=its,itf
219 muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j))
220 ENDDO
221 ENDDO
222 ELSE IF ( ( its .EQ. ids ) .AND. ( ite .NE. ide ) ) THEN
223 DO j=jts,jtf
224 DO i=its+1,itf
225 muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j))
226 ENDDO
227 ENDDO
228 i=its
229 im = its
230 if(config_flags%periodic_x) im = its-1
231 DO j=jts,jtf
232 muu(i,j) = 0.5*(mu(i,j)+mu(im,j))
233 ENDDO
234 ELSE IF ( ( its .NE. ids ) .AND. ( ite .EQ. ide ) ) THEN
235 DO j=jts,jtf
236 DO i=its,itf-1
237 muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j))
238 ENDDO
239 ENDDO
240 i=ite
241 im = ite-1
242 if(config_flags%periodic_x) im = ite
243 DO j=jts,jtf
244 muu(i,j) = 0.5*(mu(i-1,j)+mu(im,j))
245 ENDDO
246 ELSE IF ( ( its .EQ. ids ) .AND. ( ite .EQ. ide ) ) THEN
247 DO j=jts,jtf
248 DO i=its+1,itf-1
249 muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j))
250 ENDDO
251 ENDDO
252 i=its
253 im = its
254 if(config_flags%periodic_x) im = its-1
255 DO j=jts,jtf
256 muu(i,j) = 0.5*(mu(i,j)+mu(im,j))
257 ENDDO
258 i=ite
259 im = ite-1
260 if(config_flags%periodic_x) im = ite
261 DO j=jts,jtf
262 muu(i,j) = 0.5*(mu(i-1,j)+mu(im,j))
263 ENDDO
264 END IF
266 itf=MIN(ite,ide-1)
267 jtf=jte
269 IF ( ( jts .NE. jds ) .AND. ( jte .NE. jde ) ) THEN
270 DO j=jts,jtf
271 DO i=its,itf
272 muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1))
273 ENDDO
274 ENDDO
275 ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .NE. jde ) ) THEN
276 DO j=jts+1,jtf
277 DO i=its,itf
278 muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1))
279 ENDDO
280 ENDDO
281 j=jts
282 jm = jts
283 if(config_flags%periodic_y) jm = jts-1
284 DO i=its,itf
285 muv(i,j) = 0.5*(mu(i,j)+mu(i,jm))
286 ENDDO
287 ELSE IF ( ( jts .NE. jds ) .AND. ( jte .EQ. jde ) ) THEN
288 DO j=jts,jtf-1
289 DO i=its,itf
290 muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1))
291 ENDDO
292 ENDDO
293 j=jte
294 jm = jte-1
295 if(config_flags%periodic_y) jm = jte
296 DO i=its,itf
297 muv(i,j) = 0.5*(mu(i,j-1)+mu(i,jm))
298 ENDDO
299 ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .EQ. jde ) ) THEN
300 DO j=jts+1,jtf-1
301 DO i=its,itf
302 muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1))
303 ENDDO
304 ENDDO
305 j=jts
306 jm = jts
307 if(config_flags%periodic_y) jm = jts-1
308 DO i=its,itf
309 muv(i,j) = 0.5*(mu(i,j)+mu(i,jm))
310 ENDDO
311 j=jte
312 jm = jte-1
313 if(config_flags%periodic_y) jm = jte
314 DO i=its,itf
315 muv(i,j) = 0.5*(mu(i,j-1)+mu(i,jm))
316 ENDDO
317 END IF
319 END SUBROUTINE calc_mu_uv_1
321 !-------------------------------------------------------------------------------
323 ! Map scale factor comments for this routine:
324 ! Locally not changed, but sent the correct map scale factors
325 ! from module_em (msfuy, msfvx, msfty)
327 SUBROUTINE couple_momentum ( muu, ru, u, msfu, &
328 muv, rv, v, msfv, msfv_inv, &
329 mut, rw, w, msft, &
330 ids, ide, jds, jde, kds, kde, &
331 ims, ime, jms, jme, kms, kme, &
332 its, ite, jts, jte, kts, kte )
334 IMPLICIT NONE
336 ! Input data
338 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
339 ims, ime, jms, jme, kms, kme, &
340 its, ite, jts, jte, kts, kte
342 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT( OUT) :: ru, rv, rw
344 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muu, muv, mut
345 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfu, msfv, msft, msfv_inv
347 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: u, v, w
349 ! Local data
351 INTEGER :: i, j, k, itf, jtf, ktf
353 !<DESCRIPTION>
354 !
355 ! couple_momentum couples the velocities to the full column mass and
356 ! the map factors.
357 !
358 !</DESCRIPTION>
360 ktf=MIN(kte,kde-1)
362 itf=ite
363 jtf=MIN(jte,jde-1)
365 DO j=jts,jtf
366 DO k=kts,ktf
367 DO i=its,itf
368 ru(i,k,j)=u(i,k,j)*muu(i,j)/msfu(i,j)
369 ENDDO
370 ENDDO
371 ENDDO
373 itf=MIN(ite,ide-1)
374 jtf=jte
376 DO j=jts,jtf
377 DO k=kts,ktf
378 DO i=its,itf
379 rv(i,k,j)=v(i,k,j)*muv(i,j)*msfv_inv(i,j)
380 ! rv(i,k,j)=v(i,k,j)*muv(i,j)/msfv(i,j)
381 ENDDO
382 ENDDO
383 ENDDO
385 itf=MIN(ite,ide-1)
386 jtf=MIN(jte,jde-1)
388 DO j=jts,jtf
389 DO k=kts,kte
390 DO i=its,itf
391 rw(i,k,j)=w(i,k,j)*mut(i,j)/msft(i,j)
392 ENDDO
393 ENDDO
394 ENDDO
396 END SUBROUTINE couple_momentum
398 !-------------------------------------------------------------------
400 SUBROUTINE calc_mu_staggered ( mu, mub, muu, muv, &
401 ids, ide, jds, jde, kds, kde, &
402 ims, ime, jms, jme, kms, kme, &
403 its, ite, jts, jte, kts, kte )
405 IMPLICIT NONE
407 ! Input data
409 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
410 ims, ime, jms, jme, kms, kme, &
411 its, ite, jts, jte, kts, kte
413 REAL, DIMENSION( ims:ime , jms:jme ) , INTENT( OUT) :: muu, muv
414 REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu, mub
416 ! local stuff
418 INTEGER :: i, j, itf, jtf
420 !<DESCRIPTION>
421 !
422 ! calc_mu_staggered calculates the full dry air mass at the staggered
423 ! velocity points (u,v).
424 !
425 !</DESCRIPTION>
427 itf=ite
428 jtf=MIN(jte,jde-1)
430 IF ( ( its .NE. ids ) .AND. ( ite .NE. ide ) ) THEN
431 DO j=jts,jtf
432 DO i=its,itf
433 muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
434 ENDDO
435 ENDDO
436 ELSE IF ( ( its .EQ. ids ) .AND. ( ite .NE. ide ) ) THEN
437 DO j=jts,jtf
438 DO i=its+1,itf
439 muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
440 ENDDO
441 ENDDO
442 i=its
443 DO j=jts,jtf
444 muu(i,j) = mu(i,j) +mub(i,j)
445 ENDDO
446 ELSE IF ( ( its .NE. ids ) .AND. ( ite .EQ. ide ) ) THEN
447 DO j=jts,jtf
448 DO i=its,itf-1
449 muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
450 ENDDO
451 ENDDO
452 i=ite
453 DO j=jts,jtf
454 muu(i,j) = mu(i-1,j) +mub(i-1,j)
455 ENDDO
456 ELSE IF ( ( its .EQ. ids ) .AND. ( ite .EQ. ide ) ) THEN
457 DO j=jts,jtf
458 DO i=its+1,itf-1
459 muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
460 ENDDO
461 ENDDO
462 i=its
463 DO j=jts,jtf
464 muu(i,j) = mu(i,j) +mub(i,j)
465 ENDDO
466 i=ite
467 DO j=jts,jtf
468 muu(i,j) = mu(i-1,j) +mub(i-1,j)
469 ENDDO
470 END IF
472 itf=MIN(ite,ide-1)
473 jtf=jte
475 IF ( ( jts .NE. jds ) .AND. ( jte .NE. jde ) ) THEN
476 DO j=jts,jtf
477 DO i=its,itf
478 muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
479 ENDDO
480 ENDDO
481 ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .NE. jde ) ) THEN
482 DO j=jts+1,jtf
483 DO i=its,itf
484 muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
485 ENDDO
486 ENDDO
487 j=jts
488 DO i=its,itf
489 muv(i,j) = mu(i,j) +mub(i,j)
490 ENDDO
491 ELSE IF ( ( jts .NE. jds ) .AND. ( jte .EQ. jde ) ) THEN
492 DO j=jts,jtf-1
493 DO i=its,itf
494 muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
495 ENDDO
496 ENDDO
497 j=jte
498 DO i=its,itf
499 muv(i,j) = mu(i,j-1) +mub(i,j-1)
500 ENDDO
501 ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .EQ. jde ) ) THEN
502 DO j=jts+1,jtf-1
503 DO i=its,itf
504 muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
505 ENDDO
506 ENDDO
507 j=jts
508 DO i=its,itf
509 muv(i,j) = mu(i,j) +mub(i,j)
510 ENDDO
511 j=jte
512 DO i=its,itf
513 muv(i,j) = mu(i,j-1) +mub(i,j-1)
514 ENDDO
515 END IF
517 END SUBROUTINE calc_mu_staggered
519 !-------------------------------------------------------------------------------
521 SUBROUTINE couple ( mu, mub, rfield, field, name, &
522 msf, &
523 ids, ide, jds, jde, kds, kde, &
524 ims, ime, jms, jme, kms, kme, &
525 its, ite, jts, jte, kts, kte )
527 IMPLICIT NONE
529 ! Input data
531 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
532 ims, ime, jms, jme, kms, kme, &
533 its, ite, jts, jte, kts, kte
535 CHARACTER(LEN=1) , INTENT(IN ) :: name
537 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT( OUT) :: rfield
539 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu, mub, msf
541 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field
543 ! Local data
545 INTEGER :: i, j, k, itf, jtf, ktf
546 REAL , DIMENSION(ims:ime,jms:jme) :: muu , muv
548 !<DESCRIPTION>
549 !
550 ! subroutine couple couples the input variable with the dry-air
551 ! column mass (mu).
552 !
553 !</DESCRIPTION>
556 ktf=MIN(kte,kde-1)
558 IF (name .EQ. 'u')THEN
560 CALL calc_mu_staggered ( mu, mub, muu, muv, &
561 ids, ide, jds, jde, kds, kde, &
562 ims, ime, jms, jme, kms, kme, &
563 its, ite, jts, jte, kts, kte )
565 itf=ite
566 jtf=MIN(jte,jde-1)
568 DO j=jts,jtf
569 DO k=kts,ktf
570 DO i=its,itf
571 rfield(i,k,j)=field(i,k,j)*muu(i,j)/msf(i,j)
572 ENDDO
573 ENDDO
574 ENDDO
576 ELSE IF (name .EQ. 'v')THEN
578 CALL calc_mu_staggered ( mu, mub, muu, muv, &
579 ids, ide, jds, jde, kds, kde, &
580 ims, ime, jms, jme, kms, kme, &
581 its, ite, jts, jte, kts, kte )
583 itf=ite
584 itf=MIN(ite,ide-1)
585 jtf=jte
587 DO j=jts,jtf
588 DO k=kts,ktf
589 DO i=its,itf
590 rfield(i,k,j)=field(i,k,j)*muv(i,j)/msf(i,j)
591 ENDDO
592 ENDDO
593 ENDDO
595 ELSE IF (name .EQ. 'w')THEN
596 itf=MIN(ite,ide-1)
597 jtf=MIN(jte,jde-1)
598 DO j=jts,jtf
599 DO k=kts,kte
600 DO i=its,itf
601 rfield(i,k,j)=field(i,k,j)*(mu(i,j)+mub(i,j))/msf(i,j)
602 ENDDO
603 ENDDO
604 ENDDO
606 ELSE IF (name .EQ. 'h')THEN
607 itf=MIN(ite,ide-1)
608 jtf=MIN(jte,jde-1)
609 DO j=jts,jtf
610 DO k=kts,kte
611 DO i=its,itf
612 rfield(i,k,j)=field(i,k,j)*(mu(i,j)+mub(i,j))
613 ENDDO
614 ENDDO
615 ENDDO
617 ELSE
618 itf=MIN(ite,ide-1)
619 jtf=MIN(jte,jde-1)
620 DO j=jts,jtf
621 DO k=kts,ktf
622 DO i=its,itf
623 rfield(i,k,j)=field(i,k,j)*(mu(i,j)+mub(i,j))
624 ENDDO
625 ENDDO
626 ENDDO
628 ENDIF
630 END SUBROUTINE couple
633 !-------------------------------------------------------------------------------
635 SUBROUTINE calc_ww_cp ( u, v, mup, mub, ww, &
636 rdx, rdy, msftx, msfty, &
637 msfux, msfuy, msfvx, msfvx_inv, &
638 msfvy, dnw, &
639 ids, ide, jds, jde, kds, kde, &
640 ims, ime, jms, jme, kms, kme, &
641 its, ite, jts, jte, kts, kte )
643 IMPLICIT NONE
645 ! Input data
648 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
649 ims, ime, jms, jme, kms, kme, &
650 its, ite, jts, jte, kts, kte
652 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: u, v
653 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mup, mub, &
654 msftx, msfty, &
655 msfux, msfuy, &
656 msfvx, msfvy, &
657 msfvx_inv
658 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
660 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: ww
661 REAL , INTENT(IN ) :: rdx, rdy
663 ! Local data
665 INTEGER :: i, j, k, itf, jtf, ktf
666 REAL , DIMENSION( its:ite ) :: dmdt
667 REAL , DIMENSION( its:ite, kts:kte ) :: divv
668 REAL , DIMENSION( its:ite+1, jts:jte+1 ) :: muu, muv
670 !<DESCRIPTION>
671 !
672 ! calc_ww calculates omega using the velocities (u,v) and the dry-air
673 ! column mass (mup+mub).
674 ! The algorithm integrates the continuity equation through the column
675 ! followed by a diagnosis of omega.
676 !
677 !</DESCRIPTION>
679 !<DESCRIPTION>
680 !
681 ! calc_ww_cp calculates omega using the velocities (u,v) and the
682 ! column mass mu.
683 !
684 !</DESCRIPTION>
686 jtf=MIN(jte,jde-1)
687 ktf=MIN(kte,kde-1)
688 itf=MIN(ite,ide-1)
690 ! mu coupled with the appropriate map factor
692 DO j=jts,jtf
693 DO i=its,min(ite+1,ide)
694 ! u is always coupled with my
695 muu(i,j) = 0.5*(mup(i,j)+mub(i,j)+mup(i-1,j)+mub(i-1,j))/msfuy(i,j)
696 ENDDO
697 ENDDO
699 DO j=jts,min(jte+1,jde)
700 DO i=its,itf
701 ! v is always coupled with mx
702 ! muv(i,j) = 0.5*(mup(i,j)+mub(i,j)+mup(i,j-1)+mub(i,j-1))/msfvx(i,j)
703 muv(i,j) = 0.5*(mup(i,j)+mub(i,j)+mup(i,j-1)+mub(i,j-1))*msfvx_inv(i,j)
704 ENDDO
705 ENDDO
707 DO j=jts,jtf
709 DO i=its,ite
710 dmdt(i) = 0.
711 ww(i,1,j) = 0.
712 ww(i,kte,j) = 0.
713 ENDDO
715 ! Comments on the modifications for map scale factors
716 ! ADT eqn 47 / my (putting rho -> 'mu') is:
717 ! (1/my) partial d mu/dt = -mx partial d/dx(mu u/my)
718 ! -mx partial d/dy(mu v/mx)
719 ! -partial d/dz(mu w/my)
720 !
721 ! Using nu instead of z the last term becomes:
722 ! -partial d/dnu(mu (dnu/dt)/my)
723 !
724 ! Integrating with respect to nu over ALL levels, with dnu/dt=0 at top
725 ! and bottom, the last term becomes = 0
726 !
727 ! Integral|bot->top[(1/my) partial d mu/dt]dnu =
728 ! Integral|bot->top[-mx partial d/dx(mu u/my)
729 ! -mx partial d/dy(mu v/mx)]dnu
730 !
731 ! muu='mu'[on u]/my, muv='mu'[on v]/mx
732 ! (1/my) partial d mu/dt is independent of nu
733 ! => LHS = Integral|bot->top[con]dnu = conservation*(-1) = -dmdt
734 !
735 ! => dmdt = mx*Integral|bot->top[partial d/dx(mu u/my) +
736 ! partial d/dy(mu v/mx)]dnu
737 ! => dmdt = sum_bot->top[divv]
738 ! where
739 ! divv=mx*[partial d/dx(mu u/my) + partial d/dy(mu v/mx)]*delta nu
741 DO k=kts,ktf
742 DO i=its,itf
744 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)) &
745 +rdy*(muv(i,j+1)*v(i,k,j+1)-muv(i,j)*v(i,k,j)) )
747 ! dmdt(i) = dmdt(i) + dnw(k)* ( rdx*(ru(i+1,k,j)-ru(i,k,j)) &
748 ! +rdy*(rv(i,k,j+1)-rv(i,k,j)) )
750 dmdt(i) = dmdt(i) + divv(i,k)
753 ENDDO
754 ENDDO
756 ! Further map scale factor notes:
757 ! Now integrate from bottom to top, level by level:
758 ! mu dnu/dt/my [k+1] = mu dnu/dt/my [k] + [-(1/my) partial d mu/dt
759 ! -mx partial d/dx(mu u/my)
760 ! -mx partial d/dy(mu v/mx)]*dnu[k->k+1]
761 ! ww [k+1] = ww [k] -(1/my) partial d mu/dt * dnu[k->k+1] - divv[k]
762 ! = ww [k] -dmdt * dnw[k] - divv[k]
764 DO k=2,ktf
765 DO i=its,itf
767 ! ww(i,k,j)=ww(i,k-1,j) &
768 ! - dnw(k-1)* ( dmdt(i) &
769 ! +rdx*(ru(i+1,k-1,j)-ru(i,k-1,j)) &
770 ! +rdy*(rv(i,k-1,j+1)-rv(i,k-1,j)) )
772 ww(i,k,j)=ww(i,k-1,j) - dnw(k-1)*dmdt(i) - divv(i,k-1)
774 ENDDO
775 ENDDO
776 ENDDO
779 END SUBROUTINE calc_ww_cp
782 !-------------------------------------------------------------------------------
784 SUBROUTINE calc_cq ( moist, cqu, cqv, cqw, n_moist, &
785 ids, ide, jds, jde, kds, kde, &
786 ims, ime, jms, jme, kms, kme, &
787 its, ite, jts, jte, kts, kte )
789 IMPLICIT NONE
791 ! Input data
793 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
794 ims, ime, jms, jme, kms, kme, &
795 its, ite, jts, jte, kts, kte
797 INTEGER , INTENT(IN ) :: n_moist
800 REAL, DIMENSION( ims:ime, kms:kme , jms:jme , n_moist ), INTENT(IN ) :: moist
802 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT( OUT) :: cqu, cqv, cqw
804 ! Local stuff
806 REAL :: qtot
808 INTEGER :: i, j, k, itf, jtf, ktf, ispe
810 !<DESCRIPTION>
811 !
812 ! calc_cq calculates moist coefficients for the momentum equations.
813 !
814 !</DESCRIPTION>
816 itf=ite
817 jtf=MIN(jte,jde-1)
818 ktf=MIN(kte,kde-1)
820 IF( n_moist >= PARAM_FIRST_SCALAR ) THEN
822 DO j=jts,jtf
823 DO k=kts,ktf
824 DO i=its,itf
825 qtot = 0.
826 !DEC$ loop count(3)
827 DO ispe=PARAM_FIRST_SCALAR,n_moist
828 qtot = qtot + moist(i,k,j,ispe) + moist(i-1,k,j,ispe)
829 ENDDO
830 ! qtot = 0.5*( moist(i ,k,j,1)+moist(i ,k,j,2)+moist(i ,k,j,3)+ &
831 ! & moist(i-1,k,j,1)+moist(i-1,k,j,2)+moist(i-1,k,j,3) )
832 ! cqu(i,k,j) = 1./(1.+qtot)
833 cqu(i,k,j) = 1./(1.+0.5*qtot)
834 ENDDO
835 ENDDO
836 ENDDO
838 itf=MIN(ite,ide-1)
839 jtf=jte
841 DO j=jts,jtf
842 DO k=kts,ktf
843 DO i=its,itf
844 qtot = 0.
845 !DEC$ loop count(3)
846 DO ispe=PARAM_FIRST_SCALAR,n_moist
847 qtot = qtot + moist(i,k,j,ispe) + moist(i,k,j-1,ispe)
848 ENDDO
849 ! qtot = 0.5*( moist(i,k,j ,1)+moist(i,k,j ,2)+moist(i,k,j ,3)+ &
850 ! & moist(i,k,j-1,1)+moist(i,k,j-1,2)+moist(i,k,j-1,3) )
851 ! cqv(i,k,j) = 1./(1.+qtot)
852 cqv(i,k,j) = 1./(1.+0.5*qtot)
853 ENDDO
854 ENDDO
855 ENDDO
857 itf=MIN(ite,ide-1)
858 jtf=MIN(jte,jde-1)
859 DO j=jts,jtf
860 DO k=kts+1,ktf
861 DO i=its,itf
862 qtot = 0.
863 !DEC$ loop count(3)
864 DO ispe=PARAM_FIRST_SCALAR,n_moist
865 qtot = qtot + moist(i,k,j,ispe) + moist(i,k-1,j,ispe)
866 ENDDO
867 ! qtot = 0.5*( moist(i,k ,j,1)+moist(i,k ,j,2)+moist(i,k-1,j,3)+ &
868 ! & moist(i,k-1,j,1)+moist(i,k-1,j,2)+moist(i,k ,j,3) )
869 ! cqw(i,k,j) = qtot
870 cqw(i,k,j) = 0.5*qtot
871 ENDDO
872 ENDDO
873 ENDDO
875 ELSE
877 DO j=jts,jtf
878 DO k=kts,ktf
879 DO i=its,itf
880 cqu(i,k,j) = 1.
881 ENDDO
882 ENDDO
883 ENDDO
885 itf=MIN(ite,ide-1)
886 jtf=jte
888 DO j=jts,jtf
889 DO k=kts,ktf
890 DO i=its,itf
891 cqv(i,k,j) = 1.
892 ENDDO
893 ENDDO
894 ENDDO
896 itf=MIN(ite,ide-1)
897 jtf=MIN(jte,jde-1)
898 DO j=jts,jtf
899 DO k=kts+1,ktf
900 DO i=its,itf
901 cqw(i,k,j) = 0.
902 ENDDO
903 ENDDO
904 ENDDO
906 END IF
908 END SUBROUTINE calc_cq
910 !----------------------------------------------------------------------
912 SUBROUTINE calc_alt ( alt, al, alb, &
913 ids, ide, jds, jde, kds, kde, &
914 ims, ime, jms, jme, kms, kme, &
915 its, ite, jts, jte, kts, kte )
917 IMPLICIT NONE
919 ! Input data
921 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
922 ims, ime, jms, jme, kms, kme, &
923 its, ite, jts, jte, kts, kte
925 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: alb, al
926 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT( OUT) :: alt
928 ! Local stuff
930 INTEGER :: i, j, k, itf, jtf, ktf
932 !<DESCRIPTION>
933 !
934 ! calc_alt computes the full inverse density
935 !
936 !</DESCRIPTION>
938 itf=MIN(ite,ide-1)
939 jtf=MIN(jte,jde-1)
940 ktf=MIN(kte,kde-1)
942 DO j=jts,jtf
943 DO k=kts,ktf
944 DO i=its,itf
945 alt(i,k,j) = al(i,k,j)+alb(i,k,j)
946 ENDDO
947 ENDDO
948 ENDDO
951 END SUBROUTINE calc_alt
953 !----------------------------------------------------------------------
955 SUBROUTINE calc_p_rho_phi ( moist, n_moist, hypsometric_opt, &
956 al, alb, mu, muts, ph, phb, p, pb, &
957 t, p0, t0, ptop, znu, znw, dnw, rdnw, &
958 rdn, non_hydrostatic, &
959 ids, ide, jds, jde, kds, kde, &
960 ims, ime, jms, jme, kms, kme, &
961 its, ite, jts, jte, kts, kte )
963 IMPLICIT NONE
965 ! Input data
967 LOGICAL , INTENT(IN ) :: non_hydrostatic
969 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
970 ims, ime, jms, jme, kms, kme, &
971 its, ite, jts, jte, kts, kte
973 INTEGER , INTENT(IN ) :: n_moist
974 INTEGER , INTENT(IN ) :: hypsometric_opt
976 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: alb, &
977 pb, &
978 t
980 REAL, DIMENSION( ims:ime , kms:kme , jms:jme, n_moist ), INTENT(IN ) :: moist
982 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT( OUT) :: al, p
984 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: ph, phb
986 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN ) :: mu, muts
988 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: znu, znw, dnw, rdnw, rdn
990 REAL, INTENT(IN ) :: t0, p0, ptop
992 ! Local stuff
994 INTEGER :: i, j, k, itf, jtf, ktf, ispe
995 REAL :: qvf, qtot, qf1, qf2
996 REAL, DIMENSION( its:ite) :: temp,cpovcv_v
997 REAL :: pfu, phm, pfd
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>
1012 itf=MIN(ite,ide-1)
1013 jtf=MIN(jte,jde-1)
1014 ktf=MIN(kte,kde-1)
1016 #ifndef INTELMKL
1017 cpovcv_v = cpovcv
1018 #endif
1020 IF (non_hydrostatic) THEN
1022 IF (hypsometric_opt == 1) THEN
1023 DO j=jts,jtf
1024 DO k=kts,ktf
1025 DO i=its,itf
1026 al(i,k,j)=-1./muts(i,j)*(alb(i,k,j)*mu(i,j) + rdnw(k)*(ph(i,k+1,j)-ph(i,k,j)))
1027 END DO
1028 END DO
1029 END DO
1030 ELSE IF (hypsometric_opt == 2) THEN
1032 ! The relation used to get specific volume, al, is: al = -dZ/dp,
1033 ! where dp = mut * d(eta). The pressure depth, dp, is replaced with
1034 ! p*(dp/p) ~ p*LOG((p+0.5dp)/(p-0.5dp)). Difference between dp and p*dLOG(p)
1035 ! is as follows: p*dLOG(p) - dp = 1/12*(dp/p)**3 + 1/90*(dp/p)**5 + ...
1036 ! Therefore, p*dLOG(p) is always larger than dp and the difference is
1037 ! in proportion to dp/p. TKW, 02/16/2010
1039 DO j=jts,jtf
1040 DO k=kts,ktf
1041 DO i=its,itf
1042 pfu = muts(i,j)*znw(k+1)+ptop
1043 pfd = muts(i,j)*znw(k) +ptop
1044 phm = muts(i,j)*znu(k) +ptop
1045 al(i,k,j) = (ph(i,k+1,j)-ph(i,k,j)+phb(i,k+1,j)-phb(i,k,j))/phm/LOG(pfd/pfu)-alb(i,k,j)
1046 END DO
1047 END DO
1048 END DO
1049 ELSE
1050 CALL wrf_error_fatal ( 'calc_p_rho_phi: hypsometric_opt should be 1 or 2' )
1051 END IF
1053 IF (n_moist >= PARAM_FIRST_SCALAR ) THEN
1055 DO j=jts,jtf
1056 DO k=kts,ktf
1057 DO i=its,itf
1058 qvf = 1.+rvovrd*moist(i,k,j,P_QV)
1059 temp(i)=(r_d*(t0+t(i,k,j))*qvf)/(p0*(al(i,k,j)+alb(i,k,j)))
1060 ENDDO
1061 #ifdef INTELMKL
1062 CALL VPOWX ( itf-its+1, temp(its), cpovcv, p(its,k,j) )
1063 #else
1064 ! use vector version from libmassv or from compat lib in frame/libmassv.F
1065 CALL VPOW ( p(its,k,j), temp(its), cpovcv_v(its), itf-its+1 )
1066 #endif
1067 DO i=its,itf
1068 p(i,k,j)= p(i,k,j)*p0-pb(i,k,j)
1069 ENDDO
1070 ENDDO
1071 ENDDO
1073 ELSE
1075 DO j=jts,jtf
1076 DO k=kts,ktf
1077 DO i=its,itf
1078 p(i,k,j)=p0*( (r_d*(t0+t(i,k,j)))/ &
1079 (p0*(al(i,k,j)+alb(i,k,j))) )**cpovcv &
1080 -pb(i,k,j)
1081 ENDDO
1082 ENDDO
1083 ENDDO
1085 END IF
1087 ELSE
1089 ! hydrostatic pressure, al, and ph1 calc; WCS, 5 sept 2001
1092 IF (n_moist >= PARAM_FIRST_SCALAR ) THEN
1094 DO j=jts,jtf
1096 k=ktf ! top layer
1097 DO i=its,itf
1099 qtot = 0.
1100 DO ispe=PARAM_FIRST_SCALAR,n_moist
1101 qtot = qtot + moist(i,k,j,ispe)
1102 ENDDO
1103 qf2 = 1.
1104 qf1 = qtot*qf2
1106 p(i,k,j) = - 0.5*(mu(i,j)+qf1*muts(i,j))/rdnw(k)/qf2
1107 qvf = 1.+rvovrd*moist(i,k,j,P_QV)
1108 al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* &
1109 (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1111 ENDDO
1113 DO k=ktf-1,kts,-1 ! remaining layers, integrate down
1114 DO i=its,itf
1116 qtot = 0.
1117 DO ispe=PARAM_FIRST_SCALAR,n_moist
1118 qtot = qtot + 0.5*( moist(i,k ,j,ispe) + moist(i,k+1,j,ispe) )
1119 ENDDO
1120 qf2 = 1.
1121 qf1 = qtot*qf2
1123 p(i,k,j) = p(i,k+1,j) - (mu(i,j) + qf1*muts(i,j))/qf2/rdn(k+1)
1124 qvf = 1.+rvovrd*moist(i,k,j,P_QV)
1125 al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* &
1126 (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1127 ENDDO
1128 ENDDO
1130 ENDDO
1132 ELSE
1134 DO j=jts,jtf
1136 k=ktf ! top layer
1137 DO i=its,itf
1139 qtot = 0.
1140 qf2 = 1.
1141 qf1 = qtot*qf2
1143 p(i,k,j) = - 0.5*(mu(i,j)+qf1*muts(i,j))/rdnw(k)/qf2
1144 qvf = 1.
1145 al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* &
1146 (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1148 ENDDO
1150 DO k=ktf-1,kts,-1 ! remaining layers, integrate down
1151 DO i=its,itf
1153 qtot = 0.
1154 qf2 = 1.
1155 qf1 = qtot*qf2
1157 p(i,k,j) = p(i,k+1,j) - (mu(i,j) + qf1*muts(i,j))/qf2/rdn(k+1)
1158 qvf = 1.
1159 al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* &
1160 (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1161 ENDDO
1162 ENDDO
1164 ENDDO
1166 END IF
1168 IF (hypsometric_opt == 1) THEN
1169 DO j=jts,jtf
1170 DO k=2,ktf+1 ! integrate hydrostatic equation for geopotential
1171 DO i=its,itf
1172 ph(i,k,j) = ph(i,k-1,j) - (dnw(k-1))*( &
1173 (muts(i,j))*al(i,k-1,j)+ &
1174 mu(i,j)*alb(i,k-1,j) )
1175 ENDDO
1176 ENDDO
1177 ENDDO
1178 ELSE
1180 ! Revised hypsometric eq.: dZ=-al*p*dLOG(p), where p is dry pressure
1182 DO j=jts,jtf
1183 DO i=its,itf
1184 ph(i,kts,j) = phb(i,kts,j)
1185 END DO
1187 DO k=kts+1,ktf+1
1188 DO i=its,itf
1189 pfu = muts(i,j)*znw(k) +ptop
1190 pfd = muts(i,j)*znw(k-1)+ptop
1191 phm = muts(i,j)*znu(k-1)+ptop
1192 ph(i,k,j) = ph(i,k-1,j) + (al(i,k-1,j)+alb(i,k-1,j))*phm*LOG(pfd/pfu)
1193 ENDDO
1194 ENDDO
1196 DO k=kts,ktf+1
1197 DO i=its,itf
1198 ph(i,k,j) = ph(i,k,j) - phb(i,k,j)
1199 END DO
1200 END DO
1201 END DO
1203 END IF
1205 END IF
1207 END SUBROUTINE calc_p_rho_phi
1209 !----------------------------------------------------------------------
1211 SUBROUTINE calc_php ( php, ph, phb, &
1212 ids, ide, jds, jde, kds, kde, &
1213 ims, ime, jms, jme, kms, kme, &
1214 its, ite, jts, jte, kts, kte )
1216 IMPLICIT NONE
1218 ! Input data
1220 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1221 ims, ime, jms, jme, kms, kme, &
1222 its, ite, jts, jte, kts, kte
1224 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: phb, ph
1225 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT( OUT) :: php
1227 ! Local stuff
1229 INTEGER :: i, j, k, itf, jtf, ktf
1231 !<DESCRIPTION>
1232 !
1233 ! calc_php calculates the full geopotential from the reference state
1234 ! geopotential and the perturbation geopotential (phb_ph).
1235 !
1236 !</DESCRIPTION>
1238 itf=MIN(ite,ide-1)
1239 jtf=MIN(jte,jde-1)
1240 ktf=MIN(kte,kde-1)
1242 DO j=jts,jtf
1243 DO k=kts,ktf
1244 DO i=its,itf
1245 php(i,k,j) = 0.5*(phb(i,k,j)+phb(i,k+1,j)+ph(i,k,j)+ph(i,k+1,j))
1246 ENDDO
1247 ENDDO
1248 ENDDO
1250 END SUBROUTINE calc_php
1252 !-------------------------------------------------------------------------------
1254 SUBROUTINE diagnose_w( ph_tend, ph_new, ph_old, w, mu, dt, &
1255 u, v, ht, &
1256 cf1, cf2, cf3, rdx, rdy, &
1257 msftx, msfty, &
1258 ids, ide, jds, jde, kds, kde, &
1259 ims, ime, jms, jme, kms, kme, &
1260 its, ite, jts, jte, kts, kte )
1262 IMPLICIT NONE
1264 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1265 ims, ime, jms, jme, kms, kme, &
1266 its, ite, jts, jte, kts, kte
1268 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: ph_tend, &
1269 ph_new, &
1270 ph_old, &
1271 u, &
1272 v
1275 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT( OUT) :: w
1277 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mu, ht, msftx, msfty
1279 REAL, INTENT(IN ) :: dt, cf1, cf2, cf3, rdx, rdy
1281 INTEGER :: i, j, k, itf, jtf
1283 itf=MIN(ite,ide-1)
1284 jtf=MIN(jte,jde-1)
1286 !<DESCRIPTION>
1287 !
1288 ! diagnose_w diagnoses the vertical velocity from the geopoential equation.
1289 ! Used with the hydrostatic option.
1290 !
1291 !</DESCRIPTION>
1293 DO j = jts, jtf
1295 ! lower b.c. on w
1297 ! Notes on map scale factors:
1298 ! Chain rule: if Z=Z(X,Y) [true at the surface] then
1299 ! dZ/dt = dZ/dX * dX/dt + dZ/dY * dY/dt, U=dX/dt, V=dY/dt
1300 ! Using capitals to denote actual values
1301 ! In mapped values, u=U, v=V, z=Z, 1/dX=mx/dx, 1/dY=my/dy
1302 ! => w = dz/dt = mx u dz/dx + my v dz/dy
1303 ! [where dz/dx is just the surface height change between x
1304 ! gridpoints, and dz/dy is the change between y gridpoints]
1305 ! [NB: cf1, cf2 and cf3 do vertical weighting of u or v values
1306 ! nearest the surface]
1308 ! Previously msft multiplied by rdy and rdx terms.
1309 ! Now msfty multiplies rdy term, and msftx multiplies msftx term
1310 DO i = its, itf
1311 w(i,1,j)= msfty(i,j)*.5*rdy*( &
1312 (ht(i,j+1)-ht(i,j )) &
1313 *(cf1*v(i,1,j+1)+cf2*v(i,2,j+1)+cf3*v(i,3,j+1)) &
1314 +(ht(i,j )-ht(i,j-1)) &
1315 *(cf1*v(i,1,j )+cf2*v(i,2,j )+cf3*v(i,3,j )) ) &
1316 +msftx(i,j)*.5*rdx*( &
1317 (ht(i+1,j)-ht(i,j )) &
1318 *(cf1*u(i+1,1,j)+cf2*u(i+1,2,j)+cf3*u(i+1,3,j)) &
1319 +(ht(i,j )-ht(i-1,j)) &
1320 *(cf1*u(i ,1,j)+cf2*u(i ,2,j)+cf3*u(i ,3,j)) )
1321 ENDDO
1323 ! use geopotential equation to diagnose w
1325 ! Further notes on map scale factors
1326 ! If ph_tend contains: -mx partial d/dx(mu rho u/my)
1327 ! -mx partial d/dy(phi mu v/mx)
1328 ! -partial d/dz(phi mu w/my)
1329 ! then phi eqn is: partial d/dt(mu phi/my) = ph_tend + mu g w/my
1330 ! => w = [my/(mu*g)]*[partial d/dt(mu phi/my) - ph_tend]
1332 DO k = 2, kte
1333 DO i = its, itf
1334 w(i,k,j) = msfty(i,j)*( (ph_new(i,k,j)-ph_old(i,k,j))/dt &
1335 - ph_tend(i,k,j)/mu(i,j) )/g
1337 ENDDO
1338 ENDDO
1340 ENDDO
1342 END SUBROUTINE diagnose_w
1344 !-------------------------------------------------------------------------------
1346 SUBROUTINE rhs_ph( ph_tend, u, v, ww, &
1347 ph, ph_old, phb, w, &
1348 mut, muu, muv, &
1349 fnm, fnp, &
1350 rdnw, cfn, cfn1, rdx, rdy, &
1351 msfux, msfuy, msfvx, &
1352 msfvx_inv, msfvy, &
1353 msftx, msfty, &
1354 non_hydrostatic, &
1355 config_flags, &
1356 ids, ide, jds, jde, kds, kde, &
1357 ims, ime, jms, jme, kms, kme, &
1358 its, ite, jts, jte, kts, kte )
1359 IMPLICIT NONE
1361 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
1363 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1364 ims, ime, jms, jme, kms, kme, &
1365 its, ite, jts, jte, kts, kte
1367 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: &
1368 u, &
1369 v, &
1370 ww, &
1371 ph, &
1372 ph_old, &
1373 phb, &
1374 w
1376 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: ph_tend
1378 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: muu, muv, mut, &
1379 msfux, msfuy, &
1380 msfvx, msfvy, &
1381 msftx, msfty, &
1382 msfvx_inv
1384 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, fnm, fnp
1386 REAL, INTENT(IN ) :: cfn, cfn1, rdx, rdy
1388 LOGICAL, INTENT(IN ) :: non_hydrostatic
1390 ! Local stuff
1392 INTEGER :: i, j, k, itf, jtf, ktf, kz, i_start, j_start
1393 REAL :: ur, ul, ub, vr, vl, vb
1394 REAL, DIMENSION(its:ite,kts:kte) :: wdwn
1396 INTEGER :: advective_order
1398 LOGICAL :: specified
1400 !<DESCRIPTION>
1401 !
1402 ! rhs_ph calculates the large-timestep tendency terms for the geopotential
1403 ! equation. These terms include the advection and "gw". The geopotential
1404 ! equation is cast in advective form, so we don't use the flux form advection
1405 ! algorithms here.
1406 !
1407 !</DESCRIPTION>
1409 specified = .false.
1410 if(config_flags%specified .or. config_flags%nested) specified = .true.
1412 advective_order = config_flags%h_sca_adv_order
1414 itf=MIN(ite,ide-1)
1415 jtf=MIN(jte,jde-1)
1416 ktf=MIN(kte,kde-1)
1418 ! Notes on map scale factors (WCS, 2 march 2008)
1419 ! phi equation is: mu/my d/dt(phi) = -(1/my) mx mu u d/dx(phi)
1420 ! -(1/my) my mu v d/dy(phi)
1421 ! - omega d/d_eta(phi)
1422 ! +mu g w/my
1423 !
1424 ! A little further explanation...
1425 ! The tendency term we are computing here is for mu/my d/dt(phi). It is advective form
1426 ! but it is multiplied be mu/my. It will be decoupled from (mu/my) when the implicit w-phi
1427 ! solution is computed in subourine advance_w. The formulation dates from the early
1428 ! days of the mass coordinate model when we were testing both a flux and an advective formulation
1429 ! for the geopotential equation and different forms of the vertical momentum equation and the
1430 ! vertically implicit solver.
1432 ! advective form for the geopotential equation
1434 DO j = jts, jtf
1436 DO k = 2, kte
1437 DO i = its, itf
1438 wdwn(i,k) = .5*(ww(i,k,j)+ww(i,k-1,j))*rdnw(k-1) &
1439 *(ph(i,k,j)-ph(i,k-1,j)+phb(i,k,j)-phb(i,k-1,j))
1440 ENDDO
1441 ENDDO
1443 ! RHS term 3 is: - omega partial d/dnu(phi)
1445 DO k = 2, kte-1
1446 DO i = its, itf
1447 ph_tend(i,k,j) = ph_tend(i,k,j) &
1448 - (fnm(k)*wdwn(i,k+1)+fnp(k)*wdwn(i,k))
1449 ENDDO
1450 ENDDO
1452 ENDDO
1454 IF (non_hydrostatic) THEN ! add in "gw" term.
1455 DO j = jts, jtf ! in hydrostatic mode, "gw" will be diagnosed
1456 ! after the timestep to give us "w"
1457 DO i = its, itf
1458 ph_tend(i,kde,j) = 0.
1459 ENDDO
1461 DO k = 2, kte
1462 DO i = its, itf
1463 ! phi equation RHS term 4
1464 ph_tend(i,k,j) = ph_tend(i,k,j) + mut(i,j)*g*w(i,k,j)/msfty(i,j)
1465 ENDDO
1466 ENDDO
1468 ENDDO
1470 END IF
1472 ! Notes on map scale factors:
1473 ! RHS terms 1 and 2 are: -(1/my) mx u mu partial d/dx(phi)
1474 ! -(1/my) my v mu partial d/dy(phi)
1476 IF (advective_order <= 2) THEN
1478 ! y (v) advection
1480 i_start = its
1481 j_start = jts
1482 itf=MIN(ite,ide-1)
1483 jtf=MIN(jte,jde-1)
1485 IF ( (config_flags%open_ys .or. specified) .and. jts == jds ) j_start = jts+1
1486 IF ( (config_flags%open_ye .or. specified) .and. jte == jde ) jtf = jtf-2
1488 DO j = j_start, jtf
1490 DO k = 2, kte-1
1491 DO i = i_start, itf
1492 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* &
1493 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)* &
1494 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1495 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j )* &
1496 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1497 ENDDO
1498 ENDDO
1500 k = kte
1501 DO i = i_start, itf
1502 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* &
1503 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)* &
1504 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1505 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j )* &
1506 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1507 ENDDO
1509 ENDDO
1511 ! x (u) advection
1513 i_start = its
1514 j_start = jts
1515 itf=MIN(ite,ide-1)
1516 jtf=MIN(jte,jde-1)
1518 IF ( (config_flags%open_xs .or. specified) .and. its == ids ) i_start = its+1
1519 IF ( (config_flags%open_xe .or. specified) .and. ite == ide ) itf = itf-2
1521 DO j = j_start, jtf
1523 DO k = 2, kte-1
1524 DO i = i_start, itf
1525 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))* &
1526 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)* &
1527 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1528 +muu(i ,j)*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j)* &
1529 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1530 ENDDO
1531 ENDDO
1533 k = kte
1534 DO i = i_start, itf
1535 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))* &
1536 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)* &
1537 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1538 +muu(i ,j)*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux( i,j)* &
1539 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1540 ENDDO
1542 ENDDO
1544 ELSE IF (advective_order <= 4) THEN
1546 ! y (v) advection
1548 i_start = its
1549 j_start = jts
1550 itf=MIN(ite,ide-1)
1551 jtf=MIN(jte,jde-1)
1553 IF ( (config_flags%open_ys .or. specified) .and. jts == jds ) j_start = jts+2
1554 IF ( (config_flags%open_ye .or. specified) .and. jte == jde ) jtf = jtf-3
1556 DO j = j_start, jtf
1558 DO k = 2, kte-1
1559 DO i = i_start, itf
1560 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))*( &
1561 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1) &
1562 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j ))* (1./12.)*( &
1563 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1564 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1565 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1566 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1569 ENDDO
1570 ENDDO
1572 k = kte
1573 DO i = i_start, itf
1574 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))*( &
1575 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1) &
1576 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j ))* (1./12.)*( &
1577 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1578 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1579 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1580 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1582 ENDDO
1584 ENDDO
1586 ! pick up near boundary rows using 2nd order stencil for open and specified conditions
1588 IF ( (config_flags%open_ys .or. specified) .and. jts <= jds+1 ) THEN
1590 j = jds+1
1591 DO k = 2, kte-1
1592 DO i = i_start, itf
1593 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* &
1594 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)* &
1595 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1596 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j )* &
1597 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1598 ENDDO
1599 ENDDO
1601 k = kte
1602 DO i = i_start, itf
1603 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* &
1604 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)* &
1605 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1606 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j )* &
1607 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1608 ENDDO
1610 END IF
1612 IF ( (config_flags%open_ye .or. specified) .and. jte >= jde-2 ) THEN
1614 j = jde-2
1615 DO k = 2, kte-1
1616 DO i = i_start, itf
1617 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* &
1618 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)* &
1619 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1620 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j )* &
1621 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1622 ENDDO
1623 ENDDO
1625 k = kte
1626 DO i = i_start, itf
1627 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* &
1628 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)* &
1629 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1630 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j )* &
1631 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1632 ENDDO
1634 END IF
1636 ! x (u) advection
1638 i_start = its
1639 j_start = jts
1640 itf=MIN(ite,ide-1)
1641 jtf=MIN(jte,jde-1)
1643 IF ( (config_flags%open_xs) .and. its == ids ) i_start = its+2
1644 IF ( (config_flags%open_xe) .and. ite == ide ) itf = itf-3
1646 DO j = j_start, jtf
1648 DO k = 2, kte-1
1649 DO i = i_start, itf
1650 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*( &
1651 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j) &
1652 +muu(i,j )*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j) )* (1./12.)*( &
1653 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1654 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1655 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1656 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1657 ENDDO
1658 ENDDO
1660 k = kte
1661 DO i = i_start, itf
1662 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*( &
1663 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j) &
1664 +muu(i,j )*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux(i ,j) )* (1./12.)*( &
1665 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1666 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1667 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1668 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1669 ENDDO
1671 ENDDO
1673 ! pick up near boundary rows using 2nd order stencil for open and specified conditions
1675 IF ( (config_flags%open_xs .or. specified) .and. its <= ids+1 ) THEN
1677 i = ids + 1
1679 DO j = j_start, jtf
1680 DO k = 2, kte-1
1681 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))* &
1682 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)* &
1683 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1684 +muu(i ,j)*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j)* &
1685 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1686 ENDDO
1687 ENDDO
1689 k = kte
1690 DO j = j_start, jtf
1691 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))* &
1692 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)* &
1693 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1694 +muu(i ,j)*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux( i,j)* &
1695 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1696 ENDDO
1698 END IF
1700 IF ( (config_flags%open_xe .or. specified) .and. ite >= ide-2 ) THEN
1702 i = ide-2
1703 DO j = j_start, jtf
1704 DO k = 2, kte-1
1705 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))* &
1706 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)* &
1707 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1708 +muu(i ,j)*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j)* &
1709 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1710 ENDDO
1711 ENDDO
1713 k = kte
1714 DO j = j_start, jtf
1715 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))* &
1716 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)* &
1717 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1718 +muu(i ,j)*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux( i,j)* &
1719 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1720 ENDDO
1722 END IF
1724 !--------------------------
1726 ELSE IF (advective_order <= 6) THEN
1728 !!! NOTE: this branch has been looked at and fixed with changes for overdecomposition
1729 !!! the branches covering the other advective_order cases have not. 20090923. JM
1731 ! y (v) advection
1733 i_start = its
1734 j_start = jts
1735 itf=MIN(ite,ide-1)
1736 jtf=MIN(jte,jde-1)
1738 IF (config_flags%open_ys .or. specified ) j_start = max(jts,jds+3)
1739 IF (config_flags%open_ye .or. specified ) jtf = min(jtf,jde-4)
1741 DO j = j_start, jtf
1743 DO k = 2, kte-1
1744 DO i = i_start, itf
1745 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* ( &
1746 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1) &
1747 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j ) )* (1./60.)*( &
1748 45.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1749 -9.*(ph(i,k,j+2)-ph(i,k,j-2)) &
1750 +(ph(i,k,j+3)-ph(i,k,j-3)) &
1751 +45.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1752 -9.*(phb(i,k,j+2)-phb(i,k,j-2)) &
1753 +(phb(i,k,j+3)-phb(i,k,j-3)) ) )
1756 ENDDO
1757 ENDDO
1759 k = kte
1760 DO i = i_start, itf
1761 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* ( &
1762 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1) &
1763 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j ) )* (1./60.)*( &
1764 45.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1765 -9.*(ph(i,k,j+2)-ph(i,k,j-2)) &
1766 +(ph(i,k,j+3)-ph(i,k,j-3)) &
1767 +45.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1768 -9.*(phb(i,k,j+2)-phb(i,k,j-2)) &
1769 +(phb(i,k,j+3)-phb(i,k,j-3)) ) )
1771 ENDDO
1773 ENDDO
1775 ! 4th order stencil for open or specified conditions two in form the boundary
1777 IF ( (config_flags%open_ys .or. specified) .and. jts <= jds+2 .and. jds+2 <= jte ) THEN
1779 j = jds+2
1780 DO k = 2, kte-1
1781 DO i = i_start, itf
1782 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* ( &
1783 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1) &
1784 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j ) )* (1./12.)*( &
1785 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1786 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1787 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1788 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1791 ENDDO
1792 ENDDO
1794 k = kte
1795 DO i = i_start, itf
1796 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* ( &
1797 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1) &
1798 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j) )* (1./12.)*( &
1799 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1800 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1801 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1802 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1804 ENDDO
1806 END IF
1808 IF ( (config_flags%open_ye .or. specified) .and. jts <= jde-3 .and. jde-3 <= jte ) THEN
1809 j = jde-3
1810 DO k = 2, kte-1
1811 DO i = i_start, itf
1812 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* ( &
1813 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1) &
1814 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j) )* (1./12.)*( &
1815 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1816 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1817 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1818 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1821 ENDDO
1822 ENDDO
1824 k = kte
1825 DO i = i_start, itf
1826 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* ( &
1827 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1) &
1828 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j) )* (1./12.)*( &
1829 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1830 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1831 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1832 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1834 ENDDO
1836 END IF
1838 ! 2nd order stencil for open and specified conditions one row in from boundary
1840 IF ( (config_flags%open_ys .or. specified) .and. jts <= jds+1 .and. jds+1 <= jte ) THEN
1842 j = jds+1
1843 DO k = 2, kte-1
1844 DO i = i_start, itf
1845 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* &
1846 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)* &
1847 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1848 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j )* &
1849 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1850 ENDDO
1851 ENDDO
1853 k = kte
1854 DO i = i_start, itf
1855 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* &
1856 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)* &
1857 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1858 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j )* &
1859 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1860 ENDDO
1862 END IF
1864 IF ( (config_flags%open_ye .or. specified) .and. jts <= jde-2 .and. jde-2 <= jte ) THEN
1866 j = jde-2
1867 DO k = 2, kte-1
1868 DO i = i_start, itf
1869 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* &
1870 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)* &
1871 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1872 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j )* &
1873 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1874 ENDDO
1875 ENDDO
1877 k = kte
1878 DO i = i_start, itf
1879 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* &
1880 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)* &
1881 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1882 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j )* &
1883 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1884 ENDDO
1886 END IF
1888 ! x (u) advection
1890 i_start = its
1891 j_start = jts
1892 itf=MIN(ite,ide-1)
1893 jtf=MIN(jte,jde-1)
1895 IF (config_flags%open_xs .or. specified ) i_start = max(its,ids+3)
1896 IF (config_flags%open_xe .or. specified ) itf = min(itf,ide-4)
1898 DO j = j_start, jtf
1900 DO k = 2, kte-1
1901 DO i = i_start, itf
1902 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*( &
1903 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j) &
1904 +muu(i,j )*(u(i,k,j )+u(i,k-1,j ))*msfux(i,j) )* (1./60.)*( &
1905 45.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1906 -9.*(ph(i+2,k,j)-ph(i-2,k,j)) &
1907 +(ph(i+3,k,j)-ph(i-3,k,j)) &
1908 +45.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1909 -9.*(phb(i+2,k,j)-phb(i-2,k,j)) &
1910 +(phb(i+3,k,j)-phb(i-3,k,j)) ) )
1911 ENDDO
1912 ENDDO
1914 k = kte
1915 DO i = i_start, itf
1916 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*( &
1917 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j) &
1918 +muu(i,j )*(cfn*u(i ,k-1,j)+cfn1*u(i,k-2,j))*msfux(i,j) )* (1./60.)*( &
1919 45.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1920 -9.*(ph(i+2,k,j)-ph(i-2,k,j)) &
1921 +(ph(i+3,k,j)-ph(i-3,k,j)) &
1922 +45.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1923 -9.*(phb(i+2,k,j)-phb(i-2,k,j)) &
1924 +(phb(i+3,k,j)-phb(i-3,k,j)) ) )
1925 ENDDO
1927 ENDDO
1929 ! 4th order stencil two in from the boundary for open and specified conditions
1931 IF ( (config_flags%open_xs) .and. its <= ids+2 .and. ids+2 <= ite ) THEN
1932 i = ids + 2
1933 DO j = j_start, jtf
1934 DO k = 2, kte-1
1935 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*( &
1936 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j) &
1937 +muu(i,j )*(u(i,k,j )+u(i,k-1,j ))*msfux(i,j) )* (1./12.)*( &
1938 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1939 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1940 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1941 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1942 ENDDO
1943 k = kte
1944 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*( &
1945 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j) &
1946 +muu(i,j )*(cfn*u(i ,k-1,j)+cfn1*u(i,k-2,j))*msfux(i,j) )* (1./12.)*( &
1947 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1948 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1949 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1950 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1952 ENDDO
1953 END IF
1955 IF ( (config_flags%open_xe) .and. its <= ide-3 .and. ide-3 <= ite ) THEN
1956 i = ide-3
1957 DO j = j_start, jtf
1958 DO k = 2, kte-1
1959 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*( &
1960 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j) &
1961 +muu(i,j )*(u(i,k,j )+u(i,k-1,j ))*msfux(i,j) )* (1./12.)*( &
1962 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1963 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1964 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1965 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1966 ENDDO
1967 k = kte
1968 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*( &
1969 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j) &
1970 +muu(i,j )*(cfn*u(i ,k-1,j)+cfn1*u(i,k-2,j))*msfux(i,j) )* (1./12.)*( &
1971 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1972 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1973 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1974 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1976 ENDDO
1977 END IF
1979 ! 2nd order stencil for open and specified conditions one in from the boundary
1981 IF ( (config_flags%open_xs .or. specified) .and. its <= ids+1 .and. ids+1 <= ite ) THEN
1983 i = ids + 1
1985 DO j = j_start, jtf
1986 DO k = 2, kte-1
1987 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))* &
1988 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)* &
1989 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1990 +muu(i ,j)*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j)* &
1991 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1992 ENDDO
1993 ENDDO
1995 k = kte
1996 DO j = j_start, jtf
1997 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))* &
1998 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)* &
1999 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
2000 +muu(i ,j)*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux( i,j)* &
2001 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
2002 ENDDO
2004 END IF
2006 IF ( (config_flags%open_xe .or. specified) .and. its <= ide-2 .and. ide-2 <= ite ) THEN
2008 i = ide-2
2009 DO j = j_start, jtf
2010 DO k = 2, kte-1
2011 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))* &
2012 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)* &
2013 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
2014 +muu(i ,j)*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j)* &
2015 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
2016 ENDDO
2017 ENDDO
2019 k = kte
2020 DO j = j_start, jtf
2021 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))* &
2022 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)* &
2023 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
2024 +muu(i ,j)*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux( i,j)* &
2025 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
2026 ENDDO
2028 END IF
2030 END IF ! 6th order advection
2032 ! lateral open boundary conditions,
2033 ! start with north and south (y) boundaries
2035 i_start = its
2036 itf=MIN(ite,ide-1)
2038 ! south
2040 IF ( (config_flags%open_ys) .and. jts == jds ) THEN
2042 j=jts
2044 DO k=2,kde
2045 kz = min(k,kde-1)
2046 DO i = its,itf
2047 vb =.5*( fnm(kz)*(v(i,kz ,j+1)+v(i,kz ,j )) &
2048 +fnp(kz)*(v(i,kz-1,j+1)+v(i,kz-1,j )) )
2049 vl=amin1(vb,0.)
2050 ph_tend(i,k,j)=ph_tend(i,k,j)-rdy*mut(i,j)*( &
2051 +vl*(ph_old(i,k,j+1)-ph_old(i,k,j)))
2052 ENDDO
2053 ENDDO
2055 END IF
2057 ! north
2059 IF ( (config_flags%open_ye) .and. jte == jde ) THEN
2061 j=jte-1
2063 DO k=2,kde
2064 kz = min(k,kde-1)
2065 DO i = its,itf
2066 vb=.5*( fnm(kz)*(v(i,kz ,j+1)+v(i,kz ,j)) &
2067 +fnp(kz)*(v(i,kz-1,j+1)+v(i,kz-1,j)) )
2068 vr=amax1(vb,0.)
2069 ph_tend(i,k,j)=ph_tend(i,k,j)-rdy*mut(i,j)*( &
2070 +vr*(ph_old(i,k,j)-ph_old(i,k,j-1)))
2071 ENDDO
2072 ENDDO
2074 END IF
2076 ! now the east and west (y) boundaries
2078 j_start = its
2079 jtf=MIN(jte,jde-1)
2081 ! west
2083 IF ( (config_flags%open_xs) .and. its == ids ) THEN
2085 i=its
2087 DO j = jts,jtf
2088 DO k=2,kde-1
2089 kz = k
2090 ub =.5*( fnm(kz)*(u(i+1,kz ,j)+u(i ,kz ,j)) &
2091 +fnp(kz)*(u(i+1,kz-1,j)+u(i ,kz-1,j)) )
2092 ul=amin1(ub,0.)
2093 ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*( &
2094 +ul*(ph_old(i+1,k,j)-ph_old(i,k,j)))
2095 ENDDO
2097 k = kde
2098 kz = k
2099 ub =.5*( fnm(kz)*(u(i+1,kz ,j)+u(i ,kz ,j)) &
2100 +fnp(kz)*(u(i+1,kz-1,j)+u(i ,kz-1,j)) )
2101 ul=amin1(ub,0.)
2102 ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*( &
2103 +ul*(ph_old(i+1,k,j)-ph_old(i,k,j)))
2104 ENDDO
2106 END IF
2108 ! east
2110 IF ( (config_flags%open_xe) .and. ite == ide ) THEN
2112 i = ite-1
2114 DO j = jts,jtf
2115 DO k=2,kde-1
2116 kz = k
2117 ub=.5*( fnm(kz)*(u(i+1,kz ,j)+u(i,kz ,j)) &
2118 +fnp(kz)*(u(i+1,kz-1,j)+u(i,kz-1,j)) )
2119 ur=amax1(ub,0.)
2120 ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*( &
2121 +ur*(ph_old(i,k,j)-ph_old(i-1,k,j)))
2122 ENDDO
2124 k = kde
2125 kz = k-1
2126 ub=.5*( fnm(kz)*(u(i+1,kz ,j)+u(i,kz ,j)) &
2127 +fnp(kz)*(u(i+1,kz-1,j)+u(i,kz-1,j)) )
2128 ur=amax1(ub,0.)
2129 ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*( &
2130 +ur*(ph_old(i,k,j)-ph_old(i-1,k,j)))
2132 ENDDO
2134 END IF
2136 END SUBROUTINE rhs_ph
2139 !-------------------------------------------------------------------------------
2141 SUBROUTINE horizontal_pressure_gradient( ru_tend,rv_tend, &
2142 ph,alt,p,pb,al,php,cqu,cqv, &
2143 muu,muv,mu,fnm,fnp,rdnw, &
2144 cf1,cf2,cf3,rdx,rdy,msfux,msfuy,&
2145 msfvx,msfvy,msftx,msfty, &
2146 config_flags, non_hydrostatic, &
2147 top_lid, &
2148 ids, ide, jds, jde, kds, kde, &
2149 ims, ime, jms, jme, kms, kme, &
2150 its, ite, jts, jte, kts, kte )
2152 IMPLICIT NONE
2154 ! Input data
2157 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2159 LOGICAL, INTENT (IN ) :: non_hydrostatic, top_lid
2161 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2162 ims, ime, jms, jme, kms, kme, &
2163 its, ite, jts, jte, kts, kte
2165 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: &
2166 ph, &
2167 alt, &
2168 al, &
2169 p, &
2170 pb, &
2171 php, &
2172 cqu, &
2173 cqv
2176 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: &
2177 ru_tend, &
2178 rv_tend
2180 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: muu, muv, mu, &
2181 msfux, msfuy, &
2182 msfvx, msfvy, &
2183 msftx, msfty
2185 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, fnm, fnp
2187 REAL, INTENT(IN ) :: rdx, rdy, cf1, cf2, cf3
2189 INTEGER :: i,j,k, itf, jtf, ktf, i_start, j_start
2190 REAL, DIMENSION( ims:ime, kms:kme ) :: dpn
2191 REAL :: dpx, dpy
2193 LOGICAL :: specified
2195 !<DESCRIPTION>
2196 !
2197 ! horizontal_pressure_gradient calculates the
2198 ! horizontal pressure gradient terms for the large-timestep tendency
2199 ! in the horizontal momentum equations (u,v).
2200 !
2201 !</DESCRIPTION>
2203 specified = .false.
2204 if(config_flags%specified .or. config_flags%nested) specified = .true.
2206 ! Notes on map scale factors:
2207 ! Calculates the pressure gradient terms in ADT eqns 44 and 45
2208 ! With upper rho -> 'mu', these are:
2209 ! Eqn 30: -mu*(mx/my)*(1/rho)*partial dp/dx
2210 ! Eqn 31: -mu*(my/mx)*(1/rho)*partial dp/dy
2211 !
2212 ! As we are on nu, rather than height, surfaces:
2213 !
2214 ! mu dp/dx = mu alpha partial dp'/dx + (nu mu partial dmubar/dx) alpha'
2215 ! + mu partial dphi'/dx + (partial dphi/dx)*(partial dp'/dnu - mu')
2216 !
2217 ! mu dp/dy = mu alpha partial dp'/dy + (nu mu partial dmubar/dy) alpha'
2218 ! + mu partial dphi'/dy + (partial dphi/dy)*(partial dp'/dnu - mu')
2220 ! start with the north-south (y) pressure gradient
2222 itf=MIN(ite,ide-1)
2223 jtf=jte
2224 ktf=MIN(kte,kde-1)
2225 i_start = its
2226 j_start = jts
2227 IF ( (config_flags%open_ys .or. specified .or. &
2228 config_flags%nested .or. config_flags%polar ) .and. jts == jds ) j_start = jts+1
2229 IF ( (config_flags%open_ye .or. specified .or. &
2230 config_flags%nested .or. config_flags%polar ) .and. jte == jde ) jtf = jtf-1
2232 DO j = j_start, jtf
2234 IF ( non_hydrostatic ) THEN
2236 k=1
2238 DO i = i_start, itf
2239 dpn(i,k) = .5*( cf1*(p(i,k ,j-1)+p(i,k ,j)) &
2240 +cf2*(p(i,k+1,j-1)+p(i,k+1,j)) &
2241 +cf3*(p(i,k+2,j-1)+p(i,k+2,j)) )
2242 dpn(i,kde) = 0.
2243 ENDDO
2244 IF (top_lid) THEN
2245 DO i = i_start, itf
2246 dpn(i,kde) = .5*( cf1*(p(i,kde-1,j-1)+p(i,kde-1,j)) &
2247 +cf2*(p(i,kde-2,j-1)+p(i,kde-2,j)) &
2248 +cf3*(p(i,kde-3,j-1)+p(i,kde-3,j)) )
2249 ENDDO
2250 ENDIF
2252 DO k=2,ktf
2253 DO i = i_start, itf
2254 dpn(i,k) = .5*( fnm(k)*(p(i,k ,j-1)+p(i,k ,j)) &
2255 +fnp(k)*(p(i,k-1,j-1)+p(i,k-1,j)) )
2256 END DO
2257 END DO
2259 ! ADT eqn 45: -mu*(my/mx)*(1/rho)*partial dp/dy
2260 ! [alt, al are 1/rho terms; muv, mu are NOT coupled]
2261 DO K=1,ktf
2262 DO i = i_start, itf
2263 ! Here are mu dp/dy terms 1-3
2264 dpy = (msfvy(i,j)/msfvx(i,j))*.5*rdy*muv(i,j)*( &
2265 (ph (i,k+1,j)-ph (i,k+1,j-1) + ph(i,k,j)-ph(i,k,j-1)) &
2266 +(alt(i,k ,j)+alt(i,k ,j-1))*(p (i,k,j)-p (i,k,j-1)) &
2267 +(al (i,k ,j)+al (i,k ,j-1))*(pb(i,k,j)-pb(i,k,j-1)) )
2268 ! Here is mu dp/dy term 4
2269 dpy = dpy + (msfvy(i,j)/msfvx(i,j))*rdy*(php(i,k,j)-php(i,k,j-1))* &
2270 (rdnw(k)*(dpn(i,k+1)-dpn(i,k))-.5*(mu(i,j-1)+mu(i,j)))
2271 rv_tend(i,k,j) = rv_tend(i,k,j)-cqv(i,k,j)*dpy
2272 END DO
2273 END DO
2275 ELSE
2277 ! ADT eqn 45: -mu*(my/mx)*(1/rho)*partial dp/dy
2278 ! [alt, al are 1/rho terms; muv, mu are NOT coupled]
2279 DO K=1,ktf
2280 DO i = i_start, itf
2281 ! Here are mu dp/dy terms 1-3; term 4 not needed if hydrostatic
2282 dpy = (msfvy(i,j)/msfvx(i,j))*.5*rdy*muv(i,j)*( &
2283 (ph (i,k+1,j)-ph (i,k+1,j-1) + ph(i,k,j)-ph(i,k,j-1)) &
2284 +(alt(i,k ,j)+alt(i,k ,j-1))*(p (i,k,j)-p (i,k,j-1)) &
2285 +(al (i,k ,j)+al (i,k ,j-1))*(pb(i,k,j)-pb(i,k,j-1)) )
2286 rv_tend(i,k,j) = rv_tend(i,k,j)-cqv(i,k,j)*dpy
2287 END DO
2288 END DO
2290 END IF
2292 ENDDO
2294 ! now the east-west (x) pressure gradient
2296 itf=ite
2297 jtf=MIN(jte,jde-1)
2298 ktf=MIN(kte,kde-1)
2299 i_start = its
2300 j_start = jts
2301 IF ( (config_flags%open_xs .or. specified .or. &
2302 config_flags%nested ) .and. its == ids ) i_start = its+1
2303 IF ( (config_flags%open_xe .or. specified .or. &
2304 config_flags%nested ) .and. ite == ide ) itf = itf-1
2305 IF ( config_flags%periodic_x ) i_start = its
2306 IF ( config_flags%periodic_x ) itf=ite
2308 DO j = j_start, jtf
2310 IF ( non_hydrostatic ) THEN
2312 k=1
2314 DO i = i_start, itf
2315 dpn(i,k) = .5*( cf1*(p(i-1,k ,j)+p(i,k ,j)) &
2316 +cf2*(p(i-1,k+1,j)+p(i,k+1,j)) &
2317 +cf3*(p(i-1,k+2,j)+p(i,k+2,j)) )
2318 dpn(i,kde) = 0.
2319 ENDDO
2320 IF (top_lid) THEN
2321 DO i = i_start, itf
2322 dpn(i,kde) = .5*( cf1*(p(i-1,kde-1,j)+p(i,kde-1,j)) &
2323 +cf2*(p(i-1,kde-2,j)+p(i,kde-2,j)) &
2324 +cf3*(p(i-1,kde-3,j)+p(i,kde-3,j)) )
2325 ENDDO
2326 ENDIF
2328 DO k=2,ktf
2329 DO i = i_start, itf
2330 dpn(i,k) = .5*( fnm(k)*(p(i-1,k ,j)+p(i,k ,j)) &
2331 +fnp(k)*(p(i-1,k-1,j)+p(i,k-1,j)) )
2332 END DO
2333 END DO
2335 ! ADT eqn 44: -mu*(mx/my)*(1/rho)*partial dp/dx
2336 ! [alt, al are 1/rho terms; muu, mu are NOT coupled]
2337 DO K=1,ktf
2338 DO i = i_start, itf
2339 ! Here are mu dp/dy terms 1-3
2340 dpx = (msfux(i,j)/msfuy(i,j))*.5*rdx*muu(i,j)*( &
2341 (ph (i,k+1,j)-ph (i-1,k+1,j) + ph(i,k,j)-ph(i-1,k,j)) &
2342 +(alt(i,k ,j)+alt(i-1,k ,j))*(p (i,k,j)-p (i-1,k,j)) &
2343 +(al (i,k ,j)+al (i-1,k ,j))*(pb(i,k,j)-pb(i-1,k,j)) )
2344 ! Here is mu dp/dy term 4
2345 dpx = dpx + (msfux(i,j)/msfuy(i,j))*rdx*(php(i,k,j)-php(i-1,k,j))* &
2346 (rdnw(k)*(dpn(i,k+1)-dpn(i,k))-.5*(mu(i-1,j)+mu(i,j)))
2347 ru_tend(i,k,j) = ru_tend(i,k,j)-cqu(i,k,j)*dpx
2348 END DO
2349 END DO
2351 ELSE
2353 ! ADT eqn 44: -mu*(mx/my)*(1/rho)*partial dp/dx
2354 ! [alt, al are 1/rho terms; muu, mu are NOT coupled]
2355 DO K=1,ktf
2356 DO i = i_start, itf
2357 ! Here are mu dp/dy terms 1-3; term 4 not needed if hydrostatic
2358 dpx = (msfux(i,j)/msfuy(i,j))*.5*rdx*muu(i,j)*( &
2359 (ph (i,k+1,j)-ph (i-1,k+1,j) + ph(i,k,j)-ph(i-1,k,j)) &
2360 +(alt(i,k ,j)+alt(i-1,k ,j))*(p (i,k,j)-p (i-1,k,j)) &
2361 +(al (i,k ,j)+al (i-1,k ,j))*(pb(i,k,j)-pb(i-1,k,j)) )
2362 ru_tend(i,k,j) = ru_tend(i,k,j)-cqu(i,k,j)*dpx
2363 END DO
2364 END DO
2366 END IF
2368 ENDDO
2370 END SUBROUTINE horizontal_pressure_gradient
2372 !-------------------------------------------------------------------------------
2374 SUBROUTINE pg_buoy_w( rw_tend, p, cqw, mu, mub, &
2375 rdnw, rdn, g, msftx, msfty, &
2376 ids, ide, jds, jde, kds, kde, &
2377 ims, ime, jms, jme, kms, kme, &
2378 its, ite, jts, jte, kts, kte )
2380 IMPLICIT NONE
2382 ! Input data
2384 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2385 ims, ime, jms, jme, kms, kme, &
2386 its, ite, jts, jte, kts, kte
2388 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: p
2389 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: cqw
2392 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: rw_tend
2394 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mub, mu, msftx, msfty
2396 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, rdn
2398 REAL, INTENT(IN ) :: g
2400 INTEGER :: itf, jtf, i, j, k
2401 REAL :: cq1, cq2
2404 !<DESCRIPTION>
2405 !
2406 ! pg_buoy_w calculates the
2407 ! vertical pressure gradient and buoyancy terms for the large-timestep
2408 ! tendency in the vertical momentum equation.
2409 !
2410 !</DESCRIPTION>
2412 ! BUOYANCY AND PRESSURE GRADIENT TERM IN W EQUATION AT TIME T
2414 ! Map scale factor notes
2415 ! ADT eqn 46 RHS terms 6 and 7 (where 7 is "-rho g")
2416 ! Dividing by my, and using mu and nu (see Klemp et al. eqns 32, 40)
2417 ! term 6: +(g/my) partial dp'/dnu
2418 ! term 7: -(g/my) mu'
2419 !
2420 ! For moisture-free atmosphere, cq1=1, cq2=0
2421 ! => (1./msft(i,j)) * g * [rdn(k)*{p(i,k,j)-p(i,k-1,j)}-mu(i,j)]
2423 itf=MIN(ite,ide-1)
2424 jtf=MIN(jte,jde-1)
2426 DO j = jts,jtf
2428 k=kde
2429 DO i=its,itf
2430 cq1 = 1./(1.+cqw(i,k-1,j))
2431 cq2 = cqw(i,k-1,j)*cq1
2432 rw_tend(i,k,j) = rw_tend(i,k,j)+(1./msfty(i,j))*g*( &
2433 cq1*2.*rdnw(k-1)*( -p(i,k-1,j)) &
2434 -mu(i,j)-cq2*mub(i,j) )
2435 END DO
2437 DO k = 2, kde-1
2438 DO i = its,itf
2439 cq1 = 1./(1.+cqw(i,k,j))
2440 cq2 = cqw(i,k,j)*cq1
2441 cqw(i,k,j) = cq1
2442 rw_tend(i,k,j) = rw_tend(i,k,j)+(1./msfty(i,j))*g*( &
2443 cq1*rdn(k)*(p(i,k,j)-p(i,k-1,j)) &
2444 -mu(i,j)-cq2*mub(i,j) )
2445 END DO
2446 ENDDO
2449 ENDDO
2451 END SUBROUTINE pg_buoy_w
2453 !-------------------------------------------------------------------------------
2455 SUBROUTINE w_damp( rw_tend, max_vert_cfl,max_horiz_cfl, &
2456 u, v, ww, w, mut, rdnw, &
2457 rdx, rdy, msfux, msfuy, &
2458 msfvx, msfvy, dt, &
2459 config_flags, &
2460 ids, ide, jds, jde, kds, kde, &
2461 ims, ime, jms, jme, kms, kme, &
2462 its, ite, jts, jte, kts, kte )
2464 USE module_llxy
2465 IMPLICIT NONE
2467 ! Input data
2469 TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
2471 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2472 ims, ime, jms, jme, kms, kme, &
2473 its, ite, jts, jte, kts, kte
2475 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: u, v, ww, w
2477 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: rw_tend
2479 REAL, INTENT(OUT) :: max_vert_cfl
2480 REAL, INTENT(OUT) :: max_horiz_cfl
2481 REAL :: horiz_cfl
2483 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mut
2485 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw
2487 REAL, INTENT(IN) :: dt
2488 REAL, INTENT(IN) :: rdx, rdy
2489 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, msfuy
2490 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfvx, msfvy
2492 REAL :: vert_cfl, cf_n, cf_d, maxdub, maxdeta
2494 INTEGER :: itf, jtf, i, j, k, maxi, maxj, maxk
2495 INTEGER :: some
2496 CHARACTER*512 :: temp
2498 CHARACTER (LEN=256) :: time_str
2499 CHARACTER (LEN=256) :: grid_str
2501 integer :: total
2502 REAL :: msfuxt , msfxffl
2504 !<DESCRIPTION>
2505 !
2506 ! w_damp computes a damping term for the vertical velocity when the
2507 ! vertical Courant number is too large. This was found to be preferable to
2508 ! decreasing the timestep or increasing the diffusion in real-data applications
2509 ! that produced potentially-unstable large vertical velocities because of
2510 ! unphysically large heating rates coming from the cumulus parameterization
2511 ! schemes run at moderately high resolutions (dx ~ O(10) km).
2512 !
2513 ! Additionally, w_damp returns the maximum cfl values due to vertical motion and
2514 ! horizontal motion. These values are returned via the max_vert_cfl and
2515 ! max_horiz_cfl variables. (Added by T. Hutchinson, WSI, 3/5/2007)
2516 !
2517 !</DESCRIPTION>
2519 itf=MIN(ite,ide-1)
2520 jtf=MIN(jte,jde-1)
2522 some = 0
2523 max_vert_cfl = 0.
2524 max_horiz_cfl = 0.
2525 total = 0
2527 IF(config_flags%map_proj == PROJ_CASSINI ) then
2528 msfxffl = 1.0/COS(config_flags%fft_filter_lat*degrad)
2529 END IF
2531 IF ( config_flags%w_damping == 1 ) THEN
2532 DO j = jts,jtf
2534 DO k = 2, kde-1
2535 DO i = its,itf
2536 #if 1
2537 IF(config_flags%map_proj == PROJ_CASSINI ) then
2538 msfuxt = MIN(msfux(i,j), msfxffl)
2539 ELSE
2540 msfuxt = msfux(i,j)
2541 END IF
2542 vert_cfl = abs(ww(i,k,j)/mut(i,j)*rdnw(k)*dt)
2544 IF ( vert_cfl > max_vert_cfl ) THEN
2545 max_vert_cfl = vert_cfl ; maxi = i ; maxj = j ; maxk = k
2546 maxdub = w(i,k,j) ; maxdeta = -1./rdnw(k)
2547 ENDIF
2549 horiz_cfl = max( abs(u(i,k,j) * rdx * msfuxt * dt), &
2550 abs(v(i,k,j) * rdy * msfvy(i,j) * dt) )
2551 if (horiz_cfl > max_horiz_cfl) then
2552 max_horiz_cfl = horiz_cfl
2553 endif
2555 if(vert_cfl .gt. w_beta)then
2556 #else
2557 ! restructure to get rid of divide
2558 !
2559 ! This had been used for efficiency, but with the addition of returning the cfl values,
2560 ! the old version (above) was reinstated. (T. Hutchinson, 3/5/2007)
2561 !
2562 cf_n = abs(ww(i,k,j)*rdnw(k)*dt)
2563 cf_d = abs(mut(i,j))
2564 if(cf_n .gt. cf_d*w_beta )then
2565 #endif
2566 !$OMP CRITICAL(big_step_utilities)
2567 WRITE(temp,*)i,j,k,' vert_cfl,w,d(eta)=',vert_cfl,w(i,k,j),-1./rdnw(k)
2568 CALL wrf_debug ( 100 , TRIM(temp) )
2569 !$OMP END CRITICAL(big_step_utilities)
2570 if ( vert_cfl > 2. ) some = some + 1
2571 rw_tend(i,k,j) = rw_tend(i,k,j)-sign(1.,w(i,k,j))*w_alpha*(vert_cfl-w_beta)*mut(i,j)
2572 endif
2573 END DO
2574 ENDDO
2575 ENDDO
2576 ELSE
2577 ! just print
2578 DO j = jts,jtf
2580 DO k = 2, kde-1
2581 DO i = its,itf
2583 #if 1
2584 IF(config_flags%map_proj == PROJ_CASSINI ) then
2585 msfuxt = MIN(msfux(i,j), msfxffl)
2586 ELSE
2587 msfuxt = msfux(i,j)
2588 END IF
2589 vert_cfl = abs(ww(i,k,j)/mut(i,j)*rdnw(k)*dt)
2591 IF ( vert_cfl > max_vert_cfl ) THEN
2592 max_vert_cfl = vert_cfl ; maxi = i ; maxj = j ; maxk = k
2593 maxdub = w(i,k,j) ; maxdeta = -1./rdnw(k)
2594 ENDIF
2596 horiz_cfl = max( abs(u(i,k,j) * rdx * msfuxt * dt), &
2597 abs(v(i,k,j) * rdy * msfvy(i,j) * dt) )
2599 if (horiz_cfl > max_horiz_cfl) then
2600 max_horiz_cfl = horiz_cfl
2601 endif
2603 if(vert_cfl .gt. w_beta)then
2604 #else
2605 ! restructure to get rid of divide
2606 !
2607 ! This had been used for efficiency, but with the addition of returning the cfl values,
2608 ! the old version (above) was reinstated. (T. Hutchinson, 3/5/2007)
2609 !
2610 cf_n = abs(ww(i,k,j)*rdnw(k)*dt)
2611 cf_d = abs(mut(i,j))
2612 if(cf_n .gt. cf_d*w_beta )then
2613 #endif
2614 !$OMP CRITICAL(big_step_utilities)
2615 WRITE(temp,*)i,j,k,' vert_cfl,w,d(eta)=',vert_cfl,w(i,k,j),-1./rdnw(k)
2616 CALL wrf_debug ( 100 , TRIM(temp) )
2617 !$OMP END CRITICAL(big_step_utilities)
2618 if ( vert_cfl > 2. ) some = some + 1
2619 endif
2620 END DO
2621 ENDDO
2622 ENDDO
2623 ENDIF
2624 IF ( some .GT. 0 ) THEN
2625 CALL get_current_time_string( time_str )
2626 CALL get_current_grid_name( grid_str )
2627 !$OMP CRITICAL(big_step_utilities)
2628 WRITE(wrf_err_message,*)some, &
2629 ' points exceeded cfl=2 in domain '//TRIM(grid_str)//' at time '//TRIM(time_str)//' hours'
2630 CALL wrf_debug ( 0 , TRIM(wrf_err_message) )
2631 WRITE(wrf_err_message,*)'MAX AT i,j,k: ',maxi,maxj,maxk,' vert_cfl,w,d(eta)=',max_vert_cfl, &
2632 maxdub,maxdeta
2633 CALL wrf_debug ( 0 , TRIM(wrf_err_message) )
2634 !$OMP END CRITICAL(big_step_utilities)
2635 ENDIF
2637 END SUBROUTINE w_damp
2639 !-------------------------------------------------------------------------------
2641 SUBROUTINE horizontal_diffusion ( name, field, tendency, mu, &
2642 config_flags, &
2643 msfux, msfuy, msfvx, msfvx_inv, &
2644 msfvy, msftx, msfty, &
2645 khdif, xkmhd, rdx, rdy, &
2646 ids, ide, jds, jde, kds, kde, &
2647 ims, ime, jms, jme, kms, kme, &
2648 its, ite, jts, jte, kts, kte )
2650 IMPLICIT NONE
2652 ! Input data
2654 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2656 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2657 ims, ime, jms, jme, kms, kme, &
2658 its, ite, jts, jte, kts, kte
2660 CHARACTER(LEN=1) , INTENT(IN ) :: name
2662 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, xkmhd
2664 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2666 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
2668 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
2669 msfuy, &
2670 msfvx, &
2671 msfvx_inv, &
2672 msfvy, &
2673 msftx, &
2674 msfty
2676 REAL , INTENT(IN ) :: rdx, &
2677 rdy, &
2678 khdif
2680 ! Local data
2682 INTEGER :: i, j, k, itf, jtf, ktf
2684 INTEGER :: i_start, i_end, j_start, j_end
2686 REAL :: mrdx, mkrdxm, mkrdxp, &
2687 mrdy, mkrdym, mkrdyp
2689 LOGICAL :: specified
2691 !<DESCRIPTION>
2692 !
2693 ! horizontal_diffusion computes the horizontal diffusion tendency
2694 ! on model horizontal coordinate surfaces.
2695 !
2696 !</DESCRIPTION>
2698 specified = .false.
2699 if(config_flags%specified .or. config_flags%nested) specified = .true.
2701 ktf=MIN(kte,kde-1)
2703 IF (name .EQ. 'u') THEN
2705 i_start = its
2706 i_end = ite
2707 j_start = jts
2708 j_end = MIN(jte,jde-1)
2710 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2711 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite)
2712 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2713 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2714 IF ( config_flags%periodic_x ) i_start = its
2715 IF ( config_flags%periodic_x ) i_end = ite
2718 DO j = j_start, j_end
2719 DO k=kts,ktf
2720 DO i = i_start, i_end
2722 ! The interior is grad: (m_x*d/dx), the exterior is div: (m_x*m_y*d/dx(/m_y))
2723 ! setting up different averagings of m^2 partial d/dX and m^2 partial d/dY
2725 mkrdxm=(msftx(i-1,j)/msfty(i-1,j))*mu(i-1,j)*xkmhd(i-1,k,j)*rdx
2726 mkrdxp=(msftx(i,j)/msfty(i,j))*mu(i,j)*xkmhd(i,k,j)*rdx
2727 mrdx=msfux(i,j)*msfuy(i,j)*rdx
2728 mkrdym=( (msfuy(i,j)+msfuy(i,j-1))/(msfux(i,j)+msfux(i,j-1)) )* &
2729 0.25*(mu(i,j)+mu(i,j-1)+mu(i-1,j-1)+mu(i-1,j))* &
2730 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i-1,k,j-1)+xkmhd(i-1,k,j))*rdy
2731 mkrdyp=( (msfuy(i,j)+msfuy(i,j+1))/(msfux(i,j)+msfux(i,j+1)) )* &
2732 0.25*(mu(i,j)+mu(i,j+1)+mu(i-1,j+1)+mu(i-1,j))* &
2733 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j+1)+xkmhd(i-1,k,j+1)+xkmhd(i-1,k,j))*rdy
2734 ! need to do four-corners (t) for diffusion coefficient as there are
2735 ! no values at u,v points
2736 ! msfuy - has to be y as part of d/dY
2737 ! has to be u as we're at a u point
2738 mrdy=msfux(i,j)*msfuy(i,j)*rdy
2740 ! correctly averaged version of rho~ * m^2 *
2741 ! [partial d/dX(partial du^/dX) + partial d/dY(partial du^/dY)]
2742 tendency(i,k,j)=tendency(i,k,j)+( &
2743 mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) &
2744 -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2745 +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) &
2746 -mkrdym*(field(i,k,j )-field(i,k,j-1))))
2747 ENDDO
2748 ENDDO
2749 ENDDO
2751 ELSE IF (name .EQ. 'v')THEN
2753 i_start = its
2754 i_end = MIN(ite,ide-1)
2755 j_start = jts
2756 j_end = jte
2758 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2759 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2760 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2761 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte)
2762 IF ( config_flags%periodic_x ) i_start = its
2763 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2764 IF ( config_flags%polar ) j_start = MAX(jds+1,jts)
2765 IF ( config_flags%polar ) j_end = MIN(jde-1,jte)
2767 DO j = j_start, j_end
2768 DO k=kts,ktf
2769 DO i = i_start, i_end
2771 mkrdxm=( (msfvx(i,j)+msfvx(i-1,j))/(msfvy(i,j)+msfvy(i-1,j)) )* &
2772 0.25*(mu(i,j)+mu(i,j-1)+mu(i-1,j-1)+mu(i-1,j))* &
2773 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i-1,k,j-1)+xkmhd(i-1,k,j))*rdx
2774 mkrdxp=( (msfvx(i,j)+msfvx(i+1,j))/(msfvy(i,j)+msfvy(i+1,j)) )* &
2775 0.25*(mu(i,j)+mu(i,j-1)+mu(i+1,j-1)+mu(i+1,j))* &
2776 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i+1,k,j-1)+xkmhd(i+1,k,j))*rdx
2777 mrdx=msfvx(i,j)*msfvy(i,j)*rdx
2778 mkrdym=(msfty(i,j-1)/msftx(i,j-1))*xkmhd(i,k,j-1)*rdy
2779 mkrdyp=(msfty(i,j)/msftx(i,j))*xkmhd(i,k,j)*rdy
2780 mrdy=msfvx(i,j)*msfvy(i,j)*rdy
2782 tendency(i,k,j)=tendency(i,k,j)+( &
2783 mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) &
2784 -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2785 +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) &
2786 -mkrdym*(field(i,k,j )-field(i,k,j-1))))
2787 ENDDO
2788 ENDDO
2789 ENDDO
2791 ELSE IF (name .EQ. 'w')THEN
2793 i_start = its
2794 i_end = MIN(ite,ide-1)
2795 j_start = jts
2796 j_end = MIN(jte,jde-1)
2798 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2799 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2800 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2801 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2802 IF ( config_flags%periodic_x ) i_start = its
2803 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2805 DO j = j_start, j_end
2806 DO k=kts+1,ktf
2807 DO i = i_start, i_end
2809 mkrdxm=(msfux(i,j)/msfuy(i,j))* &
2810 0.25*(mu(i,j)+mu(i-1,j)+mu(i,j)+mu(i-1,j))* &
2811 0.25*(xkmhd(i,k,j)+xkmhd(i-1,k,j)+xkmhd(i,k-1,j)+xkmhd(i-1,k-1,j))*rdx
2812 mkrdxp=(msfux(i+1,j)/msfuy(i+1,j))* &
2813 0.25*(mu(i+1,j)+mu(i,j)+mu(i+1,j)+mu(i,j))* &
2814 0.25*(xkmhd(i+1,k,j)+xkmhd(i,k,j)+xkmhd(i+1,k-1,j)+xkmhd(i,k-1,j))*rdx
2815 mrdx=msftx(i,j)*msfty(i,j)*rdx
2816 ! mkrdym=(msfvy(i,j)/msfvx(i,j))* &
2817 mkrdym=(msfvy(i,j)*msfvx_inv(i,j))* &
2818 0.25*(mu(i,j)+mu(i,j-1)+mu(i,j)+mu(i,j-1))* &
2819 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i,k-1,j)+xkmhd(i,k-1,j-1))*rdy
2820 ! mkrdyp=(msfvy(i,j+1)/msfvx(i,j+1))* &
2821 mkrdyp=(msfvy(i,j+1)*msfvx_inv(i,j+1))* &
2822 0.25*(mu(i,j+1)+mu(i,j)+mu(i,j+1)+mu(i,j))* &
2823 0.25*(xkmhd(i,k,j+1)+xkmhd(i,k,j)+xkmhd(i,k-1,j+1)+xkmhd(i,k-1,j))*rdy
2824 mrdy=msftx(i,j)*msfty(i,j)*rdy
2826 tendency(i,k,j)=tendency(i,k,j)+( &
2827 mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) &
2828 -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2829 +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) &
2830 -mkrdym*(field(i,k,j )-field(i,k,j-1))))
2831 ENDDO
2832 ENDDO
2833 ENDDO
2835 ELSE
2838 i_start = its
2839 i_end = MIN(ite,ide-1)
2840 j_start = jts
2841 j_end = MIN(jte,jde-1)
2843 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2844 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2845 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2846 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2847 IF ( config_flags%periodic_x ) i_start = its
2848 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2850 DO j = j_start, j_end
2851 DO k=kts,ktf
2852 DO i = i_start, i_end
2854 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
2855 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
2856 mrdx=msftx(i,j)*msfty(i,j)*rdx
2857 ! 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
2858 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
2859 ! 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
2860 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
2861 mrdy=msftx(i,j)*msfty(i,j)*rdy
2863 tendency(i,k,j)=tendency(i,k,j)+( &
2864 mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) &
2865 -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2866 +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) &
2867 -mkrdym*(field(i,k,j )-field(i,k,j-1))))
2868 ENDDO
2869 ENDDO
2870 ENDDO
2872 ENDIF
2874 END SUBROUTINE horizontal_diffusion
2876 !-----------------------------------------------------------------------------------------
2878 SUBROUTINE horizontal_diffusion_3dmp ( name, field, tendency, mu, &
2879 config_flags, base_3d, &
2880 msfux, msfuy, msfvx, msfvx_inv, &
2881 msfvy, msftx, msfty, &
2882 khdif, xkmhd, rdx, rdy, &
2883 ids, ide, jds, jde, kds, kde, &
2884 ims, ime, jms, jme, kms, kme, &
2885 its, ite, jts, jte, kts, kte )
2887 IMPLICIT NONE
2889 ! Input data
2891 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2893 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2894 ims, ime, jms, jme, kms, kme, &
2895 its, ite, jts, jte, kts, kte
2897 CHARACTER(LEN=1) , INTENT(IN ) :: name
2899 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, &
2900 xkmhd, &
2901 base_3d
2903 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2905 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
2907 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
2908 msfuy, &
2909 msfvx, &
2910 msfvx_inv, &
2911 msfvy, &
2912 msftx, &
2913 msfty
2915 REAL , INTENT(IN ) :: rdx, &
2916 rdy, &
2917 khdif
2919 ! Local data
2921 INTEGER :: i, j, k, itf, jtf, ktf
2923 INTEGER :: i_start, i_end, j_start, j_end
2925 REAL :: mrdx, mkrdxm, mkrdxp, &
2926 mrdy, mkrdym, mkrdyp
2928 LOGICAL :: specified
2930 !<DESCRIPTION>
2931 !
2932 ! horizontal_diffusion_3dmp computes the horizontal diffusion tendency
2933 ! on model horizontal coordinate surfaces. This routine computes diffusion
2934 ! a perturbation scalar (field-base_3d).
2935 !
2936 !</DESCRIPTION>
2938 specified = .false.
2939 if(config_flags%specified .or. config_flags%nested) specified = .true.
2941 ktf=MIN(kte,kde-1)
2943 i_start = its
2944 i_end = MIN(ite,ide-1)
2945 j_start = jts
2946 j_end = MIN(jte,jde-1)
2948 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2949 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2950 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2951 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2952 IF ( config_flags%periodic_x ) i_start = its
2953 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2955 DO j = j_start, j_end
2956 DO k=kts,ktf
2957 DO i = i_start, i_end
2959 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
2960 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
2961 mrdx=msftx(i,j)*msfty(i,j)*rdx
2962 ! 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
2963 ! 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
2964 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
2965 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
2966 mrdy=msftx(i,j)*msfty(i,j)*rdy
2968 tendency(i,k,j)=tendency(i,k,j)+( &
2969 mrdx*( mkrdxp*( field(i+1,k,j) -field(i ,k,j) &
2970 -base_3d(i+1,k,j)+base_3d(i ,k,j) ) &
2971 -mkrdxm*( field(i ,k,j) -field(i-1,k,j) &
2972 -base_3d(i ,k,j)+base_3d(i-1,k,j) ) ) &
2973 +mrdy*( mkrdyp*( field(i,k,j+1) -field(i,k,j ) &
2974 -base_3d(i,k,j+1)+base_3d(i,k,j ) ) &
2975 -mkrdym*( field(i,k,j ) -field(i,k,j-1) &
2976 -base_3d(i,k,j )+base_3d(i,k,j-1) ) ) &
2977 )
2978 ENDDO
2979 ENDDO
2980 ENDDO
2982 END SUBROUTINE horizontal_diffusion_3dmp
2984 !-----------------------------------------------------------------------------------------
2986 SUBROUTINE vertical_diffusion ( name, field, tendency, &
2987 config_flags, &
2988 alt, mut, rdn, rdnw, kvdif, &
2989 ids, ide, jds, jde, kds, kde, &
2990 ims, ime, jms, jme, kms, kme, &
2991 its, ite, jts, jte, kts, kte )
2994 IMPLICIT NONE
2996 ! Input data
2998 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3000 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3001 ims, ime, jms, jme, kms, kme, &
3002 its, ite, jts, jte, kts, kte
3004 CHARACTER(LEN=1) , INTENT(IN ) :: name
3006 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3007 INTENT(IN ) :: field, &
3008 alt
3010 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3012 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
3014 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw
3016 REAL , INTENT(IN ) :: kvdif
3018 ! Local data
3020 INTEGER :: i, j, k, itf, jtf, ktf
3021 INTEGER :: i_start, i_end, j_start, j_end
3023 REAL , DIMENSION(its:ite, jts:jte) :: vfluxm, vfluxp, zz
3024 REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3026 REAL :: rdz
3028 LOGICAL :: specified
3030 !<DESCRIPTION>
3031 !
3032 ! vertical_diffusion
3033 ! computes vertical diffusion tendency.
3034 !
3035 !</DESCRIPTION>
3037 specified = .false.
3038 if(config_flags%specified .or. config_flags%nested) specified = .true.
3040 ktf=MIN(kte,kde-1)
3042 IF (name .EQ. 'w')THEN
3045 i_start = its
3046 i_end = MIN(ite,ide-1)
3047 j_start = jts
3048 j_end = MIN(jte,jde-1)
3050 j_loop_w : DO j = j_start, j_end
3052 DO k=kts,ktf-1
3053 DO i = i_start, i_end
3054 vflux(i,k)= (kvdif/alt(i,k,j))*rdnw(k)*(field(i,k+1,j)-field(i,k,j))
3055 ENDDO
3056 ENDDO
3058 DO i = i_start, i_end
3059 vflux(i,ktf)=0.
3060 ENDDO
3062 DO k=kts+1,ktf
3063 DO i = i_start, i_end
3064 tendency(i,k,j)=tendency(i,k,j) &
3065 +rdn(k)*g*g/mut(i,j)/(0.5*(alt(i,k,j)+alt(i,k-1,j))) &
3066 *(vflux(i,k)-vflux(i,k-1))
3067 ENDDO
3068 ENDDO
3070 ENDDO j_loop_w
3072 ELSE IF(name .EQ. 'm')THEN
3074 i_start = its
3075 i_end = MIN(ite,ide-1)
3076 j_start = jts
3077 j_end = MIN(jte,jde-1)
3079 j_loop_s : DO j = j_start, j_end
3081 DO k=kts,ktf-1
3082 DO i = i_start, i_end
3083 vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j))) &
3084 *(field(i,k+1,j)-field(i,k,j))
3085 ENDDO
3086 ENDDO
3088 DO i = i_start, i_end
3089 vflux(i,0)=vflux(i,1)
3090 ENDDO
3092 DO i = i_start, i_end
3093 vflux(i,ktf)=0.
3094 ENDDO
3096 DO k=kts,ktf
3097 DO i = i_start, i_end
3098 tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) &
3099 *rdnw(k)*(vflux(i,k)-vflux(i,k-1))
3100 ENDDO
3101 ENDDO
3103 ENDDO j_loop_s
3105 ENDIF
3107 END SUBROUTINE vertical_diffusion
3110 !-------------------------------------------------------------------------------
3112 SUBROUTINE vertical_diffusion_mp ( field, tendency, config_flags, &
3113 base, &
3114 alt, mut, rdn, rdnw, kvdif, &
3115 ids, ide, jds, jde, kds, kde, &
3116 ims, ime, jms, jme, kms, kme, &
3117 its, ite, jts, jte, kts, kte )
3120 IMPLICIT NONE
3122 ! Input data
3124 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3126 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3127 ims, ime, jms, jme, kms, kme, &
3128 its, ite, jts, jte, kts, kte
3130 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3131 INTENT(IN ) :: field, &
3132 alt
3134 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3136 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
3138 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, &
3139 rdnw, &
3140 base
3142 REAL , INTENT(IN ) :: kvdif
3144 ! Local data
3146 INTEGER :: i, j, k, itf, jtf, ktf
3147 INTEGER :: i_start, i_end, j_start, j_end
3149 REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3151 REAL :: rdz
3153 LOGICAL :: specified
3155 !<DESCRIPTION>
3156 !
3157 ! vertical_diffusion_mp
3158 ! computes vertical diffusion tendency of a perturbation variable
3159 ! (field-base). Note that base as a 1D (k) field.
3160 !
3161 !</DESCRIPTION>
3163 specified = .false.
3164 if(config_flags%specified .or. config_flags%nested) specified = .true.
3166 ktf=MIN(kte,kde-1)
3168 i_start = its
3169 i_end = MIN(ite,ide-1)
3170 j_start = jts
3171 j_end = MIN(jte,jde-1)
3173 j_loop_s : DO j = j_start, j_end
3175 DO k=kts,ktf-1
3176 DO i = i_start, i_end
3177 vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j))) &
3178 *(field(i,k+1,j)-field(i,k,j)-base(k+1)+base(k))
3179 ENDDO
3180 ENDDO
3182 DO i = i_start, i_end
3183 vflux(i,0)=vflux(i,1)
3184 ENDDO
3186 DO i = i_start, i_end
3187 vflux(i,ktf)=0.
3188 ENDDO
3190 DO k=kts,ktf
3191 DO i = i_start, i_end
3192 tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) &
3193 *rdnw(k)*(vflux(i,k)-vflux(i,k-1))
3194 ENDDO
3195 ENDDO
3197 ENDDO j_loop_s
3199 END SUBROUTINE vertical_diffusion_mp
3202 !-------------------------------------------------------------------------------
3204 SUBROUTINE vertical_diffusion_3dmp ( field, tendency, config_flags, &
3205 base_3d, &
3206 alt, mut, rdn, rdnw, kvdif, &
3207 ids, ide, jds, jde, kds, kde, &
3208 ims, ime, jms, jme, kms, kme, &
3209 its, ite, jts, jte, kts, kte )
3212 IMPLICIT NONE
3214 ! Input data
3216 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3218 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3219 ims, ime, jms, jme, kms, kme, &
3220 its, ite, jts, jte, kts, kte
3222 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3223 INTENT(IN ) :: field, &
3224 alt, &
3225 base_3d
3227 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3229 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
3231 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, &
3232 rdnw
3234 REAL , INTENT(IN ) :: kvdif
3236 ! Local data
3238 INTEGER :: i, j, k, itf, jtf, ktf
3239 INTEGER :: i_start, i_end, j_start, j_end
3241 REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3243 REAL :: rdz
3245 LOGICAL :: specified
3247 !<DESCRIPTION>
3248 !
3249 ! vertical_diffusion_3dmp
3250 ! computes vertical diffusion tendency of a perturbation variable
3251 ! (field-base_3d).
3252 !
3253 !</DESCRIPTION>
3255 specified = .false.
3256 if(config_flags%specified .or. config_flags%nested) specified = .true.
3258 ktf=MIN(kte,kde-1)
3260 i_start = its
3261 i_end = MIN(ite,ide-1)
3262 j_start = jts
3263 j_end = MIN(jte,jde-1)
3265 j_loop_s : DO j = j_start, j_end
3267 DO k=kts,ktf-1
3268 DO i = i_start, i_end
3269 vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j))) &
3270 *( field(i,k+1,j) -field(i,k,j) &
3271 -base_3d(i,k+1,j)+base_3d(i,k,j) )
3272 ENDDO
3273 ENDDO
3275 DO i = i_start, i_end
3276 vflux(i,0)=vflux(i,1)
3277 ENDDO
3279 DO i = i_start, i_end
3280 vflux(i,ktf)=0.
3281 ENDDO
3283 DO k=kts,ktf
3284 DO i = i_start, i_end
3285 tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) &
3286 *rdnw(k)*(vflux(i,k)-vflux(i,k-1))
3287 ENDDO
3288 ENDDO
3290 ENDDO j_loop_s
3292 END SUBROUTINE vertical_diffusion_3dmp
3295 !-------------------------------------------------------------------------------
3298 SUBROUTINE vertical_diffusion_u ( field, tendency, &
3299 config_flags, u_base, &
3300 alt, muu, rdn, rdnw, kvdif, &
3301 ids, ide, jds, jde, kds, kde, &
3302 ims, ime, jms, jme, kms, kme, &
3303 its, ite, jts, jte, kts, kte )
3306 IMPLICIT NONE
3308 ! Input data
3310 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3312 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3313 ims, ime, jms, jme, kms, kme, &
3314 its, ite, jts, jte, kts, kte
3316 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3317 INTENT(IN ) :: field, &
3318 alt
3320 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3322 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muu
3324 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw, u_base
3326 REAL , INTENT(IN ) :: kvdif
3328 ! Local data
3330 INTEGER :: i, j, k, itf, jtf, ktf
3331 INTEGER :: i_start, i_end, j_start, j_end
3333 REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3335 REAL :: rdz, zz
3337 LOGICAL :: specified
3339 !<DESCRIPTION>
3340 !
3341 ! vertical_diffusion_u computes vertical diffusion tendency for
3342 ! the u momentum equation. This routine assumes a constant eddy
3343 ! viscosity kvdif.
3344 !
3345 !</DESCRIPTION>
3347 specified = .false.
3348 if(config_flags%specified .or. config_flags%nested) specified = .true.
3350 ktf=MIN(kte,kde-1)
3352 i_start = its
3353 i_end = ite
3354 j_start = jts
3355 j_end = MIN(jte,jde-1)
3357 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
3358 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite)
3359 IF ( config_flags%periodic_x ) i_start = its
3360 IF ( config_flags%periodic_x ) i_end = ite
3363 j_loop_u : DO j = j_start, j_end
3365 DO k=kts,ktf-1
3366 DO i = i_start, i_end
3367 vflux(i,k)=kvdif*rdn(k+1)/(0.25*( alt(i ,k ,j) &
3368 +alt(i-1,k ,j) &
3369 +alt(i ,k+1,j) &
3370 +alt(i-1,k+1,j) ) ) &
3371 *(field(i,k+1,j)-field(i,k,j) &
3372 -u_base(k+1) +u_base(k) )
3373 ENDDO
3374 ENDDO
3376 DO i = i_start, i_end
3377 vflux(i,0)=vflux(i,1)
3378 ENDDO
3380 DO i = i_start, i_end
3381 vflux(i,ktf)=0.
3382 ENDDO
3384 DO k=kts,ktf-1
3385 DO i = i_start, i_end
3386 tendency(i,k,j)=tendency(i,k,j)+ &
3387 g*g*rdnw(k)/muu(i,j)/(0.5*(alt(i-1,k,j)+alt(i,k,j)))* &
3388 (vflux(i,k)-vflux(i,k-1))
3389 ENDDO
3390 ENDDO
3392 ENDDO j_loop_u
3394 END SUBROUTINE vertical_diffusion_u
3396 !-------------------------------------------------------------------------------
3399 SUBROUTINE vertical_diffusion_v ( field, tendency, &
3400 config_flags, v_base, &
3401 alt, muv, rdn, rdnw, kvdif, &
3402 ids, ide, jds, jde, kds, kde, &
3403 ims, ime, jms, jme, kms, kme, &
3404 its, ite, jts, jte, kts, kte )
3407 IMPLICIT NONE
3409 ! Input data
3411 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3413 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3414 ims, ime, jms, jme, kms, kme, &
3415 its, ite, jts, jte, kts, kte
3417 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3418 INTENT(IN ) :: field, &
3419 alt
3420 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw, v_base
3422 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3424 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muv
3426 REAL , INTENT(IN ) :: kvdif
3428 ! Local data
3430 INTEGER :: i, j, k, itf, jtf, ktf, jm1
3431 INTEGER :: i_start, i_end, j_start, j_end
3433 REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3435 REAL :: rdz, zz
3437 LOGICAL :: specified
3439 !<DESCRIPTION>
3440 !
3441 ! vertical_diffusion_v computes vertical diffusion tendency for
3442 ! the v momentum equation. This routine assumes a constant eddy
3443 ! viscosity kvdif.
3444 !
3445 !</DESCRIPTION>
3447 specified = .false.
3448 if(config_flags%specified .or. config_flags%nested) specified = .true.
3450 ktf=MIN(kte,kde-1)
3452 i_start = its
3453 i_end = MIN(ite,ide-1)
3454 j_start = jts
3455 j_end = MIN(jte,jde-1)
3457 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
3458 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte)
3460 j_loop_v : DO j = j_start, j_end
3461 ! jm1 = max(j-1,1)
3462 jm1 = j-1
3464 DO k=kts,ktf-1
3465 DO i = i_start, i_end
3466 vflux(i,k)=kvdif*rdn(k+1)/(0.25*( alt(i,k ,j ) &
3467 +alt(i,k ,jm1) &
3468 +alt(i,k+1,j ) &
3469 +alt(i,k+1,jm1) ) ) &
3470 *(field(i,k+1,j)-field(i,k,j) &
3471 -v_base(k+1) +v_base(k) )
3472 ENDDO
3473 ENDDO
3475 DO i = i_start, i_end
3476 vflux(i,0)=vflux(i,1)
3477 ENDDO
3479 DO i = i_start, i_end
3480 vflux(i,ktf)=0.
3481 ENDDO
3483 DO k=kts,ktf-1
3484 DO i = i_start, i_end
3485 tendency(i,k,j)=tendency(i,k,j)+ &
3486 g*g*rdnw(k)/muv(i,j)/(0.5*(alt(i,k,jm1)+alt(i,k,j)))* &
3487 (vflux(i,k)-vflux(i,k-1))
3488 ENDDO
3489 ENDDO
3491 ENDDO j_loop_v
3493 END SUBROUTINE vertical_diffusion_v
3495 !*************** end new mass coordinate routines
3497 !-------------------------------------------------------------------------------
3499 SUBROUTINE calculate_full ( rfield, rfieldb, rfieldp, &
3500 ids, ide, jds, jde, kds, kde, &
3501 ims, ime, jms, jme, kms, kme, &
3502 its, ite, jts, jte, kts, kte )
3504 IMPLICIT NONE
3506 ! Input data
3508 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3509 ims, ime, jms, jme, kms, kme, &
3510 its, ite, jts, jte, kts, kte
3512 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rfieldb, &
3513 rfieldp
3515 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: rfield
3517 ! Local indices.
3519 INTEGER :: i, j, k, itf, jtf, ktf
3521 !<DESCRIPTION>
3522 !
3523 ! calculate_full
3524 ! calculates full 3D field from pertubation and base field.
3525 !
3526 !</DESCRIPTION>
3528 itf=MIN(ite,ide-1)
3529 jtf=MIN(jte,jde-1)
3530 ktf=MIN(kte,kde-1)
3532 DO j=jts,jtf
3533 DO k=kts,ktf
3534 DO i=its,itf
3535 rfield(i,k,j)=rfieldb(i,k,j)+rfieldp(i,k,j)
3536 ENDDO
3537 ENDDO
3538 ENDDO
3540 END SUBROUTINE calculate_full
3542 !------------------------------------------------------------------------------
3544 SUBROUTINE coriolis ( ru, rv, rw, ru_tend, rv_tend, rw_tend, &
3545 config_flags, &
3546 msftx, msfty, msfux, msfuy, &
3547 msfvx, msfvy, &
3548 f, e, sina, cosa, fzm, fzp, &
3549 ids, ide, jds, jde, kds, kde, &
3550 ims, ime, jms, jme, kms, kme, &
3551 its, ite, jts, jte, kts, kte )
3553 IMPLICIT NONE
3555 ! Input data
3557 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
3559 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3560 ims, ime, jms, jme, kms, kme, &
3561 its, ite, jts, jte, kts, kte
3563 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: ru_tend, &
3564 rv_tend, &
3565 rw_tend
3566 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: ru, &
3567 rv, &
3568 rw
3570 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
3571 msfuy, &
3572 msfvx, &
3573 msfvy, &
3574 msftx, &
3575 msfty
3577 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: f, &
3578 e, &
3579 sina, &
3580 cosa
3582 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
3583 fzp
3585 ! Local indices.
3587 INTEGER :: i, j , k, ktf
3588 INTEGER :: i_start, i_end, j_start, j_end
3590 LOGICAL :: specified
3592 !<DESCRIPTION>
3593 !
3594 ! coriolis calculates the large timestep tendency terms in the
3595 ! u, v, and w momentum equations arise from the coriolis force.
3596 !
3597 !</DESCRIPTION>
3599 specified = .false.
3600 if(config_flags%specified .or. config_flags%nested) specified = .true.
3602 ktf=MIN(kte,kde-1)
3604 ! coriolis for u-momentum equation
3606 ! Notes on map scale factor
3607 ! cosa, sina are related to rotating the coordinate frame if desired
3608 ! generally sina=0, cosa=1
3609 ! ADT eqn 44, RHS terms 6 and 7: -2 mu w omega cos(lat)/my
3610 ! + 2 mu v omega sin(lat)/my
3611 ! Define f=2 omega sin(lat), e=2 omega cos(lat)
3612 ! => terms are: -e mu w / my + f mu v / my
3613 ! rv = mu v / mx ; rw = mu w / my
3614 ! => terms are: -e rw + f rv *mx / my
3616 i_start = its
3617 i_end = ite
3618 IF ( config_flags%open_xs .or. specified .or. &
3619 config_flags%nested) i_start = MAX(ids+1,its)
3620 IF ( config_flags%open_xe .or. specified .or. &
3621 config_flags%nested) i_end = MIN(ide-1,ite)
3622 IF ( config_flags%periodic_x ) i_start = its
3623 IF ( config_flags%periodic_x ) i_end = ite
3625 DO j = jts, MIN(jte,jde-1)
3627 DO k=kts,ktf
3628 DO i = i_start, i_end
3630 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)) &
3631 *0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
3632 - 0.5*(e(i,j)+e(i-1,j))*0.5*(cosa(i,j)+cosa(i-1,j)) &
3633 *0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
3635 ENDDO
3636 ENDDO
3638 ! boundary loops for coriolis not needed for open bdy (commented out 20100611 JD)
3639 ! IF ( (config_flags%open_xs) .and. (its == ids) ) THEN
3641 ! DO k=kts,ktf
3642 !
3643 ! ru_tend(its,k,j)=ru_tend(its,k,j) + (msfux(its,j)/msfuy(its,j))*0.5*(f(its,j)+f(its,j)) &
3644 ! *0.25*(rv(its,k,j+1)+rv(its,k,j+1)+rv(its,k,j)+rv(its,k,j)) &
3645 ! - 0.5*(e(its,j)+e(its,j))*0.5*(cosa(its,j)+cosa(its,j)) &
3646 ! *0.25*(rw(its,k+1,j)+rw(its,k,j)+rw(its,k+1,j)+rw(its,k,j))
3648 ! ENDDO
3650 ! ENDIF
3652 ! IF ( (config_flags%open_xe) .and. (ite == ide) ) THEN
3654 ! DO k=kts,ktf
3655 !
3656 ! 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)) &
3657 ! *0.25*(rv(ite-1,k,j+1)+rv(ite-1,k,j+1)+rv(ite-1,k,j)+rv(ite-1,k,j)) &
3658 ! - 0.5*(e(ite-1,j)+e(ite-1,j))*0.5*(cosa(ite-1,j)+cosa(ite-1,j)) &
3659 ! *0.25*(rw(ite-1,k+1,j)+rw(ite-1,k,j)+rw(ite-1,k+1,j)+rw(ite-1,k,j))
3661 ! ENDDO
3663 ! ENDIF
3665 ENDDO
3667 ! coriolis term for v-momentum equation
3669 ! Notes on map scale factors
3670 ! ADT eqn 45, RHS terms 6 and 6b [0 for sina=0]: -2 mu u omega sin(lat)/mx + ?
3671 ! Define f=2 omega sin(lat), e=2 omega cos(lat)
3672 ! => terms are: -f mu u / mx
3673 ! ru = mu u / my ; rw = mu w / my
3674 ! => terms are: -f ru *my / mx + ?
3676 j_start = jts
3677 j_end = jte
3679 IF ( config_flags%open_ys .or. specified .or. &
3680 config_flags%nested .or. config_flags%polar) j_start = MAX(jds+1,jts)
3681 IF ( config_flags%open_ye .or. specified .or. &
3682 config_flags%nested .or. config_flags%polar) j_end = MIN(jde-1,jte)
3684 ! boundary loops for coriolis not needed for open bdy (commented out 20100611 JD)
3685 ! IF ( (config_flags%open_ys) .and. (jts == jds) ) THEN
3687 ! DO k=kts,ktf
3688 ! DO i=its,MIN(ide-1,ite)
3689 !
3690 ! rv_tend(i,k,jts)=rv_tend(i,k,jts) - (msfvy(i,jts)/msfvx(i,jts))*0.5*(f(i,jts)+f(i,jts)) &
3691 ! *0.25*(ru(i,k,jts)+ru(i+1,k,jts)+ru(i,k,jts)+ru(i+1,k,jts)) &
3692 ! + (msfvy(i,jts)/msfvx(i,jts))*0.5*(e(i,jts)+e(i,jts))*0.5*(sina(i,jts)+sina(i,jts)) &
3693 ! *0.25*(rw(i,k+1,jts)+rw(i,k,jts)+rw(i,k+1,jts)+rw(i,k,jts))
3695 ! ENDDO
3696 ! ENDDO
3698 ! ENDIF
3700 DO j=j_start, j_end
3701 DO k=kts,ktf
3702 DO i=its,MIN(ide-1,ite)
3704 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)) &
3705 *0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
3706 + (msfvy(i,j)/msfvx(i,j))*0.5*(e(i,j)+e(i,j-1))*0.5*(sina(i,j)+sina(i,j-1)) &
3707 *0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
3709 ENDDO
3710 ENDDO
3711 ENDDO
3714 ! boundary loops for coriolis not needed for open bdy (commented out 20100611 JD)
3715 ! IF ( (config_flags%open_ye) .and. (jte == jde) ) THEN
3717 ! DO k=kts,ktf
3718 ! DO i=its,MIN(ide-1,ite)
3719 !
3720 ! 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)) &
3721 ! *0.25*(ru(i,k,jte-1)+ru(i+1,k,jte-1)+ru(i,k,jte-1)+ru(i+1,k,jte-1)) &
3722 ! + (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)) &
3723 ! *0.25*(rw(i,k+1,jte-1)+rw(i,k,jte-1)+rw(i,k+1,jte-1)+rw(i,k,jte-1))
3725 ! ENDDO
3726 ! ENDDO
3728 ! ENDIF
3730 ! coriolis term for w-mometum
3732 ! Notes on map scale factors
3733 ! ADT eqn 46/my, RHS terms 5 and 5b [0 for sina=0]: 2 mu u omega cos(lat)/my +?
3734 ! Define e=2 omega cos(lat)
3735 ! => terms are: e mu u / my + ???
3736 ! ru = mu u / my ; ru = mu v / mx
3737 ! => terms are: e ru + ???
3739 DO j=jts,MIN(jte, jde-1)
3740 DO k=kts+1,ktf
3741 DO i=its,MIN(ite, ide-1)
3743 rw_tend(i,k,j)=rw_tend(i,k,j) + e(i,j)* &
3744 (cosa(i,j)*0.5*(fzm(k)*(ru(i,k,j)+ru(i+1,k,j)) &
3745 +fzp(k)*(ru(i,k-1,j)+ru(i+1,k-1,j))) &
3746 -(msftx(i,j)/msfty(i,j))* &
3747 sina(i,j)*0.5*(fzm(k)*(rv(i,k,j)+rv(i,k,j+1)) &
3748 +fzp(k)*(rv(i,k-1,j)+rv(i,k-1,j+1))))
3750 ENDDO
3751 ENDDO
3752 ENDDO
3754 END SUBROUTINE coriolis
3756 !------------------------------------------------------------------------------
3758 SUBROUTINE perturbation_coriolis ( ru_in, rv_in, rw, ru_tend, rv_tend, rw_tend, &
3759 config_flags, &
3760 u_base, v_base, z_base, &
3761 muu, muv, phb, ph, &
3762 msftx, msfty, msfux, msfuy, msfvx, msfvy, &
3763 f, e, sina, cosa, fzm, fzp, &
3764 ids, ide, jds, jde, kds, kde, &
3765 ims, ime, jms, jme, kms, kme, &
3766 its, ite, jts, jte, kts, kte )
3768 IMPLICIT NONE
3770 ! Input data
3772 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
3774 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3775 ims, ime, jms, jme, kms, kme, &
3776 its, ite, jts, jte, kts, kte
3778 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: ru_tend, &
3779 rv_tend, &
3780 rw_tend
3781 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: ru_in, &
3782 rv_in, &
3783 rw, &
3784 ph, &
3785 phb
3788 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
3789 msfuy, &
3790 msfvx, &
3791 msfvy, &
3792 msftx, &
3793 msfty
3795 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: f, &
3796 e, &
3797 sina, &
3798 cosa
3800 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muu, &
3801 muv
3804 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
3805 fzp
3807 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: u_base, &
3808 v_base, &
3809 z_base
3811 ! Local storage
3813 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) :: ru, &
3814 rv
3816 REAL :: z_at_u, z_at_v, wkp1, wk, wkm1
3818 ! Local indices.
3820 INTEGER :: i, j , k, ktf
3821 INTEGER :: i_start, i_end, j_start, j_end
3823 LOGICAL :: specified
3825 !<DESCRIPTION>
3826 !
3827 ! perturbation_coriolis calculates the large timestep tendency terms in the
3828 ! u, v, and w momentum equations arise from the coriolis force. This version
3829 ! subtracts off the horizontal velocities from the initial sounding when
3830 ! computing the forcing terms, hence "perturbation" coriolis.
3831 !
3832 !</DESCRIPTION>
3834 specified = .false.
3835 if(config_flags%specified .or. config_flags%nested) specified = .true.
3837 ktf=MIN(kte,kde-1)
3839 ! coriolis for u-momentum equation
3841 i_start = its
3842 i_end = ite
3843 IF ( config_flags%open_xs .or. specified .or. &
3844 config_flags%nested) i_start = MAX(ids+1,its)
3845 IF ( config_flags%open_xe .or. specified .or. &
3846 config_flags%nested) i_end = MIN(ide-1,ite)
3847 IF ( config_flags%periodic_x ) i_start = its
3848 IF ( config_flags%periodic_x ) i_end = ite
3850 ! compute perturbation mu*v for use in u momentum equation
3852 DO j = jts, MIN(jte,jde-1)+1
3853 DO k=kts+1,ktf-1
3854 DO i = i_start-1, i_end
3855 z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
3856 +phb(i,k,j-1)+phb(i,k+1,j-1) &
3857 +ph(i,k,j )+ph(i,k+1,j ) &
3858 +ph(i,k,j-1)+ph(i,k+1,j-1))/g
3859 wkp1 = min(1.,max(0.,z_at_v-z_base(k))/(z_base(k+1)-z_base(k)))
3860 wkm1 = min(1.,max(0.,z_base(k)-z_at_v)/(z_base(k)-z_base(k-1)))
3861 wk = 1.-wkp1-wkm1
3862 rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*( &
3863 wkm1*v_base(k-1) &
3864 +wk *v_base(k ) &
3865 +wkp1*v_base(k+1) )
3866 ENDDO
3867 ENDDO
3868 ENDDO
3871 ! pick up top and bottom v
3873 DO j = jts, MIN(jte,jde-1)+1
3874 DO i = i_start-1, i_end
3876 k = kts
3877 z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
3878 +phb(i,k,j-1)+phb(i,k+1,j-1) &
3879 +ph(i,k,j )+ph(i,k+1,j ) &
3880 +ph(i,k,j-1)+ph(i,k+1,j-1))/g
3881 wkp1 = min(1.,max(0.,z_at_v-z_base(k))/(z_base(k+1)-z_base(k)))
3882 wk = 1.-wkp1
3883 rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*( &
3884 +wk *v_base(k ) &
3885 +wkp1*v_base(k+1) )
3887 k = ktf
3888 z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
3889 +phb(i,k,j-1)+phb(i,k+1,j-1) &
3890 +ph(i,k,j )+ph(i,k+1,j ) &
3891 +ph(i,k,j-1)+ph(i,k+1,j-1))/g
3892 wkm1 = min(1.,max(0.,z_base(k)-z_at_v)/(z_base(k)-z_base(k-1)))
3893 wk = 1.-wkm1
3894 rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*( &
3895 wkm1*v_base(k-1) &
3896 +wk *v_base(k ) )
3898 ENDDO
3899 ENDDO
3901 ! compute coriolis forcing for u
3903 ! Map scale factors: see comments above for Coriolis
3905 DO j = jts, MIN(jte,jde-1)
3907 DO k=kts,ktf
3908 DO i = i_start, i_end
3909 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)) &
3910 *0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
3911 - 0.5*(e(i,j)+e(i-1,j))*0.5*(cosa(i,j)+cosa(i-1,j)) &
3912 *0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
3913 ENDDO
3914 ENDDO
3916 ! boundary loops for perturbation coriolis is needed for open bdy (20110225 JD)
3917 IF ( (config_flags%open_xs) .and. (its == ids) ) THEN
3919 DO k=kts,ktf
3921 ru_tend(its,k,j)=ru_tend(its,k,j) + (msfux(its,j)/msfuy(its,j))*0.5*(f(its,j)+f(its,j)) &
3922 *0.25*(rv(its,k,j+1)+rv(its,k,j+1)+rv(its,k,j)+rv(its,k,j)) &
3923 - 0.5*(e(its,j)+e(its,j))*0.5*(cosa(its,j)+cosa(its,j)) &
3924 *0.25*(rw(its,k+1,j)+rw(its,k,j)+rw(its,k+1,j)+rw(its,k,j))
3926 ENDDO
3928 ENDIF
3930 IF ( (config_flags%open_xe) .and. (ite == ide) ) THEN
3932 DO k=kts,ktf
3934 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)) &
3935 *0.25*(rv(ite-1,k,j+1)+rv(ite-1,k,j+1)+rv(ite-1,k,j)+rv(ite-1,k,j)) &
3936 - 0.5*(e(ite-1,j)+e(ite-1,j))*0.5*(cosa(ite-1,j)+cosa(ite-1,j)) &
3937 *0.25*(rw(ite-1,k+1,j)+rw(ite-1,k,j)+rw(ite-1,k+1,j)+rw(ite-1,k,j))
3939 ENDDO
3941 ENDIF
3943 ENDDO
3945 ! coriolis term for v-momentum equation
3946 ! Map scale factors: see comments above for Coriolis
3948 j_start = jts
3949 j_end = jte
3951 IF ( config_flags%open_ys .or. specified .or. &
3952 config_flags%nested .or. config_flags%polar) j_start = MAX(jds+1,jts)
3953 IF ( config_flags%open_ye .or. specified .or. &
3954 config_flags%nested .or. config_flags%polar) j_end = MIN(jde-1,jte)
3956 ! compute perturbation mu*u for use in v momentum equation
3958 DO j = j_start-1,j_end
3959 DO k=kts+1,ktf-1
3960 DO i = its, MIN(ite,ide-1)+1
3961 z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
3962 +phb(i-1,k,j)+phb(i-1,k+1,j) &
3963 +ph(i ,k,j)+ph(i ,k+1,j) &
3964 +ph(i-1,k,j)+ph(i-1,k+1,j))/g
3965 wkp1 = min(1.,max(0.,z_at_u-z_base(k))/(z_base(k+1)-z_base(k)))
3966 wkm1 = min(1.,max(0.,z_base(k)-z_at_u)/(z_base(k)-z_base(k-1)))
3967 wk = 1.-wkp1-wkm1
3968 ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*( &
3969 wkm1*u_base(k-1) &
3970 +wk *u_base(k ) &
3971 +wkp1*u_base(k+1) )
3972 ENDDO
3973 ENDDO
3974 ENDDO
3976 ! pick up top and bottom u
3978 DO j = j_start-1,j_end
3979 DO i = its, MIN(ite,ide-1)+1
3981 k = kts
3982 z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
3983 +phb(i-1,k,j)+phb(i-1,k+1,j) &
3984 +ph(i ,k,j)+ph(i ,k+1,j) &
3985 +ph(i-1,k,j)+ph(i-1,k+1,j))/g
3986 wkp1 = min(1.,max(0.,z_at_u-z_base(k))/(z_base(k+1)-z_base(k)))
3987 wk = 1.-wkp1
3988 ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*( &
3989 +wk *u_base(k ) &
3990 +wkp1*u_base(k+1) )
3993 k = ktf
3994 z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
3995 +phb(i-1,k,j)+phb(i-1,k+1,j) &
3996 +ph(i ,k,j)+ph(i ,k+1,j) &
3997 +ph(i-1,k,j)+ph(i-1,k+1,j))/g
3998 wkm1 = min(1.,max(0.,z_base(k)-z_at_u)/(z_base(k)-z_base(k-1)))
3999 wk = 1.-wkm1
4000 ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*( &
4001 wkm1*u_base(k-1) &
4002 +wk *u_base(k ) )
4004 ENDDO
4005 ENDDO
4007 ! compute coriolis forcing for v momentum equation
4008 ! Map scale factors: see comments above for Coriolis
4010 ! boundary loops for perturbation coriolis is needed for open bdy (20110225 JD)
4011 IF ( (config_flags%open_ys) .and. (jts == jds) ) THEN
4013 DO k=kts,ktf
4014 DO i=its,MIN(ide-1,ite)
4016 rv_tend(i,k,jts)=rv_tend(i,k,jts) - (msfvy(i,jts)/msfvx(i,jts))*0.5*(f(i,jts)+f(i,jts)) &
4017 *0.25*(ru(i,k,jts)+ru(i+1,k,jts)+ru(i,k,jts)+ru(i+1,k,jts)) &
4018 + (msfvy(i,jts)/msfvx(i,jts))*0.5*(e(i,jts)+e(i,jts))*0.5*(sina(i,jts)+sina(i,jts)) &
4019 *0.25*(rw(i,k+1,jts)+rw(i,k,jts)+rw(i,k+1,jts)+rw(i,k,jts))
4021 ENDDO
4022 ENDDO
4024 ENDIF
4026 DO j=j_start, j_end
4027 DO k=kts,ktf
4028 DO i=its,MIN(ide-1,ite)
4030 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)) &
4031 *0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
4032 + (msfvy(i,j)/msfvx(i,j))*0.5*(e(i,j)+e(i,j-1))*0.5*(sina(i,j)+sina(i,j-1)) &
4033 *0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
4035 ENDDO
4036 ENDDO
4037 ENDDO
4040 ! boundary loops for perturbation coriolis is needed for open bdy (20110225 JD)
4041 IF ( (config_flags%open_ye) .and. (jte == jde) ) THEN
4043 DO k=kts,ktf
4044 DO i=its,MIN(ide-1,ite)
4046 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)) &
4047 *0.25*(ru(i,k,jte-1)+ru(i+1,k,jte-1)+ru(i,k,jte-1)+ru(i+1,k,jte-1)) &
4048 + (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)) &
4049 *0.25*(rw(i,k+1,jte-1)+rw(i,k,jte-1)+rw(i,k+1,jte-1)+rw(i,k,jte-1))
4051 ENDDO
4052 ENDDO
4054 ENDIF
4056 ! coriolis term for w-mometum
4057 ! Map scale factors: see comments above for Coriolis
4059 DO j=jts,MIN(jte, jde-1)
4060 DO k=kts+1,ktf
4061 DO i=its,MIN(ite, ide-1)
4063 rw_tend(i,k,j)=rw_tend(i,k,j) + e(i,j)* &
4064 (cosa(i,j)*0.5*(fzm(k)*(ru(i,k,j)+ru(i+1,k,j)) &
4065 +fzp(k)*(ru(i,k-1,j)+ru(i+1,k-1,j))) &
4066 -(msftx(i,j)/msfty(i,j))*sina(i,j)*0.5*(fzm(k)*(rv(i,k,j)+rv(i,k,j+1)) &
4067 +fzp(k)*(rv(i,k-1,j)+rv(i,k-1,j+1))))
4069 ENDDO
4070 ENDDO
4071 ENDDO
4073 END SUBROUTINE perturbation_coriolis
4075 !------------------------------------------------------------------------------
4077 SUBROUTINE curvature ( ru, rv, rw, u, v, w, ru_tend, rv_tend, rw_tend, &
4078 config_flags, &
4079 msfux, msfuy, msfvx, msfvy, msftx, msfty, &
4080 xlat, fzm, fzp, rdx, rdy, &
4081 ids, ide, jds, jde, kds, kde, &
4082 ims, ime, jms, jme, kms, kme, &
4083 its, ite, jts, jte, kts, kte )
4086 IMPLICIT NONE
4088 ! Input data
4090 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
4092 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4093 ims, ime, jms, jme, kms, kme, &
4094 its, ite, jts, jte, kts, kte
4096 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4097 INTENT(INOUT) :: ru_tend, &
4098 rv_tend, &
4099 rw_tend
4101 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4102 INTENT(IN ) :: ru, &
4103 rv, &
4104 rw, &
4105 u, &
4106 v, &
4107 w
4109 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
4110 msfuy, &
4111 msfvx, &
4112 msfvy, &
4113 msftx, &
4114 msfty, &
4115 xlat
4117 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
4118 fzp
4120 REAL , INTENT(IN ) :: rdx, &
4121 rdy
4123 ! Local data
4125 ! INTEGER :: i, j, k, itf, jtf, ktf, kp1, im, ip, jm, jp
4126 INTEGER :: i, j, k, itf, jtf, ktf
4127 INTEGER :: i_start, i_end, j_start, j_end
4128 ! INTEGER :: irmin, irmax, jrmin, jrmax
4130 REAL , DIMENSION( its-1:ite , kts:kte, jts-1:jte ) :: vxgm
4132 LOGICAL :: specified
4134 !<DESCRIPTION>
4135 !
4136 ! curvature calculates the large timestep tendency terms in the
4137 ! u, v, and w momentum equations arise from the curvature terms.
4138 !
4139 !</DESCRIPTION>
4141 specified = .false.
4142 if(config_flags%specified .or. config_flags%nested) specified = .true.
4144 itf=MIN(ite,ide-1)
4145 jtf=MIN(jte,jde-1)
4146 ktf=MIN(kte,kde-1)
4148 ! irmin = ims
4149 ! irmax = ime
4150 ! jrmin = jms
4151 ! jrmax = jme
4152 ! IF ( config_flags%open_xs ) irmin = ids
4153 ! IF ( config_flags%open_xe ) irmax = ide-1
4154 ! IF ( config_flags%open_ys ) jrmin = jds
4155 ! IF ( config_flags%open_ye ) jrmax = jde-1
4157 ! Define v cross grad m at scalar points - vxgm(i,j)
4159 i_start = its-1
4160 i_end = ite
4161 j_start = jts-1
4162 j_end = jte
4164 IF ( ( config_flags%open_xs .or. specified .or. &
4165 config_flags%nested) .and. (its == ids) ) i_start = its
4166 IF ( ( config_flags%open_xe .or. specified .or. &
4167 config_flags%nested) .and. (ite == ide) ) i_end = ite-1
4168 IF ( ( config_flags%open_ys .or. specified .or. &
4169 config_flags%nested .or. config_flags%polar) .and. (jts == jds) ) j_start = jts
4170 IF ( ( config_flags%open_ye .or. specified .or. &
4171 config_flags%nested .or. config_flags%polar) .and. (jte == jde) ) j_end = jte-1
4172 IF ( config_flags%periodic_x ) i_start = its-1
4173 IF ( config_flags%periodic_x ) i_end = ite
4175 DO j=j_start, j_end
4176 DO k=kts,ktf
4177 DO i=i_start, i_end
4178 ! Map scale factor notes:
4179 ! msf...y is constant everywhere for cylindrical map projection
4180 ! msf...x varies with y only
4181 ! But we know that this is not = 0 for cylindrical,
4182 ! therefore use msfvX in 1st line
4183 ! which => by symmetry use msfuY in 2nd line - ???
4184 vxgm(i,k,j)=0.5*(u(i,k,j)+u(i+1,k,j))*(msfvx(i,j+1)-msfvx(i,j))*rdy - &
4185 0.5*(v(i,k,j)+v(i,k,j+1))*(msfuy(i+1,j)-msfuy(i,j))*rdx
4186 ENDDO
4187 ENDDO
4188 ENDDO
4190 ! Pick up the boundary rows for open (radiation) lateral b.c.
4191 ! Rather crude at present, we are assuming there is no
4192 ! variation in this term at the boundary.
4194 IF ( ( config_flags%open_xs .or. (specified .AND. .NOT. config_flags%periodic_x) .or. &
4195 config_flags%nested) .and. (its == ids) ) THEN
4197 DO j = jts, jte-1
4198 DO k = kts, ktf
4199 vxgm(its-1,k,j) = vxgm(its,k,j)
4200 ENDDO
4201 ENDDO
4203 ENDIF
4205 IF ( ( config_flags%open_xe .or. (specified .AND. .NOT. config_flags%periodic_x) .or. &
4206 config_flags%nested) .and. (ite == ide) ) THEN
4208 DO j = jts, jte-1
4209 DO k = kts, ktf
4210 vxgm(ite,k,j) = vxgm(ite-1,k,j)
4211 ENDDO
4212 ENDDO
4214 ENDIF
4216 ! Polar boundary condition:
4217 ! The following change is needed in case one tries using the vxgm route with
4218 ! polar B.C.'s in the future, but not needed if 'tan' used
4219 IF ( ( config_flags%open_ys .or. specified .or. &
4220 config_flags%nested .or. config_flags%polar) .and. (jts == jds) ) THEN
4222 DO k = kts, ktf
4223 DO i = its-1, ite
4224 vxgm(i,k,jts-1) = vxgm(i,k,jts)
4225 ENDDO
4226 ENDDO
4228 ENDIF
4230 ! Polar boundary condition:
4231 ! The following change is needed in case one tries using the vxgm route with
4232 ! polar B.C.'s in the future, but not needed if 'tan' used
4233 IF ( ( config_flags%open_ye .or. specified .or. &
4234 config_flags%nested .or. config_flags%polar) .and. (jte == jde) ) THEN
4236 DO k = kts, ktf
4237 DO i = its-1, ite
4238 vxgm(i,k,jte) = vxgm(i,k,jte-1)
4239 ENDDO
4240 ENDDO
4242 ENDIF
4244 ! curvature term for u momentum eqn.
4246 ! Map scale factor notes:
4247 ! ADT eqn 44, RHS terms 4 and 5, in cylindrical: mu u v tan(lat)/(a my)
4248 ! - mu u w /(a my)
4249 ! ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
4250 ! => terms are:
4251 ! (mx/my)*u rv tan(lat) / a - u rw / a = (u/a)*[(mx/my) rv tan(lat) - rw]
4252 ! ru v tan(lat) / a - u rw / a
4253 ! xlat defined with end points half grid space from pole,
4254 ! hence are on u latitude points
4256 i_start = its
4257 IF ( config_flags%open_xs .or. specified .or. &
4258 config_flags%nested) i_start = MAX ( ids+1 , its )
4259 IF ( config_flags%open_xe .or. specified .or. &
4260 config_flags%nested) i_end = MIN ( ide-1 , ite )
4261 IF ( config_flags%periodic_x ) i_start = its
4262 IF ( config_flags%periodic_x ) i_end = ite
4264 ! Polar boundary condition
4265 IF ((config_flags%map_proj == 6) .OR. (config_flags%polar)) THEN
4267 DO j=jts,MIN(jde-1,jte)
4268 DO k=kts,ktf
4269 DO i=i_start,i_end
4271 ru_tend(i,k,j)=ru_tend(i,k,j) + u(i,k,j)*reradius* ( &
4272 (msfux(i,j)/msfuy(i,j))*0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+ &
4273 rv(i-1,k,j)+rv(i,k,j))*tan(xlat(i,j)*degrad) &
4274 - 0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j)) )
4275 ENDDO
4276 ENDDO
4277 ENDDO
4279 ELSE ! normal code
4282 DO j=jts,MIN(jde-1,jte)
4283 DO k=kts,ktf
4284 DO i=i_start,i_end
4286 ru_tend(i,k,j)=ru_tend(i,k,j) + 0.5*(vxgm(i,k,j)+vxgm(i-1,k,j)) &
4287 *0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
4288 - u(i,k,j)*reradius &
4289 *0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
4291 ENDDO
4292 ENDDO
4293 ENDDO
4295 END IF
4297 ! curvature term for v momentum eqn.
4299 ! Map scale factor notes
4300 ! ADT eqn 45, RHS terms 4 and 5, in cylindrical: - mu u*u tan(lat)/(a mx)
4301 ! - mu v w /(a mx)
4302 ! ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
4303 ! terms are:
4304 ! - (my/mx)*u ru tan(lat) / a - (my/mx)*v rw / a
4305 ! = - [my/(mx*a)]*[u ru tan(lat) + v rw]
4306 ! - (1/a)*[(my/mx)*u ru tan(lat) + w rv]
4307 ! xlat defined with end points half grid space from pole, hence are on
4308 ! u latitude points => av here
4309 !
4310 ! in original wrf, there was a sign error for the rw contribution
4312 j_start = jts
4313 IF ( config_flags%open_ys .or. specified .or. &
4314 config_flags%nested .or. config_flags%polar) j_start = MAX ( jds+1 , jts )
4315 IF ( config_flags%open_ye .or. specified .or. &
4316 config_flags%nested .or. config_flags%polar) j_end = MIN ( jde-1 , jte )
4318 IF ((config_flags%map_proj == 6) .OR. (config_flags%polar)) THEN
4320 DO j=j_start,j_end
4321 DO k=kts,ktf
4322 DO i=its,MIN(ite,ide-1)
4323 rv_tend(i,k,j)=rv_tend(i,k,j) - (msfvy(i,j)/msfvx(i,j))*reradius* ( &
4324 0.25*(u(i,k,j)+u(i+1,k,j)+u(i,k,j-1)+u(i+1,k,j-1))* &
4325 tan((xlat(i,j)+xlat(i,j-1))*0.5*degrad)* &
4326 0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
4327 + v(i,k,j)*0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+ &
4328 rw(i,k+1,j)+rw(i,k,j)) )
4329 ENDDO
4330 ENDDO
4331 ENDDO
4333 ELSE ! normal code
4335 DO j=j_start,j_end
4336 DO k=kts,ktf
4337 DO i=its,MIN(ite,ide-1)
4339 rv_tend(i,k,j)=rv_tend(i,k,j) - 0.5*(vxgm(i,k,j)+vxgm(i,k,j-1)) &
4340 *0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
4341 - (msfvy(i,j)/msfvx(i,j))*v(i,k,j)*reradius &
4342 *0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
4344 ENDDO
4345 ENDDO
4346 ENDDO
4348 END IF
4350 ! curvature term for vertical momentum eqn.
4352 ! Notes on map scale factors:
4353 ! ADT eqn 46, RHS term 4: [mu/(a my)]*[u*u + v*v]
4354 ! ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
4355 ! terms are: u ru / a + (mx/my)v rv / a
4357 DO j=jts,MIN(jte,jde-1)
4358 DO k=MAX(2,kts),ktf
4359 DO i=its,MIN(ite,ide-1)
4361 rw_tend(i,k,j)=rw_tend(i,k,j) + reradius* &
4362 (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))) &
4363 *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))) &
4364 +(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))) &
4365 *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))))
4367 ENDDO
4368 ENDDO
4369 ENDDO
4371 END SUBROUTINE curvature
4373 !------------------------------------------------------------------------------
4375 SUBROUTINE decouple ( rr, rfield, field, name, config_flags, &
4376 fzm, fzp, &
4377 ids, ide, jds, jde, kds, kde, &
4378 ims, ime, jms, jme, kms, kme, &
4379 its, ite, jts, jte, kts, kte )
4381 IMPLICIT NONE
4383 ! Input data
4385 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
4387 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4388 ims, ime, jms, jme, kms, kme, &
4389 its, ite, jts, jte, kts, kte
4391 CHARACTER(LEN=1) , INTENT(IN ) :: name
4393 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rfield
4395 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rr
4397 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: field
4399 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, fzp
4401 ! Local data
4403 INTEGER :: i, j, k, itf, jtf, ktf
4405 !<DESCRIPTION>
4406 !
4407 ! decouple decouples a variable from the column dry-air mass.
4408 !
4409 !</DESCRIPTION>
4411 ktf=MIN(kte,kde-1)
4413 IF (name .EQ. 'u')THEN
4414 itf=ite
4415 jtf=MIN(jte,jde-1)
4417 DO j=jts,jtf
4418 DO k=kts,ktf
4419 DO i=its,itf
4420 field(i,k,j)=rfield(i,k,j)/(0.5*(rr(i,k,j)+rr(i-1,k,j)))
4421 ENDDO
4422 ENDDO
4423 ENDDO
4425 ELSE IF (name .EQ. 'v')THEN
4426 itf=MIN(ite,ide-1)
4427 jtf=jte
4429 DO j=jts,jtf
4430 DO k=kts,ktf
4431 DO i=its,itf
4432 field(i,k,j)=rfield(i,k,j)/(0.5*(rr(i,k,j)+rr(i,k,j-1)))
4433 ENDDO
4434 ENDDO
4435 ENDDO
4437 ELSE IF (name .EQ. 'w')THEN
4438 itf=MIN(ite,ide-1)
4439 jtf=MIN(jte,jde-1)
4440 DO j=jts,jtf
4441 DO k=kts+1,ktf
4442 DO i=its,itf
4443 field(i,k,j)=rfield(i,k,j)/(fzm(k)*rr(i,k,j)+fzp(k)*rr(i,k-1,j))
4444 ENDDO
4445 ENDDO
4446 ENDDO
4448 DO j=jts,jtf
4449 DO i=its,itf
4450 field(i,kte,j) = 0.
4451 ENDDO
4452 ENDDO
4454 ELSE
4455 itf=MIN(ite,ide-1)
4456 jtf=MIN(jte,jde-1)
4457 ! For theta we will decouple tb and tp and add them to give t afterwards
4458 DO j=jts,jtf
4459 DO k=kts,ktf
4460 DO i=its,itf
4461 field(i,k,j)=rfield(i,k,j)/rr(i,k,j)
4462 ENDDO
4463 ENDDO
4464 ENDDO
4466 ENDIF
4468 END SUBROUTINE decouple
4470 !-------------------------------------------------------------------------------
4473 SUBROUTINE zero_tend ( tendency, &
4474 ids, ide, jds, jde, kds, kde, &
4475 ims, ime, jms, jme, kms, kme, &
4476 its, ite, jts, jte, kts, kte )
4479 IMPLICIT NONE
4481 ! Input data
4483 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4484 ims, ime, jms, jme, kms, kme, &
4485 its, ite, jts, jte, kts, kte
4487 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
4489 ! Local data
4491 INTEGER :: i, j, k, itf, jtf, ktf
4493 !<DESCRIPTION>
4494 !
4495 ! zero_tend sets the input tendency array to zero.
4496 !
4497 !</DESCRIPTION>
4499 DO j = jts, jte
4500 DO k = kts, kte
4501 DO i = its, ite
4502 tendency(i,k,j) = 0.
4503 ENDDO
4504 ENDDO
4505 ENDDO
4507 END SUBROUTINE zero_tend
4509 !-------------------------------------------------------------------------------
4510 ! Sets the an array on the polar v point(s) to zero
4511 SUBROUTINE zero_pole ( field, &
4512 ids, ide, jds, jde, kds, kde, &
4513 ims, ime, jms, jme, kms, kme, &
4514 its, ite, jts, jte, kts, kte )
4517 IMPLICIT NONE
4519 ! Input data
4521 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4522 ims, ime, jms, jme, kms, kme, &
4523 its, ite, jts, jte, kts, kte
4525 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: field
4527 ! Local data
4529 INTEGER :: i, k
4531 IF (jts == jds) THEN
4532 DO k = kts, kte
4533 DO i = its-1, ite+1
4534 field(i,k,jts) = 0.
4535 END DO
4536 END DO
4537 END IF
4538 IF (jte == jde) THEN
4539 DO k = kts, kte
4540 DO i = its-1, ite+1
4541 field(i,k,jte) = 0.
4542 END DO
4543 END DO
4544 END IF
4546 END SUBROUTINE zero_pole
4548 !-------------------------------------------------------------------------------
4549 ! Sets the an array on the polar v point(s)
4550 SUBROUTINE pole_point_bc ( field, &
4551 ids, ide, jds, jde, kds, kde, &
4552 ims, ime, jms, jme, kms, kme, &
4553 its, ite, jts, jte, kts, kte )
4556 IMPLICIT NONE
4558 ! Input data
4560 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4561 ims, ime, jms, jme, kms, kme, &
4562 its, ite, jts, jte, kts, kte
4564 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: field
4566 ! Local data
4568 INTEGER :: i, k
4570 IF (jts == jds) THEN
4571 DO k = kts, kte
4572 DO i = its, ite
4573 ! field(i,k,jts) = 2*field(i,k,jts+1) - field(i,k,jts+2)
4574 field(i,k,jts) = field(i,k,jts+1)
4575 END DO
4576 END DO
4577 END IF
4578 IF (jte == jde) THEN
4579 DO k = kts, kte
4580 DO i = its, ite
4581 ! field(i,k,jte) = 2*field(i,k,jte-1) - field(i,k,jte-2)
4582 field(i,k,jte) = field(i,k,jte-1)
4583 END DO
4584 END DO
4585 END IF
4587 END SUBROUTINE pole_point_bc
4589 !======================================================================
4590 ! physics prep routines
4591 !======================================================================
4593 SUBROUTINE phy_prep ( config_flags, & ! input
4594 mu, muu, muv, u, v, p, pb, alt, ph, & ! input
4595 phb, t, tsk, moist, n_moist, & ! input
4596 rho, th_phy, p_phy , pi_phy , & ! output
4597 u_phy, v_phy, p8w, t_phy, t8w, & ! output
4598 z, z_at_w, dz8w, & ! output
4599 p_hyd, p_hyd_w, dnw, & ! output
4600 fzm, fzp, znw, p_top, & ! params
4601 RTHRATEN, &
4602 RTHBLTEN, RUBLTEN, RVBLTEN, &
4603 RQVBLTEN, RQCBLTEN, RQIBLTEN, &
4604 RUCUTEN, RVCUTEN, RTHCUTEN, &
4605 RQVCUTEN, RQCCUTEN, RQRCUTEN, &
4606 RQICUTEN, RQSCUTEN, &
4607 RUSHTEN, RVSHTEN, RTHSHTEN, &
4608 RQVSHTEN, RQCSHTEN, RQRSHTEN, &
4609 RQISHTEN, RQSSHTEN, RQGSHTEN, &
4610 RTHFTEN, RQVFTEN, &
4611 RUNDGDTEN, RVNDGDTEN, RTHNDGDTEN, &
4612 RPHNDGDTEN,RQVNDGDTEN, RMUNDGDTEN, &
4613 ids, ide, jds, jde, kds, kde, &
4614 ims, ime, jms, jme, kms, kme, &
4615 its, ite, jts, jte, kts, kte )
4616 !----------------------------------------------------------------------
4617 IMPLICIT NONE
4618 !----------------------------------------------------------------------
4620 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
4622 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4623 ims, ime, jms, jme, kms, kme, &
4624 its, ite, jts, jte, kts, kte
4625 INTEGER , INTENT(IN ) :: n_moist
4627 REAL, DIMENSION( ims:ime, kms:kme , jms:jme , n_moist ), INTENT(IN) :: moist
4630 REAL , DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: TSK, mu, muu, muv
4632 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4633 INTENT( OUT) :: u_phy, &
4634 v_phy, &
4635 pi_phy, &
4636 p_phy, &
4637 p8w, &
4638 t_phy, &
4639 th_phy, &
4640 t8w, &
4641 rho, &
4642 z, &
4643 dz8w, &
4644 z_at_w
4646 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4647 INTENT( OUT) :: p_hyd, &
4648 p_hyd_w
4650 REAL , INTENT(IN ) :: p_top
4652 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4653 INTENT(IN ) :: pb, &
4654 p, &
4655 u, &
4656 v, &
4657 alt, &
4658 ph, &
4659 phb, &
4660 t
4663 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
4664 fzp
4666 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: znw, &
4667 dnw
4669 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4670 INTENT(INOUT) :: RTHRATEN
4672 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4673 INTENT(INOUT) :: RUCUTEN, &
4674 RVCUTEN, &
4675 RTHCUTEN, &
4676 RQVCUTEN, &
4677 RQCCUTEN, &
4678 RQRCUTEN, &
4679 RQICUTEN, &
4680 RQSCUTEN, &
4681 RUSHTEN, &
4682 RVSHTEN, &
4683 RTHSHTEN, &
4684 RQVSHTEN, &
4685 RQCSHTEN, &
4686 RQRSHTEN, &
4687 RQISHTEN, &
4688 RQSSHTEN, &
4689 RQGSHTEN
4691 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
4692 INTENT(INOUT) :: RUBLTEN, &
4693 RVBLTEN, &
4694 RTHBLTEN, &
4695 RQVBLTEN, &
4696 RQCBLTEN, &
4697 RQIBLTEN
4699 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
4700 INTENT(INOUT) :: RTHFTEN, &
4701 RQVFTEN
4703 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
4704 INTENT(INOUT) :: RUNDGDTEN, &
4705 RVNDGDTEN, &
4706 RTHNDGDTEN, &
4707 RPHNDGDTEN, &
4708 RQVNDGDTEN
4710 REAL, DIMENSION( ims:ime, jms:jme ) , &
4711 INTENT(INOUT) :: RMUNDGDTEN
4713 INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end, i_startu, j_startv
4714 INTEGER :: i, j, k
4715 REAL :: w1, w2, z0, z1, z2
4716 REAL :: qtot
4717 INTEGER :: n
4719 !-----------------------------------------------------------------------
4721 !<DESCRIPTION>
4722 !
4723 ! phys_prep calculates a number of diagnostic quantities needed by
4724 ! the physics routines. It also decouples the physics tendencies from
4725 ! the column dry-air mass (the physics routines expect to see/update the
4726 ! uncoupled tendencies).
4727 !
4728 !</DESCRIPTION>
4730 ! set up loop bounds for this grid's boundary conditions
4732 i_start = its
4733 i_end = min( ite,ide-1 )
4734 j_start = jts
4735 j_end = min( jte,jde-1 )
4737 k_start = kts
4738 k_end = min( kte, kde-1 )
4740 ! compute thermodynamics and velocities at pressure points (or half levels)
4742 do j = j_start,j_end
4743 do k = k_start, k_end
4744 do i = i_start, i_end
4746 th_phy(i,k,j) = t(i,k,j) + t0
4747 p_phy(i,k,j) = p(i,k,j) + pb(i,k,j)
4748 pi_phy(i,k,j) = (p_phy(i,k,j)/p1000mb)**rcp
4749 t_phy(i,k,j) = th_phy(i,k,j)*pi_phy(i,k,j)
4750 rho(i,k,j) = 1./alt(i,k,j)*(1.+moist(i,k,j,P_QV))
4751 u_phy(i,k,j) = 0.5*(u(i,k,j)+u(i+1,k,j))
4752 v_phy(i,k,j) = 0.5*(v(i,k,j)+v(i,k,j+1))
4754 enddo
4755 enddo
4756 enddo
4758 ! compute z at w points
4760 do j = j_start,j_end
4761 do k = k_start, kte
4762 do i = i_start, i_end
4763 z_at_w(i,k,j) = (phb(i,k,j)+ph(i,k,j))/g
4764 enddo
4765 enddo
4766 enddo
4768 do j = j_start,j_end
4769 do k = k_start, kte-1
4770 do i = i_start, i_end
4771 dz8w(i,k,j) = z_at_w(i,k+1,j)-z_at_w(i,k,j)
4772 enddo
4773 enddo
4774 enddo
4776 do j = j_start,j_end
4777 do i = i_start, i_end
4778 dz8w(i,kte,j) = 0.
4779 enddo
4780 enddo
4782 ! compute z at p points or half levels (average of z at full levels)
4784 do j = j_start,j_end
4785 do k = k_start, k_end
4786 do i = i_start, i_end
4787 z(i,k,j) = 0.5*(z_at_w(i,k,j) +z_at_w(i,k+1,j) )
4788 enddo
4789 enddo
4790 enddo
4792 ! interp t and p to full levels
4794 do j = j_start,j_end
4795 do k = 2, k_end
4796 do i = i_start, i_end
4797 p8w(i,k,j) = fzm(k)*p_phy(i,k,j)+fzp(k)*p_phy(i,k-1,j)
4798 t8w(i,k,j) = fzm(k)*t_phy(i,k,j)+fzp(k)*t_phy(i,k-1,j)
4799 enddo
4800 enddo
4801 enddo
4803 ! extrapolate p and t to surface and top.
4804 ! we'll use an extrapolation in z for now
4806 do j = j_start,j_end
4807 do i = i_start, i_end
4809 ! bottom
4811 z0 = z_at_w(i,1,j)
4812 z1 = z(i,1,j)
4813 z2 = z(i,2,j)
4814 w1 = (z0 - z2)/(z1 - z2)
4815 w2 = 1. - w1
4816 p8w(i,1,j) = w1*p_phy(i,1,j)+w2*p_phy(i,2,j)
4817 t8w(i,1,j) = w1*t_phy(i,1,j)+w2*t_phy(i,2,j)
4819 ! top
4821 z0 = z_at_w(i,kte,j)
4822 z1 = z(i,k_end,j)
4823 z2 = z(i,k_end-1,j)
4824 w1 = (z0 - z2)/(z1 - z2)
4825 w2 = 1. - w1
4827 ! p8w(i,kde,j) = w1*p_phy(i,kde-1,j)+w2*p_phy(i,kde-2,j)
4828 !!! bug fix extrapolate ln(p) so p is positive definite
4829 p8w(i,kde,j) = exp(w1*log(p_phy(i,kde-1,j))+w2*log(p_phy(i,kde-2,j)))
4830 t8w(i,kde,j) = w1*t_phy(i,kde-1,j)+w2*t_phy(i,kde-2,j)
4832 enddo
4833 enddo
4835 ! calculate hydrostatic pressure at both full and half levels
4836 ! first, full level p: assuming dry over model top
4838 do j = j_start,j_end
4839 do i = i_start, i_end
4840 p_hyd_w(i,kte,j) = p_top
4841 enddo
4842 enddo
4844 do j = j_start,j_end
4845 do k = kte-1, k_start, -1
4846 do i = i_start, i_end
4847 qtot = 0.
4848 do n = PARAM_FIRST_SCALAR,n_moist
4849 qtot = qtot + moist(i,k,j,n)
4850 enddo
4851 p_hyd_w(i,k,j) = p_hyd_w(i,k+1,j) - (1.+qtot)*mu(i,j)*dnw(k)
4852 ! p_hyd_w(i,k,j) = p_hyd_w(i,k+1,j)+1./alt(i,k,j)*(1.+moist(i,k,j,P_QV))*g*dz8w(i,k,j)
4853 enddo
4854 enddo
4855 enddo
4857 ! now calculate hydrostatic pressure at half levels
4859 do j = j_start,j_end
4860 do k = k_start, k_end
4861 do i = i_start, i_end
4862 p_hyd(i,k,j) = 0.5*(p_hyd_w(i,k,j)+p_hyd_w(i,k+1,j))
4863 enddo
4864 enddo
4865 enddo
4867 ! decouple all physics tendencies
4869 IF (config_flags%ra_lw_physics .gt. 0 .or. config_flags%ra_sw_physics .gt. 0) THEN
4871 DO J=j_start,j_end
4872 DO K=k_start,k_end
4873 DO I=i_start,i_end
4874 RTHRATEN(I,K,J)=RTHRATEN(I,K,J)/mu(I,J)
4875 ENDDO
4876 ENDDO
4877 ENDDO
4879 ENDIF
4881 IF (config_flags%cu_physics .gt. 0) THEN
4883 DO J=j_start,j_end
4884 DO I=i_start,i_end
4885 DO K=k_start,k_end
4886 RUCUTEN(I,K,J) =RUCUTEN(I,K,J)/mu(I,J)
4887 RVCUTEN(I,K,J) =RVCUTEN(I,K,J)/mu(I,J)
4888 RTHCUTEN(I,K,J)=RTHCUTEN(I,K,J)/mu(I,J)
4889 ENDDO
4890 ENDDO
4891 ENDDO
4893 IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
4894 DO J=j_start,j_end
4895 DO I=i_start,i_end
4896 DO K=k_start,k_end
4897 RQVCUTEN(I,K,J)=RQVCUTEN(I,K,J)/mu(I,J)
4898 ENDDO
4899 ENDDO
4900 ENDDO
4901 ENDIF
4903 IF (P_QC .ge. PARAM_FIRST_SCALAR)THEN
4904 DO J=j_start,j_end
4905 DO I=i_start,i_end
4906 DO K=k_start,k_end
4907 RQCCUTEN(I,K,J)=RQCCUTEN(I,K,J)/mu(I,J)
4908 ENDDO
4909 ENDDO
4910 ENDDO
4911 ENDIF
4913 IF (P_QR .ge. PARAM_FIRST_SCALAR)THEN
4914 DO J=j_start,j_end
4915 DO I=i_start,i_end
4916 DO K=k_start,k_end
4917 RQRCUTEN(I,K,J)=RQRCUTEN(I,K,J)/mu(I,J)
4918 ENDDO
4919 ENDDO
4920 ENDDO
4921 ENDIF
4923 IF (P_QI .ge. PARAM_FIRST_SCALAR)THEN
4924 DO J=j_start,j_end
4925 DO I=i_start,i_end
4926 DO K=k_start,k_end
4927 RQICUTEN(I,K,J)=RQICUTEN(I,K,J)/mu(I,J)
4928 ENDDO
4929 ENDDO
4930 ENDDO
4931 ENDIF
4933 IF(P_QS .ge. PARAM_FIRST_SCALAR)THEN
4934 DO J=j_start,j_end
4935 DO I=i_start,i_end
4936 DO K=k_start,k_end
4937 RQSCUTEN(I,K,J)=RQSCUTEN(I,K,J)/mu(I,J)
4938 ENDDO
4939 ENDDO
4940 ENDDO
4941 ENDIF
4943 ENDIF
4945 IF (config_flags%shcu_physics .gt. 0) THEN
4947 DO J=j_start,j_end
4948 DO I=i_start,i_end
4949 DO K=k_start,k_end
4950 RUSHTEN(I,K,J) =RUSHTEN(I,K,J)/mu(I,J)
4951 RVSHTEN(I,K,J) =RVSHTEN(I,K,J)/mu(I,J)
4952 RTHSHTEN(I,K,J)=RTHSHTEN(I,K,J)/mu(I,J)
4953 ENDDO
4954 ENDDO
4955 ENDDO
4957 IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
4958 DO J=j_start,j_end
4959 DO I=i_start,i_end
4960 DO K=k_start,k_end
4961 RQVSHTEN(I,K,J)=RQVSHTEN(I,K,J)/mu(I,J)
4962 ENDDO
4963 ENDDO
4964 ENDDO
4965 ENDIF
4967 IF (P_QC .ge. PARAM_FIRST_SCALAR)THEN
4968 DO J=j_start,j_end
4969 DO I=i_start,i_end
4970 DO K=k_start,k_end
4971 RQCSHTEN(I,K,J)=RQCSHTEN(I,K,J)/mu(I,J)
4972 ENDDO
4973 ENDDO
4974 ENDDO
4975 ENDIF
4977 IF (P_QR .ge. PARAM_FIRST_SCALAR)THEN
4978 DO J=j_start,j_end
4979 DO I=i_start,i_end
4980 DO K=k_start,k_end
4981 RQRSHTEN(I,K,J)=RQRSHTEN(I,K,J)/mu(I,J)
4982 ENDDO
4983 ENDDO
4984 ENDDO
4985 ENDIF
4987 IF (P_QI .ge. PARAM_FIRST_SCALAR)THEN
4988 DO J=j_start,j_end
4989 DO I=i_start,i_end
4990 DO K=k_start,k_end
4991 RQISHTEN(I,K,J)=RQISHTEN(I,K,J)/mu(I,J)
4992 ENDDO
4993 ENDDO
4994 ENDDO
4995 ENDIF
4997 IF(P_QS .ge. PARAM_FIRST_SCALAR)THEN
4998 DO J=j_start,j_end
4999 DO I=i_start,i_end
5000 DO K=k_start,k_end
5001 RQSSHTEN(I,K,J)=RQSSHTEN(I,K,J)/mu(I,J)
5002 ENDDO
5003 ENDDO
5004 ENDDO
5005 ENDIF
5007 IF(P_QG .ge. PARAM_FIRST_SCALAR)THEN
5008 DO J=j_start,j_end
5009 DO I=i_start,i_end
5010 DO K=k_start,k_end
5011 RQGSHTEN(I,K,J)=RQGSHTEN(I,K,J)/mu(I,J)
5012 ENDDO
5013 ENDDO
5014 ENDDO
5015 ENDIF
5017 ENDIF
5019 IF (config_flags%bl_pbl_physics .gt. 0) THEN
5021 DO J=j_start,j_end
5022 DO K=k_start,k_end
5023 DO I=i_start,i_end
5024 RUBLTEN(I,K,J) =RUBLTEN(I,K,J)/mu(I,J)
5025 RVBLTEN(I,K,J) =RVBLTEN(I,K,J)/mu(I,J)
5026 RTHBLTEN(I,K,J)=RTHBLTEN(I,K,J)/mu(I,J)
5027 ENDDO
5028 ENDDO
5029 ENDDO
5031 IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
5032 DO J=j_start,j_end
5033 DO K=k_start,k_end
5034 DO I=i_start,i_end
5035 RQVBLTEN(I,K,J)=RQVBLTEN(I,K,J)/mu(I,J)
5036 ENDDO
5037 ENDDO
5038 ENDDO
5039 ENDIF
5041 IF (P_QC .ge. PARAM_FIRST_SCALAR) THEN
5042 DO J=j_start,j_end
5043 DO K=k_start,k_end
5044 DO I=i_start,i_end
5045 RQCBLTEN(I,K,J)=RQCBLTEN(I,K,J)/mu(I,J)
5046 ENDDO
5047 ENDDO
5048 ENDDO
5049 ENDIF
5051 IF (P_QI .ge. PARAM_FIRST_SCALAR) THEN
5052 DO J=j_start,j_end
5053 DO K=k_start,k_end
5054 DO I=i_start,i_end
5055 RQIBLTEN(I,K,J)=RQIBLTEN(I,K,J)/mu(I,J)
5056 ENDDO
5057 ENDDO
5058 ENDDO
5059 ENDIF
5061 ENDIF
5063 ! decouple advective forcing required by Grell-Devenyi scheme
5065 if(( config_flags%cu_physics == GDSCHEME ) .OR. &
5066 ( config_flags%cu_physics == G3SCHEME ) .OR. &
5067 ( config_flags%cu_physics == KFETASCHEME ) .OR. &
5068 ( config_flags%cu_physics == TIEDTKESCHEME ) ) then ! Tiedtke ZCX&YQW
5070 DO J=j_start,j_end
5071 DO I=i_start,i_end
5072 DO K=k_start,k_end
5073 RTHFTEN(I,K,J)=RTHFTEN(I,K,J)/mu(I,J)
5074 ENDDO
5075 ENDDO
5076 ENDDO
5078 IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
5079 DO J=j_start,j_end
5080 DO I=i_start,i_end
5081 DO K=k_start,k_end
5082 RQVFTEN(I,K,J)=RQVFTEN(I,K,J)/mu(I,J)
5083 ENDDO
5084 ENDDO
5085 ENDDO
5086 ENDIF
5088 END IF
5090 ! fdda
5091 ! note fdda u and v tendencies are staggered, also only interior points have muu/muv,
5092 ! so only decouple those
5094 IF (config_flags%grid_fdda .gt. 0) THEN
5096 i_startu=MAX(its,ids+1)
5097 j_startv=MAX(jts,jds+1)
5099 DO J=j_start,j_end
5100 DO K=k_start,k_end
5101 DO I=i_startu,i_end
5102 RUNDGDTEN(I,K,J) =RUNDGDTEN(I,K,J)/muu(I,J)
5103 ENDDO
5104 ENDDO
5105 ENDDO
5106 DO J=j_startv,j_end
5107 DO K=k_start,k_end
5108 DO I=i_start,i_end
5109 RVNDGDTEN(I,K,J) =RVNDGDTEN(I,K,J)/muv(I,J)
5110 ENDDO
5111 ENDDO
5112 ENDDO
5113 DO J=j_start,j_end
5114 DO K=k_start,k_end
5115 DO I=i_start,i_end
5116 RTHNDGDTEN(I,K,J)=RTHNDGDTEN(I,K,J)/mu(I,J)
5117 ! RMUNDGDTEN(I,J) - no coupling
5118 ENDDO
5119 ENDDO
5120 ENDDO
5122 IF (config_flags%grid_fdda .EQ. 2) THEN
5123 DO J=j_start,j_end
5124 DO K=k_start,kte
5125 DO I=i_start,i_end
5126 RPHNDGDTEN(I,K,J)=RPHNDGDTEN(I,K,J)/mu(I,J)
5127 ENDDO
5128 ENDDO
5129 ENDDO
5131 ELSE IF (config_flags%grid_fdda .EQ. 1) THEN
5132 IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
5133 DO J=j_start,j_end
5134 DO K=k_start,k_end
5135 DO I=i_start,i_end
5136 RQVNDGDTEN(I,K,J)=RQVNDGDTEN(I,K,J)/mu(I,J)
5137 ENDDO
5138 ENDDO
5139 ENDDO
5140 ENDIF
5141 ENDIF
5143 ENDIF
5145 END SUBROUTINE phy_prep
5147 !------------------------------------------------------------
5149 SUBROUTINE moist_physics_prep_em( t_new, t_old, t0, rho, al, alb, &
5150 p, p8w, p0, pb, ph, phb, &
5151 th_phy, pii, pf, &
5152 z, z_at_w, dz8w, &
5153 dt,h_diabatic, &
5154 config_flags,fzm, fzp, &
5155 ids,ide, jds,jde, kds,kde, &
5156 ims,ime, jms,jme, kms,kme, &
5157 its,ite, jts,jte, kts,kte )
5159 IMPLICIT NONE
5161 ! Here we construct full fields
5162 ! needed by the microphysics
5164 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
5166 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
5167 INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
5168 INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
5170 REAL, INTENT(IN ) :: dt
5172 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5173 INTENT(IN ) :: al, &
5174 alb, &
5175 p, &
5176 pb, &
5177 ph, &
5178 phb
5181 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
5182 fzp
5184 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5185 INTENT( OUT) :: rho, &
5186 th_phy, &
5187 pii, &
5188 pf, &
5189 z, &
5190 z_at_w, &
5191 dz8w, &
5192 p8w
5194 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5195 INTENT(INOUT) :: h_diabatic
5197 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5198 INTENT(INOUT) :: t_new, &
5199 t_old
5201 REAL, INTENT(IN ) :: t0, p0
5202 REAL :: z0,z1,z2,w1,w2
5204 INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
5205 INTEGER :: i, j, k
5207 !--------------------------------------------------------------------
5209 !<DESCRIPTION>
5210 !
5211 ! moist_phys_prep_em calculates a number of diagnostic quantities needed by
5212 ! the microphysics routines.
5213 !
5214 !</DESCRIPTION>
5216 ! set up loop bounds for this grid's boundary conditions
5218 i_start = its
5219 i_end = min( ite,ide-1 )
5220 j_start = jts
5221 j_end = min( jte,jde-1 )
5223 k_start = kts
5224 k_end = min( kte, kde-1 )
5226 DO j = j_start, j_end
5227 DO k = k_start, kte
5228 DO i = i_start, i_end
5229 z_at_w(i,k,j) = (ph(i,k,j)+phb(i,k,j))/g
5230 ENDDO
5231 ENDDO
5232 ENDDO
5234 do j = j_start,j_end
5235 do k = k_start, kte-1
5236 do i = i_start, i_end
5237 dz8w(i,k,j) = z_at_w(i,k+1,j)-z_at_w(i,k,j)
5238 enddo
5239 enddo
5240 enddo
5242 do j = j_start,j_end
5243 do i = i_start, i_end
5244 dz8w(i,kte,j) = 0.
5245 enddo
5246 enddo
5249 ! compute full pii, rho, and z at the new time-level
5250 ! (needed for physics).
5251 ! convert perturbation theta to full theta (th_phy)
5252 ! use h_diabatic to temporarily save pre-microphysics full theta
5254 DO j = j_start, j_end
5255 DO k = k_start, k_end
5256 DO i = i_start, i_end
5258 #ifdef REVERT
5259 t_new(i,k,j) = t_new(i,k,j)-h_diabatic(i,k,j)*dt
5260 #endif
5261 th_phy(i,k,j) = t_new(i,k,j) + t0
5262 h_diabatic(i,k,j) = th_phy(i,k,j)
5263 rho(i,k,j) = 1./(al(i,k,j)+alb(i,k,j))
5264 pii(i,k,j) = ((p(i,k,j)+pb(i,k,j))/p0)**rcp
5265 z(i,k,j) = 0.5*(z_at_w(i,k,j) +z_at_w(i,k+1,j) )
5266 pf(i,k,j) = p(i,k,j)+pb(i,k,j)
5268 ENDDO
5269 ENDDO
5270 ENDDO
5272 ! interp t and p at w points
5274 do j = j_start,j_end
5275 do k = 2, k_end
5276 do i = i_start, i_end
5277 p8w(i,k,j) = fzm(k)*pf(i,k,j)+fzp(k)*pf(i,k-1,j)
5278 enddo
5279 enddo
5280 enddo
5282 ! extrapolate p and t to surface and top.
5283 ! we'll use an extrapolation in z for now
5285 do j = j_start,j_end
5286 do i = i_start, i_end
5288 ! bottom
5290 z0 = z_at_w(i,1,j)
5291 z1 = z(i,1,j)
5292 z2 = z(i,2,j)
5293 w1 = (z0 - z2)/(z1 - z2)
5294 w2 = 1. - w1
5295 p8w(i,1,j) = w1*pf(i,1,j)+w2*pf(i,2,j)
5297 ! top
5299 z0 = z_at_w(i,kte,j)
5300 z1 = z(i,k_end,j)
5301 z2 = z(i,k_end-1,j)
5302 w1 = (z0 - z2)/(z1 - z2)
5303 w2 = 1. - w1
5304 ! p8w(i,kde,j) = w1*pf(i,kde-1,j)+w2*pf(i,kde-2,j)
5305 p8w(i,kde,j) = exp(w1*log(pf(i,kde-1,j))+w2*log(pf(i,kde-2,j)))
5307 enddo
5308 enddo
5310 END SUBROUTINE moist_physics_prep_em
5312 !------------------------------------------------------------------------------
5314 SUBROUTINE moist_physics_finish_em( t_new, t_old, t0, mut, &
5315 th_phy, h_diabatic, dt, &
5316 config_flags, &
5317 #if ( WRF_DFI_RADAR == 1 )
5318 dfi_tten_rad,dfi_stage, &
5319 #endif
5320 ids,ide, jds,jde, kds,kde, &
5321 ims,ime, jms,jme, kms,kme, &
5322 its,ite, jts,jte, kts,kte )
5324 IMPLICIT NONE
5326 ! Here we construct full fields
5327 ! needed by the microphysics
5329 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
5331 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
5332 INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
5333 INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
5335 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5336 INTENT(INOUT) :: t_new, &
5337 t_old, &
5338 th_phy, &
5339 h_diabatic
5340 #if ( WRF_DFI_RADAR == 1 )
5341 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5342 INTENT(IN), OPTIONAL :: dfi_tten_rad
5343 INTEGER, INTENT(IN ) ,OPTIONAL :: dfi_stage
5344 REAL :: dfi_tten_max, old_max
5345 #endif
5347 REAL mpten, mptenmax, mptenmin
5349 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: mut
5352 REAL, INTENT(IN ) :: t0, dt
5354 INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
5355 INTEGER :: i, j, k, imax, jmax, imin, jmin
5357 !--------------------------------------------------------------------
5359 !<DESCRIPTION>
5360 !
5361 ! moist_phys_finish_em resets theta to its perturbation value and
5362 ! computes and stores the microphysics diabatic heating term.
5363 !
5364 !</DESCRIPTION>
5366 ! set up loop bounds for this grid's boundary conditions
5369 i_start = its
5370 i_end = min( ite,ide-1 )
5371 j_start = jts
5372 j_end = min( jte,jde-1 )
5373 ! i_start=max(its,ids+4)
5374 ! i_end=min(ite,ide-5)
5375 ! j_start=max(jts,jds+4)
5376 ! j_end=min(jte,jde-5)
5378 k_start = kts
5379 k_end = min( kte, kde-1 )
5381 #if ( WRF_DFI_RADAR == 1 )
5382 IF ( PRESENT(dfi_stage) .and. PRESENT(dfi_tten_rad) ) THEN
5383 IF ( dfi_stage ==DFI_FWD ) THEN
5384 !$OMP CRITICAL(big_step_utilities)
5385 WRITE(wrf_err_message,*)'Add radar tendency: i_start,j_start: ', i_start, j_start
5386 CALL wrf_debug ( 100 , TRIM(wrf_err_message) )
5387 !$OMP END CRITICAL(big_step_utilities)
5388 ENDIF
5389 ENDIF
5390 dfi_tten_max=-999
5391 old_max=-999
5392 #endif
5394 ! add microphysics theta diff to perturbation theta, set h_diabatic
5396 IF ( config_flags%no_mp_heating .eq. 0 ) THEN
5397 mptenmax = 0.
5398 mptenmin = 999.
5399 DO j = j_start, j_end
5400 DO k = k_start, k_end
5401 DO i = i_start, i_end
5402 mpten = th_phy(i,k,j)-h_diabatic(i,k,j)
5403 #if ( WRF_DFI_RADAR == 1 )
5404 if(mpten.gt.mptenmax) then
5405 mptenmax=mpten
5406 imax=i
5407 jmax=j
5408 endif
5409 if(mpten.lt.mptenmin) then
5410 mptenmin=mpten
5411 imin=i
5412 jmin=j
5413 endif
5414 mpten=min(config_flags%mp_tend_lim*dt, mpten)
5415 mpten=max(-config_flags%mp_tend_lim*dt, mpten)
5417 if(k < k_end ) then
5418 if(dfi_tten_max < dfi_tten_rad(i,k,j) ) dfi_tten_max = dfi_tten_rad(i,k,j)
5419 if(old_max < (th_phy(i,k,j)-h_diabatic(i,k,j)) ) old_max=th_phy(i,k,j)-h_diabatic(i,k,j)
5420 endif
5422 IF ( PRESENT(dfi_stage) .and. PRESENT(dfi_tten_rad) ) THEN
5423 IF ( dfi_stage == DFI_FWD .and. dfi_tten_rad(i,k,j) >= -0.1 .and. &
5424 dfi_tten_rad(i,k,j) <= 0.1 .and. k < k_end ) THEN
5425 ! add radar temp tendency
5426 ! there is radar coverage
5427 t_new(i,k,j) = t_new(i,k,j) + (dfi_tten_rad(i,k,j))*dt
5428 ELSE
5429 ! no radar coverage
5430 t_new(i,k,j) = t_new(i,k,j) + mpten
5431 ENDIF
5432 ENDIF
5433 #else
5434 t_new(i,k,j) = t_new(i,k,j) + mpten
5435 #endif
5436 h_diabatic(i,k,j) = mpten/dt
5437 ENDDO
5438 ENDDO
5439 ENDDO
5441 ELSE
5443 DO j = j_start, j_end
5444 DO k = k_start, k_end
5445 DO i = i_start, i_end
5446 ! t_new(i,k,j) = t_new(i,k,j)
5447 h_diabatic(i,k,j) = 0.
5448 ENDDO
5449 ENDDO
5450 ENDDO
5451 ENDIF
5453 END SUBROUTINE moist_physics_finish_em
5455 !----------------------------------------------------------------
5458 SUBROUTINE init_module_big_step
5459 END SUBROUTINE init_module_big_step
5461 SUBROUTINE set_tend ( field, field_adv_tend, msf, &
5462 ids, ide, jds, jde, kds, kde, &
5463 ims, ime, jms, jme, kms, kme, &
5464 its, ite, jts, jte, kts, kte )
5466 IMPLICIT NONE
5468 ! Input data
5470 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
5471 ims, ime, jms, jme, kms, kme, &
5472 its, ite, jts, jte, kts, kte
5474 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: field
5476 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN) :: field_adv_tend
5478 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: msf
5480 ! Local data
5482 INTEGER :: i, j, k, itf, jtf, ktf
5484 !<DESCRIPTION>
5485 !
5486 ! set_tend copies the advective tendency array into the tendency array.
5487 !
5488 !</DESCRIPTION>
5490 jtf = MIN(jte,jde-1)
5491 ktf = MIN(kte,kde-1)
5492 itf = MIN(ite,ide-1)
5493 DO j = jts, jtf
5494 DO k = kts, ktf
5495 DO i = its, itf
5496 field(i,k,j) = field_adv_tend(i,k,j)*msf(i,j)
5497 ENDDO
5498 ENDDO
5499 ENDDO
5501 END SUBROUTINE set_tend
5503 !------------------------------------------------------------------------------
5505 SUBROUTINE rk_rayleigh_damp( ru_tendf, rv_tendf, &
5506 rw_tendf, t_tendf, &
5507 u, v, w, t, t_init, &
5508 mut, muu, muv, ph, phb, &
5509 u_base, v_base, t_base, z_base, &
5510 dampcoef, zdamp, &
5511 ids, ide, jds, jde, kds, kde, &
5512 ims, ime, jms, jme, kms, kme, &
5513 its, ite, jts, jte, kts, kte )
5515 ! History: Apr 2005 Modifications by George Bryan, NCAR:
5516 ! - Generalized the code in a way that allows for
5517 ! simulations with steep terrain.
5518 !
5519 ! Jul 2004 Modifications by George Bryan, NCAR:
5520 ! - Modified the code to use u_base, v_base, and t_base
5521 ! arrays for the background state. Removed the hard-wired
5522 ! base-state values.
5523 ! - Modified the code to use dampcoef, zdamp, and damp_opt,
5524 ! i.e., the upper-level damper variables in namelist.input.
5525 ! Removed the hard-wired variables in the older version.
5526 ! This damper is used when damp_opt = 2.
5527 ! - Modified the code to account for the movement of the
5528 ! model surfaces with time. The code now obtains a base-
5529 ! state value by interpolation using the "_base" arrays.
5531 ! Nov 2003 Bug fix by Jason Knievel, NCAR
5533 ! Aug 2003 Meridional dimension, some comments, and
5534 ! changes in layout of the code added by
5535 ! Jason Knievel, NCAR
5537 ! Jul 2003 Original code by Bill Skamarock, NCAR
5539 ! Purpose: This routine applies Rayleigh damping to a layer at top
5540 ! of the model domain.
5542 !-----------------------------------------------------------------------
5543 ! Begin declarations.
5545 IMPLICIT NONE
5547 INTEGER, INTENT( IN ) &
5548 :: ids, ide, jds, jde, kds, kde, &
5549 ims, ime, jms, jme, kms, kme, &
5550 its, ite, jts, jte, kts, kte
5552 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
5553 :: ru_tendf, rv_tendf, rw_tendf, t_tendf
5555 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
5556 :: u, v, w, t, t_init, ph, phb
5558 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
5559 :: mut, muu, muv
5561 REAL, DIMENSION( kms:kme ) , INTENT(IN ) &
5562 :: u_base, v_base, t_base, z_base
5564 REAL, INTENT(IN ) &
5565 :: dampcoef, zdamp
5567 ! Local variables.
5569 INTEGER &
5570 :: i_start, i_end, j_start, j_end, k_start, k_end, i, j, k, ktf, k1, k2
5572 REAL &
5573 :: pii, dcoef, z, ztop
5575 REAL :: wkp1, wk, wkm1
5577 REAL, DIMENSION( kms:kme ) :: z00, u00, v00, t00
5579 ! End declarations.
5580 !-----------------------------------------------------------------------
5582 pii = 2.0 * asin(1.0)
5584 ktf = MIN( kte, kde-1 )
5586 !-----------------------------------------------------------------------
5587 ! Adjust u to base state.
5589 DO j = jts, MIN( jte, jde-1 )
5590 DO i = its, MIN( ite, ide )
5592 ! Get height at top of model
5593 ztop = 0.5*( phb(i ,kde,j)+phb(i-1,kde,j) &
5594 +ph(i ,kde,j)+ph(i-1,kde,j) )/g
5596 ! Find bottom of damping layer
5597 k1 = ktf
5598 z = ztop
5599 DO WHILE( z >= (ztop-zdamp) )
5600 z = 0.25*( phb(i ,k1,j)+phb(i ,k1+1,j) &
5601 +phb(i-1,k1,j)+phb(i-1,k1+1,j) &
5602 +ph(i ,k1,j)+ph(i ,k1+1,j) &
5603 +ph(i-1,k1,j)+ph(i-1,k1+1,j))/g
5604 z00(k1) = z
5605 k1 = k1 - 1
5606 ENDDO
5607 k1 = k1 + 2
5609 ! Get reference state at model levels
5610 DO k = k1, ktf
5611 k2 = ktf
5612 DO WHILE( z_base(k2) .gt. z00(k) )
5613 k2 = k2 - 1
5614 ENDDO
5615 if(k2+1.gt.ktf)then
5616 u00(k) = u_base(k2) + ( u_base(k2) - u_base(k2-1) ) &
5617 * ( z00(k) - z_base(k2) ) &
5618 / ( z_base(k2) - z_base(k2-1) )
5619 else
5620 u00(k) = u_base(k2) + ( u_base(k2+1) - u_base(k2) ) &
5621 * ( z00(k) - z_base(k2) ) &
5622 / ( z_base(k2+1) - z_base(k2) )
5623 endif
5624 ENDDO
5626 ! Apply the Rayleigh damper
5627 DO k = k1, ktf
5628 dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
5629 dcoef = (SIN( 0.5 * pii * dcoef ) )**2
5630 ru_tendf(i,k,j) = ru_tendf(i,k,j) - &
5631 muu(i,j) * ( dcoef * dampcoef ) * &
5632 ( u(i,k,j) - u00(k) )
5633 END DO
5635 END DO
5636 END DO
5638 ! End adjustment of u.
5639 !-----------------------------------------------------------------------
5641 !-----------------------------------------------------------------------
5642 ! Adjust v to base state.
5644 DO j = jts, MIN( jte, jde )
5645 DO i = its, MIN( ite, ide-1 )
5647 ! Get height at top of model
5648 ztop = 0.5*( phb(i,kde,j )+phb(i,kde,j-1) &
5649 +ph(i,kde,j )+ph(i,kde,j-1) )/g
5651 ! Find bottom of damping layer
5652 k1 = ktf
5653 z = ztop
5654 DO WHILE( z >= (ztop-zdamp) )
5655 z = 0.25*( phb(i,k1,j )+phb(i,k1+1,j ) &
5656 +phb(i,k1,j-1)+phb(i,k1+1,j-1) &
5657 +ph(i,k1,j )+ph(i,k1+1,j ) &
5658 +ph(i,k1,j-1)+ph(i,k1+1,j-1))/g
5659 z00(k1) = z
5660 k1 = k1 - 1
5661 ENDDO
5662 k1 = k1 + 2
5664 ! Get reference state at model levels
5665 DO k = k1, ktf
5666 k2 = ktf
5667 DO WHILE( z_base(k2) .gt. z00(k) )
5668 k2 = k2 - 1
5669 ENDDO
5670 if(k2+1.gt.ktf)then
5671 v00(k) = v_base(k2) + ( v_base(k2) - v_base(k2-1) ) &
5672 * ( z00(k) - z_base(k2) ) &
5673 / ( z_base(k2) - z_base(k2-1) )
5674 else
5675 v00(k) = v_base(k2) + ( v_base(k2+1) - v_base(k2) ) &
5676 * ( z00(k) - z_base(k2) ) &
5677 / ( z_base(k2+1) - z_base(k2) )
5678 endif
5679 ENDDO
5681 ! Apply the Rayleigh damper
5682 DO k = k1, ktf
5683 dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
5684 dcoef = (SIN( 0.5 * pii * dcoef ) )**2
5685 rv_tendf(i,k,j) = rv_tendf(i,k,j) - &
5686 muv(i,j) * ( dcoef * dampcoef ) * &
5687 ( v(i,k,j) - v00(k) )
5688 END DO
5690 END DO
5691 END DO
5693 ! End adjustment of v.
5694 !-----------------------------------------------------------------------
5696 !-----------------------------------------------------------------------
5697 ! Adjust w to base state.
5699 DO j = jts, MIN( jte, jde-1 )
5700 DO i = its, MIN( ite, ide-1 )
5701 ztop = ( phb(i,kde,j) + ph(i,kde,j) ) / g
5702 DO k = kts, MIN( kte, kde )
5703 z = ( phb(i,k,j) + ph(i,k,j) ) / g
5704 IF ( z >= (ztop-zdamp) ) THEN
5705 dcoef = 1.0 - MIN( 1.0, ( ztop - z ) / zdamp )
5706 dcoef = ( SIN( 0.5 * pii * dcoef ) )**2
5707 rw_tendf(i,k,j) = rw_tendf(i,k,j) - &
5708 mut(i,j) * ( dcoef * dampcoef ) * w(i,k,j)
5709 END IF
5710 END DO
5711 END DO
5712 END DO
5714 ! End adjustment of w.
5715 !-----------------------------------------------------------------------
5717 !-----------------------------------------------------------------------
5718 ! Adjust potential temperature to base state.
5720 DO j = jts, MIN( jte, jde-1 )
5721 DO i = its, MIN( ite, ide-1 )
5723 ! Get height at top of model
5724 ztop = ( phb(i,kde,j) + ph(i,kde,j) ) / g
5726 ! Find bottom of damping layer
5727 k1 = ktf
5728 z = ztop
5729 DO WHILE( z >= (ztop-zdamp) )
5730 z = 0.5 * ( phb(i,k1,j) + phb(i,k1+1,j) + &
5731 ph(i,k1,j) + ph(i,k1+1,j) ) / g
5732 z00(k1) = z
5733 k1 = k1 - 1
5734 ENDDO
5735 k1 = k1 + 2
5737 ! Get reference state at model levels
5738 DO k = k1, ktf
5739 k2 = ktf
5740 DO WHILE( z_base(k2) .gt. z00(k) )
5741 k2 = k2 - 1
5742 ENDDO
5743 if(k2+1.gt.ktf)then
5744 t00(k) = t_base(k2) + ( t_base(k2) - t_base(k2-1) ) &
5745 * ( z00(k) - z_base(k2) ) &
5746 / ( z_base(k2) - z_base(k2-1) )
5747 else
5748 t00(k) = t_base(k2) + ( t_base(k2+1) - t_base(k2) ) &
5749 * ( z00(k) - z_base(k2) ) &
5750 / ( z_base(k2+1) - z_base(k2) )
5751 endif
5752 ENDDO
5754 ! Apply the Rayleigh damper
5755 DO k = k1, ktf
5756 dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
5757 dcoef = (SIN( 0.5 * pii * dcoef ) )**2
5758 t_tendf(i,k,j) = t_tendf(i,k,j) - &
5759 mut(i,j) * ( dcoef * dampcoef ) * &
5760 ( t(i,k,j) - t00(k) )
5761 END DO
5763 END DO
5764 END DO
5766 ! End adjustment of potential temperature.
5767 !-----------------------------------------------------------------------
5769 END SUBROUTINE rk_rayleigh_damp
5771 !==============================================================================
5772 !==============================================================================
5774 SUBROUTINE theta_relaxation( t_tendf, t, t_init, &
5775 mut, ph, phb, &
5776 t_base, z_base, &
5777 ids, ide, jds, jde, kds, kde, &
5778 ims, ime, jms, jme, kms, kme, &
5779 its, ite, jts, jte, kts, kte )
5781 ! Purpose: Newtonian relaxation on potential temperature. Serves two
5782 ! purposes: 1) to mimic atmospheric radiation in a simple
5783 ! manner, and 2) to keep the vertical profile of temperature
5784 ! close to the initial (base-state) profile, which is useful
5785 ! for certain idealized applications.
5787 ! Reference: Rotunno and Emanuel, 1987, JAS, p. 546
5789 !-----------------------------------------------------------------------
5790 ! Begin declarations.
5792 IMPLICIT NONE
5794 INTEGER, INTENT( IN ) &
5795 :: ids, ide, jds, jde, kds, kde, &
5796 ims, ime, jms, jme, kms, kme, &
5797 its, ite, jts, jte, kts, kte
5799 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
5800 :: t_tendf
5802 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
5803 :: t, t_init, ph, phb
5805 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
5806 :: mut
5808 REAL, DIMENSION( kms:kme ) , INTENT(IN ) &
5809 :: t_base, z_base
5811 ! Local variables.
5813 INTEGER :: i, j, k, ktf, k2
5814 REAL :: tau_r , rmax , rmin , inv_tau_r , inv_g , rterm
5815 REAL, DIMENSION( kms:kme ) :: z00,t00
5817 ! End declarations.
5818 !-----------------------------------------------------------------------
5820 ! set tau_r to 12 h, following RE87
5821 tau_r = 12.0*3600.0
5823 ! limit rterm to +/- 2 K/day
5824 rmax = 2.0/86400.0
5825 rmin = -rmax
5827 ktf = MIN( kte, kde-1 )
5828 inv_tau_r = 1.0/tau_r
5829 inv_g = 1.0/g
5831 !-----------------------------------------------------------------------
5832 ! Adjust potential temperature to base state.
5834 DO j = jts, MIN( jte, jde-1 )
5835 DO i = its, MIN( ite, ide-1 )
5837 ! Get height of model levels:
5838 DO k = kts, ktf
5839 z00(k) = 0.5 * ( phb(i,k,j) + phb(i,k+1,j) + &
5840 ph(i,k,j) + ph(i,k+1,j) ) * inv_g
5841 ENDDO
5843 ! Get reference state:
5844 DO k = kts, ktf
5845 k2 = ktf
5846 DO WHILE( z_base(k2) .gt. z00(k) .and. k2 .gt. 1 )
5847 k2 = k2 - 1
5848 ENDDO
5849 if(k2+1.gt.ktf)then
5850 t00(k) = t_base(k2) + ( t_base(k2) - t_base(k2-1) ) &
5851 * ( z00(k) - z_base(k2) ) &
5852 / ( z_base(k2) - z_base(k2-1) )
5853 else
5854 t00(k) = t_base(k2) + ( t_base(k2+1) - t_base(k2) ) &
5855 * ( z00(k) - z_base(k2) ) &
5856 / ( z_base(k2+1) - z_base(k2) )
5857 endif
5858 ENDDO
5860 ! Apply the RE87 R term:
5861 DO k = kts, ktf
5862 rterm = -( t(i,k,j) - t00(k) )*inv_tau_r
5863 ! limit rterm:
5864 rterm = min( rterm , rmax )
5865 rterm = max( rterm , rmin )
5866 t_tendf(i,k,j) = t_tendf(i,k,j) + mut(i,j)*rterm
5867 END DO
5869 END DO
5870 END DO
5872 END SUBROUTINE theta_relaxation
5874 !==============================================================================
5875 !==============================================================================
5877 SUBROUTINE sixth_order_diffusion( name, field, tendency, mu, dt, &
5878 config_flags, &
5879 diff_6th_opt, diff_6th_factor, &
5880 ids, ide, jds, jde, kds, kde, &
5881 ims, ime, jms, jme, kms, kme, &
5882 its, ite, jts, jte, kts, kte )
5884 ! History: 14 Nov 2006 Name of variable changed by Jason Knievel
5885 ! 07 Jun 2006 Revised and generalized by Jason Knievel
5886 ! 25 Apr 2005 Original code by Jason Knievel, NCAR
5888 ! Purpose: Apply 6th-order, monotonic (flux-limited), numerical
5889 ! diffusion to 3-d velocity and to scalars.
5891 ! References: Ming Xue (MWR Aug 2000)
5892 ! Durran ("Numerical Methods for Wave Equations..." 1999)
5893 ! George Bryan (personal communication)
5895 !------------------------------------------------------------------------------
5896 ! Begin: Declarations.
5898 IMPLICIT NONE
5900 INTEGER, INTENT(IN) &
5901 :: ids, ide, jds, jde, kds, kde, &
5902 ims, ime, jms, jme, kms, kme, &
5903 its, ite, jts, jte, kts, kte
5905 TYPE(grid_config_rec_type), INTENT(IN) &
5906 :: config_flags
5908 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) &
5909 :: tendency
5911 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) &
5912 :: field
5914 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN) &
5915 :: mu
5917 REAL, INTENT(IN) &
5918 :: dt
5920 REAL, INTENT(IN) &
5921 :: diff_6th_factor
5923 INTEGER, INTENT(IN) &
5924 :: diff_6th_opt
5926 CHARACTER(LEN=1) , INTENT(IN) &
5927 :: name
5929 INTEGER &
5930 :: i, j, k, &
5931 i_start, i_end, &
5932 j_start, j_end, &
5933 k_start, k_end, &
5934 ktf
5936 REAL &
5937 :: dflux_x_p0, dflux_y_p0, &
5938 dflux_x_p1, dflux_y_p1, &
5939 tendency_x, tendency_y, &
5940 mu_avg_p0, mu_avg_p1, &
5941 diff_6th_coef
5943 LOGICAL &
5944 :: specified
5946 ! End: Declarations.
5947 !------------------------------------------------------------------------------
5949 !------------------------------------------------------------------------------
5950 ! Begin: Translate the diffusion factor into a diffusion coefficient. See
5951 ! Durran's text, section 2.4.3, then adjust for sixth-order diffusion (not
5952 ! fourth) and for diffusion in two dimensions (not one). For reference, a
5953 ! factor of 1.0 would mean complete diffusion of a 2dx wave in one time step,
5954 ! although application of the flux limiter reduces somewhat the effects of
5955 ! diffusion for a given coefficient.
5957 diff_6th_coef = diff_6th_factor * 0.015625 / ( 2.0 * dt )
5959 ! End: Translate diffusion factor.
5960 !------------------------------------------------------------------------------
5962 !------------------------------------------------------------------------------
5963 ! Begin: Assign limits of spatial loops depending on variable to be diffused.
5964 ! The halo regions are already filled with values by the time this subroutine
5965 ! is called, which allows the stencil to extend beyond the domains' edges.
5967 ktf = MIN( kte, kde-1 )
5969 IF ( name .EQ. 'u' ) THEN
5971 i_start = its
5972 i_end = ite
5973 j_start = jts
5974 j_end = MIN(jde-1,jte)
5975 k_start = kts
5976 k_end = ktf
5978 ELSE IF ( name .EQ. 'v' ) THEN
5980 i_start = its
5981 i_end = MIN(ide-1,ite)
5982 j_start = jts
5983 j_end = jte
5984 k_start = kts
5985 k_end = ktf
5987 ELSE IF ( name .EQ. 'w' ) THEN
5989 i_start = its
5990 i_end = MIN(ide-1,ite)
5991 j_start = jts
5992 j_end = MIN(jde-1,jte)
5993 k_start = kts+1
5994 k_end = ktf
5996 ELSE
5998 i_start = its
5999 i_end = MIN(ide-1,ite)
6000 j_start = jts
6001 j_end = MIN(jde-1,jte)
6002 k_start = kts
6003 k_end = ktf
6005 ENDIF
6007 ! End: Assignment of limits of spatial loops.
6008 !------------------------------------------------------------------------------
6010 !------------------------------------------------------------------------------
6011 ! Begin: Loop across spatial dimensions.
6013 DO j = j_start, j_end
6014 DO k = k_start, k_end
6015 DO i = i_start, i_end
6017 !------------------------------------------------------------------------------
6018 ! Begin: Diffusion in x (i index).
6020 ! Calculate the diffusive flux in x direction (from Xue's eq. 3).
6022 dflux_x_p0 = ( 10.0 * ( field(i, k,j) - field(i-1,k,j) ) &
6023 - 5.0 * ( field(i+1,k,j) - field(i-2,k,j) ) &
6024 + ( field(i+2,k,j) - field(i-3,k,j) ) )
6026 dflux_x_p1 = ( 10.0 * ( field(i+1,k,j) - field(i ,k,j) ) &
6027 - 5.0 * ( field(i+2,k,j) - field(i-1,k,j) ) &
6028 + ( field(i+3,k,j) - field(i-2,k,j) ) )
6030 ! If requested in the namelist (diff_6th_opt=2), prohibit up-gradient diffusion
6031 ! (variation on Xue's eq. 10).
6033 IF ( diff_6th_opt .EQ. 2 ) THEN
6035 IF ( dflux_x_p0 * ( field(i ,k,j)-field(i-1,k,j) ) .LE. 0.0 ) THEN
6036 dflux_x_p0 = 0.0
6037 END IF
6039 IF ( dflux_x_p1 * ( field(i+1,k,j)-field(i ,k,j) ) .LE. 0.0 ) THEN
6040 dflux_x_p1 = 0.0
6041 END IF
6043 END IF
6045 ! Apply 6th-order diffusion in x direction.
6047 IF ( name .EQ. 'u' ) THEN
6048 mu_avg_p0 = mu(i-1,j)
6049 mu_avg_p1 = mu(i ,j)
6050 ELSE IF ( name .EQ. 'v' ) THEN
6051 mu_avg_p0 = 0.25 * ( &
6052 mu(i-1,j-1) + &
6053 mu(i ,j-1) + &
6054 mu(i-1,j ) + &
6055 mu(i ,j ) )
6056 mu_avg_p1 = 0.25 * ( &
6057 mu(i ,j-1) + &
6058 mu(i+1,j-1) + &
6059 mu(i ,j ) + &
6060 mu(i+1,j ) )
6061 ELSE
6062 mu_avg_p0 = 0.5 * ( &
6063 mu(i-1,j) + &
6064 mu(i ,j) )
6065 mu_avg_p1 = 0.5 * ( &
6066 mu(i ,j) + &
6067 mu(i+1,j) )
6068 END IF
6070 tendency_x = diff_6th_coef * &
6071 ( ( mu_avg_p1 * dflux_x_p1 ) - ( mu_avg_p0 * dflux_x_p0 ) )
6073 ! End: Diffusion in x.
6074 !------------------------------------------------------------------------------
6076 !------------------------------------------------------------------------------
6077 ! Begin: Diffusion in y (j index).
6079 ! Calculate the diffusive flux in y direction (from Xue's eq. 3).
6081 dflux_y_p0 = ( 10.0 * ( field(i,k,j ) - field(i,k,j-1) ) &
6082 - 5.0 * ( field(i,k,j+1) - field(i,k,j-2) ) &
6083 + ( field(i,k,j+2) - field(i,k,j-3) ) )
6085 dflux_y_p1 = ( 10.0 * ( field(i,k,j+1) - field(i,k,j ) ) &
6086 - 5.0 * ( field(i,k,j+2) - field(i,k,j-1) ) &
6087 + ( field(i,k,j+3) - field(i,k,j-2) ) )
6089 ! If requested in the namelist (diff_6th_opt=2), prohibit up-gradient diffusion
6090 ! (variation on Xue's eq. 10).
6092 IF ( diff_6th_opt .EQ. 2 ) THEN
6094 IF ( dflux_y_p0 * ( field(i,k,j )-field(i,k,j-1) ) .LE. 0.0 ) THEN
6095 dflux_y_p0 = 0.0
6096 END IF
6098 IF ( dflux_y_p1 * ( field(i,k,j+1)-field(i,k,j ) ) .LE. 0.0 ) THEN
6099 dflux_y_p1 = 0.0
6100 END IF
6102 END IF
6104 ! Apply 6th-order diffusion in y direction.
6106 IF ( name .EQ. 'u' ) THEN
6107 mu_avg_p0 = 0.25 * ( &
6108 mu(i-1,j-1) + &
6109 mu(i ,j-1) + &
6110 mu(i-1,j ) + &
6111 mu(i ,j ) )
6112 mu_avg_p1 = 0.25 * ( &
6113 mu(i-1,j ) + &
6114 mu(i ,j ) + &
6115 mu(i-1,j+1) + &
6116 mu(i ,j+1) )
6117 ELSE IF ( name .EQ. 'v' ) THEN
6118 mu_avg_p0 = mu(i,j-1)
6119 mu_avg_p1 = mu(i,j )
6120 ELSE
6121 mu_avg_p0 = 0.5 * ( &
6122 mu(i,j-1) + &
6123 mu(i,j ) )
6124 mu_avg_p1 = 0.5 * ( &
6125 mu(i,j ) + &
6126 mu(i,j+1) )
6127 END IF
6129 tendency_y = diff_6th_coef * &
6130 ( ( mu_avg_p1 * dflux_y_p1 ) - ( mu_avg_p0 * dflux_y_p0 ) )
6132 ! End: Diffusion in y.
6133 !------------------------------------------------------------------------------
6135 !------------------------------------------------------------------------------
6136 ! Begin: Combine diffusion in x and y.
6138 tendency(i,k,j) = tendency(i,k,j) + tendency_x + tendency_y
6140 ! End: Combine diffusion in x and y.
6141 !------------------------------------------------------------------------------
6143 ENDDO
6144 ENDDO
6145 ENDDO
6147 ! End: Loop across spatial dimensions.
6148 !------------------------------------------------------------------------------
6150 END SUBROUTINE sixth_order_diffusion
6152 !==============================================================================
6153 !==============================================================================
6155 END MODULE module_big_step_utilities_em