| blob | 1a3cb6ae30883fc71c6a7157a6a1804d2032cd8c |
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_domain, ONLY : domain
16 USE module_model_constants
17 USE module_state_description
18 USE module_configure, ONLY : grid_config_rec_type
19 USE module_wrf_error
21 CONTAINS
23 !-------------------------------------------------------------------------------
25 SUBROUTINE calc_mu_uv ( config_flags, &
26 mu, mub, muu, muv, &
27 ids, ide, jds, jde, kds, kde, &
28 ims, ime, jms, jme, kms, kme, &
29 its, ite, jts, jte, kts, kte )
31 IMPLICIT NONE
33 ! Input data
35 TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
37 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
38 ims, ime, jms, jme, kms, kme, &
39 its, ite, jts, jte, kts, kte
41 REAL, DIMENSION( ims:ime , jms:jme ) , INTENT( OUT) :: muu, muv
42 REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu, mub
44 ! local stuff
46 INTEGER :: i, j, itf, jtf, im, jm
48 !<DESCRIPTION>
49 !
50 ! calc_mu_uv calculates the full column dry-air mass at the staggered
51 ! horizontal velocity points (u,v) and places the results in muu and muv.
52 ! This routine uses the reference state (mub) and perturbation state (mu)
53 !
54 !</DESCRIPTION>
57 itf=ite
58 jtf=MIN(jte,jde-1)
60 IF ( ( its .NE. ids ) .AND. ( ite .NE. ide ) ) THEN
61 DO j=jts,jtf
62 DO i=its,itf
63 muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
64 ENDDO
65 ENDDO
66 ELSE IF ( ( its .EQ. ids ) .AND. ( ite .NE. ide ) ) THEN
67 DO j=jts,jtf
68 DO i=its+1,itf
69 muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
70 ENDDO
71 ENDDO
72 i=its
73 im = its
74 if(config_flags%periodic_x) im = its-1
75 DO j=jts,jtf
76 ! muu(i,j) = mu(i,j) +mub(i,j)
77 ! fix for periodic b.c., 13 march 2004, wcs
78 muu(i,j) = 0.5*(mu(i,j)+mu(im,j)+mub(i,j)+mub(im,j))
79 ENDDO
80 ELSE IF ( ( its .NE. ids ) .AND. ( ite .EQ. ide ) ) THEN
81 DO j=jts,jtf
82 DO i=its,itf-1
83 muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
84 ENDDO
85 ENDDO
86 i=ite
87 im = ite-1
88 if(config_flags%periodic_x) im = ite
89 DO j=jts,jtf
90 ! muu(i,j) = mu(i-1,j) +mub(i-1,j)
91 ! fix for periodic b.c., 13 march 2004, wcs
92 muu(i,j) = 0.5*(mu(i-1,j)+mu(im,j)+mub(i-1,j)+mub(im,j))
93 ENDDO
94 ELSE IF ( ( its .EQ. ids ) .AND. ( ite .EQ. ide ) ) THEN
95 DO j=jts,jtf
96 DO i=its+1,itf-1
97 muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
98 ENDDO
99 ENDDO
100 i=its
101 im = its
102 if(config_flags%periodic_x) im = its-1
103 DO j=jts,jtf
104 ! muu(i,j) = mu(i,j) +mub(i,j)
105 ! fix for periodic b.c., 13 march 2004, wcs
106 muu(i,j) = 0.5*(mu(i,j)+mu(im,j)+mub(i,j)+mub(im,j))
107 ENDDO
108 i=ite
109 im = ite-1
110 if(config_flags%periodic_x) im = ite
111 DO j=jts,jtf
112 ! muu(i,j) = mu(i-1,j) +mub(i-1,j)
113 ! fix for periodic b.c., 13 march 2004, wcs
114 muu(i,j) = 0.5*(mu(i-1,j)+mu(im,j)+mub(i-1,j)+mub(im,j))
115 ENDDO
116 END IF
118 itf=MIN(ite,ide-1)
119 jtf=jte
121 IF ( ( jts .NE. jds ) .AND. ( jte .NE. jde ) ) THEN
122 DO j=jts,jtf
123 DO i=its,itf
124 muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
125 ENDDO
126 ENDDO
127 ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .NE. jde ) ) THEN
128 DO j=jts+1,jtf
129 DO i=its,itf
130 muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
131 ENDDO
132 ENDDO
133 j=jts
134 jm = jts
135 if(config_flags%periodic_y) jm = jts-1
136 DO i=its,itf
137 ! muv(i,j) = mu(i,j) +mub(i,j)
138 ! fix for periodic b.c., 13 march 2004, wcs
139 muv(i,j) = 0.5*(mu(i,j)+mu(i,jm)+mub(i,j)+mub(i,jm))
140 ENDDO
141 ELSE IF ( ( jts .NE. jds ) .AND. ( jte .EQ. jde ) ) THEN
142 DO j=jts,jtf-1
143 DO i=its,itf
144 muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
145 ENDDO
146 ENDDO
147 j=jte
148 jm = jte-1
149 if(config_flags%periodic_y) jm = jte
150 DO i=its,itf
151 muv(i,j) = mu(i,j-1) +mub(i,j-1)
152 ! fix for periodic b.c., 13 march 2004, wcs
153 muv(i,j) = 0.5*(mu(i,j-1)+mu(i,jm)+mub(i,j-1)+mub(i,jm))
154 ENDDO
155 ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .EQ. jde ) ) THEN
156 DO j=jts+1,jtf-1
157 DO i=its,itf
158 muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
159 ENDDO
160 ENDDO
161 j=jts
162 jm = jts
163 if(config_flags%periodic_y) jm = jts-1
164 DO i=its,itf
165 ! muv(i,j) = mu(i,j) +mub(i,j)
166 ! fix for periodic b.c., 13 march 2004, wcs
167 muv(i,j) = 0.5*(mu(i,j)+mu(i,jm)+mub(i,j)+mub(i,jm))
168 ENDDO
169 j=jte
170 jm = jte-1
171 if(config_flags%periodic_y) jm = jte
172 DO i=its,itf
173 ! muv(i,j) = mu(i,j-1) +mub(i,j-1)
174 ! fix for periodic b.c., 13 march 2004, wcs
175 muv(i,j) = 0.5*(mu(i,j-1)+mu(i,jm)+mub(i,j-1)+mub(i,jm))
176 ENDDO
177 END IF
179 END SUBROUTINE calc_mu_uv
181 !-------------------------------------------------------------------------------
183 SUBROUTINE calc_mu_uv_1 ( config_flags, &
184 mu, muu, muv, &
185 ids, ide, jds, jde, kds, kde, &
186 ims, ime, jms, jme, kms, kme, &
187 its, ite, jts, jte, kts, kte )
189 IMPLICIT NONE
191 ! Input data
193 TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
195 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
196 ims, ime, jms, jme, kms, kme, &
197 its, ite, jts, jte, kts, kte
199 REAL, DIMENSION( ims:ime , jms:jme ) , INTENT( OUT) :: muu, muv
200 REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
202 ! local stuff
204 INTEGER :: i, j, itf, jtf, im, jm
206 !<DESCRIPTION>
207 !
208 ! calc_mu_uv calculates the full column dry-air mass at the staggered
209 ! horizontal velocity points (u,v) and places the results in muu and muv.
210 ! This routine uses the full state (mu)
211 !
212 !</DESCRIPTION>
214 itf=ite
215 jtf=MIN(jte,jde-1)
217 IF ( ( its .NE. ids ) .AND. ( ite .NE. ide ) ) THEN
218 DO j=jts,jtf
219 DO i=its,itf
220 muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j))
221 ENDDO
222 ENDDO
223 ELSE IF ( ( its .EQ. ids ) .AND. ( ite .NE. ide ) ) THEN
224 DO j=jts,jtf
225 DO i=its+1,itf
226 muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j))
227 ENDDO
228 ENDDO
229 i=its
230 im = its
231 if(config_flags%periodic_x) im = its-1
232 DO j=jts,jtf
233 muu(i,j) = 0.5*(mu(i,j)+mu(im,j))
234 ENDDO
235 ELSE IF ( ( its .NE. ids ) .AND. ( ite .EQ. ide ) ) THEN
236 DO j=jts,jtf
237 DO i=its,itf-1
238 muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j))
239 ENDDO
240 ENDDO
241 i=ite
242 im = ite-1
243 if(config_flags%periodic_x) im = ite
244 DO j=jts,jtf
245 muu(i,j) = 0.5*(mu(i-1,j)+mu(im,j))
246 ENDDO
247 ELSE IF ( ( its .EQ. ids ) .AND. ( ite .EQ. ide ) ) THEN
248 DO j=jts,jtf
249 DO i=its+1,itf-1
250 muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j))
251 ENDDO
252 ENDDO
253 i=its
254 im = its
255 if(config_flags%periodic_x) im = its-1
256 DO j=jts,jtf
257 muu(i,j) = 0.5*(mu(i,j)+mu(im,j))
258 ENDDO
259 i=ite
260 im = ite-1
261 if(config_flags%periodic_x) im = ite
262 DO j=jts,jtf
263 muu(i,j) = 0.5*(mu(i-1,j)+mu(im,j))
264 ENDDO
265 END IF
267 itf=MIN(ite,ide-1)
268 jtf=jte
270 IF ( ( jts .NE. jds ) .AND. ( jte .NE. jde ) ) THEN
271 DO j=jts,jtf
272 DO i=its,itf
273 muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1))
274 ENDDO
275 ENDDO
276 ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .NE. jde ) ) THEN
277 DO j=jts+1,jtf
278 DO i=its,itf
279 muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1))
280 ENDDO
281 ENDDO
282 j=jts
283 jm = jts
284 if(config_flags%periodic_y) jm = jts-1
285 DO i=its,itf
286 muv(i,j) = 0.5*(mu(i,j)+mu(i,jm))
287 ENDDO
288 ELSE IF ( ( jts .NE. jds ) .AND. ( jte .EQ. jde ) ) THEN
289 DO j=jts,jtf-1
290 DO i=its,itf
291 muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1))
292 ENDDO
293 ENDDO
294 j=jte
295 jm = jte-1
296 if(config_flags%periodic_y) jm = jte
297 DO i=its,itf
298 muv(i,j) = 0.5*(mu(i,j-1)+mu(i,jm))
299 ENDDO
300 ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .EQ. jde ) ) THEN
301 DO j=jts+1,jtf-1
302 DO i=its,itf
303 muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1))
304 ENDDO
305 ENDDO
306 j=jts
307 jm = jts
308 if(config_flags%periodic_y) jm = jts-1
309 DO i=its,itf
310 muv(i,j) = 0.5*(mu(i,j)+mu(i,jm))
311 ENDDO
312 j=jte
313 jm = jte-1
314 if(config_flags%periodic_y) jm = jte
315 DO i=its,itf
316 muv(i,j) = 0.5*(mu(i,j-1)+mu(i,jm))
317 ENDDO
318 END IF
320 END SUBROUTINE calc_mu_uv_1
322 !-------------------------------------------------------------------------------
324 ! Map scale factor comments for this routine:
325 ! Locally not changed, but sent the correct map scale factors
326 ! from module_em (msfuy, msfvx, msfty)
328 SUBROUTINE couple_momentum ( muu, ru, u, msfu, &
329 muv, rv, v, msfv, msfv_inv, &
330 mut, rw, w, msft, &
331 ids, ide, jds, jde, kds, kde, &
332 ims, ime, jms, jme, kms, kme, &
333 its, ite, jts, jte, kts, kte )
335 IMPLICIT NONE
337 ! Input data
339 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
340 ims, ime, jms, jme, kms, kme, &
341 its, ite, jts, jte, kts, kte
343 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT( OUT) :: ru, rv, rw
345 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muu, muv, mut
346 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfu, msfv, msft, msfv_inv
348 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: u, v, w
350 ! Local data
352 INTEGER :: i, j, k, itf, jtf, ktf
354 !<DESCRIPTION>
355 !
356 ! couple_momentum couples the velocities to the full column mass and
357 ! the map factors.
358 !
359 !</DESCRIPTION>
361 ktf=MIN(kte,kde-1)
363 itf=ite
364 jtf=MIN(jte,jde-1)
366 DO j=jts,jtf
367 DO k=kts,ktf
368 DO i=its,itf
369 ru(i,k,j)=u(i,k,j)*muu(i,j)/msfu(i,j)
370 ENDDO
371 ENDDO
372 ENDDO
374 itf=MIN(ite,ide-1)
375 jtf=jte
377 DO j=jts,jtf
378 DO k=kts,ktf
379 DO i=its,itf
380 rv(i,k,j)=v(i,k,j)*muv(i,j)*msfv_inv(i,j)
381 ! rv(i,k,j)=v(i,k,j)*muv(i,j)/msfv(i,j)
382 ENDDO
383 ENDDO
384 ENDDO
386 itf=MIN(ite,ide-1)
387 jtf=MIN(jte,jde-1)
389 DO j=jts,jtf
390 DO k=kts,kte
391 DO i=its,itf
392 rw(i,k,j)=w(i,k,j)*mut(i,j)/msft(i,j)
393 ENDDO
394 ENDDO
395 ENDDO
397 END SUBROUTINE couple_momentum
399 !-------------------------------------------------------------------
401 SUBROUTINE calc_mu_staggered ( mu, mub, muu, muv, &
402 ids, ide, jds, jde, kds, kde, &
403 ims, ime, jms, jme, kms, kme, &
404 its, ite, jts, jte, kts, kte )
406 IMPLICIT NONE
408 ! Input data
410 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
411 ims, ime, jms, jme, kms, kme, &
412 its, ite, jts, jte, kts, kte
414 REAL, DIMENSION( ims:ime , jms:jme ) , INTENT( OUT) :: muu, muv
415 REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu, mub
417 ! local stuff
419 INTEGER :: i, j, itf, jtf
421 !<DESCRIPTION>
422 !
423 ! calc_mu_staggered calculates the full dry air mass at the staggered
424 ! velocity points (u,v).
425 !
426 !</DESCRIPTION>
428 itf=ite
429 jtf=MIN(jte,jde-1)
431 IF ( ( its .NE. ids ) .AND. ( ite .NE. ide ) ) THEN
432 DO j=jts,jtf
433 DO i=its,itf
434 muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
435 ENDDO
436 ENDDO
437 ELSE IF ( ( its .EQ. ids ) .AND. ( ite .NE. ide ) ) THEN
438 DO j=jts,jtf
439 DO i=its+1,itf
440 muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
441 ENDDO
442 ENDDO
443 i=its
444 DO j=jts,jtf
445 muu(i,j) = mu(i,j) +mub(i,j)
446 ENDDO
447 ELSE IF ( ( its .NE. ids ) .AND. ( ite .EQ. ide ) ) THEN
448 DO j=jts,jtf
449 DO i=its,itf-1
450 muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
451 ENDDO
452 ENDDO
453 i=ite
454 DO j=jts,jtf
455 muu(i,j) = mu(i-1,j) +mub(i-1,j)
456 ENDDO
457 ELSE IF ( ( its .EQ. ids ) .AND. ( ite .EQ. ide ) ) THEN
458 DO j=jts,jtf
459 DO i=its+1,itf-1
460 muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
461 ENDDO
462 ENDDO
463 i=its
464 DO j=jts,jtf
465 muu(i,j) = mu(i,j) +mub(i,j)
466 ENDDO
467 i=ite
468 DO j=jts,jtf
469 muu(i,j) = mu(i-1,j) +mub(i-1,j)
470 ENDDO
471 END IF
473 itf=MIN(ite,ide-1)
474 jtf=jte
476 IF ( ( jts .NE. jds ) .AND. ( jte .NE. jde ) ) THEN
477 DO j=jts,jtf
478 DO i=its,itf
479 muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
480 ENDDO
481 ENDDO
482 ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .NE. jde ) ) THEN
483 DO j=jts+1,jtf
484 DO i=its,itf
485 muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
486 ENDDO
487 ENDDO
488 j=jts
489 DO i=its,itf
490 muv(i,j) = mu(i,j) +mub(i,j)
491 ENDDO
492 ELSE IF ( ( jts .NE. jds ) .AND. ( jte .EQ. jde ) ) THEN
493 DO j=jts,jtf-1
494 DO i=its,itf
495 muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
496 ENDDO
497 ENDDO
498 j=jte
499 DO i=its,itf
500 muv(i,j) = mu(i,j-1) +mub(i,j-1)
501 ENDDO
502 ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .EQ. jde ) ) THEN
503 DO j=jts+1,jtf-1
504 DO i=its,itf
505 muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
506 ENDDO
507 ENDDO
508 j=jts
509 DO i=its,itf
510 muv(i,j) = mu(i,j) +mub(i,j)
511 ENDDO
512 j=jte
513 DO i=its,itf
514 muv(i,j) = mu(i,j-1) +mub(i,j-1)
515 ENDDO
516 END IF
518 END SUBROUTINE calc_mu_staggered
520 !-------------------------------------------------------------------------------
522 SUBROUTINE couple ( mu, mub, rfield, field, name, &
523 msf, &
524 ids, ide, jds, jde, kds, kde, &
525 ims, ime, jms, jme, kms, kme, &
526 its, ite, jts, jte, kts, kte )
528 IMPLICIT NONE
530 ! Input data
532 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
533 ims, ime, jms, jme, kms, kme, &
534 its, ite, jts, jte, kts, kte
536 CHARACTER(LEN=1) , INTENT(IN ) :: name
538 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT( OUT) :: rfield
540 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu, mub, msf
542 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field
544 ! Local data
546 INTEGER :: i, j, k, itf, jtf, ktf
547 REAL , DIMENSION(ims:ime,jms:jme) :: muu , muv
549 !<DESCRIPTION>
550 !
551 ! subroutine couple couples the input variable with the dry-air
552 ! column mass (mu).
553 !
554 !</DESCRIPTION>
557 ktf=MIN(kte,kde-1)
559 IF (name .EQ. 'u')THEN
561 CALL calc_mu_staggered ( mu, mub, muu, muv, &
562 ids, ide, jds, jde, kds, kde, &
563 ims, ime, jms, jme, kms, kme, &
564 its, ite, jts, jte, kts, kte )
566 itf=ite
567 jtf=MIN(jte,jde-1)
569 DO j=jts,jtf
570 DO k=kts,ktf
571 DO i=its,itf
572 rfield(i,k,j)=field(i,k,j)*muu(i,j)/msf(i,j)
573 ENDDO
574 ENDDO
575 ENDDO
577 ELSE IF (name .EQ. 'v')THEN
579 CALL calc_mu_staggered ( mu, mub, muu, muv, &
580 ids, ide, jds, jde, kds, kde, &
581 ims, ime, jms, jme, kms, kme, &
582 its, ite, jts, jte, kts, kte )
584 itf=ite
585 itf=MIN(ite,ide-1)
586 jtf=jte
588 DO j=jts,jtf
589 DO k=kts,ktf
590 DO i=its,itf
591 rfield(i,k,j)=field(i,k,j)*muv(i,j)/msf(i,j)
592 ENDDO
593 ENDDO
594 ENDDO
596 ELSE IF (name .EQ. 'w')THEN
597 itf=MIN(ite,ide-1)
598 jtf=MIN(jte,jde-1)
599 DO j=jts,jtf
600 DO k=kts,kte
601 DO i=its,itf
602 rfield(i,k,j)=field(i,k,j)*(mu(i,j)+mub(i,j))/msf(i,j)
603 ENDDO
604 ENDDO
605 ENDDO
607 ELSE IF (name .EQ. 'h')THEN
608 itf=MIN(ite,ide-1)
609 jtf=MIN(jte,jde-1)
610 DO j=jts,jtf
611 DO k=kts,kte
612 DO i=its,itf
613 rfield(i,k,j)=field(i,k,j)*(mu(i,j)+mub(i,j))
614 ENDDO
615 ENDDO
616 ENDDO
618 ELSE
619 itf=MIN(ite,ide-1)
620 jtf=MIN(jte,jde-1)
621 DO j=jts,jtf
622 DO k=kts,ktf
623 DO i=its,itf
624 rfield(i,k,j)=field(i,k,j)*(mu(i,j)+mub(i,j))
625 ENDDO
626 ENDDO
627 ENDDO
629 ENDIF
631 END SUBROUTINE couple
634 !-------------------------------------------------------------------------------
636 SUBROUTINE calc_ww_cp ( u, v, mup, mub, ww, &
637 rdx, rdy, msftx, msfty, &
638 msfux, msfuy, msfvx, msfvx_inv, &
639 msfvy, dnw, &
640 ids, ide, jds, jde, kds, kde, &
641 ims, ime, jms, jme, kms, kme, &
642 its, ite, jts, jte, kts, kte )
644 IMPLICIT NONE
646 ! Input data
649 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
650 ims, ime, jms, jme, kms, kme, &
651 its, ite, jts, jte, kts, kte
653 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: u, v
654 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mup, mub, &
655 msftx, msfty, &
656 msfux, msfuy, &
657 msfvx, msfvy, &
658 msfvx_inv
659 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
661 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: ww
662 REAL , INTENT(IN ) :: rdx, rdy
664 ! Local data
666 INTEGER :: i, j, k, itf, jtf, ktf
667 REAL , DIMENSION( its:ite ) :: dmdt
668 REAL , DIMENSION( its:ite, kts:kte ) :: divv
669 REAL , DIMENSION( its:ite+1, jts:jte+1 ) :: muu, muv
671 !<DESCRIPTION>
672 !
673 ! calc_ww calculates omega using the velocities (u,v) and the dry-air
674 ! column mass (mup+mub).
675 ! The algorithm integrates the continuity equation through the column
676 ! followed by a diagnosis of omega.
677 !
678 !</DESCRIPTION>
680 !<DESCRIPTION>
681 !
682 ! calc_ww_cp calculates omega using the velocities (u,v) and the
683 ! column mass mu.
684 !
685 !</DESCRIPTION>
687 jtf=MIN(jte,jde-1)
688 ktf=MIN(kte,kde-1)
689 itf=MIN(ite,ide-1)
691 ! mu coupled with the appropriate map factor
693 DO j=jts,jtf
694 DO i=its,min(ite+1,ide)
695 ! u is always coupled with my
696 muu(i,j) = 0.5*(mup(i,j)+mub(i,j)+mup(i-1,j)+mub(i-1,j))/msfuy(i,j)
697 ENDDO
698 ENDDO
700 DO j=jts,min(jte+1,jde)
701 DO i=its,itf
702 ! v is always coupled with mx
703 ! muv(i,j) = 0.5*(mup(i,j)+mub(i,j)+mup(i,j-1)+mub(i,j-1))/msfvx(i,j)
704 muv(i,j) = 0.5*(mup(i,j)+mub(i,j)+mup(i,j-1)+mub(i,j-1))*msfvx_inv(i,j)
705 ENDDO
706 ENDDO
708 DO j=jts,jtf
710 DO i=its,ite
711 dmdt(i) = 0.
712 ww(i,1,j) = 0.
713 ww(i,kte,j) = 0.
714 ENDDO
716 ! Comments on the modifications for map scale factors
717 ! ADT eqn 47 / my (putting rho -> 'mu') is:
718 ! (1/my) partial d mu/dt = -mx partial d/dx(mu u/my)
719 ! -mx partial d/dy(mu v/mx)
720 ! -partial d/dz(mu w/my)
721 !
722 ! Using nu instead of z the last term becomes:
723 ! -partial d/dnu(mu (dnu/dt)/my)
724 !
725 ! Integrating with respect to nu over ALL levels, with dnu/dt=0 at top
726 ! and bottom, the last term becomes = 0
727 !
728 ! Integral|bot->top[(1/my) partial d mu/dt]dnu =
729 ! Integral|bot->top[-mx partial d/dx(mu u/my)
730 ! -mx partial d/dy(mu v/mx)]dnu
731 !
732 ! muu='mu'[on u]/my, muv='mu'[on v]/mx
733 ! (1/my) partial d mu/dt is independent of nu
734 ! => LHS = Integral|bot->top[con]dnu = conservation*(-1) = -dmdt
735 !
736 ! => dmdt = mx*Integral|bot->top[partial d/dx(mu u/my) +
737 ! partial d/dy(mu v/mx)]dnu
738 ! => dmdt = sum_bot->top[divv]
739 ! where
740 ! divv=mx*[partial d/dx(mu u/my) + partial d/dy(mu v/mx)]*delta nu
742 DO k=kts,ktf
743 DO i=its,itf
745 divv(i,k) = msftx(i,j)*dnw(k)*( rdx*(muu(i+1,j)*u(i+1,k,j)-muu(i,j)*u(i,k,j)) &
746 +rdy*(muv(i,j+1)*v(i,k,j+1)-muv(i,j)*v(i,k,j)) )
748 ! dmdt(i) = dmdt(i) + dnw(k)* ( rdx*(ru(i+1,k,j)-ru(i,k,j)) &
749 ! +rdy*(rv(i,k,j+1)-rv(i,k,j)) )
751 dmdt(i) = dmdt(i) + divv(i,k)
754 ENDDO
755 ENDDO
757 ! Further map scale factor notes:
758 ! Now integrate from bottom to top, level by level:
759 ! mu dnu/dt/my [k+1] = mu dnu/dt/my [k] + [-(1/my) partial d mu/dt
760 ! -mx partial d/dx(mu u/my)
761 ! -mx partial d/dy(mu v/mx)]*dnu[k->k+1]
762 ! ww [k+1] = ww [k] -(1/my) partial d mu/dt * dnu[k->k+1] - divv[k]
763 ! = ww [k] -dmdt * dnw[k] - divv[k]
765 DO k=2,ktf
766 DO i=its,itf
768 ! ww(i,k,j)=ww(i,k-1,j) &
769 ! - dnw(k-1)* ( dmdt(i) &
770 ! +rdx*(ru(i+1,k-1,j)-ru(i,k-1,j)) &
771 ! +rdy*(rv(i,k-1,j+1)-rv(i,k-1,j)) )
773 ww(i,k,j)=ww(i,k-1,j) - dnw(k-1)*dmdt(i) - divv(i,k-1)
775 ENDDO
776 ENDDO
777 ENDDO
780 END SUBROUTINE calc_ww_cp
783 !-------------------------------------------------------------------------------
785 SUBROUTINE calc_cq ( moist, cqu, cqv, cqw, n_moist, &
786 ids, ide, jds, jde, kds, kde, &
787 ims, ime, jms, jme, kms, kme, &
788 its, ite, jts, jte, kts, kte )
790 IMPLICIT NONE
792 ! Input data
794 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
795 ims, ime, jms, jme, kms, kme, &
796 its, ite, jts, jte, kts, kte
798 INTEGER , INTENT(IN ) :: n_moist
801 REAL, DIMENSION( ims:ime, kms:kme , jms:jme , n_moist ), INTENT(IN ) :: moist
803 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT( OUT) :: cqu, cqv, cqw
805 ! Local stuff
807 REAL :: qtot
809 INTEGER :: i, j, k, itf, jtf, ktf, ispe
811 !<DESCRIPTION>
812 !
813 ! calc_cq calculates moist coefficients for the momentum equations.
814 !
815 !</DESCRIPTION>
817 itf=ite
818 jtf=MIN(jte,jde-1)
819 ktf=MIN(kte,kde-1)
821 IF( n_moist >= PARAM_FIRST_SCALAR ) THEN
823 DO j=jts,jtf
824 DO k=kts,ktf
825 DO i=its,itf
826 qtot = 0.
827 !DEC$ loop count(3)
828 DO ispe=PARAM_FIRST_SCALAR,n_moist
829 qtot = qtot + moist(i,k,j,ispe) + moist(i-1,k,j,ispe)
830 ENDDO
831 ! qtot = 0.5*( moist(i ,k,j,1)+moist(i ,k,j,2)+moist(i ,k,j,3)+ &
832 ! & moist(i-1,k,j,1)+moist(i-1,k,j,2)+moist(i-1,k,j,3) )
833 ! cqu(i,k,j) = 1./(1.+qtot)
834 cqu(i,k,j) = 1./(1.+0.5*qtot)
835 ENDDO
836 ENDDO
837 ENDDO
839 itf=MIN(ite,ide-1)
840 jtf=jte
842 DO j=jts,jtf
843 DO k=kts,ktf
844 DO i=its,itf
845 qtot = 0.
846 !DEC$ loop count(3)
847 DO ispe=PARAM_FIRST_SCALAR,n_moist
848 qtot = qtot + moist(i,k,j,ispe) + moist(i,k,j-1,ispe)
849 ENDDO
850 ! qtot = 0.5*( moist(i,k,j ,1)+moist(i,k,j ,2)+moist(i,k,j ,3)+ &
851 ! & moist(i,k,j-1,1)+moist(i,k,j-1,2)+moist(i,k,j-1,3) )
852 ! cqv(i,k,j) = 1./(1.+qtot)
853 cqv(i,k,j) = 1./(1.+0.5*qtot)
854 ENDDO
855 ENDDO
856 ENDDO
858 itf=MIN(ite,ide-1)
859 jtf=MIN(jte,jde-1)
860 DO j=jts,jtf
861 DO k=kts+1,ktf
862 DO i=its,itf
863 qtot = 0.
864 !DEC$ loop count(3)
865 DO ispe=PARAM_FIRST_SCALAR,n_moist
866 qtot = qtot + moist(i,k,j,ispe) + moist(i,k-1,j,ispe)
867 ENDDO
868 ! qtot = 0.5*( moist(i,k ,j,1)+moist(i,k ,j,2)+moist(i,k-1,j,3)+ &
869 ! & moist(i,k-1,j,1)+moist(i,k-1,j,2)+moist(i,k ,j,3) )
870 ! cqw(i,k,j) = qtot
871 cqw(i,k,j) = 0.5*qtot
872 ENDDO
873 ENDDO
874 ENDDO
876 ELSE
878 DO j=jts,jtf
879 DO k=kts,ktf
880 DO i=its,itf
881 cqu(i,k,j) = 1.
882 ENDDO
883 ENDDO
884 ENDDO
886 itf=MIN(ite,ide-1)
887 jtf=jte
889 DO j=jts,jtf
890 DO k=kts,ktf
891 DO i=its,itf
892 cqv(i,k,j) = 1.
893 ENDDO
894 ENDDO
895 ENDDO
897 itf=MIN(ite,ide-1)
898 jtf=MIN(jte,jde-1)
899 DO j=jts,jtf
900 DO k=kts+1,ktf
901 DO i=its,itf
902 cqw(i,k,j) = 0.
903 ENDDO
904 ENDDO
905 ENDDO
907 END IF
909 END SUBROUTINE calc_cq
911 !----------------------------------------------------------------------
913 SUBROUTINE calc_alt ( alt, al, alb, &
914 ids, ide, jds, jde, kds, kde, &
915 ims, ime, jms, jme, kms, kme, &
916 its, ite, jts, jte, kts, kte )
918 IMPLICIT NONE
920 ! Input data
922 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
923 ims, ime, jms, jme, kms, kme, &
924 its, ite, jts, jte, kts, kte
926 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: alb, al
927 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT( OUT) :: alt
929 ! Local stuff
931 INTEGER :: i, j, k, itf, jtf, ktf
933 !<DESCRIPTION>
934 !
935 ! calc_alt computes the full inverse density
936 !
937 !</DESCRIPTION>
939 itf=MIN(ite,ide-1)
940 jtf=MIN(jte,jde-1)
941 ktf=MIN(kte,kde-1)
943 DO j=jts,jtf
944 DO k=kts,ktf
945 DO i=its,itf
946 alt(i,k,j) = al(i,k,j)+alb(i,k,j)
947 ENDDO
948 ENDDO
949 ENDDO
952 END SUBROUTINE calc_alt
954 !----------------------------------------------------------------------
956 SUBROUTINE calc_p_rho_phi ( moist, n_moist, hypsometric_opt, &
957 al, alb, mu, muts, ph, phb, p, pb, &
958 t, p0, t0, ptop, znu, znw, dnw, rdnw, &
959 rdn, non_hydrostatic, &
960 ids, ide, jds, jde, kds, kde, &
961 ims, ime, jms, jme, kms, kme, &
962 its, ite, jts, jte, kts, kte )
964 IMPLICIT NONE
966 ! Input data
968 LOGICAL , INTENT(IN ) :: non_hydrostatic
970 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
971 ims, ime, jms, jme, kms, kme, &
972 its, ite, jts, jte, kts, kte
974 INTEGER , INTENT(IN ) :: n_moist
975 INTEGER , INTENT(IN ) :: hypsometric_opt
977 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: alb, &
978 pb, &
979 t
981 REAL, DIMENSION( ims:ime , kms:kme , jms:jme, n_moist ), INTENT(IN ) :: moist
983 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT( OUT) :: al, p
985 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: ph, phb
987 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN ) :: mu, muts
989 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: znu, znw, dnw, rdnw, rdn
991 REAL, INTENT(IN ) :: t0, p0, ptop
993 ! Local stuff
995 INTEGER :: i, j, k, itf, jtf, ktf, ispe
996 REAL :: qvf, qtot, qf1, qf2
997 REAL, DIMENSION( its:ite) :: temp,cpovcv_v
998 REAL :: pfu, phm, pfd
1000 !<DESCRIPTION>
1001 !
1002 ! For the nonhydrostatic option, calc_p_rho_phi calculates the
1003 ! diagnostic quantities pressure and (inverse) density from the
1004 ! prognostic variables using the equation of state.
1005 !
1006 ! For the hydrostatic option, calc_p_rho_phi calculates the
1007 ! diagnostic quantities (inverse) density and geopotential from the
1008 ! prognostic variables using the equation of state and the hydrostatic
1009 ! equation.
1010 !
1011 !</DESCRIPTION>
1013 itf=MIN(ite,ide-1)
1014 jtf=MIN(jte,jde-1)
1015 ktf=MIN(kte,kde-1)
1017 #ifndef INTELMKL
1018 cpovcv_v = cpovcv
1019 #endif
1021 IF (non_hydrostatic) THEN
1023 IF (hypsometric_opt == 1) THEN
1024 DO j=jts,jtf
1025 DO k=kts,ktf
1026 DO i=its,itf
1027 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)))
1028 END DO
1029 END DO
1030 END DO
1031 ELSE IF (hypsometric_opt == 2) THEN
1033 ! The relation used to get specific volume, al, is: al = -dZ/dp,
1034 ! where dp = mut * d(eta). The pressure depth, dp, is replaced with
1035 ! p*(dp/p) ~ p*LOG((p+0.5dp)/(p-0.5dp)). Difference between dp and p*dLOG(p)
1036 ! is as follows: p*dLOG(p) - dp = 1/12*(dp/p)**3 + 1/90*(dp/p)**5 + ...
1037 ! Therefore, p*dLOG(p) is always larger than dp and the difference is
1038 ! in proportion to dp/p. TKW, 02/16/2010
1040 DO j=jts,jtf
1041 DO k=kts,ktf
1042 DO i=its,itf
1043 pfu = muts(i,j)*znw(k+1)+ptop
1044 pfd = muts(i,j)*znw(k) +ptop
1045 phm = muts(i,j)*znu(k) +ptop
1046 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)
1047 END DO
1048 END DO
1049 END DO
1050 ELSE
1051 CALL wrf_error_fatal ( 'calc_p_rho_phi: hypsometric_opt should be 1 or 2' )
1052 END IF
1054 IF (n_moist >= PARAM_FIRST_SCALAR ) THEN
1056 DO j=jts,jtf
1057 DO k=kts,ktf
1058 DO i=its,itf
1059 qvf = 1.+rvovrd*moist(i,k,j,P_QV)
1060 temp(i)=(r_d*(t0+t(i,k,j))*qvf)/(p0*(al(i,k,j)+alb(i,k,j)))
1061 ENDDO
1062 #ifdef INTELMKL
1063 CALL VPOWX ( itf-its+1, temp(its), cpovcv, p(its,k,j) )
1064 #else
1065 ! use vector version from libmassv or from compat lib in frame/libmassv.F
1066 CALL VPOW ( p(its,k,j), temp(its), cpovcv_v(its), itf-its+1 )
1067 #endif
1068 DO i=its,itf
1069 p(i,k,j)= p(i,k,j)*p0-pb(i,k,j)
1070 ENDDO
1071 ENDDO
1072 ENDDO
1074 ELSE
1076 DO j=jts,jtf
1077 DO k=kts,ktf
1078 DO i=its,itf
1079 p(i,k,j)=p0*( (r_d*(t0+t(i,k,j)))/ &
1080 (p0*(al(i,k,j)+alb(i,k,j))) )**cpovcv &
1081 -pb(i,k,j)
1082 ENDDO
1083 ENDDO
1084 ENDDO
1086 END IF
1088 ELSE
1090 ! hydrostatic pressure, al, and ph1 calc; WCS, 5 sept 2001
1093 IF (n_moist >= PARAM_FIRST_SCALAR ) THEN
1095 DO j=jts,jtf
1097 k=ktf ! top layer
1098 DO i=its,itf
1100 qtot = 0.
1101 DO ispe=PARAM_FIRST_SCALAR,n_moist
1102 qtot = qtot + moist(i,k,j,ispe)
1103 ENDDO
1104 qf2 = 1.
1105 qf1 = qtot*qf2
1107 p(i,k,j) = - 0.5*(mu(i,j)+qf1*muts(i,j))/rdnw(k)/qf2
1108 qvf = 1.+rvovrd*moist(i,k,j,P_QV)
1109 al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* &
1110 (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1112 ENDDO
1114 DO k=ktf-1,kts,-1 ! remaining layers, integrate down
1115 DO i=its,itf
1117 qtot = 0.
1118 DO ispe=PARAM_FIRST_SCALAR,n_moist
1119 qtot = qtot + 0.5*( moist(i,k ,j,ispe) + moist(i,k+1,j,ispe) )
1120 ENDDO
1121 qf2 = 1.
1122 qf1 = qtot*qf2
1124 p(i,k,j) = p(i,k+1,j) - (mu(i,j) + qf1*muts(i,j))/qf2/rdn(k+1)
1125 qvf = 1.+rvovrd*moist(i,k,j,P_QV)
1126 al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* &
1127 (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1128 ENDDO
1129 ENDDO
1131 ENDDO
1133 ELSE
1135 DO j=jts,jtf
1137 k=ktf ! top layer
1138 DO i=its,itf
1140 qtot = 0.
1141 qf2 = 1.
1142 qf1 = qtot*qf2
1144 p(i,k,j) = - 0.5*(mu(i,j)+qf1*muts(i,j))/rdnw(k)/qf2
1145 qvf = 1.
1146 al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* &
1147 (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1149 ENDDO
1151 DO k=ktf-1,kts,-1 ! remaining layers, integrate down
1152 DO i=its,itf
1154 qtot = 0.
1155 qf2 = 1.
1156 qf1 = qtot*qf2
1158 p(i,k,j) = p(i,k+1,j) - (mu(i,j) + qf1*muts(i,j))/qf2/rdn(k+1)
1159 qvf = 1.
1160 al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* &
1161 (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1162 ENDDO
1163 ENDDO
1165 ENDDO
1167 END IF
1169 IF (hypsometric_opt == 1) THEN
1170 DO j=jts,jtf
1171 DO k=2,ktf+1 ! integrate hydrostatic equation for geopotential
1172 DO i=its,itf
1173 ph(i,k,j) = ph(i,k-1,j) - (dnw(k-1))*( &
1174 (muts(i,j))*al(i,k-1,j)+ &
1175 mu(i,j)*alb(i,k-1,j) )
1176 ENDDO
1177 ENDDO
1178 ENDDO
1179 ELSE
1181 ! Revised hypsometric eq.: dZ=-al*p*dLOG(p), where p is dry pressure
1183 DO j=jts,jtf
1184 DO i=its,itf
1185 ph(i,kts,j) = phb(i,kts,j)
1186 END DO
1188 DO k=kts+1,ktf+1
1189 DO i=its,itf
1190 pfu = muts(i,j)*znw(k) +ptop
1191 pfd = muts(i,j)*znw(k-1)+ptop
1192 phm = muts(i,j)*znu(k-1)+ptop
1193 ph(i,k,j) = ph(i,k-1,j) + (al(i,k-1,j)+alb(i,k-1,j))*phm*LOG(pfd/pfu)
1194 ENDDO
1195 ENDDO
1197 DO k=kts,ktf+1
1198 DO i=its,itf
1199 ph(i,k,j) = ph(i,k,j) - phb(i,k,j)
1200 END DO
1201 END DO
1202 END DO
1204 END IF
1206 END IF
1208 END SUBROUTINE calc_p_rho_phi
1210 !----------------------------------------------------------------------
1212 SUBROUTINE calc_php ( php, ph, phb, &
1213 ids, ide, jds, jde, kds, kde, &
1214 ims, ime, jms, jme, kms, kme, &
1215 its, ite, jts, jte, kts, kte )
1217 IMPLICIT NONE
1219 ! Input data
1221 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1222 ims, ime, jms, jme, kms, kme, &
1223 its, ite, jts, jte, kts, kte
1225 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: phb, ph
1226 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT( OUT) :: php
1228 ! Local stuff
1230 INTEGER :: i, j, k, itf, jtf, ktf
1232 !<DESCRIPTION>
1233 !
1234 ! calc_php calculates the full geopotential from the reference state
1235 ! geopotential and the perturbation geopotential (phb_ph).
1236 !
1237 !</DESCRIPTION>
1239 itf=MIN(ite,ide-1)
1240 jtf=MIN(jte,jde-1)
1241 ktf=MIN(kte,kde-1)
1243 DO j=jts,jtf
1244 DO k=kts,ktf
1245 DO i=its,itf
1246 php(i,k,j) = 0.5*(phb(i,k,j)+phb(i,k+1,j)+ph(i,k,j)+ph(i,k+1,j))
1247 ENDDO
1248 ENDDO
1249 ENDDO
1251 END SUBROUTINE calc_php
1253 !-------------------------------------------------------------------------------
1255 SUBROUTINE diagnose_w( ph_tend, ph_new, ph_old, w, mu, dt, &
1256 u, v, ht, &
1257 cf1, cf2, cf3, rdx, rdy, &
1258 msftx, msfty, &
1259 ids, ide, jds, jde, kds, kde, &
1260 ims, ime, jms, jme, kms, kme, &
1261 its, ite, jts, jte, kts, kte )
1263 IMPLICIT NONE
1265 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1266 ims, ime, jms, jme, kms, kme, &
1267 its, ite, jts, jte, kts, kte
1269 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: ph_tend, &
1270 ph_new, &
1271 ph_old, &
1272 u, &
1273 v
1276 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT( OUT) :: w
1278 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mu, ht, msftx, msfty
1280 REAL, INTENT(IN ) :: dt, cf1, cf2, cf3, rdx, rdy
1282 INTEGER :: i, j, k, itf, jtf
1284 itf=MIN(ite,ide-1)
1285 jtf=MIN(jte,jde-1)
1287 !<DESCRIPTION>
1288 !
1289 ! diagnose_w diagnoses the vertical velocity from the geopoential equation.
1290 ! Used with the hydrostatic option.
1291 !
1292 !</DESCRIPTION>
1294 DO j = jts, jtf
1296 ! lower b.c. on w
1298 ! Notes on map scale factors:
1299 ! Chain rule: if Z=Z(X,Y) [true at the surface] then
1300 ! dZ/dt = dZ/dX * dX/dt + dZ/dY * dY/dt, U=dX/dt, V=dY/dt
1301 ! Using capitals to denote actual values
1302 ! In mapped values, u=U, v=V, z=Z, 1/dX=mx/dx, 1/dY=my/dy
1303 ! => w = dz/dt = mx u dz/dx + my v dz/dy
1304 ! [where dz/dx is just the surface height change between x
1305 ! gridpoints, and dz/dy is the change between y gridpoints]
1306 ! [NB: cf1, cf2 and cf3 do vertical weighting of u or v values
1307 ! nearest the surface]
1309 ! Previously msft multiplied by rdy and rdx terms.
1310 ! Now msfty multiplies rdy term, and msftx multiplies msftx term
1311 DO i = its, itf
1312 w(i,1,j)= msfty(i,j)*.5*rdy*( &
1313 (ht(i,j+1)-ht(i,j )) &
1314 *(cf1*v(i,1,j+1)+cf2*v(i,2,j+1)+cf3*v(i,3,j+1)) &
1315 +(ht(i,j )-ht(i,j-1)) &
1316 *(cf1*v(i,1,j )+cf2*v(i,2,j )+cf3*v(i,3,j )) ) &
1317 +msftx(i,j)*.5*rdx*( &
1318 (ht(i+1,j)-ht(i,j )) &
1319 *(cf1*u(i+1,1,j)+cf2*u(i+1,2,j)+cf3*u(i+1,3,j)) &
1320 +(ht(i,j )-ht(i-1,j)) &
1321 *(cf1*u(i ,1,j)+cf2*u(i ,2,j)+cf3*u(i ,3,j)) )
1322 ENDDO
1324 ! use geopotential equation to diagnose w
1326 ! Further notes on map scale factors
1327 ! If ph_tend contains: -mx partial d/dx(mu rho u/my)
1328 ! -mx partial d/dy(phi mu v/mx)
1329 ! -partial d/dz(phi mu w/my)
1330 ! then phi eqn is: partial d/dt(mu phi/my) = ph_tend + mu g w/my
1331 ! => w = [my/(mu*g)]*[partial d/dt(mu phi/my) - ph_tend]
1333 DO k = 2, kte
1334 DO i = its, itf
1335 w(i,k,j) = msfty(i,j)*( (ph_new(i,k,j)-ph_old(i,k,j))/dt &
1336 - ph_tend(i,k,j)/mu(i,j) )/g
1338 ENDDO
1339 ENDDO
1341 ENDDO
1343 END SUBROUTINE diagnose_w
1345 !-------------------------------------------------------------------------------
1347 SUBROUTINE rhs_ph( ph_tend, u, v, ww, &
1348 ph, ph_old, phb, w, &
1349 mut, muu, muv, &
1350 fnm, fnp, &
1351 rdnw, cfn, cfn1, rdx, rdy, &
1352 msfux, msfuy, msfvx, &
1353 msfvx_inv, msfvy, &
1354 msftx, msfty, &
1355 non_hydrostatic, &
1356 config_flags, &
1357 ids, ide, jds, jde, kds, kde, &
1358 ims, ime, jms, jme, kms, kme, &
1359 its, ite, jts, jte, kts, kte )
1360 IMPLICIT NONE
1362 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
1364 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1365 ims, ime, jms, jme, kms, kme, &
1366 its, ite, jts, jte, kts, kte
1368 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: &
1369 u, &
1370 v, &
1371 ww, &
1372 ph, &
1373 ph_old, &
1374 phb, &
1375 w
1377 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: ph_tend
1379 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: muu, muv, mut, &
1380 msfux, msfuy, &
1381 msfvx, msfvy, &
1382 msftx, msfty, &
1383 msfvx_inv
1385 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, fnm, fnp
1387 REAL, INTENT(IN ) :: cfn, cfn1, rdx, rdy
1389 LOGICAL, INTENT(IN ) :: non_hydrostatic
1391 ! Local stuff
1393 INTEGER :: i, j, k, itf, jtf, ktf, kz, i_start, j_start
1394 REAL :: ur, ul, ub, vr, vl, vb
1395 REAL, DIMENSION(its:ite,kts:kte) :: wdwn
1397 INTEGER :: advective_order
1399 LOGICAL :: specified
1401 !<DESCRIPTION>
1402 !
1403 ! rhs_ph calculates the large-timestep tendency terms for the geopotential
1404 ! equation. These terms include the advection and "gw". The geopotential
1405 ! equation is cast in advective form, so we don't use the flux form advection
1406 ! algorithms here.
1407 !
1408 !</DESCRIPTION>
1410 specified = .false.
1411 if(config_flags%specified .or. config_flags%nested) specified = .true.
1413 advective_order = config_flags%h_sca_adv_order
1415 itf=MIN(ite,ide-1)
1416 jtf=MIN(jte,jde-1)
1417 ktf=MIN(kte,kde-1)
1419 ! Notes on map scale factors (WCS, 2 march 2008)
1420 ! phi equation is: mu/my d/dt(phi) = -(1/my) mx mu u d/dx(phi)
1421 ! -(1/my) my mu v d/dy(phi)
1422 ! - omega d/d_eta(phi)
1423 ! +mu g w/my
1424 !
1425 ! A little further explanation...
1426 ! The tendency term we are computing here is for mu/my d/dt(phi). It is advective form
1427 ! but it is multiplied be mu/my. It will be decoupled from (mu/my) when the implicit w-phi
1428 ! solution is computed in subourine advance_w. The formulation dates from the early
1429 ! days of the mass coordinate model when we were testing both a flux and an advective formulation
1430 ! for the geopotential equation and different forms of the vertical momentum equation and the
1431 ! vertically implicit solver.
1433 ! advective form for the geopotential equation
1435 DO j = jts, jtf
1437 DO k = 2, kte
1438 DO i = its, itf
1439 wdwn(i,k) = .5*(ww(i,k,j)+ww(i,k-1,j))*rdnw(k-1) &
1440 *(ph(i,k,j)-ph(i,k-1,j)+phb(i,k,j)-phb(i,k-1,j))
1441 ENDDO
1442 ENDDO
1444 ! RHS term 3 is: - omega partial d/dnu(phi)
1446 DO k = 2, kte-1
1447 DO i = its, itf
1448 ph_tend(i,k,j) = ph_tend(i,k,j) &
1449 - (fnm(k)*wdwn(i,k+1)+fnp(k)*wdwn(i,k))
1450 ENDDO
1451 ENDDO
1453 ENDDO
1455 IF (non_hydrostatic) THEN ! add in "gw" term.
1456 DO j = jts, jtf ! in hydrostatic mode, "gw" will be diagnosed
1457 ! after the timestep to give us "w"
1458 DO i = its, itf
1459 ph_tend(i,kde,j) = 0.
1460 ENDDO
1462 DO k = 2, kte
1463 DO i = its, itf
1464 ! phi equation RHS term 4
1465 ph_tend(i,k,j) = ph_tend(i,k,j) + mut(i,j)*g*w(i,k,j)/msfty(i,j)
1466 ENDDO
1467 ENDDO
1469 ENDDO
1471 END IF
1473 ! Notes on map scale factors:
1474 ! RHS terms 1 and 2 are: -(1/my) mx u mu partial d/dx(phi)
1475 ! -(1/my) my v mu partial d/dy(phi)
1477 IF (advective_order <= 2) THEN
1479 ! y (v) advection
1481 i_start = its
1482 j_start = jts
1483 itf=MIN(ite,ide-1)
1484 jtf=MIN(jte,jde-1)
1486 IF ( (config_flags%open_ys .or. specified) .and. jts == jds ) j_start = jts+1
1487 IF ( (config_flags%open_ye .or. specified) .and. jte == jde ) jtf = jtf-2
1489 DO j = j_start, jtf
1491 DO k = 2, kte-1
1492 DO i = i_start, itf
1493 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* &
1494 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)* &
1495 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1496 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j )* &
1497 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1498 ENDDO
1499 ENDDO
1501 k = kte
1502 DO i = i_start, itf
1503 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* &
1504 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)* &
1505 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1506 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j )* &
1507 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1508 ENDDO
1510 ENDDO
1512 ! x (u) advection
1514 i_start = its
1515 j_start = jts
1516 itf=MIN(ite,ide-1)
1517 jtf=MIN(jte,jde-1)
1519 IF ( (config_flags%open_xs .or. specified) .and. its == ids ) i_start = its+1
1520 IF ( (config_flags%open_xe .or. specified) .and. ite == ide ) itf = itf-2
1522 DO j = j_start, jtf
1524 DO k = 2, kte-1
1525 DO i = i_start, itf
1526 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))* &
1527 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)* &
1528 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1529 +muu(i ,j)*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j)* &
1530 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1531 ENDDO
1532 ENDDO
1534 k = kte
1535 DO i = i_start, itf
1536 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))* &
1537 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)* &
1538 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1539 +muu(i ,j)*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux( i,j)* &
1540 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1541 ENDDO
1543 ENDDO
1545 ELSE IF (advective_order <= 4) THEN
1547 ! y (v) advection
1549 i_start = its
1550 j_start = jts
1551 itf=MIN(ite,ide-1)
1552 jtf=MIN(jte,jde-1)
1554 IF ( (config_flags%open_ys .or. specified) .and. jts == jds ) j_start = jts+2
1555 IF ( (config_flags%open_ye .or. specified) .and. jte == jde ) jtf = jtf-3
1557 DO j = j_start, jtf
1559 DO k = 2, kte-1
1560 DO i = i_start, itf
1561 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))*( &
1562 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1) &
1563 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j ))* (1./12.)*( &
1564 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1565 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1566 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1567 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1570 ENDDO
1571 ENDDO
1573 k = kte
1574 DO i = i_start, itf
1575 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))*( &
1576 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1) &
1577 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j ))* (1./12.)*( &
1578 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1579 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1580 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1581 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1583 ENDDO
1585 ENDDO
1587 ! pick up near boundary rows using 2nd order stencil for open and specified conditions
1589 IF ( (config_flags%open_ys .or. specified) .and. jts <= jds+1 ) THEN
1591 j = jds+1
1592 DO k = 2, kte-1
1593 DO i = i_start, itf
1594 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* &
1595 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)* &
1596 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1597 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j )* &
1598 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1599 ENDDO
1600 ENDDO
1602 k = kte
1603 DO i = i_start, itf
1604 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* &
1605 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)* &
1606 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1607 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j )* &
1608 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1609 ENDDO
1611 END IF
1613 IF ( (config_flags%open_ye .or. specified) .and. jte >= jde-2 ) THEN
1615 j = jde-2
1616 DO k = 2, kte-1
1617 DO i = i_start, itf
1618 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* &
1619 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)* &
1620 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1621 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j )* &
1622 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1623 ENDDO
1624 ENDDO
1626 k = kte
1627 DO i = i_start, itf
1628 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* &
1629 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)* &
1630 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1631 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j )* &
1632 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1633 ENDDO
1635 END IF
1637 ! x (u) advection
1639 i_start = its
1640 j_start = jts
1641 itf=MIN(ite,ide-1)
1642 jtf=MIN(jte,jde-1)
1644 IF ( (config_flags%open_xs) .and. its == ids ) i_start = its+2
1645 IF ( (config_flags%open_xe) .and. ite == ide ) itf = itf-3
1647 DO j = j_start, jtf
1649 DO k = 2, kte-1
1650 DO i = i_start, itf
1651 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*( &
1652 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j) &
1653 +muu(i,j )*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j) )* (1./12.)*( &
1654 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1655 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1656 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1657 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1658 ENDDO
1659 ENDDO
1661 k = kte
1662 DO i = i_start, itf
1663 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*( &
1664 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j) &
1665 +muu(i,j )*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux(i ,j) )* (1./12.)*( &
1666 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1667 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1668 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1669 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1670 ENDDO
1672 ENDDO
1674 ! pick up near boundary rows using 2nd order stencil for open and specified conditions
1676 IF ( (config_flags%open_xs .or. specified) .and. its <= ids+1 ) THEN
1678 i = ids + 1
1680 DO j = j_start, jtf
1681 DO k = 2, kte-1
1682 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))* &
1683 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)* &
1684 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1685 +muu(i ,j)*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j)* &
1686 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1687 ENDDO
1688 ENDDO
1690 k = kte
1691 DO j = j_start, jtf
1692 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))* &
1693 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)* &
1694 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1695 +muu(i ,j)*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux( i,j)* &
1696 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1697 ENDDO
1699 END IF
1701 IF ( (config_flags%open_xe .or. specified) .and. ite >= ide-2 ) THEN
1703 i = ide-2
1704 DO j = j_start, jtf
1705 DO k = 2, kte-1
1706 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))* &
1707 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)* &
1708 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1709 +muu(i ,j)*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j)* &
1710 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1711 ENDDO
1712 ENDDO
1714 k = kte
1715 DO j = j_start, jtf
1716 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))* &
1717 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)* &
1718 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1719 +muu(i ,j)*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux( i,j)* &
1720 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1721 ENDDO
1723 END IF
1725 !--------------------------
1727 ELSE IF (advective_order <= 6) THEN
1729 !!! NOTE: this branch has been looked at and fixed with changes for overdecomposition
1730 !!! the branches covering the other advective_order cases have not. 20090923. JM
1732 ! y (v) advection
1734 i_start = its
1735 j_start = jts
1736 itf=MIN(ite,ide-1)
1737 jtf=MIN(jte,jde-1)
1739 IF (config_flags%open_ys .or. specified ) j_start = max(jts,jds+3)
1740 IF (config_flags%open_ye .or. specified ) jtf = min(jtf,jde-4)
1742 DO j = j_start, jtf
1744 DO k = 2, kte-1
1745 DO i = i_start, itf
1746 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* ( &
1747 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1) &
1748 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j ) )* (1./60.)*( &
1749 45.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1750 -9.*(ph(i,k,j+2)-ph(i,k,j-2)) &
1751 +(ph(i,k,j+3)-ph(i,k,j-3)) &
1752 +45.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1753 -9.*(phb(i,k,j+2)-phb(i,k,j-2)) &
1754 +(phb(i,k,j+3)-phb(i,k,j-3)) ) )
1757 ENDDO
1758 ENDDO
1760 k = kte
1761 DO i = i_start, itf
1762 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* ( &
1763 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1) &
1764 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j ) )* (1./60.)*( &
1765 45.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1766 -9.*(ph(i,k,j+2)-ph(i,k,j-2)) &
1767 +(ph(i,k,j+3)-ph(i,k,j-3)) &
1768 +45.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1769 -9.*(phb(i,k,j+2)-phb(i,k,j-2)) &
1770 +(phb(i,k,j+3)-phb(i,k,j-3)) ) )
1772 ENDDO
1774 ENDDO
1776 ! 4th order stencil for open or specified conditions two in form the boundary
1778 IF ( (config_flags%open_ys .or. specified) .and. jts <= jds+2 .and. jds+2 <= jte ) THEN
1780 j = jds+2
1781 DO k = 2, kte-1
1782 DO i = i_start, itf
1783 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* ( &
1784 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1) &
1785 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j ) )* (1./12.)*( &
1786 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1787 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1788 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1789 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1792 ENDDO
1793 ENDDO
1795 k = kte
1796 DO i = i_start, itf
1797 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* ( &
1798 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1) &
1799 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j) )* (1./12.)*( &
1800 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1801 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1802 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1803 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1805 ENDDO
1807 END IF
1809 IF ( (config_flags%open_ye .or. specified) .and. jts <= jde-3 .and. jde-3 <= jte ) THEN
1810 j = jde-3
1811 DO k = 2, kte-1
1812 DO i = i_start, itf
1813 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* ( &
1814 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1) &
1815 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j) )* (1./12.)*( &
1816 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1817 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1818 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1819 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1822 ENDDO
1823 ENDDO
1825 k = kte
1826 DO i = i_start, itf
1827 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* ( &
1828 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1) &
1829 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j) )* (1./12.)*( &
1830 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1831 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1832 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1833 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1835 ENDDO
1837 END IF
1839 ! 2nd order stencil for open and specified conditions one row in from boundary
1841 IF ( (config_flags%open_ys .or. specified) .and. jts <= jds+1 .and. jds+1 <= jte ) THEN
1843 j = jds+1
1844 DO k = 2, kte-1
1845 DO i = i_start, itf
1846 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* &
1847 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)* &
1848 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1849 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j )* &
1850 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1851 ENDDO
1852 ENDDO
1854 k = kte
1855 DO i = i_start, itf
1856 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* &
1857 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)* &
1858 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1859 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j )* &
1860 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1861 ENDDO
1863 END IF
1865 IF ( (config_flags%open_ye .or. specified) .and. jts <= jde-2 .and. jde-2 <= jte ) THEN
1867 j = jde-2
1868 DO k = 2, kte-1
1869 DO i = i_start, itf
1870 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* &
1871 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)* &
1872 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1873 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j )* &
1874 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1875 ENDDO
1876 ENDDO
1878 k = kte
1879 DO i = i_start, itf
1880 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* &
1881 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)* &
1882 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1883 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j )* &
1884 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1885 ENDDO
1887 END IF
1889 ! x (u) advection
1891 i_start = its
1892 j_start = jts
1893 itf=MIN(ite,ide-1)
1894 jtf=MIN(jte,jde-1)
1896 IF (config_flags%open_xs .or. specified ) i_start = max(its,ids+3)
1897 IF (config_flags%open_xe .or. specified ) itf = min(itf,ide-4)
1899 DO j = j_start, jtf
1901 DO k = 2, kte-1
1902 DO i = i_start, itf
1903 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*( &
1904 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j) &
1905 +muu(i,j )*(u(i,k,j )+u(i,k-1,j ))*msfux(i,j) )* (1./60.)*( &
1906 45.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1907 -9.*(ph(i+2,k,j)-ph(i-2,k,j)) &
1908 +(ph(i+3,k,j)-ph(i-3,k,j)) &
1909 +45.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1910 -9.*(phb(i+2,k,j)-phb(i-2,k,j)) &
1911 +(phb(i+3,k,j)-phb(i-3,k,j)) ) )
1912 ENDDO
1913 ENDDO
1915 k = kte
1916 DO i = i_start, itf
1917 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*( &
1918 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j) &
1919 +muu(i,j )*(cfn*u(i ,k-1,j)+cfn1*u(i,k-2,j))*msfux(i,j) )* (1./60.)*( &
1920 45.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1921 -9.*(ph(i+2,k,j)-ph(i-2,k,j)) &
1922 +(ph(i+3,k,j)-ph(i-3,k,j)) &
1923 +45.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1924 -9.*(phb(i+2,k,j)-phb(i-2,k,j)) &
1925 +(phb(i+3,k,j)-phb(i-3,k,j)) ) )
1926 ENDDO
1928 ENDDO
1930 ! 4th order stencil two in from the boundary for open and specified conditions
1932 IF ( (config_flags%open_xs) .and. its <= ids+2 .and. ids+2 <= ite ) THEN
1933 i = ids + 2
1934 DO j = j_start, jtf
1935 DO k = 2, kte-1
1936 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*( &
1937 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j) &
1938 +muu(i,j )*(u(i,k,j )+u(i,k-1,j ))*msfux(i,j) )* (1./12.)*( &
1939 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1940 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1941 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1942 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1943 ENDDO
1944 k = kte
1945 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*( &
1946 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j) &
1947 +muu(i,j )*(cfn*u(i ,k-1,j)+cfn1*u(i,k-2,j))*msfux(i,j) )* (1./12.)*( &
1948 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1949 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1950 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1951 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1953 ENDDO
1954 END IF
1956 IF ( (config_flags%open_xe) .and. its <= ide-3 .and. ide-3 <= ite ) THEN
1957 i = ide-3
1958 DO j = j_start, jtf
1959 DO k = 2, kte-1
1960 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*( &
1961 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j) &
1962 +muu(i,j )*(u(i,k,j )+u(i,k-1,j ))*msfux(i,j) )* (1./12.)*( &
1963 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1964 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1965 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1966 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1967 ENDDO
1968 k = kte
1969 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*( &
1970 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j) &
1971 +muu(i,j )*(cfn*u(i ,k-1,j)+cfn1*u(i,k-2,j))*msfux(i,j) )* (1./12.)*( &
1972 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1973 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1974 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1975 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1977 ENDDO
1978 END IF
1980 ! 2nd order stencil for open and specified conditions one in from the boundary
1982 IF ( (config_flags%open_xs .or. specified) .and. its <= ids+1 .and. ids+1 <= ite ) THEN
1984 i = ids + 1
1986 DO j = j_start, jtf
1987 DO k = 2, kte-1
1988 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))* &
1989 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)* &
1990 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1991 +muu(i ,j)*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j)* &
1992 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1993 ENDDO
1994 ENDDO
1996 k = kte
1997 DO j = j_start, jtf
1998 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))* &
1999 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)* &
2000 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
2001 +muu(i ,j)*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux( i,j)* &
2002 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
2003 ENDDO
2005 END IF
2007 IF ( (config_flags%open_xe .or. specified) .and. its <= ide-2 .and. ide-2 <= ite ) THEN
2009 i = ide-2
2010 DO j = j_start, jtf
2011 DO k = 2, kte-1
2012 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))* &
2013 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)* &
2014 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
2015 +muu(i ,j)*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j)* &
2016 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
2017 ENDDO
2018 ENDDO
2020 k = kte
2021 DO j = j_start, jtf
2022 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))* &
2023 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)* &
2024 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
2025 +muu(i ,j)*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux( i,j)* &
2026 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
2027 ENDDO
2029 END IF
2031 END IF ! 6th order advection
2033 ! lateral open boundary conditions,
2034 ! start with north and south (y) boundaries
2036 i_start = its
2037 itf=MIN(ite,ide-1)
2039 ! south
2041 IF ( (config_flags%open_ys) .and. jts == jds ) THEN
2043 j=jts
2045 DO k=2,kde
2046 kz = min(k,kde-1)
2047 DO i = its,itf
2048 vb =.5*( fnm(kz)*(v(i,kz ,j+1)+v(i,kz ,j )) &
2049 +fnp(kz)*(v(i,kz-1,j+1)+v(i,kz-1,j )) )
2050 vl=amin1(vb,0.)
2051 ph_tend(i,k,j)=ph_tend(i,k,j)-rdy*mut(i,j)*( &
2052 +vl*(ph_old(i,k,j+1)-ph_old(i,k,j)))
2053 ENDDO
2054 ENDDO
2056 END IF
2058 ! north
2060 IF ( (config_flags%open_ye) .and. jte == jde ) THEN
2062 j=jte-1
2064 DO k=2,kde
2065 kz = min(k,kde-1)
2066 DO i = its,itf
2067 vb=.5*( fnm(kz)*(v(i,kz ,j+1)+v(i,kz ,j)) &
2068 +fnp(kz)*(v(i,kz-1,j+1)+v(i,kz-1,j)) )
2069 vr=amax1(vb,0.)
2070 ph_tend(i,k,j)=ph_tend(i,k,j)-rdy*mut(i,j)*( &
2071 +vr*(ph_old(i,k,j)-ph_old(i,k,j-1)))
2072 ENDDO
2073 ENDDO
2075 END IF
2077 ! now the east and west (y) boundaries
2079 j_start = its
2080 jtf=MIN(jte,jde-1)
2082 ! west
2084 IF ( (config_flags%open_xs) .and. its == ids ) THEN
2086 i=its
2088 DO j = jts,jtf
2089 DO k=2,kde-1
2090 kz = k
2091 ub =.5*( fnm(kz)*(u(i+1,kz ,j)+u(i ,kz ,j)) &
2092 +fnp(kz)*(u(i+1,kz-1,j)+u(i ,kz-1,j)) )
2093 ul=amin1(ub,0.)
2094 ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*( &
2095 +ul*(ph_old(i+1,k,j)-ph_old(i,k,j)))
2096 ENDDO
2098 k = kde
2099 kz = k
2100 ub =.5*( fnm(kz)*(u(i+1,kz ,j)+u(i ,kz ,j)) &
2101 +fnp(kz)*(u(i+1,kz-1,j)+u(i ,kz-1,j)) )
2102 ul=amin1(ub,0.)
2103 ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*( &
2104 +ul*(ph_old(i+1,k,j)-ph_old(i,k,j)))
2105 ENDDO
2107 END IF
2109 ! east
2111 IF ( (config_flags%open_xe) .and. ite == ide ) THEN
2113 i = ite-1
2115 DO j = jts,jtf
2116 DO k=2,kde-1
2117 kz = k
2118 ub=.5*( fnm(kz)*(u(i+1,kz ,j)+u(i,kz ,j)) &
2119 +fnp(kz)*(u(i+1,kz-1,j)+u(i,kz-1,j)) )
2120 ur=amax1(ub,0.)
2121 ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*( &
2122 +ur*(ph_old(i,k,j)-ph_old(i-1,k,j)))
2123 ENDDO
2125 k = kde
2126 kz = k-1
2127 ub=.5*( fnm(kz)*(u(i+1,kz ,j)+u(i,kz ,j)) &
2128 +fnp(kz)*(u(i+1,kz-1,j)+u(i,kz-1,j)) )
2129 ur=amax1(ub,0.)
2130 ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*( &
2131 +ur*(ph_old(i,k,j)-ph_old(i-1,k,j)))
2133 ENDDO
2135 END IF
2137 END SUBROUTINE rhs_ph
2140 !-------------------------------------------------------------------------------
2142 SUBROUTINE horizontal_pressure_gradient( ru_tend,rv_tend, &
2143 ph,alt,p,pb,al,php,cqu,cqv, &
2144 muu,muv,mu,fnm,fnp,rdnw, &
2145 cf1,cf2,cf3,rdx,rdy,msfux,msfuy,&
2146 msfvx,msfvy,msftx,msfty, &
2147 config_flags, non_hydrostatic, &
2148 top_lid, &
2149 ids, ide, jds, jde, kds, kde, &
2150 ims, ime, jms, jme, kms, kme, &
2151 its, ite, jts, jte, kts, kte )
2153 IMPLICIT NONE
2155 ! Input data
2158 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2160 LOGICAL, INTENT (IN ) :: non_hydrostatic, top_lid
2162 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2163 ims, ime, jms, jme, kms, kme, &
2164 its, ite, jts, jte, kts, kte
2166 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: &
2167 ph, &
2168 alt, &
2169 al, &
2170 p, &
2171 pb, &
2172 php, &
2173 cqu, &
2174 cqv
2177 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: &
2178 ru_tend, &
2179 rv_tend
2181 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: muu, muv, mu, &
2182 msfux, msfuy, &
2183 msfvx, msfvy, &
2184 msftx, msfty
2186 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, fnm, fnp
2188 REAL, INTENT(IN ) :: rdx, rdy, cf1, cf2, cf3
2190 INTEGER :: i,j,k, itf, jtf, ktf, i_start, j_start
2191 REAL, DIMENSION( ims:ime, kms:kme ) :: dpn
2192 REAL :: dpx, dpy
2194 LOGICAL :: specified
2196 !<DESCRIPTION>
2197 !
2198 ! horizontal_pressure_gradient calculates the
2199 ! horizontal pressure gradient terms for the large-timestep tendency
2200 ! in the horizontal momentum equations (u,v).
2201 !
2202 !</DESCRIPTION>
2204 specified = .false.
2205 if(config_flags%specified .or. config_flags%nested) specified = .true.
2207 ! Notes on map scale factors:
2208 ! Calculates the pressure gradient terms in ADT eqns 44 and 45
2209 ! With upper rho -> 'mu', these are:
2210 ! Eqn 30: -mu*(mx/my)*(1/rho)*partial dp/dx
2211 ! Eqn 31: -mu*(my/mx)*(1/rho)*partial dp/dy
2212 !
2213 ! As we are on nu, rather than height, surfaces:
2214 !
2215 ! mu dp/dx = mu alpha partial dp'/dx + (nu mu partial dmubar/dx) alpha'
2216 ! + mu partial dphi'/dx + (partial dphi/dx)*(partial dp'/dnu - mu')
2217 !
2218 ! mu dp/dy = mu alpha partial dp'/dy + (nu mu partial dmubar/dy) alpha'
2219 ! + mu partial dphi'/dy + (partial dphi/dy)*(partial dp'/dnu - mu')
2221 ! start with the north-south (y) pressure gradient
2223 itf=MIN(ite,ide-1)
2224 jtf=jte
2225 ktf=MIN(kte,kde-1)
2226 i_start = its
2227 j_start = jts
2228 IF ( (config_flags%open_ys .or. specified .or. &
2229 config_flags%nested .or. config_flags%polar ) .and. jts == jds ) j_start = jts+1
2230 IF ( (config_flags%open_ye .or. specified .or. &
2231 config_flags%nested .or. config_flags%polar ) .and. jte == jde ) jtf = jtf-1
2233 DO j = j_start, jtf
2235 IF ( non_hydrostatic ) THEN
2237 k=1
2239 DO i = i_start, itf
2240 dpn(i,k) = .5*( cf1*(p(i,k ,j-1)+p(i,k ,j)) &
2241 +cf2*(p(i,k+1,j-1)+p(i,k+1,j)) &
2242 +cf3*(p(i,k+2,j-1)+p(i,k+2,j)) )
2243 dpn(i,kde) = 0.
2244 ENDDO
2245 IF (top_lid) THEN
2246 DO i = i_start, itf
2247 dpn(i,kde) = .5*( cf1*(p(i,kde-1,j-1)+p(i,kde-1,j)) &
2248 +cf2*(p(i,kde-2,j-1)+p(i,kde-2,j)) &
2249 +cf3*(p(i,kde-3,j-1)+p(i,kde-3,j)) )
2250 ENDDO
2251 ENDIF
2253 DO k=2,ktf
2254 DO i = i_start, itf
2255 dpn(i,k) = .5*( fnm(k)*(p(i,k ,j-1)+p(i,k ,j)) &
2256 +fnp(k)*(p(i,k-1,j-1)+p(i,k-1,j)) )
2257 END DO
2258 END DO
2260 ! ADT eqn 45: -mu*(my/mx)*(1/rho)*partial dp/dy
2261 ! [alt, al are 1/rho terms; muv, mu are NOT coupled]
2262 DO K=1,ktf
2263 DO i = i_start, itf
2264 ! Here are mu dp/dy terms 1-3
2265 dpy = (msfvy(i,j)/msfvx(i,j))*.5*rdy*muv(i,j)*( &
2266 (ph (i,k+1,j)-ph (i,k+1,j-1) + ph(i,k,j)-ph(i,k,j-1)) &
2267 +(alt(i,k ,j)+alt(i,k ,j-1))*(p (i,k,j)-p (i,k,j-1)) &
2268 +(al (i,k ,j)+al (i,k ,j-1))*(pb(i,k,j)-pb(i,k,j-1)) )
2269 ! Here is mu dp/dy term 4
2270 dpy = dpy + (msfvy(i,j)/msfvx(i,j))*rdy*(php(i,k,j)-php(i,k,j-1))* &
2271 (rdnw(k)*(dpn(i,k+1)-dpn(i,k))-.5*(mu(i,j-1)+mu(i,j)))
2272 rv_tend(i,k,j) = rv_tend(i,k,j)-cqv(i,k,j)*dpy
2273 END DO
2274 END DO
2276 ELSE
2278 ! ADT eqn 45: -mu*(my/mx)*(1/rho)*partial dp/dy
2279 ! [alt, al are 1/rho terms; muv, mu are NOT coupled]
2280 DO K=1,ktf
2281 DO i = i_start, itf
2282 ! Here are mu dp/dy terms 1-3; term 4 not needed if hydrostatic
2283 dpy = (msfvy(i,j)/msfvx(i,j))*.5*rdy*muv(i,j)*( &
2284 (ph (i,k+1,j)-ph (i,k+1,j-1) + ph(i,k,j)-ph(i,k,j-1)) &
2285 +(alt(i,k ,j)+alt(i,k ,j-1))*(p (i,k,j)-p (i,k,j-1)) &
2286 +(al (i,k ,j)+al (i,k ,j-1))*(pb(i,k,j)-pb(i,k,j-1)) )
2287 rv_tend(i,k,j) = rv_tend(i,k,j)-cqv(i,k,j)*dpy
2288 END DO
2289 END DO
2291 END IF
2293 ENDDO
2295 ! now the east-west (x) pressure gradient
2297 itf=ite
2298 jtf=MIN(jte,jde-1)
2299 ktf=MIN(kte,kde-1)
2300 i_start = its
2301 j_start = jts
2302 IF ( (config_flags%open_xs .or. specified .or. &
2303 config_flags%nested ) .and. its == ids ) i_start = its+1
2304 IF ( (config_flags%open_xe .or. specified .or. &
2305 config_flags%nested ) .and. ite == ide ) itf = itf-1
2306 IF ( config_flags%periodic_x ) i_start = its
2307 IF ( config_flags%periodic_x ) itf=ite
2309 DO j = j_start, jtf
2311 IF ( non_hydrostatic ) THEN
2313 k=1
2315 DO i = i_start, itf
2316 dpn(i,k) = .5*( cf1*(p(i-1,k ,j)+p(i,k ,j)) &
2317 +cf2*(p(i-1,k+1,j)+p(i,k+1,j)) &
2318 +cf3*(p(i-1,k+2,j)+p(i,k+2,j)) )
2319 dpn(i,kde) = 0.
2320 ENDDO
2321 IF (top_lid) THEN
2322 DO i = i_start, itf
2323 dpn(i,kde) = .5*( cf1*(p(i-1,kde-1,j)+p(i,kde-1,j)) &
2324 +cf2*(p(i-1,kde-2,j)+p(i,kde-2,j)) &
2325 +cf3*(p(i-1,kde-3,j)+p(i,kde-3,j)) )
2326 ENDDO
2327 ENDIF
2329 DO k=2,ktf
2330 DO i = i_start, itf
2331 dpn(i,k) = .5*( fnm(k)*(p(i-1,k ,j)+p(i,k ,j)) &
2332 +fnp(k)*(p(i-1,k-1,j)+p(i,k-1,j)) )
2333 END DO
2334 END DO
2336 ! ADT eqn 44: -mu*(mx/my)*(1/rho)*partial dp/dx
2337 ! [alt, al are 1/rho terms; muu, mu are NOT coupled]
2338 DO K=1,ktf
2339 DO i = i_start, itf
2340 ! Here are mu dp/dy terms 1-3
2341 dpx = (msfux(i,j)/msfuy(i,j))*.5*rdx*muu(i,j)*( &
2342 (ph (i,k+1,j)-ph (i-1,k+1,j) + ph(i,k,j)-ph(i-1,k,j)) &
2343 +(alt(i,k ,j)+alt(i-1,k ,j))*(p (i,k,j)-p (i-1,k,j)) &
2344 +(al (i,k ,j)+al (i-1,k ,j))*(pb(i,k,j)-pb(i-1,k,j)) )
2345 ! Here is mu dp/dy term 4
2346 dpx = dpx + (msfux(i,j)/msfuy(i,j))*rdx*(php(i,k,j)-php(i-1,k,j))* &
2347 (rdnw(k)*(dpn(i,k+1)-dpn(i,k))-.5*(mu(i-1,j)+mu(i,j)))
2348 ru_tend(i,k,j) = ru_tend(i,k,j)-cqu(i,k,j)*dpx
2349 END DO
2350 END DO
2352 ELSE
2354 ! ADT eqn 44: -mu*(mx/my)*(1/rho)*partial dp/dx
2355 ! [alt, al are 1/rho terms; muu, mu are NOT coupled]
2356 DO K=1,ktf
2357 DO i = i_start, itf
2358 ! Here are mu dp/dy terms 1-3; term 4 not needed if hydrostatic
2359 dpx = (msfux(i,j)/msfuy(i,j))*.5*rdx*muu(i,j)*( &
2360 (ph (i,k+1,j)-ph (i-1,k+1,j) + ph(i,k,j)-ph(i-1,k,j)) &
2361 +(alt(i,k ,j)+alt(i-1,k ,j))*(p (i,k,j)-p (i-1,k,j)) &
2362 +(al (i,k ,j)+al (i-1,k ,j))*(pb(i,k,j)-pb(i-1,k,j)) )
2363 ru_tend(i,k,j) = ru_tend(i,k,j)-cqu(i,k,j)*dpx
2364 END DO
2365 END DO
2367 END IF
2369 ENDDO
2371 END SUBROUTINE horizontal_pressure_gradient
2373 !-------------------------------------------------------------------------------
2375 SUBROUTINE pg_buoy_w( rw_tend, p, cqw, mu, mub, &
2376 rdnw, rdn, g, msftx, msfty, &
2377 ids, ide, jds, jde, kds, kde, &
2378 ims, ime, jms, jme, kms, kme, &
2379 its, ite, jts, jte, kts, kte )
2381 IMPLICIT NONE
2383 ! Input data
2385 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2386 ims, ime, jms, jme, kms, kme, &
2387 its, ite, jts, jte, kts, kte
2389 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: p
2390 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: cqw
2393 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: rw_tend
2395 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mub, mu, msftx, msfty
2397 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, rdn
2399 REAL, INTENT(IN ) :: g
2401 INTEGER :: itf, jtf, i, j, k
2402 REAL :: cq1, cq2
2405 !<DESCRIPTION>
2406 !
2407 ! pg_buoy_w calculates the
2408 ! vertical pressure gradient and buoyancy terms for the large-timestep
2409 ! tendency in the vertical momentum equation.
2410 !
2411 !</DESCRIPTION>
2413 ! BUOYANCY AND PRESSURE GRADIENT TERM IN W EQUATION AT TIME T
2415 ! Map scale factor notes
2416 ! ADT eqn 46 RHS terms 6 and 7 (where 7 is "-rho g")
2417 ! Dividing by my, and using mu and nu (see Klemp et al. eqns 32, 40)
2418 ! term 6: +(g/my) partial dp'/dnu
2419 ! term 7: -(g/my) mu'
2420 !
2421 ! For moisture-free atmosphere, cq1=1, cq2=0
2422 ! => (1./msft(i,j)) * g * [rdn(k)*{p(i,k,j)-p(i,k-1,j)}-mu(i,j)]
2424 itf=MIN(ite,ide-1)
2425 jtf=MIN(jte,jde-1)
2427 DO j = jts,jtf
2429 k=kde
2430 DO i=its,itf
2431 cq1 = 1./(1.+cqw(i,k-1,j))
2432 cq2 = cqw(i,k-1,j)*cq1
2433 rw_tend(i,k,j) = rw_tend(i,k,j)+(1./msfty(i,j))*g*( &
2434 cq1*2.*rdnw(k-1)*( -p(i,k-1,j)) &
2435 -mu(i,j)-cq2*mub(i,j) )
2436 END DO
2438 DO k = 2, kde-1
2439 DO i = its,itf
2440 cq1 = 1./(1.+cqw(i,k,j))
2441 cq2 = cqw(i,k,j)*cq1
2442 cqw(i,k,j) = cq1
2443 rw_tend(i,k,j) = rw_tend(i,k,j)+(1./msfty(i,j))*g*( &
2444 cq1*rdn(k)*(p(i,k,j)-p(i,k-1,j)) &
2445 -mu(i,j)-cq2*mub(i,j) )
2446 END DO
2447 ENDDO
2450 ENDDO
2452 END SUBROUTINE pg_buoy_w
2454 !-------------------------------------------------------------------------------
2456 SUBROUTINE w_damp( rw_tend, max_vert_cfl,max_horiz_cfl, &
2457 u, v, ww, w, mut, rdnw, &
2458 rdx, rdy, msfux, msfuy, &
2459 msfvx, msfvy, dt, &
2460 config_flags, &
2461 ids, ide, jds, jde, kds, kde, &
2462 ims, ime, jms, jme, kms, kme, &
2463 its, ite, jts, jte, kts, kte )
2465 USE module_llxy
2466 IMPLICIT NONE
2468 ! Input data
2470 TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
2472 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2473 ims, ime, jms, jme, kms, kme, &
2474 its, ite, jts, jte, kts, kte
2476 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: u, v, ww, w
2478 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: rw_tend
2480 REAL, INTENT(OUT) :: max_vert_cfl
2481 REAL, INTENT(OUT) :: max_horiz_cfl
2482 REAL :: horiz_cfl
2484 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mut
2486 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw
2488 REAL, INTENT(IN) :: dt
2489 REAL, INTENT(IN) :: rdx, rdy
2490 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, msfuy
2491 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfvx, msfvy
2493 REAL :: vert_cfl, cf_n, cf_d, maxdub, maxdeta
2495 INTEGER :: itf, jtf, i, j, k, maxi, maxj, maxk
2496 INTEGER :: some
2497 CHARACTER*512 :: temp
2499 CHARACTER (LEN=256) :: time_str
2500 CHARACTER (LEN=256) :: grid_str
2502 integer :: total
2503 REAL :: msfuxt , msfxffl
2505 !<DESCRIPTION>
2506 !
2507 ! w_damp computes a damping term for the vertical velocity when the
2508 ! vertical Courant number is too large. This was found to be preferable to
2509 ! decreasing the timestep or increasing the diffusion in real-data applications
2510 ! that produced potentially-unstable large vertical velocities because of
2511 ! unphysically large heating rates coming from the cumulus parameterization
2512 ! schemes run at moderately high resolutions (dx ~ O(10) km).
2513 !
2514 ! Additionally, w_damp returns the maximum cfl values due to vertical motion and
2515 ! horizontal motion. These values are returned via the max_vert_cfl and
2516 ! max_horiz_cfl variables. (Added by T. Hutchinson, WSI, 3/5/2007)
2517 !
2518 !</DESCRIPTION>
2520 itf=MIN(ite,ide-1)
2521 jtf=MIN(jte,jde-1)
2523 some = 0
2524 max_vert_cfl = 0.
2525 max_horiz_cfl = 0.
2526 total = 0
2528 IF(config_flags%map_proj == PROJ_CASSINI ) then
2529 msfxffl = 1.0/COS(config_flags%fft_filter_lat*degrad)
2530 END IF
2532 IF ( config_flags%w_damping == 1 ) THEN
2533 DO j = jts,jtf
2535 DO k = 2, kde-1
2536 DO i = its,itf
2537 #if 1
2538 IF(config_flags%map_proj == PROJ_CASSINI ) then
2539 msfuxt = MIN(msfux(i,j), msfxffl)
2540 ELSE
2541 msfuxt = msfux(i,j)
2542 END IF
2543 vert_cfl = abs(ww(i,k,j)/mut(i,j)*rdnw(k)*dt)
2545 IF ( vert_cfl > max_vert_cfl ) THEN
2546 max_vert_cfl = vert_cfl ; maxi = i ; maxj = j ; maxk = k
2547 maxdub = w(i,k,j) ; maxdeta = -1./rdnw(k)
2548 ENDIF
2550 horiz_cfl = max( abs(u(i,k,j) * rdx * msfuxt * dt), &
2551 abs(v(i,k,j) * rdy * msfvy(i,j) * dt) )
2552 if (horiz_cfl > max_horiz_cfl) then
2553 max_horiz_cfl = horiz_cfl
2554 endif
2556 if(vert_cfl .gt. w_beta)then
2557 #else
2558 ! restructure to get rid of divide
2559 !
2560 ! This had been used for efficiency, but with the addition of returning the cfl values,
2561 ! the old version (above) was reinstated. (T. Hutchinson, 3/5/2007)
2562 !
2563 cf_n = abs(ww(i,k,j)*rdnw(k)*dt)
2564 cf_d = abs(mut(i,j))
2565 if(cf_n .gt. cf_d*w_beta )then
2566 #endif
2568 WRITE(temp,*)i,j,k,' vert_cfl,w,d(eta)=',vert_cfl,w(i,k,j),-1./rdnw(k)
2569 CALL wrf_debug ( 100 , TRIM(temp) )
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 WRITE(temp,*)i,j,k,' vert_cfl,w,d(eta)=',vert_cfl,w(i,k,j),-1./rdnw(k)
2615 CALL wrf_debug ( 100 , TRIM(temp) )
2616 if ( vert_cfl > 2. ) some = some + 1
2617 endif
2618 END DO
2619 ENDDO
2620 ENDDO
2621 ENDIF
2622 IF ( some .GT. 0 ) THEN
2623 CALL get_current_time_string( time_str )
2624 CALL get_current_grid_name( grid_str )
2625 WRITE(temp,*)some, &
2626 ' points exceeded cfl=2 in domain '//TRIM(grid_str)//' at time '//TRIM(time_str)//' hours'
2627 CALL wrf_debug ( 0 , TRIM(temp) )
2628 WRITE(temp,*)'MAX AT i,j,k: ',maxi,maxj,maxk,' vert_cfl,w,d(eta)=',max_vert_cfl, &
2629 maxdub,maxdeta
2630 CALL wrf_debug ( 0 , TRIM(temp) )
2631 ENDIF
2633 END SUBROUTINE w_damp
2635 !-------------------------------------------------------------------------------
2637 SUBROUTINE horizontal_diffusion ( name, field, tendency, mu, &
2638 config_flags, &
2639 msfux, msfuy, msfvx, msfvx_inv, &
2640 msfvy, msftx, msfty, &
2641 khdif, xkmhd, rdx, rdy, &
2642 ids, ide, jds, jde, kds, kde, &
2643 ims, ime, jms, jme, kms, kme, &
2644 its, ite, jts, jte, kts, kte )
2646 IMPLICIT NONE
2648 ! Input data
2650 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2652 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2653 ims, ime, jms, jme, kms, kme, &
2654 its, ite, jts, jte, kts, kte
2656 CHARACTER(LEN=1) , INTENT(IN ) :: name
2658 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, xkmhd
2660 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2662 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
2664 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
2665 msfuy, &
2666 msfvx, &
2667 msfvx_inv, &
2668 msfvy, &
2669 msftx, &
2670 msfty
2672 REAL , INTENT(IN ) :: rdx, &
2673 rdy, &
2674 khdif
2676 ! Local data
2678 INTEGER :: i, j, k, itf, jtf, ktf
2680 INTEGER :: i_start, i_end, j_start, j_end
2682 REAL :: mrdx, mkrdxm, mkrdxp, &
2683 mrdy, mkrdym, mkrdyp
2685 LOGICAL :: specified
2687 !<DESCRIPTION>
2688 !
2689 ! horizontal_diffusion computes the horizontal diffusion tendency
2690 ! on model horizontal coordinate surfaces.
2691 !
2692 !</DESCRIPTION>
2694 specified = .false.
2695 if(config_flags%specified .or. config_flags%nested) specified = .true.
2697 ktf=MIN(kte,kde-1)
2699 IF (name .EQ. 'u') THEN
2701 i_start = its
2702 i_end = ite
2703 j_start = jts
2704 j_end = MIN(jte,jde-1)
2706 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2707 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite)
2708 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2709 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2710 IF ( config_flags%periodic_x ) i_start = its
2711 IF ( config_flags%periodic_x ) i_end = ite
2714 DO j = j_start, j_end
2715 DO k=kts,ktf
2716 DO i = i_start, i_end
2718 ! The interior is grad: (m_x*d/dx), the exterior is div: (m_x*m_y*d/dx(/m_y))
2719 ! setting up different averagings of m^2 partial d/dX and m^2 partial d/dY
2721 mkrdxm=(msftx(i-1,j)/msfty(i-1,j))*mu(i-1,j)*xkmhd(i-1,k,j)*rdx
2722 mkrdxp=(msftx(i,j)/msfty(i,j))*mu(i,j)*xkmhd(i,k,j)*rdx
2723 mrdx=msfux(i,j)*msfuy(i,j)*rdx
2724 mkrdym=( (msfuy(i,j)+msfuy(i,j-1))/(msfux(i,j)+msfux(i,j-1)) )* &
2725 0.25*(mu(i,j)+mu(i,j-1)+mu(i-1,j-1)+mu(i-1,j))* &
2726 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i-1,k,j-1)+xkmhd(i-1,k,j))*rdy
2727 mkrdyp=( (msfuy(i,j)+msfuy(i,j+1))/(msfux(i,j)+msfux(i,j+1)) )* &
2728 0.25*(mu(i,j)+mu(i,j+1)+mu(i-1,j+1)+mu(i-1,j))* &
2729 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j+1)+xkmhd(i-1,k,j+1)+xkmhd(i-1,k,j))*rdy
2730 ! need to do four-corners (t) for diffusion coefficient as there are
2731 ! no values at u,v points
2732 ! msfuy - has to be y as part of d/dY
2733 ! has to be u as we're at a u point
2734 mrdy=msfux(i,j)*msfuy(i,j)*rdy
2736 ! correctly averaged version of rho~ * m^2 *
2737 ! [partial d/dX(partial du^/dX) + partial d/dY(partial du^/dY)]
2738 tendency(i,k,j)=tendency(i,k,j)+( &
2739 mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) &
2740 -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2741 +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) &
2742 -mkrdym*(field(i,k,j )-field(i,k,j-1))))
2743 ENDDO
2744 ENDDO
2745 ENDDO
2747 ELSE IF (name .EQ. 'v')THEN
2749 i_start = its
2750 i_end = MIN(ite,ide-1)
2751 j_start = jts
2752 j_end = jte
2754 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2755 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2756 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2757 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte)
2758 IF ( config_flags%periodic_x ) i_start = its
2759 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2760 IF ( config_flags%polar ) j_start = MAX(jds+1,jts)
2761 IF ( config_flags%polar ) j_end = MIN(jde-1,jte)
2763 DO j = j_start, j_end
2764 DO k=kts,ktf
2765 DO i = i_start, i_end
2767 mkrdxm=( (msfvx(i,j)+msfvx(i-1,j))/(msfvy(i,j)+msfvy(i-1,j)) )* &
2768 0.25*(mu(i,j)+mu(i,j-1)+mu(i-1,j-1)+mu(i-1,j))* &
2769 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i-1,k,j-1)+xkmhd(i-1,k,j))*rdx
2770 mkrdxp=( (msfvx(i,j)+msfvx(i+1,j))/(msfvy(i,j)+msfvy(i+1,j)) )* &
2771 0.25*(mu(i,j)+mu(i,j-1)+mu(i+1,j-1)+mu(i+1,j))* &
2772 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i+1,k,j-1)+xkmhd(i+1,k,j))*rdx
2773 mrdx=msfvx(i,j)*msfvy(i,j)*rdx
2774 mkrdym=(msfty(i,j-1)/msftx(i,j-1))*xkmhd(i,k,j-1)*rdy
2775 mkrdyp=(msfty(i,j)/msftx(i,j))*xkmhd(i,k,j)*rdy
2776 mrdy=msfvx(i,j)*msfvy(i,j)*rdy
2778 tendency(i,k,j)=tendency(i,k,j)+( &
2779 mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) &
2780 -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2781 +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) &
2782 -mkrdym*(field(i,k,j )-field(i,k,j-1))))
2783 ENDDO
2784 ENDDO
2785 ENDDO
2787 ELSE IF (name .EQ. 'w')THEN
2789 i_start = its
2790 i_end = MIN(ite,ide-1)
2791 j_start = jts
2792 j_end = MIN(jte,jde-1)
2794 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2795 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2796 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2797 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2798 IF ( config_flags%periodic_x ) i_start = its
2799 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2801 DO j = j_start, j_end
2802 DO k=kts+1,ktf
2803 DO i = i_start, i_end
2805 mkrdxm=(msfux(i,j)/msfuy(i,j))* &
2806 0.25*(mu(i,j)+mu(i-1,j)+mu(i,j)+mu(i-1,j))* &
2807 0.25*(xkmhd(i,k,j)+xkmhd(i-1,k,j)+xkmhd(i,k-1,j)+xkmhd(i-1,k-1,j))*rdx
2808 mkrdxp=(msfux(i+1,j)/msfuy(i+1,j))* &
2809 0.25*(mu(i+1,j)+mu(i,j)+mu(i+1,j)+mu(i,j))* &
2810 0.25*(xkmhd(i+1,k,j)+xkmhd(i,k,j)+xkmhd(i+1,k-1,j)+xkmhd(i,k-1,j))*rdx
2811 mrdx=msftx(i,j)*msfty(i,j)*rdx
2812 ! mkrdym=(msfvy(i,j)/msfvx(i,j))* &
2813 mkrdym=(msfvy(i,j)*msfvx_inv(i,j))* &
2814 0.25*(mu(i,j)+mu(i,j-1)+mu(i,j)+mu(i,j-1))* &
2815 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i,k-1,j)+xkmhd(i,k-1,j-1))*rdy
2816 ! mkrdyp=(msfvy(i,j+1)/msfvx(i,j+1))* &
2817 mkrdyp=(msfvy(i,j+1)*msfvx_inv(i,j+1))* &
2818 0.25*(mu(i,j+1)+mu(i,j)+mu(i,j+1)+mu(i,j))* &
2819 0.25*(xkmhd(i,k,j+1)+xkmhd(i,k,j)+xkmhd(i,k-1,j+1)+xkmhd(i,k-1,j))*rdy
2820 mrdy=msftx(i,j)*msfty(i,j)*rdy
2822 tendency(i,k,j)=tendency(i,k,j)+( &
2823 mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) &
2824 -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2825 +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) &
2826 -mkrdym*(field(i,k,j )-field(i,k,j-1))))
2827 ENDDO
2828 ENDDO
2829 ENDDO
2831 ELSE
2834 i_start = its
2835 i_end = MIN(ite,ide-1)
2836 j_start = jts
2837 j_end = MIN(jte,jde-1)
2839 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2840 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2841 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2842 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2843 IF ( config_flags%periodic_x ) i_start = its
2844 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2846 DO j = j_start, j_end
2847 DO k=kts,ktf
2848 DO i = i_start, i_end
2850 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
2851 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
2852 mrdx=msftx(i,j)*msfty(i,j)*rdx
2853 ! 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
2854 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
2855 ! 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
2856 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
2857 mrdy=msftx(i,j)*msfty(i,j)*rdy
2859 tendency(i,k,j)=tendency(i,k,j)+( &
2860 mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) &
2861 -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2862 +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) &
2863 -mkrdym*(field(i,k,j )-field(i,k,j-1))))
2864 ENDDO
2865 ENDDO
2866 ENDDO
2868 ENDIF
2870 END SUBROUTINE horizontal_diffusion
2872 !-----------------------------------------------------------------------------------------
2874 SUBROUTINE horizontal_diffusion_3dmp ( name, field, tendency, mu, &
2875 config_flags, base_3d, &
2876 msfux, msfuy, msfvx, msfvx_inv, &
2877 msfvy, msftx, msfty, &
2878 khdif, xkmhd, rdx, rdy, &
2879 ids, ide, jds, jde, kds, kde, &
2880 ims, ime, jms, jme, kms, kme, &
2881 its, ite, jts, jte, kts, kte )
2883 IMPLICIT NONE
2885 ! Input data
2887 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2889 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2890 ims, ime, jms, jme, kms, kme, &
2891 its, ite, jts, jte, kts, kte
2893 CHARACTER(LEN=1) , INTENT(IN ) :: name
2895 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, &
2896 xkmhd, &
2897 base_3d
2899 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2901 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
2903 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
2904 msfuy, &
2905 msfvx, &
2906 msfvx_inv, &
2907 msfvy, &
2908 msftx, &
2909 msfty
2911 REAL , INTENT(IN ) :: rdx, &
2912 rdy, &
2913 khdif
2915 ! Local data
2917 INTEGER :: i, j, k, itf, jtf, ktf
2919 INTEGER :: i_start, i_end, j_start, j_end
2921 REAL :: mrdx, mkrdxm, mkrdxp, &
2922 mrdy, mkrdym, mkrdyp
2924 LOGICAL :: specified
2926 !<DESCRIPTION>
2927 !
2928 ! horizontal_diffusion_3dmp computes the horizontal diffusion tendency
2929 ! on model horizontal coordinate surfaces. This routine computes diffusion
2930 ! a perturbation scalar (field-base_3d).
2931 !
2932 !</DESCRIPTION>
2934 specified = .false.
2935 if(config_flags%specified .or. config_flags%nested) specified = .true.
2937 ktf=MIN(kte,kde-1)
2939 i_start = its
2940 i_end = MIN(ite,ide-1)
2941 j_start = jts
2942 j_end = MIN(jte,jde-1)
2944 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2945 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2946 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2947 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2948 IF ( config_flags%periodic_x ) i_start = its
2949 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2951 DO j = j_start, j_end
2952 DO k=kts,ktf
2953 DO i = i_start, i_end
2955 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
2956 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
2957 mrdx=msftx(i,j)*msfty(i,j)*rdx
2958 ! 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
2959 ! 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
2960 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
2961 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
2962 mrdy=msftx(i,j)*msfty(i,j)*rdy
2964 tendency(i,k,j)=tendency(i,k,j)+( &
2965 mrdx*( mkrdxp*( field(i+1,k,j) -field(i ,k,j) &
2966 -base_3d(i+1,k,j)+base_3d(i ,k,j) ) &
2967 -mkrdxm*( field(i ,k,j) -field(i-1,k,j) &
2968 -base_3d(i ,k,j)+base_3d(i-1,k,j) ) ) &
2969 +mrdy*( mkrdyp*( field(i,k,j+1) -field(i,k,j ) &
2970 -base_3d(i,k,j+1)+base_3d(i,k,j ) ) &
2971 -mkrdym*( field(i,k,j ) -field(i,k,j-1) &
2972 -base_3d(i,k,j )+base_3d(i,k,j-1) ) ) &
2973 )
2974 ENDDO
2975 ENDDO
2976 ENDDO
2978 END SUBROUTINE horizontal_diffusion_3dmp
2980 !-----------------------------------------------------------------------------------------
2982 SUBROUTINE vertical_diffusion ( name, field, tendency, &
2983 config_flags, &
2984 alt, mut, rdn, rdnw, kvdif, &
2985 ids, ide, jds, jde, kds, kde, &
2986 ims, ime, jms, jme, kms, kme, &
2987 its, ite, jts, jte, kts, kte )
2990 IMPLICIT NONE
2992 ! Input data
2994 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2996 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2997 ims, ime, jms, jme, kms, kme, &
2998 its, ite, jts, jte, kts, kte
3000 CHARACTER(LEN=1) , INTENT(IN ) :: name
3002 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3003 INTENT(IN ) :: field, &
3004 alt
3006 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3008 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
3010 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw
3012 REAL , INTENT(IN ) :: kvdif
3014 ! Local data
3016 INTEGER :: i, j, k, itf, jtf, ktf
3017 INTEGER :: i_start, i_end, j_start, j_end
3019 REAL , DIMENSION(its:ite, jts:jte) :: vfluxm, vfluxp, zz
3020 REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3022 REAL :: rdz
3024 LOGICAL :: specified
3026 !<DESCRIPTION>
3027 !
3028 ! vertical_diffusion
3029 ! computes vertical diffusion tendency.
3030 !
3031 !</DESCRIPTION>
3033 specified = .false.
3034 if(config_flags%specified .or. config_flags%nested) specified = .true.
3036 ktf=MIN(kte,kde-1)
3038 IF (name .EQ. 'w')THEN
3041 i_start = its
3042 i_end = MIN(ite,ide-1)
3043 j_start = jts
3044 j_end = MIN(jte,jde-1)
3046 j_loop_w : DO j = j_start, j_end
3048 DO k=kts,ktf-1
3049 DO i = i_start, i_end
3050 vflux(i,k)= (kvdif/alt(i,k,j))*rdnw(k)*(field(i,k+1,j)-field(i,k,j))
3051 ENDDO
3052 ENDDO
3054 DO i = i_start, i_end
3055 vflux(i,ktf)=0.
3056 ENDDO
3058 DO k=kts+1,ktf
3059 DO i = i_start, i_end
3060 tendency(i,k,j)=tendency(i,k,j) &
3061 +rdn(k)*g*g/mut(i,j)/(0.5*(alt(i,k,j)+alt(i,k-1,j))) &
3062 *(vflux(i,k)-vflux(i,k-1))
3063 ENDDO
3064 ENDDO
3066 ENDDO j_loop_w
3068 ELSE IF(name .EQ. 'm')THEN
3070 i_start = its
3071 i_end = MIN(ite,ide-1)
3072 j_start = jts
3073 j_end = MIN(jte,jde-1)
3075 j_loop_s : DO j = j_start, j_end
3077 DO k=kts,ktf-1
3078 DO i = i_start, i_end
3079 vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j))) &
3080 *(field(i,k+1,j)-field(i,k,j))
3081 ENDDO
3082 ENDDO
3084 DO i = i_start, i_end
3085 vflux(i,0)=vflux(i,1)
3086 ENDDO
3088 DO i = i_start, i_end
3089 vflux(i,ktf)=0.
3090 ENDDO
3092 DO k=kts,ktf
3093 DO i = i_start, i_end
3094 tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) &
3095 *rdnw(k)*(vflux(i,k)-vflux(i,k-1))
3096 ENDDO
3097 ENDDO
3099 ENDDO j_loop_s
3101 ENDIF
3103 END SUBROUTINE vertical_diffusion
3106 !-------------------------------------------------------------------------------
3108 SUBROUTINE vertical_diffusion_mp ( field, tendency, config_flags, &
3109 base, &
3110 alt, mut, rdn, rdnw, kvdif, &
3111 ids, ide, jds, jde, kds, kde, &
3112 ims, ime, jms, jme, kms, kme, &
3113 its, ite, jts, jte, kts, kte )
3116 IMPLICIT NONE
3118 ! Input data
3120 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3122 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3123 ims, ime, jms, jme, kms, kme, &
3124 its, ite, jts, jte, kts, kte
3126 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3127 INTENT(IN ) :: field, &
3128 alt
3130 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3132 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
3134 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, &
3135 rdnw, &
3136 base
3138 REAL , INTENT(IN ) :: kvdif
3140 ! Local data
3142 INTEGER :: i, j, k, itf, jtf, ktf
3143 INTEGER :: i_start, i_end, j_start, j_end
3145 REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3147 REAL :: rdz
3149 LOGICAL :: specified
3151 !<DESCRIPTION>
3152 !
3153 ! vertical_diffusion_mp
3154 ! computes vertical diffusion tendency of a perturbation variable
3155 ! (field-base). Note that base as a 1D (k) field.
3156 !
3157 !</DESCRIPTION>
3159 specified = .false.
3160 if(config_flags%specified .or. config_flags%nested) specified = .true.
3162 ktf=MIN(kte,kde-1)
3164 i_start = its
3165 i_end = MIN(ite,ide-1)
3166 j_start = jts
3167 j_end = MIN(jte,jde-1)
3169 j_loop_s : DO j = j_start, j_end
3171 DO k=kts,ktf-1
3172 DO i = i_start, i_end
3173 vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j))) &
3174 *(field(i,k+1,j)-field(i,k,j)-base(k+1)+base(k))
3175 ENDDO
3176 ENDDO
3178 DO i = i_start, i_end
3179 vflux(i,0)=vflux(i,1)
3180 ENDDO
3182 DO i = i_start, i_end
3183 vflux(i,ktf)=0.
3184 ENDDO
3186 DO k=kts,ktf
3187 DO i = i_start, i_end
3188 tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) &
3189 *rdnw(k)*(vflux(i,k)-vflux(i,k-1))
3190 ENDDO
3191 ENDDO
3193 ENDDO j_loop_s
3195 END SUBROUTINE vertical_diffusion_mp
3198 !-------------------------------------------------------------------------------
3200 SUBROUTINE vertical_diffusion_3dmp ( field, tendency, config_flags, &
3201 base_3d, &
3202 alt, mut, rdn, rdnw, kvdif, &
3203 ids, ide, jds, jde, kds, kde, &
3204 ims, ime, jms, jme, kms, kme, &
3205 its, ite, jts, jte, kts, kte )
3208 IMPLICIT NONE
3210 ! Input data
3212 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3214 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3215 ims, ime, jms, jme, kms, kme, &
3216 its, ite, jts, jte, kts, kte
3218 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3219 INTENT(IN ) :: field, &
3220 alt, &
3221 base_3d
3223 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3225 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
3227 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, &
3228 rdnw
3230 REAL , INTENT(IN ) :: kvdif
3232 ! Local data
3234 INTEGER :: i, j, k, itf, jtf, ktf
3235 INTEGER :: i_start, i_end, j_start, j_end
3237 REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3239 REAL :: rdz
3241 LOGICAL :: specified
3243 !<DESCRIPTION>
3244 !
3245 ! vertical_diffusion_3dmp
3246 ! computes vertical diffusion tendency of a perturbation variable
3247 ! (field-base_3d).
3248 !
3249 !</DESCRIPTION>
3251 specified = .false.
3252 if(config_flags%specified .or. config_flags%nested) specified = .true.
3254 ktf=MIN(kte,kde-1)
3256 i_start = its
3257 i_end = MIN(ite,ide-1)
3258 j_start = jts
3259 j_end = MIN(jte,jde-1)
3261 j_loop_s : DO j = j_start, j_end
3263 DO k=kts,ktf-1
3264 DO i = i_start, i_end
3265 vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j))) &
3266 *( field(i,k+1,j) -field(i,k,j) &
3267 -base_3d(i,k+1,j)+base_3d(i,k,j) )
3268 ENDDO
3269 ENDDO
3271 DO i = i_start, i_end
3272 vflux(i,0)=vflux(i,1)
3273 ENDDO
3275 DO i = i_start, i_end
3276 vflux(i,ktf)=0.
3277 ENDDO
3279 DO k=kts,ktf
3280 DO i = i_start, i_end
3281 tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) &
3282 *rdnw(k)*(vflux(i,k)-vflux(i,k-1))
3283 ENDDO
3284 ENDDO
3286 ENDDO j_loop_s
3288 END SUBROUTINE vertical_diffusion_3dmp
3291 !-------------------------------------------------------------------------------
3294 SUBROUTINE vertical_diffusion_u ( field, tendency, &
3295 config_flags, u_base, &
3296 alt, muu, rdn, rdnw, kvdif, &
3297 ids, ide, jds, jde, kds, kde, &
3298 ims, ime, jms, jme, kms, kme, &
3299 its, ite, jts, jte, kts, kte )
3302 IMPLICIT NONE
3304 ! Input data
3306 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3308 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3309 ims, ime, jms, jme, kms, kme, &
3310 its, ite, jts, jte, kts, kte
3312 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3313 INTENT(IN ) :: field, &
3314 alt
3316 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3318 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muu
3320 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw, u_base
3322 REAL , INTENT(IN ) :: kvdif
3324 ! Local data
3326 INTEGER :: i, j, k, itf, jtf, ktf
3327 INTEGER :: i_start, i_end, j_start, j_end
3329 REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3331 REAL :: rdz, zz
3333 LOGICAL :: specified
3335 !<DESCRIPTION>
3336 !
3337 ! vertical_diffusion_u computes vertical diffusion tendency for
3338 ! the u momentum equation. This routine assumes a constant eddy
3339 ! viscosity kvdif.
3340 !
3341 !</DESCRIPTION>
3343 specified = .false.
3344 if(config_flags%specified .or. config_flags%nested) specified = .true.
3346 ktf=MIN(kte,kde-1)
3348 i_start = its
3349 i_end = ite
3350 j_start = jts
3351 j_end = MIN(jte,jde-1)
3353 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
3354 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite)
3355 IF ( config_flags%periodic_x ) i_start = its
3356 IF ( config_flags%periodic_x ) i_end = ite
3359 j_loop_u : DO j = j_start, j_end
3361 DO k=kts,ktf-1
3362 DO i = i_start, i_end
3363 vflux(i,k)=kvdif*rdn(k+1)/(0.25*( alt(i ,k ,j) &
3364 +alt(i-1,k ,j) &
3365 +alt(i ,k+1,j) &
3366 +alt(i-1,k+1,j) ) ) &
3367 *(field(i,k+1,j)-field(i,k,j) &
3368 -u_base(k+1) +u_base(k) )
3369 ENDDO
3370 ENDDO
3372 DO i = i_start, i_end
3373 vflux(i,0)=vflux(i,1)
3374 ENDDO
3376 DO i = i_start, i_end
3377 vflux(i,ktf)=0.
3378 ENDDO
3380 DO k=kts,ktf-1
3381 DO i = i_start, i_end
3382 tendency(i,k,j)=tendency(i,k,j)+ &
3383 g*g*rdnw(k)/muu(i,j)/(0.5*(alt(i-1,k,j)+alt(i,k,j)))* &
3384 (vflux(i,k)-vflux(i,k-1))
3385 ENDDO
3386 ENDDO
3388 ENDDO j_loop_u
3390 END SUBROUTINE vertical_diffusion_u
3392 !-------------------------------------------------------------------------------
3395 SUBROUTINE vertical_diffusion_v ( field, tendency, &
3396 config_flags, v_base, &
3397 alt, muv, rdn, rdnw, kvdif, &
3398 ids, ide, jds, jde, kds, kde, &
3399 ims, ime, jms, jme, kms, kme, &
3400 its, ite, jts, jte, kts, kte )
3403 IMPLICIT NONE
3405 ! Input data
3407 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3409 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3410 ims, ime, jms, jme, kms, kme, &
3411 its, ite, jts, jte, kts, kte
3413 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3414 INTENT(IN ) :: field, &
3415 alt
3416 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw, v_base
3418 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3420 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muv
3422 REAL , INTENT(IN ) :: kvdif
3424 ! Local data
3426 INTEGER :: i, j, k, itf, jtf, ktf, jm1
3427 INTEGER :: i_start, i_end, j_start, j_end
3429 REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3431 REAL :: rdz, zz
3433 LOGICAL :: specified
3435 !<DESCRIPTION>
3436 !
3437 ! vertical_diffusion_v computes vertical diffusion tendency for
3438 ! the v momentum equation. This routine assumes a constant eddy
3439 ! viscosity kvdif.
3440 !
3441 !</DESCRIPTION>
3443 specified = .false.
3444 if(config_flags%specified .or. config_flags%nested) specified = .true.
3446 ktf=MIN(kte,kde-1)
3448 i_start = its
3449 i_end = MIN(ite,ide-1)
3450 j_start = jts
3451 j_end = MIN(jte,jde-1)
3453 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
3454 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte)
3456 j_loop_v : DO j = j_start, j_end
3457 ! jm1 = max(j-1,1)
3458 jm1 = j-1
3460 DO k=kts,ktf-1
3461 DO i = i_start, i_end
3462 vflux(i,k)=kvdif*rdn(k+1)/(0.25*( alt(i,k ,j ) &
3463 +alt(i,k ,jm1) &
3464 +alt(i,k+1,j ) &
3465 +alt(i,k+1,jm1) ) ) &
3466 *(field(i,k+1,j)-field(i,k,j) &
3467 -v_base(k+1) +v_base(k) )
3468 ENDDO
3469 ENDDO
3471 DO i = i_start, i_end
3472 vflux(i,0)=vflux(i,1)
3473 ENDDO
3475 DO i = i_start, i_end
3476 vflux(i,ktf)=0.
3477 ENDDO
3479 DO k=kts,ktf-1
3480 DO i = i_start, i_end
3481 tendency(i,k,j)=tendency(i,k,j)+ &
3482 g*g*rdnw(k)/muv(i,j)/(0.5*(alt(i,k,jm1)+alt(i,k,j)))* &
3483 (vflux(i,k)-vflux(i,k-1))
3484 ENDDO
3485 ENDDO
3487 ENDDO j_loop_v
3489 END SUBROUTINE vertical_diffusion_v
3491 !*************** end new mass coordinate routines
3493 !-------------------------------------------------------------------------------
3495 SUBROUTINE calculate_full ( rfield, rfieldb, rfieldp, &
3496 ids, ide, jds, jde, kds, kde, &
3497 ims, ime, jms, jme, kms, kme, &
3498 its, ite, jts, jte, kts, kte )
3500 IMPLICIT NONE
3502 ! Input data
3504 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3505 ims, ime, jms, jme, kms, kme, &
3506 its, ite, jts, jte, kts, kte
3508 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rfieldb, &
3509 rfieldp
3511 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: rfield
3513 ! Local indices.
3515 INTEGER :: i, j, k, itf, jtf, ktf
3517 !<DESCRIPTION>
3518 !
3519 ! calculate_full
3520 ! calculates full 3D field from pertubation and base field.
3521 !
3522 !</DESCRIPTION>
3524 itf=MIN(ite,ide-1)
3525 jtf=MIN(jte,jde-1)
3526 ktf=MIN(kte,kde-1)
3528 DO j=jts,jtf
3529 DO k=kts,ktf
3530 DO i=its,itf
3531 rfield(i,k,j)=rfieldb(i,k,j)+rfieldp(i,k,j)
3532 ENDDO
3533 ENDDO
3534 ENDDO
3536 END SUBROUTINE calculate_full
3538 !------------------------------------------------------------------------------
3540 SUBROUTINE coriolis ( ru, rv, rw, ru_tend, rv_tend, rw_tend, &
3541 config_flags, &
3542 msftx, msfty, msfux, msfuy, &
3543 msfvx, msfvy, &
3544 f, e, sina, cosa, fzm, fzp, &
3545 ids, ide, jds, jde, kds, kde, &
3546 ims, ime, jms, jme, kms, kme, &
3547 its, ite, jts, jte, kts, kte )
3549 IMPLICIT NONE
3551 ! Input data
3553 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
3555 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3556 ims, ime, jms, jme, kms, kme, &
3557 its, ite, jts, jte, kts, kte
3559 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: ru_tend, &
3560 rv_tend, &
3561 rw_tend
3562 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: ru, &
3563 rv, &
3564 rw
3566 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
3567 msfuy, &
3568 msfvx, &
3569 msfvy, &
3570 msftx, &
3571 msfty
3573 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: f, &
3574 e, &
3575 sina, &
3576 cosa
3578 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
3579 fzp
3581 ! Local indices.
3583 INTEGER :: i, j , k, ktf
3584 INTEGER :: i_start, i_end, j_start, j_end
3586 LOGICAL :: specified
3588 !<DESCRIPTION>
3589 !
3590 ! coriolis calculates the large timestep tendency terms in the
3591 ! u, v, and w momentum equations arise from the coriolis force.
3592 !
3593 !</DESCRIPTION>
3595 specified = .false.
3596 if(config_flags%specified .or. config_flags%nested) specified = .true.
3598 ktf=MIN(kte,kde-1)
3600 ! coriolis for u-momentum equation
3602 ! Notes on map scale factor
3603 ! cosa, sina are related to rotating the coordinate frame if desired
3604 ! generally sina=0, cosa=1
3605 ! ADT eqn 44, RHS terms 6 and 7: -2 mu w omega cos(lat)/my
3606 ! + 2 mu v omega sin(lat)/my
3607 ! Define f=2 omega sin(lat), e=2 omega cos(lat)
3608 ! => terms are: -e mu w / my + f mu v / my
3609 ! rv = mu v / mx ; rw = mu w / my
3610 ! => terms are: -e rw + f rv *mx / my
3612 i_start = its
3613 i_end = ite
3614 IF ( config_flags%open_xs .or. specified .or. &
3615 config_flags%nested) i_start = MAX(ids+1,its)
3616 IF ( config_flags%open_xe .or. specified .or. &
3617 config_flags%nested) i_end = MIN(ide-1,ite)
3618 IF ( config_flags%periodic_x ) i_start = its
3619 IF ( config_flags%periodic_x ) i_end = ite
3621 DO j = jts, MIN(jte,jde-1)
3623 DO k=kts,ktf
3624 DO i = i_start, i_end
3626 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)) &
3627 *0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
3628 - 0.5*(e(i,j)+e(i-1,j))*0.5*(cosa(i,j)+cosa(i-1,j)) &
3629 *0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
3631 ENDDO
3632 ENDDO
3634 ! boundary loops for coriolis not needed for open bdy (commented out 20100611 JD)
3635 ! IF ( (config_flags%open_xs) .and. (its == ids) ) THEN
3637 ! DO k=kts,ktf
3638 !
3639 ! ru_tend(its,k,j)=ru_tend(its,k,j) + (msfux(its,j)/msfuy(its,j))*0.5*(f(its,j)+f(its,j)) &
3640 ! *0.25*(rv(its,k,j+1)+rv(its,k,j+1)+rv(its,k,j)+rv(its,k,j)) &
3641 ! - 0.5*(e(its,j)+e(its,j))*0.5*(cosa(its,j)+cosa(its,j)) &
3642 ! *0.25*(rw(its,k+1,j)+rw(its,k,j)+rw(its,k+1,j)+rw(its,k,j))
3644 ! ENDDO
3646 ! ENDIF
3648 ! IF ( (config_flags%open_xe) .and. (ite == ide) ) THEN
3650 ! DO k=kts,ktf
3651 !
3652 ! 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)) &
3653 ! *0.25*(rv(ite-1,k,j+1)+rv(ite-1,k,j+1)+rv(ite-1,k,j)+rv(ite-1,k,j)) &
3654 ! - 0.5*(e(ite-1,j)+e(ite-1,j))*0.5*(cosa(ite-1,j)+cosa(ite-1,j)) &
3655 ! *0.25*(rw(ite-1,k+1,j)+rw(ite-1,k,j)+rw(ite-1,k+1,j)+rw(ite-1,k,j))
3657 ! ENDDO
3659 ! ENDIF
3661 ENDDO
3663 ! coriolis term for v-momentum equation
3665 ! Notes on map scale factors
3666 ! ADT eqn 45, RHS terms 6 and 6b [0 for sina=0]: -2 mu u omega sin(lat)/mx + ?
3667 ! Define f=2 omega sin(lat), e=2 omega cos(lat)
3668 ! => terms are: -f mu u / mx
3669 ! ru = mu u / my ; rw = mu w / my
3670 ! => terms are: -f ru *my / mx + ?
3672 j_start = jts
3673 j_end = jte
3675 IF ( config_flags%open_ys .or. specified .or. &
3676 config_flags%nested .or. config_flags%polar) j_start = MAX(jds+1,jts)
3677 IF ( config_flags%open_ye .or. specified .or. &
3678 config_flags%nested .or. config_flags%polar) j_end = MIN(jde-1,jte)
3680 ! boundary loops for coriolis not needed for open bdy (commented out 20100611 JD)
3681 ! IF ( (config_flags%open_ys) .and. (jts == jds) ) THEN
3683 ! DO k=kts,ktf
3684 ! DO i=its,MIN(ide-1,ite)
3685 !
3686 ! rv_tend(i,k,jts)=rv_tend(i,k,jts) - (msfvy(i,jts)/msfvx(i,jts))*0.5*(f(i,jts)+f(i,jts)) &
3687 ! *0.25*(ru(i,k,jts)+ru(i+1,k,jts)+ru(i,k,jts)+ru(i+1,k,jts)) &
3688 ! + (msfvy(i,jts)/msfvx(i,jts))*0.5*(e(i,jts)+e(i,jts))*0.5*(sina(i,jts)+sina(i,jts)) &
3689 ! *0.25*(rw(i,k+1,jts)+rw(i,k,jts)+rw(i,k+1,jts)+rw(i,k,jts))
3691 ! ENDDO
3692 ! ENDDO
3694 ! ENDIF
3696 DO j=j_start, j_end
3697 DO k=kts,ktf
3698 DO i=its,MIN(ide-1,ite)
3700 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)) &
3701 *0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
3702 + (msfvy(i,j)/msfvx(i,j))*0.5*(e(i,j)+e(i,j-1))*0.5*(sina(i,j)+sina(i,j-1)) &
3703 *0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
3705 ENDDO
3706 ENDDO
3707 ENDDO
3710 ! boundary loops for coriolis not needed for open bdy (commented out 20100611 JD)
3711 ! IF ( (config_flags%open_ye) .and. (jte == jde) ) THEN
3713 ! DO k=kts,ktf
3714 ! DO i=its,MIN(ide-1,ite)
3715 !
3716 ! 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)) &
3717 ! *0.25*(ru(i,k,jte-1)+ru(i+1,k,jte-1)+ru(i,k,jte-1)+ru(i+1,k,jte-1)) &
3718 ! + (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)) &
3719 ! *0.25*(rw(i,k+1,jte-1)+rw(i,k,jte-1)+rw(i,k+1,jte-1)+rw(i,k,jte-1))
3721 ! ENDDO
3722 ! ENDDO
3724 ! ENDIF
3726 ! coriolis term for w-mometum
3728 ! Notes on map scale factors
3729 ! ADT eqn 46/my, RHS terms 5 and 5b [0 for sina=0]: 2 mu u omega cos(lat)/my +?
3730 ! Define e=2 omega cos(lat)
3731 ! => terms are: e mu u / my + ???
3732 ! ru = mu u / my ; ru = mu v / mx
3733 ! => terms are: e ru + ???
3735 DO j=jts,MIN(jte, jde-1)
3736 DO k=kts+1,ktf
3737 DO i=its,MIN(ite, ide-1)
3739 rw_tend(i,k,j)=rw_tend(i,k,j) + e(i,j)* &
3740 (cosa(i,j)*0.5*(fzm(k)*(ru(i,k,j)+ru(i+1,k,j)) &
3741 +fzp(k)*(ru(i,k-1,j)+ru(i+1,k-1,j))) &
3742 -(msftx(i,j)/msfty(i,j))* &
3743 sina(i,j)*0.5*(fzm(k)*(rv(i,k,j)+rv(i,k,j+1)) &
3744 +fzp(k)*(rv(i,k-1,j)+rv(i,k-1,j+1))))
3746 ENDDO
3747 ENDDO
3748 ENDDO
3750 END SUBROUTINE coriolis
3752 !------------------------------------------------------------------------------
3754 SUBROUTINE perturbation_coriolis ( ru_in, rv_in, rw, ru_tend, rv_tend, rw_tend, &
3755 config_flags, &
3756 u_base, v_base, z_base, &
3757 muu, muv, phb, ph, &
3758 msftx, msfty, msfux, msfuy, msfvx, msfvy, &
3759 f, e, sina, cosa, fzm, fzp, &
3760 ids, ide, jds, jde, kds, kde, &
3761 ims, ime, jms, jme, kms, kme, &
3762 its, ite, jts, jte, kts, kte )
3764 IMPLICIT NONE
3766 ! Input data
3768 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
3770 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3771 ims, ime, jms, jme, kms, kme, &
3772 its, ite, jts, jte, kts, kte
3774 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: ru_tend, &
3775 rv_tend, &
3776 rw_tend
3777 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: ru_in, &
3778 rv_in, &
3779 rw, &
3780 ph, &
3781 phb
3784 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
3785 msfuy, &
3786 msfvx, &
3787 msfvy, &
3788 msftx, &
3789 msfty
3791 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: f, &
3792 e, &
3793 sina, &
3794 cosa
3796 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muu, &
3797 muv
3800 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
3801 fzp
3803 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: u_base, &
3804 v_base, &
3805 z_base
3807 ! Local storage
3809 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) :: ru, &
3810 rv
3812 REAL :: z_at_u, z_at_v, wkp1, wk, wkm1
3814 ! Local indices.
3816 INTEGER :: i, j , k, ktf
3817 INTEGER :: i_start, i_end, j_start, j_end
3819 LOGICAL :: specified
3821 !<DESCRIPTION>
3822 !
3823 ! perturbation_coriolis calculates the large timestep tendency terms in the
3824 ! u, v, and w momentum equations arise from the coriolis force. This version
3825 ! subtracts off the horizontal velocities from the initial sounding when
3826 ! computing the forcing terms, hence "perturbation" coriolis.
3827 !
3828 !</DESCRIPTION>
3830 specified = .false.
3831 if(config_flags%specified .or. config_flags%nested) specified = .true.
3833 ktf=MIN(kte,kde-1)
3835 ! coriolis for u-momentum equation
3837 i_start = its
3838 i_end = ite
3839 IF ( config_flags%open_xs .or. specified .or. &
3840 config_flags%nested) i_start = MAX(ids+1,its)
3841 IF ( config_flags%open_xe .or. specified .or. &
3842 config_flags%nested) i_end = MIN(ide-1,ite)
3843 IF ( config_flags%periodic_x ) i_start = its
3844 IF ( config_flags%periodic_x ) i_end = ite
3846 ! compute perturbation mu*v for use in u momentum equation
3848 DO j = jts, MIN(jte,jde-1)+1
3849 DO k=kts+1,ktf-1
3850 DO i = i_start-1, i_end
3851 z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
3852 +phb(i,k,j-1)+phb(i,k+1,j-1) &
3853 +ph(i,k,j )+ph(i,k+1,j ) &
3854 +ph(i,k,j-1)+ph(i,k+1,j-1))/g
3855 wkp1 = min(1.,max(0.,z_at_v-z_base(k))/(z_base(k+1)-z_base(k)))
3856 wkm1 = min(1.,max(0.,z_base(k)-z_at_v)/(z_base(k)-z_base(k-1)))
3857 wk = 1.-wkp1-wkm1
3858 rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*( &
3859 wkm1*v_base(k-1) &
3860 +wk *v_base(k ) &
3861 +wkp1*v_base(k+1) )
3862 ENDDO
3863 ENDDO
3864 ENDDO
3867 ! pick up top and bottom v
3869 DO j = jts, MIN(jte,jde-1)+1
3870 DO i = i_start-1, i_end
3872 k = kts
3873 z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
3874 +phb(i,k,j-1)+phb(i,k+1,j-1) &
3875 +ph(i,k,j )+ph(i,k+1,j ) &
3876 +ph(i,k,j-1)+ph(i,k+1,j-1))/g
3877 wkp1 = min(1.,max(0.,z_at_v-z_base(k))/(z_base(k+1)-z_base(k)))
3878 wk = 1.-wkp1
3879 rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*( &
3880 +wk *v_base(k ) &
3881 +wkp1*v_base(k+1) )
3883 k = ktf
3884 z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
3885 +phb(i,k,j-1)+phb(i,k+1,j-1) &
3886 +ph(i,k,j )+ph(i,k+1,j ) &
3887 +ph(i,k,j-1)+ph(i,k+1,j-1))/g
3888 wkm1 = min(1.,max(0.,z_base(k)-z_at_v)/(z_base(k)-z_base(k-1)))
3889 wk = 1.-wkm1
3890 rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*( &
3891 wkm1*v_base(k-1) &
3892 +wk *v_base(k ) )
3894 ENDDO
3895 ENDDO
3897 ! compute coriolis forcing for u
3899 ! Map scale factors: see comments above for Coriolis
3901 DO j = jts, MIN(jte,jde-1)
3903 DO k=kts,ktf
3904 DO i = i_start, i_end
3905 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)) &
3906 *0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
3907 - 0.5*(e(i,j)+e(i-1,j))*0.5*(cosa(i,j)+cosa(i-1,j)) &
3908 *0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
3909 ENDDO
3910 ENDDO
3912 ! boundary loops for perturbation coriolis is needed for open bdy (20110225 JD)
3913 IF ( (config_flags%open_xs) .and. (its == ids) ) THEN
3915 DO k=kts,ktf
3917 ru_tend(its,k,j)=ru_tend(its,k,j) + (msfux(its,j)/msfuy(its,j))*0.5*(f(its,j)+f(its,j)) &
3918 *0.25*(rv(its,k,j+1)+rv(its,k,j+1)+rv(its,k,j)+rv(its,k,j)) &
3919 - 0.5*(e(its,j)+e(its,j))*0.5*(cosa(its,j)+cosa(its,j)) &
3920 *0.25*(rw(its,k+1,j)+rw(its,k,j)+rw(its,k+1,j)+rw(its,k,j))
3922 ENDDO
3924 ENDIF
3926 IF ( (config_flags%open_xe) .and. (ite == ide) ) THEN
3928 DO k=kts,ktf
3930 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)) &
3931 *0.25*(rv(ite-1,k,j+1)+rv(ite-1,k,j+1)+rv(ite-1,k,j)+rv(ite-1,k,j)) &
3932 - 0.5*(e(ite-1,j)+e(ite-1,j))*0.5*(cosa(ite-1,j)+cosa(ite-1,j)) &
3933 *0.25*(rw(ite-1,k+1,j)+rw(ite-1,k,j)+rw(ite-1,k+1,j)+rw(ite-1,k,j))
3935 ENDDO
3937 ENDIF
3939 ENDDO
3941 ! coriolis term for v-momentum equation
3942 ! Map scale factors: see comments above for Coriolis
3944 j_start = jts
3945 j_end = jte
3947 IF ( config_flags%open_ys .or. specified .or. &
3948 config_flags%nested .or. config_flags%polar) j_start = MAX(jds+1,jts)
3949 IF ( config_flags%open_ye .or. specified .or. &
3950 config_flags%nested .or. config_flags%polar) j_end = MIN(jde-1,jte)
3952 ! compute perturbation mu*u for use in v momentum equation
3954 DO j = j_start-1,j_end
3955 DO k=kts+1,ktf-1
3956 DO i = its, MIN(ite,ide-1)+1
3957 z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
3958 +phb(i-1,k,j)+phb(i-1,k+1,j) &
3959 +ph(i ,k,j)+ph(i ,k+1,j) &
3960 +ph(i-1,k,j)+ph(i-1,k+1,j))/g
3961 wkp1 = min(1.,max(0.,z_at_u-z_base(k))/(z_base(k+1)-z_base(k)))
3962 wkm1 = min(1.,max(0.,z_base(k)-z_at_u)/(z_base(k)-z_base(k-1)))
3963 wk = 1.-wkp1-wkm1
3964 ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*( &
3965 wkm1*u_base(k-1) &
3966 +wk *u_base(k ) &
3967 +wkp1*u_base(k+1) )
3968 ENDDO
3969 ENDDO
3970 ENDDO
3972 ! pick up top and bottom u
3974 DO j = j_start-1,j_end
3975 DO i = its, MIN(ite,ide-1)+1
3977 k = kts
3978 z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
3979 +phb(i-1,k,j)+phb(i-1,k+1,j) &
3980 +ph(i ,k,j)+ph(i ,k+1,j) &
3981 +ph(i-1,k,j)+ph(i-1,k+1,j))/g
3982 wkp1 = min(1.,max(0.,z_at_u-z_base(k))/(z_base(k+1)-z_base(k)))
3983 wk = 1.-wkp1
3984 ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*( &
3985 +wk *u_base(k ) &
3986 +wkp1*u_base(k+1) )
3989 k = ktf
3990 z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
3991 +phb(i-1,k,j)+phb(i-1,k+1,j) &
3992 +ph(i ,k,j)+ph(i ,k+1,j) &
3993 +ph(i-1,k,j)+ph(i-1,k+1,j))/g
3994 wkm1 = min(1.,max(0.,z_base(k)-z_at_u)/(z_base(k)-z_base(k-1)))
3995 wk = 1.-wkm1
3996 ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*( &
3997 wkm1*u_base(k-1) &
3998 +wk *u_base(k ) )
4000 ENDDO
4001 ENDDO
4003 ! compute coriolis forcing for v momentum equation
4004 ! Map scale factors: see comments above for Coriolis
4006 ! boundary loops for perturbation coriolis is needed for open bdy (20110225 JD)
4007 IF ( (config_flags%open_ys) .and. (jts == jds) ) THEN
4009 DO k=kts,ktf
4010 DO i=its,MIN(ide-1,ite)
4012 rv_tend(i,k,jts)=rv_tend(i,k,jts) - (msfvy(i,jts)/msfvx(i,jts))*0.5*(f(i,jts)+f(i,jts)) &
4013 *0.25*(ru(i,k,jts)+ru(i+1,k,jts)+ru(i,k,jts)+ru(i+1,k,jts)) &
4014 + (msfvy(i,jts)/msfvx(i,jts))*0.5*(e(i,jts)+e(i,jts))*0.5*(sina(i,jts)+sina(i,jts)) &
4015 *0.25*(rw(i,k+1,jts)+rw(i,k,jts)+rw(i,k+1,jts)+rw(i,k,jts))
4017 ENDDO
4018 ENDDO
4020 ENDIF
4022 DO j=j_start, j_end
4023 DO k=kts,ktf
4024 DO i=its,MIN(ide-1,ite)
4026 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)) &
4027 *0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
4028 + (msfvy(i,j)/msfvx(i,j))*0.5*(e(i,j)+e(i,j-1))*0.5*(sina(i,j)+sina(i,j-1)) &
4029 *0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
4031 ENDDO
4032 ENDDO
4033 ENDDO
4036 ! boundary loops for perturbation coriolis is needed for open bdy (20110225 JD)
4037 IF ( (config_flags%open_ye) .and. (jte == jde) ) THEN
4039 DO k=kts,ktf
4040 DO i=its,MIN(ide-1,ite)
4042 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)) &
4043 *0.25*(ru(i,k,jte-1)+ru(i+1,k,jte-1)+ru(i,k,jte-1)+ru(i+1,k,jte-1)) &
4044 + (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)) &
4045 *0.25*(rw(i,k+1,jte-1)+rw(i,k,jte-1)+rw(i,k+1,jte-1)+rw(i,k,jte-1))
4047 ENDDO
4048 ENDDO
4050 ENDIF
4052 ! coriolis term for w-mometum
4053 ! Map scale factors: see comments above for Coriolis
4055 DO j=jts,MIN(jte, jde-1)
4056 DO k=kts+1,ktf
4057 DO i=its,MIN(ite, ide-1)
4059 rw_tend(i,k,j)=rw_tend(i,k,j) + e(i,j)* &
4060 (cosa(i,j)*0.5*(fzm(k)*(ru(i,k,j)+ru(i+1,k,j)) &
4061 +fzp(k)*(ru(i,k-1,j)+ru(i+1,k-1,j))) &
4062 -(msftx(i,j)/msfty(i,j))*sina(i,j)*0.5*(fzm(k)*(rv(i,k,j)+rv(i,k,j+1)) &
4063 +fzp(k)*(rv(i,k-1,j)+rv(i,k-1,j+1))))
4065 ENDDO
4066 ENDDO
4067 ENDDO
4069 END SUBROUTINE perturbation_coriolis
4071 !------------------------------------------------------------------------------
4073 SUBROUTINE curvature ( ru, rv, rw, u, v, w, ru_tend, rv_tend, rw_tend, &
4074 config_flags, &
4075 msfux, msfuy, msfvx, msfvy, msftx, msfty, &
4076 xlat, fzm, fzp, rdx, rdy, &
4077 ids, ide, jds, jde, kds, kde, &
4078 ims, ime, jms, jme, kms, kme, &
4079 its, ite, jts, jte, kts, kte )
4082 IMPLICIT NONE
4084 ! Input data
4086 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
4088 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4089 ims, ime, jms, jme, kms, kme, &
4090 its, ite, jts, jte, kts, kte
4092 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4093 INTENT(INOUT) :: ru_tend, &
4094 rv_tend, &
4095 rw_tend
4097 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4098 INTENT(IN ) :: ru, &
4099 rv, &
4100 rw, &
4101 u, &
4102 v, &
4103 w
4105 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
4106 msfuy, &
4107 msfvx, &
4108 msfvy, &
4109 msftx, &
4110 msfty, &
4111 xlat
4113 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
4114 fzp
4116 REAL , INTENT(IN ) :: rdx, &
4117 rdy
4119 ! Local data
4121 ! INTEGER :: i, j, k, itf, jtf, ktf, kp1, im, ip, jm, jp
4122 INTEGER :: i, j, k, itf, jtf, ktf
4123 INTEGER :: i_start, i_end, j_start, j_end
4124 ! INTEGER :: irmin, irmax, jrmin, jrmax
4126 REAL , DIMENSION( its-1:ite , kts:kte, jts-1:jte ) :: vxgm
4128 LOGICAL :: specified
4130 !<DESCRIPTION>
4131 !
4132 ! curvature calculates the large timestep tendency terms in the
4133 ! u, v, and w momentum equations arise from the curvature terms.
4134 !
4135 !</DESCRIPTION>
4137 specified = .false.
4138 if(config_flags%specified .or. config_flags%nested) specified = .true.
4140 itf=MIN(ite,ide-1)
4141 jtf=MIN(jte,jde-1)
4142 ktf=MIN(kte,kde-1)
4144 ! irmin = ims
4145 ! irmax = ime
4146 ! jrmin = jms
4147 ! jrmax = jme
4148 ! IF ( config_flags%open_xs ) irmin = ids
4149 ! IF ( config_flags%open_xe ) irmax = ide-1
4150 ! IF ( config_flags%open_ys ) jrmin = jds
4151 ! IF ( config_flags%open_ye ) jrmax = jde-1
4153 ! Define v cross grad m at scalar points - vxgm(i,j)
4155 i_start = its-1
4156 i_end = ite
4157 j_start = jts-1
4158 j_end = jte
4160 IF ( ( config_flags%open_xs .or. specified .or. &
4161 config_flags%nested) .and. (its == ids) ) i_start = its
4162 IF ( ( config_flags%open_xe .or. specified .or. &
4163 config_flags%nested) .and. (ite == ide) ) i_end = ite-1
4164 IF ( ( config_flags%open_ys .or. specified .or. &
4165 config_flags%nested .or. config_flags%polar) .and. (jts == jds) ) j_start = jts
4166 IF ( ( config_flags%open_ye .or. specified .or. &
4167 config_flags%nested .or. config_flags%polar) .and. (jte == jde) ) j_end = jte-1
4168 IF ( config_flags%periodic_x ) i_start = its-1
4169 IF ( config_flags%periodic_x ) i_end = ite
4171 DO j=j_start, j_end
4172 DO k=kts,ktf
4173 DO i=i_start, i_end
4174 ! Map scale factor notes:
4175 ! msf...y is constant everywhere for cylindrical map projection
4176 ! msf...x varies with y only
4177 ! But we know that this is not = 0 for cylindrical,
4178 ! therefore use msfvX in 1st line
4179 ! which => by symmetry use msfuY in 2nd line - ???
4180 vxgm(i,k,j)=0.5*(u(i,k,j)+u(i+1,k,j))*(msfvx(i,j+1)-msfvx(i,j))*rdy - &
4181 0.5*(v(i,k,j)+v(i,k,j+1))*(msfuy(i+1,j)-msfuy(i,j))*rdx
4182 ENDDO
4183 ENDDO
4184 ENDDO
4186 ! Pick up the boundary rows for open (radiation) lateral b.c.
4187 ! Rather crude at present, we are assuming there is no
4188 ! variation in this term at the boundary.
4190 IF ( ( config_flags%open_xs .or. (specified .AND. .NOT. config_flags%periodic_x) .or. &
4191 config_flags%nested) .and. (its == ids) ) THEN
4193 DO j = jts, jte-1
4194 DO k = kts, ktf
4195 vxgm(its-1,k,j) = vxgm(its,k,j)
4196 ENDDO
4197 ENDDO
4199 ENDIF
4201 IF ( ( config_flags%open_xe .or. (specified .AND. .NOT. config_flags%periodic_x) .or. &
4202 config_flags%nested) .and. (ite == ide) ) THEN
4204 DO j = jts, jte-1
4205 DO k = kts, ktf
4206 vxgm(ite,k,j) = vxgm(ite-1,k,j)
4207 ENDDO
4208 ENDDO
4210 ENDIF
4212 ! Polar boundary condition:
4213 ! The following change is needed in case one tries using the vxgm route with
4214 ! polar B.C.'s in the future, but not needed if 'tan' used
4215 IF ( ( config_flags%open_ys .or. specified .or. &
4216 config_flags%nested .or. config_flags%polar) .and. (jts == jds) ) THEN
4218 DO k = kts, ktf
4219 DO i = its-1, ite
4220 vxgm(i,k,jts-1) = vxgm(i,k,jts)
4221 ENDDO
4222 ENDDO
4224 ENDIF
4226 ! Polar boundary condition:
4227 ! The following change is needed in case one tries using the vxgm route with
4228 ! polar B.C.'s in the future, but not needed if 'tan' used
4229 IF ( ( config_flags%open_ye .or. specified .or. &
4230 config_flags%nested .or. config_flags%polar) .and. (jte == jde) ) THEN
4232 DO k = kts, ktf
4233 DO i = its-1, ite
4234 vxgm(i,k,jte) = vxgm(i,k,jte-1)
4235 ENDDO
4236 ENDDO
4238 ENDIF
4240 ! curvature term for u momentum eqn.
4242 ! Map scale factor notes:
4243 ! ADT eqn 44, RHS terms 4 and 5, in cylindrical: mu u v tan(lat)/(a my)
4244 ! - mu u w /(a my)
4245 ! ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
4246 ! => terms are:
4247 ! (mx/my)*u rv tan(lat) / a - u rw / a = (u/a)*[(mx/my) rv tan(lat) - rw]
4248 ! ru v tan(lat) / a - u rw / a
4249 ! xlat defined with end points half grid space from pole,
4250 ! hence are on u latitude points
4252 i_start = its
4253 IF ( config_flags%open_xs .or. specified .or. &
4254 config_flags%nested) i_start = MAX ( ids+1 , its )
4255 IF ( config_flags%open_xe .or. specified .or. &
4256 config_flags%nested) i_end = MIN ( ide-1 , ite )
4257 IF ( config_flags%periodic_x ) i_start = its
4258 IF ( config_flags%periodic_x ) i_end = ite
4260 ! Polar boundary condition
4261 IF ((config_flags%map_proj == 6) .OR. (config_flags%polar)) THEN
4263 DO j=jts,MIN(jde-1,jte)
4264 DO k=kts,ktf
4265 DO i=i_start,i_end
4267 ru_tend(i,k,j)=ru_tend(i,k,j) + u(i,k,j)*reradius* ( &
4268 (msfux(i,j)/msfuy(i,j))*0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+ &
4269 rv(i-1,k,j)+rv(i,k,j))*tan(xlat(i,j)*degrad) &
4270 - 0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j)) )
4271 ENDDO
4272 ENDDO
4273 ENDDO
4275 ELSE ! normal code
4278 DO j=jts,MIN(jde-1,jte)
4279 DO k=kts,ktf
4280 DO i=i_start,i_end
4282 ru_tend(i,k,j)=ru_tend(i,k,j) + 0.5*(vxgm(i,k,j)+vxgm(i-1,k,j)) &
4283 *0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
4284 - u(i,k,j)*reradius &
4285 *0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
4287 ENDDO
4288 ENDDO
4289 ENDDO
4291 END IF
4293 ! curvature term for v momentum eqn.
4295 ! Map scale factor notes
4296 ! ADT eqn 45, RHS terms 4 and 5, in cylindrical: - mu u*u tan(lat)/(a mx)
4297 ! - mu v w /(a mx)
4298 ! ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
4299 ! terms are:
4300 ! - (my/mx)*u ru tan(lat) / a - (my/mx)*v rw / a
4301 ! = - [my/(mx*a)]*[u ru tan(lat) + v rw]
4302 ! - (1/a)*[(my/mx)*u ru tan(lat) + w rv]
4303 ! xlat defined with end points half grid space from pole, hence are on
4304 ! u latitude points => av here
4305 !
4306 ! in original wrf, there was a sign error for the rw contribution
4308 j_start = jts
4309 IF ( config_flags%open_ys .or. specified .or. &
4310 config_flags%nested .or. config_flags%polar) j_start = MAX ( jds+1 , jts )
4311 IF ( config_flags%open_ye .or. specified .or. &
4312 config_flags%nested .or. config_flags%polar) j_end = MIN ( jde-1 , jte )
4314 IF ((config_flags%map_proj == 6) .OR. (config_flags%polar)) THEN
4316 DO j=j_start,j_end
4317 DO k=kts,ktf
4318 DO i=its,MIN(ite,ide-1)
4319 rv_tend(i,k,j)=rv_tend(i,k,j) - (msfvy(i,j)/msfvx(i,j))*reradius* ( &
4320 0.25*(u(i,k,j)+u(i+1,k,j)+u(i,k,j-1)+u(i+1,k,j-1))* &
4321 tan((xlat(i,j)+xlat(i,j-1))*0.5*degrad)* &
4322 0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
4323 + v(i,k,j)*0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+ &
4324 rw(i,k+1,j)+rw(i,k,j)) )
4325 ENDDO
4326 ENDDO
4327 ENDDO
4329 ELSE ! normal code
4331 DO j=j_start,j_end
4332 DO k=kts,ktf
4333 DO i=its,MIN(ite,ide-1)
4335 rv_tend(i,k,j)=rv_tend(i,k,j) - 0.5*(vxgm(i,k,j)+vxgm(i,k,j-1)) &
4336 *0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
4337 - (msfvy(i,j)/msfvx(i,j))*v(i,k,j)*reradius &
4338 *0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
4340 ENDDO
4341 ENDDO
4342 ENDDO
4344 END IF
4346 ! curvature term for vertical momentum eqn.
4348 ! Notes on map scale factors:
4349 ! ADT eqn 46, RHS term 4: [mu/(a my)]*[u*u + v*v]
4350 ! ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
4351 ! terms are: u ru / a + (mx/my)v rv / a
4353 DO j=jts,MIN(jte,jde-1)
4354 DO k=MAX(2,kts),ktf
4355 DO i=its,MIN(ite,ide-1)
4357 rw_tend(i,k,j)=rw_tend(i,k,j) + reradius* &
4358 (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))) &
4359 *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))) &
4360 +(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))) &
4361 *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))))
4363 ENDDO
4364 ENDDO
4365 ENDDO
4367 END SUBROUTINE curvature
4369 !------------------------------------------------------------------------------
4371 SUBROUTINE decouple ( rr, rfield, field, name, config_flags, &
4372 fzm, fzp, &
4373 ids, ide, jds, jde, kds, kde, &
4374 ims, ime, jms, jme, kms, kme, &
4375 its, ite, jts, jte, kts, kte )
4377 IMPLICIT NONE
4379 ! Input data
4381 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
4383 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4384 ims, ime, jms, jme, kms, kme, &
4385 its, ite, jts, jte, kts, kte
4387 CHARACTER(LEN=1) , INTENT(IN ) :: name
4389 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rfield
4391 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rr
4393 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: field
4395 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, fzp
4397 ! Local data
4399 INTEGER :: i, j, k, itf, jtf, ktf
4401 !<DESCRIPTION>
4402 !
4403 ! decouple decouples a variable from the column dry-air mass.
4404 !
4405 !</DESCRIPTION>
4407 ktf=MIN(kte,kde-1)
4409 IF (name .EQ. 'u')THEN
4410 itf=ite
4411 jtf=MIN(jte,jde-1)
4413 DO j=jts,jtf
4414 DO k=kts,ktf
4415 DO i=its,itf
4416 field(i,k,j)=rfield(i,k,j)/(0.5*(rr(i,k,j)+rr(i-1,k,j)))
4417 ENDDO
4418 ENDDO
4419 ENDDO
4421 ELSE IF (name .EQ. 'v')THEN
4422 itf=MIN(ite,ide-1)
4423 jtf=jte
4425 DO j=jts,jtf
4426 DO k=kts,ktf
4427 DO i=its,itf
4428 field(i,k,j)=rfield(i,k,j)/(0.5*(rr(i,k,j)+rr(i,k,j-1)))
4429 ENDDO
4430 ENDDO
4431 ENDDO
4433 ELSE IF (name .EQ. 'w')THEN
4434 itf=MIN(ite,ide-1)
4435 jtf=MIN(jte,jde-1)
4436 DO j=jts,jtf
4437 DO k=kts+1,ktf
4438 DO i=its,itf
4439 field(i,k,j)=rfield(i,k,j)/(fzm(k)*rr(i,k,j)+fzp(k)*rr(i,k-1,j))
4440 ENDDO
4441 ENDDO
4442 ENDDO
4444 DO j=jts,jtf
4445 DO i=its,itf
4446 field(i,kte,j) = 0.
4447 ENDDO
4448 ENDDO
4450 ELSE
4451 itf=MIN(ite,ide-1)
4452 jtf=MIN(jte,jde-1)
4453 ! For theta we will decouple tb and tp and add them to give t afterwards
4454 DO j=jts,jtf
4455 DO k=kts,ktf
4456 DO i=its,itf
4457 field(i,k,j)=rfield(i,k,j)/rr(i,k,j)
4458 ENDDO
4459 ENDDO
4460 ENDDO
4462 ENDIF
4464 END SUBROUTINE decouple
4466 !-------------------------------------------------------------------------------
4469 SUBROUTINE zero_tend ( tendency, &
4470 ids, ide, jds, jde, kds, kde, &
4471 ims, ime, jms, jme, kms, kme, &
4472 its, ite, jts, jte, kts, kte )
4475 IMPLICIT NONE
4477 ! Input data
4479 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4480 ims, ime, jms, jme, kms, kme, &
4481 its, ite, jts, jte, kts, kte
4483 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
4485 ! Local data
4487 INTEGER :: i, j, k, itf, jtf, ktf
4489 !<DESCRIPTION>
4490 !
4491 ! zero_tend sets the input tendency array to zero.
4492 !
4493 !</DESCRIPTION>
4495 DO j = jts, jte
4496 DO k = kts, kte
4497 DO i = its, ite
4498 tendency(i,k,j) = 0.
4499 ENDDO
4500 ENDDO
4501 ENDDO
4503 END SUBROUTINE zero_tend
4505 !-------------------------------------------------------------------------------
4506 ! Sets the an array on the polar v point(s) to zero
4507 SUBROUTINE zero_pole ( field, &
4508 ids, ide, jds, jde, kds, kde, &
4509 ims, ime, jms, jme, kms, kme, &
4510 its, ite, jts, jte, kts, kte )
4513 IMPLICIT NONE
4515 ! Input data
4517 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4518 ims, ime, jms, jme, kms, kme, &
4519 its, ite, jts, jte, kts, kte
4521 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: field
4523 ! Local data
4525 INTEGER :: i, k
4527 IF (jts == jds) THEN
4528 DO k = kts, kte
4529 DO i = its-1, ite+1
4530 field(i,k,jts) = 0.
4531 END DO
4532 END DO
4533 END IF
4534 IF (jte == jde) THEN
4535 DO k = kts, kte
4536 DO i = its-1, ite+1
4537 field(i,k,jte) = 0.
4538 END DO
4539 END DO
4540 END IF
4542 END SUBROUTINE zero_pole
4544 !-------------------------------------------------------------------------------
4545 ! Sets the an array on the polar v point(s)
4546 SUBROUTINE pole_point_bc ( field, &
4547 ids, ide, jds, jde, kds, kde, &
4548 ims, ime, jms, jme, kms, kme, &
4549 its, ite, jts, jte, kts, kte )
4552 IMPLICIT NONE
4554 ! Input data
4556 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4557 ims, ime, jms, jme, kms, kme, &
4558 its, ite, jts, jte, kts, kte
4560 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: field
4562 ! Local data
4564 INTEGER :: i, k
4566 IF (jts == jds) THEN
4567 DO k = kts, kte
4568 DO i = its, ite
4569 ! field(i,k,jts) = 2*field(i,k,jts+1) - field(i,k,jts+2)
4570 field(i,k,jts) = field(i,k,jts+1)
4571 END DO
4572 END DO
4573 END IF
4574 IF (jte == jde) THEN
4575 DO k = kts, kte
4576 DO i = its, ite
4577 ! field(i,k,jte) = 2*field(i,k,jte-1) - field(i,k,jte-2)
4578 field(i,k,jte) = field(i,k,jte-1)
4579 END DO
4580 END DO
4581 END IF
4583 END SUBROUTINE pole_point_bc
4585 !======================================================================
4586 ! physics prep routines
4587 !======================================================================
4589 SUBROUTINE phy_prep ( config_flags, & ! input
4590 mu, muu, muv, u, v, p, pb, alt, ph, & ! input
4591 phb, t, tsk, moist, n_moist, & ! input
4592 rho, th_phy, p_phy , pi_phy , & ! output
4593 u_phy, v_phy, p8w, t_phy, t8w, & ! output
4594 z, z_at_w, dz8w, & ! output
4595 p_hyd, p_hyd_w, dnw, & ! output
4596 fzm, fzp, znw, p_top, & ! params
4597 RTHRATEN, &
4598 RTHBLTEN, RUBLTEN, RVBLTEN, &
4599 RQVBLTEN, RQCBLTEN, RQIBLTEN, &
4600 RUCUTEN, RVCUTEN, RTHCUTEN, &
4601 RQVCUTEN, RQCCUTEN, RQRCUTEN, &
4602 RQICUTEN, RQSCUTEN, &
4603 RUSHTEN, RVSHTEN, RTHSHTEN, &
4604 RQVSHTEN, RQCSHTEN, RQRSHTEN, &
4605 RQISHTEN, RQSSHTEN, RQGSHTEN, &
4606 RTHFTEN, RQVFTEN, &
4607 RUNDGDTEN, RVNDGDTEN, RTHNDGDTEN, &
4608 RPHNDGDTEN,RQVNDGDTEN, RMUNDGDTEN, &
4609 ids, ide, jds, jde, kds, kde, &
4610 ims, ime, jms, jme, kms, kme, &
4611 its, ite, jts, jte, kts, kte )
4612 !----------------------------------------------------------------------
4613 IMPLICIT NONE
4614 !----------------------------------------------------------------------
4616 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
4618 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4619 ims, ime, jms, jme, kms, kme, &
4620 its, ite, jts, jte, kts, kte
4621 INTEGER , INTENT(IN ) :: n_moist
4623 REAL, DIMENSION( ims:ime, kms:kme , jms:jme , n_moist ), INTENT(IN) :: moist
4626 REAL , DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: TSK, mu, muu, muv
4628 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4629 INTENT( OUT) :: u_phy, &
4630 v_phy, &
4631 pi_phy, &
4632 p_phy, &
4633 p8w, &
4634 t_phy, &
4635 th_phy, &
4636 t8w, &
4637 rho, &
4638 z, &
4639 dz8w, &
4640 z_at_w
4642 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4643 INTENT( OUT) :: p_hyd, &
4644 p_hyd_w
4646 REAL , INTENT(IN ) :: p_top
4648 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4649 INTENT(IN ) :: pb, &
4650 p, &
4651 u, &
4652 v, &
4653 alt, &
4654 ph, &
4655 phb, &
4656 t
4659 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
4660 fzp
4662 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: znw, &
4663 dnw
4665 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4666 INTENT(INOUT) :: RTHRATEN
4668 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4669 INTENT(INOUT) :: RUCUTEN, &
4670 RVCUTEN, &
4671 RTHCUTEN, &
4672 RQVCUTEN, &
4673 RQCCUTEN, &
4674 RQRCUTEN, &
4675 RQICUTEN, &
4676 RQSCUTEN, &
4677 RUSHTEN, &
4678 RVSHTEN, &
4679 RTHSHTEN, &
4680 RQVSHTEN, &
4681 RQCSHTEN, &
4682 RQRSHTEN, &
4683 RQISHTEN, &
4684 RQSSHTEN, &
4685 RQGSHTEN
4687 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
4688 INTENT(INOUT) :: RUBLTEN, &
4689 RVBLTEN, &
4690 RTHBLTEN, &
4691 RQVBLTEN, &
4692 RQCBLTEN, &
4693 RQIBLTEN
4695 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
4696 INTENT(INOUT) :: RTHFTEN, &
4697 RQVFTEN
4699 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
4700 INTENT(INOUT) :: RUNDGDTEN, &
4701 RVNDGDTEN, &
4702 RTHNDGDTEN, &
4703 RPHNDGDTEN, &
4704 RQVNDGDTEN
4706 REAL, DIMENSION( ims:ime, jms:jme ) , &
4707 INTENT(INOUT) :: RMUNDGDTEN
4709 INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end, i_startu, j_startv
4710 INTEGER :: i, j, k
4711 REAL :: w1, w2, z0, z1, z2
4712 REAL :: qtot
4713 INTEGER :: n
4715 !-----------------------------------------------------------------------
4717 !<DESCRIPTION>
4718 !
4719 ! phys_prep calculates a number of diagnostic quantities needed by
4720 ! the physics routines. It also decouples the physics tendencies from
4721 ! the column dry-air mass (the physics routines expect to see/update the
4722 ! uncoupled tendencies).
4723 !
4724 !</DESCRIPTION>
4726 ! set up loop bounds for this grid's boundary conditions
4728 i_start = its
4729 i_end = min( ite,ide-1 )
4730 j_start = jts
4731 j_end = min( jte,jde-1 )
4733 k_start = kts
4734 k_end = min( kte, kde-1 )
4736 ! compute thermodynamics and velocities at pressure points (or half levels)
4738 do j = j_start,j_end
4739 do k = k_start, k_end
4740 do i = i_start, i_end
4742 th_phy(i,k,j) = t(i,k,j) + t0
4743 p_phy(i,k,j) = p(i,k,j) + pb(i,k,j)
4744 pi_phy(i,k,j) = (p_phy(i,k,j)/p1000mb)**rcp
4745 t_phy(i,k,j) = th_phy(i,k,j)*pi_phy(i,k,j)
4746 rho(i,k,j) = 1./alt(i,k,j)*(1.+moist(i,k,j,P_QV))
4747 u_phy(i,k,j) = 0.5*(u(i,k,j)+u(i+1,k,j))
4748 v_phy(i,k,j) = 0.5*(v(i,k,j)+v(i,k,j+1))
4750 enddo
4751 enddo
4752 enddo
4754 ! compute z at w points
4756 do j = j_start,j_end
4757 do k = k_start, kte
4758 do i = i_start, i_end
4759 z_at_w(i,k,j) = (phb(i,k,j)+ph(i,k,j))/g
4760 enddo
4761 enddo
4762 enddo
4764 do j = j_start,j_end
4765 do k = k_start, kte-1
4766 do i = i_start, i_end
4767 dz8w(i,k,j) = z_at_w(i,k+1,j)-z_at_w(i,k,j)
4768 enddo
4769 enddo
4770 enddo
4772 do j = j_start,j_end
4773 do i = i_start, i_end
4774 dz8w(i,kte,j) = 0.
4775 enddo
4776 enddo
4778 ! compute z at p points or half levels (average of z at full levels)
4780 do j = j_start,j_end
4781 do k = k_start, k_end
4782 do i = i_start, i_end
4783 z(i,k,j) = 0.5*(z_at_w(i,k,j) +z_at_w(i,k+1,j) )
4784 enddo
4785 enddo
4786 enddo
4788 ! interp t and p to full levels
4790 do j = j_start,j_end
4791 do k = 2, k_end
4792 do i = i_start, i_end
4793 p8w(i,k,j) = fzm(k)*p_phy(i,k,j)+fzp(k)*p_phy(i,k-1,j)
4794 t8w(i,k,j) = fzm(k)*t_phy(i,k,j)+fzp(k)*t_phy(i,k-1,j)
4795 enddo
4796 enddo
4797 enddo
4799 ! extrapolate p and t to surface and top.
4800 ! we'll use an extrapolation in z for now
4802 do j = j_start,j_end
4803 do i = i_start, i_end
4805 ! bottom
4807 z0 = z_at_w(i,1,j)
4808 z1 = z(i,1,j)
4809 z2 = z(i,2,j)
4810 w1 = (z0 - z2)/(z1 - z2)
4811 w2 = 1. - w1
4812 p8w(i,1,j) = w1*p_phy(i,1,j)+w2*p_phy(i,2,j)
4813 t8w(i,1,j) = w1*t_phy(i,1,j)+w2*t_phy(i,2,j)
4815 ! top
4817 z0 = z_at_w(i,kte,j)
4818 z1 = z(i,k_end,j)
4819 z2 = z(i,k_end-1,j)
4820 w1 = (z0 - z2)/(z1 - z2)
4821 w2 = 1. - w1
4823 ! p8w(i,kde,j) = w1*p_phy(i,kde-1,j)+w2*p_phy(i,kde-2,j)
4824 !!! bug fix extrapolate ln(p) so p is positive definite
4825 p8w(i,kde,j) = exp(w1*log(p_phy(i,kde-1,j))+w2*log(p_phy(i,kde-2,j)))
4826 t8w(i,kde,j) = w1*t_phy(i,kde-1,j)+w2*t_phy(i,kde-2,j)
4828 enddo
4829 enddo
4831 ! calculate hydrostatic pressure at both full and half levels
4832 ! first, full level p: assuming dry over model top
4834 do j = j_start,j_end
4835 do i = i_start, i_end
4836 p_hyd_w(i,kte,j) = p_top
4837 enddo
4838 enddo
4840 do j = j_start,j_end
4841 do k = kte-1, k_start, -1
4842 do i = i_start, i_end
4843 qtot = 0.
4844 do n = PARAM_FIRST_SCALAR,n_moist
4845 qtot = qtot + moist(i,k,j,n)
4846 enddo
4847 p_hyd_w(i,k,j) = p_hyd_w(i,k+1,j) - (1.+qtot)*mu(i,j)*dnw(k)
4848 ! 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)
4849 enddo
4850 enddo
4851 enddo
4853 ! now calculate hydrostatic pressure at half levels
4855 do j = j_start,j_end
4856 do k = k_start, k_end
4857 do i = i_start, i_end
4858 p_hyd(i,k,j) = 0.5*(p_hyd_w(i,k,j)+p_hyd_w(i,k+1,j))
4859 enddo
4860 enddo
4861 enddo
4863 ! decouple all physics tendencies
4865 IF (config_flags%ra_lw_physics .gt. 0 .or. config_flags%ra_sw_physics .gt. 0) THEN
4867 DO J=j_start,j_end
4868 DO K=k_start,k_end
4869 DO I=i_start,i_end
4870 RTHRATEN(I,K,J)=RTHRATEN(I,K,J)/mu(I,J)
4871 ENDDO
4872 ENDDO
4873 ENDDO
4875 ENDIF
4877 IF (config_flags%cu_physics .gt. 0) THEN
4879 DO J=j_start,j_end
4880 DO I=i_start,i_end
4881 DO K=k_start,k_end
4882 RUCUTEN(I,K,J) =RUCUTEN(I,K,J)/mu(I,J)
4883 RVCUTEN(I,K,J) =RVCUTEN(I,K,J)/mu(I,J)
4884 RTHCUTEN(I,K,J)=RTHCUTEN(I,K,J)/mu(I,J)
4885 ENDDO
4886 ENDDO
4887 ENDDO
4889 IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
4890 DO J=j_start,j_end
4891 DO I=i_start,i_end
4892 DO K=k_start,k_end
4893 RQVCUTEN(I,K,J)=RQVCUTEN(I,K,J)/mu(I,J)
4894 ENDDO
4895 ENDDO
4896 ENDDO
4897 ENDIF
4899 IF (P_QC .ge. PARAM_FIRST_SCALAR)THEN
4900 DO J=j_start,j_end
4901 DO I=i_start,i_end
4902 DO K=k_start,k_end
4903 RQCCUTEN(I,K,J)=RQCCUTEN(I,K,J)/mu(I,J)
4904 ENDDO
4905 ENDDO
4906 ENDDO
4907 ENDIF
4909 IF (P_QR .ge. PARAM_FIRST_SCALAR)THEN
4910 DO J=j_start,j_end
4911 DO I=i_start,i_end
4912 DO K=k_start,k_end
4913 RQRCUTEN(I,K,J)=RQRCUTEN(I,K,J)/mu(I,J)
4914 ENDDO
4915 ENDDO
4916 ENDDO
4917 ENDIF
4919 IF (P_QI .ge. PARAM_FIRST_SCALAR)THEN
4920 DO J=j_start,j_end
4921 DO I=i_start,i_end
4922 DO K=k_start,k_end
4923 RQICUTEN(I,K,J)=RQICUTEN(I,K,J)/mu(I,J)
4924 ENDDO
4925 ENDDO
4926 ENDDO
4927 ENDIF
4929 IF(P_QS .ge. PARAM_FIRST_SCALAR)THEN
4930 DO J=j_start,j_end
4931 DO I=i_start,i_end
4932 DO K=k_start,k_end
4933 RQSCUTEN(I,K,J)=RQSCUTEN(I,K,J)/mu(I,J)
4934 ENDDO
4935 ENDDO
4936 ENDDO
4937 ENDIF
4939 ENDIF
4941 IF (config_flags%shcu_physics .gt. 0) THEN
4943 DO J=j_start,j_end
4944 DO I=i_start,i_end
4945 DO K=k_start,k_end
4946 RUSHTEN(I,K,J) =RUSHTEN(I,K,J)/mu(I,J)
4947 RVSHTEN(I,K,J) =RVSHTEN(I,K,J)/mu(I,J)
4948 RTHSHTEN(I,K,J)=RTHSHTEN(I,K,J)/mu(I,J)
4949 ENDDO
4950 ENDDO
4951 ENDDO
4953 IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
4954 DO J=j_start,j_end
4955 DO I=i_start,i_end
4956 DO K=k_start,k_end
4957 RQVSHTEN(I,K,J)=RQVSHTEN(I,K,J)/mu(I,J)
4958 ENDDO
4959 ENDDO
4960 ENDDO
4961 ENDIF
4963 IF (P_QC .ge. PARAM_FIRST_SCALAR)THEN
4964 DO J=j_start,j_end
4965 DO I=i_start,i_end
4966 DO K=k_start,k_end
4967 RQCSHTEN(I,K,J)=RQCSHTEN(I,K,J)/mu(I,J)
4968 ENDDO
4969 ENDDO
4970 ENDDO
4971 ENDIF
4973 IF (P_QR .ge. PARAM_FIRST_SCALAR)THEN
4974 DO J=j_start,j_end
4975 DO I=i_start,i_end
4976 DO K=k_start,k_end
4977 RQRSHTEN(I,K,J)=RQRSHTEN(I,K,J)/mu(I,J)
4978 ENDDO
4979 ENDDO
4980 ENDDO
4981 ENDIF
4983 IF (P_QI .ge. PARAM_FIRST_SCALAR)THEN
4984 DO J=j_start,j_end
4985 DO I=i_start,i_end
4986 DO K=k_start,k_end
4987 RQISHTEN(I,K,J)=RQISHTEN(I,K,J)/mu(I,J)
4988 ENDDO
4989 ENDDO
4990 ENDDO
4991 ENDIF
4993 IF(P_QS .ge. PARAM_FIRST_SCALAR)THEN
4994 DO J=j_start,j_end
4995 DO I=i_start,i_end
4996 DO K=k_start,k_end
4997 RQSSHTEN(I,K,J)=RQSSHTEN(I,K,J)/mu(I,J)
4998 ENDDO
4999 ENDDO
5000 ENDDO
5001 ENDIF
5003 IF(P_QG .ge. PARAM_FIRST_SCALAR)THEN
5004 DO J=j_start,j_end
5005 DO I=i_start,i_end
5006 DO K=k_start,k_end
5007 RQGSHTEN(I,K,J)=RQGSHTEN(I,K,J)/mu(I,J)
5008 ENDDO
5009 ENDDO
5010 ENDDO
5011 ENDIF
5013 ENDIF
5015 IF (config_flags%bl_pbl_physics .gt. 0) THEN
5017 DO J=j_start,j_end
5018 DO K=k_start,k_end
5019 DO I=i_start,i_end
5020 RUBLTEN(I,K,J) =RUBLTEN(I,K,J)/mu(I,J)
5021 RVBLTEN(I,K,J) =RVBLTEN(I,K,J)/mu(I,J)
5022 RTHBLTEN(I,K,J)=RTHBLTEN(I,K,J)/mu(I,J)
5023 ENDDO
5024 ENDDO
5025 ENDDO
5027 IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
5028 DO J=j_start,j_end
5029 DO K=k_start,k_end
5030 DO I=i_start,i_end
5031 RQVBLTEN(I,K,J)=RQVBLTEN(I,K,J)/mu(I,J)
5032 ENDDO
5033 ENDDO
5034 ENDDO
5035 ENDIF
5037 IF (P_QC .ge. PARAM_FIRST_SCALAR) THEN
5038 DO J=j_start,j_end
5039 DO K=k_start,k_end
5040 DO I=i_start,i_end
5041 RQCBLTEN(I,K,J)=RQCBLTEN(I,K,J)/mu(I,J)
5042 ENDDO
5043 ENDDO
5044 ENDDO
5045 ENDIF
5047 IF (P_QI .ge. PARAM_FIRST_SCALAR) THEN
5048 DO J=j_start,j_end
5049 DO K=k_start,k_end
5050 DO I=i_start,i_end
5051 RQIBLTEN(I,K,J)=RQIBLTEN(I,K,J)/mu(I,J)
5052 ENDDO
5053 ENDDO
5054 ENDDO
5055 ENDIF
5057 ENDIF
5059 ! decouple advective forcing required by Grell-Devenyi scheme
5061 if(( config_flags%cu_physics == GDSCHEME ) .OR. &
5062 ( config_flags%cu_physics == G3SCHEME ) .OR. &
5063 ( config_flags%cu_physics == KFETASCHEME ) .OR. &
5064 ( config_flags%cu_physics == TIEDTKESCHEME ) ) then ! Tiedtke ZCX&YQW
5066 DO J=j_start,j_end
5067 DO I=i_start,i_end
5068 DO K=k_start,k_end
5069 RTHFTEN(I,K,J)=RTHFTEN(I,K,J)/mu(I,J)
5070 ENDDO
5071 ENDDO
5072 ENDDO
5074 IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
5075 DO J=j_start,j_end
5076 DO I=i_start,i_end
5077 DO K=k_start,k_end
5078 RQVFTEN(I,K,J)=RQVFTEN(I,K,J)/mu(I,J)
5079 ENDDO
5080 ENDDO
5081 ENDDO
5082 ENDIF
5084 END IF
5086 ! fdda
5087 ! note fdda u and v tendencies are staggered, also only interior points have muu/muv,
5088 ! so only decouple those
5090 IF (config_flags%grid_fdda .gt. 0) THEN
5092 i_startu=MAX(its,ids+1)
5093 j_startv=MAX(jts,jds+1)
5095 DO J=j_start,j_end
5096 DO K=k_start,k_end
5097 DO I=i_startu,i_end
5098 RUNDGDTEN(I,K,J) =RUNDGDTEN(I,K,J)/muu(I,J)
5099 ENDDO
5100 ENDDO
5101 ENDDO
5102 DO J=j_startv,j_end
5103 DO K=k_start,k_end
5104 DO I=i_start,i_end
5105 RVNDGDTEN(I,K,J) =RVNDGDTEN(I,K,J)/muv(I,J)
5106 ENDDO
5107 ENDDO
5108 ENDDO
5109 DO J=j_start,j_end
5110 DO K=k_start,k_end
5111 DO I=i_start,i_end
5112 RTHNDGDTEN(I,K,J)=RTHNDGDTEN(I,K,J)/mu(I,J)
5113 ! RMUNDGDTEN(I,J) - no coupling
5114 ENDDO
5115 ENDDO
5116 ENDDO
5118 IF (config_flags%grid_fdda .EQ. 2) THEN
5119 DO J=j_start,j_end
5120 DO K=k_start,kte
5121 DO I=i_start,i_end
5122 RPHNDGDTEN(I,K,J)=RPHNDGDTEN(I,K,J)/mu(I,J)
5123 ENDDO
5124 ENDDO
5125 ENDDO
5127 ELSE IF (config_flags%grid_fdda .EQ. 1) THEN
5128 IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
5129 DO J=j_start,j_end
5130 DO K=k_start,k_end
5131 DO I=i_start,i_end
5132 RQVNDGDTEN(I,K,J)=RQVNDGDTEN(I,K,J)/mu(I,J)
5133 ENDDO
5134 ENDDO
5135 ENDDO
5136 ENDIF
5137 ENDIF
5139 ENDIF
5141 END SUBROUTINE phy_prep
5143 !------------------------------------------------------------
5145 SUBROUTINE moist_physics_prep_em( t_new, t_old, t0, rho, al, alb, &
5146 p, p8w, p0, pb, ph, phb, &
5147 th_phy, pii, pf, &
5148 z, z_at_w, dz8w, &
5149 dt,h_diabatic, &
5150 config_flags,fzm, fzp, &
5151 ids,ide, jds,jde, kds,kde, &
5152 ims,ime, jms,jme, kms,kme, &
5153 its,ite, jts,jte, kts,kte )
5155 IMPLICIT NONE
5157 ! Here we construct full fields
5158 ! needed by the microphysics
5160 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
5162 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
5163 INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
5164 INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
5166 REAL, INTENT(IN ) :: dt
5168 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5169 INTENT(IN ) :: al, &
5170 alb, &
5171 p, &
5172 pb, &
5173 ph, &
5174 phb
5177 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
5178 fzp
5180 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5181 INTENT( OUT) :: rho, &
5182 th_phy, &
5183 pii, &
5184 pf, &
5185 z, &
5186 z_at_w, &
5187 dz8w, &
5188 p8w
5190 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5191 INTENT(INOUT) :: h_diabatic
5193 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5194 INTENT(INOUT) :: t_new, &
5195 t_old
5197 REAL, INTENT(IN ) :: t0, p0
5198 REAL :: z0,z1,z2,w1,w2
5200 INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
5201 INTEGER :: i, j, k
5203 !--------------------------------------------------------------------
5205 !<DESCRIPTION>
5206 !
5207 ! moist_phys_prep_em calculates a number of diagnostic quantities needed by
5208 ! the microphysics routines.
5209 !
5210 !</DESCRIPTION>
5212 ! set up loop bounds for this grid's boundary conditions
5214 i_start = its
5215 i_end = min( ite,ide-1 )
5216 j_start = jts
5217 j_end = min( jte,jde-1 )
5219 k_start = kts
5220 k_end = min( kte, kde-1 )
5222 DO j = j_start, j_end
5223 DO k = k_start, kte
5224 DO i = i_start, i_end
5225 z_at_w(i,k,j) = (ph(i,k,j)+phb(i,k,j))/g
5226 ENDDO
5227 ENDDO
5228 ENDDO
5230 do j = j_start,j_end
5231 do k = k_start, kte-1
5232 do i = i_start, i_end
5233 dz8w(i,k,j) = z_at_w(i,k+1,j)-z_at_w(i,k,j)
5234 enddo
5235 enddo
5236 enddo
5238 do j = j_start,j_end
5239 do i = i_start, i_end
5240 dz8w(i,kte,j) = 0.
5241 enddo
5242 enddo
5245 ! compute full pii, rho, and z at the new time-level
5246 ! (needed for physics).
5247 ! convert perturbation theta to full theta (th_phy)
5248 ! use h_diabatic to temporarily save pre-microphysics full theta
5250 DO j = j_start, j_end
5251 DO k = k_start, k_end
5252 DO i = i_start, i_end
5254 #ifdef REVERT
5255 t_new(i,k,j) = t_new(i,k,j)-h_diabatic(i,k,j)*dt
5256 #endif
5257 th_phy(i,k,j) = t_new(i,k,j) + t0
5258 h_diabatic(i,k,j) = th_phy(i,k,j)
5259 rho(i,k,j) = 1./(al(i,k,j)+alb(i,k,j))
5260 pii(i,k,j) = ((p(i,k,j)+pb(i,k,j))/p0)**rcp
5261 z(i,k,j) = 0.5*(z_at_w(i,k,j) +z_at_w(i,k+1,j) )
5262 pf(i,k,j) = p(i,k,j)+pb(i,k,j)
5264 ENDDO
5265 ENDDO
5266 ENDDO
5268 ! interp t and p at w points
5270 do j = j_start,j_end
5271 do k = 2, k_end
5272 do i = i_start, i_end
5273 p8w(i,k,j) = fzm(k)*pf(i,k,j)+fzp(k)*pf(i,k-1,j)
5274 enddo
5275 enddo
5276 enddo
5278 ! extrapolate p and t to surface and top.
5279 ! we'll use an extrapolation in z for now
5281 do j = j_start,j_end
5282 do i = i_start, i_end
5284 ! bottom
5286 z0 = z_at_w(i,1,j)
5287 z1 = z(i,1,j)
5288 z2 = z(i,2,j)
5289 w1 = (z0 - z2)/(z1 - z2)
5290 w2 = 1. - w1
5291 p8w(i,1,j) = w1*pf(i,1,j)+w2*pf(i,2,j)
5293 ! top
5295 z0 = z_at_w(i,kte,j)
5296 z1 = z(i,k_end,j)
5297 z2 = z(i,k_end-1,j)
5298 w1 = (z0 - z2)/(z1 - z2)
5299 w2 = 1. - w1
5300 ! p8w(i,kde,j) = w1*pf(i,kde-1,j)+w2*pf(i,kde-2,j)
5301 p8w(i,kde,j) = exp(w1*log(pf(i,kde-1,j))+w2*log(pf(i,kde-2,j)))
5303 enddo
5304 enddo
5306 END SUBROUTINE moist_physics_prep_em
5308 !------------------------------------------------------------------------------
5310 SUBROUTINE moist_physics_finish_em( t_new, t_old, t0, mut, &
5311 th_phy, h_diabatic, dt, &
5312 config_flags, &
5313 #if ( WRF_DFI_RADAR == 1 )
5314 dfi_tten_rad,dfi_stage, &
5315 #endif
5316 ids,ide, jds,jde, kds,kde, &
5317 ims,ime, jms,jme, kms,kme, &
5318 its,ite, jts,jte, kts,kte )
5320 IMPLICIT NONE
5322 ! Here we construct full fields
5323 ! needed by the microphysics
5325 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
5327 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
5328 INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
5329 INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
5331 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5332 INTENT(INOUT) :: t_new, &
5333 t_old, &
5334 th_phy, &
5335 h_diabatic
5336 #if ( WRF_DFI_RADAR == 1 )
5337 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5338 INTENT(IN), OPTIONAL :: dfi_tten_rad
5339 INTEGER, INTENT(IN ) ,OPTIONAL :: dfi_stage
5340 REAL :: dfi_tten_max, old_max
5341 #endif
5343 REAL mpten, mptenmax, mptenmin
5345 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: mut
5348 REAL, INTENT(IN ) :: t0, dt
5350 INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
5351 INTEGER :: i, j, k, imax, jmax, imin, jmin
5353 !--------------------------------------------------------------------
5355 !<DESCRIPTION>
5356 !
5357 ! moist_phys_finish_em resets theta to its perturbation value and
5358 ! computes and stores the microphysics diabatic heating term.
5359 !
5360 !</DESCRIPTION>
5362 ! set up loop bounds for this grid's boundary conditions
5365 i_start = its
5366 i_end = min( ite,ide-1 )
5367 j_start = jts
5368 j_end = min( jte,jde-1 )
5369 ! i_start=max(its,ids+4)
5370 ! i_end=min(ite,ide-5)
5371 ! j_start=max(jts,jds+4)
5372 ! j_end=min(jte,jde-5)
5374 k_start = kts
5375 k_end = min( kte, kde-1 )
5377 #if ( WRF_DFI_RADAR == 1 )
5378 IF ( PRESENT(dfi_stage) .and. PRESENT(dfi_tten_rad) ) THEN
5379 IF ( dfi_stage ==DFI_FWD ) THEN
5380 WRITE(wrf_err_message,*)'Add radar tendency: i_start,j_start: ', i_start, j_start
5381 CALL wrf_debug ( 100 , TRIM(wrf_err_message) )
5382 ENDIF
5383 ENDIF
5384 dfi_tten_max=-999
5385 old_max=-999
5386 #endif
5388 ! add microphysics theta diff to perturbation theta, set h_diabatic
5390 IF ( config_flags%no_mp_heating .eq. 0 ) THEN
5391 mptenmax = 0.
5392 mptenmin = 999.
5393 DO j = j_start, j_end
5394 DO k = k_start, k_end
5395 DO i = i_start, i_end
5396 mpten = th_phy(i,k,j)-h_diabatic(i,k,j)
5397 #if ( WRF_DFI_RADAR == 1 )
5398 if(mpten.gt.mptenmax) then
5399 mptenmax=mpten
5400 imax=i
5401 jmax=j
5402 endif
5403 if(mpten.lt.mptenmin) then
5404 mptenmin=mpten
5405 imin=i
5406 jmin=j
5407 endif
5408 mpten=min(config_flags%mp_tend_lim*dt, mpten)
5409 mpten=max(-config_flags%mp_tend_lim*dt, mpten)
5411 if(k < k_end ) then
5412 if(dfi_tten_max < dfi_tten_rad(i,k,j) ) dfi_tten_max = dfi_tten_rad(i,k,j)
5413 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)
5414 endif
5416 IF ( PRESENT(dfi_stage) .and. PRESENT(dfi_tten_rad) ) THEN
5417 IF ( dfi_stage == DFI_FWD .and. dfi_tten_rad(i,k,j) >= -0.1 .and. &
5418 dfi_tten_rad(i,k,j) <= 0.1 .and. k < k_end ) THEN
5419 ! add radar temp tendency
5420 ! there is radar coverage
5421 t_new(i,k,j) = t_new(i,k,j) + (dfi_tten_rad(i,k,j))*dt
5422 ELSE
5423 ! no radar coverage
5424 t_new(i,k,j) = t_new(i,k,j) + mpten
5425 ENDIF
5426 ENDIF
5427 #else
5428 t_new(i,k,j) = t_new(i,k,j) + mpten
5429 #endif
5430 h_diabatic(i,k,j) = mpten/dt
5431 ENDDO
5432 ENDDO
5433 ENDDO
5435 ELSE
5437 DO j = j_start, j_end
5438 DO k = k_start, k_end
5439 DO i = i_start, i_end
5440 ! t_new(i,k,j) = t_new(i,k,j)
5441 h_diabatic(i,k,j) = 0.
5442 ENDDO
5443 ENDDO
5444 ENDDO
5445 ENDIF
5447 END SUBROUTINE moist_physics_finish_em
5449 !----------------------------------------------------------------
5452 SUBROUTINE init_module_big_step
5453 END SUBROUTINE init_module_big_step
5455 SUBROUTINE set_tend ( field, field_adv_tend, msf, &
5456 ids, ide, jds, jde, kds, kde, &
5457 ims, ime, jms, jme, kms, kme, &
5458 its, ite, jts, jte, kts, kte )
5460 IMPLICIT NONE
5462 ! Input data
5464 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
5465 ims, ime, jms, jme, kms, kme, &
5466 its, ite, jts, jte, kts, kte
5468 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: field
5470 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN) :: field_adv_tend
5472 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: msf
5474 ! Local data
5476 INTEGER :: i, j, k, itf, jtf, ktf
5478 !<DESCRIPTION>
5479 !
5480 ! set_tend copies the advective tendency array into the tendency array.
5481 !
5482 !</DESCRIPTION>
5484 jtf = MIN(jte,jde-1)
5485 ktf = MIN(kte,kde-1)
5486 itf = MIN(ite,ide-1)
5487 DO j = jts, jtf
5488 DO k = kts, ktf
5489 DO i = its, itf
5490 field(i,k,j) = field_adv_tend(i,k,j)*msf(i,j)
5491 ENDDO
5492 ENDDO
5493 ENDDO
5495 END SUBROUTINE set_tend
5497 !------------------------------------------------------------------------------
5499 SUBROUTINE rk_rayleigh_damp( ru_tendf, rv_tendf, &
5500 rw_tendf, t_tendf, &
5501 u, v, w, t, t_init, &
5502 mut, muu, muv, ph, phb, &
5503 u_base, v_base, t_base, z_base, &
5504 dampcoef, zdamp, &
5505 ids, ide, jds, jde, kds, kde, &
5506 ims, ime, jms, jme, kms, kme, &
5507 its, ite, jts, jte, kts, kte )
5509 ! History: Apr 2005 Modifications by George Bryan, NCAR:
5510 ! - Generalized the code in a way that allows for
5511 ! simulations with steep terrain.
5512 !
5513 ! Jul 2004 Modifications by George Bryan, NCAR:
5514 ! - Modified the code to use u_base, v_base, and t_base
5515 ! arrays for the background state. Removed the hard-wired
5516 ! base-state values.
5517 ! - Modified the code to use dampcoef, zdamp, and damp_opt,
5518 ! i.e., the upper-level damper variables in namelist.input.
5519 ! Removed the hard-wired variables in the older version.
5520 ! This damper is used when damp_opt = 2.
5521 ! - Modified the code to account for the movement of the
5522 ! model surfaces with time. The code now obtains a base-
5523 ! state value by interpolation using the "_base" arrays.
5525 ! Nov 2003 Bug fix by Jason Knievel, NCAR
5527 ! Aug 2003 Meridional dimension, some comments, and
5528 ! changes in layout of the code added by
5529 ! Jason Knievel, NCAR
5531 ! Jul 2003 Original code by Bill Skamarock, NCAR
5533 ! Purpose: This routine applies Rayleigh damping to a layer at top
5534 ! of the model domain.
5536 !-----------------------------------------------------------------------
5537 ! Begin declarations.
5539 IMPLICIT NONE
5541 INTEGER, INTENT( IN ) &
5542 :: ids, ide, jds, jde, kds, kde, &
5543 ims, ime, jms, jme, kms, kme, &
5544 its, ite, jts, jte, kts, kte
5546 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
5547 :: ru_tendf, rv_tendf, rw_tendf, t_tendf
5549 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
5550 :: u, v, w, t, t_init, ph, phb
5552 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
5553 :: mut, muu, muv
5555 REAL, DIMENSION( kms:kme ) , INTENT(IN ) &
5556 :: u_base, v_base, t_base, z_base
5558 REAL, INTENT(IN ) &
5559 :: dampcoef, zdamp
5561 ! Local variables.
5563 INTEGER &
5564 :: i_start, i_end, j_start, j_end, k_start, k_end, i, j, k, ktf, k1, k2
5566 REAL &
5567 :: pii, dcoef, z, ztop
5569 REAL :: wkp1, wk, wkm1
5571 REAL, DIMENSION( kms:kme ) :: z00, u00, v00, t00
5573 ! End declarations.
5574 !-----------------------------------------------------------------------
5576 pii = 2.0 * asin(1.0)
5578 ktf = MIN( kte, kde-1 )
5580 !-----------------------------------------------------------------------
5581 ! Adjust u to base state.
5583 DO j = jts, MIN( jte, jde-1 )
5584 DO i = its, MIN( ite, ide )
5586 ! Get height at top of model
5587 ztop = 0.5*( phb(i ,kde,j)+phb(i-1,kde,j) &
5588 +ph(i ,kde,j)+ph(i-1,kde,j) )/g
5590 ! Find bottom of damping layer
5591 k1 = ktf
5592 z = ztop
5593 DO WHILE( z >= (ztop-zdamp) )
5594 z = 0.25*( phb(i ,k1,j)+phb(i ,k1+1,j) &
5595 +phb(i-1,k1,j)+phb(i-1,k1+1,j) &
5596 +ph(i ,k1,j)+ph(i ,k1+1,j) &
5597 +ph(i-1,k1,j)+ph(i-1,k1+1,j))/g
5598 z00(k1) = z
5599 k1 = k1 - 1
5600 ENDDO
5601 k1 = k1 + 2
5603 ! Get reference state at model levels
5604 DO k = k1, ktf
5605 k2 = ktf
5606 DO WHILE( z_base(k2) .gt. z00(k) )
5607 k2 = k2 - 1
5608 ENDDO
5609 if(k2+1.gt.ktf)then
5610 u00(k) = u_base(k2) + ( u_base(k2) - u_base(k2-1) ) &
5611 * ( z00(k) - z_base(k2) ) &
5612 / ( z_base(k2) - z_base(k2-1) )
5613 else
5614 u00(k) = u_base(k2) + ( u_base(k2+1) - u_base(k2) ) &
5615 * ( z00(k) - z_base(k2) ) &
5616 / ( z_base(k2+1) - z_base(k2) )
5617 endif
5618 ENDDO
5620 ! Apply the Rayleigh damper
5621 DO k = k1, ktf
5622 dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
5623 dcoef = (SIN( 0.5 * pii * dcoef ) )**2
5624 ru_tendf(i,k,j) = ru_tendf(i,k,j) - &
5625 muu(i,j) * ( dcoef * dampcoef ) * &
5626 ( u(i,k,j) - u00(k) )
5627 END DO
5629 END DO
5630 END DO
5632 ! End adjustment of u.
5633 !-----------------------------------------------------------------------
5635 !-----------------------------------------------------------------------
5636 ! Adjust v to base state.
5638 DO j = jts, MIN( jte, jde )
5639 DO i = its, MIN( ite, ide-1 )
5641 ! Get height at top of model
5642 ztop = 0.5*( phb(i,kde,j )+phb(i,kde,j-1) &
5643 +ph(i,kde,j )+ph(i,kde,j-1) )/g
5645 ! Find bottom of damping layer
5646 k1 = ktf
5647 z = ztop
5648 DO WHILE( z >= (ztop-zdamp) )
5649 z = 0.25*( phb(i,k1,j )+phb(i,k1+1,j ) &
5650 +phb(i,k1,j-1)+phb(i,k1+1,j-1) &
5651 +ph(i,k1,j )+ph(i,k1+1,j ) &
5652 +ph(i,k1,j-1)+ph(i,k1+1,j-1))/g
5653 z00(k1) = z
5654 k1 = k1 - 1
5655 ENDDO
5656 k1 = k1 + 2
5658 ! Get reference state at model levels
5659 DO k = k1, ktf
5660 k2 = ktf
5661 DO WHILE( z_base(k2) .gt. z00(k) )
5662 k2 = k2 - 1
5663 ENDDO
5664 if(k2+1.gt.ktf)then
5665 v00(k) = v_base(k2) + ( v_base(k2) - v_base(k2-1) ) &
5666 * ( z00(k) - z_base(k2) ) &
5667 / ( z_base(k2) - z_base(k2-1) )
5668 else
5669 v00(k) = v_base(k2) + ( v_base(k2+1) - v_base(k2) ) &
5670 * ( z00(k) - z_base(k2) ) &
5671 / ( z_base(k2+1) - z_base(k2) )
5672 endif
5673 ENDDO
5675 ! Apply the Rayleigh damper
5676 DO k = k1, ktf
5677 dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
5678 dcoef = (SIN( 0.5 * pii * dcoef ) )**2
5679 rv_tendf(i,k,j) = rv_tendf(i,k,j) - &
5680 muv(i,j) * ( dcoef * dampcoef ) * &
5681 ( v(i,k,j) - v00(k) )
5682 END DO
5684 END DO
5685 END DO
5687 ! End adjustment of v.
5688 !-----------------------------------------------------------------------
5690 !-----------------------------------------------------------------------
5691 ! Adjust w to base state.
5693 DO j = jts, MIN( jte, jde-1 )
5694 DO i = its, MIN( ite, ide-1 )
5695 ztop = ( phb(i,kde,j) + ph(i,kde,j) ) / g
5696 DO k = kts, MIN( kte, kde )
5697 z = ( phb(i,k,j) + ph(i,k,j) ) / g
5698 IF ( z >= (ztop-zdamp) ) THEN
5699 dcoef = 1.0 - MIN( 1.0, ( ztop - z ) / zdamp )
5700 dcoef = ( SIN( 0.5 * pii * dcoef ) )**2
5701 rw_tendf(i,k,j) = rw_tendf(i,k,j) - &
5702 mut(i,j) * ( dcoef * dampcoef ) * w(i,k,j)
5703 END IF
5704 END DO
5705 END DO
5706 END DO
5708 ! End adjustment of w.
5709 !-----------------------------------------------------------------------
5711 !-----------------------------------------------------------------------
5712 ! Adjust potential temperature to base state.
5714 DO j = jts, MIN( jte, jde-1 )
5715 DO i = its, MIN( ite, ide-1 )
5717 ! Get height at top of model
5718 ztop = ( phb(i,kde,j) + ph(i,kde,j) ) / g
5720 ! Find bottom of damping layer
5721 k1 = ktf
5722 z = ztop
5723 DO WHILE( z >= (ztop-zdamp) )
5724 z = 0.5 * ( phb(i,k1,j) + phb(i,k1+1,j) + &
5725 ph(i,k1,j) + ph(i,k1+1,j) ) / g
5726 z00(k1) = z
5727 k1 = k1 - 1
5728 ENDDO
5729 k1 = k1 + 2
5731 ! Get reference state at model levels
5732 DO k = k1, ktf
5733 k2 = ktf
5734 DO WHILE( z_base(k2) .gt. z00(k) )
5735 k2 = k2 - 1
5736 ENDDO
5737 if(k2+1.gt.ktf)then
5738 t00(k) = t_base(k2) + ( t_base(k2) - t_base(k2-1) ) &
5739 * ( z00(k) - z_base(k2) ) &
5740 / ( z_base(k2) - z_base(k2-1) )
5741 else
5742 t00(k) = t_base(k2) + ( t_base(k2+1) - t_base(k2) ) &
5743 * ( z00(k) - z_base(k2) ) &
5744 / ( z_base(k2+1) - z_base(k2) )
5745 endif
5746 ENDDO
5748 ! Apply the Rayleigh damper
5749 DO k = k1, ktf
5750 dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
5751 dcoef = (SIN( 0.5 * pii * dcoef ) )**2
5752 t_tendf(i,k,j) = t_tendf(i,k,j) - &
5753 mut(i,j) * ( dcoef * dampcoef ) * &
5754 ( t(i,k,j) - t00(k) )
5755 END DO
5757 END DO
5758 END DO
5760 ! End adjustment of potential temperature.
5761 !-----------------------------------------------------------------------
5763 END SUBROUTINE rk_rayleigh_damp
5765 !==============================================================================
5766 !==============================================================================
5768 SUBROUTINE theta_relaxation( t_tendf, t, t_init, &
5769 mut, ph, phb, &
5770 t_base, z_base, &
5771 ids, ide, jds, jde, kds, kde, &
5772 ims, ime, jms, jme, kms, kme, &
5773 its, ite, jts, jte, kts, kte )
5775 ! Purpose: Newtonian relaxation on potential temperature. Serves two
5776 ! purposes: 1) to mimic atmospheric radiation in a simple
5777 ! manner, and 2) to keep the vertical profile of temperature
5778 ! close to the initial (base-state) profile, which is useful
5779 ! for certain idealized applications.
5781 ! Reference: Rotunno and Emanuel, 1987, JAS, p. 546
5783 !-----------------------------------------------------------------------
5784 ! Begin declarations.
5786 IMPLICIT NONE
5788 INTEGER, INTENT( IN ) &
5789 :: ids, ide, jds, jde, kds, kde, &
5790 ims, ime, jms, jme, kms, kme, &
5791 its, ite, jts, jte, kts, kte
5793 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
5794 :: t_tendf
5796 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
5797 :: t, t_init, ph, phb
5799 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
5800 :: mut
5802 REAL, DIMENSION( kms:kme ) , INTENT(IN ) &
5803 :: t_base, z_base
5805 ! Local variables.
5807 INTEGER :: i, j, k, ktf, k2
5808 REAL :: tau_r , rmax , rmin , inv_tau_r , inv_g , rterm
5809 REAL, DIMENSION( kms:kme ) :: z00,t00
5811 ! End declarations.
5812 !-----------------------------------------------------------------------
5814 ! set tau_r to 12 h, following RE87
5815 tau_r = 12.0*3600.0
5817 ! limit rterm to +/- 2 K/day
5818 rmax = 2.0/86400.0
5819 rmin = -rmax
5821 ktf = MIN( kte, kde-1 )
5822 inv_tau_r = 1.0/tau_r
5823 inv_g = 1.0/g
5825 !-----------------------------------------------------------------------
5826 ! Adjust potential temperature to base state.
5828 DO j = jts, MIN( jte, jde-1 )
5829 DO i = its, MIN( ite, ide-1 )
5831 ! Get height of model levels:
5832 DO k = kts, ktf
5833 z00(k) = 0.5 * ( phb(i,k,j) + phb(i,k+1,j) + &
5834 ph(i,k,j) + ph(i,k+1,j) ) * inv_g
5835 ENDDO
5837 ! Get reference state:
5838 DO k = kts, ktf
5839 k2 = ktf
5840 DO WHILE( z_base(k2) .gt. z00(k) .and. k2 .gt. 1 )
5841 k2 = k2 - 1
5842 ENDDO
5843 if(k2+1.gt.ktf)then
5844 t00(k) = t_base(k2) + ( t_base(k2) - t_base(k2-1) ) &
5845 * ( z00(k) - z_base(k2) ) &
5846 / ( z_base(k2) - z_base(k2-1) )
5847 else
5848 t00(k) = t_base(k2) + ( t_base(k2+1) - t_base(k2) ) &
5849 * ( z00(k) - z_base(k2) ) &
5850 / ( z_base(k2+1) - z_base(k2) )
5851 endif
5852 ENDDO
5854 ! Apply the RE87 R term:
5855 DO k = kts, ktf
5856 rterm = -( t(i,k,j) - t00(k) )*inv_tau_r
5857 ! limit rterm:
5858 rterm = min( rterm , rmax )
5859 rterm = max( rterm , rmin )
5860 t_tendf(i,k,j) = t_tendf(i,k,j) + mut(i,j)*rterm
5861 END DO
5863 END DO
5864 END DO
5866 END SUBROUTINE theta_relaxation
5868 !==============================================================================
5869 !==============================================================================
5871 SUBROUTINE sixth_order_diffusion( name, field, tendency, mu, dt, &
5872 config_flags, &
5873 diff_6th_opt, diff_6th_factor, &
5874 ids, ide, jds, jde, kds, kde, &
5875 ims, ime, jms, jme, kms, kme, &
5876 its, ite, jts, jte, kts, kte )
5878 ! History: 14 Nov 2006 Name of variable changed by Jason Knievel
5879 ! 07 Jun 2006 Revised and generalized by Jason Knievel
5880 ! 25 Apr 2005 Original code by Jason Knievel, NCAR
5882 ! Purpose: Apply 6th-order, monotonic (flux-limited), numerical
5883 ! diffusion to 3-d velocity and to scalars.
5885 ! References: Ming Xue (MWR Aug 2000)
5886 ! Durran ("Numerical Methods for Wave Equations..." 1999)
5887 ! George Bryan (personal communication)
5889 !------------------------------------------------------------------------------
5890 ! Begin: Declarations.
5892 IMPLICIT NONE
5894 INTEGER, INTENT(IN) &
5895 :: ids, ide, jds, jde, kds, kde, &
5896 ims, ime, jms, jme, kms, kme, &
5897 its, ite, jts, jte, kts, kte
5899 TYPE(grid_config_rec_type), INTENT(IN) &
5900 :: config_flags
5902 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) &
5903 :: tendency
5905 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) &
5906 :: field
5908 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN) &
5909 :: mu
5911 REAL, INTENT(IN) &
5912 :: dt
5914 REAL, INTENT(IN) &
5915 :: diff_6th_factor
5917 INTEGER, INTENT(IN) &
5918 :: diff_6th_opt
5920 CHARACTER(LEN=1) , INTENT(IN) &
5921 :: name
5923 INTEGER &
5924 :: i, j, k, &
5925 i_start, i_end, &
5926 j_start, j_end, &
5927 k_start, k_end, &
5928 ktf
5930 REAL &
5931 :: dflux_x_p0, dflux_y_p0, &
5932 dflux_x_p1, dflux_y_p1, &
5933 tendency_x, tendency_y, &
5934 mu_avg_p0, mu_avg_p1, &
5935 diff_6th_coef
5937 LOGICAL &
5938 :: specified
5940 ! End: Declarations.
5941 !------------------------------------------------------------------------------
5943 !------------------------------------------------------------------------------
5944 ! Begin: Translate the diffusion factor into a diffusion coefficient. See
5945 ! Durran's text, section 2.4.3, then adjust for sixth-order diffusion (not
5946 ! fourth) and for diffusion in two dimensions (not one). For reference, a
5947 ! factor of 1.0 would mean complete diffusion of a 2dx wave in one time step,
5948 ! although application of the flux limiter reduces somewhat the effects of
5949 ! diffusion for a given coefficient.
5951 diff_6th_coef = diff_6th_factor * 0.015625 / ( 2.0 * dt )
5953 ! End: Translate diffusion factor.
5954 !------------------------------------------------------------------------------
5956 !------------------------------------------------------------------------------
5957 ! Begin: Assign limits of spatial loops depending on variable to be diffused.
5958 ! The halo regions are already filled with values by the time this subroutine
5959 ! is called, which allows the stencil to extend beyond the domains' edges.
5961 ktf = MIN( kte, kde-1 )
5963 IF ( name .EQ. 'u' ) THEN
5965 i_start = its
5966 i_end = ite
5967 j_start = jts
5968 j_end = MIN(jde-1,jte)
5969 k_start = kts
5970 k_end = ktf
5972 ELSE IF ( name .EQ. 'v' ) THEN
5974 i_start = its
5975 i_end = MIN(ide-1,ite)
5976 j_start = jts
5977 j_end = jte
5978 k_start = kts
5979 k_end = ktf
5981 ELSE IF ( name .EQ. 'w' ) THEN
5983 i_start = its
5984 i_end = MIN(ide-1,ite)
5985 j_start = jts
5986 j_end = MIN(jde-1,jte)
5987 k_start = kts+1
5988 k_end = ktf
5990 ELSE
5992 i_start = its
5993 i_end = MIN(ide-1,ite)
5994 j_start = jts
5995 j_end = MIN(jde-1,jte)
5996 k_start = kts
5997 k_end = ktf
5999 ENDIF
6001 ! End: Assignment of limits of spatial loops.
6002 !------------------------------------------------------------------------------
6004 !------------------------------------------------------------------------------
6005 ! Begin: Loop across spatial dimensions.
6007 DO j = j_start, j_end
6008 DO k = k_start, k_end
6009 DO i = i_start, i_end
6011 !------------------------------------------------------------------------------
6012 ! Begin: Diffusion in x (i index).
6014 ! Calculate the diffusive flux in x direction (from Xue's eq. 3).
6016 dflux_x_p0 = ( 10.0 * ( field(i, k,j) - field(i-1,k,j) ) &
6017 - 5.0 * ( field(i+1,k,j) - field(i-2,k,j) ) &
6018 + ( field(i+2,k,j) - field(i-3,k,j) ) )
6020 dflux_x_p1 = ( 10.0 * ( field(i+1,k,j) - field(i ,k,j) ) &
6021 - 5.0 * ( field(i+2,k,j) - field(i-1,k,j) ) &
6022 + ( field(i+3,k,j) - field(i-2,k,j) ) )
6024 ! If requested in the namelist (diff_6th_opt=2), prohibit up-gradient diffusion
6025 ! (variation on Xue's eq. 10).
6027 IF ( diff_6th_opt .EQ. 2 ) THEN
6029 IF ( dflux_x_p0 * ( field(i ,k,j)-field(i-1,k,j) ) .LE. 0.0 ) THEN
6030 dflux_x_p0 = 0.0
6031 END IF
6033 IF ( dflux_x_p1 * ( field(i+1,k,j)-field(i ,k,j) ) .LE. 0.0 ) THEN
6034 dflux_x_p1 = 0.0
6035 END IF
6037 END IF
6039 ! Apply 6th-order diffusion in x direction.
6041 IF ( name .EQ. 'u' ) THEN
6042 mu_avg_p0 = mu(i-1,j)
6043 mu_avg_p1 = mu(i ,j)
6044 ELSE IF ( name .EQ. 'v' ) THEN
6045 mu_avg_p0 = 0.25 * ( &
6046 mu(i-1,j-1) + &
6047 mu(i ,j-1) + &
6048 mu(i-1,j ) + &
6049 mu(i ,j ) )
6050 mu_avg_p1 = 0.25 * ( &
6051 mu(i ,j-1) + &
6052 mu(i+1,j-1) + &
6053 mu(i ,j ) + &
6054 mu(i+1,j ) )
6055 ELSE
6056 mu_avg_p0 = 0.5 * ( &
6057 mu(i-1,j) + &
6058 mu(i ,j) )
6059 mu_avg_p1 = 0.5 * ( &
6060 mu(i ,j) + &
6061 mu(i+1,j) )
6062 END IF
6064 tendency_x = diff_6th_coef * &
6065 ( ( mu_avg_p1 * dflux_x_p1 ) - ( mu_avg_p0 * dflux_x_p0 ) )
6067 ! End: Diffusion in x.
6068 !------------------------------------------------------------------------------
6070 !------------------------------------------------------------------------------
6071 ! Begin: Diffusion in y (j index).
6073 ! Calculate the diffusive flux in y direction (from Xue's eq. 3).
6075 dflux_y_p0 = ( 10.0 * ( field(i,k,j ) - field(i,k,j-1) ) &
6076 - 5.0 * ( field(i,k,j+1) - field(i,k,j-2) ) &
6077 + ( field(i,k,j+2) - field(i,k,j-3) ) )
6079 dflux_y_p1 = ( 10.0 * ( field(i,k,j+1) - field(i,k,j ) ) &
6080 - 5.0 * ( field(i,k,j+2) - field(i,k,j-1) ) &
6081 + ( field(i,k,j+3) - field(i,k,j-2) ) )
6083 ! If requested in the namelist (diff_6th_opt=2), prohibit up-gradient diffusion
6084 ! (variation on Xue's eq. 10).
6086 IF ( diff_6th_opt .EQ. 2 ) THEN
6088 IF ( dflux_y_p0 * ( field(i,k,j )-field(i,k,j-1) ) .LE. 0.0 ) THEN
6089 dflux_y_p0 = 0.0
6090 END IF
6092 IF ( dflux_y_p1 * ( field(i,k,j+1)-field(i,k,j ) ) .LE. 0.0 ) THEN
6093 dflux_y_p1 = 0.0
6094 END IF
6096 END IF
6098 ! Apply 6th-order diffusion in y direction.
6100 IF ( name .EQ. 'u' ) THEN
6101 mu_avg_p0 = 0.25 * ( &
6102 mu(i-1,j-1) + &
6103 mu(i ,j-1) + &
6104 mu(i-1,j ) + &
6105 mu(i ,j ) )
6106 mu_avg_p1 = 0.25 * ( &
6107 mu(i-1,j ) + &
6108 mu(i ,j ) + &
6109 mu(i-1,j+1) + &
6110 mu(i ,j+1) )
6111 ELSE IF ( name .EQ. 'v' ) THEN
6112 mu_avg_p0 = mu(i,j-1)
6113 mu_avg_p1 = mu(i,j )
6114 ELSE
6115 mu_avg_p0 = 0.5 * ( &
6116 mu(i,j-1) + &
6117 mu(i,j ) )
6118 mu_avg_p1 = 0.5 * ( &
6119 mu(i,j ) + &
6120 mu(i,j+1) )
6121 END IF
6123 tendency_y = diff_6th_coef * &
6124 ( ( mu_avg_p1 * dflux_y_p1 ) - ( mu_avg_p0 * dflux_y_p0 ) )
6126 ! End: Diffusion in y.
6127 !------------------------------------------------------------------------------
6129 !------------------------------------------------------------------------------
6130 ! Begin: Combine diffusion in x and y.
6132 tendency(i,k,j) = tendency(i,k,j) + tendency_x + tendency_y
6134 ! End: Combine diffusion in x and y.
6135 !------------------------------------------------------------------------------
6137 ENDDO
6138 ENDDO
6139 ENDDO
6141 ! End: Loop across spatial dimensions.
6142 !------------------------------------------------------------------------------
6144 END SUBROUTINE sixth_order_diffusion
6146 !==============================================================================
6147 !==============================================================================
6149 END MODULE module_big_step_utilities_em