| blob | 523d94337de866dded56d67b7e0a520259f7b111 |
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, &
957 al, alb, mu, muts, ph, p, pb, &
958 t, p0, t0, znu, dnw, rdnw, &
959 rdn, non_hydrostatic, &
960 ids, ide, jds, jde, kds, kde, &
961 ims, ime, jms, jme, kms, kme, &
962 its, ite, jts, jte, kts, kte )
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
976 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: alb, &
977 pb, &
978 t
980 REAL, DIMENSION( ims:ime , kms:kme , jms:jme, n_moist ), INTENT(IN ) :: moist
982 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT( OUT) :: al, p
984 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: ph
986 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN ) :: mu, muts
988 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: znu, dnw, rdnw, rdn
990 REAL, INTENT(IN ) :: t0, p0
992 ! Local stuff
994 INTEGER :: i, j, k, itf, jtf, ktf, ispe
995 REAL :: qvf, qtot, qf1, qf2
996 REAL, DIMENSION( its:ite) :: temp,cpovcv_v
999 !<DESCRIPTION>
1000 !
1001 ! For the nonhydrostatic option, calc_p_rho_phi calculates the
1002 ! diagnostic quantities pressure and (inverse) density from the
1003 ! prognostic variables using the equation of state.
1004 !
1005 ! For the hydrostatic option, calc_p_rho_phi calculates the
1006 ! diagnostic quantities (inverse) density and geopotential from the
1007 ! prognostic variables using the equation of state and the hydrostatic
1008 ! equation.
1009 !
1010 !</DESCRIPTION>
1012 itf=MIN(ite,ide-1)
1013 jtf=MIN(jte,jde-1)
1014 ktf=MIN(kte,kde-1)
1016 #ifndef INTELMKL
1017 cpovcv_v = cpovcv
1018 #endif
1020 IF (non_hydrostatic) THEN
1022 IF (n_moist >= PARAM_FIRST_SCALAR ) THEN
1024 DO j=jts,jtf
1025 DO k=kts,ktf
1026 DO i=its,itf
1027 qvf = 1.+rvovrd*moist(i,k,j,P_QV)
1028 al(i,k,j)=-1./muts(i,j)*(alb(i,k,j)*mu(i,j) &
1029 +rdnw(k)*(ph(i,k+1,j)-ph(i,k,j)))
1030 temp(i)=(r_d*(t0+t(i,k,j))*qvf)/ &
1031 (p0*(al(i,k,j)+alb(i,k,j)))
1032 ENDDO
1033 #ifdef INTELMKL
1034 CALL VPOWX ( itf-its+1, temp(its), cpovcv, p(its,k,j) )
1035 #else
1036 ! use vector version from libmassv or from compat lib in frame/libmassv.F
1037 CALL VPOW ( p(its,k,j), temp(its), cpovcv_v(its), itf-its+1 )
1038 #endif
1039 DO i=its,itf
1040 p(i,k,j)= p(i,k,j)*p0-pb(i,k,j)
1041 ENDDO
1042 ENDDO
1043 ENDDO
1045 ELSE
1047 DO j=jts,jtf
1048 DO k=kts,ktf
1049 DO i=its,itf
1050 al(i,k,j)=-1./muts(i,j)*(alb(i,k,j)*mu(i,j) &
1051 +rdnw(k)*(ph(i,k+1,j)-ph(i,k,j)))
1052 p(i,k,j)=p0*( (r_d*(t0+t(i,k,j)))/ &
1053 (p0*(al(i,k,j)+alb(i,k,j))) )**cpovcv &
1054 -pb(i,k,j)
1055 ENDDO
1056 ENDDO
1057 ENDDO
1059 END IF
1061 ELSE
1063 ! hydrostatic pressure, al, and ph1 calc; WCS, 5 sept 2001
1066 IF (n_moist >= PARAM_FIRST_SCALAR ) THEN
1068 DO j=jts,jtf
1070 k=ktf ! top layer
1071 DO i=its,itf
1073 qtot = 0.
1074 DO ispe=PARAM_FIRST_SCALAR,n_moist
1075 qtot = qtot + moist(i,k,j,ispe)
1076 ENDDO
1077 qf2 = 1./(1.+qtot)
1078 qf1 = qtot*qf2
1080 p(i,k,j) = - 0.5*(mu(i,j)+qf1*muts(i,j))/rdnw(k)/qf2
1081 qvf = 1.+rvovrd*moist(i,k,j,P_QV)
1082 al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* &
1083 (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1085 ENDDO
1087 DO k=ktf-1,kts,-1 ! remaining layers, integrate down
1088 DO i=its,itf
1090 qtot = 0.
1091 DO ispe=PARAM_FIRST_SCALAR,n_moist
1092 qtot = qtot + 0.5*( moist(i,k ,j,ispe) + moist(i,k+1,j,ispe) )
1093 ENDDO
1094 qf2 = 1./(1.+qtot)
1095 qf1 = qtot*qf2
1097 p(i,k,j) = p(i,k+1,j) - (mu(i,j) + qf1*muts(i,j))/qf2/rdn(k+1)
1098 qvf = 1.+rvovrd*moist(i,k,j,P_QV)
1099 al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* &
1100 (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1101 ENDDO
1102 ENDDO
1104 DO k=2,ktf+1 ! integrate hydrostatic equation for geopotential
1105 DO i=its,itf
1107 ! ph(i,k,j) = ph(i,k-1,j) - (1./rdnw(k-1))*( &
1108 ! (muts(i,j)+mu(i,j))*al(i,k-1,j)+ &
1109 ! mu(i,j)*alb(i,k-1,j) )
1110 ph(i,k,j) = ph(i,k-1,j) - (dnw(k-1))*( &
1111 (muts(i,j))*al(i,k-1,j)+ &
1112 mu(i,j)*alb(i,k-1,j) )
1115 ENDDO
1116 ENDDO
1118 ENDDO
1120 ELSE
1122 DO j=jts,jtf
1124 k=ktf ! top layer
1125 DO i=its,itf
1127 qtot = 0.
1128 qf2 = 1./(1.+qtot)
1129 qf1 = qtot*qf2
1131 p(i,k,j) = - 0.5*(mu(i,j)+qf1*muts(i,j))/rdnw(k)/qf2
1132 qvf = 1.
1133 al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* &
1134 (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1136 ENDDO
1138 DO k=ktf-1,kts,-1 ! remaining layers, integrate down
1139 DO i=its,itf
1141 qtot = 0.
1142 qf2 = 1./(1.+qtot)
1143 qf1 = qtot*qf2
1145 p(i,k,j) = p(i,k+1,j) - (mu(i,j) + qf1*muts(i,j))/qf2/rdn(k+1)
1146 qvf = 1.
1147 al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* &
1148 (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1149 ENDDO
1150 ENDDO
1152 DO k=2,ktf+1 ! integrate hydrostatic equation for geopotential
1153 DO i=its,itf
1155 ! ph(i,k,j) = ph(i,k-1,j) - (1./rdnw(k-1))*( &
1156 ! (muts(i,j)+mu(i,j))*al(i,k-1,j)+ &
1157 ! mu(i,j)*alb(i,k-1,j) )
1158 ph(i,k,j) = ph(i,k-1,j) - (dnw(k-1))*( &
1159 (muts(i,j))*al(i,k-1,j)+ &
1160 mu(i,j)*alb(i,k-1,j) )
1163 ENDDO
1164 ENDDO
1166 ENDDO
1168 END IF
1170 END IF
1172 END SUBROUTINE calc_p_rho_phi
1174 !----------------------------------------------------------------------
1176 SUBROUTINE calc_php ( php, ph, phb, &
1177 ids, ide, jds, jde, kds, kde, &
1178 ims, ime, jms, jme, kms, kme, &
1179 its, ite, jts, jte, kts, kte )
1181 IMPLICIT NONE
1183 ! Input data
1185 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1186 ims, ime, jms, jme, kms, kme, &
1187 its, ite, jts, jte, kts, kte
1189 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: phb, ph
1190 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT( OUT) :: php
1192 ! Local stuff
1194 INTEGER :: i, j, k, itf, jtf, ktf
1196 !<DESCRIPTION>
1197 !
1198 ! calc_php calculates the full geopotential from the reference state
1199 ! geopotential and the perturbation geopotential (phb_ph).
1200 !
1201 !</DESCRIPTION>
1203 itf=MIN(ite,ide-1)
1204 jtf=MIN(jte,jde-1)
1205 ktf=MIN(kte,kde-1)
1207 DO j=jts,jtf
1208 DO k=kts,ktf
1209 DO i=its,itf
1210 php(i,k,j) = 0.5*(phb(i,k,j)+phb(i,k+1,j)+ph(i,k,j)+ph(i,k+1,j))
1211 ENDDO
1212 ENDDO
1213 ENDDO
1215 END SUBROUTINE calc_php
1217 !-------------------------------------------------------------------------------
1219 SUBROUTINE diagnose_w( ph_tend, ph_new, ph_old, w, mu, dt, &
1220 u, v, ht, &
1221 cf1, cf2, cf3, rdx, rdy, &
1222 msftx, msfty, &
1223 ids, ide, jds, jde, kds, kde, &
1224 ims, ime, jms, jme, kms, kme, &
1225 its, ite, jts, jte, kts, kte )
1227 IMPLICIT NONE
1229 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1230 ims, ime, jms, jme, kms, kme, &
1231 its, ite, jts, jte, kts, kte
1233 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: ph_tend, &
1234 ph_new, &
1235 ph_old, &
1236 u, &
1237 v
1240 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT( OUT) :: w
1242 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mu, ht, msftx, msfty
1244 REAL, INTENT(IN ) :: dt, cf1, cf2, cf3, rdx, rdy
1246 INTEGER :: i, j, k, itf, jtf
1248 itf=MIN(ite,ide-1)
1249 jtf=MIN(jte,jde-1)
1251 !<DESCRIPTION>
1252 !
1253 ! diagnose_w diagnoses the vertical velocity from the geopoential equation.
1254 ! Used with the hydrostatic option.
1255 !
1256 !</DESCRIPTION>
1258 DO j = jts, jtf
1260 ! lower b.c. on w
1262 ! Notes on map scale factors:
1263 ! Chain rule: if Z=Z(X,Y) [true at the surface] then
1264 ! dZ/dt = dZ/dX * dX/dt + dZ/dY * dY/dt, U=dX/dt, V=dY/dt
1265 ! Using capitals to denote actual values
1266 ! In mapped values, u=U, v=V, z=Z, 1/dX=mx/dx, 1/dY=my/dy
1267 ! => w = dz/dt = mx u dz/dx + my v dz/dy
1268 ! [where dz/dx is just the surface height change between x
1269 ! gridpoints, and dz/dy is the change between y gridpoints]
1270 ! [NB: cf1, cf2 and cf3 do vertical weighting of u or v values
1271 ! nearest the surface]
1273 ! Previously msft multiplied by rdy and rdx terms.
1274 ! Now msfty multiplies rdy term, and msftx multiplies msftx term
1275 DO i = its, itf
1276 w(i,1,j)= msfty(i,j)*.5*rdy*( &
1277 (ht(i,j+1)-ht(i,j )) &
1278 *(cf1*v(i,1,j+1)+cf2*v(i,2,j+1)+cf3*v(i,3,j+1)) &
1279 +(ht(i,j )-ht(i,j-1)) &
1280 *(cf1*v(i,1,j )+cf2*v(i,2,j )+cf3*v(i,3,j )) ) &
1281 +msftx(i,j)*.5*rdx*( &
1282 (ht(i+1,j)-ht(i,j )) &
1283 *(cf1*u(i+1,1,j)+cf2*u(i+1,2,j)+cf3*u(i+1,3,j)) &
1284 +(ht(i,j )-ht(i-1,j)) &
1285 *(cf1*u(i ,1,j)+cf2*u(i ,2,j)+cf3*u(i ,3,j)) )
1286 ENDDO
1288 ! use geopotential equation to diagnose w
1290 ! Further notes on map scale factors
1291 ! If ph_tend contains: -mx partial d/dx(mu rho u/my)
1292 ! -mx partial d/dy(phi mu v/mx)
1293 ! -partial d/dz(phi mu w/my)
1294 ! then phi eqn is: partial d/dt(mu phi/my) = ph_tend + mu g w/my
1295 ! => w = [my/(mu*g)]*[partial d/dt(mu phi/my) - ph_tend]
1297 DO k = 2, kte
1298 DO i = its, itf
1299 w(i,k,j) = msfty(i,j)*( (ph_new(i,k,j)-ph_old(i,k,j))/dt &
1300 - ph_tend(i,k,j)/mu(i,j) )/g
1302 ENDDO
1303 ENDDO
1305 ENDDO
1307 END SUBROUTINE diagnose_w
1309 !-------------------------------------------------------------------------------
1311 SUBROUTINE rhs_ph( ph_tend, u, v, ww, &
1312 ph, ph_old, phb, w, &
1313 mut, muu, muv, &
1314 fnm, fnp, &
1315 rdnw, cfn, cfn1, rdx, rdy, &
1316 msfux, msfuy, msfvx, &
1317 msfvx_inv, msfvy, &
1318 msftx, msfty, &
1319 non_hydrostatic, &
1320 config_flags, &
1321 ids, ide, jds, jde, kds, kde, &
1322 ims, ime, jms, jme, kms, kme, &
1323 its, ite, jts, jte, kts, kte )
1324 IMPLICIT NONE
1326 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
1328 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1329 ims, ime, jms, jme, kms, kme, &
1330 its, ite, jts, jte, kts, kte
1332 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: &
1333 u, &
1334 v, &
1335 ww, &
1336 ph, &
1337 ph_old, &
1338 phb, &
1339 w
1341 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: ph_tend
1343 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: muu, muv, mut, &
1344 msfux, msfuy, &
1345 msfvx, msfvy, &
1346 msftx, msfty, &
1347 msfvx_inv
1349 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, fnm, fnp
1351 REAL, INTENT(IN ) :: cfn, cfn1, rdx, rdy
1353 LOGICAL, INTENT(IN ) :: non_hydrostatic
1355 ! Local stuff
1357 INTEGER :: i, j, k, itf, jtf, ktf, kz, i_start, j_start
1358 REAL :: ur, ul, ub, vr, vl, vb
1359 REAL, DIMENSION(its:ite,kts:kte) :: wdwn
1361 INTEGER :: advective_order
1363 LOGICAL :: specified
1365 !<DESCRIPTION>
1366 !
1367 ! rhs_ph calculates the large-timestep tendency terms for the geopotential
1368 ! equation. These terms include the advection and "gw". The geopotential
1369 ! equation is cast in advective form, so we don't use the flux form advection
1370 ! algorithms here.
1371 !
1372 !</DESCRIPTION>
1374 specified = .false.
1375 if(config_flags%specified .or. config_flags%nested) specified = .true.
1377 advective_order = config_flags%h_sca_adv_order
1378 ! advective_order = 2 ! original configuration (pre Oct 2001)
1380 itf=MIN(ite,ide-1)
1381 jtf=MIN(jte,jde-1)
1382 ktf=MIN(kte,kde-1)
1384 ! Notes on map scale factors (WCS, 2 march 2008)
1385 ! phi equation is: mu/my d/dt(phi) = -(1/my) mx mu u d/dx(phi)
1386 ! -(1/my) my mu v d/dy(phi)
1387 ! - omega d/d_eta(phi)
1388 ! +mu g w/my
1389 !
1390 ! A little further explanation...
1391 ! The tendency term we are computing here is for mu/my d/dt(phi). It is advective form
1392 ! but it is multiplied be mu/my. It will be decoupled from (mu/my) when the implicit w-phi
1393 ! solution is computed in subourine advance_w. The formulation dates from the early
1394 ! days of the mass coordinate model when we were testing both a flux and an advective formulation
1395 ! for the geopotential equation and different forms of the vertical momentum equation and the
1396 ! vertically implicit solver.
1398 ! advective form for the geopotential equation
1400 DO j = jts, jtf
1402 DO k = 2, kte
1403 DO i = its, itf
1404 wdwn(i,k) = .5*(ww(i,k,j)+ww(i,k-1,j))*rdnw(k-1) &
1405 *(ph(i,k,j)-ph(i,k-1,j)+phb(i,k,j)-phb(i,k-1,j))
1406 ENDDO
1407 ENDDO
1409 ! RHS term 3 is: - omega partial d/dnu(phi)
1411 DO k = 2, kte-1
1412 DO i = its, itf
1413 ph_tend(i,k,j) = ph_tend(i,k,j) &
1414 - (fnm(k)*wdwn(i,k+1)+fnp(k)*wdwn(i,k))
1415 ENDDO
1416 ENDDO
1418 ENDDO
1420 IF (non_hydrostatic) THEN ! add in "gw" term.
1421 DO j = jts, jtf ! in hydrostatic mode, "gw" will be diagnosed
1422 ! after the timestep to give us "w"
1423 DO i = its, itf
1424 ph_tend(i,kde,j) = 0.
1425 ENDDO
1427 DO k = 2, kte
1428 DO i = its, itf
1429 ! phi equation RHS term 4
1430 ph_tend(i,k,j) = ph_tend(i,k,j) + mut(i,j)*g*w(i,k,j)/msfty(i,j)
1431 ENDDO
1432 ENDDO
1434 ENDDO
1436 END IF
1438 ! Notes on map scale factors:
1439 ! RHS terms 1 and 2 are: -(1/my) mx u mu partial d/dx(phi)
1440 ! -(1/my) my v mu partial d/dy(phi)
1442 IF (advective_order <= 2) THEN
1444 ! y (v) advection
1446 i_start = its
1447 j_start = jts
1448 itf=MIN(ite,ide-1)
1449 jtf=MIN(jte,jde-1)
1451 IF ( (config_flags%open_ys) .and. jts == jds ) j_start = jts+1
1452 IF ( (config_flags%open_ye) .and. jte == jde ) jtf = jtf-1
1454 DO j = j_start, jtf
1456 DO k = 2, kte-1
1457 DO i = i_start, itf
1458 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* &
1459 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)* &
1460 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1461 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j )* &
1462 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1463 ENDDO
1464 ENDDO
1466 k = kte
1467 DO i = i_start, itf
1468 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* &
1469 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)* &
1470 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1471 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j )* &
1472 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1473 ENDDO
1475 ENDDO
1477 ! x (u) advection
1479 i_start = its
1480 j_start = jts
1481 itf=MIN(ite,ide-1)
1482 jtf=MIN(jte,jde-1)
1484 IF ( (config_flags%open_xs) .and. its == ids ) i_start = its+1
1485 IF ( (config_flags%open_xe) .and. ite == ide ) itf = itf-1
1487 DO j = j_start, jtf
1489 DO k = 2, kte-1
1490 DO i = i_start, itf
1491 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))* &
1492 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)* &
1493 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1494 +muu(i ,j)*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j)* &
1495 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1496 ENDDO
1497 ENDDO
1499 k = kte
1500 DO i = i_start, itf
1501 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))* &
1502 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)* &
1503 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1504 +muu(i ,j)*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux( i,j)* &
1505 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1506 ENDDO
1508 ENDDO
1510 ELSE IF (advective_order <= 4) THEN
1512 ! y (v) 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_ys) .and. jts == jds ) j_start = jts+1
1520 IF ( (config_flags%open_ye) .and. jte == jde ) jtf = jtf-1
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*rdy/msfty(i,j))*( &
1527 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1) &
1528 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j ))* (1./12.)*( &
1529 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1530 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1531 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1532 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1535 ENDDO
1536 ENDDO
1538 k = kte
1539 DO i = i_start, itf
1540 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))*( &
1541 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1) &
1542 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j ))* (1./12.)*( &
1543 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1544 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1545 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1546 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1548 ENDDO
1550 ENDDO
1553 ! x (u) advection
1555 i_start = its
1556 j_start = jts
1557 itf=MIN(ite,ide-1)
1558 jtf=MIN(jte,jde-1)
1560 IF ( (config_flags%open_xs) .and. its == ids ) i_start = its+1
1561 IF ( (config_flags%open_xe) .and. ite == ide ) itf = itf-1
1563 DO j = j_start, jtf
1565 DO k = 2, kte-1
1566 DO i = i_start, itf
1567 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*( &
1568 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j) &
1569 +muu(i,j )*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j) )* (1./12.)*( &
1570 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1571 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1572 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1573 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1574 ENDDO
1575 ENDDO
1577 k = kte
1578 DO i = i_start, itf
1579 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*( &
1580 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j) &
1581 +muu(i,j )*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux(i ,j) )* (1./12.)*( &
1582 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1583 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1584 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1585 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1586 ENDDO
1588 ENDDO
1590 ELSE IF (advective_order <= 6) THEN
1592 ! y (v) advection
1594 i_start = its
1595 j_start = jts
1596 itf=MIN(ite,ide-1)
1597 jtf=MIN(jte,jde-1)
1599 ! IF ( (config_flags%open_ys) .and. jts == jds ) j_start = jts+1
1600 ! IF ( (config_flags%open_ye) .and. jte == jde ) jtf = jtf-1
1602 IF (config_flags%open_ys .or. specified ) j_start = max(jts,jds+2)
1603 IF (config_flags%open_ye .or. specified ) jtf = min(jtf,jde-3)
1605 DO j = j_start, jtf
1607 DO k = 2, kte-1
1608 DO i = i_start, itf
1609 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* ( &
1610 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1) &
1611 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j ) )* (1./60.)*( &
1612 45.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1613 -9.*(ph(i,k,j+2)-ph(i,k,j-2)) &
1614 +(ph(i,k,j+3)-ph(i,k,j-3)) &
1615 +45.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1616 -9.*(phb(i,k,j+2)-phb(i,k,j-2)) &
1617 +(phb(i,k,j+3)-phb(i,k,j-3)) ) )
1620 ENDDO
1621 ENDDO
1623 k = kte
1624 DO i = i_start, itf
1625 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* ( &
1626 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1) &
1627 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j ) )* (1./60.)*( &
1628 45.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1629 -9.*(ph(i,k,j+2)-ph(i,k,j-2)) &
1630 +(ph(i,k,j+3)-ph(i,k,j-3)) &
1631 +45.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1632 -9.*(phb(i,k,j+2)-phb(i,k,j-2)) &
1633 +(phb(i,k,j+3)-phb(i,k,j-3)) ) )
1635 ENDDO
1637 ENDDO
1640 ! pick up near boundary rows using 4th order stencil
1641 ! (open bc copy only goes out to jds-1 and jde, hence 4rth is ok but 6th is too big)
1643 IF ( (config_flags%open_ys) .and. jts <= jds+1 ) THEN
1645 j = jds+1
1646 DO k = 2, kte-1
1647 DO i = i_start, itf
1648 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* ( &
1649 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1) &
1650 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j ) )* (1./12.)*( &
1651 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1652 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1653 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1654 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1657 ENDDO
1658 ENDDO
1660 k = kte
1661 DO i = i_start, itf
1662 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* ( &
1663 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1) &
1664 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j) )* (1./12.)*( &
1665 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1666 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1667 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1668 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1670 ENDDO
1672 END IF
1674 IF ( (config_flags%open_ye) .and. jte >= jde-2 ) THEN
1676 j = jde-2
1677 DO k = 2, kte-1
1678 DO i = i_start, itf
1679 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* ( &
1680 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1) &
1681 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j) )* (1./12.)*( &
1682 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1683 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1684 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1685 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1688 ENDDO
1689 ENDDO
1691 k = kte
1692 DO i = i_start, itf
1693 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* ( &
1694 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1) &
1695 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j) )* (1./12.)*( &
1696 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1697 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1698 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1699 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1701 ENDDO
1703 END IF
1705 ! x (u) advection
1707 i_start = its
1708 j_start = jts
1709 itf=MIN(ite,ide-1)
1710 jtf=MIN(jte,jde-1)
1712 IF (config_flags%open_xs .or. specified ) i_start = max(its,ids+2)
1713 IF (config_flags%open_xe .or. specified ) itf = min(itf,ide-3)
1714 IF ( config_flags%periodic_x ) i_start = its
1715 IF ( config_flags%periodic_x ) itf=MIN(ite,ide-1)
1717 DO j = j_start, jtf
1719 DO k = 2, kte-1
1720 DO i = i_start, itf
1721 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*( &
1722 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j) &
1723 +muu(i,j )*(u(i,k,j )+u(i,k-1,j ))*msfux(i,j) )* (1./60.)*( &
1724 45.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1725 -9.*(ph(i+2,k,j)-ph(i-2,k,j)) &
1726 +(ph(i+3,k,j)-ph(i-3,k,j)) &
1727 +45.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1728 -9.*(phb(i+2,k,j)-phb(i-2,k,j)) &
1729 +(phb(i+3,k,j)-phb(i-3,k,j)) ) )
1730 ENDDO
1731 ENDDO
1733 k = kte
1734 DO i = i_start, itf
1735 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*( &
1736 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j) &
1737 +muu(i,j )*(cfn*u(i ,k-1,j)+cfn1*u(i,k-2,j))*msfux(i,j) )* (1./60.)*( &
1738 45.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1739 -9.*(ph(i+2,k,j)-ph(i-2,k,j)) &
1740 +(ph(i+3,k,j)-ph(i-3,k,j)) &
1741 +45.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1742 -9.*(phb(i+2,k,j)-phb(i-2,k,j)) &
1743 +(phb(i+3,k,j)-phb(i-3,k,j)) ) )
1744 ENDDO
1746 ENDDO
1748 IF ( (config_flags%open_xs) .and. its <= ids+1 ) THEN
1749 i = ids + 1
1750 DO j = j_start, jtf
1751 DO k = 2, kte-1
1752 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*( &
1753 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j) &
1754 +muu(i,j )*(u(i,k,j )+u(i,k-1,j ))*msfux(i,j) )* (1./12.)*( &
1755 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1756 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1757 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1758 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1759 ENDDO
1760 k = kte
1761 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*( &
1762 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j) &
1763 +muu(i,j )*(cfn*u(i ,k-1,j)+cfn1*u(i,k-2,j))*msfux(i,j) )* (1./12.)*( &
1764 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1765 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1766 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1767 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1769 ENDDO
1770 END IF
1772 IF ( (config_flags%open_xe) .and. ite >= ide-2 ) THEN
1773 i = ide-2
1774 DO j = j_start, jtf
1775 DO k = 2, kte-1
1776 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*( &
1777 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j) &
1778 +muu(i,j )*(u(i,k,j )+u(i,k-1,j ))*msfux(i,j) )* (1./12.)*( &
1779 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1780 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1781 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1782 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1783 ENDDO
1784 k = kte
1785 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*( &
1786 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j) &
1787 +muu(i,j )*(cfn*u(i ,k-1,j)+cfn1*u(i,k-2,j))*msfux(i,j) )* (1./12.)*( &
1788 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1789 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1790 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1791 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1793 ENDDO
1794 END IF
1796 END IF
1798 ! lateral open boundary conditions,
1799 ! start with north and south (y) boundaries
1801 i_start = its
1802 itf=MIN(ite,ide-1)
1804 ! south
1806 IF ( (config_flags%open_ys) .and. jts == jds ) THEN
1808 j=jts
1810 DO k=2,kde
1811 kz = min(k,kde-1)
1812 DO i = its,itf
1813 vb =.5*( fnm(kz)*(v(i,kz ,j+1)+v(i,kz ,j )) &
1814 +fnp(kz)*(v(i,kz-1,j+1)+v(i,kz-1,j )) )
1815 vl=amin1(vb,0.)
1816 ph_tend(i,k,j)=ph_tend(i,k,j)-rdy*mut(i,j)*( &
1817 +vl*(ph_old(i,k,j+1)-ph_old(i,k,j)))
1818 ENDDO
1819 ENDDO
1821 END IF
1823 ! north
1825 IF ( (config_flags%open_ye) .and. jte == jde ) THEN
1827 j=jte-1
1829 DO k=2,kde
1830 kz = min(k,kde-1)
1831 DO i = its,itf
1832 vb=.5*( fnm(kz)*(v(i,kz ,j+1)+v(i,kz ,j)) &
1833 +fnp(kz)*(v(i,kz-1,j+1)+v(i,kz-1,j)) )
1834 vr=amax1(vb,0.)
1835 ph_tend(i,k,j)=ph_tend(i,k,j)-rdy*mut(i,j)*( &
1836 +vr*(ph_old(i,k,j)-ph_old(i,k,j-1)))
1837 ENDDO
1838 ENDDO
1840 END IF
1842 ! now the east and west (y) boundaries
1844 j_start = its
1845 jtf=MIN(jte,jde-1)
1847 ! west
1849 IF ( (config_flags%open_xs) .and. its == ids ) THEN
1851 i=its
1853 DO j = jts,jtf
1854 DO k=2,kde-1
1855 kz = k
1856 ub =.5*( fnm(kz)*(u(i+1,kz ,j)+u(i ,kz ,j)) &
1857 +fnp(kz)*(u(i+1,kz-1,j)+u(i ,kz-1,j)) )
1858 ul=amin1(ub,0.)
1859 ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*( &
1860 +ul*(ph_old(i+1,k,j)-ph_old(i,k,j)))
1861 ENDDO
1863 k = kde
1864 kz = k
1865 ub =.5*( fnm(kz)*(u(i+1,kz ,j)+u(i ,kz ,j)) &
1866 +fnp(kz)*(u(i+1,kz-1,j)+u(i ,kz-1,j)) )
1867 ul=amin1(ub,0.)
1868 ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*( &
1869 +ul*(ph_old(i+1,k,j)-ph_old(i,k,j)))
1870 ENDDO
1872 END IF
1874 ! east
1876 IF ( (config_flags%open_xe) .and. ite == ide ) THEN
1878 i = ite-1
1880 DO j = jts,jtf
1881 DO k=2,kde-1
1882 kz = k
1883 ub=.5*( fnm(kz)*(u(i+1,kz ,j)+u(i,kz ,j)) &
1884 +fnp(kz)*(u(i+1,kz-1,j)+u(i,kz-1,j)) )
1885 ur=amax1(ub,0.)
1886 ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*( &
1887 +ur*(ph_old(i,k,j)-ph_old(i-1,k,j)))
1888 ENDDO
1890 k = kde
1891 kz = k-1
1892 ub=.5*( fnm(kz)*(u(i+1,kz ,j)+u(i,kz ,j)) &
1893 +fnp(kz)*(u(i+1,kz-1,j)+u(i,kz-1,j)) )
1894 ur=amax1(ub,0.)
1895 ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*( &
1896 +ur*(ph_old(i,k,j)-ph_old(i-1,k,j)))
1898 ENDDO
1900 END IF
1902 END SUBROUTINE rhs_ph
1905 !-------------------------------------------------------------------------------
1907 SUBROUTINE horizontal_pressure_gradient( ru_tend,rv_tend, &
1908 ph,alt,p,pb,al,php,cqu,cqv, &
1909 muu,muv,mu,fnm,fnp,rdnw, &
1910 cf1,cf2,cf3,rdx,rdy,msfux,msfuy,&
1911 msfvx,msfvy,msftx,msfty, &
1912 config_flags, non_hydrostatic, &
1913 top_lid, &
1914 ids, ide, jds, jde, kds, kde, &
1915 ims, ime, jms, jme, kms, kme, &
1916 its, ite, jts, jte, kts, kte )
1918 IMPLICIT NONE
1920 ! Input data
1923 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
1925 LOGICAL, INTENT (IN ) :: non_hydrostatic, top_lid
1927 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1928 ims, ime, jms, jme, kms, kme, &
1929 its, ite, jts, jte, kts, kte
1931 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: &
1932 ph, &
1933 alt, &
1934 al, &
1935 p, &
1936 pb, &
1937 php, &
1938 cqu, &
1939 cqv
1942 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: &
1943 ru_tend, &
1944 rv_tend
1946 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: muu, muv, mu, &
1947 msfux, msfuy, &
1948 msfvx, msfvy, &
1949 msftx, msfty
1951 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, fnm, fnp
1953 REAL, INTENT(IN ) :: rdx, rdy, cf1, cf2, cf3
1955 INTEGER :: i,j,k, itf, jtf, ktf, i_start, j_start
1956 REAL, DIMENSION( ims:ime, kms:kme ) :: dpn
1957 REAL :: dpx, dpy
1959 LOGICAL :: specified
1961 !<DESCRIPTION>
1962 !
1963 ! horizontal_pressure_gradient calculates the
1964 ! horizontal pressure gradient terms for the large-timestep tendency
1965 ! in the horizontal momentum equations (u,v).
1966 !
1967 !</DESCRIPTION>
1969 specified = .false.
1970 if(config_flags%specified .or. config_flags%nested) specified = .true.
1972 ! Notes on map scale factors:
1973 ! Calculates the pressure gradient terms in ADT eqns 44 and 45
1974 ! With upper rho -> 'mu', these are:
1975 ! Eqn 30: -mu*(mx/my)*(1/rho)*partial dp/dx
1976 ! Eqn 31: -mu*(my/mx)*(1/rho)*partial dp/dy
1977 !
1978 ! As we are on nu, rather than height, surfaces:
1979 !
1980 ! mu dp/dx = mu alpha partial dp'/dx + (nu mu partial dmubar/dx) alpha'
1981 ! + mu partial dphi'/dx + (partial dphi/dx)*(partial dp'/dnu - mu')
1982 !
1983 ! mu dp/dy = mu alpha partial dp'/dy + (nu mu partial dmubar/dy) alpha'
1984 ! + mu partial dphi'/dy + (partial dphi/dy)*(partial dp'/dnu - mu')
1986 ! start with the north-south (y) pressure gradient
1988 itf=MIN(ite,ide-1)
1989 jtf=jte
1990 ktf=MIN(kte,kde-1)
1991 i_start = its
1992 j_start = jts
1993 IF ( (config_flags%open_ys .or. specified .or. &
1994 config_flags%nested .or. config_flags%polar ) .and. jts == jds ) j_start = jts+1
1995 IF ( (config_flags%open_ye .or. specified .or. &
1996 config_flags%nested .or. config_flags%polar ) .and. jte == jde ) jtf = jtf-1
1998 DO j = j_start, jtf
2000 IF ( non_hydrostatic ) THEN
2002 k=1
2004 DO i = i_start, itf
2005 dpn(i,k) = .5*( cf1*(p(i,k ,j-1)+p(i,k ,j)) &
2006 +cf2*(p(i,k+1,j-1)+p(i,k+1,j)) &
2007 +cf3*(p(i,k+2,j-1)+p(i,k+2,j)) )
2008 dpn(i,kde) = 0.
2009 ENDDO
2010 IF (top_lid) THEN
2011 DO i = i_start, itf
2012 dpn(i,kde) = .5*( cf1*(p(i,kde-1,j-1)+p(i,kde-1,j)) &
2013 +cf2*(p(i,kde-2,j-1)+p(i,kde-2,j)) &
2014 +cf3*(p(i,kde-3,j-1)+p(i,kde-3,j)) )
2015 ENDDO
2016 ENDIF
2018 DO k=2,ktf
2019 DO i = i_start, itf
2020 dpn(i,k) = .5*( fnm(k)*(p(i,k ,j-1)+p(i,k ,j)) &
2021 +fnp(k)*(p(i,k-1,j-1)+p(i,k-1,j)) )
2022 END DO
2023 END DO
2025 ! ADT eqn 45: -mu*(my/mx)*(1/rho)*partial dp/dy
2026 ! [alt, al are 1/rho terms; muv, mu are NOT coupled]
2027 DO K=1,ktf
2028 DO i = i_start, itf
2029 ! Here are mu dp/dy terms 1-3
2030 dpy = (msfvy(i,j)/msfvx(i,j))*.5*rdy*muv(i,j)*( &
2031 (ph (i,k+1,j)-ph (i,k+1,j-1) + ph(i,k,j)-ph(i,k,j-1)) &
2032 +(alt(i,k ,j)+alt(i,k ,j-1))*(p (i,k,j)-p (i,k,j-1)) &
2033 +(al (i,k ,j)+al (i,k ,j-1))*(pb(i,k,j)-pb(i,k,j-1)) )
2034 ! Here is mu dp/dy term 4
2035 dpy = dpy + (msfvy(i,j)/msfvx(i,j))*rdy*(php(i,k,j)-php(i,k,j-1))* &
2036 (rdnw(k)*(dpn(i,k+1)-dpn(i,k))-.5*(mu(i,j-1)+mu(i,j)))
2037 rv_tend(i,k,j) = rv_tend(i,k,j)-cqv(i,k,j)*dpy
2038 END DO
2039 END DO
2041 ELSE
2043 ! ADT eqn 45: -mu*(my/mx)*(1/rho)*partial dp/dy
2044 ! [alt, al are 1/rho terms; muv, mu are NOT coupled]
2045 DO K=1,ktf
2046 DO i = i_start, itf
2047 ! Here are mu dp/dy terms 1-3; term 4 not needed if hydrostatic
2048 dpy = (msfvy(i,j)/msfvx(i,j))*.5*rdy*muv(i,j)*( &
2049 (ph (i,k+1,j)-ph (i,k+1,j-1) + ph(i,k,j)-ph(i,k,j-1)) &
2050 +(alt(i,k ,j)+alt(i,k ,j-1))*(p (i,k,j)-p (i,k,j-1)) &
2051 +(al (i,k ,j)+al (i,k ,j-1))*(pb(i,k,j)-pb(i,k,j-1)) )
2052 rv_tend(i,k,j) = rv_tend(i,k,j)-cqv(i,k,j)*dpy
2053 END DO
2054 END DO
2056 END IF
2058 ENDDO
2060 ! now the east-west (x) pressure gradient
2062 itf=ite
2063 jtf=MIN(jte,jde-1)
2064 ktf=MIN(kte,kde-1)
2065 i_start = its
2066 j_start = jts
2067 IF ( (config_flags%open_xs .or. specified .or. &
2068 config_flags%nested ) .and. its == ids ) i_start = its+1
2069 IF ( (config_flags%open_xe .or. specified .or. &
2070 config_flags%nested ) .and. ite == ide ) itf = itf-1
2071 IF ( config_flags%periodic_x ) i_start = its
2072 IF ( config_flags%periodic_x ) itf=ite
2074 DO j = j_start, jtf
2076 IF ( non_hydrostatic ) THEN
2078 k=1
2080 DO i = i_start, itf
2081 dpn(i,k) = .5*( cf1*(p(i-1,k ,j)+p(i,k ,j)) &
2082 +cf2*(p(i-1,k+1,j)+p(i,k+1,j)) &
2083 +cf3*(p(i-1,k+2,j)+p(i,k+2,j)) )
2084 dpn(i,kde) = 0.
2085 ENDDO
2086 IF (top_lid) THEN
2087 DO i = i_start, itf
2088 dpn(i,kde) = .5*( cf1*(p(i-1,kde-1,j)+p(i,kde-1,j)) &
2089 +cf2*(p(i-1,kde-2,j)+p(i,kde-2,j)) &
2090 +cf3*(p(i-1,kde-3,j)+p(i,kde-3,j)) )
2091 ENDDO
2092 ENDIF
2094 DO k=2,ktf
2095 DO i = i_start, itf
2096 dpn(i,k) = .5*( fnm(k)*(p(i-1,k ,j)+p(i,k ,j)) &
2097 +fnp(k)*(p(i-1,k-1,j)+p(i,k-1,j)) )
2098 END DO
2099 END DO
2101 ! ADT eqn 44: -mu*(mx/my)*(1/rho)*partial dp/dx
2102 ! [alt, al are 1/rho terms; muu, mu are NOT coupled]
2103 DO K=1,ktf
2104 DO i = i_start, itf
2105 ! Here are mu dp/dy terms 1-3
2106 dpx = (msfux(i,j)/msfuy(i,j))*.5*rdx*muu(i,j)*( &
2107 (ph (i,k+1,j)-ph (i-1,k+1,j) + ph(i,k,j)-ph(i-1,k,j)) &
2108 +(alt(i,k ,j)+alt(i-1,k ,j))*(p (i,k,j)-p (i-1,k,j)) &
2109 +(al (i,k ,j)+al (i-1,k ,j))*(pb(i,k,j)-pb(i-1,k,j)) )
2110 ! Here is mu dp/dy term 4
2111 dpx = dpx + (msfux(i,j)/msfuy(i,j))*rdx*(php(i,k,j)-php(i-1,k,j))* &
2112 (rdnw(k)*(dpn(i,k+1)-dpn(i,k))-.5*(mu(i-1,j)+mu(i,j)))
2113 ru_tend(i,k,j) = ru_tend(i,k,j)-cqu(i,k,j)*dpx
2114 END DO
2115 END DO
2117 ELSE
2119 ! ADT eqn 44: -mu*(mx/my)*(1/rho)*partial dp/dx
2120 ! [alt, al are 1/rho terms; muu, mu are NOT coupled]
2121 DO K=1,ktf
2122 DO i = i_start, itf
2123 ! Here are mu dp/dy terms 1-3; term 4 not needed if hydrostatic
2124 dpx = (msfux(i,j)/msfuy(i,j))*.5*rdx*muu(i,j)*( &
2125 (ph (i,k+1,j)-ph (i-1,k+1,j) + ph(i,k,j)-ph(i-1,k,j)) &
2126 +(alt(i,k ,j)+alt(i-1,k ,j))*(p (i,k,j)-p (i-1,k,j)) &
2127 +(al (i,k ,j)+al (i-1,k ,j))*(pb(i,k,j)-pb(i-1,k,j)) )
2128 ru_tend(i,k,j) = ru_tend(i,k,j)-cqu(i,k,j)*dpx
2129 END DO
2130 END DO
2132 END IF
2134 ENDDO
2136 END SUBROUTINE horizontal_pressure_gradient
2138 !-------------------------------------------------------------------------------
2140 SUBROUTINE pg_buoy_w( rw_tend, p, cqw, mu, mub, &
2141 rdnw, rdn, g, msftx, msfty, &
2142 ids, ide, jds, jde, kds, kde, &
2143 ims, ime, jms, jme, kms, kme, &
2144 its, ite, jts, jte, kts, kte )
2146 IMPLICIT NONE
2148 ! Input data
2150 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2151 ims, ime, jms, jme, kms, kme, &
2152 its, ite, jts, jte, kts, kte
2154 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: p
2155 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: cqw
2158 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: rw_tend
2160 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mub, mu, msftx, msfty
2162 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, rdn
2164 REAL, INTENT(IN ) :: g
2166 INTEGER :: itf, jtf, i, j, k
2167 REAL :: cq1, cq2
2170 !<DESCRIPTION>
2171 !
2172 ! pg_buoy_w calculates the
2173 ! vertical pressure gradient and buoyancy terms for the large-timestep
2174 ! tendency in the vertical momentum equation.
2175 !
2176 !</DESCRIPTION>
2178 ! BUOYANCY AND PRESSURE GRADIENT TERM IN W EQUATION AT TIME T
2180 ! Map scale factor notes
2181 ! ADT eqn 46 RHS terms 6 and 7 (where 7 is "-rho g")
2182 ! Dividing by my, and using mu and nu (see Klemp et al. eqns 32, 40)
2183 ! term 6: +(g/my) partial dp'/dnu
2184 ! term 7: -(g/my) mu'
2185 !
2186 ! For moisture-free atmosphere, cq1=1, cq2=0
2187 ! => (1./msft(i,j)) * g * [rdn(k)*{p(i,k,j)-p(i,k-1,j)}-mu(i,j)]
2189 itf=MIN(ite,ide-1)
2190 jtf=MIN(jte,jde-1)
2192 DO j = jts,jtf
2194 k=kde
2195 DO i=its,itf
2196 cq1 = 1./(1.+cqw(i,k-1,j))
2197 cq2 = cqw(i,k-1,j)*cq1
2198 rw_tend(i,k,j) = rw_tend(i,k,j)+(1./msfty(i,j))*g*( &
2199 cq1*2.*rdnw(k-1)*( -p(i,k-1,j)) &
2200 -mu(i,j)-cq2*mub(i,j) )
2201 END DO
2203 DO k = 2, kde-1
2204 DO i = its,itf
2205 cq1 = 1./(1.+cqw(i,k,j))
2206 cq2 = cqw(i,k,j)*cq1
2207 cqw(i,k,j) = cq1
2208 rw_tend(i,k,j) = rw_tend(i,k,j)+(1./msfty(i,j))*g*( &
2209 cq1*rdn(k)*(p(i,k,j)-p(i,k-1,j)) &
2210 -mu(i,j)-cq2*mub(i,j) )
2211 END DO
2212 ENDDO
2215 ENDDO
2217 END SUBROUTINE pg_buoy_w
2219 !-------------------------------------------------------------------------------
2221 SUBROUTINE w_damp( rw_tend, max_vert_cfl,max_horiz_cfl, &
2222 u, v, ww, w, mut, rdnw, &
2223 rdx, rdy, msfux, msfuy, &
2224 msfvx, msfvy, dt, &
2225 config_flags, &
2226 ids, ide, jds, jde, kds, kde, &
2227 ims, ime, jms, jme, kms, kme, &
2228 its, ite, jts, jte, kts, kte )
2230 USE module_llxy
2231 IMPLICIT NONE
2233 ! Input data
2235 TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
2237 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2238 ims, ime, jms, jme, kms, kme, &
2239 its, ite, jts, jte, kts, kte
2241 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: u, v, ww, w
2243 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: rw_tend
2245 REAL, INTENT(OUT) :: max_vert_cfl
2246 REAL, INTENT(OUT) :: max_horiz_cfl
2247 REAL :: horiz_cfl
2249 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mut
2251 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw
2253 REAL, INTENT(IN) :: dt
2254 REAL, INTENT(IN) :: rdx, rdy
2255 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, msfuy
2256 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfvx, msfvy
2258 REAL :: vert_cfl, cf_n, cf_d, maxdub, maxdeta
2260 INTEGER :: itf, jtf, i, j, k, maxi, maxj, maxk
2261 INTEGER :: some
2262 CHARACTER*512 :: temp
2264 CHARACTER (LEN=256) :: time_str
2265 CHARACTER (LEN=256) :: grid_str
2267 integer :: total
2268 REAL :: msfuxt , msfxffl
2270 !<DESCRIPTION>
2271 !
2272 ! w_damp computes a damping term for the vertical velocity when the
2273 ! vertical Courant number is too large. This was found to be preferable to
2274 ! decreasing the timestep or increasing the diffusion in real-data applications
2275 ! that produced potentially-unstable large vertical velocities because of
2276 ! unphysically large heating rates coming from the cumulus parameterization
2277 ! schemes run at moderately high resolutions (dx ~ O(10) km).
2278 !
2279 ! Additionally, w_damp returns the maximum cfl values due to vertical motion and
2280 ! horizontal motion. These values are returned via the max_vert_cfl and
2281 ! max_horiz_cfl variables. (Added by T. Hutchinson, WSI, 3/5/2007)
2282 !
2283 !</DESCRIPTION>
2285 itf=MIN(ite,ide-1)
2286 jtf=MIN(jte,jde-1)
2288 some = 0
2289 max_vert_cfl = 0.
2290 max_horiz_cfl = 0.
2291 total = 0
2293 IF(config_flags%map_proj == PROJ_CASSINI ) then
2294 msfxffl = 1.0/COS(config_flags%fft_filter_lat*degrad)
2295 END IF
2297 IF ( config_flags%w_damping == 1 ) THEN
2298 DO j = jts,jtf
2300 DO k = 2, kde-1
2301 DO i = its,itf
2302 #if 1
2303 IF(config_flags%map_proj == PROJ_CASSINI ) then
2304 msfuxt = MIN(msfux(i,j), msfxffl)
2305 ELSE
2306 msfuxt = msfux(i,j)
2307 END IF
2308 vert_cfl = abs(ww(i,k,j)/mut(i,j)*rdnw(k)*dt)
2310 IF ( vert_cfl > max_vert_cfl ) THEN
2311 max_vert_cfl = vert_cfl ; maxi = i ; maxj = j ; maxk = k
2312 maxdub = w(i,k,j) ; maxdeta = -1./rdnw(k)
2313 ENDIF
2315 horiz_cfl = max( abs(u(i,k,j) * rdx * msfuxt * dt), &
2316 abs(v(i,k,j) * rdy * msfvy(i,j) * dt) )
2317 if (horiz_cfl > max_horiz_cfl) then
2318 max_horiz_cfl = horiz_cfl
2319 endif
2321 if(vert_cfl .gt. w_beta)then
2322 #else
2323 ! restructure to get rid of divide
2324 !
2325 ! This had been used for efficiency, but with the addition of returning the cfl values,
2326 ! the old version (above) was reinstated. (T. Hutchinson, 3/5/2007)
2327 !
2328 cf_n = abs(ww(i,k,j)*rdnw(k)*dt)
2329 cf_d = abs(mut(i,j))
2330 if(cf_n .gt. cf_d*w_beta )then
2331 #endif
2333 WRITE(temp,*)i,j,k,' vert_cfl,w,d(eta)=',vert_cfl,w(i,k,j),-1./rdnw(k)
2334 CALL wrf_debug ( 100 , TRIM(temp) )
2335 if ( vert_cfl > 2. ) some = some + 1
2336 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)
2337 endif
2338 END DO
2339 ENDDO
2340 ENDDO
2341 ELSE
2342 ! just print
2343 DO j = jts,jtf
2345 DO k = 2, kde-1
2346 DO i = its,itf
2348 #if 1
2349 IF(config_flags%map_proj == PROJ_CASSINI ) then
2350 msfuxt = MIN(msfux(i,j), msfxffl)
2351 ELSE
2352 msfuxt = msfux(i,j)
2353 END IF
2354 vert_cfl = abs(ww(i,k,j)/mut(i,j)*rdnw(k)*dt)
2356 IF ( vert_cfl > max_vert_cfl ) THEN
2357 max_vert_cfl = vert_cfl ; maxi = i ; maxj = j ; maxk = k
2358 maxdub = w(i,k,j) ; maxdeta = -1./rdnw(k)
2359 ENDIF
2361 horiz_cfl = max( abs(u(i,k,j) * rdx * msfuxt * dt), &
2362 abs(v(i,k,j) * rdy * msfvy(i,j) * dt) )
2364 if (horiz_cfl > max_horiz_cfl) then
2365 max_horiz_cfl = horiz_cfl
2366 endif
2368 if(vert_cfl .gt. w_beta)then
2369 #else
2370 ! restructure to get rid of divide
2371 !
2372 ! This had been used for efficiency, but with the addition of returning the cfl values,
2373 ! the old version (above) was reinstated. (T. Hutchinson, 3/5/2007)
2374 !
2375 cf_n = abs(ww(i,k,j)*rdnw(k)*dt)
2376 cf_d = abs(mut(i,j))
2377 if(cf_n .gt. cf_d*w_beta )then
2378 #endif
2379 WRITE(temp,*)i,j,k,' vert_cfl,w,d(eta)=',vert_cfl,w(i,k,j),-1./rdnw(k)
2380 CALL wrf_debug ( 100 , TRIM(temp) )
2381 if ( vert_cfl > 2. ) some = some + 1
2382 endif
2383 END DO
2384 ENDDO
2385 ENDDO
2386 ENDIF
2387 IF ( some .GT. 0 ) THEN
2388 CALL get_current_time_string( time_str )
2389 CALL get_current_grid_name( grid_str )
2390 WRITE(wrf_err_message,*)some, &
2391 ' points exceeded cfl=2 in domain '//TRIM(grid_str)//' at time '//TRIM(time_str)//' hours'
2392 CALL wrf_debug ( 0 , TRIM(wrf_err_message) )
2393 WRITE(wrf_err_message,*)'MAX AT i,j,k: ',maxi,maxj,maxk,' vert_cfl,w,d(eta)=',max_vert_cfl, &
2394 maxdub,maxdeta
2395 CALL wrf_debug ( 0 , TRIM(wrf_err_message) )
2396 ENDIF
2398 END SUBROUTINE w_damp
2400 !-------------------------------------------------------------------------------
2402 SUBROUTINE horizontal_diffusion ( name, field, tendency, mu, &
2403 config_flags, &
2404 msfux, msfuy, msfvx, msfvx_inv, &
2405 msfvy, msftx, msfty, &
2406 khdif, xkmhd, rdx, rdy, &
2407 ids, ide, jds, jde, kds, kde, &
2408 ims, ime, jms, jme, kms, kme, &
2409 its, ite, jts, jte, kts, kte )
2411 IMPLICIT NONE
2413 ! Input data
2415 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2417 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2418 ims, ime, jms, jme, kms, kme, &
2419 its, ite, jts, jte, kts, kte
2421 CHARACTER(LEN=1) , INTENT(IN ) :: name
2423 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, xkmhd
2425 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2427 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
2429 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
2430 msfuy, &
2431 msfvx, &
2432 msfvx_inv, &
2433 msfvy, &
2434 msftx, &
2435 msfty
2437 REAL , INTENT(IN ) :: rdx, &
2438 rdy, &
2439 khdif
2441 ! Local data
2443 INTEGER :: i, j, k, itf, jtf, ktf
2445 INTEGER :: i_start, i_end, j_start, j_end
2447 REAL :: mrdx, mkrdxm, mkrdxp, &
2448 mrdy, mkrdym, mkrdyp
2450 LOGICAL :: specified
2452 !<DESCRIPTION>
2453 !
2454 ! horizontal_diffusion computes the horizontal diffusion tendency
2455 ! on model horizontal coordinate surfaces.
2456 !
2457 !</DESCRIPTION>
2459 specified = .false.
2460 if(config_flags%specified .or. config_flags%nested) specified = .true.
2462 ktf=MIN(kte,kde-1)
2464 IF (name .EQ. 'u') THEN
2466 i_start = its
2467 i_end = ite
2468 j_start = jts
2469 j_end = MIN(jte,jde-1)
2471 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2472 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite)
2473 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2474 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2475 IF ( config_flags%periodic_x ) i_start = its
2476 IF ( config_flags%periodic_x ) i_end = ite
2479 DO j = j_start, j_end
2480 DO k=kts,ktf
2481 DO i = i_start, i_end
2483 ! The interior is grad: (m_x*d/dx), the exterior is div: (m_x*m_y*d/dx(/m_y))
2484 ! setting up different averagings of m^2 partial d/dX and m^2 partial d/dY
2486 mkrdxm=(msftx(i-1,j)/msfty(i-1,j))*mu(i-1,j)*xkmhd(i-1,k,j)*rdx
2487 mkrdxp=(msftx(i,j)/msfty(i,j))*mu(i,j)*xkmhd(i,k,j)*rdx
2488 mrdx=msfux(i,j)*msfuy(i,j)*rdx
2489 mkrdym=( (msfuy(i,j)+msfuy(i,j-1))/(msfux(i,j)+msfux(i,j-1)) )* &
2490 0.25*(mu(i,j)+mu(i,j-1)+mu(i-1,j-1)+mu(i-1,j))* &
2491 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i-1,k,j-1)+xkmhd(i-1,k,j))*rdy
2492 mkrdyp=( (msfuy(i,j)+msfuy(i,j+1))/(msfux(i,j)+msfux(i,j+1)) )* &
2493 0.25*(mu(i,j)+mu(i,j+1)+mu(i-1,j+1)+mu(i-1,j))* &
2494 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j+1)+xkmhd(i-1,k,j+1)+xkmhd(i-1,k,j))*rdy
2495 ! need to do four-corners (t) for diffusion coefficient as there are
2496 ! no values at u,v points
2497 ! msfuy - has to be y as part of d/dY
2498 ! has to be u as we're at a u point
2499 mrdy=msfux(i,j)*msfuy(i,j)*rdy
2501 ! correctly averaged version of rho~ * m^2 *
2502 ! [partial d/dX(partial du^/dX) + partial d/dY(partial du^/dY)]
2503 tendency(i,k,j)=tendency(i,k,j)+( &
2504 mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) &
2505 -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2506 +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) &
2507 -mkrdym*(field(i,k,j )-field(i,k,j-1))))
2508 ENDDO
2509 ENDDO
2510 ENDDO
2512 ELSE IF (name .EQ. 'v')THEN
2514 i_start = its
2515 i_end = MIN(ite,ide-1)
2516 j_start = jts
2517 j_end = jte
2519 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2520 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2521 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2522 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte)
2523 IF ( config_flags%periodic_x ) i_start = its
2524 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2525 IF ( config_flags%polar ) j_start = MAX(jds+1,jts)
2526 IF ( config_flags%polar ) j_end = MIN(jde-1,jte)
2528 DO j = j_start, j_end
2529 DO k=kts,ktf
2530 DO i = i_start, i_end
2532 mkrdxm=( (msfvx(i,j)+msfvx(i-1,j))/(msfvy(i,j)+msfvy(i-1,j)) )* &
2533 0.25*(mu(i,j)+mu(i,j-1)+mu(i-1,j-1)+mu(i-1,j))* &
2534 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i-1,k,j-1)+xkmhd(i-1,k,j))*rdx
2535 mkrdxp=( (msfvx(i,j)+msfvx(i+1,j))/(msfvy(i,j)+msfvy(i+1,j)) )* &
2536 0.25*(mu(i,j)+mu(i,j-1)+mu(i+1,j-1)+mu(i+1,j))* &
2537 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i+1,k,j-1)+xkmhd(i+1,k,j))*rdx
2538 mrdx=msfvx(i,j)*msfvy(i,j)*rdx
2539 mkrdym=(msfty(i,j-1)/msftx(i,j-1))*xkmhd(i,k,j-1)*rdy
2540 mkrdyp=(msfty(i,j)/msftx(i,j))*xkmhd(i,k,j)*rdy
2541 mrdy=msfvx(i,j)*msfvy(i,j)*rdy
2543 tendency(i,k,j)=tendency(i,k,j)+( &
2544 mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) &
2545 -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2546 +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) &
2547 -mkrdym*(field(i,k,j )-field(i,k,j-1))))
2548 ENDDO
2549 ENDDO
2550 ENDDO
2552 ELSE IF (name .EQ. 'w')THEN
2554 i_start = its
2555 i_end = MIN(ite,ide-1)
2556 j_start = jts
2557 j_end = MIN(jte,jde-1)
2559 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2560 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2561 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2562 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2563 IF ( config_flags%periodic_x ) i_start = its
2564 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2566 DO j = j_start, j_end
2567 DO k=kts+1,ktf
2568 DO i = i_start, i_end
2570 mkrdxm=(msfux(i,j)/msfuy(i,j))* &
2571 0.25*(mu(i,j)+mu(i-1,j)+mu(i,j)+mu(i-1,j))* &
2572 0.25*(xkmhd(i,k,j)+xkmhd(i-1,k,j)+xkmhd(i,k-1,j)+xkmhd(i-1,k-1,j))*rdx
2573 mkrdxp=(msfux(i+1,j)/msfuy(i+1,j))* &
2574 0.25*(mu(i+1,j)+mu(i,j)+mu(i+1,j)+mu(i,j))* &
2575 0.25*(xkmhd(i+1,k,j)+xkmhd(i,k,j)+xkmhd(i+1,k-1,j)+xkmhd(i,k-1,j))*rdx
2576 mrdx=msftx(i,j)*msfty(i,j)*rdx
2577 ! mkrdym=(msfvy(i,j)/msfvx(i,j))* &
2578 mkrdym=(msfvy(i,j)*msfvx_inv(i,j))* &
2579 0.25*(mu(i,j)+mu(i,j-1)+mu(i,j)+mu(i,j-1))* &
2580 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i,k-1,j)+xkmhd(i,k-1,j-1))*rdy
2581 ! mkrdyp=(msfvy(i,j+1)/msfvx(i,j+1))* &
2582 mkrdyp=(msfvy(i,j+1)*msfvx_inv(i,j+1))* &
2583 0.25*(mu(i,j+1)+mu(i,j)+mu(i,j+1)+mu(i,j))* &
2584 0.25*(xkmhd(i,k,j+1)+xkmhd(i,k,j)+xkmhd(i,k-1,j+1)+xkmhd(i,k-1,j))*rdy
2585 mrdy=msftx(i,j)*msfty(i,j)*rdy
2587 tendency(i,k,j)=tendency(i,k,j)+( &
2588 mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) &
2589 -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2590 +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) &
2591 -mkrdym*(field(i,k,j )-field(i,k,j-1))))
2592 ENDDO
2593 ENDDO
2594 ENDDO
2596 ELSE
2599 i_start = its
2600 i_end = MIN(ite,ide-1)
2601 j_start = jts
2602 j_end = MIN(jte,jde-1)
2604 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2605 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2606 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2607 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2608 IF ( config_flags%periodic_x ) i_start = its
2609 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2611 DO j = j_start, j_end
2612 DO k=kts,ktf
2613 DO i = i_start, i_end
2615 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
2616 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
2617 mrdx=msftx(i,j)*msfty(i,j)*rdx
2618 ! 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
2619 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
2620 ! 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
2621 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
2622 mrdy=msftx(i,j)*msfty(i,j)*rdy
2624 tendency(i,k,j)=tendency(i,k,j)+( &
2625 mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) &
2626 -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2627 +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) &
2628 -mkrdym*(field(i,k,j )-field(i,k,j-1))))
2629 ENDDO
2630 ENDDO
2631 ENDDO
2633 ENDIF
2635 END SUBROUTINE horizontal_diffusion
2637 !-----------------------------------------------------------------------------------------
2639 SUBROUTINE horizontal_diffusion_3dmp ( name, field, tendency, mu, &
2640 config_flags, base_3d, &
2641 msfux, msfuy, msfvx, msfvx_inv, &
2642 msfvy, msftx, msfty, &
2643 khdif, xkmhd, rdx, rdy, &
2644 ids, ide, jds, jde, kds, kde, &
2645 ims, ime, jms, jme, kms, kme, &
2646 its, ite, jts, jte, kts, kte )
2648 IMPLICIT NONE
2650 ! Input data
2652 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2654 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2655 ims, ime, jms, jme, kms, kme, &
2656 its, ite, jts, jte, kts, kte
2658 CHARACTER(LEN=1) , INTENT(IN ) :: name
2660 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, &
2661 xkmhd, &
2662 base_3d
2664 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2666 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
2668 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
2669 msfuy, &
2670 msfvx, &
2671 msfvx_inv, &
2672 msfvy, &
2673 msftx, &
2674 msfty
2676 REAL , INTENT(IN ) :: rdx, &
2677 rdy, &
2678 khdif
2680 ! Local data
2682 INTEGER :: i, j, k, itf, jtf, ktf
2684 INTEGER :: i_start, i_end, j_start, j_end
2686 REAL :: mrdx, mkrdxm, mkrdxp, &
2687 mrdy, mkrdym, mkrdyp
2689 LOGICAL :: specified
2691 !<DESCRIPTION>
2692 !
2693 ! horizontal_diffusion_3dmp computes the horizontal diffusion tendency
2694 ! on model horizontal coordinate surfaces. This routine computes diffusion
2695 ! a perturbation scalar (field-base_3d).
2696 !
2697 !</DESCRIPTION>
2699 specified = .false.
2700 if(config_flags%specified .or. config_flags%nested) specified = .true.
2702 ktf=MIN(kte,kde-1)
2704 i_start = its
2705 i_end = MIN(ite,ide-1)
2706 j_start = jts
2707 j_end = MIN(jte,jde-1)
2709 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2710 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2711 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2712 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2713 IF ( config_flags%periodic_x ) i_start = its
2714 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2716 DO j = j_start, j_end
2717 DO k=kts,ktf
2718 DO i = i_start, i_end
2720 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
2721 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
2722 mrdx=msftx(i,j)*msfty(i,j)*rdx
2723 ! 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
2724 ! 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
2725 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
2726 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
2727 mrdy=msftx(i,j)*msfty(i,j)*rdy
2729 tendency(i,k,j)=tendency(i,k,j)+( &
2730 mrdx*( mkrdxp*( field(i+1,k,j) -field(i ,k,j) &
2731 -base_3d(i+1,k,j)+base_3d(i ,k,j) ) &
2732 -mkrdxm*( field(i ,k,j) -field(i-1,k,j) &
2733 -base_3d(i ,k,j)+base_3d(i-1,k,j) ) ) &
2734 +mrdy*( mkrdyp*( field(i,k,j+1) -field(i,k,j ) &
2735 -base_3d(i,k,j+1)+base_3d(i,k,j ) ) &
2736 -mkrdym*( field(i,k,j ) -field(i,k,j-1) &
2737 -base_3d(i,k,j )+base_3d(i,k,j-1) ) ) &
2738 )
2739 ENDDO
2740 ENDDO
2741 ENDDO
2743 END SUBROUTINE horizontal_diffusion_3dmp
2745 !-----------------------------------------------------------------------------------------
2747 SUBROUTINE vertical_diffusion ( name, field, tendency, &
2748 config_flags, &
2749 alt, mut, rdn, rdnw, kvdif, &
2750 ids, ide, jds, jde, kds, kde, &
2751 ims, ime, jms, jme, kms, kme, &
2752 its, ite, jts, jte, kts, kte )
2755 IMPLICIT NONE
2757 ! Input data
2759 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2761 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2762 ims, ime, jms, jme, kms, kme, &
2763 its, ite, jts, jte, kts, kte
2765 CHARACTER(LEN=1) , INTENT(IN ) :: name
2767 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
2768 INTENT(IN ) :: field, &
2769 alt
2771 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2773 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
2775 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw
2777 REAL , INTENT(IN ) :: kvdif
2779 ! Local data
2781 INTEGER :: i, j, k, itf, jtf, ktf
2782 INTEGER :: i_start, i_end, j_start, j_end
2784 REAL , DIMENSION(its:ite, jts:jte) :: vfluxm, vfluxp, zz
2785 REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
2787 REAL :: rdz
2789 LOGICAL :: specified
2791 !<DESCRIPTION>
2792 !
2793 ! vertical_diffusion
2794 ! computes vertical diffusion tendency.
2795 !
2796 !</DESCRIPTION>
2798 specified = .false.
2799 if(config_flags%specified .or. config_flags%nested) specified = .true.
2801 ktf=MIN(kte,kde-1)
2803 IF (name .EQ. 'w')THEN
2806 i_start = its
2807 i_end = MIN(ite,ide-1)
2808 j_start = jts
2809 j_end = MIN(jte,jde-1)
2811 j_loop_w : DO j = j_start, j_end
2813 DO k=kts,ktf-1
2814 DO i = i_start, i_end
2815 vflux(i,k)= (kvdif/alt(i,k,j))*rdnw(k)*(field(i,k+1,j)-field(i,k,j))
2816 ENDDO
2817 ENDDO
2819 DO i = i_start, i_end
2820 vflux(i,ktf)=0.
2821 ENDDO
2823 DO k=kts+1,ktf
2824 DO i = i_start, i_end
2825 tendency(i,k,j)=tendency(i,k,j) &
2826 +rdn(k)*g*g/mut(i,j)/(0.5*(alt(i,k,j)+alt(i,k-1,j))) &
2827 *(vflux(i,k)-vflux(i,k-1))
2828 ENDDO
2829 ENDDO
2831 ENDDO j_loop_w
2833 ELSE IF(name .EQ. 'm')THEN
2835 i_start = its
2836 i_end = MIN(ite,ide-1)
2837 j_start = jts
2838 j_end = MIN(jte,jde-1)
2840 j_loop_s : DO j = j_start, j_end
2842 DO k=kts,ktf-1
2843 DO i = i_start, i_end
2844 vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j))) &
2845 *(field(i,k+1,j)-field(i,k,j))
2846 ENDDO
2847 ENDDO
2849 DO i = i_start, i_end
2850 vflux(i,0)=vflux(i,1)
2851 ENDDO
2853 DO i = i_start, i_end
2854 vflux(i,ktf)=0.
2855 ENDDO
2857 DO k=kts,ktf
2858 DO i = i_start, i_end
2859 tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) &
2860 *rdnw(k)*(vflux(i,k)-vflux(i,k-1))
2861 ENDDO
2862 ENDDO
2864 ENDDO j_loop_s
2866 ENDIF
2868 END SUBROUTINE vertical_diffusion
2871 !-------------------------------------------------------------------------------
2873 SUBROUTINE vertical_diffusion_mp ( field, tendency, config_flags, &
2874 base, &
2875 alt, mut, rdn, rdnw, kvdif, &
2876 ids, ide, jds, jde, kds, kde, &
2877 ims, ime, jms, jme, kms, kme, &
2878 its, ite, jts, jte, kts, kte )
2881 IMPLICIT NONE
2883 ! Input data
2885 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2887 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2888 ims, ime, jms, jme, kms, kme, &
2889 its, ite, jts, jte, kts, kte
2891 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
2892 INTENT(IN ) :: field, &
2893 alt
2895 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2897 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
2899 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, &
2900 rdnw, &
2901 base
2903 REAL , INTENT(IN ) :: kvdif
2905 ! Local data
2907 INTEGER :: i, j, k, itf, jtf, ktf
2908 INTEGER :: i_start, i_end, j_start, j_end
2910 REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
2912 REAL :: rdz
2914 LOGICAL :: specified
2916 !<DESCRIPTION>
2917 !
2918 ! vertical_diffusion_mp
2919 ! computes vertical diffusion tendency of a perturbation variable
2920 ! (field-base). Note that base as a 1D (k) field.
2921 !
2922 !</DESCRIPTION>
2924 specified = .false.
2925 if(config_flags%specified .or. config_flags%nested) specified = .true.
2927 ktf=MIN(kte,kde-1)
2929 i_start = its
2930 i_end = MIN(ite,ide-1)
2931 j_start = jts
2932 j_end = MIN(jte,jde-1)
2934 j_loop_s : DO j = j_start, j_end
2936 DO k=kts,ktf-1
2937 DO i = i_start, i_end
2938 vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j))) &
2939 *(field(i,k+1,j)-field(i,k,j)-base(k+1)+base(k))
2940 ENDDO
2941 ENDDO
2943 DO i = i_start, i_end
2944 vflux(i,0)=vflux(i,1)
2945 ENDDO
2947 DO i = i_start, i_end
2948 vflux(i,ktf)=0.
2949 ENDDO
2951 DO k=kts,ktf
2952 DO i = i_start, i_end
2953 tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) &
2954 *rdnw(k)*(vflux(i,k)-vflux(i,k-1))
2955 ENDDO
2956 ENDDO
2958 ENDDO j_loop_s
2960 END SUBROUTINE vertical_diffusion_mp
2963 !-------------------------------------------------------------------------------
2965 SUBROUTINE vertical_diffusion_3dmp ( field, tendency, config_flags, &
2966 base_3d, &
2967 alt, mut, rdn, rdnw, kvdif, &
2968 ids, ide, jds, jde, kds, kde, &
2969 ims, ime, jms, jme, kms, kme, &
2970 its, ite, jts, jte, kts, kte )
2973 IMPLICIT NONE
2975 ! Input data
2977 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2979 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2980 ims, ime, jms, jme, kms, kme, &
2981 its, ite, jts, jte, kts, kte
2983 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
2984 INTENT(IN ) :: field, &
2985 alt, &
2986 base_3d
2988 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2990 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
2992 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, &
2993 rdnw
2995 REAL , INTENT(IN ) :: kvdif
2997 ! Local data
2999 INTEGER :: i, j, k, itf, jtf, ktf
3000 INTEGER :: i_start, i_end, j_start, j_end
3002 REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3004 REAL :: rdz
3006 LOGICAL :: specified
3008 !<DESCRIPTION>
3009 !
3010 ! vertical_diffusion_3dmp
3011 ! computes vertical diffusion tendency of a perturbation variable
3012 ! (field-base_3d).
3013 !
3014 !</DESCRIPTION>
3016 specified = .false.
3017 if(config_flags%specified .or. config_flags%nested) specified = .true.
3019 ktf=MIN(kte,kde-1)
3021 i_start = its
3022 i_end = MIN(ite,ide-1)
3023 j_start = jts
3024 j_end = MIN(jte,jde-1)
3026 j_loop_s : DO j = j_start, j_end
3028 DO k=kts,ktf-1
3029 DO i = i_start, i_end
3030 vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j))) &
3031 *( field(i,k+1,j) -field(i,k,j) &
3032 -base_3d(i,k+1,j)+base_3d(i,k,j) )
3033 ENDDO
3034 ENDDO
3036 DO i = i_start, i_end
3037 vflux(i,0)=vflux(i,1)
3038 ENDDO
3040 DO i = i_start, i_end
3041 vflux(i,ktf)=0.
3042 ENDDO
3044 DO k=kts,ktf
3045 DO i = i_start, i_end
3046 tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) &
3047 *rdnw(k)*(vflux(i,k)-vflux(i,k-1))
3048 ENDDO
3049 ENDDO
3051 ENDDO j_loop_s
3053 END SUBROUTINE vertical_diffusion_3dmp
3056 !-------------------------------------------------------------------------------
3059 SUBROUTINE vertical_diffusion_u ( field, tendency, &
3060 config_flags, u_base, &
3061 alt, muu, rdn, rdnw, kvdif, &
3062 ids, ide, jds, jde, kds, kde, &
3063 ims, ime, jms, jme, kms, kme, &
3064 its, ite, jts, jte, kts, kte )
3067 IMPLICIT NONE
3069 ! Input data
3071 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3073 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3074 ims, ime, jms, jme, kms, kme, &
3075 its, ite, jts, jte, kts, kte
3077 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3078 INTENT(IN ) :: field, &
3079 alt
3081 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3083 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muu
3085 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw, u_base
3087 REAL , INTENT(IN ) :: kvdif
3089 ! Local data
3091 INTEGER :: i, j, k, itf, jtf, ktf
3092 INTEGER :: i_start, i_end, j_start, j_end
3094 REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3096 REAL :: rdz, zz
3098 LOGICAL :: specified
3100 !<DESCRIPTION>
3101 !
3102 ! vertical_diffusion_u computes vertical diffusion tendency for
3103 ! the u momentum equation. This routine assumes a constant eddy
3104 ! viscosity kvdif.
3105 !
3106 !</DESCRIPTION>
3108 specified = .false.
3109 if(config_flags%specified .or. config_flags%nested) specified = .true.
3111 ktf=MIN(kte,kde-1)
3113 i_start = its
3114 i_end = ite
3115 j_start = jts
3116 j_end = MIN(jte,jde-1)
3118 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
3119 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite)
3120 IF ( config_flags%periodic_x ) i_start = its
3121 IF ( config_flags%periodic_x ) i_end = ite
3124 j_loop_u : DO j = j_start, j_end
3126 DO k=kts,ktf-1
3127 DO i = i_start, i_end
3128 vflux(i,k)=kvdif*rdn(k+1)/(0.25*( alt(i ,k ,j) &
3129 +alt(i-1,k ,j) &
3130 +alt(i ,k+1,j) &
3131 +alt(i-1,k+1,j) ) ) &
3132 *(field(i,k+1,j)-field(i,k,j) &
3133 -u_base(k+1) +u_base(k) )
3134 ENDDO
3135 ENDDO
3137 DO i = i_start, i_end
3138 vflux(i,0)=vflux(i,1)
3139 ENDDO
3141 DO i = i_start, i_end
3142 vflux(i,ktf)=0.
3143 ENDDO
3145 DO k=kts,ktf-1
3146 DO i = i_start, i_end
3147 tendency(i,k,j)=tendency(i,k,j)+ &
3148 g*g*rdnw(k)/muu(i,j)/(0.5*(alt(i-1,k,j)+alt(i,k,j)))* &
3149 (vflux(i,k)-vflux(i,k-1))
3150 ENDDO
3151 ENDDO
3153 ENDDO j_loop_u
3155 END SUBROUTINE vertical_diffusion_u
3157 !-------------------------------------------------------------------------------
3160 SUBROUTINE vertical_diffusion_v ( field, tendency, &
3161 config_flags, v_base, &
3162 alt, muv, rdn, rdnw, kvdif, &
3163 ids, ide, jds, jde, kds, kde, &
3164 ims, ime, jms, jme, kms, kme, &
3165 its, ite, jts, jte, kts, kte )
3168 IMPLICIT NONE
3170 ! Input data
3172 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3174 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3175 ims, ime, jms, jme, kms, kme, &
3176 its, ite, jts, jte, kts, kte
3178 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3179 INTENT(IN ) :: field, &
3180 alt
3181 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw, v_base
3183 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3185 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muv
3187 REAL , INTENT(IN ) :: kvdif
3189 ! Local data
3191 INTEGER :: i, j, k, itf, jtf, ktf, jm1
3192 INTEGER :: i_start, i_end, j_start, j_end
3194 REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3196 REAL :: rdz, zz
3198 LOGICAL :: specified
3200 !<DESCRIPTION>
3201 !
3202 ! vertical_diffusion_v computes vertical diffusion tendency for
3203 ! the v momentum equation. This routine assumes a constant eddy
3204 ! viscosity kvdif.
3205 !
3206 !</DESCRIPTION>
3208 specified = .false.
3209 if(config_flags%specified .or. config_flags%nested) specified = .true.
3211 ktf=MIN(kte,kde-1)
3213 i_start = its
3214 i_end = MIN(ite,ide-1)
3215 j_start = jts
3216 j_end = MIN(jte,jde-1)
3218 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
3219 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte)
3221 j_loop_v : DO j = j_start, j_end
3222 ! jm1 = max(j-1,1)
3223 jm1 = j-1
3225 DO k=kts,ktf-1
3226 DO i = i_start, i_end
3227 vflux(i,k)=kvdif*rdn(k+1)/(0.25*( alt(i,k ,j ) &
3228 +alt(i,k ,jm1) &
3229 +alt(i,k+1,j ) &
3230 +alt(i,k+1,jm1) ) ) &
3231 *(field(i,k+1,j)-field(i,k,j) &
3232 -v_base(k+1) +v_base(k) )
3233 ENDDO
3234 ENDDO
3236 DO i = i_start, i_end
3237 vflux(i,0)=vflux(i,1)
3238 ENDDO
3240 DO i = i_start, i_end
3241 vflux(i,ktf)=0.
3242 ENDDO
3244 DO k=kts,ktf-1
3245 DO i = i_start, i_end
3246 tendency(i,k,j)=tendency(i,k,j)+ &
3247 g*g*rdnw(k)/muv(i,j)/(0.5*(alt(i,k,jm1)+alt(i,k,j)))* &
3248 (vflux(i,k)-vflux(i,k-1))
3249 ENDDO
3250 ENDDO
3252 ENDDO j_loop_v
3254 END SUBROUTINE vertical_diffusion_v
3256 !*************** end new mass coordinate routines
3258 !-------------------------------------------------------------------------------
3260 SUBROUTINE calculate_full ( rfield, rfieldb, rfieldp, &
3261 ids, ide, jds, jde, kds, kde, &
3262 ims, ime, jms, jme, kms, kme, &
3263 its, ite, jts, jte, kts, kte )
3265 IMPLICIT NONE
3267 ! Input data
3269 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3270 ims, ime, jms, jme, kms, kme, &
3271 its, ite, jts, jte, kts, kte
3273 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rfieldb, &
3274 rfieldp
3276 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: rfield
3278 ! Local indices.
3280 INTEGER :: i, j, k, itf, jtf, ktf
3282 !<DESCRIPTION>
3283 !
3284 ! calculate_full
3285 ! calculates full 3D field from pertubation and base field.
3286 !
3287 !</DESCRIPTION>
3289 itf=MIN(ite,ide-1)
3290 jtf=MIN(jte,jde-1)
3291 ktf=MIN(kte,kde-1)
3293 DO j=jts,jtf
3294 DO k=kts,ktf
3295 DO i=its,itf
3296 rfield(i,k,j)=rfieldb(i,k,j)+rfieldp(i,k,j)
3297 ENDDO
3298 ENDDO
3299 ENDDO
3301 END SUBROUTINE calculate_full
3303 !------------------------------------------------------------------------------
3305 SUBROUTINE coriolis ( ru, rv, rw, ru_tend, rv_tend, rw_tend, &
3306 config_flags, &
3307 msftx, msfty, msfux, msfuy, &
3308 msfvx, msfvy, &
3309 f, e, sina, cosa, fzm, fzp, &
3310 ids, ide, jds, jde, kds, kde, &
3311 ims, ime, jms, jme, kms, kme, &
3312 its, ite, jts, jte, kts, kte )
3314 IMPLICIT NONE
3316 ! Input data
3318 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
3320 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3321 ims, ime, jms, jme, kms, kme, &
3322 its, ite, jts, jte, kts, kte
3324 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: ru_tend, &
3325 rv_tend, &
3326 rw_tend
3327 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: ru, &
3328 rv, &
3329 rw
3331 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
3332 msfuy, &
3333 msfvx, &
3334 msfvy, &
3335 msftx, &
3336 msfty
3338 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: f, &
3339 e, &
3340 sina, &
3341 cosa
3343 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
3344 fzp
3346 ! Local indices.
3348 INTEGER :: i, j , k, ktf
3349 INTEGER :: i_start, i_end, j_start, j_end
3351 LOGICAL :: specified
3353 !<DESCRIPTION>
3354 !
3355 ! coriolis calculates the large timestep tendency terms in the
3356 ! u, v, and w momentum equations arise from the coriolis force.
3357 !
3358 !</DESCRIPTION>
3360 specified = .false.
3361 if(config_flags%specified .or. config_flags%nested) specified = .true.
3363 ktf=MIN(kte,kde-1)
3365 ! coriolis for u-momentum equation
3367 ! Notes on map scale factor
3368 ! cosa, sina are related to rotating the coordinate frame if desired
3369 ! generally sina=0, cosa=1
3370 ! ADT eqn 44, RHS terms 6 and 7: -2 mu w omega cos(lat)/my
3371 ! + 2 mu v omega sin(lat)/my
3372 ! Define f=2 omega sin(lat), e=2 omega cos(lat)
3373 ! => terms are: -e mu w / my + f mu v / my
3374 ! rv = mu v / mx ; rw = mu w / my
3375 ! => terms are: -e rw + f rv *mx / my
3377 i_start = its
3378 i_end = ite
3379 IF ( config_flags%open_xs .or. specified .or. &
3380 config_flags%nested) i_start = MAX(ids+1,its)
3381 IF ( config_flags%open_xe .or. specified .or. &
3382 config_flags%nested) i_end = MIN(ide-1,ite)
3383 IF ( config_flags%periodic_x ) i_start = its
3384 IF ( config_flags%periodic_x ) i_end = ite
3386 DO j = jts, MIN(jte,jde-1)
3388 DO k=kts,ktf
3389 DO i = i_start, i_end
3391 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)) &
3392 *0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
3393 - 0.5*(e(i,j)+e(i-1,j))*0.5*(cosa(i,j)+cosa(i-1,j)) &
3394 *0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
3396 ENDDO
3397 ENDDO
3399 IF ( (config_flags%open_xs) .and. (its == ids) ) THEN
3401 DO k=kts,ktf
3403 ru_tend(its,k,j)=ru_tend(its,k,j) + (msfux(its,j)/msfuy(its,j))*0.5*(f(its,j)+f(its,j)) &
3404 *0.25*(rv(its,k,j+1)+rv(its,k,j+1)+rv(its,k,j)+rv(its,k,j)) &
3405 - 0.5*(e(its,j)+e(its,j))*0.5*(cosa(its,j)+cosa(its,j)) &
3406 *0.25*(rw(its,k+1,j)+rw(its,k,j)+rw(its,k+1,j)+rw(its,k,j))
3408 ENDDO
3410 ENDIF
3412 IF ( (config_flags%open_xe) .and. (ite == ide) ) THEN
3414 DO k=kts,ktf
3416 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)) &
3417 *0.25*(rv(ite-1,k,j+1)+rv(ite-1,k,j+1)+rv(ite-1,k,j)+rv(ite-1,k,j)) &
3418 - 0.5*(e(ite-1,j)+e(ite-1,j))*0.5*(cosa(ite-1,j)+cosa(ite-1,j)) &
3419 *0.25*(rw(ite-1,k+1,j)+rw(ite-1,k,j)+rw(ite-1,k+1,j)+rw(ite-1,k,j))
3421 ENDDO
3423 ENDIF
3425 ENDDO
3427 ! coriolis term for v-momentum equation
3429 ! Notes on map scale factors
3430 ! ADT eqn 45, RHS terms 6 and 6b [0 for sina=0]: -2 mu u omega sin(lat)/mx + ?
3431 ! Define f=2 omega sin(lat), e=2 omega cos(lat)
3432 ! => terms are: -f mu u / mx
3433 ! ru = mu u / my ; rw = mu w / my
3434 ! => terms are: -f ru *my / mx + ?
3436 j_start = jts
3437 j_end = jte
3439 IF ( config_flags%open_ys .or. specified .or. &
3440 config_flags%nested .or. config_flags%polar) j_start = MAX(jds+1,jts)
3441 IF ( config_flags%open_ye .or. specified .or. &
3442 config_flags%nested .or. config_flags%polar) j_end = MIN(jde-1,jte)
3444 IF ( (config_flags%open_ys) .and. (jts == jds) ) THEN
3446 DO k=kts,ktf
3447 DO i=its,MIN(ide-1,ite)
3449 rv_tend(i,k,jts)=rv_tend(i,k,jts) - (msfvy(i,jts)/msfvx(i,jts))*0.5*(f(i,jts)+f(i,jts)) &
3450 *0.25*(ru(i,k,jts)+ru(i+1,k,jts)+ru(i,k,jts)+ru(i+1,k,jts)) &
3451 + (msfvy(i,jts)/msfvx(i,jts))*0.5*(e(i,jts)+e(i,jts))*0.5*(sina(i,jts)+sina(i,jts)) &
3452 *0.25*(rw(i,k+1,jts)+rw(i,k,jts)+rw(i,k+1,jts)+rw(i,k,jts))
3454 ENDDO
3455 ENDDO
3457 ENDIF
3459 DO j=j_start, j_end
3460 DO k=kts,ktf
3461 DO i=its,MIN(ide-1,ite)
3463 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)) &
3464 *0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
3465 + (msfvy(i,j)/msfvx(i,j))*0.5*(e(i,j)+e(i,j-1))*0.5*(sina(i,j)+sina(i,j-1)) &
3466 *0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
3468 ENDDO
3469 ENDDO
3470 ENDDO
3473 IF ( (config_flags%open_ye) .and. (jte == jde) ) THEN
3475 DO k=kts,ktf
3476 DO i=its,MIN(ide-1,ite)
3478 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)) &
3479 *0.25*(ru(i,k,jte-1)+ru(i+1,k,jte-1)+ru(i,k,jte-1)+ru(i+1,k,jte-1)) &
3480 + (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)) &
3481 *0.25*(rw(i,k+1,jte-1)+rw(i,k,jte-1)+rw(i,k+1,jte-1)+rw(i,k,jte-1))
3483 ENDDO
3484 ENDDO
3486 ENDIF
3488 ! coriolis term for w-mometum
3490 ! Notes on map scale factors
3491 ! ADT eqn 46/my, RHS terms 5 and 5b [0 for sina=0]: 2 mu u omega cos(lat)/my +?
3492 ! Define e=2 omega cos(lat)
3493 ! => terms are: e mu u / my + ???
3494 ! ru = mu u / my ; ru = mu v / mx
3495 ! => terms are: e ru + ???
3497 DO j=jts,MIN(jte, jde-1)
3498 DO k=kts+1,ktf
3499 DO i=its,MIN(ite, ide-1)
3501 rw_tend(i,k,j)=rw_tend(i,k,j) + e(i,j)* &
3502 (cosa(i,j)*0.5*(fzm(k)*(ru(i,k,j)+ru(i+1,k,j)) &
3503 +fzp(k)*(ru(i,k-1,j)+ru(i+1,k-1,j))) &
3504 -(msftx(i,j)/msfty(i,j))* &
3505 sina(i,j)*0.5*(fzm(k)*(rv(i,k,j)+rv(i,k,j+1)) &
3506 +fzp(k)*(rv(i,k-1,j)+rv(i,k-1,j+1))))
3508 ENDDO
3509 ENDDO
3510 ENDDO
3512 END SUBROUTINE coriolis
3514 !------------------------------------------------------------------------------
3516 SUBROUTINE perturbation_coriolis ( ru_in, rv_in, rw, ru_tend, rv_tend, rw_tend, &
3517 config_flags, &
3518 u_base, v_base, z_base, &
3519 muu, muv, phb, ph, &
3520 msftx, msfty, msfux, msfuy, msfvx, msfvy, &
3521 f, e, sina, cosa, fzm, fzp, &
3522 ids, ide, jds, jde, kds, kde, &
3523 ims, ime, jms, jme, kms, kme, &
3524 its, ite, jts, jte, kts, kte )
3526 IMPLICIT NONE
3528 ! Input data
3530 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
3532 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3533 ims, ime, jms, jme, kms, kme, &
3534 its, ite, jts, jte, kts, kte
3536 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: ru_tend, &
3537 rv_tend, &
3538 rw_tend
3539 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: ru_in, &
3540 rv_in, &
3541 rw, &
3542 ph, &
3543 phb
3546 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
3547 msfuy, &
3548 msfvx, &
3549 msfvy, &
3550 msftx, &
3551 msfty
3553 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: f, &
3554 e, &
3555 sina, &
3556 cosa
3558 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muu, &
3559 muv
3562 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
3563 fzp
3565 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: u_base, &
3566 v_base, &
3567 z_base
3569 ! Local storage
3571 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) :: ru, &
3572 rv
3574 REAL :: z_at_u, z_at_v, wkp1, wk, wkm1
3576 ! Local indices.
3578 INTEGER :: i, j , k, ktf
3579 INTEGER :: i_start, i_end, j_start, j_end
3581 LOGICAL :: specified
3583 !<DESCRIPTION>
3584 !
3585 ! perturbation_coriolis calculates the large timestep tendency terms in the
3586 ! u, v, and w momentum equations arise from the coriolis force. This version
3587 ! subtracts off the horizontal velocities from the initial sounding when
3588 ! computing the forcing terms, hence "perturbation" coriolis.
3589 !
3590 !</DESCRIPTION>
3592 specified = .false.
3593 if(config_flags%specified .or. config_flags%nested) specified = .true.
3595 ktf=MIN(kte,kde-1)
3597 ! coriolis for u-momentum equation
3599 i_start = its
3600 i_end = ite
3601 IF ( config_flags%open_xs .or. specified .or. &
3602 config_flags%nested) i_start = MAX(ids+1,its)
3603 IF ( config_flags%open_xe .or. specified .or. &
3604 config_flags%nested) i_end = MIN(ide-1,ite)
3605 IF ( config_flags%periodic_x ) i_start = its
3606 IF ( config_flags%periodic_x ) i_end = ite
3608 ! compute perturbation mu*v for use in u momentum equation
3610 DO j = jts, MIN(jte,jde-1)+1
3611 DO k=kts+1,ktf-1
3612 DO i = i_start-1, i_end
3613 z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
3614 +phb(i,k,j-1)+phb(i,k+1,j-1) &
3615 +ph(i,k,j )+ph(i,k+1,j ) &
3616 +ph(i,k,j-1)+ph(i,k+1,j-1))/g
3617 wkp1 = min(1.,max(0.,z_at_v-z_base(k))/(z_base(k+1)-z_base(k)))
3618 wkm1 = min(1.,max(0.,z_base(k)-z_at_v)/(z_base(k)-z_base(k-1)))
3619 wk = 1.-wkp1-wkm1
3620 rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*( &
3621 wkm1*v_base(k-1) &
3622 +wk *v_base(k ) &
3623 +wkp1*v_base(k+1) )
3624 ENDDO
3625 ENDDO
3626 ENDDO
3629 ! pick up top and bottom v
3631 DO j = jts, MIN(jte,jde-1)+1
3632 DO i = i_start-1, i_end
3634 k = kts
3635 z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
3636 +phb(i,k,j-1)+phb(i,k+1,j-1) &
3637 +ph(i,k,j )+ph(i,k+1,j ) &
3638 +ph(i,k,j-1)+ph(i,k+1,j-1))/g
3639 wkp1 = min(1.,max(0.,z_at_v-z_base(k))/(z_base(k+1)-z_base(k)))
3640 wk = 1.-wkp1
3641 rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*( &
3642 +wk *v_base(k ) &
3643 +wkp1*v_base(k+1) )
3645 k = ktf
3646 z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
3647 +phb(i,k,j-1)+phb(i,k+1,j-1) &
3648 +ph(i,k,j )+ph(i,k+1,j ) &
3649 +ph(i,k,j-1)+ph(i,k+1,j-1))/g
3650 wkm1 = min(1.,max(0.,z_base(k)-z_at_v)/(z_base(k)-z_base(k-1)))
3651 wk = 1.-wkm1
3652 rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*( &
3653 wkm1*v_base(k-1) &
3654 +wk *v_base(k ) )
3656 ENDDO
3657 ENDDO
3659 ! compute coriolis forcing for u
3661 ! Map scale factors: see comments above for Coriolis
3663 DO j = jts, MIN(jte,jde-1)
3665 DO k=kts,ktf
3666 DO i = i_start, i_end
3667 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)) &
3668 *0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
3669 - 0.5*(e(i,j)+e(i-1,j))*0.5*(cosa(i,j)+cosa(i-1,j)) &
3670 *0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
3671 ENDDO
3672 ENDDO
3674 IF ( (config_flags%open_xs) .and. (its == ids) ) THEN
3676 DO k=kts,ktf
3678 ru_tend(its,k,j)=ru_tend(its,k,j) + (msfux(its,j)/msfuy(its,j))*0.5*(f(its,j)+f(its,j)) &
3679 *0.25*(rv(its,k,j+1)+rv(its,k,j+1)+rv(its,k,j)+rv(its,k,j)) &
3680 - 0.5*(e(its,j)+e(its,j))*0.5*(cosa(its,j)+cosa(its,j)) &
3681 *0.25*(rw(its,k+1,j)+rw(its,k,j)+rw(its,k+1,j)+rw(its,k,j))
3683 ENDDO
3685 ENDIF
3687 IF ( (config_flags%open_xe) .and. (ite == ide) ) THEN
3689 DO k=kts,ktf
3691 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)) &
3692 *0.25*(rv(ite-1,k,j+1)+rv(ite-1,k,j+1)+rv(ite-1,k,j)+rv(ite-1,k,j)) &
3693 - 0.5*(e(ite-1,j)+e(ite-1,j))*0.5*(cosa(ite-1,j)+cosa(ite-1,j)) &
3694 *0.25*(rw(ite-1,k+1,j)+rw(ite-1,k,j)+rw(ite-1,k+1,j)+rw(ite-1,k,j))
3696 ENDDO
3698 ENDIF
3700 ENDDO
3702 ! coriolis term for v-momentum equation
3703 ! Map scale factors: see comments above for Coriolis
3705 j_start = jts
3706 j_end = jte
3708 IF ( config_flags%open_ys .or. specified .or. &
3709 config_flags%nested .or. config_flags%polar) j_start = MAX(jds+1,jts)
3710 IF ( config_flags%open_ye .or. specified .or. &
3711 config_flags%nested .or. config_flags%polar) j_end = MIN(jde-1,jte)
3713 ! compute perturbation mu*u for use in v momentum equation
3715 DO j = j_start-1,j_end
3716 DO k=kts+1,ktf-1
3717 DO i = its, MIN(ite,ide-1)+1
3718 z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
3719 +phb(i-1,k,j)+phb(i-1,k+1,j) &
3720 +ph(i ,k,j)+ph(i ,k+1,j) &
3721 +ph(i-1,k,j)+ph(i-1,k+1,j))/g
3722 wkp1 = min(1.,max(0.,z_at_u-z_base(k))/(z_base(k+1)-z_base(k)))
3723 wkm1 = min(1.,max(0.,z_base(k)-z_at_u)/(z_base(k)-z_base(k-1)))
3724 wk = 1.-wkp1-wkm1
3725 ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*( &
3726 wkm1*u_base(k-1) &
3727 +wk *u_base(k ) &
3728 +wkp1*u_base(k+1) )
3729 ENDDO
3730 ENDDO
3731 ENDDO
3733 ! pick up top and bottom u
3735 DO j = j_start-1,j_end
3736 DO i = its, MIN(ite,ide-1)+1
3738 k = kts
3739 z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
3740 +phb(i-1,k,j)+phb(i-1,k+1,j) &
3741 +ph(i ,k,j)+ph(i ,k+1,j) &
3742 +ph(i-1,k,j)+ph(i-1,k+1,j))/g
3743 wkp1 = min(1.,max(0.,z_at_u-z_base(k))/(z_base(k+1)-z_base(k)))
3744 wk = 1.-wkp1
3745 ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*( &
3746 +wk *u_base(k ) &
3747 +wkp1*u_base(k+1) )
3750 k = ktf
3751 z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
3752 +phb(i-1,k,j)+phb(i-1,k+1,j) &
3753 +ph(i ,k,j)+ph(i ,k+1,j) &
3754 +ph(i-1,k,j)+ph(i-1,k+1,j))/g
3755 wkm1 = min(1.,max(0.,z_base(k)-z_at_u)/(z_base(k)-z_base(k-1)))
3756 wk = 1.-wkm1
3757 ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*( &
3758 wkm1*u_base(k-1) &
3759 +wk *u_base(k ) )
3761 ENDDO
3762 ENDDO
3764 ! compute coriolis forcing for v momentum equation
3765 ! Map scale factors: see comments above for Coriolis
3767 IF ( (config_flags%open_ys) .and. (jts == jds) ) THEN
3769 DO k=kts,ktf
3770 DO i=its,MIN(ide-1,ite)
3772 rv_tend(i,k,jts)=rv_tend(i,k,jts) - (msfvy(i,jts)/msfvx(i,jts))*0.5*(f(i,jts)+f(i,jts)) &
3773 *0.25*(ru(i,k,jts)+ru(i+1,k,jts)+ru(i,k,jts)+ru(i+1,k,jts)) &
3774 + (msfvy(i,jts)/msfvx(i,jts))*0.5*(e(i,jts)+e(i,jts))*0.5*(sina(i,jts)+sina(i,jts)) &
3775 *0.25*(rw(i,k+1,jts)+rw(i,k,jts)+rw(i,k+1,jts)+rw(i,k,jts))
3777 ENDDO
3778 ENDDO
3780 ENDIF
3782 DO j=j_start, j_end
3783 DO k=kts,ktf
3784 DO i=its,MIN(ide-1,ite)
3786 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)) &
3787 *0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
3788 + (msfvy(i,j)/msfvx(i,j))*0.5*(e(i,j)+e(i,j-1))*0.5*(sina(i,j)+sina(i,j-1)) &
3789 *0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
3791 ENDDO
3792 ENDDO
3793 ENDDO
3796 IF ( (config_flags%open_ye) .and. (jte == jde) ) THEN
3798 DO k=kts,ktf
3799 DO i=its,MIN(ide-1,ite)
3801 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)) &
3802 *0.25*(ru(i,k,jte-1)+ru(i+1,k,jte-1)+ru(i,k,jte-1)+ru(i+1,k,jte-1)) &
3803 + (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)) &
3804 *0.25*(rw(i,k+1,jte-1)+rw(i,k,jte-1)+rw(i,k+1,jte-1)+rw(i,k,jte-1))
3806 ENDDO
3807 ENDDO
3809 ENDIF
3811 ! coriolis term for w-mometum
3812 ! Map scale factors: see comments above for Coriolis
3814 DO j=jts,MIN(jte, jde-1)
3815 DO k=kts+1,ktf
3816 DO i=its,MIN(ite, ide-1)
3818 rw_tend(i,k,j)=rw_tend(i,k,j) + e(i,j)* &
3819 (cosa(i,j)*0.5*(fzm(k)*(ru(i,k,j)+ru(i+1,k,j)) &
3820 +fzp(k)*(ru(i,k-1,j)+ru(i+1,k-1,j))) &
3821 -(msftx(i,j)/msfty(i,j))*sina(i,j)*0.5*(fzm(k)*(rv(i,k,j)+rv(i,k,j+1)) &
3822 +fzp(k)*(rv(i,k-1,j)+rv(i,k-1,j+1))))
3824 ENDDO
3825 ENDDO
3826 ENDDO
3828 END SUBROUTINE perturbation_coriolis
3830 !------------------------------------------------------------------------------
3832 SUBROUTINE curvature ( ru, rv, rw, u, v, w, ru_tend, rv_tend, rw_tend, &
3833 config_flags, &
3834 msfux, msfuy, msfvx, msfvy, msftx, msfty, &
3835 xlat, fzm, fzp, rdx, rdy, &
3836 ids, ide, jds, jde, kds, kde, &
3837 ims, ime, jms, jme, kms, kme, &
3838 its, ite, jts, jte, kts, kte )
3841 IMPLICIT NONE
3843 ! Input data
3845 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
3847 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3848 ims, ime, jms, jme, kms, kme, &
3849 its, ite, jts, jte, kts, kte
3851 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3852 INTENT(INOUT) :: ru_tend, &
3853 rv_tend, &
3854 rw_tend
3856 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3857 INTENT(IN ) :: ru, &
3858 rv, &
3859 rw, &
3860 u, &
3861 v, &
3862 w
3864 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
3865 msfuy, &
3866 msfvx, &
3867 msfvy, &
3868 msftx, &
3869 msfty, &
3870 xlat
3872 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
3873 fzp
3875 REAL , INTENT(IN ) :: rdx, &
3876 rdy
3878 ! Local data
3880 ! INTEGER :: i, j, k, itf, jtf, ktf, kp1, im, ip, jm, jp
3881 INTEGER :: i, j, k, itf, jtf, ktf
3882 INTEGER :: i_start, i_end, j_start, j_end
3883 ! INTEGER :: irmin, irmax, jrmin, jrmax
3885 REAL , DIMENSION( its-1:ite , kts:kte, jts-1:jte ) :: vxgm
3887 LOGICAL :: specified
3889 !<DESCRIPTION>
3890 !
3891 ! curvature calculates the large timestep tendency terms in the
3892 ! u, v, and w momentum equations arise from the curvature terms.
3893 !
3894 !</DESCRIPTION>
3896 specified = .false.
3897 if(config_flags%specified .or. config_flags%nested) specified = .true.
3899 itf=MIN(ite,ide-1)
3900 jtf=MIN(jte,jde-1)
3901 ktf=MIN(kte,kde-1)
3903 ! irmin = ims
3904 ! irmax = ime
3905 ! jrmin = jms
3906 ! jrmax = jme
3907 ! IF ( config_flags%open_xs ) irmin = ids
3908 ! IF ( config_flags%open_xe ) irmax = ide-1
3909 ! IF ( config_flags%open_ys ) jrmin = jds
3910 ! IF ( config_flags%open_ye ) jrmax = jde-1
3912 ! Define v cross grad m at scalar points - vxgm(i,j)
3914 i_start = its-1
3915 i_end = ite
3916 j_start = jts-1
3917 j_end = jte
3919 IF ( ( config_flags%open_xs .or. specified .or. &
3920 config_flags%nested) .and. (its == ids) ) i_start = its
3921 IF ( ( config_flags%open_xe .or. specified .or. &
3922 config_flags%nested) .and. (ite == ide) ) i_end = ite-1
3923 IF ( ( config_flags%open_ys .or. specified .or. &
3924 config_flags%nested .or. config_flags%polar) .and. (jts == jds) ) j_start = jts
3925 IF ( ( config_flags%open_ye .or. specified .or. &
3926 config_flags%nested .or. config_flags%polar) .and. (jte == jde) ) j_end = jte-1
3927 IF ( config_flags%periodic_x ) i_start = its-1
3928 IF ( config_flags%periodic_x ) i_end = ite
3930 DO j=j_start, j_end
3931 DO k=kts,ktf
3932 DO i=i_start, i_end
3933 ! Map scale factor notes:
3934 ! msf...y is constant everywhere for cylindrical map projection
3935 ! msf...x varies with y only
3936 ! But we know that this is not = 0 for cylindrical,
3937 ! therefore use msfvX in 1st line
3938 ! which => by symmetry use msfuY in 2nd line - ???
3939 vxgm(i,k,j)=0.5*(u(i,k,j)+u(i+1,k,j))*(msfvx(i,j+1)-msfvx(i,j))*rdy - &
3940 0.5*(v(i,k,j)+v(i,k,j+1))*(msfuy(i+1,j)-msfuy(i,j))*rdx
3941 ENDDO
3942 ENDDO
3943 ENDDO
3945 ! Pick up the boundary rows for open (radiation) lateral b.c.
3946 ! Rather crude at present, we are assuming there is no
3947 ! variation in this term at the boundary.
3949 IF ( ( config_flags%open_xs .or. (specified .AND. .NOT. config_flags%periodic_x) .or. &
3950 config_flags%nested) .and. (its == ids) ) THEN
3952 DO j = jts, jte-1
3953 DO k = kts, ktf
3954 vxgm(its-1,k,j) = vxgm(its,k,j)
3955 ENDDO
3956 ENDDO
3958 ENDIF
3960 IF ( ( config_flags%open_xe .or. (specified .AND. .NOT. config_flags%periodic_x) .or. &
3961 config_flags%nested) .and. (ite == ide) ) THEN
3963 DO j = jts, jte-1
3964 DO k = kts, ktf
3965 vxgm(ite,k,j) = vxgm(ite-1,k,j)
3966 ENDDO
3967 ENDDO
3969 ENDIF
3971 ! Polar boundary condition:
3972 ! The following change is needed in case one tries using the vxgm route with
3973 ! polar B.C.'s in the future, but not needed if 'tan' used
3974 IF ( ( config_flags%open_ys .or. specified .or. &
3975 config_flags%nested .or. config_flags%polar) .and. (jts == jds) ) THEN
3977 DO k = kts, ktf
3978 DO i = its-1, ite
3979 vxgm(i,k,jts-1) = vxgm(i,k,jts)
3980 ENDDO
3981 ENDDO
3983 ENDIF
3985 ! Polar boundary condition:
3986 ! The following change is needed in case one tries using the vxgm route with
3987 ! polar B.C.'s in the future, but not needed if 'tan' used
3988 IF ( ( config_flags%open_ye .or. specified .or. &
3989 config_flags%nested .or. config_flags%polar) .and. (jte == jde) ) THEN
3991 DO k = kts, ktf
3992 DO i = its-1, ite
3993 vxgm(i,k,jte) = vxgm(i,k,jte-1)
3994 ENDDO
3995 ENDDO
3997 ENDIF
3999 ! curvature term for u momentum eqn.
4001 ! Map scale factor notes:
4002 ! ADT eqn 44, RHS terms 4 and 5, in cylindrical: mu u v tan(lat)/(a my)
4003 ! - mu u w /(a my)
4004 ! ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
4005 ! => terms are:
4006 ! (mx/my)*u rv tan(lat) / a - u rw / a = (u/a)*[(mx/my) rv tan(lat) - rw]
4007 ! ru v tan(lat) / a - u rw / a
4008 ! xlat defined with end points half grid space from pole,
4009 ! hence are on u latitude points
4011 i_start = its
4012 IF ( config_flags%open_xs .or. specified .or. &
4013 config_flags%nested) i_start = MAX ( ids+1 , its )
4014 IF ( config_flags%open_xe .or. specified .or. &
4015 config_flags%nested) i_end = MIN ( ide-1 , ite )
4016 IF ( config_flags%periodic_x ) i_start = its
4017 IF ( config_flags%periodic_x ) i_end = ite
4019 ! Polar boundary condition
4020 IF ((config_flags%map_proj == 6) .OR. (config_flags%polar)) THEN
4022 DO j=jts,MIN(jde-1,jte)
4023 DO k=kts,ktf
4024 DO i=i_start,i_end
4026 ru_tend(i,k,j)=ru_tend(i,k,j) + u(i,k,j)*reradius* ( &
4027 (msfux(i,j)/msfuy(i,j))*0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+ &
4028 rv(i-1,k,j)+rv(i,k,j))*tan(xlat(i,j)*degrad) &
4029 - 0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j)) )
4030 ENDDO
4031 ENDDO
4032 ENDDO
4034 ELSE ! normal code
4037 DO j=jts,MIN(jde-1,jte)
4038 DO k=kts,ktf
4039 DO i=i_start,i_end
4041 ru_tend(i,k,j)=ru_tend(i,k,j) + 0.5*(vxgm(i,k,j)+vxgm(i-1,k,j)) &
4042 *0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
4043 - u(i,k,j)*reradius &
4044 *0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
4046 ENDDO
4047 ENDDO
4048 ENDDO
4050 END IF
4052 ! curvature term for v momentum eqn.
4054 ! Map scale factor notes
4055 ! ADT eqn 45, RHS terms 4 and 5, in cylindrical: mu u*u tan(lat)/(a mx)
4056 ! - mu v w /(a mx)
4057 ! ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
4058 ! terms are:
4059 ! (my/mx)*u ru tan(lat) / a - (my/mx)*v rw / a
4060 ! = [my/(mx*a)]*[u ru tan(lat) - v rw]
4061 ! (1/a)*[(my/mx)*u ru tan(lat) - w rv]
4062 ! xlat defined with end points half grid space from pole, hence are on
4063 ! u latitude points => av here
4064 !
4065 ! in original wrf, there was a sign error for the rw contribution
4067 j_start = jts
4068 IF ( config_flags%open_ys .or. specified .or. &
4069 config_flags%nested .or. config_flags%polar) j_start = MAX ( jds+1 , jts )
4070 IF ( config_flags%open_ye .or. specified .or. &
4071 config_flags%nested .or. config_flags%polar) j_end = MIN ( jde-1 , jte )
4073 IF ((config_flags%map_proj == 6) .OR. (config_flags%polar)) THEN
4075 DO j=j_start,j_end
4076 DO k=kts,ktf
4077 DO i=its,MIN(ite,ide-1)
4078 rv_tend(i,k,j)=rv_tend(i,k,j) - (msfvy(i,j)/msfvx(i,j))*reradius* ( &
4079 0.25*(u(i,k,j)+u(i+1,k,j)+u(i,k,j-1)+u(i+1,k,j-1))* &
4080 tan((xlat(i,j)+xlat(i,j-1))*0.5*degrad)* &
4081 0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
4082 - v(i,k,j)*0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+ &
4083 rw(i,k+1,j)+rw(i,k,j)) )
4084 ENDDO
4085 ENDDO
4086 ENDDO
4088 ELSE ! normal code
4090 DO j=j_start,j_end
4091 DO k=kts,ktf
4092 DO i=its,MIN(ite,ide-1)
4094 rv_tend(i,k,j)=rv_tend(i,k,j) - 0.5*(vxgm(i,k,j)+vxgm(i,k,j-1)) &
4095 *0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
4096 - (msfvy(i,j)/msfvx(i,j))*v(i,k,j)*reradius &
4097 *0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
4099 ENDDO
4100 ENDDO
4101 ENDDO
4103 END IF
4105 ! curvature term for vertical momentum eqn.
4107 ! Notes on map scale factors:
4108 ! ADT eqn 46, RHS term 4: [mu/(a my)]*[u*u + v*v]
4109 ! ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
4110 ! terms are: u ru / a + (mx/my)v rv / a
4112 DO j=jts,MIN(jte,jde-1)
4113 DO k=MAX(2,kts),ktf
4114 DO i=its,MIN(ite,ide-1)
4116 rw_tend(i,k,j)=rw_tend(i,k,j) + reradius* &
4117 (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))) &
4118 *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))) &
4119 +(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))) &
4120 *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))))
4122 ENDDO
4123 ENDDO
4124 ENDDO
4126 END SUBROUTINE curvature
4128 !------------------------------------------------------------------------------
4130 SUBROUTINE decouple ( rr, rfield, field, name, config_flags, &
4131 fzm, fzp, &
4132 ids, ide, jds, jde, kds, kde, &
4133 ims, ime, jms, jme, kms, kme, &
4134 its, ite, jts, jte, kts, kte )
4136 IMPLICIT NONE
4138 ! Input data
4140 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
4142 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4143 ims, ime, jms, jme, kms, kme, &
4144 its, ite, jts, jte, kts, kte
4146 CHARACTER(LEN=1) , INTENT(IN ) :: name
4148 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rfield
4150 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rr
4152 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: field
4154 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, fzp
4156 ! Local data
4158 INTEGER :: i, j, k, itf, jtf, ktf
4160 !<DESCRIPTION>
4161 !
4162 ! decouple decouples a variable from the column dry-air mass.
4163 !
4164 !</DESCRIPTION>
4166 ktf=MIN(kte,kde-1)
4168 IF (name .EQ. 'u')THEN
4169 itf=ite
4170 jtf=MIN(jte,jde-1)
4172 DO j=jts,jtf
4173 DO k=kts,ktf
4174 DO i=its,itf
4175 field(i,k,j)=rfield(i,k,j)/(0.5*(rr(i,k,j)+rr(i-1,k,j)))
4176 ENDDO
4177 ENDDO
4178 ENDDO
4180 ELSE IF (name .EQ. 'v')THEN
4181 itf=MIN(ite,ide-1)
4182 jtf=jte
4184 DO j=jts,jtf
4185 DO k=kts,ktf
4186 DO i=its,itf
4187 field(i,k,j)=rfield(i,k,j)/(0.5*(rr(i,k,j)+rr(i,k,j-1)))
4188 ENDDO
4189 ENDDO
4190 ENDDO
4192 ELSE IF (name .EQ. 'w')THEN
4193 itf=MIN(ite,ide-1)
4194 jtf=MIN(jte,jde-1)
4195 DO j=jts,jtf
4196 DO k=kts+1,ktf
4197 DO i=its,itf
4198 field(i,k,j)=rfield(i,k,j)/(fzm(k)*rr(i,k,j)+fzp(k)*rr(i,k-1,j))
4199 ENDDO
4200 ENDDO
4201 ENDDO
4203 DO j=jts,jtf
4204 DO i=its,itf
4205 field(i,kte,j) = 0.
4206 ENDDO
4207 ENDDO
4209 ELSE
4210 itf=MIN(ite,ide-1)
4211 jtf=MIN(jte,jde-1)
4212 ! For theta we will decouple tb and tp and add them to give t afterwards
4213 DO j=jts,jtf
4214 DO k=kts,ktf
4215 DO i=its,itf
4216 field(i,k,j)=rfield(i,k,j)/rr(i,k,j)
4217 ENDDO
4218 ENDDO
4219 ENDDO
4221 ENDIF
4223 END SUBROUTINE decouple
4225 !-------------------------------------------------------------------------------
4228 SUBROUTINE zero_tend ( tendency, &
4229 ids, ide, jds, jde, kds, kde, &
4230 ims, ime, jms, jme, kms, kme, &
4231 its, ite, jts, jte, kts, kte )
4234 IMPLICIT NONE
4236 ! Input data
4238 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4239 ims, ime, jms, jme, kms, kme, &
4240 its, ite, jts, jte, kts, kte
4242 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
4244 ! Local data
4246 INTEGER :: i, j, k, itf, jtf, ktf
4248 !<DESCRIPTION>
4249 !
4250 ! zero_tend sets the input tendency array to zero.
4251 !
4252 !</DESCRIPTION>
4254 DO j = jts, jte
4255 DO k = kts, kte
4256 DO i = its, ite
4257 tendency(i,k,j) = 0.
4258 ENDDO
4259 ENDDO
4260 ENDDO
4262 END SUBROUTINE zero_tend
4264 !-------------------------------------------------------------------------------
4265 ! Sets the an array on the polar v point(s) to zero
4266 SUBROUTINE zero_pole ( field, &
4267 ids, ide, jds, jde, kds, kde, &
4268 ims, ime, jms, jme, kms, kme, &
4269 its, ite, jts, jte, kts, kte )
4272 IMPLICIT NONE
4274 ! Input data
4276 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4277 ims, ime, jms, jme, kms, kme, &
4278 its, ite, jts, jte, kts, kte
4280 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: field
4282 ! Local data
4284 INTEGER :: i, k
4286 IF (jts == jds) THEN
4287 DO k = kts, kte
4288 DO i = its-1, ite+1
4289 field(i,k,jts) = 0.
4290 END DO
4291 END DO
4292 END IF
4293 IF (jte == jde) THEN
4294 DO k = kts, kte
4295 DO i = its-1, ite+1
4296 field(i,k,jte) = 0.
4297 END DO
4298 END DO
4299 END IF
4301 END SUBROUTINE zero_pole
4303 !-------------------------------------------------------------------------------
4304 ! Sets the an array on the polar v point(s)
4305 SUBROUTINE pole_point_bc ( field, &
4306 ids, ide, jds, jde, kds, kde, &
4307 ims, ime, jms, jme, kms, kme, &
4308 its, ite, jts, jte, kts, kte )
4311 IMPLICIT NONE
4313 ! Input data
4315 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4316 ims, ime, jms, jme, kms, kme, &
4317 its, ite, jts, jte, kts, kte
4319 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: field
4321 ! Local data
4323 INTEGER :: i, k
4325 IF (jts == jds) THEN
4326 DO k = kts, kte
4327 DO i = its, ite
4328 ! field(i,k,jts) = 2*field(i,k,jts+1) - field(i,k,jts+2)
4329 field(i,k,jts) = field(i,k,jts+1)
4330 END DO
4331 END DO
4332 END IF
4333 IF (jte == jde) THEN
4334 DO k = kts, kte
4335 DO i = its, ite
4336 ! field(i,k,jte) = 2*field(i,k,jte-1) - field(i,k,jte-2)
4337 field(i,k,jte) = field(i,k,jte-1)
4338 END DO
4339 END DO
4340 END IF
4342 END SUBROUTINE pole_point_bc
4344 !======================================================================
4345 ! physics prep routines
4346 !======================================================================
4348 SUBROUTINE phy_prep ( config_flags, & ! input
4349 mu, muu, muv, u, v, p, pb, alt, ph, & ! input
4350 phb, t, tsk, moist, n_moist, & ! input
4351 mu_3d, rho, th_phy, p_phy , pi_phy , & ! output
4352 u_phy, v_phy, p8w, t_phy, t8w, & ! output
4353 z, z_at_w, dz8w, & ! output
4354 fzm, fzp, & ! params
4355 RTHRATEN, &
4356 RTHBLTEN, RUBLTEN, RVBLTEN, &
4357 RQVBLTEN, RQCBLTEN, RQIBLTEN, &
4358 RTHCUTEN, RQVCUTEN, RQCCUTEN, &
4359 RQRCUTEN, RQICUTEN, RQSCUTEN, &
4360 RTHFTEN, RQVFTEN, &
4361 RUNDGDTEN, RVNDGDTEN, RTHNDGDTEN, &
4362 RQVNDGDTEN, RMUNDGDTEN, &
4363 ids, ide, jds, jde, kds, kde, &
4364 ims, ime, jms, jme, kms, kme, &
4365 its, ite, jts, jte, kts, kte )
4366 !----------------------------------------------------------------------
4367 IMPLICIT NONE
4368 !----------------------------------------------------------------------
4370 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
4372 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4373 ims, ime, jms, jme, kms, kme, &
4374 its, ite, jts, jte, kts, kte
4375 INTEGER , INTENT(IN ) :: n_moist
4377 REAL, DIMENSION( ims:ime, kms:kme , jms:jme , n_moist ), INTENT(IN) :: moist
4380 REAL , DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: TSK, mu, muu, muv
4382 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4383 INTENT( OUT) :: u_phy, &
4384 v_phy, &
4385 pi_phy, &
4386 p_phy, &
4387 p8w, &
4388 t_phy, &
4389 th_phy, &
4390 t8w, &
4391 mu_3d, &
4392 rho, &
4393 z, &
4394 dz8w, &
4395 z_at_w
4397 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4398 INTENT(IN ) :: pb, &
4399 p, &
4400 u, &
4401 v, &
4402 alt, &
4403 ph, &
4404 phb, &
4405 t
4408 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
4409 fzp
4411 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4412 INTENT(INOUT) :: RTHRATEN
4414 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4415 INTENT(INOUT) :: RTHCUTEN, &
4416 RQVCUTEN, &
4417 RQCCUTEN, &
4418 RQRCUTEN, &
4419 RQICUTEN, &
4420 RQSCUTEN
4422 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
4423 INTENT(INOUT) :: RUBLTEN, &
4424 RVBLTEN, &
4425 RTHBLTEN, &
4426 RQVBLTEN, &
4427 RQCBLTEN, &
4428 RQIBLTEN
4430 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
4431 INTENT(INOUT) :: RTHFTEN, &
4432 RQVFTEN
4434 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
4435 INTENT(INOUT) :: RUNDGDTEN, &
4436 RVNDGDTEN, &
4437 RTHNDGDTEN, &
4438 RQVNDGDTEN, &
4439 RMUNDGDTEN
4441 INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end, i_startu, j_startv
4442 INTEGER :: i, j, k
4443 REAL :: w1, w2, z0, z1, z2
4445 !-----------------------------------------------------------------------
4447 !<DESCRIPTION>
4448 !
4449 ! phys_prep calculates a number of diagnostic quantities needed by
4450 ! the physics routines. It also decouples the physics tendencies from
4451 ! the column dry-air mass (the physics routines expect to see/update the
4452 ! uncoupled tendencies).
4453 !
4454 !</DESCRIPTION>
4456 ! set up loop bounds for this grid's boundary conditions
4458 i_start = its
4459 i_end = min( ite,ide-1 )
4460 j_start = jts
4461 j_end = min( jte,jde-1 )
4463 k_start = kts
4464 k_end = min( kte, kde-1 )
4466 ! compute thermodynamics and velocities at pressure points
4468 do j = j_start,j_end
4469 do k = k_start, k_end
4470 do i = i_start, i_end
4472 th_phy(i,k,j) = t(i,k,j) + t0
4473 p_phy(i,k,j) = p(i,k,j) + pb(i,k,j)
4474 pi_phy(i,k,j) = (p_phy(i,k,j)/p1000mb)**rcp
4475 t_phy(i,k,j) = th_phy(i,k,j)*pi_phy(i,k,j)
4476 rho(i,k,j) = 1./alt(i,k,j)*(1.+moist(i,k,j,P_QV))
4477 mu_3d(i,k,j) = mu(i,j)
4478 u_phy(i,k,j) = 0.5*(u(i,k,j)+u(i+1,k,j))
4479 v_phy(i,k,j) = 0.5*(v(i,k,j)+v(i,k,j+1))
4481 enddo
4482 enddo
4483 enddo
4485 ! compute z at w points
4487 do j = j_start,j_end
4488 do k = k_start, kte
4489 do i = i_start, i_end
4490 z_at_w(i,k,j) = (phb(i,k,j)+ph(i,k,j))/g
4491 enddo
4492 enddo
4493 enddo
4495 do j = j_start,j_end
4496 do k = k_start, kte-1
4497 do i = i_start, i_end
4498 dz8w(i,k,j) = z_at_w(i,k+1,j)-z_at_w(i,k,j)
4499 enddo
4500 enddo
4501 enddo
4503 do j = j_start,j_end
4504 do i = i_start, i_end
4505 dz8w(i,kte,j) = 0.
4506 enddo
4507 enddo
4509 ! compute z at p points (average of z at w points)
4511 do j = j_start,j_end
4512 do k = k_start, k_end
4513 do i = i_start, i_end
4514 z(i,k,j) = 0.5*(z_at_w(i,k,j) +z_at_w(i,k+1,j) )
4515 enddo
4516 enddo
4517 enddo
4519 ! interp t and p at w points
4521 do j = j_start,j_end
4522 do k = 2, k_end
4523 do i = i_start, i_end
4524 p8w(i,k,j) = fzm(k)*p_phy(i,k,j)+fzp(k)*p_phy(i,k-1,j)
4525 t8w(i,k,j) = fzm(k)*t_phy(i,k,j)+fzp(k)*t_phy(i,k-1,j)
4526 enddo
4527 enddo
4528 enddo
4530 ! extrapolate p and t to surface and top.
4531 ! we'll use an extrapolation in z for now
4533 do j = j_start,j_end
4534 do i = i_start, i_end
4536 ! bottom
4538 z0 = z_at_w(i,1,j)
4539 z1 = z(i,1,j)
4540 z2 = z(i,2,j)
4541 w1 = (z0 - z2)/(z1 - z2)
4542 w2 = 1. - w1
4543 p8w(i,1,j) = w1*p_phy(i,1,j)+w2*p_phy(i,2,j)
4544 t8w(i,1,j) = w1*t_phy(i,1,j)+w2*t_phy(i,2,j)
4546 ! top
4548 z0 = z_at_w(i,kte,j)
4549 z1 = z(i,k_end,j)
4550 z2 = z(i,k_end-1,j)
4551 w1 = (z0 - z2)/(z1 - z2)
4552 w2 = 1. - w1
4554 ! p8w(i,kde,j) = w1*p_phy(i,kde-1,j)+w2*p_phy(i,kde-2,j)
4555 !!! bug fix extrapolate ln(p) so p is positive definite
4556 p8w(i,kde,j) = exp(w1*log(p_phy(i,kde-1,j))+w2*log(p_phy(i,kde-2,j)))
4557 t8w(i,kde,j) = w1*t_phy(i,kde-1,j)+w2*t_phy(i,kde-2,j)
4559 enddo
4560 enddo
4562 ! decouple all physics tendencies
4564 IF (config_flags%ra_lw_physics .gt. 0 .or. config_flags%ra_sw_physics .gt. 0) THEN
4566 DO J=j_start,j_end
4567 DO K=k_start,k_end
4568 DO I=i_start,i_end
4569 RTHRATEN(I,K,J)=RTHRATEN(I,K,J)/mu(I,J)
4570 ENDDO
4571 ENDDO
4572 ENDDO
4574 ENDIF
4576 IF (config_flags%cu_physics .gt. 0) THEN
4578 DO J=j_start,j_end
4579 DO I=i_start,i_end
4580 DO K=k_start,k_end
4581 RTHCUTEN(I,K,J)=RTHCUTEN(I,K,J)/mu(I,J)
4582 ENDDO
4583 ENDDO
4584 ENDDO
4586 IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
4587 DO J=j_start,j_end
4588 DO I=i_start,i_end
4589 DO K=k_start,k_end
4590 RQVCUTEN(I,K,J)=RQVCUTEN(I,K,J)/mu(I,J)
4591 ENDDO
4592 ENDDO
4593 ENDDO
4594 ENDIF
4596 IF (P_QC .ge. PARAM_FIRST_SCALAR)THEN
4597 DO J=j_start,j_end
4598 DO I=i_start,i_end
4599 DO K=k_start,k_end
4600 RQCCUTEN(I,K,J)=RQCCUTEN(I,K,J)/mu(I,J)
4601 ENDDO
4602 ENDDO
4603 ENDDO
4604 ENDIF
4606 IF (P_QR .ge. PARAM_FIRST_SCALAR)THEN
4607 DO J=j_start,j_end
4608 DO I=i_start,i_end
4609 DO K=k_start,k_end
4610 RQRCUTEN(I,K,J)=RQRCUTEN(I,K,J)/mu(I,J)
4611 ENDDO
4612 ENDDO
4613 ENDDO
4614 ENDIF
4616 IF (P_QI .ge. PARAM_FIRST_SCALAR)THEN
4617 DO J=j_start,j_end
4618 DO I=i_start,i_end
4619 DO K=k_start,k_end
4620 RQICUTEN(I,K,J)=RQICUTEN(I,K,J)/mu(I,J)
4621 ENDDO
4622 ENDDO
4623 ENDDO
4624 ENDIF
4626 IF(P_QS .ge. PARAM_FIRST_SCALAR)THEN
4627 DO J=j_start,j_end
4628 DO I=i_start,i_end
4629 DO K=k_start,k_end
4630 RQSCUTEN(I,K,J)=RQSCUTEN(I,K,J)/mu(I,J)
4631 ENDDO
4632 ENDDO
4633 ENDDO
4634 ENDIF
4636 ENDIF
4638 IF (config_flags%bl_pbl_physics .gt. 0) THEN
4640 DO J=j_start,j_end
4641 DO K=k_start,k_end
4642 DO I=i_start,i_end
4643 RUBLTEN(I,K,J) =RUBLTEN(I,K,J)/mu(I,J)
4644 RVBLTEN(I,K,J) =RVBLTEN(I,K,J)/mu(I,J)
4645 RTHBLTEN(I,K,J)=RTHBLTEN(I,K,J)/mu(I,J)
4646 ENDDO
4647 ENDDO
4648 ENDDO
4650 IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
4651 DO J=j_start,j_end
4652 DO K=k_start,k_end
4653 DO I=i_start,i_end
4654 RQVBLTEN(I,K,J)=RQVBLTEN(I,K,J)/mu(I,J)
4655 ENDDO
4656 ENDDO
4657 ENDDO
4658 ENDIF
4660 IF (P_QC .ge. PARAM_FIRST_SCALAR) THEN
4661 DO J=j_start,j_end
4662 DO K=k_start,k_end
4663 DO I=i_start,i_end
4664 RQCBLTEN(I,K,J)=RQCBLTEN(I,K,J)/mu(I,J)
4665 ENDDO
4666 ENDDO
4667 ENDDO
4668 ENDIF
4670 IF (P_QI .ge. PARAM_FIRST_SCALAR) THEN
4671 DO J=j_start,j_end
4672 DO K=k_start,k_end
4673 DO I=i_start,i_end
4674 RQIBLTEN(I,K,J)=RQIBLTEN(I,K,J)/mu(I,J)
4675 ENDDO
4676 ENDDO
4677 ENDDO
4678 ENDIF
4680 ENDIF
4682 ! decouple advective forcing required by Grell-Devenyi scheme
4684 if(( config_flags%cu_physics == GDSCHEME ) .OR. &
4685 ( config_flags%cu_physics == G3SCHEME )) then
4687 DO J=j_start,j_end
4688 DO I=i_start,i_end
4689 DO K=k_start,k_end
4690 RTHFTEN(I,K,J)=RTHFTEN(I,K,J)/mu(I,J)
4691 ENDDO
4692 ENDDO
4693 ENDDO
4695 IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
4696 DO J=j_start,j_end
4697 DO I=i_start,i_end
4698 DO K=k_start,k_end
4699 RQVFTEN(I,K,J)=RQVFTEN(I,K,J)/mu(I,J)
4700 ENDDO
4701 ENDDO
4702 ENDDO
4703 ENDIF
4705 END IF
4707 ! fdda
4708 ! note fdda u and v tendencies are staggered, also only interior points have muu/muv,
4709 ! so only decouple those
4711 IF (config_flags%grid_fdda .gt. 0) THEN
4713 i_startu=MAX(its,ids+1)
4714 j_startv=MAX(jts,jds+1)
4716 DO J=j_start,j_end
4717 DO K=k_start,k_end
4718 DO I=i_startu,i_end
4719 RUNDGDTEN(I,K,J) =RUNDGDTEN(I,K,J)/muu(I,J)
4720 ENDDO
4721 ENDDO
4722 ENDDO
4723 DO J=j_startv,j_end
4724 DO K=k_start,k_end
4725 DO I=i_start,i_end
4726 RVNDGDTEN(I,K,J) =RVNDGDTEN(I,K,J)/muv(I,J)
4727 ENDDO
4728 ENDDO
4729 ENDDO
4730 DO J=j_start,j_end
4731 DO K=k_start,k_end
4732 DO I=i_start,i_end
4733 RTHNDGDTEN(I,K,J)=RTHNDGDTEN(I,K,J)/mu(I,J)
4734 ! RMUNDGDTEN(I,J) - no coupling
4735 ENDDO
4736 ENDDO
4737 ENDDO
4738 IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
4739 DO J=j_start,j_end
4740 DO K=k_start,k_end
4741 DO I=i_start,i_end
4742 RQVNDGDTEN(I,K,J)=RQVNDGDTEN(I,K,J)/mu(I,J)
4743 ENDDO
4744 ENDDO
4745 ENDDO
4746 ENDIF
4748 ENDIF
4750 END SUBROUTINE phy_prep
4752 !------------------------------------------------------------
4754 SUBROUTINE moist_physics_prep_em( t_new, t_old, t0, rho, al, alb, &
4755 p, p8w, p0, pb, ph, phb, &
4756 th_phy, pii, pf, &
4757 z, z_at_w, dz8w, &
4758 dt,h_diabatic, &
4759 config_flags,fzm, fzp, &
4760 ids,ide, jds,jde, kds,kde, &
4761 ims,ime, jms,jme, kms,kme, &
4762 its,ite, jts,jte, kts,kte )
4764 IMPLICIT NONE
4766 ! Here we construct full fields
4767 ! needed by the microphysics
4769 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
4771 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
4772 INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
4773 INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
4775 REAL, INTENT(IN ) :: dt
4777 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4778 INTENT(IN ) :: al, &
4779 alb, &
4780 p, &
4781 pb, &
4782 ph, &
4783 phb
4786 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
4787 fzp
4789 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4790 INTENT( OUT) :: rho, &
4791 th_phy, &
4792 pii, &
4793 pf, &
4794 z, &
4795 z_at_w, &
4796 dz8w, &
4797 p8w
4799 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4800 INTENT(INOUT) :: h_diabatic
4802 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4803 INTENT(INOUT) :: t_new, &
4804 t_old
4806 REAL, INTENT(IN ) :: t0, p0
4807 REAL :: z0,z1,z2,w1,w2
4809 INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
4810 INTEGER :: i, j, k
4812 !--------------------------------------------------------------------
4814 !<DESCRIPTION>
4815 !
4816 ! moist_phys_prep_em calculates a number of diagnostic quantities needed by
4817 ! the microphysics routines.
4818 !
4819 !</DESCRIPTION>
4821 ! set up loop bounds for this grid's boundary conditions
4823 i_start = its
4824 i_end = min( ite,ide-1 )
4825 j_start = jts
4826 j_end = min( jte,jde-1 )
4828 k_start = kts
4829 k_end = min( kte, kde-1 )
4831 DO j = j_start, j_end
4832 DO k = k_start, kte
4833 DO i = i_start, i_end
4834 z_at_w(i,k,j) = (ph(i,k,j)+phb(i,k,j))/g
4835 ENDDO
4836 ENDDO
4837 ENDDO
4839 do j = j_start,j_end
4840 do k = k_start, kte-1
4841 do i = i_start, i_end
4842 dz8w(i,k,j) = z_at_w(i,k+1,j)-z_at_w(i,k,j)
4843 enddo
4844 enddo
4845 enddo
4847 do j = j_start,j_end
4848 do i = i_start, i_end
4849 dz8w(i,kte,j) = 0.
4850 enddo
4851 enddo
4854 ! compute full pii, rho, and z at the new time-level
4855 ! (needed for physics).
4856 ! convert perturbation theta to full theta (th_phy)
4857 ! use h_diabatic to temporarily save pre-microphysics full theta
4859 DO j = j_start, j_end
4860 DO k = k_start, k_end
4861 DO i = i_start, i_end
4863 #ifdef REVERT
4864 t_new(i,k,j) = t_new(i,k,j)-h_diabatic(i,k,j)*dt
4865 #endif
4866 th_phy(i,k,j) = t_new(i,k,j) + t0
4867 h_diabatic(i,k,j) = th_phy(i,k,j)
4868 rho(i,k,j) = 1./(al(i,k,j)+alb(i,k,j))
4869 pii(i,k,j) = ((p(i,k,j)+pb(i,k,j))/p0)**rcp
4870 z(i,k,j) = 0.5*(z_at_w(i,k,j) +z_at_w(i,k+1,j) )
4871 pf(i,k,j) = p(i,k,j)+pb(i,k,j)
4873 ENDDO
4874 ENDDO
4875 ENDDO
4877 ! interp t and p at w points
4879 do j = j_start,j_end
4880 do k = 2, k_end
4881 do i = i_start, i_end
4882 p8w(i,k,j) = fzm(k)*pf(i,k,j)+fzp(k)*pf(i,k-1,j)
4883 enddo
4884 enddo
4885 enddo
4887 ! extrapolate p and t to surface and top.
4888 ! we'll use an extrapolation in z for now
4890 do j = j_start,j_end
4891 do i = i_start, i_end
4893 ! bottom
4895 z0 = z_at_w(i,1,j)
4896 z1 = z(i,1,j)
4897 z2 = z(i,2,j)
4898 w1 = (z0 - z2)/(z1 - z2)
4899 w2 = 1. - w1
4900 p8w(i,1,j) = w1*pf(i,1,j)+w2*pf(i,2,j)
4902 ! top
4904 z0 = z_at_w(i,kte,j)
4905 z1 = z(i,k_end,j)
4906 z2 = z(i,k_end-1,j)
4907 w1 = (z0 - z2)/(z1 - z2)
4908 w2 = 1. - w1
4909 ! p8w(i,kde,j) = w1*pf(i,kde-1,j)+w2*pf(i,kde-2,j)
4910 p8w(i,kde,j) = exp(w1*log(pf(i,kde-1,j))+w2*log(pf(i,kde-2,j)))
4912 enddo
4913 enddo
4915 END SUBROUTINE moist_physics_prep_em
4917 !------------------------------------------------------------------------------
4919 SUBROUTINE moist_physics_finish_em( t_new, t_old, t0, mut, &
4920 th_phy, h_diabatic, dt, &
4921 config_flags, &
4922 ids,ide, jds,jde, kds,kde, &
4923 ims,ime, jms,jme, kms,kme, &
4924 its,ite, jts,jte, kts,kte )
4926 IMPLICIT NONE
4928 ! Here we construct full fields
4929 ! needed by the microphysics
4931 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
4933 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
4934 INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
4935 INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
4937 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4938 INTENT(INOUT) :: t_new, &
4939 t_old, &
4940 th_phy, &
4941 h_diabatic
4943 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: mut
4946 REAL, INTENT(IN ) :: t0, dt
4948 INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
4949 INTEGER :: i, j, k
4951 !--------------------------------------------------------------------
4953 !<DESCRIPTION>
4954 !
4955 ! moist_phys_finish_em resets theta to its perturbation value and
4956 ! computes and stores the microphysics diabatic heating term.
4957 !
4958 !</DESCRIPTION>
4960 ! set up loop bounds for this grid's boundary conditions
4963 i_start = its
4964 i_end = min( ite,ide-1 )
4965 j_start = jts
4966 j_end = min( jte,jde-1 )
4968 k_start = kts
4969 k_end = min( kte, kde-1 )
4971 ! add microphysics theta diff to perturbation theta, set h_diabatic
4973 IF ( config_flags%no_mp_heating .eq. 0 ) THEN
4974 DO j = j_start, j_end
4975 DO k = k_start, k_end
4976 DO i = i_start, i_end
4977 t_new(i,k,j) = t_new(i,k,j) + (th_phy(i,k,j)-h_diabatic(i,k,j))
4978 h_diabatic(i,k,j) = (th_phy(i,k,j)-h_diabatic(i,k,j))/dt
4979 ENDDO
4980 ENDDO
4981 ENDDO
4983 ELSE
4985 DO j = j_start, j_end
4986 DO k = k_start, k_end
4987 DO i = i_start, i_end
4988 ! t_new(i,k,j) = t_new(i,k,j)
4989 h_diabatic(i,k,j) = 0.
4990 ENDDO
4991 ENDDO
4992 ENDDO
4993 ENDIF
4995 END SUBROUTINE moist_physics_finish_em
4997 !----------------------------------------------------------------
5000 SUBROUTINE init_module_big_step
5001 END SUBROUTINE init_module_big_step
5003 SUBROUTINE set_tend ( field, field_adv_tend, msf, &
5004 ids, ide, jds, jde, kds, kde, &
5005 ims, ime, jms, jme, kms, kme, &
5006 its, ite, jts, jte, kts, kte )
5008 IMPLICIT NONE
5010 ! Input data
5012 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
5013 ims, ime, jms, jme, kms, kme, &
5014 its, ite, jts, jte, kts, kte
5016 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: field
5018 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN) :: field_adv_tend
5020 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: msf
5022 ! Local data
5024 INTEGER :: i, j, k, itf, jtf, ktf
5026 !<DESCRIPTION>
5027 !
5028 ! set_tend copies the advective tendency array into the tendency array.
5029 !
5030 !</DESCRIPTION>
5032 jtf = MIN(jte,jde-1)
5033 ktf = MIN(kte,kde-1)
5034 itf = MIN(ite,ide-1)
5035 DO j = jts, jtf
5036 DO k = kts, ktf
5037 DO i = its, itf
5038 field(i,k,j) = field_adv_tend(i,k,j)*msf(i,j)
5039 ENDDO
5040 ENDDO
5041 ENDDO
5043 END SUBROUTINE set_tend
5045 !------------------------------------------------------------------------------
5047 SUBROUTINE rk_rayleigh_damp( ru_tendf, rv_tendf, &
5048 rw_tendf, t_tendf, &
5049 u, v, w, t, t_init, &
5050 mut, muu, muv, ph, phb, &
5051 u_base, v_base, t_base, z_base, &
5052 dampcoef, zdamp, &
5053 ids, ide, jds, jde, kds, kde, &
5054 ims, ime, jms, jme, kms, kme, &
5055 its, ite, jts, jte, kts, kte )
5057 ! History: Apr 2005 Modifications by George Bryan, NCAR:
5058 ! - Generalized the code in a way that allows for
5059 ! simulations with steep terrain.
5060 !
5061 ! Jul 2004 Modifications by George Bryan, NCAR:
5062 ! - Modified the code to use u_base, v_base, and t_base
5063 ! arrays for the background state. Removed the hard-wired
5064 ! base-state values.
5065 ! - Modified the code to use dampcoef, zdamp, and damp_opt,
5066 ! i.e., the upper-level damper variables in namelist.input.
5067 ! Removed the hard-wired variables in the older version.
5068 ! This damper is used when damp_opt = 2.
5069 ! - Modified the code to account for the movement of the
5070 ! model surfaces with time. The code now obtains a base-
5071 ! state value by interpolation using the "_base" arrays.
5073 ! Nov 2003 Bug fix by Jason Knievel, NCAR
5075 ! Aug 2003 Meridional dimension, some comments, and
5076 ! changes in layout of the code added by
5077 ! Jason Knievel, NCAR
5079 ! Jul 2003 Original code by Bill Skamarock, NCAR
5081 ! Purpose: This routine applies Rayleigh damping to a layer at top
5082 ! of the model domain.
5084 !-----------------------------------------------------------------------
5085 ! Begin declarations.
5087 IMPLICIT NONE
5089 INTEGER, INTENT( IN ) &
5090 :: ids, ide, jds, jde, kds, kde, &
5091 ims, ime, jms, jme, kms, kme, &
5092 its, ite, jts, jte, kts, kte
5094 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
5095 :: ru_tendf, rv_tendf, rw_tendf, t_tendf
5097 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
5098 :: u, v, w, t, t_init, ph, phb
5100 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
5101 :: mut, muu, muv
5103 REAL, DIMENSION( kms:kme ) , INTENT(IN ) &
5104 :: u_base, v_base, t_base, z_base
5106 REAL, INTENT(IN ) &
5107 :: dampcoef, zdamp
5109 ! Local variables.
5111 INTEGER &
5112 :: i_start, i_end, j_start, j_end, k_start, k_end, i, j, k, ktf, k1, k2
5114 REAL &
5115 :: pii, dcoef, z, ztop
5117 REAL :: wkp1, wk, wkm1
5119 REAL, DIMENSION( kms:kme ) :: z00, u00, v00, t00
5121 ! End declarations.
5122 !-----------------------------------------------------------------------
5124 pii = 2.0 * asin(1.0)
5126 ktf = MIN( kte, kde-1 )
5128 !-----------------------------------------------------------------------
5129 ! Adjust u to base state.
5131 DO j = jts, MIN( jte, jde-1 )
5132 DO i = its, MIN( ite, ide )
5134 ! Get height at top of model
5135 ztop = 0.5*( phb(i ,kde,j)+phb(i-1,kde,j) &
5136 +ph(i ,kde,j)+ph(i-1,kde,j) )/g
5138 ! Find bottom of damping layer
5139 k1 = ktf
5140 z = ztop
5141 DO WHILE( z >= (ztop-zdamp) )
5142 z = 0.25*( phb(i ,k1,j)+phb(i ,k1+1,j) &
5143 +phb(i-1,k1,j)+phb(i-1,k1+1,j) &
5144 +ph(i ,k1,j)+ph(i ,k1+1,j) &
5145 +ph(i-1,k1,j)+ph(i-1,k1+1,j))/g
5146 z00(k1) = z
5147 k1 = k1 - 1
5148 ENDDO
5149 k1 = k1 + 2
5151 ! Get reference state at model levels
5152 DO k = k1, ktf
5153 k2 = ktf
5154 DO WHILE( z_base(k2) .gt. z00(k) )
5155 k2 = k2 - 1
5156 ENDDO
5157 if(k2+1.gt.ktf)then
5158 u00(k) = u_base(k2) + ( u_base(k2) - u_base(k2-1) ) &
5159 * ( z00(k) - z_base(k2) ) &
5160 / ( z_base(k2) - z_base(k2-1) )
5161 else
5162 u00(k) = u_base(k2) + ( u_base(k2+1) - u_base(k2) ) &
5163 * ( z00(k) - z_base(k2) ) &
5164 / ( z_base(k2+1) - z_base(k2) )
5165 endif
5166 ENDDO
5168 ! Apply the Rayleigh damper
5169 DO k = k1, ktf
5170 dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
5171 dcoef = (SIN( 0.5 * pii * dcoef ) )**2
5172 ru_tendf(i,k,j) = ru_tendf(i,k,j) - &
5173 muu(i,j) * ( dcoef * dampcoef ) * &
5174 ( u(i,k,j) - u00(k) )
5175 END DO
5177 END DO
5178 END DO
5180 ! End adjustment of u.
5181 !-----------------------------------------------------------------------
5183 !-----------------------------------------------------------------------
5184 ! Adjust v to base state.
5186 DO j = jts, MIN( jte, jde )
5187 DO i = its, MIN( ite, ide-1 )
5189 ! Get height at top of model
5190 ztop = 0.5*( phb(i,kde,j )+phb(i,kde,j-1) &
5191 +ph(i,kde,j )+ph(i,kde,j-1) )/g
5193 ! Find bottom of damping layer
5194 k1 = ktf
5195 z = ztop
5196 DO WHILE( z >= (ztop-zdamp) )
5197 z = 0.25*( phb(i,k1,j )+phb(i,k1+1,j ) &
5198 +phb(i,k1,j-1)+phb(i,k1+1,j-1) &
5199 +ph(i,k1,j )+ph(i,k1+1,j ) &
5200 +ph(i,k1,j-1)+ph(i,k1+1,j-1))/g
5201 z00(k1) = z
5202 k1 = k1 - 1
5203 ENDDO
5204 k1 = k1 + 2
5206 ! Get reference state at model levels
5207 DO k = k1, ktf
5208 k2 = ktf
5209 DO WHILE( z_base(k2) .gt. z00(k) )
5210 k2 = k2 - 1
5211 ENDDO
5212 if(k2+1.gt.ktf)then
5213 v00(k) = v_base(k2) + ( v_base(k2) - v_base(k2-1) ) &
5214 * ( z00(k) - z_base(k2) ) &
5215 / ( z_base(k2) - z_base(k2-1) )
5216 else
5217 v00(k) = v_base(k2) + ( v_base(k2+1) - v_base(k2) ) &
5218 * ( z00(k) - z_base(k2) ) &
5219 / ( z_base(k2+1) - z_base(k2) )
5220 endif
5221 ENDDO
5223 ! Apply the Rayleigh damper
5224 DO k = k1, ktf
5225 dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
5226 dcoef = (SIN( 0.5 * pii * dcoef ) )**2
5227 rv_tendf(i,k,j) = rv_tendf(i,k,j) - &
5228 muv(i,j) * ( dcoef * dampcoef ) * &
5229 ( v(i,k,j) - v00(k) )
5230 END DO
5232 END DO
5233 END DO
5235 ! End adjustment of v.
5236 !-----------------------------------------------------------------------
5238 !-----------------------------------------------------------------------
5239 ! Adjust w to base state.
5241 DO j = jts, MIN( jte, jde-1 )
5242 DO i = its, MIN( ite, ide-1 )
5243 ztop = ( phb(i,kde,j) + ph(i,kde,j) ) / g
5244 DO k = kts, MIN( kte, kde )
5245 z = ( phb(i,k,j) + ph(i,k,j) ) / g
5246 IF ( z >= (ztop-zdamp) ) THEN
5247 dcoef = 1.0 - MIN( 1.0, ( ztop - z ) / zdamp )
5248 dcoef = ( SIN( 0.5 * pii * dcoef ) )**2
5249 rw_tendf(i,k,j) = rw_tendf(i,k,j) - &
5250 mut(i,j) * ( dcoef * dampcoef ) * w(i,k,j)
5251 END IF
5252 END DO
5253 END DO
5254 END DO
5256 ! End adjustment of w.
5257 !-----------------------------------------------------------------------
5259 !-----------------------------------------------------------------------
5260 ! Adjust potential temperature to base state.
5262 DO j = jts, MIN( jte, jde-1 )
5263 DO i = its, MIN( ite, ide-1 )
5265 ! Get height at top of model
5266 ztop = ( phb(i,kde,j) + ph(i,kde,j) ) / g
5268 ! Find bottom of damping layer
5269 k1 = ktf
5270 z = ztop
5271 DO WHILE( z >= (ztop-zdamp) )
5272 z = 0.5 * ( phb(i,k1,j) + phb(i,k1+1,j) + &
5273 ph(i,k1,j) + ph(i,k1+1,j) ) / g
5274 z00(k1) = z
5275 k1 = k1 - 1
5276 ENDDO
5277 k1 = k1 + 2
5279 ! Get reference state at model levels
5280 DO k = k1, ktf
5281 k2 = ktf
5282 DO WHILE( z_base(k2) .gt. z00(k) )
5283 k2 = k2 - 1
5284 ENDDO
5285 if(k2+1.gt.ktf)then
5286 t00(k) = t_base(k2) + ( t_base(k2) - t_base(k2-1) ) &
5287 * ( z00(k) - z_base(k2) ) &
5288 / ( z_base(k2) - z_base(k2-1) )
5289 else
5290 t00(k) = t_base(k2) + ( t_base(k2+1) - t_base(k2) ) &
5291 * ( z00(k) - z_base(k2) ) &
5292 / ( z_base(k2+1) - z_base(k2) )
5293 endif
5294 ENDDO
5296 ! Apply the Rayleigh damper
5297 DO k = k1, ktf
5298 dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
5299 dcoef = (SIN( 0.5 * pii * dcoef ) )**2
5300 t_tendf(i,k,j) = t_tendf(i,k,j) - &
5301 mut(i,j) * ( dcoef * dampcoef ) * &
5302 ( t(i,k,j) - t00(k) )
5303 END DO
5305 END DO
5306 END DO
5308 ! End adjustment of potential temperature.
5309 !-----------------------------------------------------------------------
5311 END SUBROUTINE rk_rayleigh_damp
5313 !==============================================================================
5314 !==============================================================================
5316 SUBROUTINE sixth_order_diffusion( name, field, tendency, mu, dt, &
5317 config_flags, &
5318 diff_6th_opt, diff_6th_factor, &
5319 ids, ide, jds, jde, kds, kde, &
5320 ims, ime, jms, jme, kms, kme, &
5321 its, ite, jts, jte, kts, kte )
5323 ! History: 14 Nov 2006 Name of variable changed by Jason Knievel
5324 ! 07 Jun 2006 Revised and generalized by Jason Knievel
5325 ! 25 Apr 2005 Original code by Jason Knievel, NCAR
5327 ! Purpose: Apply 6th-order, monotonic (flux-limited), numerical
5328 ! diffusion to 3-d velocity and to scalars.
5330 ! References: Ming Xue (MWR Aug 2000)
5331 ! Durran ("Numerical Methods for Wave Equations..." 1999)
5332 ! George Bryan (personal communication)
5334 !------------------------------------------------------------------------------
5335 ! Begin: Declarations.
5337 IMPLICIT NONE
5339 INTEGER, INTENT(IN) &
5340 :: ids, ide, jds, jde, kds, kde, &
5341 ims, ime, jms, jme, kms, kme, &
5342 its, ite, jts, jte, kts, kte
5344 TYPE(grid_config_rec_type), INTENT(IN) &
5345 :: config_flags
5347 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) &
5348 :: tendency
5350 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) &
5351 :: field
5353 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN) &
5354 :: mu
5356 REAL, INTENT(IN) &
5357 :: dt
5359 REAL, INTENT(IN) &
5360 :: diff_6th_factor
5362 INTEGER, INTENT(IN) &
5363 :: diff_6th_opt
5365 CHARACTER(LEN=1) , INTENT(IN) &
5366 :: name
5368 INTEGER &
5369 :: i, j, k, &
5370 i_start, i_end, &
5371 j_start, j_end, &
5372 k_start, k_end, &
5373 ktf
5375 REAL &
5376 :: dflux_x_p0, dflux_y_p0, &
5377 dflux_x_p1, dflux_y_p1, &
5378 tendency_x, tendency_y, &
5379 mu_avg_p0, mu_avg_p1, &
5380 diff_6th_coef
5382 LOGICAL &
5383 :: specified
5385 ! End: Declarations.
5386 !------------------------------------------------------------------------------
5388 !------------------------------------------------------------------------------
5389 ! Begin: Translate the diffusion factor into a diffusion coefficient. See
5390 ! Durran's text, section 2.4.3, then adjust for sixth-order diffusion (not
5391 ! fourth) and for diffusion in two dimensions (not one). For reference, a
5392 ! factor of 1.0 would mean complete diffusion of a 2dx wave in one time step,
5393 ! although application of the flux limiter reduces somewhat the effects of
5394 ! diffusion for a given coefficient.
5396 diff_6th_coef = diff_6th_factor * 0.015625 / ( 2.0 * dt )
5398 ! End: Translate diffusion factor.
5399 !------------------------------------------------------------------------------
5401 !------------------------------------------------------------------------------
5402 ! Begin: Assign limits of spatial loops depending on variable to be diffused.
5403 ! The halo regions are already filled with values by the time this subroutine
5404 ! is called, which allows the stencil to extend beyond the domains' edges.
5406 ktf = MIN( kte, kde-1 )
5408 IF ( name .EQ. 'u' ) THEN
5410 i_start = its
5411 i_end = ite
5412 j_start = jts
5413 j_end = MIN(jde-1,jte)
5414 k_start = kts
5415 k_end = ktf
5417 ELSE IF ( name .EQ. 'v' ) THEN
5419 i_start = its
5420 i_end = MIN(ide-1,ite)
5421 j_start = jts
5422 j_end = jte
5423 k_start = kts
5424 k_end = ktf
5426 ELSE IF ( name .EQ. 'w' ) THEN
5428 i_start = its
5429 i_end = MIN(ide-1,ite)
5430 j_start = jts
5431 j_end = MIN(jde-1,jte)
5432 k_start = kts+1
5433 k_end = ktf
5435 ELSE
5437 i_start = its
5438 i_end = MIN(ide-1,ite)
5439 j_start = jts
5440 j_end = MIN(jde-1,jte)
5441 k_start = kts
5442 k_end = ktf
5444 ENDIF
5446 ! End: Assignment of limits of spatial loops.
5447 !------------------------------------------------------------------------------
5449 !------------------------------------------------------------------------------
5450 ! Begin: Loop across spatial dimensions.
5452 DO j = j_start, j_end
5453 DO k = k_start, k_end
5454 DO i = i_start, i_end
5456 !------------------------------------------------------------------------------
5457 ! Begin: Diffusion in x (i index).
5459 ! Calculate the diffusive flux in x direction (from Xue's eq. 3).
5461 dflux_x_p0 = ( 10.0 * ( field(i, k,j) - field(i-1,k,j) ) &
5462 - 5.0 * ( field(i+1,k,j) - field(i-2,k,j) ) &
5463 + ( field(i+2,k,j) - field(i-3,k,j) ) )
5465 dflux_x_p1 = ( 10.0 * ( field(i+1,k,j) - field(i ,k,j) ) &
5466 - 5.0 * ( field(i+2,k,j) - field(i-1,k,j) ) &
5467 + ( field(i+3,k,j) - field(i-2,k,j) ) )
5469 ! If requested in the namelist (diff_6th_opt=2), prohibit up-gradient diffusion
5470 ! (variation on Xue's eq. 10).
5472 IF ( diff_6th_opt .EQ. 2 ) THEN
5474 IF ( dflux_x_p0 * ( field(i ,k,j)-field(i-1,k,j) ) .LE. 0.0 ) THEN
5475 dflux_x_p0 = 0.0
5476 END IF
5478 IF ( dflux_x_p1 * ( field(i+1,k,j)-field(i ,k,j) ) .LE. 0.0 ) THEN
5479 dflux_x_p1 = 0.0
5480 END IF
5482 END IF
5484 ! Apply 6th-order diffusion in x direction.
5486 IF ( name .EQ. 'u' ) THEN
5487 mu_avg_p0 = mu(i-1,j)
5488 mu_avg_p1 = mu(i ,j)
5489 ELSE IF ( name .EQ. 'v' ) THEN
5490 mu_avg_p0 = 0.25 * ( &
5491 mu(i-1,j-1) + &
5492 mu(i ,j-1) + &
5493 mu(i-1,j ) + &
5494 mu(i ,j ) )
5495 mu_avg_p1 = 0.25 * ( &
5496 mu(i ,j-1) + &
5497 mu(i+1,j-1) + &
5498 mu(i ,j ) + &
5499 mu(i+1,j ) )
5500 ELSE
5501 mu_avg_p0 = 0.5 * ( &
5502 mu(i-1,j) + &
5503 mu(i ,j) )
5504 mu_avg_p1 = 0.5 * ( &
5505 mu(i ,j) + &
5506 mu(i+1,j) )
5507 END IF
5509 tendency_x = diff_6th_coef * &
5510 ( ( mu_avg_p1 * dflux_x_p1 ) - ( mu_avg_p0 * dflux_x_p0 ) )
5512 ! End: Diffusion in x.
5513 !------------------------------------------------------------------------------
5515 !------------------------------------------------------------------------------
5516 ! Begin: Diffusion in y (j index).
5518 ! Calculate the diffusive flux in y direction (from Xue's eq. 3).
5520 dflux_y_p0 = ( 10.0 * ( field(i,k,j ) - field(i,k,j-1) ) &
5521 - 5.0 * ( field(i,k,j+1) - field(i,k,j-2) ) &
5522 + ( field(i,k,j+2) - field(i,k,j-3) ) )
5524 dflux_y_p1 = ( 10.0 * ( field(i,k,j+1) - field(i,k,j ) ) &
5525 - 5.0 * ( field(i,k,j+2) - field(i,k,j-1) ) &
5526 + ( field(i,k,j+3) - field(i,k,j-2) ) )
5528 ! If requested in the namelist (diff_6th_opt=2), prohibit up-gradient diffusion
5529 ! (variation on Xue's eq. 10).
5531 IF ( diff_6th_opt .EQ. 2 ) THEN
5533 IF ( dflux_y_p0 * ( field(i,k,j )-field(i,k,j-1) ) .LE. 0.0 ) THEN
5534 dflux_y_p0 = 0.0
5535 END IF
5537 IF ( dflux_y_p1 * ( field(i,k,j+1)-field(i,k,j ) ) .LE. 0.0 ) THEN
5538 dflux_y_p1 = 0.0
5539 END IF
5541 END IF
5543 ! Apply 6th-order diffusion in y direction.
5545 IF ( name .EQ. 'u' ) THEN
5546 mu_avg_p0 = 0.25 * ( &
5547 mu(i-1,j-1) + &
5548 mu(i ,j-1) + &
5549 mu(i-1,j ) + &
5550 mu(i ,j ) )
5551 mu_avg_p1 = 0.25 * ( &
5552 mu(i-1,j ) + &
5553 mu(i ,j ) + &
5554 mu(i-1,j+1) + &
5555 mu(i ,j+1) )
5556 ELSE IF ( name .EQ. 'v' ) THEN
5557 mu_avg_p0 = mu(i,j-1)
5558 mu_avg_p1 = mu(i,j )
5559 ELSE
5560 mu_avg_p0 = 0.5 * ( &
5561 mu(i,j-1) + &
5562 mu(i,j ) )
5563 mu_avg_p1 = 0.5 * ( &
5564 mu(i,j ) + &
5565 mu(i,j+1) )
5566 END IF
5568 tendency_y = diff_6th_coef * &
5569 ( ( mu_avg_p1 * dflux_y_p1 ) - ( mu_avg_p0 * dflux_y_p0 ) )
5571 ! End: Diffusion in y.
5572 !------------------------------------------------------------------------------
5574 !------------------------------------------------------------------------------
5575 ! Begin: Combine diffusion in x and y.
5577 tendency(i,k,j) = tendency(i,k,j) + tendency_x + tendency_y
5579 ! End: Combine diffusion in x and y.
5580 !------------------------------------------------------------------------------
5582 ENDDO
5583 ENDDO
5584 ENDDO
5586 ! End: Loop across spatial dimensions.
5587 !------------------------------------------------------------------------------
5589 END SUBROUTINE sixth_order_diffusion
5591 !==============================================================================
5592 !==============================================================================
5594 END MODULE module_big_step_utilities_em